stream_gribapi.c

#ifdef  HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef  HAVE_LIBGRIB_API
#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
#include "file.h"
#include "gribapi_utilities.h"
#include "stream_scan.h"
#include "stream_grb.h"
#include "stream_gribapi.h"
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
#include "calendar.h"
#include "subtype.h"

#include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
#include "gribapi.h"

#include <grib_api.h>

extern int CDI_inventory_mode;

static const var_tile_t dummy_tiles = { 0, -1, -1, -1, -1, -1 };


typedef struct {
  int param;
  int level1;
  int level2;
  int ltype;
  int tsteptype;
  size_t gridsize;
  char name[32];
  VarScanKeys scanKeys;
  var_tile_t tiles;
} compvar2_t;


static
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
{
  return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype);
}

static
int getTimeunits(const long unitsOfTime)
{
  switch (unitsOfTime)
    {
    case 13:  return TUNIT_SECOND;
    case  0:  return TUNIT_MINUTE;
    case  1:  return TUNIT_HOUR;
    case 10:  return TUNIT_3HOURS;
    case 11:  return TUNIT_6HOURS;
    case 12:  return TUNIT_12HOURS;
    case  2:  return TUNIT_DAY;
    }

  return TUNIT_HOUR;
}

static
double timeunit_factor(const int tu1, const int tu2)
{
  if (tu2 == TUNIT_HOUR)
    {
      switch (tu1)
        {
        case TUNIT_SECOND:  return 3600;
        case TUNIT_MINUTE:  return 60;
        case TUNIT_HOUR:    return 1;
        case TUNIT_3HOURS:  return 1./3;
        case TUNIT_6HOURS:  return 1./6;
        case TUNIT_12HOURS: return 1./12;
        case TUNIT_DAY:     return 1./24;
        }
    }

  return 1;
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  long unitsOfTime = -1;
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);

  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);

  return getTimeunits(unitsOfTime);
}

static
void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
{
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
  const int timeunits2 = (status == 0) ? getTimeunits(unitsOfTime) : timeunits;
  //timeunits2 = gribapiGetTimeUnits(gh);

  long lpar;
  status = grib_get_long(gh, "forecastTime", &lpar);
  if ( status == 0 ) *startStep = (int) lpar;
  else
    {
      status = grib_get_long(gh, "startStep", &lpar);
      if ( status == 0 )
        *startStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
    }

  *endStep = *startStep;
  status = grib_get_long(gh, "endStep", &lpar);
  if ( status == 0 )
    *endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
  // printf("%d %d %d %d %d %g\n", *startStep, *endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));
}

static
void gribapiGetDataDateTime(grib_handle *gh, int64_t *datadate, int *datatime)
{
  long date;
  GRIB_CHECK(grib_get_long(gh, "dataDate", &date), 0);
  *datadate = (int64_t) date;
  long hour, minute, second;
  GRIB_CHECK(grib_get_long(gh, "hour", &hour), 0);
  GRIB_CHECK(grib_get_long(gh, "minute", &minute), 0);
  GRIB_CHECK(grib_get_long(gh, "second", &second), 0);
  *datatime = cdiEncodeTime((int)hour, (int)minute, (int)second);
}

static
void gribapiSetDataDateTime(grib_handle *gh, int64_t datadate, int datatime)
{
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", (long)datadate), 0);
  int hour, minute, second;
  cdiDecodeTime(datatime, &hour, &minute, &second);
  GRIB_CHECK(my_grib_set_long(gh, "hour", hour), 0);
  GRIB_CHECK(my_grib_set_long(gh, "minute", minute), 0);
  GRIB_CHECK(my_grib_set_long(gh, "second", second), 0);
}

static
int gribapiGetTimeUnitFactor(const int timeUnits)
{
  static bool lprint = true;
  switch (timeUnits)
    {
    case TUNIT_SECOND:  return     1;
    case TUNIT_MINUTE:  return    60;
    case TUNIT_HOUR:    return  3600;
    case TUNIT_3HOURS:  return 10800;
    case TUNIT_6HOURS:  return 21600;
    case TUNIT_12HOURS: return 43200;
    case TUNIT_DAY:     return 86400;
    default:
      if ( lprint )
        {
          Warning("Time unit %d unsupported", timeUnits);
          lprint = false;
        }
    }

  return 0;
}

static
void gribapiGetValidityDateTime(grib_handle *gh, int64_t *vdate, int *vtime, int64_t *sDate, int *sTime)
{
  *sDate = 0;
  *sTime = 0;

  long sigofrtime = 3;
  if ( gribEditionNumber(gh) > 1 )
    GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
  else
    GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);

  if ( sigofrtime == 3 ) //XXX: This looks like a bug to me, because timeRangeIndicator == 3 does not seem to have the same meaning as significanceOfReferenceTime == 3. I would recommend replacing this condition with `if(!gribapiTimeIsFC())`.
    {
      gribapiGetDataDateTime(gh, vdate, vtime);
    }
  else
    {
      int64_t rdate;
      int rtime;
      gribapiGetDataDateTime(gh, &rdate, &rtime);

      const int timeUnits = gribapiGetTimeUnits(gh);
      int startStep = 0, endStep = 0;
      gribapiGetSteps(gh, timeUnits, &startStep, &endStep);

      {
	int ryear, rmonth, rday, rhour, rminute, rsecond;
	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

        if ( rday > 0 )
          {
            extern int CGRIBEX_grib_calendar;
            int64_t julday;
            int secofday;
            encode_caldaysec(CGRIBEX_grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);

            int64_t time_unit_factor = gribapiGetTimeUnitFactor(timeUnits);

            //if (startStep > 0)
              {
                int64_t julday_x = julday;
                int secofday_x = secofday;
                julday_add_seconds(time_unit_factor * startStep, &julday_x, &secofday_x);
                decode_caldaysec(CGRIBEX_grib_calendar, julday_x, secofday_x, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
                *sDate = cdiEncodeDate(ryear, rmonth, rday);
                *sTime = cdiEncodeTime(rhour, rminute, 0);
              }

            julday_add_seconds(time_unit_factor * endStep, &julday, &secofday);
            decode_caldaysec(CGRIBEX_grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
          }

	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
    }
}

static
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
{
  *leveltype = 0;
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;

  long lpar;
  if ( !grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar) )       //1 byte
    {
      *leveltype = (int) lpar;

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

      if ( *lbounds )
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);  //1 byte
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
	  *level1 = (int)lpar;
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);       //1 byte
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
 	  *level2 = (int)lpar;
	}
      else
	{
          double dlevel;
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
	  if ( *leveltype == GRIB1_LTYPE_ISOBARIC ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
	  if ( *leveltype == GRIB1_LTYPE_99 || *leveltype == GRIB1_LTYPE_ISOBARIC_PA ) *leveltype = GRIB1_LTYPE_ISOBARIC;

	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
    }
}

static
double grib2ScaleFactor(const long factor)
{
  switch(factor)
    {
      case GRIB_MISSING_LONG: return 1;
      case -9: return 1000000000;
      case -8: return 100000000;
      case -7: return 10000000;
      case -6: return 1000000;
      case -5: return 100000;
      case -4: return 10000;
      case -3: return 1000;
      case -2: return 100;
      case -1: return 10;
      case  0: return 1;
      case  1: return 0.1;
      case  2: return 0.01;
      case  3: return 0.001;
      case  4: return 0.0001;
      case  5: return 0.00001;
      case  6: return 0.000001;
      case  7: return 0.0000001;
      case  8: return 0.00000001;
      case  9: return 0.000000001;
      default: return 0;
    }
}

static
int calcLevel(int level_sf, long factor, long level)
{
  double result = 0;
  if (level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
  if (level_sf) result *= level_sf;
  return (int)result;
}

static
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
                   int *level2, int *level_sf, int *level_unit)
{
  *leveltype1 = 0;
  *leveltype2 = -1;
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;

  long lpar;
  int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
  if ( status == 0 )
    {
      *leveltype1 = (int) lpar;

      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
      if ( status == 0 ) *leveltype2 = (int)lpar;

      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
      switch (*leveltype1)
        {
          case GRIB2_LTYPE_REFERENCE:
            if(*leveltype2 == 1) *lbounds = 0;
            break;

          case GRIB2_LTYPE_LANDDEPTH:
            *level_sf = 1000;
            *level_unit = CDI_UNIT_M;
            break;

          case GRIB2_LTYPE_ISOBARIC:
            *level_sf = 1000;
            *level_unit = CDI_UNIT_PA;
            break;

          case GRIB2_LTYPE_SIGMA:
            *level_sf = 1000;
            *level_unit = 0;
            break;
        }

      long factor, llevel;
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);      //1 byte
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);      //4 byte
      *level1 = calcLevel(*level_sf, factor, llevel);

      if ( *lbounds )
        {
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0); //1 byte
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0); //4 byte
          *level2 = calcLevel(*level_sf, factor, llevel);
        }
    }
}

static
void gribGetLevel(grib_handle *gh, int* leveltype1, int* leveltype2, int* lbounds, int* level1, int* level2, int* level_sf, int* level_unit, var_tile_t* tiles)
{
  if ( gribEditionNumber(gh) <= 1 )
    {
      grib1GetLevel(gh, leveltype1, lbounds, level1, level2);
      *leveltype2 = -1;
      *level_sf = 0;
      *level_unit = 0;
    }
  else
    {
      grib2GetLevel(gh, leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);

      /* read in tiles attributes (if there are any) */
      tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], 0);
      tiles->totalno_of_tileattr_pairs = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TOTALNO_OF_TILEATTR_PAIRS], -1);
      tiles->tileClassification = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILE_CLASSIFICATION], -1);
      tiles->numberOfTiles = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_TILES], -1);
      tiles->numberOfAttributes = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_ATTR], -1);
      tiles->attribute = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEATTRIBUTE], -1);
    }
}

static
void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length)
{
  string[0] = 0;

  int ret = grib_get_string(gh, key, string, &length);
  if (ret != 0)
    {
      fprintf(stderr, "grib_get_string(gh, \"%s\", ...) failed!\n", key);
      GRIB_CHECK(ret, 0);
    }
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

static
int gribapiGetEnsembleInfo(grib_handle *gh, long *typeOfEnsembleForecast, long *numberOfForecastsInEnsemble, long *perturbationNumber)
{
  int status = 0;
  if (grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast) == 0)
    {
      GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble), 0);
      GRIB_CHECK(grib_get_long(gh, "perturbationNumber", perturbationNumber), 0);
      if (*perturbationNumber > 0) status = 1;
    }

  if (status == 0)
    {
      *typeOfEnsembleForecast = 0;
      *perturbationNumber = 0;
      *numberOfForecastsInEnsemble = 0;
    }

  return status;
}

static
VarScanKeys gribapiGetScanKeys(grib_handle *gh)
{
  VarScanKeys scanKeys;
  varScanKeysInit(&scanKeys);

  long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
  gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
  scanKeys.perturbationNumber = (short)perturbationNumber;

  long typeOfGeneratingProcess = 0;
  if (grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess) == 0)
    scanKeys.typeOfGeneratingProcess = (short)typeOfGeneratingProcess;

  return scanKeys;
}

static
void gribapiGetKeys(grib_handle *gh, int varID)
{
  long tablesVersion = 0;
  if ( grib_get_long(gh, "tablesVersion", &tablesVersion) == 0 )
    varDefKeyInt(varID, CDI_KEY_TABLESVERSION, (int) tablesVersion);

  long localTablesVersion = 0;
  if ( grib_get_long(gh, "localTablesVersion", &localTablesVersion) == 0 )
    varDefKeyInt(varID, CDI_KEY_LOCALTABLESVERSION, (int) localTablesVersion);

  long typeOfGeneratingProcess = 0;
  if ( grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFGENERATINGPROCESS, (int) typeOfGeneratingProcess);

  long productDefinitionTemplate = 0;
  if ( grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate) == 0 )
    varDefKeyInt(varID, CDI_KEY_PRODUCTDEFINITIONTEMPLATE, (int) productDefinitionTemplate);

  long typeOfProcessedData = 0;
  if ( grib_get_long(gh, "typeOfProcessedData", &typeOfProcessedData) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFPROCESSEDDATA, (int) typeOfProcessedData);

  long shapeOfTheEarth = 0;
  if ( grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth) == 0 )
    varDefKeyInt(varID, CDI_KEY_SHAPEOFTHEEARTH, (int) shapeOfTheEarth);

  long backgroundProcess = 0;
  if ( grib_get_long(gh, "backgroundProcess", &backgroundProcess) == 0 )
    varDefKeyInt(varID, CDI_KEY_BACKGROUNDPROCESS, (int) backgroundProcess);

  long typeOfTimeIncrement = 0;
  if ( grib_get_long(gh, "typeOfTimeIncrement", &typeOfTimeIncrement) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFTIMEINCREMENT, (int) typeOfTimeIncrement);
  /*
  long constituentType = 0;
  if ( grib_get_long(gh, "constituentType", &constituentType) == 0 )
    varDefKeyInt(varID, CDI_KEY_CONSTITUENTTYPE, (int) constituentType);
  */

  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
  */
  long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
  gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
  if ( perturbationNumber > 0 )
    {
      varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
      varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
      varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
    }

  long section2Length = 0;
  int status = grib_get_long(gh, "section2Length", &section2Length);
  if (status == 0 && section2Length > 0)
    {
      long grib2LocalSectionNumber;
      long mpimType, mpimClass, mpimUser;
      status = grib_get_long(gh, "grib2LocalSectionNumber", &grib2LocalSectionNumber);
      if (status == 0)
        {
          size_t section2PaddingLength = 0;
          status = grib_get_size(gh, "section2Padding", &section2PaddingLength);
          if (status == 0 && section2PaddingLength > 0)
            {
              varDefKeyInt(varID, CDI_KEY_GRIB2LOCALSECTIONNUMBER, (int) grib2LocalSectionNumber);
              varDefKeyInt(varID, CDI_KEY_SECTION2PADDINGLENGTH, (int) section2PaddingLength);
              unsigned char *section2Padding = (unsigned char*) Malloc(section2PaddingLength);
              grib_get_bytes(gh, "section2Padding", section2Padding, &section2PaddingLength);
              varDefKeyBytes(varID, CDI_KEY_SECTION2PADDING, section2Padding, (int)section2PaddingLength);
              Free(section2Padding);
            }
          else if ( grib_get_long(gh, "mpimType", &mpimType) == 0 &&
                    grib_get_long(gh, "mpimClass", &mpimClass) == 0 &&
                    grib_get_long(gh, "mpimUser", &mpimUser) == 0)
            {
              varDefKeyInt(varID, CDI_KEY_MPIMTYPE, mpimType);
              varDefKeyInt(varID, CDI_KEY_MPIMCLASS, mpimClass);
              varDefKeyInt(varID, CDI_KEY_MPIMUSER, mpimUser);

              size_t revNumLen = 20;
              unsigned char revNumber[revNumLen];
              if ( grib_get_bytes(gh, "revNumber", revNumber, &revNumLen) == 0 )
                varDefKeyBytes(varID, CDI_KEY_REVNUMBER, revNumber, (int)revNumLen);

              long revStatus;
              grib_get_long(gh, "revStatus", &revStatus);
              varDefKeyInt(varID, CDI_KEY_REVSTATUS, revStatus);
            }
        }
    }
}

static
void gribapiDefProjRLL(grib_handle *gh, int gridID)
{
  double xpole = 0, ypole = 0, angle = 0;
  grib_get_double(gh, "latitudeOfSouthernPoleInDegrees",  &ypole);
  grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
  grib_get_double(gh, "angleOfRotation", &angle);
  xpole -= 180;
  if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
  if ( fabs(angle) > 0 ) angle = -angle;

  gridDefParamRLL(gridID, xpole, ypole, angle);
}

static
double shapeOfTheEarthToRadius(long shapeOfTheEarth)
{
  switch (shapeOfTheEarth)
    {
    case 0:  return 6367470.;
    case 2:  return 6378160.;
    case 6:  return 6371229.;
    case 8:  return 6371200.;
    }
  return 6367470.;
}

static
void gribapiDefProjLCC(grib_handle *gh, int gridID)
{
  long shapeOfTheEarth;
  grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
  const double a = shapeOfTheEarthToRadius(shapeOfTheEarth);
  const double rf = (shapeOfTheEarth == 2) ? 297.0 : 0;
  double lon_0, lat_1, lat_2, xval_0, yval_0;
  long projflag = 0;
  grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
  grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
  grib_get_double(gh, "LoVInDegrees", &lon_0);
  grib_get_double(gh, "Latin1InDegrees", &lat_1);
  grib_get_double(gh, "Latin2InDegrees", &lat_2);
  grib_get_long(gh, "projectionCentreFlag", &projflag);
  bool lsouth = gribbyte_get_bit((int)projflag, 1);
  if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }

  double lat_0 = lat_2;
  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;

  if ( proj_lonlat_to_lcc_func )
    {
      x_0 = xval_0; y_0 = yval_0;
      proj_lonlat_to_lcc_func(CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
      if ( IS_NOT_EQUAL(x_0, CDI_grid_missval) && IS_NOT_EQUAL(y_0, CDI_grid_missval) )
        { x_0 = -x_0; y_0 = -y_0; }
    }

  gridDefParamLCC(gridID, CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
}


static
void gribapiDefProjSTERE(grib_handle *gh, int gridID)
{
  long shapeOfTheEarth;
  grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
  const double a = shapeOfTheEarthToRadius(shapeOfTheEarth);

  double lon_0, lat_ts, xval_0, yval_0;
  grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
  grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
  grib_get_double(gh, "LaDInDegrees", &lat_ts);
  grib_get_double(gh, "orientationOfTheGridInDegrees", &lon_0);

  long southPoleOnProjectionPlane;
  grib_get_long(gh, "southPoleOnProjectionPlane", &southPoleOnProjectionPlane);
  double lat_0 = southPoleOnProjectionPlane ? -90 : 90;

  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;

  if ( proj_lonlat_to_stere_func )
    {
      x_0 = xval_0; y_0 = yval_0;
      proj_lonlat_to_stere_func(CDI_grid_missval, lon_0, lat_ts, lat_0, a, (size_t)1, &x_0, &y_0);
      if ( IS_NOT_EQUAL(x_0, CDI_grid_missval) && IS_NOT_EQUAL(y_0, CDI_grid_missval) )
        { x_0 = -x_0; y_0 = -y_0; }
    }

  gridDefParamSTERE(gridID, CDI_grid_missval, lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0);
}

static
void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
                      VarScanKeys *scanKeys, const var_tile_t *tiles, int lread_additional_keys)
{
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];

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

  const int tsteptype = gribapiGetTsteptype(gh);
  // numavg  = ISEC1_AvgNum;
  int numavg = 0;

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

  record->size      = recsize;
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
  record->ltype     = leveltype1;
  record->tsteptype = (short)tsteptype;
  record->gridsize  = gribapiGetGridsize(gh);
  record->scanKeys  = *scanKeys;
  record->tiles     = tiles ? *tiles : dummy_tiles;

  strncpy(record->varname, varname, sizeof(record->varname)-1);
  record->varname[sizeof(record->varname) - 1] = 0;

  grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
  const bool uvRelativeToGrid = gribapiGetGrid(gh, grid);

  struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
  const int gridID = gridAdded.Id;
  if      ( !gridAdded.isNew ) Free(grid);
  else if ( grid->projtype == CDI_PROJ_RLL ) gribapiDefProjRLL(gh, gridID);
  else if ( grid->projtype == CDI_PROJ_LCC ) gribapiDefProjLCC(gh, gridID);
  else if ( grid->projtype == CDI_PROJ_STERE ) gribapiDefProjSTERE(gh, gridID);

  const int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);

  switch (zaxistype)
    {
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        long lpar;
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
        /* FIXME: assert(lpar >= 0) */
        const size_t vctsize = (size_t)lpar;
        if ( vctsize > 0 )
          {
            double *vctptr = (double *) Malloc(vctsize*sizeof(double));
            size_t dummy = vctsize;
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
            varDefVCT(vctsize, vctptr);
            Free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        unsigned char uuid[CDI_UUID_SIZE];
        long lpar;
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
        if ( lpar != 6 ) fprintf(stderr, "Warning ...\n");
        GRIB_CHECK(grib_get_long(gh, "nlev", &lpar), 0);
        int nhlev = (int)lpar;
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &lpar), 0);
        int nvgrid = (int)lpar;
        size_t len = (size_t)CDI_UUID_SIZE;
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
        varDefZAxisReference(nhlev, nvgrid, uuid);
        break;
      }
    }

  // if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = CDI_DATATYPE_PACK;

  stdname[0] = 0;
  longname[0] = 0;
  units[0] = 0;

  if ( varname[0] != 0 )
    {
      size_t vlen = CDI_MAX_NAME;
      gribapiGetString(gh, "name", longname, vlen); // longname
      gribapiGetString(gh, "units", units, vlen);
      const int status = grib_get_string(gh, "cfName", stdname, &vlen);
      if ( status != 0 || vlen <= 1 || strncmp(stdname, "unknown", 7) == 0 )
        stdname[0] = 0;
    }
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);

  // add the previously read record data to the (intermediate) list of records
  int tile_index = 0;
  int varID = 0, levelID = 0;
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
	       varname, stdname, longname, units, scanKeys, tiles, &tile_index);

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

  varDefCompType(varID, comptype);

  if ( uvRelativeToGrid ) varDefKeyInt(varID, CDI_KEY_UVRELATIVETOGRID, 1);

  gribapiGetKeys(gh, varID);

  if (lread_additional_keys)
    {
      long   lval;
      double dval;
      for ( int i = 0; i < cdiNAdditionalGRIBKeys; i++ )
        {
          /* note: if the key is not defined, we do not throw an error! */
          if ( grib_get_long(gh, cdiAdditionalGRIBKeys[i], &lval) == 0 )
            varDefOptGribInt(varID, tile_index, lval, cdiAdditionalGRIBKeys[i]);
          if ( grib_get_double(gh, cdiAdditionalGRIBKeys[i], &dval) == 0 )
            varDefOptGribDbl(varID, tile_index, dval, cdiAdditionalGRIBKeys[i]);
        }
    }

  if ( varInqInst(varID) == CDI_UNDEFID )
    {
      long center, subcenter;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
      int instID = institutInq((int)center, (int)subcenter, NULL, NULL);
      if ( instID == CDI_UNDEFID )
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      long processID;
      if ( grib_get_long(gh, "generatingProcessIdentifier", &processID) == 0 )
	{
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
	  int modelID = modelInq(varInqInst(varID), (int)processID, NULL);
	  if ( modelID == CDI_UNDEFID )
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
	  varDefModel(varID, modelID);
	}
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
      int pdis, pcat, pnum;
      cdiDecodeParam(param, &pnum, &pcat, &pdis);

      if ( pdis == 255 )
	{
	  int tabnum = pcat;
	  int tableID = tableInq(varInqModel(varID), tabnum, NULL);
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
	}
    }

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

  if ( CDI_Debug )
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
	    varID, param, zaxistype, gridID, levelID);
}

static
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype,
                         size_t gridsize, char *name, VarScanKeys scanKeys, var_tile_t tiles_data)
{
  compvar2_t compVar;
  memset(&compVar, 0, sizeof(compvar2_t));
  const size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;

  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
  compVar.gridsize  = gridsize;
  //memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
  compVar.scanKeys = scanKeys;
  compVar.tiles = tiles_data;

  return compVar;
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
{
  compvar2_t compVar0;
  memset(&compVar0, 0, sizeof(compvar2_t));
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
  compVar0.gridsize  = record.gridsize;
  memcpy(compVar0.name, record.varname, sizeof(compVar.name));

  if ( flag == 0 )
    {
      if ( compVar0.tsteptype == TSTEP_INSTANT  && compVar.tsteptype == TSTEP_INSTANT3 ) compVar0.tsteptype = TSTEP_INSTANT3;
      if ( compVar0.tsteptype == TSTEP_INSTANT3 && compVar.tsteptype == TSTEP_INSTANT  ) compVar0.tsteptype = TSTEP_INSTANT;
    }

  compVar0.scanKeys = record.scanKeys;
  compVar0.tiles = record.tiles;

  return memcmp(&compVar0, &compVar, sizeof(compvar2_t));
}

static
grib_handle *gribapiGetDiskRepresentation(size_t recsize, size_t *buffersize, void **gribbuffer, int *outDatatype, int *outCompressionType)
{
  const int gribversion = (int)((char*)*gribbuffer)[7];

  if ( gribversion <= 1 ) *outCompressionType = grbDecompress(recsize, buffersize, gribbuffer);

  grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);

  bool lieee = false;

  if ( gribversion > 1 )
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = CDI_COMPRESS_JPEG;
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = CDI_COMPRESS_AEC;
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = true;
        }
    }

  if ( lieee )
    {
      *outDatatype = CDI_DATATYPE_FLT64;
      long precision;
      const int status = grib_get_long(gh, "precision", &precision);
      if ( status == 0 && precision == 1 ) *outDatatype = CDI_DATATYPE_FLT32;
    }
  else
    {
      *outDatatype = CDI_DATATYPE_PACK;
      long bitsPerValue;
      if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
        {
          if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
        }
    }

  return gh;
}

typedef enum { CHECKTIME_OK, CHECKTIME_SKIP, CHECKTIME_STOP, CHECKTIME_INCONSISTENT } checkTimeResult;

static checkTimeResult checkTime(stream_t* streamptr, compvar2_t compVar, const DateTime* verificationTime, const DateTime* expectedVTime) {
  // First determine whether the current record exists already.
  int recID = 0;
  for ( ; recID < streamptr->nrecs; recID++ )
    {
      if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
    }
  const bool recordExists =