stream_cgribex.c 53.9 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <stdio.h>
#include <float.h>  /* FLT_EPSILON */

#include "dmemory.h"
#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
15
#include "stream_grb.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16

Uwe Schulzweida's avatar
Uwe Schulzweida committed
17
18
19
#if  defined  (HAVE_LIBCGRIBEX)
#  include "cgribex.h"
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
extern int cdiInventoryMode;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
23

typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
  int param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
  int level1;
  int level2;
  int ltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
} compvar_t; 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
30


Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
int cgribexGetGridType(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
34
{
35
  int gridtype = GRID_GENERIC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
38

  switch (ISEC2_GridType)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
    case  GRIB1_GTYPE_LATLON:
40
41
42
43
44
45
46
47
48
49
    case  GRIB1_GTYPE_LATLON_ROT: { gridtype = GRID_LONLAT;   break; }
    case  GRIB1_GTYPE_LCC:        { gridtype = GRID_LCC;      break; }
    case  GRIB1_GTYPE_GAUSSIAN:   { if ( ISEC2_Reduced )
	                              gridtype = GRID_GAUSSIAN_REDUCED;
                         	    else
				      gridtype = GRID_GAUSSIAN;
          	                    break;
                                  }
    case  GRIB1_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL; break; }
    case  GRIB1_GTYPE_GME:        { gridtype = GRID_GME;      break; }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52
53
    }

  return (gridtype);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55

Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
int cgribexGetIsRotated(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
61
{
  int isRotated = 0;

62
  if ( ISEC2_GridType == GRIB1_GTYPE_LATLON_ROT )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
65
66
67
68
    {
      isRotated = 1;
    }

  return (isRotated);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70

Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
int cgribexGetZaxisHasBounds(int grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
76
{
  int lbounds = 0;

77
  switch (grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
    case GRIB1_LTYPE_HYBRID_LAYER:
    case GRIB1_LTYPE_LANDDEPTH_LAYER:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
82
83
84
85
86
87
88
      {
	lbounds = 1;
	break;
      }
    }

  return (lbounds);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
int cgribexGetTimeUnit(int *isec1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
95
96
97
98
99
{
  int timeunit = -1;
  static int lprint = TRUE;

  switch ( ISEC1_TimeUnit )
    {
100
101
102
    case ISEC1_TABLE4_MINUTE:  timeunit = TUNIT_MINUTE;  break;
    case ISEC1_TABLE4_QUARTER: timeunit = TUNIT_QUARTER; break;
    case ISEC1_TABLE4_HOUR:    timeunit = TUNIT_HOUR;    break;
103
104
105
    case ISEC1_TABLE4_3HOURS:  timeunit = TUNIT_3HOURS;  break;
    case ISEC1_TABLE4_6HOURS:  timeunit = TUNIT_6HOURS;  break;
    case ISEC1_TABLE4_12HOURS: timeunit = TUNIT_12HOURS; break;
106
    case ISEC1_TABLE4_DAY:     timeunit = TUNIT_DAY;     break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
109
    default:
      if ( lprint )
	{
110
	  Message("Time unit %d unsupported", ISEC1_TimeUnit);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
112
113
114
115
116
	  lprint = FALSE;
	}
    }

  return (timeunit);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
#if  defined  (HAVE_LIBCGRIBEX)
120
121
122
123
124
125
126
127
128
129
static
int cgribexTimeIsFC(int *isec1)
{
  int isFC = TRUE;

  if ( ISEC1_TimeRange == 10 && ISEC1_TimePeriod1 == 0 && ISEC1_TimePeriod2 == 0 )
    isFC = FALSE;

  return (isFC);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
#endif
131

Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
#if  defined  (HAVE_LIBCGRIBEX)
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
static
int cgribexGetTsteptype(int timerange)
{
  int tsteptype = 0;
  static int lprint = TRUE;

  switch ( timerange )
    {
    case  0:  tsteptype = TSTEP_INSTANT;  break;
    case  1:  tsteptype = TSTEP_INSTANT2; break;
    case  2:  tsteptype = TSTEP_RANGE;    break;
    case  3:  tsteptype = TSTEP_AVG;      break;
    case  4:  tsteptype = TSTEP_ACCUM;    break;
    case  5:  tsteptype = TSTEP_DIFF;     break;
    case 10:  tsteptype = TSTEP_INSTANT3; break;
    default:
      if ( lprint )
	{
151
	  Message("Time range %d unsupported", timerange);
152
153
154
155
156
157
	  lprint = FALSE;
	}
    }

  return (tsteptype);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
#endif
159

Uwe Schulzweida's avatar
Uwe Schulzweida committed
160
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
static
162
163
void cgribexAddRecord(int streamID, int param, int *isec1, int *isec2, double *fsec2, double *fsec3,
		      int *isec4, long recsize, off_t position, int prec, int ztype, int lmv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
165
{
  int gridtype;
166
  int zaxistype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
  int gridID = CDI_UNDEFID, varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
169
170
171
  int levelID = 0;
  int tsID, recID;
  int level1, level2;
  int numavg;
172
  int tsteptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
  int lbounds = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
  record_t *record;
175
  grid_t grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
  int vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
179

  streamptr = stream_to_pointer(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181

  vlistID = streamInqVlist(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
  tsID    = streamptr->curTsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
  recID   = recordNewEntry(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
  record  = &streamptr->tsteps[tsID].records[recID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185

186
187
  tsteptype = cgribexGetTsteptype(ISEC1_TimeRange);
  numavg    = ISEC1_AvgNum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
189
190
191

  level1  = ISEC1_Level1;
  level2  = ISEC1_Level2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
  /* fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, ISEC1_LevelType); */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
194
195

  (*record).size     = recsize;
  (*record).position = position;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
  (*record).param    = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
198
199
200
  (*record).ilevel   = level1;
  (*record).ilevel2  = level2;
  (*record).ltype    = ISEC1_LevelType;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
  gridtype = cgribexGetGridType(isec2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202

Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
205
206
207
208
209
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }

210
  memset(&grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
212
213
214
215
216
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
217
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
220
221
222
223
224
225
226
227
228
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
	grid.size  = ISEC4_NumValues;
	grid.xsize = ISEC2_NumLon;
	grid.ysize = ISEC2_NumLat;
	grid.xinc  = 0;
	grid.yinc  = 0;
	grid.xdef  = 0;
	/* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
	  {
	    if ( grid.xsize > 1 )
	      {
229
230
231
		if ( ISEC2_ResFlag && ISEC2_LonIncr > 0 )
		  if ( ISEC2_LonIncr*(grid.xsize-1) > 360000 ) ISEC2_LonIncr = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
233
234
235
236
237
238
239
240
241
242
243
244
		if ( ISEC2_ResFlag && ISEC2_LonIncr > 0 )
		  grid.xinc = ISEC2_LonIncr * 0.001;
		else
		  grid.xinc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid.xsize - 1);

		/* correct xinc if necessary */
		if ( ISEC2_FirstLon == 0 && ISEC2_LastLon > 354000 )
		  {
		    double xinc = 360. / grid.xsize;

		    if ( fabs(grid.xinc-xinc) > 0.0 )
		      {
			grid.xinc = xinc;
245
			if ( CDI_Debug ) Message("set xinc to %g", grid.xinc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
247
248
249
250
251
252
253
254
255
256
257
		      }
		  }
	      }
	    grid.xfirst = ISEC2_FirstLon * 0.001;
	    grid.xlast  = ISEC2_LastLon  * 0.001;
	    grid.xdef   = 2;	    
	  }
	grid.ydef  = 0;
	/* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
	  {
	    if ( grid.ysize > 1 )
	      {
258
259
260
		if ( ISEC2_ResFlag && ISEC2_LatIncr > 0 )
		  if ( ISEC2_LatIncr*(grid.ysize-1) > 180000 ) ISEC2_LatIncr = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
		if ( ISEC2_ResFlag && ISEC2_LatIncr > 0 )
		  grid.yinc = ISEC2_LatIncr * 0.001;
		else
		  grid.yinc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid.ysize - 1);
	      }
	    grid.yfirst = ISEC2_FirstLat * 0.001;
	    grid.ylast  = ISEC2_LastLat  * 0.001;
	    grid.ydef   = 2;	    
	  }
	break;
      }
    case GRID_GAUSSIAN_REDUCED:
      {
	grid.size   = ISEC4_NumValues;
        grid.rowlon = ISEC2_RowLonPtr;
	grid.ysize  = ISEC2_NumLat;
	grid.xinc  = 0;
	grid.yinc  = 0;
	grid.xdef  = 0;
	/* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
	  {
	    if ( grid.xsize > 1 )
	      {
		if ( ISEC2_ResFlag && ISEC2_LonIncr > 0 )
		  grid.xinc = ISEC2_LonIncr * 0.001;
		else
		  grid.xinc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid.xsize - 1);
	      }
	    grid.xfirst = ISEC2_FirstLon * 0.001;
	    grid.xlast  = ISEC2_LastLon  * 0.001;
	    grid.xdef   = 2;	    
	  }
	grid.ydef  = 0;
	/* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
	  {
	    if ( grid.ysize > 1 )
	      {
		if ( ISEC2_ResFlag && ISEC2_LatIncr > 0 )
		  grid.yinc = ISEC2_LatIncr * 0.001;
		else
		  grid.yinc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid.ysize - 1);
	      }
	    grid.yfirst = ISEC2_FirstLat * 0.001;
	    grid.ylast  = ISEC2_LastLat  * 0.001;
	    grid.ydef   = 2;	    
	  }
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
    case GRID_LCC:
310
311
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
312
	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
313
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314

315
316
317
318
	grid.size  = ISEC4_NumValues;
	grid.xsize = ISEC2_NumLon;
	grid.ysize = ISEC2_NumLat;

319
320
	grid.lcc_xinc      = ISEC2_Lambert_dx;
	grid.lcc_yinc      = ISEC2_Lambert_dy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
322
	grid.lcc_originLon = ISEC2_FirstLon * 0.001;
	grid.lcc_originLat = ISEC2_FirstLat * 0.001;
323
324
325
326
327
	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;
328

329
	grid.xdef   = 0;	    
330
331
332
333
	grid.ydef   = 0;

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
336
337
    case GRID_SPECTRAL:
      {
	grid.size  = ISEC4_NumValues;
	grid.trunc = ISEC2_PentaJ;
338
339
340
341
342
	if ( ISEC2_RepMode == 2 )
	  grid.lcomplex = 1;
	else
	  grid.lcomplex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
	break;
      }
    case GRID_GME:
      {
	grid.size  = ISEC4_NumValues;
	grid.nd    = ISEC2_GME_ND;
	grid.ni    = ISEC2_GME_NI;
	grid.ni2   = ISEC2_GME_NI2;
	grid.ni3   = ISEC2_GME_NI3;
	break;
      }
    case GRID_GENERIC:
      {
	grid.size  = ISEC4_NumValues;
	grid.xsize = 0;
	grid.ysize = 0;
	break;
      }
    default:
      {
363
	Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
364
365
366
367
368
	break;
      }
    }

  grid.isRotated = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
  if ( cgribexGetIsRotated(isec2) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
371
372
373
374
375
376
377
378
379
380
381
382
    {
      grid.isRotated = TRUE;
      grid.ypole     = - ISEC2_LatSP * 0.001;
      grid.xpole     =   ISEC2_LonSP * 0.001 - 180;
      grid.angle     = 0;
    }

  grid.xvals = NULL;
  grid.yvals = NULL;
  grid.type  = gridtype;

  gridID = varDefGrid(vlistID, grid, 0);

383
  zaxistype = grib1ltypeToZaxisType(ISEC1_LevelType);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384

385
  if ( zaxistype == ZAXIS_HYBRID || zaxistype == ZAXIS_HYBRID_HALF )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
387
388
389
390
391
392
    {
      int vctsize = ISEC2_NumVCP;
      double *vctptr = &fsec2[10];

      varDefVCT(vctsize, vctptr);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
393
  lbounds = cgribexGetZaxisHasBounds(ISEC1_LevelType);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394

395
396
  if ( prec > 32 ) prec = DATATYPE_PACK32;
  if ( prec <  0 ) prec = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397

Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2,
399
	       prec, &varID, &levelID, tsteptype, numavg, ISEC1_LevelType, NULL, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
401
402
403

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

404
  varDefZtype(varID, ztype);
405

406
407
  if ( lmv ) varDefMissval(varID, FSEC3_MissVal);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
408
  if ( varInqInst(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
410
411
412
413
    {
      int center, subcenter, instID;
      center    = ISEC1_CenterID;
      subcenter = ISEC1_SubCenterID;
      instID    = institutInq(center, subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
416
417
418
	instID = institutDef(center, subcenter, NULL, NULL);
      varDefInst(varID, instID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
  if ( varInqModel(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
421
422
    {
      int modelID;
      modelID = modelInq(varInqInst(varID), ISEC1_ModelID, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423
      if ( modelID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
424
425
426
427
	modelID = modelDef(varInqInst(varID), ISEC1_ModelID, NULL);
      varDefModel(varID, modelID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
  if ( varInqTable(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
430
    {
      int tableID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431

Uwe Schulzweida's avatar
Uwe Schulzweida committed
432
      tableID = tableInq(varInqModel(varID), ISEC1_CodeTable, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433

Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
      if ( tableID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
436
437
438
	tableID = tableDef(varInqModel(varID), ISEC1_CodeTable, NULL);
      varDefTable(varID, tableID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
440
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
442

  if ( CDI_Debug )
443
    Message("varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
445
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447

Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
#if  defined  (HAVE_LIBCGRIBEX)
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
static
void MCH_get_undef(int *isec1, double *undef_pds, double *undef_eps)
{
  /* 2010-01-13: Oliver Fuhrer */
  if ( ISEC1_CenterID == 215 ) {
    if (isec1[34] != 0 && isec1[34] != 255) {
      if (isec1[34] & 2) {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,-isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,-isec1[35]);
        }
        *undef_eps = pow(10.0,-isec1[35]-1);
      } else {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,+isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,+isec1[35]);
        }
        *undef_eps = pow(10.0,isec1[35]-1);
      }
    }
  }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
#endif
474

Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
#if  defined  (HAVE_LIBCGRIBEX)
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
static
void cgribexDecodeHeader(int *isec0, int *isec1, int *isec2, double *fsec2,
			 int *isec3, double *fsec3, int *isec4, double *fsec4, 
			 int *gribbuffer, int recsize, int *lmv)
{
  int iret = 0, ipunp = 0, iword = 0;

  gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
	   ipunp, (int *) gribbuffer, recsize, &iword, "J", &iret);

  *lmv = 0;

  if ( ISEC1_CenterID == 215 && (isec1[34] != 0 && isec1[34] != 255) )
    {
      double undef_pds, undef_eps;

      MCH_get_undef(isec1, &undef_pds, &undef_eps);
      FSEC3_MissVal = undef_pds;
      *lmv = 1;
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
498

499
int cgribexScanTimestep1(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
503
  int *isec0, *isec1, *isec2, *isec3, *isec4;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
504
  int lmv = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505
506
507
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508
  int rstatus;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
509
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
510
511
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
512
513
514
515
516
517
518
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
  int prec;
  long recsize = 0;
519
  int warn_time = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
521
522
  int warn_numavg = TRUE;
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
523
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524
  int vlistID;
525
  int ztype;
526
  long unzipsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
528
  compvar_t compVar, compVar0;
  stream_t *streamptr;
529
530
  extern int cdiSkipRecords;
  int nskip = cdiSkipRecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
531
532

  streamptr = stream_to_pointer(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533

534
  stream_check_ptr(__func__, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
535

Uwe Schulzweida's avatar
Uwe Schulzweida committed
536
  streamptr->curTsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537

Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
539
540
541
542
  isec0 = streamptr->record->sec0;
  isec1 = streamptr->record->sec1;
  isec2 = streamptr->record->sec2;
  isec3 = streamptr->record->sec3;
  isec4 = streamptr->record->sec4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
544

  tsID  = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547

  if ( tsID != 0 )
548
    Error("Internal problem! tstepsNewEntry returns %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
551

  fileID = streamInqFileID(streamID);

552
553
554
  while ( nskip-- > 0 )
    {
      recsize = gribGetSize(fileID);
555
      if ( recsize == 0 )
556
	Error("Skipping of %d records failed!", cdiSkipRecords);
557

558
559
560
561
      recpos  = fileGetPos(fileID);
      fileSetPos(fileID, recsize, SEEK_CUR);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
563
564
565
566
567
568
569
  nrecs = 0;
  while ( TRUE )
    {
      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
570
	  if ( nrecs == 0 )
571
	    Error("No GRIB records found!");
572

Uwe Schulzweida's avatar
Uwe Schulzweida committed
573
	  streamptr->ntsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574
575
576
577
578
579
580
581
582
	  break;
	}
      if ( recsize > buffersize )
	{
	  buffersize = recsize;
	  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	}

      readsize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
583
584
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585

586
      ztype = COMPRESS_NONE;
587
588
      if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	{
589
	  ztype = COMPRESS_SZIP;
590
591
592
593
594
595
596
597
	  unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	  if ( (long) buffersize < unzipsize )
	    {
	      buffersize = unzipsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
	}

598
599
      cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
			  (int *) gribbuffer, recsize, &lmv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
600

Uwe Schulzweida's avatar
Uwe Schulzweida committed
601
      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
602
603
      if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
      if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604
605
606
      level1   = ISEC1_Level1;
      level2   = ISEC1_Level2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
      gribDateTime(isec1, &vdate, &vtime);

609
610
611
612
      if ( ISEC4_NumBits > 0 && ISEC4_NumBits <= 32 )
	prec = ISEC4_NumBits;
      else
        prec = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
614
615
616
617
618
619

      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	  rdate = gribRefDate(isec1);
	  rtime = gribRefTime(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
620
	  tunit = cgribexGetTimeUnit(isec1);
621
	  fcast = cgribexTimeIsFC(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
623
624
625
626
	}
      else
	{
	  datetime.date  = vdate;
	  datetime.time  = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
627
	  compVar.param  = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
628
629
630
631
632
          compVar.level1 = level1;
          compVar.level2 = level2;
          compVar.ltype  = ISEC1_LevelType;
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
	      compVar0.param  = streamptr->tsteps[0].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
634
635
636
	      compVar0.level1 = streamptr->tsteps[0].records[recID].ilevel;
	      compVar0.level2 = streamptr->tsteps[0].records[recID].ilevel2;
	      compVar0.ltype  = streamptr->tsteps[0].records[recID].ltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637

Uwe Schulzweida's avatar
Uwe Schulzweida committed
638
	      if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
639
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
640
641
642
643
644
645
646

	  if ( cdiInventoryMode == 1 )
	    {
	      if ( recID < nrecs ) break;
	      if ( warn_time )
		if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 )
		  {
647
648
		    char paramstr[32];
		    cdiParamToString(param, paramstr, sizeof(paramstr));
649
		    Warning("Inconsistent verification time (param=%s level=%d)", paramstr, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
651
652
653
654
655
656
657
658
		    warn_time = FALSE;
		  }
	    }
	  else
	    {
	      if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;

	      if ( recID < nrecs )
		{
659
660
		  char paramstr[32];
		  cdiParamToString(param, paramstr, sizeof(paramstr));
661
		  Warning("Param=%s level=%d already exist, skipped!", paramstr, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
663
664
		  continue;
		}
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
666
667
668
669
670
	}

      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
	    {
671
	      Warning("Changing numavg from %d to %d not supported!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
674
675
676
677
678
679
680
681
682
683
		      taxis->numavg, ISEC1_AvgNum);
	      warn_numavg = FALSE;
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}

      nrecs++;

      if ( CDI_Debug )
684
	Message("%4d %8d %4d  %8d %8d %6d", nrecs, (int)recpos, param, level1, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
685

686
687
      cgribexAddRecord(streamID, param, isec1, isec2, fsec2, fsec3,
		       isec4, recsize, recpos, prec, ztype, lmv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
688
689
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
690
  streamptr->rtsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
691

692
693
  if ( nrecs == 0 ) return (CDI_EUFSTRUCT);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
695
696
697
  cdiGenVars(streamID);

  if ( fcast )
    {
698
      taxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699
700
701
702
703
704
705
      taxis->type  = TAXIS_RELATIVE;
      taxis->rdate = rdate;
      taxis->rtime = rtime;
      taxis->unit  = tunit;
    }
  else
    {
706
      taxisID = taxisCreate(TAXIS_ABSOLUTE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
707
      taxis->type  = TAXIS_ABSOLUTE;
708
      taxis->unit  = tunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
709
710
711
712
713
714
715
716
    }

  taxis->vdate = datetime0.date;
  taxis->vtime = datetime0.time;

  vlistID = streamInqVlist(streamID);
  vlistDefTaxis(vlistID, taxisID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
717
718
  nrecords = streamptr->tsteps[0].nallrecs;
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
719
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
720
721
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
722
      (record_t *) realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
723
724
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
725
726
  streamptr->tsteps[0].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[0].nrecs = nrecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
    streamptr->tsteps[0].recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729

Uwe Schulzweida's avatar
Uwe Schulzweida committed
730
731
  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
732

Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
734
735
    {
      tsID = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
736
      if ( tsID != streamptr->rtsteps )
737
	Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
738

Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
740
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
742
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
743
  if ( streamptr->ntsteps == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
744
745
746
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
748
	  streamptr->ntsteps = 0;
	  for ( varID = 0; varID < streamptr->nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
751
752
753
	    {
	      vlistDefVarTime(vlistID, varID, TIME_CONSTANT);
	    }
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754
755
756
#else
  Error("CGRIBEX support not compiled in!");
#endif
757
758

  return (0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759
760
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
761
762

int cgribexScanTimestep2(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
763
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
765
  int rstatus = 0;
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
766
767
  int *isec0, *isec1, *isec2, *isec3, *isec4;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
768
  int lmv = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
769
770
771
772
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
773
774
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
775
  DateTime datetime, datetime0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
776
777
778
779
780
781
  int tsID;
  int varID, gridID;
  size_t readsize;
  int nrecords, nrecs, recID, rindex;
  long recsize = 0;
  int warn_numavg = TRUE;
782
  int tsteptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
783
  int taxisID = -1;
784
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
  int vlistID;
786
  long unzipsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
787
788
  compvar_t compVar, compVar0;
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
789

Uwe Schulzweida's avatar
Uwe Schulzweida committed
790
  streamptr = stream_to_pointer(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
791

792
  stream_check_ptr(__func__, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793

Uwe Schulzweida's avatar
Uwe Schulzweida committed
794
795
796
797
798
799
800
  streamptr->curTsID = 1;

  isec0 = streamptr->record->sec0;
  isec1 = streamptr->record->sec1;
  isec2 = streamptr->record->sec2;
  isec3 = streamptr->record->sec3;
  isec4 = streamptr->record->sec4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801
802
803
804
805

  fileID  = streamInqFileID(streamID);
  vlistID = streamInqVlist(streamID);
  taxisID = vlistInqTaxis(vlistID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
806
807
  gribbuffer = (unsigned char *) streamptr->record->buffer;
  buffersize = streamptr->record->buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
808
                                          
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
  tsID = streamptr->rtsteps;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810
  if ( tsID != 1 )
811
    Error("Internal problem! unexpeceted timestep %d", tsID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812

Uwe Schulzweida's avatar
Uwe Schulzweida committed
813
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
814

Uwe Schulzweida's avatar
Uwe Schulzweida committed
815
  fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
817
818

  cdiCreateRecords(streamID, tsID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
819
820
821
  nrecords = streamptr->tsteps[tsID].nallrecs;
  streamptr->tsteps[1].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[1].nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
822
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
823
    streamptr->tsteps[1].recIDs[recID] = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824
825
826
      
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827
828
829
830
831
      varID = streamptr->tsteps[0].records[recID].varID;
      streamptr->tsteps[tsID].records[recID].position = 
	streamptr->tsteps[0].records[recID].position;
      streamptr->tsteps[tsID].records[recID].size     = 
	streamptr->tsteps[0].records[recID].size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
832
833
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
834
835
  rindex = 0;
  while ( TRUE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
837
838
      if ( rindex > nrecords ) break;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
840
841
842
      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);
      if ( recsize == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843
	  streamptr->ntsteps = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
844
845
846
847
848
	  break;
	}
      if ( recsize > buffersize )
	{
	  buffersize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
849
	  gribbuffer = (unsigned char *) realloc(gribbuffer, (size_t)buffersize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
850
851
852
	}

      readsize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
853
854
      rstatus = gribRead(fileID, gribbuffer, &readsize);
      if ( rstatus ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855

856
857
858
859
860
861
862
863
864
865
      if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	{
	  unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	  if ( (long) buffersize < unzipsize )
	    {
	      buffersize = unzipsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
	}

866
867
      cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
			  (int *) gribbuffer, recsize, &lmv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
868

Uwe Schulzweida's avatar
Uwe Schulzweida committed
869
      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
870
871
      if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
      if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
      level1    = ISEC1_Level1;
      level2    = ISEC1_Level2;

      gribDateTime(isec1, &vdate, &vtime);

      if ( rindex == 0 )
	{
	  if ( taxisInqType(taxisID) == TAXIS_RELATIVE )
	    {
	      taxis->type  = TAXIS_RELATIVE;
	      taxis->rdate = gribRefDate(isec1);
	      taxis->rtime = gribRefTime(isec1);
	    }
	  else
	    {
	      taxis->type  = TAXIS_ABSOLUTE;
	    }
889
	  taxis->unit  = cgribexGetTimeUnit(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
890
891
	  taxis->vdate = vdate;
	  taxis->vtime = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
892
893
894

	  datetime0.date = vdate;
	  datetime0.time = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
895
896
	}

897
898
      tsteptype = cgribexGetTsteptype(ISEC1_TimeRange);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
900
901
902
903
904
      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg &&
		(taxis->numavg != ISEC1_AvgNum) )
	    {
	  /*
905
	      Warning("Changing numavg from %d to %d not supported!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
906
907
908
909
910
911
912
913
914
915
		      taxis->numavg, ISEC1_AvgNum);
	  */
	      warn_numavg = FALSE;
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
916
917
      datetime.date  = vdate;
      datetime.time  = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
918
      compVar.param  = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
919
920
921
922
923
      compVar.level1 = level1;
      compVar.level2 = level2;
      compVar.ltype  = ISEC1_LevelType;
      for ( recID = 0; recID < nrecords; recID++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
924
	  compVar0.param  = streamptr->tsteps[tsID].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
925
926
927
	  compVar0.level1 = streamptr->tsteps[tsID].records[recID].ilevel;
	  compVar0.level2 = streamptr->tsteps[tsID].records[recID].ilevel2;
	  compVar0.ltype  = streamptr->tsteps[tsID].records[recID].ltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928

Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
	  if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
930
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
931

Uwe Schulzweida's avatar
Uwe Schulzweida committed
932
933
      if ( recID == nrecords )
	{
934
935
	  char paramstr[32];
	  cdiParamToString(param, paramstr, sizeof(paramstr));
936
	  Warning("Param=%s level=%d not defined at timestep 1!", paramstr, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
937
938
939
	  return (CDI_EUFSTRUCT);
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
      if ( cdiInventoryMode == 1 )
	{
	  if ( streamptr->tsteps[tsID].records[recID].used )
	    {
	      break;
	    }
	  else
	    {
	      streamptr->tsteps[tsID].records[recID].used = TRUE;
	      streamptr->tsteps[tsID].recIDs[rindex] = recID;
	    }    
	}
      else
	{
	  if ( streamptr->tsteps[tsID].records[recID].used )
	    {
956
957
958
	      char paramstr[32];
	      cdiParamToString(param, paramstr, sizeof(paramstr));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
959
960
	      if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;

961
	      Warning("Param=%s level=%d already exist, skipped!", paramstr, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
963
964
965
966
967
968
969
	      continue;
	    }
	  else
	    {
	      streamptr->tsteps[tsID].records[recID].used = TRUE;
	      streamptr->tsteps[tsID].recIDs[rindex] = recID;
	    }    
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
970
971

      if ( CDI_Debug )
972
	Message("%4d %8d %4d %8d %8d %6d", rindex+1, (int)recpos, param, level1, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
973

Uwe Schulzweida's avatar
Uwe Schulzweida committed
974
      streamptr->tsteps[tsID].records[recID].size = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
975

Uwe Schulzweida's avatar
Uwe Schulzweida committed
976
      compVar0.param  = streamptr->tsteps[tsID].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
977
978
979
      compVar0.level1 = streamptr->tsteps[tsID].records[recID].ilevel;
      compVar0.level2 = streamptr->tsteps[tsID].records[recID].ilevel2;
      compVar0.ltype  = streamptr->tsteps[tsID].records[recID].ltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
980

Uwe Schulzweida's avatar
Uwe Schulzweida committed
981
      if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
982
	{
983
	  Message("tsID = %d recID = %d param = %3d new %3d  level = %3d new %3d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
984
		  tsID, recID,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
		  streamptr->tsteps[tsID].records[recID].param, param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
986
		  streamptr->tsteps[tsID].records[recID].ilevel, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
987
988
989
	  return (CDI_EUFSTRUCT);
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
990
991
      streamptr->tsteps[1].records[recID].position = recpos;
      varID = streamptr->tsteps[tsID].records[recID].varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
992
993
994
      gridID = vlistInqVarGrid(vlistID, varID);
      if ( gridInqSize(gridID) == 1 && gridInqType(gridID) == GRID_LONLAT )
	{
995
996
	  if ( IS_NOT_EQUAL(gridInqXval(gridID, 0),ISEC2_FirstLon*0.001) ||
	       IS_NOT_EQUAL(gridInqYval(gridID, 0),ISEC2_FirstLat*0.001) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
997
998
	    gridChangeType(gridID, GRID_TRAJECTORY);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
999

1000
1001
1002
      if ( tsteptype != vlistInqVarTsteptype(vlistID, varID) )
	vlistDefVarTsteptype(vlistID, varID, tsteptype);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1003
      rindex++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1004
1005
1006
1007
1008
    }

  nrecs = 0;
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
      if ( ! streamptr->tsteps[tsID].records[recID].used )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1010
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1011
	  varID = streamptr->tsteps[tsID].records[recID].varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1012
1013
1014
1015
1016
1017
1018
	  vlistDefVarTime(vlistID, varID, TIME_CONSTANT);
	}
      else
	{
	  nrecs++;
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1019
  streamptr->tsteps[tsID].nrecs = nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1020

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1021
  streamptr->rtsteps = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1022

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1023
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
    {
      tsID = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1026
      if ( tsID != streamptr->rtsteps )
1027
	Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1028

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1029
1030
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1031
1032
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
1034
  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1036

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1037
  return (rstatus);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1038
1039
1040
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
1041
int cgribexScanTimestep(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1042
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1043
1044
  int rstatus = 0;
#if  defined  (HAVE_LIBCGRIBEX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1045
1046
  int *isec0, *isec1, *isec2, *isec3, *isec4;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
1047
  int lmv = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1048
1049
1050
1051
  long recsize = 0;
  off_t recpos = 0;
  unsigned char *gribbuffer;
  long buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1052
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1053
1054
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1055
1056
1057
1058
1059
1060
  DateTime datetime, datetime0;
  int tsID;
  int vrecID, recID;
  int warn_numavg = TRUE;
  size_t readsize;
  int taxisID = -1;
1061
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1062
1063
1064
1065
  int vlistID;
  int rindex, nrecs = 0;
  long unzipsize;
  compvar_t compVar, compVar0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1066
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1067
1068

  streamptr = stream_to_pointer(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1069

1070
  stream_check_ptr(__func__, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1071

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1072
  vlistID = streamInqVlist(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1073

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1074
1075
  if ( CDI_Debug )
    {
1076
1077
1078
1079
      Message("streamID = %d", streamID);
      Message("cts = %d", streamptr->curTsID);
      Message("rts = %d", streamptr->rtsteps);
      Message("nts = %d", streamptr->ntsteps);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1080
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1081

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1082
1083
1084
1085
1086
  isec0 = streamptr->record->sec0;
  isec1 = streamptr->record->sec1;
  isec2 = streamptr->record->sec2;
  isec3 = streamptr->record->sec3;
  isec4 = streamptr->record->sec4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1087

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1088
1089
  tsID  = streamptr->rtsteps;
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1090

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1091
1092
1093
1094
  if ( streamptr->tsteps[tsID].recordSize == 0 )
    {
      gribbuffer = (unsigned char *) streamptr->record->buffer;
      buffersize = streamptr->record->buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1095

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1096
      cdiCreateRecords(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1097

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1098
      nrecs = streamptr->tsteps[1].nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1099

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1100
1101
1102
1103
      streamptr->tsteps[tsID].nrecs = nrecs;
      streamptr->tsteps[tsID].recIDs = (int *) malloc(nrecs*sizeof(int));
      for ( recID = 0; recID < nrecs; recID++ )
	streamptr->tsteps[tsID].recIDs[recID] = streamptr->tsteps[1].recIDs[recID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1104

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1105
      fileID = streamInqFileID(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1106

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1107
      fileSetPos(fileID, streamptr->tsteps[tsID].position, SEEK_SET);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1108

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1109
1110
      rindex = 0;
      while ( TRUE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1111
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1112
	  if ( rindex > nrecs ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1113

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
	  recsize = gribGetSize(fileID);
	  recpos  = fileGetPos(fileID);
	  if ( recsize == 0 )
	    {
	      streamptr->ntsteps = streamptr->rtsteps + 1;
	      break;
	    }
	  if ( recsize > buffersize )
	    {
	      buffersize = recsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1126

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1127
	  readsize = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1128
1129
	  rstatus = gribRead(fileID, gribbuffer, &readsize);
	  if ( rstatus )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1130
	    {
1131
	      Error("Inconsistent timestep %d (GRIB record %d/%d)!\n", tsID+1, rindex+1,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
		    streamptr->tsteps[tsID].recordSize);
	      break;
	    }

	  if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	    {
	      unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	      if ( (long) buffersize < unzipsize )
		{
		  buffersize = unzipsize;
		  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
		}
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1145

1146
1147
	  cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
			      (int *) gribbuffer, recsize, &lmv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1148