grid.c

/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2016 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

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

#if defined(_OPENMP)
#  include <omp.h>
#endif

#include <stdio.h>
#include <stdarg.h> /* va_list */

#if defined(HAVE_LIBPROJ)
#  include "proj_api.h"
#endif

#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "error.h"
#include "grid.h"


static
void scale_vec(double scalefactor, long nvals, double *restrict values)
{
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(nvals, scalefactor, values)
#endif
  for ( long n = 0; n < nvals; ++n )
    {
      values[n] *= scalefactor;
    }
}


void grid_to_radian(const char *units, long nvals, double *restrict values, const char *description)
{
  if ( cmpstr(units, "degree") == 0 )
    {
      scale_vec(DEG2RAD, nvals, values);
    }
  else if ( cmpstr(units, "radian") == 0 )
    {
      /* No conversion necessary */
    }
  else
    {
      cdoWarning("Unknown units [%s] supplied for %s; proceeding assuming radians!", units, description);
    }
}


void grid_to_degree(const char *units, long nvals, double *restrict values, const char *description)
{
  if ( cmpstr(units, "radian") == 0 )
    {
      for ( long n = 0; n < nvals; ++n ) values[n] *= RAD2DEG;
    }
  else if ( cmpstr(units, "degree") == 0 )
    {
      /* No conversion necessary */
    }
  else
    {
      cdoWarning("Unknown units [%s] supplied for %s; proceeding assuming degress!", units, description);
    }
}


int gridToZonal(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
  double  xval = 0;
  double *yvals;

  gridtype = gridInqType(gridID1);
  gridsize = gridInqYsize(gridID1);
  gridID2  = gridCreate(gridtype, gridsize);
	  
  if ( gridtype == GRID_LONLAT   ||
       gridtype == GRID_GAUSSIAN ||
       gridtype == GRID_GENERIC )
    {
      gridDefXsize(gridID2, 1);
      gridDefYsize(gridID2, gridsize);

      gridDefXvals(gridID2, &xval);

      if ( gridInqYvals(gridID1, NULL) )
	{
	  yvals = (double*) Malloc(gridsize*sizeof(double));

	  gridInqYvals(gridID1, yvals);
	  gridDefYvals(gridID2, yvals);

	  Free(yvals);
	}
    }
  else
    {
      Error("Gridtype %s unsupported!", gridNamePtr(gridtype));
    }

  return gridID2;
}


int gridToMeridional(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
  double *xvals;
  double  yval = 0;

  gridtype = gridInqType(gridID1);
  gridsize = gridInqXsize(gridID1);
  gridID2  = gridCreate(gridtype, gridsize);
	  
  if ( gridtype == GRID_LONLAT   ||
       gridtype == GRID_GAUSSIAN ||
       gridtype == GRID_GENERIC )
    {
      gridDefXsize(gridID2, gridsize);
      gridDefYsize(gridID2, 1);

      if ( gridInqXvals(gridID1, NULL) )
	{
	  xvals = (double*) Malloc(gridsize*sizeof(double));

	  gridInqXvals(gridID1, xvals);
	  gridDefXvals(gridID2, xvals);

	  Free(xvals);
	}

      gridDefYvals(gridID2, &yval);
    }
  else
    {
      Error("Gridtype %s unsupported!", gridNamePtr(gridtype));
    }

  return gridID2;
}


void grid_gen_corners(int n, const double* restrict vals, double* restrict corners)
{
  int i;

  if ( n == 1 )
    {
      corners[0] = vals[0];
      corners[1] = vals[0];
    }
  else
    {
      for ( i = 0; i < n-1; ++i )
	corners[i+1] = 0.5*(vals[i] + vals[i+1]);

      corners[0] = 2*vals[0] - corners[1];
      corners[n] = 2*vals[n-1] - corners[n-1];
    }
}


void grid_gen_bounds(int n, const double *restrict vals, double *restrict bounds)
{
  for ( int i = 0; i < n-1; ++i )
    {
      bounds[2*i+1]   = 0.5*(vals[i] + vals[i+1]);
      bounds[2*(i+1)] = 0.5*(vals[i] + vals[i+1]);
    }

  bounds[0]     = 2*vals[0] - bounds[1];
  bounds[2*n-1] = 2*vals[n-1] - bounds[2*(n-1)];
}


void grid_check_lat_borders(int n, double *ybounds)
{
  if ( ybounds[0] > ybounds[n-1] )
    {
      if ( ybounds[0]   >  88 ) ybounds[0]   =  90;
      if ( ybounds[n-1] < -88 ) ybounds[n-1] = -90;
    }
  else
    {
      if ( ybounds[0]   < -88 ) ybounds[0]   = -90;
      if ( ybounds[n-1] >  88 ) ybounds[n-1] =  90;
    }
}


void grid_cell_center_to_bounds_X2D(const char* xunitstr, long xsize, long ysize, const double* restrict grid_center_lon, 
				    double* restrict grid_corner_lon, double dlon)
{
  long i, j, index;
  double minlon, maxlon;

  if ( ! (dlon > 0) ) dlon = 360./xsize;
  /*
  if ( xsize == 1 || (grid_center_lon[xsize-1]-grid_center_lon[0]+dlon) < 359 )
    cdoAbort("Cannot calculate Xbounds for %d vals with dlon = %g", xsize, dlon);
  */
  for ( i = 0; i < xsize; ++i )
    {
      minlon = grid_center_lon[i] - 0.5*dlon;
      maxlon = grid_center_lon[i] + 0.5*dlon;
      for ( j = 0; j < ysize; ++j )
	{
	  index = (j<<2)*xsize + (i<<2);
	  grid_corner_lon[index  ] = minlon;
	  grid_corner_lon[index+1] = maxlon;
	  grid_corner_lon[index+2] = maxlon;
	  grid_corner_lon[index+3] = minlon;
	}
    }
}

static
double genYmin(double y1, double y2)
{
  double ymin, dy;

  dy = y2 - y1;
  ymin = y1 - dy/2;

  if ( y1 < -85 && ymin < -87.5 ) ymin = -90;

  if ( cdoVerbose )
    cdoPrint("genYmin: y1 = %g  y2 = %g  dy = %g  ymin = %g", y1, y2, dy, ymin);

  return ymin;
}

static
double genYmax(double y1, double y2)
{
  double ymax, dy;

  dy = y1 - y2;
  ymax = y1 + dy/2;

  if ( y1 > 85 && ymax > 87.5 ) ymax = 90;

  if ( cdoVerbose )
    cdoPrint("genYmax: y1 = %g  y2 = %g  dy = %g  ymax = %g", y1, y2, dy, ymax);

  return ymax;
}



/*****************************************************************************/

void grid_cell_center_to_bounds_Y2D(const char* yunitstr, long xsize, long ysize, const double* restrict grid_center_lat, double* restrict grid_corner_lat)
{
  long i, j, index;
  double minlat, maxlat;
  double firstlat, lastlat;

  firstlat = grid_center_lat[0];
  lastlat  = grid_center_lat[xsize*ysize-1];

  // if ( ysize == 1 ) cdoAbort("Cannot calculate Ybounds for 1 value!");

  for ( j = 0; j < ysize; ++j )
    {
      if ( ysize == 1 )
	{
	  minlat = grid_center_lat[0] - 360./ysize;
	  maxlat = grid_center_lat[0] + 360./ysize;
	}
      else
	{
	  index = j*xsize;
	  if ( firstlat > lastlat )
	    {
	      if ( j == 0 )
		maxlat = genYmax(grid_center_lat[index], grid_center_lat[index+xsize]);
	      else
		maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);

	      if ( j == (ysize-1) )
		minlat = genYmin(grid_center_lat[index], grid_center_lat[index-xsize]);
	      else
		minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
	    }
	  else
	    {
	      if ( j == 0 )
		minlat = genYmin(grid_center_lat[index], grid_center_lat[index+xsize]);
	      else
		minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);

	      if ( j == (ysize-1) )
		maxlat = genYmax(grid_center_lat[index], grid_center_lat[index-xsize]);
	      else
		maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
	    }
	}

      for ( i = 0; i < xsize; ++i )
	{
	  index = (j<<2)*xsize + (i<<2);
	  grid_corner_lat[index  ] = minlat;
	  grid_corner_lat[index+1] = minlat;
	  grid_corner_lat[index+2] = maxlat;
	  grid_corner_lat[index+3] = maxlat;
	}
    }
}


void gridGenRotBounds(int gridID, int nx, int ny,
		      double *xbounds, double *ybounds, double *xbounds2D, double *ybounds2D)
{
  long i, j, index;
  double minlon, maxlon;
  double minlat, maxlat;
  double xpole, ypole, angle;

  xpole = gridInqXpole(gridID);
  ypole = gridInqYpole(gridID);
  angle = gridInqAngle(gridID);

  for ( j = 0; j < ny; j++ )
    {
      if ( ybounds[0] > ybounds[1] )
	{
	  maxlat = ybounds[2*j];
	  minlat = ybounds[2*j+1];
	}
      else
	{
	  maxlat = ybounds[2*j+1];
	  minlat = ybounds[2*j];
	}

      for ( i = 0; i < nx; i++ )
	{
	  minlon = xbounds[2*i];
	  maxlon = xbounds[2*i+1];

	  index = j*4*nx + 4*i;
	  xbounds2D[index+0] = lamrot_to_lam(minlat, minlon, ypole, xpole, angle);
	  xbounds2D[index+1] = lamrot_to_lam(minlat, maxlon, ypole, xpole, angle);
	  xbounds2D[index+2] = lamrot_to_lam(maxlat, maxlon, ypole, xpole, angle);
	  xbounds2D[index+3] = lamrot_to_lam(maxlat, minlon, ypole, xpole, angle);

	  ybounds2D[index+0] = phirot_to_phi(minlat, minlon, ypole, angle);
	  ybounds2D[index+1] = phirot_to_phi(minlat, maxlon, ypole, angle);
	  ybounds2D[index+2] = phirot_to_phi(maxlat, maxlon, ypole, angle);
	  ybounds2D[index+3] = phirot_to_phi(maxlat, minlon, ypole, angle);
	}
    }
}


void grid_gen_xbounds2D(int nx, int ny, const double *restrict xbounds, double *restrict xbounds2D)
{
  int index;
  double minlon, maxlon;

#if defined(_OPENMP)
#pragma omp parallel for default(none)        \
                          shared(nx, ny, xbounds, xbounds2D) \
                         private(minlon, maxlon, index)
#endif
  for ( int i = 0; i < nx; ++i )
    {
      minlon = xbounds[2*i  ];
      maxlon = xbounds[2*i+1];

      for ( int j = 0; j < ny; ++j )
	{
	  index = j*4*nx + 4*i;
	  xbounds2D[index  ] = minlon;
	  xbounds2D[index+1] = maxlon;
	  xbounds2D[index+2] = maxlon;
	  xbounds2D[index+3] = minlon;
	}
    }
}


void grid_gen_ybounds2D(int nx, int ny, const double *restrict ybounds, double *restrict ybounds2D)
{
  int index;
  double minlat, maxlat;

#if defined(_OPENMP)
#pragma omp parallel for default(none)        \
                          shared(nx, ny, ybounds, ybounds2D) \
                         private(minlat, maxlat, index)
#endif
  for ( int j = 0; j < ny; ++j )
    {
      if ( ybounds[0] > ybounds[1] )
	{
	  maxlat = ybounds[2*j  ];
	  minlat = ybounds[2*j+1];
	}
      else
	{
	  maxlat = ybounds[2*j+1];
	  minlat = ybounds[2*j  ];
	}

      for ( int i = 0; i < nx; ++i )
	{
	  index = j*4*nx + 4*i;
	  ybounds2D[index  ] = minlat;
	  ybounds2D[index+1] = minlat;
	  ybounds2D[index+2] = maxlat;
	  ybounds2D[index+3] = maxlat;
	}
    }
}

static
char *gen_param(const char *fmt, ...)
{
  va_list args;
  char str[256];
  char *rstr;
  int len;

  va_start(args, fmt);

  len = vsprintf(str, fmt, args);

  va_end(args);

  len++;
  rstr = (char*) Malloc(len*sizeof(char));
  memcpy(rstr, str, len*sizeof(char));

  return rstr;
}

static
void lcc_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
  double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
  double zlat, zlon;
  double xi, xj;
  int projflag, scanflag;
  long i;
  proj_info_t proj;

  gridInqLCC(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm, &projflag, &scanflag);
  /*
    while ( originLon < 0 ) originLon += 360;
    while ( lonParY   < 0 ) lonParY   += 360;
  */
  if ( IS_NOT_EQUAL(xincm, yincm) )
    Warning("X and Y increment must be equal on Lambert Conformal grid (Xinc = %g, Yinc = %g)\n", 
	    xincm, yincm);
  /*
  if ( IS_NOT_EQUAL(lat1, lat2) )
    Warning("Lat1 and Lat2 must be equal on Lambert Conformal grid (Lat1 = %g, Lat2 = %g)\n", 
	    lat1, lat2);
  */
  map_set(PROJ_LC, originLat, originLon, xincm, lonParY, lat1, lat2, &proj);

  for ( i = 0; i < gridsize; i++ )
    {
      xi = xvals[i];
      xj = yvals[i];
      // status = W3FB12(xi, xj, originLat, originLon, xincm, lonParY, lat1, &zlat, &zlon);
      ijll_lc(xi, xj, proj, &zlat, &zlon);
      xvals[i] = zlon;
      yvals[i] = zlat;
    }
}

static
void sinusoidal_to_geo(int gridsize, double *xvals, double *yvals)
{
#if defined(HAVE_LIBPROJ)
  char *params[20];
  projUV data, res;
  long i;

  int nbpar = 0;
  params[nbpar++] = gen_param("proj=sinu");
  params[nbpar++] = gen_param("ellps=WGS84");

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%ld] = %s", i+1, params[i]);

  projPJ proj = pj_init(nbpar, params);
  if ( !proj )
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

  for ( i = 0; i < nbpar; ++i ) Free(params[i]);

  /* proj->over = 1; */		/* allow longitude > 180 */

  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
      res = pj_inv(data, proj);
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
      if ( xvals[i] < -9000. || xvals[i] > 9000. ) xvals[i] = -9999.;
      if ( yvals[i] < -9000. || yvals[i] > 9000. ) yvals[i] = -9999.;
    }

  pj_free(proj);
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

static
void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
#if defined(HAVE_LIBPROJ)
  char *params[20];
  projUV data, res;
  double a, lon_0, lat_0;
  long i;

  gridInqLaea(gridID, &a , &lon_0, &lat_0);

  int nbpar = 0;
  params[nbpar++] = gen_param("proj=laea");
  if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
  params[nbpar++] = gen_param("lon_0=%g", lon_0);
  params[nbpar++] = gen_param("lat_0=%g", lat_0);

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%d] = %s", i+1, params[i]);

  projPJ proj = pj_init(nbpar, &params[0]);
  if ( !proj )
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

  for ( i = 0; i < nbpar; ++i ) Free(params[i]);

  /* proj->over = 1; */		/* allow longitude > 180 */

  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
      res = pj_inv(data, proj);
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
      if ( xvals[i] < -9000. || xvals[i] > 9000. ) xvals[i] = -9999.;
      if ( yvals[i] < -9000. || yvals[i] > 9000. ) yvals[i] = -9999.;
    }

  pj_free(proj);
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

static
void lcc2_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
#if defined(HAVE_LIBPROJ)
  char *params[20];
  projUV data, res;
  double a, lon_0, lat_0, lat_1, lat_2;
  long i;

  gridInqLcc2(gridID, &a , &lon_0, &lat_0, &lat_1, &lat_2);

  int nbpar = 0;
  params[nbpar++] = gen_param("proj=lcc");
  if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
  params[nbpar++] = gen_param("lon_0=%g", lon_0);
  params[nbpar++] = gen_param("lat_0=%g", lat_0);
  params[nbpar++] = gen_param("lat_1=%g", lat_1);
  params[nbpar++] = gen_param("lat_2=%g", lat_2);

  if ( cdoVerbose )
    for ( i = 0; i < nbpar; ++i )
      cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
  
  projPJ proj = pj_init(nbpar, &params[0]);
  if ( !proj )
    cdoAbort("proj error: %s", pj_strerrno(pj_errno));

  for ( i = 0; i < nbpar; ++i ) Free(params[i]);

  /* proj->over = 1; */		/* allow longitude > 180 */
  
  for ( i = 0; i < gridsize; i++ )
    {
      data.u = xvals[i];
      data.v = yvals[i];
      res = pj_inv(data, proj);
      xvals[i] = res.u*RAD2DEG;
      yvals[i] = res.v*RAD2DEG;
    }

  pj_free(proj);
#else
  cdoAbort("proj4 support not compiled in!");
#endif
}

/*
 * grib_get_reduced_row: code from GRIB_API 1.10.4
 *
 * Description:
 *   computes the number of points within the range lon_first->lon_last and the zero based indexes
 *   ilon_first,ilon_last of the first and last point given the number of points along a parallel (pl)
 *
 */
static
void grib_get_reduced_row(long pl,double lon_first,double lon_last,long* npoints,long* ilon_first, long* ilon_last )
{
  double range=0,dlon_first=0,dlon_last=0;
  long irange;

  range=lon_last-lon_first;
  if (range<0) {range+=360;lon_first-=360;}

  /* computing integer number of points and coordinates without using floating point resolution*/
  *npoints=(range*pl)/360.0+1;
  *ilon_first=(lon_first*pl)/360.0;
  *ilon_last=(lon_last*pl)/360.0;

  irange=*ilon_last-*ilon_first+1;

  if (irange != *npoints) {
    if (irange > *npoints) {
      /* checking if the first point is out of range*/
      dlon_first=((*ilon_first)*360.0)/pl;
      if (dlon_first < lon_first) {(*ilon_first)++;irange--;
      }

      /* checking if the last point is out of range*/
      dlon_last=((*ilon_last)*360.0)/pl;
      if (dlon_last > lon_last) {(*ilon_last)--;irange--;
      }
    } else {
      int ok=0;
      /* checking if the point before the first is in the range*/
      dlon_first=((*ilon_first-1)*360.0)/pl;
      if (dlon_first > lon_first) {(*ilon_first)--;irange++;ok=1;
      }

      /* checking if the point after the last is in the range*/
      dlon_last=((*ilon_last+1)*360.0)/pl;
      if (dlon_last < lon_last) {(*ilon_last)++;irange++;ok=1;
      }

      /* if neither of the two are triggered then npoints is too large */
      if (!ok) {(*npoints)--;
      }
    }

    //   assert(*npoints==irange);
  } else {
	  /* checking if the first point is out of range*/
	  dlon_first=((*ilon_first)*360.0)/pl;
	  if (dlon_first < lon_first) {
		  (*ilon_first)++;(*ilon_last)++;
	  }
  }

  if (*ilon_first<0) *ilon_first+=pl;
}

static
int qu2reg_subarea(int gridsize, int np, double xfirst, double xlast, 
		   double *array, int *rowlon, int ny, double missval, int *iret, int lmiss, int lperio, int lveggy)
{
  int nx = 0;
  /* sub area (longitudes) */
  long ilon_firstx;
  long ilon_first, ilon_last;
  int i, j;
  long row_count;
  int rlon, nwork = 0;
  int np4 = np*4;
  int size = 0;
  int wlen;
  double *work, **pwork;
  int ii;
  double *parray;

  if ( np <= 0 ) cdoAbort("Number of values between pole and equator missing!");

  grib_get_reduced_row(np4, xfirst, xlast, &row_count, &ilon_firstx, &ilon_last);
  nx = row_count;
  // printf("nx %d  %ld %ld lon1 %g lon2 %g\n", nx, ilon_firstx, ilon_last, (ilon_firstx*360.)/np4, (ilon_last*360.)/np4);

  for ( j = 0; j < ny; ++j ) nwork += rowlon[j];

  pwork = (double **) Malloc(ny*sizeof(double *));
  work  = (double*) Malloc(ny*np4*sizeof(double));
  wlen = 0;
  pwork[0] = work;
  for ( j = 1; j < ny; ++j )
    {
      wlen += rowlon[j-1];
      pwork[j] = work + wlen;
    } 
  // printf(" ny, np4, nwork %d %d %d wlen %d\n", ny, np4, nwork, wlen);

  for ( j = 0; j < ny; ++j )
    {
      rlon = rowlon[j];
      for ( i = 0; i < rlon; ++i ) pwork[j][i] = missval;
    } 

  parray = array;
  for ( j = 0; j < ny; ++j )
    {
      rlon = rowlon[j];
      row_count = 0;
      grib_get_reduced_row(rlon, xfirst, xlast, &row_count, &ilon_first, &ilon_last);
      // printf("j %d xfirst %g xlast %g rowlon %d %ld %ld %ld %g %g\n", j, xfirst, xlast, rlon, row_count, ilon_first, ilon_last, (ilon_first*360.)/rlon, (ilon_last*360.)/rlon);

      for ( i = ilon_first; i < (ilon_first+row_count); ++i )
	{
	  ii = i;
	  if ( ii >= rlon ) ii -= rlon; 
	  pwork[j][ii] = *parray;
	  parray++;
	}
      size += row_count;
    }

  if ( gridsize != size ) cdoAbort("gridsize1 inconsistent!");

  (void) qu2reg3_double(work, rowlon, ny, np4, missval, iret, lmiss, lperio, lveggy);

  wlen = 0;
  pwork[0] = work;
  for ( j = 1; j < ny; ++j )
    {
      wlen += np4;
      pwork[j] = work + wlen;
    } 

  // printf("nx, ilon_firstx %d %ld\n", nx, ilon_firstx);
  parray = array;
  for ( j = 0; j < ny; ++j )
    {
      for ( i = ilon_firstx; i < (ilon_firstx+nx); ++i )
	{
	  ii = i;
	  if ( ii >= np4 ) ii -= np4; 
	  *parray = pwork[j][ii];
	  parray++;
	}
    }

  Free(work);
  Free(pwork);

  return nx;
}


void field2regular(int gridID1, int gridID2, double missval, double *array, int nmiss)
{
  int nx = 0, ny, np;
  int gridtype;
  int lmiss, lperio, lveggy;
  int iret;
  int *rowlon;
  double xfirstandlast[2];
  double xfirst, xlast;

  gridtype = gridInqType(gridID1);

  if ( gridtype != GRID_GAUSSIAN_REDUCED ) Error("Not a reduced Gaussian grid!");

  lmiss = nmiss > 0;
  lperio = 1;
  lveggy = 0;

  ny = gridInqYsize(gridID1);
  np = gridInqNP(gridID1);

  rowlon = (int*) Malloc(ny*sizeof(int));
  gridInqRowlon(gridID1, rowlon);

  xfirstandlast[0] = 0.;
  xfirstandlast[1] = 0.;
  gridInqXvals(gridID1, xfirstandlast);
  xfirst = xfirstandlast[0];
  xlast  = xfirstandlast[1];

  if ( fabs(xfirst) > 0 || (np > 0 && fabs(xlast - (360.0-90.0/np)) > 90.0/np) )
    {
      nx = qu2reg_subarea(gridInqSize(gridID1), np, xfirst, xlast, array, rowlon, ny, missval, &iret, lmiss, lperio, lveggy);
    }
  else
    {
      nx = 2*ny;
      (void) qu2reg3_double(array, rowlon, ny, nx, missval, &iret, lmiss, lperio, lveggy);
    }

  if ( gridInqSize(gridID2) != nx*ny ) Error("Gridsize differ!");

  Free(rowlon);
}


int gridToRegular(int gridID1)
{
  int gridID2;
  int gridtype, gridsize;
  int nx = 0, ny, np;
  long i;
  double *xvals = NULL, *yvals = NULL;
  double xfirstandlast[2];
  double xfirst, xlast;

  gridtype = gridInqType(gridID1);

  if ( gridtype != GRID_GAUSSIAN_REDUCED ) Error("Not a reduced Gaussian grid!");

  ny = gridInqYsize(gridID1);
  np = gridInqNP(gridID1);

  yvals = (double*) Malloc(ny*sizeof(double));
  gridInqYvals(gridID1, yvals);

  xfirstandlast[0] = 0.;
  xfirstandlast[1] = 0.;
  gridInqXvals(gridID1, xfirstandlast);
  xfirst = xfirstandlast[0];
  xlast  = xfirstandlast[1];

  if ( fabs(xfirst) > 0 || (np > 0 && fabs(xlast - (360.0-90.0/np)) > 90.0/np) )
    {
      /* sub area (longitudes) */
      long ilon_first, ilon_last;
      long row_count;
      int np4 = np*4;
      int *rowlon = NULL;

      if ( np <= 0 ) cdoAbort("Number of values between pole and equator missing!");

      grib_get_reduced_row(np4, xfirst, xlast, &row_count, &ilon_first, &ilon_last);

      nx = row_count;
      xvals = (double*) Malloc(nx*sizeof(double));
      for ( i = 0; i < nx; ++i )
	{
	  xvals[i] = ((ilon_first+i)*360.)/np4;
	  if ( xfirst > xlast ) xvals[i] -= 360.;
	}

      Free(rowlon);
    }
  else
    {
      nx = 2*ny;
      xvals = (double*) Malloc(nx*sizeof(double));
      for ( i = 0; i < nx; ++i ) xvals[i] = i * 360./nx;
    }

  gridsize = nx*ny;

  gridID2  = gridCreate(GRID_GAUSSIAN, gridsize);
	  
  gridDefXsize(gridID2, nx);
  gridDefYsize(gridID2, ny);
  
  gridDefXvals(gridID2, xvals);
  gridDefYvals(gridID2, yvals);
  gridDefNP(gridID2, np);

  Free(xvals);
  Free(yvals);

  return gridID2;
}

static
void gridCopyMask(int gridID1, int gridID2, long gridsize)
{
  if ( gridInqMask(gridID1, NULL) )
    {
      int *mask = (int*) Malloc(gridsize*sizeof(int));
      gridInqMask(gridID1, mask);
      gridDefMask(gridID2, mask);
      Free(mask);
    }
}

static
int check_range(long n, double *vals, double valid_min, double valid_max)
{
  int status = 0;
  long i;

  for ( i = 0; i < n; ++i )
    {
      if ( vals[i] < valid_min || vals[i] > valid_max )
	{
	  status = 1;
	  break;
	}
    }

  return status;
}


int gridToCurvilinear(int gridID1, int lbounds)
{
  long index;
  int gridtype = gridInqType(gridID1);
  size_t gridsize = (size_t) gridInqSize(gridID1);
  int gridID2  = gridCreate(GRID_CURVILINEAR, (int) gridsize);
  gridDefPrec(gridID2, DATATYPE_FLT32);
	  
  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
    case GRID_LCC:
    case GRID_LCC2:
    case GRID_LAEA:
    case GRID_SINUSOIDAL:
      {
	double xscale = 1, yscale = 1;
	double *xvals = NULL, *yvals = NULL;
	double *xbounds = NULL, *ybounds = NULL;
	char xunits[CDI_MAX_NAME], yunits[CDI_MAX_NAME];

	int nx = gridInqXsize(gridID1);
	int ny = gridInqYsize(gridID1);

	gridInqXunits(gridID1, xunits);
	gridInqYunits(gridID1, yunits);

	if ( gridtype == GRID_LAEA )
	  {
	    int lvalid_xunits = FALSE;
	    int lvalid_yunits = FALSE;
	    int len;
	    len = (int) strlen(xunits);
	    if ( len == 1 && memcmp(xunits, "m",  1) == 0 ) lvalid_xunits = TRUE;
	    if ( len == 2 && memcmp(xunits, "km", 2) == 0 ) lvalid_xunits = TRUE;
	    len = (int) strlen(yunits);
	    if ( len == 1 && memcmp(yunits, "m",  1) == 0 ) lvalid_yunits = TRUE;
	    if ( len == 2 && memcmp(yunits, "km", 2) == 0 ) lvalid_yunits = TRUE;

	    if ( lvalid_xunits == FALSE )
	      cdoWarning("Possibly wrong result! Invalid x-coordinate units: \"%s\" (expected \"m\" or \"km\")", xunits);
	    if ( lvalid_yunits == FALSE )
	      cdoWarning("Possibly wrong result! Invalid y-coordinate units: \"%s\" (expected \"m\" or \"km\")", yunits);
	  }

	if ( memcmp(xunits, "km", 2) == 0 ) xscale = 1000;
	if ( memcmp(yunits, "km", 2) == 0 ) yscale = 1000;

	gridDefXsize(gridID2, nx);
	gridDefYsize(gridID2, ny);

	double *xvals2D = (double*) Malloc(gridsize*sizeof(double));
	double *yvals2D = (double*) Malloc(gridsize*sizeof(double));


	if ( gridtype == GRID_LCC )
	  {
	    for ( int j = 0; j < ny; j++ )
	      for ( int i = 0; i < nx; i++ )
		{
		  xvals2D[j*nx+i] = i+1;
		  yvals2D[j*nx+i] = j+1;
		}

	    lcc_to_geo(gridID1, gridsize, xvals2D, yvals2D);