stream_cdf.c

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


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

#include "dmemory.h"

#include "cdi.h"
#include "basetime.h"
#include "util.h"  /* gaussaw */
#include "stream_int.h"
#include "stream_cdf.h"
#include "cdf_int.h"
#include "varscan.h"
#include "vlist.h"


#if  defined  (HAVE_LIBNETCDF)
#  include "netcdf.h"
#endif

#undef  UNDEFID
#define UNDEFID  CDI_UNDEFID


void cdfDefTime(int streamID);
void cdfDefGlobalAtts(int streamID);
void cdfDefLocalAtts(int streamID);

#define  MAXNAMELEN  256

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

typedef struct {
  int     ncvarid;
  int     dimtype;
  size_t  len;
  char    name[MAXNAMELEN];
}
NCDIM;

typedef struct {
  int      vartype;
  int      isvar;
  int      islon;
  int      islat;
  int      islev;
  int      warn;
  int      timeID;
  int      code;
  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      varids[4];
  int      cellarea;
  int      calendar;
  int      tableID;
  int      truncation;
  int      defmiss;
  int      xtype;
  int      ndims;
  int      gmapid;
  int      positive;
  int      dimids[8];
  int      dimtype[8];
  int      natts;
  int     *atts;
  size_t   vlen;
  double  *vdata;
  double   missval;
  double   addoffset;
  double   scalefactor;
  char     name[MAXNAMELEN];
  char     longname[MAXNAMELEN];
  char     stdname[MAXNAMELEN];
  char     units[MAXNAMELEN];
}
NCVAR;


static int splitBasetime(char *timeunits, TAXIS *taxis)
{
  static char func[] = "splitBasetime";
  int len, i;
  char *ptu;
  int year, month, day;
  int hour = 0, minute = 0, second = 0;
  int timetype = TAXIS_ABSOLUTE;
  int rdate = -1, rtime = -1;
  int timeunit = -1;

  ptu = timeunits;
  len = (int) strlen(timeunits);

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

  (*taxis).type = timetype;

  if      ( strncmp(ptu, "sec",    3) == 0 ) timeunit = TUNIT_SECOND;
  else if ( strncmp(ptu, "minute", 6) == 0 ) timeunit = TUNIT_MINUTE;
  else if ( strncmp(ptu, "hour",   4) == 0 ) timeunit = TUNIT_HOUR;
  else if ( strncmp(ptu, "day",    3) == 0 ) timeunit = TUNIT_DAY;
  else if ( strncmp(ptu, "month",  5) == 0 ) timeunit = TUNIT_MONTH;
  else if ( strncmp(ptu, "calendar_month", 14) == 0 ) timeunit = TUNIT_MONTH;
  else if ( strncmp(ptu, "year",   4) == 0 ) timeunit = TUNIT_YEAR;
  else
    {
      Message(func, "Unsupported TIMEUNIT: %s!", ptu);
      return (1);
    }

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

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

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

	  if ( timetype == TAXIS_ABSOLUTE )
	    {
	      if ( strncmp(ptu, "%y%m%d.%f", 9) != 0 && timeunit == TUNIT_DAY )
		{
		  Message(func, "Unsupported format %s for TIMEUNIT day!", ptu);
		  timeunit = -1;
		}
	      else if ( strncmp(ptu, "%y%m.%f", 7) != 0 && timeunit == TUNIT_MONTH )
		{
		  Message(func, "Unsupported format %s for TIMEUNIT month!", ptu);
		  timeunit = -1;
		}
	    }
	  else if ( timetype == TAXIS_RELATIVE )
	    {
	      year = atoi(ptu);
	      if ( year < 0 ) ptu++;
	      while ( isdigit((int) *ptu) ) ptu++;
	      month = atoi(++ptu);
	      while ( isdigit((int) *ptu) ) ptu++;
	      day = atoi(++ptu);
	      while ( isdigit((int) *ptu) ) ptu++;

	      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);
			  if ( second != 0 )
			    Message(func, "Seconds not supported in time units!");
			}
		    }
		}

	      rdate = encode_date(year, month, day);
	      rtime = encode_time(hour, minute);
	      (*taxis).rdate = rdate;
	      (*taxis).rtime = rtime;

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

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

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

  return (0);
}


#if  defined  (HAVE_LIBNETCDF)
static void cdfGetAttInt(int fileID, int ncvarid, char *attname, int *attint)
{
  static char func[] = "cdfGetAttInt";
  size_t attlen;
  int intatt;
  int *pintatt;

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

  if ( attlen > 1 )
    pintatt = (int *) malloc(attlen*sizeof(int));
  else
    pintatt = &intatt;

  cdf_get_att_int(fileID, ncvarid, attname, pintatt);

  *attint = *pintatt;

  if ( attlen > 1 )
    free(pintatt);
}
#endif

#if  defined  (HAVE_LIBNETCDF)
static void cdfGetAttDouble(int fileID, int ncvarid, char *attname, double *attdouble)
{
  static char func[] = "cdfGetAttInt";
  size_t attlen;
  double doubleatt;
  double *pdoubleatt;

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

  if ( attlen > 1 )
    pdoubleatt = (double *) malloc(attlen*sizeof(double));
  else
    pdoubleatt = &doubleatt;

  cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);

  *attdouble = *pdoubleatt;

  if ( attlen > 1 )
    free(pdoubleatt);
}
#endif

#if  defined  (HAVE_LIBNETCDF)
static void cdfGetAttText(int fileID, int ncvarid, char *attname, char *atttext, int maxlen)
{
  size_t attlen;
  char attbuf[65636];

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

  if ( attlen < sizeof(attbuf) )
    {
      cdf_get_att_text(fileID, ncvarid, attname, attbuf);

      attbuf[attlen++] = 0;

      if ( (int) attlen > maxlen ) attlen = maxlen;
      memcpy(atttext, attbuf, attlen);
    }
  else
    {
      atttext[0] = 0;
    }
}
#endif

#if  defined  (HAVE_LIBNETCDF)
int cdfInqDatatype(int xtype)
{
  int datatype;

  if      ( xtype == NC_BYTE  )  datatype = DATATYPE_INT8;
  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                           datatype = DATATYPE_FLT64;

  return (datatype);
}
#endif

#if  defined  (HAVE_LIBNETCDF)
int cdfDefDatatype(int datatype)
{
  int xtype;

  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_FLT64 )  xtype = NC_DOUBLE;
  else                                    xtype = NC_FLOAT;

  return (xtype);
}
#endif


#if  defined  (HAVE_LIBNETCDF)
void defineAttributes(int vlistID, int varID, int fileID, int ncvarID)
{
  int natts, iatt;
  int attrtype, attrlen;
  size_t len;
  int attrint;
  double attrflt;
  char attrname[1024];
  char attrtxt[8192];

  vlistInqNatts(vlistID, varID, &natts);

  for ( iatt = 0; iatt < natts; iatt++ )
    {
      vlistInqAttr(vlistID, varID, iatt, attrname, &attrtype, &attrlen);

      if ( attrtype == DATATYPE_TXT )
	{
	  vlistInqAttrTxt(vlistID, varID, attrname, sizeof(attrtxt), attrtxt);
	  len = strlen(attrtxt);
	  cdf_put_att_text(fileID, ncvarID, attrname, len+1, attrtxt);
	}
      else if ( attrtype == DATATYPE_INT )
	{
	  vlistInqAttrInt(vlistID, varID, attrname, 1, &attrint);
	  cdf_put_att_int(fileID, ncvarID, attrname, NC_INT, 1L, &attrint);
	}
      else if ( attrtype == DATATYPE_FLT )
	{
	  vlistInqAttrFlt(vlistID, varID, attrname, 1, &attrflt);
	  cdf_put_att_double(fileID, ncvarID, attrname, NC_DOUBLE, 1L, &attrflt);
	}
    }
}
#endif


int cdfCopyRecord(int streamID2, int streamID1)
{
  static char func[] = "cdfCopyRecord";
  double *data;
  int datasize;
  int tsID1, tsID2, recID1, recID2;
  int ivarID, gridID;
  int nmiss;
  int ierr = 0;
  int vlistID1, vlistID2;
  STREAM *streamptr1;
  STREAM *streamptr2;

  streamptr1 = stream_to_pointer(streamID1);
  streamptr2 = stream_to_pointer(streamID2);

  stream_check_ptr(func, streamptr1);
  stream_check_ptr(func, streamptr2);

  vlistID1 = streamptr1->vlistID;
  vlistID2 = streamptr2->vlistID;

  tsID1 = streamptr1->curTsID;
  tsID2 = streamptr2->curTsID;

  recID1 = streamptr1->tsteps[tsID1].curRecID;
  recID2 = streamptr2->tsteps[tsID2].curRecID;

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

  gridID = vlistInqVarGrid(vlistID1, ivarID);

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

  data = (double *) malloc(datasize*sizeof(double));

  streamReadRecord(streamID1, data, &nmiss);
  streamWriteRecord(streamID2, data, nmiss);

  free(data);

  return (ierr);
}

/* not used
int cdfInqRecord(int streamID, int *varID, int *levelID)
{
  static char func[] = "cdfInqRecord";
  int tsID, recID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  stream_check_ptr(func, streamptr);

  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(func, "recID = %d  varID = %d  levelID = %d", recID, *varID, *levelID);
  
  return (recID+1);
}
*/
int cdfDefRecord(int streamID)
{
  static char func[] = "cdfDefRecord";
  int ierr = 0;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  if ( CDI_Debug )
    Message(func, "streamID = %d", streamID);

  stream_check_ptr(func, streamptr);

  return (ierr);
}

void cdfWriteGridTraj(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  int tsID, fileID;
  int lonID, latID, gridindex;
  size_t index;
  double xlon, xlat;
  int vlistID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  vlistID = streamInqVlist(streamID);
  fileID  = streamInqFileID(streamID);

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

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

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


void cdfReadGridTraj(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  int tsID, fileID;
  int lonID, latID, gridindex;
  size_t index;
  double xlon, xlat;
  int vlistID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  vlistID = streamInqVlist(streamID);
  fileID  = streamInqFileID(streamID);

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

  tsID = streamptr->curTsID;
  index = tsID;

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

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


void cdfDefMissval(int streamID, int varID, int dtype)
{
#if  defined  (HAVE_LIBNETCDF)
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  if ( streamptr->vars[varID].defmiss == FALSE )
    {
      int fileID;
      int ncvarid;
      double missval;
      int vlistID;
      int xtype;

      vlistID = streamInqVlist(streamID);
      fileID  = streamInqFileID(streamID);
      ncvarid = streamptr->vars[varID].ncvarid;
      missval = vlistInqVarMissval(vlistID, varID);
      if ( streamptr->ncmode == 2 ) cdf_redef(fileID);

      xtype = cdfDefDatatype(dtype);

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

      if ( cdiNcMissingValue == 1 )
	cdf_put_att_double(fileID, ncvarid, "missing_value", (nc_type) xtype, 1L, &missval);

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

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

void cdfWriteRecord(int streamID, double *data, int nmiss)
{
#if  defined  (HAVE_LIBNETCDF)
  static char func[] = "cdfWriteRecord";
  int varID;
  int levelID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

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

  if ( CDI_Debug )
    Message(func, "streamID = %d  varID = %d", streamID, varID);

  cdfWriteVarSliceDP(streamID, varID, levelID, data, nmiss);
#endif
}

int cdfReadRecord(int streamID, double *data, int *nmiss)
{
  static char func[] = "cdfReadRecord";
  int ierr = 0;
  int levelID, varID, tsID, recID, vrecID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  if ( CDI_Debug ) Message(func, "streamID = %d", streamID);

  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(streamID, varID, levelID, data, nmiss);

  return (ierr);
}


void cdfDefTimeValue(int streamID, int tsID)
{
#if  defined  (HAVE_LIBNETCDF)
  static char func[] = "cdfDefTimeValue";
  int fileID;
  double timevalue;
  int ncvarid;
  size_t index;
  TAXIS *taxis;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  fileID = streamInqFileID(streamID);

  if ( CDI_Debug )
    Message(func, "streamID = %d, fileID = %d", streamID, fileID);

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

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

  index = tsID;

  timevalue = encode_timeval(taxis->vdate, taxis->vtime, &streamptr->tsteps[0].taxis);
  if ( CDI_Debug ) Message(func, "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 = encode_timeval(taxis->vdate_lb, taxis->vtime_lb, &streamptr->tsteps[0].taxis);
      start[0] = tsID; count[0] = 1; start[1] = 0; count[1] = 1;
      cdf_put_vara_double(fileID, ncvarid, start, count, &timevalue);

      timevalue = encode_timeval(taxis->vdate_ub, taxis->vtime_ub, &streamptr->tsteps[0].taxis);
      start[0] = tsID; count[0] = 1; start[1] = 1; count[1] = 1;
      cdf_put_vara_double(fileID, ncvarid, start, count, &timevalue);
    }
  /*
printf("fileID = %d %d %d %f\n", fileID, time_varid, index, timevalue);
  */
#endif
}


void cdfDefTimestep(int streamID, int tsID)
{
  int vlistID;

  vlistID = streamInqVlist(streamID);

  if ( vlistHasTime(vlistID) ) cdfDefTime(streamID);

  cdfDefTimeValue(streamID, tsID);
}


void cdfDefTime(int streamID)
{
#if  defined  (HAVE_LIBNETCDF)
  int fileID;
  int time_varid;
  int time_bndsid;
  int dims[2];
  int year, month, day, hour, minute;
  char unitstr[80];
  char calstr[80];
  size_t len;
  TAXIS *taxis;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

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

  fileID = streamInqFileID(streamID);

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

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

  cdf_def_dim(fileID, "time", NC_UNLIMITED, &streamptr->basetime.ncdimid);

  dims[0] = streamptr->basetime.ncdimid;
  cdf_def_var(fileID, "time", NC_DOUBLE, 1, dims, &time_varid);

  streamptr->basetime.ncvarid = time_varid;

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

  if ( taxis->has_bounds )
    {
      /* fprintf(stderr, "time has bounds\n"); */

      cdf_def_dim(fileID, "tbnds", 2, &dims[1]);
      cdf_def_var(fileID, "time_bnds", NC_DOUBLE, 2, dims, &time_bndsid);

      streamptr->basetime.ncvarboundsid = time_bndsid;

      cdf_put_att_text(fileID, time_varid, "bounds", 9, "time_bnds");
    }

  unitstr[0] = 0;
  if ( streamptr->tsteps[0].taxis.type == TAXIS_ABSOLUTE )
    {
      if ( streamptr->tsteps[0].taxis.unit == TUNIT_MONTH )
	sprintf(unitstr, "month as %s", "%Y%m.%f");
      else
	sprintf(unitstr, "day as %s", "%Y%m%d.%f");
    }
  else
    {
      int rdate, rtime;
      int timeunit;

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

      decode_date(rdate, &year, &month, &day);
      decode_time(rtime, &hour, &minute);

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

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

  if ( taxis->has_bounds )
    if ( len )
      cdf_put_att_text(fileID, time_bndsid, "units", len, unitstr);

  if ( taxis->calendar != -1 )
    {
      calstr[0] = 0;

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

      len = strlen(calstr);
      if ( len )
	cdf_put_att_text(fileID, time_varid, "calendar", len, calstr);

      if ( taxis->has_bounds )
	if ( len )
	  cdf_put_att_text(fileID, time_bndsid, "calendar", len, calstr);
    }

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


void cdfDefComplex(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  char axisname[] = "complex";
  int index;
  int dimID = UNDEFID;
  int gridID0, gridtype0, gridindex;
  int ngrids;
  int fileID;
  int dimlen;
  int vlistID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  vlistID = streamInqVlist(streamID);
  fileID  = streamInqFileID(streamID);

  ngrids = vlistNgrids(vlistID);

  for ( index = 0; index < ngrids; index++ )
    {
      if ( streamptr->xdimID[index] != UNDEFID )
	{
	  gridID0 = vlistGrid(vlistID, index);
	  gridtype0 = gridInqType(gridID0);
	  if ( gridtype0 == GRID_SPECTRAL )
	    {
	      dimID = streamptr->xdimID[index];
	      break;
	    }
	}
    }

  if ( dimID == UNDEFID )
    {
      dimlen = 2;

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

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

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

  gridindex = vlistGridIndex(vlistID, gridID);
  streamptr->xdimID[gridindex] = dimID;

#endif
}

void cdfDefSpc(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  /*
  char longname[] = "Spherical harmonic coefficient";
  */
  char axisname[5] = "spcX";
  int index, iz = 0;
  int gridID0, gridtype0, gridindex;
  int dimID = UNDEFID;
  int ngrids;
  int fileID;
  int dimlen, dimlen0;
  int vlistID;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  vlistID = streamInqVlist(streamID);
  fileID  = streamInqFileID(streamID);

  ngrids = vlistNgrids(vlistID);

  dimlen = gridInqSize(gridID)/2;

  for ( index = 0; index < ngrids; index++ )
    {
      if ( streamptr->ydimID[index] != UNDEFID )
	{
	  gridID0 = vlistGrid(vlistID, index);
	  gridtype0 = gridInqType(gridID0);
	  if ( gridtype0 == GRID_SPECTRAL )
	    {
	      dimlen0 = gridInqSize(gridID0)/2;
	      if ( dimlen == dimlen0 )
		{
		  dimID = streamptr->ydimID[index];
		  break;
		}
	      else
		iz++;	
	    }
	}
    }

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

  gridindex = vlistGridIndex(vlistID, gridID);
  streamptr->ydimID[gridindex] = dimID;

#endif
}

void cdfDefTrajLon(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  static char func[] = "cdfDefTrajLon";
  char units[256];
  char longname[256];
  char stdname[256];
  char axisname[256];
  int gridtype, gridindex;
  int dimID = UNDEFID;
  int fileID;
  int dimlen;
  size_t len;
  int ncvarid;
  int vlistID;
  int xtype = NC_DOUBLE;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  if ( gridInqPrec(gridID) == DATATYPE_FLT32 ) xtype = NC_FLOAT;

  vlistID = streamInqVlist(streamID);
  fileID  = streamInqFileID(streamID);

  gridtype = gridInqType(gridID);
  dimlen = gridInqXsize(gridID);
  if ( dimlen != 1 ) Error(func, "Xsize isn't 1 for %s grid!", gridNamePtr(gridtype));

  gridindex = vlistGridIndex(vlistID, gridID);
  ncvarid = streamptr->xdimID[gridindex];

  gridInqXname(gridID, axisname);
  gridInqXlongname(gridID, longname);
  gridInqXstdname(gridID, stdname);
  gridInqXunits(gridID, units);

  if ( ncvarid == UNDEFID )
    {
      dimID = streamptr->basetime.ncvarid;

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

      cdf_def_var(fileID, axisname, (nc_type) xtype, 1, &dimID, &ncvarid);
      if ( (len = strlen(longname)) )
	cdf_put_att_text(fileID, ncvarid, "long_name", len, longname);
      if ( (len = strlen(units)) )
	cdf_put_att_text(fileID, ncvarid, "units", len, units);
      if ( (len = strlen(stdname)) )
	cdf_put_att_text(fileID, ncvarid, "standard_name", len, stdname);

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

  streamptr->xdimID[gridindex] = ncvarid; /* var ID for trajectory !!! */

#endif
}

void cdfDefTrajLat(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  static char func[] = "cdfDefTrajLat";
  char units[] = "degrees_north";
  char longname[] = "latitude";
  char stdname[] = "latitude";
  char axisname[] = "tlat";
  int gridtype, gridindex;
  int dimID = UNDEFID;
  int fileID;
  int dimlen;
  size_t len;
  int ncvarid;
  int vlistID;
  int xtype = NC_DOUBLE;
  STREAM *streamptr;

  streamptr = stream_to_pointer(streamID);

  if ( gridInqPrec(gridID) == DATATYPE_FLT32 ) xtype = NC_FLOAT;

  vlistID = streamInqVlist(streamID);
  fileID = streamInqFileID(streamID);

  gridtype = gridInqType(gridID);
  dimlen = gridInqYsize(gridID);
  if ( dimlen != 1 ) Error(func, "Ysize isn't 1 for %s grid!", gridNamePtr(gridtype));

  gridindex = vlistGridIndex(vlistID, gridID);
  ncvarid = streamptr->ydimID[gridindex];

  gridInqYname(gridID, axisname);
  gridInqYlongname(gridID, longname);
  gridInqYstdname(gridID, stdname);
  gridInqYunits(gridID, units);

  if ( ncvarid == UNDEFID )
    {
      dimID = streamptr->basetime.ncvarid;

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

      cdf_def_var(fileID, axisname, (nc_type) xtype, 1, &dimID, &ncvarid);
      if ( (len = strlen(longname)) )
	cdf_put_att_text(fileID, ncvarid, "long_name", len, longname);
      if ( (len = strlen(units)) )
	cdf_put_att_text(fileID, ncvarid, "units", len, units);
      if ( (len = strlen(stdname)) )
	cdf_put_att_text(fileID, ncvarid, "standard_name", len, stdname);

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

  streamptr->ydimID[gridindex] = ncvarid; /* var ID for trajectory !!! */

#endif
}


void cdfDefXaxis(int streamID, int gridID)
{
#if  defined  (HAVE_LIBNETCDF)
  static char func[] = "cdfDefXaxis";