stream_cdf.c

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

//#define TEST_GROUPS 1

#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <float.h>

#if  defined  (HAVE_LIBNETCDF)
#  include <netcdf.h>
#endif


#include "dmemory.h"

#include "cdi.h"
#include "basetime.h"
#include "gaussgrid.h"
#include "cdi_int.h"
#include "stream_cdf.h"
#include "cdf_int.h"
#include "varscan.h"
#include "vlist.h"

//#define PROJECTION_TEST

#undef  UNDEFID
#define UNDEFID  CDI_UNDEFID


void    vlistInqVarDimorder(int vlistID, int varID, int (*outDimorder)[3]);

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


#define  X_AXIS  1
#define  Y_AXIS  2
#define  Z_AXIS  3
#define  T_AXIS  4

#define  POSITIVE_UP    1
#define  POSITIVE_DOWN  2

typedef struct {
  int     ncvarid;
  int     dimtype;
  size_t  len;
  char    name[CDI_MAX_NAME];
}
ncdim_t;

#define  MAX_COORDVARS  4
#define  MAX_AUXVARS    4

typedef struct {
  int      ncid;
  int      ignore;
  int      isvar;
  int      islon;
  int      islat;
  int      islev;
  int      istime;
  int      warn;
  int      tsteptype;
  int      param;
  int      code;
  int      tabnum;
  int      climatology;
  int      bounds;
  int      gridID;
  int      zaxisID;
  int      gridtype;
  int      zaxistype;
  int      xdim;
  int      ydim;
  int      zdim;
  int      xvarid;
  int      yvarid;
  int      zvarid;
  int      tvarid;
  int      ncoordvars;
  int      coordvarids[MAX_COORDVARS];
  int      nauxvars;
  int      auxvarids[MAX_AUXVARS];
  int      cellarea;
  int      calendar;
  int      tableID;
  int      truncation;
  int      position;
  int      defmissval;
  int      deffillval;
  int      xtype;
  int      ndims;
  int      gmapid;
  int      positive;
  int      dimids[8];
  int      dimtype[8];
  int      chunks[8];
  int      chunked;
  int      chunktype;
  int      natts;
  int     *atts;
  int      deflate;
  int      lunsigned;
  int      lvalidrange;
  size_t   vlen;
  double  *vdata;
  double   missval;
  double   fillval;
  double   addoffset;
  double   scalefactor;
  double   validrange[2];
  char     name[CDI_MAX_NAME];
  char     longname[CDI_MAX_NAME];
  char     stdname[CDI_MAX_NAME];
  char     units[CDI_MAX_NAME];
  char     extra[CDI_MAX_NAME];
  ensinfo_t   *ensdata;    /* Ensemble information */
}
ncvar_t;

#ifdef HAVE_LIBNETCDF
static
void strtolower(char *str)
{
  if ( str )
    for (size_t i = 0; str[i]; ++i)
      str[i] = (char)tolower((int)str[i]);
}

static
int get_timeunit(size_t len, const char *ptu)
{
  int timeunit = -1;

  if ( len > 2 )
    {
      if      ( memcmp(ptu, "sec",    3) == 0 )          timeunit = TUNIT_SECOND;
      else if ( memcmp(ptu, "minute", 6) == 0 )          timeunit = TUNIT_MINUTE;
      else if ( memcmp(ptu, "hour",   4) == 0 )          timeunit = TUNIT_HOUR;
      else if ( memcmp(ptu, "day",    3) == 0 )          timeunit = TUNIT_DAY;
      else if ( memcmp(ptu, "month",  5) == 0 )          timeunit = TUNIT_MONTH;
      else if ( memcmp(ptu, "calendar_month", 14) == 0 ) timeunit = TUNIT_MONTH;
      else if ( memcmp(ptu, "year",   4) == 0 )          timeunit = TUNIT_YEAR;
    }
  else if ( len == 1 )
    {
      if ( ptu[0] == 's' ) timeunit = TUNIT_SECOND;
    }

  return (timeunit);
}

static
int isTimeUnits(const char *timeunits)
{
  int status = 0;

  if ( strncmp(timeunits, "sec",    3) == 0 ||
       strncmp(timeunits, "minute", 6) == 0 ||
       strncmp(timeunits, "hour",   4) == 0 ||
       strncmp(timeunits, "day",    3) == 0 ||
       strncmp(timeunits, "month",  5) == 0 ) status = 1;

  return (status);
}

static
int isTimeAxisUnits(const char *timeunits)
{
  char *ptu, *tu;
  int timetype = -1;
  int timeunit;
  int status = FALSE;

  size_t len = strlen(timeunits);
  tu = xmalloc((len+1)*sizeof(char));
  memcpy(tu, timeunits, (len+1) * sizeof(char));
  ptu = tu;

  for (size_t i = 0; i < len; i++ ) ptu[i] = (char)tolower((int)ptu[i]);

  timeunit = get_timeunit(len, ptu);
  if ( timeunit != -1 )
    {

      while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
      if ( *ptu )
        {
          while ( isspace(*ptu) ) ptu++;

          if ( memcmp(ptu, "as", 2) == 0 )
            timetype = TAXIS_ABSOLUTE;
          else if ( memcmp(ptu, "since", 5) == 0 )
            timetype = TAXIS_RELATIVE;

          if ( timetype != -1 ) status = TRUE;
        }
    }

  free(tu);

  return (status);
}

static
void scanTimeString(const char *ptu, int *rdate, int *rtime)
{
  int year, month, day;
  int hour = 0, minute = 0, second = 0;
  int v1, v2, v3;

  *rdate = 0;
  *rtime = 0;

  v1 = atoi(ptu);
  if ( v1 < 0 ) ptu++;
  while ( isdigit((int) *ptu) ) ptu++;
  v2 = atoi(++ptu);
  while ( isdigit((int) *ptu) ) ptu++;
  v3 = atoi(++ptu);
  while ( isdigit((int) *ptu) ) ptu++;

  if ( v3 > 999 && v1 < 32 )
    { year = v3; month = v2; day = v1; }
  else
    { year = v1; month = v2; day = v3; }

  while ( isspace((int) *ptu) ) ptu++;

  if ( *ptu )
    {
      while ( ! isdigit((int) *ptu) ) ptu++;

      hour = atoi(ptu);
      while ( isdigit((int) *ptu) ) ptu++;
      if ( *ptu == ':' )
        {
          ptu++;
          minute = atoi(ptu);
          while ( isdigit((int) *ptu) ) ptu++;
          if ( *ptu == ':' )
            {
              ptu++;
              second = atoi(ptu);
            }
        }
    }

  *rdate = cdiEncodeDate(year, month, day);
  *rtime = cdiEncodeTime(hour, minute, second);
}

static
int scanTimeUnit(const char *unitstr)
{
  int timeunit = -1;
  int len;

  len = (int) strlen(unitstr);
  timeunit = get_timeunit(len, unitstr);
  if ( timeunit == -1 )
    Message("Unsupported TIMEUNIT: %s!", unitstr);

  return (timeunit);
}

static
void setForecastTime(const char *timestr, taxis_t *taxis)
{
  int len;

  (*taxis).fdate = 0;
  (*taxis).ftime = 0;

  len = (int) strlen(timestr);
  if ( len == 0 ) return;

  int fdate = 0, ftime = 0;
  scanTimeString(timestr, &fdate, &ftime);

  (*taxis).fdate = fdate;
  (*taxis).ftime = ftime;
}

static
int setBaseTime(const char *timeunits, taxis_t *taxis)
{
  char *ptu, *tu;
  int timetype = TAXIS_ABSOLUTE;
  int rdate = -1, rtime = -1;
  int timeunit;

  size_t len = strlen(timeunits);
  tu = xmalloc((len+1) * sizeof (char));
  memcpy(tu, timeunits, (len+1) * sizeof (char));
  ptu = tu;

  for ( size_t i = 0; i < len; i++ ) ptu[i] = (char)tolower((int) ptu[i]);

  timeunit = get_timeunit(len, ptu);
  if ( timeunit == -1 )
    {
      Message("Unsupported TIMEUNIT: %s!", timeunits);
      return (1);
    }

  while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
  if ( *ptu )
    {
      while ( isspace(*ptu) ) ptu++;

      if ( memcmp(ptu, "as", 2) == 0 )
        timetype = TAXIS_ABSOLUTE;
      else if ( memcmp(ptu, "since", 5) == 0 )
        timetype = TAXIS_RELATIVE;

      while ( ! isspace(*ptu) && *ptu != 0 ) ptu++;
      if ( *ptu )
        {
          while ( isspace(*ptu) ) ptu++;

          if ( timetype == TAXIS_ABSOLUTE )
            {
              if ( memcmp(ptu, "%y%m%d.%f", 9) != 0 && timeunit == TUNIT_DAY )
                {
                  Message("Unsupported format %s for TIMEUNIT day!", ptu);
                  timeunit = -1;
                }
              else if ( memcmp(ptu, "%y%m.%f", 7) != 0 && timeunit == TUNIT_MONTH )
                {
                  Message("Unsupported format %s for TIMEUNIT month!", ptu);
                  timeunit = -1;
                }
            }
          else if ( timetype == TAXIS_RELATIVE )
            {
              scanTimeString(ptu, &rdate, &rtime);

              (*taxis).rdate = rdate;
              (*taxis).rtime = rtime;

              if ( CDI_Debug )
                Message("rdate = %d  rtime = %d", rdate, rtime);
            }
        }
    }

  (*taxis).type = timetype;
  (*taxis).unit = timeunit;

  free(tu);

  if ( CDI_Debug )
    Message("timetype = %d  unit = %d", timetype, timeunit);

  return (0);
}

static
void cdfGetAttInt(int fileID, int ncvarid, char *attname, int attlen, int *attint)
{
  size_t nc_attlen;
  int *pintatt;

  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);

  if ( (int)nc_attlen > attlen )
    pintatt = (int *) malloc(nc_attlen * sizeof (int));
  else
    pintatt = attint;

  cdf_get_att_int(fileID, ncvarid, attname, pintatt);

  if ( (int)nc_attlen > attlen )
    {
      memcpy(attint, pintatt, (size_t)attlen * sizeof (int));
      free(pintatt);
    }
}

static
void cdfGetAttDouble(int fileID, int ncvarid, char *attname, int attlen, double *attdouble)
{
  size_t nc_attlen;
  double *pdoubleatt;

  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);

  if ( (int)nc_attlen > attlen )
    pdoubleatt = (double *) malloc(nc_attlen * sizeof (double));
  else
    pdoubleatt = attdouble;

  cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);

  if ( (int)nc_attlen > attlen )
    {
      memcpy(attdouble, pdoubleatt, (size_t)attlen * sizeof (double));
      free(pdoubleatt);
    }
}

static
void cdfGetAttText(int fileID, int ncvarid, char *attname, int attlen, char *atttext)
{
  nc_type atttype;
  size_t nc_attlen;

  cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
  cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);

  if ( atttype == NC_CHAR )
    {
      char attbuf[65636];
      if ( nc_attlen < sizeof(attbuf) )
        {
          cdf_get_att_text(fileID, ncvarid, attname, attbuf);

          attbuf[nc_attlen++] = 0;

          if ( (int) nc_attlen > attlen ) nc_attlen = (size_t)attlen;
          memcpy(atttext, attbuf, nc_attlen);
        }
      else
        {
          atttext[0] = 0;
        }
    }
#if  defined  (HAVE_NETCDF4)
  else if ( atttype == NC_STRING )
    {
      if ( nc_attlen == 1 )
        {
          char *attbuf = NULL;
          cdf_get_att_string(fileID, ncvarid, attname, &attbuf);

          int slen = (int) strlen(attbuf);
          slen++;

          if ( slen > attlen ) slen = attlen;
          memcpy(atttext, attbuf, slen);

          free(attbuf);
        }
      else
        {
          atttext[0] = 0;
        }
    }
#endif
}

static
int xtypeIsText(int xtype)
{
  int isText = FALSE;

  if ( xtype == NC_CHAR )
    isText = TRUE;
#if  defined  (HAVE_NETCDF4)
  else if ( xtype == NC_STRING )
    isText = TRUE;
#endif

  return (isText);
}

static
int xtypeIsFloat(int xtype)
{
  int isFloat = FALSE;

  if ( xtype == NC_FLOAT || xtype == NC_DOUBLE ) isFloat = TRUE;

  return isFloat;
}

static
int cdfInqDatatype(int xtype, int lunsigned)
{
  int datatype = -1;

#if  defined  (HAVE_NETCDF4)
  if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
#endif

  if      ( xtype == NC_BYTE   )  datatype = DATATYPE_INT8;
  /* else if ( xtype == NC_CHAR   )  datatype = DATATYPE_UINT8; */
  else if ( xtype == NC_SHORT  )  datatype = DATATYPE_INT16;
  else if ( xtype == NC_INT    )  datatype = DATATYPE_INT32;
  else if ( xtype == NC_FLOAT  )  datatype = DATATYPE_FLT32;
  else if ( xtype == NC_DOUBLE )  datatype = DATATYPE_FLT64;
#if  defined  (HAVE_NETCDF4)
  else if ( xtype == NC_UBYTE  )  datatype = DATATYPE_UINT8;
  else if ( xtype == NC_LONG   )  datatype = DATATYPE_INT32;
  else if ( xtype == NC_USHORT )  datatype = DATATYPE_UINT16;
  else if ( xtype == NC_UINT   )  datatype = DATATYPE_UINT32;
  else if ( xtype == NC_INT64  )  datatype = DATATYPE_FLT64;
  else if ( xtype == NC_UINT64 )  datatype = DATATYPE_FLT64;
#endif

  return (datatype);
}

static
int cdfDefDatatype(int datatype, int filetype)
{
  int xtype;

  if ( datatype == DATATYPE_CPX32 || datatype == DATATYPE_CPX64 )
    Error("CDI/netCDF library does not support complex numbers!");

  if ( filetype == FILETYPE_NC4 )
    {
      if      ( datatype == DATATYPE_INT8   ) xtype = NC_BYTE;
      else if ( datatype == DATATYPE_INT16  ) xtype = NC_SHORT;
      else if ( datatype == DATATYPE_INT32  ) xtype = NC_INT;
#if  defined  (HAVE_NETCDF4)
      else if ( datatype == DATATYPE_UINT8  ) xtype = NC_UBYTE;
      else if ( datatype == DATATYPE_UINT16 ) xtype = NC_USHORT;
      else if ( datatype == DATATYPE_UINT32 ) xtype = NC_UINT;
#else
      else if ( datatype == DATATYPE_UINT8  ) xtype = NC_SHORT;
      else if ( datatype == DATATYPE_UINT16 ) xtype = NC_INT;
      else if ( datatype == DATATYPE_UINT32 ) xtype = NC_INT;
#endif
      else if ( datatype == DATATYPE_FLT64  ) xtype = NC_DOUBLE;
      else                                    xtype = NC_FLOAT;
    }
  else
    {
      if      ( datatype == DATATYPE_INT8   ) xtype = NC_BYTE;
      else if ( datatype == DATATYPE_INT16  ) xtype = NC_SHORT;
      else if ( datatype == DATATYPE_INT32  ) xtype = NC_INT;
      else if ( datatype == DATATYPE_UINT8  ) xtype = NC_SHORT;
      else if ( datatype == DATATYPE_UINT16 ) xtype = NC_INT;
      else if ( datatype == DATATYPE_UINT32 ) xtype = NC_INT;
      else if ( datatype == DATATYPE_FLT64  ) xtype = NC_DOUBLE;
      else                                    xtype = NC_FLOAT;
    }

  return (xtype);
}

static
void defineAttributes(int vlistID, int varID, int fileID, int ncvarID)
{
  int natts, iatt;
  int atttype, attlen;
  size_t len;
  char attname[1024];

  vlistInqNatts(vlistID, varID, &natts);

  for ( iatt = 0; iatt < natts; iatt++ )
    {
      vlistInqAtt(vlistID, varID, iatt, attname, &atttype, &attlen);

      if ( attlen == 0 ) continue;

      if ( atttype == DATATYPE_TXT )
        {
          char *atttxt = (char *)xmalloc((size_t)attlen*sizeof(char));
          vlistInqAttTxt(vlistID, varID, attname, attlen, atttxt);
          len = (size_t)attlen;
          cdf_put_att_text(fileID, ncvarID, attname, len, atttxt);
          free(atttxt);
        }
      else if ( atttype == DATATYPE_INT16 || atttype == DATATYPE_INT32 )
        {
          int *attint = (int *)xmalloc((size_t)attlen*sizeof(int));
          vlistInqAttInt(vlistID, varID, attname, attlen, &attint[0]);
          len = (size_t)attlen;
          if ( atttype == DATATYPE_INT16 )
            cdf_put_att_int(fileID, ncvarID, attname, NC_SHORT, len, attint);
          else
            cdf_put_att_int(fileID, ncvarID, attname, NC_INT, len, attint);
          free(attint);
        }
      else if ( atttype == DATATYPE_FLT32 || atttype == DATATYPE_FLT64 )
        {
          double *attflt = (double *)xmalloc((size_t)attlen*sizeof(double));
          vlistInqAttFlt(vlistID, varID, attname, attlen, attflt);
          len = (size_t)attlen;
          if ( atttype == DATATYPE_FLT32 )
            cdf_put_att_double(fileID, ncvarID, attname, NC_FLOAT, len, attflt);
          else
            cdf_put_att_double(fileID, ncvarID, attname, NC_DOUBLE, len, attflt);
          free(attflt);
        }
    }
}
#endif

void cdfCopyRecord(stream_t *streamptr2, stream_t *streamptr1)
{
  int memtype = MEMTYPE_DOUBLE;
  int vlistID1 = streamptr1->vlistID;

  int tsID1 = streamptr1->curTsID;

  int recID1 = streamptr1->tsteps[tsID1].curRecID;

  int ivarID = streamptr1->tsteps[tsID1].records[recID1].varID;

  int gridID = vlistInqVarGrid(vlistID1, ivarID);

  int datasize = gridInqSize(gridID);
  /* bug fix for constant netCDF fields */
  if ( datasize < 1048576 ) datasize = 1048576;

  double *data = xmalloc((size_t)datasize * sizeof (double));

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

  free(data);
}

/* not used
int cdfInqRecord(stream_t *streamptr, int *varID, int *levelID)
{
  int tsID, recID;

  recID = streamptr->tsteps[0].curRecID++;
  printf("cdfInqRecord recID %d %d\n", recID, streamptr->tsteps[0].curRecID);
  printf("cdfInqRecord tsID %d\n", streamptr->curTsID);

  if ( streamptr->tsteps[0].curRecID >= streamptr->tsteps[0].nrecs )
    {
      streamptr->tsteps[0].curRecID = 0;
    }

  *varID   = streamptr->tsteps[0].records[recID].varID;
  *levelID = streamptr->tsteps[0].records[recID].levelID;

  streamptr->record->varID   = *varID;
  streamptr->record->levelID = *levelID;

  if ( CDI_Debug )
    Message("recID = %d  varID = %d  levelID = %d", recID, *varID, *levelID);

  return (recID+1);
}
*/


void cdfDefRecord(stream_t *streamptr)
{
}

#if  defined  (HAVE_LIBNETCDF)
static
void cdfWriteGridTraj(stream_t *streamptr, int gridID)
{
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;

  int gridindex = vlistGridIndex(vlistID, gridID);
  int lonID = streamptr->xdimID[gridindex],
    latID = streamptr->ydimID[gridindex];

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

  cdf_put_var1_double(fileID, lonID, &index, &xlon);
  cdf_put_var1_double(fileID, latID, &index, &xlat);
}

static
void cdfReadGridTraj(stream_t *streamptr, int gridID)
{
  int vlistID = streamptr->vlistID;
  int fileID  = streamptr->fileID;

  int gridindex = vlistGridIndex(vlistID, gridID);
  int lonID = streamptr->xdimID[gridindex];
  int latID = streamptr->ydimID[gridindex];

  int tsID = streamptr->curTsID;
  size_t index = (size_t)tsID;

  double xlon, xlat;
  cdf_get_var1_double(fileID, lonID, &index, &xlon);
  cdf_get_var1_double(fileID, latID, &index, &xlat);

  gridDefXvals(gridID, &xlon);
  gridDefYvals(gridID, &xlat);
}
#endif

#if  defined  (HAVE_LIBNETCDF)
static
void cdfDefVarDeflate(int ncid, int ncvarid, int deflate_level)
{
#if  defined  (HAVE_NETCDF4)
  int retval;
  /* Set chunking, shuffle, and deflate. */
  int shuffle = 1;
  int deflate = 1;

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

  if ((retval = nc_def_var_deflate(ncid, ncvarid, shuffle, deflate, deflate_level)))
    {
      Error("nc_def_var_deflate failed, status = %d", retval);
    }
#else
  static int lwarn = TRUE;

  if ( lwarn )
    {
      lwarn = FALSE;
      Warning("Deflate compression failed, netCDF4 not available!");
    }
#endif
}
#endif

#if  defined(HAVE_LIBNETCDF) && defined(NC_SZIP_NN_OPTION_MASK)
static
void cdfDefVarSzip(int ncid, int ncvarid)
{
  int retval;
  /* Set options_mask and bits_per_pixel. */
  int options_mask = NC_SZIP_NN_OPTION_MASK;
  int bits_per_pixel = 16;

  if ((retval = nc_def_var_szip(ncid, ncvarid, options_mask, bits_per_pixel)))
    {
      if ( retval == NC_EINVAL )
        {
          static int lwarn = TRUE;

          if ( lwarn )
            {
              lwarn = FALSE;
              Warning("netCDF4/Szip compression not compiled in!");
            }
        }
      else
        Error("nc_def_var_szip failed, status = %d", retval);
    }
}
#endif

#if  defined  (HAVE_LIBNETCDF)
static
void cdfDefVarMissval(stream_t *streamptr, int varID, int dtype, int lcheck)
{
  if ( streamptr->vars[varID].defmiss == FALSE )
    {
      int fileID;
      int ncvarid;
      double missval;
      int vlistID;
      int xtype;

      vlistID = streamptr->vlistID;
      fileID  = streamptr->fileID;
      ncvarid = streamptr->vars[varID].ncvarid;
      missval = vlistInqVarMissval(vlistID, varID);

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

      xtype = cdfDefDatatype(dtype, streamptr->filetype);

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

      cdf_put_att_double(fileID, ncvarid, "_FillValue", (nc_type) xtype, 1, &missval);
      cdf_put_att_double(fileID, ncvarid, "missing_value", (nc_type) xtype, 1, &missval);

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

      streamptr->vars[varID].defmiss = TRUE;
    }
}
#endif

void cdf_write_record(stream_t *streamptr, int memtype, const void *data, int nmiss)
{
#if  defined  (HAVE_LIBNETCDF)
  int varID;
  int levelID;

  varID   = streamptr->record->varID;
  levelID = streamptr->record->levelID;

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

  cdf_write_var_slice(streamptr, varID, levelID, memtype, data, nmiss);
#endif
}

void cdfReadRecord(stream_t *streamptr, double *data, int *nmiss)
{
  int levelID, varID, tsID, recID, vrecID;

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

  tsID    = streamptr->curTsID;
  vrecID  = streamptr->tsteps[tsID].curRecID;
  recID   = streamptr->tsteps[tsID].recIDs[vrecID];
  varID   = streamptr->tsteps[tsID].records[recID].varID;
  levelID = streamptr->tsteps[tsID].records[recID].levelID;

  cdfReadVarSliceDP(streamptr, varID, levelID, data, nmiss);
}

#if  defined  (HAVE_LIBNETCDF)
static
void cdfDefTimeValue(stream_t *streamptr, int tsID)
{
  int fileID;
  int ncvarid;
  taxis_t *taxis;

  fileID = streamptr->fileID;

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

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

  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 != 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, char* taxis_name, taxis_t* taxis)
{
  int time_bndsid = -1;
  int dims[2];
  char tmpstr[CDI_MAX_NAME];

  dims[0] = nctimedimid;

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

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

  if ( taxis->climatology )
    {
      strcpy(tmpstr, "climatology");
      strcat(tmpstr, "_bnds");
      cdf_def_var(fileID, tmpstr, NC_DOUBLE, 2, dims, &time_bndsid);

      cdf_put_att_text(fileID, nctimevarid, "climatology", strlen(tmpstr), tmpstr);
    }
  else
    {
      strcpy(tmpstr, taxis_name);
      strcat(tmpstr, "_bnds");
      cdf_def_var(fileID, tmpstr, NC_DOUBLE, 2, dims, &time_bndsid);

      cdf_put_att_text(fileID, nctimevarid, "bounds", strlen(tmpstr), tmpstr);
    }

  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 year, month, day, hour, minute, second;
      int rdate, rtime;
      int timeunit;

      timeunit = taxis->unit;
      if ( timeunit == -1 ) timeunit = TUNIT_HOUR;
      rdate    = taxis->rdate;
      rtime    = taxis->rtime;
      if ( rdate == -1 )
        {
          rdate  = taxis->vdate;
          rtime  = taxis->vtime;
        }

      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-%02d-%02d %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;

  sprintf(unitstr, "%s", tunitNamePtr(timeunit));
}
#endif

static
void cdfDefCalendar(int fileID, int ncvarid, int calendar)
{
  size_t len;
  char calstr[80];

  calstr[0] = 0;

  if      ( calendar == CALENDAR_STANDARD )  strcpy(calstr, "standard");
  else if ( calendar == CALENDAR_PROLEPTIC ) strcpy(calstr, "proleptic_gregorian");
  else if ( calendar == CALENDAR_NONE )      strcpy(calstr, "none");
  else if ( calendar == CALENDAR_360DAYS )   strcpy(calstr, "360_day");
  else if ( calendar == CALENDAR_365DAYS )   strcpy(calstr, "365_day");
  else if ( calendar == CALENDAR_366DAYS )   strcpy(calstr, "366_day");

  len = strlen(calstr);

#if  defined  (HAVE_LIBNETCDF)
  if ( len ) cdf_put_att_text(fileID, ncvarid, "calendar", len, calstr);
#endif
}

static
void cdfDefTime(stream_t* streamptr)
{
#if  defined  (HAVE_LIBNETCDF)
  int fileID;
  int time_varid;
  int time_dimid;
  int time_bndsid = -1;
  char unitstr[CDI_MAX_NAME];
  char tmpstr[CDI_MAX_NAME];
  char default_name[] = "time";
  char* taxis_name = default_name;
  size_t len;
  taxis_t* taxis;

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

  fileID = streamptr->fileID;

  if ( streamptr->ncmode == 0 ) streamptr->ncmode = 1;

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

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

  if ( taxis->name && taxis->name[0] ) taxis_name = taxis->name;

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

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

  streamptr->basetime.ncvarid = time_varid;

  strcpy(tmpstr, "time");
  cdf_put_att_text(fileID, time_varid, "standard_name", strlen(tmpstr), tmpstr);

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

  cdfDefTimeUnits(unitstr, &streamptr->tsteps[0].taxis, taxis);

  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", NC_DOUBLE, 1, &time_dimid, &leadtimeid);

      streamptr->basetime.leadtimeid = leadtimeid;

      strcpy(tmpstr, "forecast_period");
      cdf_put_att_text(fileID, leadtimeid, "standard_name", strlen(tmpstr), tmpstr);

      strcpy(tmpstr, "Time elapsed since the start of the forecast");
      cdf_put_att_text(fileID, leadtimeid, "long_name", strlen(tmpstr), tmpstr);

      cdfDefForecastTimeUnits(unitstr, taxis->fc_unit);

      len = strlen(unitstr);
      if ( len ) cdf_put_att_text(fileID, leadtimeid, "units", len, unitstr);
    }

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


void cdfDefTimestep(stream_t *streamptr, int