grid.c

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

#include <string.h>

#include "dmemory.h"
#include "cdi.h"
#include "stream_int.h"
#include "grid.h"
#include "util.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 */  "cell",
  /* 10 */  "curvilinear",
  /* 11 */  "lcc",
  /* 12 */  "lcc2",
  /* 13 */  "laea",
  /* 14 */  "sinusoidal",
};


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

static int _grid_max = MAX_GRIDS;

static void grid_initialize(void);

static int _grid_init = FALSE;

#if  defined  (HAVE_LIBPTHREAD)
#  include <pthread.h>

static pthread_once_t _grid_init_thread = PTHREAD_ONCE_INIT;
static pthread_mutex_t _grid_mutex;

#  define GRID_LOCK           pthread_mutex_lock(&_grid_mutex);
#  define GRID_UNLOCK         pthread_mutex_unlock(&_grid_mutex);
#  define GRID_INIT                               \
   if ( _grid_init == FALSE ) pthread_once(&_grid_init_thread, grid_initialize);

#else

#  define GRID_LOCK
#  define GRID_UNLOCK
#  define GRID_INIT                               \
   if ( _grid_init == FALSE ) grid_initialize();

#endif


typedef struct _gridPtrToIdx {
  int idx;
  grid_t *ptr;
  struct _gridPtrToIdx *next;
} gridPtrToIdx;


static gridPtrToIdx *_gridList  = NULL;
static gridPtrToIdx *_gridAvail = NULL;


static void grid_list_new(void)
{
  static char func[] = "grid_list_new";

  assert(_gridList == NULL);

  _gridList = (gridPtrToIdx *) malloc(_grid_max*sizeof(gridPtrToIdx));
}


static void grid_list_delete(void)
{
  static char func[] = "grid_list_delete";

  if ( _gridList ) free(_gridList);
}


static void grid_init_pointer(void)
{
  int  i;
  
  for ( i = 0; i < _grid_max; i++ )
    {
      _gridList[i].next = _gridList + i + 1;
      _gridList[i].idx  = i;
      _gridList[i].ptr  = 0;
    }

  _gridList[_grid_max-1].next = 0;

  _gridAvail = _gridList;
}


grid_t *grid_to_pointer(int idx)
{
  static char func[] = "grid_to_pointer";
  grid_t *gridptr = NULL;

  GRID_INIT

  if ( idx >= 0 && idx < _grid_max )
    {
      GRID_LOCK

      gridptr = _gridList[idx].ptr;

      GRID_UNLOCK
    }
  else
    Error(func, "grid index %d undefined!", idx);

  return (gridptr);
}


/* Create an index from a pointer */
static int grid_from_pointer(grid_t *ptr)
{
  static char func[] = "grid_from_pointer";
  int      idx = -1;
  gridPtrToIdx *newptr;

  if ( ptr )
    {
      GRID_LOCK

      if ( _gridAvail )
	{
	  newptr       = _gridAvail;
	  _gridAvail   = _gridAvail->next;
	  newptr->next = 0;
	  idx	       = newptr->idx;
	  newptr->ptr  = ptr;
      
	  if ( GRID_Debug )
	    Message(func, "Pointer %p has idx %d from grid list", ptr, idx);
	}
      else
	Warning(func, "Too many open grids (limit is %d)!", _grid_max);

      GRID_UNLOCK
    }
  else
    Error(func, "Internal problem (pointer %p undefined)", ptr);

  return (idx);
}


static void grid_init_entry(grid_t *gridptr)
{
  gridptr->self         = grid_from_pointer(gridptr);

  gridptr->type         = CDI_UNDEFID;
  gridptr->mask         = 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->prec         = 0;
  gridptr->size         = 0;
  gridptr->xsize        = 0;
  gridptr->ysize        = 0;
  gridptr->isCyclic     = CDI_UNDEFID;
  gridptr->isRotated    = FALSE;
  gridptr->xpole        = 0.0;
  gridptr->ypole        = 0.0;
  gridptr->angle        = 0.0;
  gridptr->locked       = FALSE;
  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;
}


static grid_t *grid_new_entry(void)
{
  static char func[] = "grid_new_entry";
  grid_t *gridptr;

  gridptr = (grid_t *) malloc(sizeof(grid_t));

  if ( gridptr ) grid_init_entry(gridptr);

  return (gridptr);
}


static void grid_delete_entry(grid_t *gridptr)
{
  static char func[] = "grid_delete_entry";
  int idx;

  idx = gridptr->self;

  GRID_LOCK

  free(gridptr);

  _gridList[idx].next = _gridAvail;
  _gridList[idx].ptr  = 0;
  _gridAvail          = &_gridList[idx];

  GRID_UNLOCK

  if ( GRID_Debug )
    Message(func, "Removed idx %d from grid list", idx);
}


static void grid_initialize(void)
{
  char *env;

#if  defined  (HAVE_LIBPTHREAD)
  /* initialize global API mutex lock */
  pthread_mutex_init(&_grid_mutex, NULL);
#endif

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

  grid_list_new();
  atexit(grid_list_delete);

  GRID_LOCK

  grid_init_pointer();

  GRID_UNLOCK

  _grid_init = TRUE;
}


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 grid_check_ptr(const char *func, grid_t *gridptr)
{
  if ( gridptr == NULL )
    Error(func, "grid undefined!");
}


int gridSize(void)
{
  int gridsize = 0;
  int i;
  
  GRID_INIT

  GRID_LOCK

  for ( i = 0; i < _grid_max; i++ )
    if ( _gridList[i].ptr ) gridsize++;

  GRID_UNLOCK

  return (gridsize);
}


void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *xvals)
{
  int i;

  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 ( i = 0; i < xsize; i++ )
    xvals[i] = xfirst + i*xinc;
}


void calc_gaussaw(double *yvals, int ysize, double yfirst, double ylast)
{
  static char func[] = "calc_gaussaw";
  double *yw;
  int yhsize;
  int 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)
{
  static char func[] = "gridGenYvals";
  int i;

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

	  if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
	    if ( fabs(yvals[0] - yfirst) > 0.001 || fabs(yvals[ysize-1] - ylast) > 0.001 )
	      {
		double yinc = fabs(ylast-yfirst)/(ysize-1);
		double *ytmp;
		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); */
		ytmp   = (double *) malloc(ny*sizeof(double));
		calc_gaussaw(ytmp, ny, yfirst, ylast);
		for ( i = 0; i < (ny-ysize); i++ )
		  if ( fabs(ytmp[i] - yfirst) < 0.001 ) break;

		nstart = i;

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

		if ( lfound )
		  {
		    for ( i = 0; i < ysize; i++ ) yvals[i] = ytmp[i+nstart];
		  }
		else
		  {
		    Warning(func, "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;
		  }

		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(func, "unable to calculate values for %s grid!", gridNamePtr(gridtype));
  */
}
/*
static void defineXvals(int gridID)
{
  static char func[] = "defineXvals";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  if ( gridptr->type == GRID_CELL ) return;

  if ( gridptr->xvals == NULL )
    {
      double *xvals, xfirst, xlast, xinc;
      int i, xsize;

      xsize = gridptr->xsize;
      if ( xsize == 0 )
	Error(func, "xsize undefined for gridID = %d", gridID);

      xvals = (double *) malloc(xsize*sizeof(double));
      xfirst = gridptr->xfirst;
      xlast  = gridptr->xlast;
      xinc = gridptr->xinc;
      if ( ! fabs(xinc) )
	{
	  if ( xfirst > xlast )
	    xinc = (xfirst-xlast)/xsize;
	  else
	    xinc = (xlast-xfirst)/xsize;

	  gridptr->xinc = xinc;
	}
      if ( xfirst > xlast && xinc > 0 ) xinc = -xinc;
      for ( i = 0; i < xsize; i++ )
	xvals[i] = xfirst + i*xinc;

      gridptr->xvals = xvals;
    }
}
*/

static void defineYvals(int gridID)
{
  static char func[] = "defineYvals";
  int gridtype;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  gridtype = gridptr->type;

  if ( gridtype == GRID_CELL ) return;

  if ( gridptr->yvals == NULL )
    {
      double *yvals, yfirst, ylast, yinc;
      int i, ysize;

      ysize = gridptr->ysize;

      if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
	{
	  double *yw;
	  int yhsize;

	  yvals = (double *) malloc(ysize*sizeof(double));
	  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;

	  yfirst = gridptr->yfirst;
	  ylast  = gridptr->ylast;
	  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;
		}
	    }
	  gridptr->yvals = yvals;
	}
      /*     else if ( gridtype == GRID_LONLAT || gridtype == GRID_GENERIC ) */
      else
	{
	  if ( ysize == 0 )
	    Error(func, "ysize undefined for gridID = %d", gridID);

	  yvals = (double *) malloc(ysize*sizeof(double));
	  yfirst = gridptr->yfirst;
	  ylast  = gridptr->ylast;
	  yinc = gridptr->yinc;
	  if ( ! (fabs(yinc) > 0) )
	    {
	      if ( ysize > 1 )
		{
		  if ( yfirst > ylast )
		    yinc = (yfirst-ylast)/(ysize-1);
		  else
		    yinc = (ylast-yfirst)/(ysize-1);
		}

	      gridptr->yinc = yinc;
	    }
	  if ( yfirst > ylast && yinc > 0 ) yinc = -yinc;
	  for ( i = 0; i < ysize; i++ )
	    yvals[i] = yfirst + i*yinc;

	  gridptr->yvals = yvals;
	}
      /*
      else
	Error(func, "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_CELL}.
    @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)
{
  static char func[] = "gridCreate";
  int gridID;
  grid_t *gridptr;

  if ( CDI_Debug )
    Message(func, "gridtype: %d size: %d", gridtype, size);

  GRID_INIT

  gridptr = grid_new_entry();
  if ( ! gridptr ) Error(func, "No memory");

  gridID = gridptr->self;

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

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

  /*  if ( gridtype == GRID_GENERIC )     gridptr->xsize = size; */
  if ( gridtype == GRID_CELL )        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_CELL:
      {
	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);
}


/*
@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)
{
}


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, char *xname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  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, char *xlongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  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, char *xunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  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, char *yname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  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, char *ylongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  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, char *yunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

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


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

@Prototype void gridInqXname(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{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_to_pointer(gridID);

  strcpy(xname, gridptr->xname);
}


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

@Prototype void gridInqXlongname(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{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_to_pointer(gridID);

  strcpy(xlongname, gridptr->xlongname);
}


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

@Prototype void gridInqXunits(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{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_to_pointer(gridID);

  strcpy(xunits, gridptr->xunits);
}


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

  gridptr = grid_to_pointer(gridID);

  strcpy(xstdname, gridptr->xstdname);
}


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

@Prototype void gridInqYname(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{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_to_pointer(gridID);

  strcpy(yname, gridptr->yname);
}


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

@Prototype void gridInqXlongname(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{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_to_pointer(gridID);

  strcpy(ylongname, gridptr->ylongname);
}


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

@Prototype void gridInqYunits(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{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_to_pointer(gridID);

  strcpy(yunits, gridptr->yunits);
}

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

  gridptr = grid_to_pointer(gridID);

  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}

@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_CELL}.

@EndFunction
*/
int gridInqType(int gridID)
{
  static char func[] = "gridInqType";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, 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}

@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)
{
  static char func[] = "gridInqSize";
  int size = 0;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, 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.0 + 1) - 3) / 2;

  return (trunc);
}


int gridInqTrunc(int gridID)
{
  static char func[] = "gridInqTrunc";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, 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)
{
  static char func[] = "gridDefTrunc";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, 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)
{
  static char func[] = "gridDefXsize";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  if ( xsize > gridInqSize(gridID) )
    Error(func, "xsize %d is greater then gridsize %d", xsize, gridInqSize(gridID));

  if ( gridInqType(gridID) == GRID_CELL && xsize != gridInqSize(gridID) )
    Error(func, "xsize %d must be equal gridsize %d for gridtype CELL", xsize, gridInqSize(gridID));

  gridptr->xsize = xsize;

  if ( gridInqType(gridID) != GRID_CELL && 
       gridptr->xsize*gridptr->ysize > gridInqSize(gridID) )
    Error(func, "inconsistent grid declaration! (xsize %d ysize %d gridsize %d)",
	  gridptr->xsize, gridptr->ysize, gridInqSize(gridID));
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
void gridDefPrec(int gridID, int prec)
{
  static char func[] = "gridDefPrec";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  gridptr->prec = prec;
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
int gridInqPrec(int gridID)
{
  static char func[] = "gridInqPrec";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  return (gridptr->prec);
}


/*
@Function  gridInqXsize
@Title     Get the number of values of a X-axis

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

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

@Result
@func{gridInqXsize} returns the number of values of a X-axis.

@EndFunction
*/
int gridInqXsize(int gridID)
{
  static char func[] = "gridInqXsize";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  return (gridptr->xsize);
}


/*
@Function  gridDefYsize
@Title     Define the number of values of a Y-axis

@Prototype void gridDefYsize(int gridID, int ysize)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}
    @Item  ysize    Number of values of a Y-axis

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

@EndFunction
*/
void gridDefYsize(int gridID, int ysize)
{
  static char func[] = "gridDefYsize";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  if ( ysize > gridInqSize(gridID) )
    Error(func, "ysize %d is greater then gridsize %d", ysize, gridInqSize(gridID));

  if ( gridInqType(gridID) == GRID_CELL && ysize != gridInqSize(gridID) )
    Error(func, "ysize %d must be equal gridsize %d for gridtype CELL", ysize, gridInqSize(gridID));

  gridptr->ysize = ysize;

  if ( gridInqType(gridID) != GRID_CELL && 
       gridptr->xsize*gridptr->ysize > gridInqSize(gridID) )
    Error(func, "inconsistent grid declaration! (xsize %d ysize %d gridsize %d)",
	  gridptr->xsize, gridptr->ysize, gridInqSize(gridID));
}


/*
@Function  gridInqYsize
@Title     Get the number of values of a Y-axis

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

@Description
The function @func{gridInqYsize} returns the number of values of a Y-axis.

@Result
@func{gridInqYsize} returns the number of values of a Y-axis.

@EndFunction
*/
int gridInqYsize(int gridID)
{
  static char func[] = "gridInqYsize";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  return (gridptr->ysize);
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
void gridDefRowlon(int gridID, int nrowlon, int *rowlon)
{
  static char func[] = "gridDefRowlon";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);
  
  grid_check_ptr(func, gridptr);

  gridptr->rowlon = (int *) malloc(nrowlon*sizeof(int));
  gridptr->nrowlon = nrowlon;

  memcpy(gridptr->rowlon, rowlon, nrowlon*sizeof(int));
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
void gridInqRowlon(int gridID, int *rowlon)
{
  static char func[] = "gridInqRowlon";
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  if ( gridptr->rowlon == 0 )  Error(func, "undefined pointer");

  memcpy(rowlon, gridptr->rowlon, gridptr->nrowlon*sizeof(int));
}


int gridInqMask(int gridID, int *mask)
{
  static char func[] = "gridInqMask";
  int size;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  size = gridptr->size;

  if ( CDI_Debug && size == 0 )
    Warning(func, "Size undefined for gridID = %d", gridID);
    
  if ( mask && gridptr->mask )
    for ( i = 0; i < size; i++ )
      mask[i] =  gridptr->mask[i];

  if ( gridptr->mask == NULL ) size = 0;

  return (size);
}


void gridDefMask(int gridID, int *mask)
{
  static char func[] = "gridDefMask";
  int size;
  int gridtype;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  gridtype = gridptr->type;

  size = gridptr->size;

  if ( size == 0 )
    Error(func, "Size undefined for gridID = %d", gridID);
    
  if ( gridptr->mask == NULL )
    gridptr->mask = (int *) malloc(size*sizeof(int));
  else if ( CDI_Debug )
    Warning(func, "mask already defined!");

  for ( i = 0; i < size; i++ )
    gridptr->mask[i] = mask[i];
}


/*
@Function  gridInqXvals
@Title     Get all values of a X-axis

@Prototype int gridInqXvals(int gridID, double *xvals)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}
    @Item  xvals    X-values of the grid

@Description
The function @func{gridInqXvals} returns all values of the X-axis.

@Result
Upon successful completion @func{gridInqXvals} returns the number of values and
the values are stored in @func{xvals}.
Otherwise, 0 is returned and @func{xvals} is empty.

@EndFunction
*/
int gridInqXvals(int gridID, double *xvals)
{
  static char func[] = "gridInqXvals";
  int size;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  if ( gridptr->type == GRID_CURVILINEAR || gridptr->type == GRID_CELL )
    size = gridptr->size;
  else
    size = gridptr->xsize;

  if ( CDI_Debug && size == 0 )
    Warning(func, "Size undefined for gridID = %d", gridID);
    
  if ( xvals && gridptr->xvals )
    for ( i = 0; i < size; i++ )
      xvals[i] = gridptr->xvals[i];

  if ( gridptr->xvals == NULL ) size = 0;

  return (size);
}


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

@Prototype void gridDefXvals(int gridID, double *xvals)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}
    @Item  xvals    X-values of the grid

@Description
The function @func{gridDefXvals} defines all values of the X-axis.

@EndFunction
*/
void gridDefXvals(int gridID, double *xvals)
{
  static char func[] = "gridDefXvals";
  int size;
  int gridtype;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  gridtype = gridptr->type;

  if ( gridtype == GRID_CELL || gridtype == GRID_CURVILINEAR )
    size = gridptr->size;
  else
    size = gridptr->xsize;

  if ( size == 0 )
    Error(func, "Size undefined for gridID = %d", gridID);
    
  if ( gridptr->xvals == NULL )
    gridptr->xvals = (double *) malloc(size*sizeof(double));
  else if ( CDI_Debug )
    Warning(func, "values already defined!");

  for ( i = 0; i < size; i++ )
    gridptr->xvals[i] = xvals[i];
}


/*
@Function  gridInqYvals
@Title     Get all values of a Y-axis

@Prototype int gridInqYvals(int gridID, double *yvals)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}
    @Item  yvals    Y-values of the grid

@Description
The function @func{gridInqYvals} returns all values of the Y-axis.

@Result
Upon successful completion @func{gridInqYvals} returns the number of values and
the values are stored in @func{yvals}.
Otherwise, 0 is returned and @func{yvals} is empty.

@EndFunction
*/
int gridInqYvals(int gridID, double *yvals)
{
  static char func[] = "gridInqYvals";
  int size;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  if ( gridptr->type == GRID_CURVILINEAR || gridptr->type == GRID_CELL )
    size = gridptr->size;
  else
    size = gridptr->ysize;

  if ( CDI_Debug && size == 0 )
    Warning(func, "Size undefined for gridID = %d!", gridID);
    
  if ( yvals && gridptr->yvals )
    for ( i = 0; i < size; i++ )
      yvals[i] =  gridptr->yvals[i];

  if ( gridptr->yvals == NULL ) size = 0;

  return (size);
}


/*
@Function  gridDefYvals
@Title     Define the values of a Y-axis

@Prototype void gridDefYvals(int gridID, double *yvals)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}
    @Item  yvals    Y-values of the grid

@Description
The function @func{gridDefYvals} defines all values of the Y-axis.

@EndFunction
*/
void gridDefYvals(int gridID, double *yvals)
{
  static char func[] = "gridDefYvals";
  int size;
  int gridtype;
  int i;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(func, gridptr);

  gridtype = gridptr->type;

  if ( gridtype == GRID_CELL || gridtype == GRID_CURVILINEAR )
    size = gridptr->size;
  else
    size = gridptr->ysize;

  if ( size == 0 )
    Error(func, "Size undefined for gridID = %d!", gridID);
    
  if ( gridptr->yvals == NULL )
    gridptr->yvals = (double *) malloc(size*sizeof(double));
  else if ( CDI_Debug )
    Warning(func, "Values already defined!");
    
  for ( i = 0; i < size; i++ )
    gridptr->yvals[i] = yvals[i];
}


double gridInqXval(int gridID, int index)
{
  double xval = 0;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  if ( gridptr->xvals )
    xval = gridptr->xvals[index];

  return (xval);
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
double gridInqYval(int gridID, int index)
{
  double yval = 0;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  if ( gridptr->yvals )
    yval = gridptr->yvals[index];

  return (yval);
}


/*</