/* **************************************************************************** * * Routines for indexing, reading and writing grib files. Routines * are designed to be called by Fortran. * * All routines return 0 for success, 1 for failure, unless otherwise noted. * * Todd Hutchinson * WSI * 05/17/2005 * **************************************************************************** */ #include "grib1_routines.h" #include "gridnav.h" #include "wrf_projection.h" #include #include #include #include char *trim (char *str); int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid); int index_times(GribInfo *gribinfo, Times *times); int find_time(Times *times, char valid_time[15]); int get_gridnav_projection(int wrf_projection); int get_byte(int input_int, int bytenum); /* * Allocate space for the fileindex structure */ int ALLOC_INDEX_FILE(FileIndex *fileindex) { int status = 0; fileindex->gribinfo = (GribInfo *)malloc(sizeof(GribInfo)); if (fileindex->gribinfo == NULL) { fprintf(stderr,"Allocating fileindex->gribinfo failed.\n"); status = 1; return status; } fileindex->metadata = (MetaData *)malloc(sizeof(MetaData)); if (fileindex->metadata == NULL) { fprintf(stderr,"Allocating fileindex->metadata failed.\n"); status = 1; return status; } fileindex->metadata->elements = NULL; fileindex->times = (Times *)malloc(sizeof(Times)); if (fileindex->times == NULL) { fprintf(stderr,"Allocating fileindex->times failed.\n"); status = 1; return status; } fileindex->times->elements = NULL; return status; } void FREE_INDEX_FILE(FileIndex *fileindex) { int status = 0; rg_free_gribinfo_elements(fileindex->gribinfo); free(fileindex->gribinfo); free(fileindex->metadata->elements); free(fileindex->metadata); free(fileindex->times->elements); free(fileindex->times); } int INDEX_FILE(int *fid, FileIndex *fileindex) { int status; /* Index the grib records */ status = rg_setup_gribinfo_i(fileindex->gribinfo,*fid,1); if (status < 0) { fprintf(stderr,"Error setting up gribinfo structure.\n"); return 1; } /* Index the metadata section */ status = index_metadata(fileindex->gribinfo, fileindex->metadata, *fid); if (status != 0) { fprintf(stderr,"Error setting up metadata structure.\n"); return 1; } /* Setup a list of times based on times in grib records */ status = index_times(fileindex->gribinfo, fileindex->times); if (status != 0) { fprintf(stderr,"Error indexing times in grib file.\n"); return 1; } return 0; } int GET_FILEINDEX_SIZE(int *size) { *size = sizeof(FileIndex); return *size; } int GET_NUM_TIMES(FileIndex *fileindex, int *numtimes) { *numtimes = (fileindex->times)->num_elements; return *numtimes; } int GET_TIME(FileIndex *fileindex, int *idx, char time[]) { int num_times; int year, month, day, minute, hour, second; char time2[100]; num_times = GET_NUM_TIMES(fileindex,&num_times); if (*idx > num_times) { fprintf(stderr,"Tried to get time %d, but only %d times exist\n", *idx, num_times); return 1; } strcpy(time,fileindex->times->elements[*idx-1].valid_time); /* Reformat time to meet WRF time format */ sscanf(time, "%4d%2d%2d%2d%2d%2d", &year, &month, &day, &hour, &minute, &second); sprintf(time2, "%04d-%02d-%02d_%02d:%02d:%02d", year, month, day, hour, minute, second); strncpy(time,time2,19); return 0; } int GET_LEVEL1(FileIndex *fileindex, int *idx, int *level1) { *level1 = (fileindex->gribinfo)->elements[*idx].usHeight1; return *level1; } int GET_LEVEL2(FileIndex *fileindex, int *idx, int *level2) { *level2 = (fileindex->gribinfo)->elements[*idx].usHeight2; return *level2; } int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid) { int status=0; int end; char string[11]; int found_metadata=0; int idx; int pos; int fileend; int seekpos; int bytesread; char line[1000]; char element[100],datestr[100],varname[100]; char value[1000]; int incomment; int charidx; int elemidx=0; FILE *stream; /* Associate a FILE *stream with the file id */ stream = fdopen(fid,"r"); if (stream == NULL) { perror("Error associating stream with file descriptor"); status = -1; return status; } /* * First, set the position to end of grib data (the * metadata section comes after the grib data). */ idx = rg_num_elements(gribinfo) - 1; end = rg_get_end(gribinfo,idx); pos = end + 1; fileend = lseek(fid,0,SEEK_END); /* * Now, start searching for metadata */ while (pos < (fileend - 10)) { seekpos = pos; pos = lseek(fid,seekpos,SEEK_SET); if (pos != seekpos) { fprintf(stderr,"Error seeking %d bytes in file\n",end); perror(""); return 1; } bytesread = read(fid,string,10); if (bytesread != 10) { fprintf(stderr,"Invalid read, pos: %d :\n",pos); perror(""); pos += 1; continue; } if (strncmp(string,"",10) == 0) { /* We found it, so break out ! */ found_metadata = 1; break; } pos += 1; } /* Now, read metadata, line by line */ incomment = 0; while(fgets(line,1000,stream) != NULL) { trim(line); /* Set comment flag, if we found a comment */ if (strncmp(line,"",3) == 0) { strcpy(line,line+charidx+3); incomment = 0; break; } else { charidx++; } } } if (incomment) continue; /* Check for end of metadata */ if (strncmp(line,"",11) == 0) { /* We found end of data, so, break out */ break; } /* Skip blank lines */ if (strlen(line) == 0) continue; /* Parse line */ trim(line); strcpy(element,"none"); strcpy(datestr,"none"); strcpy(varname,"none"); strcpy(value,"none"); if (sscanf(line,"%[^;=];%[^;=];%[^;=]=%[^\n]",varname,datestr, element,value) == 4) { } else if (sscanf(line,"%[^;=];%[^;=]=%[^\n]",datestr,element,value) == 3) { strcpy(varname,"none"); } else if (sscanf(line,"%[^;=]=%[^\n]",element,value) == 2) { strcpy(varname,"none"); strcpy(datestr,"none"); } else { strcpy(varname,"none"); strcpy(datestr,"none"); strcpy(element,"none"); strcpy(value,"none"); fprintf(stderr,"Invalid line in metadata: \n%s",line); } trim(varname); trim(datestr); trim(element); trim(value); metadata->elements = (MetaData_Elements *)realloc( metadata->elements, (elemidx+1)*sizeof(MetaData_Elements) ); strcpy(metadata->elements[elemidx].VarName,varname); strcpy(metadata->elements[elemidx].DateStr,datestr); strcpy(metadata->elements[elemidx].Element,element); strcpy(metadata->elements[elemidx].Value,value); elemidx++; } metadata->num_elements = elemidx; return 0; } int index_times(GribInfo *gribinfo, Times *times) { int idx; int status; int numtimes=0; int date; char valid_time[15]; char tmp[15]; int swapped; times->num_elements = 0; /* Loop through elements, and build list of times */ for (idx=0; idx < gribinfo->num_elements; idx++) { /* Calculate valid time */ status = rg_get_valid_time(gribinfo,idx,valid_time); if (status != 0) { fprintf(stderr,"Could not retrieve valid time for index: %d\n",idx); continue; } /* * Check if this time is already contained in times * If not, allocate space for it, and add it to list */ if (find_time(times,valid_time) < 0) { times->num_elements++; times->elements = (Times_Elements *) realloc(times->elements,times->num_elements*sizeof(Times_Elements)); if (times->elements == NULL) { fprintf(stderr,"Allocating times->elements failed.\n"); status = 1; return status; } strcpy(times->elements[times->num_elements - 1].valid_time,valid_time); } } /* Sort times */ swapped = 1; while (swapped) { swapped=0; for (idx=1; idx < times->num_elements; idx++) { if (strcmp(times->elements[idx-1].valid_time, times->elements[idx].valid_time) > 0) { strcpy(tmp,times->elements[idx-1].valid_time); strcpy(times->elements[idx-1].valid_time, times->elements[idx].valid_time); strcpy(times->elements[idx].valid_time, tmp); swapped = 1; } } } return 0; } int find_time(Times *times, char valid_time[15]) { int idx; int found_elem = -1; for (idx = 0; idx < times->num_elements; idx++) { if (strcmp(times->elements[idx].valid_time,valid_time) == 0) { found_elem = idx; break; } } return found_elem; } int GET_METADATA_VALUE(FileIndex *fileindex, char ElementIn[], char DateStrIn[], char VarNameIn[], char Value[], int *stat, int strlen1, int strlen2, int strlen3, int strlen4, int strlen5) { int elemidx; int elemnum; char VarName[200]; char DateStr[200]; char Element[200]; int Value_Len; *stat = 0; strncpy(Element,ElementIn,strlen2); Element[strlen2] = '\0'; strncpy(DateStr,DateStrIn,strlen3); DateStr[strlen3] = '\0'; strncpy(VarName,VarNameIn,strlen4); VarName[strlen4] = '\0'; Value_Len = strlen5; elemnum = -1; for (elemidx = 0; elemidx < fileindex->metadata->num_elements; elemidx++) { if (strcmp(Element,fileindex->metadata->elements[elemidx].Element) == 0) { if (strcmp(DateStr,fileindex->metadata->elements[elemidx].DateStr) == 0) { if (strcmp(VarName, fileindex->metadata->elements[elemidx].VarName) == 0) { elemnum = elemidx; break; } } } } if (elemnum != -1) { strncpy(Value, fileindex->metadata->elements[elemnum].Value, Value_Len); } else { strncpy(Value, "none", Value_Len); *stat = 1; } /* * Pad end of string with one space. This allows Fortran internal * read function to work properly. */ if (strlen(Value) < Value_Len) { strcpy(Value + strlen(Value), " "); } return elemidx; } int GET_GRIB_INDEX(FileIndex *fileindex, int *center, int *subcenter, int *parmtbl, int *parmid, char DateStrIn[], int *leveltype, int *level1, int *level2, int *fcsttime1, int *fcsttime2, int *index, int strlen1, int strlen2) { char DateStr[1000]; FindGrib findgrib; strncpy(DateStr,DateStrIn,strlen2); DateStr[strlen2] = '\0'; grib_time_format(DateStr,DateStr); rg_init_findgrib(&findgrib); strncpy(findgrib.initdate,DateStrIn,strlen2); findgrib.initdate[strlen2] = '\0'; findgrib.parmid = *parmid; findgrib.leveltype = *leveltype; findgrib.level1 = *level1; findgrib.level2 = *level2; findgrib.fcsttime1 = *fcsttime1; findgrib.fcsttime2 = *fcsttime2; findgrib.center_id = *center; findgrib.subcenter_id = *subcenter; findgrib.parmtbl_version = *parmtbl; *index = rg_get_index(fileindex->gribinfo, &findgrib); return *index; } int GET_GRIB_INDEX_GUESS(FileIndex *fileindex, int *center, int *subcenter, int *parmtbl, int *parmid, char DateStrIn[], int *leveltype, int *level1, int *level2, int *fcsttime1,int *fcsttime2, int *guessidx, int *index, int strlen1, int strlen2) { char DateStr[1000]; FindGrib findgrib; strncpy(DateStr,DateStrIn,strlen2); DateStr[strlen2] = '\0'; grib_time_format(DateStr,DateStr); rg_init_findgrib(&findgrib); strncpy(findgrib.initdate,DateStrIn,strlen2); findgrib.initdate[strlen2] = '\0'; findgrib.parmid = *parmid; findgrib.leveltype = *leveltype; findgrib.level1 = *level1; findgrib.level2 = *level2; findgrib.fcsttime1 = *fcsttime1; findgrib.fcsttime2 = *fcsttime2; findgrib.center_id = *center; findgrib.subcenter_id = *subcenter; findgrib.parmtbl_version = *parmtbl; *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx); return *index; } int GET_GRIB_CENTER(FileIndex *fileindex, int *parmid, int *center) { *center = rg_get_center_id(fileindex->gribinfo,*parmid); return *center; } int GET_GRIB_SUBCENTER(FileIndex *fileindex, int *parmid, int *subcenter) { *subcenter = rg_get_subcenter_id(fileindex->gribinfo,*parmid); return *subcenter; } int GET_GRIB_TBLVERSION(FileIndex *fileindex, int *parmid, int *parmtbl) { *parmtbl = rg_get_parmtbl(fileindex->gribinfo,*parmid); return *parmtbl; } int GET_GRIB_PROCID(FileIndex *fileindex, int *parmid, int *proc_id) { *proc_id = rg_get_proc_id(fileindex->gribinfo,*parmid); return *proc_id; } int GET_GRIB_INDEX_VALIDTIME(FileIndex *fileindex, int *center, int *subcenter, int *parmtbl, int *parmid, char DateStrIn[], int *leveltype, int *level1, int *level2, int *index, int strlen1, int strlen2) { char DateStr[1000]; FindGrib findgrib; strncpy(DateStr,DateStrIn,strlen2); DateStr[strlen2] = '\0'; grib_time_format(DateStr,DateStr); rg_init_findgrib(&findgrib); strncpy(findgrib.validdate,DateStr,strlen2); findgrib.initdate[strlen2] = '\0'; findgrib.parmid = *parmid; findgrib.leveltype = *leveltype; findgrib.level1 = *level1; findgrib.level2 = *level2; findgrib.center_id = *center; findgrib.subcenter_id = *subcenter; findgrib.parmtbl_version = *parmtbl; *index = rg_get_index(fileindex->gribinfo, &findgrib); return *index; } int GET_GRIB_INDEX_VALIDTIME_GUESS(FileIndex *fileindex, int *center, int *subcenter, int *parmtbl, int *parmid, char DateStrIn[], int *leveltype, int *level1, int *level2, int *guessidx, int *index, int strlen1, int strlen2) { char DateStr[1000]; FindGrib findgrib; strncpy(DateStr,DateStrIn,strlen2); DateStr[strlen2] = '\0'; grib_time_format(DateStr,DateStr); rg_init_findgrib(&findgrib); strncpy(findgrib.validdate,DateStr,strlen2); findgrib.initdate[strlen2] = '\0'; findgrib.parmid = *parmid; findgrib.leveltype = *leveltype; findgrib.level1 = *level1; findgrib.level2 = *level2; findgrib.center_id = *center; findgrib.subcenter_id = *subcenter; findgrib.parmtbl_version = *parmtbl; *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx); return *index; } int GET_GRIB_INDICES(FileIndex *fileindex, int *center, int *subcenter, int *parmtbl,int *parmid, char DateStrIn[], int *leveltype, int *level1, int *level2, int *fcsttime1, int *fcsttime2, int *indices, int *num_indices, int strlen1, int strlen2) { char DateStr[1000]; int status; FindGrib findgrib; strncpy(DateStr,DateStrIn,strlen2); DateStr[strlen2] = '\0'; grib_time_format(DateStr,DateStr); rg_init_findgrib(&findgrib); strncpy(findgrib.initdate,DateStrIn,strlen2); findgrib.initdate[strlen2] = '\0'; trim(findgrib.initdate); findgrib.parmid = *parmid; findgrib.leveltype = *leveltype; findgrib.level1 = *level1; findgrib.level2 = *level2; findgrib.fcsttime1 = *fcsttime1; findgrib.fcsttime2 = *fcsttime2; findgrib.center_id = *center; findgrib.subcenter_id = *subcenter; findgrib.parmtbl_version = *parmtbl; *num_indices = rg_get_indices(fileindex->gribinfo, &findgrib, indices); return (*num_indices); } int GET_GRID_INFO_SIZE(int *size) { *size = sizeof(Grid_Info); return *size; } int LOAD_GRID_INFO(char *varnameIn, char *initdateIn, int *leveltype, int *level1, int *level2, float *fcst_time, int *accum_period, int *grid_id, int *projection, int *xpoints, int *ypoints, float *center_lat, float *center_lon, float *Di, float *Dj,float *central_lon, int *proj_center_flag, float *latin1, float *latin2, Grib1_Tables *grib_tables, Grid_Info *grid_info, int strlen1, int strlen2) { char varname[1000], initdate[1000]; strncpy(varname,varnameIn,strlen1); varname[strlen1] = '\0'; strncpy(initdate,initdateIn,strlen2); initdate[strlen2] = '\0'; strcpy(grid_info->varname, varname); strcpy(grid_info->initdate, initdate); grid_info->leveltype = *leveltype; grid_info->level1 = *level1 ; grid_info->level2 = *level2 ; grid_info->fcst_time = *fcst_time ; grid_info->accum_period = *accum_period ; grid_info->grid_id = *grid_id ; grid_info->projection = *projection ; grid_info->xpoints = *xpoints ; grid_info->ypoints = *ypoints ; grid_info->center_lat = *center_lat ; grid_info->center_lon = *center_lon; grid_info->Di = *Di ; grid_info->Dj = *Dj ; grid_info->central_lon = *central_lon ; grid_info->proj_center_flag = *proj_center_flag ; grid_info->latin1 = *latin1 ; grid_info->latin2 = *latin2 ; grid_info->grib_tables = copy_grib_tables(grib_tables); return 0; } int PRINT_GRID_INFO(Grid_Info *grid_info) { fprintf(stdout,"varname =%s\n",grid_info->varname); fprintf(stdout,"initdate =%s\n",grid_info->initdate); fprintf(stdout,"leveltype =%d\n",grid_info->leveltype); fprintf(stdout,"level1 =%d\n",grid_info->level1); fprintf(stdout,"level2 =%d\n",grid_info->level2); fprintf(stdout,"fcst_time =%f\n",grid_info->fcst_time); fprintf(stdout,"accum_period =%d\n",grid_info->accum_period); fprintf(stdout,"grid_id =%d\n",grid_info->grid_id); fprintf(stdout,"projection =%d\n",grid_info->projection); fprintf(stdout,"xpoints =%d\n",grid_info->xpoints); fprintf(stdout,"ypoints =%d\n",grid_info->ypoints); fprintf(stdout,"center_lat =%f\n",grid_info->center_lat); fprintf(stdout,"center_lon =%f\n",grid_info->center_lon); fprintf(stdout,"Di =%f\n",grid_info->Di); fprintf(stdout,"Dj =%f\n",grid_info->Dj); fprintf(stdout,"central_lon =%f\n",grid_info->central_lon); fprintf(stdout,"proj_center_flag =%d\n",grid_info->proj_center_flag); fprintf(stdout,"latin1 =%f\n",grid_info->latin1); fprintf(stdout,"latin2 =%f\n",grid_info->latin2); return 0; } int GET_SIZEOF_GRID(FileIndex *fileindex, int *index, int *numcols, int *numrows) { *numcols = rg_get_numcols(fileindex->gribinfo,*index); *numrows = rg_get_numrows(fileindex->gribinfo,*index); return (*numcols)*(*numrows); } void FREE_GRID_INFO(Grid_Info *grid_info) { FREE_GRIBMAP(grid_info->grib_tables); } int READ_GRIB(FileIndex *fileindex, int *fid, int *index, float *data) { int status; status = rg_get_data_1d(fileindex->gribinfo,*index,data); return status; } #define LINESIZE 300 #define SECS_IN_SEC 1 #define SECS_IN_MIN 60 #define MINS_IN_HOUR 60 #define MINS_IN_5MINS 5 #define HOURS_IN_DAY 24 #define MAX_FCST 65535 #define MAX_FCST_SECS MAX_FCST*SECS_IN_SEC #define MAX_FCST_MINS MAX_FCST*SECS_IN_MIN #define MAX_FCST_5MINS MAX_FCST*MINS_IN_5MINS*SECS_IN_MIN #define MAX_FCST_HOURS MAX_FCST*MINS_IN_HOUR*SECS_IN_MIN #define MAX_FCST_DAYS MAX_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN #define MAX1B_FCST 256 #define MAX1B_FCST_SECS MAX1B_FCST*SECS_IN_SEC #define MAX1B_FCST_MINS MAX1B_FCST*SECS_IN_MIN #define MAX1B_FCST_5MINS MAX1B_FCST*MINS_IN_5MINS*SECS_IN_MIN #define MAX1B_FCST_HOURS MAX1B_FCST*MINS_IN_HOUR*SECS_IN_MIN #define MAX1B_FCST_DAYS MAX1B_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN typedef struct { int time_range; int fcst_unit; int P1; int P2; int time_range_ext; int fcst_unit_ext_1; int fcst_unit_ext_2; int P1_ext; int P2_ext; } FcstTimeStruct; int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *fcst_time); /**************************************************************************** * * This function takes in metadata in the grid_info structure, output data in * the *data array, and calls routines to write the metadata and data * in grib version 1 format the open file descriptor filefd. * ****************************************************************************/ int WRITE_GRIB(Grid_Info *grid_info, int *filefd, float *data) { GRIB_HDR *gh=NULL; DATA_INPUT data_input; GEOM_IN geom_in; USER_INPUT user_input; int grid_projection; int status; float x_center, y_center; GridNav gridnav; float first_lat, first_lon, last_lat, last_lon; int year, month, day, hour, minute; float second; char varname2[1000]; int table_index; int fcst_unit; int time_range; int P1, P2; int fcst_unit_ext_1, fcst_unit_ext_2; int P1_ext, P2_ext; int time_range_ext; char errmsg[1000]; int center, subcenter, parmtbl; int tablenum; float dx, dy; FcstTimeStruct fcst_time; strcpy(varname2,grid_info->varname); trim(varname2); sscanf(grid_info->initdate,"%d-%d-%d_%d:%d:%f", &year,&month,&day,&hour,&minute,&second); /* Get coords of center of grid */ x_center = (grid_info->xpoints + 1)/2.; y_center = (grid_info->ypoints + 1)/2.; grid_projection = get_gridnav_projection(grid_info->projection); /* Convert grid spacing to degrees for LATLON projection */ if ( (grid_info->projection == WRF_LATLON) || (grid_info->projection == WRF_CASSINI) ) { dx = grid_info->Di * KM_TO_DEGREES; dy = grid_info->Dj * KM_TO_DEGREES; } else { dx = grid_info->Di; dy = grid_info->Dj; } /* Initialize grid structure */ status = GRID_init(grid_info->center_lat, grid_info->central_lon, grid_projection, grid_info->latin1, grid_info->latin2, grid_info->xpoints, grid_info->ypoints, dx, dy, grid_info->center_lat, grid_info->center_lon, x_center, y_center, &gridnav); if (!status) { fprintf(stderr,"write_grib: error from GRID_init\n"); } /* get lat/lon of lower left corner */ status = GRID_to_latlon(&gridnav, 1, 1, &first_lat, &first_lon); if (!status) { fprintf(stderr, "write_grib: error from GRID_to_latlon for first lat/lon\n"); } /* get lat/lon of upper right corner */ status = GRID_to_latlon(&gridnav, grid_info->xpoints, grid_info->ypoints, &last_lat, &last_lon); if (!status) { fprintf(stderr, "write_grib: error from GRID_to_latlon for last lat/lon\n"); } /* Read the grib parameter table */ status = GET_GRIB_PARAM(grid_info->grib_tables, varname2, ¢er, &subcenter, &parmtbl, &tablenum, &table_index, 1,strlen(varname2)); if (table_index < 0) { fprintf(stderr,\ "Skipping %s, Could not find parameter for %s in gribmap.txt\n",\ varname2,varname2); return 1; } /* * We skip any parameters that are listed in parameter 255 in gribmap.txt. * Parameter 255 is used to indicate that a WRF parameter should not be * output. It is useful for parameters that are requested to be output in * the WRF Registry, but are already implicitly output in grib. */ if (table_index == 255) { return 0; } /* * Setup the geom_in structure for the grib library. Here, we set * the generic parms. Below, we set the projection specific parms */ geom_in.nx = grid_info->xpoints; geom_in.ny = grid_info->ypoints; geom_in.first_lat = first_lat; geom_in.first_lon = first_lon; geom_in.last_lat = last_lat; geom_in.last_lon = last_lon; geom_in.scan = 64; switch (grid_info->projection) { case WRF_LATLON: case WRF_CASSINI: strcpy(geom_in.prjn_name,"spherical"); geom_in.parm_1 = dy; geom_in.parm_2 = dx; geom_in.parm_3 = -1; geom_in.usRes_flag = 0; /* * Set to 0 here, MEL grib library will reset * to 128 to indicate that direction * increments are given. */ break; case WRF_MERCATOR: strcpy(geom_in.prjn_name,"mercator"); geom_in.parm_1 = grid_info->latin1; geom_in.parm_2 = dx; geom_in.parm_3 = dy; geom_in.usRes_flag = 128; break; case WRF_LAMBERT: strcpy(geom_in.prjn_name,"lambert"); geom_in.usRes_flag = 0; /* Set to 0 here, MEL grib library will reset * to 128. */ geom_in.parm_3 = grid_info->central_lon; geom_in.x_int_dis = dx; geom_in.y_int_dis = dy; geom_in.parm_1 = grid_info->latin1; geom_in.parm_2 = grid_info->latin2; break; case WRF_POLAR_STEREO: strcpy(geom_in.prjn_name,"polar_stereo"); geom_in.usRes_flag = 0; /* Set to 0 here, MEL grib library will reset * to 128. */ geom_in.parm_3 = -1; geom_in.x_int_dis = dx*(1.+sin(60. * PI/180.)) / (1.+sin(abs(grid_info->latin1) * PI/180.)); geom_in.y_int_dis = dy*(1.+sin(60. * PI/180.)) / (1.+sin(abs(grid_info->latin1) * PI/180.)); geom_in.parm_1 = -1; geom_in.parm_2 = grid_info->central_lon; break; default: fprintf(stderr,"Error, invalid projection: %d\n",grid_info->projection); return 1; } /* * Setup the data_input structure. */ data_input.nDec_sc_fctr = grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index]; data_input.usProc_id = 220; data_input.usGrid_id = grid_info->grid_id; data_input.usParm_id = grid_info->grib_tables->grib_table_info[tablenum].parm_id[table_index]; data_input.usParm_sub_id = grid_info->grib_tables->grib_table_info[tablenum].subcenter; data_input.usLevel_id = grid_info->leveltype; if (grid_info->leveltype == 112) { data_input.nLvl_1 = grid_info->level2; data_input.nLvl_2 = grid_info->level1; } else { data_input.nLvl_1 = grid_info->level1; data_input.nLvl_2 = grid_info->level2; } data_input.nYear = year; data_input.nMonth = month; data_input.nDay = day; data_input.nHour = hour; data_input.nMinute = minute; data_input.nSecond = second; status = get_fcst_time(grid_info->accum_period, grid_info->fcst_time, &fcst_time); data_input.usFcst_id = fcst_time.fcst_unit; data_input.usFcst_per1 = fcst_time.P1; data_input.usFcst_per2 = fcst_time.P2; data_input.usTime_range_id = fcst_time.time_range; data_input.usTime_range_avg = 0; data_input.usTime_range_mis = 0; /* * This is for WSI's extended PDS section */ data_input.PDS_41 = fcst_time.fcst_unit_ext_1; data_input.PDS_42 = fcst_time.P1_ext; data_input.PDS_46 = fcst_time.fcst_unit_ext_2; data_input.PDS_47 = fcst_time.P2_ext; data_input.PDS_51 = fcst_time.time_range_ext; data_input.PDS_52 = 0; data_input.nDec_sc_fctr = grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index]; user_input.usCenter_id = grid_info->grib_tables->grib_table_info[tablenum].center; user_input.usParm_tbl = grid_info->grib_tables->grib_table_info[tablenum].parmtbl; user_input.chCase_id='0'; user_input.usSub_tbl = 0; user_input.usCenter_sub = grid_info->grib_tables->grib_table_info[tablenum].subcenter; user_input.usTrack_num = 0; user_input.usGds_bms_id = 128; user_input.usBDS_flag = 0; user_input.usBit_pack_num = 0; status = init_gribhdr(&gh,errmsg); if (status != 0) { fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg); return 1; } status = grib_enc(data_input,user_input,geom_in,data,gh,errmsg); if (status != 0) { fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg); fprintf (stderr,"\tCheck precision for %s in gribmap.txt.\n", varname2); return 1; } status = gribhdr2filed(gh,*filefd,errmsg); if (status != 0) { fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg); return 1; } free_gribhdr(&gh); return 0; } /*************************************************************************** * Function to set up a structure containing forecast time parameters * This encodes the standard grib forecast time parameters as well * as WSI's extended forecast time parameters. ***************************************************************************/ int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *ft) { /* * Added ability to output a "5-minute" forecast time unit for the * sake of WxProducer. This allows WxPro to ingest data beyond * 18 hours, and accumulation data beyond 255 units. */ /* * Initialize. */ ft->time_range = 0; ft->fcst_unit = 0; ft->P1 = 0; ft->P2 = 0; ft->fcst_unit_ext_1 = 0; ft->fcst_unit_ext_2 = 0; ft->P1_ext = 0; ft->P2_ext = 0; ft->time_range_ext = 0; if (accum_period == 0) { if (fcst_secs < MAX_FCST_SECS) { ft->time_range = 10; ft->fcst_unit = 254; ft->P1 = get_byte(fcst_secs,2); ft->P2 = get_byte(fcst_secs,1); } else if (((fcst_secs % SECS_IN_MIN) == 0) && (fcst_secs < MAX_FCST_MINS)) { ft->time_range = 10; ft->fcst_unit = 0; ft->P1 = get_byte(fcst_secs/SECS_IN_MIN,2); ft->P2 = get_byte(fcst_secs/SECS_IN_MIN,1); } else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR) == 0) && (fcst_secs < MAX_FCST_HOURS)) { ft->time_range = 10; ft->fcst_unit = 1; ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),2); ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),1); } /* * MAX_FCST_DAYS is causing an integer overflow, so, we'll just skip * the check here. It's very unlikely that someone would exceed this * anyway (5.6 million days!) */ /* else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0) && (fcst_secs < MAX_FCST_DAYS)) */ else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0)) { ft->time_range = 10; ft->fcst_unit = 2; ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),2); ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),1); } else if (((fcst_secs % SECS_IN_MIN*MINS_IN_5MINS) == 0) && (fcst_secs < MAX_FCST_5MINS)) { ft->time_range = 10; ft->fcst_unit = 50; ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),2); ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),1); } else { ft->time_range = 255; ft->fcst_unit = 0; ft->P1 = 0; ft->P2 = 0; ft->fcst_unit_ext_1 = 254; ft->fcst_unit_ext_2 = 254; ft->P1_ext = fcst_secs; ft->P2_ext = 0; ft->time_range_ext = 0; } } else /* Accumulation period is not 0 */ { if ((fcst_secs < MAX1B_FCST_HOURS) && (fcst_secs%(SECS_IN_MIN*MINS_IN_HOUR) == 0) && (accum_period%(SECS_IN_MIN*MINS_IN_HOUR) == 0)) { ft->time_range = 4; ft->fcst_unit = 1; ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_HOUR); ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR); } else if ((fcst_secs < MAX1B_FCST_MINS) && ((fcst_secs-accum_period)%SECS_IN_MIN == 0) && (fcst_secs%SECS_IN_MIN == 0)) { ft->time_range = 4; ft->fcst_unit = 0; ft->P1 = (fcst_secs-accum_period)/SECS_IN_MIN; ft->P2 = fcst_secs/SECS_IN_MIN; } else if (fcst_secs < MAX1B_FCST_SECS) { ft->time_range = 4; ft->fcst_unit = 254; ft->P1 = fcst_secs-accum_period; ft->P2 = fcst_secs; } else if ((fcst_secs < MAX1B_FCST_5MINS) && (fcst_secs%(SECS_IN_MIN*MINS_IN_5MINS) == 0) && (accum_period%(SECS_IN_MIN*MINS_IN_5MINS) == 0)) { ft->time_range = 4; ft->fcst_unit = 50; ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_5MINS); ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS); } else { ft->time_range = 255; ft->fcst_unit = 0; ft->P1 = 0; ft->P2 = 0; ft->fcst_unit_ext_1 = 254; ft->fcst_unit_ext_2 = 254; ft->P1_ext = fcst_secs - accum_period; ft->P2_ext = accum_period; ft->time_range_ext = 203; /* Duration */ } } return 0; } /****************************************************************************** * returns a byt from an input integer *****************************************************************************/ int get_byte(int input_int, int bytenum) { int out; out = ((input_int >> (bytenum-1)*8) & ~(~0 <<8)); return out; } /************************************************************************* * Converts from WRF time format to time format required by grib routines *************************************************************************/ int grib_time_format(char *DateStr, char *DateStrIn) { int year,month,day,hour,minute,second; trim(DateStrIn); if (DateStrIn[0] == '*') { strcpy(DateStr,"*"); } else { sscanf(DateStrIn,"%04d-%02d-%02d_%02d:%02d:%02d", &year,&month,&day,&hour,&minute,&second); sprintf(DateStr,"%04d%02d%02d%02d%02d%02d", year,month,day,hour,minute,second); } return 0; }