remaplib.c 52.3 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
/*
  This is a C library of the Fortran SCRIP version 1.4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
4
  ===>>> Please send bug reports to <http://mpimet.mpg.de/cdo> <<<===
Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

  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

*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
/*
41
  2013-11-08 Uwe Schulzweida: split remapgrid class to src_grid and tgt_grid
42
  2012-01-16 Uwe Schulzweida: alloc grid2_bound_box only for conservative remapping
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
43
  2011-01-07 Uwe Schulzweida: Changed remap weights from 2D to 1D array
44
  2009-05-25 Uwe Schulzweida: Changed restrict data type from double to int
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
  2009-01-11 Uwe Schulzweida: OpenMP parallelization
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
 */
Cedrick Ansorge's avatar
Cedrick Ansorge committed
47

Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
#if defined(HAVE_CONFIG_H)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
54
#  include "config.h"
#endif

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

Ralf Mueller's avatar
Ralf Mueller committed
55
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
#include "cdo.h"
#include "cdo_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
#include "grid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
#include "remap.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
#include "remap_store_link_cnsrv.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61

Uwe Schulzweida's avatar
Uwe Schulzweida committed
62

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

65

Uwe Schulzweida's avatar
Uwe Schulzweida committed
66

67
68
static bool remap_gen_weights     = true;
static bool remap_write_remap     = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
static int  remap_num_srch_bins   = 180;
70
#define  DEFAULT_MAX_ITER  100
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
long remap_max_iter        = DEFAULT_MAX_ITER;  /* Max iteration count for i, j iteration */
72

73
void remap_set_int(int remapvar, int value)
74
{
75
  if      ( remapvar == REMAP_STORE_LINK_FAST ) remap_store_link_fast = value;
76
  else if ( remapvar == REMAP_WRITE_REMAP     ) remap_write_remap     = value > 0;
77
  else if ( remapvar == REMAP_MAX_ITER        ) remap_max_iter        = value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
  else if ( remapvar == REMAP_NUM_SRCH_BINS   ) remap_num_srch_bins   = value;
79
  else if ( remapvar == REMAP_GENWEIGHTS      ) remap_gen_weights     = value > 0;
80
  else      cdoAbort("Unsupported remap variable (%d)!", remapvar);
81
82
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
double intlin(double x, double y1, double x1, double y2, double x2);

85
void remapGridFree(remapgrid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
{
87
88
  if ( grid->vgpm ) Free(grid->vgpm);
  if ( grid->mask ) Free(grid->mask);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89

90
91
92
93
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94

95
96
97
98
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99

100
101
  if ( grid->cell_area ) Free(grid->cell_area);
  if ( grid->cell_frac ) Free(grid->cell_frac);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102

103
  if ( grid->cell_bound_box ) Free(grid->cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104

105
106
  if ( grid->bin_addr ) Free(grid->bin_addr);
  if ( grid->bin_lats ) Free(grid->bin_lats);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
109
110
111

} /* remapGridFree */

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
void remapVarsFree(remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
{
114
  if ( rv->pinit )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
    {
116
117
      rv->pinit    = false;
      rv->sort_add = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

119
120
121
      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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122

123
      if ( rv->links.option )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
	{
125
	  rv->links.option = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126

127
	  if ( rv->links.num_blks )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
	    {
129
	      Free(rv->links.num_links);
130
131
	      long num_blks = rv->links.num_blks;
	      for ( long i = 0; i < num_blks; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
		{
133
134
135
		  Free(rv->links.src_add[i]);
		  Free(rv->links.dst_add[i]);
		  Free(rv->links.w_index[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
		}
137
138
139
	      Free(rv->links.src_add);
	      Free(rv->links.dst_add);
	      Free(rv->links.w_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
140
141
	    }
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
    }
143
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
    fprintf(stderr, "%s Warning: vars not initialized!\n", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146
147
148
149

} /* remapVarsFree */

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
150
void remapgrid_init(remapgrid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
153
154
155
  grid->remap_grid_type  = -1;
  grid->num_srch_bins    = remap_num_srch_bins; // only for source grid ?

  grid->num_cell_corners = 0;
156
157
  grid->luse_cell_corners  = false;
  grid->lneed_cell_corners = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158

159
160
161
162
  grid->nvgp             = 0;
  grid->vgpm             = NULL;

  grid->mask             = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163

164
165
  grid->reg2d_center_lon = NULL;
  grid->reg2d_center_lat = NULL;
166
167
  grid->reg2d_corner_lon = NULL;
  grid->reg2d_corner_lat = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168

169
170
171
172
173
  grid->cell_center_lon  = NULL;
  grid->cell_center_lat  = NULL;
  grid->cell_corner_lon  = NULL;
  grid->cell_corner_lat  = NULL;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
176
  grid->cell_area        = NULL;
  grid->cell_frac        = NULL;

177
178
179
180
  grid->cell_bound_box   = NULL;

  grid->bin_addr         = NULL;
  grid->bin_lats         = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
184
}

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

185
void remapgrid_alloc(int map_type, remapgrid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
{
187
  if ( grid->nvgp )
188
    grid->vgpm   = (int*) Malloc(grid->nvgp*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189

190
  grid->mask     = (int*) Malloc(grid->size*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
191

192
  if ( remap_write_remap || grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
    {
194
195
      grid->cell_center_lon = (double*) Malloc(grid->size*sizeof(double));
      grid->cell_center_lat = (double*) Malloc(grid->size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197

198
  if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
    {
200
      grid->cell_area = (double*) Malloc(grid->size*sizeof(double));
201
      memset(grid->cell_area, 0, grid->size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
203
    }

204
  grid->cell_frac = (double*) Malloc(grid->size*sizeof(double));
205
  memset(grid->cell_frac, 0, grid->size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206

207
  if ( grid->lneed_cell_corners )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
    {
209
210
      if ( grid->num_cell_corners > 0 )
        {
211
	  long nalloc = grid->num_cell_corners*grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212

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

216
	  grid->cell_corner_lat = (double*) Malloc(nalloc*sizeof(double));  
217
	  memset(grid->cell_corner_lat, 0, nalloc*sizeof(double));
218
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
221
    }
}

222
223
/*****************************************************************************/
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
225
void boundbox_from_corners(long size, long nc, const double *restrict corner_lon,
			   const double *restrict corner_lat, restr_t *restrict bound_box)
226
{
227
  long i4, inc, j;
228
  restr_t clon, clat;
229

Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
#pragma omp parallel for default(none)        \
  shared(bound_box, corner_lat, corner_lon, nc, size)	\
233
  private(i4, inc, j, clon, clat)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234
#endif
235
  for ( long i = 0; i < size; ++i )
236
    {
237
      i4 = i<<2; // *4
238
      inc = i*nc;
239
240
241
242
243
244
245
      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 )
246
	{
247
248
249
250
251
252
	  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;
253
254
255
256
257
	}
    }
}

static
258
void boundbox_from_center(bool lonIsCyclic, long size, long nx, long ny, const double *restrict center_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
			  const double *restrict center_lat, restr_t *restrict bound_box)
260
{
261
  long n4, i, j, k, ip1, jp1;
262
  long n_add, e_add, ne_add;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
  restr_t tmp_lats[4], tmp_lons[4];  /* temps for computing bounding boxes */
264

Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
#pragma omp parallel for default(none)        \
  shared(lonIsCyclic, size, nx, ny, center_lon, center_lat, bound_box)	\
268
  private(n4, i, j, k, ip1, jp1, n_add, e_add, ne_add, tmp_lats, tmp_lons)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
#endif
270
  for ( long n = 0; n < size; n++ )
271
    {
272
      n4 = n<<2;
273

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

279
      if ( i < (nx-1) )
280
281
282
	ip1 = i + 1;
      else
	{
283
	  /* 2009-01-09 Uwe Schulzweida: bug fix */
284
	  ip1 = lonIsCyclic ? 0 : i;
285
286
	}

287
      if ( j < (ny-1) )
288
289
290
	jp1 = j + 1;
      else
	{
291
292
	  /* 2008-12-17 Uwe Schulzweida: latitute cyclic ??? (bug fix) */
	  jp1 = j;
293
294
	}

295
296
297
      n_add  = jp1*nx + i;
      e_add  = j  *nx + ip1;
      ne_add = jp1*nx + ip1;
298
299
300

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
302
303
304
      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]);
305

Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
307
308
309
      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]);
310

Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
      bound_box[n4  ] = tmp_lats[0];
312
313
314
315
316
317
      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++ )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318
	  if ( tmp_lats[k] < bound_box[n4  ] ) bound_box[n4  ] = tmp_lats[k];
319
320
321
322
323
324
325
	  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];
	}
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344

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];
    }
}

345

Uwe Schulzweida's avatar
Uwe Schulzweida committed
346
347
void check_lon_range(long nlons, double *lons)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
349
  assert(lons != NULL);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
#pragma omp parallel for default(none) shared(nlons, lons)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
#endif
353
  for ( long n = 0; n < nlons; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
354
    {
355
356
357
358
      // remove missing values
      if ( lons[n] < -PI2 ) lons[n] = 0;
      if ( lons[n] > 2*PI2) lons[n] = PI2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
360
361
362
363
      if ( lons[n] > PI2  ) lons[n] -= PI2;
      if ( lons[n] < ZERO ) lons[n] += PI2;
    }
}

364

Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
366
void check_lat_range(long nlats, double *lats)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
367
368
  assert(lats != NULL);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
#pragma omp parallel for default(none) shared(nlats, lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
#endif
372
  for ( long n = 0; n < nlats; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
374
375
376
377
378
379
380
381
    {
      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)
{
382
  long n4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383

Uwe Schulzweida's avatar
Uwe Schulzweida committed
384
385
  assert(bound_box != NULL);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
#pragma omp parallel for default(none) shared(nlons, bound_box) private(n4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
388
#endif
389
  for ( long n = 0; n < nlons; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
    {
391
      n4 = n<<2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
393
      if ( RESTR_ABS(bound_box[n4+3] - bound_box[n4+2]) > RESTR_SCALE(PI) )
	{
394
	  bound_box[n4+2] = RESTR_SCALE(0.);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
397
398
399
400
	  bound_box[n4+3] = RESTR_SCALE(PI2);
	}
    }
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
401
void check_lat_boundbox_range(long nlats, restr_t *restrict bound_box, double *restrict lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
{
403
  long n4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404

Uwe Schulzweida's avatar
Uwe Schulzweida committed
405
406
  assert(bound_box != NULL);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
407
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408
#pragma omp parallel for default(none) shared(nlats, bound_box, lats) private(n4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
#endif
410
  for ( long n = 0; n < nlats; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
    {
412
      n4 = n<<2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
414
      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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
416
417
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
419
420
static
int expand_lonlat_grid(int gridID)
{
421
422
423
424
  long nx = gridInqXsize(gridID);
  long ny = gridInqYsize(gridID);
  long nxp4 = nx+4;
  long nyp4 = ny+4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425

426
427
  double *xvals = (double*) Malloc(nxp4*sizeof(double));
  double *yvals = (double*) Malloc(nyp4*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
429
430
  gridInqXvals(gridID, xvals+2);
  gridInqYvals(gridID, yvals+2);

431
432
  int gridtype = gridInqType(gridID);
  int gridIDnew = gridCreate(gridtype, nxp4*nyp4);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
434
  gridDefXsize(gridIDnew, nxp4);
  gridDefYsize(gridIDnew, nyp4);
435
436

  grid_copy_attributes(gridID, gridIDnew);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437

Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
439
440
441
442
443
444
445
446
447
448
449
450
  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);

451
452
  Free(xvals);
  Free(yvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453

454
455
456
  if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL )
    {
      double xpole, ypole, angle;
457
      gridInqParamRLL(gridID, &xpole, &ypole, &angle);
458
      gridDefParamRLL(gridIDnew, xpole, ypole, angle);
459
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460

461
  return gridIDnew;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
463
464
465
}

static
int expand_curvilinear_grid(int gridID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
466
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
  long i, j;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468

469
470
471
472
473
474
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475

476
477
  double *xvals = (double*) Malloc(gridsize_new*sizeof(double));
  double *yvals = (double*) Malloc(gridsize_new*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478
479
480
  gridInqXvals(gridID, xvals);
  gridInqYvals(gridID, yvals);

481
  int gridIDnew = gridCreate(GRID_CURVILINEAR, nxp4*nyp4);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
483
484
  gridDefXsize(gridIDnew, nxp4);
  gridDefYsize(gridIDnew, nyp4);

485
  grid_copy_attributes(gridID, gridIDnew);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533

  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);
  
534
535
  Free(xvals);
  Free(yvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536

537
  return gridIDnew;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
539
540
}

/*****************************************************************************/
541
542
543
544
545
546

static
void grid_check_lat_borders_rad(int n, double *ybounds)
{
  if ( ybounds[0] > ybounds[n-1] )
    {
547
548
      if ( RAD2DEG*ybounds[0]   >  PIH ) ybounds[0]   =  PIH;
      if ( RAD2DEG*ybounds[n-1] < -PIH ) ybounds[n-1] = -PIH;
549
550
551
    }
  else
    {
552
553
      if ( RAD2DEG*ybounds[0]   < -PIH ) ybounds[0]   = -PIH;
      if ( RAD2DEG*ybounds[n-1] >  PIH ) ybounds[n-1] =  PIH;
554
555
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556

557
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
void remap_define_reg2d(int gridID, remapgrid_t *grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559
{
560
561
  long nx = grid->dims[0];
  long ny = grid->dims[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562

563
564
  long nxp1 = nx + 1;
  long nyp1 = ny + 1;
565

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

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

571
572
  grid->reg2d_center_lon = (double*) Malloc(nxm*sizeof(double));
  grid->reg2d_center_lat = (double*) Malloc( ny*sizeof(double));
573
 
574
575
  grid->reg2d_center_lon[0] = 0;
  grid->reg2d_center_lat[0] = 0;
576
577
  gridInqXvals(gridID, grid->reg2d_center_lon);
  gridInqYvals(gridID, grid->reg2d_center_lat);
578

579
  /* Convert lat/lon units if required */
580

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

584
585
  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"); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586

587
  if ( grid->is_cyclic ) grid->reg2d_center_lon[nx] = grid->reg2d_center_lon[0] + PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588

589
590
  grid->reg2d_corner_lon = (double*) Malloc(nxp1*sizeof(double));
  grid->reg2d_corner_lat = (double*) Malloc(nyp1*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591

592
593
594
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
595

Uwe Schulzweida's avatar
Uwe Schulzweida committed
596
597
  //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]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598

599
}
600

601
static
602
void remap_define_grid(int map_type, int gridID, remapgrid_t *grid, const char *txt)
603
{
604
605
  bool lgrid_destroy = false;
  bool lgrid_gen_bounds = false;
606
  int gridID_gme = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607

608
  if ( gridInqType(grid->gridID) != GRID_UNSTRUCTURED && gridInqType(grid->gridID) != GRID_CURVILINEAR )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
609
    {
610
      if ( gridInqType(grid->gridID) == GRID_GME )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
	{
612
613
614
615
	  gridID_gme = gridToUnstructured(grid->gridID, 1);
	  grid->nvgp = gridInqSize(gridID_gme);
	  gridID = gridDuplicate(gridID_gme);
	  gridCompress(gridID);
616
	  grid->luse_cell_corners = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
	}
618
      else if ( remap_write_remap || grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
619
	{
620
	  lgrid_destroy = true;
621
	  gridID = gridToCurvilinear(grid->gridID, 1);
622
	  lgrid_gen_bounds = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
623
624
625
	}
    }

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

628
629
  grid->dims[0] = gridInqXsize(gridID);
  grid->dims[1] = gridInqYsize(gridID);
630
631
632
633
634
  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)));
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635

636
  grid->is_cyclic = (gridIsCircular(gridID) > 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637

638
  grid->rank = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? 1 : 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
639

640
  grid->num_cell_corners = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? gridInqNvertex(gridID) : 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
641

642
  remapgrid_alloc(map_type, grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643

Uwe Schulzweida's avatar
Uwe Schulzweida committed
644
645
646
647
648
  /* Initialize logical mask */

#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(gridsize, grid)
#endif
649
  for ( long i = 0; i < gridsize; ++i ) grid->mask[i] = TRUE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650

651
652
  if ( gridInqMask(gridID, NULL) )
    {
653
      int *mask = (int*) Malloc(gridsize*sizeof(int));
654
      gridInqMask(gridID, mask);
655
      for ( long i = 0; i < gridsize; ++i )
656
	if ( mask[i] == 0 ) grid->mask[i] = FALSE;
657
      Free(mask);
658
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
659

660
  if ( !remap_write_remap && grid->remap_grid_type == REMAP_GRID_TYPE_REG2D ) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
661

662
  if ( !(gridInqXvals(gridID, NULL) && gridInqYvals(gridID, NULL)) )
663
    cdoAbort("%s grid cell center coordinates missing!", txt);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
664

665
666
  gridInqXvals(gridID, grid->cell_center_lon);
  gridInqYvals(gridID, grid->cell_center_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
667

668
669
670
671
  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);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672

673
  if ( grid->lneed_cell_corners )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
    {
675
      if ( gridInqXbounds(gridID, NULL) && gridInqYbounds(gridID, NULL) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
	{
677
678
	  gridInqXbounds(gridID, grid->cell_corner_lon);
	  gridInqYbounds(gridID, grid->cell_corner_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
679
	}
680
      else if ( lgrid_gen_bounds )
681
	{
682
683
	  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);
684
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
685
686
      else
	{
687
	  cdoAbort("%s grid cell corner coordinates missing!", txt);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
688
689
690
691
	}
    }


692
  if ( gridInqType(grid->gridID) == GRID_GME ) gridInqMaskGME(gridID_gme, grid->vgpm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693

Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
  /* Convert lat/lon units if required */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
695

696
697
  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"); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
698
  /* Note: using units from cell center instead from bounds */
699
  if ( grid->num_cell_corners && grid->lneed_cell_corners )
700
    {
701
702
      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"); 
703
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
704

705
  if ( lgrid_destroy ) gridDestroy(gridID);
706

707
  /* Convert longitudes to 0,2pi interval */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
708

709
  check_lon_range(grid->size, grid->cell_center_lon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
710

711
  if ( grid->num_cell_corners && grid->lneed_cell_corners )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
712
713
    check_lon_range(grid->num_cell_corners*grid->size, grid->cell_corner_lon);

714
715
716
717
718
719
720
721
722
723
724
725
  /*  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)
{
726
  if ( remap_grid_basis == REMAP_GRID_BASIS_SRC || grid->luse_cell_corners )
727
    grid->cell_bound_box = (restr_t*) Malloc(4*grid->size*sizeof(restr_t));
728
729

  if ( grid->luse_cell_corners )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730
    {
731
      if ( grid->lneed_cell_corners )
732
	{
733
734
735
736
	  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);
737
	}
738
      else /* full grid search */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
	{
740
741
742
743
744
745
	  long gridsize;
	  long i, i4;
	  
	  gridsize = grid->size;
  
	  if ( cdoVerbose ) cdoPrint("Grid: bounds missing -> full grid search!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
746

747
748
749
750
751
	  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);
752
	      grid->cell_bound_box[i4+2] = RESTR_SCALE(0.);
753
754
	      grid->cell_bound_box[i4+3] = RESTR_SCALE(PI2);
	    }
755
756
	}
    }
757
  else if ( remap_grid_basis == REMAP_GRID_BASIS_SRC )
758
    {
759
      if ( grid->rank != 2 ) cdoAbort("Internal problem, grid rank = %d!", grid->rank);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
760

761
762
      long nx = grid->dims[0];
      long ny = grid->dims[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
763

764
      if ( cdoVerbose ) cdoPrint("Grid: boundbox_from_center");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
765

766
767
      boundbox_from_center(grid->is_cyclic, grid->size, nx, ny, 
			   grid->cell_center_lon, grid->cell_center_lat, grid->cell_bound_box);
768
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
769

770
771
  if ( remap_grid_basis == REMAP_GRID_BASIS_SRC || grid->lneed_cell_corners )
    check_lon_boundbox_range(grid->size, grid->cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
772

773
  /* Try to check for cells that overlap poles */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
774

775
776
777
  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);
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778

Uwe Schulzweida's avatar
Uwe Schulzweida committed
779

780
void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t *src_grid, int gridID2, remapgrid_t *tgt_grid)
781
782
{
  int reg2d_src_gridID = gridID1;
783
  int reg2d_tgt_gridID = gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
784

Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
786
787
  /* Initialize remapgrid structure */
  remapgrid_init(src_grid);
  remapgrid_init(tgt_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
788

Uwe Schulzweida's avatar
Uwe Schulzweida committed
789
  if ( map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_BICUBIC ||
790
       map_type == MAP_TYPE_DISTWGT  || map_type == MAP_TYPE_CONSERV_YAC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
791
    {
792
      if ( IS_REG2D_GRID(gridID1) ) src_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
793
794
      // src_grid->remap_grid_type = 0;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
795

Uwe Schulzweida's avatar
Uwe Schulzweida committed
796
  if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
797
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
798
      if ( IS_REG2D_GRID(gridID2) && map_type == MAP_TYPE_CONSERV_YAC ) tgt_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
799
800
801
      // else src_grid->remap_grid_type = -1;
    }

802
  if ( !remap_gen_weights && IS_REG2D_GRID(gridID2) && tgt_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
803
    {
804
805
      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;
806
807
    }

808
809
  if ( lextrapolate )
    src_grid->lextrapolate = true;
810
  else
811
    src_grid->lextrapolate = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812

813
  if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSERV_YAC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
814
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
815
      if ( src_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
816
	{
817
818
	  src_grid->luse_cell_corners  = true;
	  src_grid->lneed_cell_corners = true;
819
	}
820
821
822

      if ( tgt_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
	{
823
824
	  tgt_grid->luse_cell_corners  = true;
	  tgt_grid->lneed_cell_corners = true;
825
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
827
    }

828
829
  src_grid->gridID = gridID1;
  tgt_grid->gridID = gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830

831
  if ( !src_grid->lextrapolate && gridInqSize(src_grid->gridID) > 1 &&
832
       map_type == MAP_TYPE_DISTWGT &&
833
       ((gridInqType(gridID1) == GRID_PROJECTION && gridInqProjType(gridID1) == CDI_PROJ_RLL) ||
834
835
836
	(gridInqType(gridID1) == GRID_LONLAT && src_grid->non_global)) )
    {
      src_grid->gridID = gridID1 = expand_lonlat_grid(gridID1);
837
      reg2d_src_gridID = gridID1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
838
839
    }

840
  if ( gridInqType(gridID1) == GRID_UNSTRUCTURED )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
841
    {
842
      if ( gridInqYvals(gridID1, NULL) == 0 || gridInqXvals(gridID1, NULL) == 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843
	{
844
845
846
847
848
	  if ( gridInqNumber(gridID1) > 0 )
	    {
	      src_grid->gridID = gridID1 = referenceToGrid(gridID1);
	      if ( gridID1 == -1 ) cdoAbort("Reference to source grid not found!");
	    }
849
	}