stream_gribapi.c 80.1 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
29
#endif

extern int cdiInventoryMode;

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


#if  defined  (HAVE_LIBGRIB_API)
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
66
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
67
static
68
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
72
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
    {
74
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
    }
  else
    {
78
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
81
82
83
    }

  return (zaxistype);
}

84
static
85
int getTimeunits(long unitsOfTime)
86
87
{
  int timeunits = -1;
88
89
90
91
92
93
94
95
96
97
98
99
100

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

101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
  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;
130
  long unitsOfTime = -1;
131

132
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
133

134
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
135

136
  timeunits = getTimeunits(unitsOfTime);
137
138
139
140

  return (timeunits);
}

141
142
143
144
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
145
  int timeunits2 = timeunits;
146

147
148
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
149
150
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
151

152
  long lpar;
153
154
155
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
156
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
157
158
159
160
161
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

162
163
164
165
166
167
168
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
169
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
170
171
172
  *datatime = (int) lpar*100;
}

173
174
175
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
176
177
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
178
179
}

180
static
181
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
182
{
183
  int rdate, rtime;
184
  int timeUnits, startStep = 0, endStep;
185
186
  int tstepRange = 0;
  int range;
187
  int status;
188
  long lpar;
189
190
  long sigofrtime = 3;

191
  if ( gribEditionNumber(gh) > 1 )
192
193
194
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
195
196
197
198
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
199

200
  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())`.
201
    {
202
      gribapiGetDataDateTime(gh, vdate, vtime);
203
204
205
    }
  else
    {
206
207
      gribapiGetDataDateTime(gh, &rdate, &rtime);

208
209
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
210
      timeUnits = gribapiGetTimeUnits(gh);
211
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
212
213
214
215
216
217
218
219
220

      range = endStep - startStep;

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

221
222
223
224
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
225
226
227
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
228
229
230
231

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

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
        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);
          }
259

260
261
262
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
263
    }
264
265

  return (tstepRange);
266
267
}

268
static
269
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
270
{
271
  *leveltype = 0;
272
273
274
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
275

276
  long lpar;
277
  if(!grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar))       //1 byte
278
    {
279
      *leveltype = (int) lpar;
280

281
      switch (*leveltype)
282
283
284
285
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
286
	  { *lbounds = 1; break; }
287
	}
288

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

304
305
306
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
307
    }
308
309
}

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

339
static
340
341
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
342
{
343
344
345
  int status;
  long lpar;
  long factor;
346

347
348
  *leveltype1 = 0;
  *leveltype2 = -1;
349
350
351
352
353
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
354

355
  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
356
  if ( status == 0 )
357
    {
358
      long llevel;
359

360
      *leveltype1 = (int) lpar;
361

362
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
363
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
364
      if ( status == 0 ) *leveltype2 = (int)lpar;
365

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

389
390
391
      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);
392

393
394
395
396
397
      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);
398
        }
399
400
    }
}
401

402
403
404
405
406
407
408
409
410
411
412
413
414
static
void gribGetLevel(grib_handle *gh, int* leveltype1, int* leveltype2, int* lbounds, int* level1, int* level2, int* level_sf, int* level_unit)
{
  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);
415
    }
416
417
}

418
419
420
421
422
423
424
425
426
427
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;
}

428
#if  defined  (HAVE_LIBGRIB_API)
429
static
430
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
431
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
432
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit)
433
434
435
436
437
438
439
440
441
442
443
444
{
  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;
445
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
446
  size_t vlen;
447
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
448

Uwe Schulzweida's avatar
Uwe Schulzweida committed
449
  vlistID = streamptr->vlistID;
450
  tsID    = streamptr->curTsID;
451
  recID   = recordNewEntry(streamptr, tsID);
452
453
454
455
456
457
  record  = &streamptr->tsteps[tsID].records[recID];

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

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

460
461
462
463
464
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
465
  (*record).ltype     = leveltype1;
466
  (*record).tsteptype = tsteptype;
467
468
469
470
471
472
473
474
475
476
477
478

  //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));
479
480

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

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

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

486
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
487
    {
488
489
490
491
492
493
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494

495
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
496
497
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
        vctsize = (int)lpar;
498
499
500
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
501
            dummy = (size_t)vctsize;
502
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
503
            varDefVCT((size_t)vctsize, vctptr);
504
505
506
507
508
509
510
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
511
        unsigned char uuid[CDI_UUID_SIZE];
512
        long ltmp;
513
        long nhlev, nvgrid;
514
515

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
516
        if ( lpar != 6 )
517
518
519
          {
            fprintf(stderr, "Warning ...\n");
          }
520
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
521
        nhlev = ltmp;
522
523
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
524
525
526
        len = (size_t)CDI_UUID_SIZE;
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
527
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
528
529
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
531

532
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
533
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534

535
  stdname[0] = 0;
536
537
538
  longname[0] = 0;
  units[0] = 0;

539
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
    {
541
      vlen = CDI_MAX_NAME;
542
      gribapiGetString(gh, "name", longname, vlen);
543
      vlen = CDI_MAX_NAME;
544
      gribapiGetString(gh, "units", units, vlen);
545
546
547
548
549

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

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

559
560
  (*record).varID   = (short)varID;
  (*record).levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561

562
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563

564
565
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
566
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
567
  */
568
569
570
571
572
573
574
  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);
    }

575
576
577
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

578
579
580
581
582
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

583
584
585
586
587
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
  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 )
605
            varDefOptGribDbl(varID, dval, cdiAdditionalGRIBKeys[i]);
606
607
608
	}
      /* note: if the key is not defined, we do not throw an error! */
    }
609

Uwe Schulzweida's avatar
Uwe Schulzweida committed
610
611
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612
613
614
615
616
      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
617
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
620
621
622
623
624
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
625
      long processID;
626
627
628
      status = grib_get_long(gh, "generatingProcessIdentifier", &processID);
      if ( status == 0 )
	{
629
630
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
	  modelID = modelInq(varInqInst(varID), (int)processID, NULL);
631
	  if ( modelID == CDI_UNDEFID )
632
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
633
634
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635
636
637
638
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
643
644
645
646
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
647

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
651
652
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
654
655
656
657
658
659
	}
    }

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

  if ( CDI_Debug )
660
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
}
663
#endif
664

665
static
666
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype, char *name)
667
668
{
  compvar2_t compVar;
669
670
671
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
672

673
674
675
676
677
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
678
679
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
680
681
682
683
684

  return (compVar);
}

static
685
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
686
687
{
  compvar2_t compVar0;
688
  size_t maxlen = sizeof(compVar.name);
689

690
691
692
693
694
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
695
  memcpy(compVar0.name, record.varname, maxlen);
696

697
698
699
700
701
702
  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;
    }

703
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
704
705
706

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
707
708
#endif

709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
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
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
static void ensureBufferSize(size_t requiredSize, size_t* curSize, unsigned char** buffer) {
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
      *buffer = realloc(*buffer, *curSize);
    }
}

grib_handle* gribapiGetDiskRepresentation(long recsize, size_t* buffersize, unsigned char** gribbuffer, int* outDatatype, int* outCompressionType, long* outUnzipsize)
{
  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
    {
      if( gribGetZip(recsize, *gribbuffer, outUnzipsize) > 0 )
        {
          *outCompressionType = COMPRESS_SZIP;
          ensureBufferSize(*outUnzipsize + 100, buffersize, gribbuffer);
        }
      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;
}

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

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

802
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
803
804
805
806
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
807
  size_t buffersize = 0;
808
809
  DateTime datetime0 = { .date = 10101, .time = 0 };
  int nrecs_scanned = 0;        //Only used for debug output.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810
  int warn_time = TRUE;
811
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812
813
814
815
816
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
  grib_handle *gh = NULL;

  streamptr->curTsID = 0;

817
818
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
819
820

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

823
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824

825
  int nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
827
  while ( TRUE )
    {
828
829
      int level1 = 0, level2 = 0;
      size_t recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
831
832
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
833
834
835
836
837
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
838

839
840
      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
841
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
842

843
844
845
      int datatype, comptype;
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
846

847
      nrecs_scanned++;
848
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
849

850
851
852
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit);
853

854
      char varname[256];
855
856
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
857

858
859
      int tsteptype = gribapiGetTsteptype(gh);
      int vdate = 0, vtime = 0;
860
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
861
      DateTime datetime = { .date = vdate, .time = vtime };
862
      /*
863
      printf("%d %d %d\n", vdate, vtime, leveltype1);
864
      */
865
866

      if( datetime0.date == 10101 && datetime0.time == 0 )
867
        {
868
869
870
871
872
873
874
875
876
          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);
            }
877
        }
878
879

      if(nrecs)
880
        {
881
882
          checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, varname), &datetime, &datetime0);
          if(result == CHECKTIME_STOP)
883
            {
884
              break;
885
            }
886
887
888
889
890
891
892
893
894
895
896
          else if(result == CHECKTIME_SKIP)
            {
              gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              continue;
            }
          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);
897
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898
899
      /*
      if ( ISEC1_AvgNum )
900
901
902
903
904
905
906
907
908
909
910
911
        {
          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
912
913
914
915
      */
      nrecs++;

      if ( CDI_Debug )
916
917
918
919
        {
          char paramstr[32];
          cdiParamToString(param, paramstr, sizeof(paramstr));
          Message("%4d %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%d vtime=%d",
920
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
921
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
922

923
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
924
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);
925
926
927

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928
929
    }

930
931
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
932
933
  streamptr->rtsteps = 1;

934
935
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

936
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
937

938
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
939
940
941
942
943
944
945
946
947
948
949
950
951
952
  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;
    }

953
954
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
955

956
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
  vlistDefTaxis(vlistID, taxisID);

959
960
  int nrecords = streamptr->tsteps[0].nallrecs;
  if ( nrecords > streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
961
962
963
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
964
      (record_t *) realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
965
966
967
968
    }

  streamptr->tsteps[0].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[0].nrecs = nrecords;
969
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
970
971
972
973
974
975
976
    streamptr->tsteps[0].recIDs[recID] = recID;

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

  if ( streamptr->ntsteps == -1 )
    {
977
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
978
      if ( tsID != streamptr->rtsteps )
979
        Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
980
981
982
983
984
985
986
987

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

  if ( streamptr->ntsteps == 1 )
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
988
989
990
991
992
993
994
        {
          streamptr->ntsteps = 0;
          for ( int varID = 0; varID < streamptr->nvars; varID++ )
            {
              vlistDefVarTsteptype(vlistID, varID, TSTEP_CONSTANT);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
996
    }
#else
997
  (void)streamptr;
998
  Error("GRIB_API support not compiled in!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
999
#endif
1000
1001

  return (0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
1003
1004
}


1005
#ifdef HAVE_LIBGRIB_API
1006
int gribapiScanTimestep2(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1007
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1008
  int rstatus = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
1010
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
1011
  size_t buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1012
  int fileID;
1013
  DateTime datetime0;
1014
  // int gridID;
1015
  int recID;
1016
  //  int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1017
1018
1019
1020
  grib_handle *gh = NULL;

  streamptr->curTsID = 1;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1021
  fileID  = streamptr->fileID;
1022
1023
  int vlistID = streamptr->vlistID;
  int taxisID = vlistInqTaxis(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
1026

  gribbuffer = (unsigned char *) streamptr->record->buffer;
  buffersize = streamptr->record->buffersize;
1027

1028
  int tsID = streamptr->rtsteps;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1029
  if ( tsID != 1 )
1030
    Error("Internal problem! unexpected timestep %d", tsID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1031

1032
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
1034
1035

  fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);

1036
  cdi_create_records(streamptr, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1037

1038
  int nrecords = streamptr->tsteps[tsID].nallrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1039
1040
1041
1042
  streamptr->tsteps[1].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[1].nrecs = 0;
  for ( recID = 0; recID < nrecords; recID++ )
    streamptr->tsteps[1].recIDs[recID] = -1;
1043

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1044
1045
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1046
1047
      streamptr->tsteps[tsID].records[recID].position = streamptr->tsteps[0].records[recID].position;
      streamptr->tsteps[tsID].records[recID].size     = streamptr->tsteps[0].records[recID].size;