stream_gribapi.c 83.7 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
24
25
26
27
#  include "grib_api.h"
#endif


extern int cdiInventoryMode;

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


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

43
44
45
46
47
48
    {
      GRIB_CHECK(grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar), 0);
      gribgridtype = (int) lpar;

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

  return (gridtype);
}

70
static
71
int gribapiGetIsRotated(grib_handle *gh)
72
73
{
  int isRotated = 0;
74
75
  int gribgridtype;
  long lpar;
76
77

    {
78
79
80
81
      GRIB_CHECK(grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar), 0);
      gribgridtype = (int) lpar;

      if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
82
83
84
85
86
    }

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
static
88
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
91
92
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
    {
94
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
97
    }
  else
    {
98
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
102
103
    }

  return (zaxistype);
}

104
static
105
int getTimeunits(long unitsOfTime)
106
107
{
  int timeunits = -1;
108
109
110
111
112
113
114
115
116
117
118
119
120

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

121
122
123
124
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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
  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 gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
  int timeunits2;
  long unitsOfTime;
  long lpar;

  GRIB_CHECK(grib_get_long(gh, "stepUnits", &unitsOfTime), 0);

  timeunits2 = getTimeunits(unitsOfTime);

  GRIB_CHECK(grib_get_long(gh, "endStep", &lpar), 0);

  endStep = (int)  ((lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);

  return (endStep);
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  int timeunits = -1;
  long unitsOfTime;
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

  GRIB_CHECK(grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime), 0);

  timeunits = getTimeunits(unitsOfTime);

178
  /*
179
180
  GRIB_CHECK(grib_get_string(gh, "stepUnits", stepunits, &len), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
184
185
186
187
188
189
190
191
  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;
192
  else if ( lprint )
193
    {
194
      Message("Step units >%s< unsupported!", stepunits);
195
      lprint = FALSE;
196
    }
197
  */
198
199
200
201
202
203
204

  return (timeunits);
}

static
int gribapiTimeIsFC(grib_handle *gh)
{
205
  long editionNumber;
206
207
  int isFC = TRUE;

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

210
  if ( editionNumber > 1 )
211
212
213
214
215
216
217
    {
      long sigofrtime;

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

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
218
219
220
221
222
223
224

  return (isFC);
}

static
int gribapiGetTsteptype(grib_handle *gh)
{
225
  int tsteptype = TSTEP_INSTANT;
226
227
228
229
  static int lprint = TRUE;

  if ( gribapiTimeIsFC(gh) )
    {
230
      int status;
231
232
      size_t len = 256;
      char stepType[256];
233

234
235
      status = grib_get_string(gh, "stepType", stepType, &len);
      if ( status == 0 && len > 1 && len < 256 )
236
	{
237
	  if      ( strncmp("instant", stepType, len) == 0 ) tsteptype = TSTEP_INSTANT;
238
	  else if ( strncmp("avg",     stepType, len) == 0 ) tsteptype = TSTEP_AVG;
239
	  else if ( strncmp("accum",   stepType, len) == 0 ) tsteptype = TSTEP_ACCUM;
240
241
242
243
244
245
246
	  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;
247
	  else if ( lprint )
248
	    {
249
250
	      Message("stepType %s unsupported, set to instant!", stepType);
	      lprint = FALSE;
251
	    }
252
253

	  // printf("stepType: %s %ld %d\n", stepType, len, tsteptype);
254
255
256
257
258
259
	}
    }

  return (tsteptype);
}

260
static
261
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
262
{
263
  int rdate, rtime;
264
265
266
  int timeUnits, startStep, endStep;
  int tstepRange = 0;
  int range;
267
  long lpar;
268
269
270
271
272
  long sigofrtime = 3;
  long editionNumber;

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

273
  if ( editionNumber > 1 )
274
275
276
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
277
278
279
280
281
282
283
284
285
286

  if ( sigofrtime == 3 )
    {
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      *vdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
      *vtime = (int) lpar*100;
    }
  else
    {
287
288
289
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      rdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
290
      rtime = (int) lpar*100;
291
      GRIB_CHECK(grib_get_long(gh, "forecastTime", &lpar), 0);
292
      startStep = (int) lpar;
293
      timeUnits = gribapiGetTimeUnits(gh);
294
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
295
296
297
298
299
300
301
302
303

      range = endStep - startStep;

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

304
305
306
307
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
308
	int time_period = endStep;
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
	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;
	      }
	  }

	julday_add_seconds(addsec, &julday, &secofday);

	decode_caldaysec(grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
	/*
	  printf("new %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute);
	*/
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
343
    }
344
345

  return (tstepRange);
346
347
}

348
static
349
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
{
351
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
  size_t datasize;
354
355
  long numberOfPoints;
  long lpar;
356
357
358

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

359
  gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360
361
362
363
364
365
366
367
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
368
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369

Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
371
372
  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
373
374
375
376
377
378
379
380
381
382
383
384
  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;

385
386
387
388
389
390
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
	if ( numberOfPoints != nlon*nlat )
392
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393
		(int)numberOfPoints, nlon*nlat);
394

395
396
397
398
399
400
401
402
403
404
405
	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
406
	if ( gridtype == GRID_LONLAT )
407
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408

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

411
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
	  {
413
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
	      {
415
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416

417
418
419
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
420
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
		      {
422
			double xinc = 360. / grid->xsize;
423

424
			if ( fabs(grid->xinc-xinc) > 0.0 )
425
			  {
426
427
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
428
			  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
430
431
		      }
		  }
	      }
432
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
	  }
434
	grid->ydef = 0;
435
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
	  {
437
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
	      {
439
440
441
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
	      }
443
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
445
446
447
448
	  }
	break;
      }
    case GRID_GAUSSIAN_REDUCED:
      {
449
450
451
452
	int nlat, i;
	size_t dummy;
	long *pl;

453
454
455
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

456
457
458
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

459
	grid->size   = numberOfPoints;
460

461
        grid->rowlon = (int *) malloc(nlat*sizeof(int));
462
463
464
        pl          = (long *) malloc(nlat*sizeof(long));
	dummy       = nlat;
	GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
465
	for ( i = 0; i < nlat; ++i ) grid->rowlon[i] = pl[i];
466
467
	free(pl);

468
469
470
471
472
473
474
475
476
	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);
477

478
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
479

480
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
	  {
482
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
	      {
484
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
485
486
487
488

		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
489
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
490
		      {
491
			double xinc = 360. / grid->xsize;
492

493
			if ( fabs(grid->xinc-xinc) > 0.0 )
494
			  {
495
496
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
497
498
499
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
	      }
501
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
	  }
503
504
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505
	  {
506
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
	      {
508
509
510
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
	      }
512
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513
514
515
	  }
	break;
      }
516
      /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
518
519
    case GRID_LCC:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
520
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
521
522
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);

523
524
525
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526

527
528
529
530
531
532
533
534
535
	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
536

537
	grid->xdef   = 0;
538
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
539
540
541

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
544
    case GRID_SPECTRAL:
      {
545
546
	size_t len = 256;
	char typeOfPacking[256];
547
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
548
549
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
550

551
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
553
	grid->trunc = lpar;
554

Uwe Schulzweida's avatar
Uwe Schulzweida committed
555
556
557
558
	break;
      }
    case GRID_GME:
      {
559
560
561
562
563
	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;
564
565
566

	break;
      }
567
    case GRID_REFERENCE:
568
      {
569
        char uuid[17];
570
571
572
573
	char reference_link[8192];
	size_t len = sizeof(reference_link);
	reference_link[0] = 0;

574
	grid->size  = numberOfPoints;
575
576
	if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
	  {
577
578
	    grid->number   = lpar;
	    if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
579
580
	    if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
	      {
581
		if ( strncmp(reference_link, "file://", 7) == 0 )
582
		  grid->reference = strdupx(reference_link);
583
	      }
584
585
586
587
588
            len = (size_t) 16;
            if ( grib_get_string(gh, "uuidOfHGrid", uuid, &len) == 0)
              {
                strncpy(grid->uuid, uuid, 16);
              }
589
	  }
590

Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
592
593
594
	break;
      }
    case GRID_GENERIC:
      {
595
596
597
598
599
	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;

600
	grid->size  = numberOfPoints;
601

602
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
603
	  {
604
605
	    grid->xsize = nlon;
	    grid->ysize = nlat;
606
607
608
	  }
	else
	  {
609
610
	    grid->xsize = 0;
	    grid->ysize = 0;
611
612
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
614
615
616
	break;
      }
    default:
      {
617
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
619
620
621
	break;
      }
    }

622
  grid->isRotated = FALSE;
623
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
624
    {
625
626
627
628
      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);
629
      /* change from south to north pole */
630
631
632
633
634
635
636
637
638
639
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

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

static
640
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
641
{
642
643
644
  int status;
  long lpar;
  double dlevel;
645

646
647
648
649
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
650

651
652
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
653
    {
654
      *leveltype = (int) lpar;
655

656
      switch (*leveltype)
657
658
659
660
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
661
	  { *lbounds = 1; break; }
662
	}
663

664
665
666
667
668
669
      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;
670

671
672
673
674
675
676
677
678
679
680
	  *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;
	}
681
    }
682
683
684
}

static
685
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
686
{
687
688
689
690
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
691
692
  double dlevel;

693
694
695
696
697
698
699
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
700
    {
701
      *leveltype = (int) lpar;
702

703
704
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
705

706
      if ( *leveltype == leveltype2 && *leveltype != 255 ) *lbounds = 1;
707

708
      if ( *lbounds == 0 )
709
	{
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	    }
	  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;
	    }

726
	  *level1 = (int) dlevel;
727
	  *level2 = 0;
728
729
730
	}
      else
	{
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	      *level1 = (int) dlevel;
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfSecondFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	      *level2 = (int) dlevel;
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	      *level1 = lpar;
	      GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	      *level2 = lpar;
	    }
753
754
	}
    }
755
756
757
758
}

static
void gribapiAddRecord(int streamID, int param, grib_handle *gh,
759
		      long recsize, off_t position, int datatype, int comptype)
760
761
762
763
764
765
{
  long editionNumber;
  int zaxistype;
  int gridID = CDI_UNDEFID, varID;
  int levelID = 0;
  int tsID, recID;
766
  int level1 = 0, level2 = 0;
767
768
769
770
771
772
773
774
775
776
777
  int numavg;
  int tsteptype;
  int lbounds = 0;
  record_t *record;
  grid_t grid;
  int vlistID;
  stream_t *streamptr;
  int leveltype;
  long lpar;
  int status;
  char name[256], longname[256], units[256];
778
  size_t vlen;
779
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
780
781
782
783
784

  streamptr = stream_to_pointer(streamID);

  vlistID = streamInqVlist(streamID);
  tsID    = streamptr->curTsID;
785
  recID   = recordNewEntry(streamptr, tsID);
786
787
788
789
790
791
792
793
794
  record  = &streamptr->tsteps[tsID].records[recID];

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

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

  if ( editionNumber <= 1 )
795
    grib1GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
796
  else
797
    grib2GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
798

799
800
801
802
803
804
805
806
807
808
  // 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;

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
810
811

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
812
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
813

814
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
815
    {
816
817
818
819
820
821
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
822

823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
        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];
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
842
        if ( lpar != 3 )
843
844
845
846
847
848
849
850
851
852
853
854
          {
            fprintf(stderr, "Warning ...\n");
          }
        GRIB_CHECK(grib_get_long(gh, "nlev", &nlev), 0);
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &nvgrid), 0);

        len = (size_t) 16;
        uuid[16] = 0;
        GRIB_CHECK(grib_get_string(gh, "uuidOfVGrid", uuid, &len), 0);
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
856

857
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
858
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
859

860
861
862
863
  name[0] = 0;
  longname[0] = 0;
  units[0] = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
864
865
  vlen = 256;
  GRIB_CHECK(grib_get_string(gh, "shortName", name, &vlen), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
866
867
868
869
870
871
872
873
874
875
876
877
  if      ( vlen == 8 && memcmp(name, "unknown", vlen) == 0 ) name[0] = 0;
  else if ( vlen == 2 && memcmp(name, "~", vlen)       == 0 ) name[0] = 0;

  if ( name[0] != 0 )
    {
      vlen = 256;
      GRIB_CHECK(grib_get_string(gh, "name", longname, &vlen), 0);
      if ( vlen == 8 && memcmp(longname, "unknown", vlen) == 0 ) longname[0] = 0;
      vlen = 256;
      GRIB_CHECK(grib_get_string(gh, "units", units, &vlen), 0);
      if ( vlen == 8 && memcmp(units, "unknown", vlen) == 0 ) units[0] = 0;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
878
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
879

Uwe Schulzweida's avatar
Uwe Schulzweida committed
880
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2,
881
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
882
	       name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
883
884
885
886

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

887
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888

889
890
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
891
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
892
  */
893
894
895
896
897
898
899
900
901
  {
    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);
      }
  }
902
903
904
905
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);


Uwe Schulzweida's avatar
Uwe Schulzweida committed
906
907
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
908
909
910
911
912
      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
913
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
914
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
915
916
917
918
919
920
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921
      long processID;
922
923
924
925
926
927
928
929
      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
930
931
932
933
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
938
939
940
941
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
945
946
947
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
948
949
950
951
952
953
954
	}
    }

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

  if ( CDI_Debug )
955
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
959
}
#endif

960
int gribapiScanTimestep1(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
961
962
963
964
965
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
967
968
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
969
  int rtabnum = 0;
970
971
  int rcode = 0, level1 = 0, level2 = 0;
  int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
973
974
975
976
977
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
978
  int datatype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
979
980
  long recsize = 0;
  int warn_time = TRUE;
981
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
982
983
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
984
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
  int vlistID;
986
  int comptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
987
988
989
990
991
  long unzipsize;
  compvar2_t compVar, compVar0;
  stream_t *streamptr;
  grib_handle *gh = NULL;
  int leveltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
992
  int pdis = 0, pcat = 0, pnum = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
994
  long editionNumber;
  long lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
  int bitsPerValue;
996
  int lieee = FALSE;
997
  int lbounds;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
998
999
1000

  streamptr = stream_to_pointer(streamID);

1001
  stream_check_ptr(__func__, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
1003
1004

  streamptr->curTsID = 0;

1005
  tsID  = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1006
1007
1008
  taxis = &streamptr->tsteps[tsID].taxis;

  if ( tsID != 0 )
1009
    Error("Internal problem! tstepsNewEntry returns %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1010
1011
1012
1013
1014
1015

  fileID = streamInqFileID(streamID);

  nrecs = 0;
  while ( TRUE )
    {
1016
1017
      level1 = 0;
      level2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
	  streamptr->ntsteps = 1;
	  break;
	}
      if ( recsize > buffersize )
	{
	  buffersize = recsize;
	  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	}

      readsize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
1034
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035

1036
1037
      lieee = FALSE;

1038
      comptype = COMPRESS_NONE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1039
1040
      if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	{
1041
	  comptype = COMPRESS_SZIP;
1042
	  unzipsize += 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1043
1044
1045
1046
1047
1048
1049
1050
	  if ( (long) buffersize < unzipsize )
	    {
	      buffersize = unzipsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
	}

      gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
1051
      GRIB_CHECK(grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1052
1053
1054
1055
1056

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

      if ( editionNumber <= 1 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1057
1058
	  GRIB_CHECK(grib_get_long(gh, "table2Version", &lpar), 0);
	  rtabnum = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1059
1060
	  GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &lpar), 0);
	  rcode = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1061

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1062
	  param = cdiEncodeParam(rcode, rtabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1063

1064
	  grib1GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1065
1066
1067
	}
      else
	{
1068
1069
	  size_t len = 256;
	  char typeOfPacking[256];
1070

1071
	  status = grib_get_string(gh, "packingType", typeOfPacking, &len);
1072
1073
	  if ( status == 0 )
	    {
Uwe Schulzweida's avatar