stream_gribapi.c 81.6 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"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
24
25
26
27
28
#  include "grib_api.h"
#endif

extern int cdiInventoryMode;

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


#if  defined  (HAVE_LIBGRIB_API)
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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);
}

66
static
67
int gribapiGetIsRotated(grib_handle *gh)
68
69
{
  int isRotated = 0;
70
  int gribgridtype = -1;
71
  long lpar;
72
  int status;
73

74
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
75

76
77
78
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
79
80
81
82

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
static
84
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
88
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
    {
90
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
93
    }
  else
    {
94
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
97
98
99
    }

  return (zaxistype);
}

100
static
101
int getTimeunits(long unitsOfTime)
102
103
{
  int timeunits = -1;
104
105
106
107
108
109
110
111
112
113
114
115
116

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

117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
  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;
146
  long unitsOfTime = -1;
147

148
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
149

150
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
151

152
  timeunits = getTimeunits(unitsOfTime);
153
154
155
156

  return (timeunits);
}

157
158
159
160
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
161
  int timeunits2 = timeunits;
162

163
164
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
165
166
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
167

168
  long lpar;
169
170
171
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
172
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
173
174
175
176
177
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

178
179
180
181
182
183
184
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
185
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
186
187
188
  *datatime = (int) lpar*100;
}

189
190
191
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
192
193
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
194
195
}

196
static
197
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
198
{
199
  int rdate, rtime;
200
  int timeUnits, startStep = 0, endStep;
201
202
  int tstepRange = 0;
  int range;
203
  int status;
204
  long lpar;
205
206
207
208
209
  long sigofrtime = 3;
  long editionNumber;

  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

210
  if ( editionNumber > 1 )
211
212
213
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
214
215
216
217
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
218
219
220

  if ( sigofrtime == 3 )
    {
221
      gribapiGetDataDateTime(gh, vdate, vtime);
222
223
224
    }
  else
    {
225
226
      gribapiGetDataDateTime(gh, &rdate, &rtime);

227
228
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
229
      timeUnits = gribapiGetTimeUnits(gh);
230
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
231
232
233
234
235
236
237
238
239

      range = endStep - startStep;

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

240
241
242
243
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
244
245
246
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
247
248
249
250

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

251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
        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);
          }
278

279
280
281
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
282
    }
283
284

  return (tstepRange);
285
286
}

287
static
288
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
289
{
290
  *leveltype = 0;
291
292
293
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
294

295
296
  long lpar;
  int status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
297
  if ( status == 0 )
298
    {
299
      *leveltype = (int) lpar;
300

301
      switch (*leveltype)
302
303
304
305
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
306
	  { *lbounds = 1; break; }
307
	}
308

309
310
      if ( *lbounds == 0 )
	{
311
          double dlevel;
312
313
314
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
315
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
316

317
318
319
320
321
322
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
      else
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
323
	  *level1 = (int)lpar;
324
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
325
	  *level2 = (int)lpar;
326
	}
327
    }
328
329
}

330
331
332
static
double grib2ScaleFactor(long factor)
{
333
  double scaleFactor = 0;
334

335
336
337
338
339
  if      ( factor == 0 ) scaleFactor =    1;
  else if ( factor == 1 ) scaleFactor =    0.1;
  else if ( factor == 2 ) scaleFactor =    0.01;
  else if ( factor == 3 ) scaleFactor =    0.001;
  else if ( factor == 4 ) scaleFactor =    0.0001;
340
341
342
343
344
  else if ( factor == 5 ) scaleFactor =    0.00001;
  else if ( factor == 6 ) scaleFactor =    0.000001;
  else if ( factor == 7 ) scaleFactor =    0.0000001;
  else if ( factor == 8 ) scaleFactor =    0.00000001;
  else if ( factor == 9 ) scaleFactor =    0.000000001;
345
346
347
348

  return (scaleFactor);
}

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

357
358
  *leveltype1 = 0;
  *leveltype2 = -1;
359
360
361
362
363
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
364
365
366

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
367
    {
368
      long llevel;
369
370
      double dlevel1 = 0, dlevel2 = 0;

371
      *leveltype1 = (int) lpar;
372

373
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
374
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
375
      if ( status == 0 ) *leveltype2 = (int)lpar;
376

377
378
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
      if ( *leveltype1 == GRIB2_LTYPE_REFERENCE && *leveltype2 == 1 ) *lbounds = 0;
379

380
      if ( *leveltype1 == GRIB2_LTYPE_LANDDEPTH )
381
382
383
384
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_M;
        }
385
      else if ( *leveltype1 == GRIB2_LTYPE_ISOBARIC )
386
387
388
389
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_PA;
        }
390
      else if ( *leveltype1 == GRIB2_LTYPE_SIGMA )
391
392
393
394
        {
          *level_sf = 1000;
          *level_unit = 0;
        }
395

396
397
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
398
399
400
      if ( llevel != GRIB_MISSING_LONG )
        {
          if ( factor != GRIB_MISSING_LONG )
401
            dlevel1 = (double)llevel * grib2ScaleFactor(factor);
402
          else
403
            dlevel1 = (double)llevel;
404
        }
405

406
      if ( *level_sf != 0 ) dlevel1 *= (double)(*level_sf);
407
408

      if ( *lbounds == 1 )
409
	{
410
411
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
412
413
414
          if ( llevel != GRIB_MISSING_LONG )
            {
              if ( factor != GRIB_MISSING_LONG )
415
                dlevel2 = (double)llevel * grib2ScaleFactor(factor);
416
              else
417
                dlevel2 = (double)llevel;
418
            }
419
420
421
422
423
424

          if ( *level_sf != 0 ) dlevel2 *= (*level_sf);
        }

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
425
    }
426
427
}

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

  GRIB_CHECK(grib_get_string(gh, key, string, &length), 0);
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

438
#if  defined  (HAVE_LIBGRIB_API)
439
static
440
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
441
		      size_t recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
442
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit)
443
444
445
446
447
448
449
450
451
452
453
454
455
{
  long editionNumber;
  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;
456
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
457
  size_t vlen;
458
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
459

Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
  vlistID = streamptr->vlistID;
461
  tsID    = streamptr->curTsID;
462
  recID   = recordNewEntry(streamptr, tsID);
463
464
465
466
467
468
469
470
  record  = &streamptr->tsteps[tsID].records[recID];

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

  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

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

473
474
475
476
477
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
478
  (*record).ltype     = leveltype1;
479
  (*record).tsteptype = tsteptype;
480
481
482
483
484
485
486
487
488
489
490
491
492

  //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));
  record->varname[sizeof(record->varname) - 1] = 0;
493
494

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

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

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

500
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
    {
502
503
504
505
506
507
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508

509
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
510
511
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
        vctsize = (int)lpar;
512
513
514
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
515
            dummy = (size_t)vctsize;
516
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
517
            varDefVCT((size_t)vctsize, vctptr);
518
519
520
521
522
523
524
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
525
        unsigned char uuid[CDI_UUID_SIZE];
526
        long ltmp;
527
        long nhlev, nvgrid;
528
529

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

546
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
547
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
548

549
  stdname[0] = 0;
550
551
552
  longname[0] = 0;
  units[0] = 0;

553
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554
    {
555
      vlen = CDI_MAX_NAME;
556
      gribapiGetString(gh, "name", longname, vlen);
557
      vlen = CDI_MAX_NAME;
558
      gribapiGetString(gh, "units", units, vlen);
559
560
561
562
563

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

569
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
570
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
571
	       varname, stdname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572

573
574
  (*record).varID   = (short)varID;
  (*record).levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575

576
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577

578
579
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
580
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
581
  */
582
583
584
585
586
587
588
  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);
    }

589
590
591
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

592
593
594
595
596
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

597
598
599
600
601
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
  int    i;
  long   lval;
  double dval;

  /* we read the additional keys for the first variable record only. */
  int linitial_field = (varOptGribNentries(varID) == 0);

  for ( i = 0; i < cdiNAdditionalGRIBKeys; i++ )
    {
      if ( linitial_field )
	{
	  if ( grib_get_long(gh, cdiAdditionalGRIBKeys[i], &lval) == 0 )
            varDefOptGribInt(varID, lval, cdiAdditionalGRIBKeys[i]);
	}
      if ( linitial_field )
	{
	  if ( grib_get_double(gh, cdiAdditionalGRIBKeys[i], &dval) == 0 )
619
            varDefOptGribDbl(varID, dval, cdiAdditionalGRIBKeys[i]);
620
621
622
	}
      /* note: if the key is not defined, we do not throw an error! */
    }
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
static
680
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype, char *name)
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
696
697
698

  return (compVar);
}

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

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

711
712
713
714
715
716
  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;
    }

717
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
718
719
720

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
722
#endif

723
724
725
#define gribWarning(text, nrecs, timestep, varname, paramstr, level1, level2) \
            Warning("Record %2d (name=%s id=%s lev1=%d lev2=%d) timestep %d: %s", nrecs, varname, paramstr, level1, level2, timestep, text)

726
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
728
729
730
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
731
  size_t buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
732
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
734
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
735
  int rtabnum = 0;
736
737
  int rcode = 0, level1 = 0, level2 = 0;
  int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
738
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
740
741
742
743
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
744
  int nrecs_scanned = 0;
745
  int datatype;
746
  size_t recsize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
  int warn_time = TRUE;
748
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
751
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752
  int vlistID;
753
  int comptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754
  long unzipsize;
755
  compvar2_t compVar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
756
  grib_handle *gh = NULL;
757
  int leveltype1, leveltype2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
758
759
  long editionNumber;
  long lpar;
760
  size_t len;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
761
  int bitsPerValue;
762
  int lieee = FALSE;
763
  int lbounds;
764
  int level_sf, level_unit;
765
  int tsteptype;
766
767
  char paramstr[32];
  char varname[256];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
768
769
770

  streamptr->curTsID = 0;

771
  tsID  = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
772
773
774
  taxis = &streamptr->tsteps[tsID].taxis;

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
777
  fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778
779
780
781

  nrecs = 0;
  while ( TRUE )
    {
782
783
      level1 = 0;
      level2 = 0;
784
      recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
786
787
788
789
790
791
792
793
794
795
796
797
798
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
	  streamptr->ntsteps = 1;
	  break;
	}
      if ( recsize > buffersize )
	{
	  buffersize = recsize;
	  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	}

      readsize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
799
800
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801

802
803
      lieee = FALSE;

804
      comptype = COMPRESS_NONE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
805

806
      nrecs_scanned++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
807
      gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
808
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
810
811
812
813

      GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

      if ( editionNumber <= 1 )
	{
814
815
816
817
818
819
820
821
822
823
824
          if ( gribGetZip((long)recsize, gribbuffer, &unzipsize) > 0 )
            {
              comptype = COMPRESS_SZIP;
              unzipsize += 100;
              if ( buffersize < (size_t)unzipsize )
                {
                  buffersize = (size_t)unzipsize;
                  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
                }
            }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
825
826
	  GRIB_CHECK(grib_get_long(gh, "table2Version", &lpar), 0);
	  rtabnum = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827
828
	  GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &lpar), 0);
	  rcode = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829

Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
	  param = cdiEncodeParam(rcode, rtabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831

832
	  grib1GetLevel(gh, &leveltype1, &lbounds, &level1, &level2);
833
          level_sf = 0;
834
          level_unit = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
835
836
837
	}
      else
	{
838
839
	  size_t len = 256;
	  char typeOfPacking[256];
840

841
	  status = grib_get_string(gh, "packingType", typeOfPacking, &len);
842
843
	  if ( status == 0 )
	    {
844
	      // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
845
846
847
	      if      ( strncmp(typeOfPacking, "grid_jpeg", len)  == 0 ) comptype = COMPRESS_JPEG;
	      else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) comptype = COMPRESS_SZIP;
	      else if ( strncmp(typeOfPacking, "grid_ieee", len)  == 0 ) lieee = TRUE;
848
	    }
849

850
	  param = gribapiGetParam(gh);
851

852
	  grib2GetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
853
854
	}

855
856
857
858
859
860
861
      cdiParamToString(param, paramstr, sizeof(paramstr));

      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
      len = strlen(varname);
      if ( len > 32 ) len = 32;
      //printf("param = %s  name = %s   l1 = %d  l2 = %d\n", paramstr, varname, level1, level2);
862

863
      tsteptype = gribapiGetTsteptype(gh);
864
865
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
      /*
866
      printf("%d %d %d\n", vdate, vtime, leveltype1);
867
      */
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
      if ( lieee )
        {
          datatype = DATATYPE_FLT64;
          status = grib_get_long(gh, "precision", &lpar);
          if ( status == 0 && lpar == 1 ) datatype = DATATYPE_FLT32;
        }
      else
        {
          datatype = DATATYPE_PACK;
          status = grib_get_long(gh, "bitsPerValue", &lpar);
          if ( status == 0 )
            {
              bitsPerValue = (int) lpar;
              if ( bitsPerValue > 0 && bitsPerValue <= 32 )
                datatype = bitsPerValue;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885

Uwe Schulzweida's avatar
Uwe Schulzweida committed
886
887
888
889
      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
890
891
892

          gribapiGetDataDateTime(gh, &rdate, &rtime);

893
894
	  fcast = gribapiTimeIsFC(gh);
	  if ( fcast ) tunit = gribapiGetTimeUnits(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
895
896
897
898
899
	}
      else
	{
	  datetime.date  = vdate;
	  datetime.time  = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
900

901
	  compVar = gribapiVarSet(param, level1, level2, leveltype1, tsteptype, varname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
902

Uwe Schulzweida's avatar
Uwe Schulzweida committed
903
	  for ( recID = 0; recID < nrecs; recID++ )
904
            if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
905
906
907
908
909
910
911

	  if ( cdiInventoryMode == 1 )
	    {
	      if ( recID < nrecs ) break;
	      if ( warn_time )
		if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 )
		  {
912
913
914
915
916
917
918
919
920
921
922
923
                    if ( datetime0.date == 10101 && datetime0.time == 0 )
                      {
                        datetime0.date = datetime.date;
                        datetime0.time = datetime.time;

                        gribapiGetDataDateTime(gh, &rdate, &rtime);

                        fcast = gribapiTimeIsFC(gh);
                        if ( fcast ) tunit = gribapiGetTimeUnits(gh);
                      }
                    else
                      {
924
                        gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, varname, paramstr, level1, level2);
925
926
927
                        warn_time = FALSE;
                      }
                  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928
929
930
931
932
933
934
	    }
	  else
	    {
	      if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;

	      if ( recID < nrecs )
		{
935
		  gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, paramstr, level1, level2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
936
937
938
939
940
941
942
943
944
		  continue;
		}
	    }
	}
      /*
      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
	    {
945
	      Message("Change numavg from %d to %d not allowed!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
946
947
948
949
950
951
952
953
954
955
956
957
		      taxis->numavg, ISEC1_AvgNum);
	      warn_numavg = FALSE;
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}
      */
      nrecs++;

      if ( CDI_Debug )
958
	Message("%4d %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%d vtime=%d",
959
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960

961
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, len, varname,
962
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);
963
964
965

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
967
    }

968
969
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
970
971
  streamptr->rtsteps = 1;

972
973
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

974
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989

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

990
991
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
992

Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
  vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
994