stream_gribapi.c 93.9 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
623
624
        unsigned char uuid[CDI_UUID_SIZE];
    	/*
        char reference_link[8192];
625
626
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
627
         */
628

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

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

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

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

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

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

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

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

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

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

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

746
747
748
static
double grib2ScaleFactor(long factor)
{
749
  double scaleFactor = 0;
750

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

  return (scaleFactor);
}

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

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

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

787
      *leveltype1 = (int) lpar;
788

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

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

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

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

822
      if ( *level_sf != 0 ) dlevel1 *= (double)(*level_sf);
823
824

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

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
841
    }
842
843
}

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

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

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

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

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

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

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

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
900
901

  gridID = varDefGrid(vlistID, grid, 0);

902
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
903

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

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

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

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

953
  stdname[0] = 0;
954
955
956
  longname[0] = 0;
  units[0] = 0;

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