stream_gribapi.c 85.6 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
#  include "grib_api.h"
#endif

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

extern int cdiInventoryMode;

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


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

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

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

      switch (gribgridtype)
	{
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
	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; }
68
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
71
72
73
    }

  return (gridtype);
}

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

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

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

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

  return (isRotated);
}

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

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

  return (zaxistype);
}

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

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

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  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;
154
155
  long unitsOfTime = -1;
  int status;
156
157
158
159
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

160
  status = grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
161
162
163

  timeunits = getTimeunits(unitsOfTime);

164
  /*
165
166
  GRIB_CHECK(grib_get_string(gh, "stepUnits", stepunits, &len), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
168
169
170
171
172
173
174
175
176
177
  len--;

  if      ( memcmp(stepunits, "s",   len) == 0 ) timeunits = TUNIT_SECOND;
  else if ( memcmp(stepunits, "m",   len) == 0 ) timeunits = TUNIT_MINUTE;
  else if ( memcmp(stepunits, "h",   len) == 0 ) timeunits = TUNIT_HOUR;
  else if ( memcmp(stepunits, "3h",  len) == 0 ) timeunits = TUNIT_3HOURS;
  else if ( memcmp(stepunits, "6h",  len) == 0 ) timeunits = TUNIT_6HOURS;
  else if ( memcmp(stepunits, "12h", len) == 0 ) timeunits = TUNIT_12HOURS;
  else if ( memcmp(stepunits, "D",   len) == 0 ) timeunits = TUNIT_DAY;
  else if ( memcmp(stepunits, "M",   len) == 0 ) timeunits = TUNIT_MONTH;
  else if ( memcmp(stepunits, "Y",   len) == 0 ) timeunits = TUNIT_YEAR;
178
  else if ( lprint )
179
    {
180
      Message("Step units >%s< unsupported!", stepunits);
181
      lprint = FALSE;
182
    }
183
  */
184
185
186
187

  return (timeunits);
}

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

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

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

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

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

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

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

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

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

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

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

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

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

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

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

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

      range = endStep - startStep;

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

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

	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
359
	    break;
360
361
362
363
364
	  }

	julday_add_seconds(addsec, &julday, &secofday);

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

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

  return (tstepRange);
372
373
}

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

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

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

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

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

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

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

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

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

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

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

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

481
482
483
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

484
	grid->size   = numberOfPoints;
485

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

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

503
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
504

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

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

518
			if ( fabs(grid->xinc-xinc) > 0.0 )
519
			  {
520
521
			    grid->xinc = xinc;
			    if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
522
523
524
			  }
		      }
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
	      }
526
	    grid->xdef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
	  }
528
529
	grid->ydef  = 0;
        /* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
	  {
531
	    if ( grid->ysize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
	      {
533
534
535
		if ( editionNumber <= 1 )
		  {
		  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536
	      }
537
	    grid->ydef = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
539
540
	  }
	break;
      }
541
      /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
543
544
    case GRID_LCC:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
545
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);

548
549
550
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551

552
553
554
555
556
557
558
559
560
	grid->lcc_xinc      = ISEC2_Lambert_dx;
	grid->lcc_yinc      = ISEC2_Lambert_dy;
	grid->lcc_originLon = ISEC2_FirstLon * 0.001;
	grid->lcc_originLat = ISEC2_FirstLat * 0.001;
	grid->lcc_lonParY   = ISEC2_Lambert_Lov * 0.001;
	grid->lcc_lat1      = ISEC2_Lambert_LatS1 * 0.001;
	grid->lcc_lat2      = ISEC2_Lambert_LatS2 * 0.001;
	grid->lcc_projflag  = ISEC2_Lambert_ProjFlag;
	grid->lcc_scanflag  = ISEC2_ScanFlag;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561

562
	grid->xdef   = 0;
563
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
565
566

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
568
569
    case GRID_SPECTRAL:
      {
570
571
	size_t len = 256;
	char typeOfPacking[256];
572
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
573
574
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
575

576
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
578
	grid->trunc = lpar;
579

Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
581
582
583
	break;
      }
    case GRID_GME:
      {
584
585
586
587
588
	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;
589
590
591

	break;
      }
592
    case GRID_UNSTRUCTURED:
593
      {
594
        char uuid[17];
595
    	char reference_link[8192];
596
597
        size_t len = sizeof(reference_link);
        reference_link[0] = 0;
598

599
    	grid->size  = numberOfPoints;
600
601
602
603
        if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
          {
            grid->number   = lpar;
            if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
604
            /*
605
606
607
608
609
            if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
              {
                if ( strncmp(reference_link, "file://", 7) == 0 )
                  grid->reference = strdupx(reference_link);
              }
610
            */
611
612
613
            len = (size_t) 16;
            if ( grib_get_bytes(gh, "uuidOfHGrid", (unsigned char *) uuid, &len) == 0)
              {
614
                memcpy(grid->uuid, uuid, 16);
615
616
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
618
619
620
	break;
      }
    case GRID_GENERIC:
      {
621
622
623
624
625
	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;

626
	grid->size  = numberOfPoints;
627

628
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
629
	  {
630
631
	    grid->xsize = nlon;
	    grid->ysize = nlat;
632
633
634
	  }
	else
	  {
635
636
	    grid->xsize = 0;
	    grid->ysize = 0;
637
638
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
639
640
641
642
	break;
      }
    default:
      {
643
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
645
646
647
	break;
      }
    }

648
  grid->isRotated = FALSE;
649
  if ( gribapiGetIsRotated(gh) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
    {
651
652
653
654
      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);
655
      /* change from south to north pole */
656
657
658
659
660
661
662
663
664
665
      grid->ypole = -grid->ypole;
      grid->xpole =  grid->xpole - 180;
    }

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

static
666
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
667
{
668
669
670
  int status;
  long lpar;
  double dlevel;
671

672
673
674
675
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
676

677
678
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
679
    {
680
      *leveltype = (int) lpar;
681

682
      switch (*leveltype)
683
684
685
686
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
687
	  { *lbounds = 1; break; }
688
	}
689

690
691
692
693
694
      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;
695
	  if ( *leveltype == GRIB1_LTYPE_99 ) *leveltype = 100;
696

697
698
699
700
701
702
703
704
705
706
	  *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;
	}
707
    }
708
709
}

710
711
712
static
double grib2ScaleFactor(long factor)
{
713
  double scaleFactor = 0;
714

715
716
717
718
719
  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;
720
721
722
723

  return (scaleFactor);
}

724
static
725
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2, int *level_sf, int *level_unit)
726
{
727
728
729
730
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
731

732
733
734
735
736
737
  *leveltype  = 0;
  *lbounds    = 0;
  *level1     = 0;
  *level2     = 0;
  *level_sf   = 0;
  *level_unit = 0;
738
739
740

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
741
    {
742
      long llevel;
743
744
      double dlevel1 = 0, dlevel2 = 0;

745
      *leveltype = (int) lpar;
746

747
748
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
749

750
      if ( *leveltype != 255 && leveltype2 != 255 && leveltype2 > 0 ) *lbounds = 1;
751

752
753
754
755
756
757
758
759
760
761
      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;
        }
762

763
764
765
766
767
768
769
770
      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
      GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);
      if ( llevel != GRIB_MISSING_LONG && factor != GRIB_MISSING_LONG )
        dlevel1 = llevel*grib2ScaleFactor(factor);

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

      if ( *lbounds == 1 )
771
	{
772
773
774
775
776
777
778
779
780
781
          GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
          GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0);
          if ( llevel != GRIB_MISSING_LONG && factor != GRIB_MISSING_LONG )
            dlevel2 = llevel*grib2ScaleFactor(factor);

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

      *level1 = (int) dlevel1;
      *level2 = (int) dlevel2;
782
    }
783
784
}

785
786
787
788
789
790
791
792
793
794
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;
}

795
#if  defined  (HAVE_LIBGRIB_API)
796
static
797
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
798
		      long recsize, off_t position, int datatype, int comptype, size_t len, const char *varname,
799
                      int leveltype, int lbounds, int level1, int level2, int level_sf, int level_unit)
800
801
802
803
804
805
806
807
808
809
810
811
812
{
  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;
813
  char longname[256], units[256];
814
  size_t vlen;
815
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
816

Uwe Schulzweida's avatar
Uwe Schulzweida committed
817
  vlistID = streamptr->vlistID;
818
  tsID    = streamptr->curTsID;
819
  recID   = recordNewEntry(streamptr, tsID);
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
  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;
836
  memcpy((*record).varname, varname, len);
837
838

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
840
841

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
842
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843

844
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
845
    {
846
847
848
849
850
851
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
852

853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
        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];
869
        double dtmp;
870
871
872
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
873
        if ( lpar != 6 )
874
875
876
          {
            fprintf(stderr, "Warning ...\n");
          }
877
878
879
880
        GRIB_CHECK(grib_get_double(gh, "nlev", &dtmp), 0);
        nlev = (int) NINT(dtmp);
        GRIB_CHECK(grib_get_double(gh, "numberOfVGridUsed", &dtmp), 0);
        nvgrid = NINT(dtmp);
881
882
        len = (size_t) 16;
        uuid[16] = 0;
883
        GRIB_CHECK(grib_get_bytes(gh, "uuidOfVGrid", (unsigned char *) uuid, &len), 0);
884
885
886
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888

889
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
890
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
891

892
893
894
  longname[0] = 0;
  units[0] = 0;

895
  if ( varname[0] != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
896
897
    {
      vlen = 256;
898
      gribapiGetString(gh, "name", longname, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
      vlen = 256;
900
      gribapiGetString(gh, "units", units, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
901
    }
902
  // fprintf(stderr, "param %d name %s %s %s\n", param, name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
903

904
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, level_sf, level_unit,
905
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
906
	       varname, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
907
908
909
910

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

911
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
912

913
914
  /*
    Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure.
915
    Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars"
916
  */
917
918
919
920
921
922
923
924
925
  {
    int status;
    status = grib_get_long(gh, "typeOfEnsembleForecast", &ens_forecast_type );
    if ( status == 0 )
      {
        GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", &ens_count ), 0);
        GRIB_CHECK(grib_get_long(gh, "perturbationNumber", &ens_index ), 0);
      }
  }
926
927
928
  if ( ens_index > 0 )
    varDefEnsembleInfo(varID, (int)ens_index, (int)ens_count, (int)ens_forecast_type);

929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
  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 )
            varDefOptGribInt(varID, dval, cdiAdditionalGRIBKeys[i]);
	}
      /* note: if the key is not defined, we do not throw an error! */
    }
950

Uwe Schulzweida's avatar
Uwe Schulzweida committed
951
952
  if ( varInqInst(varID) == CDI_UNDEFID )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953
954
955
956
957
      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
958
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
959
	instID = institutDef((int)center, (int)subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960
961
962
963
964
965
      varDefInst(varID, instID);
    }

  if ( varInqModel(varID) == CDI_UNDEFID )
    {
      int modelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
      long processID;
967
968
969
970
971
972
973
974
      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
975
976
977
978
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
983
984
985
986
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
987

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
990
991
992
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
994
995
996
997
998
999
	}
    }

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

  if ( CDI_Debug )
1000
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1001
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
}
1003
#endif
1004
1005
1006
1007
1008
1009

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
1010
  int status;
1011
1012
1013
1014
1015
  long lpar;

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

1016
1017
  status = grib_get_long(gh, "parameterCategory", &lpar);
  if ( status == 0 ) pcat = (int) lpar;
1018

1019
1020
  status = grib_get_long(gh, "parameterNumber", &lpar);
  if ( status == 0 ) pnum = (int) lpar;
1021
1022

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

1024
1025
  return (param);
}
1026
1027

static
1028
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
1029
1030
{
  compvar2_t compVar;
1031
1032
1033
  size_t maxlen = sizeof(compVar.name);
  size_t len = strlen(name);
  if ( len > maxlen ) len = maxlen;
1034
1035
1036
1037
1038

  compVar.param  = param;
  compVar.level1 = level1;
  compVar.level2 = level2;
  compVar.ltype  = leveltype;
1039
1040
  memset(compVar.name, 0, maxlen);
  memcpy(compVar.name, name, len);
1041
1042
1043
1044
1045
1046
1047
1048
1049

  return (compVar);
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record)
{
  int rstatus;
  compvar2_t compVar0;
1050
  size_t maxlen = sizeof(compVar.name);
1051
1052
1053
1054
1055

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