stream_gribapi.c 90 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 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
#  include "gribapi.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
24
25
26
#  include "grib_api.h"
#endif

extern int cdiInventoryMode;

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


#if  defined  (HAVE_LIBGRIB_API)
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
static
int my_grib_set_double(grib_handle* h, const char* key, double val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_double(\tgrib_handle* h, \"%s\", %f)\n", key, val);

  return grib_set_double(h, key, val);
}

static
int my_grib_set_long(grib_handle* h, const char* key, long val)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_long(  \tgrib_handle* h, \"%s\", %ld)\n", key, val);

  return grib_set_long(h, key, val);
}

static
int my_grib_set_string(grib_handle* h, const char* key, const char* val, size_t* length)
{
  if ( cdiGribApiDebug )
    fprintf(stderr, "grib_set_string(\tgrib_handle* h, \"%s\", \"%s\")\n", key, val);

  return grib_set_string(h, key, val, length);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
static
64
int gribapiGetGridType(grib_handle *gh)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
{
66
  int gridtype = GRID_GENERIC;
67
  int gribgridtype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
69
  long lpar;

70
    {
71
72
73
74
      int status;
      status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);

      if ( status ==  0 ) gribgridtype = (int) lpar;
75
76
77

      switch (gribgridtype)
	{
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
	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_UNSTRUCTURED:  { gridtype = GRID_UNSTRUCTURED; break; }
93
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
95
96
97
98
    }

  return (gridtype);
}

99
static
100
int gribapiGetIsRotated(grib_handle *gh)
101
102
{
  int isRotated = 0;
103
  int gribgridtype = -1;
104
  long lpar;
105
  int status;
106

107
  status = grib_get_long(gh, "gridDefinitionTemplateNumber", &lpar);
108

109
110
111
  if ( status ==  0 ) gribgridtype = (int) lpar;

  if ( gribgridtype == GRIB2_GTYPE_LATLON_ROT ) isRotated = 1;
112
113
114
115

  return (isRotated);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
static
117
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
121
  int zaxistype = ZAXIS_GENERIC;

  if ( editionNumber <= 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
    {
123
      zaxistype = grib1ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
126
    }
  else
    {
127
      zaxistype = grib2ltypeToZaxisType(grib_ltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
129
130
131
132
    }

  return (zaxistype);
}

133
static
134
int getTimeunits(long unitsOfTime)
135
136
{
  int timeunits = -1;
137
138
139
140
141
142
143
144
145
146
147
148
149

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

150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  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;
179
  long unitsOfTime = -1;
180

181
  grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
182

183
  GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
184

185
  timeunits = getTimeunits(unitsOfTime);
186
187
188
189

  return (timeunits);
}

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

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

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

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

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

220
  if ( editionNumber > 1 )
221
222
223
224
225
226
227
    {
      long sigofrtime;

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

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
228
229
230
231
232
233
234

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

270
271
272
273
274
275
276
277
278
279
280
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;
}

281
282
283
static
void gribapiSetDataDateTime(grib_handle *gh, int datadate, int datatime)
{
284
285
  GRIB_CHECK(my_grib_set_long(gh, "dataDate", datadate), 0);
  GRIB_CHECK(my_grib_set_long(gh, "dataTime", datatime/100), 0);
286
287
}

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

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

302
  if ( editionNumber > 1 )
303
304
305
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
306
307
308
309
  else
    {
      GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &sigofrtime), 0);
    }
310
311
312

  if ( sigofrtime == 3 )
    {
313
      gribapiGetDataDateTime(gh, vdate, vtime);
314
315
316
    }
  else
    {
317
318
      gribapiGetDataDateTime(gh, &rdate, &rtime);

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

      range = endStep - startStep;

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

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

	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
361
	    break;
362
363
364
365
366
	  }

	julday_add_seconds(addsec, &julday, &secofday);

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

368
369
370
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
371
    }
372
373

  return (tstepRange);
374
375
}

376
static
377
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378
{
379
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
381
  size_t datasize;
382
383
  long numberOfPoints;
  long lpar;
384
385
386

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

387
  gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
388
389
390
391
392
393
394
395
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
396
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397

Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
399
400
  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
401
402
403
404
405
406
407
408
409
410
411
412
  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;

413
414
415
416
417
418
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

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

422
423
424
425
426
427
428
429
430
431
432
	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
433
	if ( gridtype == GRID_LONLAT )
434
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435

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

438
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
	  {
440
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
	      {
442
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443

444
445
446
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
447
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
		      {
449
			double xinc = 360. / grid->xsize;
450

451
			if ( fabs(grid->xinc-xinc) > 0.0 )
452
			  {
453
454
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
455
			  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
457
458
		      }
		  }
	      }
459
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
	  }
461
	grid->ydef = 0;
462
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463
	  {
464
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
465
	      {
466
467
468
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
	      }
470
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
472
473
474
475
	  }
	break;
      }
    case GRID_GAUSSIAN_REDUCED:
      {
476
477
478
479
	int nlat, i;
	size_t dummy;
	long *pl;

480
481
482
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

483
484
485
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

486
	grid->size   = numberOfPoints;
487

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

495
496
497
498
499
500
501
502
503
	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);
504

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

507
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508
	  {
509
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
510
	      {
511
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
512
513
514
515

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

520
			if ( fabs(grid->xinc-xinc) > 0.0 )
521
			  {
522
523
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
524
525
526
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
	      }
528
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
	  }
530
531
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
	  {
533
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
	      {
535
536
537
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
	      }
539
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
541
542
543
544
	  }
	break;
      }
    case GRID_LCC:
      {
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
	int nlon, nlat;

	GRIB_CHECK(grib_get_long(gh, "Nx", &lpar), 0);
	nlon = lpar;
	GRIB_CHECK(grib_get_long(gh, "Ny", &lpar), 0);
	nlat = lpar;

	if ( numberOfPoints != nlon*nlat )
	  Error("numberOfPoints (%d) and gridSize (%d) differ!", (int)numberOfPoints, nlon*nlat);

	grid->size  = numberOfPoints;
	grid->xsize = nlon;
	grid->ysize = nlat;

	GRIB_CHECK(grib_get_double(gh, "DxInMetres", &grid->lcc_xinc), 0);
	GRIB_CHECK(grib_get_double(gh, "DyInMetres", &grid->lcc_yinc), 0);
561
562
563
564
565
566
567
568
569
570
571
572
573
	GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->lcc_originLon), 0);
	GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid->lcc_originLat), 0);
	GRIB_CHECK(grib_get_double(gh, "LoVInDegrees", &grid->lcc_lonParY), 0);
	GRIB_CHECK(grib_get_double(gh, "Latin1InDegrees", &grid->lcc_lat1), 0);
	GRIB_CHECK(grib_get_double(gh, "Latin2InDegrees", &grid->lcc_lat2), 0);

        if ( editionNumber <= 1 )
          {
            GRIB_CHECK(grib_get_long(gh, "projectionCenterFlag", &lpar), 0);
            grid->lcc_projflag  = (int) lpar;
            GRIB_CHECK(grib_get_long(gh, "scanningMode", &lpar), 0);
            grid->lcc_scanflag  = (int) lpar;
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574

575
	grid->xdef   = 0;
576
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
579
580
581

	break;
      }
    case GRID_SPECTRAL:
      {
582
583
	size_t len = 256;
	char typeOfPacking[256];
584
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
585
586
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
587

588
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
590
	grid->trunc = lpar;
591

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

	break;
      }
604
    case GRID_UNSTRUCTURED:
605
      {
606
        char uuid[17];
607
    	char reference_link[8192];
608
609
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
610

611
    	grid->size  = numberOfPoints;
612
613
614
615
        if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
          {
            grid->number   = lpar;
            if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
616
            /*
617
618
619
620
621
            if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
              {
                if ( strncmp(reference_link, "file://", 7) == 0 )
                  grid->reference = strdupx(reference_link);
              }
622
            */
623
624
625
            len = (size_t) 16;
            if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
              {
626
                memcpy(grid->uuid, uuid, 16);
627
628
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
629
630
631
632
	break;
      }
    case GRID_GENERIC:
      {
633
634
635
636
637
	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;

638
	grid->size  = numberOfPoints;
639

640
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
641
	  {
642
643
	    grid->xsize = nlon;
	    grid->ysize = nlat;
644
645
646
	  }
	else
	  {
647
648
	    grid->xsize = 0;
	    grid->ysize = 0;
649
650
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
651
652
653
654
	break;
      }
    default:
      {
655
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
656
657
658
659
	break;
      }
    }

660
  grid->isRotated = FALSE;
661
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
    {
663
664
665
666
      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);
667
      /* change from south to north pole */
668
669
670
671
672
673
674
675
676
677
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

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

static
678
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
679
{
680
681
682
  int status;
  long lpar;
  double dlevel;
683

684
685
686
687
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
688

689
690
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
691
    {
692
      *leveltype = (int) lpar;
693

694
      switch (*leveltype)
695
696
697
698
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
699
	  { *lbounds = 1; break; }
700
	}
701

702
703
704
705
706
      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;
707
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
708

709
710
711
712
713
714
715
716
717
718
	  *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;
	}
719
    }
720
721
}

722
723
724
static
double grib2ScaleFactor(long factor)
{
725
  double scaleFactor = 0;
726

727
728
729
730
731
  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;
732
733
734
735
736
  else if ( factor == 5 ) scaleFactor =    0.00001;
  else if ( factor == 6 ) scaleFactor =    0.000001;
  else if ( factor == 7 ) scaleFactor =    0.0000001;
  else if ( factor == 8 ) scaleFactor =    0.00000001;
  else if ( factor == 9 ) scaleFactor =    0.000000001;
737
738
739
740

  return (scaleFactor);
}

741
static
742
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf, int *level_unit)
743
{
744
745
746
747
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
748

749
750
751
752
753
754
  *leveltype  = 0;
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
755
756
757

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
758
    {
759
      long llevel;
760
761
      double dlevel1 = 0, dlevel2 = 0;

762
      *leveltype = (int) lpar;
763

764
765
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
766

767
      if ( *leveltype != 255 && leveltype2 != 255 && leveltype2 > 0 ) *lbounds = 1;
768
      if ( *leveltype == GRIB2_LTYPE_REFERENCE && leveltype2 == 1 ) *lbounds = 0;
769

770
771
772
773
774
775
776
777
778
779
      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;
        }
780
781
782
783
784
      else if ( *leveltype == GRIB2_LTYPE_SIGMA )
        {
          *level_sf = 1000;
          *level_unit = 0;
        }
785

786
787
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
788
789
790
791
792
793
794
      if ( llevel != GRIB_MISSING_LONG )
        {
          if ( factor != GRIB_MISSING_LONG )
            dlevel1 = llevel*grib2ScaleFactor(factor);
          else
            dlevel1 = llevel;
        }
795
796
797
798

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

      if ( *lbounds == 1 )
799
	{
800
801
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
802
803
804
805
806
807
808
          if ( llevel != GRIB_MISSING_LONG )
            {
              if ( factor != GRIB_MISSING_LONG )
                dlevel2 = llevel*grib2ScaleFactor(factor);
              else
                dlevel2 = llevel;
            }
809
810
811
812
813
814

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
815
    }
816
817
}

818
819
820
821
822
823
824
825
826
827
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;
}

828
#if  defined  (HAVE_LIBGRIB_API)
829
static
830
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
831
		      long recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
832
                      int leveltype, int lbounds, int level1, int level2, int level_sf, int level_unit)
833
834
835
836
837
838
839
840
841
842
843
844
845
{
  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;
846
  char stdname[CDI_MAX_NAME], longname[CDI_MAX_NAME], units[CDI_MAX_NAME];
847
  size_t vlen;
848
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
849

Uwe Schulzweida's avatar
Uwe Schulzweida committed
850
  vlistID = streamptr->vlistID;
851
  tsID    = streamptr->curTsID;
852
  recID   = recordNewEntry(streamptr, tsID);
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
  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;
869
  memcpy((*record).varname, varname, len);
870
871

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872
873
874

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
875
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
876

877
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
878
    {
879
880
881
882
883
884
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885

886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
        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];
902
        long ltmp;
903
        long nhlev, nvgrid;
904
905

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
906
        if ( lpar != 6 )
907
908
909
          {
            fprintf(stderr, "Warning ...\n");
          }
910
        GRIB_CHECK(grib_get_long(gh, "nlev", &ltmp), 0);
911
        nhlev = ltmp;
912
913
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &ltmp), 0);
        nvgrid = ltmp;
914
915
        len = (size_t) 16;
        uuid[16] = 0;
916
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
917
        varDefZAxisReference((int) nhlev, (int) nvgrid, uuid);
918
919
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
920
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921

922
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
923
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
924

925
  stdname[0] = 0;
926
927
928
  longname[0] = 0;
  units[0] = 0;

929
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
930
    {
931
      vlen = CDI_MAX_NAME;
932
      gribapiGetString(gh, "name", longname, vlen);
933
      vlen = CDI_MAX_NAME;
934
      gribapiGetString(gh, "units", units, vlen);
935
936
937
938
939

      {
        vlen = CDI_MAX_NAME;
        status = grib_get_string(gh, "cfName", stdname, &vlen);
        if ( status != 0 || vlen <= 1 ) stdname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
940
        else if ( strncmp(stdname, "unknown", 7) == 0 ) stdname[0] = 0;
941
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
    }
943
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
944

945
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
946
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
947
	       varname, stdname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
948
949
950
951

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

952
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953

954
955
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
956
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
957
  */
958
959
960
961
962
963
964
  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);
    }

965
966
967
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

968
969
970
971
972
  long typeOfGeneratingProcess = 0;
  status = grib_get_long(gh, "typeOfGeneratingProcess", &typeOfGeneratingProcess);
  if ( status == 0 )
    varDefTypeOfGeneratingProcess(varID, (int) typeOfGeneratingProcess);

973
974
975
976
977
  long productDefinitionTemplate = 0;
  status = grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplate);
  if ( status == 0 )
    varDefProductDefinitionTemplate(varID, (int) productDefinitionTemplate);

978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
  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 )
995
            varDefOptGribDbl(varID, dval, cdiAdditionalGRIBKeys[i]);
996
997
998
	}
      /* note: if the key is not defined, we do not throw an error! */
    }
999

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1000
1001
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
1003
1004
1005
1006
      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
1007
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1008
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
1010
1011
1012
1013
1014
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1015
      long processID;
1016
1017
1018
1019
1020
1021
1022
1023
      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
1024
1025
1026
1027
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032
1033
1034
1035
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1036

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1039