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