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

#include <stdio.h>

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


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

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

extern int cdiInventoryMode;

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


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

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

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

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

  return (gridtype);
}

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

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

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

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

  return (isRotated);
}

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

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

  return (zaxistype);
}

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

  switch (unitsOfTime)
    {
    case 13:  timeunits = TUNIT_SECOND;  break;
    case  0:  timeunits = TUNIT_MINUTE;  break;
    case  1:  timeunits = TUNIT_HOUR;    break;
    case 10:  timeunits = TUNIT_3HOURS;  break;
    case 11:  timeunits = TUNIT_6HOURS;  break;
    case 12:  timeunits = TUNIT_12HOURS; break;
    case  2:  timeunits = TUNIT_DAY;     break;
    default:  timeunits = TUNIT_HOUR;    break;
    }

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
  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;
155
  int status;
156
157
158
  long unitsOfTime;
  long lpar;

159
  status = grib_get_long(gh, "stepUnits", &unitsOfTime);
160
161
162

  timeunits2 = getTimeunits(unitsOfTime);

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

165
166
  if ( status == 0 )
    endStep = (int) ((lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
167
168
169
170
171
172
173
174

  return (endStep);
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  int timeunits = -1;
175
176
  long unitsOfTime = -1;
  int status;
177
178
179
180
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

181
  status = grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
182
183
184

  timeunits = getTimeunits(unitsOfTime);

185
  /*
186
187
  GRIB_CHECK(grib_get_string(gh, "stepUnits", stepunits, &len), 0);

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

  return (timeunits);
}

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

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

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

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

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

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

267
static
268
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
269
{
270
  int rdate, rtime;
271
  int timeUnits, startStep = 0, endStep;
272
273
  int tstepRange = 0;
  int range;
274
  int status;
275
  long lpar;
276
277
278
279
280
  long sigofrtime = 3;
  long editionNumber;

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

281
  if ( editionNumber > 1 )
282
283
284
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
285
286
287
288
289
290
291
292
293
294

  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
    {
295
296
297
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      rdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
298
      rtime = (int) lpar*100;
299
300
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
301
      timeUnits = gribapiGetTimeUnits(gh);
302
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
303
304
305
306
307
308
309
310
311

      range = endStep - startStep;

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

312
313
314
315
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
316
	int time_period = endStep;
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
	int julday, secofday, addsec;

	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

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

	addsec = 0;
	switch ( timeUnits )
	  {
	  case TUNIT_SECOND:  addsec =         time_period; break;
	  case TUNIT_MINUTE:  addsec =    60 * time_period; break;
	  case TUNIT_HOUR:    addsec =  3600 * time_period; break;
	  case TUNIT_3HOURS:  addsec = 10800 * time_period; break;
	  case TUNIT_6HOURS:  addsec = 21600 * time_period; break;
	  case TUNIT_12HOURS: addsec = 43200 * time_period; break;
	  case TUNIT_DAY:     addsec = 86400 * time_period; break;
	  default:
	    if ( lprint )
	      {
	        Warning("Time unit %d unsupported", timeUnits);
		lprint = FALSE;
	      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
	    break;
341
342
343
344
345
346
347
348
349
350
351
	  }

	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);
      }
352
    }
353
354

  return (tstepRange);
355
356
}

357
static
358
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
{
360
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
362
  size_t datasize;
363
364
  long numberOfPoints;
  long lpar;
365
366
367

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

368
  gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
370
371
372
373
374
375
376
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
377
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378

Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
380
381
  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
382
383
384
385
386
387
388
389
390
391
392
393
  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;

394
395
396
397
398
399
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

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

403
404
405
406
407
408
409
410
411
412
413
	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
414
	if ( gridtype == GRID_LONLAT )
415
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416

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

419
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
	  {
421
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
422
	      {
423
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
424

425
426
427
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
428
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
		      {
430
			double xinc = 360. / grid->xsize;
431

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

461
462
463
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

464
465
466
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

467
	grid->size   = numberOfPoints;
468

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

476
477
478
479
480
481
482
483
484
	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);
485

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

488
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
	  {
490
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
	      {
492
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
493
494
495
496

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

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

531
532
533
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534

535
536
537
538
539
540
541
542
543
	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
544

545
	grid->xdef   = 0;
546
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547
548
549

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
552
    case GRID_SPECTRAL:
      {
553
554
	size_t len = 256;
	char typeOfPacking[256];
555
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
556
557
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
558

559
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
561
	grid->trunc = lpar;
562

Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
564
565
566
	break;
      }
    case GRID_GME:
      {
567
568
569
570
571
	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;
572
573
574

	break;
      }
575
    case GRID_REFERENCE:
576
      {
577
        char uuid[17];
578
579
580
    	char reference_link[8192];
	    size_t len = sizeof(reference_link);
	    reference_link[0] = 0;
581

582
583
    	grid->size  = numberOfPoints;
	    if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
584
	      {
585
586
587
588
589
590
591
592
593
594
595
596
	        grid->number   = lpar;
	        if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
	        if ( grib_get_bytes(gh, "gridDescriptionFile", (unsigned char *) reference_link, &len) == 0 )
	          {
		        if ( strncmp(reference_link, "file://", 7) == 0 )
		          grid->reference = strdupx(reference_link);
	            }
                len = (size_t) 16;
                if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
                  {
                    strncpy(grid->uuid, uuid, 16);
                  }
597
	      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
599
600
601
	break;
      }
    case GRID_GENERIC:
      {
602
603
604
605
606
	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;

607
	grid->size  = numberOfPoints;
608

609
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
610
	  {
611
612
	    grid->xsize = nlon;
	    grid->ysize = nlat;
613
614
615
	  }
	else
	  {
616
617
	    grid->xsize = 0;
	    grid->ysize = 0;
618
619
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
620
621
622
623
	break;
      }
    default:
      {
624
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
625
626
627
628
	break;
      }
    }

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

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

static
647
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
648
{
649
650
651
  int status;
  long lpar;
  double dlevel;
652

653
654
655
656
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
657

658
659
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
660
    {
661
      *leveltype = (int) lpar;
662

663
      switch (*leveltype)
664
665
666
667
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
668
	  { *lbounds = 1; break; }
669
	}
670

671
672
673
674
675
676
      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;
677

678
679
680
681
682
683
684
685
686
687
	  *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;
	}
688
    }
689
690
691
}

static
692
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf)
693
{
694
695
696
697
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
698
699
  double dlevel;

700
701
702
703
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
704
  *level_sf = 0;
705
706
707

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
708
    {
709
      *leveltype = (int) lpar;
710

711
712
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
713

714
      if ( *leveltype == leveltype2 && *leveltype != 255 ) *lbounds = 1;
715

716
      if ( *lbounds == 0 )
717
	{
718
719
720
721
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
722
723
724
725
726
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
              *level_sf = 77;
727
728
729
730
731
732
733
734
735
	    }
	  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;
	    }

736
	  *level1 = (int) dlevel;
737
	  *level2 = 0;
738
739
740
	}
      else
	{
741
742
743
744
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
745
746
747
748
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
749
750
751
	      *level1 = (int) dlevel;
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfSecondFixedSurface", &dlevel), 0);
752
753
754
755
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
756
	      *level2 = (int) dlevel;
757
              *level_sf = 77;
758
759
760
761
762
763
764
765
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	      *level1 = lpar;
	      GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	      *level2 = lpar;
	    }
766
767
	}
    }
768
769
}

770
771
772
773
774
775
776
777
778
779
static
void gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length)
{
  string[0] = 0;

  GRIB_CHECK(grib_get_string(gh, key, string, &length), 0);
  if      ( length == 8 && memcmp(string, "unknown", length) == 0 ) string[0] = 0;
  else if ( length == 2 && memcmp(string, "~", length)       == 0 ) string[0] = 0;
}

780
static
781
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
782
783
		      long recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
                      int leveltype, int lbounds, int level1, int level2, int level_sf)
784
785
786
787
788
789
790
791
792
793
794
795
796
{
  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;
797
  char longname[256], units[256];
798
  size_t vlen;
799
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
800

Uwe Schulzweida's avatar
Uwe Schulzweida committed
801
  vlistID = streamptr->vlistID;
802
  tsID    = streamptr->curTsID;
803
  recID   = recordNewEntry(streamptr, tsID);
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
  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;
820
  memcpy((*record).varname, varname, len);
821
822

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
823
824
825

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827

828
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
    {
830
831
832
833
834
835
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836

837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
        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];
853
        double dtmp;
854
855
856
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
857
        if ( lpar != 6 )
858
859
860
          {
            fprintf(stderr, "Warning ...\n");
          }
861
862
863
864
        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);
865
866
        len = (size_t) 16;
        uuid[16] = 0;
867
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
868
869
870
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
871
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872

873
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
874
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
875

876
877
878
  longname[0] = 0;
  units[0] = 0;

879
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
880
881
    {
      vlen = 256;
882
      gribapiGetString(gh, "name", longname, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
883
      vlen = 256;
884
      gribapiGetString(gh, "units", units, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885
    }
886
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887

888
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf,
889
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
890
	       varname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
891
892
893
894

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

895
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
896

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


Uwe Schulzweida's avatar
Uwe Schulzweida committed
914
915
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
916
917
918
919
920
      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
921
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
922
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
923
924
925
926
927
928
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
      long processID;
930
931
932
933
934
935
936
937
      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
938
939
940
941
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
946
947
948
949
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
950

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
953
954
955
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956
957
958
959
960
961
962
	}
    }

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

  if ( CDI_Debug )
963
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
964
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
965
}
966
967
968
969
970
971

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
972
  int status;
973
974
975
976
977
  long lpar;

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

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

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

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

986
987
  return (param);
}
988
989

static
990
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
991
992
{
  compvar2_t compVar;
993
994
995
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
996
997
998
999
1000

  compVar.param  = param;
  compVar.level1 = level1;
  compVar.level2 = level2;
  compVar.ltype  = leveltype;
1001
1002
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
1003
1004
1005
1006
1007
1008
1009
1010
1011

  return (compVar);
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record)
{
  int rstatus;
  compvar2_t compVar0;
1012
  size_t maxlen = sizeof(compVar.name);
1013
1014
1015
1016
1017

  compVar0.param  = record.param;
  compVar0.level1 = record.ilevel;
  compVar0.level2 = record.ilevel2;
  compVar0.ltype  = record.ltype;
1018
  memcpy(compVar0.name, record.varname, maxlen);
1019
1020
1021
1022
1023

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

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
#endif

1026
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
1028
1029
1030
1031
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
1034
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035
  int rtabnum = 0;
1036
1037
  int rcode = 0, level1 = 0, level2 = 0;
  int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1038
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1039
1040
1041
1042
1043
  DateTime datetime, datetime0;