stream_cdf_i.c 140 KB
Newer Older
1
#ifdef HAVE_CONFIG_H
2
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3
4
#endif

5
6
#ifdef HAVE_LIBNETCDF

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include <ctype.h>
8
#include <limits.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
9

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


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

27
28
29
#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2

Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
typedef struct {
31
  int     dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
34
  int     ncvarid;
  int     dimtype;
  size_t  len;
35
  char    name[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
}
37
ncdim_t;
38

39
40
41
#define  MAX_COORDVARS  4
#define  MAX_AUXVARS    4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
typedef struct {
43
  int      ncid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
45
  int      isvar;
  bool     ignore;
46
47
  bool     isx;
  bool     isy;
48
  bool     isc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
  bool     islon;
  bool     islat;
  bool     islev;
  bool     istime;
  bool     warn;
54
  bool     calendar;
55
  bool     climatology;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
  bool     lformulaterms;
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
  int      bounds;
  int      gridID;
  int      zaxisID;
63
  int      timetype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
68
69
70
71
  int      gridtype;
  int      zaxistype;
  int      xdim;
  int      ydim;
  int      zdim;
  int      xvarid;
  int      yvarid;
  int      zvarid;
72
  int      cvarids[MAX_COORDVARS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
  int      tvarid;
74
  int      psvarid;
75
  int      p0varid;
76
  int      ncoordvars;
77
78
79
  int      coordvarids[MAX_COORDVARS];
  int      nauxvars;
  int      auxvarids[MAX_AUXVARS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
  int      cellarea;
  int      tableID;
  int      truncation;
83
  int      position;
84
85
  bool     defmissval;
  bool     deffillval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
88
  int      xtype;
  int      gmapid;
  int      positive;
89
  int      ndims;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
91
  int      dimids[8];
  int      dimtype[8];
92
  size_t   chunks[8];
93
94
  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
113
114
  int      typeOfEnsembleForecast;
  int      numberOfForecastsInEnsemble;
  int      perturbationNumber;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
}
116
117
ncvar_t;

118
119

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

  *rdate = 0;
  *rtime = 0;

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
  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++;
            }
        }
    }
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
174
175

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

184
  return timeunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
}

187
188
189
static
void setForecastTime(const char *timestr, taxis_t *taxis)
{
Thomas Jahns's avatar
Thomas Jahns committed
190
191
192
193
194
  size_t len = strlen(timestr);
  if ( len != 0 )
    scanTimeString(timestr, &taxis->fdate, &taxis->ftime);
  else
    taxis->fdate = taxis->ftime = 0;
195
196
197
198
199
}

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

204
  size_t len = strlen(timeunits);
205
206
  while ( isspace(*timeunits) && len ) { timeunits++; len--; }

Thomas Jahns's avatar
Thomas Jahns committed
207
  char *restrict tu = (char *)Malloc((len+1) * sizeof(char));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208

Thomas Jahns's avatar
Thomas Jahns committed
209
210
  for ( size_t i = 0; i < len; i++ ) tu[i] = (char)tolower((int)timeunits[i]);
  tu[len] = 0;
211

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

Thomas Jahns's avatar
Thomas Jahns committed
219
  size_t pos = 0;
220
  while ( pos < len && !isspace(tu[pos]) ) ++pos;
Thomas Jahns's avatar
Thomas Jahns committed
221
  if ( tu[pos] )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
    {
Thomas Jahns's avatar
Thomas Jahns committed
223
      while ( isspace(tu[pos]) ) ++pos;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224

Thomas Jahns's avatar
Thomas Jahns committed
225
      if ( str_is_equal(tu+pos, "since") )
226
        taxistype = TAXIS_RELATIVE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227

228
      while ( pos < len && !isspace(tu[pos]) ) ++pos;
Thomas Jahns's avatar
Thomas Jahns committed
229
      if ( tu[pos] )
230
        {
Thomas Jahns's avatar
Thomas Jahns committed
231
          while ( isspace(tu[pos]) ) ++pos;
232

233
          if ( taxistype == TAXIS_ABSOLUTE )
234
            {
235
              if ( timeunit == TUNIT_DAY )
236
                {
237
238
239
240
241
                  if ( !str_is_equal(tu+pos, "%y%m%d.%f") )
                    {
                      Message("Unsupported format %s for TIMEUNIT day!", tu+pos);
                      timeunit = -1;
                    }
242
                }
243
              else if ( timeunit == TUNIT_MONTH )
244
                {
245
246
247
248
249
                  if ( !str_is_equal(tu+pos, "%y%m.%f") )
                    {
                      Message("Unsupported format %s for TIMEUNIT month!", tu+pos);
                      timeunit = -1;
                    }
250
251
                }
            }
252
          else if ( taxistype == TAXIS_RELATIVE )
253
            {
Thomas Jahns's avatar
Thomas Jahns committed
254
              scanTimeString(tu+pos, &rdate, &rtime);
255

256
257
              taxis->rdate = rdate;
              taxis->rtime = rtime;
258

259
              if ( CDI_Debug )
260
                Message("rdate = %lld  rtime = %d", rdate, rtime);
261
            }
262
        }
263
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264

265
  taxis->type = taxistype;
266
  taxis->unit = timeunit;
Deike Kleberg's avatar
Deike Kleberg committed
267

268
  Free(tu);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269

270
  if ( CDI_Debug )
271
    Message("taxistype = %d  unit = %d", taxistype, timeunit);
Deike Kleberg's avatar
Deike Kleberg committed
272

273
274
  return 0;
}
275

276
static
277
bool xtypeIsText(int xtype)
278
{
279
  bool isText = ( xtype == NC_CHAR )
280
#ifdef HAVE_NETCDF4
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    || ( xtype == NC_STRING )
#endif
    ;
  return isText;
}

static
bool xtypeIsFloat(nc_type xtype)
{
  bool isFloat = xtype == NC_FLOAT || xtype == NC_DOUBLE;

  return isFloat;
}

static
bool xtypeIsInt(nc_type xtype)
{
  bool isInt = xtype == NC_SHORT || xtype == NC_INT
            || xtype == NC_BYTE
300
#ifdef HAVE_NETCDF4
301
302
303
304
305
306
307
308
309
310
311
312
313
            || xtype == NC_USHORT || xtype == NC_UINT
            || xtype == NC_UBYTE
#endif
             ;

  return isInt;
}

static
int cdfInqDatatype(int xtype, bool lunsigned)
{
  int datatype = -1;

314
#ifdef HAVE_NETCDF4
315
316
317
318
  if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
#endif

  if      ( xtype == NC_BYTE   )  datatype = CDI_DATATYPE_INT8;
319
  else if ( xtype == NC_CHAR   )  datatype = CDI_DATATYPE_UINT8;
320
321
322
323
  else if ( xtype == NC_SHORT  )  datatype = CDI_DATATYPE_INT16;
  else if ( xtype == NC_INT    )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_FLOAT  )  datatype = CDI_DATATYPE_FLT32;
  else if ( xtype == NC_DOUBLE )  datatype = CDI_DATATYPE_FLT64;
324
#ifdef HAVE_NETCDF4
325
326
327
328
329
330
331
  else if ( xtype == NC_UBYTE  )  datatype = CDI_DATATYPE_UINT8;
  else if ( xtype == NC_LONG   )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_USHORT )  datatype = CDI_DATATYPE_UINT16;
  else if ( xtype == NC_UINT   )  datatype = CDI_DATATYPE_UINT32;
  else if ( xtype == NC_INT64  )  datatype = CDI_DATATYPE_FLT64;
  else if ( xtype == NC_UINT64 )  datatype = CDI_DATATYPE_FLT64;
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
332

333
334
335
336
337
338
  return datatype;
}

static
void cdfGetAttInt(int fileID, int ncvarid, const char *attname, size_t attlen, int *attint)
{
339
  *attint = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340

341
342
  nc_type atttype;
  size_t nc_attlen;
343
344
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Deike Kleberg's avatar
Deike Kleberg committed
345

346
  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
347
    {
348
349
      bool lalloc = nc_attlen > attlen;
      int *pintatt = lalloc ? (int *)(Malloc(nc_attlen*sizeof(int))) : attint;
350
      cdf_get_att_int(fileID, ncvarid, attname, pintatt);
351
      if ( lalloc )
352
        {
353
          memcpy(attint, pintatt, attlen*sizeof(int));
354
          Free(pintatt);
355
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
357
358
    }
}

359
static
360
void cdfGetAttDouble(int fileID, int ncvarid, char *attname, size_t attlen, double *attdouble)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
{
362
  *attdouble = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363

364
365
  nc_type atttype;
  size_t nc_attlen;
366
367
  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368

369
  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
    {
371
372
      bool lalloc = nc_attlen > attlen;
      double *pdoubleatt = lalloc ? (double*)Malloc(nc_attlen*sizeof(double)) : attdouble;
373
      cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);
374
      if ( lalloc )
375
        {
376
          memcpy(attdouble, pdoubleatt, attlen*sizeof(double));
377
378
379
380
          Free(pdoubleatt);
        }
    }
}
381

382
static
383
384
385
386
387
388
389
390
391
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
392
#ifdef HAVE_NETCDF4
393
394
395
396
397
398
399
400
401
402
403
           || atttype == NC_STRING
#endif
           ) )
    {
      status = true;
    }

  return status;
}

static
404
void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
405
{
406
407
408
409
410
  nc_type atttype;
  size_t nc_attlen;

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

412
  if ( atttype == NC_CHAR )
413
    {
414
415
416
417
      char attbuf[65636];
      if ( nc_attlen < sizeof(attbuf) )
        {
          cdf_get_att_text(fileID, ncvarid, attname, attbuf);
418

419
          if ( nc_attlen > (attlen-1) ) nc_attlen = (attlen-1);
420

421
422
423
424
          attbuf[nc_attlen++] = 0;
          memcpy(atttext, attbuf, nc_attlen);
        }
      else
425
        {
426
          atttext[0] = 0;
427
        }
428
    }
429
#ifdef HAVE_NETCDF4
430
  else if ( atttype == NC_STRING )
431
    {
432
      if ( nc_attlen == 1 )
433
        {
434
435
          char *attbuf = NULL;
          cdf_get_att_string(fileID, ncvarid, attname, &attbuf);
436

437
          size_t ssize = strlen(attbuf) + 1;
438

439
          if ( ssize > attlen ) ssize = attlen;
440
441
442
443
444
445
446
          memcpy(atttext, attbuf, ssize);
          atttext[ssize - 1] = 0;
          Free(attbuf);
        }
      else
        {
          atttext[0] = 0;
447
        }
448
    }
449
450
#endif
}
451

452

453
void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
{
455
456
  bool laddoffset   = IS_NOT_EQUAL(addoffset, 0);
  bool lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
457

458
  if ( laddoffset && lscalefactor )
459
    {
460
461
462
463
464
465
466
467
468
469
470
471
      for (size_t i = 0; i < size; ++i )
        data[i] = data[i] * scalefactor + addoffset;
    }
  else if (lscalefactor)
    {
      for (size_t i = 0; i < size; ++i )
        data[i] *= scalefactor;
    }
  else if (laddoffset)
    {
      for (size_t i = 0; i < size; ++i )
        data[i] += addoffset;
472
473
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
474

475
476
static
void cdfCreateRecords(stream_t *streamptr, int tsID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
  if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479

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

482
483
  int vlistID  = streamptr->vlistID;

484
485
486
487
488
489
490
491
  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
492
493
  if ( tsID == 0 )
    {
494
      int nvrecs = nrecs; /* use all records at first timestep */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495

Uwe Schulzweida's avatar
Uwe Schulzweida committed
496
      streamptr->nrecs += nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497

498
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
499
500
501
      destTstep->nrecs      = nrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
502
      destTstep->curRecID   = CDI_UNDEFID;
503
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs*sizeof (int));;
504
      for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505

506
      record_t *records = destTstep->records;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507

508
      for ( int varID = 0, recID = 0; varID < nvars; varID++ )
509
        {
510
511
512
          int zaxisID = vlistInqVarZaxis(vlistID, varID);
          int nlev    = zaxisInqSize(zaxisID);
          for ( int levelID = 0; levelID < nlev; levelID++ )
513
514
            {
              recordInitEntry(&records[recID]);
515
516
              records[recID].varID   = (short)varID;
              records[recID].levelID = (short)levelID;
517
518
519
              recID++;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
521
522
    }
  else if ( tsID == 1 )
    {
523
524
      int nvrecs = 0;
      for ( int varID = 0; varID < nvars; varID++ )
525
        {
526
          if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
527
            {
528
              int zaxisID = vlistInqVarZaxis(vlistID, varID);
529
530
531
              nvrecs += zaxisInqSize(zaxisID);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532

Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
      streamptr->nrecs += nvrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534

535
      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
536
537
538
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
539
      destTstep->curRecID   = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540

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

543
      if ( nvrecs )
544
        {
545
          destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
546
          for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
547
            {
548
              int varID = destTstep->records[recID].varID;
549
              if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
550
                {
551
                  destTstep->recIDs[vrecID++] = recID;
552
553
554
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555
556
557
    }
  else
    {
558
      if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);
559

560
      int nvrecs = streamptr->tsteps[1].nrecs;
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
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
568
      destTstep->curRecID   = CDI_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
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs * sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
573

574
      memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
576
577
    }
}

578
static
579
int cdf_time_dimid(int fileID, int ndims, int nvars, ncdim_t *ncdims)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
{
581
582
  char dimname[CDI_MAX_NAME];

583
  for ( int dimid = 0; dimid < ndims; ++dimid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
584
    {
585
      dimname[0] = 0;
586
      cdf_inq_dimname(fileID, ncdims[dimid].dimid, dimname);
587
      if ( str_is_equal(dimname, "time") || str_is_equal(dimname, "Time") ) return dimid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
589
    }

590
  for ( int varid = 0; varid < nvars; ++varid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
    {
592
      nc_type xtype;
593
      int nvdims, nvatts, dimids[9];
594
      cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
595
596
597
598
599
600
601
602
      for ( int i = 0; i < nvdims; ++i )
        for ( int dimid = 0; dimid < ndims; ++dimid )
          if ( ncdims[dimid].dimid == dimids[i] )
            {
              dimids[i] = dimid;
              break;
            }

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
                  strToLower(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 CDI_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
      ncdims[ncdimid].dimid        = CDI_UNDEFID;
630
631
      ncdims[ncdimid].ncvarid      = CDI_UNDEFID;
      ncdims[ncdimid].dimtype      = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
632
633
634
635
636
637
      ncdims[ncdimid].len          = 0;
      ncdims[ncdimid].name[0]      = 0;
    }
}

static
638
void init_ncvars(long nvars, ncvar_t *ncvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
639
{
640
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
Deike Kleberg's avatar
Deike Kleberg committed
641
    {
642
643
      ncvars[ncvarid].ncid            = CDI_UNDEFID;
      ncvars[ncvarid].isvar           = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
      ncvars[ncvarid].ignore          = false;
645
646
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
647
      ncvars[ncvarid].isc             = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
648
649
650
651
652
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
653
      ncvars[ncvarid].calendar        = false;
654
      ncvars[ncvarid].climatology     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
655
      ncvars[ncvarid].lformulaterms   = false;
656
      ncvars[ncvarid].timetype        = TIME_CONSTANT;
657
658
      ncvars[ncvarid].param           = CDI_UNDEFID;
      ncvars[ncvarid].code            = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
659
      ncvars[ncvarid].tabnum          = 0;
660
661
662
663
664
665
666
667
668
669
670
671
672
673
      ncvars[ncvarid].bounds          = CDI_UNDEFID;
      ncvars[ncvarid].gridID          = CDI_UNDEFID;
      ncvars[ncvarid].zaxisID         = CDI_UNDEFID;
      ncvars[ncvarid].gridtype        = CDI_UNDEFID;
      ncvars[ncvarid].zaxistype       = CDI_UNDEFID;
      ncvars[ncvarid].xdim            = CDI_UNDEFID;
      ncvars[ncvarid].ydim            = CDI_UNDEFID;
      ncvars[ncvarid].zdim            = CDI_UNDEFID;
      ncvars[ncvarid].xvarid          = CDI_UNDEFID;
      ncvars[ncvarid].yvarid          = CDI_UNDEFID;
      ncvars[ncvarid].zvarid          = CDI_UNDEFID;
      ncvars[ncvarid].tvarid          = CDI_UNDEFID;
      ncvars[ncvarid].psvarid         = CDI_UNDEFID;
      ncvars[ncvarid].p0varid         = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
      ncvars[ncvarid].ncoordvars      = 0;
675
      for ( int i = 0; i < MAX_COORDVARS; ++i )
676
677
678
679
        {
          ncvars[ncvarid].coordvarids[i]  = CDI_UNDEFID;
          ncvars[ncvarid].cvarids[i]      = CDI_UNDEFID;
        }
680
681
      ncvars[ncvarid].nauxvars      = 0;
      for ( int i = 0; i < MAX_AUXVARS; ++i )
682
683
684
        ncvars[ncvarid].auxvarids[i]  = CDI_UNDEFID;
      ncvars[ncvarid].cellarea        = CDI_UNDEFID;
      ncvars[ncvarid].tableID         = CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
685
686
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
687
      ncvars[ncvarid].gmapid          = CDI_UNDEFID;
688
689
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
690
      ncvars[ncvarid].truncation      = 0;
691
      ncvars[ncvarid].position        = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
692
      ncvars[ncvarid].positive        = 0;
693
      ncvars[ncvarid].chunked         = 0;
694
      ncvars[ncvarid].chunktype       = CDI_UNDEFID;
695
696
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
697
      ncvars[ncvarid].missval         = 0;
698
      ncvars[ncvarid].fillval         = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699
700
701
702
703
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
704
705
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
Deike Kleberg's avatar
Deike Kleberg committed
706
707
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
708
709
710
      ncvars[ncvarid].typeOfEnsembleForecast      = -1;
      ncvars[ncvarid].numberOfForecastsInEnsemble = -1;
      ncvars[ncvarid].perturbationNumber          = -1;
711
712
713
714
715
      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
716
717
718
    }
}

719
static
720
void cdf_set_var(ncvar_t *ncvars, int ncvarid, short isvar)
721
{
722
  if ( ncvars[ncvarid].isvar != CDI_UNDEFID &&
723
       ncvars[ncvarid].isvar != isvar   &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
724
       ncvars[ncvarid].warn  == false )
725
726
727
728
    {
      if ( ! ncvars[ncvarid].ignore )
        Warning("Inconsistent variable definition for %s!", ncvars[ncvarid].name);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
729
      ncvars[ncvarid].warn = true;
730
731
732
733
734
735
736
      isvar = FALSE;
    }

  ncvars[ncvarid].isvar = isvar;
}

static
737
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
738
{
739
  if ( ncvars[ncvarid].dimtype[dimid] != CDI_UNDEFID &&
740
741
742
743
744
745
746
747
748
749
       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
750
void scan_hybrid_formulaterms(int ncid, int ncfvarid, int *avarid, int *bvarid, int *psvarid, int *p0varid)
751
{
752
  *avarid  = -1;
753
754
  *bvarid  = -1;
  *psvarid = -1;
755
  *p0varid = -1;
756

757
  char attstring[1024];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
758
759
  cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
  char *pstring = attstring;
760

761
  bool lstop = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
  for ( int i = 0; i < 4; i++ )
763
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
765
766
767
768
769
770
771
772
773
774
775
776
777
      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;

778
      int dimvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
779
780
      int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
      if ( status_nc == NC_NOERR )
781
        {
782
783
          if      ( strcmp(tagname, "ap:") == 0 ) *avarid  = dimvarid;
          else if ( strcmp(tagname, "a:")  == 0 ) *avarid  = dimvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
784
785
786
          else if ( strcmp(tagname, "b:")  == 0 ) *bvarid  = dimvarid;
          else if ( strcmp(tagname, "ps:") == 0 ) *psvarid = dimvarid;
          else if ( strcmp(tagname, "p0:") == 0 ) *p0varid = dimvarid;
787
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
788
789
790
791
792
793
      else if ( strcmp(tagname, "ps:") != 0 )
        {
          Warning("%s - %s", nc_strerror(status_nc), varname);
        }

      if ( lstop ) break;
794
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
795
796
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
static
void readVCT(int ncid, int ndims2, size_t dimlen, size_t dimlen2, int avarid2, int bvarid2, double *vct)
{
  double *abuf = (double*) Malloc(dimlen*2*sizeof(double));
  double *bbuf = (double*) Malloc(dimlen*2*sizeof(double));
  cdf_get_var_double(ncid, avarid2, abuf);
  cdf_get_var_double(ncid, bvarid2, bbuf);

  if ( ndims2 == 2 )
    {
      for ( size_t i = 0; i < dimlen; ++i )
        {
          vct[i] = abuf[i*2];
          vct[i+dimlen+1] = bbuf[i*2];
        }
      vct[dimlen]     = abuf[dimlen*2-1];
      vct[dimlen*2+1] = bbuf[dimlen*2-1];
    }
  else
    {
      for ( size_t i = 0; i < dimlen2; ++i )
        {
          vct[i] = abuf[i];
          vct[i+dimlen+1] = bbuf[i];
        }
    }

  Free(abuf);
  Free(bbuf);
}

828
static
829
bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
830
{
831
  bool status = false;
832
833
834
835
  ncvar_t *ncvar = &ncvars[ncvarid];

  if ( strcmp(ncvar->stdname, "atmosphere_hybrid_sigma_pressure_coordinate") == 0 )
    {
836
      CDI_convention = CDI_CONVENTION_CF;
837

838
      status = true;
839
      ncvar->zaxistype = ZAXIS_HYBRID;
840
      //int ndims = ncvar->ndims;
841
      int dimid = ncvar->dimids[0];
Thomas Jahns's avatar
Thomas Jahns committed
842
      size_t dimlen = ncdims[dimid].len;
843
      int avarid1 = -1, bvarid1 = -1, psvarid1 = -1, p0varid1 = -1;
844
      int ncfvarid = ncvarid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
845
      if ( ncvars[ncfvarid].lformulaterms )
846
847
        scan_hybrid_formulaterms(ncid, ncfvarid, &avarid1, &bvarid1, &psvarid1, &p0varid1);
      // printf("avarid1, bvarid1, psvarid1, p0varid1 %d %d %d %d\n", avarid1, bvarid1, psvarid1, p0varid1);
848
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
849
850
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
851
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
852

853
      if ( ncvar->bounds != CDI_UNDEFID && ncvars[ncvar->bounds].lformulaterms )
854
855
        {
          ncfvarid = ncvar->bounds;
856
          int avarid2 = -1, bvarid2 = -1, psvarid2 = -1, p0varid2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857
          if ( ncvars[ncfvarid].lformulaterms )
858
859
            scan_hybrid_formulaterms(ncid, ncfvarid, &avarid2, &bvarid2, &psvarid2, &p0varid2);
          // printf("avarid2, bvarid2, psvarid2, p0varid2 %d %d %d %d\n", avarid2, bvarid2, psvarid2, p0varid2);
860
          if ( avarid2 != -1 && bvarid2 != -1 )
861
            {
862
              ncvars[avarid2].isvar = FALSE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
863
              ncvars[bvarid2].isvar = FALSE;
864

865
866
867
868
869
870
              int ndims2 = ncvars[avarid2].ndims;
              int dimid2 = ncvars[avarid2].dimids[0];
              size_t dimlen2 = ncdims[dimid2].len;

              if ( (ndims2 == 2 && dimid == ncvars[avarid2].dimids[0] ) ||
                   (ndims2 == 1 && dimlen == dimlen2-1 ) )
871
                {
872
                  double px = 1;
873
                  if ( p0varid1 != -1 && p0varid1 == p0varid2 )
874
875
                    cdf_get_var_double(ncid, p0varid2, &px);

Thomas Jahns's avatar
Thomas Jahns committed
876
                  size_t vctsize = (dimlen+1)*2;
877
                  double *vct = (double *) Malloc(vctsize*sizeof(double));
878

Uwe Schulzweida's avatar
Uwe Schulzweida committed
879
                  readVCT(ncid, ndims2, dimlen, dimlen2, avarid2, bvarid2