stream_gribapi.c 92.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
26
#  include "grib_api.h"
#endif

extern int cdiInventoryMode;

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


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

71
    {
72
73
74
75
      int status;
      status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);

      if ( status ==  0 ) gribgridtype = (int) lpar;
76
77
78

      switch (gribgridtype)
	{
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
	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; }
94
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
97
98
99
    }

  return (gridtype);
}

100
static
101
int gribapiGetIsRotated(grib_handle *gh)
102
103
{
  int isRotated = 0;
104
  int gribgridtype = -1;
105
  long lpar;
106
  int status;
107

108
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
109

110
111
112
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
113
114
115
116

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
static
118
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
122
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
    {
124
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
127
    }
  else
    {
128
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
131
132
133
    }

  return (zaxistype);
}

134
static
135
int getTimeunits(long unitsOfTime)
136
137
{
  int timeunits = -1;
138
139
140
141
142
143
144
145
146
147
148
149
150

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

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
  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;
180
  long unitsOfTime = -1;
181

182
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
183

184
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
185

186
  timeunits = getTimeunits(unitsOfTime);
187
188
189
190

  return (timeunits);
}

191
192
193
194
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
195
  int timeunits2 = timeunits;
196

197
198
  long unitsOfTime;
  int status = grib_get_long(gh, "stepUnits", &unitsOfTime);
199
200
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
201

202
  long lpar;
203
204
205
  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
206
    endStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
207
208
209
210
211
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

212
213
214
static
int gribapiTimeIsFC(grib_handle *gh)
{
215
  long editionNumber;
216
217
  int isFC = TRUE;

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

220
  if ( editionNumber > 1 )
221
222
223
224
225
226
227
    {
      long sigofrtime;

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

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
228
229
230
231
232
233
234

  return (isFC);
}

static
int gribapiGetTsteptype(grib_handle *gh)
{
235
  int tsteptype = TSTEP_INSTANT;
236
237
238
239
  static int lprint = TRUE;

  if ( gribapiTimeIsFC(gh) )
    {
240
      int status;
241
242
      size_t len = 256;
      char stepType[256];
243

244
245
      status = grib_get_string(gh, "stepType", stepType, &len);
      if ( status == 0 && len > 1 && len < 256 )
246
	{
247
	  if      ( strncmp("instant", stepType, len) == 0 ) tsteptype = TSTEP_INSTANT;
248
	  else if ( strncmp("avg",     stepType, len) == 0 ) tsteptype = TSTEP_AVG;
249
	  else if ( strncmp("accum",   stepType, len) == 0 ) tsteptype = TSTEP_ACCUM;
250
251
252
253
254
255
256
	  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;
257
	  else if ( lprint )
258
	    {
259
260
	      Message("stepType %s unsupported, set to instant!", stepType);
	      lprint = FALSE;
261
	    }
262
263

	  // printf("stepType: %s %ld %d\n", stepType, len, tsteptype);
264
265
266
267
268
269
	}
    }

  return (tsteptype);
}

270
271
272
273
274
275
276
277
278
279
280
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;
}

281
282
283
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
284
285
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
286
287
}

288
static
289
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
290
{
291
  int rdate, rtime;
292
  int timeUnits, startStep = 0, endStep;
293
294
  int tstepRange = 0;
  int range;
295
  int status;
296
  long lpar;
297
298
299
300
301
  long sigofrtime = 3;
  long editionNumber;

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

302
  if ( editionNumber > 1 )
303
304
305
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
306
307
308
309
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
310
311
312

  if ( sigofrtime == 3 )
    {
313
      gribapiGetDataDateTime(gh, vdate, vtime);
314
315
316
    }
  else
    {
317
318
      gribapiGetDataDateTime(gh, &rdate, &rtime);

319
320
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
321
      timeUnits = gribapiGetTimeUnits(gh);
322
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
323
324
325
326
327
328
329
330
331

      range = endStep - startStep;

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

332
333
334
335
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
336
337
338
	int julday, secofday;
	int64_t time_period = endStep;
        int64_t addsec;
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360

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

	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;
	      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
	    break;
362
363
364
365
366
	  }

	julday_add_seconds(addsec, &julday, &secofday);

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

368
369
370
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
371
    }
372
373

  return (tstepRange);
374
375
}

376
static
377
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378
{
379
380
381
  long editionNumber;
  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);

382
  int gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
384
385
386
387
388
389
390
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
391
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
392

393
  size_t datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
  GRIB_CHECK(grib_get_size(gh, "values", &datasize), 0);
395
  long numberOfPoints;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
397
  GRIB_CHECK(grib_get_long(gh, "numberOfPoints", &numberOfPoints), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
399
400
401
402
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
403
        long lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404
	GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
405
406
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlon = (int)lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
408
409
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlat = (int)lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410

411
412
413
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
414
            grid->np = (int)lpar;
415
416
          }

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

420
421
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
422
423
424
425
426
427
428
429
430
431
	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
432
	if ( gridtype == GRID_LONLAT )
433
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434

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

437
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
	  {
439
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
	      {
441
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442

443
444
445
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
446
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
		      {
448
			double xinc = 360. / grid->xsize;
449

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

478
        long lpar;
479
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
480
481
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
        grid->np = (int)lpar;
482

483
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
484
485
        /* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
	int nlat = (int)lpar;
486

487
488
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size   = (int)numberOfPoints;
489

490
491
492
        grid->rowlon = (int *) malloc((size_t)nlat * sizeof (int));
        pl          = (long *) malloc((size_t)nlat * sizeof (long));
	dummy       = (size_t)nlat;
493
	GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
494
495
        /* FIXME: assert(pl[i] >= INT_MIN && pl[i] <= INT_MIN) */
	for (int i = 0; i < nlat; ++i ) grid->rowlon[i] = (int)pl[i];
496
497
	free(pl);

498
499
500
501
502
503
504
505
506
	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);
507

508
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
509

510
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
	  {
512
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513
	      {
514
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
515
516
517
518

		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
519
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
520
		      {
521
			double xinc = 360. / grid->xsize;
522

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

	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);
565
566
567
568
569
570
571
572
573
574
575
576
577
	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
578

579
	grid->xdef   = 0;
580
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
581
582
583
584
585

	break;
      }
    case GRID_SPECTRAL:
      {
586
587
	size_t len = 256;
	char typeOfPacking[256];
588
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
589
590
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
591

592
593
594
        /* FIXME: assert(datasize >= INT_MIN && datasize <= INT_MAX) */
	grid->size  = (int)datasize;
        long lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
595
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
596
597
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
	grid->trunc = (int)lpar;
598

Uwe Schulzweida's avatar
Uwe Schulzweida committed
599
600
601
602
	break;
      }
    case GRID_GME:
      {
603
604
605
606
607
608
609
610
611
612
613
        /* 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;
614
615
616

	break;
      }
617
    case GRID_UNSTRUCTURED:
618
      {
619
        char uuid[17];
620
    	char reference_link[8192];
621
622
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
623

624
625
626
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
    	grid->size  = (int)numberOfPoints;
        long lpar;
627
628
        if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
          {
629
630
631
632
633
            /* 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;
634
            /*
635
636
637
638
639
            if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
              {
                if ( strncmp(reference_link, "file://", 7) == 0 )
                  grid->reference = strdupx(reference_link);
              }
640
            */
641
642
643
            len = (size_t) 16;
            if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
              {
644
                memcpy(grid->uuid, uuid, 16);
645
646
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
647
648
649
650
	break;
      }
    case GRID_GENERIC:
      {
651
	int nlon = 0, nlat = 0;
652
653
654
655
656
        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;
657

658
659
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
660

661
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
662
	  {
663
664
	    grid->xsize = nlon;
	    grid->ysize = nlat;
665
666
667
	  }
	else
	  {
668
669
	    grid->xsize = 0;
	    grid->ysize = 0;
670
671
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
674
675
	break;
      }
    default:
      {
676
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677
678
679
680
	break;
      }
    }

681
  grid->isRotated = FALSE;
682
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
683
    {
684
685
686
687
      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);
688
      /* change from south to north pole */
689
690
691
692
693
694
695
696
697
698
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

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

static
699
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
700
{
701
702
703
704
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
705

706
707
  long lpar;
  int status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
708
  if ( status == 0 )
709
    {
710
      *leveltype = (int) lpar;
711

712
      switch (*leveltype)
713
714
715
716
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
717
	  { *lbounds = 1; break; }
718
	}
719

720
721
      if ( *lbounds == 0 )
	{
722
          double dlevel;
723
724
725
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
726
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
727

728
729
730
731
732
733
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
      else
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
734
	  *level1 = (int)lpar;
735
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
736
	  *level2 = (int)lpar;
737
	}
738
    }
739
740
}

741
742
743
static
double grib2ScaleFactor(long factor)
{
744
  double scaleFactor = 0;
745

746
747
748
749
750
  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;
751
752
753
754
755
  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;
756
757
758
759

  return (scaleFactor);
}

760
static
761
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf, int *level_unit)
762
{
763
764
765
766
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
767

768
769
770
771
772
773
  *leveltype  = 0;
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
774
775
776

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
777
    {
778
      long llevel;
779
780
      double dlevel1 = 0, dlevel2 = 0;

781
      *leveltype = (int) lpar;
782

783
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
784
785
      /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
      if ( status == 0 ) leveltype2 = (int)lpar;
786

787
      if ( *leveltype != 255 && leveltype2 != 255 && leveltype2 > 0 ) *lbounds = 1;
788
      if ( *leveltype == GRIB2_LTYPE_REFERENCE && leveltype2 == 1 ) *lbounds = 0;
789

790
791
792
793
794
795
796
797
798
799
      if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_M;
        }
      else if ( *leveltype == GRIB2_LTYPE_ISOBARIC )
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_PA;
        }
800
801
802
803
804
      else if ( *leveltype == GRIB2_LTYPE_SIGMA )
        {
          *level_sf = 1000;
          *level_unit = 0;
        }
805

806
807
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
808
809
810
      if ( llevel != GRIB_MISSING_LONG )
        {
          if ( factor != GRIB_MISSING_LONG )
811
            dlevel1 = (double)llevel * grib2ScaleFactor(factor);
812
          else
813
            dlevel1 = (double)llevel;
814
        }
815

816
      if ( *level_sf != 0 ) dlevel1 *= (double)(*level_sf);
817
818

      if ( *lbounds == 1 )
819
	{
820
821
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
822
823
824
          if ( llevel != GRIB_MISSING_LONG )
            {
              if ( factor != GRIB_MISSING_LONG )
825
                dlevel2 = (double)llevel * grib2ScaleFactor(factor);
826
              else
827
                dlevel2 = (double)llevel;
828
            }
829
830
831
832
833
834

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
835
    }
836
837
}

838
839
840
841
842
843
844
845
846
847
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;
}

848
#if  defined  (HAVE_LIBGRIB_API)
849
static
850
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
851
		      size_t recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
852
                      int leveltype, int lbounds, int level1, int level2, int level_sf, int level_unit)
853
854
855
856
857
858
859
860
861
862
863
864
865
{
  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;
866
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
867
  size_t vlen;
868
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
869

Uwe Schulzweida's avatar
Uwe Schulzweida committed
870
  vlistID = streamptr->vlistID;
871
  tsID    = streamptr->curTsID;
872
  recID   = recordNewEntry(streamptr, tsID);
873
874
875
876
877
878
879
880
881
882
  record  = &streamptr->tsteps[tsID].records[recID];

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

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

  // fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, leveltype);

883
884
885
886
887
888
889
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
  (*record).ltype     = leveltype;
  (*record).tsteptype = tsteptype;
890
  memcpy((*record).varname, varname, len);
891
892

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
893
894
895

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
896
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
897

898
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
    {
900
901
902
903
904
905
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
906

907
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
908
909
        /* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
        vctsize = (int)lpar;
910
911
912
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
913
            dummy = (size_t)vctsize;
914
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
915
            varDefVCT((size_t)vctsize, vctptr);
916
917
918
919
920
921
922
923
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
        char uuid[17];
924
        long ltmp;
925
        long nhlev, nvgrid;
926
927

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
928
        if ( lpar != 6 )
929
930
931
          {
            fprintf(stderr, "Warning ...\n");
          }
932
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
933
        nhlev = ltmp;
934
935
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
936
937
        len = (size_t) 16;
        uuid[16] = 0;
938
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
939
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
940
941
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
943

944
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
945
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
946

947
  stdname[0] = 0;
948
949
950
  longname[0] = 0;
  units[0] = 0;

951
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
952
    {
953
      vlen = CDI_MAX_NAME;
954
      gribapiGetString(gh, "name", longname, vlen);
955
      vlen = CDI_MAX_NAME;
956
      gribapiGetString(gh, "units", units, vlen);
957
958
959
960
961

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

967
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
968
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
969
	       varname, stdname, longname, units);
Uwe Schulzweida's avatar