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

#include <stdio.h>

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


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

#  include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
28
#endif

extern int cdiInventoryMode;

29
30


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

  var_tile_t tiles;

41
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
43
44


#if  defined  (HAVE_LIBGRIB_API)
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
static
int my_grib_set_double(grib_handle* h, const char* key, double val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_double(\tgrib_handle* h, \"%s\", %f)\n", key, val);

  return grib_set_double(h, key, val);
}

static
int my_grib_set_long(grib_handle* h, const char* key, long val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_long(  \tgrib_handle* h, \"%s\", %ld)\n", key, val);

  return grib_set_long(h, key, val);
}

static
int my_grib_set_string(grib_handle* h, const char* key, const char* val, size_t* length)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_string(\tgrib_handle* h, \"%s\", \"%s\")\n", key, val);

  return grib_set_string(h, key, val, length);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
static
73
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
    {
79
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
    }
  else
    {
83
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
85
86
87
88
    }

  return (zaxistype);
}

89
static
90
int getTimeunits(long unitsOfTime)
91
92
{
  int timeunits = -1;
93
94
95
96
97
98
99
100
101
102
103
104
105

  switch (unitsOfTime)
    {
    case 13:  timeunits = TUNIT_SECOND;  break;
    case  0:  timeunits = TUNIT_MINUTE;  break;
    case  1:  timeunits = TUNIT_HOUR;    break;
    case 10:  timeunits = TUNIT_3HOURS;  break;
    case 11:  timeunits = TUNIT_6HOURS;  break;
    case 12:  timeunits = TUNIT_12HOURS; break;
    case  2:  timeunits = TUNIT_DAY;     break;
    default:  timeunits = TUNIT_HOUR;    break;
    }

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
  return (timeunits);
}

static
double timeunit_factor(int tu1, int tu2)
{
  double factor = 1;

  if ( tu2 == TUNIT_HOUR )
    {
      switch (tu1)
        {
        case TUNIT_SECOND:  factor = 3600;   break;
        case TUNIT_MINUTE:  factor = 60;     break;
        case TUNIT_HOUR:    factor = 1;      break;
        case TUNIT_3HOURS:  factor = 1./3;   break;
        case TUNIT_6HOURS:  factor = 1./6;   break;
        case TUNIT_12HOURS: factor = 1./12;  break;
        case TUNIT_DAY:     factor = 1./24;  break;
        }
    }

  return (factor);
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  int timeunits = -1;
135
  long unitsOfTime = -1;
136

137
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
138

139
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
140

141
  timeunits = getTimeunits(unitsOfTime);
142
143
144
145

  return (timeunits);
}

146
147
148
149
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
150
  int timeunits2 = timeunits;
151

152
153
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
154
155
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
156

157
  long lpar;
158
159
160
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
161
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
162
163
164
165
166
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

167
168
169
170
171
172
173
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
174
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
175
176
177
  *datatime = (int) lpar*100;
}

178
179
180
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
181
182
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
183
184
}

185
static
186
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
187
{
188
  int rdate, rtime;
189
  int timeUnits, startStep = 0, endStep;
190
191
  int tstepRange = 0;
  int range;
192
  int status;
193
  long lpar;
194
195
  long sigofrtime = 3;

196
  if ( gribEditionNumber(gh) > 1 )
197
198
199
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
200
201
202
203
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
204

205
  if ( sigofrtime == 3 )        //XXX: This looks like a bug to me, because timeRangeIndicator == 3 does not seem to have the same meaning as significanceOfReferenceTime == 3. I would recommend replacing this condition with `if(!gribapiTimeIsFC())`.
206
    {
207
      gribapiGetDataDateTime(gh, vdate, vtime);
208
209
210
    }
  else
    {
211
212
      gribapiGetDataDateTime(gh, &rdate, &rtime);

213
214
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
215
      timeUnits = gribapiGetTimeUnits(gh);
216
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
217
218
219
220
221
222
223
224
225

      range = endStep - startStep;

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

226
227
228
229
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
230
231
232
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
233
234
235
236

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

237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        if ( rday > 0 )
          {
            encode_caldaysec(grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);

            addsec = 0;
            switch ( timeUnits )
              {
              case TUNIT_SECOND:  addsec =         time_period; break;
              case TUNIT_MINUTE:  addsec =    60 * time_period; break;
              case TUNIT_HOUR:    addsec =  3600 * time_period; break;
              case TUNIT_3HOURS:  addsec = 10800 * time_period; break;
              case TUNIT_6HOURS:  addsec = 21600 * time_period; break;
              case TUNIT_12HOURS: addsec = 43200 * time_period; break;
              case TUNIT_DAY:     addsec = 86400 * time_period; break;
              default:
                if ( lprint )
                  {
                    Warning("Time unit %d unsupported", timeUnits);
                    lprint = FALSE;
                  }
                break;
              }

            julday_add_seconds(addsec, &julday, &secofday);

            decode_caldaysec(grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
          }
264

265
266
267
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
268
    }
269
270

  return (tstepRange);
271
272
}

273
static
274
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
275
{
276
  *leveltype = 0;
277
278
279
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
280

281
  long lpar;
282
  if(!grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar))       //1 byte
283
    {
284
      *leveltype = (int) lpar;
285

286
      switch (*leveltype)
287
288
289
290
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
291
	  { *lbounds = 1; break; }
292
	}
293

294
295
296
297
298
299
300
301
      if ( *lbounds )
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);  //1 byte
	  *level1 = (int)lpar;
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);       //1 byte
	  *level2 = (int)lpar;
	}
      else
302
	{
303
          double dlevel;
304
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
305
306
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
307
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
308

309
310
311
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
312
    }
313
314
}

315
316
317
static
double grib2ScaleFactor(long factor)
{
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
  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;
342
343
}

344
static
345
346
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
347
{
348
349
350
  int status;
  long lpar;
  long factor;
351

352
353
  *leveltype1 = 0;
  *leveltype2 = -1;
354
355
356
357
358
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
359

360
  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
361
  if ( status == 0 )
362
    {
363
      long llevel;
364

365
      *leveltype1 = (int) lpar;
366

367
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
368
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
369
      if ( status == 0 ) *leveltype2 = (int)lpar;
370

371
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
372
      switch(*leveltype1)
373
        {
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
          case GRIB2_LTYPE_REFERENCE:
            if(*leveltype2 == 1) *lbounds = 0;
            break;

          case GRIB2_LTYPE_LANDDEPTH:
            *level_sf = 1000;
            *level_unit = CDI_UNIT_M;
            break;

          case GRIB2_LTYPE_ISOBARIC:
            *level_sf = 1000;
            *level_unit = CDI_UNIT_PA;
            break;

          case GRIB2_LTYPE_SIGMA:
            *level_sf = 1000;
            *level_unit = 0;
            break;
392
        }
393

394
395
396
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);      //1 byte
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);      //4 byte
      *level1 = calcLevel(*level_sf, factor, llevel);
397

398
399
400
401
402
      if ( *lbounds )
        {
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0); //1 byte
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0); //4 byte
          *level2 = calcLevel(*level_sf, factor, llevel);
403
        }
404
405
    }
}
406

407
static
408
void gribGetLevel(grib_handle *gh, int* leveltype1, int* leveltype2, int* lbounds, int* level1, int* level2, int* level_sf, int* level_unit, var_tile_t* tiles)
409
410
411
412
413
414
415
416
417
418
419
{
  if ( gribEditionNumber(gh) <= 1 )
    {
      grib1GetLevel(gh, leveltype1, lbounds, level1, level2);
      *leveltype2 = -1;
      *level_sf = 0;
      *level_unit = 0;
    }
  else
    {
      grib2GetLevel(gh, leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit);
420
421
422
423
424
425
426
427

      /* read in tiles attributes (if there are any) */
      tiles->tileindex = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_TILEINDEX], -1);
      tiles->totalno_of_tileattr_pairs = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_TOTALNO_OF_TILEATTR_PAIRS], -1);
      tiles->tileClassification = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_TILE_CLASSIFICATION], -1);
      tiles->numberOfTiles = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_TILES], -1);
      tiles->numberOfAttributes = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_ATTR], -1);
      tiles->attribute = gribGetLongDefault(gh, subtypeAttributeName[SUBTYPE_ATT_TILEATTRIBUTE], -1);
428
    }
429
430
}

431
432
433
434
435
436
437
438
439
440
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;
}

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

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

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

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

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

  //FIXME: This may leave the variable name unterminated (which is the behavior that I found in the code).
  //       I don't know precisely how this field is used, so I did not change this behavior to avoid regressions,
  //       but I think that it would be better to at least add a line
  //
  //           record->varname[sizeof(record->varname) - 1] = 0;`
  //
  //       after the `strncpy()` call.
  //
  //       I would consider using strdup() (that requires POSIX-2008 compliance, though), or a similar homebrew approach.
  //       I. e. kick the fixed size array and allocate enough space, whatever that may be.
  strncpy(record->varname, varname, sizeof(record->varname));
493
494

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

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

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

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

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

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

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

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

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

      {
        vlen = CDI_MAX_NAME;
        status = grib_get_string(gh, "cfName", stdname, &vlen);
        if ( status != 0 || vlen <= 1 ) stdname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
        else if ( strncmp(stdname, "unknown", 7) == 0 ) stdname[0] = 0;
565
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
    }
567
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
568

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

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

578
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
579

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

591
592
593
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

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

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

604
605
606
607
  int    i;
  long   lval;
  double dval;

608
609
610
611
612
613
614
615
616
  if (lread_additional_keys)
    for ( i = 0; i < cdiNAdditionalGRIBKeys; i++ )
      {
        /* note: if the key is not defined, we do not throw an error! */
        if ( grib_get_long(gh, cdiAdditionalGRIBKeys[i], &lval) == 0 )
          varDefOptGribInt(varID, tile_index, lval, cdiAdditionalGRIBKeys[i]);
        if ( grib_get_double(gh, cdiAdditionalGRIBKeys[i], &dval) == 0 )
          varDefOptGribDbl(varID, tile_index, dval, cdiAdditionalGRIBKeys[i]);
      }
617

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

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
651
652
653
654
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
655

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
658
659
660
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
662
663
664
665
666
667
	}
    }

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

  if ( CDI_Debug )
668
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
669
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
670
}
671
#endif
672

673
674
static compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, 
                                int tsteptype, char *name, var_tile_t tiles_data)
675
676
{
  compvar2_t compVar;
677
678
679
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
680

681
682
683
684
685
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
686
687
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
688

689
  compVar.tiles = tiles_data;
690
691
692
693
  return (compVar);
}

static
694
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
695
696
{
  compvar2_t compVar0;
697
  size_t maxlen = sizeof(compVar.name);
698

699
700
701
702
703
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
704
  memcpy(compVar0.name, record.varname, maxlen);
705

706
707
708
709
710
711
  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;
    }

712
713
  compVar0.tiles = record.tiles;

714
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
715
716
717

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
718
719
#endif

720
721
722
723
724
725
726
727
static void ensureBufferSize(size_t requiredSize, size_t* curSize, unsigned char** buffer) {
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
      *buffer = realloc(*buffer, *curSize);
    }
}

728
729
#if  defined  (HAVE_LIBGRIB_API)
static
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
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;
}
779
#endif
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814

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)
815

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

  streamptr->curTsID = 0;

831
832
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
834

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

837
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
838

839
  int nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
840
841
  while ( TRUE )
    {
842
843
      int level1 = 0, level2 = 0;
      size_t recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
844
845
846
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
847
848
849
850
851
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
852

853
854
      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
855
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
856

857
      int datatype, comptype = 0;
858
859
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
860

861
      nrecs_scanned++;
862
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
863

864
865
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
866
      var_tile_t tiles = dummy_tiles;
867
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
868

869
      char varname[256];
870
871
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
872

873
874
      int tsteptype = gribapiGetTsteptype(gh);
      int vdate = 0, vtime = 0;
875
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
876
      DateTime datetime = { .date = vdate, .time = vtime };
877
      /*
878
      printf("%d %d %d\n", vdate, vtime, leveltype1);
879
      */
880
881

      if( datetime0.date == 10101 && datetime0.time == 0 )
882
        {
883
884
885
886
887
888
889
890
891
          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);
            }
892
        }
893
894

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

      if ( CDI_Debug )
931
932
933
934
        {
          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",
935
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
936
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
937

938
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
939
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, tiles, 1);
940
941
942

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
943
944
    }

945
946
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
947
948
  streamptr->rtsteps = 1;

949
950
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

951
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
952

953
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
954
955
956
957
958
959
960
961
962
963
964
965
966
967
  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;
    }

968
969
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
970

971
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
973
  vlistDefTaxis(vlistID, taxisID);

974
  int nrecords = streamptr->tsteps[0].nallrecs;
975
  if ( nrecords > streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
976
977
978
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
979
      (record_t *) realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
980
981
982
983
    }

  streamptr->tsteps[0].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[0].nrecs = nrecords;
984
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
986
987
988
989
990
991
    streamptr->tsteps[0].recIDs[recID] = recID;

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

  if ( streamptr->ntsteps == -1 )
    {
992
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
      if ( tsID != streamptr->rtsteps )
994
        Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
996
997
998
999
1000
1001
1002

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

  if ( streamptr->ntsteps == 1 )
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
1003
1004
1005
1006
1007
1008
1009
        {
          streamptr->ntsteps = 0;
          for ( int varID = 0; varID < streamptr->nvars; varID++ )
            {
              vlistDefVarTsteptype(vlistID, varID, TSTEP_CONSTANT);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1010
1011
    }
#else
1012
  (void)streamptr;
1013
  Error("GRIB_API support not compiled in!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1014
#endif
1015
1016

  return (0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1017
1018
1019
}


1020
#ifdef HAVE_LIBGRIB_API
1021
int gribapiScanTimestep2(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1022
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1023
  int rstatus = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
1026
  size_t buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
  int fileID;
1028
  DateTime datetime0;
1029
  // int gridID;
1030
  int recID;
1031
  //  int warn_numavg = TRUE;