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

5
#include <limits.h>
Thomas Jahns's avatar
Thomas Jahns committed
6
#include <stdbool.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include <stdio.h>
8
#include <stdlib.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
9
10
11
12
13
14
15
16
17
#include <string.h>
#include <float.h>
#include <math.h>

#include "dmemory.h"

#include "error.h"
#include "file.h"
#include "cdi.h"
18
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19
20
21
#include "varscan.h"
#include "datetime.h"
#include "ieg.h"
22
#include "stream_fcommon.h"
23
#include "stream_ieg.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
25
26
27
28
29
30
31
32
33
34
#include "vlist.h"


#undef  UNDEFID
#define UNDEFID  CDI_UNDEFID

#define SINGLE_PRECISION  4
#define DOUBLE_PRECISION  8

#if defined (HAVE_LIBIEG)

35
36

typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
  int param;
38
  int level;
39
} IEGCOMPVAR;
40
41


Thomas Jahns's avatar
Thomas Jahns committed
42
static int iegInqDatatype(int prec)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
{
  int datatype;

46
47
  if ( prec == DOUBLE_PRECISION ) datatype = DATATYPE_FLT64;
  else                            datatype = DATATYPE_FLT32;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
51
52

  return (datatype);
}


Thomas Jahns's avatar
Thomas Jahns committed
53
static int iegDefDatatype(int datatype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
56
{
  int prec;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
  if ( datatype == DATATYPE_CPX32 || datatype == DATATYPE_CPX64 )
58
    Error("CDI/IEG library does not support complex numbers!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59

60
61
  if ( datatype != DATATYPE_FLT32 && datatype != DATATYPE_FLT64 )
    datatype = DATATYPE_FLT32;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62

63
  if ( datatype == DATATYPE_FLT64 ) prec = DOUBLE_PRECISION;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
68
  else                              prec = SINGLE_PRECISION;

  return (prec);
}

69
/* not used
70
int iegInqRecord(stream_t *streamptr, int *varID, int *levelID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
72
73
{
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
  int icode, ilevel;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
  int zaxisID = -1;
  int vlistID;
77
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78

Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
  vlistID = streamptr->vlistID;
  fileID  = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
82
83
84
85
86
87

  *varID   = -1;
  *levelID = -1;

  status = iegRead(fileID, iegp);
  if ( status != 0 ) return (0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
  icode  = IEG_P_Parameter(iegp->ipdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
91
92
93
  if ( IEG_P_LevelType(iegp->ipdb) == IEG_LTYPE_HYBRID_LAYER )
    ilevel = IEG_P_Level1(iegp->ipdb);
  else
    ilevel = IEG_P_Level2(iegp->ipdb);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
  *varID = vlistInqVarID(vlistID, icode);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

96
  if ( *varID == UNDEFID ) Error("Code %d undefined", icode);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
99
100

  zaxisID = vlistInqVarZaxis(vlistID, *varID);

  *levelID = zaxisInqLevelID(zaxisID, (double) ilevel);
101

Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
  return (1);
}
104
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105

106
void iegReadRecord(stream_t *streamptr, double *data, int *nmiss)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
109
110
111
112
113
114
{
  int vlistID, fileID;
  int status;
  int recID, vrecID, tsID;
  off_t recpos;
  int varID, gridID;
  int i, size;
  double missval;
115
  void *iegp = streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116

Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
118
  vlistID = streamptr->vlistID;
  fileID  = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
121
122
123
  tsID    = streamptr->curTsID;
  vrecID  = streamptr->tsteps[tsID].curRecID;
  recID   = streamptr->tsteps[tsID].recIDs[vrecID];
  recpos  = streamptr->tsteps[tsID].records[recID].position;
  varID   = streamptr->tsteps[tsID].records[recID].varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
126
127

  fileSetPos(fileID, recpos, SEEK_SET);

  status = iegRead(fileID, iegp);
128
129
  if ( status != 0 )
    Error("Could not read IEG record!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132

  iegInqDataDP(iegp, data);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
  missval = vlistInqVarMissval(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
135
  gridID  = vlistInqVarGrid(vlistID, varID);
  size    = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136

Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
  streamptr->numvals += size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
140
141
142
143
144
145
146
147

  *nmiss = 0;
  for ( i = 0; i < size; i++ )
    if ( DBL_IS_EQUAL(data[i], missval) || DBL_IS_EQUAL(data[i], (float)missval) )
      {
	data[i] = missval;
	(*nmiss)++;
      }
}

148
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
int iegGetZaxisType(int iegleveltype)
{
  int leveltype = 0;

  switch ( iegleveltype )
    {
    case IEG_LTYPE_SURFACE:
      {
	leveltype = ZAXIS_SURFACE;
	break;
      }
    case IEG_LTYPE_99:
    case IEG_LTYPE_ISOBARIC:
      {
	leveltype = ZAXIS_PRESSURE;
	break;
      }
    case IEG_LTYPE_HEIGHT:
      {
	leveltype = ZAXIS_HEIGHT;
	break;
      }
    case IEG_LTYPE_ALTITUDE:
      {
	leveltype = ZAXIS_ALTITUDE;
	break;
      }
    case IEG_LTYPE_HYBRID:
    case IEG_LTYPE_HYBRID_LAYER:
      {
	leveltype = ZAXIS_HYBRID;
	break;
      }
    case IEG_LTYPE_LANDDEPTH:
    case IEG_LTYPE_LANDDEPTH_LAYER:
      {
	leveltype = ZAXIS_DEPTH_BELOW_LAND;
	break;
      }
    case IEG_LTYPE_SEADEPTH:
      {
	leveltype = ZAXIS_DEPTH_BELOW_SEA;
	break;
      }
    default:
      {
	leveltype = ZAXIS_GENERIC;
	break;
      }
    }

  return (leveltype);
}


Thomas Jahns's avatar
Thomas Jahns committed
204
static void iegDefTime(int *pdb, int date, int time, int taxisID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
{
206
  int year, month, day, hour, minute, second;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
  int timetype = -1;

209
  if ( taxisID != -1 ) timetype = taxisInqType(taxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
211
212

  if ( timetype == TAXIS_ABSOLUTE || timetype == TAXIS_RELATIVE )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
      cdiDecodeDate(date, &year, &month, &day);
      cdiDecodeTime(time, &hour, &minute, &second);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

      IEG_P_Year(pdb)     = year;
      IEG_P_Month(pdb)    = month;
      IEG_P_Day(pdb)      = day;
      IEG_P_Hour(pdb)     = hour;
      IEG_P_Minute(pdb)   = minute;

      pdb[15] = 1;
      pdb[16] = 0;
      pdb[17] = 0;
      pdb[18] = 10;
      pdb[36] = 1;
    }

  pdb[5] = 128;
}

Thomas Jahns's avatar
Thomas Jahns committed
232
233
234
235
236
/* find smallest power of 10 in [1000,10000000] that upon
 * multiplication results in fractional part close to zero for all
 * arguments */
static double
calc_resfac(double xfirst, double xlast, double xinc, double yfirst, double ylast, double yinc)
237
{
Thomas Jahns's avatar
Thomas Jahns committed
238
239
240
241
242
243
244
245
246
247
  double resfac = 1000.0;
  enum {
    nPwrOf10 = 5,
    nMultTests = 6,
  };
  static const double scaleFactors[nPwrOf10]
    = { 1000, 10000, 100000, 1000000, 10000000 };
  double vals[nMultTests] = { xfirst, xlast, xinc, yfirst, ylast, yinc };

  for (size_t j = 0; j < nPwrOf10; ++j )
248
    {
Thomas Jahns's avatar
Thomas Jahns committed
249
250
251
      double scaleBy = scaleFactors[j];
      bool fractionalScale = false;
      for (size_t i = 0; i < nMultTests; ++i )
252
        {
Thomas Jahns's avatar
Thomas Jahns committed
253
254
          fractionalScale = fractionalScale
            || fabs(vals[i]*scaleBy - round(vals[i]*scaleBy)) > FLT_EPSILON;
255
        }
Thomas Jahns's avatar
Thomas Jahns committed
256
      if ( !fractionalScale )
257
        {
Thomas Jahns's avatar
Thomas Jahns committed
258
          resfac = scaleBy;
259
260
261
262
          break;
        }
    }

Thomas Jahns's avatar
Thomas Jahns committed
263
  return (resfac);
264
265
}

266
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267
268
void iegDefGrid(int *gdb, int gridID)
{
269
  int gridtype = gridInqType(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
271
272

  if ( gridtype == GRID_GENERIC )
    {
273
274
275
276
277
278
      int xsize, ysize;

      xsize = gridInqXsize(gridID);
      ysize = gridInqYsize(gridID);

      if ( (ysize == 32  || ysize == 48 || ysize == 64 ||
279
	    ysize == 96  || ysize == 160) &&
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
	   (xsize == 2*ysize || xsize == 1) )
	{
	  gridtype = GRID_GAUSSIAN;
	  gridChangeType(gridID, gridtype);
	}
      else if ( (xsize == 1 && ysize == 1) || (xsize == 0 && ysize == 0) )
	{
	  gridtype = GRID_LONLAT;
	  gridChangeType(gridID, gridtype);
	}
      else if ( gridInqXvals(gridID, NULL) && gridInqYvals(gridID, NULL) )
	{
	  gridtype = GRID_LONLAT;
	  gridChangeType(gridID, gridtype);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
296
297
298
299
300
    }
  else if ( gridtype == GRID_CURVILINEAR )
    {
      gridtype = GRID_LONLAT;
    }

301
  if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
    {
303
304
      double xfirst = 0, xlast = 0, xinc = 0;
      double yfirst = 0, ylast = 0, yinc = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
305

306
307
      int nlon = gridInqXsize(gridID),
        nlat = gridInqYsize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308

309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
      if ( nlon == 0 )
	{
	  nlon = 1;
	}
      else
	{
	  xfirst = gridInqXval(gridID,      0);
	  xlast  = gridInqXval(gridID, nlon-1);
	  xinc   = gridInqXinc(gridID);
	}

      if ( nlat == 0 )
	{
	  nlat = 1;
	}
      else
	{
	  yfirst = gridInqYval(gridID,      0);
	  ylast  = gridInqYval(gridID, nlat-1);
	  yinc   = gridInqYinc(gridID);
	}

331
332
333
      if ( gridtype == GRID_GAUSSIAN )
	IEG_G_GridType(gdb) = 4;
      else if ( gridtype == GRID_LONLAT && gridIsRotated(gridID) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
336
337
	IEG_G_GridType(gdb) = 10;
      else
	IEG_G_GridType(gdb) = 0;

Thomas Jahns's avatar
Thomas Jahns committed
338
339
      double resfac = calc_resfac(xfirst, xlast, xinc, yfirst, ylast, yinc);
      int iresfac = (int)resfac;
340
341
342
343
      if ( iresfac == 1000 ) iresfac = 0;

      IEG_G_ResFac(gdb)   = iresfac;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
345
      IEG_G_NumLon(gdb)   = nlon;
      IEG_G_NumLat(gdb)   = nlat;
346
347
348
349
350
      IEG_G_FirstLat(gdb) = (int)lround(yfirst*resfac);
      IEG_G_LastLat(gdb)  = (int)lround(ylast*resfac);
      IEG_G_FirstLon(gdb) = (int)lround(xfirst*resfac);
      IEG_G_LastLon(gdb)  = (int)lround(xlast*resfac);
      IEG_G_LonIncr(gdb)  = (int)lround(xinc*resfac);
351
      if ( fabs(xinc*resfac - IEG_G_LonIncr(gdb)) > FLT_EPSILON )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
353
	IEG_G_LonIncr(gdb) = 0;

354
355
356
357
      if ( gridtype == GRID_GAUSSIAN )
	IEG_G_LatIncr(gdb) = nlat/2;
      else
	{
358
	  IEG_G_LatIncr(gdb) = (int)lround(yinc*resfac);
359
	  if ( fabs(yinc*resfac - IEG_G_LatIncr(gdb)) > FLT_EPSILON )
360
	    IEG_G_LatIncr(gdb) = 0;
361
362

	  if ( IEG_G_LatIncr(gdb) < 0 ) IEG_G_LatIncr(gdb) = -IEG_G_LatIncr(gdb);
363
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
364

365
      if ( IEG_G_NumLon(gdb) > 1 && IEG_G_NumLat(gdb) == 1 )
366
367
	if ( IEG_G_LonIncr(gdb) != 0 && IEG_G_LatIncr(gdb) == 0 ) IEG_G_LatIncr(gdb) = IEG_G_LonIncr(gdb);

368
      if ( IEG_G_NumLon(gdb) == 1 && IEG_G_NumLat(gdb) > 1 )
369
370
371
372
373
374
375
	if ( IEG_G_LonIncr(gdb) == 0 && IEG_G_LatIncr(gdb) != 0 ) IEG_G_LonIncr(gdb) = IEG_G_LatIncr(gdb);

      if ( IEG_G_LatIncr(gdb) == 0 || IEG_G_LonIncr(gdb) == 0 )
	IEG_G_ResFlag(gdb) = 0;
      else
	IEG_G_ResFlag(gdb) = 128;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
376
377
      if ( gridIsRotated(gridID) )
	{
378
379
	  IEG_G_LatSP(gdb) = - (int)lround(gridInqYpole(gridID) * resfac);
	  IEG_G_LonSP(gdb) =   (int)lround((gridInqXpole(gridID) + 180) * resfac);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
381
382
383
384
385
386
387
388
	  IEG_G_Size(gdb)  = 42;
	}
      else
	{
	  IEG_G_Size(gdb)  = 32;
	}
    }
  else
    {
389
      Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
391
392
393
394
    }

  IEG_G_ScanFlag(gdb) = 64;
}

395
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
397
398
399
400
401
402
403
404
405
406
void iegDefLevel(int *pdb, int *gdb, double *vct, int zaxisID, int levelID)
{
  double level;
  int ilevel, leveltype;
  static int warning = 1;
  static int vct_warning = 1;

  leveltype = zaxisInqType(zaxisID);

  if ( leveltype == ZAXIS_GENERIC )
    {
407
      Message("Changed zaxis type from %s to %s",
408
409
	      zaxisNamePtr(leveltype),
	      zaxisNamePtr(ZAXIS_PRESSURE));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
411
      leveltype = ZAXIS_PRESSURE;
      zaxisChangeType(zaxisID, leveltype);
412
      zaxisDefUnits(zaxisID, "Pa");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
414
415
416
417
418
419
420
421
422
    }

  /*  IEG_G_NumVCP(gdb) = 0; */

  switch (leveltype)
    {
    case ZAXIS_SURFACE:
      {
	IEG_P_LevelType(pdb) = IEG_LTYPE_SURFACE;
	IEG_P_Level1(pdb)    = 0;
423
	IEG_P_Level2(pdb)    = (int)(zaxisInqLevel(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
424
425
426
427
428
429
430
431
432
	break;
      }
    case ZAXIS_HYBRID:
      {
	int vctsize;

	if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_HYBRID_LAYER;
433
434
	    IEG_P_Level1(pdb)    = (int)(zaxisInqLbound(zaxisID, levelID));
	    IEG_P_Level2(pdb)    = (int)(zaxisInqUbound(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
436
437
438
439
	  }
	else
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_HYBRID;
	    IEG_P_Level1(pdb)    = 0;
440
	    IEG_P_Level2(pdb)    = (int)(zaxisInqLevel(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
442
443
444
445
	  }

	vctsize = zaxisInqVctSize(zaxisID);
	if ( vctsize == 0 && warning )
	  {
446
	    Warning("VCT missing. ( code = %d, zaxisID = %d )",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
448
449
450
451
452
453
454
		    IEG_P_Parameter(pdb), zaxisID);
	    warning = 0;
	  }
	if ( vctsize > 100 )
	  {
	    /*	    IEG_G_NumVCP(gdb) = 0; */
	    if ( vct_warning )
	      {
455
		Warning("VCT size of %d is too large (maximum is 100). Set to 0!", vctsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
457
458
459
460
461
		vct_warning = 0;
	      }
	  }
	else
	  {
	    IEG_G_Size(gdb) += (vctsize*4);
462
463
	    memcpy(vct, zaxisInqVctPtr(zaxisID), (size_t)vctsize/2*sizeof(double));
	    memcpy(vct+50, zaxisInqVctPtr(zaxisID)+vctsize/2, (size_t)vctsize/2*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
465
466
467
468
469
470
471
472
473
	  }
	break;
      }
    case ZAXIS_PRESSURE:
      {
	double dum;
	char units[128];

	level = zaxisInqLevel(zaxisID, levelID);
	if ( level < 0 )
474
	  Warning("pressure level of %f Pa is below 0.", level);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
476

	zaxisInqUnits(zaxisID, units);
477
	if ( memcmp(units, "hPa", 3) == 0 || memcmp(units, "mb",2 ) == 0 )
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
	  level = level*100;

	ilevel = (int) level;
	if ( level < 32768 && (level < 100 || modf(level/100, &dum) > 0) )
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_99;
	    IEG_P_Level1(pdb)    = 0;
	    IEG_P_Level2(pdb)    = ilevel;
	  }
	else
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_ISOBARIC;
	    IEG_P_Level1(pdb)    = 0;
	    IEG_P_Level2(pdb)    = ilevel/100;
	  }
	break;
      }
    case ZAXIS_HEIGHT:
      {
	level = zaxisInqLevel(zaxisID, levelID);

	ilevel = (int) level;
	IEG_P_LevelType(pdb) = IEG_LTYPE_HEIGHT;
	IEG_P_Level1(pdb)    = 0;
	IEG_P_Level2(pdb)    = ilevel;

	break;
      }
    case ZAXIS_ALTITUDE:
      {
	level = zaxisInqLevel(zaxisID, levelID);

	ilevel = (int) level;
	IEG_P_LevelType(pdb) = IEG_LTYPE_ALTITUDE;
	IEG_P_Level1(pdb)    = 0;
	IEG_P_Level2(pdb)    = ilevel;

	break;
      }
    case ZAXIS_DEPTH_BELOW_LAND:
      {
	if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_LANDDEPTH_LAYER;
522
523
	    IEG_P_Level1(pdb)    = (int)(zaxisInqLbound(zaxisID, levelID));
	    IEG_P_Level2(pdb)    = (int)(zaxisInqUbound(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
	  }
	else
	  {
	    level = zaxisInqLevel(zaxisID, levelID);

	    ilevel = (int) level;
	    IEG_P_LevelType(pdb) = IEG_LTYPE_LANDDEPTH;
	    IEG_P_Level1(pdb)    = 0;
	    IEG_P_Level2(pdb)    = ilevel;
	  }

	break;
      }
    case ZAXIS_DEPTH_BELOW_SEA:
      {
	level = zaxisInqLevel(zaxisID, levelID);

	ilevel = (int) level;
	IEG_P_LevelType(pdb) = IEG_LTYPE_SEADEPTH;
	IEG_P_Level1(pdb)    = 0;
	IEG_P_Level2(pdb)    = ilevel;

	break;
      }
    case ZAXIS_ISENTROPIC:
      {
	level = zaxisInqLevel(zaxisID, levelID);

	ilevel = (int) level;
	IEG_P_LevelType(pdb) = 113;
	IEG_P_Level1(pdb)    = 0;
	IEG_P_Level2(pdb)    = ilevel;

	break;
      }
    default:
      {
561
	Error("Unsupported zaxis type: %s", zaxisNamePtr(leveltype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
563
564
565
566
567
	break;
      }
    }
}


568
void iegCopyRecord(stream_t *streamptr2, stream_t *streamptr1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
569
{
570
  streamFCopyRecord(streamptr2, streamptr1, "IEG");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
573
}


574
void iegDefRecord(stream_t *streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
{
576
  int vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
579
580
  int gridID;
  int date, time;
  int datatype;
  int i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
581
  int param, pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
582
583
  int varID, levelID, tsID, zaxisID;
  int byteorder;
584
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585

Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
  vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
587
  byteorder = streamptr->byteorder;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588

Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
590
591
  varID   = streamptr->record->varID;
  levelID = streamptr->record->levelID;
  tsID    = streamptr->curTsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592
593
594

  gridID  = vlistInqVarGrid(vlistID, varID);
  zaxisID = vlistInqVarZaxis(vlistID, varID);
595

596
  iegInitMem(iegp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
597
598
599
600
  for ( i = 0; i < 37; i++ ) iegp->ipdb[i] = -1;

  iegp->byteswap = getByteswap(byteorder);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
601
  param =  vlistInqVarParam(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
  cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
603
604
  IEG_P_Parameter(iegp->ipdb) = pnum;
  if ( pdis == 255 ) IEG_P_CodeTable(iegp->ipdb) = pcat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
606
  date     = streamptr->tsteps[tsID].taxis.vdate;
  time     = streamptr->tsteps[tsID].taxis.vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
609
610
611

  iegDefTime(iegp->ipdb, date, time, vlistInqTaxis(vlistID));
  iegDefGrid(iegp->igdb, gridID);
  iegDefLevel(iegp->ipdb, iegp->igdb, iegp->vct, zaxisID, levelID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
612
  datatype = streamptr->record->prec;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
614
615
616
617

  iegp->dprec = iegDefDatatype(datatype);
}


618
void iegWriteRecord(stream_t *streamptr, const double *data)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
620
621
622
{
  int fileID;
  int i, gridsize, gridID;
  double refval;
623
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
624

Uwe Schulzweida's avatar
Uwe Schulzweida committed
625
  fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
  gridID = streamptr->record->gridID;
627

Uwe Schulzweida's avatar
Uwe Schulzweida committed
628
629
630
631
632
633
634
635
636
637
638
639
640
  gridsize = gridInqSize(gridID);

  refval = data[0];
  for ( i = 1; i < gridsize; i++ )
    if ( data[i] < refval ) refval = data[i];

  iegp->refval = refval;

  iegDefDataDP(iegp, data);

  iegWrite(fileID, iegp);
}

641
static
642
void iegAddRecord(stream_t *streamptr, int param, int *pdb, int *gdb, double *vct,
643
		  size_t recsize, off_t position, int prec)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
645
{
  int levelID = 0;
646
647
648
649
  int vlistID = streamptr->vlistID;
  int tsID    = streamptr->curTsID;
  int recID   = recordNewEntry(streamptr, tsID);
  record_t *record = &streamptr->tsteps[tsID].records[recID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650

651
  int level1, level2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
652
653
654
655
656
657
658
659
660
  if ( IEG_P_LevelType(pdb) == IEG_LTYPE_HYBRID_LAYER )
    {
      level1 = IEG_P_Level1(pdb);
      level2 = IEG_P_Level2(pdb);
    }
  else
    {
      level1 = IEG_P_Level2(pdb);
      level2 = 0;
661
      if ( IEG_P_LevelType(pdb) == 100 ) level1 *= 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
663

664
665
666
667
668
669
  record->size     = recsize;
  record->position = position;
  record->param    = param;
  record->ilevel   = level1;
  record->ilevel2  = level2;
  record->ltype    = IEG_P_LevelType(pdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
670

671
672
673
  int gridtype =
   ( IEG_G_GridType(gdb) == 0 || IEG_G_GridType(gdb) == 10 ) ? GRID_LONLAT :
    ( IEG_G_GridType(gdb) == 4 ) ? GRID_GAUSSIAN : GRID_GENERIC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674

675
676
  grid_t *grid = (grid_t *)Malloc(sizeof (*grid));
  grid_init(grid);
677
  cdiGridTypeInit(grid, gridtype, IEG_G_NumLon(gdb)*IEG_G_NumLat(gdb));
678
679
680
681
682
  grid->xsize = IEG_G_NumLon(gdb);
  grid->ysize = IEG_G_NumLat(gdb);
  grid->xinc  = 0;
  grid->yinc  = 0;
  grid->xdef  = 0;
683
684
685
686
687

  int iresfac = IEG_G_ResFac(gdb);
  if ( iresfac == 0 ) iresfac = 1000;
  double resfac = 1./(double) iresfac;

688
689
  /* if ( IEG_G_FirstLon != 0 || IEG_G_LastLon != 0 ) */
  {
690
    if ( grid->xsize > 1 )
691
692
      {
	if ( IEG_G_ResFlag(gdb) && IEG_G_LonIncr(gdb) > 0 )
693
	  grid->xinc = IEG_G_LonIncr(gdb) * resfac;
694
	else
695
	  grid->xinc = (IEG_G_LastLon(gdb) - IEG_G_FirstLon(gdb)) * resfac / (grid->xsize - 1);
696
697
698
699

	/* correct xinc if necessary */
	if ( IEG_G_FirstLon(gdb) == 0 && IEG_G_LastLon(gdb) > 354000 )
	  {
700
701
702
	    double xinc = 360. / grid->xsize;
            /* FIXME: why not use grid->xinc != xinc as condition? */
	    if ( fabs(grid->xinc-xinc) > 0.0 )
703
	      {
704
705
		grid->xinc = xinc;
		if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
706
707
708
	      }
	  }
      }
709
710
711
    grid->xfirst = IEG_G_FirstLon(gdb) * resfac;
    grid->xlast  = IEG_G_LastLon(gdb)  * resfac;
    grid->xdef   = 2;
712
  }
713
  grid->ydef  = 0;
714
715
  /* if ( IEG_G_FirstLat != 0 || IEG_G_LastLat != 0 ) */
  {
716
    if ( grid->ysize > 1 )
717
718
      {
	if ( IEG_G_ResFlag(gdb) && IEG_G_LatIncr(gdb) > 0 )
719
	  grid->yinc = IEG_G_LatIncr(gdb) * resfac;
720
	else
721
	  grid->yinc = (IEG_G_LastLat(gdb) - IEG_G_FirstLat(gdb)) * resfac / (grid->ysize - 1);
722
      }
723
724
725
    grid->yfirst = IEG_G_FirstLat(gdb) * resfac;
    grid->ylast  = IEG_G_LastLat(gdb)  * resfac;
    grid->ydef   = 2;
726
727
  }
  /*
728
729
730
731
732
733
734
735
  grid->xfirst= IEG_G_FirstLon(gdb) * resfac;
  grid->xlast = IEG_G_LastLon(gdb) * resfac;
  grid->xinc  = IEG_G_LonIncr(gdb) * resfac;
  grid->xdef  = 2;
  grid->yfirst= IEG_G_FirstLat(gdb) * resfac;
  grid->ylast = IEG_G_LastLat(gdb) * resfac;
  grid->yinc  = IEG_G_LatIncr(gdb) * resfac;
  grid->ydef  = 2;
736
  */
737
738
  grid->xvals = NULL;
  grid->yvals = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739

740
  grid->isRotated = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
742
  if ( IEG_G_GridType(gdb) == 10 )
    {
743
744
745
746
      grid->isRotated = TRUE;
      grid->ypole     = - IEG_G_LatSP(gdb) * resfac;
      grid->xpole     =   IEG_G_LonSP(gdb) * resfac - 180;
      grid->angle     = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
748
    }

749
750
751
  struct addIffNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
  int gridID = gridAdded.Id;
  if (!gridAdded.isNew) Free(grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752

753
  int leveltype = iegGetZaxisType(IEG_P_LevelType(pdb));
754

Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
756
757
  if ( leveltype == ZAXIS_HYBRID )
    {
      double tmpvct[100];
758
      size_t vctsize = (size_t)IEG_G_NumVCP(gdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759

760
761
      for (size_t i = 0; i < vctsize/2; i++ ) tmpvct[i] = vct[i];
      for (size_t i = 0; i < vctsize/2; i++ ) tmpvct[i+vctsize/2] = vct[i+50];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
763
764
765

      varDefVCT(vctsize, tmpvct);
    }

766
  int lbounds = IEG_P_LevelType(pdb) == IEG_LTYPE_HYBRID_LAYER ? 1 : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
767

768
  int datatype = iegInqDatatype(prec);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
769

770
  int varID;
771
  varAddRecord(recID, param, gridID, leveltype, lbounds, level1, level2, 0, 0,
772
773
	       datatype, &varID, &levelID, TSTEP_INSTANT, 0, 0, -1,
               NULL, NULL, NULL, NULL, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
774

775
776
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
777

Uwe Schulzweida's avatar
Uwe Schulzweida committed
778
779
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
780
781

  if ( CDI_Debug )
782
    Message("varID = %d gridID = %d levelID = %d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
783
784
785
	    varID, gridID, levelID);
}

Thomas Jahns's avatar
Thomas Jahns committed
786
#if 0
787
static
788
void iegCmpRecord(stream_t *streamptr, int tsID, int recID, off_t position, int param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
789
790
791
792
		  int level, int xsize, int ysize)
{
  int varID = 0;
  int levelID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793
  record_t *record;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
794
795

  record  = &streamptr->tsteps[tsID].records[recID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
796

Uwe Schulzweida's avatar
Uwe Schulzweida committed
797
  if ( param != (*record).param || level != (*record).ilevel )
798
    Error("inconsistent timestep");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
799
800
801
802
803
804

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
805
  streamptr->vars[varID].level[levelID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
806

Uwe Schulzweida's avatar
Uwe Schulzweida committed
807
808
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
810
  */
  if ( CDI_Debug )
811
    Message("varID = %d levelID = %d", varID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812
}
Thomas Jahns's avatar
Thomas Jahns committed
813
#endif
814

Thomas Jahns's avatar
Thomas Jahns committed
815
static void iegDateTime(int *pdb, int *date, int *time)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
817
818
819
820
821
822
823
824
825
826
827
828
{
  int ryear, rmonth, rday, rhour, rminute;

  ryear   = IEG_P_Year(pdb);

  rmonth  = IEG_P_Month(pdb);
  rday    = IEG_P_Day(pdb);

  rhour   = IEG_P_Hour(pdb);
  rminute = IEG_P_Minute(pdb);

  if ( rminute == -1 ) rminute = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
830
  *date = cdiEncodeDate(ryear, rmonth, rday);
  *time = cdiEncodeTime(rhour, rminute, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831
832
}

833
static
834
835
void iegScanTimestep1(stream_t *streamptr)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836
837
838
  int prec = 0;
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
840
  int tabnum;
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
841
  int rcode = 0, rlevel = 0, vdate = 0, vtime = 0;
842
  DateTime datetime0 = { LONG_MIN, LONG_MIN };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843
844
  int tsID;
  int varID;
845
  size_t recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
846
847
848
  off_t recpos;
  int nrecords, nrecs, recID;
  int taxisID = -1;
849
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
850
  int vlistID;
851
  IEGCOMPVAR compVar, compVar0;
852
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
853

Uwe Schulzweida's avatar
Uwe Schulzweida committed
854
  streamptr->curTsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855

856
  tsID  = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857
  taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
858
859

  if ( tsID != 0 )
860
    Error("Internal problem! tstepsNewEntry returns %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
861

Uwe Schulzweida's avatar
Uwe Schulzweida committed
862
  fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
863
864
865
866
867
868
869
870

  nrecs = 0;
  while ( TRUE )
    {
      recpos = fileGetPos(fileID);
      status = iegRead(fileID, iegp);
      if ( status != 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
871
	  streamptr->ntsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872
873
	  break;
	}
874
      recsize = (size_t)(fileGetPos(fileID) - recpos);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
875
876
877

      prec   = iegp->dprec;
      rcode  = IEG_P_Parameter(iegp->ipdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
878
      tabnum = IEG_P_CodeTable(iegp->ipdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
879
      param  = cdiEncodeParam(rcode, tabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
880

Uwe Schulzweida's avatar
Uwe Schulzweida committed
881
882
883
884
885
      if ( IEG_P_LevelType(iegp->ipdb) == IEG_LTYPE_HYBRID_LAYER )
	rlevel = IEG_P_Level1(iegp->ipdb);
      else
	rlevel = IEG_P_Level2(iegp->ipdb);

886
887
      if ( IEG_P_LevelType(iegp->ipdb) == 100 ) rlevel *= 100;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
888
889
890
891
892
893
894
895
896
      iegDateTime(iegp->ipdb, &vdate, &vtime);

      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	}
      else
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
897
	  compVar.param = param;
898
          compVar.level = rlevel;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
900
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
901
	      compVar0.param = streamptr->tsteps[0].records[recID].param;
902
903
904
	      compVar0.level = streamptr->tsteps[0].records[recID].ilevel;

	      if ( memcmp(&compVar0, &compVar, sizeof(IEGCOMPVAR)) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
905
906
	    }
	  if ( recID < nrecs ) break;
907
908
	  DateTime datetime = { .date = vdate, .time = vtime};
	  if ( datetimeCmp(datetime, datetime0) )
909
	    Warning("Inconsistent verification time for param %d level %d", param, rlevel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
910
911
912
913
914
	}

      nrecs++;

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

917
      iegAddRecord(streamptr, param, iegp->ipdb, iegp->igdb, iegp->vct, recsize, recpos, prec);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
918
919
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
920
  streamptr->rtsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921

922
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
923

924
  taxisID = taxisCreate(TAXIS_ABSOLUTE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
925
  taxis->type  = TAXIS_ABSOLUTE;
926
927
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928

Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
  vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
930
931
  vlistDefTaxis(vlistID, taxisID);

932
  vlist_check_contents(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
933

Uwe Schulzweida's avatar
Uwe Schulzweida committed
934
935
  nrecords = streamptr->tsteps[0].nallrecs;
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
936
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
937
938
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
939
	(record_t *) Realloc(streamptr->tsteps[0].records,
940
                             (size_t)nrecords * sizeof (record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
941
942
    }

943
  streamptr->tsteps[0].recIDs = (int *) Malloc((size_t)nrecords * sizeof (int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
944
  streamptr->tsteps[0].nrecs = nrecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
945
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
946
    streamptr->tsteps[0].recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
947

Uwe Schulzweida's avatar
Uwe Schulzweida committed
948
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
949
    {
950
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
951
      if ( tsID != streamptr->rtsteps )
952
	Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953

Uwe Schulzweida's avatar
Uwe Schulzweida committed
954
955
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956
957
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
958
  if ( streamptr->ntsteps == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
959
960
961
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
963
	  streamptr->ntsteps = 0;
	  for ( varID = 0; varID < streamptr->nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
964
	    {
965
	      vlistDefVarTsteptype(vlistID, varID, TSTEP_CONSTANT);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
967
968
969
970
	    }
	}
    }
}

971
static
972
int iegScanTimestep2(stream_t *streamptr)
973
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
974
975
  int status;
  int fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
976
977
  int tabnum;
  int param = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
978
979
980
  int rcode = 0, rlevel = 0, vdate = 0, vtime = 0;
  int tsID;
  int varID;
981
  size_t recsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
982
983
984
  off_t recpos = 0;
  int nrecords, nrecs, recID, rindex;
  int nextstep;
985
  taxis_t *taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
986
  int vlistID;
987
  IEGCOMPVAR compVar, compVar0;
988
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
989
990

  streamptr->curTsID = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
991

Uwe Schulzweida's avatar
Uwe Schulzweida committed
992
993
  vlistID = streamptr->vlistID;
  fileID  = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
994

Uwe Schulzweida's avatar
Uwe Schulzweida committed
995
  tsID = streamptr->rtsteps;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
996
  if ( tsID != 1 )
997
    Error("Internal problem! unexpected timestep %d", tsID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
998

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

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

1003
  cdi_create_records(streamptr, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1004

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1005
  nrecords = streamptr->tsteps[0].nallrecs;
1006
  streamptr->tsteps[1].recIDs = (int *) Malloc((size_t)nrecords * sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1007
  streamptr->tsteps[1].nrecs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1008
  for ( recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
    streamptr->tsteps[1].recIDs[recID] = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1010
1011
1012

  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1013
      varID = streamptr->tsteps[0].records[recID].varID;
1014
      streamptr->tsteps[tsID].records[recID].position =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1015
	streamptr->tsteps[0].records[recID].position;
1016
      streamptr->tsteps[tsID].records[recID].size     =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1017
	streamptr->tsteps[0].records[recID].size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1018
1019
1020
1021
1022
1023
1024
1025
    }

  for ( rindex = 0; rindex <= nrecords; rindex++ )
    {
      recpos = fileGetPos(fileID);
      status = iegRead(fileID, iegp);
      if ( status != 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1026
	  streamptr->ntsteps = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
1028
	  break;
	}
1029
      recsize = (size_t)(fileGetPos(fileID) - recpos);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1030
1031

      rcode  = IEG_P_Parameter(iegp->ipdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1032
      tabnum = IEG_P_CodeTable(iegp->ipdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1033
      param  = cdiEncodeParam(rcode, tabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1034

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1035
1036
1037
1038
1039
      if ( IEG_P_LevelType(iegp->ipdb) == IEG_LTYPE_HYBRID_LAYER )
	rlevel = IEG_P_Level1(iegp->ipdb);
      else
	rlevel = IEG_P_Level2(iegp->ipdb);

1040
1041
      if ( IEG_P_LevelType(iegp->ipdb) == 100 ) rlevel *= 100;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1042
1043
1044
1045
1046
1047
1048
1049
1050
      iegDateTime(iegp->ipdb, &vdate, &vtime);

      if ( rindex == 0 )
	{
	  taxis->type  = TAXIS_ABSOLUTE;
	  taxis->vdate = vdate;
	  taxis->vtime = vtime;
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1051
      compVar.param = param;
1052
      compVar.level = rlevel;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1053
1054
1055
      nextstep = FALSE;
      for ( recID = 0; recID < nrecords; recID++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1056
	  compVar0.param = streamptr->tsteps[tsID].records[recID].param;
1057
1058
1059
	  compVar0.level = streamptr->tsteps[tsID].records[recID].ilevel;

	  if ( memcmp(&compVar0, &compVar, sizeof(IEGCOMPVAR)) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1060
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1061
	      if ( streamptr->tsteps[tsID].records[recID].used )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1062
1063
1064
1065
1066
		{
		  nextstep = TRUE;
		}
	      else
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1067
1068
		  streamptr->tsteps[tsID].records[recID].used = TRUE;
		  streamptr->tsteps[tsID].recIDs[rindex] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1069
1070
1071
1072
1073
1074
		}
	      break;
	    }
	}
      if ( recID == nrecords )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1075
1076
	  char paramstr[32];
	  cdiParamToString(param, paramstr, sizeof(paramstr));
1077
	  Warning("param %s level %d not defined at timestep 1", paramstr, rlevel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1078
1079
1080
1081
1082
1083
	  return (CDI_EUFSTRUCT);
	}

      if ( nextstep ) break;

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1088
      compVar0.param = streamptr->tsteps[tsID].records[recID].param;
1089
      compVar0.level = streamptr->tsteps[tsID].records[recID].ilevel;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1090

1091
      if ( memcmp(&compVar0, &compVar, sizeof(IEGCOMPVAR)) != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1092
	{
1093
	  Message("tsID = %d recID = %d param = %3d new %3d  level = %3d new %3d",
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1094
		  tsID, recID,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1095
		  streamptr->tsteps[tsID].records[recID].param, param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1096
		  streamptr->tsteps[tsID].records[recID].ilevel, rlevel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1097
1098
1099
	  return (CDI_EUFSTRUCT);
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1100
      streamptr->tsteps[1].records[recID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1101
1102
1103
1104
1105
    }

  nrecs = 0;
  for ( recID = 0; recID < nrecords; recID++ )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1106
      if ( ! streamptr->tsteps[tsID].records[recID].used )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1107
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed