stream_cdf_o.c 64.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#ifdef HAVE_LIBNETCDF

#include "dmemory.h"
#include "cdi_int.h"
#include "cdi_uuid.h"
#include "stream_cdf.h"
#include "cdf_int.h"
#include "varscan.h"
#include "vlist.h"
#include "zaxis.h"


#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2


Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
static const char bndsName[] = "bnds";


24
25
26
27
28
29
30
31
32
33
void cdfCopyRecord(stream_t *streamptr2, stream_t *streamptr1)
{
  int vlistID1 = streamptr1->vlistID;
  int tsID     = streamptr1->curTsID;
  int vrecID   = streamptr1->tsteps[tsID].curRecID;
  int recID    = streamptr1->tsteps[tsID].recIDs[vrecID];
  int ivarID   = streamptr1->tsteps[tsID].records[recID].varID;
  int gridID   = vlistInqVarGrid(vlistID1, ivarID);
  int datasize = gridInqSize(gridID);
  int datatype = vlistInqVarDatatype(vlistID1, ivarID);
34
  int memtype  = datatype != CDI_DATATYPE_FLT32 ? MEMTYPE_DOUBLE : MEMTYPE_FLOAT;
35

Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
  void *data = Malloc((size_t)datasize
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
             * (memtype == MEMTYPE_DOUBLE ? sizeof(double) : sizeof(float)));

  int nmiss;
  cdf_read_record(streamptr1, memtype, data, &nmiss);
  cdf_write_record(streamptr2, memtype, data, nmiss);

  Free(data);
}


void cdfDefRecord(stream_t *streamptr)
{
  (void)streamptr;
}

static
void cdfDefTimeValue(stream_t *streamptr, int tsID)
{
  int fileID = streamptr->fileID;

  if ( CDI_Debug )
    Message("streamID = %d, fileID = %d", streamptr->self, fileID);

  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;

  if ( streamptr->ncmode == 1 )
    {
      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

  double timevalue = cdiEncodeTimeval(taxis->vdate, taxis->vtime, &streamptr->tsteps[0].taxis);
  if ( CDI_Debug ) Message("tsID = %d  timevalue = %f", tsID, timevalue);

  int ncvarid = streamptr->basetime.ncvarid;
72
  size_t index = (size_t)tsID;
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
  cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);

  if ( taxis->has_bounds )
    {
      size_t start[2], count[2];

      ncvarid = streamptr->basetime.ncvarboundsid;

      timevalue = cdiEncodeTimeval(taxis->vdate_lb, taxis->vtime_lb, &streamptr->tsteps[0].taxis);
      start[0] = (size_t)tsID; count[0] = 1; start[1] = 0; count[1] = 1;
      cdf_put_vara_double(fileID, ncvarid, start, count, &timevalue);

      timevalue = cdiEncodeTimeval(taxis->vdate_ub, taxis->vtime_ub, &streamptr->tsteps[0].taxis);
      start[0] = (size_t)tsID; count[0] = 1; start[1] = 1; count[1] = 1;
      cdf_put_vara_double(fileID, ncvarid, start, count, &timevalue);
    }

  ncvarid = streamptr->basetime.leadtimeid;
91
  if ( taxis->type == TAXIS_FORECAST && ncvarid != CDI_UNDEFID )
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    {
      timevalue = taxis->fc_period;
      cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);
    }
}

void cdfDefTimestep(stream_t *streamptr, int tsID)
{
  int vlistID = streamptr->vlistID;

  if ( vlistHasTime(vlistID) ) cdfDefTime(streamptr);

  cdfDefTimeValue(streamptr, tsID);
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
void cdfDefComplex(stream_t *streamptr, int gridID, int gridindex)
109
{
110
  int dimID;
111
  ncgrid_t *ncgrid = streamptr->ncgrid;
112

Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
  for ( int index = 0; index < gridindex; ++index )
114
    {
115
      if ( ncgrid[index].ncIDs[CDF_DIMID_X] != CDI_UNDEFID )
116
        {
117
          int gridID0 = ncgrid[index].gridID;
118
119
120
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_SPECTRAL || gridtype0 == GRID_FOURIER )
            {
121
              dimID = ncgrid[index].ncIDs[CDF_DIMID_X];
122
              goto dimIDEstablished;
123
124
125
126
            }
        }
    }

127
128
129
130
131
132
133
134
135
136
137
  {
    static const char axisname[] = "nc2";
    size_t dimlen = 2;
    int fileID  = streamptr->fileID;

    if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
    cdf_def_dim(fileID, axisname, dimlen, &dimID);
    cdf_enddef(fileID);
    streamptr->ncmode = 2;
  }
  dimIDEstablished:
138
  ncgrid[gridindex].gridID = gridID;
139
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = dimID;
140
141
}

142
struct idSearch
143
{
144
145
146
  int numNonMatching, foundID;
  size_t foundIdx;
};
147

148
149
150
151
152
153
154
155
156
static inline struct idSearch
cdfSearchIDBySize(size_t startIdx, size_t numIDs, const ncgrid_t ncgrid[numIDs],
                  int ncIDType, int searchType, int searchSize,
                  int (*typeInq)(int id), int (*sizeInq)(int id))
{
  int numNonMatching = 0,
    foundID = CDI_UNDEFID;
  size_t foundIdx = SIZE_MAX;
  for ( size_t index = startIdx; index < numIDs; index++ )
157
    {
158
      if ( ncgrid[index].ncIDs[ncIDType] != CDI_UNDEFID )
159
        {
160
161
162
          int id0 = ncgrid[index].gridID,
            id0Type = typeInq(id0);
          if ( id0Type == searchType )
163
            {
164
165
              int size0 = sizeInq(id0);
              if ( searchSize == size0 )
166
                {
167
168
                  foundID = ncgrid[index].ncIDs[ncIDType];
                  foundIdx = index;
169
170
                  break;
                }
171
              numNonMatching++;
172
173
174
            }
        }
    }
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  return (struct idSearch){ .numNonMatching = numNonMatching,
      .foundID = foundID, .foundIdx = foundIdx };
}

static int
cdfGridInqHalfSize(int gridID)
{
  return gridInqSize(gridID)/2;
}


static void
cdfDefSPorFC(stream_t *streamptr, int gridID, int gridindex,
             char *restrict axisname, int gridRefType)
{
  ncgrid_t *ncgrid = streamptr->ncgrid;

  size_t dimlen = (size_t)(gridInqSize(gridID))/2;

  int iz;
  int dimID;
  {
    struct idSearch search
      = cdfSearchIDBySize(0, (size_t)gridindex, ncgrid, CDF_DIMID_Y,
                          gridRefType, (int)dimlen,
                          gridInqType, cdfGridInqHalfSize);
    dimID = search.foundID;
    iz = search.numNonMatching;
  }
204

205
  if ( dimID == CDI_UNDEFID )
206
207
208
209
210
211
212
213
214
215
216
217
218
    {
      int fileID  = streamptr->fileID;
      if ( iz == 0 ) axisname[3] = '\0';
      else           sprintf(&axisname[3], "%1d", iz+1);

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      cdf_def_dim(fileID, axisname, dimlen, &dimID);

      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

219
  ncgrid[gridindex].gridID = gridID;
220
  ncgrid[gridindex].ncIDs[CDF_DIMID_Y] = dimID;
221
222
223
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
void cdfDefSP(stream_t *streamptr, int gridID, int gridindex)
225
226
227
228
229
{
  /*
  char longname[] = "Spherical harmonic coefficient";
  */
  char axisname[5] = "nspX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_SPECTRAL);
231
232
233
234
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
void cdfDefFC(stream_t *streamptr, int gridID, int gridindex)
236
237
{
  char axisname[5] = "nfcX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_FOURIER);
239
240
241
242
}

static const struct cdfDefGridAxisInqs {
  int (*axisSize)(int gridID);
243
  int (*axisDimname)(int cdiID, int key, int size, char *mesg);
244
245
  int (*axisName)(int cdiID, int key, int size, char *mesg);
  int (*axisLongname)(int cdiID, int key, int size, char *mesg);
246
247
  int (*axisUnits)(int cdiID, int key, int size, char *mesg);
  void (*axisStdname)(int cdiID, char *dimstdname);
248
249
250
251
252
  double (*axisVal)(int gridID, int index);
  const double *(*axisValsPtr)(int gridID);
  const double *(*axisBoundsPtr)(int gridID);
} gridInqsX = {
  .axisSize = gridInqXsize,
253
  .axisDimname = cdiGridInqKeyStr,
254
  .axisName = cdiGridInqKeyStr,
255
  .axisLongname = cdiGridInqKeyStr,
256
257
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqXstdname,
258
259
260
261
262
  .axisVal = gridInqXval,
  .axisValsPtr = gridInqXvalsPtr,
  .axisBoundsPtr = gridInqXboundsPtr,
}, gridInqsY = {
  .axisSize = gridInqYsize,
263
  .axisDimname = cdiGridInqKeyStr,
264
  .axisName = cdiGridInqKeyStr,
265
  .axisLongname = cdiGridInqKeyStr,
266
267
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqYstdname,
268
269
270
271
  .axisVal = gridInqYval,
  .axisValsPtr = gridInqYvalsPtr,
  .axisBoundsPtr = gridInqYboundsPtr,
}, gridInqsZ = {
272
  .axisLongname = cdiZaxisInqKeyStr,
273
274
  .axisUnits = cdiZaxisInqKeyStr,
  .axisStdname = zaxisInqStdname,
275
276
};

277
278
static
void cdfPutGridStdAtts(int fileID, int ncvarid, int gridID, int dimtype, const struct cdfDefGridAxisInqs *inqs)
279
280
{
  size_t len;
281
282
283
284
285
286
287
288
289

  char stdname[CDI_MAX_NAME];
  inqs->axisStdname(gridID, stdname);
  if ( (len = strlen(stdname)) )
    cdf_put_att_text(fileID, ncvarid, "standard_name", len, stdname);

  char longname[CDI_MAX_NAME]; longname[0] = 0;
  int keyname = (dimtype == 'Z') ? CDI_KEY_LONGNAME : (dimtype == 'X') ? CDI_KEY_XLONGNAME : CDI_KEY_YLONGNAME;
  inqs->axisLongname(gridID, keyname, CDI_MAX_NAME, longname);
290
  if ( longname[0] && (len = strlen(longname)) )
291
292
    cdf_put_att_text(fileID, ncvarid, "long_name", len, longname);

293
294
295
296
  char units[CDI_MAX_NAME]; units[0] = 0;
  keyname = (dimtype == 'Z') ? CDI_KEY_UNITS : (dimtype == 'X') ? CDI_KEY_XUNITS : CDI_KEY_YUNITS;
  inqs->axisUnits(gridID, keyname, CDI_MAX_NAME, units);
  if ( units[0] && (len = strlen(units)) )
297
    cdf_put_att_text(fileID, ncvarid, "units", len, units);
298
299
300
}

static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
cdfDefTrajLatLon(stream_t *streamptr, int gridID, int gridindex,
302
                 const struct cdfDefGridAxisInqs *inqs, int dimtype)
303
{
304
  nc_type xtype = (gridInqPrec(gridID) == CDI_DATATYPE_FLT32) ? NC_FLOAT : NC_DOUBLE;
305
  ncgrid_t *ncgrid = streamptr->ncgrid;
306
307
308

  int dimlen = inqs->axisSize(gridID);
  if ( dimlen != 1 )
309
    Error("%c size isn't 1 for %s grid!", dimtype, gridNamePtr(gridInqType(gridID)));
310

Thomas Jahns's avatar
Thomas Jahns committed
311
312
  int ncvarid
    = ncgrid[gridindex].ncIDs[dimtype == 'X' ? CDF_DIMID_X : CDF_DIMID_Y];
313

314
  if ( ncvarid == CDI_UNDEFID )
315
316
317
318
319
    {
      int dimNcID = streamptr->basetime.ncvarid;
      int fileID  = streamptr->fileID;
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

320
321
322
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
      int keyname = (dimtype == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
      inqs->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
323
      cdf_def_var(fileID, axisname, xtype, 1, &dimNcID, &ncvarid);
324
      cdfPutGridStdAtts(fileID, ncvarid, gridID, dimtype, inqs);
325
326
327
328
      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

329
  ncgrid[gridindex].gridID = gridID;
Thomas Jahns's avatar
Thomas Jahns committed
330
331
  /* var ID for trajectory !!! */
  ncgrid[gridindex].ncIDs[dimtype == 'X' ? CDF_DIMID_X : CDF_DIMID_Y] = ncvarid;
332
333
334
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
void cdfDefTrajLon(stream_t *streamptr, int gridID, int gridindex)
336
{
337
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsX, 'X');
338
339
340
341
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
void cdfDefTrajLat(stream_t *streamptr, int gridID, int gridindex)
343
{
344
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsY, 'Y');
345
346
347
348
349
350
351
}

static
int checkDimName(int fileID, size_t dimlen, char *dimname)
{
  /* check whether the dimenion name is already defined with the same length */
  unsigned iz = 0;
352
  int dimid = CDI_UNDEFID;
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
  char name[CDI_MAX_NAME];

  size_t len = strlen(dimname);
  memcpy(name, dimname, len + 1);

  do
    {
      if ( iz ) sprintf(name + len, "_%u", iz+1);

      int dimid0, status = nc_inq_dimid(fileID, name, &dimid0);
      if ( status != NC_NOERR )
        break;
      size_t dimlen0;
      cdf_inq_dimlen(fileID, dimid0, &dimlen0);
      if ( dimlen0 == dimlen )
        {
          dimid = dimid0;
          break;
        }
      iz++;
    }
  while ( iz <= 99 );


  if ( iz ) sprintf(dimname + len, "_%u", iz+1);

  return dimid;
}

static
383
void checkGridName(char *axisname, int fileID)
384
385
386
387
388
389
390
391
392
393
394
395
396
397
{
  int ncdimid;
  char axisname2[CDI_MAX_NAME];

  /* check that the name is not already defined */
  unsigned iz = 0;

  size_t axisnameLen = strlen(axisname);
  memcpy(axisname2, axisname, axisnameLen + 1);
  do
    {
      if ( iz ) sprintf(axisname2 + axisnameLen, "_%u", iz+1);

      int status = nc_inq_varid(fileID, axisname2, &ncdimid);
398
      if ( status != NC_NOERR ) break;
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451

      ++iz;
    }
  while ( iz <= 99 );

  if ( iz ) sprintf(axisname + axisnameLen, "_%u", iz+1);
}

static
int checkZaxisName(char *axisname, int fileID, int vlistID, int zaxisID, int nzaxis)
{
  char axisname2[CDI_MAX_NAME];

  /* check that the name is not already defined */
  unsigned iz = 0;

  size_t axisnameLen = strlen(axisname);
  memcpy(axisname2, axisname, axisnameLen + 1);
  do
    {
      if ( iz ) sprintf(axisname2 + axisnameLen, "_%u", iz+1);

      int ncdimid, status = nc_inq_varid(fileID, axisname2, &ncdimid);

      if ( status != NC_NOERR )
        {
          if ( iz )
            {
              /* check that the name does not exist for other zaxes */
              for ( int index = 0; index < nzaxis; index++ )
                {
                  int zaxisID0 = vlistZaxis(vlistID, index);
                  if ( zaxisID != zaxisID0 )
                    {
                      const char *axisname0 = zaxisInqNamePtr(zaxisID0);
                      if ( strcmp(axisname0, axisname2) == 0 ) goto nextSuffix;
                    }
                }
            }
          break;
        }
      nextSuffix:
      ++iz;
    }
  while (iz <= 99);


  if ( iz ) sprintf(axisname + axisnameLen, "_%u", iz+1);

  return (int)iz;
}

static void
452
cdfDefAxisCommon(stream_t *streamptr, int gridID, int gridindex, int ndims,
453
454
                 const struct cdfDefGridAxisInqs *gridAxisInq, int dimKey, char axisLetter,
                 void (*finishCyclicBounds)(double *pbounds, size_t dimlen, const double *pvals))
455
{
456
457
458
  int dimID = CDI_UNDEFID;
  int ncvarid = CDI_UNDEFID, ncbvarid = CDI_UNDEFID;
  int nvdimID = CDI_UNDEFID;
459
460
  int fileID  = streamptr->fileID;
  size_t dimlen = (size_t)gridAxisInq->axisSize(gridID);
461
462
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);

463
  ncgrid_t *ncgrid = streamptr->ncgrid;
464

465
466
467
468
  const double *pvals = gridAxisInq->axisValsPtr(gridID);
  char dimname[CDI_MAX_NAME+3]; dimname[0] = 0;
  if ( ndims && pvals == NULL ) cdiGridInqKeyStr(gridID, dimKey, CDI_MAX_NAME, dimname);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
  for ( int index = 0; index < gridindex; ++index )
470
    {
471
      int gridID0 = ncgrid[index].gridID;
472
      assert(gridID0 != CDI_UNDEFID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
474
475
476
477
478
      int gridtype0 = gridInqType(gridID0);
      if ( gridtype0 == GRID_GAUSSIAN    ||
           gridtype0 == GRID_LONLAT      ||
           gridtype0 == GRID_PROJECTION  ||
           gridtype0 == GRID_CURVILINEAR ||
           gridtype0 == GRID_GENERIC )
479
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
          size_t dimlen0 = (size_t)gridAxisInq->axisSize(gridID0);
481
482
483
484
          char dimname0[CDI_MAX_NAME]; dimname0[0] = 0;
          if ( dimname[0] ) cdiGridInqKeyStr(gridID0, dimKey, CDI_MAX_NAME, dimname0);
          bool lname = dimname0[0] ? strcmp(dimname, dimname0) == 0 : true;
          if ( dimlen == dimlen0 && lname )
485
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
487
488
              double (*inqVal)(int gridID, int index) = gridAxisInq->axisVal;
              if ( IS_EQUAL(inqVal(gridID0, 0), inqVal(gridID, 0)) &&
                   IS_EQUAL(inqVal(gridID0, (int)dimlen-1), inqVal(gridID, (int)dimlen-1)) )
489
                {
Thomas Jahns's avatar
Thomas Jahns committed
490
491
                  dimID = ncgrid[index].ncIDs[dimKey == CDI_KEY_XDIMNAME
                                              ? CDF_DIMID_X : CDF_DIMID_Y];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
492
                  break;
493
494
495
496
497
                }
            }
        }
    }

498
  if ( dimID == CDI_UNDEFID )
499
    {
500
501
502
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
      int keyname = (axisLetter == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
      gridAxisInq->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
503
504
505
      if ( axisname[0] == 0 ) Error("axis name undefined!");
      size_t axisnameLen = strlen(axisname);

506
507
508
      /* enough to append _ plus up to 100 decimal and trailing \0 */
      char extendedAxisname[axisnameLen + 4 + 1];
      memcpy(extendedAxisname, axisname, axisnameLen + 1);
509
      checkGridName(extendedAxisname, fileID);
510
      size_t extendedAxisnameLen = axisnameLen + strlen(extendedAxisname + axisnameLen);
511
512
513
514
515
516
517
518

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      if ( ndims )
        {
          if ( dimname[0] == 0 ) strcpy(dimname, extendedAxisname);
          dimID = checkDimName(fileID, dimlen, dimname);

519
          if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);
520
521
522
        }

      bool gen_bounds = false;
523
      bool grid_is_cyclic = gridIsCircular(gridID) > 0;
524
525
526
527
528
      double *pbounds = NULL;
      if ( pvals )
        {
          cdf_def_var(fileID, extendedAxisname, xtype, ndims, &dimID, &ncvarid);

529
          cdfPutGridStdAtts(fileID, ncvarid, gridID, axisLetter, gridAxisInq);
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
          {
            char axisStr[2] = { axisLetter, '\0' };
            cdf_put_att_text(fileID, ncvarid, "axis", 1, axisStr);
          }

          pbounds = (double *)gridAxisInq->axisBoundsPtr(gridID);

          if ( CDI_cmor_mode && grid_is_cyclic && !pbounds )
            {
              gen_bounds = true;
              pbounds = (double*) Malloc(2*dimlen*sizeof(double));
              for ( size_t i = 0; i < dimlen-1; ++i )
                {
                  pbounds[i*2+1]   = (pvals[i] + pvals[i+1])/2;
                  pbounds[(i+1)*2] = (pvals[i] + pvals[i+1])/2;
                }
              finishCyclicBounds(pbounds, dimlen, pvals);
            }
          if ( pbounds )
            {
              size_t nvertex = 2;
              if ( nc_inq_dimid(fileID, bndsName, &nvdimID) != NC_NOERR )
                cdf_def_dim(fileID, bndsName, nvertex, &nvdimID);
            }
554
          if ( pbounds && nvdimID != CDI_UNDEFID )
555
556
557
558
559
560
561
562
563
564
565
566
567
568
            {
              char boundsname[extendedAxisnameLen + 1 + sizeof (bndsName)];
              memcpy(boundsname, axisname, extendedAxisnameLen);
              boundsname[extendedAxisnameLen] = '_';
              memcpy(boundsname + extendedAxisnameLen + 1, bndsName, sizeof bndsName);
              int dimIDs[2] = { dimID, nvdimID };
              cdf_def_var(fileID, boundsname, xtype, 2, dimIDs, &ncbvarid);
              cdf_put_att_text(fileID, ncvarid, "bounds", extendedAxisnameLen + sizeof (bndsName), boundsname);
            }
        }

      cdf_enddef(fileID);
      streamptr->ncmode = 2;

569
570
      if ( ncvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncvarid, pvals);
      if ( ncbvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbvarid, pbounds);
571
572
      if ( gen_bounds ) Free(pbounds);

573
      if ( ndims == 0 )
574
575
        ncgrid[gridindex].ncIDs[dimKey == CDI_KEY_XDIMNAME
                                ? CDF_VARID_X : CDF_VARID_Y] = ncvarid;
576
577
    }

578
  ncgrid[gridindex].gridID = gridID;
Thomas Jahns's avatar
Thomas Jahns committed
579
580
  ncgrid[gridindex].ncIDs[dimKey == CDI_KEY_XDIMNAME
                          ? CDF_DIMID_X : CDF_DIMID_Y] = dimID;
581
582
}

583
584
static
void finishCyclicXBounds(double *pbounds, size_t dimlen, const double *pvals)
585
586
587
588
589
590
{
  pbounds[0] = (pvals[0] + pvals[dimlen-1]-360)*0.5;
  pbounds[2*dimlen-1] = (pvals[dimlen-1] + pvals[0]+360)*0.5;
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
void cdfDefXaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
592
{
593
594
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsX,
                   CDI_KEY_XDIMNAME, 'X', finishCyclicXBounds);
595
596
}

597
598
static
void finishCyclicYBounds(double *pbounds, size_t dimlen, const double *pvals)
599
600
601
602
603
604
{
  pbounds[0] = copysign(90.0, pvals[0]);
  pbounds[2*dimlen-1] = copysign(90.0, pvals[dimlen-1]);
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
void cdfDefYaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
606
{
607
608
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsY,
                   CDI_KEY_YDIMNAME, 'Y', finishCyclicYBounds);
609
610
611
612
613
614
}

static
void cdfGridCompress(int fileID, int ncvarid, int gridsize, int filetype, int comptype)
{
#if  defined  (HAVE_NETCDF4)
615
  if ( gridsize > 1 && comptype == CDI_COMPRESS_ZIP && (filetype == CDI_FILETYPE_NC4 || filetype == CDI_FILETYPE_NC4C) )
616
    {
617
      cdf_def_var_chunking(fileID, ncvarid, NC_CHUNKED, NULL);
618
619
620
621
622
623
      cdfDefVarDeflate(fileID, ncvarid, 1);
    }
#endif
}

static
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
void cdfDefGridReference(stream_t *streamptr, int gridID)
{
  int fileID  = streamptr->fileID;
  int number = gridInqNumber(gridID);

  if ( number > 0 )
    {
      cdf_put_att_int(fileID, NC_GLOBAL, "number_of_grid_used", NC_INT, 1, &number);
    }

  const char *gridfile = gridInqReferencePtr(gridID);
  if ( gridfile && gridfile[0] != 0 )
    cdf_put_att_text(fileID, NC_GLOBAL, "grid_file_uri", strlen(gridfile), gridfile);
}

static
void cdfDefGridUUID(stream_t *streamptr, int gridID)
{
  unsigned char uuidOfHGrid[CDI_UUID_SIZE];

  gridInqUUID(gridID, uuidOfHGrid);
  if ( !cdiUUIDIsNull(uuidOfHGrid) )
    {
      char uuidOfHGridStr[37];
      cdiUUID2Str(uuidOfHGrid, uuidOfHGridStr);
      if ( uuidOfHGridStr[0] != 0 && strlen(uuidOfHGridStr) == 36 )
        {
          int fileID  = streamptr->fileID;
          //if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
          cdf_put_att_text(fileID, NC_GLOBAL, "uuidOfHGrid", 36, uuidOfHGridStr);
          //if ( streamptr->ncmode == 2 ) cdf_enddef(fileID);
        }
    }
}

struct cdfDefIrregularGridCommonIDs
{
  int xdimID, ydimID, ncxvarid, ncyvarid, ncavarid;
};

static struct cdfDefIrregularGridCommonIDs
cdfDefIrregularGridCommon(stream_t *streamptr, int gridID,
                          size_t xdimlen, size_t ydimlen,
                          int ndims, const char *xdimname_default,
                          size_t nvertex, const char *vdimname_default,
                          bool setVdimname)
670
{
671
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
  int xdimID = CDI_UNDEFID;
  int ydimID = CDI_UNDEFID;
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
  int ncbxvarid = CDI_UNDEFID, ncbyvarid = CDI_UNDEFID;
  int fileID  = streamptr->fileID;
  if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

  {
    char xdimname[CDI_MAX_NAME+3];
    xdimname[0] = 0;
    cdiGridInqKeyStr(gridID, CDI_KEY_XDIMNAME, CDI_MAX_NAME, xdimname);
    if ( xdimname[0] == 0 ) strcpy(xdimname, xdimname_default);
    xdimID = checkDimName(fileID, xdimlen, xdimname);
    if ( xdimID == CDI_UNDEFID ) cdf_def_dim(fileID, xdimname, xdimlen, &xdimID);
  }

  if ( ndims == 3 )
    {
      char ydimname[CDI_MAX_NAME+3];
      ydimname[0] = 0;
      cdiGridInqKeyStr(gridID, CDI_KEY_YDIMNAME, CDI_MAX_NAME, ydimname);
      if ( ydimname[0] == 0 ) { ydimname[0] = 'y'; ydimname[1] = 0; }
      ydimID = checkDimName(fileID, ydimlen, ydimname);
      if ( ydimID == CDI_UNDEFID ) cdf_def_dim(fileID, ydimname, ydimlen, &ydimID);
    }

  int nvdimID = CDI_UNDEFID;
  int dimIDs[3];
  dimIDs[ndims-1] = CDI_UNDEFID;
  if ( setVdimname )
    {
      char vdimname[CDI_MAX_NAME+3]; vdimname[0] = 0;
      cdiGridInqKeyStr(gridID, CDI_KEY_VDIMNAME, CDI_MAX_NAME, vdimname);
      if ( vdimname[0] == 0 ) strcpy(vdimname, vdimname_default);
      nvdimID = dimIDs[ndims-1] = checkDimName(fileID, nvertex, vdimname);
      if ( nvdimID == CDI_UNDEFID )
        {
          cdf_def_dim(fileID, vdimname, nvertex, dimIDs+ndims-1);
          nvdimID = dimIDs[ndims-1];
        }
    }

  if ( ndims == 3 )
    {
      dimIDs[0] = ydimID;
      dimIDs[1] = xdimID;
    }
  else /* ndims == 2 */
    {
      dimIDs[0] = xdimID;
      cdfDefGridReference(streamptr, gridID);
      cdfDefGridUUID(streamptr, gridID);
    }

  const double *xvalsPtr = gridInqXvalsPtr(gridID),
    *xboundsPtr = NULL;
  if ( xvalsPtr )
    {
      char xaxisname[CDI_MAX_NAME]; xaxisname[0] = 0;
      cdiGridInqKeyStr(gridID, CDI_KEY_XNAME, CDI_MAX_NAME, xaxisname);
      checkGridName(xaxisname, fileID);
      cdf_def_var(fileID, xaxisname, xtype, ndims-1, dimIDs, &ncxvarid);
      cdfGridCompress(fileID, ncxvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

      cdfPutGridStdAtts(fileID, ncxvarid, gridID, 'X', &gridInqsX);

      /* attribute for Panoply */
      if ( ndims == 3 )
        cdf_put_att_text(fileID, ncxvarid, "_CoordinateAxisType", 3, "Lon");

      if ( (xboundsPtr = gridInqXboundsPtr(gridID)) && nvdimID != CDI_UNDEFID )
        {
          size_t xaxisnameLen = strlen(xaxisname);
          xaxisname[xaxisnameLen] = '_';
          memcpy(xaxisname + xaxisnameLen + 1, bndsName, sizeof (bndsName));
          cdf_def_var(fileID, xaxisname, xtype, ndims, dimIDs, &ncbxvarid);
          cdfGridCompress(fileID, ncbxvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

          cdf_put_att_text(fileID, ncxvarid, "bounds", xaxisnameLen + sizeof (bndsName), xaxisname);
        }
    }

  const double *yvalsPtr = gridInqYvalsPtr(gridID),
    *yboundsPtr = NULL;
  if ( yvalsPtr )
    {
      char yaxisname[CDI_MAX_NAME];
      gridInqYname(gridID, yaxisname);
      checkGridName(yaxisname, fileID);

      cdf_def_var(fileID, yaxisname, xtype, ndims - 1, dimIDs, &ncyvarid);
      cdfGridCompress(fileID, ncyvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

      cdfPutGridStdAtts(fileID, ncyvarid, gridID, 'Y', &gridInqsY);

      /* attribute for Panoply */
      if ( ndims == 3 )
        cdf_put_att_text(fileID, ncyvarid, "_CoordinateAxisType", 3, "Lat");

      if ( (yboundsPtr = gridInqYboundsPtr(gridID)) && nvdimID != CDI_UNDEFID )
        {
          size_t yaxisnameLen = strlen(yaxisname);
          yaxisname[yaxisnameLen] = '_';
          memcpy(yaxisname + yaxisnameLen + 1, bndsName, sizeof (bndsName));
          cdf_def_var(fileID, yaxisname, xtype, ndims, dimIDs, &ncbyvarid);
          cdfGridCompress(fileID, ncbyvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

          cdf_put_att_text(fileID, ncyvarid, "bounds", yaxisnameLen + sizeof (bndsName), yaxisname);
        }
    }

  const double *areaPtr = gridInqAreaPtr(gridID);
  if ( areaPtr )
    {
      static const char yaxisname_[] = "cell_area";
      static const char units[] = "m2";
      static const char longname[] = "area of grid cell";
      static const char stdname[] = "cell_area";

      cdf_def_var(fileID, yaxisname_, xtype, ndims-1, dimIDs, &ncavarid);

      cdf_put_att_text(fileID, ncavarid, "standard_name", sizeof (stdname) - 1, stdname);
      cdf_put_att_text(fileID, ncavarid, "long_name", sizeof (longname) - 1, longname);
      cdf_put_att_text(fileID, ncavarid, "units", sizeof (units) - 1, units);
    }

  cdf_enddef(fileID);
  streamptr->ncmode = 2;

  if ( ncxvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncxvarid,  xvalsPtr);
  if ( ncbxvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbxvarid, xboundsPtr);
  if ( ncyvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncyvarid,  yvalsPtr);
  if ( ncbyvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbyvarid, yboundsPtr);
  if ( ncavarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncavarid,  areaPtr);

  return (struct cdfDefIrregularGridCommonIDs) {
    .xdimID=xdimID, .ydimID = ydimID,
    .ncxvarid=ncxvarid, .ncyvarid=ncyvarid, .ncavarid=ncavarid
  };
}

static
void cdfDefCurvilinear(stream_t *streamptr, int gridID, int gridindex)
{
816
  ncgrid_t *ncgrid = streamptr->ncgrid;
817
818
819
820
821

  size_t dimlen = (size_t)gridInqSize(gridID);
  size_t xdimlen = (size_t)gridInqXsize(gridID);
  size_t ydimlen = (size_t)gridInqYsize(gridID);

822
  int xdimID = CDI_UNDEFID, ydimID = CDI_UNDEFID;
823
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
824
825
826
827
828
829
830
831
832
  {
    size_t ofs = 0;
    do {
      struct idSearch search
        = cdfSearchIDBySize(ofs, (size_t)gridindex, ncgrid, CDF_DIMID_X,
                            GRID_CURVILINEAR, (int)dimlen,
                            gridInqType, gridInqSize);
      size_t index = search.foundIdx;
      if ( index != SIZE_MAX )
833
        {
834
          int gridID0 = ncgrid[index].gridID;
835
836
837
838
839
840
          if (    IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0))
               && IS_EQUAL(gridInqXval(gridID0, (int)dimlen-1),
                           gridInqXval(gridID, (int)dimlen-1))
               && IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0))
               && IS_EQUAL(gridInqYval(gridID0, (int)dimlen-1),
                           gridInqYval(gridID, (int)dimlen-1)) )
841
            {
842
843
844
845
846
              xdimID = ncgrid[index].ncIDs[CDF_DIMID_X];
              ydimID = ncgrid[index].ncIDs[CDF_DIMID_Y];
              ncxvarid = ncgrid[index].ncIDs[CDF_VARID_X];
              ncyvarid = ncgrid[index].ncIDs[CDF_VARID_Y];
              break;
847
            }
848
849
850
          ofs = search.foundIdx;
          if ( ofs < (size_t)gridindex )
            continue;
851
        }
852
853
    } while (false);
  }
854

855
  if ( xdimID == CDI_UNDEFID || ydimID == CDI_UNDEFID )
856
    {
857
858
859
860
861
862
863
864
865
866
      struct cdfDefIrregularGridCommonIDs createdIDs
        = cdfDefIrregularGridCommon(streamptr, gridID,
                                    xdimlen, ydimlen, 3, "x", 4, "nv4",
                                    gridInqXboundsPtr(gridID)
                                    || gridInqYboundsPtr(gridID));
      xdimID = createdIDs.xdimID;
      ydimID = createdIDs.ydimID;
      ncxvarid = createdIDs.ncxvarid;
      ncyvarid = createdIDs.ncyvarid;
      ncavarid = createdIDs.ncavarid;
867
868
    }

869
  ncgrid[gridindex].gridID = gridID;
870
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = xdimID;
871
  ncgrid[gridindex].ncIDs[CDF_DIMID_Y] = ydimID;
872
  ncgrid[gridindex].ncIDs[CDF_VARID_X] = ncxvarid;
873
  ncgrid[gridindex].ncIDs[CDF_VARID_Y] = ncyvarid;
874
  ncgrid[gridindex].ncIDs[CDF_VARID_A] = ncavarid;
875
876
877
878
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
879
void cdfDefUnstructured(stream_t *streamptr, int gridID, int gridindex)
880
{
881
  ncgrid_t *ncgrid = streamptr->ncgrid;
882
883
884

  size_t dimlen = (size_t)gridInqSize(gridID);

885
  int dimID = CDI_UNDEFID;
886
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
887
888
889
890
891
892
893
894
895
  {
    size_t ofs = 0;
    do {
      struct idSearch search
        = cdfSearchIDBySize(ofs, (size_t)gridindex, ncgrid, CDF_DIMID_X,
                            GRID_UNSTRUCTURED, (int)dimlen,
                            gridInqType, gridInqSize);
      size_t index = search.foundIdx;
      if ( index != SIZE_MAX )
896
        {
897
          int gridID0 = ncgrid[index].gridID;
898
899
900
901
902
903
904
          if ( gridInqNvertex(gridID0) == gridInqNvertex(gridID) &&
               IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0)) &&
               IS_EQUAL(gridInqXval(gridID0, (int)dimlen-1),
                        gridInqXval(gridID, (int)dimlen-1)) &&
               IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0)) &&
               IS_EQUAL(gridInqYval(gridID0, (int)dimlen-1),
                        gridInqYval(gridID, (int)dimlen-1)) )
905
            {
906
907
908
909
910
              dimID = ncgrid[index].ncIDs[CDF_DIMID_X];
              ncxvarid = ncgrid[index].ncIDs[CDF_VARID_X];
              ncyvarid = ncgrid[index].ncIDs[CDF_VARID_Y];
              ncavarid = ncgrid[index].ncIDs[CDF_VARID_A];
              break;
911
            }
912
          ofs = search.foundIdx;
913
          if ( ofs < (size_t)gridindex )
914
            continue;
915
        }
916
917
    } while (false);
  }
918

919
  if ( dimID == CDI_UNDEFID )
920
921
    {
      size_t nvertex = (size_t)gridInqNvertex(gridID);
922
923
924
925
926
927
928
929
      struct cdfDefIrregularGridCommonIDs createdIDs
        = cdfDefIrregularGridCommon(streamptr, gridID,
                                    dimlen, 1, 2, "ncells",
                                    nvertex, "vertices", nvertex > 0);
      dimID = createdIDs.xdimID;
      ncxvarid = createdIDs.ncxvarid;
      ncyvarid = createdIDs.ncyvarid;
      ncavarid = createdIDs.ncavarid;
930
931
    }

932
  ncgrid[gridindex].gridID = gridID;
933
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = dimID;
934
  ncgrid[gridindex].ncIDs[CDF_VARID_X] = ncxvarid;
935
  ncgrid[gridindex].ncIDs[CDF_VARID_Y] = ncyvarid;
936
  ncgrid[gridindex].ncIDs[CDF_VARID_A] = ncavarid;
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
}

struct attTxtTab2
{
  const char *attName, *attVal;
  size_t valLen;
};

static
void cdf_def_vct_echam(stream_t *streamptr, int zaxisID)
{
  int type = zaxisInqType(zaxisID);

  if ( type == ZAXIS_HYBRID || type == ZAXIS_HYBRID_HALF )
    {
      int ilev = zaxisInqVctSize(zaxisID)/2;
      if ( ilev == 0 ) return;

      int mlev = ilev - 1;

      if ( streamptr->vct.ilev > 0 )
        {
          if ( streamptr->vct.ilev != ilev )
960
            Error("More than one VCT for each file unsupported!");
961
962
963
964
965
966
967
          return;
        }

      int fileID = streamptr->fileID;

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

Thomas Jahns's avatar
Thomas Jahns committed
968
      int ncdimid = -1, ncdimid2 = -1;
969
970
      int hyaiid, hybiid, hyamid = -1, hybmid = -1;

971
      cdf_def_dim(fileID, "nhyi", (size_t)ilev, &ncdimid2);
972
973
974
975
976
977
978
979
      cdf_def_var(fileID, "hyai", NC_DOUBLE, 1, &ncdimid2, &hyaiid);
      cdf_def_var(fileID, "hybi", NC_DOUBLE, 1, &ncdimid2, &hybiid);
      if ( mlev > 0 )
        {
          cdf_def_dim(fileID, "nhym", (size_t)mlev, &ncdimid);
          cdf_def_var(fileID, "hyam", NC_DOUBLE, 1, &ncdimid,  &hyamid);
          cdf_def_var(fileID, "hybm", NC_DOUBLE, 1, &ncdimid,  &hybmid);
        }
980
981

      streamptr->vct.ilev   = ilev;
982
      streamptr->vct.mlev   = mlev;
983
984
985
986
987
988
      streamptr->vct.mlevID = ncdimid;
      streamptr->vct.ilevID = ncdimid2;

      {
        static const char lname_n[] = "long_name",
          units_n[] = "units",
989
          lname_v_ai[] = "hybrid A coefficient at layer interfaces",
990
991
          units_v_ai[] = "Pa",
          lname_v_bi[] = "hybrid B coefficient at layer interfaces",
992
          units_v_bi[] = "1";
993
994
995
996
997
998
        static const struct attTxtTab2 tab[]
          = {
          { lname_n, lname_v_ai, sizeof (lname_v_ai) - 1 },
          { units_n, units_v_ai, sizeof (units_v_ai) - 1 },
          { lname_n, lname_v_bi, sizeof (lname_v_bi) - 1 },
          { units_n, units_v_bi, sizeof (units_v_bi) - 1 },
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
        };
        enum { tabLen = sizeof (tab) / sizeof (tab[0]) };
        int ids[tabLen] = { hyaiid, hyaiid, hybiid, hybiid };
        for ( size_t i = 0; i < tabLen; ++i )
          cdf_put_att_text(fileID, ids[i], tab[i].attName, tab[i].valLen, tab[i].attVal);
      }

      {
        static const char lname_n[] = "long_name",
          units_n[] = "units",
          lname_v_am[] = "hybrid A coefficient at layer midpoints",
          units_v_am[] = "Pa",
          lname_v_bm[] = "hybrid B coefficient at layer midpoints",
          units_v_bm[] = "1";
        static const struct attTxtTab2 tab[]
          = {
1015
1016
1017
1018
1019
1020
          { lname_n, lname_v_am, sizeof (lname_v_am) - 1 },
          { units_n, units_v_am, sizeof (units_v_am) - 1 },
          { lname_n, lname_v_bm, sizeof (lname_v_bm) - 1 },
          { units_n, units_v_bm, sizeof (units_v_bm) - 1 },
        };
        enum { tabLen = sizeof (tab) / sizeof (tab[0]) };
1021
        int ids[tabLen] = { hyamid, hyamid, hybmid, hybmid };
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
        for ( size_t i = 0; i < tabLen; ++i )
          cdf_put_att_text(fileID, ids[i], tab[i].attName, tab[i].valLen, tab[i].attVal);
      }

      cdf_enddef(fileID);
      streamptr->ncmode = 2;

      const double *vctptr = zaxisInqVctPtr(zaxisID);

      cdf_put_var_double(fileID, hyaiid, vctptr);
      cdf_put_var_double(fileID, hybiid, vctptr+ilev);

1034
1035
1036
      size_t start;
      size_t count = 1;
      double mval;
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
      for ( int i = 0; i < mlev; i++ )
        {
          start = (size_t)i;
          mval = (vctptr[i] + vctptr[i+1]) * 0.5;
          cdf_put_vara_double(fileID, hyamid, &start, &count, &mval);
          mval = (vctptr[ilev+i] + vctptr[ilev+i+1]) * 0.5;
          cdf_put_vara_double(fileID, hybmid, &start, &count, &mval);
        }
    }
}

static
1049
void cdf_def_vct_cf(stream_t *streamptr, int zaxisID, int nclevID, int ncbndsID, int p0status, double p0value)
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
{
  int type = zaxisInqType(zaxisID);

  if ( type == ZAXIS_HYBRID || type == ZAXIS_HYBRID_HALF )
    {
      int ilev = zaxisInqVctSize(zaxisID)/2;
      if ( ilev == 0 ) return;

      int mlev = ilev - 1;
      int hyaiid = 0, hybiid = 0, hyamid, hybmid;

      if ( streamptr->vct.ilev > 0 )
        {
          if ( streamptr->vct.ilev != ilev )
            Error("more than one VCT for each file unsupported!");
          return;
        }

      int fileID = streamptr->fileID;

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      int dimIDs[2];
      dimIDs[0] = nclevID;
      dimIDs[1] = ncbndsID;

      streamptr->vct.mlev   = mlev;
      streamptr->vct.ilev   = ilev;
      streamptr->vct.mlevID = nclevID;
      streamptr->vct.ilevID = nclevID;

1081
1082
1083
1084
      if ( p0status == 0 )
        cdf_def_var(fileID, "a", NC_DOUBLE, 1, dimIDs,  &hyamid);
      else
        cdf_def_var(fileID, "ap", NC_DOUBLE, 1, dimIDs,  &hyamid);
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
      cdf_def_var(fileID, "b",  NC_DOUBLE, 1, dimIDs,  &hybmid);

      {
        static const char lname[] = "vertical coordinate formula term: ap(k)";
        cdf_put_att_text(fileID, hyamid, "long_name", sizeof (lname) - 1, lname);
      }
      {
        static const char units[] = "Pa";
        cdf_put_att_text(fileID, hyamid, "units", sizeof (units) - 1, units);
      }
      {
        static const char lname[] = "vertical coordinate formula term: b(k)";
        cdf_put_att_text(fileID, hybmid, "long_name", sizeof (lname) - 1, lname);
      }
      {
        static const char units[] = "1";
        cdf_put_att_text(fileID, hybmid, "units", sizeof (units) - 1, units);
      }

      if ( ncbndsID != -1 )
        {
1106
1107
1108
1109
          if ( p0status == 0 )
            cdf_def_var(fileID, "a_bnds", NC_DOUBLE, 2, dimIDs, &hyaiid);
          else
            cdf_def_var(fileID, "ap_bnds", NC_DOUBLE, 2, dimIDs, &hyaiid);
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
          cdf_def_var(fileID, "b_bnds",  NC_DOUBLE, 2, dimIDs, &hybiid);
          {
            static const char lname[] = "vertical coordinate formula term: ap(k+1/2)";
            cdf_put_att_text(fileID, hyaiid, "long_name", sizeof (lname) - 1, lname);
          }
          {
            static const char units[] = "Pa";
            cdf_put_att_text(fileID, hyaiid, "units", sizeof (units) - 1, units);
          }
          {
            static const char lname[] = "vertical coordinate formula term: b(k+1/2)";
            cdf_put_att_text(fileID, hybiid, "long_name", sizeof (lname) - 1, lname);
          }
          {
            static const char units[] = "1";
            cdf_put_att_text(fileID, hybiid, "units", sizeof (units) - 1, units);
          }
        }

      cdf_enddef(fileID);
      streamptr->ncmode = 2;

1132
1133
1134
1135
1136
1137
1138
      int vctsize = zaxisInqVctSize(zaxisID);
      double vct[vctsize];
      zaxisInqVct(zaxisID, vct);

      if ( p0status == 0 && IS_NOT_EQUAL(p0value,0) )
        for ( int i = 0; i < vctsize/2; ++i ) vct[i] /= p0value;

1139
1140
1141
1142
1143
1144
      double tarray[ilev*2];

      if ( ncbndsID != -1 )
        {
          for ( int i = 0; i < mlev; ++i )
            {
1145
1146
              tarray[2*i  ] = vct[i];
              tarray[2*i+1] = vct[i+1];
1147
1148
1149
1150
1151
            }
          cdf_put_var_double(fileID, hyaiid, tarray);

          for ( int i = 0; i < mlev; ++i )
            {
1152
1153
              tarray[2*i  ] = vct[ilev+i];
              tarray[2*i+1] = vct[ilev+i+1];
1154
1155
1156
1157
1158
            }
          cdf_put_var_double(fileID, hybiid, tarray);
        }

      for ( int i = 0; i < mlev; ++i )
1159
        tarray[i] = (vct[i] + vct[i+1]) * 0.5;
1160
1161
1162
      cdf_put_var_double(fileID, hyamid, tarray);

      for ( int i = 0; i < mlev; ++i )
1163
        tarray[i] = (vct[ilev+i] + vct[ilev+i+1]) * 0.5;
1164
1165
1166
1167
1168
1169
1170
      cdf_put_var_double(fileID, hybmid, tarray);
    }
}

struct attTxtTab { const char *txt; size_t txtLen; };

static
1171
void cdf_def_zaxis_hybrid_echam(stream_t *streamptr, int type, int *ncvaridp, int zaxisID, int zaxisindex, int xtype, size_t dimlen, int *dimID, char *axisname)
1172
{
1173
  int fileID = streamptr->fileID;
1174
1175
1176
1177

  if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

  cdf_def_dim(fileID, axisname, dimlen, dimID);
1178
1179
  cdf_def_var(fileID, axisname, (nc_type) xtype, 1, dimID,  ncvaridp);
  int ncvarid = *ncvaridp;
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232

  {
    static const char sname[] = "hybrid_sigma_pressure";
    cdf_put_att_text(fileID, ncvarid, "standard_name", sizeof (sname) - 1, sname);
  }
  {
    static const char *attName[] = {
      "long_name",
      "formula",
      "formula_terms"
    };
    enum { nAtt = sizeof (attName) / sizeof (attName[0]) };
    static const char lname_m[] = "hybrid level at layer midpoints",
      formula_m[] = "hyam hybm (mlev=hyam+hybm*aps)",
      fterms_m[] = "ap: hyam b: hybm ps: aps",
      lname_i[] = "hybrid level at layer interfaces",
      formula_i[] = "hyai hybi (ilev=hyai+hybi*aps)",
      fterms_i[] = "ap: hyai b: hybi ps: aps";
    static const struct attTxtTab tab[2][nAtt] = {
      {
        { lname_i, sizeof (lname_i) - 1 },
        { formula_i, sizeof (formula_i) - 1 },
        { fterms_i, sizeof (fterms_i) - 1 }
      },
      {
        { lname_m, sizeof (lname_m) - 1 },
        { formula_m, sizeof (formula_m) - 1 },
        { fterms_m, sizeof (fterms_m) - 1 }
      }
    };

    size_t tabSelect = type == ZAXIS_HYBRID;
    for (size_t i = 0; i < nAtt; ++i)
      cdf_put_att_text(fileID, ncvarid, attName[i],
                       tab[tabSelect][i].txtLen, tab[tabSelect][i].txt);
  }

  {
    static const char units[] = "level";
    cdf_put_att_text(fileID, ncvarid, "units", sizeof (units) - 1, units);
  }
  {
    static const char direction[] = "down";
    cdf_put_att_text(fileID, ncvarid, "positive", sizeof (direction) - 1, direction);
  }

  cdf_enddef(fileID);
  streamptr->ncmode = 2;

  cdf_put_var_double(fileID, ncvarid, zaxisInqLevelsPtr(zaxisID));

  cdf_def_vct_echam(streamptr, zaxisID);

1233
  if ( *dimID == CDI_UNDEFID )
Thomas Jahns's avatar
Thomas Jahns committed
1234
1235
    streamptr->zaxisID[zaxisindex] = type == ZAXIS_HYBRID
      ? streamptr->vct.mlevID : streamptr->vct.ilevID;
1236
1237
1238
}

static
1239
void cdf_def_zaxis_hybrid_cf(stream_t *streamptr, int type, int *ncvaridp, int zaxisID, int zaxisindex, int xtype, size_t dimlen, int *dimID, char *axisname)
1240
{
1241
1242
1243
  int fileID = streamptr->fileID;
  if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

1244
1245
  char psname[CDI_MAX_NAME]; psname[0] = 0;
  cdiZaxisInqKeyStr(zaxisID, CDI_KEY_PSNAME, CDI_MAX_NAME, psname);
1246
1247
  if ( psname[0] == 0 ) strcpy(psname, "ps");

1248
1249
  char p0name[CDI_MAX_NAME]; p0name[0] = 0;
  double p0value = 1;
1250
  int p0varid = CDI_UNDEFID;
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
  int p0status = cdiZaxisInqKeyFlt(zaxisID, CDI_KEY_P0VALUE, &p0value);
  if ( p0status == 0 )
    {
      cdiZaxisInqKeyStr(zaxisID, CDI_KEY_P0NAME, CDI_MAX_NAME, p0name);
      if ( p0name[0] == 0 ) strcpy(p0name, "p0");
      cdf_def_var(fileID, p0name, NC_DOUBLE, 0, 0,  &p0varid);
      static const char longname[] = "reference pressure";
      cdf_put_att_text(fileID, p0varid, "long_name", strlen(longname), longname);
      static const char units[] = "Pa";
      cdf_put_att_text(fileID, p0varid, "units", strlen(units), units);
    }
1262

1263
1264
1265
1266
1267
1268
1269
1270
1271
  char zname[CDI_MAX_NAME]; zname[0] = 0;
  char zlongname[CDI_MAX_NAME]; zlongname[0] = 0;
  char zunits[CDI_MAX_NAME]; zunits[0] = 0;
  cdiZaxisInqKeyStr(zaxisID, CDI_KEY_NAME, CDI_MAX_NAME, zname);
  //cdiZaxisInqKeyStr(zaxisID, CDI_KEY_LONGNAME, CDI_MAX_NAME, zlongname);
  cdiZaxisInqKeyStr(zaxisID, CDI_KEY_UNITS, CDI_MAX_NAME, zunits);
  if ( zname[0] ) strcpy(axisname, zname);
  if ( zlongname[0] == 0 ) strcpy(zlongname, "hybrid sigma pressure coordinate");
  if ( zunits[0] == 0 ) strcpy(zunits, "1");
1272
1273

  cdf_def_dim(fileID, axisname, dimlen, dimID);
1274
1275
  cdf_def_var(fileID, axisname, (nc_type) xtype, 1, dimID, ncvaridp);
  int ncvarid = *ncvaridp;
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289

  {
    static const char sname[] = "standard_name",
      sname_v[] = "atmosphere_hybrid_sigma_pressure_coordinate",
      axis[] = "axis",
      axis_v[] = "Z",
      direction[] = "positive",
      direction_v[] = "down";
    struct attTxtTab2 tab[] = {
      { sname, sname_v, sizeof (sname_v) - 1 },
      { axis, axis_v, sizeof (axis_v) - 1 },
      { direction, direction_v, sizeof (direction_v) - 1 },
    };
    enum { nAtt = sizeof (tab) / sizeof (tab[0]) };
1290
    for ( size_t i = 0; i < nAtt; ++i )
1291
      cdf_put_att_text(fileID, ncvarid, tab[i].attName, tab[i].valLen, tab[i].attVal);
1292
1293
1294

    cdf_put_att_text(fileID, ncvarid, "long_name", strlen(zlongname), zlongname);
    cdf_put_att_text(fileID, ncvarid, "units", strlen(zunits), zunits);
1295
  }
1296
1297
1298
1299
1300
1301
1302
1303

  size_t len = 0;
  char txt[CDI_MAX_NAME];
  if ( p0status == 0 )
    len = (size_t)(sprintf(txt, "%s%s %s%s", "a: a b: b p0: ", p0name, "ps: ", psname));
  else
    len = (size_t)(sprintf(txt, "%s%s", "ap: ap b: b ps: ", psname));
  cdf_put_att_text(fileID, ncvarid, "formula_terms", len, txt);
1304

1305
1306
  int ncbvarid = CDI_UNDEFID;
  int nvdimID = CDI_UNDEFID;
1307
1308
1309

  double lbounds[dimlen], ubounds[dimlen], levels[dimlen];

1310
1311
1312
1313
  if ( zaxisInqLevels(zaxisID, NULL) )
    zaxisInqLevels(zaxisID, levels);
  else
    for ( size_t i = 0; i < dimlen; ++i ) levels[i] = i+1;
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329

  if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
    {
      zaxisInqLbounds(zaxisID, lbounds);
      zaxisInqUbounds(zaxisID, ubounds);
    }
  else
    {
      for ( size_t i = 0; i < dimlen; ++i ) lbounds[i] = levels[i];
      for ( size_t i = 0; i < dimlen-1; ++i ) ubounds[i] = levels[i+1];
      ubounds[dimlen-1] = levels[dimlen-1] + 1;
    }

  //if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
    {
      size_t nvertex = 2;
1330
      if ( dimlen > 1 && nc_inq_dimid(fileID, bndsName, &nvdimID) != NC_NOERR )
1331
1332
        cdf_def_dim(fileID, bndsName, nvertex, &nvdimID);

1333
      if ( nvdimID != CDI_UNDEFID )
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
        {
          size_t axisnameLen = strlen(axisname);
          axisname[axisnameLen] = '_';
          memcpy(axisname + axisnameLen + 1, bndsName, sizeof (bndsName));
          axisnameLen += sizeof (bndsName);
          int dimIDs[2] = { *dimID, nvdimID };
          cdf_def_var(fileID, axisname, (nc_type) xtype, 2, dimIDs, &ncbvarid);
          cdf_put_att_text(fileID, ncvarid, "bounds", axisnameLen, axisname);
          {
            static const char sname[] = "standard_name",
1344
              sname_v[] = "atmosphere_hybrid_sigma_pressure_coordinate";
1345
1346
1347
1348
            struct attTxtTab2 tab[] = {
              { sname, sname_v, sizeof (sname_v) - 1 },
            };
            enum { nAtt = sizeof (tab) / sizeof (tab[0]) };
1349
            for ( size_t i = 0; i < nAtt; ++i )
1350
              cdf_put_att_text(fileID, ncbvarid, tab[i].attName, tab[i].valLen, tab[i].attVal);
1351
            cdf_put_att_text(fileID, ncbvarid, "units", strlen(zunits), zunits);
1352
          }
1353
1354
1355
1356
1357
1358

          if ( p0status == 0 )
            len = (size_t)(sprintf(txt, "%s%s %s%s", "a: a_bnds b: b_bnds p0: ", p0name, "ps: ", psname));
          else
            len = (size_t)(sprintf(txt, "%s%s", "ap: ap_bnds b: b_bnds ps: ", psname));
          cdf_put_att_text(fileID, ncbvarid, "formula_terms", len, txt);
1359
1360
1361
1362
1363
1364
1365
1366
        }
    }

  cdf_enddef(fileID);
  streamptr->ncmode = 2;

  cdf_put_var_double(fileID, ncvarid, levels);

1367
  if ( p0varid != CDI_UNDEFID ) cdf_put_var_double(fileID, p0varid, &p0value);
1368