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

116
117

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

  *rdate = 0;
  *rtime = 0;

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
  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++;
            }
        }
    }
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

  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
174
175
176
static
int scanTimeUnit(const char *unitstr)
{
177
  size_t len = strlen(unitstr);
178
  int timeunit = get_timeunit(len, unitstr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
180
181
  if ( timeunit == -1 )
    Message("Unsupported TIMEUNIT: %s!", unitstr);

182
  return timeunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
184
}

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

191
  int len = (int) strlen(timestr);
192
193
194
195
196
197
198
199
200
201
202
203
  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
204
205
206
  int timetype = TAXIS_ABSOLUTE;
  int rdate = -1, rtime = -1;

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

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

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

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

226
      if ( str_is_equal(ptu, "as") )
227
        timetype = TAXIS_ABSOLUTE;
228
      else if ( str_is_equal(ptu, "since") )
229
        timetype = TAXIS_RELATIVE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
231
232

      while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
      if ( *ptu )
233
234
235
236
237
        {
          while ( isspace(*ptu) ) ptu++;

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

253
254
              (*taxis).rdate = rdate;
              (*taxis).rtime = rtime;
255

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

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

265
  Free(tu);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266

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

270
271
  return 0;
}
272

273
274
275
276
277
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
278

279
  *attint = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280

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

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

289
      cdf_get_att_int(fileID, ncvarid, attname, pintatt);
290

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

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

305
  *attdouble = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306

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

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

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

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

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

329
static
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
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
351
void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
352
{
353
354
355
356
357
  nc_type atttype;
  size_t nc_attlen;

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

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

366
          if ( nc_attlen > (attlen-1) ) nc_attlen = (attlen-1);
367

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

384
          size_t ssize = strlen(attbuf) + 1;
385

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

399
400
401
402
static
bool xtypeIsText(nc_type xtype)
{
  bool isText = false;
403

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

409
410
  return isText;
}
411

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

417
  return isFloat;
418
419
}

420
421
422
423
424
425
426
427
428
429
430
431
432
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;
}
433

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

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

443
444
445
446
447
448
449
450
451
452
453
454
455
456
  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
457

458
459
  return datatype;
}
460
461


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

467
468
469
  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
470

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

476
477
478
479
480
481
482
483
484
485
  *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
486
}
487
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488

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

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

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

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

511
512
  int vlistID  = streamptr->vlistID;

513
514
515
516
517
518
519
520
  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
521
522
  if ( tsID == 0 )
    {
523
      int nvrecs = nrecs; /* use all records at first timestep */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524

Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
      streamptr->nrecs += nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
526

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

535
      record_t *records = destTstep->records;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563

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

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

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

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592

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

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

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

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

607
static
608
int cdf_time_dimid(int fileID, int ndims, int nvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
609
{
610
611
  char dimname[80];
  for ( int dimid = 0; dimid < ndims; ++dimid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612
    {
613
      dimname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614
      cdf_inq_dimname(fileID, dimid, dimname);
615
      if ( str_is_equal(dimname, "time") || str_is_equal(dimname, "Time") ) 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
      ncvars[ncvarid].ignore          = false;
664
665
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
666
667
668
669
670
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
671
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
      ncvars[ncvarid].lformula        = false;
      ncvars[ncvarid].lformulaterms   = false;
674
      ncvars[ncvarid].tsteptype       = TSTEP_CONSTANT;
675
      ncvars[ncvarid].param           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
      ncvars[ncvarid].code            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677
      ncvars[ncvarid].tabnum          = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
678
679
680
681
682
683
684
685
686
687
688
689
690
      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;
691
      ncvars[ncvarid].psvarid         = UNDEFID;
692
      ncvars[ncvarid].p0varid         = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693
      ncvars[ncvarid].ncoordvars      = 0;
694
695
696
697
698
      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
699
700
701
702
703
      ncvars[ncvarid].cellarea        = UNDEFID;
      ncvars[ncvarid].tableID         = UNDEFID;
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
      ncvars[ncvarid].gmapid          = UNDEFID;
704
705
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706
      ncvars[ncvarid].truncation      = 0;
707
      ncvars[ncvarid].position        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
708
      ncvars[ncvarid].positive        = 0;
709
710
      ncvars[ncvarid].chunked         = 0;
      ncvars[ncvarid].chunktype       = UNDEFID;
711
712
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
713
      ncvars[ncvarid].missval         = 0;
714
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
715
716
717
718
719
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
720
721
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
722
723
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
724
      ncvars[ncvarid].ensdata         = NULL;
725
726
727
728
729
      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
730
731
732
    }
}

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

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

  ncvars[ncvarid].isvar = isvar;
}

static
751
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
752
753
754
755
756
757
758
759
760
761
762
763
{
  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
764
int scan_hybrid_formula(int ncid, int ncfvarid, int *apvarid, int *bvarid, int *psvarid, int *avarid, int *p0varid)
765
{
766
  int status = 0;
767
768
769
  *apvarid = -1;
  *bvarid  = -1;
  *psvarid = -1;
770
771
  *avarid  = -1;
  *p0varid = -1;
772
773
  char attstring[1024];
  cdfGetAttText(ncid, ncfvarid, "formula", sizeof(attstring), attstring);
774
775
  if ( strcmp(attstring, "p = ap + b*ps") == 0 )
    {
776
      status = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
777
      bool lstop = false;
778
      int dimvarid;
779
      cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
780
781
782
783
784
785
786
787
      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
788
          if ( *pstring == 0 ) lstop = true;
789
790
791
792
793
794
          *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
795
          if ( *pstring == 0 ) lstop = true;
796
797
          *pstring++ = 0;

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

          if ( lstop ) break;
        }
    }
813
814
815
  else if ( strcmp(attstring, "xxxp = a*p0 + b*ps") == 0 )
    {
      status = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
      bool lstop = false;
817
      int dimvarid;
818
      cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
819
820
821
822
823
824
825
826
      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
827
          if ( *pstring == 0 ) lstop = true;
828
829
830
831
832
833
          *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
834
          if ( *pstring == 0 ) lstop = true;
835
836
          *pstring++ = 0;

Thomas Jahns's avatar
Thomas Jahns committed
837
838
          int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
          if ( status_nc == NC_NOERR )
839
840
841
842
843
844
845
846
            {
              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
847
              Warning("%s - %s", nc_strerror(status_nc), varname);
848
849
850
851
852
853
854
            }

          if ( lstop ) break;
        }
    }

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855
856
}

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

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

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

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

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

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

902
                  double abuf[dimlen*2], bbuf[dimlen*2];
903
                  cdf_get_var_double(ncid, avarid2, abuf);
904
905
906
907
908
                  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
909
                  size_t vctsize = (dimlen+1)*2;
910
                  double *vct = (double *) Malloc(vctsize * sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
911
                  for ( size_t i = 0; i < dimlen; ++i )
912
913
914
915
916
                    {
                      vct[i] = abuf[i*2];
                      vct[i+dimlen+1] = bbuf[i*2];
                    }
                  vct[dimlen]     = abuf[dimlen*2-1];
917
                  vct[dimlen*2+1] = bbuf[dimlen*2-1];
918

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

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

  return status;
}

932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
static
void cdf_set_cdi_attr(int ncid, int ncvarid, int attnum, int cdiID, int varID)
{
  nc_type atttype;
  size_t attlen;
  char attname[CDI_MAX_NAME];

  cdf_inq_attname(ncid, ncvarid, attnum, attname);
  cdf_inq_attlen(ncid, ncvarid, attname, &attlen);
  cdf_inq_atttype(ncid, ncvarid, attname, &atttype);

  if ( xtypeIsInt(atttype) )
    {
      int attint[attlen];
      cdfGetAttInt(ncid, ncvarid, attname, (int)attlen, attint);
      int datatype = (atttype == NC_SHORT)  ? DATATYPE_INT16 :
948
                     (atttype == NC_BYTE)   ? DATATYPE_INT8 :
949
950
951
952
953
954
#if  defined  (HAVE_NETCDF4)
                     (atttype == NC_UBYTE)  ? DATATYPE_UINT8 :
                     (atttype == NC_USHORT) ? DATATYPE_UINT16 :
                     (atttype == NC_UINT)   ? DATATYPE_UINT32 :
#endif
                     DATATYPE_INT32;
955
      cdiDefAttInt(cdiID, varID, attname, datatype, (int)attlen, attint);