1 /**************************************************************************
2 * Todd Hutchinson 4/20/98
3 * tahutchinson@tasc.com (781) 942-2000 x3108
5 * 55 Walkers Brook Drive
8 * Functions in this file are used for decoding grib data. Please see the
9 * headers before each function for a full descrption.
11 * Routines in this file call functions in the Naval Research Lab's grib
12 * library. The grib library is freely available from
13 * http://www-mel.nrlmry.navy.mil/cgi-bin/order_grib. This library should
14 * be installed on your system prior to using the routines in this file.
15 * Documentation for this library is available from
16 * Master Environmental Grib Library user's manual
17 * http://mel.dmso.mil/docs/grib.pdf
18 * Note: the standard NRL grib library does not support
19 * "Little-Endian" platforms such as linux. There is a version of the NRL
20 * grib library within the WxPredictor project which does support linux.
22 * This file references the cfortran.h header file to ease the use of calling
23 * this function from a fortran routine. cfortran.h is a header file that
24 * allows for simple machine-independent calls between c and fortran. The
25 * package is available via anonymous ftp at zebra.desy.de.
27 * The grib document "A GUIDE TO THE CODE FORM FM 92-IX Ext. GRIB" may be
28 * useful to your understanding of this code. This document is available
29 * via anonymous ftp from nic.fb4.noaa.gov. Check the readme file in the
30 * root directory for further instructions.
32 ****************************************************************************/
38 #define EARTH_RADIUS 6371.229 /* in km */
39 #define PI 3.141592654
40 #define PI_OVER_180 PI/180.
48 #include "/usr/include/time.h"
53 #include "gribfuncs.h"
54 #include "gribsize.incl"
55 #include "read_grib.h"
57 /* Function Declarations */
59 void remove_element(int array
[],int index
, int size
);
60 int advance_time(int *century
, int year
, int month
, int day
, int hour
,
61 int amount
, int unit
);
62 char *advance_time_str(char startdatein
[], int amount
, char enddate
[]);
63 int date_diff(int date1
,int century1
,int date2
,int century2
);
64 int hours_since_1900(int date
,int century
);
65 int isLeapYear(int year
);
66 int get_factor2(int unit
);
67 int compare_record(GribInfo
*gribinfo
, FindGrib
*findgrib
, int gribnum
);
70 *These lines allow fortran routines to call the c routines. They are
71 * used by macros defined in cfortran.h
73 #define get_pressure_levels_STRV_A1 TERM_CHARS(' ',1)
75 FCALLSCFUN6(INT, get_pressure_levels,GET_PRESSURE_LEVELS,
76 get_pressure_levels,STRINGV,INTV,INTV,INTV,INT,INT)
77 #define setup_gribinfo_STRV_A1 TERM_CHARS(' ',1)
78 FCALLSCFUN2(INT,setup_gribinfo,SETUP_GRIBINFO,setup_gribinfo,STRINGV,INT)
79 #define get_msl_indices_STRV_A1 TERM_CHARS(' ',1)
80 FCALLSCFUN9(INT, get_msl_indices,GET_MSL_INDICES,get_msl_indices,
81 STRINGV,INTV,INTV,INTV,INTV,INTV,INT,INTV,INTV)
82 FCALLSCFUN5(INT, get_index,GET_INDEX,get_index,INT,INT,INT,INT,INT)
83 #define read_grib_STRV_A1 TERM_CHARS(' ',1)
84 FCALLSCFUN7(INT,get_dates,GET_DATES,get_dates,INTV,INTV,INTV,INT,INTV,
86 FCALLSCFUN7(INT, read_grib,READ_GRIB,read_grib,
87 STRINGV,INT,INT,INT,INT,FLOATVV,PVOID)
88 FCALLSCFUN8(INT, get_index_near_date,GET_INDEX_NEAR_DATE,get_index_near_date,
89 STRING,INT,INT,INT,INTV,INTV,INTV,INT)
91 /* The value for usLevel_id for isobaric levels */
92 #define ISOBARIC_LEVEL_ID 100
94 /*************************************************************************
95 * This function reads and decodes grib records in a list of input files
96 * and stores information about each grib record in the gribinfo array
97 * structure. The gribinfo structure can then be accessed by any function
102 * gribinfo - pointer to a previously allocated gribinfo structure. The
103 * gribinfo structure is filled in this function.
104 * files - a string array containing the names of the files containing
105 * the grib data. If called from a fortran routine, the
106 * fortran routine must set the character size of the array to
107 * be STRINGSIZE-1. The last filled element in the array should
109 * use_fcst - if TRUE, forecast fields will be included in the gribinfo
110 * structure, otherwise, only analysis fields will be included.
113 * 1 - successful call to setup_gribinfo
114 * -1 - call to setup_gribinfo failed
116 ***************************************************************************/
118 int rg_setup_gribinfo(GribInfo
*gribinfo
, char files
[][STRINGSIZE
],
128 /* Loop through input files */
130 while ((strcmp(files
[filenum
], "end") != 0 ) &&
131 (strcmp(files
[filenum
], "END") != 0 )) {
134 * This forces gribinfo to be fully initialized.
138 gribinfo
->num_elements
= 0;
141 start_elem
= gribinfo
->num_elements
;
143 fp
= fopen(files
[filenum
],"r");
146 fprintf(stderr
,"Could not open %s\n",files
[filenum
]);
151 status
= rg_setup_gribinfo_f(gribinfo
, fp
, use_fcst
);
155 "rg_setup_gribinfo_f returned non-zero status (%d), skipping %s\n",
156 status
,files
[filenum
]);
160 for (idx
=start_elem
; idx
< gribinfo
->num_elements
; idx
++)
162 strcpy(gribinfo
->elements
[idx
].filename
,
176 /*************************************************************************
178 * Similar to rg_setup_gribinfo, except, a unix file descriptor is passed in,
179 * rather than a list of files to open.
181 *************************************************************************/
183 int rg_setup_gribinfo_i(GribInfo
*gribinfo
, int fid
, int use_fcst
)
188 fp
= fdopen(fid
,"r");
191 fprintf(stderr
,"Could not open file descriptor %d\n",fid
);
196 /* This forces gribinfo to be initialized for the first time */
197 gribinfo
->num_elements
= 0;
199 status
= rg_setup_gribinfo_f(gribinfo
, fp
, use_fcst
);
203 "rg_setup_gribinfo_f returned non-zero status (%d)\n",
210 /*************************************************************************
212 * Similar to rg_setup_gribinfo, except, a file pointer is passed in, rather
213 * than a list of files to open.
215 * If gribinfo->num_elements is 0, gribinfo is initialized, otherwise,
216 * gribinfo is appended to.
218 *************************************************************************/
220 int rg_setup_gribinfo_f(GribInfo
*gribinfo
, FILE *fp
, int use_fcst
)
222 char errmsg
[ERRSIZE
];
235 /* Set the number of elements to be zero initially */
236 if (gribinfo
->num_elements
<= 0)
238 /* Allocate space for gribinfo */
239 gribinfo
->elements
= (Elements
*)calloc(ALLOCSIZE
,sizeof(Elements
));
240 if (gribinfo
->elements
== NULL
) {
241 sprintf(errmsg
,"Could not allocate %d bytes for gribinfo->elements\n",
242 ALLOCSIZE
*sizeof(Elements
));
247 /* Make storage for Grib Header */
248 nReturn
= init_gribhdr(&gh1
, errmsg
);
250 * The grib library is setup such that, when init_gribhdr == 0, it was
251 * successful. If it is 1, it failed.
253 if (nReturn
== 1) goto bail_out
;
255 /* Scan through message */
256 for (offset
= 0L; nReturn
== 0; offset
+= gh1
->msg_length
) {
257 if ((gribinfo
->num_elements
> 0) &&
258 (gribinfo
->num_elements
%ALLOCSIZE
== 0))
260 (Elements
*)realloc(gribinfo
->elements
,
261 (gribinfo
->num_elements
+ALLOCSIZE
)*
264 if (gribinfo
->elements
== NULL
) {
265 sprintf(errmsg
,"Could not allocate %d bytes for gribinfo\n",
266 (gribinfo
->num_elements
+ ALLOCSIZE
)*sizeof(Elements
));
270 /* Setup the File pointer */
271 gribinfo
->elements
[gribinfo
->num_elements
].fp
= fp
;
273 gribinfo
->elements
[gribinfo
->num_elements
].pds
=
274 (PDS_INPUT
*)malloc(1*sizeof(PDS_INPUT
));
275 gribinfo
->elements
[gribinfo
->num_elements
].gds
=
276 (grid_desc_sec
*)malloc(1*sizeof(grid_desc_sec
));
277 gribinfo
->elements
[gribinfo
->num_elements
].bms
=
278 (BMS_INPUT
*)malloc(1*sizeof(BMS_INPUT
));
279 gribinfo
->elements
[gribinfo
->num_elements
].bds_head
=
280 (BDS_HEAD_INPUT
*)malloc(1*sizeof(BDS_HEAD_INPUT
));
283 grib_fseek(fp
,&offset
, Rd_Indexfile
, gh1
, errmsg
);
285 if (nReturn
== 2) break; /* End of file error */
287 fprintf(stderr
, "Grib_fseek returned non zero status (%d)\n",
292 if (errmsg
[0] != '\0')
293 { /* NO errors but got a Warning msg from seek */
294 fprintf(stderr
,"%s; Skip Decoding...\n",errmsg
);
296 gh1
->msg_length
= 1L; /* set to 1 to bump offset up */
300 if (gh1
->msg_length
< 0) {
301 fprintf(stderr
, "Error: message returned had bad length (%ld)\n",
305 else if (gh1
->msg_length
== 0) {
306 fprintf(stderr
, "msg_length is Zero\n");
307 gh1
->msg_length
= 1L;
310 init_dec_struct(gribinfo
->elements
[gribinfo
->num_elements
].pds
,
311 gribinfo
->elements
[gribinfo
->num_elements
].gds
,
312 gribinfo
->elements
[gribinfo
->num_elements
].bms
,
313 gribinfo
->elements
[gribinfo
->num_elements
].bds_head
);
316 * gribgetpds is an undocumented function within the grib library.
317 * gribgetpds grabs the pds section from the grib message without
318 * decoding the entire grib message. The interface is as follows:
319 * first input param: a pointer to the beginning of the pds
321 * second input param: a pointer to a structure which will hold
322 * the pds information
323 * third param: the error message.
325 * If gribgetpds ever fails, it can be replaced with the following
326 * nReturn = grib_dec((char *)gh1->entire_msg, &pds, &gds, &bds_head,
327 * &bms, &grib_data, errmsg);
329 * This will degrade performance since this grib_dec decodes the
330 * entire grib message.
333 nReturn
= gribgetpds((char*)(gh1
->entire_msg
+ 8),
334 gribinfo
->elements
[gribinfo
->num_elements
].pds
,
336 if (nReturn
!= 0) goto bail_out
;
338 /* Get gds if present */
339 if (gribinfo
->elements
[gribinfo
->num_elements
].pds
->usGds_bms_id
>> 7
344 gribinfo
->elements
[gribinfo
->num_elements
].pds
->uslength
),
345 gribinfo
->elements
[gribinfo
->num_elements
].gds
,errmsg
);
346 if (nReturn
!= 0) goto bail_out
;
349 /* Get bms section if present */
350 if (gribinfo
->elements
[gribinfo
->num_elements
].pds
->usGds_bms_id
>> 6
353 fprintf(stderr,"grids with bms section not currently supported\n");
358 gribinfo
->elements
[gribinfo
->num_elements
].usGrid_id
=
359 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usGrid_id
;
360 gribinfo
->elements
[gribinfo
->num_elements
].usParm_id
=
361 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usParm_id
;
362 gribinfo
->elements
[gribinfo
->num_elements
].usLevel_id
=
363 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usLevel_id
;
364 gribinfo
->elements
[gribinfo
->num_elements
].usHeight1
=
365 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usHeight1
;
366 gribinfo
->elements
[gribinfo
->num_elements
].usHeight2
=
367 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usHeight2
;
368 gribinfo
->elements
[gribinfo
->num_elements
].center_id
=
369 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usCenter_id
;
370 gribinfo
->elements
[gribinfo
->num_elements
].parmtbl
=
371 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usParm_tbl
;
372 gribinfo
->elements
[gribinfo
->num_elements
].proc_id
=
373 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usProc_id
;
374 gribinfo
->elements
[gribinfo
->num_elements
].subcenter_id
=
375 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usCenter_sub
;
376 gribinfo
->elements
[gribinfo
->num_elements
].offset
= offset
;
377 gribinfo
->elements
[gribinfo
->num_elements
].end
=
378 offset
+ gh1
->msg_length
- 1;
381 century
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usCentury
;
383 if (gribinfo
->elements
[gribinfo
->num_elements
].pds
->usTime_range
== 10)
385 fcsttime1
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP1
*256 +
386 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP2
;
389 else if (gribinfo
->elements
[gribinfo
->num_elements
].pds
->usTime_range
391 /* This is the WSI extension to grib. 203 indicates "duration" */
392 fcsttime1
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP1
;
393 fcsttime2
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP1
+
394 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP2
;
396 fcsttime1
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP1
;
397 fcsttime2
= gribinfo
->elements
[gribinfo
->num_elements
].pds
->usP2
;
400 gribinfo
->elements
[gribinfo
->num_elements
].date
=
401 advance_time(¢ury
,
402 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usYear
,
403 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usMonth
,
404 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usDay
,
405 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usHour
,
407 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usFcst_unit_id
);
410 gribinfo
->elements
[gribinfo
->num_elements
].date
=
411 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usHour
*1 +
412 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usDay
*100 +
413 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usMonth
*10000 +
414 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usYear
*1000000;
416 gribinfo
->elements
[gribinfo
->num_elements
].century
=
417 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usCentury
;
420 (gribinfo
->elements
[gribinfo
->num_elements
].pds
->usCentury
- 1) * 100
421 + gribinfo
->elements
[gribinfo
->num_elements
].pds
->usYear
;
423 sprintf(gribinfo
->elements
[gribinfo
->num_elements
].initdate
,
424 "%04d%02d%02d%02d%02d%02d",
426 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usMonth
,
427 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usDay
,
428 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usHour
,
429 gribinfo
->elements
[gribinfo
->num_elements
].pds
->usMinute
,
433 get_factor2(gribinfo
->elements
[gribinfo
->num_elements
].pds
->usFcst_unit_id
);
434 gribinfo
->elements
[gribinfo
->num_elements
].fcsttime1
=
436 gribinfo
->elements
[gribinfo
->num_elements
].fcsttime2
=
439 advance_time_str(gribinfo
->elements
[gribinfo
->num_elements
].initdate
,
440 gribinfo
->elements
[gribinfo
->num_elements
].fcsttime1
,
441 gribinfo
->elements
[gribinfo
->num_elements
].valid_time
);
443 gribinfo
->num_elements
++;
449 /* The error condition */
451 if (errmsg
[0] != '\0') fprintf(stderr
,"\n***ERROR: %s: %s\n",
452 "setup_grib",errmsg
);
453 if (gribinfo
->elements
!= NULL
) free(gribinfo
->elements
);
454 perror("System Error ");
458 /*****************************************************************************
460 * Retrieve pressure levels from grib data. This function will pass the
461 * pressure levels for which the input parameter is available at all input
462 * times back to the calling routine.
466 * gribinfo - pointer to a previously allocated gribinfo structure.
467 * The gribinfo structure is filled in this function.
468 * dates: an array of dates to check for data
470 * If called from a fortran routine, the fortran routine must
471 * set the character size of the array to be STRINGSIZE-1
472 * centuries: an array holding the centuries for each of the
473 * dates in the array dates.
474 * parm_id: the input parameter id. From table 2 of the grib manual.
476 * finallevels: an array of pressure levels which are contained in
477 * the grib data at all input times.
479 * the number of levels in the levels array. The levels are listing
480 * in descending (by value) order, i.e., the value with the highest
481 * pressure (lowest vertical level) is the first element.
483 ****************************************************************************/
485 int rg_get_pressure_levels(GribInfo
*gribinfo
, int dates
[], int centuries
[],
486 int parm_id
[], int finallevels
[],int min_pres
,
496 int numfinallevels
= 0;
501 levelnum
= (int *)calloc(numparms
,sizeof(int));
502 tmplevels
= (int **)calloc(numparms
,sizeof(int *));
503 for (j
= 0; j
< numparms
; j
++) {
504 tmplevels
[j
] = (int *)calloc(1000,sizeof(int));
505 if (tmplevels
[j
] == NULL
) {
510 if ((levelnum
== NULL
) || (tmplevels
== NULL
)) {
512 "get_pressure_levels: Allocation of space failed, returning\n");
516 /* Loop through all parameters */
517 for (parmnum
= 0; parmnum
< numparms
; parmnum
++) {
519 levelnum
[parmnum
] = 0;
521 /* Get the array of pressure levels available at the first input time */
523 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
524 if (gribinfo
->elements
[gribnum
].date
== dates
[datenum
]) {
525 if (gribinfo
->elements
[gribnum
].century
== centuries
[datenum
]) {
526 if (gribinfo
->elements
[gribnum
].usLevel_id
== ISOBARIC_LEVEL_ID
) {
527 if (gribinfo
->elements
[gribnum
].usParm_id
== parm_id
[parmnum
]) {
528 if (gribinfo
->elements
[gribnum
].usHeight1
>= min_pres
) {
530 for (j
=0; j
< levelnum
[parmnum
]; j
++) {
531 if (tmplevels
[parmnum
][j
] ==
532 gribinfo
->elements
[gribnum
].usHeight1
) {
537 if (levelincluded
== 0) {
538 tmplevels
[parmnum
][levelnum
[parmnum
]] =
539 gribinfo
->elements
[gribnum
].usHeight1
;
549 /* Remove levels that are not contained at all subsequent times */
551 while (dates
[datenum
] != -99){
552 for (j
= 0; j
< levelnum
[parmnum
]; j
++) {
554 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
555 if (gribinfo
->elements
[gribnum
].date
== dates
[datenum
]) {
556 if (gribinfo
->elements
[gribnum
].century
== centuries
[datenum
]) {
557 if (gribinfo
->elements
[gribnum
].usLevel_id
== ISOBARIC_LEVEL_ID
)
559 if (gribinfo
->elements
[gribnum
].usParm_id
==
561 if (tmplevels
[parmnum
][j
] ==
562 gribinfo
->elements
[gribnum
].usHeight1
)
569 if (!(contains_level
)) {
570 remove_element(tmplevels
[parmnum
],j
,levelnum
[parmnum
]);
579 * Put the values for levels into an array. Remove any levels that
580 * were not found at all other levels
583 for (j
= 0; j
< levelnum
[parmnum
]; j
++) {
584 finallevels
[j
] = tmplevels
[parmnum
][j
];
588 for (i
=0; i
<numfinallevels
; i
++) {
590 for (j
=0; j
<levelnum
[parmnum
]; j
++) {
591 if (finallevels
[i
] == tmplevels
[parmnum
][j
]) {
596 if (!contains_level
) {
597 remove_element(finallevels
,i
,numfinallevels
);
607 * Sort the numfinallevels array into descending order. Use straight
610 for (j
=1; j
<numfinallevels
; j
++) {
611 tmpval
= finallevels
[j
];
612 for (i
=j
-1; i
>= 0; i
--) {
613 if (finallevels
[i
] >= tmpval
) break;
614 finallevels
[i
+1] = finallevels
[i
];
616 finallevels
[i
+1] = tmpval
;
619 return numfinallevels
;
622 /****************************************************************************
624 * Returns an array of grib indices that correspond to particular grib fields
625 * to use as sea level pressure. There will be exactly one element for each
626 * input time. If a field was not found, then this function returns NULL
630 * gribinfo - pointer to a previously allocated gribinfo structure. The
631 * gribinfo structure is filled in this function.
632 * dates: a string array of dates to check for data.
634 * If called from a fortran routine, the fortran routine must
635 * set the character size of the array to be STRINGSIZE-1
636 * centuries: an array holding the centuries for each of the
637 * dates in the array dates.
638 * usParm_id: an array of parameter identifiers that could be
639 * used as a sea level pressure field (From table 2 of
640 * grib documentation)
641 * usLevel_id: the level id that could be used as a sea level pressure
642 * field (from table 3 of the grib documentation)
643 * usHeight1: the height for the particular parameter and level
644 * (in units described by the parameter index)
645 * numparms: the number of parameters in each of the usParm_id,
646 * usLevel_id, and usHeight1 arrays.
648 * grib_index: an array of grib indices to use for the sea level
649 * pressure. The index to grib_index corresponds to
650 * the time, i.e., the first array element of grib_index
651 * corresponds to the first time, the second element to
652 * the second time, etc.
654 * Note: Values in the input arrays, usParm_id, usLevel_id, and
655 * usHeight with the same array index must correspond.
659 * -1 if no field was found.
660 ***************************************************************************/
662 int rg_get_msl_indices(GribInfo
*gribinfo
, char dates
[][STRINGSIZE
],
663 int centuries
[], int usParm_id
[],int usLevel_id
[],
664 int usHeight1
[],int infactor
[],int numparms
,
665 int grib_index
[],int outfactor
[])
672 for (parmindex
= 0; parmindex
< numparms
; parmindex
++) {
675 while ((strcmp(dates
[datenum
], "end") != 0 ) &&
676 (strcmp(dates
[datenum
], "END") != 0 )) {
678 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
679 if (gribinfo
->elements
[gribnum
].date
== atoi(dates
[datenum
])) {
680 if (gribinfo
->elements
[gribnum
].century
== centuries
[datenum
]) {
681 if ((gribinfo
->elements
[gribnum
].usParm_id
==
682 usParm_id
[parmindex
]) &&
683 (gribinfo
->elements
[gribnum
].usLevel_id
==
684 usLevel_id
[parmindex
]) &&
685 (gribinfo
->elements
[gribnum
].usHeight1
==
686 usHeight1
[parmindex
])) {
687 grib_index
[datenum
] = gribnum
;
688 outfactor
[datenum
] = infactor
[parmindex
];
699 * Break out of loop and continue on to next parameter if the current
700 * parameter was missing from a date.
703 if (foundfield
!= datenum
) break;
707 * Break out of the parameter loop once we've found a field available at all
710 if (foundfield
== datenum
) {
716 if (foundfield
== datenum
)
724 /***************************************************************************
726 * This function takes an index as input and returns a 2d grib data field
730 * gribinfo - pointer to a previously allocated gribinfo structure. The
731 * gribinfo structure is filled in this function.
732 * index - the index of gribinfo for which data is to be retrieved
733 * scale - the scale factor to multiply data by, i.e., if -2,
734 * data will be multiplied by 10^-2.
736 * grib_out - the 2 dimensional output grib data
737 * Warning: This 2d array is setup with i being the vertical
738 * dimension and j being the horizontal dimension. This
739 * is the convention used in mesoscale numerical modeling
740 * (the MM5 in particular), so it is used here.
744 ***************************************************************************/
746 int rg_get_grib(GribInfo
*gribinfo
, int index
,int scale
,
747 float **grib_out
,int *vect_comp_flag
,
748 GRIB_PROJECTION_INFO_DEF
*Proj
, BDS_HEAD_INPUT
*bds_head
)
750 char errmsg
[ERRSIZE
];
756 BDS_HEAD_INPUT dummy
;
763 /* Initialize Variables */
766 grib_data
= (float *)NULL
;
768 /* Make storage for Grib Header */
769 nReturn
= init_gribhdr (&gh1
, errmsg
);
770 if (nReturn
== 1) goto bail_out
;
772 /* Seek to the position in the grib data */
773 offset
= gribinfo
->elements
[index
].offset
;
774 nReturn
= grib_fseek(gribinfo
->elements
[index
].fp
,&offset
,
775 Rd_Indexfile
,gh1
,errmsg
);
777 fprintf(stderr
,"Grib_fseek returned error status (%d)\n",nReturn
);
780 if (errmsg
[0] != '\0')
781 { /* NO errors but got a Warning msg from seek */
782 fprintf(stderr
,"%s: Skip Decoding...\n",errmsg
);
785 if (gh1
->msg_length
<= 0) {
786 fprintf(stderr
,"Error: message returned had bad length (%ld)\n",
790 init_dec_struct(&pds
, &gds
, &bms
, &dummy
);
792 nReturn
= grib_dec((char *)gh1
->entire_msg
, &pds
, &gds
,
794 &bms
, &grib_data
, errmsg
);
796 if (nReturn
!= 0) goto bail_out
;
798 if (bms
.uslength
> 0) {
799 nReturn
= apply_bitmap(&bms
, &grib_data
, FILL_VALUE
, bds_head
,
801 if (nReturn
!= 0) goto bail_out
;
804 switch(gds
.head
.usData_type
) {
807 strcpy(Proj
->prjnmm
,"latlon");
808 Proj
->colcnt
= gds
.llg
.usNi
;
809 Proj
->rowcnt
= gds
.llg
.usNj
;
810 Proj
->origlat
= gds
.llg
.lLat1
/1000.;
811 Proj
->origlon
= gds
.llg
.lLon1
/1000.;
812 Proj
->xintdis
= (gds
.llg
.iDi
/1000.)*EARTH_RADIUS
*PI_OVER_180
;
813 Proj
->yintdis
= (gds
.llg
.iDj
/1000.)*EARTH_RADIUS
*PI_OVER_180
;
816 if ((gds
.llg
.usRes_flag
>> 3) & 1) *vect_comp_flag
= 1;
817 else *vect_comp_flag
= 0;
819 /* If the grid is a global grid, we want to set the expandlon flag
820 * so that the number of columns in the array is expanded by one and
821 * the first column of data is copied to the last column. This
822 * allows calling routines to interpolate between first and last columns
826 if (gds
.llg
.usNi
*gds
.llg
.iDi
/1000. == 360)
833 strcpy(Proj
->prjnmm
,"mercator");
834 Proj
->colcnt
= gds
.merc
.cols
;
835 Proj
->rowcnt
= gds
.merc
.rows
;
836 Proj
->origlat
= gds
.merc
.first_lat
/1000.;
837 Proj
->origlon
= gds
.merc
.first_lon
/1000.;
838 Proj
->xintdis
= gds
.merc
.lon_inc
/1000.;
839 Proj
->yintdis
= gds
.merc
.lat_inc
/1000.;
840 Proj
->parm1
= gds
.merc
.latin
/1000.;
841 Proj
->parm2
= (gds
.merc
.Lo2
/1000. - Proj
->origlon
)/gds
.merc
.cols
;
842 if ((gds
.merc
.usRes_flag
>> 3) & 1) *vect_comp_flag
= 1;
843 else *vect_comp_flag
= 0;
846 strcpy(Proj
->prjnmm
,"lambert");
847 Proj
->colcnt
= gds
.lam
.iNx
;
848 Proj
->rowcnt
= gds
.lam
.iNy
;
849 Proj
->origlat
= gds
.lam
.lLat1
/1000.;
850 Proj
->origlon
= gds
.lam
.lLon1
/1000.;
851 Proj
->xintdis
= gds
.lam
.ulDx
/1000.;
852 Proj
->yintdis
= gds
.lam
.ulDy
/1000.;
853 Proj
->parm1
= gds
.lam
.lLat_cut1
/1000.;
854 Proj
->parm2
= gds
.lam
.lLat_cut2
/1000.;
855 Proj
->parm3
= gds
.lam
.lLon_orient
/1000.;
856 if ((gds
.lam
.usRes_flag
>> 3) & 1) *vect_comp_flag
= 1;
857 else *vect_comp_flag
= 0;
860 strcpy(Proj
->prjnmm
,"polar_stereo");
861 Proj
->colcnt
= gds
.pol
.usNx
;
862 Proj
->rowcnt
= gds
.pol
.usNy
;
863 Proj
->origlat
= gds
.pol
.lLat1
/1000.;
864 Proj
->origlon
= gds
.pol
.lLon1
/1000.;
865 Proj
->xintdis
= gds
.pol
.ulDx
/1000.;
866 Proj
->yintdis
= gds
.pol
.ulDy
/1000.;
868 Proj
->parm2
= gds
.pol
.lLon_orient
/1000.;
869 if ((gds
.pol
.usRes_flag
>> 3) & 1) *vect_comp_flag
= 1;
870 else *vect_comp_flag
= 0;
873 fprintf(stderr
,"Grid not supported, gds.head.usData_type = %d\n",
874 gds
.head
.usData_type
);
875 fprintf(stderr
,"Exiting\n");
880 strcpy(Proj
->stordsc
,"+y_in_+x");
884 for (j
=0; j
< (Proj
->rowcnt
); j
++) {
885 for (i
=0; i
<(Proj
->colcnt
); i
++) {
886 grib_out
[j
][i
] = grib_data
[i
+j
*Proj
->colcnt
]*pow(10,scale
);
892 for (j
= 0; j
< Proj
->rowcnt
; j
++) {
893 grib_out
[j
][Proj
->colcnt
-1] = grib_out
[j
][0];
898 * You only reach here when there is no error, so return successfully.
903 if (grib_data
!= NULL
) {
910 /* The error condition */
912 if (errmsg
[0] != '\0') fprintf(stderr
,"\n***ERROR: %s %s\n",
914 if (grib_data
!= NULL
)
920 /***************************************************************************
922 * This function takes an index as input and returns a 2d grib data field
926 * gribinfo - pointer to a previously allocated gribinfo structure. The
927 * gribinfo structure is filled in this function.
928 * index - the index of gribinfo for which data is to be retrieved
930 * data - the 2 dimensional output grib data
931 * Warning: This 2d array is setup with i being the vertical
932 * dimension and j being the horizontal dimension. This
933 * is the convention used in mesoscale numerical modeling
934 * (the MM5 in particular), so it is used here.
938 ***************************************************************************/
940 int rg_get_data(GribInfo
*gribinfo
, int index
, float **data
)
947 numrows
= rg_get_numrows(gribinfo
,index
);
948 numcols
= rg_get_numcols(gribinfo
,index
);
950 data_1d
= (float *)calloc(numrows
*numcols
,sizeof(float));
953 fprintf(stderr
,"Allocating space for data_1d failed, index: %d\n",index
);
957 status
= rg_get_data_1d(gribinfo
, index
, data_1d
);
963 for (j
=0; j
< numrows
; j
++) {
964 for (i
=0; i
< numcols
; i
++) {
965 data
[j
][i
] = data_1d
[i
+j
*numcols
];
975 /***************************************************************************
977 * This function takes an index as input and returns a 1d grib data field
981 * gribinfo - pointer to a previously allocated gribinfo structure. The
982 * gribinfo structure is filled in this function.
983 * index - the index of gribinfo for which data is to be retrieved
985 * data - 1 dimensional output grib data
986 * Warning: This 2d array is setup with i being the vertical
987 * dimension and j being the horizontal dimension. This
988 * is the convention used in mesoscale numerical modeling
989 * (the MM5 in particular), so it is used here.
993 ***************************************************************************/
995 int rg_get_data_1d(GribInfo
*gribinfo
, int index
, float *data
)
997 char errmsg
[ERRSIZE
];
1003 BDS_HEAD_INPUT bds_head
;
1007 int numcols
, numrows
;
1010 /* Initialize Variables */
1013 grib_data
= (float *)NULL
;
1015 /* Make storage for Grib Header */
1016 nReturn
= init_gribhdr (&gh1
, errmsg
);
1017 if (nReturn
== 1) goto bail_out
;
1019 /* Seek to the position in the grib data */
1020 offset
= gribinfo
->elements
[index
].offset
;
1021 nReturn
= grib_fseek(gribinfo
->elements
[index
].fp
,&offset
,
1022 Rd_Indexfile
,gh1
,errmsg
);
1024 fprintf(stderr
,"Grib_fseek returned non-zero status (%d)\n",nReturn
);
1027 if (errmsg
[0] != '\0')
1028 { /* NO errors but got a Warning msg from seek */
1029 fprintf(stderr
,"%s: Skip Decoding...\n",errmsg
);
1032 if (gh1
->msg_length
<= 0) {
1033 fprintf(stderr
,"Error: message returned had bad length (%ld)\n",
1038 init_dec_struct(&pds
, &gds
, &bms
, &bds_head
);
1040 nReturn
= grib_dec((char *)gh1
->entire_msg
, &pds
, &gds
,
1042 &bms
, &grib_data
, errmsg
);
1044 if (nReturn
!= 0) goto bail_out
;
1046 if (bms
.uslength
> 0) {
1047 nReturn
= apply_bitmap(&bms
, &grib_data
, FILL_VALUE
, &bds_head
,
1049 if (nReturn
!= 0) goto bail_out
;
1053 * Copy the data into the permanent array
1055 numcols
= rg_get_numcols(gribinfo
,index
);
1056 numrows
= rg_get_numrows(gribinfo
,index
);
1057 memcpy(data
,grib_data
,numcols
*numrows
*sizeof(float));
1060 * You only reach here when there is no error, so return successfully.
1065 if (grib_data
!= NULL
) {
1072 /* The error condition */
1074 if (errmsg
[0] != '\0') fprintf(stderr
,"\n***ERROR: %s %s\n",
1076 if (grib_data
!= NULL
)
1082 /****************************************************************************
1083 * Returns the index of gribinfo corresponding to the input date, level,
1084 * height, and parameter.
1088 * gribinfo - pointer to a previously populated gribinfo structure.
1089 * initdate - initialization date in the form yyyymmdd[HHMMSS]. If any
1090 * part of HHMMSS is not specified, it will be set to 0.
1091 * parmid - the parameter id in the grib file
1092 * leveltype - the leveltype id from table 3/3a of the grib document.
1093 * level1 - First level of the data in units described by leveltype.
1094 * level2 - Second level of the data in units described by leveltype.
1095 * fcsttime1 - First forecast time in seconds.
1096 * fcsttime2 - Second forecast time in seconds.
1097 * Note: If an input variable is set set to -INT_MAX, then any value
1098 * will be considered a match.
1100 * if >= 0 The index of the gribinfo data that corresponds to the
1102 * if < 0 No field corresponding to the input parms was found.
1104 ***************************************************************************/
1106 int rg_get_index(GribInfo
*gribinfo
, FindGrib
*findgrib
)
1111 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
1112 if (compare_record(gribinfo
, findgrib
, gribnum
) == 1)
1114 grib_index
= gribnum
;
1122 /****************************************************************************
1123 * Same as rg_get_index, except that a guess for the record number is given.
1124 * This "guess" record is first checked to see if it matches, if so,
1125 * that grib record number is just returned. If it does not match,
1126 * full searching ensues.
1127 * Returns the index of gribinfo corresponding to the input date, level,
1128 * height, and parameter.
1132 * Same is rg_get_index, except:
1133 * guess_index - The index to check first.
1135 * Same as rg_get_index
1137 ***************************************************************************/
1139 int rg_get_index_guess(GribInfo
*gribinfo
, FindGrib
*findgrib
, int guess_index
)
1143 if (compare_record(gribinfo
, findgrib
, guess_index
) == 1) {
1144 retval
= guess_index
;
1146 retval
= rg_get_index(gribinfo
, findgrib
);
1153 /****************************************************************************
1154 * Sets all values in FindGrib to missing.
1158 * findgrib - pointer to a previously allocated findgrib structure.
1164 ***************************************************************************/
1165 int rg_init_findgrib(FindGrib
*findgrib
)
1167 strcpy(findgrib
->initdate
,"*");
1168 strcpy(findgrib
->validdate
,"*");
1169 findgrib
->parmid
= -INT_MAX
;
1170 findgrib
->parmid
= -INT_MAX
;
1171 findgrib
->leveltype
= -INT_MAX
;
1172 findgrib
->level1
= -INT_MAX
;
1173 findgrib
->level2
= -INT_MAX
;
1174 findgrib
->fcsttime1
= -INT_MAX
;
1175 findgrib
->fcsttime2
= -INT_MAX
;
1176 findgrib
->center_id
= -INT_MAX
;
1177 findgrib
->subcenter_id
= -INT_MAX
;
1178 findgrib
->parmtbl_version
= -INT_MAX
;
1183 /****************************************************************************
1184 * Returns the indices of all gribinfo entries that match the input date,
1185 * level, height, and parameter.
1189 * gribinfo - pointer to a previously populated gribinfo structure.
1190 * initdate - initialization date in the form yyyymmdd[HHMMSS]. If any
1191 * part of HHMMSS is not specified, it will be set to 0.
1192 * parmid - the parameter id in the grib file
1193 * leveltype - the leveltype id from table 3/3a of the grib document.
1194 * level1 - First level of the data in units described by leveltype.
1195 * level2 - Second level of the data in units described by leveltype.
1196 * fcsttime1 - First forecast time in seconds.
1197 * fcsttime2 - Second forecast time in seconds.
1198 * indices - an array of indices that match
1199 * num_indices - the number of matches and output indices
1201 * Note: If an input variable is set set to -INT_MAX, then any value
1202 * will be considered a match.
1204 * The number of matching indices.
1206 ***************************************************************************/
1208 int rg_get_indices(GribInfo
*gribinfo
, FindGrib
*findgrib
, int indices
[])
1213 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
1214 if (compare_record(gribinfo
, findgrib
, gribnum
) == 1) {
1215 indices
[matchnum
] = gribnum
;
1222 /*************************************************************************
1224 * Returns an array of dates that correspond to particular input grib fields.
1225 * The dates will be sorted so that the dates increase as the index increases.
1229 * gribinfo - pointer to a previously allocated gribinfo structure.
1230 * The gribinfo structure is filled in this function.
1231 * usParm_id: an array of parameter identifiers that could be
1232 * used as a sea level pressure field (From table 2 of
1233 * grib documentation)
1234 * usLevel_id: the level id that could be used as a sea level pressure
1235 * field (from table 3 of the grib documentation)
1236 * usHeight1: the height for the particular parameter and level
1237 * (in units described by the parameter index)
1238 * numparms: the number of parameters in each of the usParm_id,
1239 * usLevel_id, and usHeight1 arrays.
1241 * dates: the dates for which the input fields are available.
1243 * Note: Values in the input arrays, usParm_id, usLevel_id, and
1244 * usHeight with the same array index must correspond.
1247 * The number of dates found.
1248 *************************************************************************/
1250 int rg_get_dates(GribInfo
*gribinfo
,int usParm_id
[],int usLevel_id
[],
1251 int usHeight1
[],int numparms
,int dates
[],int century
[],
1257 int already_included
;
1259 int tmpval
,tmpval2
,tmpval3
;
1261 /* Get the dates for the given parameters */
1263 for (parmindex
= 0; parmindex
< numparms
; parmindex
++) {
1264 for (gribnum
= 0; gribnum
< gribinfo
->num_elements
; gribnum
++) {
1265 if ((gribinfo
->elements
[gribnum
].usParm_id
== usParm_id
[parmindex
]) &&
1266 (gribinfo
->elements
[gribnum
].usLevel_id
== usLevel_id
[parmindex
]) &&
1267 (gribinfo
->elements
[gribnum
].usHeight1
== usHeight1
[parmindex
])) {
1268 already_included
= 0;
1269 for (i
= 0; i
< datenum
; i
++){
1270 if ((dates
[datenum
] == gribinfo
->elements
[gribnum
].date
) &&
1271 (century
[datenum
] == gribinfo
->elements
[gribnum
].century
)) {
1272 already_included
= 1;
1276 if (!already_included
) {
1277 dates
[datenum
] = gribinfo
->elements
[gribnum
].date
;
1278 century
[datenum
] = gribinfo
->elements
[gribnum
].century
;
1279 indices
[datenum
] = gribnum
;
1286 /* Sort the dates into increasing order */
1287 for (j
= 1; j
< datenum
; j
++) {
1289 tmpval2
= indices
[j
];
1290 tmpval3
= century
[j
];
1291 for (i
=j
-1; i
>= 0; i
--) {
1292 if (dates
[i
] <= tmpval
) break;
1293 dates
[i
+1] = dates
[i
];
1294 indices
[i
+1] = indices
[i
];
1295 century
[i
+1] = century
[i
];
1297 dates
[i
+1] = tmpval
;
1298 indices
[i
+1] = tmpval2
;
1299 century
[i
+1] = tmpval3
;
1305 /****************************************************************************
1306 * This function returns the pds, gds, bms, and bms_head section of the
1310 * gribinfo - pointer to a previously allocated gribinfo structure. The
1311 * gribinfo structure is filled in this function.
1312 * index - the index of the grib record to access as indexed by
1316 * *pds - a pointer to a structure holding the pds information
1317 * *gds - a pointer to a structure holding the gds information
1318 * *bms - a pointer to a structure holding the bms information
1319 * *bds_head - a pointer to a structure holding the binary data section
1320 * header information
1322 ***************************************************************************
1324 int rg_get_grib_header(GribInfo
*gribinfo
, int index
, PDS_INPUT
*pds
,
1325 grid_desc_sec
*gds
,BMS_INPUT
*bms
)
1329 memcpy(pds
,gribinfo
->elements
[index
].pds
,sizeof(PDS_INPUT
));
1330 memcpy(gds
,gribinfo
->elements
[index
].gds
,sizeof(grid_desc_sec
));
1331 memcpy(bms
,gribinfo
->elements
[index
].bms
,sizeof(BMS_INPUT
));
1333 /* Reset the dimensions for thinned grids */
1334 if (gribinfo
->elements
[index
].gds
->head
.thin
!= NULL
) {
1335 if (gds
->head
.thin
!= NULL
) {
1336 if ((gds
->head
.usData_type
== LATLON_PRJ
) ||
1337 (gds
->head
.usData_type
== GAUSS_PRJ
) ||
1338 (gds
->head
.usData_type
== ROT_LATLON_PRJ
) ||
1339 (gds
->head
.usData_type
== ROT_GAUSS_PRJ
) ||
1340 (gds
->head
.usData_type
== STR_LATLON_PRJ
) ||
1341 (gds
->head
.usData_type
== STR_GAUSS_PRJ
) ||
1342 (gds
->head
.usData_type
== STR_ROT_LATLON_PRJ
) ||
1343 (gds
->head
.usData_type
== STR_ROT_GAUSS_PRJ
)) {
1344 ysize
= gds
->llg
.usNj
;
1345 } else if (gds
->head
.usData_type
== MERC_PRJ
) {
1346 ysize
= gds
->merc
.rows
;
1347 } else if (gds
->head
.usData_type
== POLAR_PRJ
) {
1348 ysize
= gds
->pol
.usNy
;
1349 } else if ((gds
->head
.usData_type
== LAMB_PRJ
) ||
1350 (gds
->head
.usData_type
== ALBERS_PRJ
) ||
1351 (gds
->head
.usData_type
== OBLIQ_LAMB_PRJ
)) {
1352 ysize
= gds
->lam
.iNy
;
1356 for (j
= 0; j
<ysize
; j
++) {
1357 if (gds
->head
.thin
[j
] > xsize
) {
1358 xsize
= gds
->head
.thin
[j
];
1363 if ((gds
->head
.usData_type
== LATLON_PRJ
) ||
1364 (gds
->head
.usData_type
== GAUSS_PRJ
) ||
1365 (gds
->head
.usData_type
== ROT_LATLON_PRJ
) ||
1366 (gds
->head
.usData_type
== ROT_GAUSS_PRJ
) ||
1367 (gds
->head
.usData_type
== STR_LATLON_PRJ
) ||
1368 (gds
->head
.usData_type
== STR_GAUSS_PRJ
) ||
1369 (gds
->head
.usData_type
== STR_ROT_LATLON_PRJ
) ||
1370 (gds
->head
.usData_type
== STR_ROT_GAUSS_PRJ
)) {
1371 gds
->llg
.usNi
= xsize
;
1372 gds
->llg
.iDi
= abs(gds
->llg
.lLat2
- gds
->llg
.lLat1
)/(xsize
-1);
1373 } else if (gds
->head
.usData_type
== MERC_PRJ
) {
1374 gds
->merc
.cols
= xsize
;
1375 } else if (gds
->head
.usData_type
== POLAR_PRJ
) {
1376 gds
->pol
.usNx
= xsize
;
1377 } else if ((gds
->head
.usData_type
== LAMB_PRJ
) ||
1378 (gds
->head
.usData_type
== ALBERS_PRJ
) ||
1379 (gds
->head
.usData_type
== OBLIQ_LAMB_PRJ
)) {
1380 gds
->lam
.iNx
= xsize
;
1388 /****************************************************************************
1389 * This returns the index of the gribdata for paramaters which match the input
1390 * parameters and for the date closest to the input targetdate. If dates are
1391 * not found either within hours_before or hours_after the time, then a missing
1392 * value is returned.
1396 * gribinfo - pointer to a previously allocated gribinfo structure. The
1397 * gribinfo structure is filled in this function.
1398 * targetdate: This is the date which dates in the grib data will be
1399 * compared to. (format: integer yymmddhh)
1400 * hours_before: The maximum difference in time prior to the targetdate
1401 * for which data should be searched for.
1402 * hours_after: The maximum difference in time after the targetdate for
1403 * which data should be searched for.
1404 * usParm_id: an array of parameter identifiers that could be
1405 * used as a sea level pressure field (From table 2 of
1406 * grib documentation)
1407 * usLevel_id: the level id that could be used as a sea level pressure
1408 * field (from table 3 of the grib documentation)
1409 * usHeight1: the height for the particular parameter and level
1410 * (in units described by the parameter index)
1411 * numparms: the number of parameters in each of the usParm_id,
1412 * usLevel_id, and usHeight1 arrays.
1414 * the index of the gribdata with a time closest to the target date.
1415 * -1 if there is no time within the input time limits.
1417 ****************************************************************************/
1418 int rg_get_index_near_date(GribInfo
*gribinfo
,char targetdate
[STRINGSIZE
],
1419 int century
,int hours_before
,int hours_after
,
1420 int usParm_id
[],int usLevel_id
[],int usHeight1
[],
1423 int dates
[500],indices
[500],centuries
[500];
1424 int date_before
= MISSING
;
1425 int date_after
= MISSING
;
1426 int century_before
,century_after
;
1427 int date_diff_before
= MISSING
;
1428 int date_diff_after
= MISSING
;
1429 int index_before
,index_after
;
1430 int numdates
,datenum
;
1434 itargetdate
= atoi(targetdate
);
1436 numdates
= rg_get_dates(gribinfo
,usParm_id
,usLevel_id
,usHeight1
,numparms
,
1437 dates
,centuries
,indices
);
1438 if (numdates
<= 0) {
1439 fprintf(stderr
,"get_index_near_date: No dates were found\n");
1443 for (datenum
= 0; datenum
< numdates
; datenum
++) {
1444 if ((dates
[datenum
] > itargetdate
) && (centuries
[datenum
] >= century
)) {
1445 century_after
= centuries
[datenum
];
1446 date_after
= dates
[datenum
];
1447 index_after
= indices
[datenum
];
1450 century_before
= centuries
[datenum
];
1451 date_before
= dates
[datenum
];
1452 index_before
= indices
[datenum
];
1456 if (date_after
!= MISSING
)
1457 date_diff_after
= date_diff(date_after
,century_after
,itargetdate
,century
);
1458 if (date_before
!= MISSING
)
1460 date_diff(itargetdate
,century
,date_before
,century_before
);
1462 if ((date_after
!= MISSING
) && (date_before
!= MISSING
)) {
1463 if ((date_diff_after
<= hours_after
) &&
1464 (date_diff_before
<= hours_before
)) {
1465 if (date_diff_after
< date_diff_before
)
1466 index
= index_before
;
1468 index
= index_after
;
1469 } else if (date_diff_after
<= hours_after
) {
1470 index
= index_after
;
1471 } else if (date_diff_before
<= hours_before
) {
1472 index
= index_before
;
1476 } else if (date_after
!= MISSING
) {
1477 if (date_diff_after
<= hours_after
)
1478 index
= index_after
;
1481 } else if (date_before
!= MISSING
) {
1482 if (date_diff_before
<= hours_before
)
1483 index
= index_before
;
1494 /*****************************************************************************
1496 * returns valid time ( = init time + forecast time)
1499 * gribinfo - pointer to a previously allocated gribinfo structure. The
1500 * gribinfo structure is filled in this function.
1501 * index - index number of record to get valid time from
1504 * valid_time - yyyymmddhhmmss
1510 *****************************************************************************/
1511 int rg_get_valid_time(GribInfo
*gribinfo
, int index
, char valid_time
[])
1513 strcpy(valid_time
, gribinfo
->elements
[index
].valid_time
);
1517 /*****************************************************************************
1519 * returns generating center id
1522 * gribinfo - pointer to a previously allocated gribinfo structure. The
1523 * gribinfo structure is filled in this function.
1524 * index - index number of record to get valid time from
1527 * generating center id
1530 *****************************************************************************/
1531 int rg_get_center_id(GribInfo
*gribinfo
, int index
)
1533 return gribinfo
->elements
[index
].center_id
;
1536 /*****************************************************************************
1538 * returns parameter table version number
1541 * gribinfo - pointer to a previously allocated gribinfo structure. The
1542 * gribinfo structure is filled in this function.
1543 * index - index number of record to get valid time from
1546 * parameter table version number
1549 *****************************************************************************/
1550 int rg_get_parmtbl(GribInfo
*gribinfo
, int index
)
1552 return gribinfo
->elements
[index
].parmtbl
;
1555 /*****************************************************************************
1557 * returns generating process id
1560 * gribinfo - pointer to a previously allocated gribinfo structure. The
1561 * gribinfo structure is filled in this function.
1562 * index - index number of record to get valid time from
1565 * generating process id
1568 *****************************************************************************/
1569 int rg_get_proc_id(GribInfo
*gribinfo
, int index
)
1571 return gribinfo
->elements
[index
].proc_id
;
1574 /*****************************************************************************
1576 * returns sub center id
1579 * gribinfo - pointer to a previously allocated gribinfo structure. The
1580 * gribinfo structure is filled in this function.
1581 * index - index number of record to get valid time from
1587 *****************************************************************************/
1588 int rg_get_subcenter_id(GribInfo
*gribinfo
, int index
)
1590 return gribinfo
->elements
[index
].subcenter_id
;
1593 /**************************************************************************
1595 * Interpolates grib grid data to a point location.
1599 * gribinfo - pointer to a previously allocated gribinfo structure. The
1600 * gribinfo structure is filled in this function.
1601 * index - the index of gribinfo for which data is to be retrieved.
1602 * the first grib record is number 1.
1603 * column - the column of the point in grid coordinates (can be
1604 * floating point number). leftmost column is 1.
1605 * row - the row of the point in grid coordinates (can be
1606 * floating point number). bottommost row is 1.
1609 * on success - the interpolated value at the column,row location.
1610 * on failure - -99999
1612 ***************************************************************************/
1614 float rg_get_point(GribInfo
*gribinfo
, int index
, float column
, float row
)
1617 GRIB_PROJECTION_INFO_DEF Proj
;
1618 BDS_HEAD_INPUT bds_head
;
1622 int numrows
, numcols
;
1623 int top
, left
, right
, bottom
;
1626 numrows
= rg_get_numrows(gribinfo
, index
);
1627 numcols
= rg_get_numcols(gribinfo
, index
);
1629 grib_out
= (float **)alloc_float_2d(numrows
,numcols
);
1630 if (grib_out
== NULL
) {
1631 fprintf(stderr
,"rg_get_point: Could not allocate space for grib_out\n");
1635 status
= rg_get_data(gribinfo
, index
, grib_out
);
1637 fprintf(stderr
,"rg_get_point: rg_get_data failed\n");
1641 /* Do the interpolation here */
1642 bottom
= floor(row
);
1644 left
= floor(column
);
1645 right
= floor(column
+1);
1647 y1
= (row
- bottom
) * (grib_out
[top
][left
] - grib_out
[bottom
][left
]) +
1648 grib_out
[bottom
][left
];
1649 y2
= (row
- bottom
) * (grib_out
[top
][right
] - grib_out
[bottom
][right
]) +
1650 grib_out
[bottom
][right
];
1651 outval
= (y2
- y1
) * (column
- left
) + y1
;
1653 free_float_2d(grib_out
,numrows
,numcols
);
1659 /**************************************************************************
1661 * Interpolates grib grid data to a point location.
1665 * gribinfo - pointer to a previously allocated gribinfo structure. The
1666 * gribinfo structure is filled in this function.
1667 * index - the index of gribinfo for which data is to be retrieved.
1668 * the first grib record is number 1.
1670 * pointdata- array of pointdata structures. Only the column and
1671 * row values in the structures need to be filled. On
1672 * output, the 'value' member of pointdata is filled.
1674 * numpoints- number of pointdata structures in the array.
1677 * on success - the interpolated value at the column,row location.
1678 * on failure - -99999
1680 ***************************************************************************/
1681 int rg_get_points(GribInfo
*gribinfo
, int index
, PointData pointdata
[],
1687 int numrows
, numcols
;
1688 int top
, left
, right
, bottom
;
1692 numrows
= rg_get_numrows(gribinfo
, index
);
1693 numcols
= rg_get_numcols(gribinfo
, index
);
1695 grib_out
= (float **)alloc_float_2d(numrows
,numcols
);
1696 if (grib_out
== NULL
) {
1697 fprintf(stderr
,"rg_get_points: Could not allocate space for grib_out\n");
1701 status
= rg_get_data(gribinfo
, index
, grib_out
);
1703 fprintf(stderr
,"rg_get_points: rg_get_data failed\n");
1707 for (idx
= 0; idx
< numpoints
; idx
++) {
1709 /* Change from 1 based to 0 based col/row */
1710 row
= pointdata
[idx
].row
;
1711 column
= pointdata
[idx
].column
;
1713 /* Do the interpolation here */
1714 bottom
= floor(row
);
1716 left
= floor(column
);
1717 right
= floor(column
+1);
1719 y1
= (row
- bottom
) * (grib_out
[top
][left
] - grib_out
[bottom
][left
]) +
1720 grib_out
[bottom
][left
];
1721 y2
= (row
- bottom
) * (grib_out
[top
][right
] - grib_out
[bottom
][right
]) +
1722 grib_out
[bottom
][right
];
1723 pointdata
[idx
].value
= (y2
- y1
) * (column
- left
) + y1
;
1727 free_float_2d(grib_out
,numrows
,numcols
);
1732 /**************************************************************************
1734 * Remove an element from an array and decrease, by one, indices of all
1735 * elements with an index greater than the index of the element to remove.
1739 * array - the integer array to manipulate
1740 * index - the index of the element to remove
1741 * size - the number of elements in the array
1743 ***************************************************************************/
1744 void remove_element(int array
[],int index
, int size
)
1748 for (j
= index
; j
< size
-1; j
++) {
1749 array
[j
] = array
[j
+1];
1754 /****************************************************************************
1755 * Advance the time by the input amount
1759 * century - an integer value for the century (20 for 1900's)
1760 * If the century is advanced, this value is advanced
1761 * and output to the calling routine.
1762 * year - a 2 digit value for the year.
1763 * month - a 2 digit value for the month.
1764 * day - the day of the month
1765 * hour - the hour of the day
1766 * amount - the amount to advance the time by.
1767 * unit - the units for the amount. These are values from table 4
1768 * of the grib manual.
1770 * a date in the form yymmddhh
1771 ****************************************************************************/
1773 int advance_time(int *century
, int year
, int month
, int day
, int hour
,
1774 int amount
, int unit
)
1776 int daysinmonth
[] = {31,28,31,30,31,30,31,31,30,31,30,31};
1781 hour
+= (int)((amount
/60.)+0.5);
1814 hour
+= (int)((amount
/12.)+0.5);
1816 hour
+= (int)((amount
/(60.*60.))+0.5);
1819 fprintf(stderr
,"WARNING: Could not advance time, incorrect unit: %d\n",
1824 while (hour
>= 24) {
1828 while (month
> 12) {
1833 /* if it is a leap year, change days in month for Feb now. */
1834 if (isLeapYear(year
)) daysinmonth
[1] = 29;
1836 while (day
> daysinmonth
[month
-1]) {
1837 day
-= daysinmonth
[month
-1];
1842 if (isLeapYear(year
))
1843 daysinmonth
[1] = 29;
1845 daysinmonth
[1] = 28;
1857 date
= hour
*1 + day
*100 + month
*10000 + year
*1000000;
1862 /****************************************************************************
1863 * Advance the time by the input amount
1867 * startdate - initialization date in the form yyyymmdd[HHMMSS]. If any
1868 * part of HHMMSS is not specified, it will be set to 0.
1869 * amount - the amount (in seconds) to advance the time by.
1872 * enddate[] - the time advanced to: yyyymmddHHMMSS format.
1878 ****************************************************************************/
1879 char *advance_time_str(char startdatein
[], int amount
, char enddate
[])
1886 strcpy(startdate
,startdatein
);
1887 while (strlen(startdate
) < 14) {
1888 strcpy(startdate
+(strlen(startdate
)),"0");
1891 /* This forces all calculations to use GMT time */
1895 sscanf(startdate
,"%4d%2d%2d%2d%2d%2d",&(starttp
.tm_year
),&(starttp
.tm_mon
),
1896 &(starttp
.tm_mday
),&(starttp
.tm_hour
),&(starttp
.tm_min
),
1898 starttp
.tm_mon
-= 1;
1899 starttp
.tm_year
-= 1900;
1900 time
= mktime(&starttp
);
1903 localtime_s(&endtp
, &time
);
1905 localtime_r(&time
, &endtp
);
1907 strftime(enddate
,15,"%Y%m%d%H%M%S",&endtp
);
1912 /****************************************************************************
1913 * Returns the difference in time in hours between date1 and date2
1918 * date1,date2: dates in yymmddhh format (integers)
1919 * century1,century2: centuries for each date (20 for 1900's).
1921 * the difference in time between the first and second dates in hours.
1922 ****************************************************************************/
1923 int date_diff(int date1
,int century1
,int date2
,int century2
)
1925 return (hours_since_1900(date1
,century1
) -
1926 hours_since_1900(date2
,century2
));
1929 /****************************************************************************
1930 * Returns the number of hours since Jan 1, at 00:00 1900.
1934 * date: integer in form yymmddhh
1935 * century: 2 digit century (20 for 1900's)
1937 * the number of hours since 00:00 Jan1, 1900.
1939 ****************************************************************************/
1940 int hours_since_1900(int date
,int century
)
1942 int daysinmonth
[] = {31,28,31,30,31,30,31,31,30,31,30,31};
1943 int hour
,day
,month
,year
;
1944 int days_since_1900
= 0;
1948 day
= (date
%10000)/100;
1949 month
= (date
%1000000)/10000;
1950 year
= (date
%100000000)/1000000;
1952 days_since_1900
+= day
;
1954 if (isLeapYear((century
-1)*100 + year
))
1955 daysinmonth
[1] = 29;
1957 daysinmonth
[1] = 28;
1959 for (i
= 0; i
< (month
- 1); i
++)
1960 days_since_1900
+= daysinmonth
[i
];
1962 for (i
=0; i
< (year
+ ((century
- 20)*100) - 1); i
++) {
1963 if (isLeapYear((century
- 1)*100 + year
))
1964 days_since_1900
+= 366;
1966 days_since_1900
+= 365;
1969 return days_since_1900
*24 + hour
;
1973 /****************************************************************************
1975 * Returns true if the input year is a leap year, otherwise returns false
1977 ****************************************************************************/
1978 int isLeapYear(int year
)
1980 if ( (((year
% 4) == 0) && ((year
% 100) != 0))
1981 || ((year
% 400) == 0) )
1988 /*****************************************************************************
1990 * Returns the number of grib elements (gribinfo->num_elements) processsed
1992 * gribinfo - pointer to a previously allocated gribinfo structure. The
1993 * gribinfo structure is filled in this function.
1996 * the number of elements in the gribinfo structure
1997 ****************************************************************************/
1999 int rg_num_elements(GribInfo
*gribinfo
){
2001 return gribinfo
->num_elements
;
2005 /*****************************************************************************
2007 * Deallocates the elements in the gribinfo structure and closes the files.
2010 * gribinfo - pointer to a previously allocated gribinfo structure. The
2011 * gribinfo structure is filled in this function.
2013 *****************************************************************************/
2014 void rg_free_gribinfo_elements(GribInfo
*gribinfo
)
2018 for (i
=0; i
<gribinfo
->num_elements
; i
++) {
2019 free(gribinfo
->elements
[i
].pds
);
2020 free(gribinfo
->elements
[i
].gds
);
2021 free(gribinfo
->elements
[i
].bms
);
2022 free(gribinfo
->elements
[i
].bds_head
);
2023 fclose(gribinfo
->elements
[i
].fp
);
2025 free(gribinfo
->elements
);
2028 /*****************************************************************************
2030 * Returns the value for level1 (gribinfo->usHeight1)
2032 * gribinfo - pointer to a previously allocated gribinfo structure. The
2033 * gribinfo structure is filled in this function.
2037 ****************************************************************************/
2039 int rg_get_level1(GribInfo
*gribinfo
, int index
)
2041 return gribinfo
->elements
[index
].usHeight1
;
2044 /*****************************************************************************
2046 * Returns the value for level2 (gribinfo->usHeight2)
2048 * gribinfo - pointer to a previously allocated gribinfo structure. The
2049 * gribinfo structure is filled in this function.
2053 ****************************************************************************/
2055 int rg_get_level2(GribInfo
*gribinfo
, int index
)
2057 return gribinfo
->elements
[index
].usHeight2
;
2060 /*****************************************************************************
2062 * returns number of rows in grid
2065 * gribinfo - pointer to a previously allocated gribinfo structure. The
2066 * gribinfo structure is filled in this function.
2069 * number of rows in grid
2070 *****************************************************************************/
2071 int rg_get_numrows(GribInfo
*gribinfo
,int index
)
2073 if ((gribinfo
->elements
[index
].gds
->head
.usData_type
== LATLON_PRJ
) ||
2074 (gribinfo
->elements
[index
].gds
->head
.usData_type
== GAUSS_PRJ
) ||
2075 (gribinfo
->elements
[index
].gds
->head
.usData_type
== ROT_LATLON_PRJ
) ||
2076 (gribinfo
->elements
[index
].gds
->head
.usData_type
== ROT_GAUSS_PRJ
) ||
2077 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_LATLON_PRJ
) ||
2078 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_GAUSS_PRJ
) ||
2079 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_ROT_LATLON_PRJ
)||
2080 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_ROT_GAUSS_PRJ
))
2082 return gribinfo
->elements
[index
].gds
->llg
.usNj
;
2083 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== MERC_PRJ
) {
2084 return gribinfo
->elements
[index
].gds
->merc
.rows
;
2085 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== LAMB_PRJ
) {
2086 return gribinfo
->elements
[index
].gds
->lam
.iNy
;
2087 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== POLAR_PRJ
) {
2088 return gribinfo
->elements
[index
].gds
->pol
.usNy
;
2092 /*****************************************************************************
2094 * returns number of columns in grid
2097 * gribinfo - pointer to a previously allocated gribinfo structure. The
2098 * gribinfo structure is filled in this function.
2101 * number of columns in grid
2102 *****************************************************************************/
2103 int rg_get_numcols(GribInfo
*gribinfo
,int index
)
2105 if ((gribinfo
->elements
[index
].gds
->head
.usData_type
== LATLON_PRJ
) ||
2106 (gribinfo
->elements
[index
].gds
->head
.usData_type
== GAUSS_PRJ
) ||
2107 (gribinfo
->elements
[index
].gds
->head
.usData_type
== ROT_LATLON_PRJ
) ||
2108 (gribinfo
->elements
[index
].gds
->head
.usData_type
== ROT_GAUSS_PRJ
) ||
2109 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_LATLON_PRJ
) ||
2110 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_GAUSS_PRJ
) ||
2111 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_ROT_LATLON_PRJ
)||
2112 (gribinfo
->elements
[index
].gds
->head
.usData_type
== STR_ROT_GAUSS_PRJ
))
2114 return gribinfo
->elements
[index
].gds
->llg
.usNi
;
2115 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== MERC_PRJ
) {
2116 return gribinfo
->elements
[index
].gds
->merc
.cols
;
2117 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== LAMB_PRJ
) {
2118 return gribinfo
->elements
[index
].gds
->lam
.iNx
;
2119 } else if (gribinfo
->elements
[index
].gds
->head
.usData_type
== POLAR_PRJ
) {
2120 return gribinfo
->elements
[index
].gds
->pol
.usNx
;
2124 /*****************************************************************************
2126 * returns the offset (in bytes) from the beginning of the file.
2129 * gribinfo - pointer to a filled gribinfo structure.
2132 * offset (in bytes) from beginning of file
2133 *****************************************************************************/
2134 int rg_get_offset(GribInfo
*gribinfo
,int index
)
2136 return gribinfo
->elements
[index
].offset
;
2138 /*****************************************************************************
2140 * returns the grib record ending position (in bytes) from the beginning of
2144 * gribinfo - pointer to a filled gribinfo structure.
2147 * position (in bytes) of the end of the grib record within the file.
2148 *****************************************************************************/
2149 int rg_get_end(GribInfo
*gribinfo
,int index
)
2151 return gribinfo
->elements
[index
].end
;
2153 /*****************************************************************************
2155 * returns grib id of input grid
2158 * gribinfo - pointer to a previously allocated gribinfo structure. The
2159 * gribinfo structure is filled in this function.
2162 * grib id of input grid
2163 *****************************************************************************/
2164 int rg_get_gridnum(GribInfo
*gribinfo
,int index
)
2166 return gribinfo
->elements
[index
].pds
->usGrid_id
;
2168 /*****************************************************************************
2173 * gribinfo - pointer to a previously allocated gribinfo structure. The
2174 * gribinfo structure is filled in this function.
2177 * date (yymmddhh) in integer type
2178 *****************************************************************************/
2179 int rg_get_date(GribInfo
*gribinfo
,int index
)
2181 return gribinfo
->elements
[index
].date
;
2183 /*****************************************************************************
2188 * gribinfo - pointer to a previously allocated gribinfo structure. The
2189 * gribinfo structure is filled in this function.
2192 * century in integer type
2193 *****************************************************************************/
2194 int rg_get_century(GribInfo
*gribinfo
,int index
)
2196 return gribinfo
->elements
[index
].century
;
2198 /*****************************************************************************
2200 * returns forecast time
2203 * gribinfo - pointer to a previously allocated gribinfo structure. The
2204 * gribinfo structure is filled in this function.
2207 * forecast time in units described by usFcst_unit_id
2208 *****************************************************************************/
2209 int rg_get_forecast_time(GribInfo
*gribinfo
,int index
)
2211 return gribinfo
->elements
[index
].pds
->usP1
;
2214 /*****************************************************************************
2216 * reads the gribmap file, and stores the information in the GribParameters
2220 * gribmap - pointer to a previously allocated GribParameters structure.
2221 * The gribmap structure is filled in this function.
2222 * file - the name of the gribmap file to read.
2225 * 1 - successful call to setup_gribinfo
2226 * -1 - call to setup_gribinfo failed
2228 *****************************************************************************/
2229 int rg_setup_gribmap(GribParameters
*gribmap
, char filename
[])
2233 int id
, center
, subcenter
, table
;
2236 fid
= fopen(filename
,"r");
2239 fprintf(stderr
,"Could not open %s\n",filename
);
2243 gribmap
->parms
= (GribTableEntry
*)malloc(sizeof(GribTableEntry
));
2246 while (fgets(line
,500,fid
) != NULL
)
2248 /* Skip over comments at begining of gribmap file */
2249 if (line
[0] == '#') continue;
2251 sscanf(line
,"%d:",&id
);
2254 sscanf(line
,"%d:%d:%d:%d",&id
,¢er
,&subcenter
,&table
);
2259 (GribTableEntry
*)realloc(gribmap
->parms
,
2260 (idx
+1)*sizeof(GribTableEntry
));
2261 gribmap
->parms
[idx
].center
= center
;
2262 gribmap
->parms
[idx
].subcenter
= subcenter
;
2263 gribmap
->parms
[idx
].table
= table
;
2264 sscanf(line
,"%d:%[^:]:%[^:]",&(gribmap
->parms
[idx
].parmid
),
2265 gribmap
->parms
[idx
].name
,
2266 gribmap
->parms
[idx
].comment
);
2272 gribmap
->num_entries
= idx
;
2278 /*****************************************************************************
2280 * finds the gribmap entry described by the gribmap file, and stores the information in the GribParameters
2284 * gribmap - pointer to structure that was filled by a call to
2286 * table - if set to -1, the first table the valid name will be used.
2287 * Otherwise, the table id must match as well.
2288 * name - name of the parameter to find.
2290 * gribmap_parms - pointer to GribTableEntry structure containing
2291 * information about the parameter that was found.
2294 * 1 - successful call to setup_gribinfo
2295 * -1 - call to setup_gribinfo failed
2297 *****************************************************************************/
2298 int rg_gribmap_parameter(GribParameters
*gribmap
, char name
[], int table
,
2299 GribTableEntry
*gribmap_parms
)
2305 for (idx
= 0; idx
< gribmap
->num_entries
; idx
++)
2308 if (strcmp(gribmap
->parms
[idx
].name
,name
) == 0)
2310 if ((table
== -1) || (table
== gribmap
->parms
[idx
].table
))
2312 /* We found a match! */
2321 memcpy(gribmap_parms
,&(gribmap
->parms
[idx
]),sizeof(GribTableEntry
));
2330 /*****************************************************************************
2332 * Deallocates the elements in the gribmap structure.
2335 * gribmap - pointer to a previously allocated gribmap structure. The
2336 * gribmap structure is filled in this function.
2338 *****************************************************************************/
2339 void rg_free_gribmap_elements(GribParameters
*gribmap
)
2341 free(gribmap
->parms
);
2344 /*****************************************************************************
2346 * Compares the elements in a findgrib structure with the elements in the
2347 * gribinfo structure for the input index. If they match, returns 1,
2348 * otherwise, returns 0.
2353 * index - the index of the grib record in gribinfo to compare to.
2355 *****************************************************************************/
2356 int compare_record(GribInfo
*gribinfo
, FindGrib
*findgrib
, int gribnum
)
2360 * Note (6/20/05): This searching is very inefficient. We may need to
2361 * improve this, since, for WRF, when searching through boundary data,
2362 * each search is slower that the previous, since the record to be
2363 * found turns out to be farther into the list.
2368 if ((findgrib
->center_id
== -INT_MAX
) ||
2369 findgrib
->center_id
== gribinfo
->elements
[gribnum
].center_id
) {
2370 if ((findgrib
->subcenter_id
== -INT_MAX
) ||
2371 findgrib
->subcenter_id
== gribinfo
->elements
[gribnum
].subcenter_id
) {
2372 if ((findgrib
->parmtbl_version
== -INT_MAX
) ||
2373 findgrib
->parmtbl_version
== gribinfo
->elements
[gribnum
].parmtbl
) {
2374 if ((strcmp(findgrib
->initdate
,"*") == 0) ||
2375 (strncmp(gribinfo
->elements
[gribnum
].initdate
,findgrib
->initdate
,
2376 strlen(findgrib
->initdate
)) == 0)) {
2377 if ((strcmp(findgrib
->validdate
,"*") == 0) ||
2378 (strncmp(gribinfo
->elements
[gribnum
].valid_time
,
2379 findgrib
->validdate
,
2380 strlen(findgrib
->validdate
)) == 0)) {
2381 if ((findgrib
->parmid
== -INT_MAX
) ||
2382 (findgrib
->parmid
==
2383 gribinfo
->elements
[gribnum
].usParm_id
)) {
2384 if ((findgrib
->leveltype
== -INT_MAX
) ||
2385 (findgrib
->leveltype
==
2386 gribinfo
->elements
[gribnum
].usLevel_id
)) {
2387 if ((findgrib
->level1
== -INT_MAX
) ||
2388 (findgrib
->level1
==
2389 gribinfo
->elements
[gribnum
].usHeight1
)) {
2390 if ((findgrib
->level2
== -INT_MAX
) ||
2391 (findgrib
->level2
==
2392 gribinfo
->elements
[gribnum
].usHeight2
)) {
2393 if ((findgrib
->fcsttime1
== -INT_MAX
) ||
2394 (findgrib
->fcsttime1
==
2395 gribinfo
->elements
[gribnum
].fcsttime1
)) {
2396 if ((findgrib
->fcsttime2
== -INT_MAX
) ||
2397 (findgrib
->fcsttime2
==
2398 gribinfo
->elements
[gribnum
].fcsttime2
)) {
2416 /*****************************************************************************
2418 * returns the multiplication factor to convert grib forecast times to
2422 * unit_id - grib forecast unit id, from Table 4.
2426 *****************************************************************************/
2427 int get_factor2(int unit
)
2451 /* This is a WSI (non-standard) time unit of 5 minutes */
2458 fprintf(stderr
,"Invalid unit for forecast time: %d\n",unit
);