stream_cdf_i.c 132 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 ( memcmp(ptu, "as", 2) == 0 )
227
        timetype = TAXIS_ABSOLUTE;
228
      else if ( memcmp(ptu, "since", 5) == 0 )
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
238
239
240
241
242
243
        {
          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 )
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
351
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)
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 ( (int) nc_attlen > (attlen-1) ) nc_attlen = (size_t)(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
387
388
389
390
391
392
393
          if ( ssize > (size_t)attlen ) ssize = (size_t)attlen;
          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
616
      size_t len = strlen(dimname);
      if ( len >= 4 && memcmp(dimname, "time", 4) == 0 ) return dimid;
617
      if ( len >= 4 && memcmp(dimname, "Time", 4) == 0 ) return dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
619
    }

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

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

643
  return UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
645
}

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

static
659
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
660
{
661
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
662
    {
663
      ncvars[ncvarid].ncid            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
664
      ncvars[ncvarid].isvar           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
      ncvars[ncvarid].ignore          = false;
666
667
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
668
669
670
671
672
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
673
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
675
      ncvars[ncvarid].lformula        = false;
      ncvars[ncvarid].lformulaterms   = false;
676
      ncvars[ncvarid].tsteptype       = TSTEP_CONSTANT;
677
      ncvars[ncvarid].param           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
678
      ncvars[ncvarid].code            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
679
      ncvars[ncvarid].tabnum          = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
680
681
682
683
684
685
686
687
688
689
690
691
692
      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;
693
      ncvars[ncvarid].psvarid         = UNDEFID;
694
      ncvars[ncvarid].p0varid         = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
695
      ncvars[ncvarid].ncoordvars      = 0;
696
697
698
699
700
      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
701
702
703
704
705
      ncvars[ncvarid].cellarea        = UNDEFID;
      ncvars[ncvarid].tableID         = UNDEFID;
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
      ncvars[ncvarid].gmapid          = UNDEFID;
706
707
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
708
      ncvars[ncvarid].truncation      = 0;
709
      ncvars[ncvarid].position        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
710
      ncvars[ncvarid].positive        = 0;
711
712
      ncvars[ncvarid].chunked         = 0;
      ncvars[ncvarid].chunktype       = UNDEFID;
713
714
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
715
      ncvars[ncvarid].missval         = 0;
716
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
717
718
719
720
721
722
      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;
723
      ncvars[ncvarid].extra[0]        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
724
725
726
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
727
728
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
729
730
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
731
      ncvars[ncvarid].ensdata         = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
732
733
734
    }
}

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

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

  ncvars[ncvarid].isvar = isvar;
}

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

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

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

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

          if ( lstop ) break;
        }
    }

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857
858
}

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

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

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

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

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

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

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

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

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

  return status;
}

934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
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 :
950
                     (atttype == NC_BYTE)   ? DATATYPE_INT8 :
951
952
953
954
955
956
#if  defined  (HAVE_NETCDF4)
                     (atttype == NC_UBYTE)  ? DATATYPE_UINT8 :
                     (atttype == NC_USHORT) ? DATATYPE_UINT16 :
                     (atttype == NC_UINT)   ? DATATYPE_UINT32 :
#endif
                     DATATYPE_INT32;
957
      cdiDefAttInt(cdiID, varID, attname, datatype, (int)attlen, attint);
958
959
960
961
962
963
    }
  else if ( xtypeIsFloat(atttype) )
    {
      double attflt[attlen];
      cdfGetAttDouble(ncid, ncvarid, attname, (int)attlen, attflt);
      int datatype = (atttype == NC_FLOAT) ? DATATYPE_FLT32 : DATATYPE_FLT64;
964
      cdiDefAttFlt(cdiID, varID, attname, datatype, (int)attlen, attflt);