stream_gribapi.c 110 KB
Newer Older
1
#ifdef  HAVE_CONFIG_H
2
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3
4
#endif

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


21
22
#include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
#include "gribapi.h"
23

24
#include <grib_api.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25

26
extern int CDI_inventory_mode;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27

28
static const var_tile_t dummy_tiles = { 0, -1, -1, -1, -1, -1 };
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
  int perturbationNumber;
38
  size_t gridsize;
39
  char name[32];
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
{
48
  return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
}

51
static
52
int getTimeunits(long unitsOfTime)
53
54
{
  int timeunits = -1;
55
56
57
58
59
60
61
62
63
64
65
66
67

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

68
  return timeunits;
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
}

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

90
  return factor;
91
92
93
94
95
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
96
  long unitsOfTime = -1;
97
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
98

99
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
100

101
  return getTimeunits(unitsOfTime);
102
103
}

104
static
105
void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep)
106
{
107
  int timeunits2 = timeunits;
108
109
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
110
111
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
112

113
  long lpar;
114
115
116
117
118
119
120
121
  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);
    }
122

123
124
  *endStep = *startStep;
  status = grib_get_long(gh, "endStep", &lpar);
125
  if ( status == 0 )
126
127
    *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));
128
129
}

130
static
131
void gribapiGetDataDateTime(grib_handle *gh, int64_t *datadate, int *datatime)
132
133
134
{
  long lpar;
  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
135
  *datadate = (int64_t) lpar;
136
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);  //FIXME: This looses the seconds in GRIB2 files.
137
138
139
  *datatime = (int) lpar*100;
}

140
static
141
void gribapiSetDataDateTime(grib_handle *gh, int64_t datadate, int datatime)
142
{
143
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", (long)datadate), 0);
144
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
145
146
}

147
static
148
int gribapiGetValidityDateTime(grib_handle *gh, int64_t *vdate, int *vtime)
149
{
150
  int tstepRange = 0;
151

152
  long sigofrtime = 3;
153
  if ( gribEditionNumber(gh) > 1 )
154
    GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
155
  else
156
    GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
157

158
  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())`.
159
    {
160
      gribapiGetDataDateTime(gh, vdate, vtime);
161
162
163
    }
  else
    {
164
165
      int64_t rdate;
      int rtime;
166
167
      gribapiGetDataDateTime(gh, &rdate, &rtime);

168
169
      int timeUnits = gribapiGetTimeUnits(gh);
      int startStep = 0, endStep = 0;
170
      gribapiGetSteps(gh, timeUnits, &startStep, &endStep);
171

172
      int range = endStep - startStep;
173
174
175
176
177
178
      if ( range > 0 )
	{
	  if ( startStep == 0 ) tstepRange = -1;
	  else                  tstepRange =  1;
	}

179
180
181
182
183
      {
	int ryear, rmonth, rday, rhour, rminute, rsecond;
	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

184
185
        if ( rday > 0 )
          {
186
187
            static bool lprint = true;
            extern int grib_calendar;
188
189
            int64_t julday;
            int secofday;
190
191
            encode_caldaysec(grib_calendar, ryear, rmonth, rday, rhour, rminute, rsecond, &julday, &secofday);

192
193
            int64_t time_period = endStep;
            int64_t addsec = 0;
194
195
196
197
198
199
200
201
202
203
204
205
206
            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);
207
                    lprint = false;
208
209
210
211
212
213
214
215
                  }
                break;
              }

            julday_add_seconds(addsec, &julday, &secofday);

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

217
218
219
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
220
    }
221

222
  return tstepRange;
223
224
}

225
static
226
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
227
{
228
  *leveltype = 0;
229
230
231
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
232

233
  long lpar;
234
  if ( !grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar) )       //1 byte
235
    {
236
      *leveltype = (int) lpar;
237

238
      switch (*leveltype)
239
240
241
242
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
243
	  { *lbounds = 1; break; }
244
	}
245

246
247
248
      if ( *lbounds )
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);  //1 byte
249
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
250
251
	  *level1 = (int)lpar;
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);       //1 byte
252
253
          if (lpar == GRIB_MISSING_LONG) lpar = 0;
 	  *level2 = (int)lpar;
254
255
	}
      else
256
	{
257
          double dlevel;
258
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0); //2 byte
259
	  if ( *leveltype == GRIB1_LTYPE_ISOBARIC ) dlevel *= 100;
260
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
261
	  if ( *leveltype == GRIB1_LTYPE_99 || *leveltype == GRIB1_LTYPE_ISOBARIC_PA ) *leveltype = GRIB1_LTYPE_ISOBARIC;
262

263
264
265
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
266
    }
267
268
}

269
270
271
static
double grib2ScaleFactor(long factor)
{
272
273
274
  switch(factor)
    {
      case GRIB_MISSING_LONG: return 1;
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
      case -9: return 1000000000;
      case -8: return 100000000;
      case -7: return 10000000;
      case -6: return 1000000;
      case -5: return 100000;
      case -4: return 10000;
      case -3: return 1000;
      case -2: return 100;
      case -1: return 10;
      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;
294
295
296
297
298
299
300
301
      default: return 0;
    }
}

static
int calcLevel(int level_sf, long factor, long level)
{
  double result = 0;
302
303
  if (level != GRIB_MISSING_LONG) result = (double)level*grib2ScaleFactor(factor);
  if (level_sf) result *= level_sf;
304
  return (int)result;
305
306
}

307
static
308
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1,
309
                   int *level2, int *level_sf, int *level_unit)
310
{
311
312
  *leveltype1 = 0;
  *leveltype2 = -1;
313
314
315
316
317
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
318

319
320
  long lpar;
  int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
321
  if ( status == 0 )
322
    {
323
      *leveltype1 = (int) lpar;
324

325
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
326
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
327
      if ( status == 0 ) *leveltype2 = (int)lpar;
328

329
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
330
      switch(*leveltype1)
331
        {
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
          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;
350
        }
351

352
      long factor, llevel;
353
354
355
      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);
356

357
358
359
360
361
      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);
362
        }
363
364
    }
}
365

366
static
367
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)
368
369
370
371
372
373
374
375
376
377
378
{
  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);
379
380

      /* read in tiles attributes (if there are any) */
381
      tiles->tileindex = (int)gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], 0);
Thomas Jahns's avatar
Thomas Jahns committed
382
383
384
385
386
      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);
387
    }
388
389
}

390
391
392
393
394
static
void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length)
{
  string[0] = 0;

395
396
397
398
399
400
  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);
    }
401
402
403
404
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
static
int gribapiGetEnsembleInfo(grib_handle *gh, long *typeOfEnsembleForecast, long *numberOfForecastsInEnsemble, long *perturbationNumber)
{
  int status = 0;
  if (grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast) == 0)
    {
      GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble), 0);
      GRIB_CHECK(grib_get_long(gh, "perturbationNumber", perturbationNumber), 0);
      if (*perturbationNumber > 0) status = 1;
    }

  if (status == 0)
    {
      *typeOfEnsembleForecast = 0;
      *perturbationNumber = 0;
      *numberOfForecastsInEnsemble = 0;
    }

  return status;
}

426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
static
void gribapiGetKeys(grib_handle *gh, int varID)
{
  long tablesVersion = 0;
  if ( grib_get_long(gh, "tablesVersion", &tablesVersion) == 0 )
    varDefKeyInt(varID, CDI_KEY_TABLESVERSION, (int) tablesVersion);

  long localTablesVersion = 0;
  if ( grib_get_long(gh, "localTablesVersion", &localTablesVersion) == 0 )
    varDefKeyInt(varID, CDI_KEY_LOCALTABLESVERSION, (int) localTablesVersion);

  long typeOfGeneratingProcess = 0;
  if ( grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFGENERATINGPROCESS, (int) typeOfGeneratingProcess);

  long productDefinitionTemplate = 0;
  if ( grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate) == 0 )
    varDefKeyInt(varID, CDI_KEY_PRODUCTDEFINITIONTEMPLATE, (int) productDefinitionTemplate);

  long typeOfProcessedData = 0;
  if ( grib_get_long(gh, "typeOfProcessedData", &typeOfProcessedData) == 0 )
    varDefKeyInt(varID, CDI_KEY_TYPEOFPROCESSEDDATA, (int) typeOfProcessedData);

  long shapeOfTheEarth = 0;
  if ( grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth) == 0 )
    varDefKeyInt(varID, CDI_KEY_SHAPEOFTHEEARTH, (int) shapeOfTheEarth);

  long backgroundProcess = 0;
  if ( grib_get_long(gh, "backgroundProcess", &backgroundProcess) == 0 )
    varDefKeyInt(varID, CDI_KEY_BACKGROUNDPROCESS, (int) backgroundProcess);
456
457
458
459
460

  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
  */
461
462
463
  long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0;
  gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber);
  if ( perturbationNumber > 0 )
464
    {
465
466
467
      varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast);
      varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble);
      varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber);
468
    }
469

470
  long section2Length;
471
472
  int status = grib_get_long(gh, "section2Length", &section2Length);
  if (status == 0 && section2Length > 0)
473
    {
474
      long grib2LocalSectionNumber;
475
      long mpimType, mpimClass, mpimUser;
476
477
      status = grib_get_long(gh, "grib2LocalSectionNumber", &grib2LocalSectionNumber);
      if (status == 0)
478
        {
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
          size_t section2PaddingLength = 0;
          status = grib_get_size(gh, "section2Padding", &section2PaddingLength);
          if (status == 0 && section2PaddingLength > 0)
            {
              varDefKeyInt(varID, CDI_KEY_GRIB2LOCALSECTIONNUMBER, (int) grib2LocalSectionNumber);
              varDefKeyInt(varID, CDI_KEY_SECTION2PADDINGLENGTH, (int) section2PaddingLength);
              unsigned char *section2Padding = (unsigned char*) Malloc(section2PaddingLength);
              grib_get_bytes(gh, "section2Padding", section2Padding, &section2PaddingLength);
              varDefKeyBytes(varID, CDI_KEY_SECTION2PADDING, section2Padding, (int)section2PaddingLength);
              Free(section2Padding);
            }
          else if ( grib_get_long(gh, "mpimType", &mpimType) == 0 &&
                    grib_get_long(gh, "mpimClass", &mpimClass) == 0 &&
                    grib_get_long(gh, "mpimUser", &mpimUser) == 0)
            {
              varDefKeyInt(varID, CDI_KEY_MPIMTYPE, mpimType);
              varDefKeyInt(varID, CDI_KEY_MPIMCLASS, mpimClass);
              varDefKeyInt(varID, CDI_KEY_MPIMUSER, mpimUser);

              size_t revNumLen = 20;
              unsigned char revNumber[revNumLen];
              if ( grib_get_bytes(gh, "revNumber", revNumber, &revNumLen) == 0 )
                varDefKeyBytes(varID, CDI_KEY_REVNUMBER, revNumber, (int)revNumLen);

              long revStatus;
              grib_get_long(gh, "revStatus", &revStatus);
              varDefKeyInt(varID, CDI_KEY_REVSTATUS, revStatus);
            }
507
        }
508
    }
509
510
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
static
void gribapiDefProjRLL(grib_handle *gh, int gridID)
{
  double xpole = 0, ypole = 0, angle = 0;
  grib_get_double(gh, "latitudeOfSouthernPoleInDegrees",  &ypole);
  grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
  grib_get_double(gh, "angleOfRotation", &angle);
  xpole -= 180;
  if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
  if ( fabs(angle) > 0 ) angle = -angle;

  gridDefParamRLL(gridID, xpole, ypole, angle);
}

static
void gribapiDefProjLCC(grib_handle *gh, int gridID)
{
  double a = 6367470., rf = 0;
  long earthIsOblate;
  grib_get_long(gh, "earthIsOblate", &earthIsOblate);
  if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
  double lon_0, lat_1, lat_2, xval_0, yval_0;
  long projflag = 0;
  grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
  grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
  grib_get_double(gh, "LoVInDegrees", &lon_0);
  grib_get_double(gh, "Latin1InDegrees", &lat_1);
  grib_get_double(gh, "Latin2InDegrees", &lat_2);
  grib_get_long(gh, "projectionCentreFlag", &projflag);
  bool lsouth = gribbyte_get_bit((int)projflag, 1);
  if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }

  double lat_0 = lat_2;
544
545
  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547
548
549

  if ( proj_lonlat_to_lcc_func )
    {
      x_0 = xval_0; y_0 = yval_0;
550
551
      proj_lonlat_to_lcc_func(CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
      if ( IS_NOT_EQUAL(x_0, CDI_grid_missval) && IS_NOT_EQUAL(y_0, CDI_grid_missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
553
554
        { x_0 = -x_0; y_0 = -y_0; }
    }

555
  gridDefParamLCC(gridID, CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
557
558
559
560
561
}


static
void gribapiDefProjSTERE(grib_handle *gh, int gridID)
{
562
563
564
565
566
  long shapeOfTheEarth;
  grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
  double a = (shapeOfTheEarth == 6) ? 6371229 : 6367470;

  double lon_0, lat_ts, xval_0, yval_0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567
568
  grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
  grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
569
570
571
572
  grib_get_double(gh, "LaDInDegrees", &lat_ts);
  grib_get_double(gh, "orientationOfTheGridInDegrees", &lon_0);

  long southPoleOnProjectionPlane;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
573
  grib_get_long(gh, "southPoleOnProjectionPlane", &southPoleOnProjectionPlane);
574
  double lat_0 = southPoleOnProjectionPlane ? -90 : 90;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575

576
577
  double x_0 = CDI_grid_missval;
  double y_0 = CDI_grid_missval;
578
579

  if ( proj_lonlat_to_stere_func )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
581
    {
      x_0 = xval_0; y_0 = yval_0;
582
583
      proj_lonlat_to_stere_func(CDI_grid_missval, lon_0, lat_ts, lat_0, a, (size_t)1, &x_0, &y_0);
      if ( IS_NOT_EQUAL(x_0, CDI_grid_missval) && IS_NOT_EQUAL(y_0, CDI_grid_missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
584
585
        { x_0 = -x_0; y_0 = -y_0; }
    }
586

587
  gridDefParamSTERE(gridID, CDI_grid_missval, lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
589
}

590
static
591
void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
592
                      size_t recsize, off_t position, int datatype, int comptype, const char *varname,
593
                      int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, int level_unit,
594
                      const var_tile_t *tiles, int perturbationNumber, int lread_additional_keys)
595
{
596
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
597

598
599
600
  int vlistID = streamptr->vlistID;
  int tsID    = streamptr->curTsID;
  int recID   = recordNewEntry(streamptr, tsID);
601
  record_t *record = &streamptr->tsteps[tsID].records[recID];
602

603
  int tsteptype = gribapiGetTsteptype(gh);
604
  // numavg  = ISEC1_AvgNum;
605
  int numavg = 0;
606

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

Thomas Jahns's avatar
Thomas Jahns committed
609
610
611
612
613
614
  record->size      = recsize;
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
  record->ltype     = leveltype1;
Thomas Jahns's avatar
Thomas Jahns committed
615
  record->tsteptype = (short)tsteptype;
616
  record->gridsize  = gribapiGetGridsize(gh);
Thomas Jahns's avatar
Thomas Jahns committed
617
  record->tiles = tiles ? *tiles : dummy_tiles;
618

619
620
  record->perturbationNumber = perturbationNumber;

621
622
  strncpy(record->varname, varname, sizeof(record->varname)-1);
  record->varname[sizeof(record->varname) - 1] = 0;
623

624
  grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
625
  gribapiGetGrid(gh, grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
626

Uwe Schulzweida's avatar
Uwe Schulzweida committed
627
  struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
628
  int gridID = gridAdded.Id;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
629
630
631
632
  if      ( !gridAdded.isNew ) Free(grid);
  else if ( grid->projtype == CDI_PROJ_RLL ) gribapiDefProjRLL(gh, gridID);
  else if ( grid->projtype == CDI_PROJ_LCC ) gribapiDefProjLCC(gh, gridID);
  else if ( grid->projtype == CDI_PROJ_STERE ) gribapiDefProjSTERE(gh, gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633

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

636
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637
    {
638
639
640
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
641
        long lpar;
642
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
Thomas Jahns's avatar
Thomas Jahns committed
643
        /* FIXME: assert(lpar >= 0) */
Thomas Jahns's avatar
Thomas Jahns committed
644
        size_t vctsize = (size_t)lpar;
645
646
        if ( vctsize > 0 )
          {
Thomas Jahns's avatar
Thomas Jahns committed
647
648
            double *vctptr = (double *) Malloc(vctsize*sizeof(double));
            size_t dummy = vctsize;
649
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
Thomas Jahns's avatar
Thomas Jahns committed
650
            varDefVCT(vctsize, vctptr);
651
            Free(vctptr);
652
653
654
655
656
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
657
        unsigned char uuid[CDI_UUID_SIZE];
658
        long lpar;
659
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
660
        if ( lpar != 6 ) fprintf(stderr, "Warning ...\n");
661
662
663
664
        GRIB_CHECK(grib_get_long(gh, "nlev", &lpar), 0);
        int nhlev = (int)lpar;
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &lpar), 0);
        int nvgrid = (int)lpar;
665
        size_t len = (size_t)CDI_UUID_SIZE;
666
667
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
668
        varDefZAxisReference(nhlev, nvgrid, uuid);
669
670
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672

673
674
  // if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = CDI_DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
675

676
  stdname[0] = 0;
677
678
679
  longname[0] = 0;
  units[0] = 0;

680
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
681
    {
682
      size_t vlen = CDI_MAX_NAME;
683
      gribapiGetString(gh, "name", longname, vlen);
684
      vlen = CDI_MAX_NAME;
685
      gribapiGetString(gh, "units", units, vlen);
Thomas Jahns's avatar
Thomas Jahns committed
686
687
688
689
      vlen = CDI_MAX_NAME;
      int status = grib_get_string(gh, "cfName", stdname, &vlen);
      if ( status != 0 || vlen <= 1 || strncmp(stdname, "unknown", 7) == 0 )
        stdname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
690
    }
691
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
692

693
  /* add the previously read record data to the (intermediate) list of records */
694
  int levelID = 0;
695
  int tile_index = 0, varID;
696
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
697
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype1, leveltype2,
698
	       varname, stdname, longname, units, tiles, &tile_index, perturbationNumber);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699

Thomas Jahns's avatar
Thomas Jahns committed
700
701
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
702

703
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
704

705
  gribapiGetKeys(gh, varID);
706

707
  if (lread_additional_keys)
708
709
710
711
712
713
714
715
716
717
718
719
    {
      long   lval;
      double dval;
      for ( int 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]);
        }
    }
720

Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
722
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
723
724
725
      long center, subcenter;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
726
      int instID = institutInq((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729
730
731
732
733
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
734
      long processID;
735
      if ( grib_get_long(gh, "generatingProcessIdentifier", &processID) == 0 )
736
	{
737
          /* FIXME: assert(processID >= INT_MIN && processID <= INT_MAX) */
738
	  int modelID = modelInq(varInqInst(varID), (int)processID, NULL);
739
	  if ( modelID == CDI_UNDEFID )
740
	    modelID = modelDef(varInqInst(varID), (int)processID, NULL);
741
742
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
743
744
745
746
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
      int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
748
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749

Uwe Schulzweida's avatar
Uwe Schulzweida committed
750
751
752
      if ( pdis == 255 )
	{
	  int tabnum = pcat;
753
	  int tableID = tableInq(varInqModel(varID), tabnum, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754
755
756
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757
758
759
760
761
762
763
	}
    }

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

  if ( CDI_Debug )
764
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
765
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
766
}
767

768
static
769
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, int tsteptype, long perturbationNumber,
770
                         size_t gridsize, char *name, var_tile_t tiles_data)
771
772
{
  compvar2_t compVar;
773
  memset(&compVar, 0, sizeof(compvar2_t));
774
775
776
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
777

778
779
780
781
782
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
783
  compVar.perturbationNumber = (int)perturbationNumber;
784
  compVar.gridsize  = gridsize;
785
  //memset(compVar.name, 0, maxlen);
786
  memcpy(compVar.name, name, len);
787
  compVar.tiles = tiles_data;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
788

789
  return compVar;
790
791
792
}

static
793
int gribapiVarCompare(compvar2_t compVar, record_t record, int flag)
794
795
{
  compvar2_t compVar0;
796
  memset(&compVar0, 0, sizeof(compvar2_t));
797
798
799
800
801
  compVar0.param     = record.param;
  compVar0.level1    = record.ilevel;
  compVar0.level2    = record.ilevel2;
  compVar0.ltype     = record.ltype;
  compVar0.tsteptype = record.tsteptype;
802
  compVar0.perturbationNumber = record.perturbationNumber;
803
  compVar0.gridsize  = record.gridsize;
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
804
  memcpy(compVar0.name, record.varname, sizeof(compVar.name));
805

806
807
808
809
810
811
  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;
    }

812
813
  compVar0.tiles = record.tiles;

Thomas Jahns's avatar
Thomas Jahns committed
814
  return memcmp(&compVar0, &compVar, sizeof(compvar2_t));
815
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816

817
static
818
grib_handle *gribapiGetDiskRepresentation(size_t recsize, size_t *buffersize, void **gribbuffer, int *outDatatype, int *outCompressionType)
819
{
820
  const int gribversion = (int)((char*)*gribbuffer)[7];
821

822
  if ( gribversion <= 1 ) *outCompressionType = grbDecompress(recsize, buffersize, gribbuffer);
823

824
  grib_handle *gh = grib_handle_new_from_message(NULL, *gribbuffer, recsize);
825
826
827
828

  bool lieee = false;

  if ( gribversion > 1 )
829
830
831
832
833
834
835
    {
      size_t len = 256;
      char typeOfPacking[256];

      if ( grib_get_string(gh, "packingType", typeOfPacking, &len) == 0 )
        {
          // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
836
          if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) *outCompressionType = CDI_COMPRESS_JPEG;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
837
          else if ( strncmp(typeOfPacking, "grid_ccsds", len) == 0 ) *outCompressionType = CDI_COMPRESS_AEC;
838
          else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = true;
839
840
841
842
843
        }
    }

  if ( lieee )
    {
844
      *outDatatype = CDI_DATATYPE_FLT64;
845
846
      long precision;
      int status = grib_get_long(gh, "precision", &precision);
847
      if ( status == 0 && precision == 1 ) *outDatatype = CDI_DATATYPE_FLT32;
848
849
850
    }
  else
    {
851
      *outDatatype = CDI_DATATYPE_PACK;
852
853
854
855
856
857
      long bitsPerValue;
      if ( grib_get_long(gh, "bitsPerValue", &bitsPerValue) == 0 )
        {
          if ( bitsPerValue > 0 && bitsPerValue <= 32 ) *outDatatype = (int)bitsPerValue;
        }
    }
858

859
860
861
862
863
  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
864
  // First determine whether the current record exists already.
865
866
867
868
869
870
871
  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
872
  // Then we need to know whether the verification time is consistent.
873
874
  int consistentTime = !memcmp(verificationTime, expectedVTime, sizeof(*verificationTime));

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
875
  // Finally, we make a decision.
876
  if ( CDI_inventory_mode == 1 )
877
878
879
880
881
882
883
884
885
    {
      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
886

887
888
889
890
891
892
893
894
895
896
  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)
897

898
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
900
{
  off_t recpos = 0;
901
  void *gribbuffer = NULL;
902
  size_t buffersize = 0;
903
904
  DateTime datetime0 = { .date = 10101, .time = 0 };
  int nrecs_scanned = 0;        //Only used for debug output.
905
906
  bool warn_time = true;
  // bool warn_numavg = true;
907
908
  int64_t rdate = 0;
  int rtime = 0, tunit = 0, fcast = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
909
910
911
912
  grib_handle *gh = NULL;

  streamptr->curTsID = 0;

913
914
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
915
916

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

919
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
920

Thomas Jahns's avatar
Thomas Jahns committed
921
  unsigned nrecs = 0;
922
  while ( true )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
923
    {
924
      int level1 = 0, level2 = 0;
925
926
      size_t recsize = gribGetSize(fileID);
      recpos = fileGetPos(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
927
928

      if ( recsize == 0 )