stream_gribapi.c 84.2 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
5
6
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <stdio.h>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
#include "cdi.h"
9
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
11
12
13
#include "file.h"
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
14
#include "stream_grb.h"
15
#include "calendar.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
17
18


#if  defined  (HAVE_LIBGRIB_API)
19
#  include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip, GRIB1_LTYPE_99 */
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
34
  char name[32];
} 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
45
  long lpar;

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

      if ( status ==  0 ) gribgridtype = (int) lpar;
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

      switch (gribgridtype)
	{
	case  GRIB2_GTYPE_LATLON:     { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
	                                if ( lpar == (long) GRIB_MISSING_LONG ) break;
                                      }
	case  GRIB2_GTYPE_LATLON_ROT: { gridtype = GRID_LONLAT;    break; }
	case  GRIB2_GTYPE_LCC:        { gridtype = GRID_LCC;       break; }
	case  GRIB2_GTYPE_GAUSSIAN:   { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
	                                if ( lpar == (long) GRIB_MISSING_LONG )
					  gridtype = GRID_GAUSSIAN_REDUCED;
					else
					  gridtype = GRID_GAUSSIAN;
					break;
                                      }
	case  GRIB2_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL;  break; }
	case  GRIB2_GTYPE_GME:        { gridtype = GRID_GME;       break; }
	case  GRIB2_GTYPE_NUMBER:     { gridtype = GRID_REFERENCE; break; }
	}
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
{
94
  int zaxistype = ZAXIS_GENERIC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

96
  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
    {
98
99
100
101
102
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
    }
  else
    {
      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
154
155
156
157
158
159
160
161
162
163
  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;
  long unitsOfTime = -1;
  int status;
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

  status = grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);

  timeunits = getTimeunits(unitsOfTime);

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

  len--;
168

169
170
171
172
173
174
175
176
177
178
  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;
  else if ( lprint )
179
    {
180
181
      Message("Step units >%s< unsupported!", stepunits);
      lprint = FALSE;
182
    }
183
  */
184
185
186
187

  return (timeunits);
}

188
189
190
191
static
int gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
192
  int timeunits2 = timeunits;
193
194
195
196
  int status;
  long unitsOfTime;
  long lpar;

197
198
199
  status = grib_get_long(gh, "stepUnits", &unitsOfTime);
  if ( status == 0 ) timeunits2 = getTimeunits(unitsOfTime);
  //timeunits2 = gribapiGetTimeUnits(gh);
200
201
202
203
204
205
206
207
208
209

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

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

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

218
  if ( editionNumber > 1 )
219
220
221
222
    {
      long sigofrtime;

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

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

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
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;
}

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

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

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

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

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

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

      range = endStep - startStep;

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
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;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378

379
  GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380

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
	      {
Deike Kleberg's avatar
Deike Kleberg committed
436
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437

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

445
446
447
448
449
			if ( fabs(grid->xinc-xinc) > 0.0 )
			  {
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
			  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
450
451
452
		      }
		  }
	      }
Deike Kleberg's avatar
Deike Kleberg committed
453
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
	  }
Deike Kleberg's avatar
Deike Kleberg committed
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
	      }
Deike Kleberg's avatar
Deike Kleberg committed
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
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

	grid->size   = numberOfPoints;

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

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

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

	/* 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
506
507
508
509
510
511
512
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;

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

514
515
516
517
518
519
520
			if ( fabs(grid->xinc-xinc) > 0.0 )
			  {
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
			  }
		      }
		  }
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
568
569
570
571
572
	size_t len = 256;
	char typeOfPacking[256];
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
577
578
579
	break;
      }
    case GRID_GME:
      {
580
581
582
583
584
585
586
587
588
589
	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;

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

595
596
597
598
599
600
601
602
603
604
    	grid->size  = numberOfPoints;
        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);
              }
605
            len = (size_t) 16;
606
            if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
607
              {
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
621
622
	grid->size  = numberOfPoints;

	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
643
  grid->isRotated = FALSE;
  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
  int status;
  long lpar;
664
665
  double dlevel;

666
667
668
669
670
671
672
673
674
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;

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

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

      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;
689
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
690
691
692
693
694
695
696
697
698
699
700
701

	  *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;
	}
    }
702
703
704
}

static
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
double grib2ScaleFactor(long factor)
{
  double scaleFactor = 0;

  if      ( factor == 0 ) scaleFactor =    1;
  else if ( factor == 1 ) scaleFactor =    0.1;
  else if ( factor == 2 ) scaleFactor =    0.01;
  else if ( factor == 3 ) scaleFactor =    0.001;
  else if ( factor == 4 ) scaleFactor =    0.0001;

  return (scaleFactor);
}

static
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf, int *level_unit)
720
{
721
722
723
724
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
725

726
727
728
729
730
731
  *leveltype  = 0;
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
732
733
734
735

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
    {
736
737
738
      long llevel;
      double dlevel1 = 0, dlevel2 = 0;

739
740
741
742
      *leveltype = (int) lpar;

      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
743

744
      if ( *leveltype != 255 && leveltype2 != 255 && leveltype2 > 0 ) *lbounds = 1;
745

746
747
748
749
750
751
752
753
754
755
      if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_M;
        }
      else if ( *leveltype == GRIB2_LTYPE_ISOBARIC )
        {
          *level_sf = 1000;
          *level_unit = CDI_UNIT_PA;
        }
756

757
758
759
760
761
762
763
764
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
      if ( llevel != GRIB_MISSING_LONG && factor != GRIB_MISSING_LONG )
        dlevel1 = llevel*grib2ScaleFactor(factor);

      if ( *level_sf != 0 ) dlevel1 *= (*level_sf);

      if ( *lbounds == 1 )
765
	{
766
767
768
769
770
771
772
773
774
775
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
          if ( llevel != GRIB_MISSING_LONG && factor != GRIB_MISSING_LONG )
            dlevel2 = llevel*grib2ScaleFactor(factor);

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
776
    }
777
778
779
}

static
780
781
782
783
784
785
786
787
788
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;
}

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

811
  vlistID = streamptr->vlistID;
812
  tsID    = streamptr->curTsID;
813
  recID   = recordNewEntry(streamptr, tsID);
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
  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;
830
  memcpy((*record).varname, varname, len);
831
832

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
834
835

  gridID = varDefGrid(vlistID, grid, 0);

836
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
837

838
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
    {
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;

        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];
863
        double dtmp;
864
865
866
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
867
        if ( lpar != 6 )
868
869
870
          {
            fprintf(stderr, "Warning ...\n");
          }
871
872
873
874
        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);
875
876
        len = (size_t) 16;
        uuid[16] = 0;
877
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
878
879
880
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
881
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
882

883
884
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885

886
887
  longname[0] = 0;
  units[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888

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

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

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

Deike Kleberg's avatar
Deike Kleberg committed
905
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
906

907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
  */
  {
    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);
      }
  }
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
  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! */
    }
944

Uwe Schulzweida's avatar
Uwe Schulzweida committed
945
946
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
947
948
949
950
951
      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
952
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
954
955
956
957
958
959
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960
      long processID;
961
962
963
964
965
966
967
968
      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
969
970
971
972
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
977
978
979
980
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
981

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

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

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

  if ( CDI_Debug )
994
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
996
}
997
#endif
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
  int status;
  long lpar;

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

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

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

  param = cdiEncodeParam(pnum, pcat, pdis);

  return (param);
}

static
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
{
  compvar2_t compVar;
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;

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

  return (compVar);
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record)
{
  int rstatus;
  compvar2_t compVar0;
  size_t maxlen = sizeof(compVar.name);

  compVar0.param  = record.param;
  compVar0.level1 = record.ilevel;
  compVar0.level2 = record.ilevel2;
  compVar0.ltype  = record.ltype;
  memcpy(compVar0.name, record.varname, maxlen);

  rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1056
1057
#endif

1058
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1059
1060
1061
1062
1063
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
1064
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1065
1066
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1067
  int rtabnum = 0;
1068
1069
  int rcode = 0, level1 = 0, level2 = 0;
  int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1070
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1071
1072
1073
1074
1075
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
1076
  int nrecs_scanned;
1077
  int datatype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1078
1079
  long recsize = 0;
  int warn_time = TRUE;
1080
  // int warn_numavg = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1081
1082
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
1083
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1084
  int vlistID;
Deike Kleberg's avatar
Deike Kleberg committed
1085
  int comptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1086
  long unzipsize;
1087
  compvar2_t compVar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1088
1089
1090
1091
  grib_handle *gh = NULL;
  int leveltype;
  long editionNumber;
  long lpar;
1092
  size_t len;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1093
  int bitsPerValue;
1094
1095
  int lieee = FALSE;
  int lbounds;
1096
  int level_sf, level_unit;
1097
1098
  char paramstr[32];
  char varname[256];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1099
1100
1101

  streamptr->curTsID = 0;

1102
  tsID  = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1103
1104
1105
  taxis = &streamptr->tsteps[tsID].taxis;

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

1108
  fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1109

1110
  nrecs_scanned = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1111
1112
1113
  nrecs = 0;
  while ( TRUE )
    {
1114
1115
      level1 = 0;
      level2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
      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;
1131
1132
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1133

1134
1135
      lieee = FALSE;

Deike Kleberg's avatar
Deike Kleberg committed
1136
      comptype = COMPRESS_NONE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1137
1138
      if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	{
Deike Kleberg's avatar
Deike Kleberg committed
1139
	  comptype = COMPRESS_SZIP;
1140
	  unzipsize += 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1141
1142
1143
1144
1145
1146
1147
	  if ( (long) buffersize < unzipsize )
	    {
	      buffersize = unzipsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
	}

1148
      nrecs_scanned++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1149
      gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
1150
      GRIB_CHECK(grib_set_double(gh, "missingValue", cdiDefaultMissval), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1151
1152
1153
1154
1155

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

      if ( editionNumber <= 1 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1156
1157
	  GRIB_CHECK(grib_get_long(gh, "table2Version", &lpar), 0);
	  rtabnum = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1158
1159
	  GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &lpar), 0);
	  rcode = (int) lpar;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1160

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

1163
	  grib1GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
1164
          level_sf = 0;
1165
          level_unit = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1166
1167
1168
	}
      else
	{
1169
1170
	  size_t len = 256;
	  char typeOfPacking[256];
1171
1172
1173
1174
1175

	  status = grib_get_string(gh, "packingType", typeOfPacking, &len);
	  if ( status == 0 )
	    {
	      // fprintf(stderr, "packingType %d %s\n", len, typeOfPacking);
1176
1177
	      if      ( strncmp(typeOfPacking, "grid_jpeg", len) == 0 ) comptype = COMPRESS_JPEG;
	      else if ( strncmp(typeOfPacking, "grid_ieee", len) == 0 ) lieee = TRUE;
1178
	    }
1179

1180
	  param = gribapiGetParam(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1181

1182
	  grib2GetLevel(gh, &leveltype, &lbounds, &level1, &level2, &level_sf, &level_unit);
1183
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1184

1185
      cdiParamToString(param, paramstr, sizeof(paramstr));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1186

1187
1188
1189
1190
1191
      varname[0] = 0;
      gribapiGetString(gh, "shortName", varname, sizeof(varname));
      len = strlen(varname);
      if ( len > 32 ) len = 32;
      //printf("param = %s  name = %s   l1 = %d  l2 = %d\n", paramstr, varname, level1, level2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1192

1193
1194
      gribapiGetValidityDateTime(gh, &vdate, &vtime);
      /*
1195
      printf("%d %d %d.%d.%d %d\n", vdate, vtime, pnum, pcat, pdis, leveltype);
1196
      */
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207