stream_cdf_o.c

#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


static const char bndsName[] = "bnds";


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);
  int memtype  = datatype != CDI_DATATYPE_FLT32 ? MEMTYPE_DOUBLE : MEMTYPE_FLOAT;

  void *data = Malloc((size_t)datasize
             * (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;
  if ( taxis->type == TAXIS_FORECAST && ncvarid != CDI_UNDEFID )
    {
      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] };

  for ( size_t i = 0; i < calTabSize; ++i )
    if ( calTab[i].calCode == calendar )
      {
        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";

  if ( streamptr->basetime.ncvarid != CDI_UNDEFID ) return;

  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;

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

  cdf_def_var(fileID, taxis_name, xtype, 1, &time_dimid, &time_varid);

  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;

      cdf_def_var(fileID, "leadtime", xtype, 1, &time_dimid, &leadtimeid);

      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
void cdfDefComplex(stream_t *streamptr, int gridID, int gridindex)
{
  int dimID = CDI_UNDEFID;
  int fileID  = streamptr->fileID;
  ncgrid_t *ncgrid = streamptr->ncgrid;

  for ( int index = 0; index < gridindex; ++index )
    {
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
        {
          int gridID0 = ncgrid[index].gridID;
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_SPECTRAL || gridtype0 == GRID_FOURIER )
            {
              dimID = ncgrid[index].xdimID;
              break;
            }
        }
    }

  if ( dimID == CDI_UNDEFID )
    {
      static const char axisname[] = "nc2";
      size_t dimlen = 2;

      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
      cdf_def_dim(fileID, axisname, dimlen, &dimID);
      cdf_enddef(fileID);

      streamptr->ncmode = 2;
    }

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
}

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

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

  for ( int index = 0; index < gridindex; index++ )
    {
      if ( ncgrid[index].ydimID != CDI_UNDEFID )
        {
          int gridID0 = ncgrid[index].gridID;
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == gridRefType )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0)/2;
              if ( dimlen == dimlen0 )
                {
                  dimID = ncgrid[index].ydimID;
                  break;
                }
              else
                iz++;
            }
        }
    }

  if ( dimID == CDI_UNDEFID )
    {
      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;
    }

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].ydimID = dimID;
}

static
void cdfDefSP(stream_t *streamptr, int gridID, int gridindex)
{
  /*
  char longname[] = "Spherical harmonic coefficient";
  */
  char axisname[5] = "nspX";
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_SPECTRAL);
}


static
void cdfDefFC(stream_t *streamptr, int gridID, int gridindex)
{
  char axisname[5] = "nfcX";
  cdfDefSPorFC(streamptr, gridID, gridindex, axisname, GRID_FOURIER);
}

static const struct cdfDefGridAxisInqs {
  int (*axisSize)(int gridID);
  int (*axisName)(int cdiID, int key, int size, char *mesg);
  int (*axisLongname)(int cdiID, int key, int size, char *mesg);
  int (*axisUnits)(int cdiID, int key, int size, char *mesg);
  void (*axisStdname)(int cdiID, char *dimstdname);
  double (*axisVal)(int gridID, int index);
  const double *(*axisValsPtr)(int gridID);
  const double *(*axisBoundsPtr)(int gridID);
} gridInqsX = {
  .axisSize = gridInqXsize,
  .axisName = cdiGridInqKeyStr,
  .axisLongname = cdiGridInqKeyStr,
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqXstdname,
  .axisVal = gridInqXval,
  .axisValsPtr = gridInqXvalsPtr,
  .axisBoundsPtr = gridInqXboundsPtr,
}, gridInqsY = {
  .axisSize = gridInqYsize,
  .axisName = cdiGridInqKeyStr,
  .axisLongname = cdiGridInqKeyStr,
  .axisUnits = cdiGridInqKeyStr,
  .axisStdname = gridInqYstdname,
  .axisVal = gridInqYval,
  .axisValsPtr = gridInqYvalsPtr,
  .axisBoundsPtr = gridInqYboundsPtr,
}, gridInqsZ = {
  .axisLongname = cdiZaxisInqKeyStr,
  .axisUnits = cdiZaxisInqKeyStr,
  .axisStdname = zaxisInqStdname,
};

static
void cdfPutGridStdAtts(int fileID, int ncvarid, int gridID, int dimtype, const struct cdfDefGridAxisInqs *inqs)
{
  size_t len;

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

  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)) )
    cdf_put_att_text(fileID, ncvarid, "units", len, units);
}

static void
cdfDefTrajLatLon(stream_t *streamptr, int gridID, int gridindex,
                 const struct cdfDefGridAxisInqs *inqs, int dimtype)
{
  nc_type xtype = (gridInqPrec(gridID) == CDI_DATATYPE_FLT32) ? NC_FLOAT : NC_DOUBLE;
  ncgrid_t *ncgrid = streamptr->ncgrid;

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

  int ncvarid = CDI_UNDEFID;
  if ( dimtype == 'X' )
    ncvarid = ncgrid[gridindex].xdimID;
  else
    ncvarid = ncgrid[gridindex].ydimID;

  if ( ncvarid == CDI_UNDEFID )
    {
      int dimNcID = streamptr->basetime.ncvarid;
      int fileID  = streamptr->fileID;
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      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);
      cdf_def_var(fileID, axisname, xtype, 1, &dimNcID, &ncvarid);
      cdfPutGridStdAtts(fileID, ncvarid, gridID, dimtype, inqs);
      cdf_enddef(fileID);
      streamptr->ncmode = 2;
    }

  ncgrid[gridindex].gridID = gridID;
  if ( dimtype == 'X' )
    ncgrid[gridindex].xdimID = ncvarid; /* var ID for trajectory !!! */
  else
    ncgrid[gridindex].ydimID = ncvarid; /* var ID for trajectory !!! */
}

static
void cdfDefTrajLon(stream_t *streamptr, int gridID, int gridindex)
{
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsX, 'X');
}


static
void cdfDefTrajLat(stream_t *streamptr, int gridID, int gridindex)
{
  cdfDefTrajLatLon(streamptr, gridID, gridindex, &gridInqsY, 'Y');
}

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;
  int dimid = CDI_UNDEFID;
  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
void checkGridName(char *axisname, int fileID)
{
  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);
      if ( status != NC_NOERR ) break;

      ++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
cdfDefAxisCommon(stream_t *streamptr, int gridID, int gridindex, int ndims,
                 const struct cdfDefGridAxisInqs *gridAxisInq, int dimKey, char axisLetter,
                 void (*finishCyclicBounds)(double *pbounds, size_t dimlen, const double *pvals))
{
  int dimID = CDI_UNDEFID;
  int ncvarid = CDI_UNDEFID, ncbvarid = CDI_UNDEFID;
  int nvdimID = CDI_UNDEFID;
  int fileID  = streamptr->fileID;
  size_t dimlen = (size_t)gridAxisInq->axisSize(gridID);
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);

  ncgrid_t *ncgrid = streamptr->ncgrid;

  for ( int index = 0; index < gridindex; ++index )
    {
      int gridID0 = ncgrid[index].gridID;
      assert(gridID0 != CDI_UNDEFID);
      int gridtype0 = gridInqType(gridID0);
      if ( gridtype0 == GRID_GAUSSIAN    ||
           gridtype0 == GRID_LONLAT      ||
           gridtype0 == GRID_PROJECTION  ||
           gridtype0 == GRID_CURVILINEAR ||
           gridtype0 == GRID_GENERIC )
        {
          size_t dimlen0 = (size_t)gridAxisInq->axisSize(gridID0);
          if ( dimlen == dimlen0 )
            {
              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)) )
                {
                  dimID = (dimKey == CDI_KEY_XDIMNAME) ? ncgrid[index].xdimID : ncgrid[index].ydimID;
                  break;
                }
            }
        }
    }

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

      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);
      if ( axisname[0] == 0 ) Error("axis name undefined!");
      size_t axisnameLen = strlen(axisname);

      /* enough to append _ plus up to 100 decimal and trailing \0 */
      char extendedAxisname[axisnameLen + 4 + 1];
      memcpy(extendedAxisname, axisname, axisnameLen + 1);
      checkGridName(extendedAxisname, fileID);
      size_t extendedAxisnameLen = axisnameLen + strlen(extendedAxisname + axisnameLen);

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

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

          if ( pvals == NULL )
            cdiGridInqKeyStr(gridID, dimKey, CDI_MAX_NAME, dimname);

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

          if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);
        }

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

          cdfPutGridStdAtts(fileID, ncvarid, gridID, axisLetter, gridAxisInq);
          {
            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);
            }
          if ( pbounds && nvdimID != CDI_UNDEFID )
            {
              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;

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

      if ( ndims == 0 )
        {
          if ( dimKey == CDI_KEY_XDIMNAME )
            ncgrid[gridindex].xvarID = ncvarid;
          else
            ncgrid[gridindex].yvarID = ncvarid;
        }
    }

  ncgrid[gridindex].gridID = gridID;
  if ( dimKey == CDI_KEY_XDIMNAME )
    ncgrid[gridindex].xdimID = dimID;
  else
    ncgrid[gridindex].ydimID = dimID;
}

static
void finishCyclicXBounds(double *pbounds, size_t dimlen, const double *pvals)
{
  pbounds[0] = (pvals[0] + pvals[dimlen-1]-360)*0.5;
  pbounds[2*dimlen-1] = (pvals[dimlen-1] + pvals[0]+360)*0.5;
}

static
void cdfDefXaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
{
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsX,
                   CDI_KEY_XDIMNAME, 'X', finishCyclicXBounds);
}

static
void finishCyclicYBounds(double *pbounds, size_t dimlen, const double *pvals)
{
  pbounds[0] = copysign(90.0, pvals[0]);
  pbounds[2*dimlen-1] = copysign(90.0, pvals[dimlen-1]);
}

static
void cdfDefYaxis(stream_t *streamptr, int gridID, int gridindex, int ndims)
{
  cdfDefAxisCommon(streamptr, gridID, gridindex, ndims, &gridInqsY,
                   CDI_KEY_YDIMNAME, 'Y', finishCyclicYBounds);
}

static
void cdfGridCompress(int fileID, int ncvarid, int gridsize, int filetype, int comptype)
{
#if  defined  (HAVE_NETCDF4)
  if ( gridsize > 1 && comptype == CDI_COMPRESS_ZIP && (filetype == CDI_FILETYPE_NC4 || filetype == CDI_FILETYPE_NC4C) )
    {
      nc_def_var_chunking(fileID, ncvarid, NC_CHUNKED, NULL);
      cdfDefVarDeflate(fileID, ncvarid, 1);
    }
#endif
}

static
void cdfDefCurvilinear(stream_t *streamptr, int gridID, int gridindex)
{
  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;
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);
  ncgrid_t *ncgrid = streamptr->ncgrid;

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

  for ( int index = 0; index < gridindex; index++ )
    {
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
        {
          int gridID0 = ncgrid[index].gridID;
          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)) )
                  {
                    xdimID = ncgrid[index].xdimID;
                    ydimID = ncgrid[index].ydimID;
                    ncxvarid = ncgrid[index].xvarID;
                    ncyvarid = ncgrid[index].yvarID;
                    break;
                  }
            }
        }
    }

  if ( xdimID == CDI_UNDEFID || ydimID == CDI_UNDEFID )
    {
      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 ) { xdimname[0] = 'x'; xdimname[1] = 0; }
        xdimID = checkDimName(fileID, xdimlen, xdimname);
        if ( xdimID == CDI_UNDEFID ) cdf_def_dim(fileID, xdimname, xdimlen, &xdimID);
      }
      {
        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];
      if ( gridInqXboundsPtr(gridID) || gridInqYboundsPtr(gridID) )
        {
          char vdimname[CDI_MAX_NAME+3]; vdimname[0] = 0;
          cdiGridInqKeyStr(gridID, CDI_KEY_VDIMNAME, CDI_MAX_NAME, vdimname);
          if ( vdimname[0] == 0 ) strcpy(vdimname, "nv4");
          size_t nvertex = 4;
          nvdimID = checkDimName(fileID, nvertex, vdimname);
          if ( nvdimID == CDI_UNDEFID ) cdf_def_dim(fileID, vdimname, nvertex, &nvdimID);
        }

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

      if ( gridInqXvalsPtr(gridID) )
        {
          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, 2, dimIDs, &ncxvarid);
          cdfGridCompress(fileID, ncxvarid, (int)(xdimlen*ydimlen), streamptr->filetype, streamptr->comptype);

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

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

          if ( gridInqXboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
            {
              size_t xaxisnameLen = strlen(xaxisname);
              xaxisname[xaxisnameLen] = '_';
              memcpy(xaxisname + xaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, xaxisname, xtype, 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);
          checkGridName(yaxisname, fileID);

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

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

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

          if ( gridInqYboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
            {
              size_t yaxisnameLen = strlen(yaxisname);
              yaxisname[yaxisnameLen] = '_';
              memcpy(yaxisname + yaxisnameLen + 1, bndsName, sizeof (bndsName));
              cdf_def_var(fileID, yaxisname, xtype, 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;

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

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = xdimID;
  ncgrid[gridindex].ydimID = ydimID;
  ncgrid[gridindex].xvarID = ncxvarid;
  ncgrid[gridindex].yvarID = ncyvarid;
  ncgrid[gridindex].avarID = ncavarid;
}

static
void cdfDefRgrid(stream_t *streamptr, int gridID, int gridindex)
{
  ncgrid_t *ncgrid = streamptr->ncgrid;
  int dimID = CDI_UNDEFID;

  size_t dimlen = (size_t)gridInqSize(gridID);

  int iz = 0;
  for ( int index = 0; index < gridindex; index++ )
    {
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
        {
          int gridID0 = ncgrid[index].gridID;
          int gridtype0 = gridInqType(gridID0);
          if ( gridtype0 == GRID_GAUSSIAN_REDUCED )
            {
              size_t dimlen0 = (size_t)gridInqSize(gridID0);

              if ( dimlen == dimlen0 )
                {
                  dimID = ncgrid[index].xdimID;
                  break;
                }
              iz++;
            }
        }
    }

  if ( dimID == CDI_UNDEFID )
    {
      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;
    }

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
}

static
void cdfDefGdim(stream_t *streamptr, int gridID, int gridindex)
{
  ncgrid_t *ncgrid = streamptr->ncgrid;
  int iz = 0;
  int dimID = CDI_UNDEFID;

  size_t dimlen = (size_t)gridInqSize(gridID);

  if ( gridInqYsize(gridID) == 0 )
    for ( int index = 0; index < gridindex; index++ )
      {
        if ( ncgrid[index].xdimID != CDI_UNDEFID )
          {
            int gridID0 = ncgrid[index].gridID;
            int gridtype0 = gridInqType(gridID0);
            if ( gridtype0 == GRID_GENERIC )
              {
                size_t dimlen0 = (size_t)gridInqSize(gridID0);
                if ( dimlen == dimlen0 )
                  {
                    dimID = ncgrid[index].xdimID;
                    break;
                  }
                else
                  iz++;
              }
          }
      }

  if ( gridInqXsize(gridID) == 0 )
    for ( int index = 0; index < gridindex; index++ )
      {
        if ( ncgrid[index].ydimID != CDI_UNDEFID )
          {
            int gridID0 = ncgrid[index].gridID;
            int gridtype0 = gridInqType(gridID0);
            if ( gridtype0 == GRID_GENERIC )
              {
                size_t dimlen0 = (size_t)gridInqSize(gridID0);
                if ( dimlen == dimlen0 )
                  {
                    dimID = ncgrid[index].ydimID;
                    break;
                  }
                else
                  iz++;
              }
          }
      }

  if ( dimID == CDI_UNDEFID )
    {
      int fileID  = streamptr->fileID;
      char dimname[CDI_MAX_NAME];
      strcpy(dimname, "gsize");

      dimID = checkDimName(fileID, dimlen, dimname);

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

      if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, dimname, dimlen, &dimID);

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

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
}

static
void cdfDefGridReference(stream_t *streamptr, int gridID)
{
  int fileID  = streamptr->fileID;
  int number = gridInqNumber(gridID);
  if ( number > 0 )
    cdf_put_att_int(fileID, NC_GLOBAL, "number_of_grid_used", NC_INT, 1, &number);

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

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

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

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
void cdfDefUnstructured(stream_t *streamptr, int gridID, int gridindex)
{
  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;
  nc_type xtype = (nc_type)cdfDefDatatype(gridInqPrec(gridID), streamptr->filetype);
  ncgrid_t *ncgrid = streamptr->ncgrid;

  int fileID  = streamptr->fileID;

  size_t dimlen = (size_t)gridInqSize(gridID);

  for ( int index = 0; index < gridindex; index++ )
    {
      if ( ncgrid[index].xdimID != CDI_UNDEFID )
        {
          int gridID0 = ncgrid[index].gridID;
          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)) )
		  {
		    dimID = ncgrid[index].xdimID;
                    ncxvarid = ncgrid[index].xvarID;
                    ncyvarid = ncgrid[index].yvarID;
                    ncavarid = ncgrid[index].avarID;
		    break;
		  }
            }
        }
    }

  if ( dimID == CDI_UNDEFID )
    {
      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, "ncells");
        dimID = checkDimName(fileID, dimlen, xdimname);
        if ( dimID == CDI_UNDEFID ) cdf_def_dim(fileID, xdimname, dimlen, &dimID);
      }

      size_t nvertex = (size_t)gridInqNvertex(gridID);
      if ( nvertex > 0 )
        {
          char vdimname[CDI_MAX_NAME+3];
          vdimname[0] = 0;
          cdiGridInqKeyStr(gridID, CDI_KEY_VDIMNAME, CDI_MAX_NAME, vdimname);
          if ( vdimname[0] == 0 ) strcpy(vdimname, "vertices");
          nvdimID = checkDimName(fileID, nvertex, vdimname);
          if ( nvdimID == CDI_UNDEFID ) cdf_def_dim(fileID, vdimname, nvertex, &nvdimID);
        }

      cdfDefGridReference(streamptr, gridID);

      cdfDefGridUUID(streamptr, gridID);

      if ( gridInqXvalsPtr(gridID) )
        {
          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, 1, &dimID, &ncxvarid);
          cdfGridCompress(fileID, ncxvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

          if ( gridInqXboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
            {
              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);
          checkGridName(yaxisname, fileID);
          cdf_def_var(fileID, yaxisname, xtype, 1, &dimID, &ncyvarid);
          cdfGridCompress(fileID, ncyvarid, (int)dimlen, streamptr->filetype, streamptr->comptype);

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

          if ( gridInqYboundsPtr(gridID) && nvdimID != CDI_UNDEFID )
            {
              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;

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

  ncgrid[gridindex].gridID = gridID;
  ncgrid[gridindex].xdimID = dimID;
  ncgrid[gridindex].xvarID = ncxvarid;
  ncgrid[gridindex].yvarID = ncyvarid;
  ncgrid[gridindex].avarID = ncavarid;
}

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 )
            Error("More than one VCT for each file unsupported!");
          return;
        }

      int fileID = streamptr->fileID;

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

      int ncdimid, ncdimid2;
      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;

      int hyaiid, hybiid, hyamid, hybmid;
      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);

      size_t start;
      size_t count = 1;
      double mval;