stream_gribapi.c 88.4 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
  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)
{
178
  int status;
179
  int timeunits = -1;
180
  long unitsOfTime = -1;
181

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

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

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

  return (timeunits);
}

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

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

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

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

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

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

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

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

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

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

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

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

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

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

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

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

      range = endStep - startStep;

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

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

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

	julday_add_seconds(addsec, &julday, &secofday);

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

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

  return (tstepRange);
375
376
}

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

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

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

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

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

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

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

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

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

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

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

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

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

487
	grid->size   = numberOfPoints;
488

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

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

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

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

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

521
			if ( fabs(grid->xinc-xinc) > 0.0 )
522
			  {
523
524
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
525
526
527
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
528
	      }
529
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
	  }
531
532
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
	  {
534
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
535
	      {
536
537
538
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
539
	      }
540
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
541
542
543
544
545
	  }
	break;
      }
    case GRID_LCC:
      {
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
	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);
562
563
564
565
566
567
568
569
570
571
572
573
574
	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
575

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

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

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

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

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

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

639
	grid->size  = numberOfPoints;
640

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

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

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

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

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

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

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

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

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

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

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

728
729
730
731
732
  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;
733
734
735
736
737
  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;
738
739
740
741

  return (scaleFactor);
}

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

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

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

763
      *leveltype = (int) lpar;
764

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

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

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

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

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

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

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

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

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

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

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

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

  gridID = varDefGrid(vlistID, grid, 0);

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

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

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

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

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

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

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

      {
        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
941
        else if ( strncmp(stdname, "unknown", 7) == 0 ) stdname[0] = 0;
942
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
943
    }
944
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
945

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

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

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

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

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

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

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

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

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1028
1029
1030
1031
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035
1036
1037