stream_gribapi.c 85.1 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
268
269
270
271
272
273
274
275
276
277
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;
}

278
static
279
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
280
{
281
  int rdate, rtime;
282
  int timeUnits, startStep = 0, endStep;
283
284
  int tstepRange = 0;
  int range;
285
  int status;
286
  long lpar;
287
288
289
290
291
  long sigofrtime = 3;
  long editionNumber;

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

292
  if ( editionNumber > 1 )
293
294
295
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
296
297
298

  if ( sigofrtime == 3 )
    {
299
      gribapiGetDataDateTime(gh, vdate, vtime);
300
301
302
    }
  else
    {
303
304
      gribapiGetDataDateTime(gh, &rdate, &rtime);

305
306
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
307
      timeUnits = gribapiGetTimeUnits(gh);
308
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
309
310
311
312
313
314
315
316
317

      range = endStep - startStep;

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

318
319
320
321
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
322
	int time_period = endStep;
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
	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
346
	    break;
347
348
349
350
351
352
353
354
355
356
357
	  }

	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);
      }
358
    }
359
360

  return (tstepRange);
361
362
}

363
static
364
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
{
366
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
367
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368
  size_t datasize;
369
370
  long numberOfPoints;
  long lpar;
371
372
373

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

374
  gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
375
376
377
378
379
380
381
382
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
383
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384

Uwe Schulzweida's avatar
Uwe Schulzweida committed
385
386
387
  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
388
389
390
391
392
393
394
395
396
397
398
399
  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;

400
401
402
403
404
405
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

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

409
410
411
412
413
414
415
416
417
418
419
	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
420
	if ( gridtype == GRID_LONLAT )
421
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
422

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

425
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
	  {
427
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
	      {
429
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430

431
432
433
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
434
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
		      {
436
			double xinc = 360. / grid->xsize;
437

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

467
468
469
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

470
471
472
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

473
	grid->size   = numberOfPoints;
474

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

482
483
484
485
486
487
488
489
490
	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);
491

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

494
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
	  {
496
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
	      {
498
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
499
500
501
502

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

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

537
538
539
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540

541
542
543
544
545
546
547
548
549
	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
550

551
	grid->xdef   = 0;
552
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
553
554
555

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557
558
    case GRID_SPECTRAL:
      {
559
560
	size_t len = 256;
	char typeOfPacking[256];
561
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
562
563
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
564

565
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
567
	grid->trunc = lpar;
568

Uwe Schulzweida's avatar
Uwe Schulzweida committed
569
570
571
572
	break;
      }
    case GRID_GME:
      {
573
574
575
576
577
	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;
578
579
580

	break;
      }
581
    case GRID_REFERENCE:
582
      {
583
        char uuid[17];
584
    	char reference_link[8192];
585
586
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
587

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

613
	grid->size  = numberOfPoints;
614

615
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
616
	  {
617
618
	    grid->xsize = nlon;
	    grid->ysize = nlat;
619
620
621
	  }
	else
	  {
622
623
	    grid->xsize = 0;
	    grid->ysize = 0;
624
625
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
627
628
629
	break;
      }
    default:
      {
630
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631
632
633
634
	break;
      }
    }

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

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

static
653
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
654
{
655
656
657
  int status;
  long lpar;
  double dlevel;
658

659
660
661
662
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
663

664
665
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
666
    {
667
      *leveltype = (int) lpar;
668

669
      switch (*leveltype)
670
671
672
673
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
674
	  { *lbounds = 1; break; }
675
	}
676

677
678
679
680
681
682
      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;
683

684
685
686
687
688
689
690
691
692
693
	  *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;
	}
694
    }
695
696
697
}

static
698
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf)
699
{
700
701
702
703
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
704
705
  double dlevel;

706
707
708
709
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
710
  *level_sf = 0;
711
712
713

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
714
    {
715
      *leveltype = (int) lpar;
716

717
718
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
719

720
      if ( *leveltype == leveltype2 && *leveltype != 255 ) *lbounds = 1;
721

722
      if ( *lbounds == 0 )
723
	{
724
725
726
727
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
728
729
730
731
732
	      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;
733
734
735
736
737
738
739
740
741
	    }
	  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;
	    }

742
	  *level1 = (int) dlevel;
743
	  *level2 = 0;
744
745
746
	}
      else
	{
747
748
749
750
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
751
752
753
754
	      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
755
756
757
	      *level1 = (int) dlevel;
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfSecondFixedSurface", &dlevel), 0);
758
759
760
761
	      if      ( factor == 0 ) dlevel *= 1000;  //  m to mm
	      else if ( factor == 1 ) dlevel *=  100;  // dm to mm
	      else if ( factor == 2 ) dlevel *=   10;  // cm to mm
	      else if ( factor == 3 ) dlevel *=    1;  // mm to mm
762
	      *level2 = (int) dlevel;
763
              *level_sf = 77;
764
765
766
767
768
769
770
771
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	      *level1 = lpar;
	      GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	      *level2 = lpar;
	    }
772
773
	}
    }
774
775
}

776
777
778
779
780
781
782
783
784
785
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;
}

786
static
787
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
788
789
		      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)
790
791
792
793
794
795
796
797
798
799
800
801
802
{
  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;
803
  char longname[256], units[256];
804
  size_t vlen;
805
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
806

Uwe Schulzweida's avatar
Uwe Schulzweida committed
807
  vlistID = streamptr->vlistID;
808
  tsID    = streamptr->curTsID;
809
  recID   = recordNewEntry(streamptr, tsID);
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
  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;
826
  memcpy((*record).varname, varname, len);
827
828

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
830
831

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
832
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
833

834
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
835
    {
836
837
838
839
840
841
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
842

843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
        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];
859
        double dtmp;
860
861
862
        long nlev, nvgrid;

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

879
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
880
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
881

882
883
884
  longname[0] = 0;
  units[0] = 0;

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

894
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf,
895
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
896
	       varname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
897
898
899
900

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

901
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
902

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


Uwe Schulzweida's avatar
Uwe Schulzweida committed
920
921
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
922
923
924
925
926
      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
927
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
930
931
932
933
934
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
935
      long processID;
936
937
938
939
940
941
942
943
      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
944
945
946
947
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
952
953
954
955
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
959
960
961
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
963
964
965
966
967
968
	}
    }

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

  if ( CDI_Debug )
969
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
970
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
971
}
972
973
974
975
976
977

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
978
  int status;
979
980
981
982
983
  long lpar;

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

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

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

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

992
993
  return (param);
}
994
995

static
996
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
997
998
{
  compvar2_t compVar;
999
1000
1001
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
1002
1003
1004
1005
1006

  compVar.param  = param;
  compVar.level1 = level1;
  compVar.level2 = level2;
  compVar.ltype  = leveltype;
1007
1008
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
1009
1010
1011
1012
1013
1014
1015
1016
1017

  return (compVar);
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record)
{
  int rstatus;
  compvar2_t compVar0;
1018
  size_t maxlen = sizeof(compVar.name);
1019
1020
1021
1022
1023

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