stream_gribapi.c 99.4 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
11
12
13
#include "file.h"
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
14
#include "stream_grb.h"
15
#include "calendar.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
17
18


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

extern int cdiInventoryMode;

26
27


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

  var_tile_t tiles;

38
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
41


#if  defined  (HAVE_LIBGRIB_API)
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
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
69
static
70
int gribapiGetGridType(grib_handle *gh)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
{
72
  int gridtype = GRID_GENERIC;
73
  int gribgridtype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
  long lpar;

76
    {
77
78
79
80
      int status;
      status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);

      if ( status ==  0 ) gribgridtype = (int) lpar;
81
82
83

      switch (gribgridtype)
	{
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
	case  GRIB2_GTYPE_LATLON:        { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
	                                   if ( lpar == (long) GRIB_MISSING_LONG ) break;
                                         }
	case  GRIB2_GTYPE_LATLON_ROT:    { gridtype = GRID_LONLAT;    break; }
	case  GRIB2_GTYPE_LCC:           { gridtype = GRID_LCC;       break; }
	case  GRIB2_GTYPE_GAUSSIAN:      { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
	                                   if ( lpar == (long) GRIB_MISSING_LONG )
					     gridtype = GRID_GAUSSIAN_REDUCED;
					   else
					     gridtype = GRID_GAUSSIAN;
				  	   break;
                                         }
	case  GRIB2_GTYPE_SPECTRAL:      { gridtype = GRID_SPECTRAL;  break; }
	case  GRIB2_GTYPE_GME:           { gridtype = GRID_GME;       break; }
	case  GRIB2_GTYPE_UNSTRUCTURED:  { gridtype = GRID_UNSTRUCTURED; break; }
99
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
102
103
104
    }

  return (gridtype);
}

105
static
106
int gribapiGetIsRotated(grib_handle *gh)
107
108
{
  int isRotated = 0;
109
  int gribgridtype = -1;
110
  long lpar;
111
  int status;
112

113
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
114

115
116
117
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
118
119
120
121

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
static
123
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
127
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
    {
129
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132
    }
  else
    {
133
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
135
136
137
138
    }

  return (zaxistype);
}

139
static
140
int getTimeunits(long unitsOfTime)
141
142
{
  int timeunits = -1;
143
144
145
146
147
148
149
150
151
152
153
154
155

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

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  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;
185
  long unitsOfTime = -1;
186

187
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
188

189
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
190

191
  timeunits = getTimeunits(unitsOfTime);
192
193
194
195

  return (timeunits);
}

196
197
198
199
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
200
  int timeunits2 = timeunits;
201

202
203
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
204
205
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
206

207
  long lpar;
208
209
210
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
211
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
212
213
214
215
216
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

217
218
219
static
int gribapiTimeIsFC(grib_handle *gh)
{
220
  long editionNumber;
221
222
  int isFC = TRUE;

223
  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
224

225
  if ( editionNumber > 1 )
226
227
228
229
230
231
232
    {
      long sigofrtime;

      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
233
234
235
236
237
238
239

  return (isFC);
}

static
int gribapiGetTsteptype(grib_handle *gh)
{
240
  int tsteptype = TSTEP_INSTANT;
241
242
243
244
  static int lprint = TRUE;

  if ( gribapiTimeIsFC(gh) )
    {
245
      int status;
246
247
      size_t len = 256;
      char stepType[256];
248

249
250
      status = grib_get_string(gh, "stepType", stepType, &len);
      if ( status == 0 && len > 1 && len < 256 )
251
	{
252
	  if      ( strncmp("instant", stepType, len) == 0 ) tsteptype = TSTEP_INSTANT;
253
	  else if ( strncmp("avg",     stepType, len) == 0 ) tsteptype = TSTEP_AVG;
254
	  else if ( strncmp("accum",   stepType, len) == 0 ) tsteptype = TSTEP_ACCUM;
255
256
257
258
259
260
261
	  else if ( strncmp("max",     stepType, len) == 0 ) tsteptype = TSTEP_MAX;
	  else if ( strncmp("min",     stepType, len) == 0 ) tsteptype = TSTEP_MIN;
	  else if ( strncmp("diff",    stepType, len) == 0 ) tsteptype = TSTEP_DIFF;
	  else if ( strncmp("rms",     stepType, len) == 0 ) tsteptype = TSTEP_RMS;
	  else if ( strncmp("sd",      stepType, len) == 0 ) tsteptype = TSTEP_SD;
	  else if ( strncmp("cov",     stepType, len) == 0 ) tsteptype = TSTEP_COV;
	  else if ( strncmp("ratio",   stepType, len) == 0 ) tsteptype = TSTEP_RATIO;
262
	  else if ( lprint )
263
	    {
264
	      Message("Time stepType %s unsupported, set to instant!", stepType);
265
	      lprint = FALSE;
266
	    }
267
268

	  // printf("stepType: %s %ld %d\n", stepType, len, tsteptype);
269
270
271
272
273
274
	}
    }

  return (tsteptype);
}

275
276
277
278
279
280
281
282
283
284
285
static
void gribapiGetDataDateTime(grib_handle *gh, int *datadate, int *datatime)
{
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
  *datadate = (int) lpar;
  GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
  *datatime = (int) lpar*100;
}

286
287
288
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
289
290
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
291
292
}

293
static
294
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
295
{
296
  int rdate, rtime;
297
  int timeUnits, startStep = 0, endStep;
298
299
  int tstepRange = 0;
  int range;
300
  int status;
301
  long lpar;
302
303
304
305
306
  long sigofrtime = 3;
  long editionNumber;

  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

307
  if ( editionNumber > 1 )
308
309
310
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
311
312
313
314
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
315
316
317

  if ( sigofrtime == 3 )
    {
318
      gribapiGetDataDateTime(gh, vdate, vtime);
319
320
321
    }
  else
    {
322
323
      gribapiGetDataDateTime(gh, &rdate, &rtime);

324
325
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
326
      timeUnits = gribapiGetTimeUnits(gh);
327
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
328
329
330
331
332
333
334
335
336

      range = endStep - startStep;

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

337
338
339
340
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
341
342
343
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
344
345
346
347

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

348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
        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);
          }
375

376
377
378
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
379
    }
380
381

  return (tstepRange);
382
383
}

384
static
385
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
{
387
388
389
  long editionNumber;
  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

390
  int gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
392
393
394
395
396
397
398
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
399
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400

401
  size_t datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
  GRIB_CHECK(grib_get_size(gh, "values", &datasize), 0);
403
  long numberOfPoints;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404
405
  GRIB_CHECK(grib_get_long(gh, "numberOfPoints", &numberOfPoints), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
407
408
409
410
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
411
        long lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
	GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
413
414
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlon = (int)lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
416
417
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlat = (int)lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418

419
420
421
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
422
            grid->np = (int)lpar;
423
424
          }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
	if ( numberOfPoints != nlon*nlat )
426
	  Error("numberOfPoints (%ld) and gridSize (%d) differ!", numberOfPoints, nlon*nlat);
427

428
429
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
430
431
432
433
434
435
436
437
438
439
	grid->xsize = nlon;
	grid->ysize = nlat;
	grid->xinc  = 0;
	grid->yinc  = 0;
	grid->xdef  = 0;
	GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->xfirst), 0);
	GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees",  &grid->xlast), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees",  &grid->yfirst), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees",   &grid->ylast), 0);
	GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid->xinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
	if ( gridtype == GRID_LONLAT )
441
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442

443
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444

445
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446
	  {
447
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
	      {
449
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
450

451
452
453
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
454
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
		      {
456
			double xinc = 360. / grid->xsize;
457

458
			if ( fabs(grid->xinc-xinc) > 0.0 )
459
			  {
460
461
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
462
			  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463
464
465
		      }
		  }
	      }
466
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
	  }
468
	grid->ydef = 0;
469
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470
	  {
471
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
	      {
473
474
475
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476
	      }
477
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
479
480
481
482
	  }
	break;
      }
    case GRID_GAUSSIAN_REDUCED:
      {
483
484
485
	size_t dummy;
	long *pl;

486
        long lpar;
487
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
488
489
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
        grid->np = (int)lpar;
490

491
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
492
493
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlat = (int)lpar;
494

495
496
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size   = (int)numberOfPoints;
497

498
499
500
        grid->rowlon = (int *) malloc((size_t)nlat * sizeof (int));
        pl          = (long *) malloc((size_t)nlat * sizeof (long));
	dummy       = (size_t)nlat;
501
	GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
502
503
        /* FIXME: assert(pl[i] >= INT_MIN && pl[i] <= INT_MIN) */
	for (int i = 0; i < nlat; ++i ) grid->rowlon[i] = (int)pl[i];
504
505
	free(pl);

506
507
508
509
510
511
512
513
514
	grid->ysize  = nlat;
	grid->xinc   = 0;
	grid->yinc   = 0;
	grid->xdef   = 0;
	GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->xfirst), 0);
	GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees",  &grid->xlast), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees",  &grid->yfirst), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees",   &grid->ylast), 0);
	GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid->xinc), 0);
515

516
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
517

518
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
	  {
520
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
521
	      {
522
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
523
524
525
526

		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
527
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
528
		      {
529
			double xinc = 360. / grid->xsize;
530

531
			if ( fabs(grid->xinc-xinc) > 0.0 )
532
			  {
533
534
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
535
536
537
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
	      }
539
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
	  }
541
542
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
	  {
544
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
	      {
546
547
548
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
	      }
550
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
552
553
554
555
	  }
	break;
      }
    case GRID_LCC:
      {
556
	int nlon, nlat;
557
        long lpar;
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572

	GRIB_CHECK(grib_get_long(gh, "Nx", &lpar), 0);
	nlon = lpar;
	GRIB_CHECK(grib_get_long(gh, "Ny", &lpar), 0);
	nlat = lpar;

	if ( numberOfPoints != nlon*nlat )
	  Error("numberOfPoints (%d) and gridSize (%d) differ!", (int)numberOfPoints, nlon*nlat);

	grid->size  = numberOfPoints;
	grid->xsize = nlon;
	grid->ysize = nlat;

	GRIB_CHECK(grib_get_double(gh, "DxInMetres", &grid->lcc_xinc), 0);
	GRIB_CHECK(grib_get_double(gh, "DyInMetres", &grid->lcc_yinc), 0);
573
574
575
576
577
578
579
580
581
582
583
584
585
	GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->lcc_originLon), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid->lcc_originLat), 0);
	GRIB_CHECK(grib_get_double(gh, "LoVInDegrees", &grid->lcc_lonParY), 0);
	GRIB_CHECK(grib_get_double(gh, "Latin1InDegrees", &grid->lcc_lat1), 0);
	GRIB_CHECK(grib_get_double(gh, "Latin2InDegrees", &grid->lcc_lat2), 0);

        if ( editionNumber <= 1 )
          {
            GRIB_CHECK(grib_get_long(gh, "projectionCenterFlag", &lpar), 0);
            grid->lcc_projflag  = (int) lpar;
            GRIB_CHECK(grib_get_long(gh, "scanningMode", &lpar), 0);
            grid->lcc_scanflag  = (int) lpar;
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586

587
	grid->xdef   = 0;
588
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
590
591
592
593

	break;
      }
    case GRID_SPECTRAL:
      {
594
595
	size_t len = 256;
	char typeOfPacking[256];
596
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
597
598
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
599

600
601
602
        /* FIXME: assert(datasize >= INT_MIN && datasize <= INT_MAX) */
	grid->size  = (int)datasize;
        long lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
603
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
604
605
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	grid->trunc = (int)lpar;
606

Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
609
610
	break;
      }
    case GRID_GME:
      {
611
612
613
614
615
616
617
618
619
620
621
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
        long lpar;
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	if ( grib_get_long(gh, "nd", &lpar) == 0 ) grid->nd  = (int)lpar;
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	if ( grib_get_long(gh, "Ni", &lpar) == 0 ) grid->ni  = (int)lpar;
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	if ( grib_get_long(gh, "n2", &lpar) == 0 ) grid->ni2 = (int)lpar;
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	if ( grib_get_long(gh, "n3", &lpar) == 0 ) grid->ni3 = (int)lpar;
622
623
624

	break;
      }
625
    case GRID_UNSTRUCTURED:
626
      {
627
628
629
        unsigned char uuid[CDI_UUID_SIZE];
    	/*
        char reference_link[8192];
630
631
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
632
         */
633

634
635
636
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
    	grid->size  = (int)numberOfPoints;
        long lpar;
637
638
        if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
          {
639
640
641
642
643
            /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
            grid->number   = (int)lpar;
            /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
            if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 )
              grid->position = (int)lpar;
644
            /*
645
646
647
648
649
            if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
              {
                if ( strncmp(reference_link, "file://", 7) == 0 )
                  grid->reference = strdupx(reference_link);
              }
650
            */
651
652
            size_t len = (size_t)CDI_UUID_SIZE;
            if ( grib_get_bytes(gh, "uuidOfHGrid", uuid, &len) == 0)
653
              {
654
                memcpy(grid->uuid, uuid, CDI_UUID_SIZE);
655
656
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
658
659
660
	break;
      }
    case GRID_GENERIC:
      {
661
	int nlon = 0, nlat = 0;
662
663
664
665
666
        long lpar;
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	if ( grib_get_long(gh, "Ni", &lpar) == 0 ) nlon = (int)lpar;
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	if ( grib_get_long(gh, "Nj", &lpar) == 0 ) nlat = (int)lpar;
667

668
669
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
670

671
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
672
	  {
673
674
	    grid->xsize = nlon;
	    grid->ysize = nlat;
675
676
677
	  }
	else
	  {
678
679
	    grid->xsize = 0;
	    grid->ysize = 0;
680
681
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
682
683
684
685
	break;
      }
    default:
      {
686
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
687
688
689
690
	break;
      }
    }

691
  grid->isRotated = FALSE;
692
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693
    {
694
695
696
697
      grid->isRotated = TRUE;
      GRIB_CHECK(grib_get_double(gh, "latitudeOfSouthernPoleInDegrees",  &grid->ypole), 0);
      GRIB_CHECK(grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &grid->xpole), 0);
      GRIB_CHECK(grib_get_double(gh, "angleOfRotation", &grid->angle), 0);
698
      /* change from south to north pole */
699
700
701
702
703
704
705
706
707
708
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

  grid->xvals = NULL;
  grid->yvals = NULL;
  grid->type  = gridtype;
}

static
709
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
710
{
711
  *leveltype = 0;
712
713
714
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
715

716
717
  long lpar;
  int status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
718
  if ( status == 0 )
719
    {
720
      *leveltype = (int) lpar;
721

722
      switch (*leveltype)
723
724
725
726
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
727
	  { *lbounds = 1; break; }
728
	}
729

730
731
      if ( *lbounds == 0 )
	{
732
          double dlevel;
733
734
735
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
736
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
737

738
739
740
741
742
743
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
      else
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
744
	  *level1 = (int)lpar;
745
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
746
	  *level2 = (int)lpar;
747
	}
748
    }
749
750
}

751
752
753
static
double grib2ScaleFactor(long factor)
{
754
  double scaleFactor = 0;
755

756
757
758
759
760
  if      ( factor == 0 ) scaleFactor =    1;
  else if ( factor == 1 ) scaleFactor =    0.1;
  else if ( factor == 2 ) scaleFactor =    0.01;
  else if ( factor == 3 ) scaleFactor =    0.001;
  else if ( factor == 4 ) scaleFactor =    0.0001;
761
762
763
764
765
  else if ( factor == 5 ) scaleFactor =    0.00001;
  else if ( factor == 6 ) scaleFactor =    0.000001;
  else if ( factor == 7 ) scaleFactor =    0.0000001;
  else if ( factor == 8 ) scaleFactor =    0.00000001;
  else if ( factor == 9 ) scaleFactor =    0.000000001;
766
767
768
769

  return (scaleFactor);
}

770
static
771
772
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
773
{
774
775
776
  int status;
  long lpar;
  long factor;
777

778
779
  *leveltype1 = 0;
  *leveltype2 = -1;
780
781
782
783
784
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
785
786
787

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
788
    {
789
      long llevel;
790
791
      double dlevel1 = 0, dlevel2 = 0;

792
      *leveltype1 = (int) lpar;
793

794
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
795
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
796
      if ( status == 0 ) *leveltype2 = (int)lpar;
797

798
799
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
      if ( *leveltype1 == GRIB2_LTYPE_REFERENCE && *leveltype2 == 1 ) *lbounds = 0;
800

801
      if ( *leveltype1 == GRIB2_LTYPE_LANDDEPTH )
802
803
804
805
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_M;
        }
806
      else if ( *leveltype1 == GRIB2_LTYPE_ISOBARIC )
807
808
809
810
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_PA;
        }
811
      else if ( *leveltype1 == GRIB2_LTYPE_SIGMA )
812
813
814
815
        {
          *level_sf = 1000;
          *level_unit = 0;
        }
816

817
818
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
819
820
821
      if ( llevel != GRIB_MISSING_LONG )
        {
          if ( factor != GRIB_MISSING_LONG )
822
            dlevel1 = (double)llevel * grib2ScaleFactor(factor);
823
          else
824
            dlevel1 = (double)llevel;
825
        }
826

827
      if ( *level_sf != 0 ) dlevel1 *= (double)(*level_sf);
828
829

      if ( *lbounds == 1 )
830
	{
831
832
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
833
834
835
          if ( llevel != GRIB_MISSING_LONG )
            {
              if ( factor != GRIB_MISSING_LONG )
836
                dlevel2 = (double)llevel * grib2ScaleFactor(factor);
837
              else
838
                dlevel2 = (double)llevel;
839
            }
840
841
842
843
844
845

          if ( *level_sf != 0 ) dlevel2 *= (*level_sf);
        }

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
846
    }
847
848
}

849
850
851
852
853
854
855
856
857
858
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;
}

859
#if  defined  (HAVE_LIBGRIB_API)
860
static
861
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
862
863
864
865
		      size_t recsize, off_t position, int datatype, int comptype, 
                      size_t len, const char *varname, int leveltype1, int leveltype2, 
                      int lbounds, int level1, int level2, int level_sf, int level_unit,
                      var_tile_t tiles, int lread_additional_keys)
866
867
868
869
870
871
872
873
874
875
876
877
878
{
  long editionNumber;
  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;
879
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
880
  size_t vlen;
881
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
882

Uwe Schulzweida's avatar
Uwe Schulzweida committed
883
  vlistID = streamptr->vlistID;
884
  tsID    = streamptr->curTsID;
885
  recID   = recordNewEntry(streamptr, tsID);
886
887
888
889
890
891
892
893
  record  = &streamptr->tsteps[tsID].records[recID];

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

  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

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

896
897
898
899
900
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
901
  (*record).ltype     = leveltype1;
902
  (*record).tsteptype = tsteptype;
903
  memcpy((*record).varname, varname, len);
904
905

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

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

909
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
910

911
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
912
    {
913
914
915
916
917
918
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
919

920
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
921
922
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
        vctsize = (int)lpar;
923
924
925
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
926
            dummy = (size_t)vctsize;
927
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
928
            varDefVCT((size_t)vctsize, vctptr);
929
930
931
932
933
934
935
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
936
        unsigned char uuid[CDI_UUID_SIZE];
937
        long ltmp;
938
        long nhlev, nvgrid;
939
940

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
941
        if ( lpar != 6 )
942
943
944
          {
            fprintf(stderr, "Warning ...\n");
          }
945
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
946
        nhlev = ltmp;
947
948
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
949
950
951
        len = (size_t)CDI_UUID_SIZE;
        memset(uuid, 0, CDI_UUID_SIZE);
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", uuid, &len), 0);
952
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
953
954
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956

957
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
958
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
959

960
  stdname[0] = 0;
961
962
963
  longname[0] = 0;
  units[0] = 0;

964
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
965
    {