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
  // numavg  = ISEC1_AvgNum;
665
  int numavg = 0;
666

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

Thomas Jahns's avatar
Thomas Jahns committed
669
670
671
672
673
674
  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
675
  record->tsteptype = (short)tsteptype;
676
  record->gridsize  = gribapiGetGridsize(gh);
677
678
  record->scanKeys  = *scanKeys;
  record->tiles     = tiles ? *tiles : dummy_tiles;
679

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

823
824
825
826
827
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
828
  compVar.gridsize  = gridsize;
829
  //memset(compVar.name, 0, maxlen);
830
  memcpy(compVar.name, name, len);
831
  compVar.scanKeys = scanKeys;
832
  compVar.tiles = tiles_data;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
833

834
  return compVar;
835
836
837
}

static
838
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
839
840
{
  compvar2_t compVar0;
841
  memset(&compVar0, 0, sizeof(compvar2_t));
842
843
844
845
846
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
847
  compVar0.gridsize  = record.gridsize;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
848
  memcpy(compVar0.name, record.varname, sizeof(compVar.name));
849

850
851
852
853
854
855
  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;
    }

856
  compVar0.scanKeys = record.scanKeys;
857
858
  compVar0.tiles = record.tiles;

Thomas Jahns's avatar
Thomas Jahns committed
859
  return memcmp(&compVar0, &compVar, sizeof(compvar2_t));
860
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
861

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

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

869
  grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);
870
871
872
873

  bool lieee = false;

  if ( gribversion > 1 )
874
875
876
877
878
879
880
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
881
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = CDI_COMPRESS_JPEG;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
882
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = CDI_COMPRESS_AEC;
883
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = true;
884
885
886
887
888
        }
    }

  if ( lieee )
    {
889
      *outDatatype = CDI_DATATYPE_FLT64;
890
      long precision;
891
      const int status = grib_get_long(gh, "precision", &precision);
892
      if ( status == 0 && precision == 1 ) *outDatatype = CDI_DATATYPE_FLT32;
893
894
895
    }
  else
    {
896
      *outDatatype = CDI_DATATYPE_PACK;
897
898
899
900
901
902
      long bitsPerValue;
      if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
        {
          if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
        }
    }
903

904
905
906
907
  return gh;
}

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

909
static checkTimeResult checkTime(stream_t* streamptr, compvar2_t compVar, const DateTime* verificationTime, const DateTime* expectedVTime) {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
910
  // First determine whether the current record exists already.
911
912
913
914
915
  int recID = 0;
  for ( ; recID < streamptr->nrecs; recID++ )
    {
      if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
    }
916
  const bool recordExists = recID < streamptr->nrecs;
917

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
918
  // Then we need to know whether the verification time is consistent.
919
  const bool consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));
920

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
921
  // Finally, we make a decision.
922
  if ( CDI_inventory_mode == 1 )
923
924
925
926
927
928