README_dkrz_readgrib 640 624 33261 2720 6251741771 7760 dkrz_readgrib is a program that either lists the contents of a GRIB file or converts GRIB records to binary or ASCII. Compile with gcc -o dkrz_readgrib dkrz_readgrib.c QU2REG2.c -lm Usage dkrz_readgrib -i input_file_name \ -o output_file_name \ -t output_type \ [-d debug_flag ] where input_file_name is the GRIB input file name output_file_name is the output file name output_type is either list - ASCII list of contents of input_file_name cona - ASCII output conb - binary output conb and cona give following structure of the output file header (1 record with 8 4-byte-integer) field (1 record with n 4-byte-reals) header field and so on. the header consists of 1. code number 2. level 3. date [yymmdd] 4. hour 5. nlon [or dim of spectral field, if encountered] 6. nlat [or 1, if spectral field] 7. 0 8. 0 so that the dimension of the data field is always given by header(5) * header(6). In case of trouble feel free to call or send mail. Good luck ------------- Deutsches Klimarechenzentrum GmbH --------------- Arno Hellbach email: hellbach@dkrz.de Deutsches Klimarechenzentrum Tel +49 40 41173363 Bundesstr. 55 Fax +49 40 41173298 D-20146 Hamburg --------------------------------------------------------------- dkrz_readgrib.c 750 624 33261 53045 6227427422 7211 #if defined(__convex__) #include #else #include #endif #include #if defined (sun) #include #endif #include #include #include #include #include #define INT long #define REAL float #define CHAR unsigned char #define TRUE 1 #define FALSE 0 #define SEEK_CUR 1 #define MAX_LEVELS 80 #define MAX_Longitudes 320 #define MAX_Latitudes 242 #define MAX_DimGP (MAX_Longitudes * MAX_Latitudes) #define MAX_NVCT (MAX_LEVELS * 2 + 2) INT Grib1Offset; INT nvct; INT Debug = 0; extern FILE *fp = NULL; extern char *ofile = NULL; /* usage: readgrib -t type -i infile -o outfile -d debug */ int main(int argc, char *argv[]) { extern INT gribdate_(); extern int gribin_(); FILE *cfile; static char yifile[80], yofile[80], yptype[80], ydebug[80]; static INT kunit = 11, keof, kret, ksec1[1200], ksec2[32], lens4; static INT i, b, nrec, offset= 0; static INT lt = 0, li = 0, lo = 0, ld = 0; static INT nvar, nlev, hour, ndat, nlat, nlon; static REAL undef = -9999.9, zero = 0.0; static INT one = 1; static INT mask[MAX_DimGP]; static REAL Sec4[MAX_DimGP]; static INT head[8]; if (argc == 1) { fprintf(stdout,"Type output type [list|cona|conb|conp] : "); fscanf (stdin ,"%s",yptype); fprintf(stdout,"Type infile name : "); fscanf (stdin ,"%s",yifile); fprintf(stdout,"Type outfile name : "); fscanf (stdin ,"%s",yofile); fprintf(stdout,"Debug info wanted ? [yes|no] : "); fscanf (stdin ,"%s",ydebug); if (strcmp(ydebug,"yes") == 0) Debug = 1; else Debug = 0; } else { for (i = 1; i < argc; i += 2) { if (strcmp(argv[i],"-t") == 0) { lt = 1; strcpy(yptype,argv[i+1]); } if (strcmp(argv[i],"-i") == 0) { li = 1; strcpy(yifile,argv[i+1]); } if (strcmp(argv[i],"-o") == 0) { lo = 1; strcpy(yofile,argv[i+1]); } if (strcmp(argv[i],"-d") == 0) { ld = 1; strcpy(ydebug,argv[i+1]); if (strcmp(ydebug,"yes") == 0) Debug = 1; else Debug = 0; } } if (li == 0 || lo == 0 || lt == 0) printf("readgrib : interactive usage\n" "readgrib -t type -i infile -o outfile -d debug\n"); if (lt == 0) { fprintf(stdout,"Type output type [list|cona|conb|conp] : "); fscanf (stdin ,"%s",yptype); } if (li == 0) { fprintf(stdout,"Type infile name : "); fscanf (stdin ,"%s",yifile); li = 1; } if (lo == 0) { fprintf(stdout,"Type outfile name : "); fscanf (stdin ,"%s",yofile); lo = 1; } } cfile = fopen(yofile,"w"); if (strcmp(yptype,"list") == 0) { fprintf(cfile,"List of contents of %s\n", yifile); fprintf(cfile," code level yymmdd hh lon lat\n"); fprintf(cfile," ---- ----- ------ -- --- ---\n"); } keof = 0; while (1) { if (strcmp(yptype,"list") == 0) gribin_(&kunit, (REAL *)0, ksec1, ksec2, mask, &lens4, &keof, &kret, yifile, &undef); if (strcmp(yptype,"cona") == 0 || strcmp(yptype,"conb") == 0 || strcmp(yptype,"conp") == 0 ) gribin_(&kunit, Sec4, ksec1, ksec2, mask, &lens4, &keof, &kret, yifile, &undef); if (keof != 0 || kret != 0) break; nvar = ksec1[6-1]; nlev = ksec1[8-1]; nlon = ksec2[2-1]; nlat = ksec2[3-1]; ndat = gribdate_(ksec1); hour = ndat % 100; ndat = ndat / 100; if (strcmp(yptype,"cona") == 0 || strcmp(yptype,"list") == 0 ) { if (ksec2[1-1] == 50) fprintf(cfile,"%8d%8d%8d%8d%8d%8d%8d%8d\n", nvar,nlev,ndat,hour,lens4,(INT)one,(INT)zero,(INT)zero); else fprintf(cfile,"%8d%8d%8d%8d%8d%8d%8d%8d\n", nvar,nlev,ndat,hour,nlon,nlat,(INT)zero,(INT)zero); } if (strcmp(yptype,"cona") == 0) for (i = 0; i < lens4; i++) { fprintf(cfile,"%13.6e",Sec4[i]); if ((i+1) % 8 == 0) fprintf(cfile,"\n"); } if (strcmp(yptype,"conp") == 0) { fwrite((char *)Sec4,sizeof(REAL),lens4,cfile); } if (strcmp(yptype,"conb") == 0) { head[0] = nvar; head[1] = nlev; head[2] = ndat; head[3] = hour; head[4] = nlon; head[5] = nlat; head[6] = zero; head[7] = zero; if (ksec2[1-1] == 50) { head[4] = lens4; head[5] = one; } fwrite(head ,sizeof(INT) , 8,cfile); fwrite((char *)Sec4,sizeof(REAL),lens4,cfile); } } if (kret == 0) printf("readgrib : NORMAL EXIT\n"); else printf("readgrib : ERROR EXIT with return code %d\n",kret); exit(kret); } INT gribdate_(INT *ksec1) { static INT mtinc, mtave, myref, mtuni, mtran; static INT mc, my, mm, md, mh, ay, am; if (ksec1[1-1] == 1) mc = ksec1[21-1] - 1; else mc = 0; if (mc < 0) mc = 0; my = ksec1[10-1]; myref = my - 1; mm = ksec1[11-1]; md = ksec1[12-1]; mh = ksec1[13-1]; if (Debug) printf("gribdate : cc %3d yy %3d mm %3d dd %3d hh %3d\n", mc, my, mm, md, mh); mtuni = ksec1[15-1]; mtinc = ksec1[16-1]; mtran = ksec1[18-1]; mtave = ksec1[19-1]; am = mtave % 100; ay = mtave / 100; if (mtran == 0) { if (mtuni == 1) {mh += mtinc;} if (mtuni == 2) {md += mtinc;} if (mtuni == 3) {mm += mtinc;} if (mtuni == 4) {my += mtinc;} } else { if (mtinc != 0) { if (mtuni == 0) {mh += mtinc / 60;} if (mtuni == 1) {mh += mtinc;} if (mtuni == 2) {mh += mtinc * 24;} if (mtuni == 3) {mh += mtinc * 720;} if (mh >= 24) {mh += -24;} md += mh / 24; mh %= 24; mm += (md - 1) / 30; md = (md - 1) % 30 + 1; my += (mm - 1) / 12; mm = (mm - 1) % 12 + 1; my -= myref; } } return (mh + md * 100 + mm * 10000 + (my + (mc * 100)) * 1000000); } /* gribdate_ */ void Abort(char *errtext) { fprintf(stderr,errtext); printf( errtext); exit(1); } REAL decfp(INT iexp, INT imant) { INT Negative; REAL Result; Negative = iexp > 127; Result = imant * pow(16.0,(REAL)((iexp & 127)-64)) * (1.0 /(1 << 24)); if (Negative) return(-Result); else return( Result); } void ungrib(CHAR *igrib, INT sec1[], INT sec2[], REAL *sec4, INT mask[], REAL *rmiss) { INT dim; INT jlenv,jlenf,lblk,i,m,j; INT ilat1,ilat2; INT ibits,imiss3; INT iexp,imant,iflag,irep,lnil,jscale,imiss; INT masklen,offset; INT ReducedGrid = FALSE, NNV, NPV; INT iwest, ieast, *nvcpr; INT kret, One = 1, nolat, nolng; REAL fmin,zscale; extern int qu2reg2_(); igrib += 4; igrib += Grib1Offset; lblk = ((*igrib) << 16) + ((*(igrib+1)) << 8) + *(igrib+2); igrib += 3; offset = 4; offset += Grib1Offset; offset += lblk; /***************************************/ /* block 1 - product definition block. */ /***************************************/ for (i = 0 ; i < lblk - 3 ; i++) sec1[i] = *igrib++ ; if (Grib1Offset) sec1[0] = 1; else sec1[0] = 0; sec1[18] = (sec1[18] << 8) + sec1[19]; sec1[19] = 0; if (sec1[6] == 100 || sec1[6] == 103 || sec1[6] == 105 || sec1[6] == 107 || sec1[6] == 109 || sec1[6] == 160 ) { sec1[7] = (sec1[7] << 8) + sec1[8]; sec1[8] = 0; } /* one time range can occupy two octets. */ if (sec1[17] == 10) { sec1[15] = (sec1[15] << 8) + sec1[16]; sec1[16] = 0 ; } if (lblk > 28) sec1[43] = (sec1[43] << ((sizeof(INT) - 1) * 8)) + (sec1[44] << ((sizeof(INT) - 2) * 8)) + (sec1[45] << ((sizeof(INT) - 3) * 8)) + (sec1[46] << ((sizeof(INT) - 4) * 8)); sec1[44] = sec1[45] = sec1[46] = 0; /*************************************/ /* block 2 - grid description block. */ /*************************************/ if (sec1[4] == 128 || sec1[4] == 192) { lblk = *igrib++ << 16; lblk += *igrib++ << 8; lblk += *igrib++ ; offset += lblk; /* length is 32 octets for lat/long, gaussian and spherical */ /* harmonics . for any data on hybrid levels the */ /* vertical coordinates are added. */ /* For reduced grid data NNV = 0 and NPV != 255 */ /* number of unused bits at end of block. */ NNV = *igrib++; /* next octet is reserved. */ /* or is the location of either VCT or PIR */ NPV = *igrib++; if (NNV == 0 && NPV != 255 && NPV != 0) ReducedGrid = TRUE; /* data representation type. */ sec2[0] = *igrib++; /* lat/longitude grid, gaussian grid and spherical harmonics */ /* are the only data representations handled. */ if (sec2[0] != 0 && sec2[0] != 4 && sec2[0] != 50) Abort(" grid description block not yet defined\n"); /* lat/long or gaussian grid. */ if (sec2[0] == 0 || sec2[0] == 4) { /* number of lat/long points. */ sec2[1] = *igrib++ << 8; sec2[1] += *igrib++; sec2[2] = *igrib++ << 8; sec2[2] += *igrib++; /* lat/long of origin. */ ilat1 = *igrib++ << 16; ilat1 += *igrib++ << 8; ilat1 += *igrib++ ; ilat2 = *igrib++ << 16; ilat2 += *igrib++ << 8; ilat2 += *igrib++ ; /* if sign bit set to 1 , values are negative. */ if (ilat1 < 8388608) sec2[3] = ilat1; else sec2[3] = 8388608 - ilat1; if (ilat2 < 8388608) sec2[4] = ilat2; else sec2[4] = 8388608 - ilat2; /* resolution flag. */ sec2[5] = *igrib++; /* lat/long of extreme points. */ ilat1 = *igrib++ << 16; ilat1 += *igrib++ << 8; ilat1 += *igrib++ ; ilat2 = *igrib++ << 16; ilat2 += *igrib++ << 8; ilat2 += *igrib++ ; /* if sign bit set to 1 , values are negative. */ if (ilat1 < 8388608) sec2[6] = ilat1; else sec2[6] = 8388608 - ilat1; if (ilat2 < 8388608) sec2[7] = ilat2; else sec2[7] = 8388608 - ilat2; /* direction increments / number of latitude lines. */ sec2[8] = *igrib++ << 8; sec2[8] += *igrib++; sec2[9] = *igrib++ << 8; sec2[9] += *igrib++; /* scanning mode flags. */ sec2[10] = *igrib++; /* 4 reserved octets. */ igrib += 4; } /* spherical harmonic data. */ if (sec2[0] == 50) { /* pentagonal resolution parameters. */ sec2[1] = *igrib++ << 8; sec2[1] += *igrib++ ; sec2[2] = *igrib++ << 8; sec2[2] += *igrib++ ; sec2[3] = *igrib++ << 8; sec2[3] += *igrib++ ; /* representation type and mode. */ sec2[4] = *igrib++; sec2[5] = *igrib++; /* 18 reserved octets. */ igrib += 18; } /* get dimension of 2d-field */ if (sec2[0] == 0) { nolat = (sec2[3] - sec2[6]); iwest = sec2[4]; ieast = sec2[7]; if (iwest <= ieast) iwest += 360000; nolng = iwest - ieast; /* if Di or Dj are missing */ if (sec2[8] == 65535) { nolat = 1 + nolat / sec2[9]; nolng = 1 + nolng / sec2[9]; } else { nolat = 1 + nolat / sec2[9]; nolng = 1 + nolng / sec2[8]; } if (ReducedGrid) sec2[1] = nolng; if (ReducedGrid) sec2[2] = nolat; dim = sec2[1] * sec2[2]; } if (sec2[0] == 4) { nolat = sec2[9] * 2; nolng = sec2[9] * 4; if (ReducedGrid) sec2[1] = nolng; if (ReducedGrid) sec2[2] = nolat; dim = sec2[1] * sec2[2]; } if (sec2[0] == 50) dim = (sec2[1]+1) * (sec2[2]+2); /* vertical coordinate parameters for hybrid levels. */ /* get number of vertical coordinate parameters, if any. */ if (ReducedGrid) jlenv = (lblk - 32) >> 1; else jlenv = (lblk - 32) >> 2; if (jlenv > MAX_NVCT) { fprintf(stderr," No. of vertical coordinates: %d\n",jlenv); fprintf(stderr," Array size: %d\n",MAX_NVCT); exit(1); } if (jlenv) { nvcpr = (INT *) malloc(jlenv * sizeof(INT)); if (ReducedGrid) for (i = 0 ; i < jlenv ; i++) { nvcpr[i] = *igrib++ << 8; nvcpr[i] += *igrib++ ; } else for (i = 0 ; i < jlenv ; i++) igrib += 4; } } /* Leave ungrib (decode only definition sections */ if (sec4 == NULL) return; /********************************************************************/ /* block 3 - bit map block. */ /********************************************************************/ if (sec1[4] == 64 || sec1[4] == 192) { lblk = ((*igrib) << 16) + ((*(igrib+1)) << 8) + *(igrib+2); offset += lblk; /* exit(1); */ imiss = *(igrib+3); imiss3 = imiss >> 3; igrib += 6; masklen = lblk - 6 - imiss3; if (masklen*8 - imiss != sec2[1]*sec2[2]) { printf("gribin : Mask dimension conflict !\n"); exit(1); } for (m = 0; m < masklen; m++) for (j = 7; j >= 0; j--) { mask[(m*8)+7-j] = ((igrib[m]) >> j) & 1; } igrib += masklen; /* for (m = 0; m < masklen*8; m++) { printf("%1.1d",mask[m]); if ((m+1) % 128 == 0) printf("\n"); } */ } else /* preset mask */ for (i = 0; i < dim; i++) mask[i] = 1; /********************************************************************/ /* block 4 - binary data block. */ /********************************************************************/ /* get length of binary data block. */ lblk = *igrib++ << 16; lblk += *igrib++ << 8; lblk += *igrib++ ; /* 4 bit flag / 4 bit count of unused bits at end of block octet. */ iflag = *igrib++; /* 0000---- grid point */ /* 1000---- spherical harmonics */ irep = iflag >> 7; lnil = iflag & 127; /* get scale factor */; jscale = *igrib++ << 8; jscale += *igrib++ ; if (jscale > 32767) jscale = 32768 - jscale; /* get reference value (fmin) in grib format (iexp,imant) */ iexp = *igrib++; imant = *igrib++ << 16; imant += *igrib++ << 8; imant += *igrib++ ; /* check for missing data indicators. */ imiss = (jscale == 65535 && iexp == 255 && imant == 16777215); /* convert reference value and scale factor. */ if (imiss == 0) { fmin = decfp(iexp,imant); zscale = pow(2.0,(REAL)jscale); } /* get number of bits in each data value. */ ibits = *igrib++; /* if data is in spherical harmonic form, next 4 octets */ /* contain the real (0,0) coefficient. */ if (irep == 1) { /* get real (0,0) coefficient in grib format and */ /* convert to floating point. */ iexp = *igrib++; imant = *igrib++ << 16; imant += *igrib++ << 8; imant += *igrib++ ; if (imiss) *sec4++ = 0.0; else *sec4++ = decfp(iexp,imant); } /* decode data values to floating point and store in sec4. */ /* first calculate the number of data values. */ jlenf = lblk - 11 - irep * 4 ; jlenf = (jlenf * 8 - lnil) / ibits ; /* check length of output array. */ if (jlenf+irep > MAX_DimGP) { fprintf(stderr," values to be decoded are %d\n",(jlenf)+1); fprintf(stderr," array size: %d\n",MAX_DimGP); exit(1); } if (imiss) memset(sec4,0.0,dim * sizeof(REAL)); else { switch (ibits) { case 8: for (i = 0; i < jlenf; i++) { if (mask[i] == 1) { *sec4 = fmin + zscale * *igrib ; igrib++; } else *sec4 = *rmiss; sec4++; } break; case 16: for (i = 0; i < jlenf; i++) { if (mask[i] == 1) { *sec4 = fmin + zscale * ((igrib[0] << 8) + igrib[1]); igrib += 2; } else *sec4 = *rmiss; sec4++; } break; case 32: for (i = 0; i < jlenf; i++) { if (mask[i] == 1) { *sec4 = fmin + zscale * ((igrib[0] << 24) + (igrib[1] << 16) + (igrib[2] << 8) + igrib[3]) ; igrib += 4; } else *sec4 = *rmiss; sec4++; } break; default: fprintf(stderr, " Unimplemented packing factor %d\n",ibits); exit(1); } } if (ReducedGrid) { qu2reg2_(sec4-jlenf,nvcpr,&nolat,&nolng,&One,&rmiss,&kret); if (Debug) printf("ungrib : free(nvcpr) at %8x\n",nvcpr); if (nvcpr) free(nvcpr); } } gribin_(INT *unit, REAL *sec4, INT *sec1, INT *sec2, INT *mask, INT *dlen, INT *keof, INT *kret, char *fname, REAL *rmiss) { INT Count ; INT BlockCount ; INT BlockLength; INT RecordBytes; INT GribEdition; INT maxbytes; int flen; char *mfile ; CHAR *grib, *pgrib; flen = strlen(fname); maxbytes = MAX_DimGP * 5; grib = (CHAR *) malloc(maxbytes); pgrib = grib; if (Debug) { printf("gribin : malloc grib at %8x\n",grib); } if (flen > 0) { mfile = malloc (flen+1); strncpy (mfile,fname,flen); *(mfile + flen) = '\040'; *(strchr(mfile,'\040')) = '\0'; if (fp == NULL) { ofile = malloc (128); strncpy (ofile,mfile,flen); printf("gribin : Open %s\n",mfile); } } else { printf("gribin : Invalid filename %s\n",mfile); *kret = 1; return;} if (strcmp(ofile,mfile) != 0) { fclose (fp); fp = NULL; printf("gribin : Input file changed from %s to %s\n",ofile,mfile); strcpy (ofile,mfile); } if (fp == NULL) fp = fopen (mfile,"rb"); if (fp == NULL) { printf("gribin : Open %s failed\n",mfile); *kret = -1; return;} if (*keof == 1) { printf("gribin : Rewound %s\n",mfile); rewind(fp); } if (fread(grib,1,4,fp) < 4) { if (feof(fp)) { printf("gribin : EOF encountered in %s\n",mfile); *keof = 1; *kret = 0; return;} printf("gribin : Read failed\n"); *kret = 2; return;} if (strncmp((char *)grib,"GRIB",4)) { /* Try to synchronize with next GRIB record */ Count = 0; while (++Count < 10240) { grib[0] = grib[1]; grib[1] = grib[2]; grib[2] = grib[3]; if (fread(grib+3,1,1,fp) < 1) if (feof(fp)) { printf("gribin : EOF encountered\n"); *keof = 1; *kret = 0; return;} else { printf("gribin : Read failed\n"); *keof = 1; *kret = 3; return;} if (strncmp((char *)grib,"GRIB",4) == 0) break; } if (Count > 10239) { printf("gribin : Synchronization failed\n"); *kret = 4; return;} } grib += 4; fread(grib,1,4,fp); GribEdition = grib[3]; if (GribEdition == 0) { Grib1Offset = 0; RecordBytes = 8; fseek(fp,(long) -4,SEEK_CUR); for (BlockCount = 1 ; BlockCount < 7 ; BlockCount++) { if (fread(grib,1,4,fp) < 1) { printf("gribin : Incomplete record\n"); *kret = 6; return;} if (strncmp(grib,"7777",4) == 0) break; BlockLength = (grib[0] << 16) + (grib[1] << 8) + grib[2]; RecordBytes += BlockLength; if (RecordBytes > maxbytes) { printf("gribin : GRIB record too long\n"); *kret = 7; return;} if (fread(grib+4,1,BlockLength-4,fp) < BlockLength-4) { printf("gribin : Block read error\n"); *kret = 8; return;} grib += BlockLength; } RecordBytes += 4; } else if (GribEdition == 1) { Grib1Offset = 4; RecordBytes = (grib[0] << 16) + (grib[1] << 8) + grib[2]; grib += 4; if (fread(grib,1,RecordBytes-8,fp) != RecordBytes-8) { printf("gribin : Block read error\n"); *kret = 8; return;} } else { printf("gribin : GRIB Version not known\n"); *kret = 5; return;} ungrib(pgrib,sec1,sec2,sec4,mask,rmiss); if (sec2[0] == 0 || sec2[0] == 4) *dlen = sec2[1]*sec2[2]; if (sec2[0] == 50) *dlen = (sec2[1]+1) * (sec2[1]+2); if (grib) { if (Debug) printf("gribin : free(pgrib) at %8x\n",pgrib); free(pgrib); } *keof = 0; *kret = 0; return; } masklen = lblk - 6 - imiss3; if (masklen*8 - imiss != sec2[1]*sec2[2]) { printf("gribin : Mask dimension conflict !\n"); exit(1); } for (m = 0; m < masklen; m++) for (j = 7; j >= 0; j--) { mask[(m*8)+7-j] = ((igrib[m]) >> j) & 1; } igrib += masklen; /* for (m = 0; m < masklen*8; m++) { printf("%1.1d",mask[m]); if ((m+1) % 128 == 0) printf("\n"); } */ } else /* prQU2REG2.c 640 624 33261 26133 6123256767 5432 #include #if defined (sun) #include #endif #include #include #include #include #include #define INT long #define REAL float #define CHAR unsigned char #define TRUE 1 #define FALSE 0 #define SEEK_CUR 1 #define abs(x) ((x) >= 0 ? (x) : -(x)) int scm0_(REAL *pdl, REAL *pdr, REAL *pfl, REAL *pfr, INT *klg) { /* System generated locals */ INT i_1; REAL r_1; /* Local variables */ static REAL zfac, zeps, zbeta; static INT jl; static REAL zalpha; /* **** SCM0 - Apply SCM0 limiter to derivative estimates. */ /* output: */ /* pdl = the limited derivative at the left edge of the interval */ /* pdr = the limited derivative at the right edge of the interval */ /* inputs */ /* pdl = the original derivative at the left edge */ /* pdr = the original derivative at the right edge */ /* pfl = function value at the left edge of the interval */ /* pfr = function value at the right edge of the interval */ /* klg = number of intervals where the derivatives are limited */ /* define constants */ zeps = (REAL)1e-12; zfac = ((REAL)1. - zeps) * (REAL)3.; i_1 = *klg; for (jl = 0; jl < i_1; ++jl) { if ((r_1 = pfr[jl] - pfl[jl], abs(r_1)) > zeps) { zalpha = pdl[jl] / (pfr[jl] - pfl[jl]); zbeta = pdr[jl] / (pfr[jl] - pfl[jl]); if (zalpha <= (REAL)0.) pdl[jl] = (REAL)0.; if (zbeta <= (REAL)0.) pdr[jl] = (REAL)0.; if (zalpha > zfac) pdl[jl] = zfac * (pfr[jl] - pfl[jl]); if (zbeta > zfac) pdr[jl] = zfac * (pfr[jl] - pfl[jl]); } else { pdl[jl] = (REAL)0.; pdr[jl] = (REAL)0.; } } } /* scm0_ */ int rowina2_(REAL *p, INT *ko, INT *ki, REAL *pw, INT *kcode, REAL *pmsval, INT *kret) { /* System generated locals */ INT pw_dim1, pw_offset, i_1; /* Local variables */ static REAL zwt1, zrdi, zpos; static INT jl, ip; static REAL zdo, zwt; extern int scm0_(); /* Parameter adjustments */ --p; pw_dim1 = *ko + 3; pw_offset = pw_dim1; pw -= pw_offset; /* **** ROWINA2 - Interpolation of row of values. */ /* Input Parameters. */ /* ----------------- */ /* P - Row of values to be interpolated. */ /* Dimension must be at least KO. */ /* KO - Number of values required. */ /* KI - Number of values in P on input. */ /* PW - Working array. */ /* Dimension must be at least (0:KO+2,3). */ /* KCODE - Interpolation required. */ /* 1 , linear. */ /* 3 , cubic. */ /* PMSVAL - Value used for missing data indicator. */ /* Output Parameters. */ /* ------------------ */ /* P - Now contains KO values. */ /* KRET - Return code */ /* 0, OK */ /* Non-zero, error */ /* Author. */ /* ------- */ /* J.D.Chambers ECMWF 22.07.94 */ /* ******************************** */ /* Section 1. Linear interpolation .. */ /* ******************************** */ *kret = 0; if (*kcode == 1) { /* Move input values to work array */ i_1 = *ki; for (jl = 1; jl <= i_1; ++jl) pw[jl + pw_dim1] = p[jl]; /* Arrange wrap-around value in work array */ pw[*ki + 1 + pw_dim1] = p[1]; /* Set up constants to be used to figure out weighting for */ /* values in interpolation. */ zrdi = (REAL) (*ki); zdo = (REAL)1. / (REAL) (*ko); /* Loop through the output points */ i_1 = *ko; for (jl = 1; jl <= i_1; ++jl) { /* Calculate weight from the start of row */ zpos = (jl - 1) * zdo; zwt = zpos * zrdi; /* Get the current array position(minus 1) from the weight - */ /* note the implicit truncation. */ ip = zwt; /* If the left value is missing, use the right value */ if (pw[ip + 1 + pw_dim1] == *pmsval) { p[jl] = pw[ip + 2 + pw_dim1]; /* If the right value is missing, use the left value */ } else if (pw[ip + 2 + pw_dim1] == *pmsval) { p[jl] = pw[ip + 1 + pw_dim1]; /* If neither missing, interpolate ... */ } else { /* Adjust the weight to range (0.0 to 1.0) */ zwt -= ip; /* Interpolate using the weighted values on either side */ /* of the output point position */ p[jl] = ((REAL)1. - zwt) * pw[ip + 1 + pw_dim1] + zwt * pw[ip + 2 + pw_dim1]; } } /* ******************************* */ /* Section 2. Cubic interpolation .. */ /* ******************************* */ } else if (*kcode == 3) { i_1 = *ki; for (jl = 1; jl <= i_1; ++jl) { if (p[jl] == *pmsval) { fprintf(stderr," ROWINA2: "); fprintf(stderr," Cubic interpolation not supported"); fprintf(stderr," for fields containing missing data.\n"); *kret = 1; goto L900; } pw[jl + pw_dim1] = p[jl]; } pw[pw_dim1] = p[*ki]; pw[*ki + 1 + pw_dim1] = p[1]; pw[*ki + 2 + pw_dim1] = p[2]; i_1 = *ki; for (jl = 1; jl <= i_1; ++jl) { pw[jl + (pw_dim1 << 1)] = - pw[jl - 1 + pw_dim1] / (REAL)3. - pw[jl + pw_dim1] * (REAL).5 + pw[jl + 1 +pw_dim1] - pw[jl + 2 + pw_dim1] / (REAL)6.; pw[jl + 1 + pw_dim1 * 3] = pw[jl - 1 + pw_dim1] / (REAL)6. - pw[jl + pw_dim1] + pw[jl + 1 + pw_dim1] * (REAL).5 + pw[jl + 2 + pw_dim1] / (REAL)3.; } scm0_(&pw[(pw_dim1 << 1) + 1], &pw[pw_dim1 * 3 + 2], &pw[pw_dim1 + 1], &pw[pw_dim1 + 2], ki); zrdi = (REAL) (*ki); zdo = (REAL)1. / (REAL) (*ko); i_1 = *ko; for (jl = 1; jl <= i_1; ++jl) { zpos = (jl - 1) * zdo; zwt = zpos * zrdi; ip = zwt + 1; zwt = zwt + (REAL)1. - ip; zwt1 = (REAL)1. - zwt; p[jl] = (((REAL)3. - zwt1 * (REAL)2.) * pw[ip + pw_dim1] + zwt * pw[ip + (pw_dim1 << 1)]) * zwt1 * zwt1 + (((REAL)3. - zwt * (REAL)2.) * pw[ip + 1 + pw_dim1] - zwt1 * pw[ip + 1 + pw_dim1 * 3]) * zwt * zwt; } } else { /* ************************************** */ /* Section 3. Invalid interpolation code .. */ /* ************************************** */ fprintf(stderr," ROWINA2:"); fprintf(stderr," Invalid interpolation code = %2d\n",*kcode); *kret = 2; } L900: return 0; } /* rowina2_ */ int qu2reg2_(REAL *pfield, INT *kpoint, INT *klat, INT *klon, INT *kcode, REAL *pmsval, INT *kret) { /* System generated locals */ INT i_1, i_2; /* Local variables */ static INT ilii, ilio, icode; static REAL zline[640], ztemp[204800], zwork[1929] /* was [643][3] */; static INT iregno, iquano, j210, j220, j230, j240, j225; extern int rowina2_(); /* Parameter adjustments */ --pfield; --kpoint; /* **** QU2REG - Convert quasi-regular grid data to regular. */ /* Input Parameters. */ /* ----------------- */ /* PFIELD - Array containing quasi-regular grid */ /* data. */ /* KPOINT - Array containing list of the number of */ /* points on each latitude (or longitude) of */ /* the quasi-regular grid. */ /* KLAT - Number of latitude lines */ /* KLON - Number of longitude lines */ /* KCODE - Interpolation required. */ /* 1 , linear - data quasi-regular on */ /* latitude lines. */ /* 3 , cubic - data quasi-regular on */ /* latitude lines. */ /* 11, linear - data quasi-regular on */ /* longitude lines. */ /* 13, cubic - data quasi-regular on */ /* longitude lines. */ /* PMSVAL - Value used for missing data indicator. */ /* Output Parameters. */ /* ------------------ */ /* KRET - return code */ /* 0 = OK */ /* non-zero indicates fatal error */ /* PFIELD - Array containing regular grid data. */ /* Author. */ /* ------- */ /* J.D.Chambers ECMWF 22.07.94 */ /* J.D.Chambers ECMWF 13.09.94 */ /* Add return code KRET and remove calls to ABORT. */ /* ------------------------------ */ /* Section 1. Set initial values. */ /* ------------------------------ */ *kret = 0; /* Check input parameters. */ if (*kcode != 1 && *kcode != 3 && *kcode != 11 && *kcode != 13) { fprintf(stderr," QU2REG :"); fprintf(stderr," Invalid interpolation type code = %2d\n",*kcode); *kret = 1; goto L900; } if (*klat > 320) { fprintf(stderr," QU2REG :"); fprintf(stderr," Number of latitudes is %3d.",*klat); fprintf(stderr," Maximum is 320\n"); *kret = 2; goto L900; } if (*klon > 640) { fprintf(stderr," QU2REG :"); fprintf(stderr," Number of longitudes is %3d.",*klon); fprintf(stderr," Maximum is 640\n"); *kret = 3; goto L900; } /* Set array indices to 0. */ ilii = 0; ilio = 0; /* Establish values of loop parameters. */ if (*kcode > 10) { /* Quasi-regular along longitude lines. */ iquano = *klon; iregno = *klat; icode = *kcode - 10; } else { /* Quasi-regular along latitude lines. */ iquano = *klat; iregno = *klon; icode = *kcode; } /* -------------------------------------------------------- */ /** Section 2. Interpolate field from quasi to regular grid. */ /* -------------------------------------------------------- */ i_1 = iquano; for (j230 = 1; j230 <= i_1; ++j230) { if (iregno != kpoint[j230]) { /* Line contains less values than required,so */ /* extract quasi-regular grid values for a line */ i_2 = kpoint[j230]; for (j210 = 1; j210 <= i_2; ++j210) { ++ilii; zline[j210 - 1] = pfield[ilii]; } /* and interpolate this line. */ rowina2_(zline, &iregno, &kpoint[j230], zwork, &icode, pmsval, kret); if (*kret != 0) goto L900; /* Add regular grid values for this line to the temporary array. */ i_2 = iregno; for (j220 = 1; j220 <= i_2; ++j220) { ++ilio; ztemp[ilio - 1] = zline[j220 - 1]; } } else { /* Line contains the required number of values, so add */ /* this line to the temporary array. */ i_2 = iregno; for (j225 = 1; j225 <= i_2; ++j225) { ++ilio; ++ilii; ztemp[ilio - 1] = pfield[ilii]; } } } /* Copy temporary array to user array. */ i_1 = *klon * *klat; for (j240 = 1; j240 <= i_1; ++j240) { pfield[j240] = ztemp[j240 - 1]; } /* -------------------------------------------------------- */ /* Section 9. Return to calling routine. Format statements. */ /* -------------------------------------------------------- */ L900: return 0; } /* qu2reg2_ */ late ... */ } else { /* Adjust the weight to range (0.0 to 1.0) */ zwt -= ip; /* Interpolate using the weighted values on either side */ /* of the output point position */ p[jl] = ((REAL)1. - zwt) * pw[ip + 1 + pw_dim1] + zwt * pw[ip + 2 + pw_dim1]; } } /* ******************************* */ /* Section 2. Cubic interpolation .. */ /* ******************************* */ } else if (*kcode == 3) { i_1 = *ki; for (jl = 1; jl <= i_1; ++jl) { if (p[jl] == *pmsval) { fprintf(stderr," ROWINA2: "); fprintf(stderr," Cubic interpolation not supported"); fprintf(stderr," for fields containing missing data.\n"); *kret = 1; goto L900; } pw[jl + pw_dim1] = p[jl]; } pw[pw_dim1] = p[*ki]; pw[*ki + 1 + pw_dim1] = p[1]; pw[*ki + 2 + pw_dim1] = p[2]; i_1 = *ki; for (jl = 1; jl <= i_1; ++jl) { pw[jl + (pw_dim1 << 1)] = - pw[jl - 1 + pw_dim1] / (REAL)3. - pw[jl + pw_dim1] * (REAL).5 + pw[jl + 1 +pw_dim1] - pw[jl + 2 + pw_dim1] / (REAL)6.; pw[jl + 1 + pw_dim1 * 3] = pw[jl - 1 + pw_dim1] / (REAL)6. - pw[jl + pw_dim1] + pw[jl + 1 + pw_dim1] * (REAL).5 + pw[jl + 2 + pw_dim1] / (REAL)3.; } scm0_(&pw[(pw_dim1 << 1) + 1], &pw[pw_dim1 * 3 + 2], &pw[pw_dim1 + 1], &pw[pw_dim1 + 2], ki); zrdi = (REAL) (*ki); zdo = (REAL)1. / (REAL) (*ko); i_1 = *ko; for (jl = 1; jl <= i_1; ++jl) { zpos = (jl - 1) * zdo; zwt = zpos * zrdi; ip = zwt + 1; zwt = zwt + (REAL)1. - ip; zwt1 = (REAL)1. - zwt; p[jl] = (((REAL)3. - zwt1 * (REAL)2.) * pw[ip + pw_dim1] + zwt * pw[ip + (pw_dim1 << 1)]) * zwt1 * zwt1 + (((REAL)3. - zwt * (REAL)2.) * pw[ip + 1 + pw_dim1] - zwt1 * pw[ip + 1 + pw_dim1 * 3]) * zwt * zwt; } } else { /* ************************************** */ /* Section 3. Invalid interpolation code .. */ /* ************************************** */ fprintf(stderr," ROWINA2:"); fprintf(stderr," Invalid interpolation code = %2d\n",*kcode); *kret = 2; } L900: return 0; } /* rowina2_ */ int qu2reg2_(REAL *pfield, INT *kpoint, INT *klat, INT *klon, INT *kcode, REAL *pmsval, INT *kret) { /* System generated locals */ INT i_1, i_2; /* Local variables */ static INT ilii, ilio, icode; static REAL zline[640], ztemp[204800], zwork[1929] /* was [643][3] */; static INT iregno, iquano, j210, j220, j230, j240, j225; extern int rowina2_(); /* Parameter adjustments */ --pfield; --kpoint; /* **** QU2REG - Convert quasi-regular grid data to regular. */ /* Inpu