stream_cdf_o.c 65.9 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, bool addVarToGrid,
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 || addVarToGrid)
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
568
void finishCyclicYBounds(double *pbounds, size_t dimlen, const double *pvals)
569
{
570
571
  pbounds[0] = copysign(90.0, pvals[0]);
  pbounds[2*dimlen-1] = copysign(90.0, pvals[dimlen-1]);
572
573
}

574
static
575
void cdfDefXaxis(stream_t *streamptr, int gridID, int gridindex, int ndims, bool addVarToGrid)
576
{
577
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, addVarToGrid, &gridInqsX, CDI_KEY_XDIMNAME, 'X', finishCyclicXBounds);
578
579
580
}

static
581
void cdfDefYaxis(stream_t *streamptr, int gridID, int gridindex, int ndims, bool addVarToGrid)
582
{
583
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, addVarToGrid, &gridInqsY, CDI_KEY_YDIMNAME, 'Y', finishCyclicYBounds);
584
585
586
}

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

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

603
604
  const int number = gridInqNumber(gridID);
  if ( number > 0 ) cdf_put_att_int(fileID, NC_GLOBAL, "number_of_grid_used", NC_INT, 1, &number);
605
606
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

  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)
642
{
643
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqDatatype(gridID), streamptr);
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
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
  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);
706
      cdfGridCompress(fileID, ncxvarid, xdimlen*ydimlen, streamptr->filetype, streamptr->comptype);
707
708
709
710
711
712
713
714
715

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

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

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

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

790
791
792
  const size_t dimlen = gridInqSize(gridID);
  const size_t xdimlen = gridInqXsize(gridID);
  const size_t ydimlen = gridInqYsize(gridID);
793

794
  int xdimID = CDI_UNDEFID, ydimID = CDI_UNDEFID;
795
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823

  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);
824

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

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


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

853
  const size_t dimlen = gridInqSize(gridID);
854

855
  int dimID = CDI_UNDEFID;
856
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
857
858
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

  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);
886

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

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

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

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

922
      const int mlev = ilev - 1;
923
924
925
926

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

931
      const int fileID = streamptr->fileID;
932
933
934

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

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

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

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

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

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

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

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

      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;

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

1057
1058
1059
1060
1061
1062
        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);
1063

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

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

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

1080
1081
1082
1083
1084
1085
            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);
1086

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

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

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

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

1102
      double *tarray = (double*) malloc(ilev*2*sizeof(double));
1103
1104
1105
1106
1107

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

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

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

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

      free(tarray);
      free(vct);