From don@cavalieri.gsfc.nasa.gov Fri Jun 8 16:28:46 2001 Return-Path: Received: from mx1.ncep.noaa.gov (mx.wwb.noaa.gov [140.90.193.27]) by polar.ncep.noaa.gov (8.8.8/8.8.8) with ESMTP id QAA04250 for ; Fri, 8 Jun 2001 16:28:45 -0400 Received: from cavalieri.gsfc.nasa.gov (cavalieri.gsfc.nasa.gov [128.183.106.122]) by mx1.ncep.noaa.gov (8.8.8/8.8.8) with ESMTP id QAA20160 for ; Fri, 8 Jun 2001 16:21:20 -0400 Received: from cavalieri.gsfc.nasa.gov (localhost [127.0.0.1]) by cavalieri.gsfc.nasa.gov (980427.SGI.8.8.8/980728.SGI.AUTOCF) via ESMTP id QAA85528 for ; Fri, 8 Jun 2001 16:29:07 -0400 (EDT) Sender: don@cavalieri.gsfc.nasa.gov Message-ID: <3B213592.2FA07FA1@cavalieri.gsfc.nasa.gov> Date: Fri, 08 Jun 2001 16:29:06 -0400 From: "Donald J. Cavalieri" X-Mailer: Mozilla 4.51C-SGI [en] (X11; I; IRIX64 6.5 IP30) X-Accept-Language: en MIME-Version: 1.0 To: "Grumbine, Robert W." Subject: [Fwd: [Fwd: Re: NT2 code]] Content-Type: multipart/mixed; boundary="------------3F5AA8459436A8835CA1FAA2" Status: R This is a multi-part message in MIME format. --------------3F5AA8459436A8835CA1FAA2 Content-Type: multipart/alternative; boundary="------------21030D878ED1A6D7C0CE4137" --------------21030D878ED1A6D7C0CE4137 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Bob, Here is the NT2 code with some documentation from Thorsten. Please give me or Thorsten a call if you have any questions. Also, I put a reprint of the paper in the mail to you this morning. Have a good weekend! Don -- Donald J. Cavalieri, Ph.D. Code 971/NASA Goddard Space Flight Center Greenbelt, Maryland 20771-0001 Phone: (301) 614-5901 FAX: (301) 614-5644 --------------21030D878ED1A6D7C0CE4137 Content-Type: text/html; charset=us-ascii Content-Transfer-Encoding: 7bit Bob,
Here is the NT2 code with some documentation from Thorsten.  Please give me or Thorsten a call
if you have any questions.  Also, I put a reprint of the paper in the mail to you this morning.
Have a good weekend!
Don
-- 
Donald J. Cavalieri, Ph.D.
Code 971/NASA Goddard Space Flight Center
Greenbelt, Maryland 20771-0001
Phone:  (301) 614-5901     FAX: (301) 614-5644
  --------------21030D878ED1A6D7C0CE4137-- --------------3F5AA8459436A8835CA1FAA2 Content-Type: message/rfc822 Content-Transfer-Encoding: 7bit Content-Disposition: inline Received: from beaufort.gsfc.nasa.gov (beaufort.gsfc.nasa.gov [128.183.105.29]) by cavalieri.gsfc.nasa.gov (980427.SGI.8.8.8/980728.SGI.AUTOCF) via ESMTP id PAA87084 for ; Fri, 8 Jun 2001 15:35:36 -0400 (EDT) Received: from beaufort.gsfc.nasa.gov (localhost [127.0.0.1]) by beaufort.gsfc.nasa.gov (980427.SGI.8.8.8/980728.SGI.AUTOCF) via ESMTP id PAA58238 for ; Fri, 8 Jun 2001 15:27:58 -0400 (EDT) Sender: thorsten@beaufort Message-ID: <3B21273D.B2897B94@beaufort.gsfc.nasa.gov> Date: Fri, 08 Jun 2001 15:27:57 -0400 From: Thorsten Markus X-Mailer: Mozilla 4.7C-SGI [en] (X11; I; IRIX64 6.5 IP30) X-Accept-Language: en MIME-Version: 1.0 To: don Subject: [Fwd: Re: NT2 code] Content-Type: multipart/mixed; boundary="------------FDA98DD8FBAA2EA5A32FA467" X-Mozilla-Status2: 00000000 This is a multi-part message in MIME format. --------------FDA98DD8FBAA2EA5A32FA467 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Don, here's the NT2 code and all the other files needed Thorsten ------------------------------ I attach the C-file together with the atmospheric data tables. So you should receive a total of 10 files: nt2.c 7 files *.tab 2 landmask files (I know you don't need them, they simply have the name as used in the program) The c-code is written to read the NSIDC CD-data, i.e full hemisphere with the 294 byte hdf-header (I know that is not the way hdf data are supposed to be read). Althgough it should be rather simple to change it for different data formats. Basically, all that has to be changed would be the xsize and ysize in the beginning of the main-routine, and the subroutines get_tbs and get_msk. To compile it type cc nt2.c -n32 -lm -o nt2.out and then something like nt2.out 990307 n where here 990307 is the cddate and n the hemisphere. The output is, in this case: n990307.nt2 which is a short integer array (just intarr for IDL) Please note that in the subroutine get_tbs no paths are set, so the program assumes that the CD data are your hard drive. --------------FDA98DD8FBAA2EA5A32FA467 Content-Type: text/plain; charset=us-ascii; name="nt2.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="nt2.c" #include #include #include #define n_atm 12 /*** Number of atmospheres ***/ /******************************************* Initialize short int (2byte) matrix *******************************************/ short int **imatrix(int nrl,int nrh,int ncl,int nch) { short int **m; int i; m=(short int **)calloc((nrh-nrl+1),sizeof(short int*)); if (!m) printf("Error in malloc 1"); m -= nrl; for(i=nrl;i<=nrh;i++){ m[i]=(short int *)calloc((nch-ncl+1),sizeof(short int)); if (!m[i]) printf("Error in malloc 2"); m[i]-=ncl; } return m; } /******************************************* Initialize float (4byte) matrix *******************************************/ float **fmatrix(int nrl,int nrh,int ncl,int nch) { float **m; int i; m=(float **)calloc((nrh-nrl+1),sizeof(float*)); if (!m) printf("Error in malloc 1"); m -= nrl; for(i=nrl;i<=nrh;i++){ m[i]=(float *)calloc((nch-ncl+1),sizeof(float)); if (!m[i]) printf("Error in malloc 2"); m[i]-=ncl; } return m; } /******************************************* Initialize float 3D matrix *******************************************/ float ***f3dmatrix(int nrl,int nrh,int ncl,int nch,int nzl,int nzh) { float ***m; int i,j; m=(float ***)calloc((nzh-nzl+1),sizeof(float**)); if (!m) printf("Error in malloc 1"); m -= nzl; for(j=nzl;j<=nzh;j++){ m[j]=(float **)calloc((nrh-nrl+1),sizeof(float*)); if (!m[j]) printf("Error in malloc 2"); for(i=nrl;i<=nrh;i++){ m[j][i]=(float *)calloc((nch-ncl+1),sizeof(float)); if (!m[j][i]) printf("Error in malloc 2"); m[j][i]-=ncl; } m[j]-=nrl; } return m; } /****************************************************** Read land mask ******************************************************/ void get_msk(short **msk, int xs, int ys, char *hemi) { int j; FILE *file; char hdr[300]; if (*hemi == 'n') file = fopen("landmask_north_25", "rb"); else file = fopen("landmask_south_25", "rb"); fread(hdr, sizeof(char), 300, file); for(j=0;j 500)&&(v85[y][x] > 500)&&(msk[y][x] == 0)){ /*** If data are valid and no land ***/ v19i=(double) v19[y][x]/10.; h19i=(double) h19[y][x]/10.; v22i=(double) v22[y][x]/10.; v37i=(double) v37[y][x]/10.; v85i=(double) v85[y][x]/10.; h85i=(double) h85[y][x]/10.; gr3719=(v37i-v19i)/(v37i+v19i); gr2219=(v22i-v19i)/(v22i+v19i); if ((gr3719 < 0.05) && (gr2219 < 0.045)){ /*** if passed the weather filters ***/ pr19=(v19i-h19i)/(v19i+h19i); pr85=(v85i-h85i)/(v85i+h85i); gr8519v=(v85i-v19i)/(v85i+v19i); gr8519h=(h85i-h19i)/(h85i+h19i); pr19r=-gr3719*sinphi19 + pr19*cosphi19; pr85r=-gr3719*sinphi85 + pr85*cosphi85; dgr=gr8519h-gr8519v; dmin=10000.; for (k=0;k=0)&&(ccj >=0)&& ((cai+ccj) >=0)&&((cai+ccj) < 101)){ if((*hemi == 'n')&&(gr3719 > -0.01)&& (((x > 0)&&(x < 80)&&(y >100)&&(y <190))|| ((x > 80)&&(x <160)&&(y >0)&&(y <130)))){ /** Only pixels of Bering Sea (top) and Sea of Okhotsk (bottom **/ dpr19=pr19r-LUT19thin[k][cai][ccj]; dpr85=pr85r-LUT85thin[k][cai][ccj]; ddgr=gr3719-LUTGR37[k][cai][ccj]; } else { dpr19=pr19r-LUT19[k][cai][ccj]; dpr85=pr85r-LUT85[k][cai][ccj]; ddgr=dgr-LUTDGR[k][cai][ccj]; } d=w19*dpr19*dpr19+w85*dpr85*dpr85+wgr*ddgr*ddgr; if (d < dmin){ dmin=d; imin=i; jmin=j; } } /*endif*/ } /*endfor*/ }/*endfor*/ ca=ca+imin; cc=cc+jmin; } /*do */ while ((imin !=0)||(jmin != 0)); camina[k]=ca; ccmina[k]=cc; dmina[k]=dmin; } /*k*/ bestk=20; dmin=1000; for (k=0;k