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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
static
67
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
71
  int zaxistype = ZAXIS_GENERIC;

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

  return (zaxistype);
}

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

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

100
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
  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;
129
  long unitsOfTime = -1;
130

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

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

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

  return (timeunits);
}

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

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

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

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

  return (endStep);
}

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

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

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

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

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

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

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

      range = endStep - startStep;

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

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

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

231
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
        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);
          }
258

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

  return (tstepRange);
265
266
}

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

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

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

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

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

309
310
311
static
double grib2ScaleFactor(long factor)
{
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
  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;
336
337
}

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

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

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

359
      *leveltype1 = (int) lpar;
360

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

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

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

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

401
402
403
404
405
406
407
408
409
410
411
412
413
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);
414
    }
415
416
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  return (compVar);
}

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

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

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

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

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

708
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
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)
800

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

  streamptr->curTsID = 0;

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

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

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

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

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

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

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

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

849
850
851
      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);
852

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

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

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

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

      if ( CDI_Debug )
915
916
917
918
        {
          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",
919
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
920
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921

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

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

929
930
  if ( gh ) grib_handle_delete(gh);

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

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

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

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

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

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

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

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

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

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

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

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

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


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

  streamptr->curTsID = 1;

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

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

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

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

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

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

1037
  int nrecords = streamptr->tsteps[tsID].nallrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1038
1039
1040
1041
  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;
1042

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