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

#ifdef HAVE_LIBNETCDF

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


#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2


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


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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
  void *data = Malloc((size_t)datasize
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
             * (memtype == MEMTYPE_DOUBLE ? sizeof(double) : sizeof(float)));

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

  Free(data);
}


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

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

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

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

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

  size_t index = (size_t)tsID;

  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;
  cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);

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

      ncvarid = streamptr->basetime.ncvarboundsid;

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

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

  ncvarid = streamptr->basetime.leadtimeid;
92
  if ( taxis->type == TAXIS_FORECAST && ncvarid != CDI_UNDEFID )
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    {
      timevalue = taxis->fc_period;
      cdf_put_var1_double(fileID, ncvarid, &index, &timevalue);
    }

  /*
printf("fileID = %d %d %d %f\n", fileID, time_varid, index, timevalue);
  */
}

static
int cdfDefTimeBounds(int fileID, int nctimevarid, int nctimedimid, const char *taxis_name, taxis_t* taxis)
{
  int time_bndsid = -1;
  int dims[2];

  dims[0] = nctimedimid;

  /* fprintf(stderr, "time has bounds\n"); */

  if ( nc_inq_dimid(fileID, bndsName, &dims[1]) != NC_NOERR )
    cdf_def_dim(fileID, bndsName, 2, &dims[1]);

  const char *bndsAttName, *bndsAttVal;
  size_t bndsAttValLen;
  char tmpstr[CDI_MAX_NAME];
  if ( taxis->climatology )
    {
      static const char climatology_bndsName[] = "climatology_bnds",
        climatology_bndsAttName[] = "climatology";
      bndsAttName = climatology_bndsAttName;
      bndsAttValLen = sizeof (climatology_bndsName) - 1;
      bndsAttVal = climatology_bndsName;
    }
  else
    {
      size_t taxisnameLen = strlen(taxis_name);
      memcpy(tmpstr, taxis_name, taxisnameLen);
      tmpstr[taxisnameLen] = '_';
      memcpy(tmpstr + taxisnameLen + 1, bndsName, sizeof (bndsName));
      size_t tmpstrLen = taxisnameLen + sizeof (bndsName);
      static const char generic_bndsAttName[] = "bounds";
      bndsAttName = generic_bndsAttName;
      bndsAttValLen = tmpstrLen;
      bndsAttVal = tmpstr;
    }
  cdf_def_var(fileID, bndsAttVal, NC_DOUBLE, 2, dims, &time_bndsid);
  cdf_put_att_text(fileID, nctimevarid, bndsAttName, bndsAttValLen, bndsAttVal);

  return time_bndsid;
}

static
void cdfDefTimeUnits(char *unitstr, taxis_t* taxis0, taxis_t* taxis)
{
  unitstr[0] = 0;

  if ( taxis0->type == TAXIS_ABSOLUTE )
    {
      if ( taxis0->unit == TUNIT_YEAR )
        sprintf(unitstr, "year as %s", "%Y.%f");
      else if ( taxis0->unit == TUNIT_MONTH )
        sprintf(unitstr, "month as %s", "%Y%m.%f");
      else
        sprintf(unitstr, "day as %s", "%Y%m%d.%f");
    }
  else
    {
      int timeunit = taxis->unit != -1 ? taxis->unit : TUNIT_HOUR;
      int rdate    = taxis->rdate;
      int rtime    = taxis->rtime;
      if ( rdate == -1 )
        {
          rdate  = taxis->vdate;
          rtime  = taxis->vtime;
        }

      int year, month, day, hour, minute, second;
      cdiDecodeDate(rdate, &year, &month, &day);
      cdiDecodeTime(rtime, &hour, &minute, &second);

      if ( timeunit == TUNIT_QUARTER   ) timeunit = TUNIT_MINUTE;
      if ( timeunit == TUNIT_30MINUTES ) timeunit = TUNIT_MINUTE;
      if ( timeunit == TUNIT_3HOURS  ||
	   timeunit == TUNIT_6HOURS  ||
	   timeunit == TUNIT_12HOURS ) timeunit = TUNIT_HOUR;

      sprintf(unitstr, "%s since %d-%d-%d %02d:%02d:%02d",
              tunitNamePtr(timeunit), year, month, day, hour, minute, second);
    }
}

static
void cdfDefForecastTimeUnits(char *unitstr, int timeunit)
{
  unitstr[0] = 0;

  if ( timeunit == -1 ) timeunit = TUNIT_HOUR;

  if ( timeunit == TUNIT_QUARTER   ) timeunit = TUNIT_MINUTE;
  if ( timeunit == TUNIT_30MINUTES ) timeunit = TUNIT_MINUTE;
  if ( timeunit == TUNIT_3HOURS  ||
       timeunit == TUNIT_6HOURS  ||
       timeunit == TUNIT_12HOURS ) timeunit = TUNIT_HOUR;

  strcpy(unitstr, tunitNamePtr(timeunit));
}

static
void cdfDefCalendar(int fileID, int ncvarid, int calendar)
{
  static const struct { int calCode; const char *calStr; } calTab[] = {
    { CALENDAR_STANDARD, "standard" },
    { CALENDAR_PROLEPTIC, "proleptic_gregorian" },
    { CALENDAR_NONE, "none" },
    { CALENDAR_360DAYS, "360_day" },
    { CALENDAR_365DAYS, "365_day" },
    { CALENDAR_366DAYS, "366_day" },
  };
  enum { calTabSize = sizeof calTab / sizeof calTab[0] };

214
215
  for ( size_t i = 0; i < calTabSize; ++i )
    if ( calTab[i].calCode == calendar )
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
      {
        const char *calstr = calTab[i].calStr;
        size_t len = strlen(calstr);
        cdf_put_att_text(fileID, ncvarid, "calendar", len, calstr);
        break;
      }
}


void cdfDefTime(stream_t* streamptr)
{
  int time_varid;
  int time_dimid;
  int time_bndsid = -1;
  static const char default_name[] = "time";

232
  if ( streamptr->basetime.ncvarid != CDI_UNDEFID ) return;
233
234
235
236
237
238
239
240
241
242
243
244
245

  int fileID = streamptr->fileID;

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

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

  const char *taxis_name = (taxis->name && taxis->name[0]) ? taxis->name : default_name ;

  cdf_def_dim(fileID, taxis_name, NC_UNLIMITED, &time_dimid);
  streamptr->basetime.ncdimid = time_dimid;

246
  nc_type xtype = (taxis->datatype == CDI_DATATYPE_FLT32) ? NC_FLOAT : NC_DOUBLE;
247
248

  cdf_def_var(fileID, taxis_name, xtype, 1, &time_dimid, &time_varid);
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292

  streamptr->basetime.ncvarid = time_varid;

  {
    static const char timeStr[] = "time";
    cdf_put_att_text(fileID, time_varid, "standard_name", sizeof(timeStr) - 1, timeStr);
  }

  if ( taxis->longname && taxis->longname[0] )
    cdf_put_att_text(fileID, time_varid, "long_name", strlen(taxis->longname), taxis->longname);

  if ( taxis->has_bounds )
    {
      time_bndsid = cdfDefTimeBounds(fileID, time_varid, time_dimid, taxis_name, taxis);
      streamptr->basetime.ncvarboundsid = time_bndsid;
    }

  {
    char unitstr[CDI_MAX_NAME];
    cdfDefTimeUnits(unitstr, &streamptr->tsteps[0].taxis, taxis);
    size_t len = strlen(unitstr);
    if ( len )
      {
        cdf_put_att_text(fileID, time_varid, "units", len, unitstr);
        /*
          if ( taxis->has_bounds )
          cdf_put_att_text(fileID, time_bndsid, "units", len, unitstr);
        */
      }
  }

  if ( taxis->calendar != -1 )
    {
      cdfDefCalendar(fileID, time_varid, taxis->calendar);
      /*
      if ( taxis->has_bounds )
        cdfDefCalendar(fileID, time_bndsid, taxis->calendar);
      */
    }

  if ( taxis->type == TAXIS_FORECAST )
    {
      int leadtimeid;

293
      cdf_def_var(fileID, "leadtime", xtype, 1, &time_dimid, &leadtimeid);
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331

      streamptr->basetime.leadtimeid = leadtimeid;

      {
        static const char stdname[] = "forecast_period";
        cdf_put_att_text(fileID, leadtimeid, "standard_name", sizeof(stdname) - 1, stdname);
      }

      {
        static const char lname[] = "Time elapsed since the start of the forecast";
        cdf_put_att_text(fileID, leadtimeid, "long_name", sizeof(lname) - 1, lname);
      }

      {
          char unitstr[CDI_MAX_NAME];
          cdfDefForecastTimeUnits(unitstr, taxis->fc_unit);
          size_t len = strlen(unitstr);
          if ( len )
            cdf_put_att_text(fileID, leadtimeid, "units", len, unitstr);
      }
    }

  cdf_put_att_text(fileID, time_varid, "axis", 1, "T");

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


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

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

  cdfDefTimeValue(streamptr, tsID);
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
332
void cdfDefComplex(stream_t *streamptr, int gridID, int gridindex)
333
{
334
  int dimID = CDI_UNDEFID;
335
  int fileID  = streamptr->fileID;
336
  ncgrid_t *ncgrid = streamptr->ncgrid;
337

Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
  for ( int index = 0; index < gridindex; ++index )
339
    {
340
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
341
        {
342
          int gridID0 = ncgrid[index].gridID;
343
344
345
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_SPECTRAL || gridtype0 == GRID_FOURIER )
            {
346
              dimID = ncgrid[index].xdimID;
347
348
349
350
351
              break;
            }
        }
    }

352
  if ( dimID == CDI_UNDEFID )
353
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
354
      static const char axisname[] = "nc2";
355
356
357
358
359
      size_t dimlen = 2;

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
      cdf_def_dim(fileID, axisname, dimlen, &dimID);
      cdf_enddef(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360

361
362
363
      streamptr->ncmode = 2;
    }

364
365
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
366
367
368
}

static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
cdfDefSPorFC(stream_t *streamptr, int gridID, int gridindex,
370
371
             char *restrict axisname, int gridRefType)
{
372
  int iz = 0;
373
  int dimID = CDI_UNDEFID;
374
  ncgrid_t *ncgrid = streamptr->ncgrid;
375
376
377

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

378
  for ( int index = 0; index < gridindex; index++ )
379
    {
380
      if ( ncgrid[index].ydimID != CDI_UNDEFID )
381
        {
382
          int gridID0 = ncgrid[index].gridID;
383
384
385
386
387
388
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == gridRefType )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0)/2;
              if ( dimlen == dimlen0 )
                {
389
                  dimID = ncgrid[index].ydimID;
390
391
392
393
394
395
396
397
                  break;
                }
              else
                iz++;
            }
        }
    }

398
  if ( dimID == CDI_UNDEFID )
399
400
401
402
403
404
405
406
407
408
409
410
411
    {
      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;
    }

412
413
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].ydimID = dimID;
414
415
416
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
void cdfDefSP(stream_t *streamptr, int gridID, int gridindex)
418
419
420
421
422
{
  /*
  char longname[] = "Spherical harmonic coefficient";
  */
  char axisname[5] = "nspX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_SPECTRAL);
424
425
426
427
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
void cdfDefFC(stream_t *streamptr, int gridID, int gridindex)
429
430
{
  char axisname[5] = "nfcX";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_FOURIER);
432
433
434
435
}

static const struct cdfDefGridAxisInqs {
  int (*axisSize)(int gridID);
436
437
  int (*axisName)(int cdiID, int key, int size, char *mesg);
  int (*axisLongname)(int cdiID, int key, int size, char *mesg);
438
439
  int (*axisUnits)(int cdiID, int key, int size, char *mesg);
  void (*axisStdname)(int cdiID, char *dimstdname);
440
441
442
443
444
  double (*axisVal)(int gridID, int index);
  const double *(*axisValsPtr)(int gridID);
  const double *(*axisBoundsPtr)(int gridID);
} gridInqsX = {
  .axisSize = gridInqXsize,
445
  .axisName = cdiGridInqKeyStr,
446
  .axisLongname = cdiGridInqKeyStr,
447
448
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqXstdname,
449
450
451
452
453
  .axisVal = gridInqXval,
  .axisValsPtr = gridInqXvalsPtr,
  .axisBoundsPtr = gridInqXboundsPtr,
}, gridInqsY = {
  .axisSize = gridInqYsize,
454
  .axisName = cdiGridInqKeyStr,
455
  .axisLongname = cdiGridInqKeyStr,
456
457
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqYstdname,
458
459
460
461
  .axisVal = gridInqYval,
  .axisValsPtr = gridInqYvalsPtr,
  .axisBoundsPtr = gridInqYboundsPtr,
}, gridInqsZ = {
462
  .axisLongname = cdiZaxisInqKeyStr,
463
464
  .axisUnits = cdiZaxisInqKeyStr,
  .axisStdname = zaxisInqStdname,
465
466
};

467
468
static
void cdfPutGridStdAtts(int fileID, int ncvarid, int gridID, int dimtype, const struct cdfDefGridAxisInqs *inqs)
469
470
{
  size_t len;
471
472
473
474
475
476
477
478
479

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

483
484
485
486
  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)) )
487
    cdf_put_att_text(fileID, ncvarid, "units", len, units);
488
489
490
}

static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
cdfDefTrajLatLon(stream_t *streamptr, int gridID, int gridindex,
492
                 const struct cdfDefGridAxisInqs *inqs, int dimtype)
493
{
494
  nc_type xtype = (gridInqPrec(gridID) == CDI_DATATYPE_FLT32) ? NC_FLOAT : NC_DOUBLE;
495
  ncgrid_t *ncgrid = streamptr->ncgrid;
496
497
498

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

501
  int ncvarid = CDI_UNDEFID;
502
503
504
505
  if ( dimtype == 'X' )
    ncvarid = ncgrid[gridindex].xdimID;
  else
    ncvarid = ncgrid[gridindex].ydimID;
506

507
  if ( ncvarid == CDI_UNDEFID )
508
509
510
511
512
    {
      int dimNcID = streamptr->basetime.ncvarid;
      int fileID  = streamptr->fileID;
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

513
514
515
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
      int keyname = (dimtype == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
      inqs->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
516
      cdf_def_var(fileID, axisname, xtype, 1, &dimNcID, &ncvarid);
517
      cdfPutGridStdAtts(fileID, ncvarid, gridID, dimtype, inqs);
518
519
520
521
      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

522
523
524
525
526
  ncgrid[gridindex].gridID = gridID;
  if ( dimtype == 'X' )
    ncgrid[gridindex].xdimID = ncvarid; /* var ID for trajectory !!! */
  else
    ncgrid[gridindex].ydimID = ncvarid; /* var ID for trajectory !!! */
527
528
529
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
void cdfDefTrajLon(stream_t *streamptr, int gridID, int gridindex)
531
{
532
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsX, 'X');
533
534
535
536
}


static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
void cdfDefTrajLat(stream_t *streamptr, int gridID, int gridindex)
538
{
539
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsY, 'Y');
540
541
542
543
544
545
546
}

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;
547
  int dimid = CDI_UNDEFID;
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
  char name[CDI_MAX_NAME];

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

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

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


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

  return dimid;
}

static
578
void checkGridName(char *axisname, int fileID)
579
580
581
582
583
584
585
586
587
588
589
590
591
592
{
  int ncdimid;
  char axisname2[CDI_MAX_NAME];

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

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

      int status = nc_inq_varid(fileID, axisname2, &ncdimid);
593
      if ( status != NC_NOERR ) break;
594
595
596
597
598
599
600
601
602
603
604
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
642
643
644
645
646

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

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

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

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

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

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

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


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

  return (int)iz;
}

static void
647
cdfDefAxisCommon(stream_t *streamptr, int gridID, int gridindex, int ndims,
648
649
                 const struct cdfDefGridAxisInqs *gridAxisInq, int dimKey, char axisLetter,
                 void (*finishCyclicBounds)(double *pbounds, size_t dimlen, const double *pvals))
650
{
651
652
653
  int dimID = CDI_UNDEFID;
  int ncvarid = CDI_UNDEFID, ncbvarid = CDI_UNDEFID;
  int nvdimID = CDI_UNDEFID;
654
655
  int fileID  = streamptr->fileID;
  size_t dimlen = (size_t)gridAxisInq->axisSize(gridID);
656
657
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);

658
  ncgrid_t *ncgrid = streamptr->ncgrid;
659

Uwe Schulzweida's avatar
Uwe Schulzweida committed
660
  for ( int index = 0; index < gridindex; ++index )
661
    {
662
      int gridID0 = ncgrid[index].gridID;
663
      assert(gridID0 != CDI_UNDEFID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
664
665
666
667
668
669
      int gridtype0 = gridInqType(gridID0);
      if ( gridtype0 == GRID_GAUSSIAN    ||
           gridtype0 == GRID_LONLAT      ||
           gridtype0 == GRID_PROJECTION  ||
           gridtype0 == GRID_CURVILINEAR ||
           gridtype0 == GRID_GENERIC )
670
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
672
          size_t dimlen0 = (size_t)gridAxisInq->axisSize(gridID0);
          if ( dimlen == dimlen0 )
673
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
675
676
              double (*inqVal)(int gridID, int index) = gridAxisInq->axisVal;
              if ( IS_EQUAL(inqVal(gridID0, 0), inqVal(gridID, 0)) &&
                   IS_EQUAL(inqVal(gridID0, (int)dimlen-1), inqVal(gridID, (int)dimlen-1)) )
677
                {
678
                  dimID = (dimKey == CDI_KEY_XDIMNAME) ? ncgrid[index].xdimID : ncgrid[index].ydimID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
679
                  break;
680
681
682
683
684
                }
            }
        }
    }

685
  if ( dimID == CDI_UNDEFID )
686
687
688
    {
      const double *pvals = gridAxisInq->axisValsPtr(gridID);

689
690
691
      char axisname[CDI_MAX_NAME]; axisname[0] = 0;
      int keyname = (axisLetter == 'X') ? CDI_KEY_XNAME : CDI_KEY_YNAME;
      gridAxisInq->axisName(gridID, keyname, CDI_MAX_NAME, axisname);
692
693
694
      if ( axisname[0] == 0 ) Error("axis name undefined!");
      size_t axisnameLen = strlen(axisname);

695
696
697
      /* enough to append _ plus up to 100 decimal and trailing \0 */
      char extendedAxisname[axisnameLen + 4 + 1];
      memcpy(extendedAxisname, axisname, axisnameLen + 1);
698
      checkGridName(extendedAxisname, fileID);
699
      size_t extendedAxisnameLen = axisnameLen + strlen(extendedAxisname + axisnameLen);
700
701
702
703
704

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

      if ( ndims )
        {
705
          char dimname[CDI_MAX_NAME+3]; dimname[0] = 0;
706
707

          if ( pvals == NULL )
708
            cdiGridInqKeyStr(gridID, dimKey, CDI_MAX_NAME, dimname);
709
710
711
712

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

713
          if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);
714
715
716
        }

      bool gen_bounds = false;
717
      bool grid_is_cyclic = gridIsCircular(gridID) > 0;
718
719
720
721
722
      double *pbounds = NULL;
      if ( pvals )
        {
          cdf_def_var(fileID, extendedAxisname, xtype, ndims, &dimID, &ncvarid);

723
          cdfPutGridStdAtts(fileID, ncvarid, gridID, axisLetter, gridAxisInq);
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
          {
            char axisStr[2] = { axisLetter, '\0' };
            cdf_put_att_text(fileID, ncvarid, "axis", 1, axisStr);
          }

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

          if ( CDI_cmor_mode && grid_is_cyclic && !pbounds )
            {
              gen_bounds = true;
              pbounds = (double*) Malloc(2*dimlen*sizeof(double));
              for ( size_t i = 0; i < dimlen-1; ++i )
                {
                  pbounds[i*2+1]   = (pvals[i] + pvals[i+1])/2;
                  pbounds[(i+1)*2] = (pvals[i] + pvals[i+1])/2;
                }
              finishCyclicBounds(pbounds, dimlen, pvals);
            }
          if ( pbounds )
            {
              size_t nvertex = 2;
              if ( nc_inq_dimid(fileID, bndsName, &nvdimID) != NC_NOERR )
                cdf_def_dim(fileID, bndsName, nvertex, &nvdimID);
            }
748
          if ( pbounds && nvdimID != CDI_UNDEFID )
749
750
751
752
753
754
755
756
757
758
759
760
761
762
            {
              char boundsname[extendedAxisnameLen + 1 + sizeof (bndsName)];
              memcpy(boundsname, axisname, extendedAxisnameLen);
              boundsname[extendedAxisnameLen] = '_';
              memcpy(boundsname + extendedAxisnameLen + 1, bndsName, sizeof bndsName);
              int dimIDs[2] = { dimID, nvdimID };
              cdf_def_var(fileID, boundsname, xtype, 2, dimIDs, &ncbvarid);
              cdf_put_att_text(fileID, ncvarid, "bounds", extendedAxisnameLen + sizeof (bndsName), boundsname);
            }
        }

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

763
764
      if ( ncvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncvarid, pvals);
      if ( ncbvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbvarid, pbounds);
765
766
      if ( gen_bounds ) Free(pbounds);

767
768
769
770
771
772
773
      if ( ndims == 0 )
        {
          if ( dimKey == CDI_KEY_XDIMNAME )
            ncgrid[gridindex].xvarID = ncvarid;
          else
            ncgrid[gridindex].yvarID = ncvarid;
        }
774
775
    }

776
777
778
779
780
  ncgrid[gridindex].gridID = gridID;
  if ( dimKey == CDI_KEY_XDIMNAME )
    ncgrid[gridindex].xdimID = dimID;
  else
    ncgrid[gridindex].ydimID = dimID;
781
782
}

783
784
static
void finishCyclicXBounds(double *pbounds, size_t dimlen, const double *pvals)
785
786
787
788
789
790
{
  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
791
void cdfDefXaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
792
{
793
794
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsX,
                   CDI_KEY_XDIMNAME, 'X', finishCyclicXBounds);
795
796
}

797
798
static
void finishCyclicYBounds(double *pbounds, size_t dimlen, const double *pvals)
799
800
801
802
803
804
{
  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
805
void cdfDefYaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
806
{
807
808
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsY,
                   CDI_KEY_YDIMNAME, 'Y', finishCyclicYBounds);
809
810
811
812
813
814
}

static
void cdfGridCompress(int fileID, int ncvarid, int gridsize, int filetype, int comptype)
{
#if  defined  (HAVE_NETCDF4)
815
  if ( gridsize > 1 && comptype == CDI_COMPRESS_ZIP && (filetype == CDI_FILETYPE_NC4 || filetype == CDI_FILETYPE_NC4C) )
816
817
818
819
820
821
822
823
    {
      nc_def_var_chunking(fileID, ncvarid, NC_CHUNKED, NULL);
      cdfDefVarDeflate(fileID, ncvarid, 1);
    }
#endif
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824
void cdfDefCurvilinear(stream_t *streamptr, int gridID, int gridindex)
825
{
826
827
828
829
  int xdimID = CDI_UNDEFID;
  int ydimID = CDI_UNDEFID;
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID;
  int ncbxvarid = CDI_UNDEFID, ncbyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
830
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);
831
  ncgrid_t *ncgrid = streamptr->ncgrid;
832
833
834
835
836
837
838

  int fileID  = streamptr->fileID;

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
  for ( int index = 0; index < gridindex; index++ )
840
    {
841
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
842
        {
843
          int gridID0 = ncgrid[index].gridID;
844
845
846
847
848
849
850
851
852
853
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_CURVILINEAR )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0);
              if ( dimlen == dimlen0 )
                if ( IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0)) &&
                     IS_EQUAL(gridInqXval(gridID0, (int)dimlen-1), gridInqXval(gridID, (int)dimlen-1)) &&
                     IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0)) &&
                     IS_EQUAL(gridInqYval(gridID0, (int)dimlen-1), gridInqYval(gridID, (int)dimlen-1)) )
                  {
854
855
856
857
                    xdimID = ncgrid[index].xdimID;
                    ydimID = ncgrid[index].ydimID;
                    ncxvarid = ncgrid[index].xvarID;
                    ncyvarid = ncgrid[index].yvarID;
858
859
860
861
862
863
                    break;
                  }
            }
        }
    }

864
  if ( xdimID == CDI_UNDEFID || ydimID == CDI_UNDEFID )
865
866
867
868
869
    {
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
      {
        char xdimname[CDI_MAX_NAME+3];
        xdimname[0] = 0;
870
        cdiGridInqKeyStr(gridID, CDI_KEY_XDIMNAME, CDI_MAX_NAME, xdimname);
871
872
        if ( xdimname[0] == 0 ) { xdimname[0] = 'x'; xdimname[1] = 0; }
        xdimID = checkDimName(fileID, xdimlen, xdimname);
873
        if ( xdimID == CDI_UNDEFID ) cdf_def_dim(fileID, xdimname, xdimlen, &xdimID);
874
875
876
877
      }
      {
        char ydimname[CDI_MAX_NAME+3];
        ydimname[0] = 0;
878
        cdiGridInqKeyStr(gridID, CDI_KEY_YDIMNAME, CDI_MAX_NAME, ydimname);
879
880
        if ( ydimname[0] == 0 ) { ydimname[0] = 'y'; ydimname[1] = 0; }
        ydimID = checkDimName(fileID, ydimlen, ydimname);
881
        if ( ydimID == CDI_UNDEFID ) cdf_def_dim(fileID, ydimname, ydimlen, &ydimID);
882
883
      }

884
      int nvdimID = CDI_UNDEFID;
885
886
887
      int dimIDs[3];
      if ( gridInqXboundsPtr(gridID) || gridInqYboundsPtr(gridID) )
        {
888
          char vdimname[CDI_MAX_NAME+3]; vdimname[0] = 0;
889
          cdiGridInqKeyStr(gridID, CDI_KEY_VDIMNAME, CDI_MAX_NAME, vdimname);
890
891
892
          if ( vdimname[0] == 0 ) strcpy(vdimname, "nv4");
          size_t nvertex = 4;
          nvdimID = checkDimName(fileID, nvertex, vdimname);
893
          if ( nvdimID == CDI_UNDEFID ) cdf_def_dim(fileID, vdimname, nvertex, &nvdimID);
894
895
896
897
898
899
900
901
        }

      dimIDs[0] = ydimID;
      dimIDs[1] = xdimID;
      dimIDs[2] = nvdimID;

      if ( gridInqXvalsPtr(gridID) )
        {
902
903
          char xaxisname[CDI_MAX_NAME]; xaxisname[0] = 0;
          cdiGridInqKeyStr(gridID, CDI_KEY_XNAME, CDI_MAX_NAME, xaxisname);
904
          checkGridName(xaxisname, fileID);
905
906
907
908

          cdf_def_var(fileID, xaxisname, xtype, 2, dimIDs, &ncxvarid);
          cdfGridCompress(fileID, ncxvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

909
          cdfPutGridStdAtts(fileID, ncxvarid, gridID, 'X', &gridInqsX);
910
911
912
913

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

914
          if ( gridInqXboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
            {
              size_t xaxisnameLen = strlen(xaxisname);
              xaxisname[xaxisnameLen] = '_';
              memcpy(xaxisname + xaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, xaxisname, xtype, 3, dimIDs, &ncbxvarid);
              cdfGridCompress(fileID, ncbxvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

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

      if ( gridInqYvalsPtr(gridID) )
        {
          char yaxisname[CDI_MAX_NAME];
          gridInqYname(gridID, yaxisname);
930
          checkGridName(yaxisname, fileID);
931
932
933
934

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

935
          cdfPutGridStdAtts(fileID, ncyvarid, gridID, 'Y', &gridInqsY);
936
937
938
939

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

940
          if ( gridInqYboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
            {
              size_t yaxisnameLen = strlen(yaxisname);
              yaxisname[yaxisnameLen] = '_';
              memcpy(yaxisname + yaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, yaxisname, xtype, 3, dimIDs, &ncbyvarid);
              cdfGridCompress(fileID, ncbyvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

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

      if ( gridInqAreaPtr(gridID) )
        {
          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, 2, 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;

969
970
971
972
973
      if ( ncxvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncxvarid,  gridInqXvalsPtr(gridID));
      if ( ncbxvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbxvarid, gridInqXboundsPtr(gridID));
      if ( ncyvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncyvarid,  gridInqYvalsPtr(gridID));
      if ( ncbyvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbyvarid, gridInqYboundsPtr(gridID));
      if ( ncavarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncavarid,  gridInqAreaPtr(gridID));
974
975
    }

976
977
978
979
980
981
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = xdimID;
  ncgrid[gridindex].ydimID = ydimID;
  ncgrid[gridindex].xvarID = ncxvarid;
  ncgrid[gridindex].yvarID = ncyvarid;
  ncgrid[gridindex].avarID = ncavarid;
982
983
984
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
985
void cdfDefRgrid(stream_t *streamptr, int gridID, int gridindex)
986
{
987
  ncgrid_t *ncgrid = streamptr->ncgrid;
988
  int dimID = CDI_UNDEFID;
989
990
991
992

  size_t dimlen = (size_t)gridInqSize(gridID);

  int iz = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
  for ( int index = 0; index < gridindex; index++ )
994
    {
995
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
996
        {
997
          int gridID0 = ncgrid[index].gridID;
998
999
1000
1001
1002
1003
1004
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_GAUSSIAN_REDUCED )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0);

              if ( dimlen == dimlen0 )
                {
1005
                  dimID = ncgrid[index].xdimID;
1006
1007
1008
1009
1010
1011
1012
                  break;
                }
              iz++;
            }
        }
    }

1013
  if ( dimID == CDI_UNDEFID )
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
    {
      int fileID  = streamptr->fileID;
      static bool lwarn = true;
      if ( lwarn )
        {
          Warning("Creating a NetCDF file with data on a gaussian reduced grid.");
          Warning("The further processing of the resulting file is unsupported!");
          lwarn = false;
        }

      char axisname[7] = "rgridX";
      if ( iz == 0 ) axisname[5] = '\0';
      else           sprintf(&axisname[5], "%1d", iz+1);

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

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

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

1036
1037
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
1038
1039
1040
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1041
void cdfDefGdim(stream_t *streamptr, int gridID, int gridindex)
1042
{
1043
  ncgrid_t *ncgrid = streamptr->ncgrid;
1044
  int iz = 0;
1045
  int dimID = CDI_UNDEFID;
1046
1047
1048
1049

  size_t dimlen = (size_t)gridInqSize(gridID);

  if ( gridInqYsize(gridID) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1050
    for ( int index = 0; index < gridindex; index++ )
1051
      {
1052
        if ( ncgrid[index].xdimID != CDI_UNDEFID )
1053
          {
1054
            int gridID0 = ncgrid[index].gridID;
1055
1056
1057
1058
1059
1060
            int gridtype0 = gridInqType(gridID0);
            if ( gridtype0 == GRID_GENERIC )
              {
                size_t dimlen0 = (size_t)gridInqSize(gridID0);
                if ( dimlen == dimlen0 )
                  {
1061
                    dimID = ncgrid[index].xdimID;
1062
1063
1064
1065
1066
1067
1068
1069
1070
                    break;
                  }
                else
                  iz++;
              }
          }
      }

  if ( gridInqXsize(gridID) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1071
    for ( int index = 0; index < gridindex; index++ )
1072
      {
1073
        if ( ncgrid[index].ydimID != CDI_UNDEFID )
1074
          {
1075
            int gridID0 = ncgrid[index].gridID;
1076
1077
1078
1079
1080
1081
            int gridtype0 = gridInqType(gridID0);
            if ( gridtype0 == GRID_GENERIC )
              {
                size_t dimlen0 = (size_t)gridInqSize(gridID0);
                if ( dimlen == dimlen0 )
                  {
1082
                    dimID = ncgrid[index].ydimID;
1083
1084
1085
1086
1087
1088
1089
1090
                    break;
                  }
                else
                  iz++;
              }
          }
      }

1091
  if ( dimID == CDI_UNDEFID )
1092
1093
1094
1095
1096
1097
1098
1099
1100
    {
      int fileID  = streamptr->fileID;
      char dimname[CDI_MAX_NAME];
      strcpy(dimname, "gsize");

      dimID = checkDimName(fileID, dimlen, dimname);

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

1101
      if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);
1102
1103
1104
1105
1106

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

1107
1108
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
1109
1110
1111
1112
1113
1114
1115
1116
}

static
void cdfDefGridReference(stream_t *streamptr, int gridID)
{
  int fileID  = streamptr->fileID;
  int number = gridInqNumber(gridID);
  if ( number > 0 )
1117
    cdf_put_att_int(fileID, NC_GLOBAL, "number_of_grid_used", NC_INT, 1, &number);
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164

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

static
void cdfDefZaxisUUID(stream_t *streamptr, int zaxisID)
{
  unsigned char uuidOfVGrid[CDI_UUID_SIZE];
  zaxisInqUUID(zaxisID, uuidOfVGrid);

  if ( uuidOfVGrid[0] != 0 )
    {
      char uuidOfVGridStr[37];
      cdiUUID2Str(uuidOfVGrid, uuidOfVGridStr);
      if ( uuidOfVGridStr[0] != 0 && strlen(uuidOfVGridStr) == 36 )
        {
          int fileID  = streamptr->fileID;
          if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
          cdf_put_att_text(fileID, NC_GLOBAL, "uuidOfVGrid", 36, uuidOfVGridStr);
          if ( streamptr->ncmode == 2 ) cdf_enddef(fileID);
        }
    }
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1165
void cdfDefUnstructured(stream_t *streamptr, int gridID, int gridindex)
1166
{
1167
1168
1169
1170
  int dimID = CDI_UNDEFID;
  int ncxvarid = CDI_UNDEFID, ncyvarid = CDI_UNDEFID;
  int ncbxvarid = CDI_UNDEFID, ncbyvarid = CDI_UNDEFID, ncavarid = CDI_UNDEFID;
  int nvdimID = CDI_UNDEFID;
1171
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);
1172
  ncgrid_t *ncgrid = streamptr->ncgrid;
1173
1174
1175
1176
1177

  int fileID  = streamptr->fileID;

  size_t dimlen = (size_t)gridInqSize(gridID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1178
  for ( int index = 0; index < gridindex; index++ )
1179
    {
1180
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
1181
        {
1182
          int gridID0 = ncgrid[index].gridID;
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_UNSTRUCTURED )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0);
              if ( dimlen == dimlen0 )
		if ( gridInqNvertex(gridID0) == gridInqNvertex(gridID) &&
		     IS_EQUAL(gridInqXval(gridID0, 0), gridInqXval(gridID, 0)) &&
                     IS_EQUAL(gridInqXval(gridID0, (int)dimlen-1), gridInqXval(gridID, (int)dimlen-1)) &&
		     IS_EQUAL(gridInqYval(gridID0, 0), gridInqYval(gridID, 0)) &&
                     IS_EQUAL(gridInqYval(gridID0, (int)dimlen-1), gridInqYval(gridID, (int)dimlen-1)) )
		  {
1194
1195
1196
1197
		    dimID = ncgrid[index].xdimID;
                    ncxvarid = ncgrid[index].xvarID;
                    ncyvarid = ncgrid[index].yvarID;
                    ncavarid = ncgrid[index].avarID;
1198
1199
1200
1201
1202
1203
		    break;
		  }
            }
        }
    }

1204
  if ( dimID == CDI_UNDEFID )
1205
1206
1207
1208
1209
    {
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
      {
        char xdimname[CDI_MAX_NAME+3];
        xdimname[0] = 0;
1210
        cdiGridInqKeyStr(gridID, CDI_KEY_XDIMNAME, CDI_MAX_NAME, xdimname);
1211
1212
        if ( xdimname[0] == 0 ) strcpy(xdimname, "ncells");
        dimID = checkDimName(fileID, dimlen, xdimname);
1213
        if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, xdimname, dimlen, &dimID);
1214
1215
1216
1217
1218
1219
1220
      }

      size_t nvertex = (size_t)gridInqNvertex(gridID);
      if ( nvertex > 0 )
        {
          char vdimname[CDI_MAX_NAME+3];
          vdimname[0] = 0;
1221
          cdiGridInqKeyStr(gridID, CDI_KEY_VDIMNAME, CDI_MAX_NAME, vdimname);
1222
1223
          if ( vdimname[0] == 0 ) strcpy(vdimname, "vertices");
          nvdimID = checkDimName(fileID, nvertex, vdimname);
1224
          if ( nvdimID == CDI_UNDEFID ) cdf_def_dim(fileID, vdimname, nvertex, &nvdimID);
1225
1226
1227
1228
1229
1230
1231
1232
        }

      cdfDefGridReference(streamptr, gridID);

      cdfDefGridUUID(streamptr, gridID);

      if ( gridInqXvalsPtr(gridID) )
        {
1233
1234
          char xaxisname[CDI_MAX_NAME]; xaxisname[0] = 0;
          cdiGridInqKeyStr(gridID, CDI_KEY_XNAME, CDI_MAX_NAME, xaxisname);
1235
          checkGridName(xaxisname, fileID);
1236
1237
1238
          cdf_def_var(fileID, xaxisname, xtype, 1, &dimID, &ncxvarid);
          cdfGridCompress(fileID, ncxvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

1241
          if ( gridInqXboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
            {
              int dimIDs[2] = { dimID, nvdimID };
              size_t xaxisnameLen = strlen(xaxisname);
              xaxisname[xaxisnameLen] = '_';
              memcpy(xaxisname + xaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, xaxisname, xtype, 2, dimIDs, &ncbxvarid);
              cdfGridCompress(fileID, ncbxvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

      if ( gridInqYvalsPtr(gridID) )
        {
          char yaxisname[CDI_MAX_NAME];
          gridInqYname(gridID, yaxisname);
1258
          checkGridName(yaxisname, fileID);
1259
1260
1261
          cdf_def_var(fileID, yaxisname, xtype, 1, &dimID, &ncyvarid);
          cdfGridCompress(fileID, ncyvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

1264
          if ( gridInqYboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
            {
              int dimIDs[2] = { dimID, nvdimID };
              size_t yaxisnameLen = strlen(yaxisname);
              yaxisname[yaxisnameLen] = '_';
              memcpy(yaxisname + yaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, yaxisname, xtype, 2, dimIDs, &ncbyvarid);
              cdfGridCompress(fileID, ncbyvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

      if ( gridInqAreaPtr(gridID) )
        {
          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, 1, &dimID, &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;

1294
1295
1296
1297
1298
      if ( ncxvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncxvarid,  gridInqXvalsPtr(gridID));
      if ( ncbxvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbxvarid, gridInqXboundsPtr(gridID));
      if ( ncyvarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncyvarid,  gridInqYvalsPtr(gridID));
      if ( ncbyvarid != CDI_UNDEFID ) cdf_put_var_double(fileID, ncbyvarid, gridInqYboundsPtr(gridID));
      if ( ncavarid  != CDI_UNDEFID ) cdf_put_var_double(fileID, ncavarid,  gridInqAreaPtr(gridID));
1299
1300
    }

1301
1302
1303
1304
1305
  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
  ncgrid[gridindex].xvarID = ncxvarid;
  ncgrid[gridindex].yvarID = ncyvarid;
  ncgrid[gridindex].avarID = ncavarid;
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
}

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

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

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

      int mlev = ilev - 1;

      if ( streamptr->vct.ilev > 0 )
        {
          if ( streamptr->vct.ilev != ilev )
1329
            Error("More than one VCT for each file unsupported!");
1330
1331
1332
1333
1334
1335
1336
          return;
        }

      int fileID = streamptr->fileID;

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

1337
      int ncdimid, ncdimid2;
1338
1339
1340
1341
1342
1343
1344
1345
      cdf_def_dim(fileID, "nhym", (size_t)mlev, &ncdimid);
      cdf_def_dim(fileID, "nhyi", (size_t)ilev, &ncdimid2);

      streamptr->vct.mlev   = mlev;
      streamptr->vct.ilev   = ilev;
      streamptr->vct.mlevID = ncdimid;
      streamptr->vct.ilevID = ncdimid2;

1346
      int hyaiid, hybiid, hyamid, hybmid;
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
      cdf_def_var(fileID, "hyai", NC_DOUBLE, 1, &ncdimid2, &hyaiid);
      cdf_def_var(fileID, "hybi", NC_DOUBLE, 1, &ncdimid2, &hybiid);
      cdf_def_var(fileID, "hyam", NC_DOUBLE, 1, &ncdimid,  &hyamid);
      cdf_def_var(fileID, "hybm", NC_DOUBLE, 1, &ncdimid,  &hybmid);

      {
        static const char lname_n[] = "long_name",
          lname_v_ai[] = "hybrid A coefficient at layer interfaces",
          units_n[] = "units",
          units_v_ai[] = "Pa",
          lname_v_bi[] = "hybrid B coefficient at layer interfaces",
          units_v_bi[] = "1",
          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[]
          = {
          { 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 },
          { 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]) };
        int ids[tabLen] = { hyaiid, hyaiid, hybiid, hybiid,
                            hyamid, hyamid, hybmid, hybmid };
        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);

1389
1390
1391
      size_t start;
      size_t count = 1;
      double mval;
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402