1 /**
2 * Copyright 1981-2016 ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence
5 * Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities
8 * granted to it by virtue of its status as an intergovernmental organisation
9 * nor does it submit to any jurisdiction.
10 */
11
12 #include "gdecode.h"
13 #include "sencode.h"
14
15 #define DEBUGOFF 1
16 #define DEBUG1 (debugSet > DEBUGOFF )
17 #define DEBUG2 (debugSet > (DEBUGOFF + 1) )
18 static char * debugLevel;
19 static int debugSet = 0;
20
21 #define ABS(a) ((a)<0?-(a):(a))
22 #define SMALL (1E-12)
23 #define PRACTICALLYZERO(a) (ABS((a))<(SMALL))
24 #define QUITESMALL (0.05)
25 #define ALMOSTZERO(a) (ABS((a))<(QUITESMALL))
26 #define TRUE 1
27 #define FALSE 0
28 #define CHECK 1
29 #define NOCHECK 0
30
31 #define NORTH 0
32 #define SOUTH 1
33 #define WEST 2
34 #define EAST 3
35 #define W_E_INC 4
36 #define N_S_INC 5
37 #define W_E_PTS 6
38 #define N_S_PTS 7
39
40 fortint initialiseNewGrib(gribProduct**);
41 fortint copyExistingGrib(gribProduct**,gribProduct**);
42 fortint getIntegerValue(gribProduct**,unsigned char*);
43 void adjustGridAreaDefinition(gribProduct*,int);
44
SENCODE(gribProduct ** newGrib,gribProduct ** oldGrib)45 fortint SENCODE(gribProduct ** newGrib, gribProduct ** oldGrib) {
46 fortint status;
47 /*
48 // See if DEBUG switched on.
49 */
50 if( ! debugSet ) {
51 debugLevel = getenv("GDECODE_DEBUG");
52 if( debugLevel == NULL )
53 debugSet = DEBUGOFF; /* off */
54 else {
55 int loop;
56 for( loop = 0; loop < strlen(debugLevel) ; loop++ ) {
57 if( ! isdigit(debugLevel[loop]) ) {
58 printf("Invalid number string in GDECODE_DEBUG: %s\n", debugLevel);
59 printf("GDECODE_DEBUG must comprise only digits [0-9].\n");
60 debugSet = DEBUGOFF;
61 }
62 }
63 debugSet = DEBUGOFF + atol( debugLevel );
64 }
65 if( DEBUG1 ) printf("SENCODE: GDECODE_DEBUG switched on, level = %s\n",
66 debugLevel);
67 }
68
69 if( *oldGrib == NULL ) {
70 status = initialiseNewGrib(newGrib);
71 if( DEBUG1 ) {
72 if( status )
73 printf("SENCODE: initialiseNewGrib failed.\n");
74 else
75 printf("SENCODE: initialiseNewGrib ran OK\n");
76 }
77 }
78 else {
79 status = copyExistingGrib(newGrib,oldGrib);
80 if( DEBUG1 ) {
81 if( status )
82 printf("SENCODE: copyExistingGrib failed.\n");
83 else
84 printf("SENCODE: copyExistingGrib ran OK\n");
85 }
86 }
87
88 return status;
89
90 }
91
initialiseNewGrib(gribProduct ** newGrib)92 fortint initialiseNewGrib(gribProduct** newGrib) {
93
94 if( DEBUG1 ) printf("initialiseNewGrib\n");
95
96 *newGrib = (gribProduct *) allocateMemory(sizeof(gribProduct));
97 (*newGrib)->g0 = (gribSection0 *) allocateMemory(sizeof(gribSection0));
98 (*newGrib)->g1 = (gribSection1 *) allocateMemory(sizeof(gribSection1));
99 (*newGrib)->g2 = (gribSection2 *) allocateMemory(sizeof(gribSection2));
100 (*newGrib)->g3 = NULL;
101 (*newGrib)->g4 = (gribSection4 *) allocateMemory(20);
102 (*newGrib)->g5 = (gribSection5 *) allocateMemory(sizeof(gribSection5));
103 memcpy((*((*newGrib)->g5)).end7777,"7777",4);
104
105 (*newGrib)->currentPointIndex = 0;
106 (*newGrib)->numberOfValues = 0;
107 (*newGrib)->value = NULL;
108 (*newGrib)->bitStart = NULL;
109 (*newGrib)->bitsPerValue = 0;
110 (*newGrib)->bitmapped = 0;
111 (*newGrib)->nextValueFirstBit = 0;
112 (*newGrib)->nextBit = 0;
113 (*newGrib)->scale = (fortdouble) 0.0;
114 (*newGrib)->minimum = (fortdouble) 0.0;
115 (*newGrib)->missingValue = (fortdouble) 0.0;
116 (*newGrib)->latitudeOffsets = NULL;
117 (*newGrib)->expandedValues = NULL;
118 (*newGrib)->currentPoint.latitude = NULL;
119 (*newGrib)->currentPoint.longitudeIncrement = NULL;
120 (*newGrib)->currentPoint.gridPointValue = NULL;
121 (*newGrib)->northSet =
122 (*newGrib)->southSet =
123 (*newGrib)->westSet =
124 (*newGrib)->eastSet =
125 (*newGrib)->northSouthIncrementSet =
126 (*newGrib)->westEastIncrementSet =
127 (*newGrib)->northSouthNumberOfPointsSet =
128 (*newGrib)->westEastNumberOfPointsSet = TRUE;
129
130 return (fortint) 0;
131 }
132
copyExistingGrib(gribProduct ** newGrib,gribProduct ** oldGrib)133 fortint copyExistingGrib(gribProduct** newGrib, gribProduct** oldGrib) {
134 fortint lengthOfSectionToCopy;
135 gribProduct * oldG = *oldGrib;
136 gribProduct * newG;
137
138 if( DEBUG1 ) printf("copyExistingGrib\n");
139
140 *newGrib = (gribProduct *) allocateMemory(sizeof(gribProduct));
141 newG = *newGrib;
142 /*
143 // Create section 0
144 */
145 newG->g0 = (gribSection0 *) allocateMemory(sizeof(gribSection0));
146 memcpy((newG->g0),(oldG->g0),sizeof(gribSection0));
147 /*
148 // Copy section 1
149 */
150 lengthOfSectionToCopy = g1_length(oldG);
151 newG->g1 = (gribSection1 *) allocateMemory(lengthOfSectionToCopy);
152 memcpy((newG->g1),(oldG->g1),lengthOfSectionToCopy);
153 /*
154 // Copy section 2
155 */
156 lengthOfSectionToCopy = g2_length(oldG);
157 newG->g2 = (gribSection2 *) allocateMemory(lengthOfSectionToCopy);
158 memcpy((newG->g2),(oldG->g2),lengthOfSectionToCopy);
159
160 /*
161 if( anyLatLonGrid(oldG) ) {
162 */
163 if( generalLatLonGrid(oldG) ) {
164 newG->northSet =
165 newG->southSet =
166 newG->westSet =
167 newG->eastSet =
168 newG->northSouthIncrementSet =
169 newG->westEastIncrementSet =
170 newG->northSouthNumberOfPointsSet =
171 newG->westEastNumberOfPointsSet = TRUE;
172 }
173 /*
174 // Section 3 may or may not have to be created.
175 */
176 newG->g3 = NULL;
177 /*
178 // Copy standard part of section 4
179 */
180 lengthOfSectionToCopy = 20;
181 newG->g4 = (gribSection4 *) allocateMemory(lengthOfSectionToCopy);
182 memcpy((newG->g4),(oldG->g4),lengthOfSectionToCopy);
183 /*
184 // Create section 5.
185 */
186 newG->g5 = (gribSection5 *) allocateMemory(sizeof(gribSection5));
187 memcpy((*(newG->g5)).end7777,"7777",4);
188
189
190 newG->currentPointIndex = 0;
191 newG->numberOfValues = getIntegerValue(oldGrib,(unsigned char*)"numberOfFieldValues");
192 newG->value = NULL;
193 newG->bitStart = NULL;
194 newG->bitsPerValue = oldG->bitsPerValue;
195 newG->bitmapped = oldG->bitmapped;
196 newG->nextValueFirstBit = 0;
197 newG->nextBit = 0;
198 newG->scale = (fortdouble) oldG->scale;
199 newG->minimum = (fortdouble) oldG->minimum;
200 newG->missingValue = (fortdouble) oldG->missingValue;
201 newG->latitudeOffsets = NULL;
202 newG->currentPoint.latitude = NULL;
203 newG->currentPoint.longitudeIncrement = NULL;
204 newG->expandedValues = NULL;
205 newG->currentPoint.gridPointValue = NULL;
206
207 return (fortint) 0;
208 }
209
SENPACK(gribProduct ** newGrib,unsigned char * buffer,fortint * bufferLength)210 fortint SENPACK(
211 gribProduct** newGrib, unsigned char * buffer, fortint * bufferLength)
212 {
213 fortint totalLength;
214 gribProduct* g = *newGrib;
215 fortint lengthOfSectionToMove, next;
216
217 totalLength = 8 + g1_length(g) + g2_length(g) + g4_length(g) + 4;
218 if( primaryBitmapPresent(g) ) totalLength += g3_length(g);
219
220 if( DEBUG1 ) printf("GRIB totalLength = %d bytes\n",totalLength);
221
222 if( totalLength > *bufferLength ) {
223 if( DEBUG1 ) printf("GRIB length (%d) greater than buffer length (%d)\n",
224 totalLength, *bufferLength );
225 return (fortint) -1;
226 }
227 else {
228 memcpy(buffer,(g->g0),sizeof(gribSection0));
229 MOVE3BYTES((buffer+4),&totalLength);
230 next = sizeof(gribSection0);
231
232 lengthOfSectionToMove = g1_length(g);
233 memcpy((buffer+next),(g->g1),lengthOfSectionToMove);
234 next += g1_length(g);
235
236 lengthOfSectionToMove = g2_length(g);
237 memcpy((buffer+next),(g->g2),lengthOfSectionToMove);
238 next += g2_length(g);
239
240 if( primaryBitmapPresent(g) > 0 ) {
241 lengthOfSectionToMove = g3_length(g);
242 memcpy((buffer+next),(g->g3),lengthOfSectionToMove);
243 next += g3_length(g);
244 }
245
246 lengthOfSectionToMove = g4_length(g);
247 memcpy((buffer+next),(g->g4),lengthOfSectionToMove);
248 next += g4_length(g);
249
250 lengthOfSectionToMove = 4;
251 memcpy((buffer+next),"7777",lengthOfSectionToMove);
252
253 return (fortint) totalLength;
254 }
255 }
256
ISTIME(gribProduct ** grib,fortint * HHMM)257 fortint ISTIME(gribProduct ** grib, fortint * HHMM ) {
258 gribProduct * g = *grib;
259 fortint hhmm = (fortint) (*HHMM);
260 fortint hh = hhmm/100;
261 fortint mm = MOD(hhmm,100);
262
263 if( DEBUG1 ) printf("ISTIME: HHMM = %d\n", *HHMM);
264
265 MOVE1BYTE(((g->g1)->hour),&hh);
266 MOVE1BYTE(((g)->g1->minute),&mm);
267
268 return 0;
269 }
270
RSTIME(gribProduct ** grib,double * HHMM)271 fortint RSTIME(gribProduct ** grib, double * HHMM ) {
272 fortint hhmm = (fortint) (*HHMM);
273
274 if( DEBUG1 ) printf("RSTIME: HHMM = %f\n", *HHMM);
275
276 return ISTIME(grib,&hhmm);
277 }
278
ISDATE(gribProduct ** grib,fortint * YYYYMMDD)279 fortint ISDATE(gribProduct ** grib, fortint * YYYYMMDD ) {
280 gribProduct * g = *grib;
281 fortint yyyymmdd = *YYYYMMDD;
282 fortint yyyy, mmdd, year, month, day, century;
283
284 if( DEBUG1 ) printf("ISDATE: YYYYMMDD = %d\n", yyyymmdd);
285
286 yyyy = yyyymmdd/10000;
287 mmdd = MOD(yyyymmdd,10000);
288 year = MOD(yyyy,100);
289 month = mmdd/100;
290 day = MOD(mmdd,100);
291 century = 1 + (yyyy/100);
292 if( year == 0 ) century -= 1;
293
294 MOVE1BYTE(((g->g1)->year),&year);
295 MOVE1BYTE(((g)->g1->month),&month);
296 MOVE1BYTE(((g)->g1->day),&day);
297 MOVE1BYTE(((g)->g1->century),¢ury);
298
299 return 0;
300 }
301
RSDATE(gribProduct ** grib,double * YYYYMMDD)302 fortint RSDATE(gribProduct ** grib, double * YYYYMMDD ) {
303 fortint yyyymmdd = (fortint) *YYYYMMDD;
304
305 if( DEBUG1 ) printf("RSDATE: YYYYMMDD = %f\n", *YYYYMMDD);
306
307 return ISDATE(grib,&yyyymmdd);
308 }
309
ISTABLE(gribProduct ** grib,fortint * table)310 fortint ISTABLE(gribProduct ** grib, fortint * table) {
311 gribProduct * g = *grib;
312
313 if( DEBUG1 ) printf("ISTABLE: table = %d\n", *table);
314
315 MOVE1BYTE(((g->g1)->tableVersionNumber),table);
316
317 return 0;
318 }
319
RSTABLE(gribProduct ** grib,fortdouble * table)320 fortint RSTABLE(gribProduct ** grib, fortdouble * table) {
321 gribProduct * g = *grib;
322 fortint Table = (fortint) *table;
323
324 if( DEBUG1 ) printf("RSTABLE: table = %f\n", *table);
325
326 return ISTABLE(grib,&Table);
327 }
328
ISCENTR(gribProduct ** grib,fortint * centre)329 fortint ISCENTR(gribProduct ** grib, fortint * centre) {
330 gribProduct * g = *grib;
331
332 if( DEBUG1 ) printf("ISCENTR: centre = %d\n", *centre);
333
334 MOVE1BYTE(((g->g1)->originatingCentre),centre);
335
336 return 0;
337 }
338
RSCENTR(gribProduct ** grib,fortdouble * centre)339 fortint RSCENTR(gribProduct ** grib, fortdouble * centre) {
340 gribProduct * g = *grib;
341 fortint Centre = (fortint) *centre;
342
343 if( DEBUG1 ) printf("RSCENTR: centre = %f\n", *centre);
344
345 return ISCENTR(grib,&Centre);
346 }
347
ISPARAM(gribProduct ** grib,fortint * parameter)348 fortint ISPARAM(gribProduct ** grib, fortint * parameter) {
349 gribProduct * g = *grib;
350
351 if( DEBUG1 ) printf("ISPARAM: parameter = %d\n", *parameter);
352
353 MOVE1BYTE(((g->g1)->parameter),parameter);
354
355 return 0;
356 }
357
RSPARAM(gribProduct ** grib,fortdouble * parameter)358 fortint RSPARAM(gribProduct ** grib, fortdouble * parameter) {
359 gribProduct * g = *grib;
360 fortint Parameter = (fortint) *parameter;
361
362 if( DEBUG1 ) printf("RSPARAM: parameter = %f\n", *parameter);
363
364 return ISPARAM(grib,&Parameter);
365 }
366
ISLEVTY(gribProduct ** grib,fortint * typeOfLevel)367 fortint ISLEVTY(gribProduct ** grib, fortint * typeOfLevel) {
368 gribProduct * g = *grib;
369 fortint TypeOfLevel = *typeOfLevel;
370 fortint zero =0;
371
372 if( DEBUG1 ) printf("ISLEVTY: typeOfLevel = %d\n", TypeOfLevel);
373
374 MOVE1BYTE(((g->g1)->typeOfLevel),typeOfLevel);
375
376 if( (TypeOfLevel < 20) ||
377 ((TypeOfLevel > 20) && (TypeOfLevel < 100)) ||
378 (TypeOfLevel == 102) ||
379 (TypeOfLevel == 118) ||
380 ((TypeOfLevel > 121) && (TypeOfLevel < 125)) ||
381 (TypeOfLevel == 126) ||
382 ((TypeOfLevel > 128) && (TypeOfLevel < 141)) ||
383 ((TypeOfLevel > 141) && (TypeOfLevel < 160)) ||
384 ((TypeOfLevel > 160) && (TypeOfLevel < 200)) ||
385 ((TypeOfLevel > 201) && (TypeOfLevel < 210)) ||
386 (TypeOfLevel > 210) ) {
387 MOVE2BYTES(((g->g1)->level1),&zero);
388 }
389
390 return 0;
391 }
392
RSLEVTY(gribProduct ** grib,fortdouble * typeOfLevel)393 fortint RSLEVTY(gribProduct ** grib, fortdouble * typeOfLevel) {
394 gribProduct * g = *grib;
395 fortint TypeOfLevel = (fortint) *typeOfLevel;
396
397 if( DEBUG1 ) printf("RSLEVTY: typeOfLevel = %f\n", *typeOfLevel);
398
399 return ISLEVTY(grib,&TypeOfLevel);
400 }
401
ISLEVEL(gribProduct ** grib,fortint * level)402 fortint ISLEVEL(gribProduct ** grib, fortint * level) {
403 gribProduct * g = *grib;
404 fortint zero = 0, bottomLayer, topLayer;
405
406 if( DEBUG1 ) printf("ISLEVEL: level = %d\n", *level);
407
408 switch( (int) g1_typeOfLevel(g) ) {
409
410 case 20:
411 case 100:
412 case 103:
413 case 105:
414 case 107:
415 case 109:
416 case 111:
417 case 113:
418 case 115:
419 case 117:
420 case 119:
421 case 125:
422 case 127:
423 case 160:
424 case 210:
425 if( DEBUG1 ) printf("ISLEVEL: two-byte level value\n");
426 MOVE2BYTES(((g->g1)->level1),level);
427 break;
428
429 case 101:
430 case 104:
431 case 106:
432 case 108:
433 case 110:
434 case 112:
435 case 114:
436 case 116:
437 case 120:
438 case 121:
439 case 128:
440 case 141:
441 if( DEBUG1 ) printf("ISLEVEL: top and bottom level values\n");
442 topLayer = (*level)/1000;
443 MOVE1BYTE(((g->g1)->level1),&topLayer);
444 bottomLayer = MOD((*level),1000);
445 MOVE1BYTE(((g->g1)->level2),&bottomLayer);
446 break;
447
448 default:
449 if( DEBUG1 ) printf("ISLEVEL: level value set to zero\n");
450 MOVE2BYTES(((g->g1)->level1),&zero);
451 break;
452 }
453
454
455 return 0;
456 }
457
RSLEVEL(gribProduct ** grib,fortdouble * level)458 fortint RSLEVEL(gribProduct ** grib, fortdouble * level) {
459 gribProduct * g = *grib;
460 fortint Level = (fortint) *level;
461
462 if( DEBUG1 ) printf("RSLEVEL: level = %f\n", *level);
463
464 return ISLEVEL(grib,&Level);
465 }
466
ISTUNIT(gribProduct ** grib,fortint * timeUnit)467 fortint ISTUNIT(gribProduct ** grib, fortint * timeUnit) {
468 gribProduct * g = *grib;
469
470 if( DEBUG1 ) printf("ISTUNIT: timeUnit = %d\n", *timeUnit);
471
472 switch( (int) *timeUnit ) {
473
474 case 0:
475 case 1:
476 case 2:
477 case 3:
478 case 4:
479 case 5:
480 case 6:
481 case 7:
482 case 10:
483 case 11:
484 case 12:
485 case 254:
486 MOVE1BYTE(((g->g1)->unitOfTimeRange),timeUnit);
487 break;
488 return 0;
489
490 default:
491 if( DEBUG1 ) printf("ISTUNIT: invalid time unit\n");
492 }
493 return -1;
494 }
495
RSTUNIT(gribProduct ** grib,fortdouble * timeUnit)496 fortint RSTUNIT(gribProduct ** grib, fortdouble * timeUnit) {
497 gribProduct * g = *grib;
498 fortint TimeUnit = (fortint) *timeUnit;
499
500 if( DEBUG1 ) printf("RSTUNIT: timeUnit = %f\n", *timeUnit);
501
502 return ISTUNIT(grib,&TimeUnit);
503 }
504
ISTRIND(gribProduct ** grib,fortint * timeRangeIndicator)505 fortint ISTRIND(gribProduct ** grib, fortint * timeRangeIndicator) {
506 gribProduct * g = *grib;
507
508 if( DEBUG1 ) printf("ISTRIND: timeRangeIndicator = %d\n",*timeRangeIndicator);
509 MOVE1BYTE(((g->g1)->timeRangeIndicator),timeRangeIndicator);
510
511 return 0;
512 }
513
RSTRIND(gribProduct ** grib,fortdouble * timeRangeIndicator)514 fortint RSTRIND(gribProduct ** grib, fortdouble * timeRangeIndicator) {
515 gribProduct * g = *grib;
516 fortint TimeRangeIndicator = (fortint) *timeRangeIndicator;
517
518 if( DEBUG1 ) printf("RSTRIND: timeRangeIndicator = %f\n",*timeRangeIndicator);
519
520 return ISTRIND(grib,&TimeRangeIndicator);
521 }
522
ISSTEP(gribProduct ** grib,fortint * timeStep)523 fortint ISSTEP(gribProduct ** grib, fortint * timeStep) {
524 gribProduct * g = *grib;
525 fortint zero = 0;
526
527 if( DEBUG1 ) printf("ISSTEP: timeStep = %d\n",*timeStep);
528
529 switch( (int) g1_timerange(g) ) {
530
531 case 0:
532 MOVE1BYTE(((g->g1)->P1),timeStep);
533 MOVE1BYTE(((g->g1)->P2),&zero);
534 return 0;
535
536 case 1:
537 MOVE1BYTE(((g->g1)->P1),&zero);
538 MOVE1BYTE(((g->g1)->P2),&zero);
539 return 0;
540
541 case 2:
542 case 3:
543 case 4:
544 case 5:
545 case 113:
546 case 114:
547 case 115:
548 case 116:
549 case 117:
550 case 118:
551 case 119:
552 printf("ISSTEP: time range indicator %d requires separate values for P1 and P2\n", g1_timerange(g));
553 return -1;
554
555 case 10:
556 MOVE2BYTES(((g->g1)->P1),timeStep);
557 return 0;
558
559 case 123:
560 case 124:
561 MOVE1BYTE(((g->g1)->P1),&zero);
562 MOVE1BYTE(((g->g1)->P2),timeStep);
563 return 0;
564
565 default:
566 printf("ISSTEP: unable to set step for reserved time range indicator %d\n", g1_timerange(g));
567 return -1;
568 }
569
570 }
571
RSSTEP(gribProduct ** grib,fortdouble * timeStep)572 fortint RSSTEP(gribProduct ** grib, fortdouble * timeStep) {
573 gribProduct * g = *grib;
574 fortint TimeStep = (fortint) *timeStep;
575
576 if( DEBUG1 ) printf("RSSTEP: timeStep = %f\n",*timeStep);
577
578 return ISSTEP(grib,&TimeStep);
579 }
580
ISSTEP1(gribProduct ** grib,fortint * timeStepP1)581 fortint ISSTEP1(gribProduct ** grib, fortint * timeStepP1) {
582 gribProduct * g = *grib;
583 fortint zero = 0;
584
585 if( DEBUG1 ) printf("ISSTEP1: timeStepP1 = %d\n",*timeStepP1);
586
587 switch( (int) g1_timerange(g) ) {
588
589 case 0:
590 case 2:
591 case 3:
592 case 4:
593 case 5:
594 case 113:
595 case 114:
596 case 115:
597 case 116:
598 case 117:
599 case 118:
600 case 119:
601 MOVE1BYTE(((g->g1)->P1),timeStepP1);
602 return 0;
603
604 case 1:
605 case 123:
606 case 124:
607 MOVE1BYTE(((g->g1)->P1),&zero);
608 return 0;
609
610 case 10:
611 MOVE2BYTES(((g->g1)->P1),timeStepP1);
612 return 0;
613
614 default:
615 printf("ISSTEP1: unable to set P1 for reserved time range indicator %d\n", g1_timerange(g));
616 return -1;
617 }
618
619 }
620
RSSTEP1(gribProduct ** grib,fortdouble * timeStepP1)621 fortint RSSTEP1(gribProduct ** grib, fortdouble * timeStepP1) {
622 gribProduct * g = *grib;
623 fortint TimeStepP1 = (fortint) *timeStepP1;
624
625 if( DEBUG1 ) printf("RSSTEP1: timeStepP1 = %f\n",*timeStepP1);
626
627 return ISSTEP1(grib,&TimeStepP1);
628 }
629
ISSTEP2(gribProduct ** grib,fortint * timeStepP2)630 fortint ISSTEP2(gribProduct ** grib, fortint * timeStepP2) {
631 gribProduct * g = *grib;
632 fortint zero = 0;
633
634 if( DEBUG1 ) printf("ISSTEP2: timeStepP2 = %d\n",*timeStepP2);
635
636 switch( (int) g1_timerange(g) ) {
637
638 case 0:
639 case 1:
640 MOVE1BYTE(((g->g1)->P2),&zero);
641 return 0;
642
643 case 2:
644 case 3:
645 case 4:
646 case 5:
647 case 113:
648 case 114:
649 case 115:
650 case 116:
651 case 117:
652 case 118:
653 case 119:
654 case 123:
655 case 124:
656 MOVE1BYTE(((g->g1)->P2),timeStepP2);
657 return 0;
658
659 case 10:
660 printf("ISSTEP2: unable to set P2 for time range indicator %d\n", g1_timerange(g));
661 return -1;
662
663 default:
664 printf("ISSTEP2: unable to set P2 for reserved time range indicator %d\n", g1_timerange(g));
665 return -1;
666 }
667
668 }
669
RSSTEP2(gribProduct ** grib,fortdouble * timeStepP2)670 fortint RSSTEP2(gribProduct ** grib, fortdouble * timeStepP2) {
671 gribProduct * g = *grib;
672 fortint TimeStepP2 = (fortint) *timeStepP2;
673
674 if( DEBUG1 ) printf("RSSTEP2: timeStepP2 = %f\n",*timeStepP2);
675
676 return ISSTEP2(grib,&TimeStepP2);
677 }
678
ISNUMAV(gribProduct ** grib,fortint * numberInAverage)679 fortint ISNUMAV(gribProduct ** grib, fortint * numberInAverage) {
680 gribProduct * g = *grib;
681
682 if( DEBUG1 ) printf("ISNUMAV: numberInAverage = %d\n",*numberInAverage);
683
684 MOVE2BYTES(((g->g1)->numberInAverage),numberInAverage);
685 return 0;
686 }
687
RSNUMAV(gribProduct ** grib,fortdouble * numberInAverage)688 fortint RSNUMAV(gribProduct ** grib, fortdouble * numberInAverage) {
689 gribProduct * g = *grib;
690 fortint NumberInAverage = (fortint) *numberInAverage;
691
692 if( DEBUG1 ) printf("RSNUMAV: numberInAverage = %f\n",*numberInAverage);
693
694 return ISNUMAV(grib,&NumberInAverage);
695
696 }
697
ISNUMMS(gribProduct ** grib,fortint * numberMissing)698 fortint ISNUMMS(gribProduct ** grib, fortint * numberMissing) {
699 gribProduct * g = *grib;
700
701 if( DEBUG1 ) printf("ISNUMMS: numberMissing = %d\n",*numberMissing);
702
703 MOVE1BYTE(((g->g1)->numberMissing),numberMissing);
704 return 0;
705 }
706
RSNUMMS(gribProduct ** grib,fortdouble * numberMissing)707 fortint RSNUMMS(gribProduct ** grib, fortdouble * numberMissing) {
708 gribProduct * g = *grib;
709 fortint NumberMissing = (fortint) *numberMissing;
710
711 if( DEBUG1 ) printf("RSNUMMS: numberMissing = %f\n",*numberMissing);
712
713 return ISNUMMS(grib,&NumberMissing);
714 }
715
ISSUBID(gribProduct ** grib,fortint * subCentreId)716 fortint ISSUBID(gribProduct ** grib, fortint * subCentreId) {
717 gribProduct * g = *grib;
718
719 if( DEBUG1 ) printf("ISSUBID: subCentreId = %d\n",*subCentreId);
720
721 MOVE1BYTE(((g->g1)->subCentreId),subCentreId);
722 return 0;
723 }
724
RSSUBID(gribProduct ** grib,fortdouble * subCentreId)725 fortint RSSUBID(gribProduct ** grib, fortdouble * subCentreId) {
726 gribProduct * g = *grib;
727 fortint SubCentreId = (fortint) *subCentreId;
728
729 if( DEBUG1 ) printf("RSSUBID: subCentreId = %f\n",*subCentreId);
730
731 return ISSUBID(grib,&SubCentreId);
732 }
733
ISUDECF(gribProduct ** grib,fortint * decimalScale)734 fortint ISUDECF(gribProduct ** grib, fortint * decimalScale) {
735 gribProduct * g = *grib;
736 fortint scale = *decimalScale;
737
738 if( DEBUG1 ) printf("ISUDECF: decimalScale = %d\n",*decimalScale);
739
740 if( scale < 0 ) scale = (-scale) | 0x8000;
741 MOVE2BYTES(((g->g1)->unitsDecimalScaleFactor),&scale);
742 return 0;
743 }
744
RSUDECF(gribProduct ** grib,fortdouble * decimalScale)745 fortint RSUDECF(gribProduct ** grib, fortdouble * decimalScale) {
746 gribProduct * g = *grib;
747 fortint DecimalScale = (fortint) *decimalScale;
748
749 if( DEBUG1 ) printf("RSUDECF: decimalScale = %f\n",*decimalScale);
750
751 return ISUDECF(grib,&DecimalScale);
752 }
753
ISTYPE(gribProduct ** grib,fortint * ecmwfType)754 fortint ISTYPE(gribProduct ** grib, fortint * ecmwfType) {
755 gribProduct * g = *grib;
756
757 if( DEBUG1 ) printf("ISTYPE: ecmwfType = %d\n",*ecmwfType);
758
759 if( ecmwfLocalDefinitionPresent(g) ) {
760 MOVE1BYTE(((g->g1)->local.mars.type),ecmwfType);
761 return 0;
762 }
763
764 if( centreUsingECMWFLocalDefinition((*grib)) ) {
765 MOVE1BYTE(((g->g1)->local.mars.type),ecmwfType);
766 return 0;
767 }
768 else {
769 if( DEBUG1 ) printf("ISTYPE: no ECMWF local definition present\n");
770 return (fortint) -1;
771 }
772
773 }
774
RSTYPE(gribProduct ** grib,fortdouble * ecmwfType)775 fortint RSTYPE(gribProduct ** grib, fortdouble * ecmwfType) {
776 gribProduct * g = *grib;
777 fortint EcmwfType = (fortint) *ecmwfType;
778
779 if( DEBUG1 ) printf("RSTYPE: ecmwfType = %f\n",*ecmwfType);
780
781 return ISTYPE(grib,&EcmwfType);
782 }
783
ISCLASS(gribProduct ** grib,fortint * ecmwfClass)784 fortint ISCLASS(gribProduct ** grib, fortint * ecmwfClass) {
785 gribProduct * g = *grib;
786
787 if( DEBUG1 ) printf("ISCLASS: ecmwfClass = %d\n",*ecmwfClass);
788
789 if( ecmwfLocalDefinitionPresent(g) ) {
790 MOVE1BYTE(((g->g1)->local.mars.ecmwfClass),ecmwfClass);
791 return 0;
792 }
793
794 if( centreUsingECMWFLocalDefinition((*grib)) ) {
795 MOVE1BYTE(((g->g1)->local.mars.ecmwfClass),ecmwfClass);
796 return 0;
797 }
798 else {
799 if( DEBUG1 ) printf("ISCLASS: no ECMWF local definition present\n");
800 return (fortint) -1;
801 }
802
803 }
804
RSCLASS(gribProduct ** grib,fortdouble * ecmwfClass)805 fortint RSCLASS(gribProduct ** grib, fortdouble * ecmwfClass) {
806 gribProduct * g = *grib;
807 fortint EcmwfClass = (fortint) *ecmwfClass;
808
809 if( DEBUG1 ) printf("RSCLASS: ecmwfClass = %f\n",*ecmwfClass);
810
811 return ISCLASS(grib,&EcmwfClass);
812 }
813
ISSTREM(gribProduct ** grib,fortint * ecmwfStream)814 fortint ISSTREM(gribProduct ** grib, fortint * ecmwfStream) {
815 gribProduct * g = *grib;
816
817 if( DEBUG1 ) printf("ISSTREM: ecmwfStream = %d\n",*ecmwfStream);
818
819 if( ecmwfLocalDefinitionPresent(g) ) {
820 MOVE2BYTES(((g->g1)->local.mars.stream),ecmwfStream);
821 return 0;
822 }
823
824 if( centreUsingECMWFLocalDefinition((*grib)) ) {
825 MOVE2BYTES(((g->g1)->local.mars.stream),ecmwfStream);
826 return 0;
827 }
828 else {
829 if( DEBUG1 ) printf("ISSTREM: no ECMWF local definition present\n");
830 return (fortint) -1;
831 }
832
833 }
834
RSSTREM(gribProduct ** grib,fortdouble * ecmwfStream)835 fortint RSSTREM(gribProduct ** grib, fortdouble * ecmwfStream) {
836 gribProduct * g = *grib;
837 fortint EcmwfStream = (fortint) *ecmwfStream;
838
839 if( DEBUG1 ) printf("RSSTREM: ecmwfStream = %f\n",*ecmwfStream);
840
841 return ISSTREM(grib,&EcmwfStream);
842 }
843
ISEXPVR(gribProduct ** grib,fortint * ecmwfExpver)844 fortint ISEXPVR(gribProduct ** grib, fortint * ecmwfExpver) {
845 gribProduct * g = *grib;
846
847 if( DEBUG1 ) printf("ISEXPVR: ecmwfExpver = %d\n",*ecmwfExpver);
848
849 if( ecmwfLocalDefinitionPresent(g) ) {
850 MOVE4BYTES(((g->g1)->local.mars.experimentVersionNumber),ecmwfExpver);
851 return 0;
852 }
853
854 if( centreUsingECMWFLocalDefinition((*grib)) ) {
855 MOVE4BYTES(((g->g1)->local.mars.experimentVersionNumber),ecmwfExpver);
856 return 0;
857 }
858 else {
859 if( DEBUG1 ) printf("ISEXPVR: no ECMWF local definition present\n");
860 return (fortint) -1;
861 }
862
863 }
864
RSEXPVR(gribProduct ** grib,fortdouble * ecmwfExpver)865 fortint RSEXPVR(gribProduct ** grib, fortdouble * ecmwfExpver) {
866 gribProduct * g = *grib;
867 fortint EcmwfExpver = (fortint) *ecmwfExpver;
868 gribSection1 * newG1;
869
870 if( DEBUG1 ) printf("RSEXPVR: ecmwfExpver = %f\n",*ecmwfExpver);
871
872 return ISEXPVR(grib,&EcmwfExpver);
873 }
874
ISDEFIN(gribProduct ** grib,fortint * definitionNumber)875 fortint ISDEFIN(gribProduct ** grib, fortint * definitionNumber) {
876 gribProduct * g = *grib;
877 localDefinition * local;
878 fortint length, oldCopyLength, loop;
879 gribSection1 * newG1, * oldG1;
880
881 if( DEBUG1 ) printf("ISDEFIN: definitionNumber = %d\n",*definitionNumber);
882
883 if( DEBUG1 ) {
884 if(ecmwfLocalDefinitionPresent(g) || centreUsingECMWFLocalDefinition(g)) {
885 if( *definitionNumber == g1_definition(g) ) {
886 printf("ISDEFIN: definitionNumber already has the given value\n");
887 }
888 else {
889 printf("ISDEFIN: changing the ECMWF local definition from %d to %d\n",
890 g1_definition(g), *definitionNumber);
891 }
892 }
893 else
894 printf("ISDEFIN: creating a new ECMWF local definition\n");
895 }
896
897 length = (40 + sizeof(marsHeader));
898
899 switch( (int) *definitionNumber ) {
900
901 case 1:
902 length += sizeof(ECMWFdefinition1);
903 break;
904
905 case 2:
906 length += sizeof(ECMWFdefinition2);
907 break;
908
909 case 3:
910 length += sizeof(ECMWFdefinition3);
911 break;
912
913 case 5:
914 length += sizeof(ECMWFdefinition5);
915 break;
916
917 case 6:
918 length += sizeof(ECMWFdefinition6);
919 break;
920
921 case 7:
922 length += sizeof(ECMWFdefinition7);
923 break;
924
925 case 8:
926 length += sizeof(ECMWFdefinition8);
927 break;
928
929 case 9:
930 length += sizeof(ECMWFdefinition9);
931 break;
932
933 case 10:
934 length += sizeof(ECMWFdefinition10);
935 break;
936
937 case 11:
938 length += sizeof(ECMWFdefinition11);
939 break;
940
941 case 14:
942 length = 1080;
943 break;
944
945 case 15:
946 length = 60;
947 break;
948
949 case 16:
950 length = 80;
951 break;
952
953 case 18:
954 case 19:
955 length = 120;
956 break;
957
958 case 20:
959 length += sizeof(ECMWFdefinition20);
960 break;
961
962 case 50:
963 length = 300;
964 break;
965
966 default:
967 printf("ISDEFIN: ECMWF local definition from %d not yet handled\n",
968 *definitionNumber);
969 return -1;
970 }
971
972 newG1 = (gribSection1 *) allocateMemory(length);
973 {
974 unsigned char * p = (unsigned char * ) newG1;
975 for( loop = 40; loop < length; loop++ ) p[loop] = '\0';
976 }
977
978 if( ecmwfLocalDefinitionPresent(g) || centreUsingECMWFLocalDefinition(g) )
979 oldCopyLength = (40 + sizeof(marsHeader));
980 else
981 oldCopyLength = 40;
982
983 oldG1 = (gribSection1 *) (g->g1);
984 memcpy(newG1,oldG1,oldCopyLength);
985
986 /*freeMemory(oldG1);*/
987 freeMemory(oldG1);
988
989 (g->g1) = newG1;
990 MOVE1BYTE(((g->g1)->local.mars.definition),definitionNumber);
991 MOVE3BYTES(((g->g1)->sectionLength),&length);
992 return 0;
993
994 }
995
RSDEFIN(gribProduct ** grib,fortdouble * definitionNumber)996 fortint RSDEFIN(gribProduct ** grib, fortdouble * definitionNumber) {
997 gribProduct * g = *grib;
998 fortint DefinitionNumber = (fortint) *definitionNumber;
999
1000 if( DEBUG1 ) printf("RSDEFIN: definitionNumber = %f\n",*definitionNumber);
1001
1002 return ISDEFIN(grib,&DefinitionNumber);
1003 }
1004
ISBTSPV(gribProduct ** grib,fortint * value)1005 fortint ISBTSPV(gribProduct ** grib, fortint * value) {
1006 gribProduct * g = *grib;
1007 if( DEBUG2) printf("ISBTSPV: set number of bits per value to %d\n",*value);
1008 MOVE1BYTE(((g->g4)->numberOfBitsPerValue),value);
1009 return 0;
1010 }
1011
RSBTSPV(gribProduct ** grib,fortdouble * value)1012 fortint RSBTSPV(gribProduct ** grib, fortdouble * value) {
1013 fortint Ivalue = (fortint) *value;
1014
1015 if( DEBUG2) printf("RSBTSPV\n");
1016 return ISBTSPV(grib,&Ivalue);
1017 }
1018
ISREPRS(gribProduct ** grib,fortint * value)1019 fortint ISREPRS(gribProduct ** grib, fortint * value) {
1020 gribProduct * g = *grib;
1021 fortint Value = *value, defaultValue, loop;
1022 fortint zero = 0, ninety = 90000;
1023 unsigned char * p;
1024
1025 if( DEBUG1) printf("ISREPRS: set data representation type to %d\n",Value);
1026
1027 /*
1028 //if( g2_datatype(g) == Value ) {
1029 // if( DEBUG1) printf("ISREPRS: data representation type remains unchanged\n");
1030 // return 0;
1031 //}
1032 */
1033 MOVE1BYTE(((g->g2)->NV),&zero);
1034 defaultValue = 0xff;
1035 MOVE1BYTE(((g->g2)->PV_PL),&defaultValue);
1036
1037 switch( (int) Value ) {
1038
1039 case 0:
1040 case 10:
1041 case 20:
1042 case 30:
1043 if( DEBUG1) printf("ISREPRS: setup for latitude/longitude grid\n");
1044 MOVE1BYTE(((g->g2)->dataRepresentationType),value);
1045
1046 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongParallel),&zero);
1047 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongMeridian),&zero);
1048 g->northSouthNumberOfPointsSet = g->westEastNumberOfPointsSet = FALSE;
1049
1050 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfFirstPoint),&ninety);
1051 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfFirstPoint),&zero);
1052 defaultValue = ninety | 0x800000;
1053 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfLastPoint),&defaultValue);
1054 defaultValue = 360000;
1055 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfLastPoint),&defaultValue);
1056 g->northSet = g->southSet = g->westSet = g->eastSet = TRUE;
1057
1058 defaultValue = 0x80;
1059 MOVE1BYTE(((g->g2)->grid.latlon.resolutionAndComponentsFlag),&defaultValue);
1060 MOVE2BYTES(((g->g2)->grid.latlon.iDirectionIncrement),&zero);
1061 MOVE2BYTES(((g->g2)->grid.latlon.jDirectionIncrement),&zero);
1062 g->northSouthIncrementSet = g->westEastIncrementSet = FALSE;
1063
1064 MOVE1BYTE(((g->g2)->grid.latlon.scanningMode),&zero);
1065 MOVE4BYTES(((g->g2)->grid.latlon.setToZero),&zero);
1066
1067 defaultValue = ninety | 0x800000;
1068 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfSouthPole),&defaultValue);
1069 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfSouthPole),&zero);
1070 MOVE4BYTES(((g->g2)->grid.latlon.angleOfRotationOrStretchingFactor),&zero);
1071 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfPoleOfStretching),&zero);
1072 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfPoleOfStretching),&zero);
1073 MOVE4BYTES(((g->g2)->grid.latlon.stretchingFactor),&zero);
1074 break;
1075
1076 case 50:
1077 case 60:
1078 case 70:
1079 case 80:
1080 if( DEBUG1) printf("ISREPRS: setup for spherical harmonics\n");
1081 MOVE1BYTE(((g->g2)->dataRepresentationType),value);
1082 MOVE2BYTES(((g->g2)->grid.spectral.J),&zero);
1083 MOVE2BYTES(((g->g2)->grid.spectral.K),&zero);
1084 MOVE2BYTES(((g->g2)->grid.spectral.M),&zero);
1085 defaultValue = 1;
1086 MOVE1BYTE(((g->g2)->grid.spectral.representationType),&defaultValue);
1087 MOVE1BYTE(((g->g2)->grid.spectral.representationMode),&defaultValue);
1088 p = (unsigned char *) &((g->g2)->grid.spectral.setToZero);
1089 for( loop = 0; loop < 18; loop++ ) { MOVE1BYTE(p,&zero); p++; }
1090 defaultValue = ninety | 0x800000;
1091 MOVE3BYTES(((g->g2)->grid.spectral.latitudeOfSouthPole),&defaultValue);
1092 MOVE3BYTES(((g->g2)->grid.spectral.longitudeOfSouthPole),&zero);
1093 MOVE4BYTES(((g->g2)->grid.spectral.angleOfRotationOrStretchingFactor),&zero);
1094 MOVE3BYTES(((g->g2)->grid.spectral.latitudeOfPoleOfStretching),&zero);
1095 MOVE3BYTES(((g->g2)->grid.spectral.longitudeOfPoleOfStretching),&zero);
1096 MOVE4BYTES(((g->g2)->grid.spectral.stretchingFactor),&zero);
1097 break;
1098
1099 case 4:
1100 case 14:
1101 case 24:
1102 case 34:
1103 if( DEBUG1) printf("ISREPRS: setup for gaussian grids\n");
1104 MOVE1BYTE(((g->g2)->dataRepresentationType),value);
1105
1106 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfPointsAlongParallel),&zero);
1107 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfPointsAlongMeridian),&zero);
1108 g->northSouthNumberOfPointsSet = g->westEastNumberOfPointsSet = FALSE;
1109
1110 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfFirstPoint),&ninety);
1111 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfFirstPoint),&zero);
1112 defaultValue = ninety | 0x800000;
1113 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfLastPoint),&defaultValue);
1114 defaultValue = 360000;
1115 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfLastPoint),&defaultValue);
1116 g->northSet = g->southSet = g->westSet = g->eastSet = TRUE;
1117
1118 defaultValue = 0x80;
1119 MOVE1BYTE(((g->g2)->grid.gaussian.resolutionAndComponentsFlag),&defaultValue);
1120 MOVE2BYTES(((g->g2)->grid.gaussian.iDirectionIncrement),&zero);
1121 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfParallelsBetweenPoleAndEquator),&zero);
1122 g->northSouthIncrementSet = g->westEastIncrementSet = FALSE;
1123
1124 MOVE1BYTE(((g->g2)->grid.gaussian.scanningMode),&zero);
1125 MOVE4BYTES(((g->g2)->grid.gaussian.setToZero),&zero);
1126
1127 defaultValue = ninety | 0x800000;
1128 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfSouthPole),&defaultValue);
1129 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfSouthPole),&zero);
1130 MOVE4BYTES(((g->g2)->grid.gaussian.angleOfRotationOrStretchingFactor),&zero);
1131 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfPoleOfStretching),&zero);
1132 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfPoleOfStretching),&zero);
1133 MOVE4BYTES(((g->g2)->grid.gaussian.stretchingFactor),&zero);
1134 break;
1135
1136 default:
1137 if( DEBUG1) printf("ISREPRS: data representation type not yet handled\n");
1138 break;
1139 }
1140
1141 return 0;
1142 }
1143
RSREPRS(gribProduct ** grib,fortdouble * value)1144 fortint RSREPRS(gribProduct ** grib, fortdouble * value) {
1145 fortint Ivalue = (fortint) *value;
1146 if( DEBUG2) printf("RSREPRS\n");
1147 return ISREPRS(grib,&Ivalue);
1148 }
1149
ISNWLAT(gribProduct ** grib,fortint * value)1150 fortint ISNWLAT(gribProduct ** grib, fortint * value) {
1151 gribProduct * g = *grib;
1152 fortint Value = *value;
1153 if( DEBUG2) printf("ISNWLAT: value = %d\n", Value);
1154 if( !generalLatLonGrid(g) ) return 0;
1155
1156 if( Value < 0 ) Value = 0x800000 | (-Value);
1157 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfFirstPoint),&Value);
1158 g->northSet = TRUE;
1159 adjustGridAreaDefinition(g,(int)NORTH);
1160 return 0;
1161 }
1162
RSNWLAT(gribProduct ** grib,fortdouble * value)1163 fortint RSNWLAT(gribProduct ** grib, fortdouble * value) {
1164 fortint Ivalue = (fortint) (*value * 1000.0);
1165 if( DEBUG2) printf("ISNWLAT\n");
1166 return ISNWLAT(grib,&Ivalue);
1167 }
1168
ISNWLON(gribProduct ** grib,fortint * value)1169 fortint ISNWLON(gribProduct ** grib, fortint * value) {
1170 gribProduct * g = *grib;
1171 fortint Value = *value;
1172 if( DEBUG2) printf("ISNWLON: value = %d\n", Value);
1173 if( !generalLatLonGrid(g) ) return 0;
1174
1175 if( Value < 0 ) Value = 0x800000 | (-Value);
1176 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfFirstPoint),&Value);
1177 g->westSet = TRUE;
1178 adjustGridAreaDefinition(g,(int)WEST);
1179 return 0;
1180 }
1181
RSNWLON(gribProduct ** grib,fortdouble * value)1182 fortint RSNWLON(gribProduct ** grib, fortdouble * value) {
1183 fortint Ivalue = (fortint) (*value * 1000.0);
1184 if( DEBUG2) printf("ISNWLON\n");
1185 return ISNWLON(grib,&Ivalue);
1186 }
1187
ISSELAT(gribProduct ** grib,fortint * value)1188 fortint ISSELAT(gribProduct ** grib, fortint * value) {
1189 gribProduct * g = *grib;
1190 fortint Value = *value;
1191 if( DEBUG2) printf("ISSELAT: value = %d\n", Value);
1192 if( !generalLatLonGrid(g) ) return 0;
1193
1194 if( Value < 0 ) Value = 0x800000 | (-Value);
1195 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfLastPoint),&Value);
1196 g->southSet = TRUE;
1197 adjustGridAreaDefinition(g,(int)SOUTH);
1198 return 0;
1199 }
1200
RSSELAT(gribProduct ** grib,fortdouble * value)1201 fortint RSSELAT(gribProduct ** grib, fortdouble * value) {
1202 fortint Ivalue = (fortint) (*value * 1000.0);
1203 if( DEBUG2) printf("ISSELAT\n");
1204 return ISSELAT(grib,&Ivalue);
1205 }
1206
ISSELON(gribProduct ** grib,fortint * value)1207 fortint ISSELON(gribProduct ** grib, fortint * value) {
1208 gribProduct * g = *grib;
1209 fortint Value = *value;
1210 if( DEBUG2) printf("ISSELON: value = %d\n", Value);
1211 if( !generalLatLonGrid(g) ) return 0;
1212
1213 if( Value < 0 ) Value = 0x800000 | (-Value);
1214 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfLastPoint),&Value);
1215 g->eastSet = TRUE;
1216 adjustGridAreaDefinition(g,(int)EAST);
1217 return 0;
1218 }
1219
RSSELON(gribProduct ** grib,fortdouble * value)1220 fortint RSSELON(gribProduct ** grib, fortdouble * value) {
1221 fortint Ivalue = (fortint) (*value * 1000.0);
1222 if( DEBUG2) printf("ISSELON\n");
1223 return ISSELON(grib,&Ivalue);
1224 }
1225
ISDI(gribProduct ** grib,fortint * value)1226 fortint ISDI(gribProduct ** grib, fortint * value) {
1227 gribProduct * g = *grib;
1228 if( DEBUG2) printf("ISDI: value = %d\n", *value);
1229 if( !generalLatLonGrid(g) ) return 0;
1230
1231 MOVE2BYTES(((g->g2)->grid.latlon.iDirectionIncrement), value);
1232 g->westEastIncrementSet = TRUE;
1233 adjustGridAreaDefinition(g,(int)W_E_INC);
1234 return 0;
1235 }
1236
RSDI(gribProduct ** grib,fortdouble * value)1237 fortint RSDI(gribProduct ** grib, fortdouble * value) {
1238 fortint Ivalue = (fortint) (*value * 1000.0);
1239 if( DEBUG2) printf("RSDI\n");
1240 return ISDI(grib,&Ivalue);
1241 }
1242
ISDJ(gribProduct ** grib,fortint * value)1243 fortint ISDJ(gribProduct ** grib, fortint * value) {
1244 gribProduct * g = *grib;
1245 if( DEBUG2) printf("ISDJ: value = %d\n", *value);
1246 if( !generalLatLonGrid(g) ) return 0;
1247
1248 MOVE2BYTES(((g->g2)->grid.latlon.jDirectionIncrement), value);
1249 g->northSouthIncrementSet = TRUE;
1250 adjustGridAreaDefinition(g,(int)N_S_INC);
1251 return 0;
1252 }
1253
RSDJ(gribProduct ** grib,fortdouble * value)1254 fortint RSDJ(gribProduct ** grib, fortdouble * value) {
1255 fortint Ivalue = (fortint) (*value * 1000.0);
1256 if( DEBUG2) printf("RSDJ\n");
1257 return ISDJ(grib,&Ivalue);
1258 }
1259
ISDIJ(gribProduct ** grib,fortint * value)1260 fortint ISDIJ(gribProduct ** grib, fortint * value) {
1261 gribProduct * g = *grib;
1262 if( DEBUG2) printf("ISDIJ: value = %d\n", *value);
1263 if( !generalLatLonGrid(g) ) return 0;
1264
1265 ISDI(grib,value);
1266 ISDJ(grib,value);
1267 return 0;
1268 }
1269
RSDIJ(gribProduct ** grib,fortdouble * value)1270 fortint RSDIJ(gribProduct ** grib, fortdouble * value) {
1271 fortint Ivalue = (fortint) (*value * 1000.0);
1272 if( DEBUG2) printf("RSDIJ\n");
1273 return ISDIJ(grib,&Ivalue);
1274 }
1275
ISNJ(gribProduct ** grib,fortint * value)1276 fortint ISNJ(gribProduct ** grib, fortint * value) {
1277 gribProduct * g = *grib;
1278 if( DEBUG2) printf("ISNJ: value = %d\n", *value);
1279 if( !generalLatLonGrid(g) ) return 0;
1280
1281 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongMeridian), value);
1282 g->northSouthNumberOfPointsSet = TRUE;
1283 adjustGridAreaDefinition(g,(int)N_S_PTS);
1284 return 0;
1285 }
1286
RSNJ(gribProduct ** grib,fortdouble * value)1287 fortint RSNJ(gribProduct ** grib, fortdouble * value) {
1288 fortint Ivalue = (fortint) (*value);
1289 if( DEBUG2) printf("ISNJ\n");
1290 return ISNJ(grib,&Ivalue);
1291 }
1292
ISNI(gribProduct ** grib,fortint * value)1293 fortint ISNI(gribProduct ** grib, fortint * value) {
1294 gribProduct * g = *grib;
1295 if( DEBUG2) printf("ISNI: value = %d\n", *value);
1296 if( !generalLatLonGrid(g) ) return 0;
1297
1298 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongParallel), value);
1299 g->westEastNumberOfPointsSet = TRUE;
1300 adjustGridAreaDefinition(g,(int)W_E_PTS);
1301 return 0;
1302 }
1303
RSNI(gribProduct ** grib,fortdouble * value)1304 fortint RSNI(gribProduct ** grib, fortdouble * value) {
1305 fortint Ivalue = (fortint) (*value);
1306 if( DEBUG2) printf("ISNI\n");
1307 return ISNI(grib,&Ivalue);
1308 }
1309
ISJ(gribProduct ** grib,fortint * value)1310 fortint ISJ(gribProduct ** grib, fortint * value) {
1311 gribProduct * g = *grib;
1312 fortint truncation = *value;
1313
1314 if( DEBUG2) printf("ISJ value = %d\n", truncation);
1315 if( anySpectralField(g) ) {
1316 MOVE2BYTES(((g->g2)->grid.spectral.J),value);
1317 g->numberOfValues = (truncation+1) * (truncation+2);
1318 }
1319 return 0;
1320 }
1321
RSJ(gribProduct ** grib,fortdouble * value)1322 fortint RSJ(gribProduct ** grib, fortdouble * value) {
1323 fortint Ivalue = (fortint) (*value);
1324 if( DEBUG2) printf("RSJ value = %f\n", *value);
1325 return ISJ(grib,&Ivalue);
1326 }
1327
ISK(gribProduct ** grib,fortint * value)1328 fortint ISK(gribProduct ** grib, fortint * value) {
1329 gribProduct * g = *grib;
1330 if( DEBUG2) printf("ISK value = %d\n", *value);
1331 if( anySpectralField(g) ) MOVE2BYTES(((g->g2)->grid.spectral.K),value);
1332 return 0;
1333 }
1334
RSK(gribProduct ** grib,fortdouble * value)1335 fortint RSK(gribProduct ** grib, fortdouble * value) {
1336 fortint Ivalue = (fortint) (*value);
1337 if( DEBUG2) printf("RSK value = %f\n", *value);
1338 return ISK(grib,&Ivalue);
1339 }
1340
ISM(gribProduct ** grib,fortint * value)1341 fortint ISM(gribProduct ** grib, fortint * value) {
1342 gribProduct * g = *grib;
1343 if( DEBUG2) printf("ISM value = %d\n", *value);
1344 if( anySpectralField(g) ) MOVE2BYTES(((g->g2)->grid.spectral.M),value);
1345 return 0;
1346 }
1347
RSM(gribProduct ** grib,fortdouble * value)1348 fortint RSM(gribProduct ** grib, fortdouble * value) {
1349 fortint Ivalue = (fortint) (*value);
1350 if( DEBUG2) printf("RSM value = %f\n", *value);
1351 return ISM(grib,&Ivalue);
1352 }
1353
ISJKM(gribProduct ** grib,fortint * value)1354 fortint ISJKM(gribProduct ** grib, fortint * value) {
1355 gribProduct * g = *grib;
1356 if( DEBUG2) printf("ISJKM value = %d\n", *value);
1357 ISJ(grib,value);
1358 ISK(grib,value);
1359 ISM(grib,value);
1360 return 0;
1361 }
1362
RSJKM(gribProduct ** grib,fortdouble * value)1363 fortint RSJKM(gribProduct ** grib, fortdouble * value) {
1364 fortint Ivalue = (fortint) (*value);
1365 if( DEBUG2) printf("RSJKM value = %f\n", *value);
1366 return ISJKM(grib,&Ivalue);
1367 }
1368
ISTJ(gribProduct ** grib,fortint * value)1369 fortint ISTJ(gribProduct ** grib, fortint * value) {
1370 gribProduct * g = *grib;
1371 if( DEBUG2) printf("ISTJ value = %d\n", *value);
1372 if( anyComplexPackedSpectralField(g) )
1373 MOVE1BYTE(((g->g4)->data.complexSpectral.J),value);
1374 return 0;
1375 }
1376
RSTJ(gribProduct ** grib,fortdouble * value)1377 fortint RSTJ(gribProduct ** grib, fortdouble * value) {
1378 fortint Ivalue = (fortint) (*value);
1379 if( DEBUG2) printf("RSTJ value = %f\n", *value);
1380 return ISTJ(grib,&Ivalue);
1381 }
1382
ISTK(gribProduct ** grib,fortint * value)1383 fortint ISTK(gribProduct ** grib, fortint * value) {
1384 gribProduct * g = *grib;
1385 if( DEBUG2) printf("ISTK value = %d\n", *value);
1386 if( anyComplexPackedSpectralField(g) )
1387 MOVE1BYTE(((g->g4)->data.complexSpectral.K),value);
1388 return 0;
1389 }
1390
RSTK(gribProduct ** grib,fortdouble * value)1391 fortint RSTK(gribProduct ** grib, fortdouble * value) {
1392 fortint Ivalue = (fortint) (*value);
1393 if( DEBUG2) printf("RSTK value = %f\n", *value);
1394 return ISTK(grib,&Ivalue);
1395 }
1396
ISTM(gribProduct ** grib,fortint * value)1397 fortint ISTM(gribProduct ** grib, fortint * value) {
1398 gribProduct * g = *grib;
1399 if( DEBUG2) printf("ISTM value = %d\n", *value);
1400 if( anyComplexPackedSpectralField(g) )
1401 MOVE1BYTE(((g->g4)->data.complexSpectral.M),value);
1402 return 0;
1403 }
1404
RSTM(gribProduct ** grib,fortdouble * value)1405 fortint RSTM(gribProduct ** grib, fortdouble * value) {
1406 fortint Ivalue = (fortint) (*value);
1407 if( DEBUG2) printf("RSTM value = %f\n", *value);
1408 return ISTM(grib,&Ivalue);
1409 }
1410
ISTJKM(gribProduct ** grib,fortint * value)1411 fortint ISTJKM(gribProduct ** grib, fortint * value) {
1412 gribProduct * g = *grib;
1413 if( DEBUG2) printf("ISTJKM value = %d\n", *value);
1414 ISTJ(grib,value);
1415 ISTK(grib,value);
1416 ISTM(grib,value);
1417 return 0;
1418 }
1419
RSTJKM(gribProduct ** grib,fortdouble * value)1420 fortint RSTJKM(gribProduct ** grib, fortdouble * value) {
1421 fortint Ivalue = (fortint) (*value);
1422 if( DEBUG2) printf("RSTJKM value = %f\n", *value);
1423 return ISTJKM(grib,&Ivalue);
1424 }
1425
ISLATRP(gribProduct ** grib,fortint * value)1426 fortint ISLATRP(gribProduct ** grib, fortint * value) {
1427 gribProduct * g = *grib;
1428 fortint Value = *value;
1429 if( DEBUG2) printf("ISLATRP: value = %d\n", Value);
1430 if( !generalRotatedGrid(g) ) return 0;
1431
1432 if( Value < 0 ) Value = 0x800000 | (-Value);
1433 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfSouthPole),&Value);
1434 return 0;
1435 }
1436
RSLATRP(gribProduct ** grib,fortdouble * value)1437 fortint RSLATRP(gribProduct ** grib, fortdouble * value) {
1438 fortint Ivalue = (fortint) (*value * 1000.0);
1439 if( DEBUG2) printf("RSLATRP\n");
1440 return ISLATRP(grib,&Ivalue);
1441 }
1442
ISLONRP(gribProduct ** grib,fortint * value)1443 fortint ISLONRP(gribProduct ** grib, fortint * value) {
1444 gribProduct * g = *grib;
1445 fortint Value = *value;
1446 if( DEBUG2) printf("ISLONRP: value = %d\n", Value);
1447 if( !generalRotatedGrid(g) ) return 0;
1448
1449 if( Value < 0 ) Value = 0x800000 | (-Value);
1450 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfSouthPole),&Value);
1451 return 0;
1452 }
1453
RSLONRP(gribProduct ** grib,fortdouble * value)1454 fortint RSLONRP(gribProduct ** grib, fortdouble * value) {
1455 fortint Ivalue = (fortint) (*value * 1000.0);
1456 if( DEBUG2) printf("RSLONRP\n");
1457 return ISLONRP(grib,&Ivalue);
1458 }
1459
RSROTAT(gribProduct ** grib,fortdouble * value)1460 fortint RSROTAT(gribProduct ** grib, fortdouble * value) {
1461 gribProduct * g = *grib;
1462 fortint valueExponent, valueMantissa, numberOfBitsPerInteger;
1463 fortint status;
1464
1465 if( DEBUG2) printf("RSROTAT: value = %f\n", *value);
1466 if( !generalRotatedGrid(g) ) return 0;
1467
1468 numberOfBitsPerInteger = sizeof(fortint) * 8;
1469 status = REF2GRB(value,&valueExponent,
1470 &valueMantissa,&numberOfBitsPerInteger);
1471 if( status ) {
1472 printf("RSROTAT: call to REF2GRB failed\n");
1473 exit(1);
1474 }
1475 MOVE1BYTE(((g->g2)->grid.latlon.angleOfRotationOrStretchingFactor),&valueExponent);
1476 MOVE3BYTES((((g->g2)->grid.latlon.angleOfRotationOrStretchingFactor)+1),&valueMantissa);
1477 return 0;
1478 }
1479
ISROTAT(gribProduct ** grib,fortint * value)1480 fortint ISROTAT(gribProduct ** grib, fortint * value) {
1481 fortdouble Value = ((fortdouble) *value) / 1000.0;
1482 if( DEBUG2) printf("ISROTAT: value = %d\n", *value);
1483 return RSROTAT(grib,&Value);
1484 }
1485
ISLATSP(gribProduct ** grib,fortint * value)1486 fortint ISLATSP(gribProduct ** grib, fortint * value) {
1487 gribProduct * g = *grib;
1488 fortint Value = *value;
1489 if( DEBUG2) printf("ISLATSP: value = %d\n", Value);
1490 if( !generalStretchedGrid(g) ) return 0;
1491
1492 if( Value < 0 ) Value = 0x800000 | (-Value);
1493
1494 if( generalStretchedAndRotatedGrid(g) )
1495 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfPoleOfStretching),&Value);
1496 else
1497 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfSouthPole),&Value);
1498
1499 return 0;
1500 }
1501
RSLATSP(gribProduct ** grib,fortdouble * value)1502 fortint RSLATSP(gribProduct ** grib, fortdouble * value) {
1503 fortint Ivalue = (fortint) (*value * 1000.0);
1504 if( DEBUG2) printf("RSLATSP\n");
1505 return ISLATSP(grib,&Ivalue);
1506 }
1507
ISLONSP(gribProduct ** grib,fortint * value)1508 fortint ISLONSP(gribProduct ** grib, fortint * value) {
1509 gribProduct * g = *grib;
1510 fortint Value = *value;
1511 if( DEBUG2) printf("ISLONSP: value = %d\n", Value);
1512 if( !generalStretchedGrid(g) ) return 0;
1513
1514 if( Value < 0 ) Value = 0x800000 | (-Value);
1515
1516 if( generalStretchedAndRotatedGrid(g) )
1517 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfPoleOfStretching),&Value);
1518 else
1519 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfSouthPole),&Value);
1520
1521 return 0;
1522 }
1523
RSLONSP(gribProduct ** grib,fortdouble * value)1524 fortint RSLONSP(gribProduct ** grib, fortdouble * value) {
1525 fortint Ivalue = (fortint) (*value * 1000.0);
1526 if( DEBUG2) printf("RSLONSP\n");
1527 return ISLONSP(grib,&Ivalue);
1528 }
1529
RSSFACT(gribProduct ** grib,fortdouble * value)1530 fortint RSSFACT(gribProduct ** grib, fortdouble * value) {
1531 gribProduct * g = *grib;
1532 fortint valueExponent, valueMantissa, numberOfBitsPerInteger;
1533 fortint status;
1534
1535 if( DEBUG2) printf("RSSFACT: value = %f\n", *value);
1536 if( !generalStretchedGrid(g) ) return 0;
1537
1538 numberOfBitsPerInteger = sizeof(fortint) * 8;
1539 status = REF2GRB(value,&valueExponent,
1540 &valueMantissa,&numberOfBitsPerInteger);
1541 if( status ) {
1542 printf("RSSFACT: call to REF2GRB failed\n");
1543 exit(1);
1544 }
1545
1546 if( generalStretchedAndRotatedGrid(g) ) {
1547 MOVE1BYTE(((g->g2)->grid.latlon.stretchingFactor),&valueExponent);
1548 MOVE3BYTES((((g->g2)->grid.latlon.stretchingFactor)+1),&valueMantissa);
1549 }
1550 else {
1551 MOVE1BYTE(((g->g2)->grid.latlon.angleOfRotationOrStretchingFactor),&valueExponent);
1552 MOVE3BYTES((((g->g2)->grid.latlon.angleOfRotationOrStretchingFactor)+1),&valueMantissa);
1553 }
1554 return 0;
1555 }
1556
ISSFACT(gribProduct ** grib,fortint * value)1557 fortint ISSFACT(gribProduct ** grib, fortint * value) {
1558 fortdouble Value = ((fortdouble) *value) / 1000.0;
1559 if( DEBUG2) printf("ISSFACT: value = %d\n", *value);
1560 return RSSFACT(grib,&Value);
1561 }
1562
ISSETQG(gribProduct ** grib,fortint * value)1563 fortint ISSETQG(gribProduct ** grib, fortint * value) {
1564 gribProduct * g = *grib;
1565 fortint northToSouth;
1566 fortint number = (*value), missing = 0xffff, zero = 0;
1567 fortdouble * gaussianLatitudes;
1568 fortint * countOfPointsAtLatitudes;
1569 fortint total, loop, status, one = 1, coordinate;
1570 unsigned char htype[2] = "R";
1571 fortdouble lastLongitude;
1572 gribSection2 * newSection2;
1573 fortint newSection2Size, listOffset, count;
1574 unsigned char * p;
1575
1576 if( DEBUG2) printf("ISSETQG: reduced gaussian number = %d\n", *value);
1577
1578 if( !anyGaussianGrid(g) ) return 0;
1579 /*
1580 //if( number == g2_gaussNumber(g) ) return 0;
1581 */
1582
1583 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfPointsAlongParallel),&missing);
1584 northToSouth = number * 2;
1585 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfPointsAlongMeridian),&northToSouth);
1586
1587 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfFirstPoint),&zero);
1588
1589 MOVE1BYTE(((g->g2)->grid.gaussian.resolutionAndComponentsFlag),&zero);
1590
1591 MOVE2BYTES(((g->g2)->grid.gaussian.iDirectionIncrement),&missing);
1592 MOVE2BYTES(((g->g2)->grid.gaussian.numberOfParallelsBetweenPoleAndEquator),
1593 &number);
1594 MOVE1BYTE(((g->g2)->grid.gaussian.scanningMode),&zero);
1595
1596 MOVE4BYTES(((g->g2)->grid.gaussian.setToZero),&zero);
1597 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfSouthPole),&zero);
1598 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfSouthPole),&zero);
1599 MOVE4BYTES(((g->g2)->grid.gaussian.angleOfRotationOrStretchingFactor),&zero);
1600 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfPoleOfStretching),&zero);
1601 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfPoleOfStretching),&zero);
1602 MOVE4BYTES(((g->g2)->grid.gaussian.stretchingFactor),&zero);
1603
1604 gaussianLatitudes =
1605 (fortdouble *) allocateMemory(northToSouth*sizeof(fortdouble));
1606 countOfPointsAtLatitudes =
1607 (fortint *) allocateMemory(northToSouth*sizeof(fortint));
1608 JGETGG(&number,htype,gaussianLatitudes,countOfPointsAtLatitudes,&status,
1609 (long)one);
1610 if( status ) {
1611 printf("ISSETQG: JGETGG failed for gaussian number %d\n",number);
1612 exit(1);
1613 }
1614
1615 coordinate = (fortint) ((0.0005 + gaussianLatitudes[0])*1000.0);
1616 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfFirstPoint),&coordinate);
1617
1618 coordinate = 0x800000 | coordinate;
1619 MOVE3BYTES(((g->g2)->grid.gaussian.latitudeOfLastPoint),&coordinate);
1620
1621 freeMemory(gaussianLatitudes);
1622
1623 lastLongitude = 360.0 - (360/(fortdouble)countOfPointsAtLatitudes[number-1]);
1624 coordinate = (fortint) (lastLongitude * 1000.0);
1625 if( coordinate < 0 ) coordinate = 0x800000 | (-coordinate);
1626 MOVE3BYTES(((g->g2)->grid.gaussian.longitudeOfLastPoint),&coordinate);
1627
1628 switch( (int) g2_datatype(g) ) {
1629
1630 case 4:
1631 listOffset = 32;
1632 break;
1633
1634 case 14:
1635 case 24:
1636 listOffset = 42;
1637 break;
1638
1639 case 34:
1640 listOffset = 52;
1641 break;
1642
1643 default:
1644 printf("ISSETQG: unexpected gaussian representation type = %d\n",
1645 g2_datatype(g));
1646 exit(1);
1647 }
1648 listOffset += g2_NV(g) * 4;
1649
1650 newSection2Size = listOffset + (northToSouth * 2);
1651 newSection2 = (gribSection2 *) allocateMemory(newSection2Size);
1652 memcpy(newSection2,(g->g2),sizeof(gribSection2));
1653
1654 total = 0;
1655 p = (unsigned char *) newSection2;
1656 p += listOffset;
1657 for(loop=0;loop<northToSouth;loop++) {
1658 count = countOfPointsAtLatitudes[loop];
1659 total += count;
1660 MOVE2BYTES(p,&count);
1661 p +=2;
1662 }
1663 g->numberOfValues = total;
1664 freeMemory((g->g2));
1665 g->g2 = newSection2;
1666
1667 MOVE3BYTES(((g->g2)->sectionLength),&newSection2Size);
1668 listOffset++;
1669 MOVE1BYTE(((g->g2)->PV_PL),&listOffset);
1670
1671 return 0;
1672 }
1673
RSSETQG(gribProduct ** grib,fortdouble * value)1674 fortint RSSETQG(gribProduct ** grib, fortdouble * value) {
1675 fortint Value = (fortint) *value;
1676 if( DEBUG2) printf("RSSETQG\n");
1677 return ISSETQG(grib,&Value);;
1678 }
1679
ISSETRG(gribProduct ** grib,fortint * value)1680 fortint ISSETRG(gribProduct ** grib, fortint * value) {
1681 if( DEBUG2) printf("ISSETRG: regular gaussian number = %d\n", *value);
1682 return 0;
1683 }
1684
RSSETRG(gribProduct ** grib,fortdouble * value)1685 fortint RSSETRG(gribProduct ** grib, fortdouble * value) {
1686 if( DEBUG2) printf("RSSETRG\n");
1687 return 0;
1688 }
1689
ISDUMMY(gribProduct ** grib,fortint * value)1690 fortint ISDUMMY(gribProduct ** grib, fortint * value) {
1691 if( DEBUG2) printf("ISDUMMY\n");
1692 return 0;
1693 }
1694
RSDUMMY(gribProduct ** grib,fortdouble * value)1695 fortint RSDUMMY(gribProduct ** grib, fortdouble * value) {
1696 if( DEBUG2) printf("RSDUMMY\n");
1697 return 0;
1698 }
1699
SVALUES(gribProduct ** grib,fortdouble * arrayOfUnpackedValues,fortint * arraySize,fortdouble * userSuppliedMissingValue)1700 fortint SVALUES(
1701 gribProduct ** grib,
1702 fortdouble * arrayOfUnpackedValues,
1703 fortint * arraySize,
1704 fortdouble * userSuppliedMissingValue ) {
1705 gribProduct * g = *grib;
1706 fortint loop, startValue, lengthOfSection3;
1707 fortint N, numberOfUnpackedValues, nextNonMissingValue = 0;
1708 fortint numberOfValues, numberOfMissingValues = 0;
1709 fortdouble missingValue = *userSuppliedMissingValue;
1710 fortdouble value, maximum, minimum, truncation;
1711 unsigned char * bitmap, * section3;
1712
1713 if( DEBUG1 ) printf("SVALUES\n");
1714 /*
1715 // Calculate the number of points from the GRIB headers
1716 */
1717 numberOfUnpackedValues = g->numberOfValues;
1718 if( DEBUG1 )
1719 printf("SVALUES: numberOfUnpackedValues = %d\n",numberOfUnpackedValues);
1720
1721 if( numberOfUnpackedValues <= 0 ) {
1722 printf("SVALUES: grib headers not sufficiently configured for packing\n");
1723 return -1;
1724 }
1725
1726 if( numberOfUnpackedValues > *arraySize ) {
1727 printf(
1728 "SVALUES: calculated number of values (%d) greater than array size(%d)\n",
1729 numberOfUnpackedValues, *arraySize);
1730 return -1;
1731 }
1732
1733 N = numberOfUnpackedValues;
1734 /*
1735 // For gridpoint fields, it may be necessary to create a bitmap if
1736 // there are missing data values
1737 */
1738 if( !anySpectralField(g) ) {
1739 if( DEBUG1 ) printf("SVALUES: gridpoint field\n");
1740 lengthOfSection3 = (6+(N+7)/8);
1741 if( MOD(lengthOfSection3,2) == 1 ) lengthOfSection3++;
1742 section3 = (unsigned char *) allocateMemory(lengthOfSection3);
1743 bitmap = section3 + 6;
1744 for( loop = 0; loop < (N+7)/8; loop++) bitmap[loop] = 0;
1745
1746 maximum = minimum = arrayOfUnpackedValues[0];
1747 for( loop = 0; loop < N; loop++) {
1748 value = arrayOfUnpackedValues[loop];
1749 if( value == missingValue ) {
1750 numberOfMissingValues++;
1751 }
1752 else {
1753 setBitMap(bitmap,loop);
1754 if( (value < minimum) || (minimum == missingValue) ) {
1755 minimum = value;
1756 }
1757 if( (value > maximum) || (maximum == missingValue) ) {
1758 maximum = value;
1759 }
1760 }
1761 }
1762 /*
1763 // Create section 3 if there is a bitmap
1764 */
1765 if( numberOfMissingValues > 0 ) {
1766 fortint numberOfUnusedBits = (lengthOfSection3-6)*8 - N;
1767 fortint zero = 0;
1768
1769 if( DEBUG1 ) printf("SVALUES: create section 3 bitmap\n");
1770
1771 g->bitmapped = 1;
1772 freeMemory(g->g3);
1773 g->g3 = (gribSection3 *) section3;
1774 MOVE3BYTES((g->g3)->sectionLength,&lengthOfSection3);
1775 MOVE1BYTE((g->g3)->numberOfUnusedBits,&numberOfUnusedBits);
1776 MOVE2BYTES((g->g3)->tableReference,&zero);
1777 *((g->g1)->section2and3PresentFlag) |= 0x40;
1778 /*
1779 // Repack the input array, leaving out the missing data
1780 */
1781 nextNonMissingValue = 0;
1782 for( loop = 0; loop < N; loop++) {
1783 if( arrayOfUnpackedValues[loop] != missingValue )
1784 arrayOfUnpackedValues[nextNonMissingValue++] =
1785 arrayOfUnpackedValues[loop];
1786 }
1787 }
1788 else {
1789 g->bitmapped = 0;
1790 *((g->g1)->section2and3PresentFlag) &= 0xbf;
1791 freeMemory(section3);
1792 }
1793 numberOfValues = N - numberOfMissingValues;
1794
1795 if( simplePacking(g) ) {
1796 fortint estimatedSizeOfSection4;
1797 unsigned char * section4;
1798 fortint status, lengthOfSection4, numberOfUnusedBits;
1799 fortint referenceValueExponent, referenceValueMantissa;
1800 fortint numberOfBitsPerPackedValue, numberOfBitsPerInteger;
1801 fortint firstCoefficientExponent, firstCoefficientMantissa;
1802 fortint integerScale = 0;
1803 fortdouble range = maximum - minimum;
1804 fortdouble realScale = 1.0;
1805 fortint numberOfValuesToScale, startOfPackedValueBits;
1806 fortint * arrayOfScaledValues;
1807
1808 if( DEBUG1 ) printf("SVALUES: gridpoint field with simple packing\n");
1809
1810 estimatedSizeOfSection4 = 512+(N*sizeof(fortint)*8+7)/8;
1811 section4 = (unsigned char *) allocateMemory(estimatedSizeOfSection4);
1812
1813
1814 numberOfBitsPerPackedValue = g4_bits(g);
1815 numberOfBitsPerInteger = sizeof(fortint) * 8;
1816 if( numberOfBitsPerPackedValue == numberOfBitsPerInteger)
1817 numberOfBitsPerPackedValue--;
1818
1819 if( PRACTICALLYZERO(range) )
1820 integerScale = 0;
1821 else
1822 if( (minimum != 0) && PRACTICALLYZERO(range/minimum) ){
1823 integerScale = 0;
1824 }
1825 else
1826 if( PRACTICALLYZERO(range-1.0) )
1827 integerScale = 1 - numberOfBitsPerPackedValue;
1828 else
1829 if( range > 1.0 ) {
1830 for( loop = 1; loop < 128; loop++ ) {
1831 realScale *= 2.0;
1832 if( realScale > (range + SMALL) ) {
1833 integerScale = loop - numberOfBitsPerPackedValue;
1834 break;
1835 }
1836 }
1837 }
1838 else {
1839 for( loop = 1; loop < 127; loop++ ) {
1840 realScale /= 2.0;
1841 if( realScale < (range - SMALL) ) {
1842 integerScale = 1 - loop - numberOfBitsPerPackedValue;
1843 break;
1844 }
1845 }
1846 }
1847
1848 realScale = (fortdouble) pow((double)2.0,(double)integerScale);
1849 if( integerScale < 0 ) integerScale = 0x8000 | (-integerScale);
1850 MOVE2BYTES((section4+4),&integerScale);
1851
1852 status = REF2GRB(&minimum,&referenceValueExponent,
1853 &referenceValueMantissa,&numberOfBitsPerInteger);
1854 if( status ) {
1855 printf("SVALUES: call to REF2GRB failed\n");
1856 return -1;
1857 }
1858 MOVE1BYTE((section4+6),&referenceValueExponent);
1859 MOVE3BYTES((section4+7),&referenceValueMantissa);
1860
1861 MOVE1BYTE((section4+10),&numberOfBitsPerPackedValue);
1862
1863 status = REF2GRB(&arrayOfUnpackedValues[0],&firstCoefficientExponent,
1864 &firstCoefficientMantissa,&numberOfBitsPerInteger);
1865 if( status ) {
1866 printf("SVALUES: call to REF2GRB failed\n");
1867 return -1;
1868 }
1869 MOVE1BYTE((section4+11),&firstCoefficientExponent);
1870 MOVE3BYTES((section4+12),&firstCoefficientMantissa);
1871
1872 numberOfValuesToScale = numberOfValues;
1873 arrayOfScaledValues =
1874 (fortint *) allocateMemory(numberOfValuesToScale*sizeof(fortint));
1875
1876 INSCAL(&arrayOfUnpackedValues[0],arrayOfScaledValues,
1877 &numberOfValuesToScale,&minimum,&realScale,
1878 &numberOfBitsPerPackedValue);
1879
1880 startOfPackedValueBits = 88;
1881 lengthOfSection4 = estimatedSizeOfSection4/sizeof(fortint);
1882 INXBIT((fortint*)section4,&lengthOfSection4,
1883 &startOfPackedValueBits,arrayOfScaledValues,
1884 &numberOfValuesToScale,&numberOfBitsPerInteger,
1885 &numberOfBitsPerPackedValue,"C",&status,(long)1);
1886 if( status ) {
1887 printf("SVALUES: call to INXBIT failed\n");
1888 return -1;
1889 }
1890 freeMemory(arrayOfScaledValues);
1891
1892 lengthOfSection4 = (startOfPackedValueBits + 7)/8;
1893 if( MOD(lengthOfSection4,2) == 1 ) lengthOfSection4++;
1894 MOVE3BYTES(section4,&lengthOfSection4);
1895
1896 numberOfUnusedBits = lengthOfSection4*8 - startOfPackedValueBits;
1897 if( !floatingPoint(g) ) numberOfUnusedBits |= 0x20;
1898 MOVE1BYTE((section4+3),&numberOfUnusedBits);
1899
1900 freeMemory(g->g4);
1901 g->g4 = (gribSection4 *) section4;
1902 }
1903 else {
1904 printf("Second order packing not yet handled!\n");
1905 return -1;
1906 }
1907 return 0;
1908 }
1909
1910 /*
1911 // Spherical harmonics
1912 */
1913 if( anySpectralField(g) && floatingPoint(g) ) {
1914 fortint estimatedSizeOfSection4;
1915 unsigned char * section4;
1916
1917 estimatedSizeOfSection4 = 512+(N*sizeof(fortint)*8+7)/8;
1918 section4 = (unsigned char *) allocateMemory(estimatedSizeOfSection4);
1919
1920 if( simplePacking(g) )
1921 startValue = 1;
1922 else
1923 startValue = 0;
1924 maximum = minimum = arrayOfUnpackedValues[startValue];
1925 for( loop = startValue; loop < N; loop++) {
1926 value = arrayOfUnpackedValues[loop];
1927 if( (value < minimum) || (minimum == missingValue) ) {
1928 minimum = value;
1929 }
1930 if( (value > maximum) || (maximum == missingValue) ) {
1931 maximum = value;
1932 }
1933 }
1934
1935 /*
1936 // If complex packing required, use GRIBEX routines GRSMKP and CSECT4
1937 */
1938 if( !simplePacking(g) ) {
1939 fortint * ksec1, * ksec4;
1940 fortint ktrunc, kleng, knspt, status, adjustment;
1941 fortint numberOfBitsPerInteger, numberOfBitsPerPackedValue;
1942
1943 if( DEBUG1 ) printf("SVALUES: spectral field with complex packing\n");
1944
1945 ksec1 = (fortint *) allocateMemory(1024*sizeof(fortint));
1946 for( loop = 0; loop < 1024; loop++) ksec1[loop] = 0;
1947
1948 ksec4 = (fortint *) allocateMemory(512*sizeof(fortint));
1949 for( loop = 0; loop < 512; loop++) ksec4[loop] = 0;
1950
1951 ksec1[22] = g1_scale(g);
1952
1953 ksec4[ 1] = g4_bits(g);
1954 ksec4[16] = g4_ip(g);
1955 ksec4[17] = g4_j(g);
1956 ksec4[18] = g4_k(g);
1957 ksec4[19] = g4_m(g);
1958
1959 ktrunc = g2_J(g);
1960 kleng = (estimatedSizeOfSection4+sizeof(fortint)-1)/sizeof(fortint);
1961 knspt = 0;
1962 numberOfBitsPerInteger = sizeof(fortint) * 8;
1963 numberOfBitsPerPackedValue = g4_bits(g);
1964
1965 {
1966 fortint forceComputation = 1;
1967 GRSMKP(&forceComputation);
1968 }
1969 status = CSECT4(arrayOfUnpackedValues,&ktrunc,ksec1,ksec4,
1970 (fortint *)section4,&kleng,&knspt,
1971 &numberOfBitsPerInteger,&numberOfBitsPerPackedValue);
1972 if( status ) {
1973 printf("SVALUES: call to CSECT4 failed\n");
1974 return -1;
1975 }
1976 /*
1977 // NB. Adjust the octet number at which first-order packed data begins
1978 // to be consistent with an ECMWF bug.
1979 */
1980 adjustment = TWOBYTEINT(section4+11);
1981 adjustment += 8 + g1_length(g) + g2_length(g);
1982 MOVE2BYTES((section4+11),&adjustment);
1983
1984 freeMemory(g->g4);
1985 g->g4 = (gribSection4 *) section4;
1986 }
1987 /*
1988 // Simple packing required
1989 */
1990 else {
1991 fortint status, lengthOfSection4, numberOfUnusedBits;
1992 fortint referenceValueExponent, referenceValueMantissa;
1993 fortint numberOfBitsPerPackedValue, numberOfBitsPerInteger;
1994 fortint firstCoefficientExponent, firstCoefficientMantissa;
1995 fortint integerScale = 0;
1996 fortdouble range = maximum - minimum;
1997 fortdouble realScale = 1.0;
1998 fortint numberOfValuesToScale, startOfPackedValueBits;
1999 fortint * arrayOfScaledValues;
2000
2001 if( DEBUG1 ) printf("SVALUES: spectral field with simple packing\n");
2002
2003 numberOfBitsPerPackedValue = g4_bits(g);
2004 numberOfBitsPerInteger = sizeof(fortint) * 8;
2005 if( numberOfBitsPerPackedValue == numberOfBitsPerInteger)
2006 numberOfBitsPerPackedValue--;
2007
2008 if( PRACTICALLYZERO(range) )
2009 integerScale = 0;
2010 else
2011 if( (minimum != 0) && PRACTICALLYZERO(range/minimum) )
2012 integerScale = 0;
2013 else
2014 if( PRACTICALLYZERO(range-1.0) )
2015 integerScale = 1 - numberOfBitsPerPackedValue;
2016 else
2017 if( range > 1.0 ) {
2018 for( loop = 1; loop < 128; loop++ ) {
2019 realScale *= 2.0;
2020 if( realScale > (range + SMALL) ) {
2021 integerScale = loop - numberOfBitsPerPackedValue;
2022 break;
2023 }
2024 }
2025 }
2026 else {
2027 for( loop = 1; loop < 127; loop++ ) {
2028 realScale /= 2.0;
2029 if( realScale < (range - SMALL) ) {
2030 integerScale = 1 - loop - numberOfBitsPerPackedValue;
2031 break;
2032 }
2033 }
2034 }
2035
2036 realScale = (fortdouble) pow((double)2.0,(double)integerScale);
2037 if( integerScale < 0 ) integerScale = 0x8000 | (-integerScale);
2038 MOVE2BYTES((section4+4),&integerScale);
2039
2040 status = REF2GRB(&minimum,&referenceValueExponent,
2041 &referenceValueMantissa,&numberOfBitsPerInteger);
2042 if( status ) {
2043 printf("SVALUES: call to REF2GRB failed\n");
2044 return -1;
2045 }
2046 MOVE1BYTE((section4+6),&referenceValueExponent);
2047 MOVE3BYTES((section4+7),&referenceValueMantissa);
2048
2049 MOVE1BYTE((section4+10),&numberOfBitsPerPackedValue);
2050
2051 status = REF2GRB(&arrayOfUnpackedValues[0],&firstCoefficientExponent,
2052 &firstCoefficientMantissa,&numberOfBitsPerInteger);
2053 if( status ) {
2054 printf("SVALUES: call to REF2GRB failed\n");
2055 return -1;
2056 }
2057 MOVE1BYTE((section4+11),&firstCoefficientExponent);
2058 MOVE3BYTES((section4+12),&firstCoefficientMantissa);
2059
2060 numberOfValuesToScale = N - 1;
2061 arrayOfScaledValues =
2062 (fortint *) allocateMemory(numberOfValuesToScale*sizeof(fortint));
2063
2064 INSCAL(&arrayOfUnpackedValues[1],arrayOfScaledValues,
2065 &numberOfValuesToScale,&minimum,&realScale,
2066 &numberOfBitsPerPackedValue);
2067
2068 startOfPackedValueBits = 120;
2069 lengthOfSection4 = estimatedSizeOfSection4/sizeof(fortint);
2070 INXBIT((fortint*)section4,&lengthOfSection4,
2071 &startOfPackedValueBits,arrayOfScaledValues,
2072 &numberOfValuesToScale,&numberOfBitsPerInteger,
2073 &numberOfBitsPerPackedValue,"C",&status,(long)1);
2074 if( status ) {
2075 printf("SVALUES: call to INXBIT failed\n");
2076 return -1;
2077 }
2078 freeMemory(arrayOfScaledValues);
2079
2080 lengthOfSection4 = (startOfPackedValueBits + 7)/8;
2081 if( MOD(lengthOfSection4,2) == 1 ) lengthOfSection4++;
2082 MOVE3BYTES(section4,&lengthOfSection4);
2083
2084 numberOfUnusedBits = lengthOfSection4*8 - startOfPackedValueBits;
2085 numberOfUnusedBits |= 0x80;
2086 MOVE1BYTE((section4+3),&numberOfUnusedBits);
2087
2088 freeMemory(g->g4);
2089 g->g4 = (gribSection4 *) section4;
2090 }
2091 }
2092 else {
2093 if( DEBUG1 ) printf("SVALUES: gridpoint field\n");
2094 }
2095 return 0;
2096 }
2097
setBitMap(unsigned char * bitmap,fortint bitNumber)2098 void setBitMap(unsigned char * bitmap, fortint bitNumber ) {
2099 fortint byteOffset = (bitNumber>>3);
2100 fortint bitOffset = bitNumber - (byteOffset<<3);
2101
2102 bitmap[byteOffset] |= (0x80>>bitOffset);
2103 }
2104
adjustGridAreaDefinition(gribProduct * g,int type)2105 void adjustGridAreaDefinition(gribProduct* g, int type) {
2106 fortdouble north, south, west, east, increment;
2107 fortint value, numberOfPoints ;
2108
2109 switch( type) {
2110
2111 case NORTH:
2112 if( DEBUG2 ) printf("adjustGridAreaDefinition: new north given\n");
2113 north = RGNWLAT(&g);
2114 if( g->southSet ) {
2115 south = RGSELAT(&g);
2116 if( g->northSouthIncrementSet ) {
2117 increment = RGDJ(&g);
2118 numberOfPoints = (fortint) (0.5 + (north-south)/increment) + 1;
2119 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongMeridian),
2120 &numberOfPoints);
2121 g->northSouthNumberOfPointsSet = TRUE;
2122 if( DEBUG2 )
2123 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2124 numberOfPoints);
2125 }
2126 else if( g->northSouthNumberOfPointsSet ) {
2127 numberOfPoints = IGNJ(&g) - 1;
2128 increment = (north-south)/(fortdouble) numberOfPoints;
2129 value = (fortint) (increment * 1000.0);
2130 MOVE2BYTES(((g->g2)->grid.latlon.jDirectionIncrement),&value);
2131 g->northSouthIncrementSet = TRUE;
2132 if( DEBUG2 )
2133 printf("adjustGridAreaDefinition: north-south increment set = %f\n",
2134 increment);
2135 }
2136 }
2137 else {
2138 if( g->northSouthIncrementSet && g->northSouthNumberOfPointsSet ) {
2139 increment = RGDJ(&g);
2140 numberOfPoints = IGNJ(&g) - 1;
2141 south = north - increment * (fortdouble) numberOfPoints;
2142 value = (fortint) (south * 1000.0);
2143 if( value < 0 ) value = 0x800000 | (-value);
2144 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfLastPoint),&value);
2145 g->southSet = TRUE;
2146 if( DEBUG2 )
2147 printf("adjustGridAreaDefinition: south set = %f\n", south);
2148 }
2149 }
2150 break;
2151
2152 case SOUTH:
2153 if( DEBUG2 ) printf("adjustGridAreaDefinition: new south given\n");
2154 south = RGSELAT(&g);
2155 if( g->northSet ) {
2156 north = RGNWLAT(&g);
2157 if( g->northSouthIncrementSet ) {
2158 increment = RGDJ(&g);
2159 numberOfPoints = (fortint) (0.5 + (north-south)/increment) + 1;
2160 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongMeridian),
2161 &numberOfPoints);
2162 if( DEBUG2 )
2163 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2164 numberOfPoints);
2165 }
2166 else if( g->northSouthNumberOfPointsSet ) {
2167 numberOfPoints = IGNJ(&g) - 1;
2168 increment = (north-south)/(fortdouble) numberOfPoints;
2169 value = (fortint) (increment * 1000.0);
2170 MOVE2BYTES(((g->g2)->grid.latlon.jDirectionIncrement),&value);
2171 g->northSouthIncrementSet = TRUE;
2172 if( DEBUG2 )
2173 printf("adjustGridAreaDefinition: north-south increment set = %f\n",
2174 increment);
2175 }
2176 }
2177 else {
2178 if( g->northSouthIncrementSet && g->northSouthNumberOfPointsSet ) {
2179 increment = RGDJ(&g);
2180 numberOfPoints = IGNJ(&g) - 1;
2181 north = south + increment * (fortdouble) numberOfPoints;
2182 value = (fortint) (north * 1000.0);
2183 if( value < 0 ) value = 0x800000 | (-value);
2184 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfFirstPoint),&value);
2185 g->northSet = TRUE;
2186 if( DEBUG2 )
2187 printf("adjustGridAreaDefinition: north set = %f\n", north);
2188 }
2189 }
2190 break;
2191
2192 case WEST:
2193 if( DEBUG2 ) printf("adjustGridAreaDefinition: new west given\n");
2194 west = RGNWLON(&g);
2195 if( g->eastSet ) {
2196 east = RGSELON(&g);
2197 if( g->westEastIncrementSet ) {
2198 increment = RGDI(&g);
2199 if( ALMOSTZERO(east - west - 360.0) )
2200 numberOfPoints = (fortint) (0.5 + (east-west)/increment);
2201 else
2202 numberOfPoints = (fortint) (0.5 + (east-west)/increment) + 1;
2203 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongParallel),
2204 &numberOfPoints);
2205 g->westEastNumberOfPointsSet = TRUE;
2206 if( DEBUG2 )
2207 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2208 numberOfPoints);
2209 }
2210 else if( g->westEastNumberOfPointsSet ) {
2211 numberOfPoints = IGNI(&g) - 1;
2212 /*
2213 // Test for longitude wraparound
2214 */
2215 if( ALMOSTZERO(east - west - 360.0) ) numberOfPoints++;
2216 increment = (east-west)/(fortdouble) numberOfPoints;
2217 value = (fortint) (increment * 1000.0);
2218 MOVE2BYTES(((g->g2)->grid.latlon.iDirectionIncrement),&value);
2219 g->westEastIncrementSet = TRUE;
2220 if( DEBUG2 )
2221 printf("adjustGridAreaDefinition: west-east increment set = %f\n",
2222 increment);
2223 }
2224 }
2225 else {
2226 if( g->westEastIncrementSet && g->westEastNumberOfPointsSet ) {
2227 increment = RGDI(&g);
2228 numberOfPoints = IGNI(&g) - 1;
2229 east = west + increment * (fortdouble) numberOfPoints;
2230 value = (fortint) (east * 1000.0);
2231 if( value < 0 ) value = 0x800000 | (-value);
2232 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfLastPoint),&value);
2233 g->eastSet = TRUE;
2234 if( DEBUG2 )
2235 printf("adjustGridAreaDefinition: east set = %f\n", east);
2236 }
2237 }
2238 break;
2239
2240 case EAST:
2241 if( DEBUG2 ) printf("adjustGridAreaDefinition: new east given\n");
2242 east = RGSELON(&g);
2243 if( g->westSet ) {
2244 west = RGNWLON(&g);
2245 if( g->westEastIncrementSet ) {
2246 increment = RGDI(&g);
2247 if( ALMOSTZERO(east - west - 360.0) )
2248 numberOfPoints = (fortint) (0.5 + (east-west)/increment);
2249 else
2250 numberOfPoints = (fortint) (0.5 + (east-west)/increment) + 1;
2251 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongParallel),
2252 &numberOfPoints);
2253 g->westEastNumberOfPointsSet = TRUE;
2254 if( DEBUG2 )
2255 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2256 numberOfPoints);
2257 }
2258 else if( g->westEastNumberOfPointsSet ) {
2259 numberOfPoints = IGNI(&g) - 1;
2260 /*
2261 // Test for longitude wraparound
2262 */
2263 if( ALMOSTZERO(east - west - 360.0) ) numberOfPoints++;
2264 increment = (east-west)/(fortdouble) numberOfPoints;
2265 value = (fortint) (increment * 1000.0);
2266 MOVE2BYTES(((g->g2)->grid.latlon.iDirectionIncrement),&value);
2267 g->westEastIncrementSet = TRUE;
2268 if( DEBUG2 )
2269 printf("adjustGridAreaDefinition: west-east increment set = %f\n",
2270 increment);
2271 }
2272 }
2273 else {
2274 if( g->westEastIncrementSet && g->westEastNumberOfPointsSet ) {
2275 increment = RGDI(&g);
2276 numberOfPoints = IGNI(&g);
2277 west = east - increment * (fortdouble) numberOfPoints;
2278 value = (fortint) (west * 1000.0);
2279 if( value < 0 ) value = 0x800000 | (-value);
2280 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfFirstPoint),&value);
2281 g->westSet = TRUE;
2282 if( DEBUG2 )
2283 printf("adjustGridAreaDefinition: west set = %f\n", west);
2284 }
2285 }
2286 break;
2287
2288 case W_E_INC:
2289 if( DEBUG2 ) printf("adjustGridAreaDefinition: new west-east increment given\n");
2290 increment = RGDI(&g);
2291 if( g->westSet ) {
2292 west = RGNWLON(&g);
2293 if( g->eastSet ) {
2294 east = RGSELON(&g);
2295 if( ALMOSTZERO(east - west - 360.0) )
2296 numberOfPoints = (fortint) (0.5 + (east-west)/increment);
2297 else
2298 numberOfPoints = (fortint) (0.5 + (east-west)/increment) + 1;
2299 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongParallel),
2300 &numberOfPoints);
2301 g->westEastNumberOfPointsSet = TRUE;
2302 if( DEBUG2 )
2303 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2304 numberOfPoints);
2305 }
2306 else if( g->westEastNumberOfPointsSet ) {
2307 numberOfPoints = IGNI(&g);
2308 east = west + increment *(fortdouble) numberOfPoints;
2309 value = (fortint) (east * 1000.0);
2310 if( value < 0 ) value = 0x800000 | (-value);
2311 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfLastPoint),&value);
2312 g->eastSet = TRUE;
2313 if( DEBUG2 )
2314 printf("adjustGridAreaDefinition: east set = %f\n", east);
2315 }
2316 }
2317 else {
2318 if( g->eastSet && g->westEastNumberOfPointsSet ) {
2319 east = RGSELON(&g);
2320 numberOfPoints = IGNI(&g);
2321 west = east - increment * (fortdouble) numberOfPoints;
2322 value = (fortint) (west * 1000.0);
2323 if( value < 0 ) value = 0x800000 | (-value);
2324 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfFirstPoint),&value);
2325 g->westSet = TRUE;
2326 if( DEBUG2 )
2327 printf("adjustGridAreaDefinition: west set = %f\n", west);
2328 }
2329 }
2330 break;
2331
2332 case N_S_INC:
2333 if( DEBUG2 ) printf("adjustGridAreaDefinition: new north-south increment given\n");
2334 increment = RGDJ(&g);
2335 if( g->northSet ) {
2336 north = RGNWLAT(&g);
2337 if( g->southSet ) {
2338 south = RGSELAT(&g);
2339 numberOfPoints = (fortint) (0.5 + (north-south)/increment) + 1;
2340 MOVE2BYTES(((g->g2)->grid.latlon.numberOfPointsAlongMeridian),
2341 &numberOfPoints);
2342 if( DEBUG2 )
2343 printf("adjustGridAreaDefinition: numberOfPoints set = %d\n",
2344 numberOfPoints);
2345 }
2346 else if( g->northSouthNumberOfPointsSet ) {
2347 numberOfPoints = IGNJ(&g) - 1;
2348 south = north - increment *(fortdouble) numberOfPoints;
2349 value = (fortint) (south * 1000.0);
2350 if( value < 0 ) value = 0x800000 | (-value);
2351 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfLastPoint),&value);
2352 g->southSet = TRUE;
2353 if( DEBUG2 )
2354 printf("adjustGridAreaDefinition: south set = %f\n", south);
2355 }
2356 }
2357 else {
2358 if( g->southSet && g->northSouthNumberOfPointsSet ) {
2359 south = RGSELAT(&g);
2360 numberOfPoints = IGNJ(&g);
2361 north = south + increment *(fortdouble) numberOfPoints;
2362 value = (fortint) (north * 1000.0);
2363 if( value < 0 ) value = 0x800000 | (-value);
2364 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfFirstPoint),&value);
2365 g->northSet = TRUE;
2366 if( DEBUG2 )
2367 printf("adjustGridAreaDefinition: north set = %f\n", north);
2368 }
2369 }
2370 break;
2371
2372 case W_E_PTS:
2373 if( DEBUG2 ) printf("adjustGridAreaDefinition: new west-east number of points given\n");
2374 numberOfPoints = IGNI(&g);
2375 if( g->westSet ) {
2376 west = RGNWLON(&g);
2377 if( g->eastSet ) {
2378 east = RGSELON(&g);
2379 if( ALMOSTZERO(east - west - 360.0) ) numberOfPoints--;
2380 increment = (east-west)/(fortdouble)numberOfPoints;
2381 value = (fortint) (increment * 1000.0);
2382 MOVE2BYTES(((g->g2)->grid.latlon.iDirectionIncrement),&value);
2383 g->westEastIncrementSet = TRUE;
2384 if( DEBUG2 )
2385 printf("adjustGridAreaDefinition: west-east increment set = %f\n",
2386 increment);
2387 }
2388 else if( g->westEastIncrementSet ) {
2389 increment = RGDI(&g);
2390 east = west + increment *(fortdouble) numberOfPoints;
2391 if( (east - west) > 360.0 ) east -= increment;
2392 value = (fortint) (east * 1000.0);
2393 if( value < 0 ) value = 0x800000 | (-value);
2394 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfLastPoint),&value);
2395 g->eastSet = TRUE;
2396 if( DEBUG2 )
2397 printf("adjustGridAreaDefinition: east set = %f\n", east);
2398 }
2399 }
2400 else {
2401 if( g->eastSet && g->westEastIncrementSet ) {
2402 east = RGSELON(&g);
2403 increment = RGDI(&g);
2404 west = east - increment *(fortdouble) numberOfPoints;
2405 if( (east - west) > 360.0 ) east += increment;
2406 value = (fortint) (west * 1000.0);
2407 if( value < 0 ) value = 0x800000 | (-value);
2408 MOVE3BYTES(((g->g2)->grid.latlon.longitudeOfFirstPoint),&value);
2409 g->westSet = TRUE;
2410 if( DEBUG2 )
2411 printf("adjustGridAreaDefinition: west set = %f\n", west);
2412 }
2413 }
2414 break;
2415
2416 case N_S_PTS:
2417 if( DEBUG2 ) printf("adjustGridAreaDefinition: new north-south number of points given\n");
2418 numberOfPoints = IGNJ(&g) - 1;
2419 if( g->northSet ) {
2420 north = RGNWLAT(&g);
2421 if( g->southSet ) {
2422 south = RGSELAT(&g);
2423 increment = (north-south)/(fortdouble)numberOfPoints;
2424 value = (fortint) (increment * 1000.0);
2425 MOVE2BYTES(((g->g2)->grid.latlon.jDirectionIncrement),&value);
2426 g->northSouthIncrementSet = TRUE;
2427 if( DEBUG2 )
2428 printf("adjustGridAreaDefinition: north-south increment set = %f\n", increment);
2429 }
2430 else if( g->northSouthIncrementSet ) {
2431 increment = RGDJ(&g);
2432 south = north - increment *(fortdouble) numberOfPoints;
2433 value = (fortint) (south * 1000.0);
2434 if( value < 0 ) value = 0x800000 | (-value);
2435 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfLastPoint),&value);
2436 g->southSet = TRUE;
2437 if( DEBUG2 )
2438 printf("adjustGridAreaDefinition: south set = %f\n", south);
2439 }
2440 }
2441 else {
2442 if( g->southSet && g->northSouthIncrementSet ) {
2443 south = RGSELAT(&g);
2444 increment = RGDJ(&g);
2445 north = south + increment *(fortdouble) numberOfPoints;
2446 value = (fortint) (north * 1000.0);
2447 if( value < 0 ) value = 0x800000 | (-value);
2448 MOVE3BYTES(((g->g2)->grid.latlon.latitudeOfFirstPoint),&value);
2449 g->northSet = TRUE;
2450 if( DEBUG2 )
2451 printf("adjustGridAreaDefinition: north set = %f\n", north);
2452 }
2453 }
2454 break;
2455
2456 default:
2457 printf("adjustGridAreaDefinition: not yet implemented\n");
2458 break;
2459
2460 }
2461
2462 g->numberOfValues = g2_ni(g) * g2_nj(g);
2463 return;
2464 }
2465
SPV(gribProduct ** grib,fortdouble * list,fortint * sizeList)2466 fortint SPV(gribProduct** grib, fortdouble* list, fortint* sizeList) {
2467 gribProduct * g = *grib;
2468 fortint section2Size, requiredMemory, listOffset, numberOfRows;
2469 fortint numberOfNewCoordinates = *sizeList;
2470 gribSection2 * newG2;
2471 fortint numberOfOldCoordinates, number, loop;
2472 unsigned char * nextValue, * newLocation, * newCoordinateLocation;
2473 fortint exponent, mantissa, numberOfBitsPerInteger;
2474 fortdouble nextCoordinate;
2475
2476 /*
2477 // Ensure there is enough memory to accomodate the list of vertical coordinates
2478 */
2479 section2Size = g2_length(g);
2480 requiredMemory = section2Size - g2_NV(g)*4 + numberOfNewCoordinates*4;
2481
2482 if( requiredMemory > section2Size ) {
2483 newG2 = (gribSection2 *) allocateMemory(requiredMemory);
2484 memcpy(newG2,(g->g2),section2Size);
2485 freeMemory((g->g2));
2486 (g->g2) = newG2;
2487 }
2488 MOVE3BYTES(((g->g2)->sectionLength),&requiredMemory);
2489
2490 switch( (int) g2_datatype(g) ) {
2491
2492 case 0:
2493 case 4:
2494 case 5:
2495 case 50:
2496 listOffset = 32;
2497 break;
2498
2499 case 1:
2500 case 3:
2501 case 6:
2502 case 8:
2503 case 10:
2504 case 13:
2505 case 14:
2506 case 20:
2507 case 24:
2508 case 60:
2509 case 70:
2510 listOffset = 42;
2511 break;
2512
2513 case 30:
2514 case 34:
2515 case 80:
2516 listOffset = 52;
2517 break;
2518
2519 case 90:
2520 listOffset = 44;
2521 break;
2522
2523 default:
2524 printf("SPV: does not handle representation type = %d\n",
2525 g2_datatype(g));
2526 exit(1);
2527 }
2528 /*
2529 // Move the list of the number of points at each latitude if it is a
2530 // reduced grid.
2531 */
2532 numberOfRows = g2_nj(g);
2533
2534 if( !directionIncrementsGiven(g) && (numberOfRows != 0) ) {
2535
2536 numberOfOldCoordinates = g2_NV(g);
2537
2538 if( numberOfOldCoordinates < numberOfNewCoordinates ) {
2539 nextValue = numberOfRows*2 +
2540 (unsigned char*)(g->g2) + g2_PV_PL(g) + numberOfOldCoordinates*4 - 1;
2541 newLocation = numberOfRows*2 +
2542 (unsigned char*)(g->g2) + listOffset + numberOfNewCoordinates*4;
2543 for( loop = numberOfRows; loop >= 0; loop-- ) {
2544 number = TWOBYTEINT(nextValue);
2545 nextValue -= 2;
2546 MOVE2BYTES(newLocation,&number);
2547 newLocation -= 2;
2548 }
2549 }
2550 else {
2551 {
2552 fortint temp;
2553 temp = g2_PV_PL(g);
2554 printf("********************* g2_PV_PL(g) = %d\n",temp);
2555 }
2556 nextValue =
2557 (unsigned char*)(g->g2) + g2_PV_PL(g) + numberOfOldCoordinates*4 - 1;
2558 newLocation =
2559 (unsigned char*)(g->g2) + listOffset + numberOfNewCoordinates*4;
2560 for( loop = 0; loop < numberOfRows; loop++ ) {
2561 number = TWOBYTEINT(nextValue);
2562 nextValue += 2;
2563 MOVE2BYTES(newLocation,&number);
2564 newLocation += 2;
2565 }
2566 }
2567 }
2568 /*
2569 // Add the vertical coordinates
2570 */
2571 newCoordinateLocation = (unsigned char*)(g->g2) + listOffset;
2572 numberOfBitsPerInteger = sizeof(fortint) * 8;
2573
2574 for( loop = 0; loop < numberOfNewCoordinates; loop++) {
2575 nextCoordinate = list[loop];
2576 if( REF2GRB(&nextCoordinate,&exponent,&mantissa,&numberOfBitsPerInteger) ) {
2577 printf("SPV: call to REF2GRB failed\n");
2578 exit(1);
2579 }
2580 MOVE1BYTE(newCoordinateLocation,&exponent);
2581 MOVE3BYTES((newCoordinateLocation+1),&mantissa);
2582 newCoordinateLocation += 4;
2583 }
2584 MOVE1BYTE(((g->g2)->NV),&numberOfNewCoordinates);
2585 listOffset++;
2586 MOVE1BYTE(((g->g2)->PV_PL),&listOffset);
2587 return 0;
2588 }
2589