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

#include <stdio.h>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "dmemory.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
9
10
11
12
13
#include "cdi.h"
#include "stream_int.h"
#include "file.h"
#include "varscan.h"
#include "datetime.h"
#include "vlist.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
14
#include "stream_grb.h"
15
#include "calendar.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
17
18
19


#if  defined  (HAVE_LIBGRIB_API)
#  include "cgribex.h"      /* gribGetSize, gribRead, gribGetZip */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
#  include "gribapi.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
24
25
26
27
#  include "grib_api.h"
#endif


extern int cdiInventoryMode;

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


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

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

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

      switch (gribgridtype)
	{
52
	case  GRIB2_GTYPE_LATLON:     { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
	                                if ( lpar == (long) GRIB_MISSING_LONG ) break;
54
                                      }
55
56
	case  GRIB2_GTYPE_LATLON_ROT: { gridtype = GRID_LONLAT;    break; }
	case  GRIB2_GTYPE_LCC:        { gridtype = GRID_LCC;       break; }
57
	case  GRIB2_GTYPE_GAUSSIAN:   { GRIB_CHECK(grib_get_long(gh, "Ni", &lpar), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
	                                if ( lpar == (long) GRIB_MISSING_LONG )
59
60
61
62
63
					  gridtype = GRID_GAUSSIAN_REDUCED;
					else
					  gridtype = GRID_GAUSSIAN;
					break;
                                      }
64
65
66
	case  GRIB2_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL;  break; }
	case  GRIB2_GTYPE_GME:        { gridtype = GRID_GME;       break; }
	case  GRIB2_GTYPE_NUMBER:     { gridtype = GRID_REFERENCE; break; }
67
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
69
70
71
72
    }

  return (gridtype);
}

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

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

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

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

  return (isRotated);
}

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

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

  return (zaxistype);
}

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

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

124
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 gribapiGetEndStep(grib_handle *gh, int startStep, int timeunits)
{
  int endStep = startStep;
  int timeunits2;
154
  int status;
155
156
157
  long unitsOfTime;
  long lpar;

158
  status = grib_get_long(gh, "stepUnits", &unitsOfTime);
159
160
161

  timeunits2 = getTimeunits(unitsOfTime);

162
  status = grib_get_long(gh, "endStep", &lpar);
163

164
165
  if ( status == 0 )
    endStep = (int) ((lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
166
167
168
169
170
171
172
173

  return (endStep);
}

static
int gribapiGetTimeUnits(grib_handle *gh)
{
  int timeunits = -1;
174
175
  long unitsOfTime = -1;
  int status;
176
177
178
179
  // size_t len = 8;
  //char stepunits[8];
  //static int lprint = TRUE;

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

  timeunits = getTimeunits(unitsOfTime);

184
  /*
185
186
  GRIB_CHECK(grib_get_string(gh, "stepUnits", stepunits, &len), 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
189
190
191
192
193
194
195
196
197
  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;
198
  else if ( lprint )
199
    {
200
      Message("Step units >%s< unsupported!", stepunits);
201
      lprint = FALSE;
202
    }
203
  */
204
205
206
207
208
209
210

  return (timeunits);
}

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

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

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

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

      if ( sigofrtime == 3 ) isFC = FALSE;
    }
224
225
226
227
228
229
230

  return (isFC);
}

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

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

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

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

  return (tsteptype);
}

266
static
267
int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
268
{
269
  int rdate, rtime;
270
  int timeUnits, startStep = 0, endStep;
271
272
  int tstepRange = 0;
  int range;
273
  int status;
274
  long lpar;
275
276
277
278
279
  long sigofrtime = 3;
  long editionNumber;

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

280
  if ( editionNumber > 1 )
281
282
283
    {
      GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
    }
284
285
286
287
288
289
290
291
292
293

  if ( sigofrtime == 3 )
    {
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      *vdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
      *vtime = (int) lpar*100;
    }
  else
    {
294
295
296
      GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
      rdate = (int) lpar;
      GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
297
      rtime = (int) lpar*100;
298
299
      status = grib_get_long(gh, "forecastTime", &lpar);
      if ( status == 0 ) startStep = (int) lpar;
300
      timeUnits = gribapiGetTimeUnits(gh);
301
      endStep = gribapiGetEndStep(gh, startStep, timeUnits);
302
303
304
305
306
307
308
309
310

      range = endStep - startStep;

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

311
312
313
314
      {
	static int lprint = TRUE;
	extern int grib_calendar;
	int ryear, rmonth, rday, rhour, rminute, rsecond;
315
	int time_period = endStep;
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
	int julday, secofday, addsec;

	cdiDecodeDate(rdate, &ryear, &rmonth, &rday);
	cdiDecodeTime(rtime, &rhour, &rminute, &rsecond);

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

	addsec = 0;
	switch ( timeUnits )
	  {
	  case TUNIT_SECOND:  addsec =         time_period; break;
	  case TUNIT_MINUTE:  addsec =    60 * time_period; break;
	  case TUNIT_HOUR:    addsec =  3600 * time_period; break;
	  case TUNIT_3HOURS:  addsec = 10800 * time_period; break;
	  case TUNIT_6HOURS:  addsec = 21600 * time_period; break;
	  case TUNIT_12HOURS: addsec = 43200 * time_period; break;
	  case TUNIT_DAY:     addsec = 86400 * time_period; break;
	  default:
	    if ( lprint )
	      {
	        Warning("Time unit %d unsupported", timeUnits);
		lprint = FALSE;
	      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
	    break;
340
341
342
343
344
345
346
347
348
349
350
	  }

	julday_add_seconds(addsec, &julday, &secofday);

	decode_caldaysec(grib_calendar, julday, secofday, &ryear, &rmonth, &rday, &rhour, &rminute, &rsecond);
	/*
	  printf("new %d/%d/%d %d:%d\n", ryear, rmonth, rday, rhour, rminute);
	*/
	*vdate = cdiEncodeDate(ryear, rmonth, rday);
	*vtime = cdiEncodeTime(rhour, rminute, rsecond);
      }
351
    }
352
353

  return (tstepRange);
354
355
}

356
static
357
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
358
{
359
  long editionNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360
  int gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
  size_t datasize;
362
363
  long numberOfPoints;
  long lpar;
364
365
366

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

367
  gridtype = gribapiGetGridType(gh);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368
369
370
371
372
373
374
375
  /*
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }
  */
376
  memset(grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
377

Uwe Schulzweida's avatar
Uwe Schulzweida committed
378
379
380
  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
381
382
383
384
385
386
387
388
389
390
391
392
  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;

393
394
395
396
397
398
	if ( gridtype == GRID_GAUSSIAN )
          {
            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
            grid->np = lpar;
          }

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

402
403
404
405
406
407
408
409
410
411
412
	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
413
	if ( gridtype == GRID_LONLAT )
414
	  GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415

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

418
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
	  {
420
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
	      {
422
		if ( (grid->xfirst >= grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423

424
425
426
		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
427
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
		      {
429
			double xinc = 360. / grid->xsize;
430

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

460
461
462
        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
        grid->np = lpar;

463
464
465
	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
	nlat = lpar;

466
	grid->size   = numberOfPoints;
467

468
        grid->rowlon = (int *) malloc(nlat*sizeof(int));
469
470
471
        pl          = (long *) malloc(nlat*sizeof(long));
	dummy       = nlat;
	GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
472
	for ( i = 0; i < nlat; ++i ) grid->rowlon[i] = pl[i];
473
474
	free(pl);

475
476
477
478
479
480
481
482
483
	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);
484

485
	if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
486

487
	/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488
	  {
489
	    if ( grid->xsize > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490
	      {
491
		if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
492
493
494
495

		if ( editionNumber <= 1 )
		  {
		    /* correct xinc if necessary */
496
		    if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
497
		      {
498
			double xinc = 360. / grid->xsize;
499

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

530
531
532
	grid->size  = ISEC4_NumValues;
	grid->xsize = ISEC2_NumLon;
	grid->ysize = ISEC2_NumLat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533

534
535
536
537
538
539
540
541
542
	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
543

544
	grid->xdef   = 0;
545
	grid->ydef   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547
548

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
551
    case GRID_SPECTRAL:
      {
552
553
	size_t len = 256;
	char typeOfPacking[256];
554
	GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
555
556
	grid->lcomplex = 0;
	if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
557

558
	grid->size  = datasize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559
	GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
560
	grid->trunc = lpar;
561

Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
563
564
565
	break;
      }
    case GRID_GME:
      {
566
567
568
569
570
	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;
571
572
573

	break;
      }
574
    case GRID_REFERENCE:
575
      {
576
        char uuid[17];
577
578
579
580
	char reference_link[8192];
	size_t len = sizeof(reference_link);
	reference_link[0] = 0;

581
	grid->size  = numberOfPoints;
582
583
	if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
	  {
584
585
	    grid->number   = lpar;
	    if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
586
587
	    if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
	      {
588
		if ( strncmp(reference_link, "file://", 7) == 0 )
589
		  grid->reference = strdupx(reference_link);
590
	      }
591
592
593
594
595
            len = (size_t) 16;
            if ( grib_get_string(gh, "uuidOfHGrid", uuid, &len) == 0)
              {
                strncpy(grid->uuid, uuid, 16);
              }
596
	  }
597

Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
599
600
601
	break;
      }
    case GRID_GENERIC:
      {
602
603
604
605
606
	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;

607
	grid->size  = numberOfPoints;
608

609
	if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
610
	  {
611
612
	    grid->xsize = nlon;
	    grid->ysize = nlat;
613
614
615
	  }
	else
	  {
616
617
	    grid->xsize = 0;
	    grid->ysize = 0;
618
619
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
620
621
622
623
	break;
      }
    default:
      {
624
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
625
626
627
628
	break;
      }
    }

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

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

static
647
void grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
648
{
649
650
651
  int status;
  long lpar;
  double dlevel;
652

653
654
655
656
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;
657

658
659
  status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
  if ( status == 0 )
660
    {
661
      *leveltype = (int) lpar;
662

663
      switch (*leveltype)
664
665
666
667
	{
	case GRIB1_LTYPE_SIGMA_LAYER:
	case GRIB1_LTYPE_HYBRID_LAYER:
	case GRIB1_LTYPE_LANDDEPTH_LAYER:
668
	  { *lbounds = 1; break; }
669
	}
670

671
672
673
674
675
676
      if ( *lbounds == 0 )
	{
	  GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	  if ( *leveltype == 100 ) dlevel *= 100;
	  if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
	  if ( *leveltype == 99 ) *leveltype = 100;
677

678
679
680
681
682
683
684
685
686
687
	  *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;
	}
688
    }
689
690
691
}

static
692
void grib2GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2)
693
{
694
695
696
697
  int status;
  int leveltype2 = -1;
  long lpar;
  long factor;
698
699
  double dlevel;

700
701
702
703
704
705
706
  *leveltype = 0;
  *lbounds = 0;
  *level1  = 0;
  *level2  = 0;

  status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
  if ( status == 0 )
707
    {
708
      *leveltype = (int) lpar;
709

710
711
      status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar);
      if ( status == 0 ) leveltype2 = lpar;
712

713
      if ( *leveltype == leveltype2 && *leveltype != 255 ) *lbounds = 1;
714

715
      if ( *lbounds == 0 )
716
	{
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
	      if ( *leveltype == GRIB2_LTYPE_ISOBARIC ) dlevel *= 100;
	      if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
	      if ( *leveltype == 99 ) *leveltype = 100;
	    }

733
	  *level1 = (int) dlevel;
734
	  *level2 = 0;
735
736
737
	}
      else
	{
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
	  if ( *leveltype == GRIB2_LTYPE_LANDDEPTH )
	    {
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfFirstFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	      *level1 = (int) dlevel;
	      GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0);
	      GRIB_CHECK(grib_get_double(gh, "scaledValueOfSecondFixedSurface", &dlevel), 0);
	      if      ( factor == 0 ) dlevel *= 100;   //  m to cm
	      else if ( factor == 1 ) dlevel *=  10;   // dm to cm
	      else if ( factor == 3 ) dlevel *=   0.1; // mm to cm
	      *level2 = (int) dlevel;
	    }
	  else
	    {
	      GRIB_CHECK(grib_get_long(gh, "topLevel", &lpar), 0);
	      *level1 = lpar;
	      GRIB_CHECK(grib_get_long(gh, "bottomLevel", &lpar), 0);
	      *level2 = lpar;
	    }
760
761
	}
    }
762
763
}

764
765
766
767
768
769
770
771
772
773
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;
}

774
static
775
void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
776
		      long recsize, off_t position, int datatype, int comptype)
777
778
779
780
781
782
{
  long editionNumber;
  int zaxistype;
  int gridID = CDI_UNDEFID, varID;
  int levelID = 0;
  int tsID, recID;
783
  int level1 = 0, level2 = 0;
784
785
786
787
788
789
790
791
792
793
  int numavg;
  int tsteptype;
  int lbounds = 0;
  record_t *record;
  grid_t grid;
  int vlistID;
  int leveltype;
  long lpar;
  int status;
  char name[256], longname[256], units[256];
794
  size_t vlen;
795
  long ens_index = 0, ens_count = 0, ens_forecast_type = 0;
796

Uwe Schulzweida's avatar
Uwe Schulzweida committed
797
  vlistID = streamptr->vlistID;
798
  tsID    = streamptr->curTsID;
799
  recID   = recordNewEntry(streamptr, tsID);
800
801
802
803
804
805
806
807
808
  record  = &streamptr->tsteps[tsID].records[recID];

  tsteptype = gribapiGetTsteptype(gh);
  // numavg  = ISEC1_AvgNum;
  numavg  = 0;

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

  if ( editionNumber <= 1 )
809
    grib1GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
810
  else
811
    grib2GetLevel(gh, &leveltype, &lbounds, &level1, &level2);
812

813
814
815
816
817
818
819
820
821
822
  // 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;

  gribapiGetGrid(gh, &grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
823
824
825

  gridID = varDefGrid(vlistID, grid, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
  zaxistype = gribapiGetZaxisType(editionNumber, leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827

828
  switch (zaxistype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
    {
830
831
832
833
834
835
    case ZAXIS_HYBRID:
    case ZAXIS_HYBRID_HALF:
      {
        int vctsize;
        size_t dummy;
        double *vctptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836

837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
        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];
        long nlev, nvgrid;

        GRIB_CHECK(grib_get_long(gh, "NV", &lpar), 0);
856
        if ( lpar != 3 )
857
858
859
860
861
862
863
864
865
866
867
868
          {
            fprintf(stderr, "Warning ...\n");
          }
        GRIB_CHECK(grib_get_long(gh, "nlev", &nlev), 0);
        GRIB_CHECK(grib_get_long(gh, "numberOfVGridUsed", &nvgrid), 0);

        len = (size_t) 16;
        uuid[16] = 0;
        GRIB_CHECK(grib_get_string(gh, "uuidOfVGrid", uuid, &len), 0);
        varDefZAxisReference((int) nlev, (int) nvgrid, uuid);
        break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
869
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
870

871
  // if ( datatype > 32 ) datatype = DATATYPE_PACK32;
872
  if ( datatype <  0 ) datatype = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
873

874
875
876
877
  name[0] = 0;
  longname[0] = 0;
  units[0] = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
878
  vlen = 256;
879
  gribapiGetString(gh, "shortName", name, vlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
880
881
882
883

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
890
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2,
891
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
892
	       name, longname, units);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
893
894
895
896

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

897
  varDefCompType(varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898

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


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

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
948
949
950
951
      if ( pdis == 255 )
	{
	  int tableID;
	  int tabnum = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
952

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
956
957
	  if ( tableID == CDI_UNDEFID )
	    tableID = tableDef(varInqModel(varID), tabnum, NULL);
	  varDefTable(varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
958
959
960
961
962
963
964
	}
    }

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

  if ( CDI_Debug )
965
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
967
}
968
969
970
971
972
973

static
int gribapiGetParam(grib_handle *gh)
{
  int pdis = 0, pcat = 0, pnum = 0;
  int param = 0;
974
  int status;
975
976
977
978
979
  long lpar;

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

980
981
  status = grib_get_long(gh, "parameterCategory", &lpar);
  if ( status == 0 ) pcat = (int) lpar;
982

983
984
  status = grib_get_long(gh, "parameterNumber", &lpar);
  if ( status == 0 ) pnum = (int) lpar;
985
986

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

988
989
  return (param);
}
990
991

static
992
compvar2_t gribapiVarSet(int param, int level1, int level2, int leveltype, char *name)
993
994
995
996
997
998
999
{
  compvar2_t compVar;

  compVar.param  = param;
  compVar.level1 = level1;
  compVar.level2 = level2;
  compVar.ltype  = leveltype;
1000
  //memset(compVar.name, 0, sizeof(compVar.name));
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014

  return (compVar);
}

static
int gribapiVarCompare(compvar2_t compVar, record_t record)
{
  int rstatus;
  compvar2_t compVar0;

  compVar0.param  = record.param;
  compVar0.level1 = record.ilevel;
  compVar0.level2 = record.ilevel2;
  compVar0.ltype  = record.ltype;
1015
  //memset(compVar0.name, 0, sizeof(compVar0.name));
1016
1017
1018
1019
1020

  rstatus = memcmp(&compVar0, &compVar, sizeof(compvar2_t));

  return (rstatus);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1021
1022
#endif

1023
int gribapiScanTimestep1(stream_t * streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
1026
1027
1028
{
#if  defined  (HAVE_LIBGRIB_API)
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1029
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1030
1031
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032
  int rtabnum = 0;
1033
1034
  int rcode = 0, level1 = 0, level2 = 0;
  int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1036
1037
1038
1039
1040
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
1041
  int nrecs_scanned;