varscan.c 38 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
#ifdef HAVE_CONFIG_H
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3
4
5
6
#endif


#include "cdi.h"
7
#include "cdi_int.h"
8
#include "cdi_uuid.h"
9
#include "cdi_key.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
#include "dmemory.h"
11
#include "resource_handle.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
12
13
#include "varscan.h"
#include "vlist.h"
14
#include "zaxis.h"
15
#include "subtype.h"
16
17


Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20
static size_t Vctsize = 0;
static double *Vct = NULL;

21
22
static int numberOfVerticalLevels = 0;
static int numberOfVerticalGrid = 0;
23
static unsigned char uuidVGrid[CDI_UUID_SIZE];
24

25

Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
27
28
29
30
31
32
typedef struct
{
  int      level1;
  int      level2;
  int      recID;
  int      lindex;
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
leveltable_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
34

35

Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
typedef struct
{
38
  int           subtypeIndex; //  corresponding tile in subtype_t structure (subtype->self)
39
  unsigned      nlevels;
40
  int           levelTableSize;
41
42
43
44
45
46
  leveltable_t* levelTable;
} subtypetable_t;


typedef struct
{
47
  int            varID;
48
49
50
  int            param;
  int            prec;
  int            tsteptype;
51
  VarScanKeys    scanKeys;
52
53
  int            gridID;
  int            zaxistype;
54
55
  int            ltype1;        // GRIB first level type
  int            ltype2;        // GRIB second level type
56
57
58
59
60
61
62
  int            lbounds;
  int            level_sf;
  int            level_unit;
  int            zaxisID;

  int            nsubtypes_alloc;
  int            nsubtypes;
63
  subtypetable_t *recordTable;   // ~ two-dimensional record list [nsubtypes_alloc][levelTableSize]
64
65
66
67
68
69

  int            instID;
  int            modelID;
  int            tableID;
  int            comptype;       // compression type
  int            complevel;      // compression level
70
  short          timave;
71
  bool           lmissval;
72
73
74
75
  double         missval;
  char          *name;

  /* meta-data for specification of tiles (currently only GRIB-API: */
76
  subtype_t     *tiles;
77

78
  cdi_keys_t     keys;
79

80
81
82
  int                 opt_grib_nentries;        // current no. key-value pairs
  int                 opt_grib_kvpair_size;     // current allocated size
  opt_key_val_pair_t *opt_grib_kvpair;          // (optional) list of keyword/value pairs
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
vartable_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85

86

87
static vartable_t *vartable;
88
static unsigned varTablesize = 0;
89
static unsigned varTableUsed = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

91
92
static
void paramInitEntry(unsigned varID, int param)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
{
Thomas Jahns's avatar
Thomas Jahns committed
94
  vartable[varID].varID          = (int)varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
  vartable[varID].param          = param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
  vartable[varID].prec           = 0;
97
  vartable[varID].tsteptype      = TSTEP_INSTANT;
98
  varScanKeysInit(&vartable[varID].scanKeys);
99
  vartable[varID].timave         = 0;
100
  vartable[varID].gridID         = CDI_UNDEFID;
101
  vartable[varID].zaxistype      = 0;
102
103
  vartable[varID].ltype1         = 0;
  vartable[varID].ltype2         = -1;
104
105
  vartable[varID].lbounds        = 0;
  vartable[varID].level_sf       = 0;
106
  vartable[varID].level_unit     = 0;
107
108
109
  vartable[varID].recordTable    = NULL;
  vartable[varID].nsubtypes_alloc= 0;
  vartable[varID].nsubtypes      = 0;
110
111
112
  vartable[varID].instID         = CDI_UNDEFID;
  vartable[varID].modelID        = CDI_UNDEFID;
  vartable[varID].tableID        = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
  cdiInitKeys(&vartable[varID].keys);
114
  vartable[varID].comptype       = CDI_COMPRESS_NONE;
115
  vartable[varID].complevel      = 1;
116
  vartable[varID].lmissval       = false;
117
  vartable[varID].missval        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
  vartable[varID].name           = NULL;
119
  vartable[varID].tiles          = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
}

122
// Test if a variable specified by the given meta-data has already been registered in "vartable".
123
static unsigned
124
varGetEntry(int param, int gridID, int zaxistype, int ltype1, int tsteptype, const char *name, const VarScanKeys *scanKeys, const var_tile_t *tiles)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
{
126
  for ( unsigned varID = 0; varID < varTablesize; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
    {
128
      // testing for "param" implicitly checks if we are beyond the current vartable size:
129
      if ( vartable[varID].param == param )
130
        {
131
132
133
          const int no_of_tiles = tiles ? tiles->numberOfTiles : -1;
          const int vt_no_of_tiles = vartable[varID].tiles ?
            subtypeGetGlobalDataP(vartable[varID].tiles, SUBTYPE_ATT_NUMBER_OF_TILES) : -1;
134
135
136
          if ( (vartable[varID].zaxistype  == zaxistype)               &&
               (vartable[varID].ltype1     == ltype1   )               &&
               (vartable[varID].tsteptype  == tsteptype)               &&
137
               (scanKeys == NULL || varScanKeysIsEqual(&vartable[varID].scanKeys, scanKeys))  &&
138
               (vartable[varID].gridID     == gridID   )               &&
139
               (vt_no_of_tiles == no_of_tiles) )
140
            {
141
142
              if ( name && name[0] && vartable[varID].name && vartable[varID].name[0] )
                {
143
                  if ( strcmp(name, vartable[varID].name) == 0 ) return varID;
144
145
146
                }
              else
                {
147
                  return varID;
148
                }
149
150
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
152
    }

153
  return (unsigned)-1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
155
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
158
void varFree(void)
{
159
160
  if ( CDI_Debug ) Message("call to varFree");

161
  for ( size_t varID = 0; varID < varTableUsed; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
    {
163
164
165
      if ( vartable[varID].recordTable )
        {
          for (int isub=0; isub<vartable[varID].nsubtypes_alloc; isub++)
166
167
            Free(vartable[varID].recordTable[isub].levelTable);
          Free(vartable[varID].recordTable);
168
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169

170
      if ( vartable[varID].name )     Free(vartable[varID].name);
171
      if ( vartable[varID].tiles )    subtypeDestroyPtr(vartable[varID].tiles);
172

173
174
175
      cdi_keys_t *keysp = &(vartable[varID].keys);
      cdiDeleteVarKeys(keysp);

176
      if ( vartable[varID].opt_grib_kvpair )
177
        {
178
179
180
181
182
          for (int i=0; i<vartable[varID].opt_grib_nentries; i++)
            {
              if ( vartable[varID].opt_grib_kvpair[i].keyword )
                Free(vartable[varID].opt_grib_kvpair[i].keyword);
            }
183
          Free(vartable[varID].opt_grib_kvpair);
184
        }
185
186
187
      vartable[varID].opt_grib_nentries    = 0;
      vartable[varID].opt_grib_kvpair_size = 0;
      vartable[varID].opt_grib_kvpair      = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
189
    }

190
  if ( vartable ) Free(vartable);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
191
192
  vartable = NULL;
  varTablesize = 0;
193
  varTableUsed = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194

195
  if ( Vct ) Free(Vct);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
197
198
199
  Vct = NULL;
  Vctsize = 0;
}

200
// Search for a tile subtype with subtypeIndex == tile_index.
201
static int tileGetEntry(unsigned varID, int tile_index)
202
203
204
205
206
207
208
209
210
{
  for (int isub=0; isub<vartable[varID].nsubtypes; isub++)
    if (vartable[varID].recordTable[isub].subtypeIndex == tile_index)
      return isub;
  return CDI_UNDEFID;
}


/* Resizes vartable:recordTable data structure, if necessary. */
211
static int tileNewEntry(int varID)
212
213
214
215
216
{
  int tileID = 0;
  if (vartable[varID].nsubtypes_alloc == 0)
    {
      /* create table for the first time. */
217
      vartable[varID].nsubtypes_alloc = 2;
218
219
      vartable[varID].nsubtypes       = 0;
      vartable[varID].recordTable     =
220
        (subtypetable_t *) Malloc((size_t)vartable[varID].nsubtypes_alloc * sizeof (subtypetable_t));
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
      if( vartable[varID].recordTable == NULL )
        SysError("Allocation of leveltable failed!");

      for (int isub = 0; isub<vartable[varID].nsubtypes_alloc; isub++) {
	vartable[varID].recordTable[isub].levelTable     = NULL;
        vartable[varID].recordTable[isub].levelTableSize = 0;
        vartable[varID].recordTable[isub].nlevels        = 0;
        vartable[varID].recordTable[isub].subtypeIndex   = CDI_UNDEFID;
      }
    }
  else
    {
      /* data structure large enough; find a free entry. */
      while(tileID <  vartable[varID].nsubtypes_alloc)
	{
	  if (vartable[varID].recordTable[tileID].levelTable == NULL) break;
	  tileID++;
	}
    }

  /* If the table overflows, double its size. */
  if (tileID == vartable[varID].nsubtypes_alloc)
    {
      tileID = vartable[varID].nsubtypes_alloc;
      vartable[varID].nsubtypes_alloc *= 2;
      vartable[varID].recordTable   =
247
        (subtypetable_t *) Realloc(vartable[varID].recordTable,
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
                                   (size_t)vartable[varID].nsubtypes_alloc * sizeof (subtypetable_t));
      if (vartable[varID].recordTable == NULL)
        SysError("Reallocation of leveltable failed");
      for(int isub=tileID; isub<vartable[varID].nsubtypes_alloc; isub++) {
	vartable[varID].recordTable[isub].levelTable     = NULL;
        vartable[varID].recordTable[isub].levelTableSize = 0;
        vartable[varID].recordTable[isub].nlevels        = 0;
        vartable[varID].recordTable[isub].subtypeIndex   = CDI_UNDEFID;
      }
    }

  return tileID;
}


263
static int levelNewEntry(unsigned varID, int level1, int level2, int tileID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
{
  int levelID = 0;
Thomas Jahns's avatar
Thomas Jahns committed
266
267
  int levelTableSize = vartable[varID].recordTable[tileID].levelTableSize;
  leveltable_t *levelTable = vartable[varID].recordTable[tileID].levelTable;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
268

269
  // Look for a free slot in levelTable. (Create the table the first time through).
Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
271
272
  if ( ! levelTableSize )
    {
      levelTableSize = 2;
273
274
      levelTable = (leveltable_t *) Malloc((size_t)levelTableSize * sizeof (leveltable_t));
      for ( int i = 0; i < levelTableSize; i++ ) levelTable[i].recID = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
276
277
    }
  else
    {
278
      while( levelID < levelTableSize && levelTable[levelID].recID != CDI_UNDEFID )
Thomas Jahns's avatar
Thomas Jahns committed
279
        ++levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280
    }
281
282

  // If the table overflows, double its size.
283
  if ( levelID == levelTableSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
    {
Thomas Jahns's avatar
Thomas Jahns committed
285
286
      levelTable = (leveltable_t *) Realloc(levelTable,
                                            (size_t)(levelTableSize *= 2)
287
                                            * sizeof (leveltable_t));
Thomas Jahns's avatar
Thomas Jahns committed
288
      for( int i = levelID; i < levelTableSize; i++ )
289
        levelTable[i].recID = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
291
    }

292
293
294
  levelTable[levelID].level1 = level1;
  levelTable[levelID].level2 = level2;
  levelTable[levelID].lindex = levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295

296
297
298
  vartable[varID].recordTable[tileID].nlevels        = (unsigned)levelID+1;
  vartable[varID].recordTable[tileID].levelTableSize = levelTableSize;
  vartable[varID].recordTable[tileID].levelTable     = levelTable;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299

300
  return levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
302
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
303
#define  UNDEF_PARAM  -4711
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304

305
306
static unsigned
paramNewEntry(int param)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
307
{
308
  unsigned varID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309

310
  // Look for a free slot in vartable. (Create the table the first time through).
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
312
313
  if ( ! varTablesize )
    {
      varTablesize = 2;
314
      vartable = (vartable_t *) Malloc((size_t)varTablesize * sizeof (vartable_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
      if( vartable == NULL )
	{
317
318
          Message("varTablesize = %d", varTablesize);
	  SysError("Allocation of vartable failed");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
320
	}

321
      for( unsigned i = 0; i < varTablesize; i++ )
322
323
        {
          vartable[i].param = UNDEF_PARAM;
324
          vartable[i].opt_grib_kvpair      = NULL;
325
          vartable[i].opt_grib_kvpair_size = 0;
326
          vartable[i].opt_grib_nentries    = 0;
327
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
329
330
331
332
    }
  else
    {
      while( varID < varTablesize )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
	  if ( vartable[varID].param == UNDEF_PARAM ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
336
	  varID++;
	}
    }
337
338

  // If the table overflows, double its size.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
  if ( varID == varTablesize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
    {
Thomas Jahns's avatar
Thomas Jahns committed
341
      vartable = (vartable_t *) Realloc(vartable, (size_t)(varTablesize *= 2)
342
                                        * sizeof (vartable_t));
Thomas Jahns's avatar
Thomas Jahns committed
343
      for ( size_t i = varID; i < varTablesize; i++ )
344
345
        {
          vartable[i].param = UNDEF_PARAM;
346
          vartable[i].opt_grib_kvpair      = NULL;
347
          vartable[i].opt_grib_kvpair_size = 0;
348
          vartable[i].opt_grib_nentries    = 0;
349
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
351
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
  paramInitEntry(varID, param);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353

354
  return varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355
356
}

357
// Append tile set to a subtype. Return index of the new tile (i.e. the "entry->self" value).
358
359
static
int varInsertTileSubtype(vartable_t *vptr, const var_tile_t *tiles)
360
{
361
  if ( tiles == NULL ) return 0;
362

363
  // first, generate a subtype based on the info in "tiles".
364
365
  subtype_t *subtype_ptr;
  subtypeAllocate(&subtype_ptr, SUBTYPE_TILES);
366
367
368
  subtypeDefGlobalDataP(subtype_ptr, SUBTYPE_ATT_TOTALNO_OF_TILEATTR_PAIRS, tiles->totalno_of_tileattr_pairs);
  subtypeDefGlobalDataP(subtype_ptr, SUBTYPE_ATT_TILE_CLASSIFICATION      , tiles->tileClassification);
  subtypeDefGlobalDataP(subtype_ptr, SUBTYPE_ATT_NUMBER_OF_TILES          , tiles->numberOfTiles);
369

370
  // Here, we create a tile set for comparison that contains only one tile/attribute pair (based on "tiles").
371
  struct subtype_entry_t *entry = subtypeEntryInsert(subtype_ptr);
372
373
374
  subtypeDefEntryDataP(entry, SUBTYPE_ATT_NUMBER_OF_ATTR, tiles->numberOfAttributes);
  subtypeDefEntryDataP(entry, SUBTYPE_ATT_TILEINDEX,      tiles->tileindex);
  subtypeDefEntryDataP(entry, SUBTYPE_ATT_TILEATTRIBUTE,  tiles->attribute);
375
376
377
378
379
380
381
382
383
384

  if (vptr->tiles == NULL) {
    vptr->tiles = subtype_ptr;
    return 0;
  }
  else {
    tilesetInsertP(vptr->tiles, subtype_ptr);
    subtypeDestroyPtr(subtype_ptr);
    return vptr->tiles->nentries - 1;
  }
385

386
387
388
389
  return CDI_UNDEFID;
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
void varAddRecord(int recID, int param, int gridID, int zaxistype, int lbounds,
391
		  int level1, int level2, int level_sf, int level_unit, int prec,
392
		  int *pvarID, int *plevelID, int tsteptype, int numavg, int ltype1, int ltype2,
393
		  const char *name, const VarScanKeys *scanKeys, const var_tile_t *tiles, int *tile_index)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
{
395
  unsigned varID = (cdiSplitLtype105 != 1 || zaxistype != ZAXIS_HEIGHT) ?
396
    varGetEntry(param, gridID, zaxistype, ltype1, tsteptype, name, scanKeys, tiles) : (unsigned) CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397

398
  if ( varID == (unsigned) CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
399
    {
400
      varTableUsed++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
401
      varID = paramNewEntry(param);
402
403
      vartable[varID].gridID     = gridID;
      vartable[varID].zaxistype  = zaxistype;
404
405
      vartable[varID].ltype1     = ltype1;
      vartable[varID].ltype2     = ltype2;
406
407
      vartable[varID].lbounds    = lbounds;
      vartable[varID].level_sf   = level_sf;
408
      vartable[varID].level_unit = level_unit;
409
      vartable[varID].tsteptype  = tsteptype;
410
      if (scanKeys) vartable[varID].scanKeys = *scanKeys;
411

412
      if ( numavg ) vartable[varID].timave = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413

Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
      if ( name && name[0] )         vartable[varID].name     = strdup(name);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
416
417
    }
  else
    {
418
419
420
      char paramstr[32];
      cdiParamToString(param, paramstr, sizeof(paramstr));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
422
      if ( vartable[varID].gridID != gridID )
	{
423
	  Message("param = %s gridID = %d", paramstr, gridID);
424
	  Error("horizontal grid must not change for same parameter!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
	}
426
      if ( vartable[varID].zaxistype != zaxistype )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
	{
428
	  Message("param = %s zaxistype = %d", paramstr, zaxistype);
429
	  Error("zaxistype must not change for same parameter!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430
431
432
	}
    }

433
434
  if ( prec > vartable[varID].prec ) vartable[varID].prec = prec;

435
436
  // append current tile to tile subtype info.
  const int this_tile = varInsertTileSubtype(&vartable[varID], tiles);
437
438
  int tileID = tileGetEntry(varID, this_tile);
  if ( tile_index ) (*tile_index) = this_tile;
439
440
441
442
443
444
  if ( tileID == CDI_UNDEFID )
    {
      tileID = tileNewEntry((int)varID);
      vartable[varID].recordTable[tileID].subtypeIndex = this_tile;
      vartable[varID].nsubtypes++;
    }
445

446
447
  // append current level to level table info
  const int levelID = levelNewEntry(varID, level1, level2, tileID);
448
  if (CDI_Debug)
449
    Message("vartable[%d].recordTable[%d].levelTable[%d].recID = %d; level1,2=%d,%d",
450
            varID, tileID, levelID, recID, level1, level2);
451
  vartable[varID].recordTable[tileID].levelTable[levelID].recID = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452

453
  *pvarID   = (int) varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
455
  *plevelID = levelID;
}
456
457


Uwe Schulzweida's avatar
Uwe Schulzweida committed
458
/*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
459
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
461
462
463
464
465
466
int dblcmp(const void *s1, const void *s2)
{
  int cmp = 0;

  if      ( *((double *) s1) < *((double *) s2) ) cmp = -1;
  else if ( *((double *) s1) > *((double *) s2) ) cmp =  1;

467
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470
static
471
int cmpLevelTable(const void* s1, const void* s2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
473
{
  int cmp = 0;
474
475
  const leveltable_t *x = (const leveltable_t*) s1;
  const leveltable_t *y = (const leveltable_t*) s2;
476
  // printf("%g %g  %d %d\n", x->leve11, y->level1, x, y);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
478
479
  if      ( x->level1 < y->level1 ) cmp = -1;
  else if ( x->level1 > y->level1 ) cmp =  1;

480
  return cmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
482
}

483
484
485
486
static
int cmpLevelTableInv(const void* s1, const void* s2)
{
  int cmp = 0;
487
488
  const leveltable_t *x = (const leveltable_t*) s1;
  const leveltable_t *y = (const leveltable_t*) s2;
489
  // printf("%g %g  %d %d\n", x->leve11, y->level1, x, y);
490
491
492
  if      ( x->level1 < y->level1 ) cmp =  1;
  else if ( x->level1 > y->level1 ) cmp = -1;

493
  return cmp;
494
495
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
496

497
498
void varCopyKeys(int vlistID, int varID)
{
499
  vlist_t *vlistptr = vlist_to_pointer(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
  cdiInitKeys(&vlistptr->vars[varID].keys);
501
  cdiCopyVarKeys(&vartable[varID].keys, &vlistptr->vars[varID].keys);
502
503
504
}


505
506
507
508
509
510
511
512
513
514
515
516
517
518
struct cdi_generate_varinfo
{
  int        varid;
  const char *name;
};

static
int cdi_generate_cmp_varname(const void *s1, const void *s2)
{
  const struct cdi_generate_varinfo *x = (const struct cdi_generate_varinfo *)s1,
                                    *y = (const struct cdi_generate_varinfo *)s2;
  return strcmp(x->name, y->name);
}

519
void cdi_generate_vars(stream_t *streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
{
521
  const int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
522

523
524
  int *varids = (int *) Malloc(varTableUsed*sizeof(int));
  for ( size_t varID = 0; varID < varTableUsed; varID++ ) varids[varID] = (int)varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
  if ( streamptr->sortname )
    {
      size_t varID;
      for ( varID = 0; varID < varTableUsed; varID++ )
        if (!vartable[varID].name) break;

      if ( varID == varTableUsed )
        {
          struct cdi_generate_varinfo *varInfo
            = (struct cdi_generate_varinfo *) Malloc((size_t)varTableUsed * sizeof(struct cdi_generate_varinfo));

          for ( size_t varID = 0; varID < varTableUsed; varID++ )
            {
              varInfo[varID].varid = varids[varID];
              varInfo[varID].name = vartable[varids[varID]].name;
            }
          qsort(varInfo, varTableUsed, sizeof(varInfo[0]), cdi_generate_cmp_varname);
          for ( size_t varID = 0; varID < varTableUsed; varID++ )
            {
              varids[varID] = varInfo[varID].varid;
            }
          Free(varInfo);
        }
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
  */
551
  for ( size_t index = 0; index < varTableUsed; index++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
    {
553
      const int varid      = varids[index];
554

555
556
557
558
      const int gridID     = vartable[varid].gridID;
      const int param      = vartable[varid].param;
      const int ltype1     = vartable[varid].ltype1;
      const int ltype2     = vartable[varid].ltype2;
559
      int zaxistype  = vartable[varid].zaxistype;
560
      if ( ltype1 == 0 && zaxistype == ZAXIS_GENERIC && cdiDefaultLeveltype != -1 )
561
	zaxistype = cdiDefaultLeveltype;
562
563
      const int lbounds    = vartable[varid].lbounds;
      const int prec       = vartable[varid].prec;
564
565
566
      int instID     = vartable[varid].instID;
      int modelID    = vartable[varid].modelID;
      int tableID    = vartable[varid].tableID;
567
568
569
      const int tsteptype  = vartable[varid].tsteptype;
      const int timave     = vartable[varid].timave;
      const int comptype   = vartable[varid].comptype;
570
571

      double level_sf = 1;
572
      if ( vartable[varid].level_sf != 0 ) level_sf = 1./vartable[varid].level_sf;
573

574
      /* consistency check: test if all subtypes have the same levels: */
575
      const unsigned nlevels = vartable[varid].recordTable[0].nlevels;
576
577
      for ( int isub = 1; isub < vartable[varid].nsubtypes; isub++ ) {
        if ( vartable[varid].recordTable[isub].nlevels != nlevels )
578
          {
579
            fprintf(stderr, "var \"%s\": isub = %d / %d :: "
580
                    "nlevels = %d, vartable[varid].recordTable[isub].nlevels = %d\n",
581
                    vartable[varid].name, isub, vartable[varid].nsubtypes,
582
583
584
585
                    nlevels, vartable[varid].recordTable[isub].nlevels);
            Error("zaxis size must not change for same parameter!");
          }

586
587
        const leveltable_t *t1 = vartable[varid].recordTable[isub-1].levelTable;
        const leveltable_t *t2 = vartable[varid].recordTable[isub  ].levelTable;
588
        for ( unsigned ilev = 0; ilev < nlevels; ilev++ )
589
590
591
592
          if ((t1[ilev].level1 != t2[ilev].level1)  ||
              (t1[ilev].level2 != t2[ilev].level2)  ||
              (t1[ilev].lindex != t2[ilev].lindex))
            {
593
              fprintf(stderr, "var \"%s\", varID=%d: isub = %d / %d :: "
594
                      "nlevels = %d, vartable[varid].recordTable[isub].nlevels = %d\n",
595
                      vartable[varid].name, varid, isub, vartable[varid].nsubtypes,
596
                      nlevels, vartable[varid].recordTable[isub].nlevels);
597
598
599
              Message("t1[ilev].level1=%d / t2[ilev].level1=%d", t1[ilev].level1, t2[ilev].level1);
              Message("t1[ilev].level2=%d / t2[ilev].level2=%d", t1[ilev].level2, t2[ilev].level2);
              Message("t1[ilev].lindex=%d / t2[ilev].lindex=%d", t1[ilev].lindex, t2[ilev].lindex);
600
601
602
603
604
              Error("zaxis type must not change for same parameter!");
            }
      }
      leveltable_t *levelTable = vartable[varid].recordTable[0].levelTable;

605
      if ( ltype1 == 0 && zaxistype == ZAXIS_GENERIC && nlevels == 1 && levelTable[0].level1 == 0 )
606
	zaxistype = ZAXIS_SURFACE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607

608
      double *dlevels = (double *) Malloc(nlevels*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
609

610
      if ( lbounds && zaxistype != ZAXIS_HYBRID && zaxistype != ZAXIS_HYBRID_HALF )
611
	for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
612
613
	  dlevels[levelID] = (level_sf*levelTable[levelID].level1 +
	                      level_sf*levelTable[levelID].level2)/2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614
      else
615
	for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
616
	  dlevels[levelID] = level_sf*levelTable[levelID].level1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
618
619

      if ( nlevels > 1 )
	{
620
          bool linc = true, ldec = true, lsort = false;
621
          for ( unsigned levelID = 1; levelID < nlevels; levelID++ )
622
            {
623
              // check increasing of levels
624
              linc &= (dlevels[levelID] > dlevels[levelID-1]);
625
              // check decreasing of levels
626
627
628
629
              ldec &= (dlevels[levelID] < dlevels[levelID-1]);
            }
          /*
           * always sort pressure z-axis to ensure
630
           * levelTable[levelID1].level1 < levelTable[levelID2].level1 <=> levelID1 > levelID2
631
632
           * unless already sorted in decreasing order
           */
633
          if ( (!linc && !ldec) && zaxistype == ZAXIS_PRESSURE )
634
            {
635
              qsort(levelTable, nlevels, sizeof(leveltable_t), cmpLevelTableInv);
636
637
638
639
              lsort = true;
            }
          /*
           * always sort hybrid and depth-below-land z-axis to ensure
640
           * levelTable[levelID1].level1 < levelTable[levelID2].level1 <=> levelID1 < levelID2
641
642
643
644
645
646
           * unless already sorted in increasing order
           */
          else if ( (!linc && !ldec) ||
                    zaxistype == ZAXIS_HYBRID ||
                    zaxistype == ZAXIS_DEPTH_BELOW_LAND )
            {
647
              qsort(levelTable, nlevels, sizeof(leveltable_t), cmpLevelTable);
648
649
              lsort = true;
            }
650
651
652
653

          if ( lsort )
            {
              if ( lbounds && zaxistype != ZAXIS_HYBRID && zaxistype != ZAXIS_HYBRID_HALF )
654
                for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
655
656
                  dlevels[levelID] = (level_sf*levelTable[levelID].level1 +
                                      level_sf*levelTable[levelID].level2)/2.;
657
              else
658
                for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
659
                  dlevels[levelID] = level_sf*levelTable[levelID].level1;
660
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
662
	}

663
664
      double *dlevels1 = NULL;
      double *dlevels2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
666
      if ( lbounds )
	{
667
	  dlevels1 = (double *) Malloc(nlevels*sizeof(double));
668
	  for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
669
	    dlevels1[levelID] = level_sf*levelTable[levelID].level1;
670
	  dlevels2 = (double *) Malloc(nlevels*sizeof(double));
671
	  for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
672
	    dlevels2[levelID] = level_sf*levelTable[levelID].level2;
673
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674

Uwe Schulzweida's avatar
Uwe Schulzweida committed
675
      const char **cvals = NULL;
Thomas Jahns's avatar
Thomas Jahns committed
676
      const char *unitptr = cdiUnitNamePtr(vartable[varid].level_unit);
677
      int zaxisID = varDefZaxis(vlistID, zaxistype, (int)nlevels, dlevels, cvals, 0, lbounds, dlevels1, dlevels2,
678
                                (int)Vctsize, Vct, NULL, NULL, unitptr, 0, 0, ltype1, ltype2);
679

680
681
      if ( CDI_cmor_mode && nlevels == 1 && zaxistype != ZAXIS_HYBRID ) zaxisDefScalar(zaxisID);

682
683
      if ( zaxisInqType(zaxisID) == ZAXIS_REFERENCE )
        {
684
          if ( numberOfVerticalLevels > 0 ) cdiDefKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NLEV, numberOfVerticalLevels);
685
          if ( numberOfVerticalGrid > 0 ) cdiDefKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NUMBEROFVGRIDUSED, numberOfVerticalGrid);
686
          if ( !cdiUUIDIsNull(uuidVGrid) ) cdiDefKeyBytes(zaxisID, CDI_GLOBAL, CDI_KEY_UUID, uuidVGrid, CDI_UUID_SIZE);
687
688
        }

689
690
691
      if ( lbounds ) Free(dlevels1);
      if ( lbounds ) Free(dlevels2);
      Free(dlevels);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
692

693
      // define new subtype for tile set
694
695
      int tilesetID = CDI_UNDEFID;
      if ( vartable[varid].tiles ) tilesetID = vlistDefTileSubtype(vlistID, vartable[varid].tiles);
696

697
      // generate new variable
698
      int varID = stream_new_var(streamptr, gridID, zaxisID, tilesetID);
699
      varID = vlistDefVarTiles(vlistID, gridID, zaxisID, TIME_VARYING, tilesetID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
700

701
      vlistDefVarTsteptype(vlistID, varID, tsteptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
702
      vlistDefVarParam(vlistID, varID, param);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703
      vlistDefVarDatatype(vlistID, varID, prec);
704
      vlistDefVarTimave(vlistID, varID, timave);
705
      vlistDefVarCompType(vlistID, varID, comptype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706

707
      varCopyKeys(vlistID, varID);
708

709
      if ( vartable[varid].lmissval ) vlistDefVarMissval(vlistID, varID, vartable[varid].missval);
710
      if ( vartable[varid].name )     vlistDefVarName(vlistID, varID, vartable[varid].name);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
711

712
713
      vlist_t *vlistptr = vlist_to_pointer(vlistID);
      for ( int i = 0; i < vartable[varid].opt_grib_nentries; i++ )
714
        {
715
          resize_opt_grib_entries(&vlistptr->vars[varID], vlistptr->vars[varID].opt_grib_nentries+1);
716
          vlistptr->vars[varID].opt_grib_nentries += 1;
717
          const int idx = vlistptr->vars[varID].opt_grib_nentries-1;
718

719
          vlistptr->vars[varID].opt_grib_kvpair[idx] = vartable[varid].opt_grib_kvpair[i];
720
          vlistptr->vars[varID].opt_grib_kvpair[idx].keyword = NULL;
721
722
	  if ( vartable[varid].opt_grib_kvpair[i].keyword )
	    vlistptr->vars[varID].opt_grib_kvpair[idx].keyword =
723
	      strdupx(vartable[varid].opt_grib_kvpair[i].keyword);
724
          vlistptr->vars[varID].opt_grib_kvpair[i].update = true;
725
        }
726
      // note: if the key is not defined, we do not throw an error!
727

728
      if ( cdiDefaultTableID != CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730
	  int pdis, pcat, pnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
731
	  cdiDecodeParam(param, &pnum, &pcat, &pdis);
732
          char name[CDI_MAX_NAME]; name[0] = 0;
733
734
735
          char longname[CDI_MAX_NAME]; longname[0] = 0;
          char units[CDI_MAX_NAME]; units[0] = 0;
          tableInqEntry(cdiDefaultTableID, pnum, -1, name, longname, units);
736
	  if ( name[0] )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
737
	    {
738
	      if ( tableID != CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
740
		{
		  vlistDefVarName(vlistID, varID, name);
741
742
                  if ( longname[0] ) vlistDefVarLongname(vlistID, varID, longname);
                  if ( units[0] ) vlistDefVarUnits(vlistID, varID, units);
743
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
744
745
746
	      else
		tableID = cdiDefaultTableID;
	    }
747
748
	  if ( cdiDefaultModelID != CDI_UNDEFID ) modelID = cdiDefaultModelID;
	  if ( cdiDefaultInstID  != CDI_UNDEFID )  instID = cdiDefaultInstID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
	}

751
752
753
      if ( instID  != CDI_UNDEFID ) vlistDefVarInstitut(vlistID, varID, instID);
      if ( modelID != CDI_UNDEFID ) vlistDefVarModel(vlistID, varID, modelID);
      if ( tableID != CDI_UNDEFID ) vlistDefVarTable(vlistID, varID, tableID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754
755
    }

756
  for ( size_t index = 0; index < varTableUsed; index++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757
    {
758
      int varid = varids[index];
759
      unsigned nlevels = vartable[varid].recordTable[0].nlevels;
760

761
762
      unsigned nsub = vartable[varid].nsubtypes >= 0 ? (unsigned)vartable[varid].nsubtypes : 0U;
      for ( size_t isub = 0; isub < nsub; isub++ )
763
        {
764
765
          sleveltable_t *restrict streamRecordTable = streamptr->vars[index].recordTable + isub;
          leveltable_t *restrict vartableLevelTable = vartable[varid].recordTable[isub].levelTable;
766
          for ( unsigned levelID = 0; levelID < nlevels; levelID++ )
767
            {
768
              streamRecordTable->recordID[levelID] = vartableLevelTable[levelID].recID;
769
              unsigned lindex;
770
              for ( lindex = 0; lindex < nlevels; lindex++ )
771
772
                if ( levelID == (unsigned)vartableLevelTable[lindex].lindex )
                  break;
773
774
              if ( lindex == nlevels )
                Error("Internal problem! lindex not found.");
775
              streamRecordTable->lindex[levelID] = (int)lindex;
776
777
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778
779
    }

780
  Free(varids);
781

Uwe Schulzweida's avatar
Uwe Schulzweida committed
782
783
784
785
786
787
788
789
790
  varFree();
}


void varDefVCT(size_t vctsize, double *vctptr)
{
  if ( Vct == NULL && vctptr != NULL && vctsize > 0 )
    {
      Vctsize = vctsize;
791
      Vct = (double *) Malloc(vctsize*sizeof(double));
792
      memcpy(Vct, vctptr, vctsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793
794
795
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
796

797
void varDefZAxisReference(int nhlev, int nvgrid, unsigned char uuid[CDI_UUID_SIZE])
798
{
799
  numberOfVerticalLevels = nhlev;
800
  numberOfVerticalGrid = nvgrid;
801
  memcpy(uuidVGrid, uuid, CDI_UUID_SIZE);
802
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
803
804


805
bool zaxisCompare(int zaxisID, int zaxistype, int nlevels, bool lbounds, const double *levels, const char *longname, const char *units, int ltype1, int ltype2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
806
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
807
  bool differ = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
808

809
810
811
812
813
  int ltype1_0 = 0, ltype2_0 = -1;
  cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_TYPEOFFIRSTFIXEDSURFACE, &ltype1_0);
  cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_TYPEOFSECONDFIXEDSURFACE, &ltype2_0);
  bool ltype1_is_equal = (ltype1 == ltype1_0);
  bool ltype2_is_equal = (ltype2 == ltype2_0);
814

815
  if ( ltype1_is_equal && ltype2_is_equal && (zaxistype == zaxisInqType(zaxisID) || zaxistype == ZAXIS_GENERIC) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
    {
817
      bool zlbounds = (zaxisInqLbounds(zaxisID, NULL) > 0);