/* =============================================== */ /* These include files should be standard on all */ /* UNIX systems. */ /* =============================================== */ #if defined(__convex__) #include #else #include #endif #include #include #include #include #include #include #include #include #include #include #ifdef _CRAY #include /* prototype for cray fortran externals */ fortran void CONVERT_IEEE(_fcd, int *, int *, int *, int *); _fcd chf771, chf772; #endif /* =============================================== */ /* int should be at least 32 bits long */ /* =============================================== */ #ifndef UINT #define UINT unsigned long long #endif #ifndef INT #define INT long #endif #ifndef REAL #define REAL float #endif #ifndef CHAR #define CHAR unsigned char #endif #ifndef SEEK_CUR #define SEEK_CUR 1 #endif #ifndef CLOCKS_PER_SEC #define CLOCKS_PER_SEC 1000000. #endif #define TRUE 1 #define FALSE 0 #ifndef _CRAY #endif #define LEV_SURFACE 1 #define LEV_MEANSEA 102 #define LEV_OCEAN 160 #define REP_REGULAR 0 #define REP_GAUSS 4 #define REP_SPECTRAL 50 #define MAX_Garbage 65535 #define MAX_Dates 65535 #define MAX_Records 65535 #define MAX_Codes 256 #define MAX_Months 12 #define MAX_Levels 100 #define MAX_Lons 320 /* Number of ECHAM T106 Lons */ #define MAX_Lats 242 /* Number of HOPE T042 Lats */ #define MAX_DimGP (MAX_Lons * MAX_Lats) #define M_PI 3.14159265358979323846 #define fpr1(fpr,a) fprintf(fpr,a) #define fpr2(fpr,a,b) fprintf(fpr,a,b) #define fpr3(fpr,a,b,c) fprintf(fpr,a,b,c) #define fpr4(fpr,a,b,c,d) fprintf(fpr,a,b,c,d) #define fpr5(fpr,a,b,c,d,e) fprintf(fpr,a,b,c,d,e) char DatHori[128]; /* GRIB input file name */ char DatPref[128]; /* GRIB input file name prefix */ char DatDate[128]; /* GRIB input file name date */ char DatSuff[128]; /* GRIB input file name suffix */ char DatHNew[128]; /* GRIB output file name */ char DefHori[128]; /* GrADS control file name */ char VariDef[128]; /* Model specific code-name-explanation list */ char ExpID[5]; /* Data Experiment ID */ char Model[5]; /* DKRZ Model name mnemonic */ /* -------------------------------- */ /* Variables related to GRIB record */ /* -------------------------------- */ CHAR Block4Info[8]; const char Release[] = "RELEASE 3.1"; INT EastWest ; INT SouthNorth ; INT TopDown ; INT GribPos ;/* No of bytes between '7777' and 'GRIB' */ INT Grib1Offset ;/* Offset in bytes before and after block 1 */ INT GribEdition ;/* grib[ 3+Offset] 1 Byte */ INT CenterID ;/* grib[ 8+Offset] Block1[ 4] 1 Byte */ INT ModelID ;/* grib[ 9+Offset] Block1[ 5] 1 Byte */ INT CodeGrib ;/* grib[12+Offset] Block1[ 8] 1 Byte only */ INT CodeCount = 0; INT TypeOfLevel ;/* grib[13+Offset] Block1[ 9] 1 Byte */ INT LevelGrib ;/* grib[14+Offset] Block1[10] 1 or 2 Bytes */ INT Average ;/* grib[26+Offset] Block1[21] 2 Bytes */ INT RepType ; INT Block1Length; INT Block2Length; INT Block3Length; INT Block4Length; INT GribLength ; REAL FMin ; REAL DScale ; REAL ZScale ; INT Verbose = 0; /* Verbose mode specifier */ INT DimGP ; /* Horizontal dimension of grid point field */ INT FileCount = 0; /* Number of data files opened */ INT EndOfProcess ; /* Flag to indicate clean up */ INT Lats ; /* Number of latitudes found in grid point data */ INT Levels = 0; /* Number of levels found in grid point data */ INT Lons ; /* Number of longitudes found in grid point data */ INT IsMask = FALSE; /* Indicates block3 included */ INT FillUp = FALSE; /* Interpolation */ INT MulFil = 0; /* Time increment of template files */ INT Process = FALSE; /* Flag for readgrib to write out */ INT TermCount = 0; /* Number of different dates encountered */ INT UseInputFileDat = 0; INT GribCount = 0; /* Number of GRIB records processed */ INT PressureLevelData = FALSE; REAL Undef = -9.99E+33 ; INT Zero = 0; REAL FirstXPoint ; REAL FirstYPoint ; REAL LastXPoint ; REAL LastYPoint ; REAL XIncrement ; REAL YIncrement ; char *MoName[13] = {"DEC","Jan","Feb","Mar","Apr","May","Jun", "Jul","Aug","Sep","Oct","Nov","Dec"}; struct gaindx { /* struct to hold GRIB->GrADS conversion info */ int type; /* Indexing file type */ int hinum; /* Number of ints in header */ int hfnum; /* Number of floats in header */ int intnum; /* Number of index ints (long) */ int fltnum; /* Number of index floats */ int *hipnt; /* Pointer to header int values */ float *hfpnt; /* Point er to header float values */ int *intpnt; /* Pointer to int index values */ float *fltpnt; /* Pointer to float index values */ }; struct Date { /* struct to hold date/time info */ INT Year; INT Month; INT Day; INT Hour; INT Minute; UINT YMDH; /* date/time group */ }; struct Date NewDate; struct Date OldDate; struct Date AllDate[MAX_Dates]; struct Control { /* struct to hold GRIB record/field info */ char name[10]; char unit[10]; char expl[60]; INT code ; INT isvec ; INT olev ; INT ltyp ; }; struct Control All[MAX_Codes]; struct Control Uni[MAX_Codes]; CHAR *grib ; REAL *GrADSLevs; REAL *GrADSLats; FILE *intape; FILE *outape; FILE *vd ; /* ======================================= */ /* fields and scalars provided for gribmap */ /* ======================================= */ UINT vars[MAX_Records]; UINT levs[MAX_Records]; UINT tims[MAX_Records]; int dpos[MAX_Records]; int bpos[MAX_Records]; int bits[MAX_Records]; float dsca[MAX_Records]; float zsca[MAX_Records]; float fmin[MAX_Records]; int fields = 0 ; int dates = 0 ; /* =============================================================== */ /* level is the predefined vertical resolution of ocean model data */ /* =============================================================== */ INT *DefLevs; /* ====================== */ /* Define grid properties */ /* ====================== */ REAL *XLons,*YLats; double bessel[160] = { 2.4048255577, 5.5200781103, 8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679, 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684, 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132, 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550, 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299, 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711, 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819, 96.6052679510, 99.7468198587, 102.8883742542, 106.0299309165, 109.1714896498, 112.3130502805, 115.4546126537, 118.5961766309, 121.7377420880, 124.8793089132, 128.0208770059, 131.1624462752, 134.3040166383, 137.4455880203, 140.5871603528, 143.7287335737, 146.8703076258, 150.0118824570, 153.1534580192, 156.2950342685, 159.4366149902, 162.5782012939, 165.7197875977, 168.8613739014, 172.0029602051, 175.1445465088, 178.2861328125, 181.4277191162, 184.5693054199, 187.7108917236, 190.8524780273, 193.9940643311, 197.1356506348, 200.2772369385, 203.4188232422, 206.5604095459, 209.7019958496, 212.8435821533, 215.9851684570, 219.1267547607, 222.2683410645, 225.4099273682, 228.5515136719, 231.6930999756, 234.8346862793, 237.9762725830, 241.1178588867, 244.2594451904, 247.4010314941, 250.5426177979, 253.6842041016, 256.8258056641, 259.9674072266, 263.1090087891, 266.2506103516, 269.3922119141, 272.5338134766, 275.6754150391, 278.8170166016, 281.9586181641, 285.1002197266, 288.2418212891, 291.3834228516, 294.5250244141, 297.6666259766, 300.8082275391, 303.9498291016, 307.0914306641, 310.2330322266, 313.3746337891, 316.5162353516, 319.6578369141, 322.7994384766, 325.9410400391, 329.0826416016, 332.2242431641, 335.3658447266, 338.5074462891, 341.6490478516, 344.7906494141, 347.9322509766, 351.0738525391, 354.2154541016, 357.3570556641, 360.4986572266, 363.6402587891, 366.7818603516, 369.9234619141, 373.0650634766, 376.2066650391, 379.3482666016, 382.4898681641, 385.6314697266, 388.7730712891, 391.9146728516, 395.0562744141, 398.1978759766, 401.3394775391, 404.4810791016, 407.6226806641, 410.7642822266, 413.9058837891, 417.0474853516, 420.1890869141, 423.3306884766, 426.4722900391, 429.6138916016, 432.7554931641, 435.8970947266, 439.0386962891, 442.1802978516, 445.3218994141, 448.4635009766, 451.6051025391, 454.7467041016, 457.8883056641, 461.0299072266, 464.1715087891, 467.3131103516, 470.4547119141, 473.5963134766, 476.7379150391, 479.8795166016, 483.0211181641, 486.1627197266, 489.3043212891, 492.4459228516, 495.5875244141, 498.7291259766, 501.8707275391 }; /* ==================================== */ /* Abort - Print error message and exit */ /* ==================================== */ void Abort(char *errtext) { printf(errtext); exit(1); } INT Get2Byte(CHAR ptr[]) { return ((ptr[0] << 8) + ptr[1]); } INT Get3Byte(CHAR ptr[]) { CHAR ptr0; ptr0 = ptr[0] & 127; if (ptr[0] > 127) return -((ptr0 << 16) + (ptr[1] << 8) + ptr[2]); else return ((ptr0 << 16) + (ptr[1] << 8) + ptr[2]); } INT Get4Byte(CHAR ptr[]) { CHAR ptr0; ptr0 = ptr[0] & 127; if (ptr[0] > 127) return -((ptr0 << 24) + (ptr[1] << 16) + (ptr[2] << 8) + ptr[3]); else return ((ptr0 << 24) + (ptr[1] << 16) + (ptr[2] << 8) + ptr[3]); } void confp(double fval,INT *iexp, INT *imant) { INT isign; isign = 0; if (fval < 0.0) { isign = 128; fval = -fval; } if (fval == 0.0) *iexp = 0; else *iexp = log(fval) * (1.0/log(16.0)) + 64.0 + (1.0/log(2.0)); if (*iexp < 0) *iexp = 0; if (*iexp > 127) *iexp = 127; *imant = fval/pow(16.0,(*iexp)-70.0); *iexp += isign; } double decfp(INT iexp, INT imant) { INT Negative; double Result; Negative = iexp > 127; Result = imant * pow(16.0,(REAL)((iexp & 127)-64)) * (1.0 /(1 << 24)); if (Negative) return(-Result); else return( Result); } /* ================================================ */ /* gauaw MUST use 64 bit real arithmetic (double) ! */ /* there is no convergence for 32 bit reals (float) */ /* ================================================ */ void gauaw(double *pa, double *pw, INT k) { double c,xz,pk,pkm1,pkm2,pkmrk,sp; INT k_half,is,iter,n; k_half = k >> 1; c = 1.0 / sqrt((0.5+k) * (0.5+k) + (1.0 - 4.0/(M_PI*M_PI)) * 0.25); memcpy(pa,bessel,k_half*sizeof(double)); pa[k_half] = 0.0; for (is = 0; is < k_half; is++) { xz = cos(pa[is] * c); iter = 0; while (++iter < 10) { pkm1 = xz; pkm2 = 1.0; for (n = 1; n < k; n++) { pk = ((n+n+1.0) * xz * pkm1 - n * pkm2) / (n+1.0); pkm2 = pkm1; pkm1 = pk; } pkm1 = pkm2; pkmrk = (k * (pkm1 - xz*pk)) / (1.-xz*xz); sp = pk / pkmrk; xz -= sp; if (sp < 1.0e-14 && sp > -1.0e-14) break; } if (iter > 9) Abort(" *** No convergence in gauaw ***\n"); pa[ is ] = xz; pa[k-is-1] = -xz; pw[ is ] = (2.0 * (1.0-xz*xz)) / (pkm1*pkm1*k*k); pw[k-is-1] = pw[is]; } } void msort(UINT *v, INT *c, INT n) { INT gap, i, j; UINT *o, vtem; if (n == 1) { c[0] = 0; return; } o = (UINT *) malloc(n * sizeof(UINT)); for (i = 0; i < n; i++) o[i] = v[i]; for (gap = n/2; gap > 0; gap /= 2) { for (i = gap; i < n; i++) { for (j = i-gap; j >= 0 && v[j] > v[j+gap]; j -= gap) { vtem = v[j]; v[j] = v[j+gap]; v[j+gap] = vtem; } } } for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { if (v[i] == o[j]) c[i] = j; } } free(o); } void gribmap(int fields, int dates, UINT *vars, UINT *levs, UINT *tims, int *dpos, int *bpos, int *bits, float *zsca, float *fmin) { INT moff, foff, i, j, t, f, products; INT *dindx; UINT *dinfo, *dtemp; FILE *mfile; char mname[128]; struct gaindx *pindx; products = GribCount; dtemp = (UINT *) malloc(products * sizeof(UINT)); dinfo = (UINT *) malloc(products * sizeof(UINT)); dindx = (INT *) malloc(products * sizeof(INT)); if (Verbose) { printf("fields %3d dates %3d\n",fields,dates); for (j = 0; j < products; j++) { printf("prevx %5d: %9d %3d %6d %2d %9d %9d %7.3f %7.3f\n", j,(INT)tims[j],(INT)vars[j],(INT)levs[j], bits[j],bpos[j],dpos[j],dsca[j],fmin[j]); } } for (i = 0; i < dates; i++) dinfo[i] = AllDate[i].YMDH; msort(dinfo,dindx,dates); for (j = 0; j < products; j++) { for (i = 0; i < dates ; i++) { if (dinfo[i] == tims[j]) dtemp[j] = (UINT)dindx[i]; } } for (j = 0; j < products; j++) tims[j] = dtemp[j]; for (i = 0; i < Levels; i++) dinfo[i] = (UINT)DefLevs[i]; for (j = 0; j < products; j++) { for (i = 0; i < Levels; i++) { if (dinfo[i] == levs[j]) dtemp[j] = (UINT)i; } } for (j = 0; j < products; j++) levs[j] = dtemp[j]; for (j = 0; j < products; j++) { dinfo[j] = levs[j] + vars[j] * 1000 + tims[j] * 1000000; } if (Verbose) { /* mfile = fopen("dinfo","wb"); */ /* fwrite(dinfo,sizeof(UINT),products,mfile); */ /* fclose(mfile); */ for (j = 0; j < products; j++) { printf("dinfo[%5.5d] %10d\n",j,(INT)dinfo[j]); } } msort(dinfo,dindx,products); sprintf(mname,"%s%s",DatHNew,".gmp"); mfile = fopen(mname,"wb"); if (mfile == NULL) Abort("\n *** Could not open map file ***\n"); pindx = (struct gaindx *)malloc(sizeof(struct gaindx)); pindx->type = 1; /* GRIB type */ pindx->hinum = 4; pindx->hfnum = 0; pindx->intnum = 1 * 3 * products; pindx->fltnum = 1 * 3 * products; pindx->hipnt = (int *)malloc(sizeof(int) *pindx->hinum); pindx->intpnt = (int *)malloc(sizeof(int) *pindx->intnum); pindx->fltpnt = (float *)malloc(sizeof(float)*pindx->fltnum); for (i=0; iintnum; i++) *(pindx->intpnt+i) = -999; *(pindx->hipnt) = 1; *(pindx->hipnt+1) = dates; *(pindx->hipnt+2) = fields; *(pindx->hipnt+3) = 255; for (j = 0; j < products; j++) { t = dindx[j]; if (Verbose) printf("pindx %5d: %9d %3d %6d %2d %9d %9d %7.3f %7.3f\n", t,(INT)tims[t],(INT)vars[t],(INT)levs[t], bits[t],bpos[t],dpos[t],dsca[t],fmin[t]); foff = t; moff = j * 3; *(pindx->intpnt+moff) = dpos[foff]; /* position of datarec */ *(pindx->intpnt+moff+1) = bpos[foff]; /* position of maskrec */ *(pindx->intpnt+moff+2) = bits[foff]; /* bits per data item */ *(pindx->fltpnt+moff) = dsca[foff]; /* decimal scale factor */ *(pindx->fltpnt+moff+1) = zsca[foff]; /* binary scale factor */ *(pindx->fltpnt+moff+2) = fmin[foff]; /* reference Value */ } fwrite (pindx,sizeof(struct gaindx),1,mfile); if (pindx->hinum >0) fwrite(pindx->hipnt ,sizeof(int) ,pindx->hinum ,mfile); if (pindx->hfnum >0) fwrite(pindx->hfpnt ,sizeof(float),pindx->hfnum ,mfile); if (pindx->intnum>0) fwrite(pindx->intpnt,sizeof(int) ,pindx->intnum,mfile); if (pindx->fltnum>0) fwrite(pindx->fltpnt,sizeof(float),pindx->fltnum,mfile); fclose(mfile); #ifdef _CRAY if (Verbose) { printf("Closed %s. Converting to IEEE\n",mname); } chf771 =_cptofcd(mname,strlen(mname)); CONVERT_IEEE(chf771, &pindx->hinum ,&pindx->hfnum , &pindx->intnum,&pindx->fltnum); #endif } INT GetBlock3and4(CHAR *igrib) { INT Section_Len, Flag; INT imiss, iskip; INT i, j, m, ibits, jlenf, locnp, BytePerValue; INT iflag, irep, lnil, cplx; INT jscale, iexp, imant; INT *lgrib; if (GribEdition == 1) igrib += 8; else igrib += 4; /***************************************/ /* block 1 - product definition block. */ /***************************************/ Section_Len = Get3Byte(igrib); Flag = igrib[ 7]; igrib += Section_Len; /*************************************/ /* block 2 - grid description block. */ /*************************************/ Section_Len = Get3Byte(igrib); if (Flag & 128) igrib += Section_Len; /****************************/ /* block 3 - bit map block. */ /****************************/ if (Flag & 64) { Block3Length = Get3Byte(igrib); igrib += Block3Length; } else { Block3Length = 0; } /********************************/ /* block 4 - binary data block. */ /********************************/ /* get length of binary data block. */ Block4Length = Get3Byte(igrib); /* 4 bit flag / 4 bit count of unused bits at end of block octet. */ iflag = igrib[3]; /* 0------- grid point */ /* 1------- spherical harmonics */ irep = iflag >> 7; if (irep == 1) Abort(" *** Found Spectral Data ***\n"); /* -0------ simple packing */ /* -1------ complex packing */ cplx = (iflag >> 6) & 1; if (cplx == 1) Abort(" *** Found Complex packed Data ***\n"); /* ----++++ number of unused bits at end of section) */ lnil = iflag & 15; /* scale factor (2 byte) */ jscale = Get2Byte(igrib+4); if (jscale > 32767) jscale = 32768 - jscale; /* reference value in GRIB format (iexp,imant) */ iexp = igrib[6]; imant = (igrib[7] << 16) + (igrib[8] << 8) + igrib[9]; FMin = decfp(iexp,imant); ZScale = pow(2.0,(REAL)jscale); /* get number of bits in each data value. */ ibits = igrib[10]; BytePerValue = ibits >> 3; return ibits; } /* ======= */ /* GetInfo */ /* ======= */ void GetInfo(CHAR *igrib) { CHAR *pgrib; INT Flag, Century, TimeUnit, Time1, TimeI; INT IScale, SMFlag; pgrib = igrib + 4 + Grib1Offset; /***************************************/ /* block 1 - product definition block. */ /***************************************/ CenterID = pgrib[ 4]; ModelID = pgrib[ 5]; Flag = pgrib[ 7]; CodeGrib = pgrib[ 8]; TypeOfLevel = pgrib[ 9]; LevelGrib = pgrib[10]; if (Flag & 64) IsMask = TRUE; else IsMask = FALSE; if (TypeOfLevel == 100 || TypeOfLevel == 103 || TypeOfLevel == 105 || TypeOfLevel == 107 || TypeOfLevel == 109 || TypeOfLevel == 160 || TypeOfLevel == 0 || TypeOfLevel == 99 ) LevelGrib = Get2Byte(pgrib+10); if (TypeOfLevel == 100) LevelGrib *= 100; if (TypeOfLevel == 100 || TypeOfLevel == 99) PressureLevelData = TRUE; if (GribEdition > 0) Century = pgrib[24] - 1; else Century = 0; NewDate.Year = pgrib[12]; NewDate.Month = pgrib[13]; NewDate.Day = pgrib[14]; NewDate.Hour = pgrib[15]; NewDate.Minute = pgrib[16]; TimeUnit = pgrib[17]; Time1 = Get2Byte(pgrib+18); TimeI = pgrib[20]; Average = Get2Byte(pgrib+21); Century = pgrib[24]; IScale = Get2Byte(pgrib+26); if (GribEdition == 1) if (IScale) DScale = pow(10.0,(REAL)IScale); else DScale = 1.0; else DScale = 1.0; if (Verbose) { printf("YY %2.2d MM %2.2d DD %2.2d HH %2.2d ", NewDate.Year,NewDate.Month,NewDate.Day, NewDate.Hour); printf("TU %2.2d T1 %6d ",TimeUnit,Time1); } /* Time Code for old ECHAM-1 runs */ if (Time1 && TimeI != 0) { /* NewDate.Year = 1; */ switch (TimeUnit) { case 1: NewDate.Hour += Time1; break; case 2: NewDate.Hour += Time1 * 24; break; case 3: NewDate.Hour += Time1 * 24 * 30; break; case 4: NewDate.Hour += Time1 * 24 * 30 * 12; break; case 5: NewDate.Hour += Time1 * 24 * 30 * 12 * 10; break; case 6: NewDate.Hour += Time1 * 24 * 30 * 12 * 30; break; case 7: NewDate.Hour += Time1 * 24 * 30 * 12 * 100; break; default: Abort("\n *** Unknown GRIB Time-Unit ***\n"); } while (NewDate.Hour >= 24) { NewDate.Hour -= 24; NewDate.Day++; } while (NewDate.Day > 30) { NewDate.Day -= 30; NewDate.Month++; } while (NewDate.Month > 12) { NewDate.Month -= 12; NewDate.Year++; } } if (GribEdition == 1 && Century > 0) NewDate.Year += 100 * (Century - 1); NewDate.YMDH = 1000000 * NewDate.Year + 10000 * NewDate.Month + 100 * NewDate.Day + NewDate.Hour; if (Verbose) { printf("%10d\n",(INT)NewDate.YMDH); } } char *LevelUnits[] = { "pressure level [hPa]", "height level [m]", }; char *GribAbortMsg[] = { "==== G R I B - E R R O R ====\n", "Could not read input file\n", "Synchronization failed\n", "Uncomplete record\n", "Record too long\n", "Block read error\n", "Too many blocks\n", "Non supported GRIB version\n" }; void GribAbort(INT errno) { printf(GribAbortMsg[0]); printf(GribAbortMsg[errno]); printf(GribAbortMsg[0]); exit(1); } INT readgrib(FILE *fp, CHAR *grib, INT maxbytes) { INT RecordBytes; INT RecordFill ; INT Flag, Bits ; INT error; GribPos = 0; if (fread(grib,1,4,fp) < 1) return(1); if (strncmp((char *)grib,"GRIB",4)) { /* Try to synchronize with next GRIB record */ while (++GribPos < MAX_Garbage) { grib[0] = grib[1]; grib[1] = grib[2]; grib[2] = grib[3]; if (fread(grib+3,1,1,fp) < 1) return(1); if (strncmp((char *)grib,"GRIB",4) == 0) break; } if (GribPos >= MAX_Garbage) return(2); } grib += 4; fread(grib,1,4,fp); GribEdition = grib[3]; if (GribEdition == 0) { RecordBytes = 4; Grib1Offset = 0; } else if (GribEdition == 1) { RecordBytes = 8; GribLength = (grib[0] << 16) + (grib[1] << 8) + grib[2]; grib += 4; Grib1Offset = 4; } if (fread(grib+(1-GribEdition)*4,1,4+GribEdition*4,fp) < 1) return(3); Block1Length = (grib[0] << 16) + (grib[1] << 8) + grib[2]; Flag = grib[7]; if (fread(grib+8,1,Block1Length-8,fp) < 1) return(5); RecordBytes += Block1Length; grib += Block1Length; if (Flag & 128) { if (fread(grib,1,4,fp) < 1) return(3); Block2Length = (grib[0] << 16) + (grib[1] << 8) + grib[2]; if (fread(grib+4,1,Block2Length-4,fp) < 1) return(5); RecordBytes += Block2Length; grib += Block2Length; } else { Block2Length = 0; } if (Flag & 64) { if (fread(grib,1,4,fp) < 1) return(3); Block3Length = (grib[0] << 16) + (grib[1] << 8) + grib[2]; if (fread(grib+4,1,Block3Length-4,fp) < 1) return(5); RecordBytes += Block3Length; grib += Block3Length; } else { Block3Length = 0; } if (fread(grib,1,4,fp) < 1) return(3); Block4Length = (grib[0] << 16) + (grib[1] << 8) + grib[2]; RecordBytes += 4; GribLength = 4 + Grib1Offset + Block1Length + Block2Length + Block3Length + Block4Length + 4; if (Process) { if (fread(grib+4,1,GribLength-RecordBytes,fp) < 1) return(3); } else { if (fseek(fp,GribLength-RecordBytes,SEEK_CUR) != 0) Abort("\n *** readgrib: fseek error *** \n"); } return(0); } void Prep_GrADS_DefHori(void) { UINT *cinfo; INT *cindx; INT i,j,code; INT xdim, ydim; INT y1st, yinc; REAL x1st, xinc; UINT dmin, dmax, ddiff; INT imin, imax; INT numofcodes = 0; INT ydigits, InFileMM, InFileYY; FILE *def_hori; sprintf(DefHori,"%s%s",DatHNew,".ctl"); def_hori = fopen(DefHori,"w"); if (def_hori == 0) { printf("could not write GrADS file %s\n",DefHori); return; } fpr1(def_hori,"* GrADS Descriptor file \n"); fpr2(def_hori,"* generated by grb2gas %s\n",Release); fpr3(def_hori,"TITLE %-4.4s/%-4.4s",Model,ExpID); fpr1(def_hori,"\n*\n"); if (EastWest || !SouthNorth || TopDown ) { fpr1(def_hori,"OPTIONS"); if (EastWest) fpr1(def_hori," XREV"); if (!SouthNorth) fpr1(def_hori," YREV"); if (TopDown) fpr1(def_hori," ZREV"); fpr1(def_hori,"\n*\n"); } fpr2(def_hori,"DSET ^%s\n",DatHNew); fpr1(def_hori,"DTYPE grib\n"); fpr2(def_hori,"INDEX ^%s.gmp\n*\n",DatHNew); fpr2(def_hori,"UNDEF %8.1e\n*\n",Undef); x1st = FirstXPoint; xinc = XIncrement; xdim = Lons; ydim = Lats; y1st = 0; yinc = 1; if (xdim == 1) fpr3(def_hori,"XDEF %3d LEVELS %8.4f\n*\n", xdim,x1st); else fpr4(def_hori,"XDEF %3d LINEAR %8.4f %8.4f\n*\n", xdim,x1st,xinc); if (RepType == 0) { if (ydim == 1) fpr3(def_hori,"YDEF %3d LEVELS %8.4f\n*\n", ydim,GrADSLats[y1st]); else fpr4(def_hori,"YDEF %3d LINEAR %8.4f %8.4f\n", ydim,GrADSLats[y1st],YIncrement); } else { if (ydim == 1) { fpr3(def_hori,"YDEF %3d LEVELS %9.4f\n", ydim,GrADSLats[y1st]); } else { fpr3(def_hori,"YDEF %3d LEVELS %9.4f\n", ydim,GrADSLats[y1st]); for (i = y1st+yinc; i < yinc*ydim; i += yinc) { fpr2(def_hori,"%9.4f",GrADSLats[i]); if ((i - y1st) % 8 == 0 || i == (yinc*ydim - 1)) fpr1(def_hori,"\n"); } } } fpr1(def_hori,"*\n"); if (PressureLevelData) fpr1(def_hori,"* You have pressure levels [hPa]\n"); if (Levels == 0) fpr1(def_hori,"ZDEF 1 LEVELS 0.00\n"); else fpr2(def_hori,"ZDEF %3d LEVELS ",Levels); for (j = 0; j < Levels; j++) { if (PressureLevelData) fpr2(def_hori,"%8.2f",(REAL)DefLevs[j]/(REAL)100.); else fpr2(def_hori,"%8d",DefLevs[j]); if ((j % 8 == 0) || (j == (Levels - 1))) fpr1(def_hori,"\n"); } fpr1(def_hori,"*\n"); imin = imax = 0; dmin = dmax = AllDate[0].YMDH; for (i = 0; i < TermCount; i++) { if (AllDate[i].YMDH < dmin) { dmin = AllDate[i].YMDH; imin = i; } if (AllDate[i].YMDH > dmax) { dmax = AllDate[i].YMDH; imax = i; } } dmin = AllDate[imin].Hour + (AllDate[imin].Day - 1) * 24 + (AllDate[imin].Month - 1) * 720 + (AllDate[imin].Year - 1) * 8640 ; dmax = AllDate[imax].Hour + (AllDate[imax].Day - 1) * 24 + (AllDate[imax].Month - 1) * 720 + (AllDate[imax].Year - 1) * 8640 ; if (TermCount > 1) ddiff = (dmax - dmin) / (TermCount - 1); else ddiff = 0; if (Verbose) { printf("dmin %8d dmax %8d ddiff %8d\n", (INT)dmin,(INT)dmax,(INT)ddiff); printf("AllDate[imin].YMDH %8d AllDate[imax].YMDH %8d\n", (INT)AllDate[imin].YMDH,(INT)AllDate[imax].YMDH); } fpr2(def_hori,"TDEF %3d LINEAR ",TermCount); fpr3(def_hori," %2.2dZ%2.2d", AllDate[0].Hour,AllDate[0].Day); fpr3(def_hori,"%3s%04d ",MoName[AllDate[0].Month],AllDate[0].Year); if (ddiff % 8640 == 0) { ddiff = ddiff / 8640; if (ddiff == 0) fpr1(def_hori,"1mo\n"); else fpr2(def_hori,"%3dyr\n",(INT)ddiff); } else { if (ddiff % 720 == 0) { ddiff = ddiff / 720; fpr2(def_hori,"%3dmo\n",(INT)ddiff); } else { if (ddiff % 24 == 0) { ddiff = ddiff / 24; fpr2(def_hori,"%6ddy\n",(INT)ddiff); } else fpr2(def_hori,"%6dhr\n",(INT)ddiff); } } /* Variable definition */ numofcodes = CodeCount; fpr2(def_hori,"VARS %3d\n",numofcodes); cindx = (INT *) malloc(numofcodes*sizeof( INT)); cinfo = (UINT *) malloc(numofcodes*sizeof(UINT)); for (i = 0; i < numofcodes; i++) cinfo[i] = (UINT)(Uni[i].code); msort(cinfo,cindx,numofcodes); for (j = 0; j < numofcodes; j++) { i = cindx[j]; if (Uni[i].olev == 0) Uni[i].olev = Zero; else Uni[i].olev = Levels; fpr5(def_hori,"%-8.8s %4d %3d,%3d,0", Uni[i].name,Uni[i].olev,Uni[i].code,Uni[i].ltyp); fpr3(def_hori," :%s:%s:::\n",Uni[i].expl,Uni[i].unit); } fpr1(def_hori,"ENDVARS\n"); if (def_hori) fclose(def_hori); free(cindx); free(cinfo); } void PostProcess(void) { if (Verbose) printf(" -> Processed Year %4d Month %3d\n", OldDate.Year, OldDate.Month); if (EndOfProcess) { if (FileCount < 1) MulFil = FALSE; EndOfProcess = FALSE; Prep_GrADS_DefHori(); fields = GribCount / TermCount; dates = TermCount; gribmap(fields, dates, vars, levs, tims, dpos, bpos, bits, zsca, fmin); } } /* ================= */ /* switch input file */ /* ================= */ void GetGrADSLats(CHAR *igrib) { INT GridListFlag; INT StagGridFlag; INT GridCoorOffset; INT NumOfLats; INT i,j; double *gmu, *gwt; GrADSLats = (REAL *)malloc(Lats*sizeof(REAL)); if (RepType == 0) { for (i = 0; i < Lats; i++) { if (FirstYPoint > LastYPoint) GrADSLats[i] = LastYPoint + i * YIncrement; else GrADSLats[i] = FirstYPoint + i * YIncrement; } } else { gmu = (double *) malloc(Lats * sizeof(double)); gwt = (double *) malloc(Lats * sizeof(double)); gauaw(gmu,gwt,Lats); for (i = 0; i < Lats; i++) { GrADSLats[i] = 180. * asin(gmu[Lats-(i+1)]) / M_PI; } free(gmu); free(gwt); } return; } /*****************/ /* Ocean Control */ /*****************/ void ProcessControl(void) { INT i,k,code,term; INT griberr, fpos = 0; INT IBits, NLons, NLats, NDim, imiss; INT NewCode, OldCode = -1, NewTerm; INT NewGribLength, NewGribOffset; while (1) { griberr = readgrib(intape,grib,DimGP*5); if (griberr == 0) { GetInfo(grib); IBits = GetBlock3and4(grib); } fpos += GribPos; if (IsMask) bpos[GribCount] = fpos + 4 + Grib1Offset + Block1Length + Block2Length + 6; else bpos[GribCount] = -999; dpos[GribCount] = fpos + 4 + Grib1Offset + Block1Length + Block2Length + Block3Length + 11; bits[GribCount] = IBits; dsca[GribCount] = DScale; zsca[GribCount] = ZScale; fmin[GribCount] = FMin; vars[GribCount] = CodeGrib; levs[GribCount] = LevelGrib; /* tims[GribCount] = NewDate.Hour + */ /* (NewDate.Day - 1) * 24 + */ /* (NewDate.Month - 1) * 720 + */ /* (NewDate.Year - 1) * 8640 ; */ tims[GribCount] = NewDate.YMDH; if (griberr == 0) { GribCount++; fpos += GribLength; } if (griberr == 1) { if (Verbose) printf(" griberr = 1;\n"); if (Verbose) printf(" EndOfProcess = TRUE;\n"); EndOfProcess = TRUE; NewTerm = TRUE; for (term = 0; term < TermCount; term++) { if (memcmp(&OldDate,&AllDate[term],sizeof(struct Date)) == 0) { NewTerm = FALSE; break; } } if (Verbose) { printf("%010.10d\n",(INT)OldDate.YMDH); } if (NewTerm) { AllDate[TermCount] = OldDate; TermCount += 1; } PostProcess(); return; } if (griberr) GribAbort(griberr); /* if (TermCount == 0) */ if (OldCode != CodeGrib) { NewCode = TRUE; for (code = 0; code < CodeCount; code++) { if (CodeGrib == Uni[code].code) { NewCode = FALSE; break; } } if (NewCode) { Uni[CodeCount].code = CodeGrib; Uni[CodeCount].olev = (LevelGrib == 0) ? 0 : Levels; Uni[CodeCount].ltyp = TypeOfLevel; strcpy(Uni[CodeCount].name,All[CodeGrib].name); strcpy(Uni[CodeCount].unit,All[CodeGrib].unit); strcpy(Uni[CodeCount].expl,All[CodeGrib].expl); if (strlen(All[CodeGrib].name) == 0) { sprintf(Uni[CodeCount].name,"var%003d",CodeGrib); sprintf(Uni[CodeCount].unit,"unknown"); sprintf(Uni[CodeCount].expl,"unknown"); } if (Verbose) { printf(" Index %3d Code %3d Levels %3d\n", CodeCount, CodeGrib, Uni[CodeCount].olev); } CodeCount++; } } OldCode = CodeGrib; if ((OldDate.YMDH > 0) && memcmp(&NewDate,&OldDate,sizeof(struct Date))) { NewTerm = TRUE; for (term = 0; term < TermCount; term++) { if (memcmp(&OldDate,&AllDate[term],sizeof(struct Date)) == 0) { NewTerm = FALSE; break; } } if (Verbose) { printf("%010.10d\n",(INT)OldDate.YMDH); } if (NewTerm) { AllDate[TermCount] = OldDate; TermCount += 1; PostProcess(); } } OldDate = NewDate; } } void dimcalc(void) { DimGP = Lons * Lats; if (RepType == 0) { printf(" -> Processing Regular Grid \n" ); } else if (RepType == 4) { printf(" -> Processing Gauss-Grid \n" ); } else { printf(" -> Processing unknown grid type \n" ); printf(" -> ERROR EXIT ...\n" ); exit(1); } printf(" Lons %3d Lats %3d Levs %4d\n", Lons,Lats,Levels); } /* ----------------------------------------------------------- */ /* Extract basic dimension information from GRIB header fields */ /* ----------------------------------------------------------- */ void precntl(void) { CHAR head[15000], *pgrib; INT dumlev[MAX_Levels]; INT headerr; INT ResIndex, RepIndex; INT ExpIndex; INT code, i; INT OldLevel; UINT OldYMDH; INT DeepScanVar, SchonDa = FALSE; INT InputFileYYMM, InputFileYY, InputFileMM; char line[81],cstr[6],name[10],unit[10],expl[60]; char *codemap_path; headerr = readgrib(intape,head,sizeof(head)); if (headerr) GribAbort(headerr); GetInfo(head); OldYMDH = NewDate.YMDH; InputFileYYMM = atol(DatDate); InputFileYY = InputFileYYMM / 100; InputFileMM = InputFileYYMM % 100; ExpIndex = 4+Grib1Offset +46; RepIndex = 4+Grib1Offset+Block1Length+ 5; ResIndex = 4+Grib1Offset+Block1Length+ 6; if (CenterID == 74) { strcpy(Model,"HAD2"); switch (ModelID) { case 1: strcpy(ExpID,"HCCI"); break;; case 2: strcpy(ExpID,"HCGG"); break;; case 3: strcpy(ExpID,"HCGS"); break;; } } else if (CenterID == 7) { strcpy(Model,"GFDL"); switch (ModelID) { case 1: strcpy(ExpID,"GFCI"); break;; case 2: strcpy(ExpID,"GFGG"); break;; case 3: strcpy(ExpID,"GFGS"); break;; } } else if (CenterID == 8) { strcpy(Model,"NCAR"); switch (ModelID) { case 1: strcpy(ExpID,"NCCI"); break;; case 2: strcpy(ExpID,"NCGG"); break;; case 3: strcpy(ExpID,"NCGS"); break;; } } else if (CenterID == 35) { strcpy(Model,"NIES"); switch (ModelID) { case 1: strcpy(ExpID,"NICI"); break;; case 2: strcpy(ExpID,"NIGG"); break;; case 3: strcpy(ExpID,"NIGS"); break;; } } else if (CenterID == 54) { strcpy(Model,"CCCM"); switch (ModelID) { case 1: strcpy(ExpID,"CCCI"); break;; case 2: strcpy(ExpID,"CCGG"); break;; case 3: strcpy(ExpID,"CCGS"); break;; } } else if (CenterID == 79) { strcpy(Model,"E4O3"); switch (ModelID) { case 1: strcpy(ExpID,"MPCI"); break;; case 2: strcpy(ExpID,"MPGG"); break;; case 3: strcpy(ExpID,"MPGS"); break;; } } else if (CenterID == 68) { strcpy(Model,"CCRP"); switch (ModelID) { case 1: strcpy(ExpID,"CSCI"); break;; case 2: strcpy(ExpID,"CSGG"); break;; case 3: strcpy(ExpID,"CSGS"); break;; } } else { strcpy(Model,"ITEM"); strcpy(ExpID,"NONE"); } Lons = Get2Byte(head+ResIndex); Lats = Get2Byte(head+ResIndex+2); RepType = *(head+RepIndex); pgrib = head+ResIndex+21; SouthNorth = (*pgrib >> 6) & 1; EastWest = (*pgrib >> 7) & 1; if (RepType == 0) { FirstYPoint = (REAL)Get3Byte(head+ResIndex+ 4) * 1.0e-3; FirstXPoint = (REAL)Get3Byte(head+ResIndex+ 7) * 1.0e-3; LastYPoint = (REAL)Get3Byte(head+ResIndex+11) * 1.0e-3; LastXPoint = (REAL)Get3Byte(head+ResIndex+14) * 1.0e-3; XIncrement = (REAL)Get2Byte(head+ResIndex+17) * 1.0e-3; YIncrement = (REAL)Get2Byte(head+ResIndex+19) * 1.0e-3; if (FirstYPoint > 90. || FirstYPoint < -90. || LastYPoint < -90. || LastYPoint > 90. || (LastYPoint == FirstYPoint && Lats > 1)) { printf("*** ------------------------------------------------\n"); printf("*** Geographical inconsistency in Latitudes !!!!!!!!\n"); printf("*** FirstYPoint = %6.2f LastYPoint = %6.2f !\n", FirstYPoint, LastYPoint ); if (LastYPoint == FirstYPoint && Lats > 1) { printf("*** ERROR EXIT ...\n"); exit(1); } else { printf("*** Warning but no severe error\n"); } } } else if (RepType == 4) { FirstYPoint = 90.0; FirstXPoint = 0.0; LastYPoint = -90.0; XIncrement = 360. / Lons; LastXPoint = FirstXPoint + (Lons - 1) * XIncrement; YIncrement = 0.0; } else Abort("\n *** Unknown Data Representation Type ***\n"); if (Verbose) { printf("FirstXPoint = %7.3f FirstYPoint = %7.3f\n", FirstXPoint, FirstYPoint); printf(" LastXPoint = %7.3f LastYPoint = %7.3f\n", LastXPoint, LastYPoint); printf(" XIncrement = %7.3f YIncrement = %7.3f\n", XIncrement, YIncrement); } GetGrADSLats(head); vd = fopen(VariDef,"r"); if (vd == NULL) { printf(" Your Variable definition file name isn't adequate\n"); printf(" Check -r option value !\n"); exit(1); } if (Verbose) printf(" Using %s for variable definitions\n",VariDef); fgets(line,80,vd); fgets(line,80,vd); while(fgets(line,80,vd) != NULL) { sscanf(line,"%s %s %s %[^\n]",cstr,name,unit,expl); code = atol(cstr); for (i = 0; i < strlen(name); i++) All[code].name[i] = name[i]; for (i = 0; i < strlen(unit); i++) All[code].unit[i] = unit[i]; for (i = 0; i < strlen(expl); i++) All[code].expl[i] = expl[i]; } if (vd) fclose(vd); while (!SchonDa) { if (LevelGrib) { if (Levels == 0) { dumlev[Levels++] = LevelGrib; } else { SchonDa = FALSE; for (i = 0; i < Levels; i++) { if (dumlev[i] == LevelGrib) SchonDa = TRUE; } if (!SchonDa) dumlev[Levels++] = LevelGrib; } } headerr = readgrib(intape,head,sizeof(head)); if (headerr) { if (!feof(intape)) { GribAbort(headerr); } else { break; } } GetInfo(head); } DefLevs = (INT *)malloc(Levels * sizeof(INT)); if (Verbose) printf("Found %3d Levels :",Levels); for (i = 0; i < Levels; i++) { DefLevs[i] = dumlev[i]; if (Verbose) { if (i % 8 == 0) printf("\n"); printf(" %4d",DefLevs[i]); } } if (Verbose) printf("\n"); /* TopDown = (TypeOfLevel == 100 || TypeOfLevel == 109 || TypeOfLevel == 110) && (DefLevs[0] < DefLevs[Levels]); */ } void Usage(void) { printf("\n%s\ngrb2gas [options] InputFile [OutputFile]\n",Release); printf(" option -r : specify your own variable def-file name\n"); printf(" option -v : verbose mode (extended output)\n"); printf(" option -h : help (this output)\n"); printf(" InputFile : GRIB file\n"); exit(1); } int main(argc,argv) int argc; char *argv[]; { INT code, i; char line[81]; /*********************/ /* print information */ /*********************/ if (argc < 2) Abort("\n *** Missing arguments ***\n"); /***********************/ /* options & filenames */ /***********************/ i = 1; while (i < argc) { if (argv[i][0] == '-') { if (argv[i][1] == 'v') Verbose = 1; else if (argv[i][1] == 'r') { i++; strcpy(VariDef,argv[i]); if (VariDef[0] == '\0') Usage(); } else Usage(); } else if (DatHori[0] == '\0') strcpy(DatHori,argv[i]); else if (DatHNew[0] == '\0') strcpy(DatHNew,argv[i]); else Usage(); i++; } if (DatHori[0] == '\0') { printf("*** Missing input filename ***\n"); Usage(); } printf("*** Input File: %-25s *\n",DatHori); /*******************/ /* open input file */ /*******************/ intape = fopen(DatHori,"rb"); if (intape == 0) { printf("could not open input file %s\n",DatHori); exit(1); } strcpy(DatHNew,DatHori); /******************/ /* pre-processing */ /******************/ precntl(); Process = TRUE; rewind(intape); /*******************/ /* initializations */ /*******************/ dimcalc(); /***********************/ /* open output file(s) */ /***********************/ grib = (CHAR *)malloc(DimGP * 5 * sizeof(CHAR)); YLats = (REAL *)malloc(Lats * sizeof(REAL)); ProcessControl(); fclose(intape); fclose(outape); printf(" -> Grads Control File: %s.ctl\n",DatHNew); printf(" -> Grads GribMap File: %s.gmp",DatHNew); #ifdef _CRAY printf(" and %s.gmp.ieee\n",DatHNew); #else printf("\n"); #endif printf(" -> NORMAL EXIT\n"); return(0); }