stream_gribapi.c 82.3 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
5
6
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <stdio.h>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
#include "cdi.h"
9
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
#include "file.h"
11
12
#include "gribapi_utilities.h"
#include "stream_grb.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


#if  defined  (HAVE_LIBGRIB_API)
22
#  include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
#  include "gribapi.h"
24
25

#  include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
27
28
29
#endif

extern int cdiInventoryMode;

Thomas Jahns's avatar
Thomas Jahns committed
30
#if  defined  (HAVE_LIBGRIB_API)
31
static const var_tile_t dummy_tiles = { -1, -1, -1, -1, -1, -1 };
Thomas Jahns's avatar
Thomas Jahns committed
32
#endif
33

Uwe Schulzweida's avatar
Uwe Schulzweida committed
34
typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
  int param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
38
  int level1;
  int level2;
  int ltype;
39
  int tsteptype;
40
  char name[32];
41
42
43

  var_tile_t tiles;

44
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
46
47


#if  defined  (HAVE_LIBGRIB_API)
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
static
int my_grib_set_double(grib_handle* h, const char* key, double val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_double(\tgrib_handle* h, \"%s\", %f)\n", key, val);

  return grib_set_double(h, key, val);
}

static
int my_grib_set_long(grib_handle* h, const char* key, long val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_long(  \tgrib_handle* h, \"%s\", %ld)\n", key, val);

  return grib_set_long(h, key, val);
}

static
int my_grib_set_string(grib_handle* h, const char* key, const char* val, size_t* length)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_string(\tgrib_handle* h, \"%s\", \"%s\")\n", key, val);

  return grib_set_string(h, key, val, length);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
static
76
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
80
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
    {
82
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
85
    }
  else
    {
86
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
89
90
91
    }

  return (zaxistype);
}

92
static
93
int getTimeunits(long unitsOfTime)
94
95
{
  int timeunits = -1;
96
97
98
99
100
101
102
103
104
105
106
107
108

  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;
    }

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
  return (timeunits);
}

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;
        }
    }

  return (factor);
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  int timeunits = -1;
138
  long unitsOfTime = -1;
139

140
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
141

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

144
  timeunits = getTimeunits(unitsOfTime);
145
146
147
148

  return (timeunits);
}

149
150
151
152
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
153
  int timeunits2 = timeunits;
154

155
156
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
157
158
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
159

160
  long lpar;
161
162
163
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
164
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
165
166
167
168
169
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

170
171
172
173
174
175
176
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
177
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
178
179
180
  *datatime = (int) lpar*100;
}

181
182
183
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
184
185
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
186
187
}

188
static
189
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
190
{
191
  int rdate, rtime;
192
  int timeUnits, startStep = 0, endStep;
193
194
  int tstepRange = 0;
  int range;
195
  int status;
196
  long lpar;
197
198
  long sigofrtime = 3;

199
  if ( gribEditionNumber(gh) > 1 )
200
201
202
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
203
204
205
206
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
207

208
  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())`.
209
    {
210
      gribapiGetDataDateTime(gh, vdate, vtime);
211
212
213
    }
  else
    {
214
215
      gribapiGetDataDateTime(gh, &rdate, &rtime);

216
217
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
218
      timeUnits = gribapiGetTimeUnits(gh);
219
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
220
221
222
223
224
225
226
227
228

      range = endStep - startStep;

      if ( range > 0 )
	{
	  if ( startStep == 0 ) tstepRange = -1;
	  else                  tstepRange =  1;
	}

229
230
231
232
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
233
234
235
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
236
237
238
239

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

240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
        if ( rday > 0 )
          {
            encode_caldaysec(grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);

            addsec = 0;
            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);
                    lprint = FALSE;
                  }
                break;
              }

            julday_add_seconds(addsec, &julday, &secofday);

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

268
269
270
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
271
    }
272
273

  return (tstepRange);
274
275
}

276
static
277
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
278
{
279
  *leveltype = 0;
280
281
282
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
283

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

289
      switch (*leveltype)
290
291
292
293
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
294
	  { *lbounds = 1; break; }
295
	}
296

297
298
299
300
301
302
303
304
      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
305
	{
306
          double dlevel;
307
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
308
309
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
310
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
311

312
313
314
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
315
    }
316
317
}

318
319
320
static
double grib2ScaleFactor(long factor)
{
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
  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;
Thomas Jahns's avatar
Thomas Jahns committed
342
  if(level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
343
344
  if(level_sf) result *= level_sf;
  return (int)result;
345
346
}

347
static
348
349
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
350
{
351
352
353
  int status;
  long lpar;
  long factor;
354

355
356
  *leveltype1 = 0;
  *leveltype2 = -1;
357
358
359
360
361
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
362

363
  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
364
  if ( status == 0 )
365
    {
366
      long llevel;
367

368
      *leveltype1 = (int) lpar;
369

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

374
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
375
      switch(*leveltype1)
376
        {
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
          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;
395
        }
396

397
398
399
      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);
400

401
402
403
404
405
      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);
406
        }
407
408
    }
}
409

410
static
411
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)
412
413
414
415
416
417
418
419
420
421
422
{
  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);
423
424

      /* read in tiles attributes (if there are any) */
Thomas Jahns's avatar
Thomas Jahns committed
425
426
427
428
429
430
      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);
431
    }
432
433
}

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

439
440
441
442
443
444
  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);
    }
445
446
447
448
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

449
#if  defined  (HAVE_LIBGRIB_API)
450
static
451
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
452
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
453
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
454
                      const var_tile_t *tiles, int lread_additional_keys)
455
456
457
458
459
460
461
462
463
464
465
466
{
  int zaxistype;
  int gridID = CDI_UNDEFID, varID;
  int levelID = 0;
  int tsID, recID;
  int numavg;
  int tsteptype;
  record_t *record;
  grid_t grid;
  int vlistID;
  long lpar;
  int status;
467
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
468
  size_t vlen;
469
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
470

Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
  vlistID = streamptr->vlistID;
472
  tsID    = streamptr->curTsID;
473
  recID   = recordNewEntry(streamptr, tsID);
474
475
476
477
478
479
  record  = &streamptr->tsteps[tsID].records[recID];

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

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

482
483
484
485
486
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
487
  (*record).ltype     = leveltype1;
488
  (*record).tsteptype = tsteptype;
489
490
  if ( tiles ) (*record).tiles = *tiles;
  else         (*record).tiles = dummy_tiles;
491
492
493
494
495
496
497
498
499
500
501
502

  //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));
503
504

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505

506
  gridID = varDefGrid(vlistID, &grid, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507

508
  zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
509

510
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
    {
512
513
514
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
Thomas Jahns's avatar
Thomas Jahns committed
515
        size_t vctsize;
516
517
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
518

519
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
520
521
        /* FIXME: assert(lpar >= 0) */
        vctsize = (size_t)lpar;
522
523
524
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
525
            dummy = vctsize;
526
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
527
            varDefVCT(vctsize, vctptr);
528
529
530
531
532
533
534
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
535
        unsigned char uuid[CDI_UUID_SIZE];
536
        long ltmp;
537
        long nhlev, nvgrid;
538
539

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
540
        if ( lpar != 6 )
541
542
543
          {
            fprintf(stderr, "Warning ...\n");
          }
544
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
545
        nhlev = ltmp;
546
547
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
548
549
550
        len = (size_t)CDI_UUID_SIZE;
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
551
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
552
553
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555

556
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
557
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
558

559
  stdname[0] = 0;
560
561
562
  longname[0] = 0;
  units[0] = 0;

563
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
    {
565
      vlen = CDI_MAX_NAME;
566
      gribapiGetString(gh, "name", longname, vlen);
567
      vlen = CDI_MAX_NAME;
568
      gribapiGetString(gh, "units", units, vlen);
569
570
571
572
573

      {
        vlen = CDI_MAX_NAME;
        status = grib_get_string(gh, "cfName", stdname, &vlen);
        if ( status != 0 || vlen <= 1 ) stdname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574
        else if ( strncmp(stdname, "unknown", 7) == 0 ) stdname[0] = 0;
575
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
    }
577
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
578

579
580
  /* add the previously read record data to the (intermediate) list of records */
  int tile_index = -1;
581
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
582
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
583
	       varname, stdname, longname, units, tiles, &tile_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
584

585
586
  (*record).varID   = (short)varID;
  (*record).levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
587

588
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
589

590
591
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
592
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
593
  */
594
595
596
597
598
599
600
  status = grib_get_long(gh, "typeOfEnsembleForecast", &ens_forecast_type );
  if ( status == 0 )
    {
      GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", &ens_count ), 0);
      GRIB_CHECK(grib_get_long(gh, "perturbationNumber", &ens_index ), 0);
    }

601
602
603
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

604
605
606
607
608
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

609
610
611
612
613
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

614
615
616
617
  int    i;
  long   lval;
  double dval;

618
619
620
621
622
623
624
625
626
  if (lread_additional_keys)
    for ( 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]);
      }
627

Uwe Schulzweida's avatar
Uwe Schulzweida committed
628
629
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
630
631
632
633
634
      long center, subcenter;
      int instID;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
      instID    = institutInq((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
636
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637
638
639
640
641
642
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643
      long processID;
644
645
646
      status = grib_get_long(gh, "generatingProcessIdentifier", &processID);
      if ( status == 0 )
	{
647
648
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
	  modelID = modelInq(varInqInst(varID), (int)processID, NULL);
649
	  if ( modelID == CDI_UNDEFID )
650
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
651
652
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
654
655
656
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
      int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
658

Uwe Schulzweida's avatar
Uwe Schulzweida committed
659
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
660

Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
662
663
664
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665

Uwe Schulzweida's avatar
Uwe Schulzweida committed
666
	  tableID = tableInq(varInqModel(varID), tabnum, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
667

Uwe Schulzweida's avatar
Uwe Schulzweida committed
668
669
670
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
672
673
674
675
676
677
	}
    }

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

  if ( CDI_Debug )
678
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
679
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
680
}
681
#endif
682

683
684
static compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, 
                                int tsteptype, char *name, var_tile_t tiles_data)
685
686
{
  compvar2_t compVar;
687
688
689
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
690

691
692
693
694
695
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
696
697
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
698

699
  compVar.tiles = tiles_data;
700
701
  return (compVar);
}
702
#endif
703

704
#ifdef HAVE_LIBGRIB_API
705
static
706
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
707
708
{
  compvar2_t compVar0;
709
  size_t maxlen = sizeof(compVar.name);
710

711
712
713
714
715
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
716
  memcpy(compVar0.name, record.varname, maxlen);
717

718
719
720
721
722
723
  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;
    }

724
725
  compVar0.tiles = record.tiles;

726
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
727
728
729

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730

731
732
733
734
735
736
737
738
static void ensureBufferSize(size_t requiredSize, size_t* curSize, unsigned char** buffer) {
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
      *buffer = realloc(*buffer, *curSize);
    }
}

739
static
Thomas Jahns's avatar
Thomas Jahns committed
740
grib_handle* gribapiGetDiskRepresentation(size_t recsize, size_t* buffersize, unsigned char** gribbuffer, int* outDatatype, int* outCompressionType, long* outUnzipsize)
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
{
  int lieee = FALSE;

  grib_handle* gh = grib_handle_new_from_message(NULL, (void *) *gribbuffer, recsize);
  if(gribEditionNumber(gh) > 1)
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = COMPRESS_JPEG;
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = COMPRESS_SZIP;
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = TRUE;
        }
    }
  else
    {
Thomas Jahns's avatar
Thomas Jahns committed
760
      if( gribGetZip((long)recsize, *gribbuffer, outUnzipsize) > 0 )
761
762
        {
          *outCompressionType = COMPRESS_SZIP;
Thomas Jahns's avatar
Thomas Jahns committed
763
          ensureBufferSize((size_t)*outUnzipsize + 100, buffersize, gribbuffer);
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
        }
      else
        {
          *outCompressionType = COMPRESS_NONE;
        }
    }

  if ( lieee )
    {
      *outDatatype = DATATYPE_FLT64;
      long precision;
      int status = grib_get_long(gh, "precision", &precision);
      if ( status == 0 && precision == 1 ) *outDatatype = DATATYPE_FLT32;
    }
  else
    {
      *outDatatype = DATATYPE_PACK;
      long bitsPerValue;
      if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
        {
          if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
        }
    }
  return gh;
}
789
#endif
790

791
#if  defined  (HAVE_LIBGRIB_API)
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
typedef enum { CHECKTIME_OK, CHECKTIME_SKIP, CHECKTIME_STOP, CHECKTIME_INCONSISTENT } checkTimeResult;
static checkTimeResult checkTime(stream_t* streamptr, compvar2_t compVar, const DateTime* verificationTime, const DateTime* expectedVTime) {
  //First determine whether the current record exists already.
  int recID = 0;
  for ( ; recID < streamptr->nrecs; recID++ )
    {
      if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
    }
  int recordExists = recID < streamptr->nrecs;

  //Then we need to know whether the verification time is consistent.
  int consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));

  //Finally, we make a decision.
  if ( cdiInventoryMode == 1 )
    {
      if ( recordExists ) return CHECKTIME_STOP;
      if ( !consistentTime ) return CHECKTIME_INCONSISTENT;
    }
  else
    {
      if ( !consistentTime ) return CHECKTIME_STOP;
      if ( recordExists ) return CHECKTIME_SKIP;
    }
  return CHECKTIME_OK;
}
818
#endif
819
820
821
822
823
824
825
826

#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)
827

828
#if  defined  (HAVE_LIBGRIB_API)
829
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
831
832
{
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
833
  size_t buffersize = 0;
834
835
  DateTime datetime0 = { .date = 10101, .time = 0 };
  int nrecs_scanned = 0;        //Only used for debug output.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836
  int warn_time = TRUE;
837
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
838
839
840
841
842
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
  grib_handle *gh = NULL;

  streamptr->curTsID = 0;

843
844
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
845
846

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

849
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
850

Thomas Jahns's avatar
Thomas Jahns committed
851
  unsigned nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
852
853
  while ( TRUE )
    {
854
855
      int level1 = 0, level2 = 0;
      size_t recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
856
857
858
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
859
860
861
862
863
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
864

865
866
      size_t readsize = recsize;
      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
867
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
868

869
      int datatype, comptype = 0;
870
871
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872

873
      nrecs_scanned++;
874
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
875

876
877
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
878
      var_tile_t tiles = dummy_tiles;
879
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
880

881
      char varname[256];
882
883
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
884

885
886
      int tsteptype = gribapiGetTsteptype(gh);
      int vdate = 0, vtime = 0;
887
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
888
      DateTime datetime = { .date = vdate, .time = vtime };
889
      /*
890
      printf("%d %d %d\n", vdate, vtime, leveltype1);
891
      */
892
893

      if( datetime0.date == 10101 && datetime0.time == 0 )
894
        {
895
896
897
898
899
900
901
902
903
          if( memcmp(&datetime, &datetime0, sizeof(datetime)) || !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.
            {
              datetime0 = datetime;

              gribapiGetDataDateTime(gh, &rdate, &rtime);

              fcast = gribapiTimeIsFC(gh);
              if ( fcast ) tunit = gribapiGetTimeUnits(gh);
            }
904
        }
905
906

      if(nrecs)
907
        {
908
          checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, varname, tiles), &datetime, &datetime0);
909
          if(result == CHECKTIME_STOP)
910
            {
911
              break;
912
            }
913
          else if(result == CHECKTIME_SKIP)
914
            {
915
916
              gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              continue;
917
            }
918
919
920
921
922
923
          else if(result == CHECKTIME_INCONSISTENT && warn_time)
            {
              gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              warn_time = FALSE;
            }
          assert(result == CHECKTIME_OK || result == CHECKTIME_INCONSISTENT);
924
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
925
926
      /*
      if ( ISEC1_AvgNum )
927
928
929
930
931
932
933
934
935
936
937
938
        {
          if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
            {
              Message("Change numavg from %d to %d not allowed!",
                      taxis->numavg, ISEC1_AvgNum);
              warn_numavg = FALSE;
            }
          else
            {
              taxis->numavg = ISEC1_AvgNum;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
939
940
941
942
      */
      nrecs++;

      if ( CDI_Debug )
943
944
945
        {
          char paramstr[32];
          cdiParamToString(param, paramstr, sizeof(paramstr));
Thomas Jahns's avatar
Thomas Jahns committed
946
          Message("%4u %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%d vtime=%d",
947
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
948
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
949

950
951
      var_tile_t *ptiles = NULL;
      if ( memcmp(&tiles, &dummy_tiles, sizeof(var_tile_t)) != 0 ) ptiles = &tiles;
952
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
953
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, ptiles, 1);
954
955
956

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
    }

959
960
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
961
962
  streamptr->rtsteps = 1;

963
964
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

965
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966

967
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
968
969
970
971
972
973
974
975
976
977
978
979
980
981
  if ( fcast )
    {
      taxisID = taxisCreate(TAXIS_RELATIVE);
      taxis->type  = TAXIS_RELATIVE;
      taxis->rdate = rdate;
      taxis->rtime = rtime;
      taxis->unit  = tunit;
    }
  else
    {
      taxisID = taxisCreate(TAXIS_ABSOLUTE);
      taxis->type  = TAXIS_ABSOLUTE;
    }

982
983
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
984

985
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
986
987
  vlistDefTaxis(vlistID, taxisID);

988
  int nrecords = streamptr->tsteps[0].nallrecs;
989
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
990
991
992
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Thomas Jahns's avatar
Thomas Jahns committed
993
        (record_t *)xrealloc(streamptr->tsteps[0].records, (size_t)nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
994
995
    }

Thomas Jahns's avatar
Thomas Jahns committed
996
  streamptr->tsteps[0].recIDs = (int *)xmalloc((size_t)nrecords*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
997
  streamptr->tsteps[0].nrecs = nrecords;