stream_gribapi.c 81.9 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
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
                      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
487
488
489
490
491
492
493
494
495
496
497

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

719
720
  compVar0.tiles = record.tiles;

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

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

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

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

786
#if  defined  (HAVE_LIBGRIB_API)
787
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
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;
}
813
#endif
814
815
816
817
818
819
820
821

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

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

  streamptr->curTsID = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

952
953
  if ( gh ) grib_handle_delete(gh);

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

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

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

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

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

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

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

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

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

  if ( streamptr->ntsteps == -1 )
    {
999
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1000
      if ( tsID != streamptr->rtsteps )
1001
        Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
1003
1004
1005
1006
1007
1008
1009

      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
    }

  if ( streamptr->ntsteps ==