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

#include <stdio.h>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
#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
  int tsteptype;
32
  char name[32];
33
} compvar2_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
34
35
36


#if  defined  (HAVE_LIBGRIB_API)
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
63
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
64
static
65
int gribapiGetGridType(grib_handle *gh)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
{
67
  int gridtype = GRID_GENERIC;
68
  int gribgridtype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
  long lpar;

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

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

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

  return (gridtype);
}

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

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

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

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

  return (isRotated);
}

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

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

  return (zaxistype);
}

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

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

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
179
  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;
180
  long unitsOfTime = -1;
181

182
  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
	  *level1 = (int) dlevel;
	  *level2 = 0;
	}
      else
	{
	  GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
716
	  *level1 = (int)lpar;
717
	  GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
718
	  *level2 = (int)lpar;
719
	}
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
  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);

864
865
866
867
868
869
870
  (*record).size      = recsize;
  (*record).position  = position;
  (*record).param     = param;
  (*record).ilevel    = level1;
  (*record).ilevel2   = level2;
  (*record).ltype     = leveltype;
  (*record).tsteptype = tsteptype;
871
  memcpy((*record).varname, varname, len);
872
873

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

  gridID = varDefGrid(vlistID, grid, 0);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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