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

5
#if  defined  (HAVE_LIBGRIB_API)
6
#include <limits.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
8
#include <stdio.h>

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


23
#  include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
#  include "gribapi.h"
25
26

#  include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
28
29

extern int cdiInventoryMode;

30
static const var_tile_t dummy_tiles = { -1, -1, -1, -1, -1, -1 };
31

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

  var_tile_t tiles;

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


static
46
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
    {
52
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
54
55
    }
  else
    {
56
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
58
59
60
61
    }

  return (zaxistype);
}

62
static
63
int getTimeunits(long unitsOfTime)
64
65
{
  int timeunits = -1;
66
67
68
69
70
71
72
73
74
75
76
77
78

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

79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  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;
108
  long unitsOfTime = -1;
109

110
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
111

112
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
113

114
  timeunits = getTimeunits(unitsOfTime);
115
116
117
118

  return (timeunits);
}

119
static
120
void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
121
{
122
  int timeunits2 = timeunits;
123
124
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
125
126
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
127

128
  long lpar;
129
130
131
132
133
134
135
136
  status = grib_get_long(gh, "forecastTime", &lpar);
  if ( status == 0 ) *startStep = (int) lpar;
  else
    {
      status = grib_get_long(gh, "startStep", &lpar);
      if ( status == 0 )
        *startStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
    }
137

138
139
  *endStep = *startStep;
  status = grib_get_long(gh, "endStep", &lpar);
140
  if ( status == 0 )
141
142
    *endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
  // printf("%d %d %d %d %d %g\n", *startStep, *endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));
143
144
}

145
146
147
148
149
150
151
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
152
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
153
154
155
  *datatime = (int) lpar*100;
}

156
157
158
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
159
160
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
161
162
}

163
static
164
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
165
{
166
  int rdate, rtime;
167
  int timeUnits, startStep = 0, endStep;
168
169
  int tstepRange = 0;
  int range;
170
171
  long sigofrtime = 3;

172
  if ( gribEditionNumber(gh) > 1 )
173
174
175
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
176
177
178
179
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
180

181
  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())`.
182
    {
183
      gribapiGetDataDateTime(gh, vdate, vtime);
184
185
186
    }
  else
    {
187
188
      gribapiGetDataDateTime(gh, &rdate, &rtime);

189
      timeUnits = gribapiGetTimeUnits(gh);
190
      gribapiGetSteps(gh, timeUnits, &startStep, &endStep);
191
192
193
194
195
196
197
198
199

      range = endStep - startStep;

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

200
201
202
203
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
204
205
206
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
207
208
209
210

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

211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
        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);
          }
238

239
240
241
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
242
    }
243
244

  return (tstepRange);
245
246
}

247
static
248
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
249
{
250
  *leveltype = 0;
251
252
253
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
254

255
  long lpar;
256
  if(!grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar))       //1 byte
257
    {
258
      *leveltype = (int) lpar;
259

260
      switch (*leveltype)
261
262
263
264
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
265
	  { *lbounds = 1; break; }
266
	}
267

268
269
270
271
272
273
274
275
      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
276
	{
277
          double dlevel;
278
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
279
280
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
281
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
282

283
284
285
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
286
    }
287
288
}

289
290
291
static
double grib2ScaleFactor(long factor)
{
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
  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
313
  if(level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
314
315
  if(level_sf) result *= level_sf;
  return (int)result;
316
317
}

318
static
319
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
320
                   int *level2, int *level_sf, int *level_unit)
321
{
322
323
324
  int status;
  long lpar;
  long factor;
325

326
327
  *leveltype1 = 0;
  *leveltype2 = -1;
328
329
330
331
332
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
333

334
  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
335
  if ( status == 0 )
336
    {
337
      long llevel;
338

339
      *leveltype1 = (int) lpar;
340

341
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
342
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
343
      if ( status == 0 ) *leveltype2 = (int)lpar;
344

345
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
346
      switch(*leveltype1)
347
        {
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
          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;
366
        }
367

368
369
370
      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);
371

372
373
374
375
376
      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);
377
        }
378
379
    }
}
380

381
static
382
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)
383
384
385
386
387
388
389
390
391
392
393
{
  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);
394
395

      /* read in tiles attributes (if there are any) */
Thomas Jahns's avatar
Thomas Jahns committed
396
397
398
399
400
401
      tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], -1);
      tiles->totalno_of_tileattr_pairs = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TOTALNO_OF_TILEATTR_PAIRS], -1);
      tiles->tileClassification = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILE_CLASSIFICATION], -1);
      tiles->numberOfTiles = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_TILES], -1);
      tiles->numberOfAttributes = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_NUMBER_OF_ATTR], -1);
      tiles->attribute = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEATTRIBUTE], -1);
402
    }
403
404
}

405
406
407
408
409
static
void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length)
{
  string[0] = 0;

410
411
412
413
414
415
  int ret = grib_get_string(gh, key, string, &length);
  if (ret != 0)
    {
      fprintf(stderr, "grib_get_string(gh, \"%s\", ...) failed!\n", key);
      GRIB_CHECK(ret, 0);
    }
416
417
418
419
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

420
static
421
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
422
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
423
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
424
                      const var_tile_t *tiles, int lread_additional_keys)
425
{
426
  int varID;
427
428
429
430
  int levelID = 0;
  grid_t grid;
  long lpar;
  int status;
431
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
432
  size_t vlen;
433
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
434

435
436
437
438
  int vlistID = streamptr->vlistID;
  int tsID    = streamptr->curTsID;
  int recID   = recordNewEntry(streamptr, tsID);
  record_t *record  = &streamptr->tsteps[tsID].records[recID];
439

440
  int tsteptype = gribapiGetTsteptype(gh);
441
  // numavg  = ISEC1_AvgNum;
442
  int numavg  = 0;
443

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

Thomas Jahns's avatar
Thomas Jahns committed
446
447
448
449
450
451
452
453
  record->size      = recsize;
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
  record->ltype     = leveltype1;
  record->tsteptype = tsteptype;
  record->tiles = tiles ? *tiles : dummy_tiles;
454
455
456
457
458
459
460
461
462
463
464
465

  //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));
466
467

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

469
  int gridID = varDefGrid(vlistID, &grid, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470

471
  int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472

473
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
474
    {
475
476
477
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
Thomas Jahns's avatar
Thomas Jahns committed
478
        size_t vctsize;
479
480
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481

482
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
483
484
        /* FIXME: assert(lpar >= 0) */
        vctsize = (size_t)lpar;
485
486
        if ( vctsize > 0 )
          {
487
            vctptr = (double *) Malloc(vctsize*sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
488
            dummy = vctsize;
489
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
490
            varDefVCT(vctsize, vctptr);
491
            Free(vctptr);
492
493
494
495
496
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
497
        unsigned char uuid[CDI_UUID_SIZE];
498
        long ltmp;
499
        long nhlev, nvgrid;
500
501

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
502
        if ( lpar != 6 )
503
504
505
          {
            fprintf(stderr, "Warning ...\n");
          }
506
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
507
        nhlev = ltmp;
508
509
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
510
        size_t len = (size_t)CDI_UUID_SIZE;
511
512
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
513
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
514
515
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
516
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517

518
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
519
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
520

521
  stdname[0] = 0;
522
523
524
  longname[0] = 0;
  units[0] = 0;

525
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526
    {
527
      vlen = CDI_MAX_NAME;
528
      gribapiGetString(gh, "name", longname, vlen);
529
      vlen = CDI_MAX_NAME;
530
      gribapiGetString(gh, "units", units, vlen);
531
532
533
534
535

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

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

Thomas Jahns's avatar
Thomas Jahns committed
547
548
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549

550
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551

552
553
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
554
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
555
  */
556
557
558
559
560
561
562
  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);
    }

563
564
565
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

566
567
568
569
570
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

571
572
573
574
575
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

576
577
578
579
  int    i;
  long   lval;
  double dval;

580
581
582
583
584
585
586
587
588
  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]);
      }
589

Uwe Schulzweida's avatar
Uwe Schulzweida committed
590
591
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592
593
594
595
596
      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
597
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
599
600
601
602
603
604
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
      long processID;
606
607
608
      status = grib_get_long(gh, "generatingProcessIdentifier", &processID);
      if ( status == 0 )
	{
609
610
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
	  modelID = modelInq(varInqInst(varID), (int)processID, NULL);
611
	  if ( modelID == CDI_UNDEFID )
612
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
613
614
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
616
617
618
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
623
624
625
626
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
627

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
630
631
632
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
634
635
636
637
638
639
	}
    }

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

  if ( CDI_Debug )
640
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
641
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
}
643

644
645
static compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, 
                                int tsteptype, char *name, var_tile_t tiles_data)
646
647
{
  compvar2_t compVar;
648
649
650
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
651

652
653
654
655
656
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
657
658
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
659
  compVar.tiles = tiles_data;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
660

661
662
663
664
  return (compVar);
}

static
665
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
666
667
{
  compvar2_t compVar0;
668
669
670
671
672
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
673
  memcpy(compVar0.name, record.varname, sizeof(compVar.name));
674

675
676
677
678
679
680
  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;
    }

681
682
  compVar0.tiles = record.tiles;

683
  int rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));
684
685
686

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

688
689
690
static
void ensureBufferSize(size_t requiredSize, size_t *curSize, void **buffer)
{
691
692
693
  if ( *curSize < requiredSize )
    {
      *curSize = requiredSize;
694
      *buffer = Realloc(*buffer, *curSize);
695
696
697
    }
}

698
static
699
grib_handle *gribapiGetDiskRepresentation(size_t recsize, size_t *buffersize, void **gribbuffer, int *outDatatype, int *outCompressionType, long *outUnzipsize)
700
701
702
{
  int lieee = FALSE;

703
  grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
  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
719
      if( gribGetZip((long)recsize, *gribbuffer, outUnzipsize) > 0 )
720
721
        {
          *outCompressionType = COMPRESS_SZIP;
Thomas Jahns's avatar
Thomas Jahns committed
722
          ensureBufferSize((size_t)*outUnzipsize + 100, buffersize, gribbuffer);
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
        }
      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) {
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
751
  // First determine whether the current record exists already.
752
753
754
755
756
757
758
  int recID = 0;
  for ( ; recID < streamptr->nrecs; recID++ )
    {
      if ( gribapiVarCompare(compVar, streamptr->tsteps[0].records[recID], 1) == 0 ) break;
    }
  int recordExists = recID < streamptr->nrecs;

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
759
  // Then we need to know whether the verification time is consistent.
760
761
  int consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
762
  // Finally, we make a decision.
763
764
765
766
767
768
769
770
771
772
  if ( cdiInventoryMode == 1 )
    {
      if ( recordExists ) return CHECKTIME_STOP;
      if ( !consistentTime ) return CHECKTIME_INCONSISTENT;
    }
  else
    {
      if ( !consistentTime ) return CHECKTIME_STOP;
      if ( recordExists ) return CHECKTIME_SKIP;
    }
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
773

774
775
776
777
778
779
780
781
782
783
  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)
784

785
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
786
787
{
  off_t recpos = 0;
788
  void *gribbuffer = NULL;
789
  size_t buffersize = 0;
790
791
  DateTime datetime0 = { .date = 10101, .time = 0 };
  int nrecs_scanned = 0;        //Only used for debug output.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
792
  int warn_time = TRUE;
793
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
794
795
796
797
798
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
  grib_handle *gh = NULL;

  streamptr->curTsID = 0;

799
800
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801
802

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

805
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
806

Thomas Jahns's avatar
Thomas Jahns committed
807
  unsigned nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
808
809
  while ( TRUE )
    {
810
811
      int level1 = 0, level2 = 0;
      size_t recsize = (size_t)gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812
813
814
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
815
816
817
818
819
        {
          streamptr->ntsteps = 1;
          break;
        }
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
820

821
      size_t readsize = recsize;
822
      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
823
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824

825
      int datatype, comptype = 0;
826
827
      long unzipsize;
      gh = gribapiGetDiskRepresentation(recsize, &buffersize, &gribbuffer, &datatype, &comptype, &unzipsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
828

829
      nrecs_scanned++;
830
      GRIB_CHECK(my_grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831

832
833
      int param = gribapiGetParam(gh);
      int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit;
834
      var_tile_t tiles = dummy_tiles;
835
      gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles);
836

837
      char varname[256];
838
839
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
840

841
      int tsteptype = gribapiGetTsteptype(gh);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
842

843
      int vdate = 0, vtime = 0;
844
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
845
      DateTime datetime = { .date = vdate, .time = vtime };
846
      /*
847
      printf("%d %d %d\n", vdate, vtime, leveltype1);
848
      */
849

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
850
      if ( datetime0.date == 10101 && datetime0.time == 0 )
851
        {
852
          if( datetimeCmp(datetime, datetime0) || !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.
853
854
855
856
857
858
859
860
            {
              datetime0 = datetime;

              gribapiGetDataDateTime(gh, &rdate, &rtime);

              fcast = gribapiTimeIsFC(gh);
              if ( fcast ) tunit = gribapiGetTimeUnits(gh);
            }
861
        }
862

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
863
      if ( nrecs )
864
        {
865
          checkTimeResult result = checkTime(streamptr, gribapiVarSet(param, level1, level2, leveltype1, tsteptype, varname, tiles), &datetime, &datetime0);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
866
          if ( result == CHECKTIME_STOP )
867
            {
868
              break;
869
            }
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
870
          else if ( result == CHECKTIME_SKIP )
871
            {
872
873
              gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              continue;
874
            }
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
875
          else if ( result == CHECKTIME_INCONSISTENT && warn_time )
876
877
878
879
880
            {
              gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, varname, param, level1, level2);
              warn_time = FALSE;
            }
          assert(result == CHECKTIME_OK || result == CHECKTIME_INCONSISTENT);
881
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
882
883
      /*
      if ( ISEC1_AvgNum )
884
885
886
887
888
889
890
891
892
893
894
895
        {
          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
896
897
898
899
      */
      nrecs++;

      if ( CDI_Debug )
900
901
902
        {
          char paramstr[32];
          cdiParamToString(param, paramstr, sizeof(paramstr));
Thomas Jahns's avatar
Thomas Jahns committed
903
          Message("%4u %8d name=%s id=%s ltype=%d lev1=%d lev2=%d vdate=%d vtime=%d",
904
                nrecs, (int)recpos, varname, paramstr, leveltype1, level1, level2, vdate, vtime);
905
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
906

907
908
      var_tile_t *ptiles = NULL;
      if ( memcmp(&tiles, &dummy_tiles, sizeof(var_tile_t)) != 0 ) ptiles = &tiles;
909
      gribapiAddRecord(streamptr, param, gh, recsize, recpos, datatype, comptype, varname,
910
                       leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit, ptiles, 1);
911
912
913

      grib_handle_delete(gh);
      gh = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
914
915
    }

916
917
  if ( gh ) grib_handle_delete(gh);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
918
919
  streamptr->rtsteps = 1;

920
921
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

922
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
923

924
  int taxisID = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
925
926
927
928
929
930
931
932
933
934
935
936
937
938
  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;
    }

939
940
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
941

942
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
943
944
  vlistDefTaxis(vlistID, taxisID);

945
  int nrecords = streamptr->tsteps[0].nallrecs;
946
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
947
948
949
    {
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
950
        (record_t *) Realloc(streamptr->tsteps[0].records, (size_t)nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
951
952
    }

953
  streamptr->tsteps[0].recIDs = (int *) Malloc((size_t)nrecords*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
954
  streamptr->tsteps[0].nrecs = nrecords;