source: trunk/WRF.COMMON/WRFV3/external/io_grib1/grib1_util/read_grib.c @ 2759

 /**************************************************************************
 * Todd Hutchinson          4/20/98
 * tahutchinson@tasc.com    (781) 942-2000 x3108
 * TASC 
 * 55 Walkers Brook Drive
 * Reading, MA  01867
 *
 * Functions in this file are used for decoding grib data.  Please see the
 * headers before each function for a full descrption.
 *
 * Routines in this file call functions in the Naval Research Lab's grib 
 * library.  The grib library is freely available from 
 * http://www-mel.nrlmry.navy.mil/cgi-bin/order_grib.  This library should
 * be installed on your system prior to using the routines in this file.
 * Documentation for this library is available from
 *  Master Environmental Grib Library user's manual
 *  http://mel.dmso.mil/docs/grib.pdf
 * Note: the standard NRL grib library does not support 
 * "Little-Endian" platforms such as linux.  There is a version of the NRL 
 * grib library within the WxPredictor project which does support linux.
 *
 * This file references the cfortran.h header file to ease the use of calling 
 * this function from a fortran routine.  cfortran.h is a header file that 
 * allows for simple machine-independent calls between c and fortran.  The 
 * package is available via anonymous ftp at zebra.desy.de.
 *
 * The grib document "A GUIDE TO THE CODE FORM FM 92-IX Ext. GRIB" may be
 * useful to your understanding of this code.  This document is available 
 * via anonymous ftp from nic.fb4.noaa.gov.  Check the readme file in the
 * root directory for further instructions.
 *
 ****************************************************************************/

#define ERRSIZE 2000
#define ALLOCSIZE 30
#define MISSING -999

#define EARTH_RADIUS 6371.229        /* in km */
#define PI           3.141592654
#define PI_OVER_180  PI/180.

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#ifdef MACOS
#include "/usr/include/time.h"
#else
#include <time.h>
#endif
#include "cfortran.h"
#include "gribfuncs.h"
#include "gribsize.incl"
#include "read_grib.h"

/* Function Declarations */

void remove_element(int array[],int index, int size);
int advance_time(int *century, int year, int month, int day, int hour,
                 int amount, int unit);
char *advance_time_str(char startdatein[], int amount, char enddate[]);
int date_diff(int date1,int century1,int date2,int century2);
int hours_since_1900(int date,int century);
int isLeapYear(int year);
int get_factor2(int unit);
int compare_record(GribInfo *gribinfo, FindGrib *findgrib, int gribnum);

/* 
 *These lines allow fortran routines to call the c routines.  They are
 * used by macros defined in cfortran.h
 */
#define get_pressure_levels_STRV_A1 TERM_CHARS(' ',1)
/*
FCALLSCFUN6(INT, get_pressure_levels,GET_PRESSURE_LEVELS,
            get_pressure_levels,STRINGV,INTV,INTV,INTV,INT,INT)
#define setup_gribinfo_STRV_A1 TERM_CHARS(' ',1)
FCALLSCFUN2(INT,setup_gribinfo,SETUP_GRIBINFO,setup_gribinfo,STRINGV,INT)
#define get_msl_indices_STRV_A1 TERM_CHARS(' ',1)
FCALLSCFUN9(INT, get_msl_indices,GET_MSL_INDICES,get_msl_indices,
            STRINGV,INTV,INTV,INTV,INTV,INTV,INT,INTV,INTV)
FCALLSCFUN5(INT, get_index,GET_INDEX,get_index,INT,INT,INT,INT,INT)
#define read_grib_STRV_A1 TERM_CHARS(' ',1)
FCALLSCFUN7(INT,get_dates,GET_DATES,get_dates,INTV,INTV,INTV,INT,INTV,
            INTV,INTV)
FCALLSCFUN7(INT, read_grib,READ_GRIB,read_grib,
            STRINGV,INT,INT,INT,INT,FLOATVV,PVOID)
FCALLSCFUN8(INT, get_index_near_date,GET_INDEX_NEAR_DATE,get_index_near_date,
            STRING,INT,INT,INT,INTV,INTV,INTV,INT)
*/
/* The value for usLevel_id for isobaric levels */
#define ISOBARIC_LEVEL_ID 100  

/*************************************************************************
 * This function reads and decodes grib records in a list of input files 
 * and stores information about each grib record in the gribinfo array 
 * structure.  The gribinfo structure can then be accessed by any function 
 * within this file.
 * 
 * Interface:
 *   Input:
 *      gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *      files - a string array containing the names of the files containing
 *              the grib data.  If called from a fortran routine, the
 *              fortran routine must set the character size of the array to
 *              be STRINGSIZE-1. The last filled element in the array should
 *              be "END".
 *      use_fcst - if TRUE, forecast fields will be included in the gribinfo
 *              structure, otherwise, only analysis fields will be included.
 *
 *    Return:
 *       1  - successful call to setup_gribinfo
 *       -1 - call to setup_gribinfo failed
 *
 ***************************************************************************/

int rg_setup_gribinfo(GribInfo *gribinfo, char files[][STRINGSIZE],
                      int use_fcst)
{
  FILE *fp;
  int filenum;
  int nReturn;
  int idx;
  int status;
  int start_elem;
  
  /* Loop through input files */
  filenum = 0;
  while ((strcmp(files[filenum], "end") != 0 ) && 
         (strcmp(files[filenum], "END") != 0 )) {
    
    /* 
     * This forces gribinfo to be fully initialized.
     */
    if (filenum == 0) 
      {
        gribinfo->num_elements = 0;
      }

    start_elem = gribinfo->num_elements;

    fp = fopen(files[filenum],"r");
    if (fp == NULL)
      {
        fprintf(stderr,"Could not open %s\n",files[filenum]);
        nReturn = -1;
        break;
      }
    
    status = rg_setup_gribinfo_f(gribinfo, fp, use_fcst);
    if (status != 1) 
      {
        fprintf(stderr, 
            "rg_setup_gribinfo_f returned non-zero status (%d), skipping %s\n",
                status,files[filenum]);
        continue;
      }

    for (idx=start_elem; idx < gribinfo->num_elements; idx++) 
      {
        strcpy(gribinfo->elements[idx].filename,
               files[filenum]);
      }

    
    filenum++;
    nReturn = 1;
  }

  return nReturn;
}



/*************************************************************************
 *
 * Similar to rg_setup_gribinfo, except, a unix file descriptor is passed in, 
 * rather than a list of files to open.
 *
 *************************************************************************/

int rg_setup_gribinfo_i(GribInfo *gribinfo, int fid, int use_fcst)
{
  FILE *fp;
  int status;
  
  fp = fdopen(fid,"r");
  if (fp == NULL)
    {
      fprintf(stderr,"Could not open file descriptor %d\n",fid);
      status = -1;
      return status;
    }
  
  /* This forces gribinfo to be initialized for the first time */
  gribinfo->num_elements = 0;
  
  status = rg_setup_gribinfo_f(gribinfo, fp, use_fcst);
  if (status != 1) 
    {
      fprintf(stderr, 
              "rg_setup_gribinfo_f returned non-zero status (%d)\n",
              status);
    }
  
  return status;
}

/*************************************************************************
 *
 * Similar to rg_setup_gribinfo, except, a file pointer is passed in, rather
 * than a list of files to open.
 * 
 * If gribinfo->num_elements is 0, gribinfo is initialized, otherwise, 
 *    gribinfo is appended to.
 *
 *************************************************************************/

int rg_setup_gribinfo_f(GribInfo *gribinfo, FILE *fp, int use_fcst)
{
  char errmsg[ERRSIZE];
  int nReturn=0;
  long offset;
  int filenum;
  int Rd_Indexfile=0;
  GRIB_HDR *gh1;
  long tmpoffset=0;
  int century;
  int year4d;
  int fcsttime1=0;
  int fcsttime2=0;
  int factor=0;

  /* Set the number of elements to be zero initially */
  if (gribinfo->num_elements <= 0)
    {
      /* Allocate space for gribinfo */
      gribinfo->elements = (Elements *)calloc(ALLOCSIZE,sizeof(Elements));
      if (gribinfo->elements == NULL) {
        sprintf(errmsg,"Could not allocate %d bytes for gribinfo->elements\n",
                ALLOCSIZE*sizeof(Elements));
        goto bail_out;
      }
    }
  
  /* Make storage for Grib Header */
  nReturn = init_gribhdr(&gh1, errmsg);
  /* 
   * The grib library is setup such that, when init_gribhdr == 0, it was 
   * successful.  If it is 1, it failed.
   */
  if (nReturn == 1) goto bail_out;
  
  /* Scan through message */
  for (offset = 0L; nReturn == 0; offset += gh1->msg_length) {
    if ((gribinfo->num_elements > 0) && 
        (gribinfo->num_elements%ALLOCSIZE == 0))
      gribinfo->elements = 
        (Elements *)realloc(gribinfo->elements,
                            (gribinfo->num_elements+ALLOCSIZE)*
                            sizeof(Elements));

    if (gribinfo->elements == NULL) {
      sprintf(errmsg,"Could not allocate %d bytes for gribinfo\n",
              (gribinfo->num_elements + ALLOCSIZE)*sizeof(Elements));
      goto bail_out;
    }

    /* Setup the File pointer */
    gribinfo->elements[gribinfo->num_elements].fp = fp;

    gribinfo->elements[gribinfo->num_elements].pds = 
      (PDS_INPUT *)malloc(1*sizeof(PDS_INPUT));
    gribinfo->elements[gribinfo->num_elements].gds = 
      (grid_desc_sec *)malloc(1*sizeof(grid_desc_sec));
    gribinfo->elements[gribinfo->num_elements].bms = 
      (BMS_INPUT *)malloc(1*sizeof(BMS_INPUT));
    gribinfo->elements[gribinfo->num_elements].bds_head = 
      (BDS_HEAD_INPUT *)malloc(1*sizeof(BDS_HEAD_INPUT));
    errmsg[0] = '\0';
    nReturn = 
      grib_fseek(fp,&offset, Rd_Indexfile, gh1, errmsg);
    if (nReturn != 0) {
      if (nReturn == 2) break;   /* End of file error */
      else {
        fprintf(stderr, "Grib_fseek returned non zero status (%d)\n",
                nReturn);
        goto bail_out;
      }
    }
    if (errmsg[0] != '\0')
      { /* NO errors but got a Warning msg from seek */
        fprintf(stderr,"%s; Skip Decoding...\n",errmsg);
        errmsg[0] = '\0';   
        gh1->msg_length = 1L;    /* set to 1 to bump offset up */
        continue;
      }
      
    if (gh1->msg_length < 0) {
      fprintf(stderr, "Error:  message returned had bad length (%ld)\n",
              gh1->msg_length);
      goto bail_out;
    }
    else if (gh1->msg_length == 0) {
      fprintf(stderr, "msg_length is Zero\n");
      gh1->msg_length = 1L;
      continue;
    }
    init_dec_struct(gribinfo->elements[gribinfo->num_elements].pds, 
                    gribinfo->elements[gribinfo->num_elements].gds, 
                    gribinfo->elements[gribinfo->num_elements].bms, 
                    gribinfo->elements[gribinfo->num_elements].bds_head);
    
    /* 
     * gribgetpds is an undocumented function within the grib library.
     * gribgetpds grabs the pds section from the grib message without
     * decoding the entire grib message.  The interface is as follows:
     *     first input param: a pointer to the beginning of the pds 
     *                        section.
     *     second input param: a pointer to a structure which will hold
     *                        the pds information
     *     third param: the error message.
     * 
     * If gribgetpds ever fails, it can be replaced with the following 
     *    nReturn = grib_dec((char *)gh1->entire_msg, &pds, &gds, &bds_head, 
     *              &bms, &grib_data, errmsg);
     * 
     * This will degrade performance since this grib_dec decodes the 
     *    entire grib message.
     */
    
    nReturn = gribgetpds((char*)(gh1->entire_msg + 8),
                         gribinfo->elements[gribinfo->num_elements].pds,
                         errmsg);
    if (nReturn != 0) goto bail_out;

    /* Get gds if present */
    if (gribinfo->elements[gribinfo->num_elements].pds->usGds_bms_id >> 7 
        & 1) {
      nReturn = 
        gribgetgds((char*)
                   (gh1->entire_msg+8+
                    gribinfo->elements[gribinfo->num_elements].pds->uslength),
                   gribinfo->elements[gribinfo->num_elements].gds,errmsg);
      if (nReturn != 0) goto bail_out;
    }

    /* Get bms section if present */
    if (gribinfo->elements[gribinfo->num_elements].pds->usGds_bms_id >> 6 
        & 1) {
      /*
        fprintf(stderr,"grids with bms section not currently supported\n");
        return -1;
      */
    }
    
    gribinfo->elements[gribinfo->num_elements].usGrid_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usGrid_id;
    gribinfo->elements[gribinfo->num_elements].usParm_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usParm_id;
    gribinfo->elements[gribinfo->num_elements].usLevel_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usLevel_id;
    gribinfo->elements[gribinfo->num_elements].usHeight1 = 
      gribinfo->elements[gribinfo->num_elements].pds->usHeight1;
    gribinfo->elements[gribinfo->num_elements].usHeight2 = 
      gribinfo->elements[gribinfo->num_elements].pds->usHeight2;
    gribinfo->elements[gribinfo->num_elements].center_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usCenter_id;
    gribinfo->elements[gribinfo->num_elements].parmtbl = 
      gribinfo->elements[gribinfo->num_elements].pds->usParm_tbl;
    gribinfo->elements[gribinfo->num_elements].proc_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usProc_id;
    gribinfo->elements[gribinfo->num_elements].subcenter_id = 
      gribinfo->elements[gribinfo->num_elements].pds->usCenter_sub;
    gribinfo->elements[gribinfo->num_elements].offset = offset;
    gribinfo->elements[gribinfo->num_elements].end = 
      offset + gh1->msg_length - 1;
  
    if (use_fcst) {
      century = gribinfo->elements[gribinfo->num_elements].pds->usCentury;
      
      if (gribinfo->elements[gribinfo->num_elements].pds->usTime_range == 10)
        {
          fcsttime1 = gribinfo->elements[gribinfo->num_elements].pds->usP1*256 + 
            gribinfo->elements[gribinfo->num_elements].pds->usP2;
          fcsttime2 = 0;
        }
      else if (gribinfo->elements[gribinfo->num_elements].pds->usTime_range 
               == 203) {
        /* This is the WSI extension to grib.  203 indicates "duration" */
        fcsttime1 = gribinfo->elements[gribinfo->num_elements].pds->usP1;
        fcsttime2 = gribinfo->elements[gribinfo->num_elements].pds->usP1 +
          gribinfo->elements[gribinfo->num_elements].pds->usP2;
      } else {
        fcsttime1 = gribinfo->elements[gribinfo->num_elements].pds->usP1;
        fcsttime2 = gribinfo->elements[gribinfo->num_elements].pds->usP2;
      }
      
      gribinfo->elements[gribinfo->num_elements].date = 
        advance_time(&century,
                     gribinfo->elements[gribinfo->num_elements].pds->usYear,
                     gribinfo->elements[gribinfo->num_elements].pds->usMonth,
                     gribinfo->elements[gribinfo->num_elements].pds->usDay,
                     gribinfo->elements[gribinfo->num_elements].pds->usHour,
                     fcsttime1,
                     gribinfo->elements[gribinfo->num_elements].pds->usFcst_unit_id);
    } 
    else {
      gribinfo->elements[gribinfo->num_elements].date = 
        gribinfo->elements[gribinfo->num_elements].pds->usHour*1 + 
        gribinfo->elements[gribinfo->num_elements].pds->usDay*100 + 
        gribinfo->elements[gribinfo->num_elements].pds->usMonth*10000 + 
        gribinfo->elements[gribinfo->num_elements].pds->usYear*1000000;
    } 
    gribinfo->elements[gribinfo->num_elements].century = 
      gribinfo->elements[gribinfo->num_elements].pds->usCentury;
    
    year4d = 
        (gribinfo->elements[gribinfo->num_elements].pds->usCentury - 1) * 100
        + gribinfo->elements[gribinfo->num_elements].pds->usYear;

    sprintf(gribinfo->elements[gribinfo->num_elements].initdate,
            "%04d%02d%02d%02d%02d%02d",
            year4d,
            gribinfo->elements[gribinfo->num_elements].pds->usMonth,
            gribinfo->elements[gribinfo->num_elements].pds->usDay,
            gribinfo->elements[gribinfo->num_elements].pds->usHour,
            gribinfo->elements[gribinfo->num_elements].pds->usMinute,
            0);       
    
    factor = 
      get_factor2(gribinfo->elements[gribinfo->num_elements].pds->usFcst_unit_id);
    gribinfo->elements[gribinfo->num_elements].fcsttime1 = 
      fcsttime1 * factor;
    gribinfo->elements[gribinfo->num_elements].fcsttime2 = 
      fcsttime2 * factor;
    
    advance_time_str(gribinfo->elements[gribinfo->num_elements].initdate,
                     gribinfo->elements[gribinfo->num_elements].fcsttime1,
                     gribinfo->elements[gribinfo->num_elements].valid_time);

    gribinfo->num_elements++;
  }
  
  free_gribhdr(&gh1);
  return 1;
  
  /* The error condition */
 bail_out:
  if (errmsg[0] != '\0') fprintf(stderr,"\n***ERROR: %s: %s\n",
                                 "setup_grib",errmsg);
  if (gribinfo->elements != NULL) free(gribinfo->elements);
  perror("System Error ");
  return -1;
}

/*****************************************************************************
 *
 * Retrieve pressure levels from grib data.  This function will pass the 
 * pressure levels for which the input parameter is available at all input 
 * times back to the calling routine.  
 *
 * Interface
 *      Input: 
 *         gribinfo - pointer to a previously allocated gribinfo structure.  
 *                 The gribinfo structure is filled in this function.
 *         dates: an array of dates to check for data
 *                format: yymmddhh
 *                If called from a fortran routine, the fortran routine must 
 *                set the character size of the array to be STRINGSIZE-1
 *         centuries: an array holding the centuries for each of the 
 *                dates in the array dates.
 *         parm_id: the input parameter id.  From table 2 of the grib manual.
 *      Output:
 *         finallevels: an array of pressure levels which are contained in 
 *                the grib data at all input times.
 *      Return:
 *         the number of levels in the levels array.  The levels are listing
 *         in descending (by value) order, i.e., the value with the highest
 *         pressure (lowest vertical level) is the first element.
 *
 ****************************************************************************/

int rg_get_pressure_levels(GribInfo *gribinfo, int dates[], int centuries[], 
                        int parm_id[], int finallevels[],int min_pres,
                        int numparms)
{
  int datenum;
  int gribnum;
  int *levelnum;
  int levelincluded;
  int i,j;
  int contains_level;
  int **tmplevels;
  int numfinallevels = 0;
  int parmnum;
  int tmpval;
  
  /* Allocate space */
  levelnum = (int *)calloc(numparms,sizeof(int));
  tmplevels = (int **)calloc(numparms,sizeof(int *));
  for (j = 0; j < numparms; j++) {
    tmplevels[j] = (int *)calloc(1000,sizeof(int));
    if (tmplevels[j] == NULL) {
      tmplevels = NULL;
      break;
    }
  }
  if ((levelnum == NULL) || (tmplevels == NULL)) {
    fprintf(stderr,
            "get_pressure_levels: Allocation of space failed, returning\n");
    return -1;
  }
  
  /* Loop through all parameters */
  for (parmnum = 0; parmnum < numparms; parmnum++) {

    levelnum[parmnum] = 0;

  /* Get the array of pressure levels available at the first input time */
    datenum = 0;
    for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
      if (gribinfo->elements[gribnum].date == dates[datenum]) {
        if (gribinfo->elements[gribnum].century == centuries[datenum]) {
          if (gribinfo->elements[gribnum].usLevel_id == ISOBARIC_LEVEL_ID) {
            if (gribinfo->elements[gribnum].usParm_id == parm_id[parmnum]) {
              if (gribinfo->elements[gribnum].usHeight1 >= min_pres) {
                levelincluded = 0;
                for (j=0; j < levelnum[parmnum]; j++) {
                  if (tmplevels[parmnum][j] == 
                      gribinfo->elements[gribnum].usHeight1) {
                    levelincluded = 1;
                    break;
                  }
                }
                if (levelincluded == 0) {
                  tmplevels[parmnum][levelnum[parmnum]] = 
                    gribinfo->elements[gribnum].usHeight1;
                  levelnum[parmnum]++;
                }
              }
            }
          }
        }
      }
    }
    
    /* Remove levels that are not contained at all subsequent times */
    datenum++;
    while (dates[datenum] != -99){
      for (j = 0; j < levelnum[parmnum]; j++) {
        contains_level = 0;
        for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
          if (gribinfo->elements[gribnum].date == dates[datenum]) {
            if (gribinfo->elements[gribnum].century == centuries[datenum]) {
              if (gribinfo->elements[gribnum].usLevel_id == ISOBARIC_LEVEL_ID)
                {
                  if (gribinfo->elements[gribnum].usParm_id == 
                      parm_id[parmnum]) {
                    if (tmplevels[parmnum][j] == 
                        gribinfo->elements[gribnum].usHeight1)
                      contains_level = 1;
                  }
                }
            }
          }
        }
        if (!(contains_level)) {
          remove_element(tmplevels[parmnum],j,levelnum[parmnum]);
          levelnum[parmnum]--;
          j--;
        }
      }
      datenum++;
    }

    /* 
     * Put the values for levels into an array.  Remove any levels that
     * were not found at all other levels
     */
    if (parmnum == 0) {
      for (j = 0; j < levelnum[parmnum]; j++) {
        finallevels[j] = tmplevels[parmnum][j];
        numfinallevels++;
      }
    } else {
      for (i=0; i<numfinallevels; i++) {
        contains_level = 0;
        for (j=0; j<levelnum[parmnum]; j++) {
          if (finallevels[i] == tmplevels[parmnum][j]) {
            contains_level = 1;
            break;
          }
        }
        if (!contains_level) {
          remove_element(finallevels,i,numfinallevels);
          numfinallevels--;
          i--;
        }
      }
    }

  }

  /* 
   * Sort the numfinallevels array into descending order. Use straight 
   * insertion.
   */
  for (j=1; j<numfinallevels; j++) {
    tmpval = finallevels[j];
    for (i=j-1; i >= 0; i--) {
      if (finallevels[i] >= tmpval) break;
      finallevels[i+1] = finallevels[i];
    }
    finallevels[i+1] = tmpval;
  }

  return numfinallevels;
}

/****************************************************************************
 *
 * Returns an array of grib indices that correspond to particular grib fields
 * to use as sea level pressure.  There will be exactly one element for each
 * input time.  If a field was not found, then this function returns NULL
 *
 * Interface:
 *    Input:
 *      gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *         dates: a string array of dates to check for data.
 *                format: yymmddhh
 *                If called from a fortran routine, the fortran routine must 
 *                set the character size of the array to be STRINGSIZE-1
 *         centuries: an array holding the centuries for each of the 
 *                dates in the array dates.
 *         usParm_id: an array of parameter identifiers that could be
 *                    used as a sea level pressure field (From table 2 of 
 *                    grib documentation)
 *         usLevel_id: the level id that could be used as a sea level pressure
 *                    field (from table 3 of the grib documentation)
 *         usHeight1: the height for the particular parameter and level
 *                    (in units described by the parameter index)
 *         numparms:  the number of parameters in each of the usParm_id, 
 *                    usLevel_id, and usHeight1 arrays.
 *    Output:
 *         grib_index: an array of grib indices to use for the sea level 
 *                    pressure.  The index to grib_index corresponds to 
 *                    the time, i.e., the first array element of grib_index 
 *                    corresponds to the first time, the second element to 
 *                    the second time, etc.
 *         
 *            Note: Values in the input arrays, usParm_id, usLevel_id, and 
 *                    usHeight with the same array index must correspond.
 *
 *    Return:
 *          1 for success
 *          -1 if no field was found.
 ***************************************************************************/

int rg_get_msl_indices(GribInfo *gribinfo, char dates[][STRINGSIZE],
                    int centuries[], int usParm_id[],int usLevel_id[],
                    int usHeight1[],int infactor[],int numparms,
                    int grib_index[],int outfactor[])
{
  int parmindex;
  int datenum = 0;
  int gribnum;
  int foundfield=0;

  for (parmindex = 0; parmindex < numparms; parmindex++) {

    datenum = 0;
    while ((strcmp(dates[datenum], "end") != 0 ) && 
           (strcmp(dates[datenum], "END") != 0 )) {

      for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
        if (gribinfo->elements[gribnum].date == atoi(dates[datenum])) {
          if (gribinfo->elements[gribnum].century == centuries[datenum]) {
            if ((gribinfo->elements[gribnum].usParm_id == 
                 usParm_id[parmindex]) &&
                (gribinfo->elements[gribnum].usLevel_id == 
                 usLevel_id[parmindex]) &&
                (gribinfo->elements[gribnum].usHeight1 == 
                 usHeight1[parmindex])) {
              grib_index[datenum] = gribnum;
              outfactor[datenum] = infactor[parmindex];
              foundfield++;
              break;
            }
          }
        }
      }

      datenum++;

      /* 
       * Break out of loop and continue on to next parameter if the current 
       * parameter was missing from a date. 
       */

      if (foundfield != datenum) break;
    }

/* 
 * Break out of the parameter loop once we've found a field available at all 
 * dates 
 */
    if (foundfield == datenum) {
      break;
    }

  }

  if (foundfield == datenum)
    return 1;
  else
    return -1;

}


/***************************************************************************
 * 
 * This function takes an index as input and returns a 2d grib data field
 *
 * Interface:
 *    input: 
 *      gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *       index - the index of gribinfo for which data is to be retrieved
 *       scale - the scale factor to multiply data by, i.e., if -2, 
 *               data will be multiplied by 10^-2.
 *    output:
 *       grib_out - the 2 dimensional output grib data
 *               Warning: This 2d array is setup with i being the vertical
 *               dimension and j being the horizontal dimension.  This
 *               is the convention used in mesoscale numerical modeling
 *               (the MM5 in particular), so it is used here.
 *    return:
 *       1 for success
 *      -1 for failure
 ***************************************************************************/

int rg_get_grib(GribInfo *gribinfo, int index,int scale,
             float **grib_out,int *vect_comp_flag, 
             GRIB_PROJECTION_INFO_DEF *Proj, BDS_HEAD_INPUT *bds_head)
{
  char errmsg[ERRSIZE];
  int nReturn=0;
  long offset;
  int Rd_Indexfile=0;
  BMS_INPUT bms;
  PDS_INPUT pds;
  BDS_HEAD_INPUT dummy;
  grid_desc_sec gds;
  GRIB_HDR *gh1;
  int i,j;
  int expandlon = 0;
  float *grib_data;

  /* Initialize Variables */
  errmsg[0] = '\0';
  offset = 0L;
  grib_data = (float *)NULL;

  /* Make storage for Grib Header */
  nReturn = init_gribhdr (&gh1, errmsg);
  if (nReturn == 1) goto bail_out;
  
  /* Seek to the position in the grib data */
  offset = gribinfo->elements[index].offset;
  nReturn = grib_fseek(gribinfo->elements[index].fp,&offset,
                      Rd_Indexfile,gh1,errmsg);
  if (nReturn != 1) {
    fprintf(stderr,"Grib_fseek returned error status (%d)\n",nReturn);
    goto bail_out;
  }
  if (errmsg[0] != '\0')
    { /* NO errors but got a Warning msg from seek */ 
      fprintf(stderr,"%s: Skip Decoding...\n",errmsg);
      errmsg[0] = '\0';
    }
  if (gh1->msg_length <= 0) {
    fprintf(stderr,"Error: message returned had bad length (%ld)\n",
            gh1->msg_length);
    goto bail_out;
  }
  init_dec_struct(&pds, &gds, &bms, &dummy);
 
  nReturn = grib_dec((char *)gh1->entire_msg, &pds, &gds, 
                     bds_head,
                     &bms, &grib_data, errmsg);
 
  if (nReturn != 0) goto bail_out;
  
  if (bms.uslength > 0) {
    nReturn = apply_bitmap(&bms, &grib_data, FILL_VALUE, bds_head, 
                           errmsg);
    if (nReturn != 0) goto bail_out;
  }
 
  switch(gds.head.usData_type) {
  case 0:
  case 4:
    strcpy(Proj->prjnmm,"latlon");
    Proj->colcnt = gds.llg.usNi;
    Proj->rowcnt = gds.llg.usNj;
    Proj->origlat = gds.llg.lLat1/1000.;
    Proj->origlon = gds.llg.lLon1/1000.;
    Proj->xintdis = (gds.llg.iDi/1000.)*EARTH_RADIUS*PI_OVER_180;
    Proj->yintdis = (gds.llg.iDj/1000.)*EARTH_RADIUS*PI_OVER_180;
    Proj->parm1 = 0.;
    Proj->parm2 = 0.;
    if ((gds.llg.usRes_flag >> 3) & 1) *vect_comp_flag = 1;
    else *vect_comp_flag = 0;

    /* If the grid is a global grid, we want to set the expandlon flag
     * so that the number of columns in the array is expanded by one and
     * the first column of data is copied to the last column.  This
     * allows calling routines to interpolate between first and last columns
     * of data.
     */

    if (gds.llg.usNi*gds.llg.iDi/1000. == 360)
      expandlon = 1;
    else
      expandlon = 0;
     
    break;
  case 1: 
    strcpy(Proj->prjnmm,"mercator");
    Proj->colcnt = gds.merc.cols;
    Proj->rowcnt = gds.merc.rows;
    Proj->origlat = gds.merc.first_lat/1000.;
    Proj->origlon = gds.merc.first_lon/1000.;
    Proj->xintdis = gds.merc.lon_inc/1000.;
    Proj->yintdis = gds.merc.lat_inc/1000.;
    Proj->parm1 = gds.merc.latin/1000.;
    Proj->parm2 = (gds.merc.Lo2/1000. - Proj->origlon)/gds.merc.cols;
    if ((gds.merc.usRes_flag >> 3) & 1) *vect_comp_flag = 1;
    else *vect_comp_flag = 0;
    break;
  case 3:
    strcpy(Proj->prjnmm,"lambert");
    Proj->colcnt = gds.lam.iNx;
    Proj->rowcnt = gds.lam.iNy;
    Proj->origlat = gds.lam.lLat1/1000.;
    Proj->origlon = gds.lam.lLon1/1000.;
    Proj->xintdis = gds.lam.ulDx/1000.;
    Proj->yintdis = gds.lam.ulDy/1000.;
    Proj->parm1 = gds.lam.lLat_cut1/1000.;
    Proj->parm2 = gds.lam.lLat_cut2/1000.;
    Proj->parm3 = gds.lam.lLon_orient/1000.;
    if ((gds.lam.usRes_flag >> 3) & 1) *vect_comp_flag = 1;
    else *vect_comp_flag = 0;
    break;
  case 5:
    strcpy(Proj->prjnmm,"polar_stereo");
    Proj->colcnt = gds.pol.usNx;
    Proj->rowcnt = gds.pol.usNy;
    Proj->origlat = gds.pol.lLat1/1000.;
    Proj->origlon = gds.pol.lLon1/1000.;
    Proj->xintdis = gds.pol.ulDx/1000.;
    Proj->yintdis = gds.pol.ulDy/1000.;
    Proj->parm1 = 60.;
    Proj->parm2 = gds.pol.lLon_orient/1000.;   
    if ((gds.pol.usRes_flag >> 3) & 1) *vect_comp_flag = 1;
    else *vect_comp_flag = 0;
    break;
  default:
    fprintf(stderr,"Grid not supported, gds.head.usData_type = %d\n",
            gds.head.usData_type);
    fprintf(stderr,"Exiting\n");
    exit(-1);
    break;
  }  

  strcpy(Proj->stordsc,"+y_in_+x");
  Proj->origx = 1;
  Proj->origy = 1;

  for (j=0; j< (Proj->rowcnt); j++) {
    for (i=0; i<(Proj->colcnt); i++) {
      grib_out[j][i] = grib_data[i+j*Proj->colcnt]*pow(10,scale);
    }
  }

  if (expandlon) {
    (Proj->colcnt)++;
    for (j = 0; j < Proj->rowcnt; j++) {
      grib_out[j][Proj->colcnt-1] = grib_out[j][0];
    }
  }

  /* 
   * You only reach here when there is no error, so return successfully.
   */
      
  nReturn = 0;

  if (grib_data != NULL) {
    free_gribhdr(&gh1);
    free(grib_data);
  } 

  return 1;
      
  /* The error condition */
 bail_out:
  if (errmsg[0] != '\0') fprintf(stderr,"\n***ERROR: %s %s\n",
                                 "get_grib",errmsg);
  if (grib_data != NULL)
    free(grib_data);
  free_gribhdr(&gh1);
  return -1;
}

/***************************************************************************
 * 
 * This function takes an index as input and returns a 2d grib data field
 *
 * Interface:
 *    input: 
 *      gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *       index - the index of gribinfo for which data is to be retrieved
 *    output:
 *       data - the 2 dimensional output grib data
 *               Warning: This 2d array is setup with i being the vertical
 *               dimension and j being the horizontal dimension.  This
 *               is the convention used in mesoscale numerical modeling
 *               (the MM5 in particular), so it is used here.
 *    return:
 *       1 for success
 *      -1 for failure
 ***************************************************************************/

int rg_get_data(GribInfo *gribinfo, int index, float **data)
{
  float *data_1d;
  int i,j;
  int numrows,numcols;
  int status;

  numrows = rg_get_numrows(gribinfo,index);
  numcols = rg_get_numcols(gribinfo,index);
  
  data_1d = (float *)calloc(numrows*numcols,sizeof(float));
  if (data_1d == 0)
    {
      fprintf(stderr,"Allocating space for data_1d failed, index: %d\n",index);
      return -1;
    }

  status = rg_get_data_1d(gribinfo, index, data_1d);
  if (status != 1)
    {
      return status;
    }

  for (j=0; j< numrows; j++) {
    for (i=0; i < numcols; i++) {
      data[j][i] = data_1d[i+j*numcols];
    }
  }

  free(data_1d);

  return 1;
  
}

/***************************************************************************
 * 
 * This function takes an index as input and returns a 1d grib data field
 *
 * Interface:
 *    input: 
 *      gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *       index - the index of gribinfo for which data is to be retrieved
 *    output:
 *       data - 1 dimensional output grib data
 *               Warning: This 2d array is setup with i being the vertical
 *               dimension and j being the horizontal dimension.  This
 *               is the convention used in mesoscale numerical modeling
 *               (the MM5 in particular), so it is used here.
 *    return:
 *       1 for success
 *      -1 for failure
 ***************************************************************************/

int rg_get_data_1d(GribInfo *gribinfo, int index, float *data)
{
  char errmsg[ERRSIZE];
  int nReturn=0;
  long offset;
  int Rd_Indexfile=0;
  BMS_INPUT bms;
  PDS_INPUT pds;
  BDS_HEAD_INPUT bds_head;
  grid_desc_sec gds;
  GRIB_HDR *gh1;
  int i,j;
  int numcols, numrows;
  float *grib_data;

  /* Initialize Variables */
  errmsg[0] = '\0';
  offset = 0L;
  grib_data = (float *)NULL;

  /* Make storage for Grib Header */
  nReturn = init_gribhdr (&gh1, errmsg);
  if (nReturn == 1) goto bail_out;
  
  /* Seek to the position in the grib data */
  offset = gribinfo->elements[index].offset;
  nReturn = grib_fseek(gribinfo->elements[index].fp,&offset,
                      Rd_Indexfile,gh1,errmsg);
  if (nReturn != 0) {
    fprintf(stderr,"Grib_fseek returned non-zero status (%d)\n",nReturn);
    goto bail_out;
  }
  if (errmsg[0] != '\0')
    { /* NO errors but got a Warning msg from seek */ 
      fprintf(stderr,"%s: Skip Decoding...\n",errmsg);
      errmsg[0] = '\0';
    }
  if (gh1->msg_length <= 0) {
    fprintf(stderr,"Error: message returned had bad length (%ld)\n",
            gh1->msg_length);
    goto bail_out;
  }

  init_dec_struct(&pds, &gds, &bms, &bds_head);
 
  nReturn = grib_dec((char *)gh1->entire_msg, &pds, &gds, 
                     &bds_head,
                     &bms, &grib_data, errmsg);
 
  if (nReturn != 0) goto bail_out;
  
  if (bms.uslength > 0) {
    nReturn = apply_bitmap(&bms, &grib_data, FILL_VALUE, &bds_head, 
                           errmsg);
    if (nReturn != 0) goto bail_out;
  }
 
  /*
   * Copy the data into the permanent array
   */ 
  numcols = rg_get_numcols(gribinfo,index);
  numrows = rg_get_numrows(gribinfo,index);
  memcpy(data,grib_data,numcols*numrows*sizeof(float));

  /* 
   * You only reach here when there is no error, so return successfully.
   */
      
  nReturn = 0;

  if (grib_data != NULL) {
    free_gribhdr(&gh1);
    free(grib_data);
  } 

  return 1;
      
  /* The error condition */
 bail_out:
  if (errmsg[0] != '\0') fprintf(stderr,"\n***ERROR: %s %s\n",
                                 "get_grib",errmsg);
  if (grib_data != NULL)
    free(grib_data);
  free_gribhdr(&gh1);
  return -1;
}

/****************************************************************************
 * Returns the index of gribinfo corresponding to the input date, level, 
 * height, and parameter.
 *
 * Interface:
 *   Input: 
 *      gribinfo  - pointer to a previously populated gribinfo structure.  
 *      initdate  -  initialization date in the form yyyymmdd[HHMMSS].  If any
 *                  part of HHMMSS is not specified, it will be set to 0.
 *      parmid    - the parameter id in the grib file
 *      leveltype - the leveltype id from table 3/3a of the grib document.
 *      level1    - First level of the data in units described by leveltype.
 *      level2    - Second level of the data in units described by leveltype.
 *      fcsttime1 - First forecast time in seconds.
 *      fcsttime2 - Second forecast time in seconds.
 *    Note: If an input variable is set set to -INT_MAX, then any value
 *                will be considered a match.
 *    Return:
 *       if >= 0    The index of the gribinfo data that corresponds to the
 *                  input parameters
 *       if < 0     No field corresponding to the input parms was found.
 *
 ***************************************************************************/

int rg_get_index(GribInfo *gribinfo, FindGrib *findgrib)
{
  int gribnum;
  int grib_index=-1;

  for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
    if (compare_record(gribinfo, findgrib, gribnum) == 1)
      {
        grib_index = gribnum;
        break;
      }
  }
  return grib_index;
}


/****************************************************************************
 * Same as rg_get_index, except that a guess for the record number is given.
 *   This "guess" record is first checked to see if it matches, if so, 
 *      that grib record number is just returned.  If it does not match, 
 *      full searching ensues.
 * Returns the index of gribinfo corresponding to the input date, level, 
 * height, and parameter.
 *
 * Interface:
 *   Input: 
 *      Same is rg_get_index, except:
 *      guess_index  - The index to check first.
 *    Return:
 *      Same as rg_get_index
 *
 ***************************************************************************/

int rg_get_index_guess(GribInfo *gribinfo, FindGrib *findgrib, int guess_index)
{
  int retval;

  if (compare_record(gribinfo, findgrib, guess_index) == 1) {
    retval = guess_index;
  } else {
    retval = rg_get_index(gribinfo, findgrib);
  }

  return retval;
}


/****************************************************************************
 * Sets all values in FindGrib to missing.
 *
 * Interface:
 *   Input: 
 *      findgrib  - pointer to a previously allocated findgrib structure.
 *
 *    Return:
 *       1  for success.
 *       -1 for failure.
 *
 ***************************************************************************/
int rg_init_findgrib(FindGrib *findgrib)
{
  strcpy(findgrib->initdate,"*");
  strcpy(findgrib->validdate,"*");
  findgrib->parmid          = -INT_MAX;
  findgrib->parmid          = -INT_MAX;
  findgrib->leveltype       = -INT_MAX;
  findgrib->level1          = -INT_MAX;
  findgrib->level2          = -INT_MAX;
  findgrib->fcsttime1       = -INT_MAX;
  findgrib->fcsttime2       = -INT_MAX;
  findgrib->center_id       = -INT_MAX;
  findgrib->subcenter_id    = -INT_MAX;
  findgrib->parmtbl_version = -INT_MAX;
  
  return 1;
}

/****************************************************************************
 * Returns the indices of all gribinfo entries that match the input date, 
 * level, height, and parameter.
 *
 * Interface:
 *   Input: 
 *      gribinfo  - pointer to a previously populated gribinfo structure.  
 *      initdate  -  initialization date in the form yyyymmdd[HHMMSS].  If any
 *                  part of HHMMSS is not specified, it will be set to 0.
 *      parmid    - the parameter id in the grib file
 *      leveltype - the leveltype id from table 3/3a of the grib document.
 *      level1    - First level of the data in units described by leveltype.
 *      level2    - Second level of the data in units described by leveltype.
 *      fcsttime1 - First forecast time in seconds.
 *      fcsttime2 - Second forecast time in seconds.
 *      indices   - an array of indices that match
 *      num_indices - the number of matches and output indices
 * 
 *    Note: If an input variable is set set to -INT_MAX, then any value
 *                will be considered a match.
 *    Return:
 *          The number of matching indices.
 *
 ***************************************************************************/

int rg_get_indices(GribInfo *gribinfo, FindGrib *findgrib, int indices[])
{
  int gribnum;
  int matchnum = 0;

  for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
    if (compare_record(gribinfo, findgrib, gribnum) == 1) {
      indices[matchnum] = gribnum;
      matchnum++;
    }    
  }
  return matchnum;
}

/*************************************************************************
 *
 * Returns an array of dates that correspond to particular input grib fields.
 * The dates will be sorted so that the dates increase as the index increases.
 *
 * Interface:
 *    Input:
 *         gribinfo - pointer to a previously allocated gribinfo structure.  
 *                    The gribinfo structure is filled in this function.
 *         usParm_id: an array of parameter identifiers that could be
 *                    used as a sea level pressure field (From table 2 of 
 *                    grib documentation)
 *         usLevel_id: the level id that could be used as a sea level pressure
 *                    field (from table 3 of the grib documentation)
 *         usHeight1: the height for the particular parameter and level
 *                    (in units described by the parameter index)
 *         numparms:  the number of parameters in each of the usParm_id, 
 *                    usLevel_id, and usHeight1 arrays.
 *    Output:
 *         dates:     the dates for which the input fields are available.
 *         
 *            Note: Values in the input arrays, usParm_id, usLevel_id, and 
 *                    usHeight with the same array index must correspond.
 *
 *    Return:
 *          The number of dates found.
 *************************************************************************/

int rg_get_dates(GribInfo *gribinfo,int usParm_id[],int usLevel_id[],
              int usHeight1[],int numparms,int dates[],int century[],
              int indices[])
{
  int datenum=0;
  int gribnum;
  int parmindex;
  int already_included;
  int i,j;
  int tmpval,tmpval2,tmpval3;

  /* Get the dates for the given parameters */

  for (parmindex = 0; parmindex < numparms; parmindex++) {
    for (gribnum = 0; gribnum < gribinfo->num_elements; gribnum++) {
      if ((gribinfo->elements[gribnum].usParm_id == usParm_id[parmindex]) &&
          (gribinfo->elements[gribnum].usLevel_id == usLevel_id[parmindex]) &&
          (gribinfo->elements[gribnum].usHeight1 == usHeight1[parmindex])) {
        already_included = 0;
        for (i = 0; i < datenum; i++){
          if ((dates[datenum] == gribinfo->elements[gribnum].date) && 
              (century[datenum] == gribinfo->elements[gribnum].century)) {
            already_included = 1;
            break;
          }
        }
        if (!already_included) {
          dates[datenum] = gribinfo->elements[gribnum].date;
          century[datenum] = gribinfo->elements[gribnum].century;
          indices[datenum] = gribnum;
          datenum++;
        }
      }
    }
  }

  /* Sort the dates into increasing order */
  for (j = 1; j < datenum; j++) {
    tmpval = dates[j];
    tmpval2 = indices[j];
    tmpval3 = century[j];
    for (i=j-1; i >= 0; i--) {
      if (dates[i] <= tmpval) break;
      dates[i+1] = dates[i];
      indices[i+1] = indices[i];
      century[i+1] = century[i];
    }
    dates[i+1] = tmpval;
    indices[i+1] = tmpval2;
    century[i+1] = tmpval3;
  }

  return datenum;
}

/****************************************************************************
 * This function returns the pds, gds, bms, and bms_head section of the
 * grib element
 *
 * Input:
 *   gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *              gribinfo structure is filled in this function.
 *   index - the index of the grib record to access as indexed by 
 *           setup_gribinfo
 *   
 * Output:
 *   *pds - a pointer to a structure holding the pds information
 *   *gds - a pointer to a structure holding the gds information
 *   *bms - a pointer to a structure holding the bms information
 *   *bds_head - a pointer to a structure holding the binary data section
 *             header information
 *
 ***************************************************************************
 */
int rg_get_grib_header(GribInfo *gribinfo, int index, PDS_INPUT *pds, 
                    grid_desc_sec *gds,BMS_INPUT *bms)
{
  int xsize,ysize,j;
  
  memcpy(pds,gribinfo->elements[index].pds,sizeof(PDS_INPUT));
  memcpy(gds,gribinfo->elements[index].gds,sizeof(grid_desc_sec));
  memcpy(bms,gribinfo->elements[index].bms,sizeof(BMS_INPUT));
  
  /* Reset the dimensions for thinned grids */
  if (gribinfo->elements[index].gds->head.thin != NULL) {
    if (gds->head.thin != NULL) {
      if ((gds->head.usData_type == LATLON_PRJ) || 
          (gds->head.usData_type == GAUSS_PRJ) ||
          (gds->head.usData_type == ROT_LATLON_PRJ) ||
          (gds->head.usData_type == ROT_GAUSS_PRJ) ||
          (gds->head.usData_type == STR_LATLON_PRJ) ||
          (gds->head.usData_type == STR_GAUSS_PRJ) ||
          (gds->head.usData_type == STR_ROT_LATLON_PRJ) ||
          (gds->head.usData_type == STR_ROT_GAUSS_PRJ)) {
        ysize = gds->llg.usNj;
      } else if (gds->head.usData_type == MERC_PRJ) {
        ysize = gds->merc.rows;
      } else if (gds->head.usData_type == POLAR_PRJ) {
        ysize = gds->pol.usNy;
      } else if ((gds->head.usData_type == LAMB_PRJ) ||
                 (gds->head.usData_type == ALBERS_PRJ) ||
                 (gds->head.usData_type == OBLIQ_LAMB_PRJ)) {
        ysize = gds->lam.iNy;
      }

      xsize = 0;
      for (j = 0; j<ysize; j++) {
        if (gds->head.thin[j] > xsize) {
          xsize = gds->head.thin[j];
        }
      }
      

      if ((gds->head.usData_type == LATLON_PRJ) || 
          (gds->head.usData_type == GAUSS_PRJ) ||
          (gds->head.usData_type == ROT_LATLON_PRJ) ||
          (gds->head.usData_type == ROT_GAUSS_PRJ) ||
          (gds->head.usData_type == STR_LATLON_PRJ) ||
          (gds->head.usData_type == STR_GAUSS_PRJ) ||
          (gds->head.usData_type == STR_ROT_LATLON_PRJ) ||
          (gds->head.usData_type == STR_ROT_GAUSS_PRJ)) {
        gds->llg.usNi = xsize;
        gds->llg.iDi = abs(gds->llg.lLat2 - gds->llg.lLat1)/(xsize-1);
      } else if (gds->head.usData_type == MERC_PRJ) {
        gds->merc.cols = xsize;
      } else if (gds->head.usData_type == POLAR_PRJ) {
        gds->pol.usNx = xsize;
      } else if ((gds->head.usData_type == LAMB_PRJ) ||
                 (gds->head.usData_type == ALBERS_PRJ) ||
                 (gds->head.usData_type == OBLIQ_LAMB_PRJ)) {
        gds->lam.iNx = xsize;
      }
      
    }
  }
  return 1;
}

/****************************************************************************
 * This returns the index of the gribdata for paramaters which match the input
 * parameters and for the date closest to the input targetdate.  If dates are
 * not found either within hours_before or hours_after the time, then a missing
 * value is returned.
 *
 * Interface:
 *   Input:
 *     gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                  gribinfo structure is filled in this function.
 *     targetdate:  This is the date which dates in the grib data will be 
 *                  compared to.  (format: integer yymmddhh)
 *     hours_before: The maximum difference in time prior to the targetdate
 *                  for which data should be searched for.
 *     hours_after: The maximum difference in time after the targetdate for
 *                  which data should be searched for.
 *     usParm_id:   an array of parameter identifiers that could be
 *                  used as a sea level pressure field (From table 2 of 
 *                  grib documentation)
 *     usLevel_id:  the level id that could be used as a sea level pressure
 *                  field (from table 3 of the grib documentation)
 *     usHeight1:   the height for the particular parameter and level
 *                  (in units described by the parameter index)
 *     numparms:    the number of parameters in each of the usParm_id, 
 *                  usLevel_id, and usHeight1 arrays.
 *   Return:
 *     the index of the gribdata with a time closest to the target date.
 *     -1 if there is no time within the input time limits.
 * 
 ****************************************************************************/
int rg_get_index_near_date(GribInfo *gribinfo,char targetdate[STRINGSIZE],
                        int century,int hours_before,int hours_after,
                        int usParm_id[],int usLevel_id[],int usHeight1[],
                        int numparms)
{
  int dates[500],indices[500],centuries[500];
  int date_before = MISSING;
  int date_after = MISSING;
  int century_before,century_after;
  int date_diff_before = MISSING;
  int date_diff_after = MISSING;
  int index_before,index_after;
  int numdates,datenum;
  int index;
  int itargetdate;

  itargetdate = atoi(targetdate);

  numdates = rg_get_dates(gribinfo,usParm_id,usLevel_id,usHeight1,numparms,
                          dates,centuries,indices);
  if (numdates <= 0) {
    fprintf(stderr,"get_index_near_date: No dates were found\n");
    return -1;
  }
  
  for (datenum = 0; datenum < numdates; datenum++) {
    if ((dates[datenum] > itargetdate) && (centuries[datenum] >= century)) {
      century_after = centuries[datenum];
      date_after = dates[datenum];
      index_after = indices[datenum];
      break;
    } else {
      century_before = centuries[datenum];
      date_before = dates[datenum];
      index_before = indices[datenum];
    }
  }
  
  if (date_after != MISSING)
    date_diff_after = date_diff(date_after,century_after,itargetdate,century);
  if (date_before != MISSING)
    date_diff_before = 
      date_diff(itargetdate,century,date_before,century_before);

  if ((date_after != MISSING) && (date_before != MISSING)) {
    if ((date_diff_after <= hours_after) && 
        (date_diff_before <= hours_before)) {
      if (date_diff_after < date_diff_before)
        index = index_before;
      else
        index = index_after;
    } else if (date_diff_after <= hours_after) {
        index = index_after;
    } else if (date_diff_before <= hours_before) {
        index = index_before;
    } else {
      index = -1;
    }
  } else if (date_after != MISSING) {
    if (date_diff_after <= hours_after)
      index = index_after;
    else 
      index = -1;
  } else if (date_before != MISSING) {
    if (date_diff_before <= hours_before)
      index = index_before;
    else
      index = -1;
  } else {
    index = -1;
  }

  return index;
  
}

/*****************************************************************************
 *
 * returns valid time ( = init time + forecast time)
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *    index    - index number of record to get valid time from
 *
 * Output:
 *    valid_time - yyyymmddhhmmss
 *
 * Return:
 *    0 for success
 *   -1 for error
 *    
 *****************************************************************************/
int rg_get_valid_time(GribInfo *gribinfo, int index, char valid_time[])
{
  strcpy(valid_time, gribinfo->elements[index].valid_time);
  return 0;
}

/*****************************************************************************
 *
 * returns generating center id
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *    index    - index number of record to get valid time from
 *
 * Return:
 *    generating center id
 *   -1 for error
 *    
 *****************************************************************************/
int rg_get_center_id(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].center_id;
}

/*****************************************************************************
 *
 * returns parameter table version number
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *    index    - index number of record to get valid time from
 *
 * Return:
 *    parameter table version number
 *   -1 for error
 *    
 *****************************************************************************/
int rg_get_parmtbl(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].parmtbl;
}

/*****************************************************************************
 *
 * returns generating process id
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *    index    - index number of record to get valid time from
 *
 * Return:
 *    generating process id
 *   -1 for error
 *    
 *****************************************************************************/
int rg_get_proc_id(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].proc_id;
}

/*****************************************************************************
 *
 * returns sub center id
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *    index    - index number of record to get valid time from
 *
 * Return:
 *    sub center id
 *   -1 for error
 *    
 *****************************************************************************/
int rg_get_subcenter_id(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].subcenter_id;
}

/**************************************************************************
 *
 * Interpolates grib grid data to a point location.
 *
 * Interface:
 *     input: 
 *       gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                  gribinfo structure is filled in this function.
 *       index    - the index of gribinfo for which data is to be retrieved.
 *                  the first grib record is number 1.
 *       column   - the column of the point in grid coordinates (can be 
 *                  floating point number).  leftmost column is 1.
 *       row      - the row of the point in grid coordinates (can be 
 *                  floating point number).  bottommost row is 1.
 *
 *    return:
 *       on success - the interpolated value at the column,row location.
 *       on failure - -99999
 * 
 ***************************************************************************/

float rg_get_point(GribInfo *gribinfo, int index, float column, float row)
{
  int status;
  GRIB_PROJECTION_INFO_DEF Proj;
  BDS_HEAD_INPUT bds_head;
  int dummy;
  float **grib_out;
  float y1, y2;
  int numrows, numcols;
  int top, left, right, bottom;
  float outval;
  
  numrows = rg_get_numrows(gribinfo, index);
  numcols = rg_get_numcols(gribinfo, index);

  grib_out = (float **)alloc_float_2d(numrows,numcols);
  if (grib_out == NULL) {
    fprintf(stderr,"rg_get_point: Could not allocate space for grib_out\n");
    return -99999;
  }

  status = rg_get_data(gribinfo, index, grib_out);
  if (status < 0) {
    fprintf(stderr,"rg_get_point: rg_get_data failed\n");
    return -99999;
  }

  /* Do the interpolation here */
  bottom = floor(row);
  top = floor(row+1);
  left = floor(column);
  right = floor(column+1);
  
  y1 = (row - bottom) * (grib_out[top][left] - grib_out[bottom][left]) + 
    grib_out[bottom][left];
  y2 = (row - bottom) * (grib_out[top][right] - grib_out[bottom][right]) + 
    grib_out[bottom][right];
  outval = (y2 - y1) * (column - left) + y1;

  free_float_2d(grib_out,numrows,numcols);

  return outval;
  
}

/**************************************************************************
 *
 * Interpolates grib grid data to a point location.
 *
 * Interface:
 *     input: 
 *       gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                  gribinfo structure is filled in this function.
 *       index    - the index of gribinfo for which data is to be retrieved.
 *                  the first grib record is number 1.
 *     input and output:
 *       pointdata- array of pointdata structures.  Only the column and
 *                  row values in the structures need to be filled.  On
 *                  output, the 'value' member of pointdata is filled.
 *     input:
 *       numpoints- number of pointdata structures in the array.
 *
 *    return:
 *       on success - the interpolated value at the column,row location.
 *       on failure - -99999
 * 
 ***************************************************************************/
int rg_get_points(GribInfo *gribinfo, int index, PointData pointdata[], 
                   int numpoints)
{
  int status;
  float **grib_out;
  float y1, y2;
  int numrows, numcols;
  int top, left, right, bottom;
  float column, row;
  int idx;

  numrows = rg_get_numrows(gribinfo, index);
  numcols = rg_get_numcols(gribinfo, index);

  grib_out = (float **)alloc_float_2d(numrows,numcols);
  if (grib_out == NULL) {
    fprintf(stderr,"rg_get_points: Could not allocate space for grib_out\n");
    return -99999;
  }

  status = rg_get_data(gribinfo, index, grib_out);
  if (status < 0) {
    fprintf(stderr,"rg_get_points: rg_get_data failed\n");
    return -99999;
  }

  for (idx = 0; idx < numpoints; idx++) {

    /* Change from 1 based to 0 based col/row */
    row = pointdata[idx].row;
    column = pointdata[idx].column;
    
    /* Do the interpolation here */
    bottom = floor(row);
    top = floor(row+1);
    left = floor(column);
    right = floor(column+1);
    
    y1 = (row - bottom) * (grib_out[top][left] - grib_out[bottom][left]) + 
      grib_out[bottom][left];
    y2 = (row - bottom) * (grib_out[top][right] - grib_out[bottom][right]) + 
      grib_out[bottom][right];
    pointdata[idx].value = (y2 - y1) * (column - left) + y1;

  }

  free_float_2d(grib_out,numrows,numcols);

  return 1;
}

/**************************************************************************
 *
 * Remove an element from an array and decrease, by one, indices of all 
 * elements with an index greater than the index of the element to remove.
 *
 * Interface:
 *     input: 
 *        array - the integer array to manipulate
 *        index - the index of the element to remove
 *        size - the number of elements in the array
 *
 ***************************************************************************/
void remove_element(int array[],int index, int size)
{
  int j;
  
  for (j = index; j < size-1; j++) {
    array[j] = array[j+1];
  }

}

/****************************************************************************
 * Advance the time by the input amount
 *
 * Interface:
 *     Input: 
 *         century - an integer value for the century (20 for 1900's)
 *                   If the century is advanced, this value is advanced
 *                   and output to the calling routine.
 *         year - a 2 digit value for the year.
 *         month - a 2 digit value for the month.
 *         day - the day of the month
 *         hour - the hour of the day
 *         amount - the amount to advance the time by.
 *         unit - the units for the amount.  These are values from table 4
 *                of the grib manual.
 *    return:
 *         a date in the form yymmddhh
 ****************************************************************************/

int advance_time(int *century, int year, int month, int day, int hour,
                 int amount, int unit)
{
  int daysinmonth[] = {31,28,31,30,31,30,31,31,30,31,30,31};
  int date;

  switch(unit) {
  case 0:
    hour += (int)((amount/60.)+0.5);
    break;
  case 1:
    hour += amount;
    break;
  case 2:
    day += amount;
    break;
  case 3:
    month += amount;
    break;
  case 4:
    year += amount;
    break;
  case 5:
    year += 10*amount;
    break;
  case 6:
    year += 30*amount;
    break;
  case 7:
    year += 100*amount;
    break;
  case 10:
    hour += 3*amount;
    break;
  case 11:
    hour += 6*amount;
    break;
  case 12:
    hour += 12*amount;
    break;
  case 50:
    hour += (int)((amount/12.)+0.5);
  case 254:
    hour += (int)((amount/(60.*60.))+0.5);
    break;
  default:
    fprintf(stderr,"WARNING: Could not advance time, incorrect unit: %d\n",
            unit);
    return -1;
  }

  while (hour >= 24) {
    day++;
    hour -= 24;
  }
  while (month > 12) {
    year++;
    month -= 12;
  }

  /* if it is a leap year, change days in month for Feb now. */ 
  if (isLeapYear(year)) daysinmonth[1] = 29;

  while (day > daysinmonth[month-1]) {
    day -= daysinmonth[month-1];
    month++;
    if (month > 12) {
      year++;
      month -= 12;
      if (isLeapYear(year))
        daysinmonth[1] = 29;
      else
        daysinmonth[1] = 28;
    }
  }

  if (year > 100) {
    (*century)++;
  }

  if (year >= 100) {
    year -= 100;
  }

  date = hour*1 + day*100 + month*10000 + year*1000000;

  return date;

}
/****************************************************************************
 * Advance the time by the input amount
 *
 * Interface:
 *     Input: 
 *         startdate  - initialization date in the form yyyymmdd[HHMMSS].  If any
 *                        part of HHMMSS is not specified, it will be set to 0.
 *         amount     - the amount (in seconds) to advance the time by.
 *
 *    Output:
 *         enddate[] - the time advanced to: yyyymmddHHMMSS format. 
 *
 *    Return:
 *          1 - success
 *         -1 - failure
 *
 ****************************************************************************/
char *advance_time_str(char startdatein[], int amount, char enddate[])
{
  struct tm starttp;
  struct tm endtp;
  char startdate[15];
  time_t time;

  strcpy(startdate,startdatein);
  while (strlen(startdate) < 14) {
    strcpy(startdate+(strlen(startdate)),"0");
  }

  /* This forces all calculations to use GMT time */
  putenv("TZ=GMT0");
  tzset();

  sscanf(startdate,"%4d%2d%2d%2d%2d%2d",&(starttp.tm_year),&(starttp.tm_mon),
         &(starttp.tm_mday),&(starttp.tm_hour),&(starttp.tm_min),
         &(starttp.tm_sec));
  starttp.tm_mon -= 1;
  starttp.tm_year -= 1900;
  time = mktime(&starttp);
  time += amount;
  localtime_r(&time, &endtp);
  strftime(enddate,15,"%Y%m%d%H%M%S",&endtp);
  
  return enddate;
}

/****************************************************************************
 * Returns the difference in time in hours between date1 and date2 
 * (date1-date2). 
 * 
 * Interface:
 *   Input:
 *     date1,date2:  dates in yymmddhh format (integers)
 *     century1,century2: centuries for each date (20 for 1900's).
 *   Return:
 *     the difference in time between the first and second dates in hours.
 ****************************************************************************/
int date_diff(int date1,int century1,int date2,int century2)
{
  return (hours_since_1900(date1,century1) - 
           hours_since_1900(date2,century2));
}

/****************************************************************************
 * Returns the number of hours since Jan 1, at 00:00 1900.
 * 
 * Interface:
 *   Input:
 *     date:    integer in form yymmddhh
 *     century: 2 digit century (20 for 1900's)
 *   Return:
 *     the number of hours since 00:00 Jan1, 1900.
 *
 ****************************************************************************/
int hours_since_1900(int date,int century)
{
  int daysinmonth[] = {31,28,31,30,31,30,31,31,30,31,30,31};
  int hour,day,month,year;
  int days_since_1900 = 0;
  int i;

  hour = date%100;
  day = (date%10000)/100;
  month = (date%1000000)/10000;
  year = (date%100000000)/1000000;

  days_since_1900 += day;
  
  if (isLeapYear((century-1)*100 + year))
    daysinmonth[1] = 29;
  else 
    daysinmonth[1] = 28;
  
  for (i = 0; i < (month - 1); i++)
    days_since_1900 += daysinmonth[i];

  for (i=0; i < (year + ((century - 20)*100) - 1); i++) {
    if (isLeapYear((century - 1)*100 + year))
      days_since_1900 += 366;
    else
      days_since_1900 += 365;
  }

  return days_since_1900*24 + hour;

}

/****************************************************************************
 *
 * Returns true if the input year is a leap year, otherwise returns false
 *
 ****************************************************************************/
int isLeapYear(int year)
{
  if ( (((year % 4) == 0) && ((year % 100) != 0)) 
       || ((year % 400) == 0) ) 
    return 1;
  else
    return 0;
                                                  
}

/*****************************************************************************
 *
 * Returns the number of grib elements (gribinfo->num_elements) processsed
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *   the number of elements in the gribinfo structure
 ****************************************************************************/

int rg_num_elements(GribInfo *gribinfo){

  return gribinfo->num_elements;

}

/*****************************************************************************
 *
 * Deallocates the elements in the gribinfo structure and closes the files.
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 *****************************************************************************/
void rg_free_gribinfo_elements(GribInfo *gribinfo)
{
  int i;
  
  for (i=0; i<gribinfo->num_elements; i++) {
    free(gribinfo->elements[i].pds);
    free(gribinfo->elements[i].gds);
    free(gribinfo->elements[i].bms);
    free(gribinfo->elements[i].bds_head);
    fclose(gribinfo->elements[i].fp);
  }
  free(gribinfo->elements);
}

/*****************************************************************************
 *
 * Returns the value for level1 (gribinfo->usHeight1)
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *   value for level1
 ****************************************************************************/

int rg_get_level1(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].usHeight1;
}

/*****************************************************************************
 *
 * Returns the value for level2 (gribinfo->usHeight2)
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *   value for level1
 ****************************************************************************/

int rg_get_level2(GribInfo *gribinfo, int index)
{
  return gribinfo->elements[index].usHeight2;
}

/*****************************************************************************
 *
 * returns number of rows in grid
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    number of rows in grid
 *****************************************************************************/
int rg_get_numrows(GribInfo *gribinfo,int index)
{
  if ((gribinfo->elements[index].gds->head.usData_type == LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == ROT_LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == ROT_GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_ROT_LATLON_PRJ)|| 
      (gribinfo->elements[index].gds->head.usData_type == STR_ROT_GAUSS_PRJ)) 
    {
    return gribinfo->elements[index].gds->llg.usNj;
  } else if (gribinfo->elements[index].gds->head.usData_type == MERC_PRJ) {
    return gribinfo->elements[index].gds->merc.rows;
  } else if (gribinfo->elements[index].gds->head.usData_type == LAMB_PRJ) {
    return gribinfo->elements[index].gds->lam.iNy;
  } else if (gribinfo->elements[index].gds->head.usData_type == POLAR_PRJ) {
    return gribinfo->elements[index].gds->pol.usNy;
  }

}
/*****************************************************************************
 *
 * returns number of columns in grid
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    number of columns in grid
 *****************************************************************************/
int rg_get_numcols(GribInfo *gribinfo,int index)
{
  if ((gribinfo->elements[index].gds->head.usData_type == LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == ROT_LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == ROT_GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_LATLON_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_GAUSS_PRJ) || 
      (gribinfo->elements[index].gds->head.usData_type == STR_ROT_LATLON_PRJ)|| 
      (gribinfo->elements[index].gds->head.usData_type == STR_ROT_GAUSS_PRJ)) 
    {
      return gribinfo->elements[index].gds->llg.usNi;
  } else if (gribinfo->elements[index].gds->head.usData_type == MERC_PRJ) {
    return gribinfo->elements[index].gds->merc.cols;
  } else if (gribinfo->elements[index].gds->head.usData_type == LAMB_PRJ) {
    return gribinfo->elements[index].gds->lam.iNx;
  } else if (gribinfo->elements[index].gds->head.usData_type == POLAR_PRJ) {
    return gribinfo->elements[index].gds->pol.usNx;
  }

}
/*****************************************************************************
 *
 * returns the offset (in bytes) from the beginning of the file.
 *
 * Input:
 *    gribinfo - pointer to a filled gribinfo structure.
 *
 * Return:
 *    offset (in bytes) from beginning of file
 *****************************************************************************/
int rg_get_offset(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].offset;
}
/*****************************************************************************
 *
 * returns the grib record ending position (in bytes) from the beginning of 
 *    the file.
 *
 * Input:
 *    gribinfo - pointer to a filled gribinfo structure.
 *
 * Return:
 *    position (in bytes) of the end of the grib record within the file.
 *****************************************************************************/
int rg_get_end(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].end;
}
/*****************************************************************************
 *
 * returns grib id of input grid
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    grib id of input grid
 *****************************************************************************/
int rg_get_gridnum(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].pds->usGrid_id;
}
/*****************************************************************************
 *
 * returns date
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    date (yymmddhh) in integer type
 *****************************************************************************/
int rg_get_date(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].date;
}
/*****************************************************************************
 *
 * returns century
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    century in integer type
 *****************************************************************************/
int rg_get_century(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].century;
}
/*****************************************************************************
 *
 * returns forecast time
 *
 * Input:
 *    gribinfo - pointer to a previously allocated gribinfo structure.  The 
 *                 gribinfo structure is filled in this function.
 *
 * Return:
 *    forecast time in units described by usFcst_unit_id
 *****************************************************************************/
int rg_get_forecast_time(GribInfo *gribinfo,int index)
{
  return gribinfo->elements[index].pds->usP1;
}

/*****************************************************************************
 *
 * reads the gribmap file, and stores the information in the GribParameters
 *    structure.
 *
 * Input:
 *   gribmap - pointer to a previously allocated GribParameters structure.  
 *              The gribmap structure is filled in this function.
 *   file - the name of the gribmap file to read.
 *
 *   Return:
 *      1  - successful call to setup_gribinfo
 *     -1 - call to setup_gribinfo failed
 *
 *****************************************************************************/
int rg_setup_gribmap(GribParameters *gribmap, char filename[])
{
  FILE *fid;
  char line[500];
  int id, center, subcenter, table;
  int idx;
  
  fid = fopen(filename,"r");
  if (fid == NULL)
    {
      fprintf(stderr,"Could not open %s\n",filename);
      return -1;
    }

  gribmap->parms = (GribTableEntry *)malloc(sizeof(GribTableEntry));

  idx = 0;
  while (fgets(line,500,fid) != NULL)
    {
      /* Skip over comments at begining of gribmap file */
      if (line[0] == '#') continue;

      sscanf(line,"%d:",&id);
      if (id == -1) 
        {
          sscanf(line,"%d:%d:%d:%d",&id,&center,&subcenter,&table);
        } 
      else
        {
          gribmap->parms = 
            (GribTableEntry *)realloc(gribmap->parms,
                                        (idx+1)*sizeof(GribTableEntry));
          gribmap->parms[idx].center = center;
          gribmap->parms[idx].subcenter = subcenter;
          gribmap->parms[idx].table = table;
          sscanf(line,"%d:%[^:]:%[^:]",&(gribmap->parms[idx].parmid),
                 gribmap->parms[idx].name,
                 gribmap->parms[idx].comment);
          idx++;
        }
    }


  gribmap->num_entries = idx;

  close(fid);
  return 1;
}

/*****************************************************************************
 *
 * finds the gribmap entry described by  the gribmap file, and stores the information in the GribParameters
 *    structure.
 *
 * Input:
 *   gribmap - pointer to structure that was filled by a call to 
 *             rg_setup_gribmap
 *   table   - if set to -1, the first table the valid name will be used.
 *             Otherwise, the table id must match as well.
 *   name    - name of the parameter to find.
 * Output
 *   gribmap_parms - pointer to GribTableEntry structure containing 
 *             information about the parameter that was found.
 *
 *   Return:
 *      1  - successful call to setup_gribinfo
 *     -1 - call to setup_gribinfo failed
 *
 *****************************************************************************/
int rg_gribmap_parameter(GribParameters *gribmap, char name[], int table,
                         GribTableEntry *gribmap_parms)
{
  int idx;
  int found;

  found = 0;
  for (idx = 0; idx < gribmap->num_entries; idx++) 
    {
      
      if (strcmp(gribmap->parms[idx].name,name) == 0)
        {
          if ((table == -1) || (table == gribmap->parms[idx].table))
            {
              /* We found a match! */
              found = 1;
              break;
            }
        }
    }
  
  if (found) 
    {
      memcpy(gribmap_parms,&(gribmap->parms[idx]),sizeof(GribTableEntry));
      return 1;
    }
  else
    {
      return -1;
    }
}

/*****************************************************************************
 *
 * Deallocates the elements in the gribmap structure.
 *
 * Input:
 *    gribmap - pointer to a previously allocated gribmap structure.  The 
 *                 gribmap structure is filled in this function.
 *
 *****************************************************************************/
void rg_free_gribmap_elements(GribParameters *gribmap)
{
  free(gribmap->parms);
}

/*****************************************************************************
 *
 * Compares the elements in a findgrib structure with the elements in the
 *   gribinfo structure for the input index.  If they match, returns 1, 
 *   otherwise, returns 0.
 *
 * Input:
 *    gribinfo
 *    findgrib
 *    index    - the index of the grib record in gribinfo to compare to.
 *
 *****************************************************************************/
int compare_record(GribInfo *gribinfo, FindGrib *findgrib, int gribnum)
{

  /* 
   *  Note (6/20/05): This searching is very inefficient. We may need to 
   *   improve this, since, for WRF, when searching through boundary data, 
   *   each search is slower that the previous, since the record to be 
   *   found turns out to be farther into the list.
   */
  
  int retval = 0;

  if ((findgrib->center_id == -INT_MAX) || 
      findgrib->center_id == gribinfo->elements[gribnum].center_id) {
    if ((findgrib->subcenter_id == -INT_MAX) ||
        findgrib->subcenter_id == gribinfo->elements[gribnum].subcenter_id) {
      if ((findgrib->parmtbl_version == -INT_MAX) || 
          findgrib->parmtbl_version == gribinfo->elements[gribnum].parmtbl) {
        if ((strcmp(findgrib->initdate,"*") == 0) || 
            (strncmp(gribinfo->elements[gribnum].initdate,findgrib->initdate,
                     strlen(findgrib->initdate)) == 0)) {
          if ((strcmp(findgrib->validdate,"*") == 0) || 
              (strncmp(gribinfo->elements[gribnum].valid_time,
                       findgrib->validdate,
                       strlen(findgrib->validdate)) == 0)) {
            if ((findgrib->parmid == -INT_MAX) || 
                (findgrib->parmid == 
                 gribinfo->elements[gribnum].usParm_id)) {
              if ((findgrib->leveltype == -INT_MAX) || 
                  (findgrib->leveltype == 
                   gribinfo->elements[gribnum].usLevel_id)) {
                if ((findgrib->level1 == -INT_MAX) || 
                    (findgrib->level1 == 
                     gribinfo->elements[gribnum].usHeight1)) {
                  if ((findgrib->level2 == -INT_MAX) || 
                      (findgrib->level2 == 
                       gribinfo->elements[gribnum].usHeight2)) {
                    if ((findgrib->fcsttime1 == -INT_MAX) || 
                        (findgrib->fcsttime1 == 
                         gribinfo->elements[gribnum].fcsttime1)) {
                      if ((findgrib->fcsttime2 == -INT_MAX) || 
                          (findgrib->fcsttime2 == 
                           gribinfo->elements[gribnum].fcsttime2)) {
                        retval = 1;
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }

  return retval;
}


/*****************************************************************************
 *
 * returns the multiplication factor to convert grib forecast times to 
 *   seconds.
 *
 * Input: 
 *    unit_id  - grib forecast unit id, from Table 4.
 *
 * Return:
 *    conversion factor
 *****************************************************************************/
int get_factor2(int unit)
{
  int factor;
  
  switch (unit) {
  case 0:
    factor = 60;
    break;
  case 1:
    factor = 60*60;
    break;
  case 2:
    factor = 60*60*24;
    break;
  case 10:
    factor = 60*60*3;
    break; 
  case 11:
    factor = 60*60*3;
    break;
  case 12:
    factor = 60*60*12;
    break;
  case 50:
    /* This is a WSI (non-standard) time unit of 5 minutes */
    factor = 5*60;
    break;
  case 254:
    factor = 1;
    break;
  default:
    fprintf(stderr,"Invalid unit for forecast time: %d\n",unit);
    factor = 0;
  }
  return factor;
}