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

5
6
#ifdef HAVE_LIBNETCDF

7
//#define TEST_GROUPS 1
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8

Uwe Schulzweida's avatar
Uwe Schulzweida committed
9
#include <ctype.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10

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

#undef  UNDEFID
#define UNDEFID  CDI_UNDEFID

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

30
31
32
#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2

Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
34
35
36
typedef struct {
  int     ncvarid;
  int     dimtype;
  size_t  len;
37
  char    name[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
}
39
ncdim_t;
40
41
42
#define  MAX_COORDVARS  4
#define  MAX_AUXVARS    4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
typedef struct {
44
  int      ncid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
46
47
48
49
50
51
  int      isvar;
  bool     ignore;
  bool     islon;
  bool     islat;
  bool     islev;
  bool     istime;
  bool     warn;
52
  bool     climatology;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
54
  bool     lformula;
  bool     lformulaterms;
55
  int      tsteptype;
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
62
63
64
65
66
67
68
69
70
  int      bounds;
  int      gridID;
  int      zaxisID;
  int      gridtype;
  int      zaxistype;
  int      xdim;
  int      ydim;
  int      zdim;
  int      xvarid;
  int      yvarid;
  int      zvarid;
  int      tvarid;
71
  int      psvarid;
72
  int      p0varid;
73
  int      ncoordvars;
74
75
76
  int      coordvarids[MAX_COORDVARS];
  int      nauxvars;
  int      auxvarids[MAX_AUXVARS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79
80
  int      cellarea;
  int      calendar;
  int      tableID;
  int      truncation;
81
  int      position;
82
83
  bool     defmissval;
  bool     deffillval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
85
86
87
88
89
  int      xtype;
  int      ndims;
  int      gmapid;
  int      positive;
  int      dimids[8];
  int      dimtype[8];
90
91
92
  int      chunks[8];
  int      chunked;
  int      chunktype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
  int      natts;
94
  int      deflate;
95
96
  bool     lunsigned;
  bool     lvalidrange;
Thomas Jahns's avatar
Thomas Jahns committed
97
  int     *atts;
98
99
  size_t   vctsize;
  double  *vct;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
  double   missval;
101
  double   fillval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
  double   addoffset;
  double   scalefactor;
Deike Kleberg's avatar
Deike Kleberg committed
104
  double   validrange[2];
105
106
107
108
  char     name[CDI_MAX_NAME];
  char     longname[CDI_MAX_NAME];
  char     stdname[CDI_MAX_NAME];
  char     units[CDI_MAX_NAME];
109
  char     extra[CDI_MAX_NAME];
110
  ensinfo_t   *ensdata;    /* Ensemble information */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
}
112
113
ncvar_t;

114
115

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

  *rdate = 0;
  *rtime = 0;

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
  if ( *ptu )
    {
      v1 = atoi(ptu);
      if ( v1 < 0 ) ptu++;
      while ( isdigit((int) *ptu) ) ptu++;
      if ( *ptu )
        {
          v2 = atoi(++ptu);
          while ( isdigit((int) *ptu) ) ptu++;
          if ( *ptu )
            {
              v3 = atoi(++ptu);
              while ( isdigit((int) *ptu) ) ptu++;
            }
        }
    }
141
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

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

  while ( isspace((int) *ptu) ) ptu++;

  if ( *ptu )
    {
      while ( ! isdigit((int) *ptu) ) ptu++;

      hour = atoi(ptu);
      while ( isdigit((int) *ptu) ) ptu++;
      if ( *ptu == ':' )
        {
          ptu++;
          minute = atoi(ptu);
          while ( isdigit((int) *ptu) ) ptu++;
          if ( *ptu == ':' )
            {
              ptu++;
              second = atoi(ptu);
            }
        }
    }

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174
static
int scanTimeUnit(const char *unitstr)
{
175
  size_t len = strlen(unitstr);
176
  int timeunit = get_timeunit(len, unitstr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
178
179
  if ( timeunit == -1 )
    Message("Unsupported TIMEUNIT: %s!", unitstr);

180
  return timeunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
}

183
184
185
static
void setForecastTime(const char *timestr, taxis_t *taxis)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
188
  (*taxis).fdate = 0;
  (*taxis).ftime = 0;

189
  int len = (int) strlen(timestr);
190
191
192
193
194
195
196
197
198
199
200
201
  if ( len == 0 ) return;

  int fdate = 0, ftime = 0;
  scanTimeString(timestr, &fdate, &ftime);

  (*taxis).fdate = fdate;
  (*taxis).ftime = ftime;
}

static
int setBaseTime(const char *timeunits, taxis_t *taxis)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
203
204
  int timetype = TAXIS_ABSOLUTE;
  int rdate = -1, rtime = -1;

205
  size_t len = strlen(timeunits);
206
  char *tu = (char *) Malloc((len+1) * sizeof (char));
207
  memcpy(tu, timeunits, (len+1) * sizeof (char));
208
  char *ptu = tu;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209

210
  for ( size_t i = 0; i < len; i++ ) ptu[i] = (char)tolower((int) ptu[i]);
211

212
  int timeunit = get_timeunit(len, ptu);
213
  if ( timeunit == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
    {
215
      Message("Unsupported TIMEUNIT: %s!", timeunits);
216
      return 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
218
219
220
221
    }

  while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
  if ( *ptu )
    {
222
      while ( isspace(*ptu) ) ptu++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223

224
      if ( memcmp(ptu, "as", 2) == 0 )
225
        timetype = TAXIS_ABSOLUTE;
226
      else if ( memcmp(ptu, "since", 5) == 0 )
227
        timetype = TAXIS_RELATIVE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
230

      while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
      if ( *ptu )
231
232
233
234
235
236
237
238
239
240
241
        {
          while ( isspace(*ptu) ) ptu++;

          if ( timetype == TAXIS_ABSOLUTE )
            {
              if ( memcmp(ptu, "%y%m%d.%f", 9) != 0 && timeunit == TUNIT_DAY )
                {
                  Message("Unsupported format %s for TIMEUNIT day!", ptu);
                  timeunit = -1;
                }
              else if ( memcmp(ptu, "%y%m.%f", 7) != 0 && timeunit == TUNIT_MONTH )
242
243
244
245
246
247
248
249
                {
                  Message("Unsupported format %s for TIMEUNIT month!", ptu);
                  timeunit = -1;
                }
            }
          else if ( timetype == TAXIS_RELATIVE )
            {
              scanTimeString(ptu, &rdate, &rtime);
250

251
252
              (*taxis).rdate = rdate;
              (*taxis).rtime = rtime;
253

254
255
256
              if ( CDI_Debug )
                Message("rdate = %d  rtime = %d", rdate, rtime);
            }
257
        }
258
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259

260
261
  (*taxis).type = timetype;
  (*taxis).unit = timeunit;
Deike Kleberg's avatar
Deike Kleberg committed
262

263
  Free(tu);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264

265
266
  if ( CDI_Debug )
    Message("timetype = %d  unit = %d", timetype, timeunit);
Deike Kleberg's avatar
Deike Kleberg committed
267

268
269
  return 0;
}
270

271
272
273
274
275
static
void cdfGetAttInt(int fileID, int ncvarid, const char *attname, int attlen, int *attint)
{
  nc_type atttype;
  size_t nc_attlen;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276

277
  *attint = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
278

279
280
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Deike Kleberg's avatar
Deike Kleberg committed
281

282
283
284
285
  if ( atttype != NC_CHAR )
    {
      int *pintatt = (int)nc_attlen > attlen
        ? (int *)(Malloc(nc_attlen * sizeof (int))) : attint;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286

287
      cdf_get_att_int(fileID, ncvarid, attname, pintatt);
288

289
290
291
292
      if ( (int)nc_attlen > attlen )
        {
          memcpy(attint, pintatt, (size_t)attlen * sizeof (int));
          Free(pintatt);
293
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
295
296
    }
}

297
static
298
void cdfGetAttDouble(int fileID, int ncvarid, char *attname, int attlen, double *attdouble)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
{
300
301
  nc_type atttype;
  size_t nc_attlen;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302

303
  *attdouble = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304

305
306
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
307

308
  if ( atttype != NC_CHAR )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
    {
310
      double *pdoubleatt = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311

312
313
314
315
      if ( (int)nc_attlen > attlen )
        pdoubleatt = (double *) Malloc(nc_attlen * sizeof (double));
      else
        pdoubleatt = attdouble;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316

317
318
319
320
321
322
323
324
325
      cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);

      if ( (int)nc_attlen > attlen )
        {
          memcpy(attdouble, pdoubleatt, (size_t)attlen * sizeof (double));
          Free(pdoubleatt);
        }
    }
}
326

327
static
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
bool cdfCheckAttText(int fileID, int ncvarid, const char *attname)
{
  bool status = false;
  nc_type atttype;

  int status_nc = nc_inq_atttype(fileID, ncvarid, attname, &atttype);

  if ( status_nc == NC_NOERR
       && (atttype == NC_CHAR
#if  defined  (HAVE_NETCDF4)
           || atttype == NC_STRING
#endif
           ) )
    {
      status = true;
    }

  return status;
}

static
void cdfGetAttText(int fileID, int ncvarid, const char *attname, int attlen, char *atttext)
350
{
351
352
353
354
355
  nc_type atttype;
  size_t nc_attlen;

  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
356

357
  if ( atttype == NC_CHAR )
358
    {
359
360
361
362
      char attbuf[65636];
      if ( nc_attlen < sizeof(attbuf) )
        {
          cdf_get_att_text(fileID, ncvarid, attname, attbuf);
363

364
          if ( (int) nc_attlen > (attlen-1) ) nc_attlen = (size_t)(attlen-1);
365

366
367
368
369
          attbuf[nc_attlen++] = 0;
          memcpy(atttext, attbuf, nc_attlen);
        }
      else
370
        {
371
          atttext[0] = 0;
372
        }
373
    }
374
375
#if  defined  (HAVE_NETCDF4)
  else if ( atttype == NC_STRING )
376
    {
377
      if ( nc_attlen == 1 )
378
        {
379
380
          char *attbuf = NULL;
          cdf_get_att_string(fileID, ncvarid, attname, &attbuf);
381

382
          size_t ssize = strlen(attbuf) + 1;
383

384
385
386
387
388
389
390
391
          if ( ssize > (size_t)attlen ) ssize = (size_t)attlen;
          memcpy(atttext, attbuf, ssize);
          atttext[ssize - 1] = 0;
          Free(attbuf);
        }
      else
        {
          atttext[0] = 0;
392
        }
393
    }
394
395
#endif
}
396

397
398
399
400
static
bool xtypeIsText(nc_type xtype)
{
  bool isText = false;
401

402
403
404
405
  if      ( xtype == NC_CHAR ) isText = true;
#if  defined  (HAVE_NETCDF4)
  else if ( xtype == NC_STRING ) isText = true;
#endif
406

407
408
  return isText;
}
409

410
411
412
413
static
bool xtypeIsFloat(nc_type xtype)
{
  bool isFloat = xtype == NC_FLOAT || xtype == NC_DOUBLE;
414

415
  return isFloat;
416
417
}

418
419
420
421
422
423
424
425
426
427
428
429
430
static
bool xtypeIsInt(nc_type xtype)
{
  bool isInt = xtype == NC_SHORT || xtype == NC_INT
            || xtype == NC_BYTE
#if  defined  (HAVE_NETCDF4)
            || xtype == NC_USHORT || xtype == NC_UINT
            || xtype == NC_UBYTE
#endif
             ;

  return isInt;
}
431

432
static
433
int cdfInqDatatype(int xtype, bool lunsigned)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
{
435
  int datatype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436

437
438
439
#if  defined  (HAVE_NETCDF4)
  if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
440

441
442
443
444
445
446
447
448
449
450
451
452
453
454
  if      ( xtype == NC_BYTE   )  datatype = DATATYPE_INT8;
  /* else if ( xtype == NC_CHAR   )  datatype = DATATYPE_UINT8; */
  else if ( xtype == NC_SHORT  )  datatype = DATATYPE_INT16;
  else if ( xtype == NC_INT    )  datatype = DATATYPE_INT32;
  else if ( xtype == NC_FLOAT  )  datatype = DATATYPE_FLT32;
  else if ( xtype == NC_DOUBLE )  datatype = DATATYPE_FLT64;
#if  defined  (HAVE_NETCDF4)
  else if ( xtype == NC_UBYTE  )  datatype = DATATYPE_UINT8;
  else if ( xtype == NC_LONG   )  datatype = DATATYPE_INT32;
  else if ( xtype == NC_USHORT )  datatype = DATATYPE_UINT16;
  else if ( xtype == NC_UINT   )  datatype = DATATYPE_UINT32;
  else if ( xtype == NC_INT64  )  datatype = DATATYPE_FLT64;
  else if ( xtype == NC_UINT64 )  datatype = DATATYPE_FLT64;
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455

456
457
  return datatype;
}
458
459


460
461
462
463
/* not used
int cdfInqRecord(stream_t *streamptr, int *varID, int *levelID)
{
  int tsID, recID;
464

465
466
467
  recID = streamptr->tsteps[0].curRecID++;
  printf("cdfInqRecord recID %d %d\n", recID, streamptr->tsteps[0].curRecID);
  printf("cdfInqRecord tsID %d\n", streamptr->curTsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468

469
  if ( streamptr->tsteps[0].curRecID >= streamptr->tsteps[0].nrecs )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470
    {
471
      streamptr->tsteps[0].curRecID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
    }
473

474
475
476
477
478
479
480
481
482
483
  *varID   = streamptr->tsteps[0].records[recID].varID;
  *levelID = streamptr->tsteps[0].records[recID].levelID;

  streamptr->record->varID   = *varID;
  streamptr->record->levelID = *levelID;

  if ( CDI_Debug )
    Message("recID = %d  varID = %d  levelID = %d", recID, *varID, *levelID);

  return (recID+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484
}
485
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486

487
void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488
{
489
490
  int laddoffset   = IS_NOT_EQUAL(addoffset, 0);
  int lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
491

492
493
  if ( laddoffset || lscalefactor )
    {
494
      for ( size_t i = 0; i < size; ++i )
495
496
497
498
499
500
        {
          if ( lscalefactor ) data[i] *= scalefactor;
          if ( laddoffset )   data[i] += addoffset;
        }
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501

502
503
static
void cdfCreateRecords(stream_t *streamptr, int tsID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505
  if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
506

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

509
510
  int vlistID  = streamptr->vlistID;

511
512
513
514
515
516
517
518
  tsteps_t* sourceTstep = streamptr->tsteps;
  tsteps_t* destTstep = sourceTstep + tsID;

  int nvars = vlistNvars(vlistID);
  int nrecs = vlistNrecs(vlistID);

  if ( nrecs <= 0 ) return;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
520
  if ( tsID == 0 )
    {
521
      int nvrecs = nrecs; /* use all records at first timestep */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
522

Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
      streamptr->nrecs += nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524

525
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
526
527
528
529
      destTstep->nrecs      = nrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
530
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs*sizeof (int));;
531
      for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532

533
      record_t *records = destTstep->records;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534

535
      for ( int varID = 0, recID = 0; varID < nvars; varID++ )
536
        {
537
538
539
          int zaxisID = vlistInqVarZaxis(vlistID, varID);
          int nlev    = zaxisInqSize(zaxisID);
          for ( int levelID = 0; levelID < nlev; levelID++ )
540
541
            {
              recordInitEntry(&records[recID]);
542
543
              records[recID].varID   = (short)varID;
              records[recID].levelID = (short)levelID;
544
545
546
              recID++;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547
548
549
    }
  else if ( tsID == 1 )
    {
550
551
      int nvrecs = 0;
      for ( int varID = 0; varID < nvars; varID++ )
552
        {
553
          if ( vlistInqVarTsteptype(vlistID, varID) != TSTEP_CONSTANT )
554
            {
555
              int zaxisID = vlistInqVarZaxis(vlistID, varID);
556
557
558
              nvrecs += zaxisInqSize(zaxisID);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559

Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561

562
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
563
564
565
566
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567

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

570
      if ( nvrecs )
571
        {
572
          destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
573
          for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
574
            {
575
              int varID = destTstep->records[recID].varID;
576
              if ( vlistInqVarTsteptype(vlistID, varID) != TSTEP_CONSTANT )
577
                {
578
                  destTstep->recIDs[vrecID++] = recID;
579
580
581
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
582
583
584
    }
  else
    {
585
      if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);
586

587
      int nvrecs = streamptr->tsteps[1].nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588

Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
590

591
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
592
593
594
595
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
596

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

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

601
      memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
603
604
    }
}

605
static
606
int cdf_time_dimid(int fileID, int ndims, int nvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
{
608
609
  char dimname[80];
  for ( int dimid = 0; dimid < ndims; ++dimid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
610
    {
611
      dimname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612
      cdf_inq_dimname(fileID, dimid, dimname);
613
614
      size_t len = strlen(dimname);
      if ( len >= 4 && memcmp(dimname, "time", 4) == 0 ) return dimid;
615
      if ( len >= 4 && memcmp(dimname, "Time", 4) == 0 ) return dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
616
617
    }

618
  for ( int varid = 0; varid < nvars; ++varid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
    {
620
      nc_type xtype;
621
      int nvdims, nvatts, dimids[9];
622
      cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
623
      if ( nvdims == 1 )
624
        {
625
626
          char sbuf[CDI_MAX_NAME];
          for ( int iatt = 0; iatt < nvatts; ++iatt )
627
            {
628
              sbuf[0] = 0;
629
630
              cdf_inq_attname(fileID, varid, iatt, sbuf);
              if ( strncmp(sbuf, "units", 5) == 0 )
631
                {
632
                  cdfGetAttText(fileID, varid, "units", sizeof(sbuf), sbuf);
633
                  str_tolower(sbuf);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
634

635
                  if ( is_time_units(sbuf) ) return dimids[0];
636
637
638
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
639
640
    }

641
  return UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
643
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
static
645
void init_ncdims(long ndims, ncdim_t *ncdims)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
646
{
647
  for ( long ncdimid = 0; ncdimid < ndims; ncdimid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
648
649
650
651
652
653
654
655
656
    {
      ncdims[ncdimid].ncvarid      = UNDEFID;
      ncdims[ncdimid].dimtype      = UNDEFID;
      ncdims[ncdimid].len          = 0;
      ncdims[ncdimid].name[0]      = 0;
    }
}

static
657
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
658
{
659
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
660
    {
661
      ncvars[ncvarid].ncid            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
      ncvars[ncvarid].isvar           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
663
664
665
666
667
668
      ncvars[ncvarid].ignore          = false;
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
669
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
670
671
      ncvars[ncvarid].lformula        = false;
      ncvars[ncvarid].lformulaterms   = false;
672
      ncvars[ncvarid].tsteptype       = TSTEP_CONSTANT;
673
      ncvars[ncvarid].param           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
      ncvars[ncvarid].code            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
675
      ncvars[ncvarid].tabnum          = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
677
678
679
680
681
682
683
684
685
686
687
688
      ncvars[ncvarid].calendar        = FALSE;
      ncvars[ncvarid].bounds          = UNDEFID;
      ncvars[ncvarid].gridID          = UNDEFID;
      ncvars[ncvarid].zaxisID         = UNDEFID;
      ncvars[ncvarid].gridtype        = UNDEFID;
      ncvars[ncvarid].zaxistype       = UNDEFID;
      ncvars[ncvarid].xdim            = UNDEFID;
      ncvars[ncvarid].ydim            = UNDEFID;
      ncvars[ncvarid].zdim            = UNDEFID;
      ncvars[ncvarid].xvarid          = UNDEFID;
      ncvars[ncvarid].yvarid          = UNDEFID;
      ncvars[ncvarid].zvarid          = UNDEFID;
      ncvars[ncvarid].tvarid          = UNDEFID;
689
      ncvars[ncvarid].psvarid         = UNDEFID;
690
      ncvars[ncvarid].p0varid         = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
691
      ncvars[ncvarid].ncoordvars      = 0;
692
693
694
695
696
      for ( int i = 0; i < MAX_COORDVARS; ++i )
        ncvars[ncvarid].coordvarids[i]  = UNDEFID;
      ncvars[ncvarid].nauxvars      = 0;
      for ( int i = 0; i < MAX_AUXVARS; ++i )
        ncvars[ncvarid].auxvarids[i]  = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
697
698
699
700
701
      ncvars[ncvarid].cellarea        = UNDEFID;
      ncvars[ncvarid].tableID         = UNDEFID;
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
      ncvars[ncvarid].gmapid          = UNDEFID;
702
703
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
704
      ncvars[ncvarid].truncation      = 0;
705
      ncvars[ncvarid].position        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706
      ncvars[ncvarid].positive        = 0;
707
708
      ncvars[ncvarid].chunked         = 0;
      ncvars[ncvarid].chunktype       = UNDEFID;
709
710
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
711
      ncvars[ncvarid].missval         = 0;
712
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
713
714
715
716
717
718
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].name[0]         = 0;
      ncvars[ncvarid].longname[0]     = 0;
      ncvars[ncvarid].stdname[0]      = 0;
      ncvars[ncvarid].units[0]        = 0;
719
      ncvars[ncvarid].extra[0]        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
720
721
722
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
723
724
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
725
726
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
727
      ncvars[ncvarid].ensdata         = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
729
730
    }
}

731
static
732
void cdf_set_var(ncvar_t *ncvars, int ncvarid, short isvar)
733
734
735
{
  if ( ncvars[ncvarid].isvar != UNDEFID &&
       ncvars[ncvarid].isvar != isvar   &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
736
       ncvars[ncvarid].warn  == false )
737
738
739
740
    {
      if ( ! ncvars[ncvarid].ignore )
        Warning("Inconsistent variable definition for %s!", ncvars[ncvarid].name);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
      ncvars[ncvarid].warn = true;
742
743
744
745
746
747
748
      isvar = FALSE;
    }

  ncvars[ncvarid].isvar = isvar;
}

static
749
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
750
751
752
753
754
755
756
757
758
759
760
761
{
  if ( ncvars[ncvarid].dimtype[dimid] != UNDEFID &&
       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
762
int scan_hybrid_formula(int ncid, int ncfvarid, int *apvarid, int *bvarid, int *psvarid, int *avarid, int *p0varid)
763
{
764
  int status = 0;
765
766
767
  *apvarid = -1;
  *bvarid  = -1;
  *psvarid = -1;
768
769
  *avarid  = -1;
  *p0varid = -1;
Thomas Jahns's avatar
Thomas Jahns committed
770
  enum { attstringlen = 8192 }; char attstring[attstringlen];
771
772
773
  cdfGetAttText(ncid, ncfvarid, "formula", attstringlen, attstring);
  if ( strcmp(attstring, "p = ap + b*ps") == 0 )
    {
774
      status = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
775
      bool lstop = false;
776
777
778
779
780
781
782
783
784
785
      int dimvarid;
      cdfGetAttText(ncid, ncfvarid, "formula_terms", attstringlen, attstring);
      char *pstring = attstring;

      for ( int i = 0; i < 3; i++ )
        {
          while ( isspace((int) *pstring) ) pstring++;
          if ( *pstring == 0 ) break;
          char *tagname = pstring;
          while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
786
          if ( *pstring == 0 ) lstop = true;
787
788
789
790
791
792
          *pstring++ = 0;

          while ( isspace((int) *pstring) ) pstring++;
          if ( *pstring == 0 ) break;
          char *varname = pstring;
          while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793
          if ( *pstring == 0 ) lstop = true;
794
795
          *pstring++ = 0;

Thomas Jahns's avatar
Thomas Jahns committed
796
797
          int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
          if ( status_nc == NC_NOERR )
798
799
800
801
802
            {
              if      ( strcmp(tagname, "ap:") == 0 ) *apvarid = dimvarid;
              else if ( strcmp(tagname, "b:")  == 0 ) *bvarid  = dimvarid;
              else if ( strcmp(tagname, "ps:") == 0 ) *psvarid = dimvarid;
            }
803
          else if ( strcmp(tagname, "ps:") != 0 )
804
            {
Thomas Jahns's avatar
Thomas Jahns committed
805
              Warning("%s - %s", nc_strerror(status_nc), varname);
806
807
808
809
810
            }

          if ( lstop ) break;
        }
    }
811
812
813
  else if ( strcmp(attstring, "xxxp = a*p0 + b*ps") == 0 )
    {
      status = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
814
      bool lstop = false;
815
816
817
818
819
820
821
822
823
824
      int dimvarid;
      cdfGetAttText(ncid, ncfvarid, "formula_terms", attstringlen, attstring);
      char *pstring = attstring;

      for ( int i = 0; i < 4; i++ )
        {
          while ( isspace((int) *pstring) ) pstring++;
          if ( *pstring == 0 ) break;
          char *tagname = pstring;
          while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
825
          if ( *pstring == 0 ) lstop = true;
826
827
828
829
830
831
          *pstring++ = 0;

          while ( isspace((int) *pstring) ) pstring++;
          if ( *pstring == 0 ) break;
          char *varname = pstring;
          while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
832
          if ( *pstring == 0 ) lstop = true;
833
834
          *pstring++ = 0;

Thomas Jahns's avatar
Thomas Jahns committed
835
836
          int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
          if ( status_nc == NC_NOERR )
837
838
839
840
841
842
843
844
            {
              if      ( strcmp(tagname, "a:")  == 0 ) *avarid  = dimvarid;
              else if ( strcmp(tagname, "b:")  == 0 ) *bvarid  = dimvarid;
              else if ( strcmp(tagname, "ps:") == 0 ) *psvarid = dimvarid;
              else if ( strcmp(tagname, "p0:") == 0 ) *p0varid = dimvarid;
            }
          else if ( strcmp(tagname, "ps:") != 0 )
            {
Thomas Jahns's avatar
Thomas Jahns committed
845
              Warning("%s - %s", nc_strerror(status_nc), varname);
846
847
848
849
850
851
852
            }

          if ( lstop ) break;
        }
    }

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
853
854
}

855
static
856
bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
857
{
858
  bool status = false;
859
860
861
862
863
  int ncfvarid = ncvarid;
  ncvar_t *ncvar = &ncvars[ncvarid];

  if ( strcmp(ncvar->stdname, "atmosphere_hybrid_sigma_pressure_coordinate") == 0 )
    {
864
865
      cdiConvention = CDI_CONVENTION_CF;

866
      status = true;
867
868
      ncvar->zaxistype = ZAXIS_HYBRID;
      int dimid = ncvar->dimids[0];
Thomas Jahns's avatar
Thomas Jahns committed
869
      size_t dimlen = ncdims[dimid].len;
870

871
872
      int ret;
      int apvarid1 = -1, bvarid1 = -1, psvarid1 = -1, avarid1 = -1, p0varid1 = -1;
873
      if ( ncvars[ncfvarid].lformula && ncvars[ncfvarid].lformulaterms )
874
        ret = scan_hybrid_formula(ncid, ncfvarid, &apvarid1, &bvarid1, &psvarid1, &avarid1, &p0varid1);
875
876
877
      if ( apvarid1 != -1 ) ncvars[apvarid1].isvar = FALSE;
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar  = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
878
879
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
880
881
882
883

      if ( ncvar->bounds != UNDEFID && ncvars[ncvar->bounds].lformula && ncvars[ncvar->bounds].lformulaterms )
        {
          ncfvarid = ncvar->bounds;
884
885
          int apvarid2 = -1, bvarid2 = -1, psvarid2 = -1, avarid2 = -1, p0varid2 = -1;
          ret = 0;
886
          if ( ncvars[ncfvarid].lformula && ncvars[ncfvarid].lformulaterms )
887
888
889
            ret = scan_hybrid_formula(ncid, ncfvarid, &apvarid2, &bvarid2, &psvarid2, &avarid2, &p0varid2);
          if ( ret == 1 ) avarid2 = apvarid2;
          if ( avarid2 != -1 && bvarid2 != -1 )
890
            {
891
              ncvars[avarid2].isvar = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
892
              ncvars[bvarid2].isvar = FALSE;
893

894
              if ( dimid == ncvars[avarid2].dimids[0] && ncdims[ncvars[avarid2].dimids[1]].len == 2 )
895
                {
896
897
898
899
                  double px = 1;
                  if ( ret == 2 && p0varid1 == p0varid2 )
                    cdf_get_var_double(ncid, p0varid2, &px);

900
                  double abuf[dimlen*2], bbuf[dimlen*2];
901
                  cdf_get_var_double(ncid, avarid2, abuf);
902
903
904
905
906
                  cdf_get_var_double(ncid, bvarid2, bbuf);
                  /*
                  for ( int i = 0; i < dimlen; ++i )
                    printf("%d  %g %g    %g %g\n", i, abuf[i*2], abuf[i*2+1], bbuf[i*2], bbuf[i*2+1]);
                  */
Thomas Jahns's avatar
Thomas Jahns committed
907
                  size_t vctsize = (dimlen+1)*2;
908
                  double *vct = (double *) Malloc(vctsize * sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
909
                  for ( size_t i = 0; i < dimlen; ++i )
910
911
912
913
914
                    {
                      vct[i] = abuf[i*2];
                      vct[i+dimlen+1] = bbuf[i*2];
                    }
                  vct[dimlen]     = abuf[dimlen*2-1];
915
                  vct[dimlen*2+1] = bbuf[dimlen*2-1];
916

917
918
919
                  if ( ret == 2 && IS_NOT_EQUAL(px, 1) )
                    for ( size_t i = 0; i < dimlen+1; ++i ) vct[i] *= px;

920
921
922
923
924
925
926
927
928
929
                  ncvar->vct = vct;
                  ncvar->vctsize = vctsize;
                }
            }
        }
    }

  return status;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
930
static
931
void cdf_print_vars(const ncvar_t *ncvars, int nvars, const char *oname)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
932
{
933
  char axis[7];
934
  static const char iaxis[] = {'t', 'z', 'y', 'x'};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
935

936
937
  fprintf(stderr, "%s:\n", oname);

938
  for ( int ncvarid = 0; ncvarid < nvars; ncvarid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
939
    {
940
      int ndim = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
941
      if ( ncvars[ncvarid].isvar )
942
943
        {
          axis[ndim++] = 'v';
944
          axis[ndim++] = ':';
945
          for ( int i = 0; i < ncvars[ncvarid].ndims; i++ )
946
947
948
949
950
951
952
            {/*
              if      ( ncvars[ncvarid].tvarid != -1 ) axis[ndim++] = iaxis[0];
              else if ( ncvars[ncvarid].zvarid != -1 ) axis[ndim++] = iaxis[1];
              else if ( ncvars[ncvarid].yvarid != -1 ) axis[ndim++] = iaxis[2];
              else if ( ncvars[ncvarid].xvarid != -1 ) axis[ndim++] = iaxis[3];
              else
             */
953
954
955
956
957
              if      ( ncvars[ncvarid].dimtype[i] == T_AXIS ) axis[ndim++