stream_cgribex.c 50 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
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
int cgribexGetGridType(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
34
  /*  static char func[] = "cgribexGetGridType"; */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
37
38
  int gridtype = 0;

  switch (ISEC2_GridType)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
    case  GTYPE_LATLON:
    case  GTYPE_LATLON_ROT:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
43
44
      {
	gridtype = GRID_LONLAT;
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
    case  GTYPE_LCC:
46
      {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
	gridtype = GRID_LCC;
48
49
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
    case  GTYPE_GAUSSIAN:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
52
53
54
55
56
57
      {
	if ( ISEC2_Reduced )
	  gridtype = GRID_GAUSSIAN_REDUCED;
	else
	  gridtype = GRID_GAUSSIAN;
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
    case  GTYPE_SPECTRAL:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
61
62
      {
	gridtype = GRID_SPECTRAL;
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
    case  GTYPE_GME:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
68
69
70
71
72
73
74
75
76
77
      {
	gridtype = GRID_GME;
	break;
      }
    default:
      {
	gridtype = GRID_GENERIC;
	break;
      }
    }

  return (gridtype);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
int cgribexGetIsRotated(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
  /*  static char func[] = "cgribexGetIsRotated"; */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
84
85
86
87
88
89
90
91
92
  int isRotated = 0;

  if ( ISEC2_GridType == 10 )
    {
      isRotated = 1;
    }

  return (isRotated);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
int cgribexGetZaxisType(int grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
{
96
  int zaxistype = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97

98
  switch ( grb_ltype )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
    {
    case LTYPE_SURFACE:
      {
102
	zaxistype = ZAXIS_SURFACE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
104
	break;
      }
105
106
107
108
109
    case LTYPE_MEANSEA:
      {
	zaxistype = ZAXIS_MEANSEA;
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
111
112
    case LTYPE_99:
    case LTYPE_ISOBARIC:
      {
113
	zaxistype = ZAXIS_PRESSURE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
116
117
	break;
      }
    case LTYPE_HEIGHT:
      {
118
	zaxistype = ZAXIS_HEIGHT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
121
122
	break;
      }
    case LTYPE_ALTITUDE:
      {
123
	zaxistype = ZAXIS_ALTITUDE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
127
    case LTYPE_SIGMA:
      {
128
	zaxistype = ZAXIS_SIGMA;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
133
    case LTYPE_HYBRID:
    case LTYPE_HYBRID_LAYER:
      {
134
	zaxistype = ZAXIS_HYBRID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
137
138
139
	break;
      }
    case LTYPE_LANDDEPTH:
    case LTYPE_LANDDEPTH_LAYER:
      {
140
	zaxistype = ZAXIS_DEPTH_BELOW_LAND;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
	break;
      }
143
144
    case LTYPE_ISENTROPIC:
      {
145
	zaxistype = ZAXIS_ISENTROPIC;
146
147
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
149
    case LTYPE_SEADEPTH:
      {
150
	zaxistype = ZAXIS_DEPTH_BELOW_SEA;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
152
153
154
	break;
      }
    default:
      {
155
	zaxistype = ZAXIS_GENERIC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
158
159
	break;
      }
    }

160
  return (zaxistype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
int cgribexGetZaxisHasBounds(int grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
167
{
  int lbounds = 0;

168
  switch (grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
171
172
173
174
175
176
177
178
179
180
    {
    case LTYPE_HYBRID_LAYER:
    case LTYPE_LANDDEPTH_LAYER:
      {
	lbounds = 1;
	break;
      }
    }

  return (lbounds);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
int cgribexGetTimeUnit(int *isec1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
  static char func[] = "cgribexGetTimeUnit";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
187
188
189
  int timeunit = -1;
  static int lprint = TRUE;

  switch ( ISEC1_TimeUnit )
    {
190
191
192
193
    case ISEC1_TABLE4_MINUTE:  timeunit = TUNIT_MINUTE;  break;
    case ISEC1_TABLE4_QUARTER: timeunit = TUNIT_QUARTER; break;
    case ISEC1_TABLE4_HOUR:    timeunit = TUNIT_HOUR;    break;
    case ISEC1_TABLE4_DAY:     timeunit = TUNIT_DAY;     break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
196
197
198
199
200
201
202
203
204
205
    default:
      if ( lprint )
	{
	  Message(func, "Time unit %d unsupported", ISEC1_TimeUnit);
	  lprint = FALSE;
	}
    }

  return (timeunit);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
static
207
208
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
209
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
  static char func[] = "cgribexAddRecord";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
  int gridtype;
212
  int zaxistype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
  int gridID = CDI_UNDEFID, varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
215
216
217
218
  int levelID = 0;
  int tsID, recID;
  int level1, level2;
  int numavg;
  int lbounds = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
  record_t *record;
220
  grid_t grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
  int vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224

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

  vlistID = streamInqVlist(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
  tsID    = streamptr->curTsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
  recID   = recordNewEntry(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
  record  = &streamptr->tsteps[tsID].records[recID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
231
232
233
234
235

  numavg  = ISEC1_AvgNum;

  level1  = ISEC1_Level1;
  level2  = ISEC1_Level2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
  /* fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, ISEC1_LevelType); */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
238
239

  (*record).size     = recsize;
  (*record).position = position;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
240
  (*record).param    = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
243
244
  (*record).ilevel   = level1;
  (*record).ilevel2  = level2;
  (*record).ltype    = ISEC1_LevelType;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
  gridtype = cgribexGetGridType(isec2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246

Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
252
253
    {
      gridtype = GRID_GAUSSIAN;
      ISEC2_NumLon = 2*ISEC2_NumLat;
      ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
    }

254
  memset(&grid, 0, sizeof(grid_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
256
257
258
259
260
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261
	  Error(func, "numberOfPoints (%d) and gridSize (%d) differ!",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
		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 )
	      {
		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;
			if ( CDI_Debug ) Message(func, "set xinc to %g", grid.xinc);
		      }
		  }
	      }
	    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;
      }
    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
347
    case GRID_LCC:
348
349
      {
	if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
	  Error(func, "numberOfPoints (%d) and gridSize (%d) differ!",
351
		ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352

353
354
355
356
	grid.size  = ISEC4_NumValues;
	grid.xsize = ISEC2_NumLon;
	grid.ysize = ISEC2_NumLat;

357
358
	grid.lcc_xinc      = ISEC2_Lambert_dx;
	grid.lcc_yinc      = ISEC2_Lambert_dy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
360
	grid.lcc_originLon = ISEC2_FirstLon * 0.001;
	grid.lcc_originLat = ISEC2_FirstLat * 0.001;
361
362
363
364
365
	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;
366

367
	grid.xdef   = 0;	    
368
369
370
371
	grid.ydef   = 0;

	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
    case GRID_SPECTRAL:
      {
	grid.size  = ISEC4_NumValues;
	grid.trunc = ISEC2_PentaJ;
	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:
      {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
	Error(func, "%s grid unsupported!", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
398
399
400
401
	break;
      }
    }

  grid.isRotated = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
  if ( cgribexGetIsRotated(isec2) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403
404
405
406
407
408
409
410
411
412
413
414
415
    {
      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);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
  zaxistype = cgribexGetZaxisType(ISEC1_LevelType);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417

418
  if ( zaxistype == ZAXIS_HYBRID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
420
421
422
423
424
425
    {
      int vctsize = ISEC2_NumVCP;
      double *vctptr = &fsec2[10];

      varDefVCT(vctsize, vctptr);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
  lbounds = cgribexGetZaxisHasBounds(ISEC1_LevelType);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427

428
429
  if ( prec > 32 ) prec = DATATYPE_PACK32;
  if ( prec <  0 ) prec = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430

Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2,
432
	       prec, &varID, &levelID, numavg, ISEC1_LevelType, NULL, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
434
435
436

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

437
  varDefZtype(varID, ztype);
438

439
440
  if ( lmv ) varDefMissval(varID, FSEC3_MissVal);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
  if ( varInqInst(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
443
444
445
446
    {
      int center, subcenter, instID;
      center    = ISEC1_CenterID;
      subcenter = ISEC1_SubCenterID;
      instID    = institutInq(center, subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
449
450
451
	instID = institutDef(center, subcenter, NULL, NULL);
      varDefInst(varID, instID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
452
  if ( varInqModel(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
454
455
    {
      int modelID;
      modelID = modelInq(varInqInst(varID), ISEC1_ModelID, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
      if ( modelID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
458
459
460
	modelID = modelDef(varInqInst(varID), ISEC1_ModelID, NULL);
      varDefModel(varID, modelID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
461
  if ( varInqTable(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
463
    {
      int tableID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
      if ( tableID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468
469
470
471
	tableID = tableDef(varInqModel(varID), ISEC1_CodeTable, NULL);
      varDefTable(varID, tableID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
473
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
474
475

  if ( CDI_Debug )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476
477
    Message(func, "varID = %d  param = %d  zaxistype = %d  gridID = %d  levelID = %d",
	    varID, param, zaxistype, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
479
}

480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
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);
      }
    }
  }
}

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
527
528

void cgribexScanTimestep1(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
  static char func[] = "cgribexScanTimestep1";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
531
532
  int *isec0, *isec1, *isec2, *isec3, *isec4;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
533
  int lmv = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
535
536
537
538
  off_t recpos = 0;
  unsigned char *gribbuffer = NULL;
  long buffersize = 0;
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
539
540
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
541
542
543
544
545
546
547
  DateTime datetime, datetime0;
  int tsID;
  int varID;
  size_t readsize;
  int nrecords, nrecs, recID;
  int prec;
  long recsize = 0;
548
  int warn_time = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
551
552
553
  int warn_numavg = TRUE;
  int taxisID = -1;
  int rdate = 0, rtime = 0, tunit = 0, fcast = 0;
  TAXIS *taxis;
  int vlistID;
554
  int ztype;
555
  long unzipsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
557
  compvar_t compVar, compVar0;
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
559

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
561
  stream_check_ptr(func, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562

Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
  streamptr->curTsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564

Uwe Schulzweida's avatar
Uwe Schulzweida committed
565
566
567
568
569
  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
570
571

  tsID  = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
573
574
575
576
577
578
579
580
581
582
583
584
585
586

  if ( tsID != 0 )
    Error(func, "Internal problem! tstepsNewEntry returns %d", tsID);

  fileID = streamInqFileID(streamID);

  nrecs = 0;
  while ( TRUE )
    {
      recsize = gribGetSize(fileID);
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
587
	  streamptr->ntsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
589
590
591
592
593
594
595
596
597
598
599
	  break;
	}
      if ( recsize > buffersize )
	{
	  buffersize = recsize;
	  gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	}

      readsize = recsize;
      status = gribRead(fileID, gribbuffer, &readsize);
      if ( status ) break;

600
      ztype = COMPRESS_NONE;
601
602
      if ( gribGetZip(recsize, gribbuffer, &unzipsize) > 0 )
	{
603
	  ztype = COMPRESS_SZIP;
604
605
606
607
608
609
610
611
	  unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
	  if ( (long) buffersize < unzipsize )
	    {
	      buffersize = unzipsize;
	      gribbuffer = (unsigned char *) realloc(gribbuffer, buffersize);
	    }
	}

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
616
617
      if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
      if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
619
620
      level1   = ISEC1_Level1;
      level2   = ISEC1_Level2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
621
622
      gribDateTime(isec1, &vdate, &vtime);

623
624
625
626
      if ( ISEC4_NumBits > 0 && ISEC4_NumBits <= 32 )
	prec = ISEC4_NumBits;
      else
        prec = DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
627
628
629
630
631
632
633

      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	  rdate = gribRefDate(isec1);
	  rtime = gribRefTime(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
634
	  tunit = cgribexGetTimeUnit(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635
636
637
638
639
640
	  fcast = gribTimeIsFC(isec1);
	}
      else
	{
	  datetime.date  = vdate;
	  datetime.time  = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
641
	  compVar.param  = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
643
644
645
646
          compVar.level1 = level1;
          compVar.level2 = level2;
          compVar.ltype  = ISEC1_LevelType;
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
647
	      compVar0.param  = streamptr->tsteps[0].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
648
649
650
	      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
651

Uwe Schulzweida's avatar
Uwe Schulzweida committed
652
	      if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
655
656
657
658
659
660

	  if ( cdiInventoryMode == 1 )
	    {
	      if ( recID < nrecs ) break;
	      if ( warn_time )
		if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 )
		  {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
		    Warning(func, "Inconsistent verification time (param %d level %d)", param, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
663
664
665
666
667
668
669
670
		    warn_time = FALSE;
		  }
	    }
	  else
	    {
	      if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;

	      if ( recID < nrecs )
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
		  Warning(func, "Param %d level %d already exist, skipped!", param, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
674
		  continue;
		}
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
	}

      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
	    {
	      Message(func, "Change numavg from %d to %d not allowed!",
		      taxis->numavg, ISEC1_AvgNum);
	      warn_numavg = FALSE;
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}

      nrecs++;

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

696
697
      cgribexAddRecord(streamID, param, isec1, isec2, fsec2, fsec3,
		       isec4, recsize, recpos, prec, ztype, lmv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
698
699
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
700
  streamptr->rtsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
701
702
703
704
705

  cdiGenVars(streamID);

  if ( fcast )
    {
706
      taxisID = taxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
707
708
709
710
711
712
713
      taxis->type  = TAXIS_RELATIVE;
      taxis->rdate = rdate;
      taxis->rtime = rtime;
      taxis->unit  = tunit;
    }
  else
    {
714
      taxisID = taxisCreate(TAXIS_ABSOLUTE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
715
      taxis->type  = TAXIS_ABSOLUTE;
716
      taxis->unit  = tunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
717
718
719
720
721
722
723
724
    }

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
725
726
  nrecords = streamptr->tsteps[0].nallrecs;
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
729
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730
      (record_t *) realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
731
732
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
734
  streamptr->tsteps[0].recIDs = (int *) malloc(nrecords*sizeof(int));
  streamptr->tsteps[0].nrecs = nrecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
735
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
736
    streamptr->tsteps[0].recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
737

Uwe Schulzweida's avatar
Uwe Schulzweida committed
738
739
  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
740

Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
742
743
    {
      tsID = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
744
      if ( tsID != streamptr->rtsteps )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
746
	Error(func, "Internal error. tsID = %d", tsID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
748
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
751
  if ( streamptr->ntsteps == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752
753
754
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
756
	  streamptr->ntsteps = 0;
	  for ( varID = 0; varID < streamptr->nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757
758
759
760
761
762
763
	    {
	      vlistDefVarTime(vlistID, varID, TIME_CONSTANT);
	    }
	}
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
765

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
794
  stream_check_ptr(func, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
795

Uwe Schulzweida's avatar
Uwe Schulzweida committed
796
797
798
799
800
801
802
  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
803
804
805
806
807

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

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

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

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

  cdiCreateRecords(streamID, tsID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
821
822
823
  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
824
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
825
    streamptr->tsteps[1].recIDs[recID] = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
827
828
      
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
830
831
832
833
      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
834
835
    }

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

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

      readsize = recsize;
      status = gribRead(fileID, gribbuffer, &readsize);
      if ( status ) break;

858
859
860
861
862
863
864
865
866
867
      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);
	    }
	}

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
871
      param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
872
873
      if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
      if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
874
875
876
877
878
879
880
881
882
883
884
885
      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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
886
	      taxis->unit  = cgribexGetTimeUnit(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887
888
889
890
891
892
893
	    }
	  else
	    {
	      taxis->type  = TAXIS_ABSOLUTE;
	    }
	  taxis->vdate = vdate;
	  taxis->vtime = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
894
895
896

	  datetime0.date = vdate;
	  datetime0.time = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
	}

      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg &&
		(taxis->numavg != ISEC1_AvgNum) )
	    {
	  /*
	      Message(func, "change numavg from %d to %d not allowed!",
		      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 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
934
	  Warning(func, "Param %d level %d not defined at timestep 1!", param, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
935
936
937
	  return (CDI_EUFSTRUCT);
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
938
939
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 )
	    {
	      if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
956
	      Warning(func, "Param %d level %d already exist, skipped!", param, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
959
960
961
962
963
964
	      continue;
	    }
	  else
	    {
	      streamptr->tsteps[tsID].records[recID].used = TRUE;
	      streamptr->tsteps[tsID].recIDs[rindex] = recID;
	    }    
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
965
966

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
971
      compVar0.param  = streamptr->tsteps[tsID].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
973
974
      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
975

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
986
      streamptr->tsteps[1].records[recID].position = recpos;
      varID = streamptr->tsteps[tsID].records[recID].varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
987
988
989
      gridID = vlistInqVarGrid(vlistID, varID);
      if ( gridInqSize(gridID) == 1 && gridInqType(gridID) == GRID_LONLAT )
	{
990
991
	  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
992
993
	    gridChangeType(gridID, GRID_TRAJECTORY);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
994
995

      rindex++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
996
997
998
999
1000
    }

  nrecs = 0;
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1001
      if ( ! streamptr->tsteps[tsID].records[recID].used )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1003
	  varID = streamptr->tsteps[tsID].records[recID].varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1004
1005
1006
1007
1008
1009
1010
	  vlistDefVarTime(vlistID, varID, TIME_CONSTANT);
	}
      else
	{
	  nrecs++;
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1011
  streamptr->tsteps[tsID].nrecs = nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1012

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1013
  streamptr->rtsteps = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1014

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1015
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1016
1017
    {
      tsID = tstepsNewEntry(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1018
      if ( tsID != streamptr->rtsteps )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1019
1020
	Error(func, "Internal error. tsID = %d", tsID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1021
1022
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1023
1024
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1025
1026
  streamptr->record->buffer     = gribbuffer;
  streamptr->record->buffersize = buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
1028
1029
1030
1031

  return (0);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032
int cgribexScanTimestep(int streamID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1034
1035
1036
  static char func[] = "cgribexScanTimestep";
  int *isec0, *isec1, *isec2, *isec3, *isec4;
  double fsec2[512], fsec3[2], *fsec4 = NULL;
1037
  int lmv = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1038
1039
1040
1041
1042
  long recsize = 0;
  off_t recpos = 0;
  unsigned char *gribbuffer;
  long buffersize = 0;
  int status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1043
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1044
1045
  int param = 0;
  int level1 = 0, level2 = 0, vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
  DateTime datetime, datetime0;
  int tsID;
  int vrecID, recID;
  int warn_numavg = TRUE;
  size_t readsize;
  int taxisID = -1;
  TAXIS *taxis;
  int vlistID;
  int rindex, nrecs = 0;
  long unzipsize;
  compvar_t compVar, compVar0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1057
  stream_t *streamptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1058
1059

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1061
  stream_check_ptr(func, streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1062

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1063
  vlistID = streamInqVlist(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1064

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1065
1066
1067
1068
1069
1070
1071
  if ( CDI_Debug )
    {
      Message(func, "streamID = %d", streamID);
      Message(func, "cts = %d", streamptr->curTsID);
      Message(func, "rts = %d", streamptr->rtsteps);
      Message(func, "nts = %d", streamptr->ntsteps);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1072

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1073
1074
1075
1076
1077
  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
1078

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1079
1080
  tsID  = streamptr->rtsteps;
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1081

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1082
1083
1084
1085
  if ( streamptr->tsteps[tsID].recordSize == 0 )
    {
      gribbuffer = (unsigned char *) streamptr->record->buffer;
      buffersize = streamptr->record->buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1086

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1087
      cdiCreateRecords(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1088

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1091
1092
1093
1094
      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
1095

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1096
      fileID = streamInqFileID(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1097

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1100
1101
      rindex = 0;
      while ( TRUE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1102
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1103
	  if ( rindex > nrecs ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1104

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
	  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
1117

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
	  readsize = recsize;
	  status = gribRead(fileID, gribbuffer, &readsize);
	  if ( status )
	    {
	      Error(func, "Inconsistent timestep %d (GRIB record %d/%d)!\n", tsID+1, rindex+1,
		    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
1136

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1140
	  param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
	  if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
	  if ( ISEC1_LevelType ==  99 ) ISEC1_LevelType = 100;
	  level1   = ISEC1_Level1;
	  level2   = ISEC1_Level2;

	  gribDateTime(isec1, &vdate, &vtime);

	  if ( rindex == nrecs ) break;

	  if ( rindex == 0 )
	    {
	      taxisID = vlistInqTaxis(vlistID);
	      if ( taxisInqType(taxisID) == TAXIS_RELATIVE )
		{
		  taxis->type  = TAXIS_RELATIVE;
		  taxis->rdate = gribRefDate(isec1);
		  taxis->rtime = gribRefTime(isec1);
		  taxis->unit  = cgribexGetTimeUnit(isec1);
		}
	      else
		{
		  taxis->type  = TAXIS_ABSOLUTE;
		}
	      taxis->vdate = vdate;
	      taxis->vtime = vtime;

	      datetime0.date = vdate;
	      datetime0.time = vtime;
	    }

	  if ( ISEC1_AvgNum )
	    {
	      if (  taxis->numavg && warn_numavg &&
		   (taxis->numavg != ISEC1_AvgNum) )
		{
	      /*
		  Message(func, "Change numavg from %d to %d not allowed!",
			  streamptr->tsteps[tsID].taxis.numavg, ISEC1_AvgNum);
	      */
		  warn_numavg = FALSE;
		}
	      else
		{
		  taxis->numavg = ISEC1_AvgNum;
		}
	    }
	  
	  datetime.date  = vdate;
	  datetime.time  = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1190
	  compVar.param  = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1191
1192
1193
1194
1195
1196
          compVar.level1 = level1;
          compVar.level2 = level2;
          compVar.ltype  = ISEC1_LevelType;
	  for ( vrecID = 0; vrecID < nrecs; vrecID++ )
	    {
	      recID   = streamptr->tsteps[1].recIDs[vrecID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1197
	      compVar0.param  = streamptr->tsteps[tsID].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1198
1199
1200
1201
1202
1203
1204
1205
1206
	      compVar0.level1 = streamptr->tsteps[tsID].records[recID].ilevel;
	      compVar0.level2 = streamptr->tsteps[tsID].records[recID].ilevel2;
	      compVar0.ltype  = streamptr->tsteps[tsID].records[recID].ltype;

	      if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) == 0 ) break;
	    }

	  if ( vrecID == nrecs )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1207
	      Warning(func, "Param %d level %d not available at timestep %d!", param, level1, tsID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226

	      if ( cdiInventoryMode == 1 )
		return (CDI_EUFSTRUCT);
	      else
		continue;
	    }

	  if ( cdiInventoryMode == 1 )
	    {
	      streamptr->tsteps[tsID].records[recID].used = TRUE;
	      streamptr->tsteps[tsID].recIDs[rindex] = recID;
	    }
	  else
	    {
	      if ( streamptr->tsteps[tsID].records[recID].used )
		{
		  if ( memcmp(&datetime, &datetime0, sizeof(DateTime)) != 0 ) break;
		  
		  if ( CDI_Debug )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1227
		    Warning(func, "Param %d level %d already exist, skipped!", param, level1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1228
1229
1230
1231
1232
1233
1234
1235
1236

		  continue;
		}
	      else
		{
		  streamptr->tsteps[tsID].records[recID].used = TRUE;
		  streamptr->tsteps[tsID].recIDs[rindex] = recID;
		}    
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1237

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1241
	  compVar0.param  = streamptr->tsteps[tsID].records[recID].param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1242
1243
1244
	  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
1245

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1246
1247
	  if ( memcmp(&compVar0, &compVar, sizeof(compvar_t)) != 0 )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1248
	      Message(func, "tsID = %d recID = %d param = %3d new %3d  level = %3d new %3d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1249
		      tsID, recID,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1250
		      streamptr->tsteps[tsID].records[recID].param, param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1251
1252
1253
1254
1255
1256
		      streamptr->tsteps[tsID].records[recID].ilevel, level1);
	      Error(func, "Invalid, unsupported or inconsistent record structure");
	    }
	  
	  streamptr->tsteps[tsID].records[recID].position = recpos;
	  streamptr->tsteps[tsID].records[recID].size = recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1257

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1258
	  if ( CDI_Debug )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1259
	    Message(func, "%4d %8d %4d %8d %8d %6d", rindex, (int)recpos, param, level1, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1260

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1261
1262
	  rindex++;
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1263

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1264
1265
1266
1267
1268
      for ( vrecID = 0; vrecID < nrecs; vrecID++ )
	{
	  recID   = streamptr->tsteps[tsID].recIDs[vrecID];
	  if ( ! streamptr->tsteps[tsID].records[recID].used ) break;
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1269

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1270
1271
      if ( vrecID < nrecs )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1272
1273
	  Warning(func, "Param %d level %d not found at timestep %d!",
		  streamptr->tsteps[tsID].records[recID].param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1274
1275