stream_gribapi.c 82 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"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20


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

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

extern int cdiInventoryMode;

29
static var_tile_t dummy_tiles = { -1, -1, -1, -1, -1, -1 };
30

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

  var_tile_t tiles;

41
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
43
44


#if  defined  (HAVE_LIBGRIB_API)
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
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
72
static
73
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
  int zaxistype = ZAXIS_GENERIC;

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

  return (zaxistype);
}

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

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

106
107
108
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
  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;
135
  long unitsOfTime = -1;
136

137
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
138

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

141
  timeunits = getTimeunits(unitsOfTime);
142
143
144
145

  return (timeunits);
}

146
147
148
149
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
150
  int timeunits2 = timeunits;
151

152
153
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
154
155
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
156

157
  long lpar;
158
159
160
  status = grib_get_long(gh, "endStep", &lpar);

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

  return (endStep);
}

167
168
169
170
171
172
173
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

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

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

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

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

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

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

      range = endStep - startStep;

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

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

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

237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        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);
          }
264

265
266
267
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
268
    }
269
270

  return (tstepRange);
271
272
}

273
static
274
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
275
{
276
  *leveltype = 0;
277
278
279
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
280

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

286
      switch (*leveltype)
287
288
289
290
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
291
	  { *lbounds = 1; break; }
292
	}
293

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

309
310
311
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
312
    }
313
314
}

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

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

352
353
  *leveltype1 = 0;
  *leveltype2 = -1;
354
355
356
357
358
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
359

360
  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
361
  if ( status == 0 )
362
    {
363
      long llevel;
364

365
      *leveltype1 = (int) lpar;
366

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

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

394
395
396
      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);
397

398
399
400
401
402
      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);
403
        }
404
405
    }
}
406

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

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

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

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

446
#if  defined  (HAVE_LIBGRIB_API)
447
static
448
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
449
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
450
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
451
                      const var_tile_t *tiles, int lread_additional_keys)
452
453
454
455
456
457
458
459
460
461
462
463
{
  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;
464
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
465
  size_t vlen;
466
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
467

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

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

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

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

  //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));
499
500

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

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

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

506
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
    {
508
509
510
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
Thomas Jahns's avatar
Thomas Jahns committed
511
        size_t vctsize;
512
513
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
514

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

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

552
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
553
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554

555
  stdname[0] = 0;
556
557
558
  longname[0] = 0;
  units[0] = 0;

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

      {
        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
570
        else if ( strncmp(stdname, "unknown", 7) == 0 ) stdname[0] = 0;
571
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572
    }
573
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574

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

581
582
  (*record).varID   = (short)varID;
  (*record).levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
583

584
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585

586
587
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
588
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
589
  */
590
591
592
593
594
595
596
  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);
    }

597
598
599
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

600
601
602
603
604
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

605
606
607
608
609
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

610
611
612
613
  int    i;
  long   lval;
  double dval;

614
615
616
617
618
619
620
621
622
  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]);
      }
623

Uwe Schulzweida's avatar
Uwe Schulzweida committed
624
625
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
627
628
629
630
      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
631
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
632
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
634
635
636
637
638
      varDefInst(varID, instID);
    }

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
658
659
660
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
664
665
666
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
667
668
669
670
671
672
673
	}
    }

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

  if ( CDI_Debug )
674
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
675
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
}
677
#endif
678

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

687
688
689
690
691
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
692
693
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
694

695
  compVar.tiles = tiles_data;
696
697
  return (compVar);
}
698
#endif
699

700
#ifdef HAVE_LIBGRIB_API
701
static
702
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
703
704
{
  compvar2_t compVar0;
705
  size_t maxlen = sizeof(compVar.name);
706

707
708
709
710
711
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
712
  memcpy(compVar0.name, record.varname, maxlen);
713

714
715
716
717
718
719
  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;
    }

720
721
  compVar0.tiles = record.tiles;

722
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
723
724
725

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

727
728
729
730
731
732
733
734
static void ensureBufferSize(size_t requiredSize, size_t* curSize, unsigned char** buffer) {
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
      *buffer = realloc(*buffer, *curSize);
    }
}

735
static
Thomas Jahns's avatar
Thomas Jahns committed
736
grib_handle* gribapiGetDiskRepresentation(size_t recsize, size_t* buffersize, unsigned char** gribbuffer, int* outDatatype, int* outCompressionType, long* outUnzipsize)
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
{
  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
756
      if( gribGetZip((long)recsize, *gribbuffer, outUnzipsize) > 0 )
757
758
        {
          *outCompressionType = COMPRESS_SZIP;
Thomas Jahns's avatar
Thomas Jahns committed
759
          ensureBufferSize((size_t)*outUnzipsize + 100, buffersize, gribbuffer);
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
        }
      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;
}
785
#endif
786

787
#if  defined  (HAVE_LIBGRIB_API)
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
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;
}
814
#endif
815
816
817
818
819
820
821
822

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

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

  streamptr->curTsID = 0;

839
840
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
841
842

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

845
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
846

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

      if ( recsize == 0 )
855
856
857
858
859
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
860

861
862
      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
863
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
864

865
      int datatype, comptype = 0;
866
867
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
868

869
      nrecs_scanned++;
870
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
871

872
873
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
874
      var_tile_t tiles = dummy_tiles;
875
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
876

877
      char varname[256];
878
879
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
880

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

      if( datetime0.date == 10101 && datetime0.time == 0 )
890
        {
891
892
893
894
895
896
897
898
899
          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);
            }
900
        }
901
902

      if(nrecs)
903
        {
904
          checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, varname, tiles), &datetime, &datetime0);
905
          if(result == CHECKTIME_STOP)
906
            {
907
              break;
908
            }
909
          else if(result == CHECKTIME_SKIP)
910
            {
911
912
              gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              continue;
913
            }
914
915
916
917
918
919
          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);
920
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921
922
      /*
      if ( ISEC1_AvgNum )
923
924
925
926
927
928
929
930
931
932
933
934
        {
          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
935
936
937
938
      */
      nrecs++;

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

946
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
947
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, &tiles, 1);
948
949
950

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
951
952
    }

953
954
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
956
  streamptr->rtsteps = 1;

957
958
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

959
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960

961
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
963
964
965
966
967
968
969
970
971
972
973
974
975
  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;
    }

976
977
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
978

979
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
980
981
  vlistDefTaxis(vlistID, taxisID);

982
  int nrecords = streamptr->tsteps[0].nallrecs;
983
  if ( nrecords > streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
984
985
986
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Thomas Jahns's avatar
Thomas Jahns committed
987
        (record_t *)xrealloc(streamptr->tsteps[0].records, (size_t)nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
988
989
    }

Thomas Jahns's avatar
Thomas Jahns committed
990
  streamptr->tsteps[0].recIDs = (int *)xmalloc((size_t)nrecords*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
991
  streamptr->tsteps[0].nrecs = nrecords;
992
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
994
995
996
997
998
999
    streamptr->tsteps[0].recIDs[recID] = recID;

  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = buffersize;

  if ( streamptr->ntsteps == -1 )
    {
1000
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1001
      if ( tsID != streamptr->rtsteps )