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 <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <ctype.h>
16
17
18 #include "PBGroutines.h"
19
20 #ifdef gribAPI
21 long date_to_Julian(long ddate);
22 long Julian_to_date(long jdate);
23 #else
24 long date_to_julian(long ddate);
25 long julian_to_date(long jdate);
26 #endif
27
28 fortint pbginitInput(char* filename, fortint filename_len);
29 fortint pbginitOutput(char* filename, fortint filename_len);
30 void dsgnbt_(fortint *, fortint *, fortint *, fortint *);
31 void copyName(char* * , char* , fortint );
32
33 fortint readgrib(FILE * file, unsigned char * buffer, fortint * prod_len);
34
35 void pbgindx(fortint);
36
37 void gribdata( unsigned char * header,
38 fortint number_of_fields,
39 fortint * parameter,
40 fortint * level,
41 fortint * date,
42 fortint * time,
43 fortint * timestep,
44 fortint * localDefinitionNumber,
45 fortint * type,
46 fortint * stream,
47 fortint * repres,
48 fortint * levtype,
49 fortint * number,
50 fortint * vdate,
51 fortint * vtime,
52 fortint * tcNumber,
53 fortint * tcTotal,
54 fortint * clusterMethod,
55 fortint * tcStep,
56 fortint * clusterStepEnd,
57 fortint * tcNorth,
58 fortint * tcWest,
59 fortint * tcSouth,
60 fortint * tcEast,
61 fortint * clusterOpFcNumber,
62 fortint * clusterControlFcNumber,
63 fortint * tcNumOfFcs,
64 fortint * probScale,
65 fortint * probIndicator,
66 fortint * probLower,
67 fortint * probUpper,
68 fortint * ni,
69 fortint * nj,
70 fortint * nlat,
71 fortint * nlon,
72 fortint * slat,
73 fortint * slon,
74 fortint * di,
75 fortint * dj,
76 fortint * splat,
77 fortint * splon,
78 fortint * quasi,
79 fortint * directionNumber,
80 fortint * frequencyNumber,
81 fortint * optimisationTime,
82 fortint * leadTime,
83 fortint * sensitiveAreaDomain,
84 fortint * sensitiveAreaMethod,
85 fortint * verifyingMonth,
86 fortint * averagingPeriod,
87 fortint * forecastMonth,
88 fortint * referenceDate,
89 fortint * climateDateFrom,
90 fortint * climateDateTo,
91 fortint * thresholdUnitsScale,
92 fortint * thresholdIndicator,
93 fortint * lowerThreshold,
94 fortint * upperThreshold
95 );
96
97 static gribfile * latestFile(collection);
98 static gribfile * currentFile(collection, fortint);
99
100 /*
101 //
102 */
103 #define DEBUGOFF 1
104 #define DEBUG1 (debugSet > DEBUGOFF )
105 #define DEBUG2 (debugSet > (DEBUGOFF + 1) )
106 static char * debugLevel;
107 static int debugSet = 0;
108
109 #define ONEBYTELONG(a) (long) ( *(a) )
110 #define TWOBYTELONG(a) (long) ( (*(a))<<8 | (*((a)+1))<<0 )
111 #define THREEBYTELONG(a) (long) (TWOBYTELONG((a))<<8 | (*((a)+2))<<0 )
112 #define FOURBYTELONG(a) (long) (THREEBYTELONG((a))<<8 | (*((a)+3))<<0 )
113
114 #define MISSING -1
115 #define NOTGIVEN(x) ((x) == MISSING )
116 #define MATCH(a,b) (NOTGIVEN((a)) || ((a) == (b)))
117
118 #define HEADLEN 10000
119
120 fortint exists(char* , fortint );
121 fortint addRead(char* , fortint );
122 fortint addWrite(char* , fortint );
123 void removeFile(char* , fortint );
124
125 collection openFiles = {exists,addRead,addWrite,removeFile,-1,-1,NULL};
126
pbgindx(fortint thisFile)127 void pbgindx(fortint thisFile) {
128 /*
129 // Fills in details of current file
130 */
131 unsigned char header[HEADLEN];
132 fortint status;
133 fortint headerlen;
134 off_t space;
135 fortint number_of_fields = 0;
136 gribfile * file = currentFile(openFiles,thisFile);
137
138 /*
139 // Loop through products in the file
140 */
141 file->offset[0] = 0;
142 do {
143 headerlen = HEADLEN;
144 status = readgrib( file->fp, header, &headerlen);
145
146 /*
147 // Accept product if OK or if only problem is 'buffer too small' or
148 // on EOF
149 */
150 if( (status!=0) && (status!=-3) && (status!=-1) ) exit(1);
151 if( status == -1 ) break;
152
153 /*
154 // Note the product byte offset in the file
155 */
156 file->length[number_of_fields] = (off_t) headerlen;
157 file->offset[number_of_fields] =
158 ftell(file->fp) - file->length[number_of_fields];
159
160
161 /*
162 // If it's a GRIB edition 1 product, pick up more information
163 // (parameter,level,..)
164 */
165 if( strncmp((char *)header,"GRIB", 4)==0 ) {
166 /*
167 // Eliminate GRIB editions 0 and -1
168 */
169 if( (header[ 7] != '\1') ) break;
170 if( (header[19] == '\0') && (header[20] == '\0') &&
171 (header[21] == '\0') && (header[22] == '\0') &&
172 (header[23] == '\0') ) break;
173
174 gribdata( header,
175 number_of_fields,
176 file->parameter,
177 file->level,
178 file->date,
179 file->time,
180 file->timestep,
181 file->localDefinitionNumber,
182 file->type,
183 file->stream,
184 file->repres,
185 file->levtype,
186 file->number,
187 file->vdate,
188 file->vtime,
189 file->tcNumber,
190 file->tcTotal,
191 file->clusterMethod,
192 file->tcStep,
193 file->clusterStepEnd,
194 file->tcNorth,
195 file->tcWest,
196 file->tcSouth,
197 file->tcEast,
198 file->clusterOpFcNumber,
199 file->clusterControlFcNumber,
200 file->tcNumOfFcs,
201 file->probScale,
202 file->probIndicator,
203 file->probLower,
204 file->probUpper,
205 file->ni,
206 file->nj,
207 file->nlat,
208 file->nlon,
209 file->slat,
210 file->slon,
211 file->di,
212 file->dj,
213 file->splat,
214 file->splon,
215 file->quasi,
216 file->directionNumber,
217 file->frequencyNumber,
218 file->optimisationTime,
219 file->leadTime,
220 file->sensitiveAreaDomain,
221 file->sensitiveAreaMethod,
222 file->verifyingMonth,
223 file->averagingPeriod,
224 file->forecastMonth,
225 file->referenceDate,
226 file->climateDateFrom,
227 file->climateDateTo,
228 file->thresholdUnitsScale,
229 file->thresholdIndicator,
230 file->lowerThreshold,
231 file->upperThreshold
232 );
233
234 #define REALLOC(a) (a) = realloc((a) , (size_t) space)
235
236 number_of_fields++;
237 if( number_of_fields == file->max ) {
238 fortint newSize = (file->max)*2;
239 file->max = newSize;
240
241 space = (off_t) (sizeof(off_t)*newSize);
242 REALLOC(file->offset);
243
244 space = (off_t) (sizeof(fortint)*newSize);
245 REALLOC(file->length);
246 REALLOC(file->parameter);
247 REALLOC(file->level);
248 REALLOC(file->date);
249 REALLOC(file->time);
250 REALLOC(file->timestep);
251 REALLOC(file->localDefinitionNumber);
252 REALLOC(file->type);
253 REALLOC(file->stream);
254 REALLOC(file->repres);
255 REALLOC(file->levtype);
256 REALLOC(file->number);
257 REALLOC(file->vdate);
258 REALLOC(file->vtime);
259 REALLOC(file->tcNumber);
260 REALLOC(file->tcTotal);
261 REALLOC(file->clusterMethod);
262 REALLOC(file->tcStep);
263 REALLOC(file->clusterStepEnd);
264 REALLOC(file->tcNorth);
265 REALLOC(file->tcWest);
266 REALLOC(file->tcSouth);
267 REALLOC(file->tcEast);
268 REALLOC(file->clusterOpFcNumber);
269 REALLOC(file->clusterControlFcNumber);
270 REALLOC(file->tcNumOfFcs);
271 REALLOC(file->probScale);
272 REALLOC(file->probIndicator);
273 REALLOC(file->probLower);
274 REALLOC(file->probUpper);
275 REALLOC(file->ni);
276 REALLOC(file->nj);
277 REALLOC(file->nlat);
278 REALLOC(file->nlon);
279 REALLOC(file->slat);
280 REALLOC(file->slon);
281 REALLOC(file->di);
282 REALLOC(file->dj);
283 REALLOC(file->splat);
284 REALLOC(file->splon);
285 REALLOC(file->quasi);
286 REALLOC(file->directionNumber);
287 REALLOC(file->frequencyNumber);
288 REALLOC(file->optimisationTime);
289 REALLOC(file->leadTime);
290 REALLOC(file->sensitiveAreaDomain);
291 REALLOC(file->sensitiveAreaMethod);
292 REALLOC(file->verifyingMonth);
293 REALLOC(file->averagingPeriod);
294 REALLOC(file->forecastMonth);
295 REALLOC(file->referenceDate);
296 REALLOC(file->climateDateFrom);
297 REALLOC(file->climateDateTo);
298 REALLOC(file->thresholdUnitsScale);
299 REALLOC(file->thresholdIndicator);
300 REALLOC(file->lowerThreshold);
301 REALLOC(file->upperThreshold);
302 }
303 }
304
305 } while( (!feof(file->fp)) &&
306 (number_of_fields < file->max)
307 );
308
309 file->count = number_of_fields;
310
311 return;
312 }
313
gribdata(unsigned char * header,fortint index,fortint * parameter,fortint * level,fortint * date,fortint * time,fortint * timestep,fortint * localDefinitionNumber,fortint * type,fortint * stream,fortint * repres,fortint * levtype,fortint * number,fortint * vdate,fortint * vtime,fortint * tcNumber,fortint * tcTotal,fortint * clusterMethod,fortint * tcStep,fortint * clusterStepEnd,fortint * tcNorth,fortint * tcWest,fortint * tcSouth,fortint * tcEast,fortint * clusterOpFcNumber,fortint * clusterControlFcNumber,fortint * tcNumOfFcs,fortint * probScale,fortint * probIndicator,fortint * probLower,fortint * probUpper,fortint * ni,fortint * nj,fortint * nlat,fortint * nlon,fortint * slat,fortint * slon,fortint * di,fortint * dj,fortint * splat,fortint * splon,fortint * quasi,fortint * directionNumber,fortint * frequencyNumber,fortint * optimisationTime,fortint * leadTime,fortint * sensitiveAreaDomain,fortint * sensitiveAreaMethod,fortint * verifyingMonth,fortint * averagingPeriod,fortint * forecastMonth,fortint * referenceDate,fortint * climateDateFrom,fortint * climateDateTo,fortint * thresholdUnitsScale,fortint * thresholdIndicator,fortint * lowerThreshold,fortint * upperThreshold)314 void gribdata( unsigned char * header,
315 fortint index,
316 fortint * parameter,
317 fortint * level,
318 fortint * date,
319 fortint * time,
320 fortint * timestep,
321 fortint * localDefinitionNumber,
322 fortint * type,
323 fortint * stream,
324 fortint * repres,
325 fortint * levtype,
326 fortint * number,
327 fortint * vdate,
328 fortint * vtime,
329 fortint * tcNumber,
330 fortint * tcTotal,
331 fortint * clusterMethod,
332 fortint * tcStep,
333 fortint * clusterStepEnd,
334 fortint * tcNorth,
335 fortint * tcWest,
336 fortint * tcSouth,
337 fortint * tcEast,
338 fortint * clusterOpFcNumber,
339 fortint * clusterControlFcNumber,
340 fortint * tcNumOfFcs,
341 fortint * probScale,
342 fortint * probIndicator,
343 fortint * probLower,
344 fortint * probUpper,
345 fortint * ni,
346 fortint * nj,
347 fortint * nlat,
348 fortint * nlon,
349 fortint * slat,
350 fortint * slon,
351 fortint * di,
352 fortint * dj,
353 fortint * splat,
354 fortint * splon,
355 fortint * quasi,
356 fortint * directionNumber,
357 fortint * frequencyNumber,
358 fortint * optimisationTime,
359 fortint * leadTime,
360 fortint * sensitiveAreaDomain,
361 fortint * sensitiveAreaMethod,
362 fortint * verifyingMonth,
363 fortint * averagingPeriod,
364 fortint * forecastMonth,
365 fortint * referenceDate,
366 fortint * climateDateFrom,
367 fortint * climateDateTo,
368 fortint * thresholdUnitsScale,
369 fortint * thresholdIndicator,
370 fortint * lowerThreshold,
371 fortint * upperThreshold
372 )
373 /*
374 // Returns parameter, level, date, time , etc information
375 // for a GRIB product.
376 */
377 {
378
379 /*
380 // Decode the GRIB product
381 */
382
383 #define POLAR_STEREO 5
384
385 #define SEC1 7
386 #define LEN1 8
387 #define TABLE 4
388 #define PARAM 9
389 #define LEVELTYPE 10
390 #define LEVEL1 11
391 #define LEVEL2 12
392 #define YY 13
393 #define MM 14
394 #define DD 15
395 #define HOUR 16
396 #define MIN 17
397 #define TIMEUNIT 18
398 #define STEP1 19
399 #define STEP2 20
400 #define TIMERANGE 21
401 #define CENTURY 25
402 #define LOCALDEF 41
403 #define CLASS 42
404 #define TYPE 43
405 #define STREAM 44
406 #define NUMBER 50
407 #define REPRES 6
408 #define NI 7
409 #define NJ 9
410 #define NLAT 11
411 #define NLON 14
412 #define SLAT 18
413 #define SLON 21
414 #define DI 24
415 #define DJ 26
416 #define DI_PS 21
417 #define DJ_PS 24
418 #define SPLAT 33
419 #define SPLON 36
420
421 #define ISOTHERMAL_LEVEL 20
422 #define ISOBARIC_LEVEL 100
423 #define ALTITUDE 103
424 #define HEIGHT 105
425 #define SIGMA_LEVEL 107
426 #define HYBRID_LEVEL 109
427 #define DEPTH_BELOW_LAND 111
428 #define ISENTROPIC_LEVEL 113
429 #define PRESSURE_DIFF 115
430 #define VORTICITY_SURF 117
431 #define HEIGHT_HIGH_PREC 125
432 #define SATELLITE_BAND 127
433 #define DEPTH_BELOW_SEA 160
434 #define ENTIRE_ATMOSPHERE 200
435 #define ENTIRE_OCEAN 201
436 #define ISOBARIC_HIGH_PREC 210
437
438 #define PROBABILITY 16
439 #define PROBSCALE 52
440 #define PROBINDIC 53
441 #define PROBLOWER 54
442 #define PROBUPPER 56
443
444 #define CLUSTERMEAN 14
445 #define CLUSTERSD 15
446 #define CLUSTERNUMBER 50
447 #define CLUSTERTOTAL 51
448 #define CLUSTERMETHOD 53
449 #define CLUSTERSTEP 54
450 #define CLUSTERSTEPEND 56
451 #define CLUSTERNORTH 58
452 #define CLUSTERWEST 61
453 #define CLUSTERSOUTH 64
454 #define CLUSTEREAST 67
455 #define CLUSTEROPFCNUM 70
456 #define CLUSTERCONTROLFCNUM 71
457 #define CLUSTERFORECASTS 72
458
459 #define TUBETYPE 24
460 #define TUBENUMBER 50
461 #define TUBETOTAL 51
462 #define TUBESTEP 71
463 #define TUBENORTH 55
464 #define TUBEWEST 58
465 #define TUBESOUTH 61
466 #define TUBEEAST 64
467 #define TUBEFORECASTS 79
468
469 #define CLUSTERMEAN 14
470 #define CLUSTERSD 15
471 #define CLUSTERNUMBER 50
472 #define CLUSTERTOTAL 51
473 #define CLUSTERMETHOD 53
474 #define CLUSTERSTEP 54
475 #define CLUSTERSTEPEND 56
476 #define CLUSTERNORTH 58
477 #define CLUSTERWEST 61
478 #define CLUSTERSOUTH 64
479 #define CLUSTEREAST 67
480 #define CLUSTEROPFCNUM 70
481 #define CLUSTERCONTROLFCNUM 71
482 #define CLUSTERFORECASTS 72
483
484 #define DIRECTIONNUMBER 52
485 #define FREQUENCYNUMBER 53
486
487 #define OPTIMISATIONTIME 92
488 #define LEADTIME 93
489 #define SENSITIVEAREADOMAIN 94
490 #define SENSITIVEAREAMETHOD 95
491
492 #define VERIFYINGMONTH 56
493 #define AVERAGINGPERIOD 60
494 #define FORECASTMONTH 61
495 #define REFERENCEDATE 63
496 #define CLIMATEDATEFROM 67
497 #define CLIMATEDATETO 71
498 #define THRESHOLDUNITSSCALE 75
499 #define THRESHOLDINDICATOR 76
500 #define LOWERTHRESHOLD 77
501 #define UPPERTHRESHOLD 79
502
503 fortint table, leveltype, timeunit, timerange, multiplier, century, year;
504 fortint twoByteNumber, ECclass, ensembleNumber;
505 fortint section2Offset;
506 fortint status, eight = 8, sixteen = 16, twentyFour = 24;
507
508 table = header[SEC1+TABLE];
509 parameter[index] = table*1000 + header[SEC1+PARAM];
510
511 leveltype = header[SEC1+LEVELTYPE];
512 levtype[index] = leveltype;
513
514 switch( leveltype ) {
515
516 case ISOTHERMAL_LEVEL :
517 case ISOBARIC_LEVEL :
518 case ALTITUDE :
519 case HEIGHT :
520 case SIGMA_LEVEL :
521 case HYBRID_LEVEL :
522 case DEPTH_BELOW_LAND :
523 case ISENTROPIC_LEVEL :
524 case PRESSURE_DIFF :
525 case VORTICITY_SURF :
526 case HEIGHT_HIGH_PREC :
527 case SATELLITE_BAND :
528 case DEPTH_BELOW_SEA :
529 case ENTIRE_ATMOSPHERE :
530 case ENTIRE_OCEAN :
531 case ISOBARIC_HIGH_PREC:
532 /* level[index] = header[SEC1+LEVEL1]*256 + header[SEC1+LEVEL2]; */
533 level[index] = TWOBYTELONG(&header[SEC1+LEVEL1]);
534 break;
535
536 default: level[index] = header[SEC1+LEVEL1];
537 }
538
539 century = header[SEC1+CENTURY];
540 year = (century-1)*100 + header[SEC1+YY];
541 date[index] = year*10000 + header[SEC1+MM]*100 + header[SEC1+DD];
542
543 time[index] = header[SEC1+HOUR]*100 + header[SEC1+MIN];
544
545 timeunit = header[SEC1+TIMEUNIT];
546
547 #define DAY 2
548 #define MONTH 3
549 #define THREE_HOURS 10
550 #define SIX_HOURS 11
551 #define TWELVE_HOURS 12
552
553 switch( timeunit ) {
554
555 case DAY :
556 multiplier = 24;
557 break;
558
559 case MONTH :
560 multiplier = 30*24;
561 break;
562
563 case THREE_HOURS :
564 multiplier = 3;
565 break;
566
567 case SIX_HOURS :
568 multiplier = 6;
569 break;
570
571 case TWELVE_HOURS :
572 multiplier = 12;
573 break;
574
575 default: multiplier = 1;
576 }
577
578 #define RANGE_P1_TO_P2 2
579 #define ACCUMULATION 4
580 #define DIFFERENCE 5
581 #define TWOBYTE 10
582
583 timerange = header[SEC1+TIMERANGE];
584
585 switch( timerange ) {
586
587 case TWOBYTE :
588 timestep[index] = TWOBYTELONG(&header[SEC1+STEP1])*multiplier;
589 break;
590
591 default: timestep[index] = header[SEC1+STEP2]*multiplier*10000 +
592 header[SEC1+STEP1]*multiplier;
593 }
594
595 type[index] = header[SEC1+TYPE];
596
597 #define ANALYSIS 2
598
599 if( type[index] == ANALYSIS ) {
600 vdate[index] = date[index];
601 vtime[index] = time[index];
602 }
603 else {
604 long iverd = (time[index] + timestep[index]*100) / 2400;
605 #ifdef gribAPI
606 long ijvday = date_to_Julian(date[index]) + iverd;
607 vdate[index] = Julian_to_date(ijvday);
608 #else
609 long ijvday = date_to_julian(date[index]) + iverd;
610 vdate[index] = julian_to_date(ijvday);
611 #endif
612 vtime[index] = (time[index] + timestep[index]*100) - (iverd*2400);
613 }
614
615 localDefinitionNumber[index] = header[SEC1+LOCALDEF];
616 ECclass = header[SEC1+CLASS];
617 stream[index] = TWOBYTELONG(&header[SEC1+STREAM]);
618
619 twoByteNumber = ( (localDefinitionNumber[index] == 4) && (stream[index] == 1090) ) ||
620 (localDefinitionNumber[index] == 9) ||
621 (localDefinitionNumber[index] == 15) ||
622 (localDefinitionNumber[index] == 16) ||
623 (localDefinitionNumber[index] == 21) ||
624 (localDefinitionNumber[index] == 23);
625
626 if( twoByteNumber )
627 ensembleNumber = TWOBYTELONG(&header[SEC1+NUMBER]);
628 else
629 ensembleNumber = header[SEC1+NUMBER];
630 number[index] = ensembleNumber;
631
632 /*
633 // Extra information for tubes
634 */
635 if( type[index] == TUBETYPE ) {
636 tcNumber[index] = header[SEC1+TUBENUMBER];
637 tcTotal[index] = header[SEC1+TUBETOTAL];
638 tcStep[index] = TWOBYTELONG(&header[SEC1+TUBESTEP]);
639 tcNorth[index] = THREEBYTELONG(&header[SEC1+TUBENORTH]);
640 dsgnbt_(&tcNorth[index],&tcNorth[index],&twentyFour,&status);
641 tcWest[index] = THREEBYTELONG(&header[SEC1+TUBEWEST]);
642 dsgnbt_(&tcWest[index],&tcWest[index],&twentyFour,&status);
643 tcSouth[index] = THREEBYTELONG(&header[SEC1+TUBESOUTH]);
644 dsgnbt_(&tcSouth[index],&tcSouth[index],&twentyFour,&status);
645 tcEast[index] = THREEBYTELONG(&header[SEC1+TUBEEAST]);
646 dsgnbt_(&tcEast[index],&tcEast[index],&twentyFour,&status);
647 tcNumOfFcs[index] = header[SEC1+TUBEFORECASTS];
648 }
649
650 /*
651 // Extra information for clusters
652 */
653 if( (type[index] == CLUSTERMEAN) || (type[index] == CLUSTERSD) ) {
654 tcNumber[index] = header[SEC1+CLUSTERNUMBER];
655 tcTotal[index] = header[SEC1+CLUSTERTOTAL];
656 clusterMethod[index] = header[SEC1+CLUSTERMETHOD];
657 tcStep[index] = TWOBYTELONG(&header[SEC1+CLUSTERSTEP]);
658 clusterStepEnd[index] = TWOBYTELONG(&header[SEC1+CLUSTERSTEPEND]);
659 tcNorth[index] = THREEBYTELONG(&header[SEC1+CLUSTERNORTH]);
660 dsgnbt_(&tcNorth[index],&tcNorth[index],&twentyFour,&status);
661 tcWest[index] = THREEBYTELONG(&header[SEC1+CLUSTERWEST]);
662 dsgnbt_(&tcWest[index],&tcWest[index],&twentyFour,&status);
663 tcSouth[index] = THREEBYTELONG(&header[SEC1+CLUSTERSOUTH]);
664 dsgnbt_(&tcSouth[index],&tcSouth[index],&twentyFour,&status);
665 tcEast[index] = THREEBYTELONG(&header[SEC1+CLUSTEREAST]);
666 dsgnbt_(&tcEast[index],&tcEast[index],&twentyFour,&status);
667 clusterOpFcNumber[index] = header[SEC1+CLUSTEROPFCNUM];
668 clusterControlFcNumber[index] = header[SEC1+CLUSTERCONTROLFCNUM];
669 tcNumOfFcs[index] = header[SEC1+CLUSTERFORECASTS];
670 }
671 /*
672 // Extra information for probabilities
673 */
674 if( type[index] == PROBABILITY ) {
675 probScale[index] = header[SEC1+PROBSCALE];
676 dsgnbt_(&probScale[index],&probScale[index],&eight,&status);
677 probIndicator[index] = header[SEC1+PROBINDIC];
678 probLower[index] = 65535;
679 if( probIndicator[index] != 2 ) {
680 probLower[index] = TWOBYTELONG(&header[SEC1+PROBLOWER]);
681 dsgnbt_(&probLower[index],&probLower[index],&sixteen,&status);
682 }
683 probUpper[index] = 65535;
684 if( probIndicator[index] != 1 ) {
685 probUpper[index] = TWOBYTELONG(&header[SEC1+PROBUPPER]);
686 dsgnbt_(&probUpper[index],&probUpper[index],&sixteen,&status);
687 }
688 }
689 /*
690 // Extra information for wave 2D spectra
691 */
692 if( localDefinitionNumber[index] == 13 ) {
693 directionNumber[index] = header[SEC1+DIRECTIONNUMBER];
694 frequencyNumber[index] = header[SEC1+FREQUENCYNUMBER];
695 }
696 /*
697 // Extra information for sensitive area forecasts
698 */
699 if( localDefinitionNumber[index] == 21 ) {
700 optimisationTime[index] = header[SEC1+OPTIMISATIONTIME];
701 leadTime[index] = header[SEC1+LEADTIME];
702 sensitiveAreaDomain[index] = header[SEC1+SENSITIVEAREADOMAIN];
703 sensitiveAreaMethod[index] =
704 TWOBYTELONG(&header[SEC1+SENSITIVEAREAMETHOD]);
705 }
706 /*
707 // Extra information for coupled atmospheric, wave and ocean models
708 */
709 if( localDefinitionNumber[index] == 23 ) {
710 unsigned char scale;
711
712 verifyingMonth[index] = FOURBYTELONG(&header[SEC1+VERIFYINGMONTH]);
713 averagingPeriod[index] = ONEBYTELONG(&header[SEC1+AVERAGINGPERIOD]);
714 forecastMonth[index] = TWOBYTELONG(&header[SEC1+FORECASTMONTH]);
715 referenceDate[index] = FOURBYTELONG(&header[SEC1+REFERENCEDATE]);
716 climateDateFrom[index] = FOURBYTELONG(&header[SEC1+CLIMATEDATEFROM]);
717 climateDateTo[index] = FOURBYTELONG(&header[SEC1+CLIMATEDATETO]);
718 scale = header[SEC1+THRESHOLDUNITSSCALE];
719 if( scale & 0x80 ) {
720 scale &= 0x7f;
721 thresholdUnitsScale[index] = - ONEBYTELONG(&scale);
722 }
723 else
724 thresholdUnitsScale[index] = ONEBYTELONG(&scale);
725 thresholdIndicator[index] = ONEBYTELONG(&header[SEC1+THRESHOLDINDICATOR]);
726 lowerThreshold[index] = TWOBYTELONG(&header[SEC1+LOWERTHRESHOLD]);
727 upperThreshold[index] = TWOBYTELONG(&header[SEC1+UPPERTHRESHOLD]);
728 }
729
730 /*
731 // Section 2
732 */
733 section2Offset = THREEBYTELONG(&header[LEN1]);
734 repres[index] = header[LEN1+section2Offset-1+REPRES];
735
736 ni[index] = TWOBYTELONG(&header[LEN1+section2Offset-1+NI]);
737 nj[index] = TWOBYTELONG(&header[LEN1+section2Offset-1+NJ]);
738
739 nlat[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+NLAT]);
740 dsgnbt_(&nlat[index],&nlat[index],&twentyFour,&status);
741
742 nlon[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+NLON]);
743 dsgnbt_(&nlon[index],&nlon[index],&twentyFour,&status);
744
745 slat[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+SLAT]);
746 dsgnbt_(&slat[index],&slat[index],&twentyFour,&status);
747
748 /*
749 // Section 2 is different for polar stereographic
750 */
751 if( repres[index] == POLAR_STEREO ) {
752 di[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+DI_PS]);
753 dj[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+DJ_PS]);
754 }
755 else {
756 slon[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+SLON]);
757 dsgnbt_(&slon[index],&slon[index],&twentyFour,&status);
758
759 di[index] = TWOBYTELONG(&header[LEN1+section2Offset-1+DI]);
760 dj[index] = TWOBYTELONG(&header[LEN1+section2Offset-1+DJ]);
761
762 splat[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+SPLAT]);
763 dsgnbt_(&splat[index],&splat[index],&twentyFour,&status);
764
765 splon[index] = THREEBYTELONG(&header[LEN1+section2Offset-1+SPLON]);
766 dsgnbt_(&splon[index],&splon[index],&twentyFour,&status);
767 }
768
769 quasi[index] = ( (ni[index] == 0xFFFF) || (nj[index] == 0xFFFF) );
770
771 return;
772
773 }
774
775 /*
776 // ****************************************************************
777 */
pbgtotl_(char * filename,fortint filename_len)778 fortint pbgtotl_(char* filename, fortint filename_len) {
779 /*
780 // Returns number of GRIB products in the file.
781 */
782 fortint thisFile;
783 gribfile * file;
784
785 thisFile = pbginitInput(filename,filename_len);
786 file = currentFile(openFiles, thisFile);
787
788 if( DEBUG1 ) {
789 char* pfile;
790 copyName(&pfile, filename, filename_len);
791 printf("PBGTOTL: Number of GRIBs in file %s = %ld\n", pfile, file->count);
792 free(pfile);
793 }
794 return ( file->count );
795 }
796
pbgtotl(char * filename,fortint filename_len)797 fortint pbgtotl(char* filename, fortint filename_len) {
798 return pbgtotl_(filename,filename_len);
799 }
800
801 /*
802 // ****************************************************************
803 */
pbgleng_(char * filename,fortint * n,fortint filename_len)804 fortint pbgleng_(char* filename, fortint * n, fortint filename_len) {
805 /*
806 // Returns length (in bytes) of GRIB product n or -1 if n is greater
807 // than the number of products in the file, or less than 1.
808 */
809 fortint index = (*n)-1;
810 fortint thisFile;
811 gribfile * file;
812
813 if( index >= 0 ) {
814 thisFile = pbginitInput(filename,filename_len);
815 file = currentFile(openFiles, thisFile);
816 if( index < file->count ) {
817 if( DEBUG1 ) {
818 char* pfile;
819 copyName(&pfile, filename, filename_len);
820 printf("PBGLENG: length of GRIB %ld in file %s = %ld\n",
821 (index+1), pfile, file->length[index]);
822 free(pfile);
823 }
824 return ( file->length[index] );
825 }
826 }
827
828 return (-1);
829 }
830
pbgleng(char * filename,fortint * n,fortint filename_len)831 fortint pbgleng(char* filename, fortint * n, fortint filename_len) {
832 return pbgleng_(filename,n,filename_len);
833 }
834
835 /*
836 // ****************************************************************
837 */
pbgoffs_(char * filename,fortint * n,fortint filename_len)838 fortint pbgoffs_(char* filename, fortint * n, fortint filename_len) {
839 /*
840 // Returns offset (in bytes) of GRIB product n or -1 if n is greater
841 // than the number of products in the file, or less than 1.
842 */
843 fortint index = (*n)-1;
844 fortint thisFile;
845 gribfile * file;
846
847 if( index >= 0 ) {
848 thisFile = pbginitInput(filename,filename_len);
849 file = currentFile(openFiles, thisFile);
850 if( index < file->count ) {
851 if( DEBUG1 ) {
852 char* pfile;
853 copyName(&pfile, filename, filename_len);
854 printf("PBGOFFS: offset of GRIB %ld in file %s = %lld\n",
855 (index+1), pfile, file->offset[index]);
856 free(pfile);
857 }
858 return ( file->offset[index] );
859 }
860 }
861
862 return (-1);
863 }
864
pbgoffs(char * filename,fortint * n,fortint filename_len)865 fortint pbgoffs(char* filename, fortint * n, fortint filename_len) {
866 return pbgoffs_(filename,n,filename_len);
867 }
868
869 /*
870 // ****************************************************************
871 */
pbgparm_(char * filename,fortint * n,fortint filename_len)872 fortint pbgparm_(char* filename, fortint * n, fortint filename_len) {
873 /*
874 // Returns parameter number of GRIB product n or -1 if n is greater
875 // than the number of products in the file, or less than 1.
876 */
877 fortint index = (*n)-1;
878 fortint thisFile;
879 gribfile * file;
880
881 if( index >= 0 ) {
882 thisFile = pbginitInput(filename,filename_len);
883 file = currentFile(openFiles, thisFile);
884 if( index < file->count ) {
885 if( DEBUG1 ) {
886 char* pfile;
887 copyName(&pfile, filename, filename_len);
888 printf("PBGPARM: parameter number of GRIB %ld in file %s = %ld\n",
889 (index+1), pfile, file->parameter[index]);
890 free(pfile);
891 }
892 return ( file->parameter[index] );
893 }
894 }
895
896 return (-1);
897 }
898
pbgparm(char * filename,fortint * n,fortint filename_len)899 fortint pbgparm(char* filename, fortint * n, fortint filename_len) {
900 return pbgparm_(filename,n,filename_len);
901 }
902
903 /*
904 // ****************************************************************
905 */
pbglevl_(char * filename,fortint * n,fortint filename_len)906 fortint pbglevl_(char* filename, fortint * n, fortint filename_len) {
907 /*
908 // Returns level of GRIB product n or -1 if n is greater
909 // than the number of products in the file, or less than 1.
910 */
911 fortint index = (*n)-1;
912 fortint thisFile;
913 gribfile * file;
914
915 if( index >= 0 ) {
916 thisFile = pbginitInput(filename,filename_len);
917 file = currentFile(openFiles, thisFile);
918 if( index < file->count ) {
919 if( DEBUG1 ) {
920 char* pfile;
921 copyName(&pfile, filename, filename_len);
922 printf("PBGLEVL: level of GRIB %ld in file %s = %ld\n",
923 (index+1), pfile, file->level[index]);
924 free(pfile);
925 }
926 return ( file->level[index] );
927 }
928 }
929
930 return (-1);
931 }
932
pbglevl(char * filename,fortint * n,fortint filename_len)933 fortint pbglevl(char* filename, fortint * n, fortint filename_len) {
934 return pbglevl_(filename,n,filename_len);
935 }
936
937 /*
938 // ****************************************************************
939 */
pbgdate_(char * filename,fortint * n,fortint filename_len)940 fortint pbgdate_(char* filename, fortint * n, fortint filename_len) {
941 /*
942 // Returns date (as YYMMDD) of GRIB product n or -1 if n is greater
943 // than the number of products in the file, or less than 1.
944 */
945 fortint index = (*n)-1;
946 fortint thisFile;
947 gribfile * file;
948
949 if( index >= 0 ) {
950 thisFile = pbginitInput(filename,filename_len);
951 file = currentFile(openFiles, thisFile);
952 if( index < file->count ) {
953 if( DEBUG1 ) {
954 char* pfile;
955 copyName(&pfile, filename, filename_len);
956 printf("PBGDATE: date of GRIB %ld in file %s = %ld\n",
957 (index+1), pfile, file->date[index]);
958 free(pfile);
959 }
960 return ( file->date[index] );
961 }
962 }
963
964 return (-1);
965 }
966
pbgdate(char * filename,fortint * n,fortint filename_len)967 fortint pbgdate(char* filename, fortint * n, fortint filename_len) {
968 return pbgdate_(filename,n,filename_len);
969 }
970
971 /*
972 // ****************************************************************
973 */
pbgtime_(char * filename,fortint * n,fortint filename_len)974 fortint pbgtime_(char* filename, fortint * n, fortint filename_len) {
975 /*
976 // Returns time (as HHMM) of GRIB product n or -1 if n is greater
977 // than the number of products in the file, or less than 1.
978 */
979 fortint index = (*n)-1;
980 fortint thisFile;
981 gribfile * file;
982
983 if( index >= 0 ) {
984 thisFile = pbginitInput(filename,filename_len);
985 file = currentFile(openFiles, thisFile);
986 if( index < file->count ) {
987 if( DEBUG1 ) {
988 char* pfile;
989 copyName(&pfile, filename, filename_len);
990 printf("PBGTIME: time of GRIB %ld in file %s = %ld\n",
991 (index+1), pfile, file->time[index]);
992 free(pfile);
993 }
994 return ( file->time[index] );
995 }
996 }
997
998 return (-1);
999 }
1000
pbgtime(char * filename,fortint * n,fortint filename_len)1001 fortint pbgtime(char* filename, fortint * n, fortint filename_len) {
1002 return pbgtime_(filename,n,filename_len);
1003 }
1004
1005 /*
1006 // ****************************************************************
1007 */
pbgstep_(char * filename,fortint * n,fortint filename_len)1008 fortint pbgstep_(char* filename, fortint * n, fortint filename_len) {
1009 /*
1010 // Returns timestep of GRIB product n or -1 if n is greater
1011 // than the number of products in the file, or less than 1.
1012 */
1013 fortint index = (*n)-1;
1014 fortint thisFile;
1015 gribfile * file;
1016
1017 if( index >= 0 ) {
1018 thisFile = pbginitInput(filename,filename_len);
1019 file = currentFile(openFiles, thisFile);
1020 if( index < file->count ) {
1021 if( DEBUG1 ) {
1022 char* pfile;
1023 copyName(&pfile, filename, filename_len);
1024 printf("PBGSTEP: timestep of GRIB %ld in file %s = %ld\n",
1025 (index+1), pfile, file->timestep[index]);
1026 free(pfile);
1027 }
1028 return ( file->timestep[index] );
1029 }
1030 }
1031
1032 return (-1);
1033 }
1034
pbgstep(char * filename,fortint * n,fortint filename_len)1035 fortint pbgstep(char* filename, fortint * n, fortint filename_len) {
1036 return pbgstep_(filename,n,filename_len);
1037 }
1038
1039 /*
1040 // ****************************************************************
1041 */
pbgtype_(char * filename,fortint * n,fortint filename_len)1042 fortint pbgtype_(char* filename, fortint * n, fortint filename_len) {
1043 /*
1044 // Returns field type of GRIB product n or -1 if n is greater
1045 // than the number of products in the file, or less than 1.
1046 */
1047 fortint index = (*n)-1;
1048 fortint thisFile;
1049 gribfile * file;
1050
1051 if( index >= 0 ) {
1052 thisFile = pbginitInput(filename,filename_len);
1053 file = currentFile(openFiles, thisFile);
1054 if( index < file->count ) {
1055 if( DEBUG1 ) {
1056 char* pfile;
1057 copyName(&pfile, filename, filename_len);
1058 printf("PBGTYPE: type of GRIB %ld in file %s = %ld\n",
1059 (index+1), pfile, file->type[index]);
1060 free(pfile);
1061 }
1062 return ( file->type[index] );
1063 }
1064 }
1065
1066 return (-1);
1067 }
1068
pbgtype(char * filename,fortint * n,fortint filename_len)1069 fortint pbgtype(char* filename, fortint * n, fortint filename_len) {
1070 return pbgtype_(filename,n,filename_len);
1071 }
1072
1073 /*
1074 // ****************************************************************
1075 */
pbgstrm_(char * filename,fortint * n,fortint filename_len)1076 fortint pbgstrm_(char* filename, fortint * n, fortint filename_len) {
1077 /*
1078 // Returns stream of GRIB product n or -1 if n is greater
1079 // than the number of products in the file, or less than 1.
1080 */
1081 fortint index = (*n)-1;
1082 fortint thisFile;
1083 gribfile * file;
1084
1085 if( index >= 0 ) {
1086 thisFile = pbginitInput(filename,filename_len);
1087 file = currentFile(openFiles, thisFile);
1088 if( index < file->count ) {
1089 if( DEBUG1 ) {
1090 char* pfile;
1091 copyName(&pfile, filename, filename_len);
1092 printf("PBGSTRM: stream of GRIB %ld in file %s = %ld\n",
1093 (index+1), pfile, file->stream[index]);
1094 free(pfile);
1095 }
1096 return ( file->stream[index] );
1097 }
1098 }
1099
1100 return (-1);
1101 }
1102
pbgstrm(char * filename,fortint * n,fortint filename_len)1103 fortint pbgstrm(char* filename, fortint * n, fortint filename_len) {
1104 return pbgstrm_(filename,n,filename_len);
1105 }
1106
1107 /*
1108 // ****************************************************************
1109 */
pbgrepr_(char * filename,fortint * n,fortint filename_len)1110 fortint pbgrepr_(char* filename, fortint * n, fortint filename_len) {
1111 /*
1112 // Returns repres of GRIB product n or -1 if n is greater
1113 // than the number of products in the file, or less than 1.
1114 */
1115 fortint index = (*n)-1;
1116 fortint thisFile;
1117 gribfile * file;
1118
1119 if( index >= 0 ) {
1120 thisFile = pbginitInput(filename,filename_len);
1121 file = currentFile(openFiles, thisFile);
1122 if( index < file->count ) {
1123 if( DEBUG1 ) {
1124 char* pfile;
1125 copyName(&pfile, filename, filename_len);
1126 printf("PBGREPR: repres of GRIB %ld in file %s = %ld\n",
1127 (index+1), pfile, file->repres[index]);
1128 free(pfile);
1129 }
1130 return ( file->repres[index] );
1131 }
1132 }
1133
1134 return (-1);
1135 }
1136
pbgrepr(char * filename,fortint * n,fortint filename_len)1137 fortint pbgrepr(char* filename, fortint * n, fortint filename_len) {
1138 return pbgrepr_(filename,n,filename_len);
1139 }
1140
1141 /*
1142 // ****************************************************************
1143 */
pbglevt_(char * filename,fortint * n,fortint filename_len)1144 fortint pbglevt_(char* filename, fortint * n, fortint filename_len) {
1145 /*
1146 // Returns levtype of GRIB product n or -1 if n is greater
1147 // than the number of products in the file, or less than 1.
1148 */
1149 fortint index = (*n)-1;
1150 fortint thisFile;
1151 gribfile * file;
1152
1153 if( index >= 0 ) {
1154 thisFile = pbginitInput(filename,filename_len);
1155 file = currentFile(openFiles, thisFile);
1156 if( index < file->count ) {
1157 if( DEBUG1 ) {
1158 char* pfile;
1159 copyName(&pfile, filename, filename_len);
1160 printf("PBGLEVT: levtype of GRIB %ld in file %s = %ld\n",
1161 (index+1), pfile, file->levtype[index]);
1162 free(pfile);
1163 }
1164 return ( file->levtype[index] );
1165 }
1166 }
1167
1168 return (-1);
1169 }
1170
pbglevt(char * filename,fortint * n,fortint filename_len)1171 fortint pbglevt(char* filename, fortint * n, fortint filename_len) {
1172 return pbglevt_(filename,n,filename_len);
1173 }
1174
1175 /*
1176 // ****************************************************************
1177 */
pbgfind_(char * filename,fortint * param,fortint * level,fortint * date,fortint * time,fortint * step,fortint * n,fortint filename_len)1178 fortint pbgfind_(
1179 char* filename,
1180 fortint * param,
1181 fortint * level,
1182 fortint * date,
1183 fortint * time,
1184 fortint * step,
1185 fortint * n,
1186 fortint filename_len) {
1187 /*
1188 // Returns the index of the GRIB product with the specified param, level, ..
1189 // Searching starts after GRIB product n.
1190 */
1191 fortint index = (*n);
1192 fortint p = *param, l = *level, d = *date, t = *time, s = *step;
1193 fortint loop;
1194 fortint thisFile;
1195 gribfile * file;
1196
1197 if( index >= 0 ) {
1198 thisFile = pbginitInput(filename,filename_len);
1199 file = currentFile(openFiles, thisFile);
1200
1201 for( loop = index; loop < file->count; loop++ ) {
1202 if( p == file->parameter[loop] )
1203 if( l == file->level[loop] )
1204 if( d == file->date[loop] )
1205 if( t == file->time[loop] )
1206 if( s == file->timestep[loop] )
1207 return (loop+1);
1208 }
1209 }
1210
1211 return (-1);
1212 }
1213
pbgfind(char * filename,fortint * param,fortint * level,fortint * date,fortint * time,fortint * step,fortint * n,fortint filename_len)1214 fortint pbgfind(
1215 char* filename,
1216 fortint * param,
1217 fortint * level,
1218 fortint * date,
1219 fortint * time,
1220 fortint * step,
1221 fortint * n,
1222 fortint filename_len) {
1223 return pbgfind_(filename,param,level,date,time,step,n,filename_len);
1224 }
1225
1226 /*
1227 // ****************************************************************
1228 */
pbgvfind_(char * filename,fortint * param,fortint * level,fortint * vdate,fortint * vtime,fortint * status,fortint * n,fortint filename_len)1229 fortint pbgvfind_(
1230 char* filename,
1231 fortint * param,
1232 fortint * level,
1233 fortint * vdate,
1234 fortint * vtime,
1235 fortint * status,
1236 fortint * n,
1237 fortint filename_len) {
1238 /*
1239 // Returns the index of the GRIB product with the specified param, level,
1240 // verifying data and time (vdate & vtime).
1241 // Searching starts after GRIB product n.
1242 //
1243 // Returns the index of the first match found, or -1 if no match is found.
1244 // If more than one match is found (eg analysis and forecast fields),
1245 // status returns the index of the second match; otherwise status = 0.
1246 */
1247 fortint index = (*n);
1248 fortint p = *param, l = *level, d = *vdate, t = *vtime;
1249 fortint loop;
1250 fortint thisFile;
1251 gribfile * file;
1252
1253 *status = 0;
1254
1255 if( index >= 0 ) {
1256 thisFile = pbginitInput(filename,filename_len);
1257 file = currentFile(openFiles, thisFile);
1258
1259 for( loop = index; loop < file->count; loop++ ) {
1260 if( p == file->parameter[loop] )
1261 if( l == file->level[loop] )
1262 if( d == file->vdate[loop] )
1263 if( t == file->vtime[loop] ) {
1264 long innerloop;
1265 for( innerloop = loop+1; innerloop < file->count; innerloop++ ) {
1266 if( p == file->parameter[innerloop] )
1267 if( l == file->level[innerloop] )
1268 if( d == file->vdate[innerloop] )
1269 if( t == file->vtime[innerloop] )
1270 *status = innerloop+1 ;
1271 }
1272 return (loop+1);
1273 }
1274 }
1275 }
1276
1277 return (-1);
1278 }
1279
pbgvfind(char * filename,fortint * param,fortint * level,fortint * vdate,fortint * vtime,fortint * status,fortint * n,fortint filename_len)1280 fortint pbgvfind(
1281 char* filename,
1282 fortint * param,
1283 fortint * level,
1284 fortint * vdate,
1285 fortint * vtime,
1286 fortint * status,
1287 fortint * n,
1288 fortint filename_len) {
1289 return pbgvfind_(filename,param,level,vdate,vtime,status,n,filename_len);
1290 }
1291
1292 /*
1293 // ****************************************************************
1294 */
1295 #define L_PARAM 0
1296 #define L_LEVEL 1
1297 #define L_DATE 2
1298 #define L_TIME 3
1299 #define L_STEP 4
1300 #define L_ENSEMBLE 5
1301 #define L_TYPE 6
1302 #define L_STREAM 7
1303 #define L_REPRES 8
1304 #define L_LEVTYPE 9
1305 #define L_VDATE 10
1306 #define L_VTIME 11
1307 #define L_NI 12
1308 #define L_NJ 13
1309 #define L_NLAT 14
1310 #define L_NLON 15
1311 #define L_SLAT 16
1312 #define L_SLON 17
1313 #define L_DI 18
1314 #define L_DJ 19
1315 #define L_SPLAT 20
1316 #define L_SPLON 21
1317 #define L_QUASI 22
1318
1319 #define L_PROBSCALE 23
1320 #define L_PROBINDIC 24
1321 #define L_PROBLOWER 25
1322 #define L_PROBUPPER 26
1323
1324 #define L_TUBENUM 23
1325 #define L_TUBETOT 24
1326 #define L_TUBEN 25
1327 #define L_TUBEW 26
1328 #define L_TUBES 27
1329 #define L_TUBEE 28
1330 #define L_TUBESTP 29
1331 #define L_TUBEFC 30
1332
1333 #define L_CLUSNUM 23
1334 #define L_CLUSTOT 24
1335 #define L_CLUSMETH 25
1336 #define L_CLUSSTP 26
1337 #define L_CLUSSTPE 27
1338 #define L_CLUSN 28
1339 #define L_CLUSW 29
1340 #define L_CLUSS 30
1341 #define L_CLUSE 31
1342 #define L_CLUSOPFC 32
1343 #define L_CLUSCFFC 33
1344 #define L_CLUSFC 34
1345
1346 #define L_DIRECTIONNUMBER 23
1347 #define L_FREQUENCYNUMBER 24
1348
1349 #define L_OPTIMISATIONTIME 23
1350 #define L_LEADTIME 24
1351 #define L_SENSITIVEAREADOMAIN 25
1352 #define L_SENSITIVEAREAMETHOD 26
1353
1354 #define L_VERIFYINGMONTH 23
1355 #define L_AVERAGINGPERIOD 24
1356 #define L_FORECASTMONTH 25
1357 #define L_REFERENCEDATE 26
1358 #define L_CLIMATEDATEFROM 27
1359 #define L_CLIMATEDATETO 28
1360 #define L_THRESHOLDUNITSSCALE 29
1361 #define L_THRESHOLDINDICATOR 30
1362 #define L_LOWERTHRESHOLD 31
1363 #define L_UPPERTHRESHOLD 32
1364
1365 #define L_LOCALDEF 35
1366
1367 #define L_ALIST_LENGTH 12
1368 #define L_BLIST_LENGTH 36
1369
pbgafind_(char * filename,fortint * list,fortint * n,fortint filename_len)1370 fortint pbgafind_(
1371 char* filename,
1372 fortint * list,
1373 fortint * n,
1374 fortint filename_len) {
1375 /*
1376 // Returns the index of the GRIB product with the specified param, level, ..
1377 // Searching starts after GRIB product n.
1378 */
1379 fortint index = (*n);
1380 fortint param = list[L_PARAM],
1381 level = list[L_LEVEL],
1382 date = list[L_DATE],
1383 time = list[L_TIME],
1384 step = list[L_STEP],
1385 ensemble = list[L_ENSEMBLE],
1386 type = list[L_TYPE],
1387 stream = list[L_STREAM],
1388 repres = list[L_REPRES],
1389 levtype = list[L_LEVTYPE],
1390 vdate = list[L_VDATE],
1391 vtime = list[L_VTIME];
1392 fortint loop;
1393 fortint thisFile;
1394 gribfile * file;
1395
1396 if( DEBUG1 ) {
1397 char* pfile;
1398 fortint loop;
1399
1400 copyName(&pfile, filename, filename_len);
1401 printf("PBGAFIND: searching file %s\n", pfile);
1402 free(pfile);
1403
1404 for( loop = 0; loop < L_ALIST_LENGTH ; loop++ ) {
1405 if( list[loop] != MISSING )
1406 printf("PBGAFIND: ilist[%ld] = %ld\n", (loop+1), list[loop] );
1407 }
1408 }
1409
1410 if( index >= 0 ) {
1411 thisFile = pbginitInput(filename,filename_len);
1412 file = currentFile(openFiles, thisFile);
1413
1414 for( loop = index; loop < file->count; loop++ ) {
1415 if( MATCH(param,file->parameter[loop]) )
1416 if( MATCH(level,file->level[loop]) )
1417 if( MATCH(date,file->date[loop]) )
1418 if( MATCH(time,file->time[loop]) )
1419 if( MATCH(step,file->timestep[loop]) )
1420 if( MATCH(ensemble,file->number[loop]) )
1421 if( MATCH(type,file->type[loop]) )
1422 if( MATCH(stream,file->stream[loop]) )
1423 if( MATCH(repres,file->repres[loop]) )
1424 if( MATCH(levtype,file->levtype[loop]) )
1425 if( MATCH(vdate,file->vdate[loop]) )
1426 if( MATCH(vtime,file->vtime[loop]) ) {
1427 if( DEBUG1 ) {
1428 char* pfile;
1429 copyName(&pfile, filename, filename_len);
1430 printf("PBGAFIND: matching GRIB found at position %ld in file %s\n",
1431 (loop+1), pfile);
1432 free(pfile);
1433 }
1434 return (loop+1);
1435 }
1436 }
1437 }
1438
1439 return (-1);
1440 }
1441
pbgafind(char * filename,fortint * list,fortint * n,fortint filename_len)1442 fortint pbgafind(char* filename,fortint* list,fortint* n,fortint filename_len) {
1443 return pbgafind_(filename,list,n,filename_len);
1444 }
1445
1446 /*
1447 // ****************************************************************
1448 */
pbgbfind_(char * filename,fortint * list,fortint * n,fortint filename_len)1449 fortint pbgbfind_(
1450 char* filename,
1451 fortint * list,
1452 fortint * n,
1453 fortint filename_len) {
1454 /*
1455 // Returns the index of the GRIB product with the specified param, level, ..
1456 // Searching starts after GRIB product n.
1457 */
1458 fortint index = (*n);
1459 fortint param = list[L_PARAM],
1460 level = list[L_LEVEL],
1461 date = list[L_DATE],
1462 time = list[L_TIME],
1463 step = list[L_STEP],
1464 localDef = list[L_LOCALDEF],
1465 ensemble = list[L_ENSEMBLE],
1466 type = list[L_TYPE],
1467 stream = list[L_STREAM],
1468 repres = list[L_REPRES],
1469 levtype = list[L_LEVTYPE],
1470 vdate = list[L_VDATE],
1471 vtime = list[L_VTIME],
1472 ni = list[L_NI],
1473 nj = list[L_NJ],
1474 nlat = list[L_NLAT],
1475 nlon = list[L_NLON],
1476 slat = list[L_SLAT],
1477 slon = list[L_SLON],
1478 di = list[L_DI],
1479 dj = list[L_DJ],
1480 splat = list[L_SPLAT],
1481 splon = list[L_SPLON],
1482 quasi = list[L_QUASI];
1483 fortint tcNumber = MISSING,
1484 tcTotal = MISSING,
1485 clusterMethod = MISSING,
1486 tcStep = MISSING,
1487 clusterStepEnd = MISSING,
1488 tcNorth = MISSING,
1489 tcWest = MISSING,
1490 tcSouth = MISSING,
1491 tcEast = MISSING,
1492 clusterOpFcNumber = MISSING,
1493 clusterControlFcNumber = MISSING,
1494 tcNumOfFcs = MISSING,
1495 probScale = MISSING,
1496 probIndicator = MISSING,
1497 probLower = MISSING,
1498 probUpper = MISSING,
1499 directionNumber = MISSING,
1500 frequencyNumber = MISSING,
1501 optimisationTime = MISSING,
1502 leadTime = MISSING,
1503 sensitiveAreaDomain = MISSING,
1504 sensitiveAreaMethod = MISSING,
1505 verifyingMonth = MISSING,
1506 averagingPeriod = MISSING,
1507 forecastMonth = MISSING,
1508 referenceDate = MISSING,
1509 climateDateFrom = MISSING,
1510 climateDateTo = MISSING,
1511 thresholdUnitsScale = MISSING,
1512 thresholdIndicator = MISSING,
1513 lowerThreshold = MISSING,
1514 upperThreshold = MISSING;
1515 fortint loop;
1516 fortint thisFile;
1517 gribfile * file;
1518
1519 if( DEBUG1 ) {
1520 char* pfile;
1521 fortint loop;
1522
1523 copyName(&pfile, filename, filename_len);
1524 printf("PBGBFIND: searching file %s\n", pfile);
1525 free(pfile);
1526
1527 for( loop = 0; loop < L_BLIST_LENGTH ; loop++ ) {
1528 if( list[loop] != MISSING )
1529 printf("PBGBFIND: ilist[%ld] = %ld\n", (loop+1), list[loop] );
1530 }
1531 }
1532
1533 if( (type == CLUSTERMEAN) || (type == CLUSTERSD) ) {
1534 tcNumber = list[L_CLUSNUM];
1535 tcTotal = list[L_CLUSTOT];
1536 clusterMethod = list[L_CLUSMETH];
1537 tcStep = list[L_CLUSSTP];
1538 clusterStepEnd = list[L_CLUSSTPE];
1539 tcNorth = list[L_CLUSN];
1540 tcWest = list[L_CLUSW];
1541 tcSouth = list[L_CLUSS];
1542 tcEast = list[L_CLUSE];
1543 clusterOpFcNumber = list[L_CLUSOPFC];
1544 clusterControlFcNumber = list[L_CLUSCFFC];
1545 tcNumOfFcs = list[L_CLUSFC];
1546 }
1547
1548 if( type == TUBETYPE) {
1549 tcNumber = list[L_TUBENUM];
1550 tcTotal = list[L_TUBETOT];
1551 tcStep = list[L_TUBESTP];
1552 tcNorth = list[L_TUBEN];
1553 tcWest = list[L_TUBEW];
1554 tcSouth = list[L_TUBES];
1555 tcEast = list[L_TUBEE];
1556 tcNumOfFcs = list[L_TUBEFC];
1557 }
1558
1559 if( type == PROBABILITY) {
1560 probScale = list[L_PROBSCALE];
1561 probIndicator = list[L_PROBINDIC];
1562 probLower = list[L_PROBLOWER];
1563 probUpper = list[L_PROBUPPER];
1564 }
1565
1566 if( localDef == 13 ) {
1567 directionNumber = list[L_DIRECTIONNUMBER];
1568 frequencyNumber = list[L_FREQUENCYNUMBER];
1569 }
1570
1571 if( localDef == 21 ) {
1572 optimisationTime = list[L_OPTIMISATIONTIME];
1573 leadTime = list[L_LEADTIME];
1574 sensitiveAreaDomain = list[L_SENSITIVEAREADOMAIN];
1575 sensitiveAreaMethod = list[L_SENSITIVEAREAMETHOD];
1576 }
1577
1578 if( localDef == 23 ) {
1579 verifyingMonth = list[L_VERIFYINGMONTH];
1580 averagingPeriod = list[L_AVERAGINGPERIOD];
1581 forecastMonth = list[L_FORECASTMONTH];
1582 referenceDate = list[L_REFERENCEDATE];
1583 climateDateFrom = list[L_CLIMATEDATEFROM];
1584 climateDateTo = list[L_CLIMATEDATETO];
1585 thresholdUnitsScale = list[L_THRESHOLDUNITSSCALE];
1586 thresholdIndicator = list[L_THRESHOLDINDICATOR];
1587 lowerThreshold = list[L_LOWERTHRESHOLD];
1588 upperThreshold = list[L_UPPERTHRESHOLD];
1589 }
1590
1591 if( index >= 0 ) {
1592 thisFile = pbginitInput(filename,filename_len);
1593 file = currentFile(openFiles, thisFile);
1594
1595 for( loop = index; loop < file->count; loop++ ) {
1596 if( MATCH(param,file->parameter[loop]) )
1597 if( MATCH(level,file->level[loop]) )
1598 if( MATCH(step,file->timestep[loop]) )
1599 if( MATCH(time,file->time[loop]) )
1600 if( MATCH(date,file->date[loop]) )
1601 if( MATCH(ensemble,file->number[loop]) )
1602 if( MATCH(type,file->type[loop]) )
1603 if( MATCH(stream,file->stream[loop]) )
1604 if( MATCH(repres,file->repres[loop]) )
1605 if( MATCH(levtype,file->levtype[loop]) )
1606 if( MATCH(vdate,file->vdate[loop]) )
1607 if( MATCH(vtime,file->vtime[loop]) )
1608 if( MATCH(ni,file->ni[loop]) )
1609 if( MATCH(nj,file->nj[loop]) )
1610 if( MATCH(nlat,file->nlat[loop]) )
1611 if( MATCH(nlon,file->nlon[loop]) )
1612 if( MATCH(slat,file->slat[loop]) )
1613 if( MATCH(slon,file->slon[loop]) )
1614 if( MATCH(di,file->di[loop]) )
1615 if( MATCH(dj,file->dj[loop]) )
1616 if( MATCH(splat,file->splat[loop]) )
1617 if( MATCH(splon,file->splon[loop]) )
1618 if( MATCH(quasi,file->quasi[loop]) )
1619 if( MATCH(tcNumber,file->tcNumber[loop]) )
1620 if( MATCH(tcTotal,file->tcTotal[loop]) )
1621 if( MATCH(clusterMethod,file->clusterMethod[loop]) )
1622 if( MATCH(tcStep,file->tcStep[loop]) )
1623 if( MATCH(clusterStepEnd,file->clusterStepEnd[loop]) )
1624 if( MATCH(tcNorth,file->tcNorth[loop]) )
1625 if( MATCH(tcWest,file->tcWest[loop]) )
1626 if( MATCH(tcSouth,file->tcSouth[loop]) )
1627 if( MATCH(tcEast,file->tcEast[loop]) )
1628 if( MATCH(clusterOpFcNumber,file->clusterOpFcNumber[loop]) )
1629 if( MATCH(clusterControlFcNumber,file->clusterControlFcNumber[loop]) )
1630 if( MATCH(tcNumOfFcs,file->tcNumOfFcs[loop]) )
1631 if( MATCH(probScale,file->probScale[loop]) )
1632 if( MATCH(probIndicator,file->probIndicator[loop]) )
1633 if( MATCH(probLower,file->probLower[loop]) )
1634 if( MATCH(probUpper,file->probUpper[loop]) )
1635 if( MATCH(directionNumber,file->directionNumber[loop]) )
1636 if( MATCH(frequencyNumber,file->frequencyNumber[loop]) )
1637 if( MATCH(optimisationTime,file->optimisationTime[loop]) )
1638 if( MATCH(leadTime,file->leadTime[loop]) )
1639 if( MATCH(sensitiveAreaDomain,file->sensitiveAreaDomain[loop]) )
1640 if( MATCH(sensitiveAreaMethod,file->sensitiveAreaMethod[loop]) )
1641 if( MATCH(verifyingMonth,file->verifyingMonth[loop]) )
1642 if( MATCH(averagingPeriod,file->averagingPeriod[loop]) )
1643 if( MATCH(forecastMonth,file->forecastMonth[loop]) )
1644 if( MATCH(referenceDate,file->referenceDate[loop]) )
1645 if( MATCH(climateDateFrom,file->climateDateFrom[loop]) )
1646 if( MATCH(climateDateTo,file->climateDateTo[loop]) )
1647 if( MATCH(thresholdUnitsScale,file->thresholdUnitsScale[loop]) )
1648 if( MATCH(thresholdIndicator,file->thresholdIndicator[loop]) )
1649 if( MATCH(lowerThreshold,file->lowerThreshold[loop]) )
1650 if( MATCH(upperThreshold,file->upperThreshold[loop]) ) {
1651 if( DEBUG1 ) {
1652 char* pfile;
1653 copyName(&pfile, filename, filename_len);
1654 printf("PBGBFIND: matching GRIB found at position %ld in file %s\n",
1655 (loop+1), pfile);
1656 free(pfile);
1657 }
1658 return (loop+1);
1659 }
1660 }
1661 }
1662
1663 return (-1);
1664 }
1665
pbgbfind(char * filename,fortint * list,fortint * n,fortint filename_len)1666 fortint pbgbfind(char* filename,fortint* list,fortint* n,fortint filename_len) {
1667 return pbgbfind_(filename,list,n,filename_len);
1668 }
1669
1670 /*
1671 // ****************************************************************
1672 */
1673
1674 #define INDEX_TOO_BIG -1
1675 #define FSEEK_ERROR -2
1676 #define FREAD_ERROR -2
1677 #define BUFF_TOO_SMALL -3
1678
pbgget_(char * filename,fortint * buffer,fortint * bufflen,fortint * n,fortint filename_len)1679 fortint pbgget_(char* filename, fortint * buffer, fortint * bufflen,
1680 fortint * n, fortint filename_len) {
1681 /*
1682 // Gets the nth GRIB product.
1683 */
1684 fortint index = (*n)-1;
1685 fortint length, status;
1686 off_t offset;
1687 fortint thisFile;
1688 FILE * fp;
1689 gribfile * file;
1690
1691 if( DEBUG1 ) {
1692 char* pfile;
1693 copyName(&pfile, filename, filename_len);
1694 printf("PBGGET: getting GRIB number %ld in file %s\n",
1695 (index+1), pfile);
1696 free(pfile);
1697 }
1698
1699 if( index >= 0 ) {
1700 thisFile = pbginitInput(filename,filename_len);
1701 file = currentFile(openFiles, thisFile);
1702 if( index < file->count ) {
1703
1704 length = file->length[index];
1705
1706 if( DEBUG1 ) printf("PBGGET: length of GRIB number %ld = %ld\n",
1707 (index+1), length);
1708
1709 if( *bufflen < length ) {
1710 fprintf(stderr,
1711 "PBGGET: user buffer too small, %ld bytes required\n", length);
1712 return BUFF_TOO_SMALL;
1713 }
1714
1715 offset = file->offset[index];
1716
1717 if( DEBUG1 ) printf("PBGGET: offset of GRIB number %ld = %lld\n",
1718 (index+1), offset);
1719
1720 fp = file->fp;
1721 #ifdef FOPEN64
1722 if( fseek( fp, offset, 0) ) {
1723 #else
1724 if( fseek( fp, offset, 0) ) {
1725 #endif
1726 perror("PBGGET: error in fseek");
1727 return FSEEK_ERROR;
1728 }
1729
1730 status = fread(buffer, 1, length, fp);
1731 if( status != length ) {
1732 fprintf(stderr,"PBGGET: error in fread\n");
1733 return FREAD_ERROR;
1734 }
1735
1736 return length;
1737 }
1738 }
1739
1740 return INDEX_TOO_BIG;
1741
1742 }
1743
1744 fortint pbgget(char* filename, fortint * buffer, fortint * bufflen,
1745 fortint * n, fortint filename_len) {
1746 return pbgget_(filename,buffer,bufflen,n,filename_len);
1747 }
1748
1749 /*
1750 // ****************************************************************
1751 */
1752
1753 #define FWRITE_ERROR -1
1754
1755 fortint pbgput_(char* filename, fortint * buffer, fortint * bufflen,
1756 fortint filename_len) {
1757 /*
1758 // Write an array to a file.
1759 */
1760 fortint length;
1761 fortint thisFile;
1762 FILE * fp;
1763 gribfile * file;
1764
1765 if( DEBUG1 ) {
1766 char* pfile;
1767 copyName(&pfile, filename, filename_len);
1768 printf("PBGPUT: putting %ld bytes to file %s\n", *bufflen, pfile);
1769 free(pfile);
1770 }
1771
1772 thisFile = pbginitOutput(filename,filename_len);
1773 file = currentFile(openFiles, thisFile);
1774
1775 fp = file->fp;
1776 length = fwrite(buffer, 1, *bufflen, fp);
1777
1778 if( DEBUG1 ) printf("PBGPUT: number of bytes written = %ld\n", length);
1779
1780 if( length!=(*bufflen) )
1781 return FWRITE_ERROR;
1782 else
1783 return (*bufflen);
1784 }
1785
1786 fortint pbgput(char* filename, fortint * buffer, fortint * bufflen,
1787 fortint filename_len) {
1788 return pbgput_(filename,buffer,bufflen,filename_len);
1789 }
1790
1791 /*
1792 // ****************************************************************
1793 */
1794 void pbgclos_(char* filename, fortint filename_len) {
1795
1796 openFiles.removeFile(filename, filename_len);
1797 return;
1798
1799 }
1800
1801 void pbgclos(char* filename, fortint filename_len) {
1802 pbgclos_(filename,filename_len);
1803 }
1804
1805 /*
1806 // ****************************************************************
1807 */
1808 fortint pbginitInput(char* filename, fortint filename_len) {
1809 /*
1810 // Initialise indices for current file if necessary.
1811 */
1812 fortint thisFile;
1813
1814 if( DEBUG2 ) {
1815 char* pfile;
1816 copyName(&pfile, filename, filename_len);
1817 printf("PBG_pbginitInput: checking if file %s already open\n", pfile);
1818 free(pfile);
1819 }
1820
1821 thisFile = openFiles.exists(filename,filename_len);
1822
1823 if( thisFile == -1 ) {
1824 if( DEBUG2 ) printf("PBG_pbginitInput: file not yet open\n");
1825 return ( openFiles.addRead(filename, filename_len) );
1826 }
1827
1828 if( DEBUG2 ) printf("PBG_pbginitInput: file has open slot %ld\n", thisFile);
1829 return thisFile;
1830
1831 }
1832
1833 fortint pbginitOutput(char* filename, fortint filename_len) {
1834 /*
1835 // Initialise output file if necessary.
1836 */
1837 fortint thisFile;
1838
1839 if( DEBUG2 ) {
1840 char* pfile;
1841 copyName(&pfile, filename, filename_len);
1842 printf("PBG_pbginitOutput: checking if file %s already open\n", pfile);
1843 free(pfile);
1844 }
1845
1846 thisFile = openFiles.exists(filename,filename_len);
1847
1848 if( thisFile == -1 ) {
1849 if( DEBUG2 ) printf("PBG_pbginitOutput: file not yet open\n");
1850 return ( openFiles.addWrite(filename, filename_len) );
1851 }
1852
1853 if( DEBUG2 ) printf("PBG_pbginitOutput: file has open slot %ld\n", thisFile);
1854 return thisFile;
1855
1856 }
1857
1858 fortint exists(char* filename , fortint filename_len) {
1859 fortint n;
1860 char* pfile;
1861 gribfile * file;
1862
1863
1864 /*
1865 // See if DEBUG switched on.
1866 */
1867 if( ! debugSet ) {
1868 debugLevel = getenv("PBG_DEBUG");
1869 if( debugLevel == NULL )
1870 debugSet = DEBUGOFF; /* off */
1871 else {
1872 int loop;
1873 for( loop = 0; loop < strlen(debugLevel) ; loop++ ) {
1874 if( ! isdigit(debugLevel[loop]) ) {
1875 printf("Invalid number string in PBG_DEBUG: %s\n", debugLevel);
1876 printf("PBG_DEBUG must comprise only digits [0-9].\n");
1877 debugSet = DEBUGOFF;
1878 }
1879 }
1880 debugSet = DEBUGOFF + atol( debugLevel );
1881 }
1882 if( DEBUG2 ) printf("PBG_exists: PBG_DEBUG switched on\n");
1883 }
1884
1885 if( openFiles.count > -1 ) {
1886
1887 copyName(&pfile, filename, filename_len);
1888
1889 if( DEBUG2 ) printf("PBG_exists: looking for filename = %s\n", pfile);
1890
1891 for( n = 0; n <= openFiles.count; n++ ) {
1892 file = currentFile(openFiles, n);
1893 if( file != 0 ) {
1894 if( strcmp( file->fname, pfile ) == 0 ) {
1895 free(pfile);
1896 if( DEBUG2 ) printf("PBG_exists: file found in slot= %ld\n", n);
1897 return n;
1898 }
1899 }
1900 }
1901
1902 free(pfile);
1903 }
1904 return -1;
1905 }
1906
1907 fortint addFile(char* filename , fortint filename_len, char readwriteflag) {
1908 char* pfile;
1909 char mode[2] = " ";
1910 FILE * in;
1911 gribfile * previousLatest, * latest;
1912
1913 openFiles.count++;
1914
1915 copyName(&pfile, filename, filename_len);
1916
1917 if( DEBUG2 ) printf("PBG_addFile: adding filename = %s\n", pfile);
1918
1919 mode[0] = readwriteflag;
1920
1921 in = fopen(pfile,mode);
1922
1923 if( in == NULL ) {
1924 perror("Error opening file");
1925 exit(1);
1926 }
1927
1928 if( openFiles.count == 0 ) {
1929 openFiles.files = (gribfile *)malloc(sizeof(gribfile));
1930 latest = openFiles.files;
1931 }
1932 else {
1933 previousLatest = latestFile(openFiles);
1934 previousLatest->next = (gribfile *)malloc(sizeof(gribfile));
1935 latest = previousLatest->next;
1936 }
1937
1938 latest->fp = in;
1939 copyName(&(latest->fname),filename,filename_len);
1940 latest->readwriteflag = readwriteflag;
1941 latest->max = MAX_NUMBER_OF_GRIBS;
1942 latest->count = 0;
1943
1944 #define MALLOC(a) (a) = (void *) malloc(space)
1945
1946 if( readwriteflag == 'r' ) {
1947 fortint space;
1948
1949 space = sizeof(off_t)*MAX_NUMBER_OF_GRIBS;
1950 MALLOC(latest->offset);
1951
1952 space = sizeof(fortint)*MAX_NUMBER_OF_GRIBS;
1953 MALLOC(latest->length);
1954 MALLOC(latest->parameter);
1955 MALLOC(latest->level);
1956 MALLOC(latest->date);
1957 MALLOC(latest->time);
1958 MALLOC(latest->timestep);
1959 MALLOC(latest->localDefinitionNumber);
1960 MALLOC(latest->type);
1961 MALLOC(latest->stream);
1962 MALLOC(latest->repres);
1963 MALLOC(latest->levtype);
1964 MALLOC(latest->number);
1965 MALLOC(latest->vdate);
1966 MALLOC(latest->vtime);
1967 MALLOC(latest->tcNumber);
1968 MALLOC(latest->tcTotal);
1969 MALLOC(latest->clusterMethod);
1970 MALLOC(latest->tcStep);
1971 MALLOC(latest->clusterStepEnd);
1972 MALLOC(latest->tcNorth);
1973 MALLOC(latest->tcWest);
1974 MALLOC(latest->tcSouth);
1975 MALLOC(latest->tcEast);
1976 MALLOC(latest->clusterOpFcNumber);
1977 MALLOC(latest->clusterControlFcNumber);
1978 MALLOC(latest->tcNumOfFcs);
1979 MALLOC(latest->probScale);
1980 MALLOC(latest->probIndicator);
1981 MALLOC(latest->probLower);
1982 MALLOC(latest->probUpper);
1983 MALLOC(latest->ni);
1984 MALLOC(latest->nj);
1985 MALLOC(latest->nlat);
1986 MALLOC(latest->nlon);
1987 MALLOC(latest->slat);
1988 MALLOC(latest->slon);
1989 MALLOC(latest->di);
1990 MALLOC(latest->dj);
1991 MALLOC(latest->splat);
1992 MALLOC(latest->splon);
1993 MALLOC(latest->quasi);
1994 MALLOC(latest->directionNumber);
1995 MALLOC(latest->frequencyNumber);
1996 MALLOC(latest->optimisationTime);
1997 MALLOC(latest->leadTime);
1998 MALLOC(latest->sensitiveAreaDomain);
1999 MALLOC(latest->sensitiveAreaMethod);
2000 MALLOC(latest->verifyingMonth);
2001 MALLOC(latest->averagingPeriod);
2002 MALLOC(latest->forecastMonth);
2003 MALLOC(latest->referenceDate);
2004 MALLOC(latest->climateDateFrom);
2005 MALLOC(latest->climateDateTo);
2006 MALLOC(latest->thresholdUnitsScale);
2007 MALLOC(latest->thresholdIndicator);
2008 MALLOC(latest->lowerThreshold);
2009 MALLOC(latest->upperThreshold);
2010 }
2011
2012 latest->next = 0;
2013
2014 if( DEBUG2 ) printf("PBG_addFile: adding file %s in slot = %ld\n",
2015 pfile, openFiles.count);
2016
2017 free(pfile);
2018 return openFiles.count;
2019
2020 }
2021
2022 void removeFile(char* filename , fortint filename_len) {
2023 fortint thisFile;
2024 int status;
2025
2026 if( DEBUG2 ) {
2027 char* pfile;
2028 copyName(&pfile, filename, filename_len);
2029 printf("PBG_removeFile: trying to remove filename = %s\n", pfile);
2030 free(pfile);
2031 }
2032
2033 thisFile = openFiles.exists(filename,filename_len);
2034
2035 if( thisFile != -1 ) {
2036 fortint index;
2037 gribfile * previous = openFiles.files;
2038 gribfile * current;
2039
2040 current = previous;
2041
2042 for( index = 0; index < thisFile; index++ ) {
2043 previous = current;
2044 current = current->next;
2045 }
2046
2047 status = fclose(current->fp);
2048 if( status != 0 ) {
2049 perror("Error closing file");
2050 exit(1);
2051 }
2052 if( DEBUG2 ) printf("PBG_removeFile: removing file %s from slot %ld\n",
2053 current->fname, thisFile);
2054 free(current->fname);
2055
2056 /*
2057 // Input files have arrays of values
2058 */
2059 if( current->readwriteflag == 'r' ) {
2060 free(current->offset);
2061 free(current->length);
2062 free(current->parameter);
2063 free(current->level);
2064 free(current->date);
2065 free(current->time);
2066 free(current->timestep);
2067 free(current->localDefinitionNumber);
2068 free(current->type);
2069 free(current->stream);
2070 free(current->repres);
2071 free(current->levtype);
2072 free(current->number);
2073 free(current->vdate);
2074 free(current->vtime);
2075 free(current->tcNumber);
2076 free(current->tcTotal);
2077 free(current->clusterMethod);
2078 free(current->tcStep);
2079 free(current->clusterStepEnd);
2080 free(current->tcNorth);
2081 free(current->tcWest);
2082 free(current->tcSouth);
2083 free(current->tcEast);
2084 free(current->clusterOpFcNumber);
2085 free(current->clusterControlFcNumber);
2086 free(current->tcNumOfFcs);
2087 free(current->probScale);
2088 free(current->probIndicator);
2089 free(current->probLower);
2090 free(current->probUpper);
2091 free(current->ni);
2092 free(current->nj);
2093 free(current->nlat);
2094 free(current->nlon);
2095 free(current->slat);
2096 free(current->slon);
2097 free(current->di);
2098 free(current->dj);
2099 free(current->splat);
2100 free(current->splon);
2101 free(current->quasi);
2102 free(current->directionNumber);
2103 free(current->frequencyNumber);
2104 free(current->optimisationTime);
2105 free(current->leadTime);
2106 free(current->sensitiveAreaDomain);
2107 free(current->sensitiveAreaMethod);
2108 free(current->verifyingMonth);
2109 free(current->averagingPeriod);
2110 free(current->forecastMonth);
2111 free(current->referenceDate);
2112 free(current->climateDateFrom);
2113 free(current->climateDateTo);
2114 free(current->thresholdUnitsScale);
2115 free(current->thresholdIndicator);
2116 free(current->lowerThreshold);
2117 free(current->upperThreshold);
2118 }
2119
2120 if( thisFile == 0 )
2121 openFiles.files = (current->next);
2122 else
2123 previous->next = (current->next);
2124 free(current);
2125
2126 openFiles.count--;
2127 }
2128 }
2129
2130 fortint addRead(char* filename , fortint filename_len) {
2131 fortint thisFile;
2132
2133 if( DEBUG2 ) {
2134 char* pfile;
2135 copyName(&pfile, filename, filename_len);
2136 printf("PBG_addRead: add for reading filename = %s\n", pfile);
2137 free(pfile);
2138 }
2139
2140 thisFile = addFile(filename,filename_len,'r');
2141
2142 pbgindx( thisFile);
2143
2144 return thisFile;
2145 }
2146
2147 fortint addWrite(char* filename , fortint filename_len) {
2148 fortint thisFile;
2149
2150 if( DEBUG2 ) {
2151 char* pfile;
2152 copyName(&pfile, filename, filename_len);
2153 printf("PBG_addWrite: add for writing filename = %s\n", pfile);
2154 free(pfile);
2155 }
2156
2157 thisFile = addFile(filename,filename_len,'w');
2158
2159 return thisFile;
2160 }
2161
2162 void copyName(char* * pfile, char* filename, fortint filename_len) {
2163 /*
2164 // Copies FORTRAN filename to C string.
2165 */
2166 char * blankCharacterPointer;
2167 fortint space;
2168
2169 *pfile = (char *) malloc((size_t) (filename_len+1) );
2170 memcpy(*pfile, filename, (size_t) filename_len);
2171 space = (fortint) filename_len;
2172 (*pfile)[space] = '\0';
2173
2174 blankCharacterPointer = strchr(filename, ' ');
2175 if( blankCharacterPointer != NULL ) {
2176 space = (fortint) (blankCharacterPointer - filename);
2177 (*pfile)[space] = '\0';
2178 }
2179
2180 return;
2181 }
2182
2183 gribfile * latestFile(collection openFiles) {
2184 fortint index;
2185 gribfile * last = openFiles.files;
2186
2187 for( index = 1; index < openFiles.count; index++ )
2188 last = last->next;
2189
2190 return last;
2191
2192 }
2193
2194 gribfile * currentFile(collection openFiles, fortint thisFile) {
2195 fortint index;
2196 gribfile * current = openFiles.files;
2197
2198 for( index = 0; index < thisFile; index++ )
2199 current = current->next;
2200
2201 return current;
2202 }
2203
2204 #ifdef gribAPI
2205 long Julian_to_date(long jdate) {
2206 #else
2207 long julian_to_date(long jdate) {
2208 #endif
2209 long x,y,d,m,e;
2210 long day,month,year;
2211
2212 x = 4 * jdate - 6884477;
2213 y = (x / 146097) * 100;
2214 e = x % 146097;
2215 d = e / 4;
2216
2217 x = 4 * d + 3;
2218 y = (x / 1461) + y;
2219 e = x % 1461;
2220 d = e / 4 + 1;
2221
2222 x = 5 * d - 3;
2223 m = x / 153 + 1;
2224 e = x % 153;
2225 d = e / 5 + 1;
2226
2227 if( m < 11 )
2228 month = m + 2;
2229 else
2230 month = m - 10;
2231
2232 day = d;
2233 year = y + m / 11;
2234
2235 return year * 10000 + month * 100 + day;
2236 }
2237
2238 #ifdef gribAPI
2239 long date_to_Julian(long ddate) {
2240 #else
2241 long date_to_julian(long ddate) {
2242 #endif
2243 long m1,y1,a,b,c,d,j1;
2244 long month,day,year;
2245
2246 year = ddate / 10000;
2247 ddate %= 10000;
2248 month = ddate / 100;
2249 ddate %= 100;
2250 day = ddate;
2251
2252
2253 if (year < 100) year = year + 1900;
2254
2255 if (month > 2)
2256 {
2257 m1 = month - 3;
2258 y1 = year;
2259 }
2260 else
2261 {
2262 m1 = month + 9;
2263 y1 = year - 1;
2264 }
2265
2266 a = 146097*(y1/100)/4;
2267 d = y1 % 100;
2268 b = 1461*d/4;
2269 c = (153*m1+2)/5+day+1721119;
2270 j1 = a+b+c;
2271
2272 return(j1);
2273 }
2274
2275 fortint pbggeth012_(char*,fortint*,fortint,fortint,fortint);
2276 fortint soffset012_(unsigned char*,fortint*,fortint*,fortint*);
2277
2278 #define XFIND_ERROR -1
2279
2280 fortint pbgxfind_(fortint * grib1, char* filename, fortint filename_len)
2281 {
2282 unsigned char* buffer1 = (unsigned char*) grib1;
2283 static fortint* ibuffer2 = NULL;
2284 static fortint buffer2Length = 0;
2285 fortint status, totalLength, headerLength, numberOfGribs, loop;
2286 fortint is0, is1, is2;
2287 fortint js0, js1, js2;
2288
2289 numberOfGribs = pbgtotl_(filename,filename_len);
2290 if( numberOfGribs < 1 ) {
2291 if( DEBUG1 ) {
2292 char* pfile;
2293 copyName(&pfile, filename, filename_len);
2294 printf("PBXFIND: No GRIBs in file %s\n", pfile);
2295 free(pfile);
2296 }
2297 return XFIND_ERROR;
2298 }
2299
2300 status = soffset012_(buffer1, &is0, &is1, &is2);
2301 if( status ) return XFIND_ERROR;
2302
2303 totalLength = THREEBYTELONG(buffer1+is0+4);
2304 headerLength = (is1 - is0);
2305 headerLength += THREEBYTELONG(buffer1+is1);
2306 if( is2 != 0 ) headerLength += THREEBYTELONG(buffer1+is2);
2307
2308 for( loop = 0; loop < numberOfGribs; loop++) {
2309 fortint next, totalLength2, headerLength2;
2310 unsigned char * buffer2;
2311
2312 next = loop + 1;
2313
2314 if( ibuffer2 != NULL ) {
2315 if( buffer2Length < headerLength ) {
2316 ibuffer2 = (fortint*) realloc(ibuffer2,headerLength);
2317 if( ibuffer2 == NULL ) {
2318 printf("pbgxfind_: realloc failed\n");
2319 return XFIND_ERROR;
2320 }
2321 buffer2Length = headerLength;
2322 }
2323 }
2324 else {
2325 ibuffer2 = (fortint*) malloc(headerLength);
2326 if( ibuffer2 == NULL ) {
2327 printf("pbgxfind_: malloc failed\n");
2328 return XFIND_ERROR;
2329 }
2330 buffer2Length = headerLength;
2331 }
2332
2333 totalLength2 =
2334 pbggeth012_(filename,ibuffer2,buffer2Length,next,filename_len);
2335
2336 if( totalLength2 > 1 ) {
2337
2338 buffer2 = (unsigned char*) ibuffer2;
2339 status = soffset012_(buffer2, &js0, &js1, &js2);
2340 if( status ) return XFIND_ERROR;
2341
2342 if( totalLength == totalLength2 ) {
2343 headerLength2 = (js1 - js0);
2344 headerLength2 += THREEBYTELONG(buffer2+js1);
2345 if( js2 != 0 ) headerLength2 += THREEBYTELONG(buffer2+js2);
2346
2347 if( headerLength == headerLength2 ) {
2348 if( (memcmp((buffer1+is0),(buffer2+js0),(size_t)headerLength) == 0) )
2349 return next;
2350 }
2351 }
2352 }
2353
2354 }
2355
2356 return XFIND_ERROR;
2357 }
2358
2359 fortint pbgxfind(fortint * grib1, char* filename, fortint filename_len) {
2360 return pbgxfind_(grib1,filename,filename_len);
2361 }
2362
2363 #define GRIB_TOO_SMALL -4
2364
2365 fortint pbggeth012_(char* filename, fortint* buffer, fortint bufflen,
2366 fortint n, fortint filename_len) {
2367 /*
2368 // Gets section 0,1 and 2 for the nth GRIB product.
2369 */
2370 fortint index = n - 1;
2371 fortint length, status;
2372 off_t offset;
2373 fortint thisFile;
2374 FILE * fp;
2375 gribfile * file;
2376
2377 if( DEBUG1 ) {
2378 char* pfile;
2379 copyName(&pfile, filename, filename_len);
2380 printf("pbggeth012: getting GRIB number %ld in file %s\n",
2381 (index+1), pfile);
2382 free(pfile);
2383 }
2384
2385 if( index >= 0 ) {
2386 thisFile = pbginitInput(filename,filename_len);
2387 file = currentFile(openFiles, thisFile);
2388 if( index < file->count ) {
2389
2390 length = file->length[index];
2391
2392 if( DEBUG1 ) printf("pbggeth012: length of GRIB number %ld = %ld\n",
2393 (index+1), length);
2394
2395 if( length < bufflen ) return GRIB_TOO_SMALL;
2396
2397 offset = file->offset[index];
2398
2399 if( DEBUG1 ) printf("pbggeth012: offset of GRIB number %ld = %lld\n",
2400 (index+1), offset);
2401
2402 fp = file->fp;
2403 #ifdef FOPEN64
2404 if( fseek( fp, offset, 0) ) {
2405 #else
2406 if( fseek( fp, offset, 0) ) {
2407 #endif
2408 perror("pbggeth012: error in fseek");
2409 return FSEEK_ERROR;
2410 }
2411
2412 status = fread(buffer, 1, bufflen, fp);
2413 if( status != bufflen ) {
2414 fprintf(stderr,"pbggeth012: error in fread\n");
2415 return FREAD_ERROR;
2416 }
2417
2418 return length;
2419 }
2420 }
2421
2422 return INDEX_TOO_BIG;
2423
2424 }
2425
2426 #define ERROR(a,b) {perror(a);return b;}
2427 #define GRIB 0x47524942
2428 #define len3oct(p) ((((long)*(p))<<16) + (((long)*(p+1))<<8) + (long)*(p+2))
2429 #define BIT1 0x80
2430 #define BIT2 0x40
2431 #define BIT3 0x20
2432 #define BIT4 0x10
2433 #define BIT5 0x08
2434 #define BIT6 0x04
2435 #define BIT7 0x02
2436 #define BIT8 0x01
2437
2438 static int grab(unsigned char * , unsigned char * , long ,long ,long * );
2439
2440 fortint soffset012_(
2441 unsigned char * buffer,
2442 fortint* is0,
2443 fortint* is1,
2444 fortint* is2 ) {
2445 long s0, s1, s2, edition;
2446 int large = 0;
2447 int found = 0;
2448 int code = 0;
2449 long bytes_read = 0, advance;
2450 unsigned char p, edit_num, flag23;
2451 unsigned char size[3];
2452 int section0 = 8, section1, section2;
2453 long total;
2454
2455 /*
2456 // Read bytes until "GRIB" found
2457 */
2458 do {
2459 if( grab(buffer, &p, 1, 1, &bytes_read) != 0) return 1;
2460 code = ( (code << 8) + p ) & 0xFFFFFFFF;
2461 if (code == GRIB ) found = 1;
2462 } while ( ! found );
2463 s0 = bytes_read - 4;
2464 bytes_read = 4;
2465 /*
2466 // Now find out which edition of GRIB is present (default is 1)
2467 */
2468 edition = 1;
2469 s1 = s0 + 8;
2470 if( (*(buffer+21-s0) == '\0') && (*(buffer+22-s0) == '\0') ) {
2471 edition = -1; /* GRIB edition -1 */
2472 s1 = s0;
2473 section1 = 20;
2474 section0 = 4;
2475 }
2476 else {
2477 if( grab(buffer, size, 3, 1, &bytes_read) != 0) return 1;
2478 total = len3oct(size);
2479 if( total == 24 ) {
2480 /*
2481 // Move past the edition number
2482 */
2483 if( grab(buffer, &edit_num, 1, 1, &bytes_read) != 0) return 1;
2484 edition = 0; /* GRIB edition 0 */
2485 section1 = 24;
2486 s1 = s0 + 4;
2487 section0 = 4;
2488 }
2489 }
2490
2491 if( edition == 1 ) {
2492 /*
2493 // See if it is an extra large (wave) product
2494 */
2495 if( total > 0x800000 ) {
2496 total = (total&0x7fffff) * 120;
2497 large = 1;
2498 }
2499 /*
2500 // Move past the edition number
2501 */
2502 if( grab(buffer, &edit_num, 1, 1, &bytes_read) != 0) return 1;
2503 /*
2504 // Read length of section 1
2505 */
2506 if( grab(buffer, size, 3, 1, &bytes_read) != 0) return 1;
2507 section1 = len3oct(size);
2508 }
2509 /*
2510 // Now figure out if section 2 is present
2511 */
2512 advance = 4;
2513 bytes_read += advance;
2514 if( grab(buffer, &flag23, 1, 1, &bytes_read) != 0) return 1;
2515 section2 = flag23 & BIT1;
2516
2517 /*
2518 // Advance to end of section 1
2519 */
2520 advance = section1 - (bytes_read - section0);
2521 bytes_read += advance;
2522 /*
2523 // Read section 2 length if it is given
2524 */
2525 if( section2 ) {
2526 s2 = s0 + bytes_read;
2527 if( grab(buffer, size, 3, 1, &bytes_read) != 0) return 1;
2528 section2 = len3oct(size);
2529 advance = section2 - (bytes_read - section0 - section1);
2530 bytes_read += advance;
2531 }
2532 else {
2533 section2 = 0;
2534 s2 = 0;
2535 }
2536 /*
2537 // Success!
2538 */
2539 *is0 = (fortint) s0;
2540 *is1 = (fortint) s1;
2541 *is2 = (fortint) s2;
2542 return 0;
2543 }
2544
2545 static int grab(unsigned char * buffer, unsigned char * where, long size,long cnt,long * num_bytes_read)
2546 {
2547 long number = size*cnt;
2548
2549 memcpy(where, (buffer+(*num_bytes_read)), number);
2550 *num_bytes_read += number;
2551
2552 return 0;
2553 }
2554