cdf_write.c 40 KB
Newer Older
1
2
#ifdef HAVE_CONFIG_H
#include "config.h"
3
4
5
6
7
8
9
10
11
12
13
#endif

#ifdef HAVE_LIBNETCDF

#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
#include "stream_cdf.h"
#include "cdf.h"
#include "cdf_int.h"
#include "vlist.h"
14
#include "vlist_var.h"
15
16
17
18


void cdfDefVarDeflate(int ncid, int ncvarid, int deflate_level)
{
19
#if  defined(HAVE_NETCDF4)
20
21
22
23
24
25
26
  int retval;
  /* Set chunking, shuffle, and deflate. */
  int shuffle = 1;
  int deflate = 1;

  if ( deflate_level < 1 || deflate_level > 9 ) deflate_level = 1;

27
  if ( (retval = nc_def_var_deflate(ncid, ncvarid, shuffle, deflate, deflate_level)) )
28
29
30
31
    {
      Error("nc_def_var_deflate failed, status = %d", retval);
    }
#else
32
  static bool lwarn = true;
33
34
  if ( lwarn )
    {
35
      lwarn = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
      Warning("Deflate compression failed, NetCDF4 not available!");
37
38
39
40
    }
#endif
}

41
int cdfDefDatatype(int datatype, stream_t *streamptr)
42
43
44
{
  int xtype = NC_FLOAT;

45
  if ( streamptr->filetype == CDI_FILETYPE_NC4 )
46
    {
47
48
49
      if      ( datatype == CDI_DATATYPE_INT8   ) xtype = NC_BYTE;
      else if ( datatype == CDI_DATATYPE_INT16  ) xtype = NC_SHORT;
      else if ( datatype == CDI_DATATYPE_INT32  ) xtype = NC_INT;
50
#ifdef  HAVE_NETCDF4
51
52
53
      else if ( datatype == CDI_DATATYPE_UINT8  ) xtype = NC_UBYTE;
      else if ( datatype == CDI_DATATYPE_UINT16 ) xtype = NC_USHORT;
      else if ( datatype == CDI_DATATYPE_UINT32 ) xtype = NC_UINT;
54
      else if ( datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64 )
55
        Error("CDI library does not support complex numbers with NetCDF4/HDF5!");
56
#else
57
58
59
      else if ( datatype == CDI_DATATYPE_UINT8  ) xtype = NC_SHORT;
      else if ( datatype == CDI_DATATYPE_UINT16 ) xtype = NC_INT;
      else if ( datatype == CDI_DATATYPE_UINT32 ) xtype = NC_INT;
60
61
      else if ( datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64 )
        Error("CDI library does not support complex numbers with NetCDF4 classic!");
62
#endif
63
      else if ( datatype == CDI_DATATYPE_FLT64  ) xtype = NC_DOUBLE;
64
      else if ( datatype == CDI_DATATYPE_FLT32  ) xtype = NC_FLOAT;
65
66
67
    }
  else
    {
68
69
70
71
72
73
74
      if      ( datatype == CDI_DATATYPE_INT8   ) xtype = NC_BYTE;
      else if ( datatype == CDI_DATATYPE_INT16  ) xtype = NC_SHORT;
      else if ( datatype == CDI_DATATYPE_INT32  ) xtype = NC_INT;
      else if ( datatype == CDI_DATATYPE_UINT8  ) xtype = NC_SHORT;
      else if ( datatype == CDI_DATATYPE_UINT16 ) xtype = NC_INT;
      else if ( datatype == CDI_DATATYPE_UINT32 ) xtype = NC_INT;
      else if ( datatype == CDI_DATATYPE_FLT64  ) xtype = NC_DOUBLE;
75
76
77
      else if ( datatype == CDI_DATATYPE_FLT32  ) xtype = NC_FLOAT;
      else if ( datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64 )
        Error("CDI library does not support complex numbers with NetCDF classic!");
78
79
80
81
82
83
84
85
    }

  return xtype;
}

static
void cdfDefVarMissval(stream_t *streamptr, int varID, int dtype, int lcheck)
{
86
  if ( streamptr->vars[varID].defmiss == false )
87
88
89
90
91
92
93
94
    {
      int vlistID = streamptr->vlistID;
      int fileID  = streamptr->fileID;
      int ncvarid = streamptr->vars[varID].ncvarid;
      double missval = vlistInqVarMissval(vlistID, varID);

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

95
      int xtype = cdfDefDatatype(dtype, streamptr);
96
97
98

      if ( xtype == NC_BYTE && missval > 127 && missval < 256 ) xtype = NC_INT;

99
100
      if ( lcheck == 0 ||
           streamptr->ncmode != 2 ||
101
           streamptr->filetype == CDI_FILETYPE_NC ||
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
           streamptr->filetype == CDI_FILETYPE_NC2||
           streamptr->filetype == CDI_FILETYPE_NC5 )
104
        cdf_put_att_double(fileID, ncvarid, "_FillValue", (nc_type) xtype, 1, &missval);
105

106
107
108
109
      cdf_put_att_double(fileID, ncvarid, "missing_value", (nc_type) xtype, 1, &missval);

      if ( lcheck && streamptr->ncmode == 2 ) cdf_enddef(fileID);

110
      streamptr->vars[varID].defmiss = true;
111
112
113
114
115
116
117
118
119
120
    }
}

static
void cdfDefInstitut(stream_t *streamptr)
{
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;
  int instID  = vlistInqInstitut(vlistID);

121
  if ( instID != CDI_UNDEFID )
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    {
      const char *longname = institutInqLongnamePtr(instID);
      if ( longname )
	{
	  size_t len = strlen(longname);
	  if ( len > 0 )
	    {
	      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
	      cdf_put_att_text(fileID, NC_GLOBAL, "institution", len, longname);
	      if ( streamptr->ncmode == 2 ) cdf_enddef(fileID);
	    }
	}
    }
}

static
void cdfDefSource(stream_t *streamptr)
{
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;
  int modelID = vlistInqModel(vlistID);

144
  if ( modelID != CDI_UNDEFID )
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    {
      const char *longname = modelInqNamePtr(modelID);
      if ( longname )
	{
          size_t len = strlen(longname);
	  if ( len > 0 )
	    {
	      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
	      cdf_put_att_text(fileID, NC_GLOBAL, "source", len, longname);
	      if ( streamptr->ncmode == 2 ) cdf_enddef(fileID);
	    }
	}
    }
}

static inline
void *resizeBuf(void **buf, size_t *bufSize, size_t reqSize)
{
163
  if ( reqSize > *bufSize )
164
165
166
167
168
169
170
    {
      *buf = Realloc(*buf, reqSize);
      *bufSize = reqSize;
    }
  return *buf;
}

171
172

void cdfDefineAttributes(int cdiID, int varID, int fileID, int ncvarID)
173
174
175
176
177
178
179
180
{
  int atttype, attlen;
  size_t len;
  char attname[CDI_MAX_NAME+1];
  void *attBuf = NULL;
  size_t attBufSize = 0;

  int natts;
181
  cdiInqNatts(cdiID, varID, &natts);
182

183
  for ( int iatt = 0; iatt < natts; ++iatt )
184
    {
185
      cdiInqAtt(cdiID, varID, iatt, attname, &atttype, &attlen);
186
187
188

      if ( attlen == 0 ) continue;

189
      if ( atttype == CDI_DATATYPE_TXT )
190
191
192
        {
          size_t attSize = (size_t)attlen*sizeof(char);
          char *atttxt = (char *)resizeBuf(&attBuf, &attBufSize, attSize);
193
          cdiInqAttTxt(cdiID, varID, attname, attlen, atttxt);
194
195
196
          len = (size_t)attlen;
          cdf_put_att_text(fileID, ncvarID, attname, len, atttxt);
        }
197
198
199
      else if ( atttype == CDI_DATATYPE_INT8  || atttype == CDI_DATATYPE_UINT8  ||
                atttype == CDI_DATATYPE_INT16 || atttype == CDI_DATATYPE_UINT16 ||
                atttype == CDI_DATATYPE_INT32 || atttype == CDI_DATATYPE_UINT32 )
200
201
202
        {
          size_t attSize = (size_t)attlen*sizeof(int);
          int *attint = (int *)resizeBuf(&attBuf, &attBufSize, attSize);
203
          cdiInqAttInt(cdiID, varID, attname, attlen, &attint[0]);
204
          len = (size_t)attlen;
205
206
          nc_type xtype = (atttype == CDI_DATATYPE_INT8)  ? NC_BYTE :
                          (atttype == CDI_DATATYPE_INT16) ? NC_SHORT :
207
#if  defined  (HAVE_NETCDF4)
208
209
210
                          (atttype == CDI_DATATYPE_UINT8)  ? NC_UBYTE :
                          (atttype == CDI_DATATYPE_UINT16) ? NC_USHORT :
                          (atttype == CDI_DATATYPE_UINT32) ? NC_UINT :
211
212
213
#endif
                          NC_INT;
          cdf_put_att_int(fileID, ncvarID, attname, xtype, len, attint);
214
        }
215
      else if ( atttype == CDI_DATATYPE_FLT32 || atttype == CDI_DATATYPE_FLT64 )
216
217
218
        {
          size_t attSize = (size_t)attlen * sizeof(double);
          double *attflt = (double *)resizeBuf(&attBuf, &attBufSize, attSize);
219
          cdiInqAttFlt(cdiID, varID, attname, attlen, attflt);
220
          len = (size_t)attlen;
221
          if ( atttype == CDI_DATATYPE_FLT32 )
222
223
224
225
226
227
228
229
230
            {
              float attflt_sp[len];
              for ( size_t i = 0; i < len; ++i ) attflt_sp[i] = (float)attflt[i];
              cdf_put_att_float(fileID, ncvarID, attname, NC_FLOAT, len, attflt_sp);
            }
          else
            cdf_put_att_double(fileID, ncvarID, attname, NC_DOUBLE, len, attflt);
        }
    }
231

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
  Free(attBuf);
}

static
void cdfDefGlobalAtts(stream_t *streamptr)
{
  if ( streamptr->globalatts ) return;

  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;

  cdfDefSource(streamptr);
  cdfDefInstitut(streamptr);

  int natts;
247
  cdiInqNatts(vlistID, CDI_GLOBAL, &natts);
248
249
250

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

251
  cdfDefineAttributes(vlistID, CDI_GLOBAL, fileID, NC_GLOBAL);
252
253
254
255
256
257
258
259
260
261
262
263
264

  if ( natts > 0 && streamptr->ncmode == 2 ) cdf_enddef(fileID);

  streamptr->globalatts = 1;
}

static
void cdfDefLocalAtts(stream_t *streamptr)
{
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;

  if ( streamptr->localatts ) return;
265
  if ( vlistInqInstitut(vlistID) != CDI_UNDEFID ) return;
266
267
268
269
270
271
272
273

  streamptr->localatts = 1;

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

  for ( int varID = 0; varID < streamptr->nvars; varID++ )
    {
      int instID = vlistInqVarInstitut(vlistID, varID);
274
      if ( instID != CDI_UNDEFID )
275
276
277
278
279
280
281
282
283
284
285
286
287
288
	{
          int ncvarid = streamptr->vars[varID].ncvarid;
  	  const char *name = institutInqNamePtr(instID);
	  if ( name )
	    {
              size_t len = strlen(name);
	      cdf_put_att_text(fileID, ncvarid, "institution", len, name);
	    }
	}
      }

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

289
290
291
static
void cdf_get_gmapvarname(int gridID, char *gmapvarname)
{
292
  int pgridID = gridID;
293
  char mapping[CDI_MAX_NAME]; mapping[0] = 0;
294
  cdiGridInqKeyStr(pgridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapping);
295
296
297
298

  if ( !mapping[0] )
    {
      int projID = gridInqProj(gridID);
299
300
301
      if ( projID != CDI_UNDEFID )
        {
          pgridID = projID;
302
          cdiGridInqKeyStr(pgridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapping);
303
        }
304
305
306
    }

  if ( mapping[0] )
307
    cdiGridInqKeyStr(pgridID, CDI_KEY_MAPPING, CDI_MAX_NAME, gmapvarname);
308
309
}

310
311
312
313
314
315
316
317
318
319
320
321
322
323
static
int nc_grid_index(stream_t *streamptr, int gridID)
{
  int index = 0;
  int vlistID = streamptr->vlistID;
  int ngrids = vlistNgrids(vlistID);
  for ( index = 0; index < ngrids; ++index )
    if ( streamptr->ncgrid[index].gridID == gridID ) break;

  assert(index < ngrids);

  return index;
}

324
325
326
327
static
int cdfDefVar(stream_t *streamptr, int varID)
{
  int ncvarid = -1;
328
  int xid = CDI_UNDEFID, yid = CDI_UNDEFID;
329
330
331
332
333
334
335
336
337
338
339
340
341
342
  size_t xsize = 0, ysize = 0;
  int dims[4];
  size_t chunks[4] = {0,0,0,0};
  int ndims = 0;
  int tablenum;
  int dimorder[3];
  size_t iax = 0;
  char axis[5];

  int fileID  = streamptr->fileID;

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

343
  if ( streamptr->vars[varID].ncvarid != CDI_UNDEFID )
344
345
    return streamptr->vars[varID].ncvarid;

346
347
348
349
350
351
  int vlistID  = streamptr->vlistID;
  int gridID   = vlistInqVarGrid(vlistID, varID);
  int zaxisID  = vlistInqVarZaxis(vlistID, varID);
  int timetype = vlistInqVarTimetype(vlistID, varID);
  int code     = vlistInqVarCode(vlistID, varID);
  int param    = vlistInqVarParam(vlistID, varID);
352
353
354
355
356
357
358
  int pnum, pcat, pdis;
  cdiDecodeParam(param, &pnum, &pcat, &pdis);

  int chunktype = vlistInqVarChunkType(vlistID, varID);

  vlistInqVarDimorder(vlistID, varID, &dimorder);

359
  size_t gridsize = gridInqSize(gridID);
360
  bool lchunk = (gridsize >= 16);
361
  int gridtype  = gridInqType(gridID);
362
  int gridindex = nc_grid_index(streamptr, gridID);
363
364
  if ( gridtype != GRID_TRAJECTORY )
    {
365
      xid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
366
      yid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_Y];
367
368
      if ( xid != CDI_UNDEFID ) cdf_inq_dimlen(fileID, xid, &xsize);
      if ( yid != CDI_UNDEFID ) cdf_inq_dimlen(fileID, yid, &ysize);
369
370
371
372
    }

  int zaxisindex = vlistZaxisIndex(vlistID, zaxisID);
  int zid = streamptr->zaxisID[zaxisindex];
373
  bool zaxis_is_scalar = false;
374
  if ( zid == CDI_UNDEFID ) zaxis_is_scalar = zaxisInqScalar(zaxisID) > 0;
375

376
  if ( dimorder[0] != 3 ) lchunk = false; /* ZYX and ZXY */
377

378
  if ( ((dimorder[0]>0)+(dimorder[1]>0)+(dimorder[2]>0)) < ((xid!=CDI_UNDEFID)+(yid!=CDI_UNDEFID)+(zid!=CDI_UNDEFID)) )
379
380
381
382
383
384
385
    {
      printf("zid=%d  yid=%d  xid=%d\n", zid, yid, xid);
      Error("Internal problem, dimension order missing!");
    }

  int tid = streamptr->basetime.ncdimid;

386
  if ( vlistHasTime(vlistID) && timetype != TIME_CONSTANT )
387
    {
388
      if ( tid == CDI_UNDEFID ) Error("Internal problem, time undefined!");
389
390
391
392
393
      chunks[ndims] = 1;
      dims[ndims++] = tid;
      axis[iax++] = 'T';
    }
  /*
394
395
396
  if ( zid != CDI_UNDEFID ) axis[iax++] = 'Z';
  if ( zid != CDI_UNDEFID ) chunks[ndims] = 1;
  if ( zid != CDI_UNDEFID ) dims[ndims++] = zid;
397

398
399
  if ( yid != CDI_UNDEFID ) chunks[ndims] = ysize;
  if ( yid != CDI_UNDEFID ) dims[ndims++] = yid;
400

401
402
  if ( xid != CDI_UNDEFID ) chunks[ndims] = xsize;
  if ( xid != CDI_UNDEFID ) dims[ndims++] = xid;
403
  */
404
  size_t chunk_size_max = 65536;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
405
406
407
408
409
410
  if ( chunktype != CDI_CHUNK_LINES && gridsize > INT32_MAX )
    {
      if ( CDI_Debug ) fprintf(stderr, "gridsize > %d, changed chunktype to CDI_CHUNK_LINES!\n", INT32_MAX);
      chunktype = CDI_CHUNK_LINES;
    }

411
412
  for ( int id = 0; id < 3; ++id )
    {
413
      if ( dimorder[id] == 3 && zid != CDI_UNDEFID )
414
415
416
417
418
419
        {
          axis[iax++] = 'Z';
          chunks[ndims] = 1;
          dims[ndims] = zid;
          ndims++;
        }
420
      else if ( dimorder[id] == 2 && yid != CDI_UNDEFID )
421
        {
422
423
424
425
          if ( chunktype == CDI_CHUNK_AUTO )
            chunks[ndims] = (chunk_size_max > gridsize) ? ysize : chunk_size_max/xsize;
          else
            chunks[ndims] = (chunktype == CDI_CHUNK_LINES) ? 1 : ysize;
426
427
428
          dims[ndims] = yid;
          ndims++;
        }
429
      else if ( dimorder[id] == 1 && xid != CDI_UNDEFID )
430
        {
431
432
433
          if ( chunktype == CDI_CHUNK_AUTO && yid == CDI_UNDEFID )
            chunks[ndims] = (chunk_size_max > xsize) ? xsize : chunk_size_max;
          else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
            chunks[ndims] = (xsize > INT32_MAX) ? INT32_MAX : xsize;
435
436
437
438
439
440
          dims[ndims] = xid;
          ndims++;
        }
    }

  if ( CDI_Debug )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
    fprintf(stderr, "chunktype %d  chunks %zu %zu %zu %zu\n", chunktype, chunks[0], chunks[1], chunks[2], chunks[3]);
442

443
444
445
446
447
448
449
450
451
  char varname[CDI_MAX_NAME];
  char name[CDI_MAX_NAME]; name[0] = 0;
  char longname[CDI_MAX_NAME]; longname[0] = 0;
  char stdname[CDI_MAX_NAME]; stdname[0] = 0;
  char units[CDI_MAX_NAME]; units[0] = 0;
  if ( vlistInqVarNamePtr(vlistID, varID) ) vlistInqVarName(vlistID, varID, name);
  vlistInqVarLongname(vlistID, varID, longname);
  vlistInqVarStdname(vlistID, varID, stdname);
  vlistInqVarUnits(vlistID, varID, units);
452

453
  int tableID  = vlistInqVarTable(vlistID, varID);
454
  if ( !name[0] ) tableInqEntry(tableID, code, -1, name, longname, units);
455
  if ( name[0] )
456
457
458
    {
      sprintf(varname, "%s", name);

459
460
      bool checkname = true;
      int iz = 0;
461
462
463
464
465

      while ( checkname )
        {
          if ( iz ) sprintf(varname, "%s_%d", name, iz+1);

466
467
          int status = nc_inq_varid(fileID, varname, &ncvarid);
          if ( status != NC_NOERR ) checkname = false;
468
469
470
471
472
473
474
475
476
477
478
479
480
481

          if ( checkname ) iz++;

          if ( iz >= CDI_MAX_NAME ) Error("Double entry of variable name '%s'!", name);
        }

      if ( strcmp(name, varname) != 0 )
        {
          if ( iz == 1 )
            Warning("Changed double entry of variable name '%s' to '%s'!", name, varname);
          else
            Warning("Changed multiple entry of variable name '%s' to '%s'!", name, varname);
        }

482
      strcpy(name, varname);
483
484
485
486
487
488
489
490
491
492
493
494
495
    }
  else
    {
      if ( code < 0 ) code = -code;
      if ( pnum < 0 ) pnum = -pnum;

      if ( pdis == 255 )
	sprintf(varname, "var%d", code);
      else
	sprintf(varname, "param%d.%d.%d", pnum, pcat, pdis);

      char *varname2 = varname+strlen(varname);

496
      bool checkname = true;
497
498
499
500
501
502
503
      int iz = 0;

      while ( checkname )
        {
          if ( iz ) sprintf(varname2, "_%d", iz+1);

          int status = nc_inq_varid(fileID, varname, &ncvarid);
504
          if ( status != NC_NOERR ) checkname = false;
505
506
507
508
509
510

          if ( checkname ) iz++;

          if ( iz >= CDI_MAX_NAME ) break;
        }

511
      strcpy(name, varname);
512
513
514
515
516
517
518
      code = 0;
      pdis = 255;
    }

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

  int dtype = vlistInqVarDatatype(vlistID, varID);
519
  int xtype = cdfDefDatatype(dtype, streamptr);
520
521
522
523

  cdf_def_var(fileID, name, (nc_type) xtype, ndims, dims, &ncvarid);

#if  defined  (HAVE_NETCDF4)
524
  if ( lchunk && (streamptr->filetype == CDI_FILETYPE_NC4 || streamptr->filetype == CDI_FILETYPE_NC4C) )
525
    cdf_def_var_chunking(fileID, ncvarid, NC_CHUNKED, chunks);
526
527
#endif

528
  if ( streamptr->comptype == CDI_COMPRESS_ZIP )
529
    {
530
      if ( lchunk && (streamptr->filetype == CDI_FILETYPE_NC4 || streamptr->filetype == CDI_FILETYPE_NC4C) )
531
532
533
534
535
536
537
        {
          cdfDefVarDeflate(fileID, ncvarid, streamptr->complevel);
        }
      else
        {
          if ( lchunk )
            {
538
              static bool lwarn = true;
539
540
              if ( lwarn )
                {
541
                  lwarn = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
                  Warning("Deflate compression is only available for NetCDF4!");
543
544
545
546
547
                }
            }
        }
    }

548
549
550
  if ( *stdname )  cdf_put_att_text(fileID, ncvarid, "standard_name", strlen(stdname), stdname);
  if ( *longname ) cdf_put_att_text(fileID, ncvarid, "long_name", strlen(longname), longname);
  if ( *units )    cdf_put_att_text(fileID, ncvarid, "units", strlen(units), units);
551
552
553
554
555
556
557
558
559
560
561

  if ( code > 0 && pdis == 255 )
    cdf_put_att_int(fileID, ncvarid, "code", NC_INT, 1, &code);

  if ( pdis != 255 )
    {
      char paramstr[32];
      cdiParamToString(param, paramstr, sizeof(paramstr));
      cdf_put_att_text(fileID, ncvarid, "param", strlen(paramstr), paramstr);
    }

562
  if ( tableID != CDI_UNDEFID )
563
564
565
566
567
568
    {
      tablenum = tableInqNum(tableID);
      if ( tablenum > 0 )
        cdf_put_att_int(fileID, ncvarid, "table", NC_INT, 1, &tablenum);
    }

569
  char coordinates[CDI_MAX_NAME]; coordinates[0] = 0;
570

571
  if ( zaxis_is_scalar || zaxisInqType(zaxisID) == ZAXIS_CHAR )
572
573
574
575
576
577
578
579
580
581
    {
      int nczvarID = streamptr->nczvarID[zaxisindex];
      if ( nczvarID != CDI_UNDEFID )
        {
          size_t len = strlen(coordinates);
          if ( len ) coordinates[len++] = ' ';
          cdf_inq_varname(fileID, nczvarID, coordinates+len);
        }
    }

582
  if ( gridtype != GRID_GENERIC && gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN &&
583
       gridtype != GRID_PROJECTION && gridtype != GRID_CURVILINEAR && gridtype != GRID_CHARXY )
584
585
586
    {
      size_t len = strlen(gridNamePtr(gridtype));
      if ( len > 0 )
587
        cdf_put_att_text(fileID, ncvarid, "CDI_grid_type", len, gridNamePtr(gridtype));
588
589
    }

590
  char gmapvarname[CDI_MAX_NAME]; gmapvarname[0] = 0;
591

592
  cdf_get_gmapvarname(gridID, gmapvarname);
593
594
595
596

  if ( gmapvarname[0] ) cdf_put_att_text(fileID, ncvarid, "grid_mapping", strlen(gmapvarname), gmapvarname);

  if ( gridtype == GRID_TRAJECTORY )
597
    {
598
      cdf_put_att_text(fileID, ncvarid, "coordinates", 9, "tlon tlat" );
599
    }
600
  else if ( gridtype == GRID_LONLAT && xid == CDI_UNDEFID && yid == CDI_UNDEFID && gridsize == 1 )
601
    {
602
      int ncxvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_X];
603
      int ncyvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_Y];
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
      if ( ncyvarID != CDI_UNDEFID )
        {
          size_t len = strlen(coordinates);
          if ( len ) coordinates[len++] = ' ';
          cdf_inq_varname(fileID, ncyvarID, coordinates+len);
        }
      if ( ncxvarID != CDI_UNDEFID )
        {
          size_t len = strlen(coordinates);
          if ( len ) coordinates[len++] = ' ';
          cdf_inq_varname(fileID, ncxvarID, coordinates+len);
        }
    }
  else if ( gridtype == GRID_UNSTRUCTURED || gridtype == GRID_CURVILINEAR )
    {
      char cellarea[CDI_MAX_NAME] = "area: ";
620
      int ncxvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_X];
621
      int ncyvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_Y];
622
      int ncavarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_A];
623
624
      // CMOR order: coordinates = "lat lon"
      if ( cdiCoordinatesLonLat )
625
        {
626
627
628
629
630
631
632
633
634
635
636
637
          if ( ncxvarID != CDI_UNDEFID )
            {
              size_t len = strlen(coordinates);
              if ( len ) coordinates[len++] = ' ';
              cdf_inq_varname(fileID, ncxvarID, coordinates+len);
            }
          if ( ncyvarID != CDI_UNDEFID )
            {
              size_t len = strlen(coordinates);
              if ( len ) coordinates[len++] = ' ';
              cdf_inq_varname(fileID, ncyvarID, coordinates+len);
            }
638
        }
639
      else
640
        {
641
642
643
644
645
646
647
648
649
650
651
652
          if ( ncyvarID != CDI_UNDEFID )
            {
              size_t len = strlen(coordinates);
              if ( len ) coordinates[len++] = ' ';
              cdf_inq_varname(fileID, ncyvarID, coordinates+len);
            }
          if ( ncxvarID != CDI_UNDEFID )
            {
              size_t len = strlen(coordinates);
              if ( len ) coordinates[len++] = ' ';
              cdf_inq_varname(fileID, ncxvarID, coordinates+len);
            }
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
        }

      if ( ncavarID != CDI_UNDEFID )
        {
          size_t len = strlen(cellarea);
          cdf_inq_varname(fileID, ncavarID, cellarea+len);
          len = strlen(cellarea);
          cdf_put_att_text(fileID, ncvarid, "cell_measures", len, cellarea);
        }

      if ( gridtype == GRID_UNSTRUCTURED )
        {
          int position = gridInqPosition(gridID);
          if ( position > 0 )
            cdf_put_att_int(fileID, ncvarid, "number_of_grid_in_reference", NC_INT, 1, &position);
        }
    }
  else if ( gridtype == GRID_SPECTRAL || gridtype == GRID_FOURIER )
    {
      int gridTruncation = gridInqTrunc(gridID);
      axis[iax++] = '-';
      axis[iax++] = '-';
      cdf_put_att_text(fileID, ncvarid, "axis", iax, axis);
      cdf_put_att_int(fileID, ncvarid, "truncation", NC_INT, 1, &gridTruncation);
    }
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
  else if ( gridtype == GRID_CHARXY )
    {
      if ( gridInqXIsc(gridID) )
        {
          int ncxvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_X];
          size_t len = strlen(coordinates);
          if ( len ) coordinates[len++] = ' ';
          cdf_inq_varname(fileID, ncxvarID, coordinates+len);
        }
      else if ( gridInqYIsc(gridID) )
        {
          int ncyvarID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_Y];
          size_t len = strlen(coordinates);
          if ( len ) coordinates[len++] = ' ';
          cdf_inq_varname(fileID, ncyvarID, coordinates+len);
        }
    }
695
696
697
698
699
700
701
702

  size_t len = strlen(coordinates);
  if ( len ) cdf_put_att_text(fileID, ncvarid, "coordinates", len, coordinates);

  /*  if ( xtype == NC_BYTE || xtype == NC_SHORT || xtype == NC_INT ) */
    {
      int astype = NC_DOUBLE;

703
704
705
706
      double addoffset   = vlistInqVarAddoffset(vlistID, varID);
      double scalefactor = vlistInqVarScalefactor(vlistID, varID);
      bool laddoffset   = IS_NOT_EQUAL(addoffset, 0);
      bool lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722

      if ( laddoffset || lscalefactor )
        {
          if ( IS_EQUAL(addoffset,   (double) ((float) addoffset)) &&
               IS_EQUAL(scalefactor, (double) ((float) scalefactor)) )
            {
              astype = NC_FLOAT;
            }

          if ( xtype == (int) NC_FLOAT ) astype = NC_FLOAT;

          cdf_put_att_double(fileID, ncvarid, "add_offset",   (nc_type) astype, 1, &addoffset);
          cdf_put_att_double(fileID, ncvarid, "scale_factor", (nc_type) astype, 1, &scalefactor);
        }
    }

723
  if ( dtype == CDI_DATATYPE_UINT8 && xtype == NC_BYTE )
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
    {
      int validrange[2] = {0, 255};
      cdf_put_att_int(fileID, ncvarid, "valid_range", NC_SHORT, 2, validrange);
      cdf_put_att_text(fileID, ncvarid, "_Unsigned", 4, "true");
    }

  streamptr->vars[varID].ncvarid = ncvarid;

  if ( vlistInqVarMissvalUsed(vlistID, varID) )
    cdfDefVarMissval(streamptr, varID, vlistInqVarDatatype(vlistID, varID), 0);

  if ( zid == -1 )
    {
      if ( zaxisInqType(zaxisID) == ZAXIS_CLOUD_BASE          ||
           zaxisInqType(zaxisID) == ZAXIS_CLOUD_TOP           ||
           zaxisInqType(zaxisID) == ZAXIS_ISOTHERM_ZERO       ||
           zaxisInqType(zaxisID) == ZAXIS_TOA                 ||
           zaxisInqType(zaxisID) == ZAXIS_SEA_BOTTOM          ||
           zaxisInqType(zaxisID) == ZAXIS_LAKE_BOTTOM         ||
           zaxisInqType(zaxisID) == ZAXIS_SEDIMENT_BOTTOM     ||
           zaxisInqType(zaxisID) == ZAXIS_SEDIMENT_BOTTOM_TA  ||
           zaxisInqType(zaxisID) == ZAXIS_SEDIMENT_BOTTOM_TW  ||
           zaxisInqType(zaxisID) == ZAXIS_MIX_LAYER           ||
           zaxisInqType(zaxisID) == ZAXIS_ATMOSPHERE )
        {
          zaxisInqName(zaxisID, varname);
          cdf_put_att_text(fileID, ncvarid, "level_type", strlen(varname), varname);
        }
    }

754
755
756
757
758
759
760
  int perturbationNumber, numberOfForecastsInEnsemble, typeOfEnsembleForecast;
  if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, &perturbationNumber) == 0 )
    cdf_put_att_int(fileID, ncvarid, "realization", NC_INT, 1, &perturbationNumber);
  if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, &numberOfForecastsInEnsemble) == 0 )
    cdf_put_att_int(fileID, ncvarid, "ensemble_members", NC_INT, 1, &numberOfForecastsInEnsemble);
  if ( cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast) == 0 )
    cdf_put_att_int(fileID, ncvarid, "forecast_init_type", NC_INT, 1, &typeOfEnsembleForecast);
761
762

  /* Attributes */
763
  cdfDefineAttributes(vlistID, varID, fileID, ncvarid);
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

  /* if ( streamptr->ncmode == 2 ) cdf_enddef(fileID); */

  return ncvarid;
}


void cdfEndDef(stream_t *streamptr)
{
  cdfDefGlobalAtts(streamptr);
  cdfDefLocalAtts(streamptr);

  if ( streamptr->accessmode == 0 )
    {
      int fileID  = streamptr->fileID;
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      int nvars =  streamptr->nvars;
      for ( int varID = 0; varID < nvars; varID++ )
	cdfDefVar(streamptr, varID);

      if ( streamptr->ncmode == 2 )
        {
          if ( CDI_netcdf_hdr_pad == 0UL )
            cdf_enddef(fileID);
          else
            cdf__enddef(fileID, CDI_netcdf_hdr_pad);
        }

      streamptr->accessmode = 1;
    }
}

static
void cdfWriteGridTraj(stream_t *streamptr, int gridID)
{
800
  int gridindex = nc_grid_index(streamptr, gridID);
801
  int lonID = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
802
  int latID = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_Y];
803
804
805
806
807
808

  double xlon = gridInqXval(gridID, 0);
  double xlat = gridInqYval(gridID, 0);
  int tsID = streamptr->curTsID;
  size_t index = (size_t)tsID;

809
  int fileID = streamptr->fileID;
810
811
812
813
814
  cdf_put_var1_double(fileID, lonID, &index, &xlon);
  cdf_put_var1_double(fileID, latID, &index, &xlat);
}

static
815
void cdf_write_var_data(int fileID, int vlistID, int varID, int ncvarid, int dtype, size_t nvals, size_t xsize, size_t ysize,
816
                        bool swapxy, size_t *start, size_t *count, int memtype, const void *data, size_t nmiss)
817
818
819
820
821
822
823
824
{
  const double *pdata_dp = (const double *) data;
  double *mdata_dp = NULL;
  double *sdata_dp = NULL;
  const float *pdata_sp = (const float *) data;
  float *mdata_sp = NULL;
  float *sdata_sp = NULL;

825
  /*  if ( dtype == CDI_DATATYPE_INT8 || dtype == CDI_DATATYPE_INT16 || dtype == CDI_DATATYPE_INT32 ) */
826
827
828
829
830
831
832
833
834
835
836
    {
      double missval      = vlistInqVarMissval(vlistID, varID);
      double addoffset    = vlistInqVarAddoffset(vlistID, varID);
      double scalefactor  = vlistInqVarScalefactor(vlistID, varID);
      bool laddoffset     = IS_NOT_EQUAL(addoffset, 0);
      bool lscalefactor   = IS_NOT_EQUAL(scalefactor, 1);

      if ( laddoffset || lscalefactor )
        {
          if ( memtype == MEMTYPE_FLOAT )
            {
837
838
              mdata_sp = (float *) Malloc(nvals*sizeof(float));
              memcpy(mdata_sp, pdata_sp, nvals*sizeof(float));
839
840
841
842
              pdata_sp = mdata_sp;

              if ( nmiss > 0 )
                {
843
                  for ( size_t i = 0; i < nvals; i++ )
844
845
846
847
848
849
850
851
852
853
854
855
                    {
                      double temp = mdata_sp[i];
                      if ( !DBL_IS_EQUAL(temp, missval) )
                        {
                          if ( laddoffset )   temp -= addoffset;
                          if ( lscalefactor ) temp /= scalefactor;
                          mdata_sp[i] = (float)temp;
                        }
                    }
                }
              else
                {
856
                  for ( size_t i = 0; i < nvals; i++ )
857
858
859
860
861
862
863
864
865
866
                    {
                      double temp = mdata_sp[i];
                      if ( laddoffset )   temp -= addoffset;
                      if ( lscalefactor ) temp /= scalefactor;
                      mdata_sp[i] = (float)temp;
                    }
                }
            }
          else
            {
867
868
              mdata_dp = (double *) Malloc(nvals*sizeof(double));
              memcpy(mdata_dp, pdata_dp, nvals*sizeof(double));
869
870
871
872
              pdata_dp = mdata_dp;

              if ( nmiss > 0 )
                {
873
                  for ( size_t i = 0; i < nvals; i++ )
874
875
876
877
878
879
880
881
882
883
                    {
                      if ( !DBL_IS_EQUAL(mdata_dp[i], missval) )
                        {
                          if ( laddoffset )   mdata_dp[i] -= addoffset;
                          if ( lscalefactor ) mdata_dp[i] /= scalefactor;
                        }
                    }
                }
              else
                {
884
                  for ( size_t i = 0; i < nvals; i++ )
885
886
887
888
889
890
891
892
                    {
                      if ( laddoffset )   mdata_dp[i] -= addoffset;
                      if ( lscalefactor ) mdata_dp[i] /= scalefactor;
                    }
                }
            }
        }

893
894
      if ( dtype == CDI_DATATYPE_UINT8 || dtype == CDI_DATATYPE_INT8 ||
           dtype == CDI_DATATYPE_INT16 || dtype == CDI_DATATYPE_INT32 )
895
896
897
898
899
        {
          if ( memtype == MEMTYPE_FLOAT )
            {
              if ( mdata_sp == NULL )
                {
900
901
                  mdata_sp = (float *) Malloc(nvals*sizeof(float));
                  memcpy(mdata_sp, pdata_sp, nvals*sizeof(float));
902
903
904
                  pdata_sp = mdata_sp;
                }

905
              for ( size_t i = 0; i < nvals; i++ ) mdata_sp[i] = roundf(mdata_sp[i]);
906

907
              if ( dtype == CDI_DATATYPE_UINT8 )
908
909
910
911
912
                {
                  nc_type xtype;
                  cdf_inq_vartype(fileID, ncvarid, &xtype);
                  if ( xtype == NC_BYTE )
                    {
913
                      for ( size_t i = 0; i < nvals; ++i )
914
915
916
917
918
919
920
921
                        if ( mdata_sp[i] > 127 ) mdata_sp[i] -= 256;
                    }
                }
            }
          else
            {
              if ( mdata_dp == NULL )
                {
922
923
                  mdata_dp = (double *) Malloc(nvals*sizeof(double));
                  memcpy(mdata_dp, pdata_dp, nvals*sizeof(double));
924
925
926
                  pdata_dp = mdata_dp;
                }

927
              for ( size_t i = 0; i < nvals; i++ ) mdata_dp[i] = round(mdata_dp[i]);
928

929
              if ( dtype == CDI_DATATYPE_UINT8 )
930
931
932
933
934
                {
                  nc_type xtype;
                  cdf_inq_vartype(fileID, ncvarid, &xtype);
                  if ( xtype == NC_BYTE )
                    {
935
                      for ( size_t i = 0; i < nvals; ++i )
936
937
938
939
940
941
942
943
                        if ( mdata_dp[i] > 127 ) mdata_dp[i] -= 256;
                    }
                }
            }
        }

      if ( CDF_Debug && memtype != MEMTYPE_FLOAT )
        {
944
945
          double fmin =  1.0e200;
          double fmax = -1.0e200;
946
          for ( size_t i = 0; i < nvals; ++i )
947
948
949
950
951
952
953
            {
              if ( !DBL_IS_EQUAL(pdata_dp[i], missval) )
                {
                  if ( pdata_dp[i] < fmin ) fmin = pdata_dp[i];
                  if ( pdata_dp[i] > fmax ) fmax = pdata_dp[i];
                }
            }
954
          Message("nvals = %zu, nmiss = %d, missval = %g, minval = %g, maxval = %g",
955
956
957
958
                  nvals, nmiss, missval, fmin, fmax);
        }
    }

959
  if ( swapxy ) // implemented only for cdf_write_var_slice()
960
    {
961
      size_t gridsize = xsize*ysize;
962
963
      if ( memtype == MEMTYPE_FLOAT )
        {
964
          sdata_sp = (float *) Malloc(gridsize*sizeof(float));
965
966
967
968
969
970
971
          for ( size_t j = 0; j < ysize; ++j )
            for ( size_t i = 0; i < xsize; ++i )
              sdata_sp[i*ysize+j] = pdata_sp[j*xsize+i];
          pdata_sp = sdata_sp;
        }
      else
        {
972
          sdata_dp = (double *) Malloc(gridsize*sizeof (double));
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
          for ( size_t j = 0; j < ysize; ++j )
            for ( size_t i = 0; i < xsize; ++i )
              sdata_dp[i*ysize+j] = pdata_dp[j*xsize+i];
          pdata_dp = sdata_dp;
        }
    }

  if ( memtype == MEMTYPE_FLOAT )
    cdf_put_vara_float(fileID, ncvarid, start, count, pdata_sp);
  else
    cdf_put_vara_double(fileID, ncvarid, start, count, pdata_dp);

  if ( mdata_dp ) Free(mdata_dp);
  if ( sdata_dp ) Free(sdata_dp);
  if ( mdata_sp ) Free(mdata_sp);
  if ( sdata_sp ) Free(sdata_sp);
}


992
void cdf_write_var(stream_t *streamptr, int varID, int memtype, const void *data, size_t nmiss)
993
{
994
995
  if ( streamptr->accessmode == 0 ) cdfEndDef(streamptr);

996
997
998
999
  size_t xsize = 0, ysize = 0;
  size_t size;
  size_t start[5];
  size_t count[5];
1000
  bool swapxy = false;
For faster browsing, not all history is shown. View entire blame