stream_cdf_i.c

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_LIBNETCDF

#include <ctype.h>
#include <limits.h>

#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 "cdf_util.h"
#include "cdf_lazy_grid.h"


#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     dimid;
  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      isvar;
  bool     ignore;
  bool     isx;
  bool     isy;
  bool     isc;
  bool     islon;
  bool     islat;
  bool     islev;
  bool     istime;
  bool     warn;
  bool     calendar;
  bool     climatology;
  bool     lformulaterms;
  int      param;
  int      code;
  int      tabnum;
  int      bounds;
  int      gridID;
  int      zaxisID;
  int      timetype;
  int      gridtype;
  int      zaxistype;
  int      xdim;
  int      ydim;
  int      zdim;
  int      xvarid;
  int      yvarid;
  int      rpvarid;
  int      zvarid;
  int      cvarids[MAX_COORDVARS];
  int      tvarid;
  int      psvarid;
  int      p0varid;
  int      ncoordvars;
  int      coordvarids[MAX_COORDVARS];
  int      nauxvars;
  int      auxvarids[MAX_AUXVARS];
  int      cellarea;
  int      tableID;
  int      truncation;
  int      position;
  int      numLPE;
  bool     defmissval;
  bool     deffillval;
  int      xtype;
  int      gmapid;
  int      positive;
  int      ndims;
  int      dimids[8];
  int      dimtype[8];
  size_t   chunks[8];
  int      chunked;
  int      chunktype;
  int      natts;
  int      deflate;
  bool     lunsigned;
  bool     lvalidrange;
  int     *atts;
  size_t   vctsize;
  double  *vct;
  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];
  int      typeOfEnsembleForecast;
  int      numberOfForecastsInEnsemble;
  int      perturbationNumber;
}
ncvar_t;


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

  if ( *ptu ) sscanf(ptu, "%d-%d-%d %d:%d:%d", &v1, &v2, &v3, &hour, &minute, &second);

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

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

static
int scanTimeUnit(const char *unitstr)
{
  const size_t len = strlen(unitstr);
  const int timeunit = get_timeunit(len, unitstr);
  if ( timeunit == -1 ) Message("Unsupported TIMEUNIT: %s!", unitstr);
  return timeunit;
}

static
void setForecastTime(const char *timestr, taxis_t *taxis)
{
  const size_t len = strlen(timestr);
  if ( len != 0 )
    scanTimeString(timestr, &taxis->fdate, &taxis->ftime);
  else
    taxis->fdate = taxis->ftime = 0;
}

static
int setBaseTime(const char *timeunits, taxis_t *taxis)
{
  int taxistype = TAXIS_ABSOLUTE;

  size_t len = strlen(timeunits);
  while ( isspace(*timeunits) && len ) { timeunits++; len--; }

  char *tu = (char *)Malloc((len+1) * sizeof(char));

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

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

  size_t pos = 0;
  while ( pos < len && !isspace(tu[pos]) ) ++pos;
  if ( tu[pos] )
    {
      while ( isspace(tu[pos]) ) ++pos;

      if ( strStartsWith(tu+pos, "since") ) taxistype = TAXIS_RELATIVE;

      while ( pos < len && !isspace(tu[pos]) ) ++pos;
      if ( tu[pos] )
        {
          while ( isspace(tu[pos]) ) ++pos;

          if ( taxistype == TAXIS_ABSOLUTE )
            {
              if ( timeunit == TUNIT_DAY )
                {
                  if ( !strStartsWith(tu+pos, "%y%m%d.%f") )
                    {
                      Warning("Unsupported format %s for TIMEUNIT day!", tu+pos);
                      timeunit = -1;
                    }
                }
              else if ( timeunit == TUNIT_MONTH )
                {
                  if ( !strStartsWith(tu+pos, "%y%m.%f") )
                    {
                      Warning("Unsupported format %s for TIMEUNIT month!", tu+pos);
                      timeunit = -1;
                    }
                }
              else if ( timeunit == TUNIT_YEAR )
                {
                  if ( !strStartsWith(tu+pos, "%y.%f") )
                    {
                      Warning("Unsupported format %s for TIMEUNIT year!", tu+pos);
                      timeunit = -1;
                    }
                }
              else
                {
                  Warning("Unsupported format for time units: %s!", tu);
                }
            }
          else if ( taxistype == TAXIS_RELATIVE )
            {
              int64_t rdate = -1;
              int rtime = -1;
              scanTimeString(tu+pos, &rdate, &rtime);
              taxis->rdate = rdate;
              taxis->rtime = rtime;

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

  taxis->type = taxistype;
  taxis->unit = timeunit;

  Free(tu);

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

  return 0;
}

static
bool xtypeIsText(int xtype)
{
  const bool isText = (xtype == NC_CHAR)
#ifdef HAVE_NETCDF4
    || (xtype == NC_STRING)
#endif
    ;
  return isText;
}

static
bool xtypeIsFloat(nc_type xtype)
{
  const bool isFloat = xtype == NC_FLOAT || xtype == NC_DOUBLE;
  return isFloat;
}

static
bool xtypeIsInt(nc_type xtype)
{
  const bool isInt = xtype == NC_SHORT || xtype == NC_INT
            || xtype == NC_BYTE
#ifdef HAVE_NETCDF4
            || xtype == NC_USHORT || xtype == NC_UINT
            || xtype == NC_UBYTE
#endif
             ;

  return isInt;
}

static
int cdfInqDatatype(stream_t *streamptr, int xtype, bool lunsigned)
{
  int datatype = -1;

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

  if      ( xtype == NC_BYTE   )  datatype = CDI_DATATYPE_INT8;
  else if ( xtype == NC_CHAR   )  datatype = CDI_DATATYPE_UINT8;
  else if ( xtype == NC_SHORT  )  datatype = CDI_DATATYPE_INT16;
  else if ( xtype == NC_INT    )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_FLOAT  )  datatype = CDI_DATATYPE_FLT32;
  else if ( xtype == NC_DOUBLE )  datatype = CDI_DATATYPE_FLT64;
#ifdef HAVE_NETCDF4
  else if ( xtype == NC_UBYTE  )  datatype = CDI_DATATYPE_UINT8;
  else if ( xtype == NC_LONG   )  datatype = CDI_DATATYPE_INT32;
  else if ( xtype == NC_USHORT )  datatype = CDI_DATATYPE_UINT16;
  else if ( xtype == NC_UINT   )  datatype = CDI_DATATYPE_UINT32;
  else if ( xtype == NC_INT64  )  datatype = CDI_DATATYPE_FLT64;
  else if ( xtype == NC_UINT64 )  datatype = CDI_DATATYPE_FLT64;
  else
    {
      if (xtype != streamptr->nc_complex_float_id && xtype != streamptr->nc_complex_double_id)
        {
          int isUserDefinedType = false;
#ifdef NC_FIRSTUSERTYPEID
          isUserDefinedType = xtype >= NC_FIRSTUSERTYPEID;
#endif
          if (isUserDefinedType)
            {
              int fileID = streamptr->fileID;
              size_t nfields = 0, compoundsize = 0;
              int status = nc_inq_compound(fileID, xtype, NULL, &compoundsize, &nfields);
              if (status == NC_NOERR && nfields == 2 && (compoundsize == 8 || compoundsize == 16))
                {
                  nc_type field_type1 = -1, field_type2 = -1;
                  int field_dims1 = 0, field_dims2 = 0;
                  nc_inq_compound_field(fileID, xtype, 0, NULL, NULL, &field_type1, &field_dims1, NULL);
                  nc_inq_compound_field(fileID, xtype, 1, NULL, NULL, &field_type2, &field_dims2, NULL);
                  if (field_type1 == field_type2 && field_dims1 == 0 && field_dims2 == 0)
                    {
                      if      (field_type1 == NC_FLOAT)  streamptr->nc_complex_float_id = xtype;
                      else if (field_type1 == NC_DOUBLE) streamptr->nc_complex_double_id = xtype;
                    }
                }
            }
        }
      if      ( xtype == streamptr->nc_complex_float_id  )  datatype = CDI_DATATYPE_CPX32;
      else if ( xtype == streamptr->nc_complex_double_id )  datatype = CDI_DATATYPE_CPX64;
    }
#endif

  return datatype;
}

static
void cdfGetAttInt(int fileID, int ncvarid, const char *attname, size_t attlen, int *attint)
{
  *attint = 0;

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

  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
    {
      const bool lalloc = nc_attlen > attlen;
      int *pintatt = lalloc ? (int *)(Malloc(nc_attlen*sizeof(int))) : attint;
      cdf_get_att_int(fileID, ncvarid, attname, pintatt);
      if ( lalloc )
        {
          memcpy(attint, pintatt, attlen*sizeof(int));
          Free(pintatt);
        }
    }
}

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

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

  if ( xtypeIsFloat(atttype) || xtypeIsInt(atttype) )
    {
      const bool lalloc = nc_attlen > attlen;
      double *pdoubleatt = lalloc ? (double*)Malloc(nc_attlen*sizeof(double)) : attdouble;
      cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);
      if ( lalloc )
        {
          memcpy(attdouble, pdoubleatt, attlen*sizeof(double));
          Free(pdoubleatt);
        }
    }
}

static
bool cdfCheckAttText(int fileID, int ncvarid, const char *attname)
{
  bool status = false;
  nc_type atttype;

  const int status_nc = nc_inq_atttype(fileID, ncvarid, attname, &atttype);
  if ( status_nc == NC_NOERR && (atttype == NC_CHAR
#ifdef HAVE_NETCDF4
           || atttype == NC_STRING
#endif
           ) )
    {
      status = true;
    }

  return status;
}

static
void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
{
  atttext[0] = 0;

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

          if ( nc_attlen > (attlen-1) ) nc_attlen = (attlen-1);

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

          size_t ssize = strlen(attbuf) + 1;
          if ( ssize > attlen ) ssize = attlen;
          memcpy(atttext, attbuf, ssize);
          atttext[ssize - 1] = 0;
          Free(attbuf);
        }
    }
#endif
}


void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
{
  const bool laddoffset   = IS_NOT_EQUAL(addoffset, 0);
  const bool lscalefactor = IS_NOT_EQUAL(scalefactor, 1);

  if ( laddoffset && lscalefactor )
    {
      for (size_t i = 0; i < size; ++i) data[i] = data[i] * scalefactor + addoffset;
    }
  else if (lscalefactor)
    {
      for (size_t i = 0; i < size; ++i) data[i] *= scalefactor;
    }
  else if (laddoffset)
    {
      for (size_t i = 0; i < size; ++i) data[i] += addoffset;
    }
}

static
void cdfCreateRecords(stream_t *streamptr, int tsID)
{
  if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;

  if ( streamptr->tsteps[tsID].nallrecs > 0 ) return;

  int vlistID  = streamptr->vlistID;

  tsteps_t* sourceTstep = streamptr->tsteps;
  tsteps_t* destTstep = sourceTstep + tsID;

  const int nvars = vlistNvars(vlistID);
  const int nrecs = vlistNrecs(vlistID);
  if ( nrecs <= 0 ) return;

  if ( tsID == 0 )
    {
      const int nvrecs = nrecs; /* use all records at first timestep */

      streamptr->nrecs += nrecs;

      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
      destTstep->nrecs      = nrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = CDI_UNDEFID;
      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs*sizeof (int));;
      for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;

      record_t *records = destTstep->records;

      for ( int varID = 0, recID = 0; varID < nvars; varID++ )
        {
          const int zaxisID = vlistInqVarZaxis(vlistID, varID);
          const int nlev    = zaxisInqSize(zaxisID);
          for ( int levelID = 0; levelID < nlev; levelID++ )
            {
              recordInitEntry(&records[recID]);
              records[recID].varID   = (short)varID;
              records[recID].levelID = (short)levelID;
              recID++;
            }
        }
    }
  else if ( tsID == 1 )
    {
      int nvrecs = 0;
      for ( int varID = 0; varID < nvars; varID++ )
        {
          if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
            {
              int zaxisID = vlistInqVarZaxis(vlistID, varID);
              nvrecs += zaxisInqSize(zaxisID);
            }
        }

      streamptr->nrecs += nvrecs;

      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = CDI_UNDEFID;

      memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));

      if ( nvrecs )
        {
          destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
          for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
            {
              int varID = destTstep->records[recID].varID;
              if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
                {
                  destTstep->recIDs[vrecID++] = recID;
                }
            }
        }
    }
  else
    {
      if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);

      const int nvrecs = streamptr->tsteps[1].nrecs;

      streamptr->nrecs += nvrecs;

      destTstep->records    = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
      destTstep->nrecs      = nvrecs;
      destTstep->nallrecs   = nrecs;
      destTstep->recordSize = nrecs;
      destTstep->curRecID   = CDI_UNDEFID;

      memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));

      destTstep->recIDs     = (int *) Malloc((size_t)nvrecs * sizeof(int));

      memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
    }
}

static
int cdf_time_dimid(int fileID, int ndims, int nvars, ncdim_t *ncdims)
{
  char dimname[CDI_MAX_NAME];

  for ( int dimid = 0; dimid < ndims; ++dimid )
    {
      dimname[0] = 0;
      cdf_inq_dimname(fileID, ncdims[dimid].dimid, dimname);
      if ( strIsEqual(dimname, "time") || strIsEqual(dimname, "Time") ) return dimid;
    }

  for ( int varid = 0; varid < nvars; ++varid )
    {
      nc_type xtype;
      int nvdims, nvatts, dimids[9];
      cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
      for ( int i = 0; i < nvdims; ++i )
        for ( int dimid = 0; dimid < ndims; ++dimid )
          if ( ncdims[dimid].dimid == dimids[i] )
            {
              dimids[i] = dimid;
              break;
            }

      if ( nvdims == 1 )
        {
          char sbuf[CDI_MAX_NAME];
          for ( int iatt = 0; iatt < nvatts; ++iatt )
            {
              sbuf[0] = 0;
              cdf_inq_attname(fileID, varid, iatt, sbuf);
              if ( strncmp(sbuf, "units", 5) == 0 )
                {
                  cdfGetAttText(fileID, varid, "units", sizeof(sbuf), sbuf);
                  strToLower(sbuf);

                  if ( is_time_units(sbuf) ) return dimids[0];
                }
            }
        }
    }

  return CDI_UNDEFID;
}

static
void init_ncdims(long ndims, ncdim_t *ncdims)
{
  for ( long ncdimid = 0; ncdimid < ndims; ncdimid++ )
    {
      ncdims[ncdimid].dimid        = CDI_UNDEFID;
      ncdims[ncdimid].ncvarid      = CDI_UNDEFID;
      ncdims[ncdimid].dimtype      = CDI_UNDEFID;
      ncdims[ncdimid].len          = 0;
      ncdims[ncdimid].name[0]      = 0;
    }
}

static
void init_ncvars(long nvars, ncvar_t *ncvars)
{
  for ( long ncvarid = 0; ncvarid < nvars; ++ncvarid )
    {
      ncvars[ncvarid].ncid            = CDI_UNDEFID;
      ncvars[ncvarid].isvar           = CDI_UNDEFID;
      ncvars[ncvarid].ignore          = false;
      ncvars[ncvarid].isx             = false;
      ncvars[ncvarid].isy             = false;
      ncvars[ncvarid].isc             = false;
      ncvars[ncvarid].islon           = false;
      ncvars[ncvarid].islat           = false;
      ncvars[ncvarid].islev           = false;
      ncvars[ncvarid].istime          = false;
      ncvars[ncvarid].warn            = false;
      ncvars[ncvarid].calendar        = false;
      ncvars[ncvarid].climatology     = false;
      ncvars[ncvarid].lformulaterms   = false;
      ncvars[ncvarid].timetype        = TIME_CONSTANT;
      ncvars[ncvarid].param           = CDI_UNDEFID;
      ncvars[ncvarid].code            = CDI_UNDEFID;
      ncvars[ncvarid].tabnum          = 0;
      ncvars[ncvarid].bounds          = CDI_UNDEFID;
      ncvars[ncvarid].gridID          = CDI_UNDEFID;
      ncvars[ncvarid].zaxisID         = CDI_UNDEFID;
      ncvars[ncvarid].gridtype        = CDI_UNDEFID;
      ncvars[ncvarid].zaxistype       = CDI_UNDEFID;
      ncvars[ncvarid].xdim            = CDI_UNDEFID;
      ncvars[ncvarid].ydim            = CDI_UNDEFID;
      ncvars[ncvarid].zdim            = CDI_UNDEFID;
      ncvars[ncvarid].xvarid          = CDI_UNDEFID;
      ncvars[ncvarid].yvarid          = CDI_UNDEFID;
      ncvars[ncvarid].rpvarid         = CDI_UNDEFID;
      ncvars[ncvarid].zvarid          = CDI_UNDEFID;
      ncvars[ncvarid].tvarid          = CDI_UNDEFID;
      ncvars[ncvarid].psvarid         = CDI_UNDEFID;
      ncvars[ncvarid].p0varid         = CDI_UNDEFID;
      ncvars[ncvarid].ncoordvars      = 0;
      for ( int i = 0; i < MAX_COORDVARS; ++i )
        {
          ncvars[ncvarid].coordvarids[i]  = CDI_UNDEFID;
          ncvars[ncvarid].cvarids[i]      = CDI_UNDEFID;
        }
      ncvars[ncvarid].nauxvars      = 0;
      for ( int i = 0; i < MAX_AUXVARS; ++i )
        ncvars[ncvarid].auxvarids[i]  = CDI_UNDEFID;
      ncvars[ncvarid].cellarea        = CDI_UNDEFID;
      ncvars[ncvarid].tableID         = CDI_UNDEFID;
      ncvars[ncvarid].xtype           = 0;
      ncvars[ncvarid].ndims           = 0;
      ncvars[ncvarid].gmapid          = CDI_UNDEFID;
      ncvars[ncvarid].vctsize         = 0;
      ncvars[ncvarid].vct             = NULL;
      ncvars[ncvarid].truncation      = 0;
      ncvars[ncvarid].position        = 0;
      ncvars[ncvarid].numLPE          = 0;
      ncvars[ncvarid].positive        = 0;
      ncvars[ncvarid].chunked         = 0;
      ncvars[ncvarid].chunktype       = CDI_UNDEFID;
      ncvars[ncvarid].defmissval      = false;
      ncvars[ncvarid].deffillval      = false;
      ncvars[ncvarid].missval         = 0;
      ncvars[ncvarid].fillval         = 0;
      ncvars[ncvarid].addoffset       = 0;
      ncvars[ncvarid].scalefactor     = 1;
      ncvars[ncvarid].natts           = 0;
      ncvars[ncvarid].atts            = NULL;
      ncvars[ncvarid].deflate         = 0;
      ncvars[ncvarid].lunsigned       = false;
      ncvars[ncvarid].lvalidrange     = false;
      ncvars[ncvarid].validrange[0]   = VALIDMISS;
      ncvars[ncvarid].validrange[1]   = VALIDMISS;
      ncvars[ncvarid].typeOfEnsembleForecast      = -1;
      ncvars[ncvarid].numberOfForecastsInEnsemble = -1;
      ncvars[ncvarid].perturbationNumber          = -1;
      memset(ncvars[ncvarid].name, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].longname, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].stdname, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].units, 0, CDI_MAX_NAME);
      memset(ncvars[ncvarid].extra, 0, CDI_MAX_NAME);
    }
}

static
void cdf_set_var(ncvar_t *ncvars, int ncvarid, short isvar)
{
  if ( ncvars[ncvarid].isvar != CDI_UNDEFID &&
       ncvars[ncvarid].isvar != isvar   &&
       ncvars[ncvarid].warn  == false )
    {
      if ( ! ncvars[ncvarid].ignore )
        Warning("Inconsistent variable definition for %s!", ncvars[ncvarid].name);

      ncvars[ncvarid].warn = true;
      isvar = FALSE;
    }

  ncvars[ncvarid].isvar = isvar;
}

static
void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
{
  if ( ncvars[ncvarid].dimtype[dimid] != CDI_UNDEFID &&
       ncvars[ncvarid].dimtype[dimid] != dimtype )
    {
      Warning("Inconsistent dimension definition for %s! dimid = %d;  type = %d;  newtype = %d",
              ncvars[ncvarid].name, dimid, ncvars[ncvarid].dimtype[dimid], dimtype);
    }

  ncvars[ncvarid].dimtype[dimid] = dimtype;
}

static
void scan_hybrid_formulaterms(int ncid, int ncfvarid, int *avarid, int *bvarid, int *psvarid, int *p0varid)
{
  *avarid  = -1;
  *bvarid  = -1;
  *psvarid = -1;
  *p0varid = -1;

  char attstring[1024];
  cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
  char *pstring = attstring;

  bool lstop = false;
  for ( int i = 0; i < 4; i++ )
    {
      while ( isspace((int) *pstring) ) pstring++;
      if ( *pstring == 0 ) break;
      char *tagname = pstring;
      while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
      if ( *pstring == 0 ) lstop = true;
      *pstring++ = 0;

      while ( isspace((int) *pstring) ) pstring++;
      if ( *pstring == 0 ) break;
      char *varname = pstring;
      while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
      if ( *pstring == 0 ) lstop = true;
      *pstring++ = 0;

      int dimvarid;
      int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
      if ( status_nc == NC_NOERR )
        {
          if      ( strIsEqual(tagname, "ap:") ) *avarid  = dimvarid;
          else if ( strIsEqual(tagname, "a:")  ) *avarid  = dimvarid;
          else if ( strIsEqual(tagname, "b:")  ) *bvarid  = dimvarid;
          else if ( strIsEqual(tagname, "ps:") ) *psvarid = dimvarid;
          else if ( strIsEqual(tagname, "p0:") ) *p0varid = dimvarid;
        }
      else if ( !strIsEqual(tagname, "ps:") )
        {
          Warning("%s - %s", nc_strerror(status_nc), varname);
        }

      if ( lstop ) break;
    }
}

static
void readVCT(int ncid, int ndims2, size_t dimlen, size_t dimlen2, int avarid2, int bvarid2, double *vct)
{
  double *abuf = (double*) Malloc(dimlen*2*sizeof(double));
  double *bbuf = (double*) Malloc(dimlen*2*sizeof(double));
  cdf_get_var_double(ncid, avarid2, abuf);
  cdf_get_var_double(ncid, bvarid2, bbuf);

  if ( ndims2 == 2 )
    {
      for ( size_t i = 0; i < dimlen; ++i )
        {
          vct[i] = abuf[i*2];
          vct[i+dimlen+1] = bbuf[i*2];
        }
      vct[dimlen]     = abuf[dimlen*2-1];
      vct[dimlen*2+1] = bbuf[dimlen*2-1];
    }
  else
    {
      for ( size_t i = 0; i < dimlen2; ++i )
        {
          vct[i] = abuf[i];
          vct[i+dimlen+1] = bbuf[i];
        }
    }

  Free(abuf);
  Free(bbuf);
}

static
bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
{
  bool status = false;
  ncvar_t *ncvar = &ncvars[ncvarid];

  if ( strIsEqual(ncvar->stdname, "atmosphere_hybrid_sigma_pressure_coordinate") )
    {
      CDI_convention = CDI_CONVENTION_CF;

      status = true;
      ncvar->zaxistype = ZAXIS_HYBRID;
      //int ndims = ncvar->ndims;
      int dimid = ncvar->dimids[0];
      size_t dimlen = ncdims[dimid].len;
      int avarid1 = -1, bvarid1 = -1, psvarid1 = -1, p0varid1 = -1;
      int ncfvarid = ncvarid;
      if ( ncvars[ncfvarid].lformulaterms )
        scan_hybrid_formulaterms(ncid, ncfvarid, &avarid1, &bvarid1, &psvarid1, &p0varid1);
      // printf("avarid1, bvarid1, psvarid1, p0varid1 %d %d %d %d\n", avarid1, bvarid1, psvarid1, p0varid1);
      if ( avarid1  != -1 ) ncvars[avarid1].isvar = FALSE;
      if ( bvarid1  != -1 ) ncvars[bvarid1].isvar = FALSE;
      if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
      if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;