stream_ieg.c 35.5 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
  return datatype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52
}


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
  else                              prec = SINGLE_PRECISION;

66
  return prec;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
}

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

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

  status = iegRead(fileID, iegp);
86
  if ( status != 0 ) return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87

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

102
  return 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
}
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
  void *iegp = streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109

110
111
112
113
114
115
116
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;
  int tsID    = streamptr->curTsID;
  int vrecID  = streamptr->tsteps[tsID].curRecID;
  int recID   = streamptr->tsteps[tsID].recIDs[vrecID];
  int varID   = streamptr->tsteps[tsID].records[recID].varID;
  off_t recpos = streamptr->tsteps[tsID].records[recID].position;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
118
119

  fileSetPos(fileID, recpos, SEEK_SET);

120
  int status = iegRead(fileID, iegp);
121
122
  if ( status != 0 )
    Error("Could not read IEG record!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
124
125

  iegInqDataDP(iegp, data);

126
127
128
  double missval = vlistInqVarMissval(vlistID, varID);
  int gridID  = vlistInqVarGrid(vlistID, varID);
  int size    = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129

Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
  streamptr->numvals += size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132

  *nmiss = 0;
133
  for ( int i = 0; i < size; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
135
136
137
138
139
140
    if ( DBL_IS_EQUAL(data[i], missval) || DBL_IS_EQUAL(data[i], (float)missval) )
      {
	data[i] = missval;
	(*nmiss)++;
      }
}

141
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
144
145
146
147
148
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
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;
      }
    }

193
  return leveltype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
196
}


Thomas Jahns's avatar
Thomas Jahns committed
197
static void iegDefTime(int *pdb, int date, int time, int taxisID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
199
{
  int timetype = -1;
200
  if ( taxisID != -1 ) timetype = taxisInqType(taxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
202
203

  if ( timetype == TAXIS_ABSOLUTE || timetype == TAXIS_RELATIVE )
    {
204
      int year, month, day, hour, minute, second;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
206
      cdiDecodeDate(date, &year, &month, &day);
      cdiDecodeTime(time, &hour, &minute, &second);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223

      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
224
225
226
227
228
/* 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)
229
{
Thomas Jahns's avatar
Thomas Jahns committed
230
231
232
233
234
235
236
237
238
239
  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 )
240
    {
Thomas Jahns's avatar
Thomas Jahns committed
241
242
243
      double scaleBy = scaleFactors[j];
      bool fractionalScale = false;
      for (size_t i = 0; i < nMultTests; ++i )
244
        {
Thomas Jahns's avatar
Thomas Jahns committed
245
246
          fractionalScale = fractionalScale
            || fabs(vals[i]*scaleBy - round(vals[i]*scaleBy)) > FLT_EPSILON;
247
        }
Thomas Jahns's avatar
Thomas Jahns committed
248
      if ( !fractionalScale )
249
        {
Thomas Jahns's avatar
Thomas Jahns committed
250
          resfac = scaleBy;
251
252
253
254
          break;
        }
    }

255
  return resfac;
256
257
}

258
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
260
void iegDefGrid(int *gdb, int gridID)
{
261
  int gridtype = gridInqType(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
262
263
264

  if ( gridtype == GRID_GENERIC )
    {
265
266
      int xsize = gridInqXsize(gridID);
      int ysize = gridInqYsize(gridID);
267
268

      if ( (ysize == 32  || ysize == 48 || ysize == 64 ||
269
	    ysize == 96  || ysize == 160) &&
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
	   (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
285
286
287
288
289
290
    }
  else if ( gridtype == GRID_CURVILINEAR )
    {
      gridtype = GRID_LONLAT;
    }

291
  if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
292
    {
293
294
      double xfirst = 0, xlast = 0, xinc = 0;
      double yfirst = 0, ylast = 0, yinc = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295

296
297
      int nlon = gridInqXsize(gridID),
        nlat = gridInqYsize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298

299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
      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);
	}

321
322
323
      if ( gridtype == GRID_GAUSSIAN )
	IEG_G_GridType(gdb) = 4;
      else if ( gridtype == GRID_LONLAT && gridIsRotated(gridID) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
325
326
327
	IEG_G_GridType(gdb) = 10;
      else
	IEG_G_GridType(gdb) = 0;

Thomas Jahns's avatar
Thomas Jahns committed
328
329
      double resfac = calc_resfac(xfirst, xlast, xinc, yfirst, ylast, yinc);
      int iresfac = (int)resfac;
330
331
332
333
      if ( iresfac == 1000 ) iresfac = 0;

      IEG_G_ResFac(gdb)   = iresfac;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
      IEG_G_NumLon(gdb)   = nlon;
      IEG_G_NumLat(gdb)   = nlat;
336
337
338
339
340
      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);
341
      if ( fabs(xinc*resfac - IEG_G_LonIncr(gdb)) > FLT_EPSILON )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
343
	IEG_G_LonIncr(gdb) = 0;

344
345
346
347
      if ( gridtype == GRID_GAUSSIAN )
	IEG_G_LatIncr(gdb) = nlat/2;
      else
	{
348
	  IEG_G_LatIncr(gdb) = (int)lround(yinc*resfac);
349
	  if ( fabs(yinc*resfac - IEG_G_LatIncr(gdb)) > FLT_EPSILON )
350
	    IEG_G_LatIncr(gdb) = 0;
351
352

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

355
      if ( IEG_G_NumLon(gdb) > 1 && IEG_G_NumLat(gdb) == 1 )
356
357
	if ( IEG_G_LonIncr(gdb) != 0 && IEG_G_LatIncr(gdb) == 0 ) IEG_G_LatIncr(gdb) = IEG_G_LonIncr(gdb);

358
      if ( IEG_G_NumLon(gdb) == 1 && IEG_G_NumLat(gdb) > 1 )
359
360
361
362
363
364
365
	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
366
367
      if ( gridIsRotated(gridID) )
	{
368
369
	  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
370
371
372
373
374
375
376
377
378
	  IEG_G_Size(gdb)  = 42;
	}
      else
	{
	  IEG_G_Size(gdb)  = 32;
	}
    }
  else
    {
379
      Error("Unsupported grid type: %s", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
381
382
383
384
    }

  IEG_G_ScanFlag(gdb) = 64;
}

385
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
387
388
void iegDefLevel(int *pdb, int *gdb, double *vct, int zaxisID, int levelID)
{
  double level;
389
  int ilevel;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390

391
  int leveltype = zaxisInqType(zaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
393
394

  if ( leveltype == ZAXIS_GENERIC )
    {
395
      Message("Changed zaxis type from %s to %s",
396
397
	      zaxisNamePtr(leveltype),
	      zaxisNamePtr(ZAXIS_PRESSURE));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
399
      leveltype = ZAXIS_PRESSURE;
      zaxisChangeType(zaxisID, leveltype);
400
      zaxisDefUnits(zaxisID, "Pa");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
401
402
403
404
405
406
407
408
409
410
    }

  /*  IEG_G_NumVCP(gdb) = 0; */

  switch (leveltype)
    {
    case ZAXIS_SURFACE:
      {
	IEG_P_LevelType(pdb) = IEG_LTYPE_SURFACE;
	IEG_P_Level1(pdb)    = 0;
411
	IEG_P_Level2(pdb)    = (int)(zaxisInqLevel(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
413
414
415
416
417
418
	break;
      }
    case ZAXIS_HYBRID:
      {
	if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_HYBRID_LAYER;
419
420
	    IEG_P_Level1(pdb)    = (int)(zaxisInqLbound(zaxisID, levelID));
	    IEG_P_Level2(pdb)    = (int)(zaxisInqUbound(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
422
423
424
425
	  }
	else
	  {
	    IEG_P_LevelType(pdb) = IEG_LTYPE_HYBRID;
	    IEG_P_Level1(pdb)    = 0;
426
	    IEG_P_Level2(pdb)    = (int)(zaxisInqLevel(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
428
	  }

429
	int vctsize = zaxisInqVctSize(zaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430
431
	if ( vctsize > 100 )
	  {
432
            static bool vct_warning = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
434
435
	    /*	    IEG_G_NumVCP(gdb) = 0; */
	    if ( vct_warning )
	      {
436
		Warning("VCT size of %d is too large (maximum is 100). Set to 0!", vctsize);
437
		vct_warning = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
439
440
441
442
	      }
	  }
	else
	  {
	    IEG_G_Size(gdb) += (vctsize*4);
443
444
	    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
445
446
447
448
449
450
451
452
453
454
	  }
	break;
      }
    case ZAXIS_PRESSURE:
      {
	double dum;
	char units[128];

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

	zaxisInqUnits(zaxisID, units);
458
	if ( memcmp(units, "hPa", 3) == 0 || memcmp(units, "mb",2 ) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
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
	  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;
503
504
	    IEG_P_Level1(pdb)    = (int)(zaxisInqLbound(zaxisID, levelID));
	    IEG_P_Level2(pdb)    = (int)(zaxisInqUbound(zaxisID, levelID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
	  }
	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:
      {
542
	Error("Unsupported zaxis type: %s", zaxisNamePtr(leveltype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
544
545
546
547
548
	break;
      }
    }
}


549
void iegCopyRecord(stream_t *streamptr2, stream_t *streamptr1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
{
551
  streamFCopyRecord(streamptr2, streamptr1, "IEG");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
553
554
}


555
void iegDefRecord(stream_t *streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
{
557
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
558

559
560
  int vlistID = streamptr->vlistID;
  int byteorder = streamptr->byteorder;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561

562
563
564
  int varID   = streamptr->record->varID;
  int levelID = streamptr->record->levelID;
  int tsID    = streamptr->curTsID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
565

566
567
  int gridID  = vlistInqVarGrid(vlistID, varID);
  int zaxisID = vlistInqVarZaxis(vlistID, varID);
568

569
  iegInitMem(iegp);
570
  for ( int i = 0; i < 37; i++ ) iegp->ipdb[i] = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
573

  iegp->byteswap = getByteswap(byteorder);

574
575
  int param =  vlistInqVarParam(vlistID, varID);
  int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
  cdiDecodeParam(param, &pnum, &pcat, &pdis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
  IEG_P_Parameter(iegp->ipdb) = pnum;
  if ( pdis == 255 ) IEG_P_CodeTable(iegp->ipdb) = pcat;
579
580
  int date     = streamptr->tsteps[tsID].taxis.vdate;
  int time     = streamptr->tsteps[tsID].taxis.vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
581
582
583
584
585

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

586
  int datatype = streamptr->record->prec;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
587
588
589
590
591

  iegp->dprec = iegDefDatatype(datatype);
}


592
void iegWriteRecord(stream_t *streamptr, const double *data)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
{
594
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
595

596
597
  int fileID = streamptr->fileID;
  int gridID = streamptr->record->gridID;
598

599
  int gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
600

601
602
  double refval = data[0];
  for ( int i = 1; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
603
604
605
606
607
608
609
610
611
    if ( data[i] < refval ) refval = data[i];

  iegp->refval = refval;

  iegDefDataDP(iegp, data);

  iegWrite(fileID, iegp);
}

612
static
613
void iegAddRecord(stream_t *streamptr, int param, int *pdb, int *gdb, double *vct,
614
		  size_t recsize, off_t position, int prec)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
{
616
617
618
619
  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
620

621
  int level1, level2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
623
624
625
626
627
628
629
630
  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;
631
      if ( IEG_P_LevelType(pdb) == 100 ) level1 *= 100;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
632
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633

634
635
636
637
638
639
  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
640

641
642
643
  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
644

645
646
  grid_t *grid = (grid_t *)Malloc(sizeof (*grid));
  grid_init(grid);
647
  cdiGridTypeInit(grid, gridtype, IEG_G_NumLon(gdb)*IEG_G_NumLat(gdb));
648
649
650
651
  int xsize = IEG_G_NumLon(gdb);
  int ysize = IEG_G_NumLat(gdb);
  grid->x.size = xsize;
  grid->y.size = ysize;
652
653
654
  grid->xinc  = 0;
  grid->yinc  = 0;
  grid->xdef  = 0;
655
656
657
658
659

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

660
661
  /* if ( IEG_G_FirstLon != 0 || IEG_G_LastLon != 0 ) */
  {
662
    if ( xsize > 1 )
663
664
      {
	if ( IEG_G_ResFlag(gdb) && IEG_G_LonIncr(gdb) > 0 )
665
	  grid->xinc = IEG_G_LonIncr(gdb) * resfac;
666
	else
667
	  grid->xinc = (IEG_G_LastLon(gdb) - IEG_G_FirstLon(gdb)) * resfac / (xsize - 1);
668
669
670
671

	/* correct xinc if necessary */
	if ( IEG_G_FirstLon(gdb) == 0 && IEG_G_LastLon(gdb) > 354000 )
	  {
672
	    double xinc = 360. / xsize;
673
674
            /* FIXME: why not use grid->xinc != xinc as condition? */
	    if ( fabs(grid->xinc-xinc) > 0.0 )
675
	      {
676
677
		grid->xinc = xinc;
		if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
678
679
680
	      }
	  }
      }
681
682
683
    grid->xfirst = IEG_G_FirstLon(gdb) * resfac;
    grid->xlast  = IEG_G_LastLon(gdb)  * resfac;
    grid->xdef   = 2;
684
  }
685
  grid->ydef  = 0;
686
687
  /* if ( IEG_G_FirstLat != 0 || IEG_G_LastLat != 0 ) */
  {
688
    if ( ysize > 1 )
689
690
      {
	if ( IEG_G_ResFlag(gdb) && IEG_G_LatIncr(gdb) > 0 )
691
	  grid->yinc = IEG_G_LatIncr(gdb) * resfac;
692
	else
693
	  grid->yinc = (IEG_G_LastLat(gdb) - IEG_G_FirstLat(gdb)) * resfac / (ysize - 1);
694
      }
695
696
697
    grid->yfirst = IEG_G_FirstLat(gdb) * resfac;
    grid->ylast  = IEG_G_LastLat(gdb)  * resfac;
    grid->ydef   = 2;
698
699
  }
  /*
700
701
702
703
704
705
706
707
  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;
708
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
709

710
  grid->isRotated = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
711
712
  if ( IEG_G_GridType(gdb) == 10 )
    {
713
714
715
716
      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
717
718
    }

719
720
721
  struct addIffNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
  int gridID = gridAdded.Id;
  if (!gridAdded.isNew) Free(grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
722

723
  int leveltype = iegGetZaxisType(IEG_P_LevelType(pdb));
724

Uwe Schulzweida's avatar
Uwe Schulzweida committed
725
726
727
  if ( leveltype == ZAXIS_HYBRID )
    {
      double tmpvct[100];
728
      size_t vctsize = (size_t)IEG_G_NumVCP(gdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729

730
731
      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
732
733
734
735

      varDefVCT(vctsize, tmpvct);
    }

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

738
  int datatype = iegInqDatatype(prec);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739

740
  int varID;
741
  int levelID = 0;
742
  varAddRecord(recID, param, gridID, leveltype, lbounds, level1, level2, 0, 0,
743
744
	       datatype, &varID, &levelID, TSTEP_INSTANT, 0, 0, -1,
               NULL, NULL, NULL, NULL, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745

746
747
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
748

Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
751
752

  if ( CDI_Debug )
753
    Message("varID = %d gridID = %d levelID = %d", varID, gridID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754
755
}

Thomas Jahns's avatar
Thomas Jahns committed
756
#if 0
757
static
758
void iegCmpRecord(stream_t *streamptr, int tsID, int recID, off_t position, int param,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759
760
761
762
		  int level, int xsize, int ysize)
{
  int varID = 0;
  int levelID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
763

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
766
  if ( param != (*record).param || level != (*record).ilevel )
767
    Error("inconsistent timestep");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
768
769
770
771
772
773

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
776
777
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778
779
  */
  if ( CDI_Debug )
780
    Message("varID = %d levelID = %d", varID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
781
}
Thomas Jahns's avatar
Thomas Jahns committed
782
#endif
783

Thomas Jahns's avatar
Thomas Jahns committed
784
static void iegDateTime(int *pdb, int *date, int *time)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
{
786
  int ryear   = IEG_P_Year(pdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
787

788
789
  int rmonth  = IEG_P_Month(pdb);
  int rday    = IEG_P_Day(pdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
790

791
792
  int rhour   = IEG_P_Hour(pdb);
  int rminute = IEG_P_Minute(pdb);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793
794
795

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
796
797
  *date = cdiEncodeDate(ryear, rmonth, rday);
  *time = cdiEncodeTime(rhour, rminute, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
798
799
}

800
static
801
802
void iegScanTimestep1(stream_t *streamptr)
{
803
  DateTime datetime0 = { LONG_MIN, LONG_MIN };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
804
  off_t recpos;
805
  int nrecs;
806
  IEGCOMPVAR compVar, compVar0;
807
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
808

Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
  streamptr->curTsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810

811
812
  int tsID  = tstepsNewEntry(streamptr);
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
813
814

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

817
  int fileID = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
818
819
820
821
822

  nrecs = 0;
  while ( TRUE )
    {
      recpos = fileGetPos(fileID);
823
      int status = iegRead(fileID, iegp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824
825
      if ( status != 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
	  streamptr->ntsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
827
828
	  break;
	}
829
      size_t recsize = (size_t)(fileGetPos(fileID) - recpos);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830

831
832
833
834
      int prec   = iegp->dprec;
      int rcode  = IEG_P_Parameter(iegp->ipdb);
      int tabnum = IEG_P_CodeTable(iegp->ipdb);
      int param  = cdiEncodeParam(rcode, tabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
835

836
      int rlevel = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
837
838
839
840
841
      if ( IEG_P_LevelType(iegp->ipdb) == IEG_LTYPE_HYBRID_LAYER )
	rlevel = IEG_P_Level1(iegp->ipdb);
      else
	rlevel = IEG_P_Level2(iegp->ipdb);

842
843
      if ( IEG_P_LevelType(iegp->ipdb) == 100 ) rlevel *= 100;

844
      int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
845
846
847
848
849
850
851
852
853
      iegDateTime(iegp->ipdb, &vdate, &vtime);

      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	}
      else
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
854
	  compVar.param = param;
855
          compVar.level = rlevel;
856
          int recID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857
858
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
859
	      compVar0.param = streamptr->tsteps[0].records[recID].param;
860
861
862
	      compVar0.level = streamptr->tsteps[0].records[recID].ilevel;

	      if ( memcmp(&compVar0, &compVar, sizeof(IEGCOMPVAR)) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
863
864
	    }
	  if ( recID < nrecs ) break;
865
866
	  DateTime datetime = { .date = vdate, .time = vtime};
	  if ( datetimeCmp(datetime, datetime0) )
867
	    Warning("Inconsistent verification time for param %d level %d", param, rlevel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
868
869
870
871
872
	}

      nrecs++;

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

875
      iegAddRecord(streamptr, param, iegp->ipdb, iegp->igdb, iegp->vct, recsize, recpos, prec);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
876
877
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
878
  streamptr->rtsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
879

880
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
881

882
  int taxisID = taxisCreate(TAXIS_ABSOLUTE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
883
  taxis->type  = TAXIS_ABSOLUTE;
884
885
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
886

887
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888
889
  vlistDefTaxis(vlistID, taxisID);

890
  vlist_check_contents(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
891

892
  int nrecords = streamptr->tsteps[0].nallrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
893
  if ( nrecords < streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
894
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
895
896
      streamptr->tsteps[0].recordSize = nrecords;
      streamptr->tsteps[0].records =
897
	(record_t *) Realloc(streamptr->tsteps[0].records,
898
                             (size_t)nrecords * sizeof (record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
899
900
    }

901
  streamptr->tsteps[0].recIDs = (int *) Malloc((size_t)nrecords * sizeof (int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
902
  streamptr->tsteps[0].nrecs = nrecords;
903
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
904
    streamptr->tsteps[0].recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
905

Uwe Schulzweida's avatar
Uwe Schulzweida committed
906
  if ( streamptr->ntsteps == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
907
    {
908
      tsID = tstepsNewEntry(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
909
      if ( tsID != streamptr->rtsteps )
910
	Error("Internal error. tsID = %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
911

Uwe Schulzweida's avatar
Uwe Schulzweida committed
912
913
      streamptr->tsteps[tsID-1].next   = TRUE;
      streamptr->tsteps[tsID].position = recpos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
914
915
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
916
  if ( streamptr->ntsteps == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
917
918
919
    {
      if ( taxis->vdate == 0 && taxis->vtime == 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
920
	  streamptr->ntsteps = 0;
921
	  for ( int varID = 0; varID < streamptr->nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
922
	    {
923
	      vlistDefVarTsteptype(vlistID, varID, TSTEP_CONSTANT);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
924
925
926
927
928
	    }
	}
    }
}

929
static
930
int iegScanTimestep2(stream_t *streamptr)
931
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
932
  off_t recpos = 0;
933
  IEGCOMPVAR compVar, compVar0;
934
  iegrec_t *iegp = (iegrec_t*) streamptr->record->exsep;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
935
936

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

938
939
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
940

941
  int tsID = streamptr->rtsteps;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
  if ( tsID != 1 )
943
    Error("Internal problem! unexpected timestep %d", tsID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
944

945
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
946

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

949
  cdi_create_records(streamptr, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
950

951
  int nrecords = streamptr->tsteps[0].nallrecs;
952
  streamptr->tsteps[1].recIDs = (int *) Malloc((size_t)nrecords * sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
953
  streamptr->tsteps[1].nrecs = 0;
954
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
    streamptr->tsteps[1].recIDs[recID] = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
956

957
  for ( int recID = 0; recID < nrecords; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
958
    {
959
      streamptr->tsteps[tsID].records[recID].position =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
960
	streamptr->tsteps[0].records[recID].position;
961
      streamptr->tsteps[tsID].records[recID].size     =
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
	streamptr->tsteps[0].records[recID].size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
963
964
    }

965
  for ( int rindex = 0; rindex <= nrecords; rindex++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
966
967
    {
      recpos = fileGetPos(fileID);
968
      int status = iegRead(fileID, iegp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
969
970
      if ( status != 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
971
	  streamptr->ntsteps = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
973
	  break;
	}
974
      size_t recsize = (size_t)(fileGetPos(fileID) - recpos);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
975

976
977
978
      int rcode  = IEG_P_Parameter(iegp->ipdb);
      int tabnum = IEG_P_CodeTable(iegp->ipdb);
      int param  = cdiEncodeParam(rcode, tabnum, 255);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
979

980
      int rlevel = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
981
982
983
984
985
      if ( IEG_P_LevelType(iegp->ipdb) == IEG_LTYPE_HYBRID_LAYER )
	rlevel = IEG_P_Level1(iegp->ipdb);
      else
	rlevel = IEG_P_Level2(iegp->ipdb);

986
987
      if ( IEG_P_LevelType(iegp->ipdb) == 100 ) rlevel *= 100;

988
      int vdate = 0, vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
989
990
991
992
993
994
995
996
997
      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
998
      compVar.param = param;
999
      compVar.level = rlevel;
1000
1001
      bool nextstep = false;
      int recID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1002
1003
      for ( recID = 0; recID < nrecords; recID++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1004
	  compVar0.param = streamptr->tsteps[tsID].records[recID].param;
1005
1006
1007
	  compVar0.level = streamptr->tsteps[tsID].records[recID].ilevel;

	  if ( memcmp(&compVar0, &compVar, sizeof(IEGCOMPVAR)) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1008
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
	      if ( streamptr->tsteps[tsID].records[recID].used )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1010
		{
1011
		  nextstep = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1012
1013
1014
		}
	      else
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1015
1016
		  streamptr->tsteps[tsID].records[recID].used = TRUE;
		  streamptr->tsteps[tsID].recIDs[rindex] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1017
1018
1019
1020
1021
1022
		}
	      break;
	    }
	}
      if ( recID == nrecords )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1023
1024
	  char paramstr[32];
	  cdiParamToString(param, paramstr, sizeof(paramstr));
1025
	  Warning("param %s level %d not defined at timestep 1", paramstr, rlevel);
1026
	  return CDI_EUFSTRUCT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
1028
1029
1030
1031
	}

      if ( nextstep ) break;

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

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

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