stream_gribapi.c 94 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
	      Message("Time stepType %s unsupported, set to instant!", stepType);
260
	      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

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

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
        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);
          }
370

371
372
373
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
374
    }
375
376

  return (tstepRange);
377
378
}

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

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

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

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

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

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

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

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

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

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

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

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

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

490
491
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size   = (int)numberOfPoints;
492

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

501
502
503
504
505
506
507
508
509
	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);
510

511
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
512

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

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

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

	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);
568
569
570
571
572
573
574
575
576
577
578
579
580
	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
581

582
	grid->xdef   = 0;
583
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
584
585
586
587
588

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

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

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

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

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

661
662
        /* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
	grid->size  = (int)numberOfPoints;
663

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

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

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

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

static
702
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
703
{
704
  *leveltype = 0;
705
706
707
  *lbounds   = 0;
  *level1    = 0;
  *level2    = 0;
708

709
710
  long lpar;
  int status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
711
  if ( status == 0 )
712
    {
713
      *leveltype = (int) lpar;
714

715
      switch (*leveltype)
716
717
718
719
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
720
	  { *lbounds = 1; break; }
721
	}
722

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

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

744
745
746
static
double grib2ScaleFactor(long factor)
{
747
  double scaleFactor = 0;
748

749
750
751
752
753
  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;
754
755
756
757
758
  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;
759
760
761
762

  return (scaleFactor);
}

763
static
764
765
void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, 
                   int *level2, int *level_sf, int *level_unit)
766
{
767
768
769
  int status;
  long lpar;
  long factor;
770

771
772
  *leveltype1 = 0;
  *leveltype2 = -1;
773
774
775
776
777
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
778
779
780

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
781
    {
782
      long llevel;
783
784
      double dlevel1 = 0, dlevel2 = 0;

785
      *leveltype1 = (int) lpar;
786

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

791
792
      if ( *leveltype1 != 255 && *leveltype2 != 255 && *leveltype2 > 0 ) *lbounds = 1;
      if ( *leveltype1 == GRIB2_LTYPE_REFERENCE && *leveltype2 == 1 ) *lbounds = 0;
793

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

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

820
      if ( *level_sf != 0 ) dlevel1 *= (double)(*level_sf);
821
822

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

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
839
    }
840
841
}

842
843
844
845
846
847
848
849
850
851
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;
}

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

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

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

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

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

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

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
897
898
899

  gridID = varDefGrid(vlistID, grid, 0);

900
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
901

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

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

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

948
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
949
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
950

951
  stdname[0] = 0;
952
953
954
  longname[0] = 0;
  units[0] = 0;

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