stream_gribapi.c 110 KB
Newer Older
1
#ifdef  HAVE_CONFIG_H
2
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3 4
#endif

5
#ifdef  HAVE_LIBGRIB_API
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "cdi.h"
8
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
9
#include "file.h"
10
#include "gribapi_utilities.h"
11
#include "stream_scan.h"
12
#include "stream_grb.h"
13
#include "stream_gribapi.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
14 15 16
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
17
#include "calendar.h"
18
#include "subtype.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19

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

23
#include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24

25
extern int CDI_inventory_mode;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26

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

29

Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
  int param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32 33 34
  int level1;
  int level2;
  int ltype;
35
  int tsteptype;
36
  size_t gridsize;
37
  char name[32];
38
  VarScanKeys scanKeys;
39
  var_tile_t tiles;
40
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41 42 43


static
44
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
{
46
  return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47 48
}

49
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
int getTimeunits(const long unitsOfTime)
51
{
52 53
  switch (unitsOfTime)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54 55 56 57 58 59 60
    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;
61 62
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
  return TUNIT_HOUR;
64 65 66
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
double timeunit_factor(const int tu1, const int tu2)
68
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
  if (tu2 == TUNIT_HOUR)
70 71 72
    {
      switch (tu1)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73 74 75 76 77 78 79
        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;
80 81 82
        }
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
  return 1;
84 85 86 87 88
}

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

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

94
  return getTimeunits(unitsOfTime);
95 96
}

97
static
98
void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
99
{
100 101
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
102
  const int timeunits2 = (status == 0) ? getTimeunits(unitsOfTime) : timeunits;
103
  //timeunits2 = gribapiGetTimeUnits(gh);
104

105
  long lpar;
106 107 108 109 110 111 112 113
  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);
    }
114

115 116
  *endStep = *startStep;
  status = grib_get_long(gh, "endStep", &lpar);
117
  if ( status == 0 )
118 119
    *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));
120 121
}

122
static
123
void gribapiGetDataDateTime(grib_handle *gh, int64_t *datadate, int *datatime)
124
{
125 126 127 128 129 130 131 132
  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);
133 134
}

135
static
136
void gribapiSetDataDateTime(grib_handle *gh, int64_t datadate, int datatime)
137
{
138
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", (long)datadate), 0);
139 140 141 142 143
  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);
144 145
}

146
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
int gribapiGetTimeUnitFactor(const int timeUnits)
148
{
149
  static bool lprint = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150 151 152 153 154 155 156 157 158
  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;
159 160 161 162 163 164 165 166 167 168 169 170
    default:
      if ( lprint )
        {
          Warning("Time unit %d unsupported", timeUnits);
          lprint = false;
        }
    }

  return 0;
}

static
171
void gribapiGetValidityDateTime(grib_handle *gh, int64_t *vdate, int *vtime, int64_t *sDate, int *sTime)
172
{
173 174
  *sDate = 0;
  *sTime = 0;
175

176
  long sigofrtime = 3;
177
  if ( gribEditionNumber(gh) > 1 )
178
    GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
179
  else
180
    GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
181

182
  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())`.
183
    {
184
      gribapiGetDataDateTime(gh, vdate, vtime);
185 186 187
    }
  else
    {
188 189
      int64_t rdate;
      int rtime;
190 191
      gribapiGetDataDateTime(gh, &rdate, &rtime);

192
      const int timeUnits = gribapiGetTimeUnits(gh);
193
      int startStep = 0, endStep = 0;
194
      gribapiGetSteps(gh, timeUnits, &startStep, &endStep);
195

196 197 198 199 200
      {
	int ryear, rmonth, rday, rhour, rminute, rsecond;
	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

201 202
        if ( rday > 0 )
          {
203
            extern int CGRIBEX_grib_calendar;
204 205
            int64_t julday;
            int secofday;
206
            encode_caldaysec(CGRIBEX_grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);
207

208 209
            int64_t time_unit_factor = gribapiGetTimeUnitFactor(timeUnits);

210
            //if (startStep > 0)
211
              {
212 213 214 215
                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);
216 217
                *sDate = cdiEncodeDate(ryear, rmonth, rday);
                *sTime = cdiEncodeTime(rhour, rminute, 0);
218 219
              }

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

224 225 226
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
227 228 229
    }
}

230
static
231
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
232
{
233
  *leveltype = 0;
234 235 236
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
237

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

243
      switch (*leveltype)
244 245 246 247
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
248
	  { *lbounds = 1; break; }
249
	}
250

251 252 253
      if ( *lbounds )
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);  //1 byte
254
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
255 256
	  *level1 = (int)lpar;
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);       //1 byte
257 258
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
 	  *level2 = (int)lpar;
259 260
	}
      else
261
	{
262
          double dlevel;
263
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
264
	  if ( *leveltype == GRIB1_LTYPE_ISOBARIC ) dlevel *= 100;
265
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
266
	  if ( *leveltype == GRIB1_LTYPE_99 || *leveltype == GRIB1_LTYPE_ISOBARIC_PA ) *leveltype = GRIB1_LTYPE_ISOBARIC;
267

268 269 270
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
271
    }
272 273
}

274
static
275
double grib2ScaleFactor(const long factor)
276
{
277 278 279
  switch(factor)
    {
      case GRIB_MISSING_LONG: return 1;
280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
      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;
299 300 301 302 303 304 305 306
      default: return 0;
    }
}

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

312
static
313
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
314
                   int *level2, int *level_sf, int *level_unit)
315
{
316 317
  *leveltype1 = 0;
  *leveltype2 = -1;
318 319 320 321 322
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
323

324 325
  long lpar;
  int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
326
  if ( status == 0 )
327
    {
328
      *leveltype1 = (int) lpar;
329

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

334
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
      switch (*leveltype1)
336
        {
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354
          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;
355
        }
356

357
      long factor, llevel;
358 359 360
      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);
361

362 363 364 365 366
      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);
367
        }
368 369
    }
}
370

371
static
372
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)
373 374 375 376 377 378 379 380 381 382 383
{
  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);
384 385

      /* read in tiles attributes (if there are any) */
386
      tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], 0);
Thomas Jahns's avatar
Thomas Jahns committed
387 388 389 390 391
      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);
392
    }
393 394
}

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

400 401 402 403 404 405
  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);
    }
406 407 408 409
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
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;
}

431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447
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;
}

448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
static
void gribapiGetNameKeys(grib_handle *gh, int varID)
{
  char string[CDI_MAX_NAME];

  size_t vlen = CDI_MAX_NAME;
  gribapiGetString(gh, "name", string, vlen); // longname
  if (string[0]) varDefKeyString(varID, CDI_KEY_LONGNAME, string);

  gribapiGetString(gh, "units", string, vlen);
  if (string[0]) varDefKeyString(varID, CDI_KEY_UNITS, string);

  string[0] = 0;
  const int status = grib_get_string(gh, "cfName", string, &vlen);
  if ( status != 0 || vlen <= 1 || strncmp(string, "unknown", 7) == 0 ) string[0] = 0;
  if (string[0]) varDefKeyString(varID, CDI_KEY_STDNAME, string);
}

466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
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);
496 497 498 499

  long typeOfTimeIncrement = 0;
  if ( grib_get_long(gh, "typeOfTimeIncrement", &typeOfTimeIncrement) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFTIMEINCREMENT, (int) typeOfTimeIncrement);
500 501 502 503 504
  /*
  long constituentType = 0;
  if ( grib_get_long(gh, "constituentType", &constituentType) == 0 )
    varDefKeyInt(varID, CDI_KEY_CONSTITUENTTYPE, (int) constituentType);
  */
505 506 507 508 509

  /*
    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"
  */
510 511 512
  long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
  gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
  if ( perturbationNumber > 0 )
513
    {
514 515 516
      varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
      varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
      varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
517
    }
518

Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
  long section2Length = 0;
520 521
  int status = grib_get_long(gh, "section2Length", &section2Length);
  if (status == 0 && section2Length > 0)
522
    {
523
      long grib2LocalSectionNumber;
524
      long mpimType, mpimClass, mpimUser;
525 526
      status = grib_get_long(gh, "grib2LocalSectionNumber", &grib2LocalSectionNumber);
      if (status == 0)
527
        {
528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555
          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);
            }
556
        }
557
    }
558 559
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
560 561 562 563 564 565 566 567 568 569 570 571 572 573
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);
}

574 575 576 577 578
static
double shapeOfTheEarthToRadius(long shapeOfTheEarth)
{
  switch (shapeOfTheEarth)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
579 580 581 582
    case 0:  return 6367470.;
    case 2:  return 6378160.;
    case 6:  return 6371229.;
    case 8:  return 6371200.;
583 584 585 586
    }
  return 6367470.;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
587 588 589
static
void gribapiDefProjLCC(grib_handle *gh, int gridID)
{
590 591 592 593
  long shapeOfTheEarth;
  grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
  const double a = shapeOfTheEarthToRadius(shapeOfTheEarth);
  const double rf = (shapeOfTheEarth == 2) ? 297.0 : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
594 595 596 597 598 599 600 601 602 603 604 605
  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;
606 607
  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
608 609 610 611

  if ( proj_lonlat_to_lcc_func )
    {
      x_0 = xval_0; y_0 = yval_0;
612 613
      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) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614 615 616
        { x_0 = -x_0; y_0 = -y_0; }
    }

617
  gridDefParamLCC(gridID, CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618 619 620 621 622 623
}


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

  double lon_0, lat_ts, xval_0, yval_0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
629 630
  grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
  grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
631 632 633 634
  grib_get_double(gh, "LaDInDegrees", &lat_ts);
  grib_get_double(gh, "orientationOfTheGridInDegrees", &lon_0);

  long southPoleOnProjectionPlane;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635
  grib_get_long(gh, "southPoleOnProjectionPlane", &southPoleOnProjectionPlane);
636
  double lat_0 = southPoleOnProjectionPlane ? -90 : 90;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637

638 639
  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;
640 641

  if ( proj_lonlat_to_stere_func )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642 643
    {
      x_0 = xval_0; y_0 = yval_0;
644 645
      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) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
646 647
        { x_0 = -x_0; y_0 = -y_0; }
    }
648

649
  gridDefParamSTERE(gridID, CDI_grid_missval, lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650 651
}

652
static
653
void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
654
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
655
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
656
                      VarScanKeys *scanKeys, const var_tile_t *tiles, int lread_additional_keys)
657
{
658 659 660
  const int vlistID = streamptr->vlistID;
  const int tsID    = streamptr->curTsID;
  const int recID   = recordNewEntry(streamptr, tsID);
661
  record_t *record = &streamptr->tsteps[tsID].records[recID];
662

663
  const int tsteptype = gribapiGetTsteptype(gh);
664
  int numavg = 0;
665

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

Thomas Jahns's avatar
Thomas Jahns committed
668 669 670 671 672 673
  record->size      = recsize;
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
  record->ltype     = leveltype1;
Thomas Jahns's avatar
Thomas Jahns committed
674
  record->tsteptype = (short)tsteptype;
675
  record->gridsize  = gribapiGetGridsize(gh);
676 677
  record->scanKeys  = *scanKeys;
  record->tiles     = tiles ? *tiles : dummy_tiles;
678

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

682
  grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
683
  const bool uvRelativeToGrid = gribapiGetGrid(gh, grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
684

Uwe Schulzweida's avatar
Uwe Schulzweida committed
685
  struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
686
  const int gridID = gridAdded.Id;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
687 688 689 690
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
691

692
  const int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693

694
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
695
    {
696 697 698
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
699
        long lpar;
700
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
701
        /* FIXME: assert(lpar >= 0) */
702
        const size_t vctsize = (size_t)lpar;
703 704
        if ( vctsize > 0 )
          {
Thomas Jahns's avatar
Thomas Jahns committed
705 706
            double *vctptr = (double *) Malloc(vctsize*sizeof(double));
            size_t dummy = vctsize;
707
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
708
            varDefVCT(vctsize, vctptr);
709
            Free(vctptr);
710 711 712 713 714
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
715
        unsigned char uuid[CDI_UUID_SIZE];
716
        long lpar;
717
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
718
        if ( lpar != 6 ) fprintf(stderr, "Warning ...\n");
719 720 721 722
        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;
723
        size_t len = (size_t)CDI_UUID_SIZE;
724 725
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
726
        varDefZAxisReference(nhlev, nvgrid, uuid);
727 728
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730

731 732
  // if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = CDI_DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
733

734 735 736
  // add the previously read record data to the (intermediate) list of records
  int tile_index = 0;
  int varID = 0, levelID = 0;
737
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
738
	       datatype, &varID, &levelID, tsteptype, leveltype1, leveltype2,
739
	       varname, scanKeys, tiles, &tile_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
740

741
  record->varID = (short)varID;
742
  record->levelID = levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
743

744
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745

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

748
  if (varname[0]) gribapiGetNameKeys(gh, varID);
749
  gribapiGetKeys(gh, varID);
750

751
  if (lread_additional_keys)
752 753 754 755 756
    {
      long   lval;
      double dval;
      for ( int i = 0; i < cdiNAdditionalGRIBKeys; i++ )
        {
757
          // note: if the key is not defined, we do not throw an error!
758 759 760 761 762 763
          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]);
        }
    }
764

Uwe Schulzweida's avatar
Uwe Schulzweida committed
765 766
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
767 768 769
      long center, subcenter;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
770
      int instID = institutInq((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
771
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
772
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
773 774 775 776 777
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778
      long processID;
779
      if ( grib_get_long(gh, "generatingProcessIdentifier", &processID) == 0 )
780
	{
781
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
782
	  int modelID = modelInq(varInqInst(varID), (int)processID, NULL);
783
	  if ( modelID == CDI_UNDEFID )
784
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
785 786
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
787 788 789 790
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
791
      int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
792
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793

Uwe Schulzweida's avatar
Uwe Schulzweida committed
794 795 796
      if ( pdis == 255 )
	{
	  int tabnum = pcat;
797
	  int tableID = tableInq(varInqModel(varID), tabnum, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
798 799 800
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801 802 803 804 805 806 807
	}
    }

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

  if ( CDI_Debug )
808
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810
}
811

812
static
813 814
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)
815 816
{
  compvar2_t compVar;
817
  memset(&compVar, 0, sizeof(compvar2_t));
818
  const size_t maxlen = sizeof(compVar.name);
819 820
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
821

822 823