1 /* Copyright (C) 1988-2005 by Brian Doty and the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 #ifdef HAVE_CONFIG_H
7 #include <config.h>
8
9 /* If autoconfed, only include malloc.h when it's presen */
10 #ifdef HAVE_MALLOC_H
11 #include <malloc.h>
12 #endif
13
14 #else /* undef HAVE_CONFIG_H */
15
16 #include <malloc.h>
17
18 #endif /* HAVE_CONFIG_H */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <unistd.h>
23 #include <math.h>
24 #include <string.h>
25
26 #include "grads.h"
27 #include "gagmap.h"
28
29 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
30 extern struct gamfcmn mfcmn;
31 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
32
33 /* Routine to scan a grib file and output a grib index (map) file. */
34
35 /*---- 960829 Modification to run outside of gribmap for LATS
36
37 Mike Fiorino to run outside gribmap
38
39 ----*/
40
41 /* Values output into the grib map file:
42
43 Header:
44
45 hipnt info: 0 - version number (1)
46 1 - number of times in file
47 2 - number of records per time
48 3 - Grid type
49 255 - user defined grid. descriptor
50 describes grid exactly; one record
51 per grid.
52 29 - Predefined grid set 29 and 30.
53 Two records per grid.
54
55 hfpnt info: None
56
57 Info:
58
59 intpnt info (for each mapped grib record) :
60 0 - position of start of data in file
61 1 - position of start of bit map in file
62 2 - number of bits per data element
63
64 fltpnt info :
65 0 - decimal scale factor for this record
66 1 - binary scale factor
67 2 - reference value
68
69 */
70
71
72 struct dt rtime; /* time based on input map */
73 struct dt ftime; /* time based on dd file */
74
75 struct gaindx *rindx;
76
77
78 /* internal prototypes */
79
80 int wtgmap(void) ;
81 void putint(int, unsigned char *,int *) ;
82
83 /* global to this file only */
84
85 int nelements=3;
86 char mrec[512];
87 static off_t flen;
88 /*---------------------------------------
89
90 main routine
91
92 ---------------------------------------*/
93
gribmap(void)94 int gribmap (void) {
95
96
97 /*---
98 declare internal vars
99 ---*/
100
101 off_t filesize;
102 unsigned char vermap;
103 int idum;
104 float fdum;
105 unsigned char urec[4];
106 int nt,i,j;
107 int update_map=0;
108
109 float map_grvals[8];
110 int rt,rtoff,toff,nxytime;
111 int didmatch=0,rcgr;
112
113 /*---
114 initialize this global variable here...
115 ---*/
116
117 mfile=NULL;
118
119 /*---
120 open the dd file
121 ---*/
122
123 if (ifile==NULL) {
124 printf ("\n");
125 cnt = nxtcmd (cmd,"Enter name of Data Descriptor file: ");
126 if (cnt==0) return(1);
127 getwrd(rec,cmd,250);
128 ifile = rec;
129 }
130
131 pfi = getpfi();
132 if (pfi==NULL) {
133 printf ("Memory allocation error \n");
134 return(1);
135 }
136
137 pindx = (struct gaindx *)malloc(sizeof(struct gaindx));
138 if (pindx==NULL) {
139 printf ("Memory allocation error for pindx\n");
140 return(1);
141 }
142
143 rindx = (struct gaindx *)malloc(sizeof(struct gaindx));
144 if (rindx==NULL) {
145 printf ("Memory allocation error for rindx\n");
146 return(1);
147 }
148
149 /* ---------------- open the dd file ------------------- */
150
151
152 rc = gaddes (ifile, pfi, 0);
153 if (rc) {
154 if(! quiet) printf ("File name is: %s\n",ifile);
155 return(1);
156 }
157 if (pfi->idxflg!=1) {
158 printf ("Data Descriptor file is not for GRIB data\n");
159 return(1);
160 }
161
162 /* save the initial time from the data descriptor file for the tau0 option and the map file*/
163
164 btimdd.yr = *(pfi->abvals[3]);
165 btimdd.mo = *(pfi->abvals[3]+1);
166 btimdd.dy = *(pfi->abvals[3]+2);
167 btimdd.hr = *(pfi->abvals[3]+3);
168 btimdd.mn = *(pfi->abvals[3]+4);
169 if (no_min) btimdd.mn = 0;
170
171 if(diag) printf(" btim %i %i %i %i %i\n",btimdd.yr,btimdd.mo,btimdd.dy,btimdd.hr,btimdd.mn);
172
173
174 /* Set up for this grid type */
175
176 if (pfi->grbgrd<-900 || pfi->grbgrd==255) {
177 nrec = 1;
178 gtype[0] = 255;
179 } else if (pfi->grbgrd>-1 && pfi->ppflag) {
180 gtype[0] = pfi->grbgrd;
181 nrec=1;
182 } else if (pfi->grbgrd==29) {
183 nrec = 2;
184 gtype[0] = 29;
185 gtype[1] = 30;
186 if (pfi->dnum[0]!=144 || pfi->dnum[1]!=73 ||
187 pfi->linear[0]!=1 || pfi->linear[1]!=1 ||
188 *(pfi->grvals[0])!= 2.5 || *(pfi->grvals[0]+1) != -2.5 ||
189 *(pfi->grvals[1])!= 2.5 || *(pfi->grvals[1]+1) != -92.5 ) {
190 printf ("Error in grid specification for GRIB grid type 29/30.\n");
191 printf (" Grid scaling must indicate a 2.5 x 2.5 grid\n");
192 printf (" Grid size must be 144 x 73\n");
193 printf (" Grid must go from 0 to 357.5 and -90 to 90\n");
194 return(1);
195 }
196 } else {
197 nrec = 1;
198 gtype[0] = pfi->grbgrd;
199 }
200
201 /* Set index stuff */
202
203 pindx->type = ver; /* gribmap version */
204 pindx->hinum = 4;
205 pindx->hfnum = 0;
206 pindx->intnum = nrec * nelements * pfi->trecs * pfi->dnum[3];
207 pindx->fltnum = nrec * nelements * pfi->trecs * pfi->dnum[3];
208 pindx->hipnt = (int *)malloc(sizeof(int)*pindx->hinum);
209 pindx->intpnt = (int *)malloc(sizeof(int)*pindx->intnum);
210 pindx->fltpnt = (float *)malloc(sizeof(float)*pindx->fltnum);
211
212 if (pindx->hipnt==NULL || pindx->intpnt==NULL || pindx->fltpnt==NULL) {
213 printf ("memory allocation error\n");
214 return(1);
215 }
216 for (i=0; i<pindx->intnum; i++) *(pindx->intpnt+i) = -999;
217 for (i=0; i<pindx->fltnum; i++) *(pindx->fltpnt+i) = -999; /* try initializing */
218 *(pindx->hipnt) = ver;
219 *(pindx->hipnt+1) = pfi->dnum[3];
220 *(pindx->hipnt+2) = pfi->trecs;
221 *(pindx->hipnt+3) = pfi->grbgrd;
222 if (pfi->grbgrd<-900) *(pindx->hipnt+3) = 255;
223
224 if(! quiet) printf ("\ngribmap: Scanning binary GRIB file(s):\n");
225
226
227 /* If there are multiple files in a time series, we need to
228 loop through them based on the info in pfi */
229
230 tcur = 0;
231 gfile = NULL;
232
233 if(pfi->tmplat && update ) {
234
235 if(!quiet) {
236 printf("UUUUU: Updating the gribmap for the TEMPLATED GRIB data described by the:\n");
237 printf("UUUUU:\n");
238 printf("UUUUU: %s\n",pfi->dnam);
239 printf("UUUUU:\n");
240 printf("UUUUU: GrADS data descriptor (.ctl) file....\n\n");
241 }
242
243 if (nrec != 1) {
244 printf ("Attempting to update a special case GRIB file nrec = %d\n",nrec);
245 printf ("No can do, bye \n");
246 return(42);
247 }
248
249 mfile = fopen(pfi->mnam,"rb+");
250
251 if (mfile==NULL) {
252
253 printf ("Could not open index file: %s for updating\n",pfi->mnam);
254 printf ("Create the map........\n");
255 update=0;
256
257 } else {
258
259 fseeko(mfile,0L,SEEK_END);
260 filesize=ftello(mfile);
261
262 if(filesize>0) {
263
264 fseeko(mfile,1,0);
265 fread(&vermap,sizeof(unsigned char),1,mfile);
266 if(diag) printf(" vermap =%d\n",vermap);
267
268 if(vermap == 2) {
269
270 fseeko(mfile,2,0);
271 fread(mrec,sizeof(unsigned char),4,mfile);
272 rindx->hinum=gagby(mrec,0,4);
273 if(diag) printf(" hinum = %d\n",rindx->hinum);
274
275 fread(mrec,sizeof(unsigned char),4,mfile);
276 rindx->hfnum=gagby(mrec,0,4);
277 if(diag) printf(" hfnum = %d\n",rindx->hfnum);
278
279 fread(mrec,sizeof(unsigned char),4,mfile);
280 rindx->intnum=gagby(mrec,0,4);
281 if(diag) printf(" intnum = %d\n",rindx->intnum);
282
283 fread(mrec,sizeof(unsigned char),4,mfile);
284 rindx->fltnum=gagby(mrec,0,4);
285 if(diag) printf(" fltnum = %d\n",rindx->fltnum);
286
287 rindx->hipnt = NULL;
288 rindx->hfpnt = NULL;
289 rindx->intpnt = NULL;
290 rindx->fltpnt = NULL;
291
292 /* skip the begining time struct info */
293
294 fread(mrec,sizeof(unsigned char),7,mfile);
295
296 if (rindx->hinum>0) {
297 rindx->hipnt = (int *)malloc(sizeof(int)*rindx->hinum);
298 if (rindx->hipnt==NULL) goto err8;
299
300 for(i=0;i<rindx->hinum;i++) {
301 fread(mrec,sizeof(unsigned char),4,mfile);
302 idum=gagby(mrec,0,4);
303 if(gagbb(mrec,0,1)) idum=-idum;
304 *(rindx->hipnt+i)=idum;
305 if(diag) printf(" 1 i = %d rindx->hipnt = %d\n",i,idum);
306 }
307
308 }
309
310 if (rindx->hfnum>0) {
311 rindx->hfpnt = (float *)malloc(sizeof(float)*rindx->hfnum);
312 if (rindx->hfpnt==NULL) goto err8;
313 fread (rindx->hfpnt,sizeof(float),rindx->hfnum,mfile);
314 }
315
316 if (rindx->intnum>0) {
317 rindx->intpnt = (int *)malloc(sizeof(int)*rindx->intnum);
318 if (rindx->intpnt==NULL) goto err8;
319 for(i=0;i<rindx->intnum;i++) {
320 fread(mrec,sizeof(unsigned char),4,mfile);
321 idum=gagby(mrec,0,4);
322 if(gagbb(mrec,0,1)) idum=-gagbb(mrec,1,31);
323 *(rindx->intpnt+i)=idum;
324 if(diag) printf(" 2 i = %d rindx->intpnt = %d\n",i,idum);
325 }
326
327 }
328
329 if (rindx->fltnum>0) {
330 rindx->fltpnt = (float *)malloc(sizeof(float)*rindx->fltnum);
331 if (rindx->fltpnt==NULL) goto err8;
332
333 for(i=0;i<rindx->fltnum;i++) {
334 fread(urec,sizeof(unsigned char),4,mfile);
335 fdum=ibm2flt(urec);
336 *(rindx->fltpnt+i)=fdum;
337 if(diag) printf(" 3 i = %d rindx->fltpnt = %g\n",i,fdum);
338 }
339
340 }
341
342 /* read in the grid to absolute time factors at the end of the regular map info */
343
344 for(i=0;i<8;i++) {
345 fread(urec,sizeof(unsigned char),4,mfile);
346 fdum=ibm2flt(urec);
347 *(map_grvals+i)=fdum;
348 if(diag) printf(" grvals i = %d map_grvals = %g\n",i,fdum);
349 }
350
351 }
352
353
354 nt=*(rindx->hipnt+1);
355 nxytime=pfi->trecs*nrec*nelements;
356
357 for(i=1;i<=nt;i++) {
358 gr2t(map_grvals,(float)i,&rtime);
359 rt=(int)(t2gr(map_grvals,&rtime)+0.01);
360 rtoff=(rt-1)*nxytime;
361 if(diag) printf("rrrr %d %d %d %d %d rt = %d rtoff = %d\n",
362 rtime.yr,rtime.mo,rtime.dy,rtime.hr,rtime.mn,rt,rtoff);
363 }
364
365
366 } else {
367
368 printf("WARNING: Map file for updating has 0 size\n");
369 printf(" will go ahead and do the map...\n\n");
370 update=0;
371
372 }
373
374 }
375
376 /*---
377 close the map file now that we've gotten all the information we need; including the map
378 ---*/
379
380 fclose(mfile);
381 mfile=NULL;
382
383 }
384
385
386 while (1) {
387
388 if (pfi->tmplat) {
389
390 if (gfile!=NULL) {
391 fclose(gfile);
392 gfile=NULL;
393 }
394
395 told = tcur;
396
397 if (tcur==0) {
398
399 tcur=1;
400
401 if(update) {
402
403 tcur=1;
404 told=tcur;
405
406 while (told==*(pfi->fnums+tcur-1) && tcur<=pfi->dnum[3]) {
407
408 gr2t(pfi->grvals[3],(float)tcur,&dtim);
409 rt=t2gr(map_grvals,&dtim);
410 rt=(int)(t2gr(map_grvals,&dtim)+0.01);
411 rtoff=(rt-1)*nxytime;
412 toff=(tcur-1)*nxytime;
413
414 /*---
415 check if a match for first time in the file
416 ---*/
417
418 update_map=0;
419 if( rt > 0 && rt <= nt ) update_map=1;
420
421 if(diag) {
422 printf("\nuuu0 told= %d tcur= %d dtim: %d %d %d %d %d\n",
423 told,tcur,dtim.yr,dtim.mo,dtim.dy,dtim.hr,dtim.mn);
424 printf("uuu0 rt = %d rtoff = %d toff = %d \n",rt,rtoff,toff);
425 printf("uuu0 told= %d tcur= %d update_map = %d\n",told,tcur,update_map);
426 }
427
428 /*---------------------
429 check if update not available in MIDDLE of a time chunk (told -> tcur)
430 ------------------*/
431
432 if(update_map==0 && (tcur != told) ) {
433 if(diag) printf("\nmmmmmmm 000 check told = %d tcur=%d\n\n",told,tcur);
434 tcur=told;
435 update=0;
436 break;
437 }
438
439 if(update_map) {
440 for(j=0;j<nxytime;j++) {
441 *(pindx->intpnt+toff+j)=*(rindx->intpnt+rtoff+j);
442 *(pindx->fltpnt+toff+j)=*(rindx->fltpnt+rtoff+j);
443 }
444 }
445
446 tcur++;
447
448 }
449
450 tcur=1;
451 }
452
453 } else {
454
455 while (told==*(pfi->fnums+tcur-1) && tcur<=pfi->dnum[3]) {
456
457 if(update) {
458
459 gr2t(pfi->grvals[3],(float)tcur,&dtim);
460 rt=(int)(t2gr(map_grvals,&dtim)+0.01);
461 rtoff=(rt-1)*nxytime;
462 toff=(tcur-1)*nxytime;
463
464 update_map=0;
465 if( rt > 0 && rt <= nt ) update_map=1;
466
467 if(diag) {
468 printf("\nuuu1 told= %d tcur= %d dtim: %d %d %d %d %d\n",told,tcur,dtim.yr,dtim.mo,dtim.dy,dtim.hr,dtim.mn);
469 printf("uuu1 rt = %d rtoff = %d toff = %d\n",rt,rtoff,toff);
470 printf("uuu1 update_map = %d\n",update_map);
471 }
472
473 /*---------------------
474 check if update not available in MIDDLE of a time chunk (told -> tcur)
475 ------------------*/
476
477 if(update_map==0 && (tcur != told) ) {
478 if(diag) printf("\nmmmmmmm check told = %d tcur=%d\n\n",told,tcur);
479 tcur=told;
480 update=0;
481 break;
482 }
483
484 if(update_map) {
485 for(j=0;j<nxytime;j++) {
486 *(pindx->intpnt+toff+j)=*(rindx->intpnt+rtoff+j);
487 *(pindx->fltpnt+toff+j)=*(rindx->fltpnt+rtoff+j);
488 }
489 }
490
491 }
492
493 tcur++;
494
495 }
496
497 /*---
498 tcur has been updated; check if we are at the end of the file...
499 ---*/
500
501 if(update) {
502 gr2t(pfi->grvals[3],(float)tcur,&dtim);
503 rt=(int)(t2gr(map_grvals,&dtim)+0.01);
504 update_map=0;
505 if( rt > 0 && rt <= nt ) update_map=1;
506 }
507
508 }
509
510 if (tcur>pfi->dnum[3]) break;
511 if (tcur != *(pfi->fnums+tcur-1)) {
512 printf ("Logic error in file template structure\n");
513 return(1);
514 }
515
516 gr2t(pfi->grvals[3], (float)tcur, &dtim);
517 gr2t(pfi->grvals[3], 1.0, &dtimi);
518
519 if(diag) printf("11111 Before open test: tcur= %d update= %d update_map=%d\n",
520 tcur,update,update_map);
521
522 if(update_map==0) {
523 ch = gafndt(pfi->name, &dtim, &dtimi, pfi->abvals[3], pfi->pchsub1, tcur);
524 if (ch==NULL) {
525 printf ("Memory allocation error\n");
526 return(1);
527 }
528 }
529
530 } else {
531
532 /*------- non template ----------*/
533
534 ch = pfi->name;
535
536 }
537
538 /*---
539 Open this GRIB file and position to start of first record
540 ---*/
541
542
543 if(diag) printf("22222 After open test: update= %d update_mape= %d\n",update,update_map);
544
545 if(update_map==0) {
546
547 if(!quiet) printf("gribmap: Opening GRIB file %s\n",ch);
548 gfile = fopen(ch,"rb");
549
550 if (gfile==NULL) {
551 if (pfi->tmplat) {
552 printf ("Warning: Could not open GRIB file: %s\n",ch);
553 continue;
554 } else {
555 printf ("Could not open GRIB file. File name is:\n");
556 printf (" %s\n",ch);
557 return(1);
558 }
559 }
560
561 if (pfi->tmplat) free(ch);
562
563 /* Get file size */
564
565 fseeko(gfile,0L,2);
566 flen = ftello(gfile);
567
568 /*---
569 Set up to skip appropriate amount and position
570 ---*/
571
572 if (skip > -1) fpos = skip;
573 else {
574 fseeko (gfile,0,0);
575 rc = fread (rec,1,100,gfile);
576 if (rc<100) {
577 printf ("I/O Error reading header\n");
578 return(1);
579 }
580 len = gagby(rec,88,4);
581 fpos = len*2 + 100;
582 }
583
584
585 /*--- main loop
586
587 1) read a grib record (gribhdr)
588 2) compare to each 2-d variable in the 5-D data volume
589 defined by the dd for a match (gribrec)
590
591 ---- main loop */
592
593 irec=1;
594 while (1) {
595 rc=gribhdr(&ghdr);
596 if (rc) break;
597 rcgr=gribrec(&ghdr,pfi,pindx);
598 if(rcgr==0) didmatch=1;
599 if(rcgr>=100) didmatch=rcgr;
600 irec++;
601 }
602
603 /*---
604 see how we did
605 ---*/
606
607 if (rc==50) {
608 printf ("I/O Error reading GRIB file\n");
609 printf ("Possible cause: premature EOF\n");
610 break;
611 }
612 if (rc>1 && rc!=98) {
613 printf ("GRIB file format error \n");
614 printf ("rc = %i\n",rc);
615 return(rc);
616 }
617
618 /*---
619 break out if not templating
620 ---*/
621
622 if (!pfi->tmplat) break;
623
624 /*---
625 end for update option
626 ---*/
627
628 }
629
630 /*---
631 return for another file if templating
632 ---*/
633
634 }
635 /*---
636 ALL DONE the last(only) GRIB file
637 ---*/
638
639 if(!quiet) printf ("gribmap: Reached EOF\n");
640
641 /*---
642 check if file closed already for case where template was set,
643 but it was not templated and the template code above closed it.
644 ---*/
645
646 if (gfile!=NULL) {
647 fclose (gfile);
648 gfile=NULL;
649 }
650
651 /*---
652 open the map file
653 ---*/
654
655 if(mfile!=NULL) {
656 if(diag) printf("-------------- file already open.....\n");
657 } else {
658 if(write_map) {
659 mfile = fopen(pfi->mnam,"wb");
660 } else {
661 printf("\nWARNING!!!! NOT writing the map !!!\n\n");
662 }
663 }
664
665 if (mfile==NULL && write_map) {
666 printf ("Could not open index file: %s\n",pfi->mnam);
667 return(1);
668 } else if(write_map) {
669 if(!quiet) printf("gribmap: Writing the map...\n\n");
670 }
671
672 /*---
673 output the map depending on version #
674 ---*/
675
676 if(ver==1) {
677 fwrite (pindx,sizeof(struct gaindx),1,mfile);
678 if (pindx->hinum>0) fwrite(pindx->hipnt,sizeof(int),pindx->hinum,mfile);
679 if (pindx->hfnum>0) fwrite(pindx->hfpnt,sizeof(float),pindx->hfnum,mfile);
680 if (pindx->intnum>0) fwrite(pindx->intpnt,sizeof(int),pindx->intnum,mfile);
681 if (pindx->fltnum>0) fwrite(pindx->fltpnt,sizeof(float),pindx->fltnum,mfile);
682 fclose (mfile);
683 }
684
685 /*---
686 new machine independent version
687 ---*/
688
689 if(ver==2) {
690
691 rc=wtgmap();
692 if(rc == 601) {
693 printf("\nGRIBMAP failed because of overflow in float -> IBM float conversion!!!\n");
694 printf("Contact developer fiorino@llnl.gov...\n\n");
695 return(601);
696 }
697 fclose (mfile);
698 }
699
700 return(didmatch);
701
702 /*--- errors ---*/
703
704 err8:
705 printf("Open Error: Memory allocation Error\n");
706 return(99);
707
708 }
709
710 /* Read a GRIB header, and get needed info out of it. */
711
gribhdr(struct grhdr * ghdr)712 int gribhdr (struct grhdr *ghdr) {
713 struct dt atim;
714 char rec[50000],*ch,*gds;
715 int i,len ,rc,sign,mant;
716 off_t cpos;
717
718 if (fpos+50>=flen) return(1);
719
720
721 /* look for data between records BY DEFAULT */
722
723 i = 0;
724 fpos += i;
725 rc = fseek(gfile,fpos,0);
726 if (rc) return(50);
727 ch=&rec[0];
728 rc = fread(ch,sizeof(char),4,gfile);
729
730 while ( (fpos < flen-4) && (i < scanlim) && !( *ch=='G' && *(ch+1)=='R' &&
731 *(ch+2)=='I' && *(ch+3)=='B' ) ) {
732 fpos++;
733 i++;
734 rc = fseeko(gfile,fpos,0);
735 if (rc) return(50);
736 rc = fread(ch,sizeof(char),4,gfile);
737 if (rc<4) return(50);
738 }
739
740 if (i == scanlim ) {
741 printf ("GRIB header not found in scanning between records\n");
742 printf ("->%c%c%c%c<-\n",*rec,*(rec+1),*(rec+2),*(rec+3));
743 if(scaneof==1 || scanEOF==1) {
744 if(scaneof) {
745 printf("...... punching out -e option set ......\n");
746 return(98);
747 }
748 if(scanEOF) {
749 printf("...... junk found, continue -E option set ......\n");
750 return(0);
751 }
752 } else {
753 return(52);
754 }
755 } else if (fpos == flen -4) {
756 if(scaneof==1 || scanEOF==1) {
757 if(scaneof) {
758 printf("...... punching out -e option set flen-4 check ......\n");
759 return(98);
760 }
761 if(scanEOF) {
762 printf("...... junk found, continue -E option set flen-4 check ......\n");
763 return(0);
764 }
765 } else {
766 return (53);
767 }
768 }
769
770 /* SUCCESS redo the initial read */
771
772 rc = fseek(gfile,fpos,0);
773 if (rc) return(50);
774 rc = fread(rec,1,8,gfile);
775 if (rc<8){
776 if(fpos+8 >= flen) return(61);
777 else return(62);
778 }
779
780
781
782 cpos = fpos;
783 ghdr->vers = gagby(rec,7,1);
784 if (ghdr->vers>1) {
785 printf ("File is not GRIB version 0 or 1, 0 or 1 is required. \n");
786 printf ("Version number is %i\n",ghdr->vers);
787 if(scaneof) {
788 printf("...... punching out -e option set junk at end ......\n");
789 return(98);
790 } else {
791 return (99);
792 }
793 }
794
795 if (ghdr->vers==0) {
796 cpos += 4;
797 rc = fseek(gfile,cpos,0);
798 if (rc) return(50);
799 } else {
800 ghdr->len = gagby(rec,4,3);
801 cpos = cpos + 8;
802 rc = fseeko(gfile,cpos,0);
803 if (rc) return(50);
804 }
805
806 /* Get PDS length, read rest of PDS */
807
808 rc = fread(rec,1,3,gfile);
809 if (rc<3) return(50);
810 len = gagby(rec,0,3);
811 ghdr->pdslen = len;
812 cpos = cpos + len;
813 rc = fread(rec+3,1,len-3,gfile);
814 if (rc<len-3) return(50);
815
816 /* Get info from PDS */
817
818 ghdr->id = gagby(rec,6,1);
819 ghdr->gdsflg = gagbb(rec+7,0,1);
820 ghdr->bmsflg = gagbb(rec+7,1,1);
821 ghdr->parm = gagby(rec,8,1);
822 ghdr->ltyp = gagby(rec,9,1);
823 ghdr->level = gagby(rec,10,2);
824 ghdr->l1 = gagby(rec,10,1);
825 ghdr->l2 = gagby(rec,11,1);
826 if (mpiflg) { /* Use initial time from the descriptor for mpiflg */
827 ghdr->btim.yr = *(pfi->abvals[3]);
828 ghdr->btim.mo = *(pfi->abvals[3]+1);
829 ghdr->btim.dy = *(pfi->abvals[3]+2);
830 ghdr->btim.hr = *(pfi->abvals[3]+3);
831 ghdr->btim.mn = *(pfi->abvals[3]+4);
832 if (no_min) ghdr->btim.mn = 0;
833 } else {
834 ghdr->btim.yr = gagby(rec,12,1);
835 ghdr->btim.mo = gagby(rec,13,1);
836 ghdr->btim.dy = gagby(rec,14,1);
837 ghdr->btim.hr = gagby(rec,15,1);
838 ghdr->btim.mn = gagby(rec,16,1);
839 if (no_min) ghdr->btim.mn = 0;
840 }
841 if (ghdr->btim.hr>23) ghdr->btim.hr = 0; /* Special for NCAR */
842 if (len>24) {
843 ghdr->cent = gagby(rec,24,1);
844 ghdr->btim.yr = ghdr->btim.yr + (ghdr->cent-1)*100;
845 } else {
846 ghdr->cent = -999;
847 if (!(mpiflg) || !(mfcmn.fullyear)) {
848 if (ghdr->btim.yr>49) ghdr->btim.yr += 1900;
849 if (ghdr->btim.yr<50) ghdr->btim.yr += 2000;
850 }
851 }
852 ghdr->ftu = gagby(rec,17,1);
853 ghdr->tri = gagby(rec,20,1);
854 if (ghdr->tri==10) {
855 ghdr->p1 = gagby(rec,18,2);
856 ghdr->p2 = 0;
857 } else {
858 ghdr->p1 = gagby(rec,18,1);
859 ghdr->p2 = gagby(rec,19,1);
860 }
861
862 ghdr->fcstt = ghdr->p1;
863 if (ghdr->tri>1 && ghdr->tri<6) ghdr->fcstt=ghdr->p2;
864 if ((tauave) && ghdr->tri==3) ghdr->fcstt=ghdr->p1;
865 atim.yr=0; atim.mo=0; atim.dy=0; atim.hr=0; atim.mn=0;
866 if (ghdr->ftu==0) atim.mn = ghdr->fcstt;
867 else if (ghdr->ftu==1) atim.hr = ghdr->fcstt;
868 else if (ghdr->ftu==10) atim.hr = ghdr->fcstt*3; /* added 3Hr incr */
869 else if (ghdr->ftu==11) atim.hr = ghdr->fcstt*6; /* added 6Hr incr */
870 else if (ghdr->ftu==12) atim.hr = ghdr->fcstt*12; /* Special NCEP 12Hr incr */
871 else if (ghdr->ftu==2) atim.dy = ghdr->fcstt;
872 else if (ghdr->ftu==3) atim.mo = ghdr->fcstt;
873 else if (ghdr->ftu==4) atim.yr = ghdr->fcstt;
874 else ghdr->fcstt = -999;
875 /* mf */
876 /* if notau != 0 then FORCE the valid DTG to be the base DTG */
877 if(notau) ghdr->fcstt = -999 ;
878 /* mf */
879
880 /* add the forecast time to the time of this grib field */
881
882 if (ghdr->fcstt>-900) {
883 timadd(&(ghdr->btim),&atim);
884 ghdr->dtim.yr = atim.yr;
885 ghdr->dtim.mo = atim.mo;
886 ghdr->dtim.dy = atim.dy;
887 ghdr->dtim.hr = atim.hr;
888 ghdr->dtim.mn = atim.mn;
889 } else {
890 ghdr->dtim.yr = ghdr->btim.yr;
891 ghdr->dtim.mo = ghdr->btim.mo;
892 ghdr->dtim.dy = ghdr->btim.dy;
893 ghdr->dtim.hr = ghdr->btim.hr;
894 ghdr->dtim.mn = ghdr->btim.mn;
895 }
896 if (len>25) {
897 ghdr->dsf = (float)gagbb(rec+26,1,15);
898 i = gagbb(rec+26,0,1);
899 if (i) ghdr->dsf = -1.0*ghdr->dsf;
900 ghdr->dsf = pow(10.0,ghdr->dsf);
901 } else ghdr->dsf = 1.0;
902
903
904 /* If it is there, get info from GDS */
905
906 if (ghdr->gdsflg) {
907 rc = fread(rec,1,3,gfile);
908 if (rc<3) return(50);
909 len = gagby(rec,0,3);
910 ghdr->gdslen = len;
911 cpos = cpos + len;
912
913 /*mf --- malloc for my generic grid where I code the lon/lats in the GDS --mf*/
914
915 gds=(char *)malloc(len+3);
916 if(gds == NULL) {
917 return(51);
918 }
919
920 rc = fread(gds+3,1,len-3,gfile);
921 if (rc<len-3) return(50);
922 ghdr->gtyp = gagby(gds,4,1);
923 ghdr->gicnt = gagby(gds,6,2);
924 ghdr->gjcnt = gagby(gds,8,2);
925 ghdr->gsf1 = gagbb(gds+27,0,1);
926 ghdr->gsf2 = gagbb(gds+27,1,1);
927 ghdr->gsf3 = gagbb(gds+27,2,1);
928
929 /*
930 printf ("gds len,typ,i,j,flags = %i %i %i %i %i %i %i\n",
931 len,ghdr->gtyp,ghdr->gicnt,ghdr->gjcnt,ghdr->gsf1,ghdr->gsf2,ghdr->gsf3);
932 */
933 free(gds);
934
935 } else ghdr->gdslen = 0;
936
937 /* Get necessary info about BMS if it is there */
938
939 if (ghdr->bmsflg) {
940 rc = fread(rec,1,6,gfile);
941 if (rc<6) return(50);
942 len = gagby(rec,0,3);
943 ghdr->bmsflg = len;
944 ghdr->bnumr = gagby(rec,4,2);
945 ghdr->bpos = cpos+6;
946 cpos = cpos + len;
947 rc = fseeko(gfile,cpos,0);
948 } else ghdr->bmslen = 0;
949
950 /* Get necessary info from data header */
951
952 rc = fread(rec,1,11,gfile);
953 if (rc<11) return(50);
954 len = gagby(rec,0,3);
955 ghdr->bdslen = len;
956 ghdr->iflg = gagbb(rec+3,0,2);
957 i = gagby(rec,4,2);
958 if (i>32767) i = 32768-i;
959 ghdr->bsf = pow(2.0,(float)i);
960
961 i = gagby(rec,6,1);
962 sign = 0;
963 if (i>127) {
964 sign = 1;
965 i = i - 128;
966 }
967 mant = gagby(rec,7,3);
968 if (sign) mant = -mant;
969 ghdr->ref = (float)mant * pow(16.0,(float)(i-70));
970
971 ghdr->bnum = gagby(rec,10,1);
972 ghdr->dpos = cpos+11;
973 /*
974 printf ("len,flg,bsf,ref,bnum %i %i %g %g %i\n",len,ghdr->iflg,
975 ghdr->bsf,ghdr->ref,ghdr->bnum);
976 */
977
978 if (ghdr->vers==0) {
979 fpos = fpos + 8 + ghdr->pdslen + ghdr->gdslen +
980 ghdr->bmslen + ghdr->bdslen;
981 } else fpos = fpos + ghdr->len;
982
983 return(0);
984
985 }
986
987 /* Routine to determine the location of the GRIB record in terms of
988 the GrADS data set, and fill in the proper values at the proper
989 slot location. */
990
gribrec(struct grhdr * ghdr,struct gafile * pfi,struct gaindx * pindx)991 int gribrec (struct grhdr *ghdr, struct gafile *pfi, struct gaindx *pindx) {
992 float (*conv) (float *, float);
993 struct gavar *pvar;
994 int i,ioff,iz,it,joff,nsiz,flag;
995 float z,t;
996
997 /* Verify that we are looking at the proper grid type */
998
999 joff =0;
1000 nsiz = nrec * nelements ;
1001 if (ghdr->iflg) {
1002 if (verb) {
1003 printf ("GRIB record contains harmonic or complex packing\n");
1004 printf (" Record is skipped.\n");
1005 printf (" Variable is %i\n",ghdr->parm);
1006 }
1007 return(10);
1008 }
1009 if (pfi->grbgrd==255 || pfi->grbgrd<-900) {
1010 if (!ghdr->gdsflg) {
1011 if (verb) {
1012 printf ("GRIB record contains pre-defined grid type: ");
1013 printf ("GrADS descriptor specifies type 255\n");
1014 gribpr(ghdr);
1015 }
1016 return(20);
1017 }
1018 if ( pfi->ppflag) {
1019 if ( ghdr->gicnt != 65535 &&
1020 ((ghdr->gicnt != pfi->ppisiz) || (ghdr->gjcnt != pfi->ppjsiz)) ) {
1021 if (verb) {
1022 printf ("GRIB grid size does not match descriptor: ");
1023 gribpr(ghdr);
1024 }
1025 return(300);
1026 }
1027 } else {
1028 if ( ghdr->gicnt != 65535 &&
1029 ((ghdr->gicnt != pfi->dnum[0]) || (ghdr->gjcnt != pfi->dnum[1])) ) {
1030 if (verb) {
1031 printf ("GRIB grid size does not match descriptor:");
1032 gribpr(ghdr);
1033 }
1034 return(301);
1035 }
1036 }
1037 } else {
1038
1039 /*---
1040 special case for
1041 GRIB grid number (dtype grib NNN) == 29
1042 ---*/
1043
1044 if (pfi->grbgrd==29) {
1045 if (ghdr->id!=29 && ghdr->id!=30) {
1046 if (verb) {
1047 printf("Record has wrong GRIB grid type: ") ;
1048 gribpr(ghdr);
1049 }
1050 return(400);
1051 }
1052 if (ghdr->id==29) joff = nelements;
1053 nsiz = 2 * nelements ;
1054 } else {
1055 if (ghdr->id != pfi->grbgrd) {
1056 if (verb) {
1057 printf("%s","Record has wrong GRIB grid type: ");
1058 gribpr(ghdr);
1059 }
1060 return(401);
1061 }
1062 }
1063 }
1064
1065 /* Calculate the grid time for this record. If it is non-integer
1066 or if it is out of bounds, just return. */
1067
1068 /* only look for tauoff hour taus */
1069
1070 if( tauflg && (ghdr->ftu==1 && ghdr->fcstt!=tauoff) ) {
1071 if (verb) {
1072 printf("%s %d","--f-- Forecast Time does not match : ",tauoff);
1073 gribpr(ghdr);
1074 }
1075 return(32);
1076 }
1077
1078 /*
1079 printf("qqq %i%i%i%i%i ... %i%i%i%i%i\n",
1080 ghdr->btim.yr,ghdr->btim.mo,ghdr->btim.dy,ghdr->btim.hr,ghdr->btim.hr,
1081 btimdd.yr,btimdd.mo,btimdd.dy,btimdd.hr,btimdd.mn);
1082 */
1083
1084 /* check if base DTG in grib field matches initial time in DD */
1085
1086 if(tau0 &&
1087 ((ghdr->btim.yr != btimdd.yr) ||
1088 (ghdr->btim.mo != btimdd.mo) ||
1089 (ghdr->btim.dy != btimdd.dy) ||
1090 (ghdr->btim.hr != btimdd.hr) ||
1091 (ghdr->btim.mn != btimdd.mn)) ) {
1092 if (verb) {
1093 printf("%s","--b-- Base Time does not match Initial Time in DD: ");
1094 gribpr(ghdr);
1095 }
1096 return(34);
1097 }
1098
1099 t = t2gr(pfi->abvals[3],&(ghdr->dtim));
1100
1101 if (t<0.99 || t>((float)(pfi->dnum[3])+0.01)) {
1102 if (verb) {
1103 printf("%s","----- Time out of bounds: ");
1104 gribpr(ghdr);
1105 }
1106 return(36);
1107 }
1108
1109 it = (int)(t+0.01);
1110
1111 if (fabs((float)it - t)>0.01) {
1112 if (verb) {
1113 printf("----- Time non-integral. %g %g: ",(float)it,t);
1114 gribpr(ghdr);
1115 }
1116 return(38);
1117 }
1118 /* it = (int)(t+0.1); */
1119 it = (it-1)*pfi->trecs;
1120
1121 /* See if we can match up this grid with a variable in
1122 the data descriptor file */
1123
1124 pvar = pfi->pvar1;
1125 i = 0;
1126 flag=0;
1127
1128
1129 while (i<pfi->vnum) {
1130
1131 if (pvar->levels>0) {
1132
1133 /*
1134 multi level data
1135 */
1136
1137 if (pvar->units[0]==ghdr->parm && pvar->units[1]==ghdr->ltyp) {
1138 /*
1139 mf 940602 - look for time range indicator match
1140 */
1141 if(diag) printf("qqq1 %d %d %d %d\n",pvar->units[3],ghdr->tri,pvar->units[0],pvar->units[1]);
1142 if(pvar->units[3] < -900 || pvar->units[3] == ghdr->tri) {
1143 conv = pfi->ab2gr[2];
1144 z = conv(pfi->abvals[2],ghdr->level);
1145 if(diag) printf("qqq2 %g %g \n",z,(float)(pvar->levels));
1146 if (z>0.99 && z<((float)(pvar->levels)+0.01)) {
1147 iz = (int)(z+0.5);
1148 /*
1149 search for level match
1150 */
1151 if(diag) printf("qqq3 %d\n",iz);
1152 if ( fabs(z-(float)iz)<0.01) {
1153 iz = (int)(z+0.5);
1154 if(diag) printf("qqq4 match\n");
1155 ioff = pvar->recoff + iz - 1;
1156 gribfill (it+ioff,joff,nsiz,ghdr,pindx);
1157 flag=1;
1158 i = pfi->vnum + 1; /* Force loop to stop */
1159 }
1160 }
1161 }
1162 }
1163 } else {
1164
1165 /*
1166 sfc data
1167 */
1168 if (pvar->units[0]==ghdr->parm && pvar->units[1]==ghdr->ltyp) {
1169 if ( (pvar->units[3]<-900 && pvar->units[2]==ghdr->level) ||
1170 (pvar->units[3]>-900 && pvar->units[2]==ghdr->l1 && pvar->units[3]==ghdr->l2 )
1171
1172 /*mf 940602
1173 test for time range indicator
1174 NB!!!! - if a variable has all the same parameters except the
1175 time range indicator then explicitly set it up!!
1176 */
1177
1178 /*mf 951129
1179
1180 took this test out.... require match in tri test for level (2 byte number) only
1181 || (pvar->units[3] == ghdr->tri && pvar->units[2]==ghdr->l1)
1182
1183 */
1184 || (pvar->units[3] == ghdr->tri && pvar->units[2]==ghdr->level) ) {
1185 ioff = pvar->recoff;
1186 gribfill (it+ioff,joff,nsiz,ghdr,pindx);
1187 i = pfi->vnum+1; /* Force loop to stop */
1188 flag=1;
1189 }
1190 }
1191 }
1192 pvar++; i++;
1193 }
1194
1195
1196 if (flag && verb) printf("!!!!! MATCH: ");
1197 if (!flag && verb) printf("..... NOOOO: ");
1198 if (verb) gribpr(ghdr);
1199
1200 return (flag ? 0 : 1);
1201
1202 }
1203
1204
1205 /* Fill in values for this record, now that we have found how it
1206 matches. We are not handling the time aspect as yet */
1207
gribfill(int ioff,int joff,int nsiz,struct grhdr * ghdr,struct gaindx * pindx)1208 void gribfill (int ioff, int joff, int nsiz, struct grhdr *ghdr,
1209 struct gaindx *pindx) {
1210 ioff = nsiz*ioff + joff;
1211 *(pindx->intpnt+ioff) = ghdr->dpos;
1212 if (ghdr->bmsflg) *(pindx->intpnt+ioff+1) = ghdr->bpos;
1213 *(pindx->intpnt+ioff+2) = ghdr->bnum;
1214 *(pindx->fltpnt+ioff) = ghdr->dsf;
1215 *(pindx->fltpnt+ioff+1) = ghdr->bsf;
1216 *(pindx->fltpnt+ioff+2) = ghdr->ref;
1217 }
1218
gribpr(struct grhdr * ghdr)1219 void gribpr(struct grhdr *ghdr) {
1220
1221 printf ("% 4i % 10i % 3i % 1i % 1i % 3i % 3i %-5i % 10i % 10i % 2i ",
1222 irec,fpos,ghdr->id,ghdr->gdsflg,ghdr->bmsflg,ghdr->parm,ghdr->ltyp,
1223 ghdr->level,ghdr->dpos,ghdr->bpos,ghdr->bnum);
1224 printf ("btim: %04i%02i%02i%02i:%02i ",ghdr->btim.yr,
1225 ghdr->btim.mo,ghdr->btim.dy,ghdr->btim.hr,ghdr->btim.mn);
1226 printf ("tau: % 6i ",ghdr->fcstt);
1227 printf ("dtim: %04i%02i%02i%02i:%02i ",ghdr->dtim.yr,
1228 ghdr->dtim.mo,ghdr->dtim.dy,ghdr->dtim.hr,ghdr->dtim.mn);
1229 printf("\n");
1230
1231 }
1232
wtgmap(void)1233 int wtgmap(void) {
1234
1235 int i,nb,bcnt,idum;
1236 float fdum;
1237 unsigned char *map;
1238 unsigned char ibmfloat[4];
1239
1240 /* calculate the size of the ver==1 index file */
1241
1242 nb=
1243 2 + (4*4) + /* version in byte 2, then 4 ints with number of each data type */
1244 pindx->hinum*sizeof(int)+
1245 pindx->hfnum*sizeof(int)+
1246 pindx->intnum*sizeof(int)+
1247 pindx->fltnum*sizeof(float) ;
1248
1249 /* add additional info */
1250
1251 nb=nb+7; /* base time (+ sec) for compatibility with earlier version 2 maps */
1252 nb=nb+8*4; /* grvals for time <-> grid conversion */
1253
1254 if(diag) {
1255 printf("\nqqq nb= %d\n",nb);
1256 printf(" btimdd.yr = %d\n",btimdd.yr);
1257 printf(" btimdd.mo = %d\n",btimdd.mo);
1258 printf(" btimdd.dy = %d\n",btimdd.dy);
1259 printf(" btimdd.hr = %d\n",btimdd.hr);
1260 printf(" btimdd.mn = %d\n",btimdd.mn);
1261 }
1262 /* allocate space for the map */
1263
1264 map = (unsigned char *)malloc(nb);
1265 if(map == NULL) {
1266 fprintf(stderr,"malloc error creating the map in gagmap\n");
1267 return(60);
1268 }
1269
1270
1271 if(diag) {
1272
1273 for(i=0;i<pindx->hinum;i++) {
1274 printf("i = %d pindx->hipnt = %d\n",i,*(pindx->hipnt+i));
1275 }
1276 for(i=0;i<pindx->hfnum;i++) {
1277 printf("i = %d pindx->hfpnt = %g\n",i,*(pindx->hfpnt+i));
1278 }
1279 for(i=0;i<pindx->intnum;i++) {
1280 printf("i = %d pindx->intpnt = %d\n",i,*(pindx->intpnt+i));
1281 }
1282 for(i=0;i<pindx->fltnum;i++) {
1283 printf("i = %d pindx->fltpnt = %g\n",i,*(pindx->fltpnt+i));
1284 }
1285
1286 }
1287 bcnt=0;
1288
1289 gapby(0,map,bcnt,1); bcnt++ ; /* set the first byte to 0 */
1290 gapby(ver,map,bcnt,1); bcnt++ ; /* set the second byte to the version number */
1291
1292 putint(pindx->hinum,map,&bcnt); /* # ints in header */
1293 putint(pindx->hfnum,map,&bcnt); /* # floats in header */
1294 putint(pindx->intnum,map,&bcnt); /* # index ints */
1295 putint(pindx->fltnum,map,&bcnt); /* # index floats */
1296
1297 /*---
1298 write out base time for consistency with earlier ver 2 maps
1299 not really needed
1300 ---*/
1301
1302 gapby(btimdd.yr,map,bcnt,2); bcnt+=2 ; /* initial year */
1303 gapby(btimdd.mo,map,bcnt,1); bcnt++ ; /* initial month */
1304 gapby(btimdd.dy,map,bcnt,1); bcnt++ ; /* initial day */
1305 gapby(btimdd.hr,map,bcnt,1); bcnt++ ; /* initial hour */
1306 gapby(btimdd.mn,map,bcnt,1); bcnt++ ; /* initial minute */
1307 gapby(0,map,bcnt,1); bcnt++ ; /* initial second */
1308
1309
1310 /*---
1311 load header
1312 ---*/
1313
1314 if(pindx->hinum) {
1315
1316 for(i=0;i<pindx->hinum;i++) {
1317 idum=*(pindx->hipnt+i);
1318 putint(idum,map,&bcnt);
1319 if(diag) printf("i = %d pindx->hipnt = %d\n",i,idum);
1320 }
1321
1322 }
1323
1324 if(pindx->hfnum) {
1325
1326 /* blank for now */
1327
1328 }
1329
1330 /*---
1331 load index
1332 ---*/
1333
1334 for(i=0;i<pindx->intnum;i++) {
1335 idum=*(pindx->intpnt+i);
1336 putint(idum,map,&bcnt);
1337 if(diag) printf("i = %d pindx->intpnt = %d\n",i,idum);
1338 }
1339
1340 for(i=0;i<pindx->fltnum;i++) {
1341 fdum=*(pindx->fltpnt+i);
1342 rc=flt2ibm(fdum, ibmfloat);
1343 if(rc<0) return(601);
1344 if(diag) printf("i = %d pindx->fltpnt = %g\n",i,fdum);
1345 memcpy(&map[bcnt],ibmfloat,4); bcnt+=4;
1346 }
1347
1348
1349 /* write out the factors for converting from grid to absolute time */
1350
1351 for(i=0;i<8;i++) {
1352 fdum=*(pfi->grvals[3]+i);
1353 rc=flt2ibm(fdum, ibmfloat);
1354 if(rc<0) return(601);
1355 if(diag) printf("i = %d grvals = %g\n",i,fdum);
1356 memcpy(&map[bcnt],ibmfloat,4); bcnt+=4;
1357 }
1358
1359
1360 if(diag) printf("qqqqq the size of the map is %d the calc size is %d\n",bcnt,nb);
1361
1362 if(write_map) {
1363 fwrite(map,1,bcnt,mfile);
1364 }
1365 free(map);
1366
1367 return(0);
1368
1369 }
1370 /*---
1371
1372 routine to dump a 4 byte int into a character stream
1373 Mike Fiorino 960829
1374
1375 ---*/
1376
putint(int dum,unsigned char * buf,int * off)1377 void putint(int dum, unsigned char *buf,int *off) {
1378
1379 int offset;
1380 offset=*off;
1381 if(dum < 0) {
1382 dum=-dum;
1383 gapby(dum,buf,offset,4);
1384 gapbb(1,buf,offset*8,1);
1385 } else {
1386 gapby(dum,buf,offset,4);
1387 }
1388 offset+=4;
1389 *off=offset;
1390
1391 }
1392