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