remaplib.c

/*
  This is a C library of the Fortran SCRIP version 1.4

  ===>>> Please send bug reports to <http://mpimet.mpg.de/cdo> <<<===

  Spherical Coordinate Remapping and Interpolation Package (SCRIP)
  ================================================================

  SCRIP is a software package which computes addresses and weights for
  remapping and interpolating fields between grids in spherical coordinates.
  It was written originally for remapping fields to other grids in a coupled
  climate model, but is sufficiently general that it can be used in other 
  applications as well. The package should work for any grid on the surface
  of a sphere. SCRIP currently supports four remapping options:

  Conservative remapping
  ----------------------
  First- and second-order conservative remapping as described in
  Jones (1999, Monthly Weather Review, 127, 2204-2210).

  Bilinear interpolation
  ----------------------
  Slightly generalized to use a local bilinear approximation
  (only logically-rectangular grids).

  Bicubic interpolation
  ----------------------
  Similarly generalized (only logically-rectangular grids).

  Distance-weighted averaging
  ---------------------------
  Distance-weighted average of a user-specified number of nearest neighbor values.

  Documentation
  =============

  http://climate.lanl.gov/Software/SCRIP/SCRIPusers.pdf

*/
/*
  2013-11-08 Uwe Schulzweida: split remapgrid class to src_grid and tgt_grid
  2012-01-16 Uwe Schulzweida: alloc grid2_bound_box only for conservative remapping
  2011-01-07 Uwe Schulzweida: Changed remap weights from 2D to 1D array
  2009-05-25 Uwe Schulzweida: Changed restrict data type from double to int
  2009-01-11 Uwe Schulzweida: OpenMP parallelization
 */

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

#include <limits.h>
#include <time.h>

#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
#include "remap_store_link_cnsrv.h"


#define IS_REG2D_GRID(gridID)  (gridInqType(gridID) == GRID_LONLAT || gridInqType(gridID) == GRID_GAUSSIAN)



static bool remap_gen_weights     = true;
static bool remap_write_remap     = false;
static int  remap_num_srch_bins   = 180;
#define  DEFAULT_MAX_ITER  100
long remap_max_iter        = DEFAULT_MAX_ITER;  /* Max iteration count for i, j iteration */

void remap_set_int(int remapvar, int value)
{
  if      ( remapvar == REMAP_STORE_LINK_FAST ) remap_store_link_fast = value;
  else if ( remapvar == REMAP_WRITE_REMAP     ) remap_write_remap     = value > 0;
  else if ( remapvar == REMAP_MAX_ITER        ) remap_max_iter        = value;
  else if ( remapvar == REMAP_NUM_SRCH_BINS   ) remap_num_srch_bins   = value;
  else if ( remapvar == REMAP_GENWEIGHTS      ) remap_gen_weights     = value > 0;
  else      cdoAbort("Unsupported remap variable (%d)!", remapvar);
}

double intlin(double x, double y1, double x1, double y2, double x2);

void remapGridFree(remapgrid_t *grid)
{
  if ( grid->vgpm ) Free(grid->vgpm);
  if ( grid->mask ) Free(grid->mask);

  if ( grid->reg2d_center_lat ) Free(grid->reg2d_center_lat);
  if ( grid->reg2d_center_lon ) Free(grid->reg2d_center_lon);
  if ( grid->reg2d_corner_lat ) Free(grid->reg2d_corner_lat);
  if ( grid->reg2d_corner_lon ) Free(grid->reg2d_corner_lon);

  if ( grid->cell_center_lat ) Free(grid->cell_center_lat);
  if ( grid->cell_center_lon ) Free(grid->cell_center_lon);
  if ( grid->cell_corner_lat ) Free(grid->cell_corner_lat);
  if ( grid->cell_corner_lon ) Free(grid->cell_corner_lon);

  if ( grid->cell_area ) Free(grid->cell_area);
  if ( grid->cell_frac ) Free(grid->cell_frac);

  if ( grid->cell_bound_box ) Free(grid->cell_bound_box);

  if ( grid->bin_addr ) Free(grid->bin_addr);
  if ( grid->bin_lats ) Free(grid->bin_lats);

} /* remapGridFree */

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

void remapVarsFree(remapvars_t *rv)
{
  if ( rv->pinit )
    {
      rv->pinit    = false;
      rv->sort_add = false;

      if ( rv->src_cell_add ) Free(rv->src_cell_add);
      if ( rv->tgt_cell_add ) Free(rv->tgt_cell_add);
      if ( rv->wts ) Free(rv->wts);

      if ( rv->links.option )
	{
	  rv->links.option = false;

	  if ( rv->links.num_blks )
	    {
	      Free(rv->links.num_links);
	      long num_blks = rv->links.num_blks;
	      for ( long i = 0; i < num_blks; ++i )
		{
		  Free(rv->links.src_add[i]);
		  Free(rv->links.dst_add[i]);
		  Free(rv->links.w_index[i]);
		}
	      Free(rv->links.src_add);
	      Free(rv->links.dst_add);
	      Free(rv->links.w_index);
	    }
	}
    }
  else
    fprintf(stderr, "%s Warning: vars not initialized!\n", __func__);

} /* remapVarsFree */

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

void remapgrid_init(remapgrid_t *grid)
{
  grid->remap_grid_type  = -1;
  grid->num_srch_bins    = remap_num_srch_bins; // only for source grid ?

  grid->num_cell_corners = 0;
  grid->luse_cell_corners  = false;
  grid->lneed_cell_corners = false;

  grid->nvgp             = 0;
  grid->vgpm             = NULL;

  grid->mask             = NULL;

  grid->reg2d_center_lon = NULL;
  grid->reg2d_center_lat = NULL;
  grid->reg2d_corner_lon = NULL;
  grid->reg2d_corner_lat = NULL;

  grid->cell_center_lon  = NULL;
  grid->cell_center_lat  = NULL;
  grid->cell_corner_lon  = NULL;
  grid->cell_corner_lat  = NULL;

  grid->cell_area        = NULL;
  grid->cell_frac        = NULL;

  grid->cell_bound_box   = NULL;

  grid->bin_addr         = NULL;
  grid->bin_lats         = NULL;
}

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

void remapgrid_alloc(int map_type, remapgrid_t *grid)
{
  if ( grid->nvgp )
    grid->vgpm   = (int*) Malloc(grid->nvgp*sizeof(int));

  grid->mask     = (int*) Malloc(grid->size*sizeof(int));

  if ( remap_write_remap || grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
    {
      grid->cell_center_lon = (double*) Malloc(grid->size*sizeof(double));
      grid->cell_center_lat = (double*) Malloc(grid->size*sizeof(double));
    }

  if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC )
    {
      grid->cell_area = (double*) Malloc(grid->size*sizeof(double));
      memset(grid->cell_area, 0, grid->size*sizeof(double));
    }

  grid->cell_frac = (double*) Malloc(grid->size*sizeof(double));
  memset(grid->cell_frac, 0, grid->size*sizeof(double));

  if ( grid->lneed_cell_corners )
    {
      if ( grid->num_cell_corners > 0 )
        {
	  long nalloc = grid->num_cell_corners*grid->size;

	  grid->cell_corner_lon = (double*) Malloc(nalloc*sizeof(double));
	  memset(grid->cell_corner_lon, 0, nalloc*sizeof(double));

	  grid->cell_corner_lat = (double*) Malloc(nalloc*sizeof(double));  
	  memset(grid->cell_corner_lat, 0, nalloc*sizeof(double));
	}
    }
}

/*****************************************************************************/
static
void boundbox_from_corners(long size, long nc, const double *restrict corner_lon,
			   const double *restrict corner_lat, restr_t *restrict bound_box)
{
  long i4, inc, j;
  restr_t clon, clat;

#if defined(_OPENMP)
#pragma omp parallel for default(none)        \
  shared(bound_box, corner_lat, corner_lon, nc, size)	\
  private(i4, inc, j, clon, clat)
#endif
  for ( long i = 0; i < size; ++i )
    {
      i4 = i<<2; // *4
      inc = i*nc;
      clat = RESTR_SCALE(corner_lat[inc]);
      clon = RESTR_SCALE(corner_lon[inc]);
      bound_box[i4  ] = clat;
      bound_box[i4+1] = clat;
      bound_box[i4+2] = clon;
      bound_box[i4+3] = clon;
      for ( j = 1; j < nc; ++j )
	{
	  clat = RESTR_SCALE(corner_lat[inc+j]);
	  clon = RESTR_SCALE(corner_lon[inc+j]);
	  if ( clat < bound_box[i4  ] ) bound_box[i4  ] = clat;
	  if ( clat > bound_box[i4+1] ) bound_box[i4+1] = clat;
	  if ( clon < bound_box[i4+2] ) bound_box[i4+2] = clon;
	  if ( clon > bound_box[i4+3] ) bound_box[i4+3] = clon;
	}
    }
}

static
void boundbox_from_center(bool lonIsCyclic, long size, long nx, long ny, const double *restrict center_lon,
			  const double *restrict center_lat, restr_t *restrict bound_box)
{
  long n4, i, j, k, ip1, jp1;
  long n_add, e_add, ne_add;
  restr_t tmp_lats[4], tmp_lons[4];  /* temps for computing bounding boxes */

#if defined(_OPENMP)
#pragma omp parallel for default(none)        \
  shared(lonIsCyclic, size, nx, ny, center_lon, center_lat, bound_box)	\
  private(n4, i, j, k, ip1, jp1, n_add, e_add, ne_add, tmp_lats, tmp_lons)
#endif
  for ( long n = 0; n < size; n++ )
    {
      n4 = n<<2;

      /* Find N,S and NE points to this grid point */
      
      j = n/nx;
      i = n - j*nx;

      if ( i < (nx-1) )
	ip1 = i + 1;
      else
	{
	  /* 2009-01-09 Uwe Schulzweida: bug fix */
	  ip1 = lonIsCyclic ? 0 : i;
	}

      if ( j < (ny-1) )
	jp1 = j + 1;
      else
	{
	  /* 2008-12-17 Uwe Schulzweida: latitute cyclic ??? (bug fix) */
	  jp1 = j;
	}

      n_add  = jp1*nx + i;
      e_add  = j  *nx + ip1;
      ne_add = jp1*nx + ip1;

      /* Find N,S and NE lat/lon coords and check bounding box */

      tmp_lats[0] = RESTR_SCALE(center_lat[n]);
      tmp_lats[1] = RESTR_SCALE(center_lat[e_add]);
      tmp_lats[2] = RESTR_SCALE(center_lat[ne_add]);
      tmp_lats[3] = RESTR_SCALE(center_lat[n_add]);

      tmp_lons[0] = RESTR_SCALE(center_lon[n]);
      tmp_lons[1] = RESTR_SCALE(center_lon[e_add]);
      tmp_lons[2] = RESTR_SCALE(center_lon[ne_add]);
      tmp_lons[3] = RESTR_SCALE(center_lon[n_add]);

      bound_box[n4  ] = tmp_lats[0];
      bound_box[n4+1] = tmp_lats[0];
      bound_box[n4+2] = tmp_lons[0];
      bound_box[n4+3] = tmp_lons[0];

      for ( k = 1; k < 4; k++ )
	{
	  if ( tmp_lats[k] < bound_box[n4  ] ) bound_box[n4  ] = tmp_lats[k];
	  if ( tmp_lats[k] > bound_box[n4+1] ) bound_box[n4+1] = tmp_lats[k];
	  if ( tmp_lons[k] < bound_box[n4+2] ) bound_box[n4+2] = tmp_lons[k];
	  if ( tmp_lons[k] > bound_box[n4+3] ) bound_box[n4+3] = tmp_lons[k];
	}
    }
}


void remapgrid_get_lonlat(remapgrid_t *grid, unsigned cell_add, double *plon, double *plat)
{
  if ( grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
    {
      unsigned nx = grid->dims[0];
      unsigned iy = cell_add/nx;
      unsigned ix = cell_add - iy*nx;
      *plat = grid->reg2d_center_lat[iy];
      *plon = grid->reg2d_center_lon[ix];
      if ( *plon < 0 ) *plon += PI2;
    }
  else
    {
      *plat = grid->cell_center_lat[cell_add];
      *plon = grid->cell_center_lon[cell_add];
    }
}


void check_lon_range(long nlons, double *lons)
{
  assert(lons != NULL);

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(nlons, lons)
#endif
  for ( long n = 0; n < nlons; ++n )
    {
      // remove missing values
      if ( lons[n] < -PI2 ) lons[n] = 0;
      if ( lons[n] > 2*PI2) lons[n] = PI2;

      if ( lons[n] > PI2  ) lons[n] -= PI2;
      if ( lons[n] < ZERO ) lons[n] += PI2;
    }
}


void check_lat_range(long nlats, double *lats)
{
  assert(lats != NULL);

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(nlats, lats)
#endif
  for ( long n = 0; n < nlats; ++n )
    {
      if ( lats[n] >  PIH ) lats[n] =  PIH;
      if ( lats[n] < -PIH ) lats[n] = -PIH;
    }
}

static
void check_lon_boundbox_range(long nlons, restr_t *bound_box)
{
  long n4;

  assert(bound_box != NULL);

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(nlons, bound_box) private(n4)
#endif
  for ( long n = 0; n < nlons; ++n )
    {
      n4 = n<<2;
      if ( RESTR_ABS(bound_box[n4+3] - bound_box[n4+2]) > RESTR_SCALE(PI) )
	{
	  bound_box[n4+2] = RESTR_SCALE(0.);
	  bound_box[n4+3] = RESTR_SCALE(PI2);
	}
    }
}

static
void check_lat_boundbox_range(long nlats, restr_t *restrict bound_box, double *restrict lats)
{
  long n4;

  assert(bound_box != NULL);

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(nlats, bound_box, lats) private(n4)
#endif
  for ( long n = 0; n < nlats; ++n )
    {
      n4 = n<<2;
      if ( RESTR_SCALE(lats[n]) < bound_box[n4  ] ) bound_box[n4  ] = RESTR_SCALE(-PIH);
      if ( RESTR_SCALE(lats[n]) > bound_box[n4+1] ) bound_box[n4+1] = RESTR_SCALE( PIH);
    }
}

static
int expand_lonlat_grid(int gridID)
{
  long nx = gridInqXsize(gridID);
  long ny = gridInqYsize(gridID);
  long nxp4 = nx+4;
  long nyp4 = ny+4;

  double *xvals = (double*) Malloc(nxp4*sizeof(double));
  double *yvals = (double*) Malloc(nyp4*sizeof(double));
  gridInqXvals(gridID, xvals+2);
  gridInqYvals(gridID, yvals+2);

  int gridtype = gridInqType(gridID);
  int gridIDnew = gridCreate(gridtype, nxp4*nyp4);
  gridDefXsize(gridIDnew, nxp4);
  gridDefYsize(gridIDnew, nyp4);

  grid_copy_attributes(gridID, gridIDnew);

  xvals[0] = xvals[2] - 2*gridInqXinc(gridID);
  xvals[1] = xvals[2] - gridInqXinc(gridID);
  xvals[nxp4-2] = xvals[nx+1] + gridInqXinc(gridID);
  xvals[nxp4-1] = xvals[nx+1] + 2*gridInqXinc(gridID);

  yvals[0] = yvals[2] - 2*gridInqYinc(gridID);
  yvals[1] = yvals[2] - gridInqYinc(gridID);
  yvals[nyp4-2] = yvals[ny+1] + gridInqYinc(gridID);
  yvals[nyp4-1] = yvals[ny+1] + 2*gridInqYinc(gridID);

  gridDefXvals(gridIDnew, xvals);
  gridDefYvals(gridIDnew, yvals);

  Free(xvals);
  Free(yvals);

  if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL )
    {
      double xpole, ypole, angle;
      gridInqParamRLL(gridID, &xpole, &ypole, &angle);
      gridDefParamRLL(gridIDnew, xpole, ypole, angle);
    }

  return gridIDnew;
}

static
int expand_curvilinear_grid(int gridID)
{
  long i, j;

  long gridsize = gridInqSize(gridID);
  long nx = gridInqXsize(gridID);
  long ny = gridInqYsize(gridID);
  long nxp4 = nx+4;
  long nyp4 = ny+4;
  int gridsize_new = gridsize + 4*(nx+2) + 4*(ny+2);

  double *xvals = (double*) Malloc(gridsize_new*sizeof(double));
  double *yvals = (double*) Malloc(gridsize_new*sizeof(double));
  gridInqXvals(gridID, xvals);
  gridInqYvals(gridID, yvals);

  int gridIDnew = gridCreate(GRID_CURVILINEAR, nxp4*nyp4);
  gridDefXsize(gridIDnew, nxp4);
  gridDefYsize(gridIDnew, nyp4);

  grid_copy_attributes(gridID, gridIDnew);

  for ( j = ny-1; j >= 0; j-- )
    for ( i = nx-1; i >= 0; i-- )
      xvals[(j+2)*(nx+4)+i+2] = xvals[j*nx+i];

  for ( j = ny-1; j >= 0; j-- )
    for ( i = nx-1; i >= 0; i-- )
      yvals[(j+2)*(nx+4)+i+2] = yvals[j*nx+i];

  for ( j = 2; j < nyp4-2; j++ )
    {
      xvals[j*nxp4  ] = intlin(3.0, xvals[j*nxp4+3], 0.0, xvals[j*nxp4+2], 1.0);
      xvals[j*nxp4+1] = intlin(2.0, xvals[j*nxp4+3], 0.0, xvals[j*nxp4+2], 1.0); 
      yvals[j*nxp4  ] = intlin(3.0, yvals[j*nxp4+3], 0.0, yvals[j*nxp4+2], 1.0); 
      yvals[j*nxp4+1] = intlin(2.0, yvals[j*nxp4+3], 0.0, yvals[j*nxp4+2], 1.0); 

      xvals[j*nxp4+nxp4-2] = intlin(2.0, xvals[j*nxp4+nxp4-4], 0.0, xvals[j*nxp4+nxp4-3], 1.0); 
      xvals[j*nxp4+nxp4-1] = intlin(3.0, xvals[j*nxp4+nxp4-4], 0.0, xvals[j*nxp4+nxp4-3], 1.0); 
      yvals[j*nxp4+nxp4-2] = intlin(2.0, yvals[j*nxp4+nxp4-4], 0.0, yvals[j*nxp4+nxp4-3], 1.0); 
      yvals[j*nxp4+nxp4-1] = intlin(3.0, yvals[j*nxp4+nxp4-4], 0.0, yvals[j*nxp4+nxp4-3], 1.0); 
    }

  for ( i = 0; i < nxp4; i++ )
    {
      xvals[0*nxp4+i] = intlin(3.0, xvals[3*nxp4+i], 0.0, xvals[2*nxp4+i], 1.0);
      xvals[1*nxp4+i] = intlin(2.0, xvals[3*nxp4+i], 0.0, xvals[2*nxp4+i], 1.0);
      yvals[0*nxp4+i] = intlin(3.0, yvals[3*nxp4+i], 0.0, yvals[2*nxp4+i], 1.0);
      yvals[1*nxp4+i] = intlin(2.0, yvals[3*nxp4+i], 0.0, yvals[2*nxp4+i], 1.0);

      xvals[(nyp4-2)*nxp4+i] = intlin(2.0, xvals[(nyp4-4)*nxp4+i], 0.0, xvals[(nyp4-3)*nxp4+i], 1.0);
      xvals[(nyp4-1)*nxp4+i] = intlin(3.0, xvals[(nyp4-4)*nxp4+i], 0.0, xvals[(nyp4-3)*nxp4+i], 1.0);
      yvals[(nyp4-2)*nxp4+i] = intlin(2.0, yvals[(nyp4-4)*nxp4+i], 0.0, yvals[(nyp4-3)*nxp4+i], 1.0);
      yvals[(nyp4-1)*nxp4+i] = intlin(3.0, yvals[(nyp4-4)*nxp4+i], 0.0, yvals[(nyp4-3)*nxp4+i], 1.0);
    }
  /*
    {
    FILE *fp;
    fp = fopen("xvals.asc", "w");
    for ( i = 0; i < gridsize_new; i++ ) fprintf(fp, "%g\n", xvals[i]);
    fclose(fp);
    fp = fopen("yvals.asc", "w");
    for ( i = 0; i < gridsize_new; i++ ) fprintf(fp, "%g\n", yvals[i]);
    fclose(fp);
    }
  */
  gridDefXvals(gridIDnew, xvals);
  gridDefYvals(gridIDnew, yvals);
  
  Free(xvals);
  Free(yvals);

  return gridIDnew;
}

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

static
void grid_check_lat_borders_rad(int n, double *ybounds)
{
  if ( ybounds[0] > ybounds[n-1] )
    {
      if ( RAD2DEG*ybounds[0]   >  PIH ) ybounds[0]   =  PIH;
      if ( RAD2DEG*ybounds[n-1] < -PIH ) ybounds[n-1] = -PIH;
    }
  else
    {
      if ( RAD2DEG*ybounds[0]   < -PIH ) ybounds[0]   = -PIH;
      if ( RAD2DEG*ybounds[n-1] >  PIH ) ybounds[n-1] =  PIH;
    }
}

static
void remap_define_reg2d(int gridID, remapgrid_t *grid)
{
  long nx = grid->dims[0];
  long ny = grid->dims[1];

  long nxp1 = nx + 1;
  long nyp1 = ny + 1;

  long nxm = nx;
  if ( grid->is_cyclic ) nxm++;

  if ( grid->size != nx*ny ) cdoAbort("Internal error, wrong dimensions!");

  grid->reg2d_center_lon = (double*) Malloc(nxm*sizeof(double));
  grid->reg2d_center_lat = (double*) Malloc( ny*sizeof(double));
 
  grid->reg2d_center_lon[0] = 0;
  grid->reg2d_center_lat[0] = 0;
  gridInqXvals(gridID, grid->reg2d_center_lon);
  gridInqYvals(gridID, grid->reg2d_center_lat);

  /* Convert lat/lon units if required */

  char yunits[CDI_MAX_NAME]; yunits[0] = 0;
  cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);

  grid_to_radian(yunits, nx, grid->reg2d_center_lon, "grid reg2d center lon"); 
  grid_to_radian(yunits, ny, grid->reg2d_center_lat, "grid reg2d center lat"); 

  if ( grid->is_cyclic ) grid->reg2d_center_lon[nx] = grid->reg2d_center_lon[0] + PI2;

  grid->reg2d_corner_lon = (double*) Malloc(nxp1*sizeof(double));
  grid->reg2d_corner_lat = (double*) Malloc(nyp1*sizeof(double));

  grid_gen_corners(nx, grid->reg2d_center_lon, grid->reg2d_corner_lon);
  grid_gen_corners(ny, grid->reg2d_center_lat, grid->reg2d_corner_lat);
  grid_check_lat_borders_rad(ny+1, grid->reg2d_corner_lat);

  //for ( long i = 0; i < nxp1; ++i ) printf("lon %ld %g\n", i, grid->reg2d_corner_lon[i]);
  //for ( long i = 0; i < nyp1; ++i ) printf("lat %ld %g\n", i, grid->reg2d_corner_lat[i]);

}

static
void remap_define_grid(int map_type, int gridID, remapgrid_t *grid, const char *txt)
{
  bool lgrid_destroy = false;
  bool lgrid_gen_bounds = false;
  int gridID_gme = -1;

  if ( gridInqType(grid->gridID) != GRID_UNSTRUCTURED && gridInqType(grid->gridID) != GRID_CURVILINEAR )
    {
      if ( gridInqType(grid->gridID) == GRID_GME )
	{
	  gridID_gme = gridToUnstructured(grid->gridID, 1);
	  grid->nvgp = gridInqSize(gridID_gme);
	  gridID = gridDuplicate(gridID_gme);
	  gridCompress(gridID);
	  grid->luse_cell_corners = true;
	}
      else if ( remap_write_remap || grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
	{
	  lgrid_destroy = true;
	  gridID = gridToCurvilinear(grid->gridID, 1);
	  lgrid_gen_bounds = true;
	}
    }

  long gridsize = grid->size = gridInqSize(gridID);

  grid->dims[0] = gridInqXsize(gridID);
  grid->dims[1] = gridInqYsize(gridID);
  if ( gridInqType(grid->gridID) != GRID_UNSTRUCTURED )
    {
      if ( grid->dims[0] == 0 ) cdoAbort("%s grid without longitude coordinates!", gridNamePtr(gridInqType(grid->gridID)));
      if ( grid->dims[1] == 0 ) cdoAbort("%s grid without latitude coordinates!", gridNamePtr(gridInqType(grid->gridID)));
    }

  grid->is_cyclic = (gridIsCircular(gridID) > 0);

  grid->rank = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? 1 : 2;

  grid->num_cell_corners = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? gridInqNvertex(gridID) : 4;

  remapgrid_alloc(map_type, grid);

  /* Initialize logical mask */

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(gridsize, grid)
#endif
  for ( long i = 0; i < gridsize; ++i ) grid->mask[i] = TRUE;

  if ( gridInqMask(gridID, NULL) )
    {
      int *mask = (int*) Malloc(gridsize*sizeof(int));
      gridInqMask(gridID, mask);
      for ( long i = 0; i < gridsize; ++i )
	if ( mask[i] == 0 ) grid->mask[i] = FALSE;
      Free(mask);
    }

  if ( !remap_write_remap && grid->remap_grid_type == REMAP_GRID_TYPE_REG2D ) return;

  if ( !(gridInqXvals(gridID, NULL) && gridInqYvals(gridID, NULL)) )
    cdoAbort("%s grid cell center coordinates missing!", txt);

  gridInqXvals(gridID, grid->cell_center_lon);
  gridInqYvals(gridID, grid->cell_center_lat);

  char xunits[CDI_MAX_NAME]; xunits[0] = 0;
  char yunits[CDI_MAX_NAME]; yunits[0] = 0;
  cdiGridInqKeyStr(gridID, CDI_KEY_XUNITS, CDI_MAX_NAME, xunits);
  cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);

  if ( grid->lneed_cell_corners )
    {
      if ( gridInqXbounds(gridID, NULL) && gridInqYbounds(gridID, NULL) )
	{
	  gridInqXbounds(gridID, grid->cell_corner_lon);
	  gridInqYbounds(gridID, grid->cell_corner_lat);
	}
      else if ( lgrid_gen_bounds )
	{
	  grid_cell_center_to_bounds_X2D(xunits, grid->dims[0], grid->dims[1], grid->cell_center_lon, grid->cell_corner_lon, 0);
	  grid_cell_center_to_bounds_Y2D(yunits, grid->dims[0], grid->dims[1], grid->cell_center_lat, grid->cell_corner_lat);
	}
      else
	{
	  cdoAbort("%s grid cell corner coordinates missing!", txt);
	}
    }


  if ( gridInqType(grid->gridID) == GRID_GME ) gridInqMaskGME(gridID_gme, grid->vgpm);

  /* Convert lat/lon units if required */

  grid_to_radian(xunits, grid->size, grid->cell_center_lon, "grid center lon"); 
  grid_to_radian(yunits, grid->size, grid->cell_center_lat, "grid center lat"); 
  /* Note: using units from cell center instead from bounds */
  if ( grid->num_cell_corners && grid->lneed_cell_corners )
    {
      grid_to_radian(xunits, grid->num_cell_corners*grid->size, grid->cell_corner_lon, "grid corner lon"); 
      grid_to_radian(yunits, grid->num_cell_corners*grid->size, grid->cell_corner_lat, "grid corner lat"); 
    }

  if ( lgrid_destroy ) gridDestroy(gridID);

  /* Convert longitudes to 0,2pi interval */

  check_lon_range(grid->size, grid->cell_center_lon);

  if ( grid->num_cell_corners && grid->lneed_cell_corners )
    check_lon_range(grid->num_cell_corners*grid->size, grid->cell_corner_lon);

  /*  Make sure input latitude range is within the machine values for +/- pi/2 */

  check_lat_range(grid->size, grid->cell_center_lat);

  if ( grid->num_cell_corners && grid->lneed_cell_corners )
    check_lat_range(grid->num_cell_corners*grid->size, grid->cell_corner_lat);
}

/*  Compute bounding boxes for restricting future grid searches */
static
void cell_bounding_boxes(remapgrid_t *grid, int remap_grid_basis)
{
  if ( remap_grid_basis == REMAP_GRID_BASIS_SRC || grid->luse_cell_corners )
    grid->cell_bound_box = (restr_t*) Malloc(4*grid->size*sizeof(restr_t));

  if ( grid->luse_cell_corners )
    {
      if ( grid->lneed_cell_corners )
	{
	  if ( cdoVerbose ) cdoPrint("Grid: boundbox_from_corners");

	  boundbox_from_corners(grid->size, grid->num_cell_corners, 
				grid->cell_corner_lon, grid->cell_corner_lat, grid->cell_bound_box);
	}
      else /* full grid search */
	{
	  long gridsize;
	  long i, i4;
	  
	  gridsize = grid->size;
  
	  if ( cdoVerbose ) cdoPrint("Grid: bounds missing -> full grid search!");

	  for ( i = 0; i < gridsize; ++i )
	    {
	      i4 = i<<2;
	      grid->cell_bound_box[i4  ] = RESTR_SCALE(-PIH);
	      grid->cell_bound_box[i4+1] = RESTR_SCALE( PIH);
	      grid->cell_bound_box[i4+2] = RESTR_SCALE(0.);
	      grid->cell_bound_box[i4+3] = RESTR_SCALE(PI2);
	    }
	}
    }
  else if ( remap_grid_basis == REMAP_GRID_BASIS_SRC )
    {
      if ( grid->rank != 2 ) cdoAbort("Internal problem, grid rank = %d!", grid->rank);

      long nx = grid->dims[0];
      long ny = grid->dims[1];

      if ( cdoVerbose ) cdoPrint("Grid: boundbox_from_center");

      boundbox_from_center(grid->is_cyclic, grid->size, nx, ny, 
			   grid->cell_center_lon, grid->cell_center_lat, grid->cell_bound_box);
    }

  if ( remap_grid_basis == REMAP_GRID_BASIS_SRC || grid->lneed_cell_corners )
    check_lon_boundbox_range(grid->size, grid->cell_bound_box);

  /* Try to check for cells that overlap poles */

  if ( remap_grid_basis == REMAP_GRID_BASIS_SRC || grid->lneed_cell_corners )
    check_lat_boundbox_range(grid->size, grid->cell_bound_box, grid->cell_center_lat);
}


void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t *src_grid, int gridID2, remapgrid_t *tgt_grid)
{
  int reg2d_src_gridID = gridID1;
  int reg2d_tgt_gridID = gridID2;

  /* Initialize remapgrid structure */
  remapgrid_init(src_grid);
  remapgrid_init(tgt_grid);

  if ( map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_BICUBIC ||
       map_type == MAP_TYPE_DISTWGT  || map_type == MAP_TYPE_CONSERV_YAC )
    {
      if ( IS_REG2D_GRID(gridID1) ) src_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
      // src_grid->remap_grid_type = 0;
    }

  if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
    {
      if ( IS_REG2D_GRID(gridID2) && map_type == MAP_TYPE_CONSERV_YAC ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
      // else src_grid->remap_grid_type = -1;
    }

  if ( !remap_gen_weights && IS_REG2D_GRID(gridID2) && tgt_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
    {
      if ( map_type == MAP_TYPE_DISTWGT ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
      if ( map_type == MAP_TYPE_BILINEAR && src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
    }

  if ( lextrapolate )
    src_grid->lextrapolate = true;
  else
    src_grid->lextrapolate = false;

  if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC )
    {
      if ( src_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
	{
	  src_grid->luse_cell_corners  = true;
	  src_grid->lneed_cell_corners = true;
	}

      if ( tgt_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
	{
	  tgt_grid->luse_cell_corners  = true;
	  tgt_grid->lneed_cell_corners = true;
	}
    }

  src_grid->gridID = gridID1;
  tgt_grid->gridID = gridID2;

  if ( !src_grid->lextrapolate && gridInqSize(src_grid->gridID) > 1 &&
       map_type == MAP_TYPE_DISTWGT &&
       ((gridInqType(gridID1) == GRID_PROJECTION && gridInqProjType(gridID1) == CDI_PROJ_RLL) ||
	(gridInqType(gridID1) == GRID_LONLAT && src_grid->non_global)) )
    {
      src_grid->gridID = gridID1 = expand_lonlat_grid(gridID1);
      reg2d_src_gridID = gridID1;
    }

  if ( gridInqType(gridID1) == GRID_UNSTRUCTURED )
    {
      if ( gridInqYvals(gridID1, NULL) == 0 || gridInqXvals(gridID1, NULL) == 0 )
	{
	  if ( gridInqNumber(gridID1) > 0 )
	    {
	      src_grid->gridID = gridID1 = referenceToGrid(gridID1);
	      if ( gridID1 == -1 ) cdoAbort("Reference to source grid not found!");
	    }
	}