stream_cdf_o.c 65.8 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
#ifdef HAVE_CONFIG_H
2
#include "config.h"
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#endif

#ifdef HAVE_LIBNETCDF

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


#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2


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


23
24
25
26
27
28
29
30
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);
31
  size_t datasize = gridInqSize(gridID);
32
  int datatype = vlistInqVarDatatype(vlistID1, ivarID);
33
  int memtype  = datatype != CDI_DATATYPE_FLT32 ? MEMTYPE_DOUBLE : MEMTYPE_FLOAT;
34

35
  void *data = Malloc(datasize * (memtype == MEMTYPE_DOUBLE ? sizeof(double) : sizeof(float)));
36

37
  size_t nmiss;
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
  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)
{
53
  const int fileID = streamptr->fileID;
54

55
  if ( CDI_Debug ) Message("streamID = %d, fileID = %d", streamptr->self, fileID);
56
57
58
59
60
61
62
63
64
65
66
67
68

  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;
69
  size_t index = (size_t)tsID;
70
71
72
73
74
  cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);

  if ( taxis->has_bounds )
    {
      ncvarid = streamptr->basetime.ncvarboundsid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
      if ( ncvarid == CDI_UNDEFID ) Error("Call to taxisWithBounds() missing!");
76
77

      timevalue = cdiEncodeTimeval(taxis->vdate_lb, taxis->vtime_lb, &streamptr->tsteps[0].taxis);
78
      size_t start[2], count[2];
79
80
81
82
83
84
85
86
87
      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;
88
  if ( taxis->type == TAXIS_FORECAST && ncvarid != CDI_UNDEFID )
89
90
91
92
93
94
95
96
97
98
99
100
    {
      timevalue = taxis->fc_period;
      cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);
    }
}

void cdfDefTimestep(stream_t *streamptr, int tsID)
{
  cdfDefTimeValue(streamptr, tsID);
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
void cdfDefComplex(stream_t *streamptr, int gridID, int gridindex)
102
{
103
  int dimID;
104
  ncgrid_t *ncgrid = streamptr->ncgrid;
105

Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
  for ( int index = 0; index < gridindex; ++index )
107
    {
108
      if ( ncgrid[index].ncIDs[CDF_DIMID_X] != CDI_UNDEFID )
109
        {
110
111
          const int gridID0 = ncgrid[index].gridID;
          const int gridtype0 = gridInqType(gridID0);
112
113
          if ( gridtype0 == GRID_SPECTRAL || gridtype0 == GRID_FOURIER )
            {
114
              dimID = ncgrid[index].ncIDs[CDF_DIMID_X];
115
              goto dimIDEstablished;
116
117
118
119
            }
        }
    }

120
121
  {
    static const char axisname[] = "nc2";
122
123
    const size_t dimlen = 2;
    const int fileID  = streamptr->fileID;
124
125
126
127
128
129
130

    if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
    cdf_def_dim(fileID, axisname, dimlen, &dimID);
    cdf_enddef(fileID);
    streamptr->ncmode = 2;
  }
  dimIDEstablished:
131
  ncgrid[gridindex].gridID = gridID;
132
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = dimID;
133
134
}

135
struct idSearch
136
{
137
138
139
  int numNonMatching, foundID;
  size_t foundIdx;
};
140

141
static inline struct idSearch
142
cdfSearchIDBySize(size_t startIdx, size_t numIDs, const ncgrid_t ncgrid[/*numIDs*/],
143
                  int ncIDType, int searchType, int searchSize,
144
                  int (*typeInq)(int id), size_t (*sizeInq)(int id))
145
146
147
148
149
{
  int numNonMatching = 0,
    foundID = CDI_UNDEFID;
  size_t foundIdx = SIZE_MAX;
  for ( size_t index = startIdx; index < numIDs; index++ )
150
    {
151
      if ( ncgrid[index].ncIDs[ncIDType] != CDI_UNDEFID )
152
        {
153
154
155
          int id0 = ncgrid[index].gridID,
            id0Type = typeInq(id0);
          if ( id0Type == searchType )
156
            {
157
158
              int size0 = sizeInq(id0);
              if ( searchSize == size0 )
159
                {
160
161
                  foundID = ncgrid[index].ncIDs[ncIDType];
                  foundIdx = index;
162
163
                  break;
                }
164
              numNonMatching++;
165
166
167
            }
        }
    }
168
169
170
171
  return (struct idSearch){ .numNonMatching = numNonMatching,
      .foundID = foundID, .foundIdx = foundIdx };
}

172
static size_t
173
174
175
176
177
178
cdfGridInqHalfSize(int gridID)
{
  return gridInqSize(gridID)/2;
}

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

183
  const size_t dimlen = gridInqSize(gridID)/2;
184

185
186
187
188
189
  struct idSearch search
    = cdfSearchIDBySize(0, (size_t)gridindex, ncgrid, CDF_DIMID_Y,
                        gridRefType, (int)dimlen, gridInqType, cdfGridInqHalfSize);
  int dimID = search.foundID;
  const int iz = search.numNonMatching;
190

191
  if ( dimID == CDI_UNDEFID )
192
193
194
195
196
197
198
199
200
201
202
203
204
    {
      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;
    }

205
  ncgrid[gridindex].gridID = gridID;
206
  ncgrid[gridindex].ncIDs[CDF_DIMID_Y] = dimID;
207
208
209
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
void cdfDefSP(stream_t *streamptr, int gridID, int gridindex)
211
{
212
  // char longname[] = "Spherical harmonic coefficient";
213
  char axisname[5] = "nspX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_SPECTRAL);
215
216
217
218
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
void cdfDefFC(stream_t *streamptr, int gridID, int gridindex)
220
221
{
  char axisname[5] = "nfcX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_FOURIER);
223
224
225
}

static const struct cdfDefGridAxisInqs {
226
  size_t (*axisSize)(int gridID);
227
  int (*axisDimname)(int cdiID, int key, int size, char *mesg);
228
229
  int (*axisName)(int cdiID, int key, int size, char *mesg);
  int (*axisLongname)(int cdiID, int key, int size, char *mesg);
230
231
  int (*axisUnits)(int cdiID, int key, int size, char *mesg);
  void (*axisStdname)(int cdiID, char *dimstdname);
232
  double (*axisVal)(int gridID, size_t index);
233
234
235
236
  const double *(*axisValsPtr)(int gridID);
  const double *(*axisBoundsPtr)(int gridID);
} gridInqsX = {
  .axisSize = gridInqXsize,
237
  .axisDimname = cdiGridInqKeyStr,
238
  .axisName = cdiGridInqKeyStr,
239
  .axisLongname = cdiGridInqKeyStr,
240
241
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqXstdname,
242
243
244
245
246
  .axisVal = gridInqXval,
  .axisValsPtr = gridInqXvalsPtr,
  .axisBoundsPtr = gridInqXboundsPtr,
}, gridInqsY = {
  .axisSize = gridInqYsize,
247
  .axisDimname = cdiGridInqKeyStr,
248
  .axisName = cdiGridInqKeyStr,
249
  .axisLongname = cdiGridInqKeyStr,
250
251
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqYstdname,
252
253
254
255
  .axisVal = gridInqYval,
  .axisValsPtr = gridInqYvalsPtr,
  .axisBoundsPtr = gridInqYboundsPtr,
}, gridInqsZ = {
256
  .axisLongname = cdiZaxisInqKeyStr,
257
258
  .axisUnits = cdiZaxisInqKeyStr,
  .axisStdname = zaxisInqStdname,
259
260
};

261
262
static
void cdfPutGridStdAtts(int fileID, int ncvarid, int gridID, int dimtype, const struct cdfDefGridAxisInqs *inqs)
263
264
{
  size_t len;
265
266
267
268
269
270
271
272
273

  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);
274
  if ( longname[0] && (len = strlen(longname)) )
275
276
    cdf_put_att_text(fileID, ncvarid, "long_name", len, longname);

277
278
279
280
  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)) )
281
    cdf_put_att_text(fileID, ncvarid, "units", len, units);
282
283
284
}

static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285
cdfDefTrajLatLon(stream_t *streamptr, int gridID, int gridindex,
286
                 const struct cdfDefGridAxisInqs *inqs, int dimtype)
287
{
288
  const nc_type xtype = (gridInqDatatype(gridID) == CDI_DATATYPE_FLT32) ? NC_FLOAT : NC_DOUBLE;
289
  ncgrid_t *ncgrid = streamptr->ncgrid;
290

291
  const size_t dimlen = inqs->axisSize(gridID);
292
  if ( dimlen != 1 )
293
    Error("%c size isn't 1 for %s grid!", dimtype, gridNamePtr(gridInqType(gridID)));
294

295
  int ncvarid = ncgrid[gridindex].ncIDs[dimtype == 'X' ? CDF_DIMID_X : CDF_DIMID_Y];
296

297
  if ( ncvarid == CDI_UNDEFID )
298
299
    {
      int dimNcID = streamptr->basetime.ncvarid;
300
      const int fileID  = streamptr->fileID;
301
302
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

303
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
304
      const int keyname = (dimtype == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
305
      inqs->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
306
      cdf_def_var(fileID, axisname, xtype, 1, &dimNcID, &ncvarid);
307
      cdfPutGridStdAtts(fileID, ncvarid, gridID, dimtype, inqs);
308
309
310
311
      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

312
  ncgrid[gridindex].gridID = gridID;
Thomas Jahns's avatar
Thomas Jahns committed
313
314
  /* var ID for trajectory !!! */
  ncgrid[gridindex].ncIDs[dimtype == 'X' ? CDF_DIMID_X : CDF_DIMID_Y] = ncvarid;
315
316
317
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318
void cdfDefTrajLon(stream_t *streamptr, int gridID, int gridindex)
319
{
320
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsX, 'X');
321
322
323
324
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
void cdfDefTrajLat(stream_t *streamptr, int gridID, int gridindex)
326
{
327
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsY, 'Y');
328
329
330
331
332
333
334
}

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;
335
  int dimid = CDI_UNDEFID;
336
337
  char name[CDI_MAX_NAME];

338
  const size_t len = strlen(dimname);
339
340
341
342
343
344
  memcpy(name, dimname, len + 1);

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

345
346
347
      int dimid0;
      const int status = nc_inq_dimid(fileID, name, &dimid0);
      if ( status != NC_NOERR ) break;
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
      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
366
void checkGridName(char *axisname, int fileID)
367
368
369
370
371
372
373
{
  int ncdimid;
  char axisname2[CDI_MAX_NAME];

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

374
  const size_t axisnameLen = strlen(axisname);
375
  memcpy(axisname2, axisname, axisnameLen + 1);
376

377
378
379
380
  do
    {
      if ( iz ) sprintf(axisname2 + axisnameLen, "_%u", iz+1);

381
      if ( nc_inq_varid(fileID, axisname2, &ncdimid) != NC_NOERR ) break;
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397

      ++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;

398
  const size_t axisnameLen = strlen(axisname);
399
400
401
402
403
  memcpy(axisname2, axisname, axisnameLen + 1);
  do
    {
      if ( iz ) sprintf(axisname2 + axisnameLen, "_%u", iz+1);

404
405
      int ncdimid;
      const int status = nc_inq_varid(fileID, axisname2, &ncdimid);
406
407
408
409
410
411
412
      if ( status != NC_NOERR )
        {
          if ( iz )
            {
              /* check that the name does not exist for other zaxes */
              for ( int index = 0; index < nzaxis; index++ )
                {
413
                  const int zaxisID0 = vlistZaxis(vlistID, index);
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
                  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
435
cdfDefAxisCommon(stream_t *streamptr, int gridID, int gridindex, int ndims,
436
437
                 const struct cdfDefGridAxisInqs *gridAxisInq, int dimKey, char axisLetter,
                 void (*finishCyclicBounds)(double *pbounds, size_t dimlen, const double *pvals))
438
{
439
  int dimID = CDI_UNDEFID;
440
441
442
  const int fileID  = streamptr->fileID;
  const size_t dimlen = gridAxisInq->axisSize(gridID);
  const nc_type xtype = (nc_type)cdfDefDatatype(gridInqDatatype(gridID), streamptr);
443

444
  ncgrid_t *ncgrid = streamptr->ncgrid;
445

446
447
448
449
  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
450
  for ( int index = 0; index < gridindex; ++index )
451
    {
452
      const int gridID0 = ncgrid[index].gridID;
453
      assert(gridID0 != CDI_UNDEFID);
454
      const int gridtype0 = gridInqType(gridID0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
456
457
458
      if ( gridtype0 == GRID_GAUSSIAN    ||
           gridtype0 == GRID_LONLAT      ||
           gridtype0 == GRID_PROJECTION  ||
           gridtype0 == GRID_GENERIC )
459
        {
460
          const size_t dimlen0 = gridAxisInq->axisSize(gridID0);
461
462
          char dimname0[CDI_MAX_NAME]; dimname0[0] = 0;
          if ( dimname[0] ) cdiGridInqKeyStr(gridID0, dimKey, CDI_MAX_NAME, dimname0);
463
          const bool lname = dimname0[0] ? strcmp(dimname, dimname0) == 0 : true;
464
          if ( dimlen == dimlen0 && lname )
465
            {
466
              double (*inqVal)(int gridID, size_t index) = gridAxisInq->axisVal;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
              if ( IS_EQUAL(inqVal(gridID0, 0), inqVal(gridID, 0)) &&
468
                   IS_EQUAL(inqVal(gridID0, dimlen-1), inqVal(gridID, dimlen-1)) )
469
                {
470
                  dimID = ncgrid[index].ncIDs[dimKey == CDI_KEY_XDIMNAME ? CDF_DIMID_X : CDF_DIMID_Y];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
                  break;
472
473
474
475
476
                }
            }
        }
    }

477
  if ( dimID == CDI_UNDEFID )
478
    {
479
      int ncvarid = CDI_UNDEFID, ncbvarid = CDI_UNDEFID;
480
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
481
      const int keyname = (axisLetter == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
482
      gridAxisInq->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
483
484
      if ( axisname[0] == 0 ) Error("axis name undefined!");

485
      checkGridName(axisname, fileID);
486
      const size_t axisnameLen = strlen(axisname);
487
488
489
490
491

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

      if ( ndims )
        {
492
          if ( dimname[0] == 0 ) strcpy(dimname, axisname);
493
494
          dimID = checkDimName(fileID, dimlen, dimname);

495
          if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);
496
497
498
        }

      bool gen_bounds = false;
499
      const bool grid_is_cyclic = gridIsCircular(gridID) > 0;
500
501
502
      double *pbounds = NULL;
      if ( pvals )
        {
503
          cdf_def_var(fileID, axisname, xtype, ndims, &dimID, &ncvarid);
504

505
          cdfPutGridStdAtts(fileID, ncvarid, gridID, axisLetter, gridAxisInq);
506
507
508
509
510
          {
            char axisStr[2] = { axisLetter, '\0' };
            cdf_put_att_text(fileID, ncvarid, "axis", 1, axisStr);
          }

511
          size_t nvertex = gridInqNvertex(gridID);
512
513
514
515
516
          pbounds = (double *)gridAxisInq->axisBoundsPtr(gridID);

          if ( CDI_cmor_mode && grid_is_cyclic && !pbounds )
            {
              gen_bounds = true;
517
              nvertex = 2;
518
519
520
521
522
523
524
525
              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);
            }
526
527

          int nvdimID = CDI_UNDEFID;
528
529
530
531
532
          if ( pbounds )
            {
              if ( nc_inq_dimid(fileID, bndsName, &nvdimID) != NC_NOERR )
                cdf_def_dim(fileID, bndsName, nvertex, &nvdimID);
            }
533
          if ( pbounds && nvdimID != CDI_UNDEFID )
534
            {
535
536
537
538
              char boundsname[CDI_MAX_NAME];
              memcpy(boundsname, axisname, axisnameLen);
              boundsname[axisnameLen] = '_';
              memcpy(boundsname + axisnameLen + 1, bndsName, sizeof(bndsName));
539
540
              int dimIDs[2] = { dimID, nvdimID };
              cdf_def_var(fileID, boundsname, xtype, 2, dimIDs, &ncbvarid);
541
              cdf_put_att_text(fileID, ncvarid, "bounds", axisnameLen + sizeof(bndsName), boundsname);
542
543
544
545
546
547
            }
        }

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

548
549
      if ( ncvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncvarid, pvals);
      if ( ncbvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbvarid, pbounds);
550
551
      if ( gen_bounds ) Free(pbounds);

552
      if ( ndims == 0 )
553
        ncgrid[gridindex].ncIDs[dimKey == CDI_KEY_XDIMNAME ? CDF_VARID_X : CDF_VARID_Y] = ncvarid;
554
555
    }

556
  ncgrid[gridindex].gridID = gridID;
557
  ncgrid[gridindex].ncIDs[dimKey == CDI_KEY_XDIMNAME ? CDF_DIMID_X : CDF_DIMID_Y] = dimID;
558
559
}

560
561
static
void finishCyclicXBounds(double *pbounds, size_t dimlen, const double *pvals)
562
563
564
565
566
567
{
  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
568
void cdfDefXaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
569
{
570
571
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsX,
                   CDI_KEY_XDIMNAME, 'X', finishCyclicXBounds);
572
573
}

574
575
static
void finishCyclicYBounds(double *pbounds, size_t dimlen, const double *pvals)
576
577
578
579
580
581
{
  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
582
void cdfDefYaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
583
{
584
585
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsY,
                   CDI_KEY_YDIMNAME, 'Y', finishCyclicYBounds);
586
587
588
}

static
589
void cdfGridCompress(int fileID, int ncvarid, size_t gridsize, int filetype, int comptype)
590
591
{
#if  defined  (HAVE_NETCDF4)
592
  if ( gridsize > 1 && comptype == CDI_COMPRESS_ZIP && (filetype == CDI_FILETYPE_NC4 || filetype == CDI_FILETYPE_NC4C) )
593
    {
594
      cdf_def_var_chunking(fileID, ncvarid, NC_CHUNKED, NULL);
595
596
597
598
599
600
      cdfDefVarDeflate(fileID, ncvarid, 1);
    }
#endif
}

static
601
602
void cdfDefGridReference(stream_t *streamptr, int gridID)
{
603
  const int fileID = streamptr->fileID;
604

605
606
  const int number = gridInqNumber(gridID);
  if ( number > 0 ) cdf_put_att_int(fileID, NC_GLOBAL, "number_of_grid_used", NC_INT, 1, &number);
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643

  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)
644
{
645
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqDatatype(gridID), streamptr);
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
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
  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);
708
      cdfGridCompress(fileID, ncxvarid, xdimlen*ydimlen, streamptr->filetype, streamptr->comptype);
709
710
711
712
713
714
715
716
717

      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 )
        {
718
          const size_t xaxisnameLen = strlen(xaxisname);
719
720
721
          xaxisname[xaxisnameLen] = '_';
          memcpy(xaxisname + xaxisnameLen + 1, bndsName, sizeof (bndsName));
          cdf_def_var(fileID, xaxisname, xtype, ndims, dimIDs, &ncbxvarid);
722
          cdfGridCompress(fileID, ncbxvarid, xdimlen*ydimlen, streamptr->filetype, streamptr->comptype);
723
724
725
726
727
728
729
730
731
732
733
734
735
736

          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);
737
      cdfGridCompress(fileID, ncyvarid, xdimlen*ydimlen, streamptr->filetype, streamptr->comptype);
738
739
740
741
742
743
744
745
746

      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 )
        {
747
          const size_t yaxisnameLen = strlen(yaxisname);
748
749
750
          yaxisname[yaxisnameLen] = '_';
          memcpy(yaxisname + yaxisnameLen + 1, bndsName, sizeof (bndsName));
          cdf_def_var(fileID, yaxisname, xtype, ndims, dimIDs, &ncbyvarid);
751
          cdfGridCompress(fileID, ncbyvarid, xdimlen*ydimlen, streamptr->filetype, streamptr->comptype);
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

          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)
{
790
  ncgrid_t *ncgrid = streamptr->ncgrid;
791

792
793
794
  const size_t dimlen = gridInqSize(gridID);
  const size_t xdimlen = gridInqXsize(gridID);
  const size_t ydimlen = gridInqYsize(gridID);
795

796
  int xdimID = CDI_UNDEFID, ydimID = CDI_UNDEFID;
797
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
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

  size_t ofs = 0;
  do {
    struct idSearch search
      = cdfSearchIDBySize(ofs, (size_t)gridindex, ncgrid, CDF_DIMID_X,
                          GRID_CURVILINEAR, (int)dimlen, gridInqType, gridInqSize);
    const size_t index = search.foundIdx;
    if ( index != SIZE_MAX )
      {
        const int gridID0 = ncgrid[index].gridID;
        if (    IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0))
                && IS_EQUAL(gridInqXval(gridID0, dimlen-1),
                            gridInqXval(gridID, dimlen-1))
                && IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0))
                && IS_EQUAL(gridInqYval(gridID0, dimlen-1),
                            gridInqYval(gridID, dimlen-1)) )
          {
            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;
          }
        ofs = search.foundIdx;
        if ( ofs < (size_t)gridindex )
          continue;
      }
  } while (false);
826

827
  if ( xdimID == CDI_UNDEFID || ydimID == CDI_UNDEFID )
828
    {
829
830
831
832
833
834
835
836
837
838
      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;
839
840
    }

841
  ncgrid[gridindex].gridID = gridID;
842
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = xdimID;
843
  ncgrid[gridindex].ncIDs[CDF_DIMID_Y] = ydimID;
844
  ncgrid[gridindex].ncIDs[CDF_VARID_X] = ncxvarid;
845
  ncgrid[gridindex].ncIDs[CDF_VARID_Y] = ncyvarid;
846
  ncgrid[gridindex].ncIDs[CDF_VARID_A] = ncavarid;
847
848
849
850
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
851
void cdfDefUnstructured(stream_t *streamptr, int gridID, int gridindex)
852
{
853
  ncgrid_t *ncgrid = streamptr->ncgrid;
854

855
  const size_t dimlen = gridInqSize(gridID);
856

857
  int dimID = CDI_UNDEFID;
858
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887

  size_t ofs = 0;
  do {
    struct idSearch search
      = cdfSearchIDBySize(ofs, (size_t)gridindex, ncgrid, CDF_DIMID_X,
                          GRID_UNSTRUCTURED, (int)dimlen, gridInqType, gridInqSize);
    const size_t index = search.foundIdx;
    if ( index != SIZE_MAX )
      {
        const int gridID0 = ncgrid[index].gridID;
        if ( gridInqNvertex(gridID0) == gridInqNvertex(gridID) &&
             IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0)) &&
             IS_EQUAL(gridInqXval(gridID0, dimlen-1),
                      gridInqXval(gridID, dimlen-1)) &&
             IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0)) &&
             IS_EQUAL(gridInqYval(gridID0, dimlen-1),
                      gridInqYval(gridID, dimlen-1)) )
          {
            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;
          }
        ofs = search.foundIdx;
        if ( ofs < (size_t)gridindex )
          continue;
      }
  } while (false);
888

889
  if ( dimID == CDI_UNDEFID )
890
    {
891
      const size_t nvertex = (size_t)gridInqNvertex(gridID);
892
893
894
895
896
897
898
899
      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;
900
901
    }

902
  ncgrid[gridindex].gridID = gridID;
903
  ncgrid[gridindex].ncIDs[CDF_DIMID_X] = dimID;
904
  ncgrid[gridindex].ncIDs[CDF_VARID_X] = ncxvarid;
905
  ncgrid[gridindex].ncIDs[CDF_VARID_Y] = ncyvarid;
906
  ncgrid[gridindex].ncIDs[CDF_VARID_A] = ncavarid;
907
908
909
910
911
912
913
914
915
916
917
}

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

static
void cdf_def_vct_echam(stream_t *streamptr, int zaxisID)
{
918
  const int type = zaxisInqType(zaxisID);
919
920
  if ( type == ZAXIS_HYBRID || type == ZAXIS_HYBRID_HALF )
    {
921
      const int ilev = zaxisInqVctSize(zaxisID)/2;
922
923
      if ( ilev == 0 ) return;

924
      const int mlev = ilev - 1;
925
926
927
928

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

933
      const int fileID = streamptr->fileID;
934
935
936

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

Thomas Jahns's avatar
Thomas Jahns committed
937
      int ncdimid = -1, ncdimid2 = -1;
938
939
      int hyaiid, hybiid, hyamid = -1, hybmid = -1;

940
      cdf_def_dim(fileID, "nhyi", (size_t)ilev, &ncdimid2);
941
942
943
944
945
946
947
948
      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);
        }
949
950

      streamptr->vct.ilev   = ilev;
951
      streamptr->vct.mlev   = mlev;
952
953
954
955
956
957
      streamptr->vct.mlevID = ncdimid;
      streamptr->vct.ilevID = ncdimid2;

      {
        static const char lname_n[] = "long_name",
          units_n[] = "units",
958
          lname_v_ai[] = "hybrid A coefficient at layer interfaces",
959
960
          units_v_ai[] = "Pa",
          lname_v_bi[] = "hybrid B coefficient at layer interfaces",
961
          units_v_bi[] = "1";
962
963
964
965
966
967
        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 },
968
969
        };
        enum { tabLen = sizeof (tab) / sizeof (tab[0]) };
970
        const int ids[tabLen] = { hyaiid, hyaiid, hybiid, hybiid };
971
972
973
974
975
976
977
978
979
980
981
982
983
        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[]
          = {
984
985
986
987
988
989
          { 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]) };
990
        const int ids[tabLen] = { hyamid, hyamid, hybmid, hybmid };
991
992
993
994
995
996
997
998
999
1000
1001
1002
        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);

1003
1004
1005
      size_t start;
      size_t count = 1;
      double mval;
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
      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
1018
void cdf_def_vct_cf(stream_t *streamptr, int zaxisID, int nclevID, int ncbndsID, int p0status, double p0value)
1019
{
1020
  const int type = zaxisInqType(zaxisID);
1021
1022
  if ( type == ZAXIS_HYBRID || type == ZAXIS_HYBRID_HALF )
    {
1023
      const int ilev = zaxisInqVctSize(zaxisID)/2;
1024
1025
      if ( ilev == 0 ) return;

1026
      const int mlev = ilev - 1;
1027
1028
1029
1030
1031
1032
1033
1034
1035
      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;
        }

1036
      const int fileID = streamptr->fileID;
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048

      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;

1049
1050
1051
1052
      if ( p0status == 0 )
        cdf_def_var(fileID, "a", NC_DOUBLE, 1, dimIDs,  &hyamid);
      else
        cdf_def_var(fileID, "ap", NC_DOUBLE, 1, dimIDs,  &hyamid);
1053
1054
1055
1056
1057
      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);
1058

1059
1060
1061
1062
1063
1064
        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);
1065

1066
1067
1068
1069
1070
1071
        static const char units[] = "1";
        cdf_put_att_text(fileID, hybmid, "units", sizeof (units) - 1, units);
      }

      if ( ncbndsID != -1 )
        {
1072
1073
1074
1075
          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);
1076
          cdf_def_var(fileID, "b_bnds",  NC_DOUBLE, 2, dimIDs, &hybiid);
1077

1078
1079
1080
          {
            static const char lname[] = "vertical coordinate formula term: ap(k+1/2)";
            cdf_put_att_text(fileID, hyaiid, "long_name", sizeof (lname) - 1, lname);
1081

1082
1083
1084
1085
1086
1087
            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);
1088

1089
1090
1091
1092
1093
1094
1095
1096
            static const char units[] = "1";
            cdf_put_att_text(fileID, hybiid, "units", sizeof (units) - 1, units);
          }
        }

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

1097
      const int vctsize = zaxisInqVctSize(zaxisID);
1098
      double *vct = (double*) malloc(vctsize*sizeof(double));;
1099
1100
1101
1102
1103
      zaxisInqVct(zaxisID, vct);

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

1104
      double *tarray = (double*) malloc(ilev*2*sizeof(double));
1105
1106
1107
1108
1109

      if ( ncbndsID != -1 )
        {
          for ( int i = 0; i < mlev; ++i )
            {
1110
1111
              tarray[2*i  ] = vct[i];
              tarray[2*i+1] = vct[i+1];
1112
1113
1114
1115
1116
            }
          cdf_put_var_double(fileID, hyaiid, tarray);

          for ( int i = 0; i < mlev; ++i )
            {
1117
1118
              tarray[2*i  ] = vct[ilev+i];
              tarray[2*i+1] = vct[ilev+i+1];
1119
1120
1121
1122
1123
            }
          cdf_put_var_double(fileID, hybiid, tarray);
        }

      for ( int i = 0; i < mlev; ++i )
1124
        tarray[i] = (vct[i] + vct[i+1]) * 0.5;
1125
1126
1127
      cdf_put_var_double(fileID, hyamid, tarray);

      for ( int i = 0; i < mlev; ++i )
1128
        tarray[i] = (vct[ilev+i] + vct[ilev+i+1]) * 0.5;
1129
      cdf_put_var_double(fileID, hybmid, tarray);
1130
1131
1132

      free(tarray);
      free(vct);
1133
1134
1135
1136
1137
1138
    }
}

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

static
1139
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)