1 : /*****************************************************************************
2 : * grib2api.c
3 : *
4 : * DESCRIPTION
5 : * This file contains the API to the GRIB2 libraries which is as close as
6 : * possible to the "official" NWS GRIB2 Library's API. The reason for this
7 : * is so we can use NCEP's unofficial GRIB2 library while having minimal
8 : * impact on the already written drivers.
9 : *
10 : * HISTORY
11 : * 12/2003 Arthur Taylor (MDL / RSIS): Created.
12 : *
13 : * NOTES
14 : *****************************************************************************
15 : */
16 : #include <stdlib.h>
17 : #include <string.h>
18 : #include <math.h>
19 : #include "grib2api.h"
20 : #include "grib2.h"
21 : #include "scan.h"
22 : #include "memendian.h"
23 : #include "myassert.h"
24 : #include "gridtemplates.h"
25 : #include "pdstemplates.h"
26 : #include "drstemplates.h"
27 :
28 : //#include "config.h" /*config.h created by configure - ADT mod*/
29 :
30 : /* Declare the external FORTRAN routine
31 : * gcc has two __ if there is one _ in the procedure name. */
32 : #ifdef _FORTRAN
33 : extern void UNPK_G2MDL
34 : (sInt4 * kfildo, sInt4 * jmin, sInt4 * lbit, sInt4 * nov, sInt4 * iwork,
35 : float * ain, sInt4 * iain, sInt4 * nd2x3, sInt4 * idat, sInt4 * nidat,
36 : float * rdat, sInt4 * nrdat, sInt4 * is0, sInt4 * ns0, sInt4 * is1,
37 : sInt4 * ns1, sInt4 * is2, sInt4 * ns2, sInt4 * is3, sInt4 * ns3,
38 : sInt4 * is4, sInt4 * ns4, sInt4 * is5, sInt4 * ns5, sInt4 * is6,
39 : sInt4 * ns6, sInt4 * is7, sInt4 * ns7, sInt4 * ib, sInt4 * ibitmap,
40 : sInt4 * ipack, sInt4 * nd5, float * xmissp, float * xmisss,
41 : sInt4 * inew, sInt4 * iclean, sInt4 * l3264b, sInt4 * iendpk, sInt4 * jer,
42 : sInt4 * ndjer, sInt4 * kjer);
43 :
44 : extern void PK_G2MDL
45 : (sInt4 * kfildo, sInt4 * jmax, sInt4 * jmin, sInt4 * lbit, sInt4 * nov,
46 : sInt4 * misslx, float * a, sInt4 * ia, sInt4 * newbox, sInt4 * newboxp,
47 : float * ain, sInt4 * iain, sInt4 * nx, sInt4 * ny, sInt4 * idat,
48 : sInt4 * nidat, float * rdat, sInt4 * nrdat, sInt4 * is0, sInt4 * ns0,
49 : sInt4 * is1, sInt4 * ns1, sInt4 * is3, sInt4 * ns3, sInt4 * is4,
50 : sInt4 * ns4, sInt4 * is5, sInt4 * ns5, sInt4 * is6, sInt4 * ns6,
51 : sInt4 * is7, sInt4 * ns7, sInt4 * ib, sInt4 * ibitmap, sInt4 * ipack,
52 : sInt4 * nd5, sInt4 * missp, float * xmissp, sInt4 * misss,
53 : float * xmisss, sInt4 * inew, sInt4 * minpk, sInt4 * iclean,
54 : sInt4 * l3264b, sInt4 * jer, sInt4 * ndjer, sInt4 * kjer);
55 : #endif
56 :
57 : /*****************************************************************************
58 : * mdl_LocalUnpack() --
59 : *
60 : * Arthur Taylor / MDL
61 : *
62 : * PURPOSE
63 : * Unpack the local use data assuming that it was packed using the MDL
64 : * encoder. This assumes "local" starts at octet 6 (ie skipping over the
65 : * length and section ID octets)
66 : *
67 : * In Section 2, GRIB2 provides for local use data. The MDL encoder packs
68 : * that data, and signifies that it has done so by setting octet 6 to 1.
69 : * This may have been inadvisable, since the GRIB2 specs did not state that
70 : * octet 6 was special.
71 : *
72 : * ARGUMENTS
73 : * local = The section 2 data prior to being unpacked. (Input)
74 : * locallen = The length of "local" (Input)
75 : * idat = Unpacked MDL local data (if it was an integer) (Output)
76 : * nidat = length of idat. (Input)
77 : * rdat = Unpacked MDL local data (if it was a float) (Output)
78 : * nrdat = length of rdat. (Input)
79 : *
80 : * FILES/DATABASES: None
81 : *
82 : * RETURNS: int
83 : * 0 = Ok.
84 : * 1 = Mixed local data types?
85 : * 2 = nrdat is not large enough.
86 : * 3 = nidat is not large enough.
87 : * 4 = Too many bits to unpack, should be a max of 32 bits.
88 : * 5 = Locallen is too small
89 : *
90 : * HISTORY
91 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
92 : *
93 : * NOTES
94 : *****************************************************************************
95 : */
96 : #define GRIB_UNSIGN_INT2(a,b) ((a<<8)+b)
97 0 : static int mdl_LocalUnpack (unsigned char *local, sInt4 locallen,
98 : sInt4 * idat, sInt4 * nidat, float * rdat,
99 : sInt4 * nrdat)
100 : {
101 0 : int BytesUsed = 0; /* How many bytes have been used. */
102 : unsigned int numGroup; /* Number of groups */
103 : int i; /* Counter over the groups. */
104 : sInt4 numVal; /* Number of values in a given group. */
105 : float refVal; /* The reference value in a given group. */
106 : unsigned int scale; /* The power of 10 scale value in the group. */
107 : sInt4 recScale10; /* 1 / 10**scale.. For faster computations. */
108 : unsigned char numBits; /* # of bits for a single element in the group. */
109 : char f_dataType; /* If the local data is a float or integer. */
110 0 : char f_firstType = 0; /* The type of the first group of local data. */
111 0 : int curIndex = 0; /* Where to store the current data. */
112 : int j; /* Counter over the number of values in a group. */
113 : size_t numUsed; /* Number of bytes used in call to memBitRead. */
114 : uChar bufLoc; /* Where to read for more bits in the data stream. */
115 : uInt4 uli_temp; /* Temporary storage to hold unpacked data. */
116 :
117 0 : if (locallen < BytesUsed + 3) {
118 : #ifdef DEBUG
119 0 : printf ("Locallen is too small.\n");
120 : #endif
121 0 : return 5;
122 : }
123 : /* The calling routine should check octet 6, which is local[0], to be 1,
124 : * so we just assert it is 1. */
125 0 : myAssert (local[0] == 1);
126 0 : numGroup = GRIB_UNSIGN_INT2 (local[1], local[2]);
127 0 : local += 3;
128 0 : BytesUsed += 3;
129 0 : myAssert (*nrdat > 1);
130 0 : myAssert (*nidat > 1);
131 0 : idat[0] = 0;
132 0 : rdat[0] = 0;
133 :
134 0 : for (i = 0; i < numGroup; i++) {
135 0 : if (locallen < BytesUsed + 12) {
136 : #ifdef DEBUG
137 0 : printf ("Locallen is too small.\n");
138 : #endif
139 0 : return 5;
140 : }
141 0 : MEMCPY_BIG (&numVal, local, sizeof (sInt4));
142 0 : MEMCPY_BIG (&refVal, local + 4, sizeof (float));
143 0 : scale = GRIB_UNSIGN_INT2 (local[8], local[9]);
144 0 : recScale10 = 1 / pow (10.0, scale);
145 0 : numBits = local[10];
146 0 : if (numBits >= 32) {
147 : #ifdef DEBUG
148 0 : printf ("Too many bits too unpack.\n");
149 : #endif
150 0 : return 4;
151 : }
152 0 : f_dataType = local[11];
153 0 : local += 12;
154 0 : BytesUsed += 12;
155 0 : if (locallen < BytesUsed + ((numBits * numVal) + 7) / 8) {
156 : #ifdef DEBUG
157 0 : printf ("Locallen is too small.\n");
158 : #endif
159 0 : return 5;
160 : }
161 0 : if (i == 0) {
162 0 : f_firstType = f_dataType;
163 0 : } else if (f_firstType != f_dataType) {
164 : #ifdef DEBUG
165 0 : printf ("Local use data has mixed float/integer type?\n");
166 : #endif
167 0 : return 1;
168 : }
169 0 : bufLoc = 8;
170 0 : if (f_dataType == 0) {
171 : /* Floating point data. */
172 0 : if (*nrdat < (curIndex + numVal + 3)) {
173 : #ifdef DEBUG
174 0 : printf ("nrdat is not large enough.\n");
175 : #endif
176 0 : return 2;
177 : }
178 0 : rdat[curIndex] = numVal;
179 0 : curIndex++;
180 0 : rdat[curIndex] = scale;
181 0 : curIndex++;
182 0 : for (j = 0; j < numVal; j++) {
183 0 : memBitRead (&uli_temp, sizeof (sInt4), local, numBits,
184 : &bufLoc, &numUsed);
185 0 : local += numUsed;
186 0 : BytesUsed += numUsed;
187 0 : rdat[curIndex] = (refVal + uli_temp) * recScale10;
188 0 : curIndex++;
189 : }
190 0 : rdat[curIndex] = 0;
191 : } else {
192 : /* Integer point data. */
193 0 : if (*nidat < (curIndex + numVal + 3)) {
194 : #ifdef DEBUG
195 0 : printf ("nidat is not large enough.\n");
196 : #endif
197 0 : return 3;
198 : }
199 0 : idat[curIndex] = numVal;
200 0 : curIndex++;
201 0 : idat[curIndex] = scale;
202 0 : curIndex++;
203 0 : for (j = 0; j < numVal; j++) {
204 0 : memBitRead (&uli_temp, sizeof (sInt4), local, numBits,
205 : &bufLoc, &numUsed);
206 0 : local += numUsed;
207 0 : BytesUsed += numUsed;
208 0 : idat[curIndex] = (refVal + uli_temp) * recScale10;
209 0 : curIndex++;
210 : }
211 0 : idat[curIndex] = 0;
212 : }
213 : }
214 0 : return 0;
215 : }
216 :
217 : /*****************************************************************************
218 : * fillOutSectLen() --
219 : *
220 : * Arthur Taylor / MDL
221 : *
222 : * PURPOSE
223 : * The GRIB2 API returns the lengths of each GRIB2 section along with the
224 : * meta data. NCEP's routines did not, so this routine fills in that part.
225 : * This routine has to consider which grid is being unpacked.
226 : *
227 : * c_ipack is passed in after section 1 (since section 1 is not repeated)
228 : *
229 : * ARGUMENTS
230 : * c_ipack = Complete GRIB2 message to look for the section lengths in. (In)
231 : * lenCpack = Length of c_ipack (Input)
232 : * subgNum = Which sub rid to find the section lengths of. (Input)
233 : * is2 = Section 2 data (Output)
234 : * is3 = Section 3 data (Output)
235 : * is4 = Section 4 data (Output)
236 : * is5 = Section 5 data (Output)
237 : * is6 = Section 6 data (Output)
238 : * is7 = Section 7 data (Output)
239 : *
240 : * FILES/DATABASES: None
241 : *
242 : * RETURNS: int
243 : * 0 = Ok.
244 : * 1 = c_ipack is not large enough
245 : * 2 = invalid section number
246 : *
247 : * HISTORY
248 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
249 : *
250 : * NOTES
251 : *****************************************************************************
252 : */
253 2 : static int fillOutSectLen (unsigned char *c_ipack, int lenCpack,
254 : int subgNum, sInt4 * is2, sInt4 * is3,
255 : sInt4 * is4, sInt4 * is5, sInt4 * is6,
256 : sInt4 * is7)
257 : {
258 2 : sInt4 offset = 0; /* How far in c_ipack we have read. */
259 : sInt4 sectLen; /* The length of the current section. */
260 : int sectId; /* The current section number. */
261 2 : unsigned int gNum = 0; /* Which sub group we are currently working with. */
262 :
263 2 : if (lenCpack < 5) {
264 : #ifdef DEBUG
265 0 : printf ("Cpack is not large enough.\n");
266 : #endif
267 0 : return 1;
268 : }
269 : /* assert that we start with data in either section 2 or 3. */
270 2 : myAssert ((c_ipack[4] == 2) || (c_ipack[4] == 3));
271 14 : while (gNum <= subgNum) {
272 10 : if (lenCpack < offset + 5) {
273 : #ifdef DEBUG
274 0 : printf ("Cpack is not large enough.\n");
275 : #endif
276 0 : return 1;
277 : }
278 10 : MEMCPY_BIG (§Len, c_ipack + offset, sizeof (sInt4));
279 : /* Check if we just read section 8. If so, then it is "7777" =
280 : * 926365495 regardless of endian'ness. */
281 10 : if (sectLen == 926365495L) {
282 : #ifdef DEBUG
283 0 : printf ("Shouldn't see sect 8. Should stop after correct sect 7\n");
284 : #endif
285 0 : return 2;
286 : }
287 10 : sectId = c_ipack[offset + 4];
288 10 : switch (sectId) {
289 : case 2:
290 0 : is2[0] = sectLen;
291 0 : break;
292 : case 3:
293 2 : is3[0] = sectLen;
294 2 : break;
295 : case 4:
296 2 : is4[0] = sectLen;
297 2 : break;
298 : case 5:
299 2 : is5[0] = sectLen;
300 2 : break;
301 : case 6:
302 2 : is6[0] = sectLen;
303 2 : break;
304 : case 7:
305 2 : is7[0] = sectLen;
306 2 : gNum++;
307 2 : break;
308 : default:
309 : #ifdef DEBUG
310 0 : printf ("Invalid section id %d.\n", sectId);
311 : #endif
312 0 : return 2;
313 : }
314 10 : offset += sectLen;
315 : }
316 2 : return 0;
317 : }
318 :
319 : /*****************************************************************************
320 : * TransferInt() --
321 : *
322 : * Arthur Taylor / MDL
323 : *
324 : * PURPOSE
325 : * To transfer the data from "fld" and "bmap" to "iain" and "ib". The API
326 : * attempts to rearrange the order, so that anything returned from it has
327 : * scan mode 0100????
328 : *
329 : * ARGUMENTS
330 : * fld = The expanded grid from NCEPs routines (Input)
331 : * ngrdpts = Length of fld (Input)
332 : * ibitmap = flag if we have a bitmap or not. (Input)
333 : * bmap = bitmap from NCEPs routines. (Input)
334 : * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input)
335 : * scan = The scan orientation of fld/bmap/iain/ib (Input/Output)
336 : * nx, ny = The dimmensions of the grid. (Input)
337 : * iclean = 1 means the user wants the unpacked data returned without
338 : * missing values in it. 0 means embed the missing values. (In)
339 : * xmissp = The primary missing value to use if iclean = 0. (Input).
340 : * iain = The grid to return. (Output)
341 : * nd2x3 = The length of iain. (Input)
342 : * ib = Bitmap if it was packed (Output)
343 : *
344 : * FILES/DATABASES: None
345 : *
346 : * RETURNS: int
347 : * 0 = Ok.
348 : * 1 = nd2x3 is too small.
349 : * 2 = nx*ny != ngrdpts
350 : *
351 : * HISTORY
352 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
353 : *
354 : * NOTES
355 : * May want to disable the scan adjustment in the future.
356 : *****************************************************************************
357 : */
358 0 : static int TransferInt (float * fld, sInt4 ngrdpts, sInt4 ibitmap,
359 : sInt4 * bmap, char f_ignoreScan, sInt4 * scan,
360 : sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp,
361 : sInt4 * iain, sInt4 nd2x3, sInt4 * ib)
362 : {
363 : int i; /* loop counter over all grid points. */
364 : sInt4 x, y; /* Where we are in a grid of scan value 0100???? */
365 : int curIndex; /* Where in iain to store the current data. */
366 :
367 0 : if (nd2x3 < ngrdpts) {
368 : #ifdef DEBUG
369 0 : printf ("nd2x3(%d) is < ngrdpts(%d)\n", nd2x3, ngrdpts);
370 : #endif
371 0 : return 1;
372 : }
373 0 : if (f_ignoreScan || ((*scan & 0xf0) == 64)) {
374 0 : if (ibitmap) {
375 0 : for (i = 0; i < ngrdpts; i++) {
376 0 : ib[i] = bmap[i];
377 : /* Check if we are supposed to insert xmissp into the field */
378 0 : if ((iclean != 0) && (ib[i] == 0)) {
379 0 : iain[i] = xmissp;
380 : } else {
381 0 : iain[i] = fld[i];
382 : }
383 : }
384 : } else {
385 0 : for (i = 0; i < ngrdpts; i++) {
386 0 : iain[i] = fld[i];
387 : }
388 : }
389 : } else {
390 0 : if (nx * ny != ngrdpts) {
391 : #ifdef DEBUG
392 0 : printf ("nx * ny (%d) != ngrdpts(%d)\n", nx * ny, ngrdpts);
393 : #endif
394 0 : return 2;
395 : }
396 0 : if (ibitmap) {
397 0 : for (i = 0; i < ngrdpts; i++) {
398 0 : ScanIndex2XY (i, &x, &y, *scan, nx, ny);
399 : /* ScanIndex returns value as if scan was 0100(0000) */
400 0 : curIndex = (x - 1) + (y - 1) * nx;
401 0 : myAssert (curIndex < nd2x3);
402 0 : ib[curIndex] = bmap[i];
403 : /* Check if we are supposed to insert xmissp into the field */
404 0 : if ((iclean != 0) && (ib[curIndex] == 0)) {
405 0 : iain[i] = xmissp;
406 : } else {
407 0 : iain[curIndex] = fld[i];
408 : }
409 : }
410 : } else {
411 0 : for (i = 0; i < ngrdpts; i++) {
412 0 : ScanIndex2XY (i, &x, &y, *scan, nx, ny);
413 : /* ScanIndex returns value as if scan was 0100(0000) */
414 0 : curIndex = (x - 1) + (y - 1) * nx;
415 0 : myAssert (curIndex < nd2x3);
416 0 : iain[curIndex] = fld[i];
417 : }
418 : }
419 0 : *scan = 64 + (*scan & 0x0f);
420 : }
421 0 : return 0;
422 : }
423 :
424 : /*****************************************************************************
425 : * TransferFloat() --
426 : *
427 : * Arthur Taylor / MDL
428 : *
429 : * PURPOSE
430 : * To transfer the data from "fld" and "bmap" to "ain" and "ib". The API
431 : * attempts to rearrange the order, so that anything returned from it has
432 : * scan mode 0100????
433 : *
434 : * ARGUMENTS
435 : * fld = The expanded grid from NCEPs routines (Input)
436 : * ngrdpts = Length of fld (Input)
437 : * ibitmap = flag if we have a bitmap or not. (Input)
438 : * bmap = bitmap from NCEPs routines. (Input)
439 : * scan = The scan orientation of fld/bmap/iain/ib (Input/Output)
440 : * f_ignoreScan = Flag to ignore the attempt at changing the scan (Input)
441 : * nx, ny = The dimmensions of the grid. (Input)
442 : * iclean = 1 means the user wants the unpacked data returned without
443 : * missing values in it. 0 means embed the missing values. (In)
444 : * xmissp = The primary missing value to use if iclean = 0. (Input).
445 : * ain = The grid to return. (Output)
446 : * nd2x3 = The length of iain. (Input)
447 : * ib = Bitmap if it was packed (Output)
448 : *
449 : * FILES/DATABASES: None
450 : *
451 : * RETURNS: int
452 : * 0 = Ok.
453 : * 1 = nd2x3 is too small.
454 : * 2 = nx*ny != ngrdpts
455 : *
456 : * HISTORY
457 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
458 : *
459 : * NOTES
460 : * May want to disable the scan adjustment in the future.
461 : *****************************************************************************
462 : */
463 2 : static int TransferFloat (float * fld, sInt4 ngrdpts, sInt4 ibitmap,
464 : sInt4 * bmap, char f_ignoreScan, sInt4 * scan,
465 : sInt4 nx, sInt4 ny, sInt4 iclean, float xmissp,
466 : float * ain, sInt4 nd2x3, sInt4 * ib)
467 : {
468 : int i; /* loop counter over all grid points. */
469 : sInt4 x, y; /* Where we are in a grid of scan value 0100???? */
470 : int curIndex; /* Where in ain to store the current data. */
471 :
472 2 : if (nd2x3 < ngrdpts) {
473 : #ifdef DEBUG
474 0 : printf ("nd2x3(%d) is < ngrdpts(%d)\n", nd2x3, ngrdpts);
475 : #endif
476 0 : return 1;
477 : }
478 2 : if (f_ignoreScan || ((*scan & 0xf0) == 64)) {
479 0 : if (ibitmap) {
480 0 : for (i = 0; i < ngrdpts; i++) {
481 0 : ib[i] = bmap[i];
482 : /* Check if we are supposed to insert xmissp into the field */
483 0 : if ((iclean != 0) && (ib[i] == 0)) {
484 0 : ain[i] = xmissp;
485 : } else {
486 0 : ain[i] = fld[i];
487 : }
488 : }
489 : } else {
490 0 : for (i = 0; i < ngrdpts; i++) {
491 0 : ain[i] = fld[i];
492 : }
493 : }
494 : } else {
495 2 : if (nx * ny != ngrdpts) {
496 : #ifdef DEBUG
497 0 : printf ("nx * ny (%d) != ngrdpts(%d)\n", nx * ny, ngrdpts);
498 : #endif
499 0 : return 2;
500 : }
501 2 : if (ibitmap) {
502 0 : for (i = 0; i < ngrdpts; i++) {
503 0 : ScanIndex2XY (i, &x, &y, *scan, nx, ny);
504 : /* ScanIndex returns value as if scan was 0100(0000) */
505 0 : curIndex = (x - 1) + (y - 1) * nx;
506 0 : myAssert (curIndex < nd2x3);
507 0 : ib[curIndex] = bmap[i];
508 : /* Check if we are supposed to insert xmissp into the field */
509 0 : if ((iclean != 0) && (ib[curIndex] == 0)) {
510 0 : ain[i] = xmissp;
511 : } else {
512 0 : ain[curIndex] = fld[i];
513 : }
514 : }
515 : } else {
516 45668 : for (i = 0; i < ngrdpts; i++) {
517 45666 : ScanIndex2XY (i, &x, &y, *scan, nx, ny);
518 : /* ScanIndex returns value as if scan was 0100(0000) */
519 45666 : curIndex = (x - 1) + (y - 1) * nx;
520 45666 : myAssert (curIndex < nd2x3);
521 45666 : ain[curIndex] = fld[i];
522 : }
523 : }
524 2 : *scan = 64 + (*scan & 0x0f);
525 : }
526 : /*
527 : if (1==1) {
528 : int i;
529 : for (i=0; i < ngrdpts; i++) {
530 : if (ain[i] != 9999) {
531 : fprintf (stderr, "grib2api.c: ain[%d] %f, fld[%d] %f\n", i, ain[i],
532 : i, fld[i]);
533 : }
534 : }
535 : }
536 : */
537 2 : return 0;
538 : }
539 :
540 : #ifdef PKNCEP
541 : int pk_g2ncep (sInt4 * kfildo, float * ain, sInt4 * iain, sInt4 * nx,
542 : sInt4 * ny, sInt4 * idat, sInt4 * nidat, float * rdat,
543 : sInt4 * nrdat, sInt4 * is0, sInt4 * ns0, sInt4 * is1,
544 : sInt4 * ns1, sInt4 * is3, sInt4 * ns3, sInt4 * is4,
545 : sInt4 * ns4, sInt4 * is5, sInt4 * ns5, sInt4 * is6,
546 : sInt4 * ns6, sInt4 * is7, sInt4 * ns7, sInt4 * ib,
547 : sInt4 * ibitmap, unsigned char *cgrib, sInt4 * nd5,
548 : sInt4 * missp, float * xmissp, sInt4 * misss,
549 : float * xmisss, sInt4 * inew, sInt4 * minpk, sInt4 * iclean,
550 : sInt4 * l3264b, sInt4 * jer, sInt4 * ndjer, sInt4 * kjer)
551 : {
552 : g2int listsec0[2];
553 : g2int listsec1[13];
554 : g2int igds[5];
555 : g2int igdstmpl[];
556 : int ierr; /* Holds the error code from a called routine. */
557 : int i;
558 :
559 : listsec0[0] = is0[6];
560 : listsec0[1] = is0[7];
561 :
562 : listsec1[0] = is1[5];
563 : listsec1[1] = is1[7];
564 : listsec1[2] = is1[9];
565 : listsec1[3] = is1[10];
566 : listsec1[4] = is1[11];
567 : listsec1[5] = is1[12];
568 : listsec1[6] = is1[14];
569 : listsec1[7] = is1[15];
570 : listsec1[8] = is1[16];
571 : listsec1[9] = is1[17];
572 : listsec1[10] = is1[18];
573 : listsec1[11] = is1[19];
574 : listsec1[12] = is1[20];
575 :
576 : ierr = g2_create (cgrib, listsec0, listsec1);
577 : printf ("Length = %d\n", ierr);
578 :
579 : if ((idat[0] != 0) || (rdat[0] != 0)) {
580 : printf ("Don't handle this yet.\n");
581 : /*
582 : ierr = g2_addlocal (cgrib, unsigned char *csec2, g2int lcsec2);
583 : */
584 : }
585 :
586 : igds[0] = is3[5];
587 : igds[1] = is3[6];
588 : igds[2] = is3[10];
589 : igds[3] = is3[11];
590 : igds[4] = is3[12];
591 :
592 : IS3(15) - IS3(nn) = Grid Definition Template, stored in bytes 15-nn (*)
593 :
594 : ierr = g2_addgrid (cgrib, igds, g2int *igdstmpl, g2int *ideflist,
595 : g2int idefnum)
596 :
597 :
598 :
599 :
600 : return 0;
601 : /*
602 :
603 : To start a new GRIB2 message, call function g2_create. G2_create
604 : encodes Sections 0 and 1 at the beginning of the message. This routine
605 : must be used to create each message.
606 :
607 : Routine g2_addlocal can be used to add a Local Use Section ( Section 2 ).
608 : Note that this section is optional and need not appear in a GRIB2 message.
609 :
610 : Function g2_addgrid is used to encode a grid definition into Section 3.
611 : This grid definition defines the geometry of the the data values in the
612 : fields that follow it. g2_addgrid can be called again to change the grid
613 : definition describing subsequent data fields.
614 :
615 : Each data field is added to the GRIB2 message using routine g2_addfield,
616 : which adds Sections 4, 5, 6, and 7 to the message.
617 :
618 : After all desired data fields have been added to the GRIB2 message, a
619 : call to function g2_gribend is needed to add the final section 8 to the
620 : message and to update the length of the message. A call to g2_gribend
621 : is required for each GRIB2 message.
622 : */
623 :
624 : }
625 : #endif
626 :
627 : /*****************************************************************************
628 : * unpk_g2ncep() --
629 : *
630 : * Arthur Taylor / MDL
631 : *
632 : * PURPOSE
633 : * This procedure is a wrapper around the NCEP GRIB2 routines to interface
634 : * their routines with the "official NWS" GRIB2 API. The reason for this is
635 : * so drivers that have been written to use the "official NWS" GRIB2 API, can
636 : * use the NCEP library with minimal disruption.
637 : *
638 : * ARGUMENTS
639 : * kfildo = Unit number for output diagnostics (C ignores this). (Input)
640 : * ain = contains the data if the original data was float data.
641 : * (size = nd2x3) (Output)
642 : * iain = contains the data if the original data was integer data.
643 : * (size = nd2x3) (Output)
644 : * nd2x3 = length of ain, iain, ib (Input)
645 : * (at least the size of num grid points)
646 : * idat = local use data if any that were unpacked from Section 2. (Output)
647 : * nidat = length of idat. (Input)
648 : * rdat = floating point local use data (Output)
649 : * nrdat = length of rdat. (Input)
650 : * is0 = Section 0 data (Output)
651 : * ns0 = length of is0 (16 is fine) (Input)
652 : * is1 = Section 1 data (Output)
653 : * ns1 = length of is1 (21 is fine) (Input)
654 : * is2 = Section 2 data (Output)
655 : * ns2 = length of is2 () (Input)
656 : * is3 = Section 3 data (Output)
657 : * ns3 = length of is3 (96 or 1600) (Input)
658 : * is4 = Section 4 data (Output)
659 : * ns4 = length of is4 (60) (Input)
660 : * is5 = Section 5 data (Output)
661 : * ns5 = length of is5 (49 is fine) (Input)
662 : * is6 = Section 6 data (Output)
663 : * ns6 = length of is6 (6 is fine) (Input)
664 : * is7 = Section 7 data (Output)
665 : * ns7 = length of is7 (8 is fine) (Input)
666 : * ib = Bitmap if user requested it, and it was packed (Output)
667 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output)
668 : * c_ipack = The message to unpack (Input)
669 : * nd5 = 1/4 the size of c_ipack (Input)
670 : * xmissp = The floating point representation for the primary missing value.
671 : * (Input/Output)
672 : * xmisss = The floating point representation for the secondary missing
673 : * value (Input/Output)
674 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
675 : * iclean = 1 means the user wants the unpacked data returned without
676 : * missing values in it. 0 means embed the missing values. (Input)
677 : * l3264b = Integer word length in bits (32 or 64) (Input)
678 : * iendpk = A 1 means no more grids in this message, a 0 means more grids.
679 : * (Output)
680 : * jer(ndjer,2) = error codes along with severity. (Output)
681 : * ndjer = 1/2 length of jer. (>= 15) (Input)
682 : * kjer = number of error messages stored in jer.
683 : *
684 : * FILES/DATABASES: None
685 : *
686 : * RETURNS: void
687 : * "jer" is used to store error messages, and kjer is used to denote how
688 : * many there were. Jer is set up as a 2 column array, with the first
689 : * column in the first half of the array, and denoting the GRIB2 section
690 : * an error occurred in. The second column denoting the severity.
691 : *
692 : * The possibilities from unpk_g2ncep() are as follows:
693 : * ker=1 jer[0,0]=0 jer[0,1]=2: Message is not formed correctly
694 : * Request for an invalid subgrid
695 : * Problems unpacking the data.
696 : * problems expanding the data.
697 : * Calling dimmensions were too small.
698 : * ker=2 jer[1,0]=100 jer[1,1]=2: Error unpacking section 1.
699 : * ker=3 jer[2,0]=200 jer[2,1]=2: Error unpacking section 2.
700 : * ker=4 jer[3,0]=300 jer[3,1]=2: Error unpacking section 3.
701 : * ker=5 jer[4,0]=400 jer[4,1]=2: Error unpacking section 4.
702 : * ker=6 jer[5,0]=500 jer[5,1]=2: Error unpacking section 5.
703 : * Data Template not implemented.
704 : * Durring Transfer, nx * ny != ngrdpts.
705 : * ker=7 jer[6,0]=600 jer[6,1]=2: Error unpacking section 6.
706 : * ker=8 jer[7,0]=700 jer[7,1]=2: Error unpacking section 7.
707 : * ker=9 jer[8,0]=2001 jer[8,1]=2: nd2x3 is not large enough.
708 : * ker=9 jer[8,0]=2003 jer[8,1]=2: undefined sect 3 template
709 : * (see gridtemplates.h).
710 : * ker=9 jer[8,0]=2004 jer[8,1]=2: undefined sect 4 template
711 : * (see pdstemplates.h).
712 : * ker=9 jer[8,0]=2005 jer[8,1]=2: undefined sect 5 template
713 : * (see drstemplates.h).
714 : * ker=9 jer[8,0]=9999 jer[8,1]=2: NCEP returns an unrecognized error.
715 : *
716 : * HISTORY
717 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
718 : *
719 : * NOTES
720 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
721 : * MDL attempts to always return grids in scan mode 0100????
722 : *
723 : * ToDo: Check length of parameters better.
724 : * ToDo: Probably shouldn't abort if they have problems expanding the data.
725 : * ToDo: Better handling of:
726 : * gfld->list_opt = (Used if gfld->numoct_opt .ne. 0) This array contains
727 : * the number of grid points contained in each row (or column). (part of
728 : * Section 3) This element is a pointer to an array that holds the data.
729 : * This pointer is nullified if gfld->numoct_opt=0.
730 : * gfld->num_opt = (Used if gfld->numoct_opt .ne. 0) The number of entries in
731 : * array ideflist. i.e. number of rows (or columns) for which optional
732 : * grid points are defined. This value is set to zero, if
733 : * gfld->numoct_opt=0.
734 : * gfld->coord_list = Real array containing floating point values intended to
735 : * document the vertical discretisation associated to model data on
736 : * hybrid coordinate vertical levels. (part of Section 4) This element
737 : * is a pointer to an array that holds the data.
738 : * gfld->num_coord = number of values in array gfld->coord_list[].
739 : *****************************************************************************
740 : */
741 2 : void unpk_g2ncep (sInt4 * kfildo, float * ain, sInt4 * iain, sInt4 * nd2x3,
742 : sInt4 * idat, sInt4 * nidat, float * rdat, sInt4 * nrdat,
743 : sInt4 * is0, sInt4 * ns0, sInt4 * is1, sInt4 * ns1,
744 : sInt4 * is2, sInt4 * ns2, sInt4 * is3, sInt4 * ns3,
745 : sInt4 * is4, sInt4 * ns4, sInt4 * is5, sInt4 * ns5,
746 : sInt4 * is6, sInt4 * ns6, sInt4 * is7, sInt4 * ns7,
747 : sInt4 * ib, sInt4 * ibitmap, unsigned char *c_ipack,
748 : sInt4 * nd5, float * xmissp, float * xmisss,
749 : sInt4 * inew, sInt4 * iclean, sInt4 * l3264b,
750 : sInt4 * iendpk, sInt4 * jer, sInt4 * ndjer, sInt4 * kjer)
751 : {
752 : int i; /* A counter used for a number of purposes. */
753 : static unsigned int subgNum = 0; /* The sub grid we read most recently.
754 : * This is primarily to help with the
755 : * inew option. */
756 : int ierr; /* Holds the error code from a called routine. */
757 : sInt4 listsec0[3];
758 : sInt4 listsec1[13];
759 : static sInt4 numfields = 1; /* Number of sub Grids in this message */
760 : sInt4 numlocal; /* Number of local sections in this message. */
761 : int unpack; /* Tell g2_getfld to unpack the message. */
762 : int expand; /* Tell g2_getflt to attempt to expand the bitmap. */
763 : gribfield *gfld; /* Holds the data after g2_getfld unpacks it. */
764 : sInt4 gridIndex; /* index in templatesgrid[] for this sect 3 templat */
765 : sInt4 pdsIndex; /* index in templatespds[] for this sect 4 template */
766 : sInt4 drsIndex; /* index in templatesdrs[] for this sect 5 template */
767 : int curIndex; /* Where in is3, is4, or is5 to store meta data */
768 : int scanIndex; /* Where in is3 to find the scan mode. */
769 : int nxIndex; /* Where in is3 to find the number of x values. */
770 : int nyIndex; /* Where in is3 to find the number of y values. */
771 : float f_temp; /* Assist with handling weird MDL behavior in is5[] */
772 : char f_ignoreScan; /* Flag to ignore the attempt at changing the scan */
773 : sInt4 dummyScan; /* Dummy place holder for call to Transfer routines
774 : * if ignoring scan. */
775 :
776 2 : myAssert (*ndjer >= 8);
777 : /* Init the error handling array. */
778 2 : memset ((void *) jer, 0, 2 * *ndjer * sizeof (sInt4));
779 18 : for (i = 0; i < 8; i++) {
780 16 : jer[i] = i * 100;
781 : }
782 2 : *kjer = 8;
783 :
784 : /* The first time in, figure out how many grids there are, and store it
785 : * in numfields for subsequent calls with inew != 1. */
786 2 : if (*inew == 1) {
787 2 : subgNum = 0;
788 2 : ierr = g2_info (c_ipack, listsec0, listsec1, &numfields, &numlocal);
789 2 : if (ierr != 0) {
790 0 : switch (ierr) {
791 : case 1: /* Beginning characters "GRIB" not found. */
792 : case 2: /* GRIB message is not Edition 2. */
793 : case 3: /* Could not find Section 1, where expected. */
794 : case 4: /* End string "7777" found, but not where expected. */
795 : case 5: /* End string "7777" not found at end of message. */
796 : case 6: /* Invalid section number found. */
797 0 : jer[0 + *ndjer] = 2;
798 0 : *kjer = 1;
799 0 : break;
800 : default:
801 0 : jer[8 + *ndjer] = 2;
802 0 : jer[8] = 9999; /* Unknown error message. */
803 0 : *kjer = 9;
804 : }
805 0 : return;
806 : }
807 : } else {
808 0 : if (subgNum + 1 >= numfields) {
809 : /* Field request error. */
810 0 : jer[0 + *ndjer] = 2;
811 0 : *kjer = 1;
812 0 : return;
813 : }
814 0 : subgNum++;
815 : }
816 :
817 : /* Expand the desired subgrid. */
818 2 : unpack = 1;
819 2 : expand = 1;
820 2 : ierr = g2_getfld (c_ipack, subgNum + 1, unpack, expand, &gfld);
821 2 : if (ierr != 0) {
822 0 : switch (ierr) {
823 : case 1: /* Beginning characters "GRIB" not found. */
824 : case 2: /* GRIB message is not Edition 2. */
825 : case 3: /* The data field request number was not positive. */
826 : case 4: /* End string "7777" found, but not where expected. */
827 : case 6: /* message did not contain requested # of fields */
828 : case 7: /* End string "7777" not found at end of message. */
829 : case 8: /* Unrecognized Section encountered. */
830 0 : jer[0 + *ndjer] = 2;
831 0 : *kjer = 1;
832 0 : break;
833 : case 15: /* Error unpacking Section 1. */
834 0 : jer[1 + *ndjer] = 2;
835 0 : *kjer = 2;
836 0 : break;
837 : case 16: /* Error unpacking Section 2. */
838 0 : jer[2 + *ndjer] = 2;
839 0 : *kjer = 3;
840 0 : break;
841 : case 10: /* Error unpacking Section 3. */
842 0 : jer[3 + *ndjer] = 2;
843 0 : *kjer = 4;
844 0 : break;
845 : case 11: /* Error unpacking Section 4. */
846 0 : jer[4 + *ndjer] = 2;
847 0 : *kjer = 5;
848 0 : break;
849 : case 9: /* Data Template 5.NN not implemented. */
850 : case 12: /* Error unpacking Section 5. */
851 0 : jer[5 + *ndjer] = 2;
852 0 : *kjer = 6;
853 0 : break;
854 : case 13: /* Error unpacking Section 6. */
855 0 : jer[6 + *ndjer] = 2;
856 0 : *kjer = 7;
857 0 : break;
858 : case 14: /* Error unpacking Section 7. */
859 0 : jer[7 + *ndjer] = 2;
860 0 : *kjer = 8;
861 0 : break;
862 : default:
863 0 : jer[8 + *ndjer] = 2;
864 0 : jer[8] = 9999; /* Unknown error message. */
865 0 : *kjer = 9;
866 : break;
867 : }
868 0 : g2_free (gfld);
869 0 : return;
870 : }
871 : /* Check if data wasn't unpacked. */
872 2 : if (!gfld->unpacked) {
873 0 : jer[0 + *ndjer] = 2;
874 0 : *kjer = 1;
875 0 : g2_free (gfld);
876 0 : return;
877 : }
878 :
879 : /* Start going through the gfld structure and converting it to the needed
880 : * data output formats. */
881 2 : myAssert (*ns0 >= 16);
882 2 : MEMCPY_BIG (&(is0[0]), c_ipack, sizeof (sInt4));
883 2 : is0[6] = gfld->discipline;
884 2 : is0[7] = gfld->version;
885 2 : MEMCPY_BIG (&(is0[8]), c_ipack + 8, sizeof (sInt4));
886 : /* The following assert fails only if the GRIB message is more that 4
887 : * giga-bytes large, which I think would break the fortran library. */
888 2 : myAssert (is0[8] == 0);
889 2 : MEMCPY_BIG (&(is0[8]), c_ipack + 12, sizeof (sInt4));
890 :
891 2 : myAssert (*ns1 >= 21);
892 2 : myAssert (gfld->idsectlen >= 13);
893 2 : MEMCPY_BIG (&(is1[0]), c_ipack + 16, sizeof (sInt4));
894 2 : is1[4] = c_ipack[20];
895 2 : is1[5] = gfld->idsect[0];
896 2 : is1[7] = gfld->idsect[1];
897 2 : is1[9] = gfld->idsect[2];
898 2 : is1[10] = gfld->idsect[3];
899 2 : is1[11] = gfld->idsect[4];
900 2 : is1[12] = gfld->idsect[5]; /* Year */
901 2 : is1[14] = gfld->idsect[6]; /* Month */
902 2 : is1[15] = gfld->idsect[7]; /* Day */
903 2 : is1[16] = gfld->idsect[8]; /* Hour */
904 2 : is1[17] = gfld->idsect[9]; /* Min */
905 2 : is1[18] = gfld->idsect[10]; /* Sec */
906 2 : is1[19] = gfld->idsect[11];
907 2 : is1[20] = gfld->idsect[12];
908 :
909 : /* Fill out section lengths (separate procedure because of possibility of
910 : * having multiple grids. Should combine fillOutSectLen g2_info, and
911 : * g2_getfld into one procedure to optimize it. */
912 2 : fillOutSectLen (c_ipack + 16 + is1[0], 4 * *nd5 - 15 - is1[0], subgNum,
913 : is2, is3, is4, is5, is6, is7);
914 :
915 : /* Check if there is section 2 data. */
916 2 : if (gfld->locallen > 0) {
917 : /* The + 1 is so we don't overwrite the section length */
918 0 : memset ((void *) (is2 + 1), 0, (*ns2 - 1) * sizeof (sInt4));
919 0 : is2[4] = 2;
920 0 : is2[5] = gfld->local[0];
921 : /* check if MDL Local use simple packed data */
922 0 : if (is2[5] == 1) {
923 0 : mdl_LocalUnpack (gfld->local, gfld->locallen, idat, nidat,
924 : rdat, nrdat);
925 : } else {
926 : /* local use section was not MDL packed, return it in is2. */
927 0 : for (i = 0; i < gfld->locallen; i++) {
928 0 : is2[i + 5] = gfld->local[i];
929 : }
930 : }
931 : } else {
932 : /* API specified that is2[0] = 0; idat[0] = 0, and rdat[0] = 0 */
933 2 : is2[0] = 0;
934 2 : idat[0] = 0;
935 2 : rdat[0] = 0;
936 : }
937 :
938 2 : is3[4] = 3;
939 2 : is3[5] = gfld->griddef;
940 2 : is3[6] = gfld->ngrdpts;
941 2 : if (*nd2x3 < gfld->ngrdpts) {
942 0 : jer[8 + *ndjer] = 2;
943 0 : jer[8] = 2001; /* nd2x3 is not large enough */
944 0 : *kjer = 9;
945 0 : g2_free (gfld);
946 0 : return;
947 : }
948 2 : is3[10] = gfld->numoct_opt;
949 2 : is3[11] = gfld->interp_opt;
950 2 : is3[12] = gfld->igdtnum;
951 2 : gridIndex = getgridindex (gfld->igdtnum);
952 2 : if (gridIndex == -1) {
953 0 : jer[8 + *ndjer] = 2;
954 0 : jer[8] = 2003; /* undefined sect 3 template */
955 0 : *kjer = 9;
956 0 : g2_free (gfld);
957 0 : return;
958 : }
959 2 : curIndex = 14;
960 40 : for (i = 0; i < gfld->igdtlen; i++) {
961 38 : const struct gridtemplate *templatesgrid = get_templatesgrid();
962 38 : is3[curIndex] = gfld->igdtmpl[i];
963 38 : curIndex += abs (templatesgrid[gridIndex].mapgrid[i]);
964 : }
965 : /* API attempts to return grid in scan mode 0100????. Find the necessary
966 : * indexes into the is3 array for the attempt. */
967 2 : switch (gfld->igdtnum) {
968 : case 0:
969 : case 1:
970 : case 2:
971 : case 3:
972 : case 40:
973 : case 41:
974 : case 42:
975 : case 43:
976 0 : scanIndex = 72 - 1;
977 0 : nxIndex = 31 - 1;
978 0 : nyIndex = 35 - 1;
979 0 : break;
980 : case 10:
981 2 : scanIndex = 60 - 1;
982 2 : nxIndex = 31 - 1;
983 2 : nyIndex = 35 - 1;
984 2 : break;
985 : case 20:
986 : case 30:
987 : case 31:
988 0 : scanIndex = 65 - 1;
989 0 : nxIndex = 31 - 1;
990 0 : nyIndex = 35 - 1;
991 0 : break;
992 : case 90:
993 0 : scanIndex = 64 - 1;
994 0 : nxIndex = 31 - 1;
995 0 : nyIndex = 35 - 1;
996 0 : break;
997 : case 110:
998 0 : scanIndex = 57 - 1;
999 0 : nxIndex = 31 - 1;
1000 0 : nyIndex = 35 - 1;
1001 0 : break;
1002 : case 50:
1003 : case 51:
1004 : case 52:
1005 : case 53:
1006 : case 100:
1007 : case 120:
1008 : case 1000:
1009 : case 1200:
1010 : default:
1011 0 : scanIndex = -1;
1012 0 : nxIndex = -1;
1013 0 : nyIndex = -1;
1014 : }
1015 :
1016 2 : is4[4] = 4;
1017 2 : is4[5] = gfld->num_coord;
1018 2 : is4[7] = gfld->ipdtnum;
1019 2 : pdsIndex = getpdsindex (gfld->ipdtnum);
1020 2 : if (pdsIndex == -1) {
1021 0 : jer[8 + *ndjer] = 2;
1022 0 : jer[8] = 2004; /* undefined sect 4 template */
1023 0 : *kjer = 9;
1024 0 : g2_free (gfld);
1025 0 : return;
1026 : }
1027 2 : curIndex = 9;
1028 60 : for (i = 0; i < gfld->ipdtlen; i++) {
1029 58 : const struct pdstemplate *templatespds = get_templatespds();
1030 58 : is4[curIndex] = gfld->ipdtmpl[i];
1031 58 : curIndex += abs (templatespds[pdsIndex].mappds[i]);
1032 : }
1033 :
1034 2 : is5[4] = 5;
1035 2 : is5[5] = gfld->ndpts;
1036 2 : is5[9] = gfld->idrtnum;
1037 2 : drsIndex = getdrsindex (gfld->idrtnum);
1038 2 : if (drsIndex == -1) {
1039 0 : jer[8 + *ndjer] = 2;
1040 0 : jer[8] = 2005; /* undefined sect 5 template */
1041 0 : *kjer = 9;
1042 0 : g2_free (gfld);
1043 0 : return;
1044 : }
1045 2 : curIndex = 11;
1046 38 : for (i = 0; i < gfld->idrtlen; i++) {
1047 36 : const struct drstemplate *templatesdrs = get_templatesdrs();
1048 36 : is5[curIndex] = gfld->idrtmpl[i];
1049 36 : curIndex += abs (templatesdrs[drsIndex].mapdrs[i]);
1050 : }
1051 : /* Mimic MDL's handling of reference value (IS5(12)) */
1052 2 : memcpy (&f_temp, &(is5[11]), sizeof (float));
1053 2 : is5[11] = (sInt4) f_temp;
1054 2 : if ((is5[9] == 2) || (is5[9] == 3)) {
1055 2 : if (is5[20] == 0) {
1056 2 : memcpy (&(f_temp), &(is5[23]), sizeof (float));
1057 2 : *xmissp = f_temp;
1058 2 : is5[23] = (sInt4) f_temp;
1059 2 : memcpy (&(f_temp), &(is5[27]), sizeof (float));
1060 2 : *xmisss = f_temp;
1061 2 : is5[27] = (sInt4) f_temp;
1062 : } else {
1063 0 : *xmissp = is5[23];
1064 0 : *xmisss = is5[27];
1065 : }
1066 : }
1067 :
1068 2 : is6[4] = 6;
1069 2 : is6[5] = gfld->ibmap;
1070 2 : is7[4] = 7;
1071 :
1072 2 : if (subgNum + 1 == numfields) {
1073 2 : *iendpk = 1;
1074 : } else {
1075 0 : *iendpk = 0;
1076 : }
1077 :
1078 2 : if ((gfld->ibmap == 0) || (gfld->ibmap == 254)) {
1079 0 : *ibitmap = 1;
1080 : } else {
1081 2 : *ibitmap = 0;
1082 : }
1083 :
1084 : /* Check type of original field, before transfering the memory. */
1085 2 : myAssert (*ns5 > 20);
1086 : /* Check if NCEP had problems expanding the data. If so we currently
1087 : * abort. May need to revisit this behavior. */
1088 2 : if (!gfld->expanded) {
1089 0 : jer[0 + *ndjer] = 2;
1090 0 : *kjer = 1;
1091 0 : g2_free (gfld);
1092 0 : return;
1093 : }
1094 2 : f_ignoreScan = 0;
1095 : /* Check if integer type... is5[20] == 1 implies integer (code table
1096 : * 5.1), but only for certain templates. (0,1,2,3,40,41,40000,40010).
1097 : * (not 50,51) */
1098 2 : if ((is5[20] == 1) && ((is5[9] != 50) && (is5[9] != 51))) {
1099 : /* integer data, use iain */
1100 0 : if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) {
1101 0 : ierr = TransferInt (gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap,
1102 : 1, &(dummyScan), 0, 0, *iclean, *xmissp,
1103 : iain, *nd2x3, ib);
1104 : } else {
1105 0 : ierr = TransferInt (gfld->fld, gfld->ngrdpts, *ibitmap, gfld->bmap,
1106 0 : f_ignoreScan, &(is3[scanIndex]), is3[nxIndex],
1107 0 : is3[nyIndex], *iclean, *xmissp, iain, *nd2x3,
1108 : ib);
1109 : }
1110 : } else {
1111 : /* float data, use ain */
1112 2 : if ((scanIndex < 0) || (nxIndex < 0) || (nyIndex < 0)) {
1113 0 : ierr = TransferFloat (gfld->fld, gfld->ngrdpts, *ibitmap,
1114 0 : gfld->bmap, 1, &(dummyScan), 0, 0, *iclean,
1115 : *xmissp, ain, *nd2x3, ib);
1116 : } else {
1117 10 : ierr = TransferFloat (gfld->fld, gfld->ngrdpts, *ibitmap,
1118 4 : gfld->bmap, f_ignoreScan, &(is3[scanIndex]),
1119 4 : is3[nxIndex], is3[nyIndex], *iclean, *xmissp,
1120 : ain, *nd2x3, ib);
1121 : }
1122 : }
1123 2 : if (ierr != 0) {
1124 0 : switch (ierr) {
1125 : case 1: /* nd2x3 is too small */
1126 0 : jer[0 + *ndjer] = 2;
1127 0 : *kjer = 1;
1128 0 : break;
1129 : case 2: /* nx * ny != ngrdpts */
1130 0 : jer[5 + *ndjer] = 2;
1131 0 : *kjer = 6;
1132 0 : break;
1133 : default:
1134 0 : jer[8 + *ndjer] = 2;
1135 0 : jer[8] = 9999; /* Unknown error message. */
1136 0 : *kjer = 9;
1137 : }
1138 0 : g2_free (gfld);
1139 0 : return;
1140 : }
1141 2 : g2_free (gfld);
1142 : }
1143 :
1144 : /*****************************************************************************
1145 : * validate() --
1146 : *
1147 : * Arthur Taylor / MDL
1148 : *
1149 : * PURPOSE
1150 : * Output all the return values from the API to a file. This is primarily
1151 : * for diagnostics purposes.
1152 : *
1153 : * ARGUMENTS
1154 : * filename = File to write the data to. (Input)
1155 : * ain = The data if the original data was float data. (nd2x3) (Input)
1156 : * iain = The data if the original data was integer data. (nd2x3) (Input)
1157 : * nd2x3 = length of ain, iain, ib (Input)
1158 : * idat = MDL local use data (if any were unpacked from Sect 2). (Input)
1159 : * nidat = length of idat. (Input)
1160 : * rdat = floating point local use data (Input)
1161 : * nrdat = length of rdat. (Input)
1162 : * is0, ns0 = Section 0 data, length of is0 (16) (Input)
1163 : * is1, ns1 = Section 1 data, length of is1 (21) (Input)
1164 : * is2, ns2 = Section 2 data, length of is2 (??) (Input)
1165 : * is3, ns3 = Section 3 data, length of is3 (96 or 1600) (Input)
1166 : * is4, ns4 = Section 4 data, length of is4 (60) (Input)
1167 : * is5, ns5 = Section 5 data, length of is5 (49) (Input)
1168 : * is6, ns6 = Section 6 data, length of is6 (6) (Input)
1169 : * is7, ns7 = Section 7 data, length of ns7 (8) (Input)
1170 : * ib = Bitmap if user requested it, and it was packed (Input)
1171 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Input)
1172 : * xmissp = Primary missing value. (Input)
1173 : * xmisss = Secondary missing value (Input)
1174 : * iendpk = flag if there is more grids in the message. (Input)
1175 : * jer(ndjer,2) = error codes along with severity. (Input)
1176 : * ndjer = length of jer. (15) (Input)
1177 : * kjer = number of error messages stored in jer. (Input)
1178 : *
1179 : * FILES/DATABASES: None
1180 : *
1181 : * RETURNS: int
1182 : * 0 = Ok.
1183 : * -1 = Problems opening the file.
1184 : *
1185 : * HISTORY
1186 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1187 : *
1188 : * NOTES
1189 : *****************************************************************************
1190 : */
1191 : #ifdef notdef
1192 : static int validate (char *filename, float * ain, sInt4 * iain,
1193 : sInt4 * nd2x3, sInt4 * idat, sInt4 * nidat,
1194 : float * rdat, sInt4 * nrdat, sInt4 * is0, sInt4 * ns0,
1195 : sInt4 * is1, sInt4 * ns1, sInt4 * is2, sInt4 * ns2,
1196 : sInt4 * is3, sInt4 * ns3, sInt4 * is4, sInt4 * ns4,
1197 : sInt4 * is5, sInt4 * ns5, sInt4 * is6, sInt4 * ns6,
1198 : sInt4 * is7, sInt4 * ns7, sInt4 * ib, sInt4 * ibitmap,
1199 : float * xmissp, float * xmisss, sInt4 * iendpk,
1200 : sInt4 * jer, sInt4 * ndjer, sInt4 * kjer)
1201 : {
1202 : FILE *fp; /* Open file to write to. */
1203 : int i; /* Counter for printing to file. */
1204 :
1205 : if ((fp = fopen (filename, "wt")) == NULL) {
1206 : #ifdef DEBUG
1207 : printf ("unable to open %s for write\n", filename);
1208 : #endif
1209 : return -1;
1210 : }
1211 : for (i = 0; i < *ns0; i++) {
1212 : fprintf (fp, "Sect 0 : %d of %d : %d\n", i, *ns0, is0[i]);
1213 : }
1214 : for (i = 0; i < *ns1; i++) {
1215 : fprintf (fp, "Sect 1 : %d of %d : %d\n", i, *ns1, is1[i]);
1216 : }
1217 : for (i = 0; i < *ns2; i++) {
1218 : fprintf (fp, "Sect 2 : %d of %d : %d\n", i, *ns2, is2[i]);
1219 : }
1220 : for (i = 0; i < idat[0]; i++) {
1221 : fprintf (fp, "idat : %d of %d : %d\n", i, idat[0], idat[i]);
1222 : }
1223 : for (i = 0; i < rdat[0]; i++) {
1224 : fprintf (fp, "rdat : %d of %f : %f\n", i, rdat[0], rdat[i]);
1225 : }
1226 : for (i = 0; i < *ns3; i++) {
1227 : fprintf (fp, "Sect 3 : %d of %d : %d\n", i, *ns3, is3[i]);
1228 : }
1229 : for (i = 0; i < *ns4; i++) {
1230 : fprintf (fp, "Sect 4 : %d of %d : %d\n", i, *ns4, is4[i]);
1231 : }
1232 : for (i = 0; i < *ns5; i++) {
1233 : fprintf (fp, "Sect 5 : %d of %d : %d\n", i, *ns5, is5[i]);
1234 : }
1235 : for (i = 0; i < *ns6; i++) {
1236 : fprintf (fp, "Sect 6 : %d of %d : %d\n", i, *ns6, is6[i]);
1237 : }
1238 : for (i = 0; i < *ns7; i++) {
1239 : fprintf (fp, "Sect 7 : %d of %d : %d\n", i, *ns7, is7[i]);
1240 : }
1241 : fprintf (fp, "Xmissp = %f\n", *xmissp);
1242 : fprintf (fp, "xmisss = %f\n", *xmisss);
1243 : if ((is5[9] == 0) || (is5[9] == 1) || (is5[9] == 2) || (is5[9] == 3)) {
1244 : if (is5[20] == 1) {
1245 : for (i = 0; i < *nd2x3; i++) {
1246 : fprintf (fp, "Int Data : %d of %d : %d\n", i, *nd2x3, iain[i]);
1247 : }
1248 : }
1249 : } else {
1250 : for (i = 0; i < *nd2x3; i++) {
1251 : fprintf (fp, "Float Data : %d of %d : %f\n", i, *nd2x3, ain[i]);
1252 : }
1253 : }
1254 : fprintf (fp, "ibitmap = %d\n", *ibitmap);
1255 : if (*ibitmap) {
1256 : for (i = 0; i < *nd2x3; i++) {
1257 : fprintf (fp, "Bitmap Data : %d of %d : %d\n", i, *nd2x3, ib[i]);
1258 : }
1259 : }
1260 : for (i = 0; i < *ndjer; i++) {
1261 : fprintf (fp, "jer(i,1) %d jer(i,2) %d\n", jer[i], jer[i + *ndjer]);
1262 : }
1263 : fprintf (fp, "kjer = %d\n", *kjer);
1264 : fprintf (fp, "iendpk = %d\n", *iendpk);
1265 : fclose (fp);
1266 : return 0;
1267 : }
1268 : #endif /* def notdef */
1269 :
1270 : /*****************************************************************************
1271 : * clear() --
1272 : *
1273 : * Arthur Taylor / MDL
1274 : *
1275 : * PURPOSE
1276 : * Clear all the return values from the API. This is primarily for
1277 : * diagnostics purposes.
1278 : *
1279 : * ARGUMENTS
1280 : * ain = The data if the original data was float data. (nd2x3) (Output)
1281 : * iain = The data if the original data was integer data. (nd2x3) (Output)
1282 : * nd2x3 = length of ain, iain, ib (Input)
1283 : * idat = MDL local use data (if any were unpacked from Sect 2). (Output)
1284 : * nidat = length of idat. (Input)
1285 : * rdat = floating point local use data (Output)
1286 : * nrdat = length of rdat. (Input)
1287 : * is0, ns0 = Section 0 data, length of is0 (16) (Output) (Input)
1288 : * is1, ns1 = Section 1 data, length of is1 (21) (Output) (Input)
1289 : * is2, ns2 = Section 2 data, length of is2 (??) (Output) (Input)
1290 : * is3, ns3 = Section 3 data, length of is3 (96 or 1600) (Output) (Input)
1291 : * is4, ns4 = Section 4 data, length of is4 (60) (Output) (Input)
1292 : * is5, ns5 = Section 5 data, length of is5 (49) (Output) (Input)
1293 : * is6, ns6 = Section 6 data, length of is6 (6) (Output) (Input)
1294 : * is7, ns7 = Section 7 data, length of ns7 (8) (Output) (Input)
1295 : * ib = Bitmap if user requested it, and it was packed (Output)
1296 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output)
1297 : * xmissp = Primary missing value. (Output)
1298 : * xmisss = Secondary missing value (Output)
1299 : * iendpk = flag if there is more grids in the message. (Output)
1300 : * jer(ndjer,2) = error codes along with severity. (Output)
1301 : * ndjer = length of jer. (15) (Input)
1302 : * kjer = number of error messages stored in jer. (Output)
1303 : *
1304 : * FILES/DATABASES: None
1305 : *
1306 : * RETURNS: int
1307 : * 0 = Ok.
1308 : * -1 = Problems opening the file.
1309 : *
1310 : * HISTORY
1311 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1312 : *
1313 : * NOTES
1314 : *****************************************************************************
1315 : */
1316 : /*
1317 : static void clear (float * ain, sInt4 * iain, sInt4 * nd2x3, sInt4 * idat,
1318 : sInt4 * nidat, float * rdat, sInt4 * nrdat, sInt4 * is0,
1319 : sInt4 * ns0, sInt4 * is1, sInt4 * ns1, sInt4 * is2,
1320 : sInt4 * ns2, sInt4 * is3, sInt4 * ns3, sInt4 * is4,
1321 : sInt4 * ns4, sInt4 * is5, sInt4 * ns5, sInt4 * is6,
1322 : sInt4 * ns6, sInt4 * is7, sInt4 * ns7, sInt4 * ib,
1323 : sInt4 * ibitmap, float * xmissp, float * xmisss,
1324 : sInt4 * iendpk, sInt4 * jer, sInt4 * ndjer, sInt4 * kjer)
1325 : {
1326 : memset ((void *) ain, 0, *nd2x3 * sizeof (float));
1327 : memset ((void *) iain, 0, *nd2x3 * sizeof (sInt4));
1328 : memset ((void *) idat, 0, *nidat * sizeof (sInt4));
1329 : memset ((void *) rdat, 0, *nrdat * sizeof (float));
1330 : memset ((void *) is0, 0, *ns0 * sizeof (sInt4));
1331 : memset ((void *) is1, 0, *ns1 * sizeof (sInt4));
1332 : memset ((void *) is2, 0, *ns2 * sizeof (sInt4));
1333 : memset ((void *) is3, 0, *ns3 * sizeof (sInt4));
1334 : memset ((void *) is4, 0, *ns4 * sizeof (sInt4));
1335 : if ((is5[9] == 2) || (is5[9] == 3)) {
1336 : *xmissp = 0;
1337 : *xmisss = 0;
1338 : }
1339 : memset ((void *) is5, 0, *ns5 * sizeof (sInt4));
1340 : memset ((void *) is6, 0, *ns6 * sizeof (sInt4));
1341 : memset ((void *) is7, 0, *ns7 * sizeof (sInt4));
1342 : memset ((void *) ib, 0, *nd2x3 * sizeof (sInt4));
1343 : *ibitmap = 0;
1344 : *iendpk = 0;
1345 : memset ((void *) jer, 0, 2 * *ndjer * sizeof (sInt4));
1346 : *kjer = 0;
1347 : }
1348 : */
1349 :
1350 : /*****************************************************************************
1351 : * BigByteCpy() --
1352 : *
1353 : * Arthur Taylor / MDL
1354 : *
1355 : * PURPOSE
1356 : * This is so we can copy upto 4 bytes from a big endian 4 byte int data
1357 : * stream.
1358 : *
1359 : * The reason this is needed is because the GRIB2 API required the GRIB2
1360 : * message to be passed in as a big endian 4 byte int data stream instead of
1361 : * something more reasonable (such as a stream of 1 byte char's) By having
1362 : * this routine we reduce the number of memswaps of the message on a little
1363 : * endian system.
1364 : *
1365 : * ARGUMENTS
1366 : * dst = Where to copy the data to. (Output)
1367 : * ipack = The 4 byte int data stream. (Input)
1368 : * nd5 = length of ipack (Input)
1369 : * startInt = Which integer to start reading in ipack. (Input)
1370 : * startByte = Which byte in the startInt to start reading from. (Input)
1371 : * numByte = How many bytes to read. (Input)
1372 : *
1373 : * FILES/DATABASES: None
1374 : *
1375 : * RETURNS: NULL
1376 : *
1377 : * HISTORY
1378 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1379 : *
1380 : * NOTES
1381 : *****************************************************************************
1382 : */
1383 32 : static void BigByteCpy (sInt4 * dst, sInt4 * ipack, sInt4 nd5,
1384 : unsigned int startInt, unsigned int startByte,
1385 : int numByte)
1386 : {
1387 : static int Lshift[] = { 0, 8, 16, 24 }; /* Amounts to shift left by. */
1388 : unsigned int intIndex; /* Where in ipack to read from. */
1389 : unsigned int byteIndex; /* Where in intIndex to read from. */
1390 : uInt4 curInt; /* An unsigned version of the current int. */
1391 : unsigned int curByte; /* The current byte we have read. */
1392 : int i; /* Loop counter over number of bytes to read. */
1393 :
1394 32 : myAssert (numByte <= 4);
1395 32 : myAssert (startByte < 4);
1396 32 : *dst = 0;
1397 32 : intIndex = startInt;
1398 32 : byteIndex = startByte;
1399 106 : for (i = 0; i < numByte; i++) {
1400 74 : curInt = (uInt4) ipack[intIndex];
1401 74 : curByte = (curInt << Lshift[byteIndex]) >> 24;
1402 74 : *dst = (*dst << 8) + curByte;
1403 74 : byteIndex++;
1404 74 : if (byteIndex == 4) {
1405 16 : byteIndex = 0;
1406 16 : intIndex++;
1407 : }
1408 : }
1409 32 : }
1410 :
1411 : /*****************************************************************************
1412 : * FindTemplateIDs() --
1413 : *
1414 : * Arthur Taylor / MDL
1415 : *
1416 : * PURPOSE
1417 : * This is so we can find out which templates are in this GRIB2 message.
1418 : * Which allows us to determine if we can call the MDL routines, or if we
1419 : * should call the NCEP routines instead.
1420 : *
1421 : * ARGUMENTS
1422 : * ipack = The 4 byte int data stream. (Input)
1423 : * nd5 = length of ipack (Input)
1424 : * subgNum = Which subgrid to look at. (Input)
1425 : * gdsTmpl = The gds template number for this subgrid. (Output)
1426 : * pdsTmpl = The pds template number for this subgrid. (Output)
1427 : * drsTmpl = The drs template number for this subgrid. (Output)
1428 : * f_noBitmap = 0 has a bitmap, else doesn't have a bitmap. (Output)
1429 : * orderDiff = The order of the differencing in template 5.3 (1st, 2nd) (out)
1430 : *
1431 : * FILES/DATABASES: None
1432 : *
1433 : * RETURNS: int
1434 : * 0 = Ok.
1435 : * 2 = invalid section number
1436 : *
1437 : * HISTORY
1438 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1439 : *
1440 : * NOTES
1441 : *****************************************************************************
1442 : */
1443 2 : static int FindTemplateIDs (sInt4 * ipack, sInt4 nd5, int subgNum,
1444 : sInt4 * gdsTmpl, sInt4 * pdsTmpl,
1445 : sInt4 * drsTmpl, sInt4 * numGrps,
1446 : uChar * f_noBitmap, sInt4 * orderDiff)
1447 : {
1448 2 : unsigned int gNum = 0; /* Which sub group we are currently working with. */
1449 : sInt4 offset; /* Where in the bytes stream we currently are. */
1450 : sInt4 sectLen; /* The length of the current section. */
1451 : sInt4 sectId; /* The current section number. */
1452 : sInt4 li_temp; /* A temporary holder for the bitmap flag. */
1453 :
1454 : /* Jump over section 0. */
1455 2 : offset = 16;
1456 14 : while (gNum <= subgNum) {
1457 10 : BigByteCpy (§Len, ipack, nd5, (offset / 4), (offset % 4), 4);
1458 : /* Check if we just read section 8. If so, then it is "7777" =
1459 : * 926365495 regardless of endian'ness. */
1460 10 : if (sectLen == 926365495L) {
1461 : #ifdef DEBUG
1462 0 : printf ("Shouldn't see sect 8. Should stop after correct sect 5\n");
1463 : #endif
1464 0 : return 2;
1465 : }
1466 10 : BigByteCpy (§Id, ipack, nd5, ((offset + 4) / 4),
1467 10 : ((offset + 4) % 4), 1);
1468 10 : switch (sectId) {
1469 : case 1:
1470 : case 2:
1471 2 : break;
1472 : case 3:
1473 2 : BigByteCpy (gdsTmpl, ipack, nd5, ((offset + 12) / 4),
1474 2 : ((offset + 12) % 4), 2);
1475 2 : break;
1476 : case 4:
1477 2 : BigByteCpy (pdsTmpl, ipack, nd5, ((offset + 7) / 4),
1478 2 : ((offset + 7) % 4), 2);
1479 2 : break;
1480 : case 5:
1481 2 : BigByteCpy (drsTmpl, ipack, nd5, ((offset + 9) / 4),
1482 2 : ((offset + 9) % 4), 2);
1483 4 : if ((*drsTmpl == 2) || (*drsTmpl == 3)) {
1484 2 : BigByteCpy (numGrps, ipack, nd5, ((offset + 31) / 4),
1485 2 : ((offset + 31) % 4), 4);
1486 : } else {
1487 0 : *numGrps = 0;
1488 : }
1489 2 : if (*drsTmpl == 3) {
1490 2 : BigByteCpy (&li_temp, ipack, nd5, ((offset + 44) / 4),
1491 2 : ((offset + 44) % 4), 1);
1492 2 : *orderDiff = li_temp;
1493 : } else {
1494 0 : *orderDiff = 0;
1495 : }
1496 2 : break;
1497 : case 6:
1498 2 : BigByteCpy (&li_temp, ipack, nd5, ((offset + 5) / 4),
1499 2 : ((offset + 5) % 4), 1);
1500 2 : if (li_temp == 255) {
1501 2 : *f_noBitmap = 1;
1502 : }
1503 2 : gNum++;
1504 2 : break;
1505 : case 7:
1506 0 : break;
1507 : default:
1508 : #ifdef DEBUG
1509 0 : printf ("Invalid section id %d.\n", sectId);
1510 : #endif
1511 0 : return 2;
1512 : }
1513 10 : offset += sectLen;
1514 : }
1515 2 : return 0;
1516 : }
1517 :
1518 : /*****************************************************************************
1519 : * unpk_grib2() --
1520 : *
1521 : * Arthur Taylor / MDL
1522 : *
1523 : * PURPOSE
1524 : * This procedure is the main API for decoding GRIB2 messages. It is
1525 : * intended to be a branching routine to call either the MDL GRIB2 library,
1526 : * or the NCEP GRIB2 library, depending on which template it sees in the
1527 : * message.
1528 : *
1529 : * ARGUMENTS
1530 : * kfildo = Unit number for output diagnostics (C ignores this). (Input)
1531 : * ain = contains the data if the original data was float data.
1532 : * (size = nd2x3) (Output)
1533 : * iain = contains the data if the original data was integer data.
1534 : * (size = nd2x3) (Output)
1535 : * nd2x3 = length of ain, iain, ib (Input)
1536 : * (at least the size of num grid points)
1537 : * idat = local use data if any that were unpacked from Section 2. (Output)
1538 : * nidat = length of idat. (Input)
1539 : * rdat = floating point local use data (Output)
1540 : * nrdat = length of rdat. (Input)
1541 : * is0 = Section 0 data (Output)
1542 : * ns0 = length of is0 (16 is fine) (Input)
1543 : * is1 = Section 1 data (Output)
1544 : * ns1 = length of is1 (21 is fine) (Input)
1545 : * is2 = Section 2 data (Output)
1546 : * ns2 = length of is2 () (Input)
1547 : * is3 = Section 3 data (Output)
1548 : * ns3 = length of is3 (96 or 1600) (Input)
1549 : * is4 = Section 4 data (Output)
1550 : * ns4 = length of is4 (60) (Input)
1551 : * is5 = Section 5 data (Output)
1552 : * ns5 = length of is5 (49 is fine) (Input)
1553 : * is6 = Section 6 data (Output)
1554 : * ns6 = length of is6 (6 is fine) (Input)
1555 : * is7 = Section 7 data (Output)
1556 : * ns7 = length of is7 (8 is fine) (Input)
1557 : * ib = Bitmap if user requested it, and it was packed (Output)
1558 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Output)
1559 : * ipack = The message to unpack (This is assumed to be Big endian) (Input)
1560 : * nd5 = The size of ipack (Input)
1561 : * xmissp = The floating point representation for the primary missing value.
1562 : * (Input/Output)
1563 : * xmisss = The floating point representation for the secondary missing
1564 : * value (Input/Output)
1565 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
1566 : * iclean = 1 means the user wants the unpacked data returned without
1567 : * missing values in it. 0 means embed the missing values. (Input)
1568 : * l3264b = Integer word length in bits (32 or 64) (Input)
1569 : * iendpk = A 1 means no more grids in this message, a 0 means more grids.
1570 : * (Output)
1571 : * jer(ndjer,2) = error codes along with severity. (Output)
1572 : * ndjer = 1/2 length of jer. (>= 15) (Input)
1573 : * kjer = number of error messages stored in jer. (Output)
1574 : *
1575 : * FILES/DATABASES: None
1576 : *
1577 : * RETURNS: void
1578 : *
1579 : * HISTORY
1580 : * 12/2003 Arthur Taylor (MDL/RSIS): Created.
1581 : *
1582 : * NOTES
1583 : * Inefficiencies: have to memswap ipack multiple times.
1584 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
1585 : * MDL attempts to always return grids in scan mode 0100????
1586 : * ToDo: Check length of parameters better.
1587 : *
1588 : * According to MDL's unpk_grib2.f, it currently supports (so for others we
1589 : * call the NCEP routines):
1590 : * TEMPLATE 3.0 EQUIDISTANT CYLINDRICAL LATITUDE/LONGITUDE
1591 : * TEMPLATE 3.10 MERCATOR
1592 : * TEMPLATE 3.20 POLAR STEREOGRAPHIC
1593 : * TEMPLATE 3.30 LAMBERT
1594 : * TEMPLATE 3.90 ORTHOGRAPHIC SPACE VIEW
1595 : * TEMPLATE 3.110 EQUATORIAL AZIMUTHAL EQUIDISTANT
1596 : * TEMPLATE 3.120 AZIMUTH-RANGE (RADAR)
1597 : *
1598 : * TEMPLATE 4.0 ANALYSIS OR FORECAST AT A LEVEL AND POINT
1599 : * TEMPLATE 4.1 INDIVIDUAL ENSEMBLE
1600 : * TEMPLATE 4.2 DERIVED FORECAST BASED ON ENSEMBLES
1601 : * TEMPLATE 4.8 AVERAGE, ACCUMULATION, EXTREMES
1602 : * TEMPLATE 4.20 RADAR
1603 : * TEMPLATE 4.30 SATELLITE
1604 : *
1605 : * TEMPLATE 5.0 SIMPLE PACKING
1606 : * TEMPLATE 5.2 COMPLEX PACKING
1607 : * TEMPLATE 5.3 COMPLEX PACKING AND SPATIAL DIFFERENCING
1608 : *
1609 : * Correction to "unpk_grib2.f" : It also supports:
1610 : * TEMPLATE 4.9 Probability forecast in a time interval
1611 : *
1612 : *****************************************************************************
1613 : */
1614 2 : void unpk_grib2 (sInt4 * kfildo, float * ain, sInt4 * iain, sInt4 * nd2x3,
1615 : sInt4 * idat, sInt4 * nidat, float * rdat, sInt4 * nrdat,
1616 : sInt4 * is0, sInt4 * ns0, sInt4 * is1, sInt4 * ns1,
1617 : sInt4 * is2, sInt4 * ns2, sInt4 * is3, sInt4 * ns3,
1618 : sInt4 * is4, sInt4 * ns4, sInt4 * is5, sInt4 * ns5,
1619 : sInt4 * is6, sInt4 * ns6, sInt4 * is7, sInt4 * ns7,
1620 : sInt4 * ib, sInt4 * ibitmap, sInt4 * ipack, sInt4 * nd5,
1621 : float * xmissp, float * xmisss, sInt4 * inew,
1622 : sInt4 * iclean, sInt4 * l3264b, sInt4 * iendpk, sInt4 * jer,
1623 : sInt4 * ndjer, sInt4 * kjer)
1624 : {
1625 : unsigned char *c_ipack; /* The compressed data as char instead of sInt4
1626 : * so it is easier to work with. */
1627 : sInt4 gdsTmpl;
1628 : sInt4 pdsTmpl;
1629 : sInt4 drsTmpl;
1630 : sInt4 numGrps;
1631 2 : char f_useMDL = 0; /* Instructed 3/8/2005 10:30 to not use MDL. */
1632 : uChar f_noBitmap; /* 0 if bitmap, else no bitmap. */
1633 : sInt4 orderDiff;
1634 :
1635 2 : if (FindTemplateIDs (ipack, *nd5, 0, &gdsTmpl, &pdsTmpl, &drsTmpl,
1636 : &numGrps, &f_noBitmap, &orderDiff) != 0) {
1637 0 : jer[0 + *ndjer] = 2;
1638 0 : jer[0] = 3000; /* Couldn't figure out which templates are used. */
1639 0 : *kjer = 1;
1640 : }
1641 2 : if ((gdsTmpl != 0) && (gdsTmpl != 10) && (gdsTmpl != 20) &&
1642 0 : (gdsTmpl != 30) && (gdsTmpl != 90) && (gdsTmpl != 110) &&
1643 0 : (gdsTmpl != 120)) {
1644 0 : f_useMDL = 0;
1645 : }
1646 4 : if ((pdsTmpl != 0) && (pdsTmpl != 1) && (pdsTmpl != 2) &&
1647 2 : (pdsTmpl != 8) && (pdsTmpl != 9) && (pdsTmpl != 20) &&
1648 0 : (pdsTmpl != 30)) {
1649 0 : f_useMDL = 0;
1650 : }
1651 2 : if ((drsTmpl != 0) && (drsTmpl != 2) && (drsTmpl != 3)) {
1652 0 : f_useMDL = 0;
1653 : }
1654 : /* MDL GRIB2 lib does not support drsTmpl 2 or 3 if there is a bitmap. */
1655 2 : if ((!f_noBitmap) && ((drsTmpl == 2) || (drsTmpl == 3))) {
1656 0 : f_useMDL = 0;
1657 : }
1658 : /* MDL GRIB2 lib does not support anything but second order differencing. */
1659 2 : if ((drsTmpl == 3) && (orderDiff != 2) && (orderDiff != 0)) {
1660 2 : f_useMDL = 0;
1661 : }
1662 :
1663 : #ifdef _FORTRAN
1664 : if (f_useMDL) {
1665 : jmin = (sInt4 *) malloc (numGrps * sizeof (sInt4));
1666 : lbit = (sInt4 *) malloc (numGrps * sizeof (sInt4));
1667 : nov = (sInt4 *) malloc (numGrps * sizeof (sInt4));
1668 : iwork = (sInt4 *) malloc ((*nd2x3) * sizeof (sInt4));
1669 :
1670 : UNPK_G2MDL (kfildo, jmin, lbit, nov, iwork, ain, iain, nd2x3, idat,
1671 : nidat, rdat, nrdat, is0, ns0, is1, ns1, is2, ns2, is3, ns3,
1672 : is4, ns4, is5, ns5, is6, ns6, is7, ns7, ib, ibitmap, ipack,
1673 : nd5, xmissp, xmisss, inew, iclean, l3264b, iendpk, jer,
1674 : ndjer, kjer);
1675 : /*
1676 : if (1==1) {
1677 : int i;
1678 : for (i=0; i < *nd2x3; i++) {
1679 : if (ain[i] != 9999) {
1680 : fprintf (stderr, "here grib2api.c: ain[%d] %f\n", i, ain[i]);
1681 : }
1682 : }
1683 : }
1684 : */
1685 : free (jmin);
1686 : free (lbit);
1687 : free (nov);
1688 : free (iwork);
1689 :
1690 : /* Check the behavior of the NCEP routines. */
1691 : /*
1692 : #ifdef DEBUG
1693 : validate ("check2.txt", ain, iain, nd2x3, idat, nidat, rdat, nrdat,
1694 : is0, ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5,
1695 : is6, ns6, is7, ns7, ib, ibitmap, xmissp, xmisss, iendpk,
1696 : jer, ndjer, kjer);
1697 : clear (ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0, ns0, is1, ns1,
1698 : is2, ns2, is3, ns3, is4, ns4, is5, ns5, is6, ns6, is7, ns7, ib,
1699 : ibitmap, xmissp, xmisss, iendpk, jer, ndjer, kjer);
1700 : #ifdef LITTLE_ENDIAN
1701 : memswp (ipack, sizeof (sInt4), *nd5);
1702 : #endif
1703 : c_ipack = (unsigned char *) ipack;
1704 : unpk_g2ncep (kfildo, ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0,
1705 : ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5,
1706 : is6, ns6, is7, ns7, ib, ibitmap, c_ipack, nd5, xmissp,
1707 : xmisss, inew, iclean, l3264b, iendpk, jer, ndjer, kjer);
1708 : #ifdef LITTLE_ENDIAN
1709 : memswp (ipack, sizeof (sInt4), *nd5);
1710 : #endif
1711 : validate ("check1.txt", ain, iain, nd2x3, idat, nidat, rdat, nrdat,
1712 : is0, ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5,
1713 : is6, ns6, is7, ns7, ib, ibitmap, xmissp, xmisss, iendpk, jer,
1714 : ndjer, kjer);
1715 : #endif
1716 : */
1717 : } else {
1718 : #endif
1719 : /* Since NCEP's code works with byte level stuff, we need to un-do the
1720 : * byte swap of the (int *) data, then cast it to an (unsigned char *) */
1721 : #ifdef LITTLE_ENDIAN
1722 : /* Can't make this dependent on inew, since they could have a sequence
1723 : * of get first message... do stuff, get second message, which
1724 : * unfortunately means they would have to get the first message again,
1725 : * causing 2 swaps, and breaking their second request for the first
1726 : * message (as well as their second message). */
1727 2 : memswp (ipack, sizeof (sInt4), *nd5);
1728 : #endif
1729 2 : c_ipack = (unsigned char *) ipack;
1730 2 : unpk_g2ncep (kfildo, ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0,
1731 : ns0, is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5,
1732 : is6, ns6, is7, ns7, ib, ibitmap, c_ipack, nd5, xmissp,
1733 : xmisss, inew, iclean, l3264b, iendpk, jer, ndjer, kjer);
1734 :
1735 : #ifdef LITTLE_ENDIAN
1736 : /* Swap back because we could be called again for the subgrid data. */
1737 2 : memswp (ipack, sizeof (sInt4), *nd5);
1738 : #endif
1739 : #ifdef _FORTRAN
1740 : }
1741 : #endif
1742 2 : }
1743 :
1744 : /* Not sure I need this... It is intended to provide a way to call it from
1745 : * FORTRAN, but I'm not sure it is needed. */
1746 : /* gcc has two __ if there is one _ in the procedure name. */
1747 0 : void unpk_grib2__ (sInt4 * kfildo, float * ain, sInt4 * iain,
1748 : sInt4 * nd2x3, sInt4 * idat, sInt4 * nidat, float * rdat,
1749 : sInt4 * nrdat, sInt4 * is0, sInt4 * ns0, sInt4 * is1,
1750 : sInt4 * ns1, sInt4 * is2, sInt4 * ns2, sInt4 * is3,
1751 : sInt4 * ns3, sInt4 * is4, sInt4 * ns4, sInt4 * is5,
1752 : sInt4 * ns5, sInt4 * is6, sInt4 * ns6, sInt4 * is7,
1753 : sInt4 * ns7, sInt4 * ib, sInt4 * ibitmap, sInt4 * ipack,
1754 : sInt4 * nd5, float * xmissp, float * xmisss,
1755 : sInt4 * inew, sInt4 * iclean, sInt4 * l3264b,
1756 : sInt4 * iendpk, sInt4 * jer, sInt4 * ndjer, sInt4 * kjer)
1757 : {
1758 0 : unpk_grib2 (kfildo, ain, iain, nd2x3, idat, nidat, rdat, nrdat, is0, ns0,
1759 : is1, ns1, is2, ns2, is3, ns3, is4, ns4, is5, ns5, is6, ns6,
1760 : is7, ns7, ib, ibitmap, ipack, nd5, xmissp, xmisss, inew,
1761 : iclean, l3264b, iendpk, jer, ndjer, kjer);
1762 0 : }
1763 :
1764 : /*****************************************************************************
1765 : * C_pkGrib2() --
1766 : *
1767 : * Arthur Taylor / MDL
1768 : *
1769 : * PURPOSE
1770 : * This procedure is the main API for encoding GRIB2 messages. It is
1771 : * intended to call NCEP's GRIB2 library.
1772 : *
1773 : * ARGUMENTS
1774 : * ain = Contains the data to pack if the original data was float data.
1775 : * (size = nd2x3) (Input)
1776 : * iain = Contains the data to pack if the original data was integer data.
1777 : * (size = nd2x3) (Input)
1778 : * nx = The number of rows in the gridded product. (Input)
1779 : * ny = The number of columns in the gridded products. (Input)
1780 : * idat = local use data if it is integer (stored in section 2). (Input)
1781 : * nidat = length of idat. (Input)
1782 : * rdat = local use data if it is a float (Input)
1783 : * nrdat = length of rdat. (Input)
1784 : * is0 = Section 0 data (element 7 should be set by caller) (Input/Output)
1785 : * ns0 = length of is0 (16 is fine) (Input)
1786 : * is1 = Section 1 data (Input/Output)
1787 : * ns1 = length of is1 (21 is fine) (Input)
1788 : * is3 = Section 3 data (Input/Output)
1789 : * ns3 = length of is3 (96 or 1600) (Input)
1790 : * is4 = Section 4 data (Input/Output)
1791 : * ns4 = length of is4 (60) (Input)
1792 : * is5 = Section 5 data (Input/Output)
1793 : * ns5 = length of is5 (49 is fine) (Input)
1794 : * is6 = Section 6 data (Input/Output)
1795 : * ns6 = length of is6 (6 is fine) (Input)
1796 : * is7 = Section 7 data (Input/Output)
1797 : * ns7 = length of is7 (8 is fine) (Input)
1798 : * ib = Bitmap if user requested it, and it was packed (Input/Output)
1799 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Input/Output)
1800 : * ipack = The message to unpack (This is assumed to be Big endian) (Output)
1801 : * nd5 = The size of ipack (250 + NX*NY + (NX*NY)/8 + num local data) (In)
1802 : * missp = The integer representation for the primary missing value. (Input)
1803 : * xmissp = The float representation for the primary missing value. (Input)
1804 : * misss = The integer representation for the secondary missing value. (In)
1805 : * xmisss = The float representation for the secondary missing value (Input)
1806 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
1807 : * minpk = minimum size of groups in complex and second order differencing
1808 : * methods. Recommended value 14 (Input)
1809 : * iclean = 1 means no primary missing values embedded in the data field,
1810 : * 0 means there are primary missing values in the data (Input/Out)
1811 : * l3264b = Integer word length in bits (32 or 64) (Input)
1812 : * jer(ndjer,2) = error codes along with severity. (Output)
1813 : * ndjer = 1/2 length of jer. (>= 15) (Input)
1814 : * kjer = number of error messages stored in jer. (Output)
1815 : *
1816 : * FILES/DATABASES: None
1817 : *
1818 : * RETURNS: void
1819 : *
1820 : * HISTORY
1821 : * 1/2004 Arthur Taylor (MDL/RSIS): Created.
1822 : *
1823 : * NOTES
1824 : * Inefficiencies: have to memswap ipack multiple times.
1825 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
1826 : * MDL attempts to always return grids in scan mode 0100????
1827 : * ToDo: Check length of parameters better.
1828 : *
1829 : * According to MDL's pk_grib2.f, it currently supports (so for others we
1830 : * call the NCEP routines):
1831 : * TEMPLATE 3.0 EQUIDISTANT CYLINDRICAL LATITUDE/LONGITUDE
1832 : * TEMPLATE 3.10 MERCATOR
1833 : * TEMPLATE 3.20 POLAR STEREOGRAPHIC
1834 : * TEMPLATE 3.30 LAMBERT
1835 : * TEMPLATE 3.90 ORTHOGRAPHIC SPACE VIEW
1836 : * TEMPLATE 3.110 EQUATORIAL AZIMUTHAL EQUIDISTANT
1837 : * TEMPLATE 3.120 AZIMUTH-RANGE (RADAR)
1838 : *
1839 : * TEMPLATE 4.0 ANALYSIS OR FORECAST AT A LEVEL AND POINT
1840 : * TEMPLATE 4.1 INDIVIDUAL ENSEMBLE
1841 : * TEMPLATE 4.2 DERIVED FORECAST BASED ON ENSEMBLES
1842 : * TEMPLATE 4.8 AVERAGE, ACCUMULATION, EXTREMES
1843 : * TEMPLATE 4.20 RADAR
1844 : * TEMPLATE 4.30 SATELLITE
1845 : *
1846 : * TEMPLATE 5.0 SIMPLE PACKING
1847 : * TEMPLATE 5.2 COMPLEX PACKING
1848 : * TEMPLATE 5.3 COMPLEX PACKING AND SPATIAL DIFFERENCING
1849 : *
1850 : * Correction to "pk_grib2.f" : It also supports:
1851 : * TEMPLATE 4.9 Probability forecast in a time interval
1852 : *
1853 : *****************************************************************************
1854 : */
1855 0 : int C_pkGrib2 (unsigned char *cgrib, sInt4 *sec0, sInt4 *sec1,
1856 : unsigned char *csec2, sInt4 lcsec2,
1857 : sInt4 *igds, sInt4 *igdstmpl, sInt4 *ideflist,
1858 : sInt4 idefnum, sInt4 ipdsnum, sInt4 *ipdstmpl,
1859 : float *coordlist, sInt4 numcoord, sInt4 idrsnum,
1860 : sInt4 *idrstmpl, float *fld, sInt4 ngrdpts,
1861 : sInt4 ibmap, sInt4 *bmap)
1862 : {
1863 : int ierr; /* error value from grib2 library. */
1864 :
1865 0 : if ((ierr = g2_create (cgrib, sec0, sec1)) == -1) {
1866 : /* Tried to use for version other than GRIB Ed 2 */
1867 0 : return -1;
1868 : }
1869 :
1870 0 : if ((ierr = g2_addlocal (cgrib, csec2, lcsec2)) < 0) {
1871 : /* Some how got a bad section 2. Should be impossible unless an
1872 : * assert was broken. */
1873 0 : return -2;
1874 : }
1875 :
1876 0 : if ((ierr = g2_addgrid (cgrib, igds, igdstmpl, ideflist, idefnum)) < 0) {
1877 : /* Some how got a bad section 3. Only way would be should be with an
1878 : * unsupported template number unless an assert was broken. */
1879 0 : return -3;
1880 : }
1881 :
1882 0 : if ((ierr = g2_addfield (cgrib, ipdsnum, ipdstmpl, coordlist, numcoord,
1883 : idrsnum, idrstmpl, fld, ngrdpts, ibmap,
1884 : bmap)) < 0) {
1885 0 : return -4;
1886 : }
1887 :
1888 0 : if ((ierr = g2_gribend (cgrib)) < 0) {
1889 0 : return -5;
1890 : }
1891 :
1892 0 : return ierr;
1893 : }
1894 :
1895 : /*****************************************************************************
1896 : * pk_grib2() --
1897 : *
1898 : * Arthur Taylor / MDL
1899 : *
1900 : * PURPOSE
1901 : * This procedure is the main API for encoding GRIB2 messages. It is
1902 : * intended to be a branching routine to call either the MDL GRIB2 library,
1903 : * or the NCEP GRIB2 library, depending on which template it sees in the
1904 : * message.
1905 : * Currently it calls both for debug purposes.
1906 : *
1907 : * ARGUMENTS
1908 : * kfildo = Unit number for output diagnostics (C ignores this). (Input)
1909 : * ain = Contains the data to pack if the original data was float data.
1910 : * (size = nd2x3) (Input)
1911 : * iain = Contains the data to pack if the original data was integer data.
1912 : * (size = nd2x3) (Input)
1913 : * nx = The number of rows in the gridded product. (Input)
1914 : * ny = The number of columns in the gridded products. (Input)
1915 : * idat = local use data if it is integer (stored in section 2). (Input)
1916 : * nidat = length of idat. (Input)
1917 : * rdat = local use data if it is a float (Input)
1918 : * nrdat = length of rdat. (Input)
1919 : * is0 = Section 0 data (element 7 should be set by caller) (Input/Output)
1920 : * ns0 = length of is0 (16 is fine) (Input)
1921 : * is1 = Section 1 data (Input/Output)
1922 : * ns1 = length of is1 (21 is fine) (Input)
1923 : * is3 = Section 3 data (Input/Output)
1924 : * ns3 = length of is3 (96 or 1600) (Input)
1925 : * is4 = Section 4 data (Input/Output)
1926 : * ns4 = length of is4 (60) (Input)
1927 : * is5 = Section 5 data (Input/Output)
1928 : * ns5 = length of is5 (49 is fine) (Input)
1929 : * is6 = Section 6 data (Input/Output)
1930 : * ns6 = length of is6 (6 is fine) (Input)
1931 : * is7 = Section 7 data (Input/Output)
1932 : * ns7 = length of is7 (8 is fine) (Input)
1933 : * ib = Bitmap if user requested it, and it was packed (Input/Output)
1934 : * ibitmap = 0 means ib is invalid, 1 means ib is valid. (Input/Output)
1935 : * ipack = The message to unpack (This is assumed to be Big endian) (Output)
1936 : * nd5 = The size of ipack (250 + NX*NY + (NX*NY)/8 + num local data) (In)
1937 : * missp = The integer representation for the primary missing value. (Input)
1938 : * xmissp = The float representation for the primary missing value. (Input)
1939 : * misss = The integer representation for the secondary missing value. (In)
1940 : * xmisss = The float representation for the secondary missing value (Input)
1941 : * new = 1 if this is the first grid to be unpacked, else 0. (Input)
1942 : * minpk = minimum size of groups in complex and second order differencing
1943 : * methods. Recommended value 14 (Input)
1944 : * iclean = 1 means no primary missing values embedded in the data field,
1945 : * 0 means there are primary missing values in the data (Input/Out)
1946 : * l3264b = Integer word length in bits (32 or 64) (Input)
1947 : * jer(ndjer,2) = error codes along with severity. (Output)
1948 : * ndjer = 1/2 length of jer. (>= 15) (Input)
1949 : * kjer = number of error messages stored in jer. (Output)
1950 : *
1951 : * FILES/DATABASES: None
1952 : *
1953 : * RETURNS: void
1954 : *
1955 : * HISTORY
1956 : * 1/2004 Arthur Taylor (MDL/RSIS): Created.
1957 : *
1958 : * NOTES
1959 : * Inefficiencies: have to memswap ipack multiple times.
1960 : * MDL handles is5[12], is5[23], and is5[27] in an "interesting" manner.
1961 : * MDL attempts to always return grids in scan mode 0100????
1962 : * ToDo: Check length of parameters better.
1963 : *
1964 : * According to MDL's pk_grib2.f, it currently supports (so for others we
1965 : * call the NCEP routines):
1966 : * TEMPLATE 3.0 EQUIDISTANT CYLINDRICAL LATITUDE/LONGITUDE
1967 : * TEMPLATE 3.10 MERCATOR
1968 : * TEMPLATE 3.20 POLAR STEREOGRAPHIC
1969 : * TEMPLATE 3.30 LAMBERT
1970 : * TEMPLATE 3.90 ORTHOGRAPHIC SPACE VIEW
1971 : * TEMPLATE 3.110 EQUATORIAL AZIMUTHAL EQUIDISTANT
1972 : * TEMPLATE 3.120 AZIMUTH-RANGE (RADAR)
1973 : *
1974 : * TEMPLATE 4.0 ANALYSIS OR FORECAST AT A LEVEL AND POINT
1975 : * TEMPLATE 4.1 INDIVIDUAL ENSEMBLE
1976 : * TEMPLATE 4.2 DERIVED FORECAST BASED ON ENSEMBLES
1977 : * TEMPLATE 4.8 AVERAGE, ACCUMULATION, EXTREMES
1978 : * TEMPLATE 4.20 RADAR
1979 : * TEMPLATE 4.30 SATELLITE
1980 : *
1981 : * TEMPLATE 5.0 SIMPLE PACKING
1982 : * TEMPLATE 5.2 COMPLEX PACKING
1983 : * TEMPLATE 5.3 COMPLEX PACKING AND SPATIAL DIFFERENCING
1984 : *
1985 : * Correction to "pk_grib2.f" : It also supports:
1986 : * TEMPLATE 4.9 Probability forecast in a time interval
1987 : *
1988 : *****************************************************************************
1989 : */
1990 0 : void pk_grib2 (sInt4 * kfildo, float * ain, sInt4 * iain, sInt4 * nx,
1991 : sInt4 * ny, sInt4 * idat, sInt4 * nidat, float * rdat,
1992 : sInt4 * nrdat, sInt4 * is0, sInt4 * ns0, sInt4 * is1,
1993 : sInt4 * ns1, sInt4 * is3, sInt4 * ns3, sInt4 * is4,
1994 : sInt4 * ns4, sInt4 * is5, sInt4 * ns5, sInt4 * is6,
1995 : sInt4 * ns6, sInt4 * is7, sInt4 * ns7, sInt4 * ib,
1996 : sInt4 * ibitmap, sInt4 * ipack, sInt4 * nd5, sInt4 * missp,
1997 : float * xmissp, sInt4 * misss, float * xmisss, sInt4 * inew,
1998 : sInt4 * minpk, sInt4 * iclean, sInt4 * l3264b, sInt4 * jer,
1999 : sInt4 * ndjer, sInt4 * kjer)
2000 : {
2001 : #ifndef _FORTRAN
2002 :
2003 0 : printf ("Can not pack things unless using FORTRAN!\n");
2004 : return;
2005 :
2006 : #else
2007 :
2008 : sInt4 gdsTmpl;
2009 : sInt4 pdsTmpl;
2010 : sInt4 drsTmpl;
2011 : sInt4 *jmax;
2012 : sInt4 *jmin;
2013 : sInt4 *lbit;
2014 : sInt4 *nov;
2015 : sInt4 *misslx;
2016 : sInt4 *newbox;
2017 : sInt4 *newboxp;
2018 : sInt4 *ia;
2019 : /* float *a;*/
2020 : char f_useMDL = 1;
2021 : int i;
2022 :
2023 : myAssert (*ndjer >= 8);
2024 : /* Init the error handling array. */
2025 : memset ((void *) jer, 0, 2 * *ndjer * sizeof (sInt4));
2026 : for (i = 0; i < 8; i++) {
2027 : jer[i] = i * 100;
2028 : }
2029 : *kjer = 8;
2030 :
2031 : gdsTmpl = is3[13 - 1];
2032 : pdsTmpl = is4[8 - 1];
2033 : drsTmpl = is5[10 - 1];
2034 : if ((gdsTmpl != 0) && (gdsTmpl != 10) && (gdsTmpl != 20) &&
2035 : (gdsTmpl != 30) && (gdsTmpl != 90) && (gdsTmpl != 110) &&
2036 : (gdsTmpl != 120)) {
2037 : f_useMDL = 0;
2038 : }
2039 : if ((pdsTmpl != 0) && (pdsTmpl != 1) && (pdsTmpl != 2) &&
2040 : (pdsTmpl != 8) && (pdsTmpl != 9) && (pdsTmpl != 20) &&
2041 : (pdsTmpl != 30)) {
2042 : f_useMDL = 0;
2043 : }
2044 : if ((drsTmpl != 0) && (drsTmpl != 2) && (drsTmpl != 3)) {
2045 : f_useMDL = 0;
2046 : }
2047 : /*
2048 : printf ("Forcing it to not use MDL encoder. \n");
2049 : f_useMDL = 0;
2050 : */
2051 : if (f_useMDL) {
2052 : /* pk_cmplx.f ==> jmax(M), jmin(M), nov(M), lbit(M) */
2053 : jmax = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2054 : jmin = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2055 : nov = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2056 : lbit = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2057 : /* pk_grib2.f ==> a(NX,NY), ia(NX,NY) */
2058 : /* a = (float *) malloc ((*nx) * (*ny) * sizeof (float));*/
2059 : ia = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2060 : /* pack_gp.f ==> misslx(NDG = NXY = nx*ny) */
2061 : misslx = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2062 : /* Able to use jmax(nxy) for iwork(nxy) in int_map.f, and pk_missp.f */
2063 : /* Able to use jmax(nxy) for work(nxy) in flt_map.f */
2064 : /* reduce.f ==> newbox(ndg=nxy), newboxp(ndg=nxy) */
2065 : newbox = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2066 : newboxp = (sInt4 *) malloc ((*nx) * (*ny) * sizeof (sInt4));
2067 :
2068 : /* other known automatic arrays... itemp(nidat) pk_sect2.f */
2069 : /* other known automatic arrays... rtemp(nrdat) pk_sect2.f */
2070 :
2071 : /* a and ia are equivalenced here. */
2072 : #ifdef _FORTRAN
2073 : PK_G2MDL (kfildo, jmax, jmin, lbit, nov, misslx, (float *) ia, ia,
2074 : newbox, newboxp, ain, iain, nx, ny, idat, nidat, rdat, nrdat,
2075 : is0, ns0, is1, ns1, is3, ns3, is4, ns4, is5, ns5, is6, ns6,
2076 : is7, ns7, ib, ibitmap, ipack, nd5, missp, xmissp, misss,
2077 : xmisss, inew, minpk, iclean, l3264b, jer, ndjer, kjer);
2078 : #endif /* _FORTRAN */
2079 :
2080 : free (jmax);
2081 : free (jmin);
2082 : free (nov);
2083 : free (lbit);
2084 : free (misslx);
2085 : /* free (a);*/
2086 : free (ia);
2087 : free (newbox);
2088 : free (newboxp);
2089 : } else {
2090 :
2091 :
2092 : #ifdef PKNCEP
2093 : #ifdef LITTLE_ENDIAN
2094 : /* Can't make this dependent on inew, since they could have a sequence
2095 : * of get first message... do stuff, get second message, which
2096 : * unfortunately means they would have to get the first message again,
2097 : * causing 2 swaps, and breaking their second request for the first
2098 : * message (as well as their second message). */
2099 : /*
2100 : memswp (ipack, sizeof (sInt4), *nd5);
2101 : */
2102 : #endif
2103 : c_ipack = (unsigned char *) ipack;
2104 : pk_g2ncep (kfildo, ain, iain, nx, ny, idat, nidat, rdat, nrdat, is0,
2105 : ns0, is1, ns1, is3, ns3, is4, ns4, is5, ns5, is6, ns6, is7,
2106 : ns7, ib, ibitmap, c_ipack, nd5, missp, xmissp, misss, xmisss,
2107 : inew, minpk, iclean, l3264b, jer, ndjer, kjer);
2108 :
2109 : #ifdef LITTLE_ENDIAN
2110 : /* Swap back because we could be called again for the subgrid data. */
2111 : memswp (ipack, sizeof (sInt4), *nd5);
2112 : #endif
2113 : #else
2114 :
2115 : printf ("Unable to use MDL Pack library?\n");
2116 : printf ("gdsTmpl : %d , pdsTmpl %d : drsTmpl %d\n", gdsTmpl,
2117 : pdsTmpl, drsTmpl);
2118 : jer[0 + *ndjer] = 31415926;
2119 : *kjer = 1;
2120 : #endif
2121 : }
2122 :
2123 : #endif /* _FORTRAN */
2124 : }
|
