grib1_routines.c
References to this file elsewhere.
1 /*
2 ****************************************************************************
3 *
4 * Routines for indexing, reading and writing grib files. Routines
5 * are designed to be called by Fortran.
6 *
7 * All routines return 0 for success, 1 for failure, unless otherwise noted.
8 *
9 * Todd Hutchinson
10 * WSI
11 * 05/17/2005
12 *
13 ****************************************************************************
14 */
15
16 #include "grib1_routines.h"
17 #include "gridnav.h"
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <ctype.h>
21 #include <limits.h>
22
23 char *trim (char *str);
24 int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid);
25 int index_times(GribInfo *gribinfo, Times *times);
26 int find_time(Times *times, char valid_time[15]);
27 int get_gridnav_projection(int wrf_projection);
28 int get_byte(int input_int, int bytenum);
29
30 /*
31 * Allocate space for the fileindex structure
32 */
33
34 int ALLOC_INDEX_FILE(FileIndex *fileindex)
35 {
36 int status = 0;
37
38 fileindex->gribinfo = (GribInfo *)malloc(sizeof(GribInfo));
39 if (fileindex->gribinfo == NULL) {
40 fprintf(stderr,"Allocating fileindex->gribinfo failed.\n");
41 status = 1;
42 return status;
43 }
44
45 fileindex->metadata = (MetaData *)malloc(sizeof(MetaData));
46 if (fileindex->metadata == NULL) {
47 fprintf(stderr,"Allocating fileindex->metadata failed.\n");
48 status = 1;
49 return status;
50 }
51 fileindex->metadata->elements = NULL;
52
53 fileindex->times = (Times *)malloc(sizeof(Times));
54 if (fileindex->times == NULL) {
55 fprintf(stderr,"Allocating fileindex->times failed.\n");
56 status = 1;
57 return status;
58 }
59 fileindex->times->elements = NULL;
60
61 return status;
62 }
63
64
65 void FREE_INDEX_FILE(FileIndex *fileindex)
66 {
67 int status = 0;
68
69 rg_free_gribinfo_elements(fileindex->gribinfo);
70 free(fileindex->gribinfo);
71
72 free(fileindex->metadata->elements);
73 free(fileindex->metadata);
74
75 free(fileindex->times->elements);
76 free(fileindex->times);
77
78 }
79
80
81 int INDEX_FILE(int *fid, FileIndex *fileindex)
82 {
83
84 int status;
85 /* Index the grib records */
86
87 status = rg_setup_gribinfo_i(fileindex->gribinfo,*fid,1);
88 if (status < 0) {
89 fprintf(stderr,"Error setting up gribinfo structure.\n");
90 return 1;
91 }
92
93 /* Index the metadata section */
94
95 status = index_metadata(fileindex->gribinfo, fileindex->metadata, *fid);
96 if (status != 0) {
97 fprintf(stderr,"Error setting up metadata structure.\n");
98 return 1;
99 }
100
101 /* Setup a list of times based on times in grib records */
102
103 status = index_times(fileindex->gribinfo, fileindex->times);
104 if (status != 0) {
105 fprintf(stderr,"Error indexing times in grib file.\n");
106 return 1;
107 }
108
109 return 0;
110 }
111
112
113 int GET_FILEINDEX_SIZE(int *size)
114 {
115 *size = sizeof(FileIndex);
116 return *size;
117 }
118
119
120 int GET_NUM_TIMES(FileIndex *fileindex, int *numtimes)
121 {
122 *numtimes = (fileindex->times)->num_elements;
123 return *numtimes;
124 }
125
126
127 int GET_TIME(FileIndex *fileindex, int *idx, char time[])
128 {
129 int num_times;
130 int year, month, day, minute, hour, second;
131 char time2[100];
132
133 num_times = GET_NUM_TIMES(fileindex,&num_times);
134 if (*idx > num_times)
135 {
136 fprintf(stderr,"Tried to get time %d, but only %d times exist\n",
137 *idx, num_times);
138 return 1;
139 }
140
141 strcpy(time,fileindex->times->elements[*idx-1].valid_time);
142
143 /* Reformat time to meet WRF time format */
144
145 sscanf(time, "%4d%2d%2d%2d%2d%2d",
146 &year, &month, &day, &hour, &minute, &second);
147 sprintf(time2, "%04d-%02d-%02d_%02d:%02d:%02d",
148 year, month, day, hour, minute, second);
149 strncpy(time,time2,19);
150
151 return 0;
152 }
153
154
155 int GET_LEVEL1(FileIndex *fileindex, int *idx, int *level1)
156 {
157
158 *level1 = (fileindex->gribinfo)->elements[*idx].usHeight1;
159
160 return *level1;
161
162 }
163
164
165 int GET_LEVEL2(FileIndex *fileindex, int *idx, int *level2)
166 {
167
168 *level2 = (fileindex->gribinfo)->elements[*idx].usHeight2;
169
170 return *level2;
171
172 }
173
174
175 int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid)
176 {
177 int status=0;
178 int end;
179 char string[11];
180 int found_metadata=0;
181 int idx;
182 int pos;
183 int fileend;
184 int seekpos;
185 int bytesread;
186 char line[1000];
187 char element[100],datestr[100],varname[100];
188 char value[1000];
189 int incomment;
190 int charidx;
191 int elemidx=0;
192 FILE *stream;
193
194
195 /* Associate a FILE *stream with the file id */
196 stream = fdopen(fid,"r");
197 if (stream == NULL)
198 {
199 perror("Error associating stream with file descriptor");
200 status = -1;
201 return status;
202 }
203
204 /*
205 * First, set the position to end of grib data (the
206 * metadata section comes after the grib data).
207 */
208 idx = rg_num_elements(gribinfo) - 1;
209 end = rg_get_end(gribinfo,idx);
210 pos = end + 1;
211 fileend = lseek(fid,0,SEEK_END);
212
213 /*
214 * Now, start searching for metadata
215 */
216 while (pos < (fileend - 10))
217 {
218 seekpos = pos;
219 pos = lseek(fid,seekpos,SEEK_SET);
220 if (pos != seekpos)
221 {
222 fprintf(stderr,"Error seeking %d bytes in file\n",end);
223 perror("");
224 return 1;
225 }
226
227 bytesread = read(fid,string,10);
228 if (bytesread != 10)
229 {
230 fprintf(stderr,"Invalid read, pos: %d :\n",pos);
231 perror("");
232 pos += 1;
233 continue;
234 }
235
236 if (strncmp(string,"<METADATA>",10) == 0)
237 {
238 /* We found it, so break out ! */
239 found_metadata = 1;
240 break;
241 }
242 pos += 1;
243 }
244
245
246 /* Now, read metadata, line by line */
247 incomment = 0;
248 while(fgets(line,1000,stream) != NULL)
249 {
250 trim(line);
251
252 /* Set comment flag, if we found a comment */
253 if (strncmp(line,"<!--",4) == 0)
254 {
255 incomment = 1;
256 strcpy(line,line+4);
257 }
258
259 /* Search for end of comment */
260 if (incomment)
261 {
262 charidx = 0;
263 while (charidx < strlen(line))
264 {
265
266 if (strncmp(line+charidx,"-->",3) == 0)
267 {
268 strcpy(line,line+charidx+3);
269 incomment = 0;
270 break;
271 }
272 else
273 {
274 charidx++;
275 }
276 }
277 }
278
279 if (incomment) continue;
280
281
282 /* Check for end of metadata */
283 if (strncmp(line,"</METADATA>",11) == 0)
284 {
285 /* We found end of data, so, break out */
286 break;
287 }
288
289 /* Skip blank lines */
290 if (strlen(line) == 0) continue;
291
292
293 /* Parse line */
294 trim(line);
295 strcpy(element,"none");
296 strcpy(datestr,"none");
297 strcpy(varname,"none");
298 strcpy(value,"none");
299 if (sscanf(line,"%[^;=];%[^;=];%[^;=]=%[^\n]",varname,datestr,
300 element,value) == 4)
301 {
302 }
303 else if (sscanf(line,"%[^;=];%[^;=]=%[^\n]",datestr,element,value) == 3)
304 {
305 strcpy(varname,"none");
306 }
307 else if (sscanf(line,"%[^;=]=%[^\n]",element,value) == 2)
308 {
309 strcpy(varname,"none");
310 strcpy(datestr,"none");
311 }
312 else
313 {
314 strcpy(varname,"none");
315 strcpy(datestr,"none");
316 strcpy(element,"none");
317 strcpy(value,"none");
318 fprintf(stderr,"Invalid line in metadata: \n%s",line);
319 }
320
321 trim(varname);
322 trim(datestr);
323 trim(element);
324 trim(value);
325
326 metadata->elements =
327 (MetaData_Elements *)realloc( metadata->elements,
328 (elemidx+1)*sizeof(MetaData_Elements) );
329 strcpy(metadata->elements[elemidx].VarName,varname);
330 strcpy(metadata->elements[elemidx].DateStr,datestr);
331 strcpy(metadata->elements[elemidx].Element,element);
332 strcpy(metadata->elements[elemidx].Value,value);
333
334 elemidx++;
335
336 }
337
338 metadata->num_elements = elemidx;
339
340 return 0;
341 }
342
343
344
345
346 int index_times(GribInfo *gribinfo, Times *times)
347 {
348 int idx;
349 int status;
350 int numtimes=0;
351 int date;
352 char valid_time[15];
353 char tmp[15];
354 int swapped;
355
356 times->num_elements = 0;
357
358 /* Loop through elements, and build list of times */
359
360 for (idx=0; idx < gribinfo->num_elements; idx++)
361 {
362 /* Calculate valid time */
363 status = rg_get_valid_time(gribinfo,idx,valid_time);
364 if (status != 0)
365 {
366 fprintf(stderr,"Could not retrieve valid time for index: %d\n",idx);
367 continue;
368 }
369
370 /*
371 * Check if this time is already contained in times
372 * If not, allocate space for it, and add it to list
373 */
374 if (find_time(times,valid_time) < 0)
375 {
376 times->num_elements++;
377 times->elements =
378 (Times_Elements *)
379 realloc(times->elements,times->num_elements*sizeof(Times_Elements));
380 if (times->elements == NULL)
381 {
382 fprintf(stderr,"Allocating times->elements failed.\n");
383 status = 1;
384 return status;
385 }
386 strcpy(times->elements[times->num_elements - 1].valid_time,valid_time);
387 }
388 }
389
390 /* Sort times */
391 swapped = 1;
392 while (swapped)
393 {
394 swapped=0;
395 for (idx=1; idx < times->num_elements; idx++)
396 {
397 if (strcmp(times->elements[idx-1].valid_time,
398 times->elements[idx].valid_time) > 0)
399 {
400 strcpy(tmp,times->elements[idx-1].valid_time);
401 strcpy(times->elements[idx-1].valid_time,
402 times->elements[idx].valid_time);
403 strcpy(times->elements[idx].valid_time, tmp);
404 swapped = 1;
405 }
406 }
407 }
408
409 return 0;
410 }
411
412
413
414 int find_time(Times *times, char valid_time[15])
415 {
416 int idx;
417 int found_elem = -1;
418
419 for (idx = 0; idx < times->num_elements; idx++)
420 {
421 if (strcmp(times->elements[idx].valid_time,valid_time) == 0)
422 {
423 found_elem = idx;
424 break;
425 }
426 }
427
428 return found_elem;
429
430 }
431
432
433 int GET_METADATA_VALUE(FileIndex *fileindex, char ElementIn[],
434 char DateStrIn[], char VarNameIn[], char Value[],
435 int *stat, int strlen1, int strlen2, int strlen3,
436 int strlen4, int strlen5)
437 {
438 int elemidx;
439 int elemnum;
440 char VarName[200];
441 char DateStr[200];
442 char Element[200];
443 int Value_Len;
444
445 *stat = 0;
446
447 strncpy(Element,ElementIn,strlen2);
448 Element[strlen2] = '\0';
449 strncpy(DateStr,DateStrIn,strlen3);
450 DateStr[strlen3] = '\0';
451 strncpy(VarName,VarNameIn,strlen4);
452 VarName[strlen4] = '\0';
453 Value_Len = strlen5;
454
455 elemnum = -1;
456 for (elemidx = 0; elemidx < fileindex->metadata->num_elements; elemidx++)
457 {
458 if (strcmp(Element,fileindex->metadata->elements[elemidx].Element) == 0)
459 {
460 if (strcmp(DateStr,fileindex->metadata->elements[elemidx].DateStr)
461 == 0)
462 {
463 if (strcmp(VarName,
464 fileindex->metadata->elements[elemidx].VarName) == 0)
465 {
466 elemnum = elemidx;
467 break;
468 }
469 }
470 }
471 }
472
473 if (elemnum != -1)
474 {
475 strncpy(Value, fileindex->metadata->elements[elemnum].Value, Value_Len);
476 }
477 else
478 {
479 strncpy(Value, "none", Value_Len);
480 *stat = 1;
481 }
482
483 /*
484 * Pad end of string with one space. This allows Fortran internal
485 * read function to work properly.
486 */
487 if (strlen(Value) < Value_Len)
488 {
489 strcpy(Value + strlen(Value), " ");
490 }
491
492 return elemidx;
493 }
494
495
496 int GET_GRIB_INDEX(FileIndex *fileindex, int *center, int *subcenter,
497 int *parmtbl, int *parmid, char DateStrIn[],
498 int *leveltype, int *level1, int *level2, int *fcsttime1,
499 int *fcsttime2, int *index, int strlen1, int strlen2)
500 {
501 char DateStr[1000];
502 FindGrib findgrib;
503
504 strncpy(DateStr,DateStrIn,strlen2);
505 DateStr[strlen2] = '\0';
506 grib_time_format(DateStr,DateStr);
507
508 rg_init_findgrib(&findgrib);
509
510 strncpy(findgrib.initdate,DateStrIn,strlen2);
511 findgrib.initdate[strlen2] = '\0';
512 findgrib.parmid = *parmid;
513 findgrib.leveltype = *leveltype;
514 findgrib.level1 = *level1;
515 findgrib.level2 = *level2;
516 findgrib.fcsttime1 = *fcsttime1;
517 findgrib.fcsttime2 = *fcsttime2;
518 findgrib.center_id = *center;
519 findgrib.subcenter_id = *subcenter;
520 findgrib.parmtbl_version = *parmtbl;
521
522 *index = rg_get_index(fileindex->gribinfo, &findgrib);
523
524 return *index;
525
526 }
527
528
529 int GET_GRIB_INDEX_GUESS(FileIndex *fileindex, int *center, int *subcenter,
530 int *parmtbl, int *parmid, char DateStrIn[],
531 int *leveltype, int *level1, int *level2,
532 int *fcsttime1,int *fcsttime2, int *guessidx,
533 int *index, int strlen1, int strlen2)
534 {
535 char DateStr[1000];
536 FindGrib findgrib;
537
538 strncpy(DateStr,DateStrIn,strlen2);
539 DateStr[strlen2] = '\0';
540 grib_time_format(DateStr,DateStr);
541
542 rg_init_findgrib(&findgrib);
543
544 strncpy(findgrib.initdate,DateStrIn,strlen2);
545 findgrib.initdate[strlen2] = '\0';
546 findgrib.parmid = *parmid;
547 findgrib.leveltype = *leveltype;
548 findgrib.level1 = *level1;
549 findgrib.level2 = *level2;
550 findgrib.fcsttime1 = *fcsttime1;
551 findgrib.fcsttime2 = *fcsttime2;
552 findgrib.center_id = *center;
553 findgrib.subcenter_id = *subcenter;
554 findgrib.parmtbl_version = *parmtbl;
555
556 *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx);
557
558 return *index;
559
560 }
561
562
563 int GET_GRIB_CENTER(FileIndex *fileindex, int *parmid, int *center)
564 {
565
566 *center = rg_get_center_id(fileindex->gribinfo,*parmid);
567
568 return *center;
569
570 }
571
572
573 int GET_GRIB_SUBCENTER(FileIndex *fileindex, int *parmid, int *subcenter)
574 {
575
576 *subcenter = rg_get_subcenter_id(fileindex->gribinfo,*parmid);
577
578 return *subcenter;
579
580 }
581
582
583 int GET_GRIB_TBLVERSION(FileIndex *fileindex, int *parmid, int *parmtbl)
584 {
585
586 *parmtbl = rg_get_parmtbl(fileindex->gribinfo,*parmid);
587
588 return *parmtbl;
589
590 }
591
592
593 int GET_GRIB_PROCID(FileIndex *fileindex, int *parmid, int *proc_id)
594 {
595
596 *proc_id = rg_get_proc_id(fileindex->gribinfo,*parmid);
597
598 return *proc_id;
599
600 }
601
602
603 int GET_GRIB_INDEX_VALIDTIME(FileIndex *fileindex, int *center,
604 int *subcenter, int *parmtbl, int *parmid,
605 char DateStrIn[], int *leveltype, int *level1, int *level2,
606 int *index, int strlen1, int strlen2)
607 {
608 char DateStr[1000];
609 FindGrib findgrib;
610
611 strncpy(DateStr,DateStrIn,strlen2);
612 DateStr[strlen2] = '\0';
613 grib_time_format(DateStr,DateStr);
614
615 rg_init_findgrib(&findgrib);
616
617 strncpy(findgrib.validdate,DateStr,strlen2);
618 findgrib.initdate[strlen2] = '\0';
619 findgrib.parmid = *parmid;
620 findgrib.leveltype = *leveltype;
621 findgrib.level1 = *level1;
622 findgrib.level2 = *level2;
623 findgrib.center_id = *center;
624 findgrib.subcenter_id = *subcenter;
625 findgrib.parmtbl_version = *parmtbl;
626
627 *index = rg_get_index(fileindex->gribinfo, &findgrib);
628
629 return *index;
630 }
631
632
633 int GET_GRIB_INDEX_VALIDTIME_GUESS(FileIndex *fileindex, int *center,
634 int *subcenter, int *parmtbl, int *parmid,
635 char DateStrIn[], int *leveltype,
636 int *level1, int *level2, int *guessidx,
637 int *index, int strlen1, int strlen2)
638 {
639 char DateStr[1000];
640 FindGrib findgrib;
641
642 strncpy(DateStr,DateStrIn,strlen2);
643 DateStr[strlen2] = '\0';
644 grib_time_format(DateStr,DateStr);
645
646 rg_init_findgrib(&findgrib);
647
648 strncpy(findgrib.validdate,DateStr,strlen2);
649 findgrib.initdate[strlen2] = '\0';
650 findgrib.parmid = *parmid;
651 findgrib.leveltype = *leveltype;
652 findgrib.level1 = *level1;
653 findgrib.level2 = *level2;
654 findgrib.center_id = *center;
655 findgrib.subcenter_id = *subcenter;
656 findgrib.parmtbl_version = *parmtbl;
657
658 *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx);
659
660 return *index;
661 }
662
663
664 int GET_GRIB_INDICES(FileIndex *fileindex, int *center, int *subcenter,
665 int *parmtbl,int *parmid, char DateStrIn[],
666 int *leveltype, int *level1, int *level2, int *fcsttime1,
667 int *fcsttime2, int *indices, int *num_indices,
668 int strlen1, int strlen2)
669 {
670 char DateStr[1000];
671 int status;
672 FindGrib findgrib;
673
674 strncpy(DateStr,DateStrIn,strlen2);
675 DateStr[strlen2] = '\0';
676 grib_time_format(DateStr,DateStr);
677
678 rg_init_findgrib(&findgrib);
679
680 strncpy(findgrib.initdate,DateStrIn,strlen2);
681 findgrib.initdate[strlen2] = '\0';
682 trim(findgrib.initdate);
683 findgrib.parmid = *parmid;
684 findgrib.leveltype = *leveltype;
685 findgrib.level1 = *level1;
686 findgrib.level2 = *level2;
687 findgrib.fcsttime1 = *fcsttime1;
688 findgrib.fcsttime2 = *fcsttime2;
689 findgrib.center_id = *center;
690 findgrib.subcenter_id = *subcenter;
691 findgrib.parmtbl_version = *parmtbl;
692
693 *num_indices = rg_get_indices(fileindex->gribinfo, &findgrib, indices);
694
695 return (*num_indices);
696
697 }
698
699
700 int GET_GRID_INFO_SIZE(int *size)
701 {
702
703 *size = sizeof(Grid_Info);
704
705 return *size;
706
707 }
708
709
710 int LOAD_GRID_INFO(char *varnameIn, char *initdateIn, int *leveltype,
711 int *level1, int *level2, float *fcst_time,
712 int *accum_period, int *grid_id, int *projection,
713 int *xpoints, int *ypoints, float *center_lat,
714 float *center_lon, float *Di, float *Dj,float *central_lon,
715 int *proj_center_flag, float *latin1,
716 float *latin2, Grib1_Tables *grib_tables,
717 Grid_Info *grid_info, int strlen1, int strlen2)
718 {
719
720 char varname[1000], initdate[1000];
721
722 strncpy(varname,varnameIn,strlen1);
723 varname[strlen1] = '\0';
724 strncpy(initdate,initdateIn,strlen2);
725 initdate[strlen2] = '\0';
726
727 strcpy(grid_info->varname, varname);
728 strcpy(grid_info->initdate, initdate);
729 grid_info->leveltype = *leveltype;
730 grid_info->level1 = *level1 ;
731 grid_info->level2 = *level2 ;
732 grid_info->fcst_time = *fcst_time ;
733 grid_info->accum_period = *accum_period ;
734 grid_info->grid_id = *grid_id ;
735 grid_info->projection = *projection ;
736 grid_info->xpoints = *xpoints ;
737 grid_info->ypoints = *ypoints ;
738 grid_info->center_lat = *center_lat ;
739 grid_info->center_lon = *center_lon;
740 grid_info->Di = *Di ;
741 grid_info->Dj = *Dj ;
742 grid_info->central_lon = *central_lon ;
743 grid_info->proj_center_flag = *proj_center_flag ;
744 grid_info->latin1 = *latin1 ;
745 grid_info->latin2 = *latin2 ;
746 grid_info->grib_tables = copy_grib_tables(grib_tables);
747
748 return 0;
749
750 }
751
752 int PRINT_GRID_INFO(Grid_Info *grid_info)
753 {
754
755 fprintf(stdout,"varname =%s\n",grid_info->varname);
756 fprintf(stdout,"initdate =%s\n",grid_info->initdate);
757 fprintf(stdout,"leveltype =%d\n",grid_info->leveltype);
758 fprintf(stdout,"level1 =%d\n",grid_info->level1);
759 fprintf(stdout,"level2 =%d\n",grid_info->level2);
760 fprintf(stdout,"fcst_time =%f\n",grid_info->fcst_time);
761 fprintf(stdout,"accum_period =%d\n",grid_info->accum_period);
762 fprintf(stdout,"grid_id =%d\n",grid_info->grid_id);
763 fprintf(stdout,"projection =%d\n",grid_info->projection);
764 fprintf(stdout,"xpoints =%d\n",grid_info->xpoints);
765 fprintf(stdout,"ypoints =%d\n",grid_info->ypoints);
766 fprintf(stdout,"center_lat =%f\n",grid_info->center_lat);
767 fprintf(stdout,"center_lon =%f\n",grid_info->center_lon);
768 fprintf(stdout,"Di =%f\n",grid_info->Di);
769 fprintf(stdout,"Dj =%f\n",grid_info->Dj);
770 fprintf(stdout,"central_lon =%f\n",grid_info->central_lon);
771 fprintf(stdout,"proj_center_flag =%d\n",grid_info->proj_center_flag);
772 fprintf(stdout,"latin1 =%f\n",grid_info->latin1);
773 fprintf(stdout,"latin2 =%f\n",grid_info->latin2);
774
775 return 0;
776
777 }
778
779
780 int GET_SIZEOF_GRID(FileIndex *fileindex, int *index, int *numcols,
781 int *numrows)
782 {
783
784 *numcols = rg_get_numcols(fileindex->gribinfo,*index);
785
786 *numrows = rg_get_numrows(fileindex->gribinfo,*index);
787
788 return (*numcols)*(*numrows);
789
790 }
791
792
793 void FREE_GRID_INFO(Grid_Info *grid_info)
794 {
795 FREE_GRIBMAP(grid_info->grib_tables);
796 }
797
798
799 int READ_GRIB(FileIndex *fileindex, int *fid, int *index, float *data)
800 {
801 int status;
802
803 status = rg_get_data_1d(fileindex->gribinfo,*index,data);
804
805 return status;
806 }
807
808 #define WRF_LATLON 0
809 #define WRF_LAMBERT 1
810 #define WRF_POLAR_STEREO 2
811 #define WRF_MERCATOR 3
812
813 #define LINESIZE 300
814
815 #define SECS_IN_SEC 1
816 #define SECS_IN_MIN 60
817 #define MINS_IN_HOUR 60
818 #define MINS_IN_5MINS 5
819 #define HOURS_IN_DAY 24
820
821 #define MAX_FCST 65535
822 #define MAX_FCST_SECS MAX_FCST*SECS_IN_SEC
823 #define MAX_FCST_MINS MAX_FCST*SECS_IN_MIN
824 #define MAX_FCST_5MINS MAX_FCST*MINS_IN_5MINS*SECS_IN_MIN
825 #define MAX_FCST_HOURS MAX_FCST*MINS_IN_HOUR*SECS_IN_MIN
826 #define MAX_FCST_DAYS MAX_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN
827
828 #define MAX1B_FCST 256
829 #define MAX1B_FCST_SECS MAX1B_FCST*SECS_IN_SEC
830 #define MAX1B_FCST_MINS MAX1B_FCST*SECS_IN_MIN
831 #define MAX1B_FCST_5MINS MAX1B_FCST*MINS_IN_5MINS*SECS_IN_MIN
832 #define MAX1B_FCST_HOURS MAX1B_FCST*MINS_IN_HOUR*SECS_IN_MIN
833 #define MAX1B_FCST_DAYS MAX1B_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN
834
835 #define PI 3.1415
836 typedef struct {
837 int time_range;
838 int fcst_unit;
839 int P1;
840 int P2;
841
842 int time_range_ext;
843 int fcst_unit_ext_1;
844 int fcst_unit_ext_2;
845 int P1_ext;
846 int P2_ext;
847 } FcstTimeStruct;
848
849 int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *fcst_time);
850
851 /****************************************************************************
852 *
853 * This function takes in metadata in the grid_info structure, output data in
854 * the *data array, and calls routines to write the metadata and data
855 * in grib version 1 format the open file descriptor filefd.
856 *
857 ****************************************************************************/
858
859 int WRITE_GRIB(Grid_Info *grid_info, int *filefd, float *data)
860 {
861
862 GRIB_HDR *gh=NULL;
863 DATA_INPUT data_input;
864 GEOM_IN geom_in;
865 USER_INPUT user_input;
866 int grid_projection;
867 int status;
868 float x_center, y_center;
869 GridNav gridnav;
870 float first_lat, first_lon, last_lat, last_lon;
871 int year, month, day, hour, minute;
872 float second;
873 char varname2[1000];
874 int table_index;
875 int fcst_unit;
876 int time_range;
877 int P1, P2;
878 int fcst_unit_ext_1, fcst_unit_ext_2;
879 int P1_ext, P2_ext;
880 int time_range_ext;
881 char errmsg[1000];
882 int center, subcenter, parmtbl;
883 int tablenum;
884 FcstTimeStruct fcst_time;
885
886 strcpy(varname2,grid_info->varname);
887 trim(varname2);
888
889 sscanf(grid_info->initdate,"%d-%d-%d_%d:%d:%f",
890 &year,&month,&day,&hour,&minute,&second);
891
892 /* Get coords of center of grid */
893 x_center = (grid_info->xpoints + 1)/2.;
894 y_center = (grid_info->ypoints + 1)/2.;
895
896 grid_projection = get_gridnav_projection(grid_info->projection);
897
898 /* Initialize grid structure */
899 status = GRID_init(grid_info->center_lat, grid_info->central_lon,
900 grid_projection,
901 grid_info->latin1, grid_info->latin2,
902 grid_info->xpoints, grid_info->ypoints, grid_info->Di,
903 grid_info->Dj,
904 grid_info->center_lat, grid_info->center_lon,
905 x_center, y_center,
906 &gridnav);
907 if (!status)
908 {
909 fprintf(stderr,"write_grib: error from GRID_init\n");
910 }
911
912 /* get lat/lon of lower left corner */
913 status = GRID_to_latlon(&gridnav, 1, 1, &first_lat, &first_lon);
914 if (!status)
915 {
916 fprintf(stderr,
917 "write_grib: error from GRID_to_latlon for first lat/lon\n");
918 }
919
920 /* get lat/lon of upper right corner */
921 status = GRID_to_latlon(&gridnav, grid_info->xpoints, grid_info->ypoints,
922 &last_lat, &last_lon);
923 if (!status)
924 {
925 fprintf(stderr,
926 "write_grib: error from GRID_to_latlon for last lat/lon\n");
927 }
928
929 /* Read the grib parameter table */
930 status = GET_GRIB_PARAM(grid_info->grib_tables, varname2, ¢er,
931 &subcenter, &parmtbl, &tablenum, &table_index,
932 1,strlen(varname2));
933 if (table_index < 0)
934 {
935 fprintf(stderr,\
936 "Skipping %s, Could not find parameter for %s in gribmap.txt\n",\
937 varname2,varname2);
938 return 1;
939 }
940
941 /*
942 * We skip any parameters that are listed in parameter 255 in gribmap.txt.
943 * Parameter 255 is used to indicate that a WRF parameter should not be
944 * output. It is useful for parameters that are requested to be output in
945 * the WRF Registry, but are already implicitly output in grib.
946 */
947
948 if (table_index == 255)
949 {
950 return 0;
951 }
952
953 /*
954 * Setup the geom_in structure for the grib library. Here, we set
955 * the generic parms. Below, we set the projection specific parms
956 */
957 geom_in.nx = grid_info->xpoints;
958 geom_in.ny = grid_info->ypoints;
959 geom_in.first_lat = first_lat;
960 geom_in.first_lon = first_lon;
961 geom_in.last_lat = last_lat;
962 geom_in.last_lon = last_lon;
963 geom_in.scan = 64;
964
965 switch (grid_info->projection)
966 {
967 case WRF_LATLON:
968 strcpy(geom_in.prjn_name,"spherical");
969 geom_in.parm_1 = grid_info->Dj;
970 geom_in.parm_2 = grid_info->Di;
971 geom_in.parm_3 = -1;
972 geom_in.usRes_flag = 0; /*
973 * Set to 0 here, MEL grib library will reset
974 * to 128 to indicate that direction
975 * increments are given.
976 */
977 break;
978 case WRF_MERCATOR:
979 strcpy(geom_in.prjn_name,"mercator");
980 geom_in.parm_1 = grid_info->latin1;
981 geom_in.parm_2 = grid_info->Di;
982 geom_in.parm_3 = grid_info->Dj;
983 geom_in.usRes_flag = 128;
984 break;
985 case WRF_LAMBERT:
986 strcpy(geom_in.prjn_name,"lambert");
987 geom_in.usRes_flag = 0; /* Set to 0 here, MEL grib library will reset
988 * to 128.
989 */
990 geom_in.parm_3 = grid_info->central_lon;
991 geom_in.x_int_dis = grid_info->Di;
992 geom_in.y_int_dis = grid_info->Dj;
993 geom_in.parm_1 = grid_info->latin1;
994 geom_in.parm_2 = grid_info->latin2;
995 break;
996 case WRF_POLAR_STEREO:
997 strcpy(geom_in.prjn_name,"polar_stereo");
998 geom_in.usRes_flag = 0; /* Set to 0 here, MEL grib library will reset
999 * to 128.
1000 */
1001
1002 geom_in.parm_3 = -1;
1003 geom_in.x_int_dis = grid_info->Di*(1.+sin(60. * PI/180.))
1004 / (1.+sin(abs(grid_info->latin1) * PI/180.));
1005 geom_in.y_int_dis = grid_info->Dj*(1.+sin(60. * PI/180.))
1006 / (1.+sin(abs(grid_info->latin1) * PI/180.));
1007 geom_in.parm_1 = -1;
1008 geom_in.parm_2 = grid_info->central_lon;
1009 break;
1010 default:
1011 fprintf(stderr,"Error, invalid projection: %d\n",grid_info->projection);
1012 return 1;
1013 }
1014
1015 /*
1016 * Setup the data_input structure.
1017 */
1018 data_input.nDec_sc_fctr =
1019 grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index];
1020 data_input.usProc_id = 220;
1021 data_input.usGrid_id = grid_info->grid_id;
1022 data_input.usParm_id =
1023 grid_info->grib_tables->grib_table_info[tablenum].parm_id[table_index];
1024 data_input.usParm_sub_id =
1025 grid_info->grib_tables->grib_table_info[tablenum].subcenter;
1026 data_input.usLevel_id = grid_info->leveltype;
1027
1028 if (grid_info->leveltype == 112) {
1029 data_input.nLvl_1 = grid_info->level2;
1030 data_input.nLvl_2 = grid_info->level1;
1031 } else {
1032 data_input.nLvl_1 = grid_info->level1;
1033 data_input.nLvl_2 = grid_info->level2;
1034 }
1035
1036 data_input.nYear = year;
1037 data_input.nMonth = month;
1038 data_input.nDay = day;
1039 data_input.nHour = hour;
1040 data_input.nMinute = minute;
1041 data_input.nSecond = second;
1042
1043 status = get_fcst_time(grid_info->accum_period, grid_info->fcst_time,
1044 &fcst_time);
1045
1046 data_input.usFcst_id = fcst_time.fcst_unit;
1047 data_input.usFcst_per1 = fcst_time.P1;
1048 data_input.usFcst_per2 = fcst_time.P2;
1049 data_input.usTime_range_id = fcst_time.time_range;
1050 data_input.usTime_range_avg = 0;
1051 data_input.usTime_range_mis = 0;
1052 /*
1053 * This is for WSI's extended PDS section
1054 */
1055 data_input.PDS_41 = fcst_time.fcst_unit_ext_1;
1056 data_input.PDS_42 = fcst_time.P1_ext;
1057 data_input.PDS_46 = fcst_time.fcst_unit_ext_2;
1058 data_input.PDS_47 = fcst_time.P2_ext;
1059 data_input.PDS_51 = fcst_time.time_range_ext;
1060 data_input.PDS_52 = 0;
1061
1062
1063 data_input.nDec_sc_fctr =
1064 grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index];
1065 user_input.usCenter_id =
1066 grid_info->grib_tables->grib_table_info[tablenum].center;
1067 user_input.usParm_tbl =
1068 grid_info->grib_tables->grib_table_info[tablenum].parmtbl;
1069 user_input.chCase_id='0';
1070 user_input.usSub_tbl = 0;
1071 user_input.usCenter_sub =
1072 grid_info->grib_tables->grib_table_info[tablenum].subcenter;
1073 user_input.usTrack_num = 0;
1074 user_input.usGds_bms_id = 128;
1075 user_input.usBDS_flag = 0;
1076 user_input.usBit_pack_num = 0;
1077
1078 status = init_gribhdr(&gh,errmsg);
1079 if (status != 0) {
1080 fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1081 return 1;
1082 }
1083
1084 status = grib_enc(data_input,user_input,geom_in,data,gh,errmsg);
1085 if (status != 0) {
1086 fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1087 fprintf (stderr,"\tCheck precision for %s in gribmap.txt.\n",
1088 varname2);
1089 return 1;
1090 }
1091
1092 status = gribhdr2filed(gh,*filefd,errmsg);
1093 if (status != 0) {
1094 fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1095 return 1;
1096 }
1097
1098 free_gribhdr(&gh);
1099
1100 return 0;
1101 }
1102
1103 /***************************************************************************
1104 * Function to set up a structure containing forecast time parameters
1105 * This encodes the standard grib forecast time parameters as well
1106 * as WSI's extended forecast time parameters.
1107 ***************************************************************************/
1108
1109 int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *ft)
1110 {
1111 /*
1112 * Added ability to output a "5-minute" forecast time unit for the
1113 * sake of WxProducer. This allows WxPro to ingest data beyond
1114 * 18 hours, and accumulation data beyond 255 units.
1115 */
1116
1117 /*
1118 * Initialize.
1119 */
1120 ft->time_range = 0;
1121 ft->fcst_unit = 0;
1122 ft->P1 = 0;
1123 ft->P2 = 0;
1124 ft->fcst_unit_ext_1 = 0;
1125 ft->fcst_unit_ext_2 = 0;
1126 ft->P1_ext = 0;
1127 ft->P2_ext = 0;
1128 ft->time_range_ext = 0;
1129
1130 if (accum_period == 0)
1131 {
1132 if (fcst_secs < MAX_FCST_SECS)
1133 {
1134 ft->time_range = 10;
1135 ft->fcst_unit = 254;
1136 ft->P1 = get_byte(fcst_secs,2);
1137 ft->P2 = get_byte(fcst_secs,1);
1138 }
1139 else if (((fcst_secs % SECS_IN_MIN) == 0) &&
1140 (fcst_secs < MAX_FCST_MINS))
1141 {
1142 ft->time_range = 10;
1143 ft->fcst_unit = 0;
1144 ft->P1 = get_byte(fcst_secs/SECS_IN_MIN,2);
1145 ft->P2 = get_byte(fcst_secs/SECS_IN_MIN,1);
1146 }
1147 else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR) == 0) &&
1148 (fcst_secs < MAX_FCST_HOURS))
1149 {
1150 ft->time_range = 10;
1151 ft->fcst_unit = 1;
1152 ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),2);
1153 ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),1);
1154 }
1155 /*
1156 * MAX_FCST_DAYS is causing an integer overflow, so, we'll just skip
1157 * the check here. It's very unlikely that someone would exceed this
1158 * anyway (5.6 million days!)
1159 */
1160 /*
1161 else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0) &&
1162 (fcst_secs < MAX_FCST_DAYS))
1163 */
1164 else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0))
1165 {
1166 ft->time_range = 10;
1167 ft->fcst_unit = 2;
1168 ft->P1 =
1169 get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),2);
1170 ft->P2 =
1171 get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),1);
1172 }
1173 else if (((fcst_secs % SECS_IN_MIN*MINS_IN_5MINS) == 0)
1174 && (fcst_secs < MAX_FCST_5MINS))
1175 {
1176 ft->time_range = 10;
1177 ft->fcst_unit = 50;
1178 ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),2);
1179 ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),1);
1180 }
1181 else
1182 {
1183 ft->time_range = 255;
1184 ft->fcst_unit = 0;
1185 ft->P1 = 0;
1186 ft->P2 = 0;
1187
1188 ft->fcst_unit_ext_1 = 254;
1189 ft->fcst_unit_ext_2 = 254;
1190 ft->P1_ext = fcst_secs;
1191 ft->P2_ext = 0;
1192 ft->time_range_ext = 0;
1193 }
1194 }
1195 else /* Accumulation period is not 0 */
1196 {
1197 if ((fcst_secs < MAX1B_FCST_HOURS) &&
1198 (fcst_secs%(SECS_IN_MIN*MINS_IN_HOUR) == 0) &&
1199 (accum_period%(SECS_IN_MIN*MINS_IN_HOUR) == 0))
1200 {
1201 ft->time_range = 4;
1202 ft->fcst_unit = 1;
1203 ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_HOUR);
1204 ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR);
1205 }
1206 else if ((fcst_secs < MAX1B_FCST_MINS) &&
1207 ((fcst_secs-accum_period)%SECS_IN_MIN == 0) &&
1208 (fcst_secs%SECS_IN_MIN == 0))
1209 {
1210 ft->time_range = 4;
1211 ft->fcst_unit = 0;
1212 ft->P1 = (fcst_secs-accum_period)/SECS_IN_MIN;
1213 ft->P2 = fcst_secs/SECS_IN_MIN;
1214 }
1215 else if (fcst_secs < MAX1B_FCST_SECS)
1216 {
1217 ft->time_range = 4;
1218 ft->fcst_unit = 254;
1219 ft->P1 = fcst_secs-accum_period;
1220 ft->P2 = fcst_secs;
1221 }
1222 else if ((fcst_secs < MAX1B_FCST_5MINS) &&
1223 (fcst_secs%(SECS_IN_MIN*MINS_IN_5MINS) == 0) &&
1224 (accum_period%(SECS_IN_MIN*MINS_IN_5MINS) == 0))
1225 {
1226 ft->time_range = 4;
1227 ft->fcst_unit = 50;
1228 ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_5MINS);
1229 ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS);
1230 }
1231 else
1232 {
1233 ft->time_range = 255;
1234 ft->fcst_unit = 0;
1235 ft->P1 = 0;
1236 ft->P2 = 0;
1237
1238 ft->fcst_unit_ext_1 = 254;
1239 ft->fcst_unit_ext_2 = 254;
1240 ft->P1_ext = fcst_secs - accum_period;
1241 ft->P2_ext = accum_period;
1242 ft->time_range_ext = 203; /* Duration */
1243 }
1244 }
1245 return 0;
1246 }
1247
1248
1249 /******************************************************************************
1250 * returns a byt from an input integer
1251 *****************************************************************************/
1252
1253 int get_byte(int input_int, int bytenum)
1254 {
1255 int out;
1256 out = ((input_int >> (bytenum-1)*8) & ~(~0 <<8));
1257 return out;
1258 }
1259
1260 /*************************************************************************
1261 * Converts from WRF time format to time format required by grib routines
1262 *************************************************************************/
1263 int grib_time_format(char *DateStr, char *DateStrIn)
1264 {
1265 int year,month,day,hour,minute,second;
1266
1267 trim(DateStrIn);
1268 if (DateStrIn[0] == '*') {
1269 strcpy(DateStr,"*");
1270 }
1271 else
1272 {
1273 sscanf(DateStrIn,"%04d-%02d-%02d_%02d:%02d:%02d",
1274 &year,&month,&day,&hour,&minute,&second);
1275 sprintf(DateStr,"%04d%02d%02d%02d%02d%02d",
1276 year,month,day,hour,minute,second);
1277 }
1278
1279 return 0;
1280 }