/***************************************************************************** * grib2api.c * * DESCRIPTION * This file contains the API to the GRIB2 libraries which is as close as * possible to the "official" NWS GRIB2 Library's API. The reason for this * is so we can use NCEP's unofficial GRIB2 library while having minimal * impact on the already written drivers. * * HISTORY * 12/2003 Arthur Taylor (MDL / RSIS): Created. * * NOTES ***************************************************************************** */ #include #include #include #include "libaat.h" #include "mdl_g2c.h" #include "grib2.h" #ifdef MEMWATCH #include "memwatch.h" #endif /* Would prefer to include "gribtemplates.h" but can't since NCEP decided to * put a const struct in it. The result is if I include it, I get 2 copies * of the data. */ extern sInt4 getgridindex(sInt4 number); /* have now put this in grib2api.h... may bring it back. */ #define MAXGRIDTEMP 23 /* maximum number of templates */ #define MAXGRIDMAPLEN 200 /* maximum template map length */ struct gridtemplate { sInt4 template_num; sInt4 mapgridlen; sInt4 needext; sInt4 mapgrid[MAXGRIDMAPLEN]; }; extern const struct gridtemplate templatesgrid[MAXGRIDTEMP]; /* Would prefer to include "pdstemplates.h" but can't since NCEP decided to * put a const struct in it. The result is if I include it, I get 2 copies * of the data. */ extern sInt4 getpdsindex(sInt4 number); /* have now put this in grib2api.h... may bring it back. */ #define MAXPDSTEMP 23 /* maximum number of templates */ #define MAXPDSMAPLEN 200 /* maximum template map length */ struct pdstemplate { sInt4 template_num; sInt4 mappdslen; sInt4 needext; sInt4 mappds[MAXPDSMAPLEN]; }; extern const struct pdstemplate templatespds[MAXPDSTEMP]; /* Would prefer to include "drstemplates.h" but can't since NCEP decided to * put a const struct in it. The result is if I include it, I get 2 copies * of the data. */ extern sInt4 getdrsindex(sInt4 number); #define MAXDRSTEMP 8 /* maximum number of templates */ #define MAXDRSMAPLEN 200 /* maximum template map length */ struct drstemplate { sInt4 template_num; sInt4 mapdrslen; sInt4 needext; sInt4 mapdrs[MAXDRSMAPLEN]; }; extern const struct drstemplate templatesdrs[MAXDRSTEMP]; /***************************************************************************** * mdl_LocalUnpack() -- * * Arthur Taylor / MDL * * PURPOSE * Unpack the local use data assuming that it was packed using the MDL * encoder. This assumes "local" starts at octet 6 (ie skipping over the * length and section ID octets) * * In Section 2, GRIB2 provides for local use data. The MDL encoder packs * that data, and signifies that it has done so by setting octet 6 to 1. * This may have been inadvisable, since the GRIB2 specs did not state that * octet 6 was special. * * ARGUMENTS * local = The section 2 data prior to being unpacked. (Input) * locallen = The length of "local" (Input) * idat = Unpacked MDL local data (if it was an integer) (Output) * nidat = length of idat. (Input) * rdat = Unpacked MDL local data (if it was a float) (Output) * nrdat = length of rdat. (Input) * * FILES/DATABASES: None * * RETURNS: int * 0 = Ok. * 1 = Mixed local data types? * 2 = nrdat is not large enough. * 3 = nidat is not large enough. * 4 = Too many bits to unpack, should be a max of 32 bits. * 5 = Locallen is too small * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES ***************************************************************************** */ #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b) static int mdl_LocalUnpack(unsigned char *local, sInt4 locallen, sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat) { int BytesUsed = 0; /* How many bytes have been used. */ unsigned int numGroup; /* Number of groups */ int i; /* Counter over the groups. */ sInt4 numVal; /* Number of values in a given group. */ float refVal; /* The reference value in a given group. */ unsigned int scale; /* The power of 10 scale value in the group. */ sInt4 recScale10; /* 1 / 10**scale.. For faster computations. */ unsigned char numBits; /* # of bits for a single element in the group. */ char f_dataType; /* If the local data is a float or integer. */ char f_firstType = 0; /* The type of the first group of local data. */ int curIndex = 0; /* Where to store the current data. */ int j; /* Counter over the number of values in a group. */ size_t numUsed; /* Number of bytes used in call to memBitRead. */ uChar bufLoc; /* Where to read for more bits in the data stream. */ uInt4 uli_temp; /* Temporary storage to hold unpacked data. */ if (locallen < BytesUsed + 3) { #ifdef DEBUG printf("Locallen is too small.\n"); #endif return 5; } /* The calling routine should check octet 6, which is local[0], to be 1, * so we just assert it is 1. */ myAssert(local[0] == 1); numGroup = GRIB_UNSIGN_INT2(local[1], local[2]); local += 3; BytesUsed += 3; myAssert(*nrdat > 1); myAssert(*nidat > 1); idat[0] = 0; rdat[0] = 0; for (i = 0; i < numGroup; i++) { if (locallen < BytesUsed + 12) { #ifdef DEBUG printf("Locallen is too small.\n"); #endif return 5; } MEMCPY_BIG(&numVal, local, sizeof(sInt4)); MEMCPY_BIG(&refVal, local + 4, sizeof(float)); scale = GRIB_UNSIGN_INT2(local[8], local[9]); recScale10 = 1 / pow(10, scale); numBits = local[10]; if (numBits >= 32) { #ifdef DEBUG printf("Too many bits too unpack.\n"); #endif return 4; } f_dataType = local[11]; local += 12; BytesUsed += 12; if (locallen < BytesUsed + ((numBits * numVal) + 7) / 8) { #ifdef DEBUG printf("Locallen is too small.\n"); #endif return 5; } if (i == 0) { f_firstType = f_dataType; } else if (f_firstType != f_dataType) { #ifdef DEBUG printf("Local use data has mixed float/integer type?\n"); #endif return 1; } bufLoc = 8; if (f_dataType == 0) { /* Floating point data. */ if (*nrdat < (curIndex + numVal + 3)) { #ifdef DEBUG printf("nrdat is not large enough.\n"); #endif return 2; } rdat[curIndex] = numVal; curIndex++; rdat[curIndex] = scale; curIndex++; for (j = 0; j < numVal; j++) { memBitRead(&uli_temp, sizeof(sInt4), local, numBits, &bufLoc, &numUsed); local += numUsed; BytesUsed += numUsed; rdat[curIndex] = (refVal + uli_temp) * recScale10; curIndex++; } rdat[curIndex] = 0; } else { /* Integer point data. */ if (*nidat < (curIndex + numVal + 3)) { #ifdef DEBUG printf("nidat is not large enough.\n"); #endif return 3; } idat[curIndex] = numVal; curIndex++; idat[curIndex] = scale; curIndex++; for (j = 0; j < numVal; j++) { memBitRead(&uli_temp, sizeof(sInt4), local, numBits, &bufLoc, &numUsed); local += numUsed; BytesUsed += numUsed; idat[curIndex] = (refVal + uli_temp) * recScale10; curIndex++; } idat[curIndex] = 0; } } return 0; } /***************************************************************************** * fillOutSectLen() -- * * Arthur Taylor / MDL * * PURPOSE * The GRIB2 API returns the lengths of each GRIB2 section along with the * meta data. NCEP's routines did not, so this routine fills in that part. * This routine has to consider which grid is being unpacked. * * c_ipack is passed in after section 1 (since section 1 is not repeated) * * ARGUMENTS * c_ipack = Complete GRIB2 message to look for the section lengths in. (In) * lenCpack = Length of c_ipack (Input) * subgNum = Which sub rid to find the section lengths of. (Input) * is2 = Section 2 data (Output) * is3 = Section 3 data (Output) * is4 = Section 4 data (Output) * is5 = Section 5 data (Output) * is6 = Section 6 data (Output) * is7 = Section 7 data (Output) * * FILES/DATABASES: None * * RETURNS: int * 0 = Ok. * 1 = c_ipack is not large enough * 2 = invalid section number * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES ***************************************************************************** */ static int fillOutSectLen(unsigned char *c_ipack, int lenCpack, int subgNum, sInt4 *is2, sInt4 *is3, sInt4 *is4, sInt4 *is5, sInt4 *is6, sInt4 *is7) { sInt4 offset = 0; /* How far in c_ipack we have read. */ sInt4 sectLen; /* The length of the current section. */ int sectId; /* The current section number. */ unsigned int gNum = 0; /* Which sub group we are currently working with. */ if (lenCpack < 5) { #ifdef DEBUG printf("Cpack is not large enough.\n"); #endif return 1; } /* assert that we start with data in either section 2 or 3. */ myAssert((c_ipack[4] == 2) || (c_ipack[4] == 3)); while (gNum <= subgNum) { if (lenCpack < offset + 5) { #ifdef DEBUG printf("Cpack is not large enough.\n"); #endif return 1; } MEMCPY_BIG(§Len, c_ipack + offset, sizeof(sInt4)); /* Check if we just read section 8. If so, then it is "7777" = * 926365495 regardless of endian'ness. */ if (sectLen == 926365495L) { #ifdef DEBUG printf("Shouldn't see sect 8. Should stop after correct sect 7\n"); #endif return 2; } sectId = c_ipack[offset + 4]; switch (sectId) { case 2: is2[0] = sectLen; break; case 3: is3[0] = sectLen; break; case 4: is4[0] = sectLen; break; case 5: is5[0] = sectLen; break; case 6: is6[0] = sectLen; break; case 7: is7[0] = sectLen; gNum++; break; default: #ifdef DEBUG printf("Invalid section id %d.\n", sectId); #endif return 2; } offset += sectLen; } return 0; } /***************************************************************************** * TransferInt() -- * * Arthur Taylor / MDL * * PURPOSE * To transfer the data from "fld" and "bmap" to "iain" and "ib". The API * attempts to rearrange the order, so that anything returned from it has * scan mode 0100???? * * ARGUMENTS * fld = The expanded grid from NCEPs routines (Input) * ngrdpts = Length of fld (Input) * ibitmap = flag if we have a bitmap or not. (Input) * bmap = bitmap from NCEPs routines. (Input) * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input) * scan = The scan orientation of fld/bmap/iain/ib (Input/Output) * nx, ny = The dimmensions of the grid. (Input) * iclean = 1 means the user wants the unpacked data returned without * missing values in it. 0 means embed the missing values. (In) * xmissp = The primary missing value to use if iclean = 0. (Input). * iain = The grid to return. (Output) * nd2x3 = The length of iain. (Input) * ib = Bitmap if it was packed (Output) * * FILES/DATABASES: None * * RETURNS: int * 0 = Ok. * 1 = nd2x3 is too small. * 2 = nx*ny != ngrdpts * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES * May want to disable the scan adjustment in the future. ***************************************************************************** */ static int TransferInt(float *fld, sInt4 ngrdpts, sInt4 ibitmap, sInt4 *bmap, char f_ignoreScan, sInt4 *scan, sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp, sInt4 *iain, sInt4 nd2x3, sInt4 *ib) { int i; /* loop counter over all grid points. */ sInt4 x, y; /* Where we are in a grid of scan value 0100???? */ int curIndex; /* Where in iain to store the current data. */ if (nd2x3 < ngrdpts) { #ifdef DEBUG printf("nd2x3(%ld) is < ngrdpts(%ld)\n", nd2x3, ngrdpts); #endif return 1; } if (f_ignoreScan || ((*scan & 0xf0) == 64)) { if (ibitmap) { for (i = 0; i < ngrdpts; i++) { ib[i] = bmap[i]; /* Check if we are supposed to insert xmissp into the field */ if ((iclean != 0) && (ib[i] == 0)) { iain[i] = xmissp; } else { iain[i] = fld[i]; } } } else { for (i = 0; i < ngrdpts; i++) { iain[i] = fld[i]; } } } else { if (nx * ny != ngrdpts) { #ifdef DEBUG printf("nx * ny (%ld) != ngrdpts(%ld)\n", nx * ny, ngrdpts); #endif return 2; } if (ibitmap) { for (i = 0; i < ngrdpts; i++) { ScanIndex2XY(i, &x, &y, *scan, nx, ny); /* ScanIndex returns value as if scan was 0100(0000) */ curIndex = (x - 1) + (y - 1) * nx; myAssert(curIndex < nd2x3); ib[curIndex] = bmap[i]; /* Check if we are supposed to insert xmissp into the field */ if ((iclean != 0) && (ib[curIndex] == 0)) { iain[i] = xmissp; } else { iain[curIndex] = fld[i]; } } } else { for (i = 0; i < ngrdpts; i++) { ScanIndex2XY(i, &x, &y, *scan, nx, ny); /* ScanIndex returns value as if scan was 0100(0000) */ curIndex = (x - 1) + (y - 1) * nx; myAssert(curIndex < nd2x3); iain[curIndex] = fld[i]; } } *scan = 64 + (*scan & 0x0f); } return 0; } /***************************************************************************** * TransferFloat() -- * * Arthur Taylor / MDL * * PURPOSE * To transfer the data from "fld" and "bmap" to "ain" and "ib". The API * attempts to rearrange the order, so that anything returned from it has * scan mode 0100???? * * ARGUMENTS * fld = The expanded grid from NCEPs routines (Input) * ngrdpts = Length of fld (Input) * ibitmap = flag if we have a bitmap or not. (Input) * bmap = bitmap from NCEPs routines. (Input) * scan = The scan orientation of fld/bmap/iain/ib (Input/Output) * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input) * nx, ny = The dimmensions of the grid. (Input) * iclean = 1 means the user wants the unpacked data returned without * missing values in it. 0 means embed the missing values. (In) * xmissp = The primary missing value to use if iclean = 0. (Input). * ain = The grid to return. (Output) * nd2x3 = The length of iain. (Input) * ib = Bitmap if it was packed (Output) * * FILES/DATABASES: None * * RETURNS: int * 0 = Ok. * 1 = nd2x3 is too small. * 2 = nx*ny != ngrdpts * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES * May want to disable the scan adjustment in the future. ***************************************************************************** */ static int TransferFloat(float *fld, sInt4 ngrdpts, sInt4 ibitmap, sInt4 *bmap, char f_ignoreScan, sInt4 *scan, sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp, float *ain, sInt4 nd2x3, sInt4 *ib) { int i; /* loop counter over all grid points. */ sInt4 x, y; /* Where we are in a grid of scan value 0100???? */ int curIndex; /* Where in ain to store the current data. */ if (nd2x3 < ngrdpts) { #ifdef DEBUG printf("nd2x3(%ld) is < ngrdpts(%ld)\n", nd2x3, ngrdpts); #endif return 1; } if (f_ignoreScan || ((*scan & 0xf0) == 64)) { if (ibitmap) { for (i = 0; i < ngrdpts; i++) { ib[i] = bmap[i]; /* Check if we are supposed to insert xmissp into the field */ if ((iclean != 0) && (ib[i] == 0)) { ain[i] = xmissp; } else { ain[i] = fld[i]; } } } else { for (i = 0; i < ngrdpts; i++) { ain[i] = fld[i]; } } } else { if (nx * ny != ngrdpts) { #ifdef DEBUG printf("nx * ny (%ld) != ngrdpts(%ld)\n", nx * ny, ngrdpts); #endif return 2; } if (ibitmap) { for (i = 0; i < ngrdpts; i++) { ScanIndex2XY(i, &x, &y, *scan, nx, ny); /* ScanIndex returns value as if scan was 0100(0000) */ curIndex = (x - 1) + (y - 1) * nx; myAssert(curIndex < nd2x3); ib[curIndex] = bmap[i]; /* Check if we are supposed to insert xmissp into the field */ if ((iclean != 0) && (ib[curIndex] == 0)) { ain[i] = xmissp; } else { ain[curIndex] = fld[i]; } } } else { for (i = 0; i < ngrdpts; i++) { ScanIndex2XY(i, &x, &y, *scan, nx, ny); /* ScanIndex returns value as if scan was 0100(0000) */ curIndex = (x - 1) + (y - 1) * nx; myAssert(curIndex < nd2x3); ain[curIndex] = fld[i]; } } *scan = 64 + (*scan & 0x0f); } /* if (1==1) { int i; for (i=0; i < ngrdpts; i++) { if (ain[i] != 9999) { fprintf (stderr, "grib2api.c: ain[%d] %f, fld[%d] %f\n", i, ain[i], i, fld[i]); } } } */ return 0; } #ifdef PKNCEP int pk_g2ncep(sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nx, sInt4 *ny, sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat, sInt4 *is0, sInt4 *ns0, sInt4 *is1, sInt4 *ns1, sInt4 *is3, sInt4 *ns3, sInt4 *is4, sInt4 *ns4, sInt4 *is5, sInt4 *ns5, sInt4 *is6, sInt4 *ns6, sInt4 *is7, sInt4 *ns7, sInt4 *ib, sInt4 *ibitmap, unsigned char *cgrib, sInt4 *nd5, sInt4 *missp, float *xmissp, sInt4 *misss, float *xmisss, sInt4 *inew, sInt4 *minpk, sInt4 *iclean, sInt4 *l3264b, sInt4 *jer, sInt4 *ndjer, sInt4 *kjer) { g2int listsec0[2]; g2int listsec1[13]; g2int igds[5]; g2int igdstmpl[]; int ierr; /* Holds the error code from a called routine. */ int i; listsec0[0] = is0[6]; listsec0[1] = is0[7]; listsec1[0] = is1[5]; listsec1[1] = is1[7]; listsec1[2] = is1[9]; listsec1[3] = is1[10]; listsec1[4] = is1[11]; listsec1[5] = is1[12]; listsec1[6] = is1[14]; listsec1[7] = is1[15]; listsec1[8] = is1[16]; listsec1[9] = is1[17]; listsec1[10] = is1[18]; listsec1[11] = is1[19]; listsec1[12] = is1[20]; ierr = g2_create(cgrib, listsec0, listsec1); printf("Length = %d\n", ierr); if ((idat[0] != 0) || (rdat[0] != 0)) { printf("Don't handle this yet.\n"); /* ierr = g2_addlocal (cgrib, unsigned char *csec2, g2int lcsec2); */ } igds[0] = is3[5]; igds[1] = is3[6]; igds[2] = is3[10]; igds[3] = is3[11]; igds[4] = is3[12]; IS3(15) - IS3(nn) = Grid Definition Template, stored in bytes 15 - nn(*) ierr = g2_addgrid(cgrib, igds, g2int * igdstmpl, g2int * ideflist, g2int idefnum) return 0; /* To start a new GRIB2 message, call function g2_create. G2_create encodes Sections 0 and 1 at the beginning of the message. This routine must be used to create each message. Routine g2_addlocal can be used to add a Local Use Section ( Section 2 ). Note that this section is optional and need not appear in a GRIB2 message. Function g2_addgrid is used to encode a grid definition into Section 3. This grid definition defines the geometry of the the data values in the fields that follow it. g2_addgrid can be called again to change the grid definition describing subsequent data fields. Each data field is added to the GRIB2 message using routine g2_addfield, which adds Sections 4, 5, 6, and 7 to the message. After all desired data fields have been added to the GRIB2 message, a call to function g2_gribend is needed to add the final section 8 to the message and to update the length of the message. A call to g2_gribend is required for each GRIB2 message. */ } #endif /***************************************************************************** * unpk_g2ncep() -- * * Arthur Taylor / MDL * * PURPOSE * This procedure is a wrapper around the NCEP GRIB2 routines to interface * their routines with the "official NWS" GRIB2 API. The reason for this is * so drivers that have been written to use the "official NWS" GRIB2 API, can * use the NCEP library with minimal disruption. * * ARGUMENTS * kfildo = Unit number for output diagnostics (C ignores this). (Input) * ain = contains the data if the original data was float data. * (size = nd2x3) (Output) * iain = contains the data if the original data was integer data. * (size = nd2x3) (Output) * nd2x3 = length of ain, iain, ib (Input) * (at least the size of num grid points) * idat = local use data if any that were unpacked from Section 2. (Output) * nidat = length of idat. (Input) * rdat = floating point local use data (Output) * nrdat = length of rdat. (Input) * is0 = Section 0 data (Output) * ns0 = length of is0 (16 is fine) (Input) * is1 = Section 1 data (Output) * ns1 = length of is1 (21 is fine) (Input) * is2 = Section 2 data (Output) * ns2 = length of is2 () (Input) * is3 = Section 3 data (Output) * ns3 = length of is3 (96 or 1600) (Input) * is4 = Section 4 data (Output) * ns4 = length of is4 (60) (Input) * is5 = Section 5 data (Output) * ns5 = length of is5 (49 is fine) (Input) * is6 = Section 6 data (Output) * ns6 = length of is6 (6 is fine) (Input) * is7 = Section 7 data (Output) * ns7 = length of is7 (8 is fine) (Input) * ib = Bitmap if user requested it, and it was packed (Output) * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output) * c_ipack = The message to unpack (Input) * nd5 = 1/4 the size of c_ipack (Input) * xmissp = The floating point representation for the primary missing value. * (Input/Output) * xmisss = The floating point representation for the secondary missing * value (Input/Output) * new = 1 if this is the first grid to be unpacked, else 0. (Input) * iclean = 1 means the user wants the unpacked data returned without * missing values in it. 0 means embed the missing values. (Input) * l3264b = Integer word length in bits (32 or 64) (Input) * iendpk = A 1 means no more grids in this message, a 0 means more grids. * (Output) * jer(ndjer,2) = error codes along with severity. (Output) * ndjer = 1/2 length of jer. (>= 15) (Input) * kjer = number of error messages stored in jer. * * FILES/DATABASES: None * * RETURNS: void * "jer" is used to store error messages, and kjer is used to denote how * many there were. Jer is set up as a 2 column array, with the first * column in the first half of the array, and denoting the GRIB2 section * an error occurred in. The second column denoting the severity. * * The possibilities from unpk_g2ncep() are as follows: * ker=1 jer[0,0]=0 jer[0,1]=2: Message is not formed correctly * Request for an invalid subgrid * Problems unpacking the data. * problems expanding the data. * Calling dimmensions were too small. * ker=2 jer[1,0]=100 jer[1,1]=2: Error unpacking section 1. * ker=3 jer[2,0]=200 jer[2,1]=2: Error unpacking section 2. * ker=4 jer[3,0]=300 jer[3,1]=2: Error unpacking section 3. * ker=5 jer[4,0]=400 jer[4,1]=2: Error unpacking section 4. * ker=6 jer[5,0]=500 jer[5,1]=2: Error unpacking section 5. * Data Template not implemented. * Durring Transfer, nx * ny != ngrdpts. * ker=7 jer[6,0]=600 jer[6,1]=2: Error unpacking section 6. * ker=8 jer[7,0]=700 jer[7,1]=2: Error unpacking section 7. * ker=9 jer[8,0]=2001 jer[8,1]=2: nd2x3 is not large enough. * ker=9 jer[8,0]=2003 jer[8,1]=2: undefined sect 3 template * (see gridtemplates.h). * ker=9 jer[8,0]=2004 jer[8,1]=2: undefined sect 4 template * (see pdstemplates.h). * ker=9 jer[8,0]=2005 jer[8,1]=2: undefined sect 5 template * (see drstemplates.h). * ker=9 jer[8,0]=9999 jer[8,1]=2: NCEP returns an unrecognized error. * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner. * MDL attempts to always return grids in scan mode 0100???? * * ToDo: Check length of parameters better. * ToDo: Probably shouldn't abort if they have problems expanding the data. * ToDo: Better handling of: * gfld->list_opt = (Used if gfld->numoct_opt .ne. 0) This array contains * the number of grid points contained in each row (or column). (part of * Section 3) This element is a pointer to an array that holds the data. * This pointer is nullified if gfld->numoct_opt=0. * gfld->num_opt = (Used if gfld->numoct_opt .ne. 0) The number of entries in * array ideflist. i.e. number of rows (or columns) for which optional * grid points are defined. This value is set to zero, if * gfld->numoct_opt=0. * gfld->coord_list = Real array containing floating point values intended to * document the vertical discretisation associated to model data on * hybrid coordinate vertical levels. (part of Section 4) This element * is a pointer to an array that holds the data. * gfld->num_coord = number of values in array gfld->coord_list[]. ***************************************************************************** */ void unpk_g2ncep(sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nd2x3, sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat, sInt4 *is0, sInt4 *ns0, sInt4 *is1, sInt4 *ns1, sInt4 *is2, sInt4 *ns2, sInt4 *is3, sInt4 *ns3, sInt4 *is4, sInt4 *ns4, sInt4 *is5, sInt4 *ns5, sInt4 *is6, sInt4 *ns6, sInt4 *is7, sInt4 *ns7, sInt4 *ib, sInt4 *ibitmap, unsigned char *c_ipack, sInt4 *nd5, float *xmissp, float *xmisss, sInt4 *inew, sInt4 *iclean, sInt4 *l3264b, sInt4 *iendpk, sInt4 *jer, sInt4 *ndjer, sInt4 *kjer) { int i; /* A counter used for a number of purposes. */ static unsigned int subgNum = 0; /* The sub grid we read most recently. * This is primarily to help with the * inew option. */ int ierr; /* Holds the error code from a called routine. */ sInt4 listsec0[3]; sInt4 listsec1[13]; static sInt4 numfields = 1; /* Number of sub Grids in this message */ sInt4 numlocal; /* Number of local sections in this message. */ int unpack; /* Tell g2_getfld to unpack the message. */ int expand; /* Tell g2_getflt to attempt to expand the bitmap. */ gribfield *gfld; /* Holds the data after g2_getfld unpacks it. */ sInt4 gridIndex; /* index in templatesgrid[] for this sect 3 templat */ sInt4 pdsIndex; /* index in templatespds[] for this sect 4 template */ sInt4 drsIndex; /* index in templatesdrs[] for this sect 5 template */ int curIndex; /* Where in is3, is4, or is5 to store meta data */ int scanIndex; /* Where in is3 to find the scan mode. */ int nxIndex; /* Where in is3 to find the number of x values. */ int nyIndex; /* Where in is3 to find the number of y values. */ float f_temp; /* Assist with handling weird MDL behavior in is5[] */ char f_ignoreScan; /* Flag to ignore the attempt at changing the scan */ sInt4 dummyScan; /* Dummy place holder for call to Transfer routines * if ignoring scan. */ myAssert(*ndjer >= 8); /* Init the error handling array. */ memset((void *)jer, 0, 2 * *ndjer * sizeof(sInt4)); for (i = 0; i < 8; i++) { jer[i] = i * 100; } *kjer = 8; /* The first time in, figure out how many grids there are, and store it in * numfields for subsequent calls with inew != 1. */ if (*inew == 1) { subgNum = 0; ierr = g2_info(c_ipack, listsec0, listsec1, &numfields, &numlocal); if (ierr != 0) { switch (ierr) { case 1: /* Beginning characters "GRIB" not found. */ case 2: /* GRIB message is not Edition 2. */ case 3: /* Could not find Section 1, where expected. */ case 4: /* End string "7777" found, but not where expected. */ case 5: /* End string "7777" not found at end of message. */ case 6: /* Invalid section number found. */ jer[0 + *ndjer] = 2; *kjer = 1; break; default: jer[8 + *ndjer] = 2; jer[8] = 9999; /* Unknown error message. */ *kjer = 9; } return; } } else { if (subgNum + 1 >= numfields) { /* Field request error. */ jer[0 + *ndjer] = 2; *kjer = 1; return; } subgNum++; } /* Expand the desired subgrid. */ unpack = 1; expand = 1; ierr = g2_getfld(c_ipack, subgNum + 1, unpack, expand, &gfld); if (ierr != 0) { switch (ierr) { case 1: /* Beginning characters "GRIB" not found. */ case 2: /* GRIB message is not Edition 2. */ case 3: /* The data field request number was not positive. */ case 4: /* End string "7777" found, but not where expected. */ case 6: /* message did not contain requested # of fields */ case 7: /* End string "7777" not found at end of message. */ case 8: /* Unrecognized Section encountered. */ jer[0 + *ndjer] = 2; *kjer = 1; break; case 15: /* Error unpacking Section 1. */ jer[1 + *ndjer] = 2; *kjer = 2; break; case 16: /* Error unpacking Section 2. */ jer[2 + *ndjer] = 2; *kjer = 3; break; case 10: /* Error unpacking Section 3. */ jer[3 + *ndjer] = 2; *kjer = 4; break; case 11: /* Error unpacking Section 4. */ jer[4 + *ndjer] = 2; *kjer = 5; break; case 9: /* Data Template 5.NN not implemented. */ case 12: /* Error unpacking Section 5. */ jer[5 + *ndjer] = 2; *kjer = 6; break; case 13: /* Error unpacking Section 6. */ jer[6 + *ndjer] = 2; *kjer = 7; break; case 14: /* Error unpacking Section 7. */ jer[7 + *ndjer] = 2; *kjer = 8; break; default: jer[8 + *ndjer] = 2; jer[8] = 9999; /* Unknown error message. */ *kjer = 9; break; } g2_free(gfld); return; } /* Check if data wasn't unpacked. */ if (!gfld->unpacked) { jer[0 + *ndjer] = 2; *kjer = 1; g2_free(gfld); return; } /* Start going through the gfld structure and converting it to the needed * data output formats. */ myAssert(*ns0 >= 16); MEMCPY_BIG(&(is0[0]), c_ipack, sizeof(sInt4)); is0[6] = gfld->discipline; is0[7] = gfld->version; MEMCPY_BIG(&(is0[8]), c_ipack + 8, sizeof(sInt4)); /* The following assert fails only if the GRIB message is more that 4 * giga-bytes large, which I think would break the fortran library. */ myAssert(is0[8] == 0); MEMCPY_BIG(&(is0[8]), c_ipack + 12, sizeof(sInt4)); myAssert(*ns1 >= 21); myAssert(gfld->idsectlen >= 13); MEMCPY_BIG(&(is1[0]), c_ipack + 16, sizeof(sInt4)); is1[4] = c_ipack[20]; is1[5] = gfld->idsect[0]; is1[7] = gfld->idsect[1]; is1[9] = gfld->idsect[2]; is1[10] = gfld->idsect[3]; is1[11] = gfld->idsect[4]; is1[12] = gfld->idsect[5]; /* Year */ is1[14] = gfld->idsect[6]; /* Month */ is1[15] = gfld->idsect[7]; /* Day */ is1[16] = gfld->idsect[8]; /* Hour */ is1[17] = gfld->idsect[9]; /* Min */ is1[18] = gfld->idsect[10]; /* Sec */ is1[19] = gfld->idsect[11]; is1[20] = gfld->idsect[12]; /* Fill out section lengths (separate procedure because of possibility of * having multiple grids. Should combine fillOutSectLen g2_info, and * g2_getfld into one procedure to optimize it. */ fillOutSectLen(c_ipack + 16 + is1[0], 4 * *nd5 - 15 - is1[0], subgNum, is2, is3, is4, is5, is6, is7); /* Check if there is section 2 data. */ if (gfld->locallen > 0) { /* The + 1 is so we don't overwrite the section length */ memset((void *)(is2 + 1), 0, (*ns2 - 1) * sizeof(sInt4)); is2[4] = 2; is2[5] = gfld->local[0]; /* check if MDL Local use simple packed data */ if (is2[5] == 1) { mdl_LocalUnpack(gfld->local, gfld->locallen, idat, nidat, rdat, nrdat); } else { /* local use section was not MDL packed, return it in is2. */ for (i = 0; i < gfld->locallen; i++) { is2[i + 5] = gfld->local[i]; } } } else { /* API specified that is2[0] = 0; idat[0] = 0, and rdat[0] = 0 */ is2[0] = 0; idat[0] = 0; rdat[0] = 0; } is3[4] = 3; is3[5] = gfld->griddef; is3[6] = gfld->ngrdpts; if (*nd2x3 < gfld->ngrdpts) { jer[8 + *ndjer] = 2; jer[8] = 2001; /* nd2x3 is not large enough */ *kjer = 9; g2_free(gfld); return; } is3[10] = gfld->numoct_opt; is3[11] = gfld->interp_opt; is3[12] = gfld->igdtnum; gridIndex = getgridindex(gfld->igdtnum); if (gridIndex == -1) { jer[8 + *ndjer] = 2; jer[8] = 2003; /* undefined sect 3 template */ *kjer = 9; g2_free(gfld); return; } curIndex = 14; for (i = 0; i < gfld->igdtlen; i++) { is3[curIndex] = gfld->igdtmpl[i]; curIndex += abs(templatesgrid[gridIndex].mapgrid[i]); } /* API attempts to return grid in scan mode 0100????. Find the necessary * indexes into the is3 array for the attempt. */ switch (gfld->igdtnum) { case 0: case 1: case 2: case 3: case 40: case 41: case 42: case 43: scanIndex = 72 - 1; nxIndex = 31 - 1; nyIndex = 35 - 1; break; case 10: scanIndex = 60 - 1; nxIndex = 31 - 1; nyIndex = 35 - 1; break; case 20: case 30: case 31: scanIndex = 65 - 1; nxIndex = 31 - 1; nyIndex = 35 - 1; break; case 90: scanIndex = 64 - 1; nxIndex = 31 - 1; nyIndex = 35 - 1; break; case 110: scanIndex = 57 - 1; nxIndex = 31 - 1; nyIndex = 35 - 1; break; case 50: case 51: case 52: case 53: case 100: case 120: case 1000: case 1200: default: scanIndex = -1; nxIndex = -1; nyIndex = -1; } is4[4] = 4; is4[5] = gfld->num_coord; is4[7] = gfld->ipdtnum; pdsIndex = getpdsindex(gfld->ipdtnum); if (pdsIndex == -1) { jer[8 + *ndjer] = 2; jer[8] = 2004; /* undefined sect 4 template */ *kjer = 9; g2_free(gfld); return; } curIndex = 9; for (i = 0; i < gfld->ipdtlen; i++) { is4[curIndex] = gfld->ipdtmpl[i]; curIndex += abs(templatespds[pdsIndex].mappds[i]); } is5[4] = 5; is5[5] = gfld->ndpts; is5[9] = gfld->idrtnum; drsIndex = getdrsindex(gfld->idrtnum); if (drsIndex == -1) { jer[8 + *ndjer] = 2; jer[8] = 2005; /* undefined sect 5 template */ *kjer = 9; g2_free(gfld); return; } curIndex = 11; for (i = 0; i < gfld->idrtlen; i++) { is5[curIndex] = gfld->idrtmpl[i]; curIndex += abs(templatesdrs[drsIndex].mapdrs[i]); } /* Mimic MDL's handling of reference value (IS5(12)) */ memcpy(&f_temp, &(is5[11]), sizeof(float)); is5[11] = (sInt4)f_temp; if ((is5[9] == 2) || (is5[9] == 3)) { if (is5[20] == 0) { memcpy(&(f_temp), &(is5[23]), sizeof(float)); *xmissp = f_temp; is5[23] = (sInt4)f_temp; memcpy(&(f_temp), &(is5[27]), sizeof(float)); *xmisss = f_temp; is5[27] = (sInt4)f_temp; } else { *xmissp = is5[23]; *xmisss = is5[27]; } } is6[4] = 6; is6[5] = gfld->ibmap; is7[4] = 7; if (subgNum + 1 == numfields) { *iendpk = 1; } else { *iendpk = 0; } if ((gfld->ibmap == 0) || (gfld->ibmap == 254)) { *ibitmap = 1; } else { *ibitmap = 0; } /* Check type of original field, before transfering the memory. */ myAssert(*ns5 > 20); /* Check if NCEP had problems expanding the data. If so we currently * abort. May need to revisit this behavior. */ if (!gfld->expanded) { jer[0 + *ndjer] = 2; *kjer = 1; g2_free(gfld); return; } f_ignoreScan = 0; /* Check if integer type... is5[20] == 1 implies integer (code table 5.1), * but only for certain templates. (0,1,2,3,40,41,40000,40010). (not 50,51) */ if ((is5[20] == 1) && ((is5[9] != 50) && (is5[9] != 51))) { /* integer data, use iain */ if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) { ierr = TransferInt(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap, 1, &(dummyScan), 0, 0, *iclean, *xmissp, iain, *nd2x3, ib); } else { ierr = TransferInt(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap, f_ignoreScan, &(is3[scanIndex]), is3[nxIndex], is3[nyIndex], *iclean, *xmissp, iain, *nd2x3, ib); } } else { /* float data, use ain */ if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) { ierr = TransferFloat(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap, 1, &(dummyScan), 0, 0, *iclean, *xmissp, ain, *nd2x3, ib); } else { ierr = TransferFloat(gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap, f_ignoreScan, &(is3[scanIndex]), is3[nxIndex], is3[nyIndex], *iclean, *xmissp, ain, *nd2x3, ib); } } if (ierr != 0) { switch (ierr) { case 1: /* nd2x3 is too small */ jer[0 + *ndjer] = 2; *kjer = 1; break; case 2: /* nx * ny != ngrdpts */ jer[5 + *ndjer] = 2; *kjer = 6; break; default: jer[8 + *ndjer] = 2; jer[8] = 9999; /* Unknown error message. */ *kjer = 9; } g2_free(gfld); return; } g2_free(gfld); } /***************************************************************************** * BigByteCpy() -- * * Arthur Taylor / MDL * * PURPOSE * This is so we can copy upto 4 bytes from a big endian 4 byte int data * stream. * * The reason this is needed is because the GRIB2 API required the GRIB2 * message to be passed in as a big endian 4 byte int data stream instead of * something more reasonable (such as a stream of 1 byte char's) By having * this routine we reduce the number of memswaps of the message on a little * endian system. * * ARGUMENTS * dst = Where to copy the data to. (Output) * ipack = The 4 byte int data stream. (Input) * nd5 = length of ipack (Input) * startInt = Which integer to start reading in ipack. (Input) * startByte = Which byte in the startInt to start reading from. (Input) * numByte = How many bytes to read. (Input) * * FILES/DATABASES: None * * RETURNS: NULL * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES ***************************************************************************** */ static void BigByteCpy(sInt4 *dst, sInt4 *ipack, sInt4 nd5, unsigned int startInt, unsigned int startByte, int numByte) { static char Lshift[] = { 0, 8, 16, 24 }; /* Amounts to shift left by. */ unsigned int intIndex; /* Where in ipack to read from. */ unsigned int byteIndex; /* Where in intIndex to read from. */ uInt4 curInt; /* An unsigned version of the current int. */ unsigned int curByte; /* The current byte we have read. */ int i; /* Loop counter over number of bytes to read. */ myAssert(numByte <= 4); myAssert(startByte < 4); *dst = 0; intIndex = startInt; byteIndex = startByte; for (i = 0; i < numByte; i++) { curInt = (uInt4)ipack[intIndex]; curByte = (curInt << Lshift[byteIndex]) >> 24; *dst = (*dst << 8) + curByte; byteIndex++; if (byteIndex == 4) { byteIndex = 0; intIndex++; } } } /***************************************************************************** * FindTemplateIDs() -- * * Arthur Taylor / MDL * * PURPOSE * This is so we can find out which templates are in this GRIB2 message. * Which allows us to determine if we can call the MDL routines, or if we * should call the NCEP routines instead. * * ARGUMENTS * ipack = The 4 byte int data stream. (Input) * nd5 = length of ipack (Input) * subgNum = Which subgrid to look at. (Input) * gdsTmpl = The gds template number for this subgrid. (Output) * pdsTmpl = The pds template number for this subgrid. (Output) * drsTmpl = The drs template number for this subgrid. (Output) * f_noBitmap = 0 has a bitmap, else doesn't have a bitmap. (Output) * orderDiff = The order of the differencing in template 5.3 (1st, 2nd) (out) * * FILES/DATABASES: None * * RETURNS: int * 0 = Ok. * 2 = invalid section number * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES ***************************************************************************** */ static int FindTemplateIDs(sInt4 *ipack, sInt4 nd5, int subgNum, sInt4 *gdsTmpl, sInt4 *pdsTmpl, sInt4 *drsTmpl, sInt4 *numGrps, uChar *f_noBitmap, sInt4 *orderDiff) { unsigned int gNum = 0; /* Which sub group we are currently working with. */ sInt4 offset; /* Where in the bytes stream we currently are. */ sInt4 sectLen; /* The length of the current section. */ sInt4 sectId; /* The current section number. */ sInt4 li_temp; /* A temporary holder for the bitmap flag. */ /* Jump over section 0. */ offset = 16; while (gNum <= subgNum) { BigByteCpy(§Len, ipack, nd5, (offset / 4), (offset % 4), 4); /* Check if we just read section 8. If so, then it is "7777" = * 926365495 regardless of endian'ness. */ if (sectLen == 926365495L) { #ifdef DEBUG printf("Shouldn't see sect 8. Should stop after correct sect 5\n"); #endif return 2; } BigByteCpy(§Id, ipack, nd5, ((offset + 4) / 4), ((offset + 4) % 4), 1); switch (sectId) { case 1: case 2: break; case 3: BigByteCpy(gdsTmpl, ipack, nd5, ((offset + 12) / 4), ((offset + 12) % 4), 2); break; case 4: BigByteCpy(pdsTmpl, ipack, nd5, ((offset + 7) / 4), ((offset + 7) % 4), 2); break; case 5: BigByteCpy(drsTmpl, ipack, nd5, ((offset + 9) / 4), ((offset + 9) % 4), 2); if ((*drsTmpl == 2) || (*drsTmpl == 3)) { BigByteCpy(numGrps, ipack, nd5, ((offset + 31) / 4), ((offset + 31) % 4), 4); } else { *numGrps = 0; } if (*drsTmpl == 3) { BigByteCpy(&li_temp, ipack, nd5, ((offset + 44) / 4), ((offset + 44) % 4), 1); *orderDiff = li_temp; } else { *orderDiff = 0; } break; case 6: BigByteCpy(&li_temp, ipack, nd5, ((offset + 5) / 4), ((offset + 5) % 4), 1); if (li_temp == 255) { *f_noBitmap = 1; } gNum++; break; case 7: break; default: #ifdef DEBUG printf("Invalid section id %ld.\n", sectId); #endif return 2; } offset += sectLen; } return 0; } /***************************************************************************** * unpk_grib2() -- * * Arthur Taylor / MDL * * PURPOSE * This procedure is the main API for decoding GRIB2 messages. It is * intended to be a branching routine to call either the MDL GRIB2 library, * or the NCEP GRIB2 library, depending on which template it sees in the * message. * * ARGUMENTS * kfildo = Unit number for output diagnostics (C ignores this). (Input) * ain = contains the data if the original data was float data. * (size = nd2x3) (Output) * iain = contains the data if the original data was integer data. * (size = nd2x3) (Output) * nd2x3 = length of ain, iain, ib (Input) * (at least the size of num grid points) * idat = local use data if any that were unpacked from Section 2. (Output) * nidat = length of idat. (Input) * rdat = floating point local use data (Output) * nrdat = length of rdat. (Input) * is0 = Section 0 data (Output) * ns0 = length of is0 (16 is fine) (Input) * is1 = Section 1 data (Output) * ns1 = length of is1 (21 is fine) (Input) * is2 = Section 2 data (Output) * ns2 = length of is2 () (Input) * is3 = Section 3 data (Output) * ns3 = length of is3 (96 or 1600) (Input) * is4 = Section 4 data (Output) * ns4 = length of is4 (60) (Input) * is5 = Section 5 data (Output) * ns5 = length of is5 (49 is fine) (Input) * is6 = Section 6 data (Output) * ns6 = length of is6 (6 is fine) (Input) * is7 = Section 7 data (Output) * ns7 = length of is7 (8 is fine) (Input) * ib = Bitmap if user requested it, and it was packed (Output) * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output) * ipack = The message to unpack (This is assumed to be Big endian) (Input) * nd5 = The size of ipack (Input) * xmissp = The floating point representation for the primary missing value. * (Input/Output) * xmisss = The floating point representation for the secondary missing * value (Input/Output) * new = 1 if this is the first grid to be unpacked, else 0. (Input) * iclean = 1 means the user wants the unpacked data returned without * missing values in it. 0 means embed the missing values. (Input) * l3264b = Integer word length in bits (32 or 64) (Input) * iendpk = A 1 means no more grids in this message, a 0 means more grids. * (Output) * jer(ndjer,2) = error codes along with severity. (Output) * ndjer = 1/2 length of jer. (>= 15) (Input) * kjer = number of error messages stored in jer. (Output) * * FILES/DATABASES: None * * RETURNS: void * * HISTORY * 12/2003 Arthur Taylor (MDL/RSIS): Created. * * NOTES * Inefficiencies: have to memswap ipack multiple times. * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner. * MDL attempts to always return grids in scan mode 0100???? * ToDo: Check length of parameters better. * * According to MDL's unpk_grib2.f, it currently supports (so for others we * call the NCEP routines): * TEMPLATE 3.0 EQUIDISTANT CYLINDRICAL LATITUDE/LONGITUDE * TEMPLATE 3.10 MERCATOR * TEMPLATE 3.20 POLAR STEREOGRAPHIC * TEMPLATE 3.30 LAMBERT * TEMPLATE 3.90 ORTHOGRAPHIC SPACE VIEW * TEMPLATE 3.110 EQUATORIAL AZIMUTHAL EQUIDISTANT * TEMPLATE 3.120 AZIMUTH-RANGE (RADAR) * * TEMPLATE 4.0 ANALYSIS OR FORECAST AT A LEVEL AND POINT * TEMPLATE 4.1 INDIVIDUAL ENSEMBLE * TEMPLATE 4.2 DERIVED FORECAST BASED ON ENSEMBLES * TEMPLATE 4.8 AVERAGE, ACCUMULATION, EXTREMES * TEMPLATE 4.20 RADAR * TEMPLATE 4.30 SATELLITE * * TEMPLATE 5.0 SIMPLE PACKING * TEMPLATE 5.2 COMPLEX PACKING * TEMPLATE 5.3 COMPLEX PACKING AND SPATIAL DIFFERENCING * * Correction to "unpk_grib2.f" : It also supports: * TEMPLATE 4.9 Probability forecast in a time interval * ***************************************************************************** */ void unpk_grib2(sInt4 *kfildo, float *ain, sInt4 *iain, sInt4 *nd2x3, sInt4 *idat, sInt4 *nidat, float *rdat, sInt4 *nrdat, sInt4 *is0, sInt4 *ns0, sInt4 *is1, sInt4 *ns1, sInt4 *is2, sInt4 *ns2, sInt4 *is3, sInt4 *ns3, sInt4 *is4, sInt4 *ns4, sInt4 *is5, sInt4 *ns5, sInt4 *is6, sInt4 *ns6, sInt4 *is7, sInt4 *ns7, sInt4 *ib, sInt4 *ibitmap, sInt4 *ipack, sInt4 *nd5, float *xmissp, float *xmisss, sInt4 *inew, sInt4 *iclean, sInt4 *l3264b, sInt4 *iendpk, sInt4 *jer, sInt4 *ndjer, sInt4 *kjer) { unsigned char *c_ipack; /* The compressed data as char instead of sInt4 so * it is easier to work with. */ char f_useMDL = 0; /* Instructed 3/8/2005 10:30 to not use MDL. */ /* Since NCEP's code works with byte level stuff, we need to un-do the * byte swap of the (int *) data, then cast it to an (unsigned char *) */ #ifndef WORDS_BIGENDIAN /* Can't make this dependent on inew, since they could have a sequence of * get first message... do stuff, get second message, which unfortunately * means they would have to get the first message again, causing 2 swaps, * and breaking their second request for the first message (as well as * their second message). */ memswp(ipack, sizeof(sInt4), *nd5); #endif c_ipack = (unsigned char *)ipack; unpk_g2ncep(kfildo, ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0, ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5, is6, ns6, is7, ns7, ib, ibitmap, c_ipack, nd5, xmissp, xmisss, inew, iclean, l3264b, iendpk, jer, ndjer, kjer); #ifndef WORDS_BIGENDIAN /* Swap back because we could be called again for the subgrid data. */ memswp(ipack, sizeof(sInt4), *nd5); #endif } /***************************************************************************** * C_pkGrib2() -- * * Arthur Taylor / MDL * * PURPOSE * This procedure is the main API for encoding GRIB2 messages. It is * intended to call NCEP's GRIB2 library. * * RETURNS: void * * HISTORY * 1/2004 Arthur Taylor (MDL/RSIS): Created. * * NOTES ***************************************************************************** */ int C_pkGrib2(unsigned char *cgrib, sInt4 *sec0, sInt4 *sec1, unsigned char *csec2, sInt4 lcsec2, sInt4 *igds, sInt4 *igdstmpl, sInt4 *ideflist, sInt4 idefnum, sInt4 ipdsnum, sInt4 *ipdstmpl, float *coordlist, sInt4 numcoord, sInt4 idrsnum, sInt4 *idrstmpl, float *fld, sInt4 ngrdpts, sInt4 ibmap, sInt4 *bmap) { int ierr; /* error value from grib2 library. */ if ((ierr = g2_create(cgrib, sec0, sec1)) == -1) { /* Tried to use for version other than GRIB Ed 2 */ return -1; } if ((ierr = g2_addlocal(cgrib, csec2, lcsec2)) < 0) { /* Some how got a bad section 2. Should be impossible unless an assert * was broken. */ return -2; } if ((ierr = g2_addgrid(cgrib, igds, igdstmpl, ideflist, idefnum)) < 0) { /* Some how got a bad section 3. Only way would be should be with an * unsupported template number unless an assert was broken. */ return -3; } if ((ierr = g2_addfield(cgrib, ipdsnum, ipdstmpl, coordlist, numcoord, idrsnum, idrstmpl, fld, ngrdpts, ibmap, bmap)) < 0) { return -4; } if ((ierr = g2_gribend(cgrib)) < 0) { return -5; } return ierr; }