stream_gribapi.c 80.9 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
5
6
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <stdio.h>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
#include "cdi.h"
9
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
#include "file.h"
11
12
#include "gribapi_utilities.h"
#include "stream_grb.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
13
14
15
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
#include "stream_grb.h"
17
#include "calendar.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20


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

#  include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
28
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
  switch(factor)
    {
      case GRIB_MISSING_LONG: return 1;
      case 0: return 1;
      case 1: return 0.1;
      case 2: return 0.01;
      case 3: return 0.001;
      case 4: return 0.0001;
      case 5: return 0.00001;
      case 6: return 0.000001;
      case 7: return 0.0000001;
      case 8: return 0.00000001;
      case 9: return 0.000000001;
      default: return 0;
    }
}

static
int calcLevel(int level_sf, long factor, long level)
{
  double result = 0;
Thomas Jahns's avatar
Thomas Jahns committed
334
  if(level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
335
336
  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
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
Thomas Jahns's avatar
Thomas Jahns committed
491
        size_t vctsize;
492
493
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494

495
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
496
497
        /* FIXME: assert(lpar >= 0) */
        vctsize = (size_t)lpar;
498
499
500
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
501
            dummy = vctsize;
502
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
503
            varDefVCT(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

  return (compVar);
}
683
#endif
684

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

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

699
700
701
702
703
704
  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;
    }

705
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
706
707
708

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

710
711
712
713
714
715
716
717
static void ensureBufferSize(size_t requiredSize, size_t* curSize, unsigned char** buffer) {
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
      *buffer = realloc(*buffer, *curSize);
    }
}

718
static
Thomas Jahns's avatar
Thomas Jahns committed
719
grib_handle* gribapiGetDiskRepresentation(size_t recsize, size_t* buffersize, unsigned char** gribbuffer, int* outDatatype, int* outCompressionType, long* outUnzipsize)
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
{
  int lieee = FALSE;

  grib_handle* gh = grib_handle_new_from_message(NULL, (void *) *gribbuffer, recsize);
  if(gribEditionNumber(gh) > 1)
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = COMPRESS_JPEG;
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = COMPRESS_SZIP;
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = TRUE;
        }
    }
  else
    {
Thomas Jahns's avatar
Thomas Jahns committed
739
      if( gribGetZip((long)recsize, *gribbuffer, outUnzipsize) > 0 )
740
741
        {
          *outCompressionType = COMPRESS_SZIP;
Thomas Jahns's avatar
Thomas Jahns committed
742
          ensureBufferSize((size_t)*outUnzipsize + 100, buffersize, gribbuffer);
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
        }
      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;
}
768
#endif
769

770
#if  defined  (HAVE_LIBGRIB_API)
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
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;
}
797
#endif
798
799
800
801
802
803
804
805

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

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

  streamptr->curTsID = 0;

822
823
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824
825

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

828
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829

Thomas Jahns's avatar
Thomas Jahns committed
830
  unsigned nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831
832
  while ( TRUE )
    {
833
834
      int level1 = 0, level2 = 0;
      size_t recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
835
836
837
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
838
839
840
841
842
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843

844
845
      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
846
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
847

848
      int datatype, comptype = 0;
849
850
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
851

852
      nrecs_scanned++;
853
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
854

855
856
857
      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);
858

859
      char varname[256];
860
861
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
862

863
864
      int tsteptype = gribapiGetTsteptype(gh);
      int vdate = 0, vtime = 0;
865
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
866
      DateTime datetime = { .date = vdate, .time = vtime };
867
      /*
868
      printf("%d %d %d\n", vdate, vtime, leveltype1);
869
      */
870
871

      if( datetime0.date == 10101 && datetime0.time == 0 )
872
        {
873
874
875
876
877
878
879
880
881
          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);
            }
882
        }
883
884

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

      if ( CDI_Debug )
921
922
923
        {
          char paramstr[32];
          cdiParamToString(param, paramstr, sizeof(paramstr));
Thomas Jahns's avatar
Thomas Jahns committed
924
          Message("%4u %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%d vtime=%d",
925
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
926
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
927

928
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
929
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);
930
931
932

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
933
934
    }

935
936
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
937
938
  streamptr->rtsteps = 1;

939
940
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

941
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942

943
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
944
945
946
947
948
949
950
951
952
953
954
955
956
957
  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;
    }

958
959
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960

961
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
963
  vlistDefTaxis(vlistID, taxisID);

964
  int nrecords = streamptr->tsteps[0].nallrecs;
965
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
967
968
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Thomas Jahns's avatar
Thomas Jahns committed
969
        (record_t *)xrealloc(streamptr->tsteps[0].records, (size_t)nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
970
971
    }

Thomas Jahns's avatar
Thomas Jahns committed
972
  streamptr->tsteps[0].recIDs = (int *)xmalloc((size_t)nrecords*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
973
  streamptr->tsteps[0].nrecs = nrecords;
974
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
975
976
977
978
979
980
981
    streamptr->tsteps[0].recIDs[recID] = recID;

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

  if ( streamptr->ntsteps == -1 )
    {
982
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
983
      if ( tsID != streamptr->rtsteps )
984
        Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
986
987
988
989
990
991
992

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

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

  return (0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1007
1008
1009
}


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