stream_gribapi.c 85.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
9
10
11
12
13
#include "cdi.h"
#include "stream_int.h"
#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
19


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

24
#define  NINT(x)  ((x) < 0 ? (int)((x)-.5) : (int)((x)+.5))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
28

extern int cdiInventoryMode;

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


#if  defined  (HAVE_LIBGRIB_API)
static
39
int gribapiGetGridType(grib_handle *gh)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
{
41
  int gridtype = GRID_GENERIC;
42
  int gribgridtype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
  long lpar;

45
    {
46
47
48
49
      int status;
      status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);

      if ( status ==  0 ) gribgridtype = (int) lpar;
50
51
52

      switch (gribgridtype)
	{
53
	case  GRIB2_GTYPE_LATLON:     { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
	                                if ( lpar == (long) GRIB_MISSING_LONG ) break;
55
                                      }
56
57
	case  GRIB2_GTYPE_LATLON_ROT: { gridtype = GRID_LONLAT;    break; }
	case  GRIB2_GTYPE_LCC:        { gridtype = GRID_LCC;       break; }
58
	case  GRIB2_GTYPE_GAUSSIAN:   { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
	                                if ( lpar == (long) GRIB_MISSING_LONG )
60
61
62
63
64
					  gridtype = GRID_GAUSSIAN_REDUCED;
					else
					  gridtype = GRID_GAUSSIAN;
					break;
                                      }
65
66
67
	case  GRIB2_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL;  break; }
	case  GRIB2_GTYPE_GME:        { gridtype = GRID_GME;       break; }
	case  GRIB2_GTYPE_NUMBER:     { gridtype = GRID_REFERENCE; break; }
68
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
71
72
73
    }

  return (gridtype);
}

74
static
75
int gribapiGetIsRotated(grib_handle *gh)
76
77
{
  int isRotated = 0;
78
  int gribgridtype = -1;
79
  long lpar;
80
  int status;
81

82
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
83

84
85
86
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
87
88
89
90

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
static
92
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
95
96
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
    {
98
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
    }
  else
    {
102
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
104
105
106
107
    }

  return (zaxistype);
}

108
static
109
int getTimeunits(long unitsOfTime)
110
111
{
  int timeunits = -1;
112
113
114
115
116
117
118
119
120
121
122
123
124

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

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  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;
154
155
  long unitsOfTime = -1;
  int status;
156
157
158
159
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

160
  status = grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
161
162
163

  timeunits = getTimeunits(unitsOfTime);

164
  /*
165
166
  GRIB_CHECK(grib_get_string(gh, "stepUnits", stepunits, &len), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
168
169
170
171
172
173
174
175
176
177
  len--;

  if      ( memcmp(stepunits, "s",   len) == 0 ) timeunits = TUNIT_SECOND;
  else if ( memcmp(stepunits, "m",   len) == 0 ) timeunits = TUNIT_MINUTE;
  else if ( memcmp(stepunits, "h",   len) == 0 ) timeunits = TUNIT_HOUR;
  else if ( memcmp(stepunits, "3h",  len) == 0 ) timeunits = TUNIT_3HOURS;
  else if ( memcmp(stepunits, "6h",  len) == 0 ) timeunits = TUNIT_6HOURS;
  else if ( memcmp(stepunits, "12h", len) == 0 ) timeunits = TUNIT_12HOURS;
  else if ( memcmp(stepunits, "D",   len) == 0 ) timeunits = TUNIT_DAY;
  else if ( memcmp(stepunits, "M",   len) == 0 ) timeunits = TUNIT_MONTH;
  else if ( memcmp(stepunits, "Y",   len) == 0 ) timeunits = TUNIT_YEAR;
178
  else if ( lprint )
179
    {
180
      Message("Step units >%s< unsupported!", stepunits);
181
      lprint = FALSE;
182
    }
183
  */
184
185
186
187

  return (timeunits);
}

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
  int timeunits2;
  int status;
  long unitsOfTime;
  long lpar;

  // status = grib_get_long(gh, "stepUnits", &unitsOfTime);

  // timeunits2 = getTimeunits(unitsOfTime);
  timeunits2 = gribapiGetTimeUnits(gh);

  status = grib_get_long(gh, "endStep", &lpar);

  if ( status == 0 )
    endStep = (int) ((lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
  // printf("%d %d %d %d %d %g\n", startStep, endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2));

  return (endStep);
}

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

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

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

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

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

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

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

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

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

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

301
  if ( editionNumber > 1 )
302
303
304
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
305
306
307

  if ( sigofrtime == 3 )
    {
308
      gribapiGetDataDateTime(gh, vdate, vtime);
309
310
311
    }
  else
    {
312
313
      gribapiGetDataDateTime(gh, &rdate, &rtime);

314
315
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
316
      timeUnits = gribapiGetTimeUnits(gh);
317
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
318
319
320
321
322
323
324
325
326

      range = endStep - startStep;

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

327
328
329
330
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
331
	int time_period = endStep;
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
	int julday, secofday, addsec;

	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
355
	    break;
356
357
358
359
360
	  }

	julday_add_seconds(addsec, &julday, &secofday);

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

362
363
364
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
365
    }
366
367

  return (tstepRange);
368
369
}

370
static
371
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
{
373
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
375
  size_t datasize;
376
377
  long numberOfPoints;
  long lpar;
378
379
380

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
393
394
  GRIB_CHECK(grib_get_size(gh, "values", &datasize), 0);
  GRIB_CHECK(grib_get_long(gh, "numberOfPoints", &numberOfPoints), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
397
398
399
400
401
402
403
404
405
406
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
	int nlon, nlat;

	GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
	nlon = lpar;
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

407
408
409
410
411
412
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
	if ( numberOfPoints != nlon*nlat )
414
	  Error("numberOfPoints (%d) and gridSize (%d) differ!", (int)numberOfPoints, nlon*nlat);
415

416
417
418
419
420
421
422
423
424
425
426
	grid->size  = numberOfPoints;
	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
427
	if ( gridtype == GRID_LONLAT )
428
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429

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

432
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
	  {
434
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
	      {
436
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437

438
439
440
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
441
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
		      {
443
			double xinc = 360. / grid->xsize;
444

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

474
475
476
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

477
478
479
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

480
	grid->size   = numberOfPoints;
481

482
        grid->rowlon = (int *) malloc(nlat*sizeof(int));
483
484
485
        pl          = (long *) malloc(nlat*sizeof(long));
	dummy       = nlat;
	GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
486
	for ( i = 0; i < nlat; ++i ) grid->rowlon[i] = pl[i];
487
488
	free(pl);

489
490
491
492
493
494
495
496
497
	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);
498

499
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
500

501
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
	  {
503
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
	      {
505
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
506
507
508
509

		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
510
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
511
		      {
512
			double xinc = 360. / grid->xsize;
513

514
			if ( fabs(grid->xinc-xinc) > 0.0 )
515
			  {
516
517
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
518
519
520
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
521
	      }
522
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
	  }
524
525
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526
	  {
527
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
528
	      {
529
530
531
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
	      }
533
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
535
536
	  }
	break;
      }
537
      /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
539
540
    case GRID_LCC:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
541
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
543
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);

544
545
546
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547

548
549
550
551
552
553
554
555
556
	grid->lcc_xinc      = ISEC2_Lambert_dx;
	grid->lcc_yinc      = ISEC2_Lambert_dy;
	grid->lcc_originLon = ISEC2_FirstLon * 0.001;
	grid->lcc_originLat = ISEC2_FirstLat * 0.001;
	grid->lcc_lonParY   = ISEC2_Lambert_Lov * 0.001;
	grid->lcc_lat1      = ISEC2_Lambert_LatS1 * 0.001;
	grid->lcc_lat2      = ISEC2_Lambert_LatS2 * 0.001;
	grid->lcc_projflag  = ISEC2_Lambert_ProjFlag;
	grid->lcc_scanflag  = ISEC2_ScanFlag;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557

558
	grid->xdef   = 0;
559
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
561
562

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
565
    case GRID_SPECTRAL:
      {
566
567
	size_t len = 256;
	char typeOfPacking[256];
568
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
569
570
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
571

572
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
573
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
574
	grid->trunc = lpar;
575

Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
577
578
579
	break;
      }
    case GRID_GME:
      {
580
581
582
583
584
	grid->size  = numberOfPoints;
	if ( grib_get_long(gh, "nd", &lpar) == 0 ) grid->nd  = lpar;
	if ( grib_get_long(gh, "Ni", &lpar) == 0 ) grid->ni  = lpar;
	if ( grib_get_long(gh, "n2", &lpar) == 0 ) grid->ni2 = lpar;
	if ( grib_get_long(gh, "n3", &lpar) == 0 ) grid->ni3 = lpar;
585
586
587

	break;
      }
588
    case GRID_REFERENCE:
589
      {
590
        char uuid[17];
591
    	char reference_link[8192];
592
593
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
594

595
    	grid->size  = numberOfPoints;
596
597
598
599
600
601
602
603
604
605
606
607
        if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
          {
            grid->number   = lpar;
            if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
            if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
              {
                if ( strncmp(reference_link, "file://", 7) == 0 )
                  grid->reference = strdupx(reference_link);
              }
            len = (size_t) 16;
            if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
              {
608
                memcpy(grid->uuid, uuid, 16);
609
610
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
612
613
614
	break;
      }
    case GRID_GENERIC:
      {
615
616
617
618
619
	int nlon = 0, nlat = 0;

	if ( grib_get_long(gh, "Ni", &lpar) == 0 ) nlon = lpar;
	if ( grib_get_long(gh, "Nj", &lpar) == 0 ) nlat = lpar;

620
	grid->size  = numberOfPoints;
621

622
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
623
	  {
624
625
	    grid->xsize = nlon;
	    grid->ysize = nlat;
626
627
628
	  }
	else
	  {
629
630
	    grid->xsize = 0;
	    grid->ysize = 0;
631
632
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
634
635
636
	break;
      }
    default:
      {
637
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
638
639
640
641
	break;
      }
    }

642
  grid->isRotated = FALSE;
643
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
    {
645
646
647
648
      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);
649
      /* change from south to north pole */
650
651
652
653
654
655
656
657
658
659
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

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

static
660
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
661
{
662
663
664
  int status;
  long lpar;
  double dlevel;
665

666
667
668
669
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
670

671
672
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
673
    {
674
      *leveltype = (int) lpar;
675

676
      switch (*leveltype)
677
678
679
680
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
681
	  { *lbounds = 1; break; }
682
	}
683

684
685
686
687
688
689
      if ( *lbounds == 0 )
	{
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
	  if ( *leveltype == 99 ) *leveltype = 100;
690

691
692
693
694
695
696
697
698
699
700
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
      else
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	  *level1 = lpar;
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	  *level2 = lpar;
	}
701
    }
702
703
704
}

static
705
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf)
706
{
707
708
709
710
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
711
712
  double dlevel;

713
714
715
716
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
717
  *level_sf = 0;
718
719
720

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
721
    {
722
      *leveltype = (int) lpar;
723

724
725
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
726

727
      if ( *leveltype != 255 && *leveltype2 != 255 &&  *leveltype2 > 0 ) *lbounds = 1;
728

729
      if ( *lbounds == 0 )
730
	{
731
732
733
734
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
735
736
737
738
739
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
              *level_sf = 77;
740
741
742
743
744
745
746
747
748
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	      if ( *leveltype == GRIB2_LTYPE_ISOBARIC ) dlevel *= 100;
	      if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
	      if ( *leveltype == 99 ) *leveltype = 100;
	    }

749
	  *level1 = (int) dlevel;
750
	  *level2 = 0;
751
752
753
	}
      else
	{
754
755
756
757
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
758
759
760
761
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
762
763
764
	      *level1 = (int) dlevel;
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfSecondFixedSurface", &dlevel), 0);
765
766
767
768
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
769
	      *level2 = (int) dlevel;
770
              *level_sf = 77;
771
772
773
774
775
776
777
778
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	      *level1 = lpar;
	      GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	      *level2 = lpar;
	    }
779
780
	}
    }
781
782
}

783
784
785
786
787
788
789
790
791
792
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;
}

793
static
794
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
795
796
		      long recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
                      int leveltype, int lbounds, int level1, int level2, int level_sf)
797
798
799
800
801
802
803
804
805
806
807
808
809
{
  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;
810
  char longname[256], units[256];
811
  size_t vlen;
812
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
813

Uwe Schulzweida's avatar
Uwe Schulzweida committed
814
  vlistID = streamptr->vlistID;
815
  tsID    = streamptr->curTsID;
816
  recID   = recordNewEntry(streamptr, tsID);
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
  record  = &streamptr->tsteps[tsID].records[recID];

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

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

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

  (*record).size     = recsize;
  (*record).position = position;
  (*record).param    = param;
  (*record).ilevel   = level1;
  (*record).ilevel2  = level2;
  (*record).ltype    = leveltype;
833
  memcpy((*record).varname, varname, len);
834
835

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836
837
838

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
840

841
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
842
    {
843
844
845
846
847
848
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
849

850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
        vctsize = lpar;
        if ( vctsize > 0 )
          {
            vctptr = (double *) malloc(vctsize*sizeof(double));
            dummy = vctsize;
            GRIB_CHECK(grib_get_double_array(gh, "pv", vctptr, &dummy), 0);
            varDefVCT(vctsize, vctptr);
            free(vctptr);
          }
        break;
      }
    case ZAXIS_REFERENCE:
      {
        size_t len;
        char uuid[17];
866
        double dtmp;
867
868
869
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
870
        if ( lpar != 6 )
871
872
873
          {
            fprintf(stderr, "Warning ...\n");
          }
874
875
876
877
        GRIB_CHECK(grib_get_double(gh, "nlev", &dtmp), 0);
        nlev = (int) NINT(dtmp);
        GRIB_CHECK(grib_get_double(gh, "numberOfVGridUsed", &dtmp), 0);
        nvgrid = NINT(dtmp);
878
879
        len = (size_t) 16;
        uuid[16] = 0;
880
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
881
882
883
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
884
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885

886
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
887
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888

889
890
891
  longname[0] = 0;
  units[0] = 0;

892
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
893
894
    {
      vlen = 256;
895
      gribapiGetString(gh, "name", longname, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
896
      vlen = 256;
897
      gribapiGetString(gh, "units", units, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898
    }
899
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
900

901
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf,
902
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
903
	       varname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
904
905
906
907

  (*record).varID   = varID;
  (*record).levelID = levelID;

908
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
909

910
911
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
912
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
913
  */
914
915
916
917
918
919
920
921
922
  {
    int status;
    status = grib_get_long(gh, "typeOfEnsembleForecast", &ens_forecast_type );
    if ( status == 0 )
      {
        GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", &ens_count ), 0);
        GRIB_CHECK(grib_get_long(gh, "perturbationNumber", &ens_index ), 0);
      }
  }
923
924
925
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
  int    i;
  long   lval;
  double dval;

  /* we read the additional keys for the first variable record only. */
  int linitial_field = (varOptGribNentries(varID) == 0);

  for ( i = 0; i < cdiNAdditionalGRIBKeys; i++ )
    {
      if ( linitial_field )
	{
	  if ( grib_get_long(gh, cdiAdditionalGRIBKeys[i], &lval) == 0 )
            varDefOptGribInt(varID, lval, cdiAdditionalGRIBKeys[i]);
	}
      if ( linitial_field )
	{
	  if ( grib_get_double(gh, cdiAdditionalGRIBKeys[i], &dval) == 0 )
            varDefOptGribInt(varID, dval, cdiAdditionalGRIBKeys[i]);
	}
      /* note: if the key is not defined, we do not throw an error! */
    }
947

Uwe Schulzweida's avatar
Uwe Schulzweida committed
948
949
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
950
951
952
953
954
      long center, subcenter;
      int instID;
      GRIB_CHECK(grib_get_long(gh, "centre", &center), 0);
      GRIB_CHECK(grib_get_long(gh, "subCentre", &subcenter), 0);
      instID    = institutInq((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
959
960
961
962
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
963
      long processID;
964
965
966
967
968
969
970
971
      status = grib_get_long(gh, "generatingProcessIdentifier", &processID);
      if ( status == 0 )
	{
	  modelID = modelInq(varInqInst(varID), processID, NULL);
	  if ( modelID == CDI_UNDEFID )
	    modelID = modelDef(varInqInst(varID), processID, NULL);
	  varDefModel(varID, modelID);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
973
974
975
    }

  if ( varInqTable(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
976
      int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
977

Uwe Schulzweida's avatar
Uwe Schulzweida committed
978
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
979

Uwe Schulzweida's avatar
Uwe Schulzweida committed
980
981
982
983
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
984

Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
	  tableID = tableInq(varInqModel(varID), tabnum, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
986

Uwe Schulzweida's avatar
Uwe Schulzweida committed
987
988
989
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
990
991
992
993
994
995
996
	}
    }

  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;

  if ( CDI_Debug )
997
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
998
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
999
}
1000
1001
1002
1003
1004
1005

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
1006
  int status;
1007
1008
1009
1010
1011
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "discipline", &lpar), 0);
  pdis = (int) lpar;

1012
1013
  status = grib_get_long(gh, "parameterCategory", &lpar);
  if ( status == 0 ) pcat = (int) lpar;
1014

1015
1016
  status = grib_get_long(gh, "parameterNumber", &lpar);
  if ( status == 0 ) pnum = (int) lpar;
1017
1018

  param = cdiEncodeParam(pnum, pcat, pdis);
1019

1020
1021
  return (param);
}
1022
1023

static
1024
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
1025
1026
{
  compvar2_t compVar;
1027
1028
1029
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
1030
1031
1032
1033
1034

  compVar.param  = param;
  compVar.level1 = level1;
  compVar.level2 = level2;
  compVar.ltype  = leveltype;
1035
1036
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);