stream_gribapi.c 105 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
11
#include "gribapi_utilities.h"
#include "stream_grb.h"
12
#include "stream_gribapi.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
13
14
15
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
#include "stream_grb.h"
17
#include "calendar.h"
18
#include "subtype.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19
20


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

24
#include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27

extern int cdiInventoryMode;

28
static const var_tile_t dummy_tiles = { 0, -1, -1, -1, -1, -1 };
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
38
39
40
41
42
43
44
#ifdef HIRLAM_EXTENSIONS
    // NOTE: tsteptype MUST be part of attributes used to compare variables!
    // Modern NWP models (HARMONIE, HIRLAM) use timeRangeIndicator to specify
    // if the field is instantanous or accumulated.
    // Both types are typically in the same GRIB-file.
    // (181; 105, 0, timeRangeIndicator=0) .. instantanous rain
    // (181; 105, 0, timeRangeIndicator=4) .. accumulated rain  .. both can be in the same grib file
#endif // HIRLAM_EXTENSIONS
45
  char name[32];
46
47
48

  var_tile_t tiles;

49
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52


static
53
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
{
55
  return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
}

58
static
59
int getTimeunits(long unitsOfTime)
60
61
{
  int timeunits = -1;
62
63
64
65
66
67
68
69
70
71
72
73
74

  switch (unitsOfTime)
    {
    case 13:  timeunits = TUNIT_SECOND;  break;
    case  0:  timeunits = TUNIT_MINUTE;  break;
    case  1:  timeunits = TUNIT_HOUR;    break;
    case 10:  timeunits = TUNIT_3HOURS;  break;
    case 11:  timeunits = TUNIT_6HOURS;  break;
    case 12:  timeunits = TUNIT_12HOURS; break;
    case  2:  timeunits = TUNIT_DAY;     break;
    default:  timeunits = TUNIT_HOUR;    break;
    }

75
  return timeunits;
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
}

static
double timeunit_factor(int tu1, int tu2)
{
  double factor = 1;

  if ( tu2 == TUNIT_HOUR )
    {
      switch (tu1)
        {
        case TUNIT_SECOND:  factor = 3600;   break;
        case TUNIT_MINUTE:  factor = 60;     break;
        case TUNIT_HOUR:    factor = 1;      break;
        case TUNIT_3HOURS:  factor = 1./3;   break;
        case TUNIT_6HOURS:  factor = 1./6;   break;
        case TUNIT_12HOURS: factor = 1./12;  break;
        case TUNIT_DAY:     factor = 1./24;  break;
        }
    }

97
  return factor;
98
99
100
101
102
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
103
  long unitsOfTime = -1;
104

105
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
106

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

109
  int timeunits = getTimeunits(unitsOfTime);
110

111
  return timeunits;
112
113
}

114
static
115
void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
116
{
117
  int timeunits2 = timeunits;
118
119
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
120
121
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
122

123
  long lpar;
124
125
126
127
128
129
130
131
  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);
    }
132

133
134
  *endStep = *startStep;
  status = grib_get_long(gh, "endStep", &lpar);
135
  if ( status == 0 )
136
137
    *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));
138
139
}

140
141
142
143
144
145
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;
  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
146
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
147
148
149
  *datatime = (int) lpar*100;
}

150
151
152
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
153
154
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
155
156
}

157
static
158
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
159
{
160
  int tstepRange = 0;
161

162
  long sigofrtime = 3;
163
  if ( gribEditionNumber(gh) > 1 )
164
    GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
165
  else
166
    GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
167

168
  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())`.
169
    {
170
      gribapiGetDataDateTime(gh, vdate, vtime);
171
172
173
    }
  else
    {
174
      int rdate, rtime;
175
176
      gribapiGetDataDateTime(gh, &rdate, &rtime);

177
178
      int timeUnits = gribapiGetTimeUnits(gh);
      int startStep = 0, endStep = 0;
179
      gribapiGetSteps(gh, timeUnits, &startStep, &endStep);
180

181
      int range = endStep - startStep;
182
183
184
185
186
187
      if ( range > 0 )
	{
	  if ( startStep == 0 ) tstepRange = -1;
	  else                  tstepRange =  1;
	}

188
189
190
191
192
      {
	int ryear, rmonth, rday, rhour, rminute, rsecond;
	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

193
194
        if ( rday > 0 )
          {
195
196
197
            static bool lprint = true;
            extern int grib_calendar;
            int julday, secofday;
198
199
            encode_caldaysec(grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);

200
201
            int64_t time_period = endStep;
            int64_t addsec = 0;
202
203
204
205
206
207
208
209
210
211
212
213
214
            switch ( timeUnits )
              {
              case TUNIT_SECOND:  addsec =         time_period; break;
              case TUNIT_MINUTE:  addsec =    60 * time_period; break;
              case TUNIT_HOUR:    addsec =  3600 * time_period; break;
              case TUNIT_3HOURS:  addsec = 10800 * time_period; break;
              case TUNIT_6HOURS:  addsec = 21600 * time_period; break;
              case TUNIT_12HOURS: addsec = 43200 * time_period; break;
              case TUNIT_DAY:     addsec = 86400 * time_period; break;
              default:
                if ( lprint )
                  {
                    Warning("Time unit %d unsupported", timeUnits);
215
                    lprint = false;
216
217
218
219
220
221
222
223
                  }
                break;
              }

            julday_add_seconds(addsec, &julday, &secofday);

            decode_caldaysec(grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
          }
224

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

230
  return tstepRange;
231
232
}

233
static
234
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
235
{
236
  *leveltype = 0;
237
238
239
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
240

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

246
      switch (*leveltype)
247
248
249
250
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
251
	  { *lbounds = 1; break; }
252
	}
253

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

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

275
276
277
static
double grib2ScaleFactor(long factor)
{
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
  switch(factor)
    {
      case GRIB_MISSING_LONG: return 1;
      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;
299
300
  if (level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
  if (level_sf) result *= level_sf;
301
  return (int)result;
302
303
}

304
static
305
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
306
                   int *level2, int *level_sf, int *level_unit)
307
{
308
309
  *leveltype1 = 0;
  *leveltype2 = -1;
310
311
312
313
314
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
315

316
317
  long lpar;
  int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
318
  if ( status == 0 )
319
    {
320
      long llevel;
321

322
      *leveltype1 = (int) lpar;
323

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

328
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
329
      switch(*leveltype1)
330
        {
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
          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;
349
        }
350

351
      long factor;
352
353
354
      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);
355

356
357
358
359
360
      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);
361
        }
362
363
    }
}
364

365
static
366
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)
367
368
369
370
371
372
373
374
375
376
377
{
  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);
378
379

      /* read in tiles attributes (if there are any) */
Thomas Jahns's avatar
Thomas Jahns committed
380
381
382
383
384
385
      tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], -1);
      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);
386
    }
387
388
}

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

394
395
396
397
398
399
  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);
    }
400
401
402
403
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
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);
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451

  /*
    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;
  if ( grib_get_long(gh, "typeOfEnsembleForecast", &typeOfEnsembleForecast) == 0 )
    {
      long perturbationNumber = 0, numberOfForecastsInEnsemble = 0;
      GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", &numberOfForecastsInEnsemble), 0);
      GRIB_CHECK(grib_get_long(gh, "perturbationNumber", &perturbationNumber), 0);
      if ( perturbationNumber > 0 )
        {
          varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
          varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
          varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
        }
    }
452

453
454
455
  long section2Length;
  grib_get_long(gh, "section2Length", &section2Length);
  if ( section2Length > 0 )
456
    {
457
458
459
460
461
462
463
464
465
466
      long grib2LocalSectionNumber;
      grib_get_long(gh, "grib2LocalSectionNumber", &grib2LocalSectionNumber);
      size_t size;
      grib_get_size(gh, "section2Padding", &size);
      if ( size > 0 )
        {
          varDefKeyInt(varID, CDI_KEY_GRIB2LOCALSECTIONNUMBER, (int) grib2LocalSectionNumber);
          varDefKeyInt(varID, CDI_KEY_SECTION2PADDINGLENGTH, (int) size);
          unsigned char *section2Padding = (unsigned char*) Malloc(size);
          grib_get_bytes(gh, "section2Padding", section2Padding, &size);
467
          varDefKeyBytes(varID, CDI_KEY_SECTION2PADDING, section2Padding, (int)size);
468
469
          Free(section2Padding);
        }
470
    }
471
472
}

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

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

486
  int tsteptype = gribapiGetTsteptype(gh);
487
  // numavg  = ISEC1_AvgNum;
488
  int numavg = 0;
489

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

Thomas Jahns's avatar
Thomas Jahns committed
492
493
494
495
496
497
  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
498
  record->tsteptype = (short)tsteptype;
499
  record->gridsize  = gribapiGetGridsize(gh);
Thomas Jahns's avatar
Thomas Jahns committed
500
  record->tiles = tiles ? *tiles : dummy_tiles;
501
502
503
504
505
506
507
508
509
510
511
512

  //FIXME: This may leave the variable name unterminated (which is the behavior that I found in the code).
  //       I don't know precisely how this field is used, so I did not change this behavior to avoid regressions,
  //       but I think that it would be better to at least add a line
  //
  //           record->varname[sizeof(record->varname) - 1] = 0;`
  //
  //       after the `strncpy()` call.
  //
  //       I would consider using strdup() (that requires POSIX-2008 compliance, though), or a similar homebrew approach.
  //       I. e. kick the fixed size array and allocate enough space, whatever that may be.
  strncpy(record->varname, varname, sizeof(record->varname));
513

514
  grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
515
  gribapiGetGrid(gh, grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
516

Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
  struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
518
  int gridID = gridAdded.Id;
519
520
521
522
523
524
525
  if ( !gridAdded.isNew ) Free(grid);
  else if ( grid->projtype == CDI_PROJ_RLL )
    {
      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);
526
      xpole -= 180;
527
528
529
      if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
      if ( fabs(angle) > 0 ) angle = -angle;

530
      gridDefParamRLL(gridID, xpole, ypole, angle);
531
    }
532
533
534
535
536
537
  else if ( grid->projtype == CDI_PROJ_LCC )
    {
      double a = 6367470., rf = 0;
      long earthIsOblate;
      grib_get_long(gh, "earthIsOblate", &earthIsOblate);
      if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
538
      double lon_0, lat_1, lat_2, xval_0, yval_0;
539
540
541
542
543
544
      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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
      grib_get_long(gh, "projectionCentreFlag", &projflag);
Thomas Jahns's avatar
Thomas Jahns committed
546
      bool lsouth = gribbyte_get_bit((int)projflag, 1);
547
548
      if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }

549
550
551
552
553
554
555
      double lat_0 = lat_2;
      double x_0 = grid_missval;
      double y_0 = grid_missval;

      if ( proj_lonlat_to_lcc_func )
        {
          x_0 = xval_0; y_0 = yval_0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
          proj_lonlat_to_lcc_func(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
557
558
559
          if ( IS_NOT_EQUAL(x_0, grid_missval) && IS_NOT_EQUAL(y_0, grid_missval) )
            { x_0 = -x_0; y_0 = -y_0; }
        }
560
      gridDefParamLCC(gridID, grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
561
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562

563
  int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564

565
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
    {
567
568
569
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
570
        long lpar;
571
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
572
        /* FIXME: assert(lpar >= 0) */
Thomas Jahns's avatar
Thomas Jahns committed
573
        size_t vctsize = (size_t)lpar;
574
575
        if ( vctsize > 0 )
          {
Thomas Jahns's avatar
Thomas Jahns committed
576
577
            double *vctptr = (double *) Malloc(vctsize*sizeof(double));
            size_t dummy = vctsize;
578
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
579
            varDefVCT(vctsize, vctptr);
580
            Free(vctptr);
581
582
583
584
585
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
586
        unsigned char uuid[CDI_UUID_SIZE];
587
        long lpar;
588
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
589
        if ( lpar != 6 ) fprintf(stderr, "Warning ...\n");
590
591
592
593
        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;
594
        size_t len = (size_t)CDI_UUID_SIZE;
595
596
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
597
        varDefZAxisReference(nhlev, nvgrid, uuid);
598
599
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
600
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
601

602
603
  // if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = CDI_DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604

605
  stdname[0] = 0;
606
607
608
  longname[0] = 0;
  units[0] = 0;

609
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
610
    {
611
      size_t vlen = CDI_MAX_NAME;
612
      gribapiGetString(gh, "name", longname, vlen);
613
      vlen = CDI_MAX_NAME;
614
      gribapiGetString(gh, "units", units, vlen);
Thomas Jahns's avatar
Thomas Jahns committed
615
616
617
618
      vlen = CDI_MAX_NAME;
      int status = grib_get_string(gh, "cfName", stdname, &vlen);
      if ( status != 0 || vlen <= 1 || strncmp(stdname, "unknown", 7) == 0 )
        stdname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
    }
620
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
621

622
  /* add the previously read record data to the (intermediate) list of records */
623
  int levelID = 0;
624
  int tile_index = 0, varID;
625
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
626
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
627
	       varname, stdname, longname, units, tiles, &tile_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
628

Thomas Jahns's avatar
Thomas Jahns committed
629
630
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631

632
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633

634
  gribapiGetKeys(gh, varID);
635

636
  if (lread_additional_keys)
637
638
639
640
641
642
643
644
645
646
647
648
    {
      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]);
        }
    }
649

Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
651
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
652
653
654
      long center, subcenter;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
655
      int instID = institutInq((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
656
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
658
659
660
661
662
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
663
      long processID;
664
      if ( grib_get_long(gh, "generatingProcessIdentifier", &processID) == 0 )
665
	{
666
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
667
	  int modelID = modelInq(varInqInst(varID), (int)processID, NULL);
668
	  if ( modelID == CDI_UNDEFID )
669
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
670
671
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
674
675
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
      int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
678

Uwe Schulzweida's avatar
Uwe Schulzweida committed
679
680
681
      if ( pdis == 255 )
	{
	  int tabnum = pcat;
682
	  int tableID = tableInq(varInqModel(varID), tabnum, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
683
684
685
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
686
687
688
689
690
691
692
	}
    }

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

  if ( CDI_Debug )
693
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
695
}
696

697
698
699
static
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype,
                         size_t gridsize, char *name, var_tile_t tiles_data)
700
701
{
  compvar2_t compVar;
702
  memset(&compVar, 0, sizeof(compvar2_t));
703
704
705
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
706

707
708
709
710
711
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
712
  compVar.gridsize  = gridsize;
713
  //memset(compVar.name, 0, maxlen);
714
  memcpy(compVar.name, name, len);
715
  compVar.tiles = tiles_data;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
716

717
  return compVar;
718
719
720
}

static
721
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
722
723
{
  compvar2_t compVar0;
724
  memset(&compVar0, 0, sizeof(compvar2_t));
725
726
727
728
729
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
730
  compVar0.gridsize  = record.gridsize;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
731
  memcpy(compVar0.name, record.varname, sizeof(compVar.name));
732

733
734
735
736
737
738
  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;
    }

739
740
  compVar0.tiles = record.tiles;

Thomas Jahns's avatar
Thomas Jahns committed
741
  return memcmp(&compVar0, &compVar, sizeof(compvar2_t));
742
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
743

744
745
746
static
void ensureBufferSize(size_t requiredSize, size_t *curSize, void **buffer)
{
747
748
749
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
750
      *buffer = Realloc(*buffer, *curSize);
751
752
753
    }
}

754
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
grib_handle *gribapiGetDiskRepresentation(size_t recsize, size_t *buffersize, void **gribbuffer, int *outDatatype, int *outCompressionType, size_t *outUnzipsize)
756
{
757
758
759
760
761
762
763
764
765
766
767
768
769
770
  int gribversion = (int)((char*)*gribbuffer)[7];

  if ( gribversion <= 1 )
    {
      if ( gribGetZip(recsize, *gribbuffer, outUnzipsize) > 0 )
        {
          *outCompressionType = CDI_COMPRESS_SZIP;
          ensureBufferSize(*outUnzipsize + 100, buffersize, gribbuffer);
        }
      else
        {
          *outCompressionType = CDI_COMPRESS_NONE;
        }
    }
771

772
  grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);
773
774
775
776

  bool lieee = false;

  if ( gribversion > 1 )
777
778
779
780
781
782
783
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
784
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = CDI_COMPRESS_JPEG;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = CDI_COMPRESS_AEC;
786
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = true;
787
788
789
790
791
        }
    }

  if ( lieee )
    {
792
      *outDatatype = CDI_DATATYPE_FLT64;
793
794
      long precision;
      int status = grib_get_long(gh, "precision", &precision);
795
      if ( status == 0 && precision == 1 ) *outDatatype = CDI_DATATYPE_FLT32;
796
797
798
    }
  else
    {
799
      *outDatatype = CDI_DATATYPE_PACK;
800
801
802
803
804
805
      long bitsPerValue;
      if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
        {
          if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
        }
    }
806

807
808
809
810
811
  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) {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
812
  // First determine whether the current record exists already.
813
814
815
816
817
818
819
  int recID = 0;
  for ( ; recID < streamptr->nrecs; recID++ )
    {
      if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
    }
  int recordExists = recID < streamptr->nrecs;

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
820
  // Then we need to know whether the verification time is consistent.
821
822
  int consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
823
  // Finally, we make a decision.
824
825
826
827
828
829
830
831
832
833
  if ( cdiInventoryMode == 1 )
    {
      if ( recordExists ) return CHECKTIME_STOP;
      if ( !consistentTime ) return CHECKTIME_INCONSISTENT;
    }
  else
    {
      if ( !consistentTime ) return CHECKTIME_STOP;
      if ( recordExists ) return CHECKTIME_SKIP;
    }
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
834

835
836
837
838
839
840
841
842
843
844
  return CHECKTIME_OK;
}

#define gribWarning(text, nrecs, timestep, varname, param, level1, level2) do \
  { \
    char paramstr[32]; \
    cdiParamToString(param, paramstr, sizeof(paramstr)); \
    Warning("Record %2d (name=%s id=%s lev1=%d lev2=%d) timestep %d: %s", nrecs, varname, paramstr, level1, level2, timestep, text); \
  } \
while(0)
845

846
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
847
848
{
  off_t recpos = 0;
849
  void *gribbuffer = NULL;
850
  size_t buffersize = 0;
851
852
  DateTime datetime0 = { .date = 10101, .time = 0 };
  int nrecs_scanned = 0;        //Only used for debug output.
853
854
  bool warn_time = true;
  // bool warn_numavg = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855
856
857
858
859
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
  grib_handle *gh = NULL;

  streamptr->curTsID = 0;

860
861
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
862
863

  if ( tsID != 0 )
864
    Error("Internal problem! tstepsNewEntry returns %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
865

866
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
867

Thomas Jahns's avatar
Thomas Jahns committed
868
  unsigned nrecs = 0;
869
  while ( true )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
870
    {
871
      int level1 = 0, level2 = 0;
872
873
      size_t recsize = gribGetSize(fileID);
      recpos = fileGetPos(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
874
875

      if ( recsize == 0 )
876
877
878
879
880
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
881

882
      size_t readsize = recsize;
883
      int rstatus = gribRead(fileID, gribbuffer, &readsize); //Search for next 'GRIB', read the following record, and position file offset after it.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
884
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885

886
      int datatype, comptype = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887
      size_t unzipsize;
888
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
889

890
      nrecs_scanned++;
891
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
892

893
894
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
895
      var_tile_t tiles = dummy_tiles;
896
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
897

898
      char varname[256];
899
900
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
901

902
      int tsteptype = gribapiGetTsteptype(gh);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
903

904
      int vdate = 0, vtime = 0;
905
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
906
907
      DateTime datetime = { .date = vdate, .time = vtime };

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
908
      if ( datetime0.date == 10101 && datetime0.time == 0 )
909
        {
910
          if ( datetimeCmp(datetime, datetime0) || !nrecs ) //Do we really need this condition? I have included it in order not to change the number of times gribapiGetDataDateTime() etc. get called. But if those are sideeffect-free, this condition should be removed.
911
912
913
914
915
916
917
918
            {
              datetime0 = datetime;

              gribapiGetDataDateTime(gh, &rdate, &rtime);

              fcast = gribapiTimeIsFC(gh);
              if ( fcast ) tunit = gribapiGetTimeUnits(gh);
            }
919
        }
920

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
921
      if ( nrecs )
922
        {
923
924
925
          size_t gridsize = gribapiGetGridsize(gh);
          checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, gridsize, varname, tiles),
                                             &datetime, &datetime0);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
926
          if ( result == CHECKTIME_STOP )
927
            {
928
              break;
929
            }