/wrfv2_fire/external/io_grib1/grib1_util/read_grib.c

Large files files are truncated, but you can click here to view the full file

 /**************************************************************************
 * 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 fil…