remap_conserv.cc 35.4 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
#include "remap_store_link.h"
6

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
extern "C" {
8
9
10
#include "clipping/clipping.h"
#include "clipping/area.h"
#include "clipping/geometry.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
11
}
12

13
//#define STIMER
Uwe Schulzweida's avatar
Uwe Schulzweida committed
14
15
16
17
18

#ifdef STIMER
#include <time.h>
#endif

19
20
typedef struct {
  enum yac_edge_type *src_edge_type;
21
  size_t srch_corners;
22
  size_t max_srch_cells;
23
24
25
26
27
28
29
  double *partial_areas;
  double *partial_weights;
  struct grid_cell *src_grid_cells;
  struct grid_cell *overlap_buffer;
} search_t;

static
30
void search_realloc(size_t num_srch_cells, search_t *search)
31
{
32
  size_t max_srch_cells = search->max_srch_cells;
33
34
35
36
37
38
39
  double *partial_areas   = search->partial_areas;
  double *partial_weights = search->partial_weights;
  struct grid_cell *overlap_buffer = search->overlap_buffer;
  struct grid_cell *src_grid_cells = search->src_grid_cells;

  if ( num_srch_cells > max_srch_cells )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
41
42
43
      partial_areas   = (double*) Realloc(partial_areas,   num_srch_cells*sizeof(double));
      partial_weights = (double*) Realloc(partial_weights, num_srch_cells*sizeof(double));
      overlap_buffer = (struct grid_cell*) Realloc(overlap_buffer, num_srch_cells*sizeof(struct grid_cell));
      src_grid_cells = (struct grid_cell*) Realloc(src_grid_cells, num_srch_cells*sizeof(struct grid_cell));
44

45
      for ( size_t n = max_srch_cells; n < num_srch_cells; ++n )
46
47
48
49
50
51
52
53
54
        {
          overlap_buffer[n].array_size      = 0;
          overlap_buffer[n].num_corners     = 0;
          overlap_buffer[n].edge_type       = NULL;
          overlap_buffer[n].coordinates_x   = NULL;
          overlap_buffer[n].coordinates_y   = NULL;
          overlap_buffer[n].coordinates_xyz = NULL;
        }

55
      for ( size_t n = max_srch_cells; n < num_srch_cells; ++n )
56
57
58
59
        {
          src_grid_cells[n].array_size      = search->srch_corners;
          src_grid_cells[n].num_corners     = search->srch_corners;
          src_grid_cells[n].edge_type       = search->src_edge_type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
61
62
          src_grid_cells[n].coordinates_x   = (double*) Malloc(search->srch_corners*sizeof(double));
          src_grid_cells[n].coordinates_y   = (double*) Malloc(search->srch_corners*sizeof(double));
          src_grid_cells[n].coordinates_xyz = (double*) Malloc(3*search->srch_corners*sizeof(double));
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
        }

      max_srch_cells = num_srch_cells;

      search->max_srch_cells  = max_srch_cells;
      search->partial_areas   = partial_areas;
      search->partial_weights = partial_weights;
      search->overlap_buffer  = overlap_buffer;
      search->src_grid_cells  = src_grid_cells;
    }
}

static
void search_free(search_t *search)
{
78
  size_t max_srch_cells  = search->max_srch_cells;
79
80
81
82
83
  double *partial_areas   = search->partial_areas;
  double *partial_weights = search->partial_weights;
  struct grid_cell *overlap_buffer = search->overlap_buffer;
  struct grid_cell *src_grid_cells = search->src_grid_cells;

84
  for ( size_t n = 0; n < max_srch_cells; n++ )
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    {
      if ( overlap_buffer[n].array_size > 0 )
        {
          Free(overlap_buffer[n].coordinates_x);
          Free(overlap_buffer[n].coordinates_y);
          if ( overlap_buffer[n].coordinates_xyz ) Free(overlap_buffer[n].coordinates_xyz);
          if ( overlap_buffer[n].edge_type ) Free(overlap_buffer[n].edge_type);
        }
      
      Free(src_grid_cells[n].coordinates_x);
      Free(src_grid_cells[n].coordinates_y);
      Free(src_grid_cells[n].coordinates_xyz);
    }

  Free(partial_areas);
  Free(partial_weights);
}

103

Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
106
int rect_grid_search2(long *imin, long *imax, double xmin, double xmax, long nxm, const double *restrict xm);

static
107
size_t get_srch_cells_reg2d(const size_t *restrict src_grid_dims, 
108
                            const double *restrict src_corner_lat, const double *restrict src_corner_lon,
109
                            const double *restrict tgt_cell_bound_box, size_t *srch_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
{
111
  bool debug = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
113
  long nx = src_grid_dims[0];
  long ny = src_grid_dims[1];
114
  size_t num_srch_cells = 0;  // num cells in restricted search arrays
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115

116
117
  long nxp1 = nx+1;
  long nyp1 = ny+1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

119
120
  double src_lon_min = src_corner_lon[0];
  double src_lon_max = src_corner_lon[nx];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
123

  long imin = nxp1, imax = -1, jmin = nyp1, jmax = -1;

124
125
126
  int lfound = rect_grid_search2(&jmin, &jmax, tgt_cell_bound_box[0], tgt_cell_bound_box[1], nyp1, src_corner_lat);
  if ( !lfound ) return 0;
  // printf("lfound, jmin, jmax %d %ld %ld\n", lfound, jmin, jmax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
  // if ( jmin > 0 ) jmin--;
  // if ( jmax < (ny-2) ) jmax++;
129
130
131

  double bound_lon1 = tgt_cell_bound_box[2];
  double bound_lon2 = tgt_cell_bound_box[3];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
133
134
135
136
137
138
139
140
141
  if ( bound_lon1 <= src_lon_max && bound_lon2 >= src_lon_min )
    {
      if ( debug ) printf("  b1 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
      if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
      if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
      lfound = rect_grid_search2(&imin, &imax, bound_lon1, bound_lon2, nxp1, src_corner_lon);
      if ( lfound )
	{
	  if ( debug )
	    printf("   %g %g imin %ld  imax %ld  jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin], RAD2DEG*src_corner_lon[imax+1], imin, imax, jmin, jmax);
142
143
	  for ( long jm = jmin; jm <= jmax; ++jm )
	    for ( long im = imin; im <= imax; ++im )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
	      srch_add[num_srch_cells++] = jm*nx + im;
	}
    }

  bound_lon1 = tgt_cell_bound_box[2];
  bound_lon2 = tgt_cell_bound_box[3];
  if ( bound_lon1 <= src_lon_min && bound_lon2 >= src_lon_min )
    {
      long imin2 = nxp1, imax2 = -1;
      bound_lon1 += 2*M_PI;
      bound_lon2 += 2*M_PI;
      if ( debug ) printf("  b2 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
      if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
      if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
      lfound = rect_grid_search2(&imin2, &imax2, bound_lon1, bound_lon2, nxp1, src_corner_lon);
      if ( lfound )
	{
	  if ( imax != -1 && imin2 <= imax ) imin2 = imax+1;
	  if ( imax != -1 && imax2 <= imax ) imax2 = imax+1;
	  if ( imin2 >= 0 && imax2 < nxp1 )
	    {
	      if ( debug )
		printf("   %g %g imin %ld  imax %ld  jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin2], RAD2DEG*src_corner_lon[imax2+1], imin2, imax2, jmin, jmax);
167
168
	      for ( long jm = jmin; jm <= jmax; ++jm )
		for ( long im = imin2; im <= imax2; ++im )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
		  srch_add[num_srch_cells++] = jm*nx + im;
	    }
	}
    }

  bound_lon1 = tgt_cell_bound_box[2];
  bound_lon2 = tgt_cell_bound_box[3];
  if ( bound_lon1 <= src_lon_max && bound_lon2 >= src_lon_max )
    {
      long imin3 = nxp1, imax3 = -1;
      bound_lon1 -= 2*M_PI;
      bound_lon2 -= 2*M_PI;
      if ( debug ) printf("  b3 %g %g\n", bound_lon1*RAD2DEG, bound_lon2*RAD2DEG);
      if ( bound_lon1 < src_lon_min && bound_lon2 > src_lon_min ) bound_lon1 = src_lon_min;
      if ( bound_lon2 > src_lon_max && bound_lon1 < src_lon_max ) bound_lon2 = src_lon_max;
      lfound = rect_grid_search2(&imin3, &imax3, bound_lon1, bound_lon2, nxp1, src_corner_lon);
      if ( lfound )
	{
	  if ( imin != nxp1 && imin3 >= imin ) imin3 = imin-1;
	  if ( imax != nxp1 && imax3 >= imin ) imax3 = imin-1;
	  if ( imin3 >= 0 && imin3 < nxp1 )
	    {
	      if ( debug )
		printf("   %g %g imin %ld  imax %ld  jmin %ld jmax %ld\n", RAD2DEG*src_corner_lon[imin3], RAD2DEG*src_corner_lon[imax3+1], imin3, imax3, jmin, jmax);
193
194
	      for ( long jm = jmin; jm <= jmax; ++jm )
		for ( long im = imin3; im <= imax3; ++im )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
196
197
198
199
		  srch_add[num_srch_cells++] = jm*nx + im;
	    }
	}
    }

200
  if ( debug ) printf(" -> num_srch_cells: %zu\n", num_srch_cells);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201

Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
  return num_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
}

static
void restrict_boundbox(const double *restrict grid_bound_box, double *restrict bound_box)
{
  if ( bound_box[0] < grid_bound_box[0] && bound_box[1] > grid_bound_box[0] ) bound_box[0] = grid_bound_box[0];
  if ( bound_box[1] > grid_bound_box[1] && bound_box[0] < grid_bound_box[1] ) bound_box[1] = grid_bound_box[1];

  if ( bound_box[2] >= grid_bound_box[3] && (bound_box[3]-2*M_PI) > grid_bound_box[2] ) { bound_box[2] -= 2*M_PI; bound_box[3] -= 2*M_PI; }
  if ( bound_box[3] <= grid_bound_box[2] && (bound_box[2]-2*M_PI) < grid_bound_box[3] ) { bound_box[2] += 2*M_PI; bound_box[3] += 2*M_PI; }
  //  if ( bound_box[2] < grid_bound_box[2] && bound_box[3] > grid_bound_box[2] ) bound_box[2] = grid_bound_box[2];
  //  if ( bound_box[3] > grid_bound_box[3] && bound_box[2] < grid_bound_box[3] ) bound_box[3] = grid_bound_box[3];
}

static
218
void boundbox_from_corners_reg2d(size_t grid_add, const size_t *restrict grid_dims, const double *restrict corner_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
				 const double *restrict corner_lat, double *restrict bound_box)
{
221
222
223
  size_t nx = grid_dims[0];
  size_t iy = grid_add/nx;
  size_t ix = grid_add - iy*nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224

225
226
  double clat1 = corner_lat[iy  ];
  double clat2 = corner_lat[iy+1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243

  if ( clat2 > clat1 )
    {
      bound_box[0] = clat1;
      bound_box[1] = clat2;
    }
  else
    {
      bound_box[0] = clat2;
      bound_box[1] = clat1;
    }

  bound_box[2] = corner_lon[ix  ];
  bound_box[3] = corner_lon[ix+1];
}

static
244
void boundbox_from_corners1(size_t ic, size_t nc, const double *restrict corner_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
246
			    const double *restrict corner_lat, double *restrict bound_box)
{
247
  size_t inc = ic*nc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248

249
250
  double clat = corner_lat[inc];
  double clon = corner_lon[inc];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
252
253
254
255
256

  bound_box[0] = clat;
  bound_box[1] = clat;
  bound_box[2] = clon;
  bound_box[3] = clon;

257
  for ( size_t j = 1; j < nc; ++j )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    {
      clat = corner_lat[inc+j];
      clon = corner_lon[inc+j];

      if ( clat < bound_box[0] ) bound_box[0] = clat;
      if ( clat > bound_box[1] ) bound_box[1] = clat;
      if ( clon < bound_box[2] ) bound_box[2] = clon;
      if ( clon > bound_box[3] ) bound_box[3] = clon;
    }

  if ( fabs(bound_box[3] - bound_box[2]) > PI )
    {
      bound_box[2] = 0;
      bound_box[3] = PI2;
    }

  /*
  double dlon = fabs(bound_box[3] - bound_box[2]);

  if ( dlon > PI )
    {
      if ( bound_box[3] > bound_box[2] && (bound_box[3]-PI2) < 0. )
	{
	  double tmp = bound_box[2];
	  bound_box[2] = bound_box[3] - PI2;
	  bound_box[3] = tmp;
	}
      else
	{
	  bound_box[2] = 0;
	  bound_box[3] = PI2;
	}
    }
  */
}

static
295
void boundbox_from_corners1r(size_t ic, size_t nc, const double *restrict corner_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
297
			     const double *restrict corner_lat, restr_t *restrict bound_box)
{
298
  size_t inc = ic*nc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299

300
301
  restr_t clat = RESTR_SCALE(corner_lat[inc]);
  restr_t clon = RESTR_SCALE(corner_lon[inc]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
303
304
305
306
307

  bound_box[0] = clat;
  bound_box[1] = clat;
  bound_box[2] = clon;
  bound_box[3] = clon;

308
  for ( size_t j = 1; j < nc; ++j )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312
313
314
315
316
317
318
    {
      clat = RESTR_SCALE(corner_lat[inc+j]);
      clon = RESTR_SCALE(corner_lon[inc+j]);

      if ( clat < bound_box[0] ) bound_box[0] = clat;
      if ( clat > bound_box[1] ) bound_box[1] = clat;
      if ( clon < bound_box[2] ) bound_box[2] = clon;
      if ( clon > bound_box[3] ) bound_box[3] = clon;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
  if ( RESTR_ABS(bound_box[3] - bound_box[2]) > RESTR_SCALE(PI) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
    {
      bound_box[2] = 0;
      bound_box[3] = RESTR_SCALE(PI2);
    }
  /*
  if ( RESTR_ABS(bound_box[3] - bound_box[2]) > RESTR_SCALE(PI) )
    {
      if ( bound_box[3] > bound_box[2] && (bound_box[3]-RESTR_SCALE(PI2)) < RESTR_SCALE(0.) )
	{
	  restr_t tmp = bound_box[2];
	  bound_box[2] = bound_box[3] - RESTR_SCALE(PI2);
	  bound_box[3] = tmp;
	}
    }
  */
}

//#if defined(HAVE_LIBYAC)

339
340
341
static
double gridcell_area(struct grid_cell cell)
{
342
  return yac_huiliers_area(cell);
343
344
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
static
346
void cdo_compute_overlap_areas(size_t N, search_t *search, struct grid_cell target_cell)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
{
348
349
350
351
  double *partial_areas   = search->partial_areas;
  struct grid_cell *overlap_buffer = search->overlap_buffer;
  struct grid_cell *source_cells = search->src_grid_cells;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
353
  /* Do the clipping and get the cell for the overlapping area */

354
  yac_cell_clipping(N, source_cells, target_cell, overlap_buffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355
356
357

  /* Get the partial areas for the overlapping regions */

358
  for ( size_t n = 0; n < N; n++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
    {
360
      partial_areas[n] = gridcell_area(overlap_buffer[n]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
362
363
    }

#ifdef VERBOSE
364
  for ( size_t n = 0; n < N; n++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
366
367
    printf("overlap area : %lf\n", partial_areas[n]);
#endif
}
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409

static double const tol = 1.0e-12;

enum cell_type {
  UNDEF_CELL,
  LON_LAT_CELL,
  LAT_CELL,
  GREAT_CIRCLE_CELL,
  MIXED_CELL
};
/*
static enum cell_type get_cell_type(struct grid_cell target_cell) {

  int count_lat_edges = 0, count_great_circle_edges = 0;

   if ((target_cell.num_corners == 4) &&
       ((target_cell.edge_type[0] == LAT_CIRCLE &&
         target_cell.edge_type[1] == LON_CIRCLE &&
         target_cell.edge_type[2] == LAT_CIRCLE &&
         target_cell.edge_type[3] == LON_CIRCLE) ||
        (target_cell.edge_type[0] == LON_CIRCLE &&
         target_cell.edge_type[1] == LAT_CIRCLE &&
         target_cell.edge_type[2] == LON_CIRCLE &&
         target_cell.edge_type[3] == LAT_CIRCLE)))
      return LON_LAT_CELL;
   else
      for (unsigned i = 0; i < target_cell.num_corners; ++i)
         if (target_cell.edge_type[i] == LON_CIRCLE ||
             target_cell.edge_type[i] == GREAT_CIRCLE)
            count_great_circle_edges++;
         else
            count_lat_edges++;

   if (count_lat_edges && count_great_circle_edges)
      return MIXED_CELL;
   else if (count_lat_edges)
      return LAT_CELL;
   else
      return GREAT_CIRCLE_CELL;
}
*/
static
410
void cdo_compute_concave_overlap_areas(size_t N, search_t *search, struct grid_cell target_cell,
411
				       double target_node_x, double target_node_y)
412
{
413
414
415
416
  double *partial_areas   = search->partial_areas;
  struct grid_cell *overlap_buffer = search->overlap_buffer;
  struct grid_cell *source_cell = search->src_grid_cells;

417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
  /*
  enum cell_type target_cell_type = UNDEF_CELL;

  if ( target_cell.num_corners > 3 )
    target_cell_type = get_cell_type(target_cell);

  if ( target_cell.num_corners < 4 || target_cell_type == LON_LAT_CELL )
    {
      cdo_compute_overlap_areas(N, overlap_buffer, source_cell, target_cell, partial_areas);
      return;
    }

  if ( target_node_x == NULL || target_node_y == NULL )
    cdoAbort("Internal problem (cdo_compute_concave_overlap_areas): missing target point coordinates!");
  */
432
  /*
433
434
435
436
  struct grid_cell target_partial_cell =
    {.coordinates_x   = (double[3]){-1, -1, -1},
     .coordinates_y   = (double[3]){-1, -1, -1},
     .coordinates_xyz = (double[3*3]){-1, -1, -1},
437
     .edge_type       = (enum yac_edge_type[3]) {GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE},
438
     .num_corners     = 3};
439
440
441
442
  */
  double coordinates_x[3] = {-1, -1, -1};
  double coordinates_y[3] = {-1, -1, -1};
  double coordinates_xyz[9] = {-1, -1, -1};
443
  enum yac_edge_type edge_types[3] = {GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
445
446
447
448
449
  struct grid_cell target_partial_cell;
  target_partial_cell.coordinates_x   = coordinates_x;
  target_partial_cell.coordinates_y   = coordinates_y;
  target_partial_cell.coordinates_xyz = coordinates_xyz;
  target_partial_cell.edge_type       = edge_types;
  target_partial_cell.num_corners     = 3;
450
451
452

  /* Do the clipping and get the cell for the overlapping area */

453
  for ( size_t n = 0; n < N; n++) partial_areas[n] = 0.0;
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504

  // common node point to all partial target cells
  target_partial_cell.coordinates_x[0] = target_node_x;
  target_partial_cell.coordinates_y[0] = target_node_y;

  LLtoXYZ ( target_node_x, target_node_y, target_partial_cell.coordinates_xyz );

  for ( unsigned num_corners = 0; num_corners < target_cell.num_corners; ++num_corners )
    {
      unsigned corner_a = num_corners;
      unsigned corner_b = (num_corners+1)%target_cell.num_corners;

      // skip clipping and area calculation for degenerated triangles
      //
      // If this is not sufficient, instead we can try something like:
      //
      //     struct point_list target_list
      //     init_point_list(&target_list);
      //     generate_point_list(&target_list, target_cell);
      //     struct grid_cell temp_target_cell;
      //     generate_overlap_cell(target_list, temp_target_cell);
      //     free_point_list(&target_list);
      //
      // and use temp_target_cell for triangulation.
      //
      // Compared to the if statement below the alternative seems
      // to be quite costly.

      if ( ( ( fabs(target_cell.coordinates_xyz[0+3*corner_a]-target_cell.coordinates_xyz[0+3*corner_b]) < tol ) &&
	     ( fabs(target_cell.coordinates_xyz[1+3*corner_a]-target_cell.coordinates_xyz[1+3*corner_b]) < tol ) &&
	     ( fabs(target_cell.coordinates_xyz[2+3*corner_a]-target_cell.coordinates_xyz[2+3*corner_b]) < tol ) ) ||
	   ( ( fabs(target_cell.coordinates_xyz[0+3*corner_a]-target_partial_cell.coordinates_xyz[0]) < tol    ) &&
	     ( fabs(target_cell.coordinates_xyz[1+3*corner_a]-target_partial_cell.coordinates_xyz[1]) < tol    ) &&
	     ( fabs(target_cell.coordinates_xyz[2+3*corner_a]-target_partial_cell.coordinates_xyz[2]) < tol    ) ) ||
	   ( ( fabs(target_cell.coordinates_xyz[0+3*corner_b]-target_partial_cell.coordinates_xyz[0]) < tol    ) &&
	     ( fabs(target_cell.coordinates_xyz[1+3*corner_b]-target_partial_cell.coordinates_xyz[1]) < tol    ) &&
	     ( fabs(target_cell.coordinates_xyz[2+3*corner_b]-target_partial_cell.coordinates_xyz[2]) < tol    ) ) )
	continue;

      target_partial_cell.coordinates_x[1] = target_cell.coordinates_x[corner_a];
      target_partial_cell.coordinates_y[1] = target_cell.coordinates_y[corner_a];
      target_partial_cell.coordinates_x[2] = target_cell.coordinates_x[corner_b];
      target_partial_cell.coordinates_y[2] = target_cell.coordinates_y[corner_b];

      target_partial_cell.coordinates_xyz[0+3*1] = target_cell.coordinates_xyz[0+3*corner_a];
      target_partial_cell.coordinates_xyz[1+3*1] = target_cell.coordinates_xyz[1+3*corner_a];
      target_partial_cell.coordinates_xyz[2+3*1] = target_cell.coordinates_xyz[2+3*corner_a];
      target_partial_cell.coordinates_xyz[0+3*2] = target_cell.coordinates_xyz[0+3*corner_b];
      target_partial_cell.coordinates_xyz[1+3*2] = target_cell.coordinates_xyz[1+3*corner_b];
      target_partial_cell.coordinates_xyz[2+3*2] = target_cell.coordinates_xyz[2+3*corner_b];

505
      yac_cell_clipping((unsigned)N, source_cell, target_partial_cell, overlap_buffer);
506
507
508

      /* Get the partial areas for the overlapping regions as sum over the partial target cells. */

509
      for (size_t n = 0; n < N; n++)
510
	{
511
	  partial_areas[n] += gridcell_area(overlap_buffer[n]);
512
513
514
515
516
517
518
	  // we cannot use pole_area because it is rather inaccurate for great circle
	  // edges that are nearly circles of longitude
	  //partial_areas[n] = pole_area (overlap_buffer[n]);
	}
    }

#ifdef VERBOSE
519
520
  for (size_t n = 0; n < N; n++)
    printf("overlap area %zu: %lf\n", n, partial_areas[n]);
521
522
523
#endif
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
//#endif

static
int get_lonlat_circle_index(remapgrid_t *remap_grid)
{
  int lonlat_circle_index = -1;

  if ( remap_grid->num_cell_corners == 4 )
    {
      if ( remap_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  lonlat_circle_index = 1;
 	}
      else
	{
	  const double* cell_corner_lon = remap_grid->cell_corner_lon;
	  const double* cell_corner_lat = remap_grid->cell_corner_lat;
541
542
543
	  size_t gridsize = remap_grid->size;
	  size_t num_i = 0, num_eq0 = 0, num_eq1 = 0;
	  long iadd = gridsize/3-1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
544
545
546

	  if ( iadd == 0 ) iadd++;

547
	  for ( size_t i = 0; i < gridsize; i += iadd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
	    {
	      num_i++;

	      if ( IS_EQUAL(cell_corner_lon[i*4+1], cell_corner_lon[i*4+2]) &&
		   IS_EQUAL(cell_corner_lon[i*4+3], cell_corner_lon[i*4+0]) &&
		   IS_EQUAL(cell_corner_lat[i*4+0], cell_corner_lat[i*4+1]) &&
		   IS_EQUAL(cell_corner_lat[i*4+2], cell_corner_lat[i*4+3]) )
		{  
		  num_eq1++;
		}
	      else if ( IS_EQUAL(cell_corner_lon[i*4+0], cell_corner_lon[i*4+1]) &&
			IS_EQUAL(cell_corner_lon[i*4+2], cell_corner_lon[i*4+3]) &&
			IS_EQUAL(cell_corner_lat[i*4+1], cell_corner_lat[i*4+2]) &&
			IS_EQUAL(cell_corner_lat[i*4+3], cell_corner_lat[i*4+0]) )
		{
		  num_eq0++;
		}
	    }

	  if ( num_i == num_eq1 ) lonlat_circle_index = 1;
	  if ( num_i == num_eq0 ) lonlat_circle_index = 0;	      
	}
    }

  //printf("lonlat_circle_index %d\n", lonlat_circle_index);

574
  return lonlat_circle_index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
576
577
}


578
static
579
void remapNormalizeWeights(remapgrid_t *tgt_grid, remapvars_t *rv)
580
{
581
  // Include centroids in weights and normalize using destination area if requested
582
583
584
  size_t num_links = rv->num_links;
  size_t num_wts = rv->num_wts;
  size_t tgt_cell_add;       // current linear address for target grid cell
585
  double norm_factor = 0; // factor for normalizing wts
586
587
588
589
590
591
592
593

  if ( rv->norm_opt == NORM_OPT_DESTAREA )
    {
#if defined(SX)
#pragma vdir nodep
#endif
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
594
  shared(num_wts, num_links, rv, tgt_grid) \
595
  private(tgt_cell_add, norm_factor)
596
#endif
597
      for ( size_t n = 0; n < num_links; ++n )
598
	{
599
	  tgt_cell_add = rv->tgt_cell_add[n];
600

601
          if ( IS_NOT_EQUAL(tgt_grid->cell_area[tgt_cell_add], 0) )
602
	    norm_factor = 1./tgt_grid->cell_area[tgt_cell_add];
603
          else
604
            norm_factor = 0.;
605
606
607
608
609
610
611
612
613
614
615

	  rv->wts[n*num_wts] *= norm_factor;
	}
    }
  else if ( rv->norm_opt == NORM_OPT_FRACAREA )
    {
#if defined(SX)
#pragma vdir nodep
#endif
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
616
  shared(num_wts, num_links, rv, tgt_grid) \
617
  private(tgt_cell_add, norm_factor)
618
#endif
619
      for ( size_t n = 0; n < num_links; ++n )
620
	{
621
	  tgt_cell_add = rv->tgt_cell_add[n];
622

623
          if ( IS_NOT_EQUAL(tgt_grid->cell_frac[tgt_cell_add], 0) )
624
	    norm_factor = 1./tgt_grid->cell_frac[tgt_cell_add];
625
          else
626
            norm_factor = 0.;
627
628
629
630
631
632
633
634
635

	  rv->wts[n*num_wts] *= norm_factor;
	}
    }
  else if ( rv->norm_opt == NORM_OPT_NONE )
    {
    }
}

636
static
637
void set_yac_coordinates(int remap_grid_type, size_t cell_add, size_t num_cell_corners, remapgrid_t *remap_grid, struct grid_cell *yac_grid_cell)
638
639
640
{
  if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
    {
641
642
643
      size_t nx = remap_grid->dims[0];
      size_t iy = cell_add/nx;
      size_t ix = cell_add - iy*nx;
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
      const double *restrict reg2d_corner_lon = remap_grid->reg2d_corner_lon;
      const double *restrict reg2d_corner_lat = remap_grid->reg2d_corner_lat;
      double *restrict coordinates_x = yac_grid_cell->coordinates_x;
      double *restrict coordinates_y = yac_grid_cell->coordinates_y;

      coordinates_x[0] = reg2d_corner_lon[ix  ];
      coordinates_y[0] = reg2d_corner_lat[iy  ];
      coordinates_x[1] = reg2d_corner_lon[ix+1];
      coordinates_y[1] = reg2d_corner_lat[iy  ];
      coordinates_x[2] = reg2d_corner_lon[ix+1];
      coordinates_y[2] = reg2d_corner_lat[iy+1];
      coordinates_x[3] = reg2d_corner_lon[ix  ];
      coordinates_y[3] = reg2d_corner_lat[iy+1];
    }
  else
    {
      const double *restrict cell_corner_lon = remap_grid->cell_corner_lon;
      const double *restrict cell_corner_lat = remap_grid->cell_corner_lat;
      double *restrict coordinates_x = yac_grid_cell->coordinates_x;
      double *restrict coordinates_y = yac_grid_cell->coordinates_y;
664
      for ( size_t ic = 0; ic < num_cell_corners; ++ic )
665
666
667
668
669
670
671
672
673
        {
          coordinates_x[ic] = cell_corner_lon[cell_add*num_cell_corners+ic];
          coordinates_y[ic] = cell_corner_lat[cell_add*num_cell_corners+ic];
        }
    }
      
  const double *restrict coordinates_x = yac_grid_cell->coordinates_x;
  const double *restrict coordinates_y = yac_grid_cell->coordinates_y;
  double *restrict coordinates_xyz = yac_grid_cell->coordinates_xyz;
674
  for ( size_t ic = 0; ic < num_cell_corners; ++ic )
675
676
677
    LLtoXYZ(coordinates_x[ic], coordinates_y[ic], coordinates_xyz+ic*3);
}

678
679
680
static
void reg2d_bound_box(remapgrid_t *remap_grid, double *grid_bound_box)
{
681
682
  size_t nx = remap_grid->dims[0];
  size_t ny = remap_grid->dims[1];
683
684
685
686
687
688
689
690
691
692
693
694
695
696
  const double *restrict reg2d_corner_lon = remap_grid->reg2d_corner_lon;
  const double *restrict reg2d_corner_lat = remap_grid->reg2d_corner_lat;

  grid_bound_box[0] = reg2d_corner_lat[0];
  grid_bound_box[1] = reg2d_corner_lat[ny];
  if ( grid_bound_box[0] > grid_bound_box[1] )
    {
      grid_bound_box[0] = reg2d_corner_lat[ny];
      grid_bound_box[1] = reg2d_corner_lat[0];
    }
  grid_bound_box[2] = reg2d_corner_lon[0];
  grid_bound_box[3] = reg2d_corner_lon[nx];
}

697

Uwe Schulzweida's avatar
Uwe Schulzweida committed
698
void remap_conserv_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699
{
700
  bool lcheck = true;
701
  size_t srch_corners;  // num of corners of srch cells
Uwe Schulzweida's avatar
Uwe Schulzweida committed
702
703
704
705
706
707
708
709
710
711
712
713

  /* Variables necessary if segment manages to hit pole */
  int src_remap_grid_type = src_grid->remap_grid_type;
  int tgt_remap_grid_type = tgt_grid->remap_grid_type;
  extern int timer_remap_con;

  if ( cdoVerbose ) cdoPrint("Called %s()", __func__);

  progressInit();

  if ( cdoTimer ) timer_start(timer_remap_con);

714
715
  size_t src_grid_size = src_grid->size;
  size_t tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
716

717
718
  size_t src_num_cell_corners = src_grid->num_cell_corners;
  size_t tgt_num_cell_corners = tgt_grid->num_cell_corners;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
719

720
  size_t max_num_cell_corners = src_num_cell_corners;
721
722
  if ( tgt_num_cell_corners > max_num_cell_corners ) max_num_cell_corners = tgt_num_cell_corners;

723
  std::vector<enum yac_edge_type> great_circle_type(max_num_cell_corners);
724
  for ( size_t i = 0; i < max_num_cell_corners; ++i ) great_circle_type[i] = GREAT_CIRCLE;
725

726
  enum yac_edge_type lonlat_circle_type[] = {LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727

728
729
  enum yac_edge_type *src_edge_type = great_circle_type.data();
  enum yac_edge_type *tgt_edge_type = great_circle_type.data();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730

731
732
  enum cell_type target_cell_type = UNDEF_CELL;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
734
735
736
737
738
739
740
741
  if ( src_num_cell_corners == 4 )
    {
      int lonlat_circle_index = get_lonlat_circle_index(src_grid);
      if ( lonlat_circle_index >= 0 ) src_edge_type = &lonlat_circle_type[lonlat_circle_index];
    }

  if ( tgt_num_cell_corners == 4 )
    {
      int lonlat_circle_index = get_lonlat_circle_index(tgt_grid);
742
743
744
745
746
747
748
749
750
751
      if ( lonlat_circle_index >= 0 )
	{
	  target_cell_type = LON_LAT_CELL;
	  tgt_edge_type = &lonlat_circle_type[lonlat_circle_index];
	}
    }

  if ( !(tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL) )
    {
      if ( tgt_grid->cell_center_lon == NULL || tgt_grid->cell_center_lat == NULL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752
	cdoAbort("Internal problem (%s): missing target point coordinates!", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
753
754
755
    }


Uwe Schulzweida's avatar
Uwe Schulzweida committed
756
  std::vector<struct grid_cell> tgt_grid_cell(ompNumThreads);  
757
  for ( int i = 0; i < ompNumThreads; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
758
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759
760
761
762
763
764
      tgt_grid_cell[i].array_size      = tgt_num_cell_corners;
      tgt_grid_cell[i].num_corners     = tgt_num_cell_corners;
      tgt_grid_cell[i].edge_type       = tgt_edge_type;
      tgt_grid_cell[i].coordinates_x   = new double[tgt_num_cell_corners];
      tgt_grid_cell[i].coordinates_y   = new double[tgt_num_cell_corners];
      tgt_grid_cell[i].coordinates_xyz = new double[3*tgt_num_cell_corners];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
765
766
    }

767
  std::vector<search_t> search(ompNumThreads);
768
  for ( int i = 0; i < ompNumThreads; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
769
    {
770
771
772
773
774
775
776
      search[i].srch_corners    = src_num_cell_corners;
      search[i].src_edge_type   = src_edge_type;
      search[i].max_srch_cells  = 0;
      search[i].partial_areas   = NULL;
      search[i].partial_weights = NULL;
      search[i].src_grid_cells  = NULL;
      search[i].overlap_buffer  = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
777
778
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
779
780
  std::vector<size_t *> srch_add(ompNumThreads);
  for ( int i = 0; i < ompNumThreads; ++i ) srch_add[i] = new size_t[src_grid_size];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
781
782
783

  srch_corners = src_num_cell_corners;

784
  double src_grid_bound_box[4];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
785
  if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
786
    reg2d_bound_box(src_grid, src_grid_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
787

788
  weightlinks_t *weightlinks = (weightlinks_t *) Malloc(tgt_grid_size*sizeof(weightlinks_t));
789
  
790
  double findex = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
791

792
793
  size_t sum_srch_cells = 0;
  size_t sum_srch_cells2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
794

Uwe Schulzweida's avatar
Uwe Schulzweida committed
795
796
797
#ifdef STIMER
  double stimer = 0;
#endif
798
799

  // Loop over destination grid
800

801
802
#ifdef  HAVE_OPENMP4
#pragma omp parallel for schedule(dynamic) default(none)  reduction(+:findex) \
803
  shared(ompNumThreads, src_remap_grid_type, tgt_remap_grid_type, src_grid_bound_box, \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
804
	 rv, cdoVerbose, tgt_num_cell_corners, target_cell_type, \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
805
         weightlinks,  srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, \
806
	 search, srch_add, tgt_grid_cell, sum_srch_cells, sum_srch_cells2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
807
#endif
808
  for ( size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
809
    {
810
      double partial_weight;
811
      size_t src_cell_add;       // current linear address for source grid cell
812
813
      size_t num_srch_cells;
      size_t n, num_weights, num_weights_old;
814
      int ompthID = cdo_omp_get_thread_num();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
815
816

      findex++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
817
      if ( ompthID == 0 ) progressStatus(0, 1, findex/tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
818

819
820
      weightlinks[tgt_cell_add].nlinks = 0;	

821
      // Get search cells
Uwe Schulzweida's avatar
Uwe Schulzweida committed
822
823
824
825
#ifdef STIMER
      clock_t start = clock();
#endif
          
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
827
828
      if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D && tgt_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  double tgt_cell_bound_box[4];
829
	  boundbox_from_corners_reg2d(tgt_cell_add, tgt_grid->dims, tgt_grid->reg2d_corner_lon, tgt_grid->reg2d_corner_lat, tgt_cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
	  restrict_boundbox(src_grid_bound_box, tgt_cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831

Uwe Schulzweida's avatar
Uwe Schulzweida committed
832
	  num_srch_cells = get_srch_cells_reg2d(src_grid->dims, src_grid->reg2d_corner_lat, src_grid->reg2d_corner_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
						tgt_cell_bound_box, srch_add[ompthID]);
834
835
836
837

	  if ( num_srch_cells == 1 && src_grid->dims[0] == 1 && src_grid->dims[1] == 1 &&
	       IS_EQUAL(src_grid->reg2d_corner_lat[0], src_grid->reg2d_corner_lat[1]) && 
	       IS_EQUAL(src_grid->reg2d_corner_lon[0], src_grid->reg2d_corner_lon[1]) ) num_srch_cells = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
838
839
840
841
	}
      else if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  double tgt_cell_bound_box[4];
842
	  boundbox_from_corners1(tgt_cell_add, tgt_num_cell_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat, tgt_cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843
	  restrict_boundbox(src_grid_bound_box, tgt_cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
844

Uwe Schulzweida's avatar
Uwe Schulzweida committed
845
	  num_srch_cells = get_srch_cells_reg2d(src_grid->dims, src_grid->reg2d_corner_lat, src_grid->reg2d_corner_lon,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
846
						tgt_cell_bound_box, srch_add[ompthID]);
847
848
849
850

	  if ( num_srch_cells == 1 && src_grid->dims[0] == 1 && src_grid->dims[1] == 1 &&
	       IS_EQUAL(src_grid->reg2d_corner_lat[0], src_grid->reg2d_corner_lat[1]) && 
	       IS_EQUAL(src_grid->reg2d_corner_lon[0], src_grid->reg2d_corner_lon[1]) ) num_srch_cells = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
851
852
853
854
	}
      else
	{
	  restr_t tgt_cell_bound_box_r[4];
855
	  boundbox_from_corners1r(tgt_cell_add, tgt_num_cell_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat, tgt_cell_bound_box_r);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
856

857
          size_t nbins = src_grid->num_srch_bins;
858
	  num_srch_cells = get_srch_cells(tgt_cell_add, nbins, tgt_grid->bin_addr, src_grid->bin_addr,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
859
					  tgt_cell_bound_box_r, src_grid->cell_bound_box, src_grid_size, srch_add[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
860
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
861
862
863
864
#ifdef STIMER
      clock_t finish = clock();
      stimer += ((double)(finish-start))/CLOCKS_PER_SEC;
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
865

Uwe Schulzweida's avatar
Uwe Schulzweida committed
866
      if ( 0 && cdoVerbose ) sum_srch_cells += num_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
867

868
      if ( 0 && cdoVerbose )
869
	printf("tgt_cell_add %zu  num_srch_cells %zu\n", tgt_cell_add, num_srch_cells);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
870
871
872

      if ( num_srch_cells == 0 ) continue;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
873
      set_yac_coordinates(tgt_remap_grid_type, tgt_cell_add, tgt_num_cell_corners, tgt_grid, &tgt_grid_cell[ompthID]);
874

875
      // Create search arrays
Uwe Schulzweida's avatar
Uwe Schulzweida committed
876

877
      search_realloc(num_srch_cells, &search[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
878

879
880
881
      double *partial_areas   = search[ompthID].partial_areas;
      double *partial_weights = search[ompthID].partial_weights;
      struct grid_cell *src_grid_cells = search[ompthID].src_grid_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
882
883
884

      for ( n = 0; n < num_srch_cells; ++n )
	{
885
	  size_t srch_corners_new = srch_corners;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
886
	  src_cell_add = srch_add[ompthID][n];
887
          set_yac_coordinates(src_remap_grid_type, src_cell_add, srch_corners_new, src_grid, &src_grid_cells[n]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
888
889
	}

890
891
      if ( tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
892
	  cdo_compute_overlap_areas(num_srch_cells, &search[ompthID], tgt_grid_cell[ompthID]);
893
894
895
	}
      else
	{
896
897
	  double cell_center_lon = tgt_grid->cell_center_lon[tgt_cell_add];
	  double cell_center_lat = tgt_grid->cell_center_lat[tgt_cell_add];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898
	  cdo_compute_concave_overlap_areas(num_srch_cells, &search[ompthID], tgt_grid_cell[ompthID], cell_center_lon, cell_center_lat);
899
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
900

Uwe Schulzweida's avatar
Uwe Schulzweida committed
901
902
      double tgt_area = gridcell_area(tgt_grid_cell[ompthID]);
      // double tgt_area = cell_area(&tgt_grid_cell[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
903
904
905
906
907
908

      for ( num_weights = 0, n = 0; n < num_srch_cells; ++n )
	{
	  if ( partial_areas[n] > 0 )
	    {
	      partial_areas[num_weights] = partial_areas[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
909
	      srch_add[ompthID][num_weights] = srch_add[ompthID][n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
910
911
912
913
	      num_weights++;
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
914
      if ( 0 && cdoVerbose ) sum_srch_cells2 += num_weights;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
915
916
917
918

      for ( n = 0; n < num_weights; ++n )
	partial_weights[n] = partial_areas[n] / tgt_area;

919
      if ( rv->norm_opt == NORM_OPT_FRACAREA )
920
	yac_correct_weights((unsigned)num_weights, partial_weights);
921
922
923

      for ( n = 0; n < num_weights; ++n )
	partial_weights[n] *= tgt_area;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
924

925
      num_weights_old = num_weights;
926
927
      num_weights = 0;
      for ( n = 0; n < num_weights_old; ++n )
928
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
	  src_cell_add = srch_add[ompthID][n];
930
931
	  if ( partial_weights[n] <= 0. ) src_cell_add = src_grid_size;
	  if ( src_cell_add != src_grid_size )
932
933
	    {
	      partial_weights[num_weights] = partial_weights[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
934
	      srch_add[ompthID][num_weights] = src_cell_add;
935
936
937
938
	      num_weights++;
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
939
940
      for ( n = 0; n < num_weights; ++n )
	{
941
	  partial_weight = partial_weights[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
	  src_cell_add = srch_add[ompthID][n];
943

944
945
946
#if defined(_OPENMP)
#pragma omp atomic
#endif
947
	  src_grid->cell_area[src_cell_add] += partial_weight;
948
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
949
950


951
952
953
      num_weights_old = num_weights;
      for ( num_weights = 0, n = 0; n < num_weights_old; ++n )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
954
	  src_cell_add = srch_add[ompthID][n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
955
956
957
958
959
960

	  /*
	    Store the appropriate addresses and weights. 
	    Also add contributions to cell areas.
	    The source grid mask is the master mask
	  */
961
	  if ( src_grid->mask[src_cell_add] )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
962
	    {
963
	      partial_weights[num_weights] = partial_weights[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
964
	      srch_add[ompthID][num_weights] = src_cell_add;
965
966
967
968
969
970
	      num_weights++;
	    }
	}

      for ( n = 0; n < num_weights; ++n )
	{
971
	  partial_weight = partial_weights[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
972
	  src_cell_add = srch_add[ompthID][n];
973

Uwe Schulzweida's avatar
Uwe Schulzweida committed
974
975
976
#if defined(_OPENMP)
#pragma omp atomic
#endif
977
	  src_grid->cell_frac[src_cell_add] += partial_weight;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
978
		  
979
	  tgt_grid->cell_frac[tgt_cell_add] += partial_weight;
980
981
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
982
      store_weightlinks(1, num_weights, srch_add[ompthID], partial_weights, tgt_cell_add, weightlinks);
983

984
985
      tgt_grid->cell_area[tgt_cell_add] = tgt_area; 
      // printf("area %d %g %g\n", tgt_cell_add, tgt_grid->cell_area[tgt_cell_add], tgt_area);
986
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
987
988
989
990

#ifdef STIMER
printf("stime = %gs\n", stimer);
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
991

992
  if ( 0 && cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
993
    {
994
995
      printf("sum_srch_cells : %zu\n", sum_srch_cells);
      printf("sum_srch_cells2: %zu\n", sum_srch_cells2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
996
997
    }

998
  // Finished with all cells: deallocate search arrays
Uwe Schulzweida's avatar
Uwe Schulzweida committed
999

1000
  for ( int ompthID = 0; ompthID < ompNumThreads; ++ompthID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1001
    {
1002
      search_free(&search[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1003

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1004
1005
1006
      delete[] tgt_grid_cell[ompthID].coordinates_x;
      delete[] tgt_grid_cell[ompthID].coordinates_y;
      delete[] tgt_grid_cell[ompthID].coordinates_xyz;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1007

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1008
      delete[] srch_add[ompthID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1009
    }
1010
  
1011
  weightlinks2remaplinks(1, tgt_grid_size, weightlinks, rv);
1012

1013
  if ( weightlinks ) Free(weightlinks);
1014

1015
1016
  // Normalize weights using destination area if requested
  remapNormalizeWeights(tgt_grid, rv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1017

1018
  if ( cdoVerbose ) cdoPrint("Total number of links = %zu", rv->num_links);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1019
  
1020
  for ( size_t n = 0; n < src_grid_size; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1021
1022
    if ( IS_NOT_EQUAL(src_grid->cell_area[n], 0) ) src_grid->cell_frac[n] /= src_grid->cell_area[n];

1023
  for ( size_t n = 0; n < tgt_grid_size; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1024
1025
    if ( IS_NOT_EQUAL(tgt_grid->cell_area[n], 0) ) tgt_grid->cell_frac[n] /= tgt_grid->cell_area[n];

1026
  // Perform some error checking on final weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1027
1028
  if ( lcheck )
    {
1029
1030
      remapCheckArea(src_grid_size, src_grid->cell_area, "Source");
      remapCheckArea(tgt_grid_size, tgt_grid->cell_area, "Target");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1031

1032
1033
      remapCheckWeights(rv->num_links, rv->num_wts, rv->norm_opt, rv->src_cell_add, rv->tgt_cell_add, rv->wts);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1034
1035
1036

  if ( cdoTimer ) timer_stop(timer_remap_con);

1037
} // remap_weights_conserv
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1038
1039


1040
1041
//void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval)
void remap_conserv(remapgrid_t *, remapgrid_t *, const double* restrict , double* restrict , double )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1042
{
1043
} // remap_conserv