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
  int gridtype = GRID_GENERIC;
40
  int gribgridtype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
  long lpar;

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

      if ( status ==  0 ) gribgridtype = (int) lpar;
48
49
50

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

  return (gridtype);
}

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

80
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
81

82
83
84
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
85
86
87
88

  return (isRotated);
}

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

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

  return (zaxistype);
}

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

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

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

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

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

  return (timeunits);
}

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

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

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

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

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
220
221
222
223
224
225
226

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

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

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

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

  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
    {
289
290
291
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      rdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
292
      rtime = (int) lpar*100;
293
      GRIB_CHECK(grib_get_long(gh, "forecastTime", &lpar), 0);
294
      startStep = (int) lpar;
295
      timeUnits = gribapiGetTimeUnits(gh);
296
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
297
298
299
300
301
302
303
304
305

      range = endStep - startStep;

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

306
307
308
309
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
310
	int time_period = endStep;
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
343
344
	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);
      }
345
    }
346
347

  return (tstepRange);
348
349
}

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

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

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

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

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

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

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

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

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

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

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

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

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

460
	grid->size   = numberOfPoints;
461

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

601
	grid->size  = numberOfPoints;
602

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

709
      if ( *lbounds == 0 )
710
	{
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
	  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;
	    }

727
	  *level1 = (int) dlevel;
728
	  *level2 = 0;
729
730
731
	}
      else
	{
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
	  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;
	    }
754
755
	}
    }
756
757
758
759
}

static
void gribapiAddRecord(int streamID, int param, grib_handle *gh,
760
		      long recsize, off_t position, int datatype, int comptype)
761
762
763
764
765
766
{
  long editionNumber;
  int zaxistype;
  int gridID = CDI_UNDEFID, varID;
  int levelID = 0;
  int tsID, recID;
767
  int level1 = 0, level2 = 0;
768
769
770
771
772
773
774
775
776
777
778
  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];
779
  size_t vlen;
780
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
781
782
783

  streamptr = stream_to_pointer(streamID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
784
  vlistID = streamptr->vlistID;
785
  tsID    = streamptr->curTsID;
786
  recID   = recordNewEntry(streamptr, tsID);
787
788
789
790
791
792
793
794
795
  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 )
796
    grib1GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
797
  else
798
    grib2GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
799

800
801
802
803
804
805
806
807
808
809
  // 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
810
811
812

  gridID = varDefGrid(vlistID, grid, 0);

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

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

824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
        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);
843
        if ( lpar != 3 )
844
845
846
847
848
849
850
851
852
853
854
855
          {
            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
856
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
865
866
  vlen = 256;
  GRIB_CHECK(grib_get_string(gh, "shortName", name, &vlen), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
867
868
869
870
871
872
873
874
875
876
877
878
  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
879
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
880

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

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

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

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


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

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

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

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

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

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

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

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

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

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

  streamptr = stream_to_pointer(streamID);

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

  streamptr->curTsID = 0;

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1012
  fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1013
1014
1015
1016

  nrecs = 0;
  while ( TRUE )
    {
1017
1018
      level1 = 0;
      level2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
      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
1034
1035
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1036

1037
1038
      lieee = FALSE;

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

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

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

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