grid.c

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

#include <string.h>
#include <float.h>  /* FLT_EPSILON */
#include <limits.h> /* INT_MAX     */

#include "dmemory.h"
#include "cdi.h"
#include "cdi_cksum.h"
#include "cdi_int.h"
#include "grid.h"
#include "gaussgrid.h"
#include "resource_handle.h"
#include "resource_unpack.h"
#include "namespace.h"
#include "serialize.h"

#ifndef  RAD2DEG
#define  RAD2DEG  (180./M_PI)   /* conversion for rad to deg */
#endif

#ifndef  DEG2RAD
#define  DEG2RAD  (M_PI/180.)   /* conversion for deg to rad */
#endif


char *Grids[] = {
  /*  0 */  "undefined",
  /*  1 */  "generic",
  /*  2 */  "gaussian",
  /*  3 */  "gaussian reduced",
  /*  4 */  "lonlat",
  /*  5 */  "spectral",
  /*  6 */  "fourier",
  /*  7 */  "gme",
  /*  8 */  "trajectory",
  /*  9 */  "unstructured",
  /* 10 */  "curvilinear",
  /* 11 */  "lcc",
  /* 12 */  "lcc2",
  /* 13 */  "laea",
  /* 14 */  "sinusoidal",
  /* 15 */  "projection",
};


static int    gridCompareP    ( void * gridptr1, void * gridptr2 );
static void   gridDestroyP    ( void * gridptr );
static void   gridPrintP      ( void * gridptr, FILE * fp );
static int    gridGetPackSize ( void * gridptr, void *context);
static void   gridPack        ( void * gridptr, void * buff, int size,
				int *position, void *context);
static int    gridTxCode      ( void );

resOps gridOps = { gridCompareP, gridDestroyP, gridPrintP
                   , gridGetPackSize, gridPack, gridTxCode
};


static int  GRID_Debug = 0;   /* If set to 1, debugging */


void grid_init(grid_t *gridptr)
{
  gridptr->self         = CDI_UNDEFID;
  gridptr->type         = CDI_UNDEFID;
  gridptr->proj         = CDI_UNDEFID;
  gridptr->mask         = NULL;
  gridptr->mask_gme     = NULL;
  gridptr->xvals        = NULL;
  gridptr->yvals        = NULL;
  gridptr->area         = NULL;
  gridptr->xbounds      = NULL;
  gridptr->ybounds      = NULL;
  gridptr->rowlon       = NULL;
  gridptr->nrowlon      = 0;
  gridptr->xinc         = 0.0;
  gridptr->yinc         = 0.0;
  gridptr->lcc_originLon = 0.0;
  gridptr->lcc_originLat = 0.0;
  gridptr->lcc_lonParY  = 0.0;
  gridptr->lcc_lat1     = 0.0;
  gridptr->lcc_lat2     = 0.0;
  gridptr->lcc_xinc     = 0.0;
  gridptr->lcc_yinc     = 0.0;
  gridptr->lcc_projflag = 0;
  gridptr->lcc_scanflag = 0;
  gridptr->lcc_defined  = FALSE;
  gridptr->lcc2_lon_0   = 0.0;
  gridptr->lcc2_lat_0   = 0.0;
  gridptr->lcc2_lat_1   = 0.0;
  gridptr->lcc2_lat_2   = 0.0;
  gridptr->lcc2_a       = 0.0;
  gridptr->lcc2_defined = FALSE;
  gridptr->laea_lon_0   = 0.0;
  gridptr->laea_lat_0   = 0.0;
  gridptr->laea_a       = 0.0;
  gridptr->laea_defined = FALSE;
  gridptr->trunc        = 0;
  gridptr->nvertex      = 0;
  gridptr->nd           = 0;
  gridptr->ni           = 0;
  gridptr->ni2          = 0;
  gridptr->ni3          = 0;
  gridptr->number       = 0;
  gridptr->position     = 0;
  gridptr->reference    = NULL;
  gridptr->prec         = 0;
  gridptr->size         = 0;
  gridptr->xsize        = 0;
  gridptr->ysize        = 0;
  gridptr->np           = 0;
  gridptr->xdef         = 0;
  gridptr->ydef         = 0;
  gridptr->isCyclic     = CDI_UNDEFID;
  gridptr->isRotated    = FALSE;
  gridptr->xpole        = 0.0;
  gridptr->ypole        = 0.0;
  gridptr->angle        = 0.0;
  gridptr->locked       = FALSE;
  gridptr->lcomplex     = 0;
  gridptr->xname[0]     = 0;
  gridptr->yname[0]     = 0;
  gridptr->xlongname[0] = 0;
  gridptr->ylongname[0] = 0;
  gridptr->xunits[0]    = 0;
  gridptr->yunits[0]    = 0;
  gridptr->xstdname[0]  = 0;
  gridptr->ystdname[0]  = 0;
  memset(gridptr->uuid, 0, 16);
  gridptr->name         = NULL;
}


void grid_free(grid_t *gridptr)
{
  if ( gridptr->mask      ) free(gridptr->mask);
  if ( gridptr->mask_gme  ) free(gridptr->mask_gme);
  if ( gridptr->xvals     ) free(gridptr->xvals);
  if ( gridptr->yvals     ) free(gridptr->yvals);
  if ( gridptr->area      ) free(gridptr->area);
  if ( gridptr->xbounds   ) free(gridptr->xbounds);
  if ( gridptr->ybounds   ) free(gridptr->ybounds);
  if ( gridptr->rowlon    ) free(gridptr->rowlon);
  if ( gridptr->reference ) free(gridptr->reference);
  if ( gridptr->name      ) free(gridptr->name);

  grid_init(gridptr);
}

static
grid_t *gridNewEntry ( void )
{
  grid_t *gridptr;

  gridptr = ( grid_t *) xmalloc ( sizeof ( grid_t ));
  grid_init ( gridptr );
  gridptr->self = reshPut (( void * ) gridptr, &gridOps );

  return gridptr;
}

static
void gridInit (void)
{
  static int gridInitialized = 0;
  char *env;

  if ( gridInitialized ) return;

  gridInitialized = 1;

  env = getenv("GRID_DEBUG");
  if ( env ) GRID_Debug = atoi(env);
}

static
void grid_copy(grid_t *gridptr2, grid_t *gridptr1)
{
  int gridID2;

  gridID2 = gridptr2->self;
  memcpy(gridptr2, gridptr1, sizeof(grid_t));
  gridptr2->self = gridID2;
}

static
void gridCheckPtr(const char *caller, int gridID, grid_t *gridptr)
{
  if ( gridptr == NULL )
    Errorc("grid %d undefined!", gridID);
}

#define  grid_check_ptr(gridID, gridptr)  gridCheckPtr(__func__, gridID, gridptr)


int gridSize(void)
{
  return reshCountType ( &gridOps );
}


void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *xvals)
{
  if ( (! (fabs(xinc) > 0)) && xsize > 1 )
    {
      if ( xfirst >= xlast )
        {
          while ( xfirst >= xlast ) xlast += 360;
          xinc = (xlast-xfirst)/(xsize);
        }
      else
        {
          xinc = (xlast-xfirst)/(xsize-1);
        }
    }

  for ( int i = 0; i < xsize; ++i )
    xvals[i] = xfirst + i*xinc;
}

static
void calc_gaussgrid(double *yvals, int ysize, double yfirst, double ylast)
{
  double *yw;
  long yhsize;
  long i;

  yw = (double *) malloc(ysize*sizeof(double));
  gaussaw(yvals, yw, ysize);
  free(yw);
  for ( i = 0; i < ysize; i++ )
    yvals[i] = asin(yvals[i])/M_PI*180.0;

  if ( yfirst < ylast && yfirst > -90.0 && ylast < 90.0 )
    {
      double ytmp;
      yhsize = ysize/2;
      for ( i = 0; i < yhsize; i++ )
        {
          ytmp = yvals[i];
          yvals[i] = yvals[ysize-i-1];
          yvals[ysize-i-1] = ytmp;
        }
    }
}


void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *yvals)
{
  long i;
  double deleps = 0.002;

  if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
    {
      if ( ysize > 2 )
	{
	  calc_gaussgrid(yvals, ysize, yfirst, ylast);

	  if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
	    if ( fabs(yvals[0] - yfirst) > deleps || fabs(yvals[ysize-1] - ylast) > deleps )
	      {
		double yinc = fabs(ylast-yfirst)/(ysize-1);
		double *ytmp = NULL;
		int nstart, lfound = 0;
		int ny = (int) (180./yinc + 0.5);
		ny -= ny%2;
		/* printf("%g %g %g %g %g %d\n", ylast, yfirst, ylast-yfirst,yinc, 180/yinc, ny); */
		if ( ny > ysize && ny < 4096 )
		  {
		    ytmp = (double *) malloc(ny*sizeof(double));
		    calc_gaussgrid(ytmp, ny, yfirst, ylast);
		    for ( i = 0; i < (ny-ysize); i++ )
		      if ( fabs(ytmp[i] - yfirst) < deleps ) break;

		    nstart = i;

		    if ( (nstart+ysize-1) < ny )
		      if ( fabs(ytmp[nstart+ysize-1] - ylast) < deleps ) lfound = 1;
		  }

		if ( lfound )
		  {
		    for ( i = 0; i < ysize; i++ ) yvals[i] = ytmp[i+nstart];
		  }
		else
		  {
		    Warning("Cannot calculate gaussian latitudes for lat1 = %g latn = %g!", yfirst, ylast);
		    for ( i = 0; i < ysize; i++ ) yvals[i] = 0;
		    yvals[0] = yfirst;
		    yvals[ysize-1] = ylast;
		  }

		if ( ytmp ) free(ytmp);
	      }
	}
      else
        {
          yvals[0] = yfirst;
          yvals[ysize-1] = ylast;
        }
    }
  /*     else if ( gridtype == GRID_LONLAT || gridtype == GRID_GENERIC ) */
  else
    {
      if ( (! (fabs(yinc) > 0)) && ysize > 1 )
        {
          if ( IS_EQUAL(yfirst, ylast) && IS_NOT_EQUAL(yfirst, 0) ) ylast *= -1;

          if ( yfirst > ylast )
            yinc = (yfirst-ylast)/(ysize-1);
          else if ( yfirst < ylast )
            yinc = (ylast-yfirst)/(ysize-1);
          else
            {
              if ( ysize%2 != 0 )
                {
                  yinc = 180.0/(ysize-1);
                  yfirst = -90;
                }
              else
                {
                  yinc = 180.0/ysize;
                  yfirst = -90 + yinc/2;
                }
            }
        }

      if ( yfirst > ylast && yinc > 0 ) yinc = -yinc;

      for ( i = 0; i < ysize; i++ )
        yvals[i] = yfirst + i*yinc;
    }
  /*
    else
    Error("unable to calculate values for %s grid!", gridNamePtr(gridtype));
  */
}

/*
@Function  gridCreate
@Title     Create a horizontal Grid

@Prototype int gridCreate(int gridtype, int size)
@Parameter
    @Item  gridtype  The type of the grid, one of the set of predefined CDI grid types.
                     The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_GAUSSIAN},
                     @func{GRID_LONLAT}, @func{GRID_LCC}, @func{GRID_SPECTRAL},
                     @func{GRID_GME}, @func{GRID_CURVILINEAR} and @func{GRID_UNSTRUCTURED} and.
    @Item  size      Number of gridpoints.

@Description
The function @func{gridCreate} creates a horizontal Grid.

@Result
@func{gridCreate} returns an identifier to the Grid.

@Example
Here is an example using @func{gridCreate} to create a regular lon/lat Grid:

@Source
#include "cdi.h"
   ...
#define  nlon  12
#define  nlat   6
   ...
double lons[nlon] = {0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330};
double lats[nlat] = {-75, -45, -15, 15, 45, 75};
int gridID;
   ...
gridID = gridCreate(GRID_LONLAT, nlon*nlat);
gridDefXsize(gridID, nlon);
gridDefYsize(gridID, nlat);
gridDefXvals(gridID, lons);
gridDefYvals(gridID, lats);
   ...
@EndSource
@EndFunction
*/
int gridCreate(int gridtype, int size)
{
  int gridID;
  grid_t *gridptr;

  if ( CDI_Debug ) Message("gridtype=%s  size=%d", gridNamePtr(gridtype), size);

  if ( size < 0 || size > INT_MAX ) Error("Grid size (%d) out of bounds (0 - %d)!", size, INT_MAX);

  gridInit();

  gridptr = gridNewEntry();
  if ( ! gridptr ) Error("No memory");

  gridID = gridptr->self;

  if ( CDI_Debug ) Message("gridID: %d", gridID);

  gridptr->type = gridtype;
  gridptr->size = size;

  /*  if ( gridtype == GRID_GENERIC )     gridptr->xsize = size; */
  if ( gridtype == GRID_UNSTRUCTURED )  gridptr->xsize = size;
  if ( gridtype == GRID_CURVILINEAR  )  gridptr->nvertex = 4;

  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
    case GRID_GAUSSIAN_REDUCED:
    case GRID_CURVILINEAR:
    case GRID_TRAJECTORY:
      {
        if ( gridtype == GRID_TRAJECTORY )
          {
            gridDefXname(gridID, "tlon");
            gridDefYname(gridID, "tlat");
          }
        else
          {
            gridDefXname(gridID, "lon");
            gridDefYname(gridID, "lat");
          }
        gridDefXlongname(gridID, "longitude");
        gridDefYlongname(gridID, "latitude");

        /*
        if ( gridtype == GRID_CURVILINEAR )
          {
            strcpy(gridptr->xstdname, "grid_longitude");
            strcpy(gridptr->ystdname, "grid_latitude");
            gridDefXunits(gridID, "degrees");
            gridDefYunits(gridID, "degrees");
          }
        else
        */
          {
            strcpy(gridptr->xstdname, "longitude");
            strcpy(gridptr->ystdname, "latitude");
            gridDefXunits(gridID, "degrees_east");
            gridDefYunits(gridID, "degrees_north");
          }

        break;
      }
    case GRID_GME:
    case GRID_UNSTRUCTURED:
      {
        gridDefXname(gridID, "lon");
        gridDefYname(gridID, "lat");
        strcpy(gridptr->xstdname, "longitude");
        strcpy(gridptr->ystdname, "latitude");
        gridDefXunits(gridID, "degrees_east");
        gridDefYunits(gridID, "degrees_north");
        break;
      }
    case GRID_GENERIC:
      {
        gridDefXname(gridID, "x");
        gridDefYname(gridID, "y");
        strcpy(gridptr->xstdname, "grid_longitude");
        strcpy(gridptr->ystdname, "grid_latitude");
        gridDefXunits(gridID, "degrees");
        gridDefYunits(gridID, "degrees");
        break;
      }
    case GRID_LCC2:
    case GRID_SINUSOIDAL:
    case GRID_LAEA:
      {
        gridDefXname(gridID, "x");
        gridDefYname(gridID, "y");
        strcpy(gridptr->xstdname, "projection_x_coordinate");
        strcpy(gridptr->ystdname, "projection_y_coordinate");
        gridDefXunits(gridID, "m");
        gridDefYunits(gridID, "m");
        break;
      }
    }

  return (gridID);
}

static
void gridDestroyKernel( grid_t * gridptr )
{
  int id;

  xassert ( gridptr );

  id = gridptr->self;

  if ( gridptr->mask      ) free(gridptr->mask);
  if ( gridptr->mask_gme  ) free(gridptr->mask_gme);
  if ( gridptr->xvals     ) free(gridptr->xvals);
  if ( gridptr->yvals     ) free(gridptr->yvals);
  if ( gridptr->area      ) free(gridptr->area);
  if ( gridptr->xbounds   ) free(gridptr->xbounds);
  if ( gridptr->ybounds   ) free(gridptr->ybounds);
  if ( gridptr->rowlon    ) free(gridptr->rowlon);
  if ( gridptr->reference ) free(gridptr->reference);

  free ( gridptr );

  reshRemove ( id, &gridOps );
}

/*
@Function  gridDestroy
@Title     Destroy a horizontal Grid

@Prototype void gridDestroy(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.

@EndFunction
*/
void gridDestroy(int gridID)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  gridDestroyKernel ( gridptr );
}

void gridDestroyP ( void * gridptr )
{
  gridDestroyKernel (( grid_t * ) gridptr );
}


char *gridNamePtr(int gridtype)
{
  char *name;
  int size = (int) (sizeof(Grids)/sizeof(char *));

  if ( gridtype >= 0 && gridtype < size )
    name = Grids[gridtype];
  else
    name = Grids[GRID_GENERIC];

  return (name);
}


void gridName(int gridtype, char *gridname)
{
  strcpy(gridname, gridNamePtr(gridtype));
}

/*
@Function  gridDefXname
@Title     Define the name of a X-axis

@Prototype void gridDefXname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the X-axis.

@Description
The function @func{gridDefXname} defines the name of a X-axis.

@EndFunction
*/
void gridDefXname(int gridID, const char *xname)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning ("%s", "Operation not executed." );
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( xname )
    strcpy(gridptr->xname, xname);
}

/*
@Function  gridDefXlongname
@Title     Define the longname of a X-axis

@Prototype void gridDefXlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the X-axis.

@Description
The function @func{gridDefXlongname} defines the longname of a X-axis.

@EndFunction
*/
void gridDefXlongname(int gridID, const char *xlongname)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning ("%s", "Operation not executed." );
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  if ( xlongname )
    strcpy(gridptr->xlongname, xlongname);
}

/*
@Function  gridDefXunits
@Title     Define the units of a X-axis

@Prototype void gridDefXunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the X-axis.

@Description
The function @func{gridDefXunits} defines the units of a X-axis.

@EndFunction
*/
void gridDefXunits(int gridID, const char *xunits)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed." );
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( xunits )
    strcpy(gridptr->xunits, xunits);
}

/*
@Function  gridDefYname
@Title     Define the name of a Y-axis

@Prototype void gridDefYname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the Y-axis.

@Description
The function @func{gridDefYname} defines the name of a Y-axis.

@EndFunction
*/
void gridDefYname(int gridID, const char *yname)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed.");
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( yname )
    strcpy(gridptr->yname, yname);
}

/*
@Function  gridDefYlongname
@Title     Define the longname of a Y-axis

@Prototype void gridDefYlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the Y-axis.

@Description
The function @func{gridDefYlongname} defines the longname of a Y-axis.

@EndFunction
*/
void gridDefYlongname(int gridID, const char *ylongname)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed.");
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( ylongname )
    strcpy(gridptr->ylongname, ylongname);
}

/*
@Function  gridDefYunits
@Title     Define the units of a Y-axis

@Prototype void gridDefYunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the Y-axis.

@Description
The function @func{gridDefYunits} defines the units of a Y-axis.

@EndFunction
*/
void gridDefYunits(int gridID, const char *yunits)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed.");
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( yunits )
    strcpy(gridptr->yunits, yunits);
}

/*
@Function  gridInqXname
@Title     Get the name of a X-axis

@Prototype void gridInqXname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  name     Name of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXname} returns the name of a X-axis.

@Result
@func{gridInqXname} returns the name of the X-axis to the parameter name.

@EndFunction
*/
void gridInqXname(int gridID, char *xname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(xname, gridptr->xname);
}

/*
@Function  gridInqXlongname
@Title     Get the longname of a X-axis

@Prototype void gridInqXlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  longname Longname of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXlongname} returns the longname of a X-axis.

@Result
@func{gridInqXlongname} returns the longname of the X-axis to the parameter longname.

@EndFunction
*/
void gridInqXlongname(int gridID, char *xlongname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(xlongname, gridptr->xlongname);
}

/*
@Function  gridInqXunits
@Title     Get the units of a X-axis

@Prototype void gridInqXunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  units    Units of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXunits} returns the units of a X-axis.

@Result
@func{gridInqXunits} returns the units of the X-axis to the parameter units.

@EndFunction
*/
void gridInqXunits(int gridID, char *xunits)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(xunits, gridptr->xunits);
}


void gridInqXstdname(int gridID, char *xstdname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(xstdname, gridptr->xstdname);
}

/*
@Function  gridInqYname
@Title     Get the name of a Y-axis

@Prototype void gridInqYname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  name     Name of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYname} returns the name of a Y-axis.

@Result
@func{gridInqYname} returns the name of the Y-axis to the parameter name.

@EndFunction
*/
void gridInqYname(int gridID, char *yname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(yname, gridptr->yname);
}

/*
@Function  gridInqYlongname
@Title     Get the longname of a Y-axis

@Prototype void gridInqXlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  longname Longname of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYlongname} returns the longname of a Y-axis.

@Result
@func{gridInqYlongname} returns the longname of the Y-axis to the parameter longname.

@EndFunction
*/
void gridInqYlongname(int gridID, char *ylongname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(ylongname, gridptr->ylongname);
}

/*
@Function  gridInqYunits
@Title     Get the units of a Y-axis

@Prototype void gridInqYunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  units    Units of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYunits} returns the units of a Y-axis.

@Result
@func{gridInqYunits} returns the units of the Y-axis to the parameter units.

@EndFunction
*/
void gridInqYunits(int gridID, char *yunits)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(yunits, gridptr->yunits);
}

void gridInqYstdname(int gridID, char *ystdname)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  strcpy(ystdname, gridptr->ystdname);
}

/*
@Function  gridInqType
@Title     Get the type of a Grid

@Prototype int gridInqType(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.

@Description
The function @func{gridInqType} returns the type of a Grid.

@Result
@func{gridInqType} returns the type of the grid,
one of the set of predefined CDI grid types.
The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_GAUSSIAN},
@func{GRID_LONLAT}, @func{GRID_LCC}, @func{GRID_SPECTRAL}, @func{GRID_GME},
@func{GRID_CURVILINEAR} and @func{GRID_UNSTRUCTURED}.

@EndFunction
*/
int gridInqType(int gridID)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  return (gridptr->type);
}


/*
@Function  gridInqSize
@Title     Get the size of a Grid

@Prototype int gridInqSize(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.

@Description
The function @func{gridInqSize} returns the size of a Grid.

@Result
@func{gridInqSize} returns the number of grid points of a Grid.

@EndFunction
*/
int gridInqSize(int gridID)
{
  int size = 0;
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  size = gridptr->size;

  if ( ! size )
    {
      int xsize, ysize;

      xsize = gridptr->xsize;
      ysize = gridptr->ysize;

      if ( ysize )
        size = xsize *ysize;
      else
        size = xsize;

      gridptr->size = size;
    }

  return (size);
}

static
int nsp2trunc(int nsp)
{
  /*  nsp = (trunc+1)*(trunc+1)              */
  /*      => trunc^2 + 3*trunc - (x-2) = 0   */
  /*                                         */
  /*  with:  y^2 + p*y + q = 0               */
  /*         y = -p/2 +- sqrt((p/2)^2 - q)   */
  /*         p = 3 and q = - (x-2)           */
  int trunc;

  trunc = (int) (sqrt(nsp*4 + 1.) - 3) / 2;

  return (trunc);
}


int gridInqTrunc(int gridID)
{
  grid_t *gridptr;

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  if ( gridptr->trunc == 0 )
    {
      if ( gridptr->type == GRID_SPECTRAL )
        gridptr->trunc = nsp2trunc(gridptr->size);
      /*
      else if      ( gridptr->type == GRID_GAUSSIAN )
        gridptr->trunc = nlat2trunc(gridptr->ysize);
      */
    }

  return (gridptr->trunc);
}


void gridDefTrunc(int gridID, int trunc)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed.");
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );

  grid_check_ptr(gridID, gridptr);

  gridptr->trunc = trunc;
}

/*
@Function  gridDefXsize
@Title     Define the number of values of a X-axis

@Prototype void gridDefXsize(int gridID, int xsize)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  xsize    Number of values of a X-axis.

@Description
The function @func{gridDefXsize} defines the number of values of a X-axis.

@EndFunction
*/
void gridDefXsize(int gridID, int xsize)
{
  grid_t *gridptr;

  if ( reshGetStatus ( gridID, &gridOps ) == CLOSED )
    {
      Warning("%s", "Operation not executed.");
      return;
    }

  gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );