stream_cdf_i.c 145 KB
Newer Older
1
#ifdef HAVE_CONFIG_H
2
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3
4
#endif

5
6
#ifdef HAVE_LIBNETCDF

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include <ctype.h>
8
#include <limits.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
9

Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
#include "dmemory.h"
11
#include "cdi_int.h"
12
#include "cdi_uuid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
13
14
15
16
#include "stream_cdf.h"
#include "cdf_int.h"
#include "varscan.h"
#include "vlist.h"
17
#include "cdf_util.h"
18
#include "cdf_lazy_grid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19
20
21
22
23
24
25


#define  X_AXIS  1
#define  Y_AXIS  2
#define  Z_AXIS  3
#define  T_AXIS  4

26
27
28
#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2

Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
typedef struct {
30
  int     dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
33
  int     ncvarid;
  int     dimtype;
  size_t  len;
34
  char    name[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
}
36
ncdim_t;
37

38
39
40
#define  MAX_COORDVARS  4
#define  MAX_AUXVARS    4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
typedef struct {
42
  int      ncid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
  int      isvar;
  bool     ignore;
45
46
  bool     isx;
  bool     isy;
47
  bool     isc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
51
52
  bool     islon;
  bool     islat;
  bool     islev;
  bool     istime;
  bool     warn;
53
  bool     calendar;
54
  bool     climatology;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
  bool     lformulaterms;
56
  int      param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
  int      code;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
  int      tabnum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
61
  int      bounds;
  int      gridID;
  int      zaxisID;
62
  int      timetype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
65
66
67
68
69
  int      gridtype;
  int      zaxistype;
  int      xdim;
  int      ydim;
  int      zdim;
  int      xvarid;
  int      yvarid;
70
  int      rpvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
  int      zvarid;
72
  int      cvarids[MAX_COORDVARS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
  int      tvarid;
74
  int      psvarid;
75
  int      p0varid;
76
  int      ncoordvars;
77
78
79
  int      coordvarids[MAX_COORDVARS];
  int      nauxvars;
  int      auxvarids[MAX_AUXVARS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
  int      cellarea;
  int      tableID;
  int      truncation;
83
  int      position;
84
  int      numLPE;
85
86
  bool     defmissval;
  bool     deffillval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
89
  int      xtype;
  int      gmapid;
  int      positive;
90
  int      ndims;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
  int      dimids[8];
  int      dimtype[8];
93
  size_t   chunks[8];
94
95
  int      chunked;
  int      chunktype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
  int      natts;
97
  int      deflate;
98
99
  bool     lunsigned;
  bool     lvalidrange;
Thomas Jahns's avatar
Thomas Jahns committed
100
  int     *atts;
101
102
  size_t   vctsize;
  double  *vct;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
  double   missval;
104
  double   fillval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
  double   addoffset;
  double   scalefactor;
Deike Kleberg's avatar
Deike Kleberg committed
107
  double   validrange[2];
108
109
110
111
  char     name[CDI_MAX_NAME];
  char     longname[CDI_MAX_NAME];
  char     stdname[CDI_MAX_NAME];
  char     units[CDI_MAX_NAME];
112
  char     extra[CDI_MAX_NAME];
113
114
115
  int      typeOfEnsembleForecast;
  int      numberOfForecastsInEnsemble;
  int      perturbationNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
}
117
118
ncvar_t;

119
120

static
121
void scanTimeString(const char *ptu, int64_t *rdate, int *rtime)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
{
123
  int year = 1, month = 1, day = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
  int hour = 0, minute = 0, second = 0;
125
  int v1 = 1, v2 = 1, v3 = 1;
126

Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
  if ( *ptu ) sscanf(ptu, "%d-%d-%d %d:%d:%d", &v1, &v2, &v3, &hour, &minute, &second);
128
129
130
131
132
133
134
135
136
137

  if ( v3 > 999 && v1 < 32 )
    { year = v3; month = v2; day = v1; }
  else
    { year = v1; month = v2; day = v3; }

  *rdate = cdiEncodeDate(year, month, day);
  *rtime = cdiEncodeTime(hour, minute, second);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
140
static
int scanTimeUnit(const char *unitstr)
{
141
142
143
  const size_t len = strlen(unitstr);
  const int timeunit = get_timeunit(len, unitstr);
  if ( timeunit == -1 ) Message("Unsupported TIMEUNIT: %s!", unitstr);
144
  return timeunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146
}

147
148
149
static
void setForecastTime(const char *timestr, taxis_t *taxis)
{
150
  const size_t len = strlen(timestr);
Thomas Jahns's avatar
Thomas Jahns committed
151
152
153
154
  if ( len != 0 )
    scanTimeString(timestr, &taxis->fdate, &taxis->ftime);
  else
    taxis->fdate = taxis->ftime = 0;
155
156
157
158
159
}

static
int setBaseTime(const char *timeunits, taxis_t *taxis)
{
160
  int taxistype = TAXIS_ABSOLUTE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161

162
  size_t len = strlen(timeunits);
163
164
  while ( isspace(*timeunits) && len ) { timeunits++; len--; }

165
  char *tu = (char *)Malloc((len+1) * sizeof(char));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166

Thomas Jahns's avatar
Thomas Jahns committed
167
168
  for ( size_t i = 0; i < len; i++ ) tu[i] = (char)tolower((int)timeunits[i]);
  tu[len] = 0;
169

Thomas Jahns's avatar
Thomas Jahns committed
170
  int timeunit = get_timeunit(len, tu);
171
  if ( timeunit == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
    {
173
      Message("Unsupported TIMEUNIT: %s!", timeunits);
174
      return 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
176
    }

Thomas Jahns's avatar
Thomas Jahns committed
177
  size_t pos = 0;
178
  while ( pos < len && !isspace(tu[pos]) ) ++pos;
Thomas Jahns's avatar
Thomas Jahns committed
179
  if ( tu[pos] )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
    {
Thomas Jahns's avatar
Thomas Jahns committed
181
      while ( isspace(tu[pos]) ) ++pos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182

183
      if ( strStartsWith(tu+pos, "since") ) taxistype = TAXIS_RELATIVE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
      while ( pos < len && !isspace(tu[pos]) ) ++pos;
Thomas Jahns's avatar
Thomas Jahns committed
186
      if ( tu[pos] )
187
        {
Thomas Jahns's avatar
Thomas Jahns committed
188
          while ( isspace(tu[pos]) ) ++pos;
189

190
          if ( taxistype == TAXIS_ABSOLUTE )
191
            {
192
              if ( timeunit == TUNIT_DAY )
193
                {
194
                  if ( !strStartsWith(tu+pos, "%y%m%d.%f") )
195
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
                      Warning("Unsupported format %s for TIMEUNIT day!", tu+pos);
197
198
                      timeunit = -1;
                    }
199
                }
200
              else if ( timeunit == TUNIT_MONTH )
201
                {
202
                  if ( !strStartsWith(tu+pos, "%y%m.%f") )
203
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
                      Warning("Unsupported format %s for TIMEUNIT month!", tu+pos);
205
206
                      timeunit = -1;
                    }
207
                }
208
209
210
211
212
213
214
215
              else if ( timeunit == TUNIT_YEAR )
                {
                  if ( !strStartsWith(tu+pos, "%y.%f") )
                    {
                      Warning("Unsupported format %s for TIMEUNIT year!", tu+pos);
                      timeunit = -1;
                    }
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
217
218
219
              else
                {
                  Warning("Unsupported format for time units: %s!", tu);
                }
220
            }
221
          else if ( taxistype == TAXIS_RELATIVE )
222
            {
223
224
              int64_t rdate = -1;
              int rtime = -1;
Thomas Jahns's avatar
Thomas Jahns committed
225
              scanTimeString(tu+pos, &rdate, &rtime);
226
227
              taxis->rdate = rdate;
              taxis->rtime = rtime;
228

229
              if ( CDI_Debug ) Message("rdate = %lld  rtime = %d", rdate, rtime);
230
            }
231
        }
232
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233

234
  taxis->type = taxistype;
235
  taxis->unit = timeunit;
Deike Kleberg's avatar
Deike Kleberg committed
236

237
  Free(tu);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238

239
  if ( CDI_Debug ) Message("taxistype = %d  unit = %d", taxistype, timeunit);
Deike Kleberg's avatar
Deike Kleberg committed
240

241
242
  return 0;
}
243

244
static
245
bool xtypeIsText(int xtype)
246
{
247
  const bool isText = (xtype == NC_CHAR)
248
#ifdef HAVE_NETCDF4
Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
    || (xtype == NC_STRING)
250
251
252
253
254
255
256
257
#endif
    ;
  return isText;
}

static
bool xtypeIsFloat(nc_type xtype)
{
258
  const bool isFloat = xtype == NC_FLOAT || xtype == NC_DOUBLE;
259
260
261
262
263
264
  return isFloat;
}

static
bool xtypeIsInt(nc_type xtype)
{
265
  const bool isInt = xtype == NC_SHORT || xtype == NC_INT
266
            || xtype == NC_BYTE
267
#ifdef HAVE_NETCDF4
268
269
270
271
272
273
274
275
276
            || xtype == NC_USHORT || xtype == NC_UINT
            || xtype == NC_UBYTE
#endif
             ;

  return isInt;
}

static
277
int cdfInqDatatype(stream_t *streamptr, int xtype, bool lunsigned)
278
279
280
{
  int datatype = -1;

281
#ifdef HAVE_NETCDF4
282
283
284
285
  if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
#endif

  if      ( xtype == NC_BYTE   )  datatype = CDI_DATATYPE_INT8;
286
  else if ( xtype == NC_CHAR   )  datatype = CDI_DATATYPE_UINT8;
287
288
289
290
  else if ( xtype == NC_SHORT  )  datatype = CDI_DATATYPE_INT16;
  else if ( xtype == NC_INT    )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_FLOAT  )  datatype = CDI_DATATYPE_FLT32;
  else if ( xtype == NC_DOUBLE )  datatype = CDI_DATATYPE_FLT64;
291
#ifdef HAVE_NETCDF4
292
293
294
295
296
297
  else if ( xtype == NC_UBYTE  )  datatype = CDI_DATATYPE_UINT8;
  else if ( xtype == NC_LONG   )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_USHORT )  datatype = CDI_DATATYPE_UINT16;
  else if ( xtype == NC_UINT   )  datatype = CDI_DATATYPE_UINT32;
  else if ( xtype == NC_INT64  )  datatype = CDI_DATATYPE_FLT64;
  else if ( xtype == NC_UINT64 )  datatype = CDI_DATATYPE_FLT64;
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
  else
    {
      if (xtype != streamptr->nc_complex_float_id && xtype != streamptr->nc_complex_double_id)
        {
          int isUserDefinedType = false;
#ifdef NC_FIRSTUSERTYPEID
          isUserDefinedType = xtype >= NC_FIRSTUSERTYPEID;
#endif
          if (isUserDefinedType)
            {
              int fileID = streamptr->fileID;
              size_t nfields = 0, compoundsize = 0;
              int status = nc_inq_compound(fileID, xtype, NULL, &compoundsize, &nfields);
              if (status == NC_NOERR && nfields == 2 && (compoundsize == 8 || compoundsize == 16))
                {
                  nc_type field_type1 = -1, field_type2 = -1;
                  int field_dims1 = 0, field_dims2 = 0;
                  nc_inq_compound_field(fileID, xtype, 0, NULL, NULL, &field_type1, &field_dims1, NULL);
                  nc_inq_compound_field(fileID, xtype, 1, NULL, NULL, &field_type2, &field_dims2, NULL);
                  if (field_type1 == field_type2 && field_dims1 == 0 && field_dims2 == 0)
                    {
                      if      (field_type1 == NC_FLOAT)  streamptr->nc_complex_float_id = xtype;
                      else if (field_type1 == NC_DOUBLE) streamptr->nc_complex_double_id = xtype;
                    }
                }
            }
        }
      if      ( xtype == streamptr->nc_complex_float_id  )  datatype = CDI_DATATYPE_CPX32;
      else if ( xtype == streamptr->nc_complex_double_id )  datatype = CDI_DATATYPE_CPX64;
    }
328
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329

330
331
332
333
334
335
  return datatype;
}

static
void cdfGetAttInt(int fileID, int ncvarid, const char *attname, size_t attlen, int *attint)
{
336
  *attint = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337

338
339
  nc_type atttype;
  size_t nc_attlen;
340
341
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Deike Kleberg's avatar
Deike Kleberg committed
342

343
  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
344
    {
345
      const bool lalloc = nc_attlen > attlen;
346
      int *pintatt = lalloc ? (int *)(Malloc(nc_attlen*sizeof(int))) : attint;
347
      cdf_get_att_int(fileID, ncvarid, attname, pintatt);
348
      if ( lalloc )
349
        {
350
          memcpy(attint, pintatt, attlen*sizeof(int));
351
          Free(pintatt);
352
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
354
355
    }
}

356
static
357
void cdfGetAttDouble(int fileID, int ncvarid, char *attname, size_t attlen, double *attdouble)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
358
{
359
  *attdouble = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360

361
362
  nc_type atttype;
  size_t nc_attlen;
363
364
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365

366
  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
367
    {
368
      const bool lalloc = nc_attlen > attlen;
369
      double *pdoubleatt = lalloc ? (double*)Malloc(nc_attlen*sizeof(double)) : attdouble;
370
      cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);
371
      if ( lalloc )
372
        {
373
          memcpy(attdouble, pdoubleatt, attlen*sizeof(double));
374
375
376
377
          Free(pdoubleatt);
        }
    }
}
378

379
static
380
381
382
383
384
bool cdfCheckAttText(int fileID, int ncvarid, const char *attname)
{
  bool status = false;
  nc_type atttype;

385
386
  const int status_nc = nc_inq_atttype(fileID, ncvarid, attname, &atttype);
  if ( status_nc == NC_NOERR && (atttype == NC_CHAR
387
#ifdef HAVE_NETCDF4
388
389
390
391
392
393
394
395
396
397
398
           || atttype == NC_STRING
#endif
           ) )
    {
      status = true;
    }

  return status;
}

static
399
void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
400
{
401
402
  atttext[0] = 0;

403
404
405
406
  nc_type atttype;
  size_t nc_attlen;
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
407

408
  if ( atttype == NC_CHAR )
409
    {
410
411
412
413
      char attbuf[65636];
      if ( nc_attlen < sizeof(attbuf) )
        {
          cdf_get_att_text(fileID, ncvarid, attname, attbuf);
414

415
          if ( nc_attlen > (attlen-1) ) nc_attlen = (attlen-1);
416

417
418
419
          attbuf[nc_attlen++] = 0;
          memcpy(atttext, attbuf, nc_attlen);
        }
420
    }
421
#ifdef HAVE_NETCDF4
422
  else if ( atttype == NC_STRING )
423
    {
424
      if ( nc_attlen == 1 )
425
        {
426
427
          char *attbuf = NULL;
          cdf_get_att_string(fileID, ncvarid, attname, &attbuf);
428

429
          size_t ssize = strlen(attbuf) + 1;
430
          if ( ssize > attlen ) ssize = attlen;
431
432
433
434
          memcpy(atttext, attbuf, ssize);
          atttext[ssize - 1] = 0;
          Free(attbuf);
        }
435
    }
436
437
#endif
}
438

439

440
void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
{
442
443
  const bool laddoffset   = IS_NOT_EQUAL(addoffset, 0);
  const bool lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
444

445
  if ( laddoffset && lscalefactor )
446
    {
447
      for (size_t i = 0; i < size; ++i) data[i] = data[i] * scalefactor + addoffset;
448
449
450
    }
  else if (lscalefactor)
    {
451
      for (size_t i = 0; i < size; ++i) data[i] *= scalefactor;
452
453
454
    }
  else if (laddoffset)
    {
455
      for (size_t i = 0; i < size; ++i) data[i] += addoffset;
456
457
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
458

459
460
static
void cdfCreateRecords(stream_t *streamptr, int tsID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
461
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
  if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463

Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
  if ( streamptr->tsteps[tsID].nallrecs > 0 ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
465

466
467
  int vlistID  = streamptr->vlistID;

468
469
470
  tsteps_t* sourceTstep = streamptr->tsteps;
  tsteps_t* destTstep = sourceTstep + tsID;

471
472
  const int nvars = vlistNvars(vlistID);
  const int nrecs = vlistNrecs(vlistID);
473
474
  if ( nrecs <= 0 ) return;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
476
  if ( tsID == 0 )
    {
477
      const int nvrecs = nrecs; /* use all records at first timestep */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478

Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
      streamptr->nrecs += nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
480

481
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
482
483
484
      destTstep->nrecs      = nrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
485
      destTstep->curRecID   = CDI_UNDEFID;
486
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs*sizeof (int));;
487
      for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488

489
      record_t *records = destTstep->records;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490

491
      for ( int varID = 0, recID = 0; varID < nvars; varID++ )
492
        {
493
494
          const int zaxisID = vlistInqVarZaxis(vlistID, varID);
          const int nlev    = zaxisInqSize(zaxisID);
495
          for ( int levelID = 0; levelID < nlev; levelID++ )
496
497
            {
              recordInitEntry(&records[recID]);
498
499
              records[recID].varID   = (short)varID;
              records[recID].levelID = (short)levelID;
500
501
502
              recID++;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
503
504
505
    }
  else if ( tsID == 1 )
    {
506
507
      int nvrecs = 0;
      for ( int varID = 0; varID < nvars; varID++ )
508
        {
509
          if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
510
            {
511
              int zaxisID = vlistInqVarZaxis(vlistID, varID);
512
513
514
              nvrecs += zaxisInqSize(zaxisID);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
515

Uwe Schulzweida's avatar
Uwe Schulzweida committed
516
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517

518
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
519
520
521
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
522
      destTstep->curRecID   = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523

524
      memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525

526
      if ( nvrecs )
527
        {
528
          destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
529
          for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
530
            {
531
              int varID = destTstep->records[recID].varID;
532
              if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
533
                {
534
                  destTstep->recIDs[vrecID++] = recID;
535
536
537
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
539
540
    }
  else
    {
541
      if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);
542

543
      const int nvrecs = streamptr->tsteps[1].nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
544

Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546

547
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
548
549
550
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
551
      destTstep->curRecID   = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552

553
      memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554

555
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs * sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556

557
      memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
559
560
    }
}

561
static
562
int cdf_time_dimid(int fileID, int ndims, int nvars, ncdim_t *ncdims)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
{
564
565
  char dimname[CDI_MAX_NAME];

566
  for ( int dimid = 0; dimid < ndims; ++dimid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567
    {
568
      dimname[0] = 0;
569
      cdf_inq_dimname(fileID, ncdims[dimid].dimid, dimname);
570
      if ( strIsEqual(dimname, "time") || strIsEqual(dimname, "Time") ) return dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
    }

573
  for ( int varid = 0; varid < nvars; ++varid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574
    {
575
      nc_type xtype;
576
      int nvdims, nvatts, dimids[9];
577
      cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
578
579
580
581
582
583
584
585
      for ( int i = 0; i < nvdims; ++i )
        for ( int dimid = 0; dimid < ndims; ++dimid )
          if ( ncdims[dimid].dimid == dimids[i] )
            {
              dimids[i] = dimid;
              break;
            }

586
      if ( nvdims == 1 )
587
        {
588
589
          char sbuf[CDI_MAX_NAME];
          for ( int iatt = 0; iatt < nvatts; ++iatt )
590
            {
591
              sbuf[0] = 0;
592
593
              cdf_inq_attname(fileID, varid, iatt, sbuf);
              if ( strncmp(sbuf, "units", 5) == 0 )
594
                {
595
                  cdfGetAttText(fileID, varid, "units", sizeof(sbuf), sbuf);
596
                  strToLower(sbuf);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
597

598
                  if ( is_time_units(sbuf) ) return dimids[0];
599
600
601
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
603
    }

604
  return CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
606
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
static
608
void init_ncdims(long ndims, ncdim_t *ncdims)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
609
{
610
  for ( long ncdimid = 0; ncdimid < ndims; ncdimid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
    {
612
      ncdims[ncdimid].dimid        = CDI_UNDEFID;
613
614
      ncdims[ncdimid].ncvarid      = CDI_UNDEFID;
      ncdims[ncdimid].dimtype      = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
616
617
618
619
620
      ncdims[ncdimid].len          = 0;
      ncdims[ncdimid].name[0]      = 0;
    }
}

static
621
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
{
623
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
624
    {
625
626
      ncvars[ncvarid].ncid            = CDI_UNDEFID;
      ncvars[ncvarid].isvar           = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
627
      ncvars[ncvarid].ignore          = false;
628
629
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
630
      ncvars[ncvarid].isc             = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631
632
633
634
635
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
636
      ncvars[ncvarid].calendar        = false;
637
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
638
      ncvars[ncvarid].lformulaterms   = false;
639
      ncvars[ncvarid].timetype        = TIME_CONSTANT;
640
641
      ncvars[ncvarid].param           = CDI_UNDEFID;
      ncvars[ncvarid].code            = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
      ncvars[ncvarid].tabnum          = 0;
643
644
645
646
647
648
649
650
651
652
      ncvars[ncvarid].bounds          = CDI_UNDEFID;
      ncvars[ncvarid].gridID          = CDI_UNDEFID;
      ncvars[ncvarid].zaxisID         = CDI_UNDEFID;
      ncvars[ncvarid].gridtype        = CDI_UNDEFID;
      ncvars[ncvarid].zaxistype       = CDI_UNDEFID;
      ncvars[ncvarid].xdim            = CDI_UNDEFID;
      ncvars[ncvarid].ydim            = CDI_UNDEFID;
      ncvars[ncvarid].zdim            = CDI_UNDEFID;
      ncvars[ncvarid].xvarid          = CDI_UNDEFID;
      ncvars[ncvarid].yvarid          = CDI_UNDEFID;
653
      ncvars[ncvarid].rpvarid         = CDI_UNDEFID;
654
655
656
657
      ncvars[ncvarid].zvarid          = CDI_UNDEFID;
      ncvars[ncvarid].tvarid          = CDI_UNDEFID;
      ncvars[ncvarid].psvarid         = CDI_UNDEFID;
      ncvars[ncvarid].p0varid         = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
658
      ncvars[ncvarid].ncoordvars      = 0;
659
      for ( int i = 0; i < MAX_COORDVARS; ++i )
660
661
662
663
        {
          ncvars[ncvarid].coordvarids[i]  = CDI_UNDEFID;
          ncvars[ncvarid].cvarids[i]      = CDI_UNDEFID;
        }
664
665
      ncvars[ncvarid].nauxvars      = 0;
      for ( int i = 0; i < MAX_AUXVARS; ++i )
666
667
668
        ncvars[ncvarid].auxvarids[i]  = CDI_UNDEFID;
      ncvars[ncvarid].cellarea        = CDI_UNDEFID;
      ncvars[ncvarid].tableID         = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
669
670
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
671
      ncvars[ncvarid].gmapid          = CDI_UNDEFID;
672
673
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
      ncvars[ncvarid].truncation      = 0;
675
      ncvars[ncvarid].position        = 0;
676
      ncvars[ncvarid].numLPE          = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677
      ncvars[ncvarid].positive        = 0;
678
      ncvars[ncvarid].chunked         = 0;
679
      ncvars[ncvarid].chunktype       = CDI_UNDEFID;
680
681
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
682
      ncvars[ncvarid].missval         = 0;
683
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
684
685
686
687
688
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
689
690
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
691
692
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
693
694
695
      ncvars[ncvarid].typeOfEnsembleForecast      = -1;
      ncvars[ncvarid].numberOfForecastsInEnsemble = -1;
      ncvars[ncvarid].perturbationNumber          = -1;
696
697
698
699
700
      memset(ncvars[ncvarid].name, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].longname, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].stdname, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].units, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].extra, 0, CDI_MAX_NAME);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
701
702
703
    }
}

704
static
705
void cdf_set_var(ncvar_t *ncvars, int ncvarid, short isvar)
706
{
707
  if ( ncvars[ncvarid].isvar != CDI_UNDEFID &&
708
       ncvars[ncvarid].isvar != isvar   &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
709
       ncvars[ncvarid].warn  == false )
710
711
712
713
    {
      if ( ! ncvars[ncvarid].ignore )
        Warning("Inconsistent variable definition for %s!", ncvars[ncvarid].name);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
714
      ncvars[ncvarid].warn = true;
715
716
717
718
719
720
721
      isvar = FALSE;
    }

  ncvars[ncvarid].isvar = isvar;
}

static
722
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
723
{
724
  if ( ncvars[ncvarid].dimtype[dimid] != CDI_UNDEFID &&
725
726
727
728
729
730
731
732
733
734
       ncvars[ncvarid].dimtype[dimid] != dimtype )
    {
      Warning("Inconsistent dimension definition for %s! dimid = %d;  type = %d;  newtype = %d",
              ncvars[ncvarid].name, dimid, ncvars[ncvarid].dimtype[dimid], dimtype);
    }

  ncvars[ncvarid].dimtype[dimid] = dimtype;
}

static
735
void scan_hybrid_formulaterms(int ncid, int ncfvarid, int *avarid, int *bvarid, int *psvarid, int *p0varid)
736
{
737
  *avarid  = -1;
738
739
  *bvarid  = -1;
  *psvarid = -1;
740
  *p0varid = -1;
741

742
  char attstring[1024];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
743
744
  cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
  char *pstring = attstring;
745

746
  bool lstop = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
  for ( int i = 0; i < 4; i++ )
748
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
750
751
752
753
754
755
756
757
758
759
760
761
762
      while ( isspace((int) *pstring) ) pstring++;
      if ( *pstring == 0 ) break;
      char *tagname = pstring;
      while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
      if ( *pstring == 0 ) lstop = true;
      *pstring++ = 0;

      while ( isspace((int) *pstring) ) pstring++;
      if ( *pstring == 0 ) break;
      char *varname = pstring;
      while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
      if ( *pstring == 0 ) lstop = true;
      *pstring++ = 0;

763
      int dimvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
765
      int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
      if ( status_nc == NC_NOERR )
766
        {
767
768
769
770
771
          if      ( strIsEqual(tagname, "ap:") ) *avarid  = dimvarid;
          else if ( strIsEqual(tagname, "a:")  ) *avarid  = dimvarid;
          else if ( strIsEqual(tagname, "b:")  ) *bvarid  = dimvarid;
          else if ( strIsEqual(tagname, "ps:") ) *psvarid = dimvarid;
          else if ( strIsEqual(tagname, "p0:") ) *p0varid = dimvarid;
772
        }
773
      else if ( !strIsEqual(tagname, "ps:") )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
774
775
776
777
778
        {
          Warning("%s - %s", nc_strerror(status_nc), varname);
        }

      if ( lstop ) break;
779
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
780
781
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
static
void readVCT(int ncid, int ndims2, size_t dimlen, size_t dimlen2, int avarid2, int bvarid2, double *vct)
{
  double *abuf = (double*) Malloc(dimlen*2*sizeof(double));
  double *bbuf = (double*) Malloc(dimlen*2*sizeof(double));
  cdf_get_var_double(ncid, avarid2, abuf);
  cdf_get_var_double(ncid, bvarid2, bbuf);

  if ( ndims2 == 2 )
    {
      for ( size_t i = 0; i < dimlen; ++i )
        {
          vct[i] = abuf[i*2];
          vct[i+dimlen+1] = bbuf[i*2];
        }
      vct[dimlen]     = abuf[dimlen*2-1];
      vct[dimlen*2+1] = bbuf[dimlen*2-1];
    }
  else
    {
      for ( size_t i = 0; i < dimlen2; ++i )
        {
          vct[i] = abuf[i];
          vct[i+dimlen+1] = bbuf[i];
        }
    }

  Free(abuf);
  Free(bbuf);
}

813
static
814
bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
815
{
816
  bool status = false;
817
818
  ncvar_t *ncvar = &ncvars[ncvarid];

819
  if ( strIsEqual(ncvar->stdname, "atmosphere_hybrid_sigma_pressure_coordinate") )
820
    {
821
      CDI_convention = CDI_CONVENTION_CF;
822

823
      status = true;
824
      ncvar->zaxistype = ZAXIS_HYBRID;
825
      //int ndims = ncvar->ndims;
826
      int dimid = ncvar->dimids[0];
Thomas Jahns's avatar
Thomas Jahns committed
827
      size_t dimlen = ncdims[dimid].len;
828
      int avarid1 = -1, bvarid1 = -1, psvarid1 = -1, p0varid1 = -1;
829
      int ncfvarid = ncvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
      if ( ncvars[ncfvarid].lformulaterms )
831
832
        scan_hybrid_formulaterms(ncid, ncfvarid, &avarid1, &bvarid1, &psvarid1, &p0varid1);
      // printf("avarid1, bvarid1, psvarid1, p0varid1 %d %d %d %d\n", avarid1, bvarid1, psvarid1, p0varid1);
833
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
834
835
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
836
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
837

838
      if ( ncvar->bounds != CDI_UNDEFID && ncvars[ncvar->bounds].lformulaterms )
839
840
        {
          ncfvarid = ncvar->bounds;
841
          int avarid2 = -1, bvarid2 = -1, psvarid2 = -1, p0varid2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
842
          if ( ncvars[ncfvarid].lformulaterms )
843
844
            scan_hybrid_formulaterms(ncid, ncfvarid, &avarid2, &bvarid2, &psvarid2, &p0varid2);
          // printf("avarid2, bvarid2, psvarid2, p0varid2 %d %d %d %d\n", avarid2, bvarid2, psvarid2, p0varid2);
845
          if ( avarid2 != -1 && bvarid2 != -1 )
846
            {
847
              ncvars[avarid2].isvar = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
848
              ncvars[bvarid2].isvar = FALSE;
849

850
851
852
              const int ndims2 = ncvars[avarid2].ndims;
              const int dimid2 = ncvars[avarid2].dimids[0];
              const size_t dimlen2 = ncdims[dimid2].len;
853
854
855

              if ( (ndims2 == 2 && dimid == ncvars[avarid2].dimids[0] ) ||
                   (ndims2 == 1 && dimlen == dimlen2-1 ) )
856
                {
857
                  double px = 1;
858
                  if ( p0varid1 != -1 && p0varid1 == p0varid2 )
859
860
                    cdf_get_var_double(ncid, p0varid2, &px);