stream_cgribex.c

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <limits.h>
#include <stdio.h>

#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
#include "file.h"
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
#include "stream_grb.h"
#include "stream_cgribex.h"

#if  defined  (HAVE_LIBCGRIBEX)
#  include "cgribex.h"
#endif

typedef struct {
  int param;
  int level1;
  int level2;
  int ltype;
  int tsteptype;
} compvar_t;


#if  defined  (HAVE_LIBCGRIBEX)
static
int cgribexGetGridType(int *isec2)
{
  int gridtype = GRID_GENERIC;

  switch (ISEC2_GridType)
    {
    case  GRIB1_GTYPE_LATLON:     { if ( ISEC2_Reduced )      break; }
    case  GRIB1_GTYPE_LATLON_ROT: { gridtype = GRID_LONLAT;   break; }
    case  GRIB1_GTYPE_LCC:        { gridtype = GRID_LCC;      break; }
    case  GRIB1_GTYPE_GAUSSIAN:   { if ( ISEC2_Reduced )
	                              gridtype = GRID_GAUSSIAN_REDUCED;
                         	    else
				      gridtype = GRID_GAUSSIAN;
          	                    break;
                                  }
    case  GRIB1_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL; break; }
    case  GRIB1_GTYPE_GME:        { gridtype = GRID_GME;      break; }
    }

  return gridtype;
}

static
bool cgribexGetIsRotated(int *isec2)
{
  bool isRotated = (ISEC2_GridType == GRIB1_GTYPE_LATLON_ROT) ? true : false;

  return isRotated;
}

static
int cgribexGetZaxisHasBounds(int grb_ltype)
{
  int lbounds = 0;

  switch (grb_ltype)
    {
    case GRIB1_LTYPE_SIGMA_LAYER:
    case GRIB1_LTYPE_HYBRID_LAYER:
    case GRIB1_LTYPE_LANDDEPTH_LAYER:
      {
	lbounds = 1;
	break;
      }
    }

  return lbounds;
}

static
int cgribexGetTimeUnit(int *isec1)
{
  int timeunit = TUNIT_HOUR;
  static bool lprint = true;

  switch ( ISEC1_TimeUnit )
    {
    case ISEC1_TABLE4_MINUTE:    timeunit = TUNIT_MINUTE;    break;
    case ISEC1_TABLE4_QUARTER:   timeunit = TUNIT_QUARTER;   break;
    case ISEC1_TABLE4_30MINUTES: timeunit = TUNIT_30MINUTES; break;
    case ISEC1_TABLE4_HOUR:      timeunit = TUNIT_HOUR;      break;
    case ISEC1_TABLE4_3HOURS:    timeunit = TUNIT_3HOURS;    break;
    case ISEC1_TABLE4_6HOURS:    timeunit = TUNIT_6HOURS;    break;
    case ISEC1_TABLE4_12HOURS:   timeunit = TUNIT_12HOURS;   break;
    case ISEC1_TABLE4_DAY:       timeunit = TUNIT_DAY;       break;
    default:
      if ( lprint )
	{
	  Message("GRIB time unit %d unsupported!", ISEC1_TimeUnit);
	  lprint = false;
	}
      break;
    }

  return timeunit;
}

static
bool cgribexTimeIsFC(int *isec1)
{
  bool isFC = (ISEC1_TimeRange == 10 && ISEC1_TimePeriod1 == 0 && ISEC1_TimePeriod2 == 0) ? false : true;

  return isFC;
}

static
int cgribexGetTsteptype(int timerange)
{
  int tsteptype = TSTEP_INSTANT;
  static bool lprint = true;

  switch ( timerange )
    {
    case  0:  tsteptype = TSTEP_INSTANT;  break;
    case  1:  tsteptype = TSTEP_INSTANT2; break;
    case  2:  tsteptype = TSTEP_RANGE;    break;
    case  3:  tsteptype = TSTEP_AVG;      break;
    case  4:  tsteptype = TSTEP_ACCUM;    break;
    case  5:  tsteptype = TSTEP_DIFF;     break;
    case 10:  tsteptype = TSTEP_INSTANT3; break;
    default:
      if ( lprint )
	{
	  Message("Time range indicator %d unsupported, set to 0!", timerange);
	  lprint = false;
	}
      break;
    }

  return tsteptype;
}

static
void cgribexGetGrid(stream_t *streamptr, int *isec2, double *fsec2, int *isec4, grid_t *grid, int iret)
{
  bool compyinc = true;
  int gridtype = cgribexGetGridType(isec2);

  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED && iret != -801 )
    {
      int ilat, nlon = 0;
      for ( ilat = 0; ilat < ISEC2_NumLat; ++ilat )
        if ( ISEC2_RowLon(ilat) > nlon ) nlon = ISEC2_RowLon(ilat);
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = nlon;
      ISEC4_NumValues = nlon*ISEC2_NumLat;
      compyinc = false;
    }

  grid_init(grid);
  cdiGridTypeInit(grid, gridtype, 0);
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
	  Error("numberOfPoints (%d) and gridSize (%d) differ!", ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
	grid->size  = ISEC4_NumValues;
	grid->x.size = ISEC2_NumLon;
	grid->y.size = ISEC2_NumLat;
        if ( gridtype == GRID_GAUSSIAN ) grid->np = ISEC2_NumPar;
	grid->xinc  = 0;
	grid->yinc  = 0;
	grid->xdef  = 0;
	/* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
	  {
	    if ( grid->x.size > 1 )
	      {
                bool recompinc = true;

                if ( ISEC2_LastLon < ISEC2_FirstLon && ISEC2_LastLon < 0 ) ISEC2_LastLon += 360000;

		if ( ISEC2_ResFlag && ISEC2_LonIncr > 0 )
                  {
                    if ( abs(ISEC2_LastLon - (ISEC2_FirstLon+ISEC2_LonIncr*(grid->x.size-1))) <= 2 )
                      {
                        recompinc = false;
                        grid->xinc = ISEC2_LonIncr * 0.001;
                      }
                  }

		/* recompute xinc if necessary */
                if ( recompinc ) grid->xinc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid->x.size-1);

		/* correct xinc if necessary */
		if ( ISEC2_FirstLon == 0 && ISEC2_LastLon > 354000 && ISEC2_LastLon < 360000 )
		  {
		    double xinc = 360. / grid->x.size;

		    if ( fabs(grid->xinc-xinc) > 0.0 )
		      {
			grid->xinc = xinc;
			if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
		      }
		  }
	      }
	    grid->xfirst = ISEC2_FirstLon * 0.001;
	    grid->xlast  = ISEC2_LastLon  * 0.001;
	    grid->xdef   = 2;
	  }
	grid->ydef  = 0;
	/* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
	  {
	    if ( grid->y.size > 1 && compyinc )
	      {
                bool recompinc = true;
		if ( ISEC2_ResFlag && ISEC2_LatIncr > 0 )
                  {
                    if ( abs(ISEC2_LastLat - (ISEC2_FirstLat+ISEC2_LatIncr*(grid->y.size-1))) <= 2 )
                      {
                        recompinc = false;
                        grid->yinc = ISEC2_LatIncr * 0.001;
                      }
                  }

		/* recompute yinc if necessary */
                if ( recompinc ) grid->yinc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid->y.size - 1);
	      }
	    grid->yfirst = ISEC2_FirstLat * 0.001;
	    grid->ylast  = ISEC2_LastLat  * 0.001;
	    grid->ydef   = 2;
	  }
	break;
      }
    case GRID_GAUSSIAN_REDUCED:
      {
        grid->np      = ISEC2_NumPar;
	grid->size    = ISEC4_NumValues;
        grid->rowlon  = ISEC2_RowLonPtr;
	grid->nrowlon = ISEC2_NumLat;
	grid->y.size   = ISEC2_NumLat;
	grid->xinc    = 0;
	grid->yinc    = 0;
	grid->xdef    = 0;
	/* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
	  {
	    if ( grid->x.size > 1 )
	      {
                if ( ISEC2_LastLon < ISEC2_FirstLon && ISEC2_LastLon < 0 ) ISEC2_LastLon += 360000;

		if ( ISEC2_ResFlag && ISEC2_LonIncr > 0 )
		  grid->xinc = ISEC2_LonIncr * 0.001;
		else
		  grid->xinc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid->x.size - 1);
	      }
	    grid->xfirst = ISEC2_FirstLon * 0.001;
	    grid->xlast  = ISEC2_LastLon  * 0.001;
	    grid->xdef   = 2;
	  }
	grid->ydef  = 0;
	/* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
	  {
	    if ( grid->y.size > 1 )
	      {
		if ( ISEC2_ResFlag && ISEC2_LatIncr > 0 )
		  grid->yinc = ISEC2_LatIncr * 0.001;
		else
		  grid->yinc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid->y.size - 1);
	      }
	    grid->yfirst = ISEC2_FirstLat * 0.001;
	    grid->ylast  = ISEC2_LastLat  * 0.001;
	    grid->ydef   = 2;
	  }
	break;
      }
    case GRID_LCC:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);

	grid->size  = ISEC4_NumValues;
	grid->x.size = ISEC2_NumLon;
	grid->y.size = ISEC2_NumLat;

	grid->lcc.xinc      = ISEC2_Lambert_dx;
	grid->lcc.yinc      = ISEC2_Lambert_dy;
	grid->lcc.originLon = ISEC2_FirstLon * 0.001;
	grid->lcc.originLat = ISEC2_FirstLat * 0.001;
	grid->lcc.lonParY   = ISEC2_Lambert_Lov * 0.001;
	grid->lcc.lat1      = ISEC2_Lambert_LatS1 * 0.001;
	grid->lcc.lat2      = ISEC2_Lambert_LatS2 * 0.001;
	grid->lcc.projflag  = ISEC2_Lambert_ProjFlag;
	grid->lcc.scanflag  = ISEC2_ScanFlag;

	grid->xdef   = 0;
	grid->ydef   = 0;

	break;
      }
    case GRID_SPECTRAL:
      {
	grid->size  = ISEC4_NumValues;
	grid->trunc = ISEC2_PentaJ;
	if ( ISEC2_RepMode == 2 )
	  grid->lcomplex = 1;
	else
	  grid->lcomplex = 0;

	break;
      }
    case GRID_GME:
      {
	grid->size  = ISEC4_NumValues;
	grid->gme.nd    = ISEC2_GME_ND;
	grid->gme.ni    = ISEC2_GME_NI;
	grid->gme.ni2   = ISEC2_GME_NI2;
	grid->gme.ni3   = ISEC2_GME_NI3;
	break;
      }
    case GRID_GENERIC:
      {
	grid->size  = ISEC4_NumValues;
	grid->x.size = 0;
	grid->y.size = 0;
	break;
      }
    default:
      {
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
	break;
      }
    }

  grid->isRotated = FALSE;
  if ( cgribexGetIsRotated(isec2) )
    {
      grid->isRotated = TRUE;
      grid->ypole     = - ISEC2_LatSP*0.001;
      grid->xpole     =   ISEC2_LonSP*0.001 - 180;
      grid->angle     = - FSEC2_RotAngle;
    }

  grid->type  = gridtype;
}

static
void cgribexAddRecord(stream_t * streamptr, int param, int *isec1, int *isec2, double *fsec2, double *fsec3,
		      int *isec4, long recsize, off_t position, int datatype, int comptype, int lmv, int iret)
{
  int varID;
  int levelID = 0;

  int vlistID = streamptr->vlistID;
  int tsID    = streamptr->curTsID;
  int recID   = recordNewEntry(streamptr, tsID);
  record_t *record = &streamptr->tsteps[tsID].records[recID];

  int tsteptype = cgribexGetTsteptype(ISEC1_TimeRange);
  int numavg    = ISEC1_AvgNum;

  int level1  = ISEC1_Level1;
  int level2  = ISEC1_Level2;

  /* fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, ISEC1_LevelType); */

  record->size      = (size_t)recsize;
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
  record->ltype     = ISEC1_LevelType;
  record->tsteptype = (short)tsteptype;

  grid_t *gridptr = (grid_t*) Malloc(sizeof(*gridptr));
  cgribexGetGrid(streamptr, isec2, fsec2, isec4, gridptr, iret);

  struct addIffNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, gridptr, 0);
  int gridID = gridAdded.Id;
  if ( gridAdded.isNew )
    {
      if ( gridptr->nrowlon )
        {
          size_t nrowlon = (size_t) gridptr->nrowlon;
          int *rowlon = gridptr->rowlon;
          gridptr->rowlon = (int*) Malloc(nrowlon * sizeof(int));
          memcpy(gridptr->rowlon, rowlon, nrowlon * sizeof(int));
        }
    }
  else
    Free(gridptr);

  int zaxistype = grib1ltypeToZaxisType(ISEC1_LevelType);

  if ( zaxistype == ZAXIS_HYBRID || zaxistype == ZAXIS_HYBRID_HALF )
    {
      size_t vctsize = (size_t)ISEC2_NumVCP;
      double *vctptr = &fsec2[10];

      varDefVCT(vctsize, vctptr);
    }

  int lbounds = cgribexGetZaxisHasBounds(ISEC1_LevelType);

  if ( datatype > 32 ) datatype = DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = DATATYPE_PACK;

  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, 0, 0,
	       datatype, &varID, &levelID, tsteptype, numavg, ISEC1_LevelType, -1,
               NULL, NULL, NULL, NULL, NULL, NULL);

  record->varID   = (short)varID;
  record->levelID = (short)levelID;

  varDefCompType(varID, comptype);

  if ( ISEC1_LocalFLag )
    {
      if      ( ISEC1_CenterID == 78  && isec1[36] == 253 ) // DWD local extension
        varDefEnsembleInfo(varID, isec1[54], isec1[53], isec1[52]);
      else if ( ISEC1_CenterID == 252 && isec1[36] ==   1 ) // MPIM local extension
        varDefEnsembleInfo(varID, isec1[38], isec1[39], isec1[37]);
    }

  if ( lmv ) varDefMissval(varID, FSEC3_MissVal);

  if ( varInqInst(varID) == CDI_UNDEFID )
    {
      int center, subcenter, instID;
      center    = ISEC1_CenterID;
      subcenter = ISEC1_SubCenterID;
      instID    = institutInq(center, subcenter, NULL, NULL);
      if ( instID == CDI_UNDEFID )
	instID = institutDef(center, subcenter, NULL, NULL);
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
      modelID = modelInq(varInqInst(varID), ISEC1_ModelID, NULL);
      if ( modelID == CDI_UNDEFID )
	modelID = modelDef(varInqInst(varID), ISEC1_ModelID, NULL);
      varDefModel(varID, modelID);
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
      int tableID;

      tableID = tableInq(varInqModel(varID), ISEC1_CodeTable, NULL);

      if ( tableID == CDI_UNDEFID )
	tableID = tableDef(varInqModel(varID), ISEC1_CodeTable, NULL);
      varDefTable(varID, tableID);
    }

  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
}

static
void MCH_get_undef(int *isec1, double *undef_pds, double *undef_eps)
{
  /* 2010-01-13: Oliver Fuhrer */
  if ( ISEC1_CenterID == 215 ) {
    if (isec1[34] != 0 && isec1[34] != 255) {
      if (isec1[34] & 2) {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,-isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,-isec1[35]);
        }
        *undef_eps = pow(10.0,-isec1[35]-1);
      } else {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,+isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,+isec1[35]);
        }
        *undef_eps = pow(10.0,isec1[35]-1);
      }
    }
  }
}

static
void cgribexDecodeHeader(int *isec0, int *isec1, int *isec2, double *fsec2,
			 int *isec3, double *fsec3, int *isec4, double *fsec4,
			 int *gribbuffer, int recsize, int *lmv, int *iret)
{
  int ipunp = 0, iword = 0;

  memset(isec1, 0, 256*sizeof(int));

  gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
	   ipunp, (int *) gribbuffer, recsize, &iword, "J", iret);

  *lmv = 0;

  if ( ISEC1_CenterID == 215 && (isec1[34] != 0 && isec1[34] != 255) )
    {
      double undef_pds, undef_eps;

      MCH_get_undef(isec1, &undef_pds, &undef_eps);
      FSEC3_MissVal = undef_pds;
      *lmv = 1;
    }
}

static
compvar_t cgribexVarSet(int param, int level1, int level2, int leveltype, int trange)
{
  compvar_t compVar;
  int tsteptype = cgribexGetTsteptype(trange);

  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;

  return compVar;
}

static inline int
cgribexVarCompare(compvar_t compVar, record_t record, int flag)
{
  int tstepDiff = (!((flag == 0) & (((compVar.tsteptype == TSTEP_INSTANT)
                                     & (record.tsteptype == TSTEP_INSTANT3))
                                    |((compVar.tsteptype == TSTEP_INSTANT3)
                                      & (record.tsteptype == TSTEP_INSTANT)))))
    & (compVar.tsteptype != record.tsteptype);
  int rstatus = (compVar.param != record.param)
    |           (compVar.level1 != record.ilevel)
    |           (compVar.level2 != record.ilevel2)
    |           (compVar.ltype != record.ltype)
    |           tstepDiff;
  return rstatus;
}
#endif

#define gribWarning(text, nrecs, timestep, paramstr, level1, level2) \
            Warning("Record %2d (id=%s lev1=%d lev2=%d) timestep %d: %s", nrecs, paramstr, level1, level2, timestep, text)

#if  defined  (HAVE_LIBCGRIBEX)

static inline void
cgribexScanTsFixNtsteps(stream_t *streamptr, off_t recpos)
{
  if ( streamptr->ntsteps == -1 )
    {
      int tsID = tstepsNewEntry(streamptr);
      if ( tsID != streamptr->rtsteps )
	Error("Internal error. tsID = %d", tsID);

      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
    }
}

static inline void
cgribexScanTsConstAdjust(stream_t *streamptr, taxis_t *taxis)
{
  int vlistID = streamptr->vlistID;
  if ( streamptr->ntsteps == 1 )
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
	{
	  streamptr->ntsteps = 0;
	  for (int varID = 0; varID < streamptr->nvars; varID++ )
	    {
	      vlistDefVarTsteptype(vlistID, varID, TSTEP_CONSTANT);
	    }
	}
    }
}


int cgribexScanTimestep1(stream_t * streamptr)
{
  double fsec2[512], fsec3[2], *fsec4 = NULL;
  int lmv = 0, iret = 0;
  off_t recpos = 0;
  void *gribbuffer = NULL;
  size_t buffersize = 0;
  int rstatus;
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
  DateTime datetime, datetime0 = { LONG_MIN, LONG_MIN };
  size_t readsize;
  unsigned nrecords, recID;
  int nrecs_scanned = 0;
  int datatype;
  long recsize = 0;
  bool warn_time = true;
  bool warn_numavg = true;
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0;
  bool fcast = false;
  int vlistID;
  int comptype;
  long unzipsize;
  char paramstr[32];
  int nskip = cdiSkipRecords;

  streamptr->curTsID = 0;

  int *isec0 = streamptr->record->sec0;
  int *isec1 = streamptr->record->sec1;
  int *isec2 = streamptr->record->sec2;
  int *isec3 = streamptr->record->sec3;
  int *isec4 = streamptr->record->sec4;

  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;

  if ( tsID != 0 )
    Error("Internal problem! tstepsNewEntry returns %d", tsID);

  int fileID = streamptr->fileID;

  while ( nskip-- > 0 )
    {
      recsize = gribGetSize(fileID);
      if ( recsize == 0 )
	Error("Skipping of %d records failed!", cdiSkipRecords);

      recpos  = fileGetPos(fileID);
      fileSetPos(fileID, (off_t)recsize, SEEK_CUR);
    }

  unsigned nrecs = 0;
  while ( TRUE )
    {
      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
	  if ( nrecs == 0 )
	    Error("No GRIB records found!");

	  streamptr->ntsteps = 1;
	  break;
	}
      if ( (size_t)recsize > buffersize )
	{
	  buffersize = (size_t)recsize;
	  gribbuffer = Realloc(gribbuffer, buffersize);
	}

      readsize = (size_t)recsize;
      rstatus = gribRead(fileID, (unsigned char *)gribbuffer, &readsize);
      if ( rstatus ) break;

      comptype = COMPRESS_NONE;
      if ( gribGetZip(recsize, (unsigned char *)gribbuffer, &unzipsize) > 0 )
	{
	  comptype = COMPRESS_SZIP;
	  unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	  if ( buffersize < (size_t)unzipsize )
	    {
	      buffersize = (size_t)unzipsize;
	      gribbuffer = Realloc(gribbuffer, buffersize);
	    }
	}

      nrecs_scanned++;
      cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
			  (int *) gribbuffer, (int)recsize, &lmv, &iret);

      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
      cdiParamToString(param, paramstr, sizeof(paramstr));

      if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
      if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
      level1   = ISEC1_Level1;
      level2   = ISEC1_Level2;

      gribDateTime(isec1, &vdate, &vtime);

      if ( ISEC4_NumBits > 0 && ISEC4_NumBits <= 32 )
	datatype = ISEC4_NumBits;
      else
        datatype = DATATYPE_PACK;

      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	  rdate = gribRefDate(isec1);
	  rtime = gribRefTime(isec1);
	  tunit = cgribexGetTimeUnit(isec1);
	  fcast = cgribexTimeIsFC(isec1);
	}
      else
	{
	  datetime.date  = vdate;
	  datetime.time  = vtime;
	  compvar_t compVar = cgribexVarSet(param, level1, level2, ISEC1_LevelType, ISEC1_TimeRange);
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
	      if ( cgribexVarCompare(compVar, streamptr->tsteps[0].records[recID], 0) == 0 ) break;
	    }

	  if ( cdiInventoryMode == 1 )
	    {
	      if ( recID < nrecs ) break;
	      if ( warn_time )
		if ( datetimeCmp(datetime, datetime0) != 0 )
		  {
                    gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, paramstr, level1, level2);
		    warn_time = false;
		  }
	    }
	  else
	    {
	      if ( datetimeCmp(datetime, datetime0) != 0 ) break;

	      if ( recID < nrecs )
		{
		  gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, paramstr, level1, level2);
		  continue;
		}
	    }
	}

      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
	    {
	      Warning("Changing numavg from %d to %d not supported!", taxis->numavg, ISEC1_AvgNum);
	      warn_numavg = false;
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}

      nrecs++;

      if ( CDI_Debug )
	Message("Read record %2d (id=%s lev1=%d lev2=%d) %8d %6d", nrecs_scanned, paramstr, level1, level2, vdate, vtime);

      cgribexAddRecord(streamptr, param, isec1, isec2, fsec2, fsec3,
		       isec4, recsize, recpos, datatype, comptype, lmv, iret);
    }

  streamptr->rtsteps = 1;

  if ( nrecs == 0 ) return CDI_EUFSTRUCT;

  cdi_generate_vars(streamptr);

  if ( fcast )
    {
      taxisID = taxisCreate(TAXIS_RELATIVE);
      taxis->type  = TAXIS_RELATIVE;
      taxis->rdate = rdate;
      taxis->rtime = rtime;
      taxis->unit  = tunit;
    }
  else
    {
      taxisID = taxisCreate(TAXIS_ABSOLUTE);
      taxis->type  = TAXIS_ABSOLUTE;
      taxis->unit  = tunit;
    }

  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;

  vlistID = streamptr->vlistID;
  vlistDefTaxis(vlistID, taxisID);

  nrecords = (unsigned)streamptr->tsteps[0].nallrecs;
  if ( nrecords < (unsigned)streamptr->tsteps[0].recordSize )
    {
      streamptr->tsteps[0].recordSize = (int)nrecords;
      streamptr->tsteps[0].records =
      (record_t *) Realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
    }

  streamptr->tsteps[0].recIDs = (int *) Malloc(nrecords*sizeof(int));
  streamptr->tsteps[0].nrecs = (int)nrecords;
  for ( recID = 0; recID < nrecords; recID++ )
    streamptr->tsteps[0].recIDs[recID] = (int)recID;

  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = (size_t)buffersize;

  cgribexScanTsFixNtsteps(streamptr, recpos);
  cgribexScanTsConstAdjust(streamptr, taxis);

  return 0;
}


int cgribexScanTimestep2(stream_t * streamptr)
{
  int rstatus = 0;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
  int lmv = 0, iret = 0;
  off_t recpos = 0;
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
  DateTime datetime, datetime0 = { LONG_MIN, LONG_MIN };
  int varID, gridID;
  size_t readsize;
  int nrecs, recID;
  long recsize = 0;
  bool warn_numavg = true;
  int tsteptype;
  long unzipsize;
  char paramstr[32];

  streamptr->curTsID = 1;

  int *isec0 = streamptr->record->sec0;
  int *isec1 = streamptr->record->sec1;
  int *isec2 = streamptr->record->sec2;
  int *isec3 = streamptr->record->sec3;
  int *isec4 = streamptr->record->sec4;

  int fileID  = streamptr->fileID;
  int vlistID = streamptr->vlistID;
  int taxisID = vlistInqTaxis(vlistID);

  void *gribbuffer = streamptr->record->buffer;
  size_t buffersize = streamptr->record->buffersize;

  int tsID = streamptr->rtsteps;
  if ( tsID != 1 )
    Error("Internal problem! unexpected timestep %d", tsID+1);

  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;

  fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);

  cdi_create_records(streamptr, tsID);

  int nrecords = streamptr->tsteps[tsID].nallrecs;
  if ( nrecords ) streamptr->tsteps[1].recIDs = (int *) Malloc((size_t)nrecords * sizeof(int));
  streamptr->tsteps[1].nrecs = 0;
  for ( recID = 0; recID < nrecords; recID++ )
    streamptr->tsteps[1].recIDs[recID] = -1;

  for ( recID = 0; recID < nrecords; recID++ )
    {
      varID = streamptr->tsteps[0].records[recID].varID;
      streamptr->tsteps[tsID].records[recID].position =	streamptr->tsteps[0].records[recID].position;
      streamptr->tsteps[tsID].records[recID].size     =	streamptr->tsteps[0].records[recID].size;
    }

  int nrecs_scanned = nrecords;
  int rindex = 0;
  while ( TRUE )
    {
      if ( rindex > nrecords ) break;

      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);
      if ( recsize == 0 )
	{
	  streamptr->ntsteps = 2;
	  break;
	}
      if ( (size_t)recsize > buffersize )
	{
	  buffersize = (size_t)recsize;
	  gribbuffer = Realloc(gribbuffer, buffersize);
	}

      readsize = (size_t)recsize;
      rstatus = gribRead(fileID, (unsigned char *)gribbuffer, &readsize);
      if ( rstatus ) break;

      if ( gribGetZip(recsize, (unsigned char *)gribbuffer, &unzipsize) > 0 )
	{
	  unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	  if ( buffersize < (size_t)unzipsize )
	    {
	      buffersize = (size_t)unzipsize;
	      gribbuffer = Realloc(gribbuffer, buffersize);
	    }
	}

      cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
			  (int *) gribbuffer, (int)recsize, &lmv, &iret);

      nrecs_scanned++;

      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
      cdiParamToString(param, paramstr, sizeof(paramstr));

      if ( ISEC1_LevelType == 100 )