stream_gribapi.c 93.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
  *leveltype = 0;
702
703
704
  *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
762
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
763
{
764
765
766
  int status;
  long lpar;
  long factor;
767

768
769
  *leveltype1 = 0;
  *leveltype2 = -1;
770
771
772
773
774
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
775
776
777

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

782
      *leveltype1 = (int) lpar;
783

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  gridID = varDefGrid(vlistID, grid, 0);

897
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898

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

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

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

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

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

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

      {
        vlen = CDI_MAX_NAME;
        status = grib_get_string(gh, "cfName", stdname, &vlen);
        if ( status != 0 || vlen <= 1 ) stdname[0] = 0;