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

115
116

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

  *rdate = 0;
  *rtime = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

269
270
  return 0;
}
271

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

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

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

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

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

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

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

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

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

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

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

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

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

328
static
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
350
void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
351
{
352
353
354
355
356
  nc_type atttype;
  size_t nc_attlen;

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

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

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

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

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

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

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

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

408
409
  return isText;
}
410

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

416
  return isFloat;
417
418
}

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

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

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

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

457
458
  return datatype;
}
459
460


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

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

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

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

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

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

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

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

510
511
  int vlistID  = streamptr->vlistID;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

static
656
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
{
658
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
659
    {
660
      ncvars[ncvarid].ncid            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
      ncvars[ncvarid].isvar           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662
      ncvars[ncvarid].ignore          = false;
663
664
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
666
667
668
669
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
670
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
      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
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
718
719
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
720
721
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
722
      ncvars[ncvarid].ensdata         = NULL;
723
724
725
726
727
      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
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
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
int scan_hybrid_formulaterms(int ncid, int ncfvarid, int *apvarid, int *bvarid, int *psvarid, int *avarid, int *p0varid)
763
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
  int status = 1;
765
766
767
  *apvarid = -1;
  *bvarid  = -1;
  *psvarid = -1;
768
769
  *avarid  = -1;
  *p0varid = -1;
770
  char attstring[1024];
771

Uwe Schulzweida's avatar
Uwe Schulzweida committed
772
773
774
775
  bool lstop = false;
  int dimvarid;
  cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
  char *pstring = attstring;
776

Uwe Schulzweida's avatar
Uwe Schulzweida committed
777
  for ( int i = 0; i < 4; i++ )
778
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
      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;

      int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
      if ( status_nc == NC_NOERR )
795
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
796
797
798
799
          if      ( strcmp(tagname, "ap:") == 0 ) *apvarid = 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;
800
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801
802
803
804
805
806
      else if ( strcmp(tagname, "ps:") != 0 )
        {
          Warning("%s - %s", nc_strerror(status_nc), varname);
        }

      if ( lstop ) break;
807
808
809
    }

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810
811
}

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

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

822
      status = true;
823
824
      ncvar->zaxistype = ZAXIS_HYBRID;
      int dimid = ncvar->dimids[0];
Thomas Jahns's avatar
Thomas Jahns committed
825
      size_t dimlen = ncdims[dimid].len;
826

Uwe Schulzweida's avatar
Uwe Schulzweida committed
827
      int ret = 0;
828
      int apvarid1 = -1, bvarid1 = -1, psvarid1 = -1, avarid1 = -1, p0varid1 = -1;
829
      int ncfvarid = ncvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
831
832
      if ( ncvars[ncfvarid].lformulaterms )
        ret = scan_hybrid_formulaterms(ncid, ncfvarid, &apvarid1, &bvarid1, &psvarid1, &avarid1, &p0varid1);
      // printf("ret %d apvarid1, bvarid1, psvarid1, avarid1, p0varid1 %d %d %d %d %d\n", ret, apvarid1, bvarid1, psvarid1, avarid1, p0varid1);
833
834
835
      if ( apvarid1 != -1 ) ncvars[apvarid1].isvar = FALSE;
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar  = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
836
837
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
838

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

853
              if ( dimid == ncvars[avarid2].dimids[0] && ncdims[ncvars[avarid2].dimids[1]].len == 2 )
854
                {
855
856
857
858
                  double px = 1;
                  if ( ret == 2 && p0varid1 == p0varid2 )
                    cdf_get_var_double(ncid, p0varid2, &px);

859
                  double abuf[dimlen*2], bbuf[dimlen*2];
860
                  cdf_get_var_double(ncid, avarid2, abuf);
861
862
863
864
865
                  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
866
                  size_t vctsize = (dimlen+1)*2;
867
                  double *vct = (double *) Malloc(vctsize * sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
868
                  for ( size_t i = 0; i < dimlen; ++i )
869
870
871
872
873
                    {
                      vct[i] = abuf[i*2];
                      vct[i+dimlen+1] = bbuf[i*2];
                    }
                  vct[dimlen]     = abuf[dimlen*2-1];
874
                  vct[dimlen*2+1] = bbuf[dimlen*2-1];
875

876
877
878
                  if ( ret == 2 && IS_NOT_EQUAL(px, 1) )
                    for ( size_t i = 0; i < dimlen+1; ++i ) vct[i] *= px;

879
880
881
882
883
884
885
886
887
888
                  ncvar->vct = vct;
                  ncvar->vctsize = vctsize;
                }
            }
        }
    }

  return status;
}

889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
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 :
905
                     (atttype == NC_BYTE)   ? DATATYPE_INT8 :
906
907
908
909
910
911
#if  defined  (HAVE_NETCDF4)
                     (atttype == NC_UBYTE)  ? DATATYPE_UINT8 :
                     (atttype == NC_USHORT) ? DATATYPE_UINT16 :
                     (atttype == NC_UINT)   ? DATATYPE_UINT32 :
#endif
                     DATATYPE_INT32;
912
      cdiDefAttInt(cdiID, varID, attname, datatype, (int)attlen, attint);
913
914
915
916
917
918
    }
  else if ( xtypeIsFloat(atttype) )
    {
      double attflt[attlen];
      cdfGetAttDouble(ncid, ncvarid, attname, (int)attlen, attflt);
      int datatype = (atttype == NC_FLOAT) ? DATATYPE_FLT32 : DATATYPE_FLT64;
919
      cdiDefAttFlt(cdiID, varID, attname, datatype, (int)attlen, attflt);
920
921
922
    }
  else if ( xtypeIsText(atttype) )
    {
923
924
      char attstring[8192];
      cdfGetAttText(ncid, ncvarid, attname, sizeof(attstring), attstring);
925
      cdiDefAttTxt(cdiID, varID, attname, (int)attlen, attstring);
926
927
928
    }
}

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

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

937
  for ( int ncvarid = 0; ncvarid < nvars; ncvarid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
938
    {
939
      int ndim = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
940
      if ( ncvars[ncvarid].isvar )
941
942
        {
          axis[ndim++] = 'v';
943
          axis[ndim++] = ':';