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
//#define PROJECTION_TEST

Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
25
26
27
28
29
30
31
#undef  UNDEFID
#define UNDEFID  CDI_UNDEFID

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

32
33
34
#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
typedef struct {
46
  int      ncid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
48
49
50
51
52
53
  int      isvar;
  bool     ignore;
  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
void cdfGetAttText(int fileID, int ncvarid,const char *attname, int attlen, char *atttext)
331
{
332
333
334
335
336
  nc_type atttype;
  size_t nc_attlen;

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

338
  if ( atttype == NC_CHAR )
339
    {
340
341
342
343
      char attbuf[65636];
      if ( nc_attlen < sizeof(attbuf) )
        {
          cdf_get_att_text(fileID, ncvarid, attname, attbuf);
344

345
          if ( (int) nc_attlen > (attlen-1) ) nc_attlen = (size_t)(attlen-1);
346

347
348
349
350
          attbuf[nc_attlen++] = 0;
          memcpy(atttext, attbuf, nc_attlen);
        }
      else
351
        {
352
          atttext[0] = 0;
353
        }
354
    }
355
356
#if  defined  (HAVE_NETCDF4)
  else if ( atttype == NC_STRING )
357
    {
358
      if ( nc_attlen == 1 )
359
        {
360
361
          char *attbuf = NULL;
          cdf_get_att_string(fileID, ncvarid, attname, &attbuf);
362

363
          size_t ssize = strlen(attbuf) + 1;
364

365
366
367
368
369
370
371
372
          if ( ssize > (size_t)attlen ) ssize = (size_t)attlen;
          memcpy(atttext, attbuf, ssize);
          atttext[ssize - 1] = 0;
          Free(attbuf);
        }
      else
        {
          atttext[0] = 0;
373
        }
374
    }
375
376
#endif
}
377

378
379
380
381
static
bool xtypeIsText(nc_type xtype)
{
  bool isText = false;
382

383
384
385
386
  if      ( xtype == NC_CHAR ) isText = true;
#if  defined  (HAVE_NETCDF4)
  else if ( xtype == NC_STRING ) isText = true;
#endif
387

388
389
  return isText;
}
390

391
392
393
394
static
bool xtypeIsFloat(nc_type xtype)
{
  bool isFloat = xtype == NC_FLOAT || xtype == NC_DOUBLE;
395

396
  return isFloat;
397
398
}

399
400
401
402
403
404
405
406
407
408
409
410
411
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;
}
412

413
static
414
int cdfInqDatatype(int xtype, bool lunsigned)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
{
416
  int datatype = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417

418
419
420
#if  defined  (HAVE_NETCDF4)
  if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421

422
423
424
425
426
427
428
429
430
431
432
433
434
435
  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
436

437
438
  return datatype;
}
439
440


441
442
443
444
/* not used
int cdfInqRecord(stream_t *streamptr, int *varID, int *levelID)
{
  int tsID, recID;
445

446
447
448
  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
449

450
  if ( streamptr->tsteps[0].curRecID >= streamptr->tsteps[0].nrecs )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
451
    {
452
      streamptr->tsteps[0].curRecID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
    }
454

455
456
457
458
459
460
461
462
463
464
  *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
465
}
466
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467

468
void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
{
470
471
  int laddoffset   = IS_NOT_EQUAL(addoffset, 0);
  int lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
472

473
474
  if ( laddoffset || lscalefactor )
    {
475
      for ( size_t i = 0; i < size; ++i )
476
477
478
479
480
481
        {
          if ( lscalefactor ) data[i] *= scalefactor;
          if ( laddoffset )   data[i] += addoffset;
        }
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482

483
484
static
void cdfCreateRecords(stream_t *streamptr, int tsID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
  if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
487

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

490
491
  int vlistID  = streamptr->vlistID;

492
493
494
495
496
497
498
499
  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
500
501
  if ( tsID == 0 )
    {
502
      int nvrecs = nrecs; /* use all records at first timestep */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
503

Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
      streamptr->nrecs += nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505

506
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
507
508
509
510
      destTstep->nrecs      = nrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
511
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs*sizeof (int));;
512
      for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513

514
      record_t *records = destTstep->records;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
515

516
      for ( int varID = 0, recID = 0; varID < nvars; varID++ )
517
        {
518
519
520
          int zaxisID = vlistInqVarZaxis(vlistID, varID);
          int nlev    = zaxisInqSize(zaxisID);
          for ( int levelID = 0; levelID < nlev; levelID++ )
521
522
            {
              recordInitEntry(&records[recID]);
523
524
              records[recID].varID   = (short)varID;
              records[recID].levelID = (short)levelID;
525
526
527
              recID++;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
528
529
530
    }
  else if ( tsID == 1 )
    {
531
532
      int nvrecs = 0;
      for ( int varID = 0; varID < nvars; varID++ )
533
        {
534
          if ( vlistInqVarTsteptype(vlistID, varID) != TSTEP_CONSTANT )
535
            {
536
              int zaxisID = vlistInqVarZaxis(vlistID, varID);
537
538
539
              nvrecs += zaxisInqSize(zaxisID);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540

Uwe Schulzweida's avatar
Uwe Schulzweida committed
541
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542

543
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
544
545
546
547
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
548

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

551
      if ( nvrecs )
552
        {
553
          destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
554
          for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
555
            {
556
              int varID = destTstep->records[recID].varID;
557
              if ( vlistInqVarTsteptype(vlistID, varID) != TSTEP_CONSTANT )
558
                {
559
                  destTstep->recIDs[vrecID++] = recID;
560
561
562
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
564
565
    }
  else
    {
566
      if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);
567

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
570
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571

572
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
573
574
575
576
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577

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

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

582
      memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
583
584
585
    }
}

586
static
587
int cdf_time_dimid(int fileID, int ndims, int nvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
{
589
590
  char dimname[80];
  for ( int dimid = 0; dimid < ndims; ++dimid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
    {
592
      dimname[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
      cdf_inq_dimname(fileID, dimid, dimname);
594
595
      size_t len = strlen(dimname);
      if ( len >= 4 && memcmp(dimname, "time", 4) == 0 ) return dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
596
597
    }

598
  for ( int varid = 0; varid < nvars; ++varid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
599
    {
600
      nc_type xtype;
601
      int nvdims, nvatts, dimids[9];
602
      cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
603
      if ( nvdims == 1 )
604
        {
605
606
          char sbuf[CDI_MAX_NAME];
          for ( int iatt = 0; iatt < nvatts; ++iatt )
607
            {
608
              sbuf[0] = 0;
609
610
              cdf_inq_attname(fileID, varid, iatt, sbuf);
              if ( strncmp(sbuf, "units", 5) == 0 )
611
                {
612
                  cdfGetAttText(fileID, varid, "units", sizeof(sbuf), sbuf);
613
                  str_tolower(sbuf);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
614

615
                  if ( is_time_units(sbuf) ) return dimids[0];
616
617
618
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
620
    }

621
  return UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622
623
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
624
static
625
void init_ncdims(long ndims, ncdim_t *ncdims)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
{
627
  for ( long ncdimid = 0; ncdimid < ndims; ncdimid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
628
629
630
631
632
633
634
635
636
    {
      ncdims[ncdimid].ncvarid      = UNDEFID;
      ncdims[ncdimid].dimtype      = UNDEFID;
      ncdims[ncdimid].len          = 0;
      ncdims[ncdimid].name[0]      = 0;
    }
}

static
637
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
638
{
639
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
640
    {
641
      ncvars[ncvarid].ncid            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
      ncvars[ncvarid].isvar           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643
644
645
646
647
648
      ncvars[ncvarid].ignore          = false;
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
649
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
651
      ncvars[ncvarid].lformula        = false;
      ncvars[ncvarid].lformulaterms   = false;
652
      ncvars[ncvarid].tsteptype       = TSTEP_CONSTANT;
653
      ncvars[ncvarid].param           = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
      ncvars[ncvarid].code            = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
655
      ncvars[ncvarid].tabnum          = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
656
657
658
659
660
661
662
663
664
665
666
667
668
      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;
669
      ncvars[ncvarid].psvarid         = UNDEFID;
670
      ncvars[ncvarid].p0varid         = UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
      ncvars[ncvarid].ncoordvars      = 0;
672
673
674
675
676
      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
677
678
679
680
681
      ncvars[ncvarid].cellarea        = UNDEFID;
      ncvars[ncvarid].tableID         = UNDEFID;
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
      ncvars[ncvarid].gmapid          = UNDEFID;
682
683
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
684
      ncvars[ncvarid].truncation      = 0;
685
      ncvars[ncvarid].position        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
686
      ncvars[ncvarid].positive        = 0;
687
688
      ncvars[ncvarid].chunked         = 0;
      ncvars[ncvarid].chunktype       = UNDEFID;
689
690
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
691
      ncvars[ncvarid].missval         = 0;
692
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693
694
695
696
697
698
      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;
699
      ncvars[ncvarid].extra[0]        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
700
701
702
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
703
704
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
705
706
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
707
      ncvars[ncvarid].ensdata         = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
708
709
710
    }
}

711
static
712
void cdf_set_var(ncvar_t *ncvars, int ncvarid, short isvar)
713
714
715
{
  if ( ncvars[ncvarid].isvar != UNDEFID &&
       ncvars[ncvarid].isvar != isvar   &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
716
       ncvars[ncvarid].warn  == false )
717
718
719
720
    {
      if ( ! ncvars[ncvarid].ignore )
        Warning("Inconsistent variable definition for %s!", ncvars[ncvarid].name);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
      ncvars[ncvarid].warn = true;
722
723
724
725
726
727
728
      isvar = FALSE;
    }

  ncvars[ncvarid].isvar = isvar;
}

static
729
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
730
731
732
733
734
735
736
737
738
739
740
741
{
  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
742
int scan_hybrid_formula(int ncid, int ncfvarid, int *apvarid, int *bvarid, int *psvarid, int *avarid, int *p0varid)
743
{
744
  int status = 0;
745
746
747
  *apvarid = -1;
  *bvarid  = -1;
  *psvarid = -1;
748
749
  *avarid  = -1;
  *p0varid = -1;
Thomas Jahns's avatar
Thomas Jahns committed
750
  enum { attstringlen = 8192 }; char attstring[attstringlen];
751
752
753
  cdfGetAttText(ncid, ncfvarid, "formula", attstringlen, attstring);
  if ( strcmp(attstring, "p = ap + b*ps") == 0 )
    {
754
      status = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
      bool lstop = false;
756
757
758
759
760
761
762
763
764
765
      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
766
          if ( *pstring == 0 ) lstop = true;
767
768
769
770
771
772
          *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
773
          if ( *pstring == 0 ) lstop = true;
774
775
          *pstring++ = 0;

Thomas Jahns's avatar
Thomas Jahns committed
776
777
          int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
          if ( status_nc == NC_NOERR )
778
779
780
781
782
            {
              if      ( strcmp(tagname, "ap:") == 0 ) *apvarid = dimvarid;
              else if ( strcmp(tagname, "b:")  == 0 ) *bvarid  = dimvarid;
              else if ( strcmp(tagname, "ps:") == 0 ) *psvarid = dimvarid;
            }
783
          else if ( strcmp(tagname, "ps:") != 0 )
784
            {
Thomas Jahns's avatar
Thomas Jahns committed
785
              Warning("%s - %s", nc_strerror(status_nc), varname);
786
787
788
789
790
            }

          if ( lstop ) break;
        }
    }
791
792
793
  else if ( strcmp(attstring, "xxxp = a*p0 + b*ps") == 0 )
    {
      status = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
794
      bool lstop = false;
795
796
797
798
799
800
801
802
803
804
      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
805
          if ( *pstring == 0 ) lstop = true;
806
807
808
809
810
811
          *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
812
          if ( *pstring == 0 ) lstop = true;
813
814
          *pstring++ = 0;

Thomas Jahns's avatar
Thomas Jahns committed
815
816
          int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
          if ( status_nc == NC_NOERR )
817
818
819
820
821
822
823
824
            {
              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
825
              Warning("%s - %s", nc_strerror(status_nc), varname);
826
827
828
829
830
831
832
            }

          if ( lstop ) break;
        }
    }

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
834
}

835
static
836
bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
837
{
838
  bool status = false;
839
840
841
842
843
  int ncfvarid = ncvarid;
  ncvar_t *ncvar = &ncvars[ncvarid];

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

846
      status = true;
847
848
      ncvar->zaxistype = ZAXIS_HYBRID;
      int dimid = ncvar->dimids[0];
Thomas Jahns's avatar
Thomas Jahns committed
849
      size_t dimlen = ncdims[dimid].len;
850

851
852
      int ret;
      int apvarid1 = -1, bvarid1 = -1, psvarid1 = -1, avarid1 = -1, p0varid1 = -1;
853
      if ( ncvars[ncfvarid].lformula && ncvars[ncfvarid].lformulaterms )
854
        ret = scan_hybrid_formula(ncid, ncfvarid, &apvarid1, &bvarid1, &psvarid1, &avarid1, &p0varid1);
855
856
857
      if ( apvarid1 != -1 ) ncvars[apvarid1].isvar = FALSE;
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar  = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
858
859
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
860
861
862
863

      if ( ncvar->bounds != UNDEFID && ncvars[ncvar->bounds].lformula && ncvars[ncvar->bounds].lformulaterms )
        {
          ncfvarid = ncvar->bounds;
864
865
          int apvarid2 = -1, bvarid2 = -1, psvarid2 = -1, avarid2 = -1, p0varid2 = -1;
          ret = 0;
866
          if ( ncvars[ncfvarid].lformula && ncvars[ncfvarid].lformulaterms )
867
868
869
            ret = scan_hybrid_formula(ncid, ncfvarid, &apvarid2, &bvarid2, &psvarid2, &avarid2, &p0varid2);
          if ( ret == 1 ) avarid2 = apvarid2;
          if ( avarid2 != -1 && bvarid2 != -1 )
870
            {
871
              ncvars[avarid2].isvar = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
872
              ncvars[bvarid2].isvar = FALSE;
873

874
              if ( dimid == ncvars[avarid2].dimids[0] && ncdims[ncvars[avarid2].dimids[1]].len == 2 )
875
                {
876
877
878
879
                  double px = 1;
                  if ( ret == 2 && p0varid1 == p0varid2 )
                    cdf_get_var_double(ncid, p0varid2, &px);

880
                  double abuf[dimlen*2], bbuf[dimlen*2];
881
                  cdf_get_var_double(ncid, avarid2, abuf);
882
883
884
885
886
                  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
887
                  size_t vctsize = (dimlen+1)*2;
888
                  double *vct = (double *) Malloc(vctsize * sizeof(double));
Thomas Jahns's avatar
Thomas Jahns committed
889
                  for ( size_t i = 0; i < dimlen; ++i )
890
891
892
893
894
                    {
                      vct[i] = abuf[i*2];
                      vct[i+dimlen+1] = bbuf[i*2];
                    }
                  vct[dimlen]     = abuf[dimlen*2-1];
895
                  vct[dimlen*2+1] = bbuf[dimlen*2-1];
896

897
898
899
                  if ( ret == 2 && IS_NOT_EQUAL(px, 1) )
                    for ( size_t i = 0; i < dimlen+1; ++i ) vct[i] *= px;

900
901
902
903
904
905
906
907
908
909
                  ncvar->vct = vct;
                  ncvar->vctsize = vctsize;
                }
            }
        }
    }

  return status;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
910
static
911
void cdf_print_vars(const ncvar_t *ncvars, int nvars, const char *oname)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
912
{
913
  char axis[7];
914
  static const char iaxis[] = {'t', 'z', 'y', 'x'};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
915

916
917
  fprintf(stderr, "%s:\n", oname);

918
  for ( int ncvarid = 0; ncvarid < nvars; ncvarid++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
919
    {
920
      int ndim = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
921
      if ( ncvars[ncvarid].isvar )
922
923
        {
          axis[ndim++] = 'v';
924
          axis[ndim++] = ':';
925
          for ( int i = 0; i < ncvars[ncvarid].ndims; i++ )
926
927
928
929
930
931
932
            {/*
              if      ( ncvars[ncvarid].tvarid != -1 ) axis[ndim++] = iaxis[0];
              else if ( ncvars[ncvarid].zvarid != -1 ) axis[ndim++] = iaxis[1];
              else if ( ncvars[ncvarid].yvarid != -1 ) axis[ndim++] = iaxis[2];
              else if ( ncvars[ncvarid].xvarid != -1 ) axis[ndim++] = iaxis[3];
              else
             */
933
934
935
936
937
              if      ( ncvars[ncvarid].dimtype[i] == T_AXIS ) axis[ndim++] = iaxis[0];
              else if ( ncvars[ncvarid].dimtype[i] == Z_AXIS ) axis[ndim++] = iaxis[1];
              else if ( ncvars[ncvarid].dimtype[i] == Y_AXIS ) axis[ndim++] = iaxis[2];
              else if ( ncvars[ncvarid].dimtype[i] == X_AXIS ) axis[ndim++] = iaxis[3];
              else                                             axis[ndim++] = '?';
938
939
            }
        }