remap_conserv.c 39.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

7

Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
int rect_grid_search2(long *imin, long *imax, double xmin, double xmax, long nxm, const double *restrict xm);

static
long get_srch_cells_reg2d(const int *restrict src_grid_dims, 
			  const double *restrict src_corner_lat, const double *restrict src_corner_lon,
			  const double *restrict tgt_cell_bound_box, int *srch_add)
{
  long nx = src_grid_dims[0];
  long ny = src_grid_dims[1];
  long num_srch_cells;  /* num cells in restricted search arrays   */
  int lfound;
  long nxp1, nyp1;
  double src_lon_min, src_lon_max;
  int debug = 0;

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

  src_lon_min = src_corner_lon[0];
  src_lon_max = src_corner_lon[nx];

  double bound_lon1, bound_lon2;

  num_srch_cells = 0;

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

  lfound = rect_grid_search2(&jmin, &jmax, tgt_cell_bound_box[0], tgt_cell_bound_box[1], nyp1, src_corner_lat);
  // if ( jmin > 0 ) jmin--;
  // if ( jmax < (ny-2) ) jmax++;
  bound_lon1 = tgt_cell_bound_box[2];
  bound_lon2 = tgt_cell_bound_box[3];
  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);
	  for ( jm = jmin; jm <= jmax; ++jm )
	    for ( im = imin; im <= imax; ++im )
	      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);
	      for ( jm = jmin; jm <= jmax; ++jm )
		for ( im = imin2; im <= imax2; ++im )
		  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);
	      for ( jm = jmin; jm <= jmax; ++jm )
		for ( im = imin3; im <= imax3; ++im )
		  srch_add[num_srch_cells++] = jm*nx + im;
	    }
	}
    }

  if ( debug ) printf(" -> num_srch_cells: %ld\n", num_srch_cells);

  return (num_srch_cells);
}

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
void boundbox_from_corners_reg2d(long grid_add, const int *restrict grid_dims, const double *restrict corner_lon,
				 const double *restrict corner_lat, double *restrict bound_box)
{
  long nx = grid_dims[0];
131
132
  long iy = grid_add/nx;
  long ix = grid_add - iy*nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133

134
135
  double clat1 = corner_lat[iy  ];
  double clat2 = corner_lat[iy+1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

  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
void boundbox_from_corners1(long ic, long nc, const double *restrict corner_lon,
			    const double *restrict corner_lat, double *restrict bound_box)
{
  long inc, j;
  double clon, clat;

  inc = ic*nc;

  clat = corner_lat[inc];
  clon = corner_lon[inc];

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

  for ( j = 1; j < nc; ++j )
    {
      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
void boundbox_from_corners1r(long ic, long nc, const double *restrict corner_lon,
			     const double *restrict corner_lat, restr_t *restrict bound_box)
{
  long inc, j;
  restr_t clon, clat;

  inc = ic*nc;

  clat = RESTR_SCALE(corner_lat[inc]);
  clon = RESTR_SCALE(corner_lon[inc]);

  bound_box[0] = clat;
  bound_box[1] = clat;
  bound_box[2] = clon;
  bound_box[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[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
234
  if ( RESTR_ABS(bound_box[3] - bound_box[2]) > RESTR_SCALE(PI) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
    {
      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)
#include "clipping/clipping.h"
#include "clipping/area.h"
#include "clipping/geometry.h"

257
258
259
static
double gridcell_area(struct grid_cell cell)
{
260
  return yac_huiliers_area(cell);
261
262
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
264
265
266
267
268
269
270
271
static
void cdo_compute_overlap_areas(unsigned N,
			       struct grid_cell *overlap_buffer,
			       struct grid_cell *source_cells,
			       struct grid_cell  target_cell,
			       double *partial_areas)
{
  /* Do the clipping and get the cell for the overlapping area */

272
  yac_cell_clipping(N, source_cells, target_cell, overlap_buffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
275
276
277

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

  for ( unsigned n = 0; n < N; n++ )
    {
278
      partial_areas[n] = gridcell_area(overlap_buffer[n]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
280
281
282
283
284
285
    }

#ifdef VERBOSE
  for ( unsigned n = 0; n < N; n++ )
    printf("overlap area : %lf\n", partial_areas[n]);
#endif
}
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350

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
void cdo_compute_concave_overlap_areas(unsigned N,
				       struct grid_cell *overlap_buffer,
				       struct grid_cell *source_cell,
				       struct grid_cell  target_cell,
				       double target_node_x,
				       double target_node_y,
				       double * partial_areas)
{
  /*
  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!");
  */
351
  /*
352
353
354
355
  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},
356
     .edge_type       = (enum yac_edge_type[3]) {GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE},
357
     .num_corners     = 3};
358
359
360
361
  */
  double coordinates_x[3] = {-1, -1, -1};
  double coordinates_y[3] = {-1, -1, -1};
  double coordinates_xyz[9] = {-1, -1, -1};
362
  enum yac_edge_type edge_types[3] = {GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE};
363
364
365
366
367
368
  struct grid_cell target_partial_cell =
    {.coordinates_x   = coordinates_x,
     .coordinates_y   = coordinates_y,
     .coordinates_xyz = coordinates_xyz,
     .edge_type       = edge_types,
     .num_corners     = 3};
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
410
411
412
413
414
415
416
417
418
419
420
421
422
423

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

  for ( unsigned n = 0; n < N; n++) partial_areas[n] = 0.0;

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

424
      yac_cell_clipping(N, source_cell, target_partial_cell, overlap_buffer);
425
426
427
428
429

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

      for (unsigned n = 0; n < N; n++)
	{
430
	  partial_areas[n] += gridcell_area(overlap_buffer[n]);
431
432
433
434
435
436
437
438
439
440
441
442
	  // 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
  for (unsigned n = 0; n < N; n++)
    printf("overlap area %i: %lf \n", n, partial_areas[n]);
#endif
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
443
444
445
446
447
448
449
450
451
452
453
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
//#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;
	  int gridsize = remap_grid->size;
	  int num_i = 0, num_eq0 = 0, num_eq1 = 0;
	  int iadd = gridsize/3-1;

	  if ( iadd == 0 ) iadd++;

	  for ( int i = 0; i < gridsize; i += iadd )
	    {
	      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);

493
  return (lonlat_circle_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
495
496
}


497
498
499
500
501
502
503
504

static
void normalize_weights(remapgrid_t *tgt_grid, remapvars_t *rv)
{
  /* Include centroids in weights and normalize using destination area if requested */
  long n;
  long num_wts = rv->num_wts;
  long num_links = rv->num_links;
505
  long tgt_cell_add;       /* current linear address for target grid cell   */
506
507
508
509
510
511
512
513
514
515
  double norm_factor = 0;  /* factor for normalizing wts */

  if ( rv->norm_opt == NORM_OPT_DESTAREA )
    {
#if defined(SX)
#pragma vdir nodep
#endif
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
  shared(num_wts, num_links, rv, tgt_grid)	\
516
  private(n, tgt_cell_add, norm_factor)
517
518
519
#endif
      for ( n = 0; n < num_links; ++n )
	{
520
	  tgt_cell_add = rv->tgt_cell_add[n];
521

522
523
          if ( IS_NOT_EQUAL(tgt_grid->cell_area[tgt_cell_add], 0) )
	    norm_factor = ONE/tgt_grid->cell_area[tgt_cell_add];
524
525
526
527
528
529
530
531
532
533
534
535
536
537
          else
            norm_factor = ZERO;

	  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) \
  shared(num_wts, num_links, rv, tgt_grid)	\
538
  private(n, tgt_cell_add, norm_factor)
539
540
541
#endif
      for ( n = 0; n < num_links; ++n )
	{
542
	  tgt_cell_add = rv->tgt_cell_add[n];
543

544
545
          if ( IS_NOT_EQUAL(tgt_grid->cell_frac[tgt_cell_add], 0) )
	    norm_factor = ONE/tgt_grid->cell_frac[tgt_cell_add];
546
547
548
549
550
551
552
553
554
555
556
          else
            norm_factor = ZERO;

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

557

Uwe Schulzweida's avatar
Uwe Schulzweida committed
558
void remap_conserv_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559
560
561
562
563
564
565
566
{
  /* local variables */

  int    lcheck = TRUE;

  long   ioffset;
  long   src_num_cell_corners;
  long   tgt_num_cell_corners;
567
568
  long   src_cell_add;       /* current linear address for source grid cell   */
  long   tgt_cell_add;       /* current linear address for target grid cell   */
569
  long   k;                  /* generic counters                        */
570
  long   nbins;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
572
573
574
575
  long   num_wts;
  long   max_srch_cells;     /* num cells in restricted search arrays  */
  long   num_srch_cells;     /* num cells in restricted search arrays  */
  long   srch_corners;       /* num of corners of srch cells           */
  int*   srch_add;           /* global address of cells in srch arrays */
576
  int    i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594

  /* Variables necessary if segment manages to hit pole */
  long nx = 0, ny = 0;
  int src_remap_grid_type = src_grid->remap_grid_type;
  int tgt_remap_grid_type = tgt_grid->remap_grid_type;
  double src_grid_bound_box[4];
  int lyac = FALSE;
  extern int timer_remap_con;

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

  progressInit();

  nbins = src_grid->num_srch_bins;
  num_wts = rv->num_wts;

  if ( cdoTimer ) timer_start(timer_remap_con);

595
596
  long src_grid_size = src_grid->size;
  long tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
597
598
599
600

  src_num_cell_corners = src_grid->num_cell_corners;
  tgt_num_cell_corners = tgt_grid->num_cell_corners;

601
602
603
  int max_num_cell_corners = src_num_cell_corners;
  if ( tgt_num_cell_corners > max_num_cell_corners ) max_num_cell_corners = tgt_num_cell_corners;

604
  enum yac_edge_type great_circle_type[32];
605
606
  for ( int i = 0; i < max_num_cell_corners; ++i ) great_circle_type[i] = GREAT_CIRCLE;

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

609
610
  enum yac_edge_type *src_edge_type = great_circle_type;
  enum yac_edge_type *tgt_edge_type = great_circle_type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611

612
613
  enum cell_type target_cell_type = UNDEF_CELL;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
614
615
616
617
618
619
620
621
622
  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);
623
624
625
626
627
628
629
630
631
632
633
      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 )
	cdoAbort("Internal problem (remap_weights_conserv): missing target point coordinates!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
634
635
636
637
638
639
640
641
    }

  double tgt_area;

  struct grid_cell* tgt_grid_cell;
  struct grid_cell* tgt_grid_cell2[ompNumThreads];  
  for ( i = 0; i < ompNumThreads; ++i )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
      tgt_grid_cell2[i] = (struct grid_cell*) malloc(sizeof(struct grid_cell));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643
644
645
      tgt_grid_cell2[i]->array_size      = tgt_num_cell_corners;
      tgt_grid_cell2[i]->num_corners     = tgt_num_cell_corners;
      tgt_grid_cell2[i]->edge_type       = tgt_edge_type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
646
647
648
      tgt_grid_cell2[i]->coordinates_x   = (double*) malloc(tgt_num_cell_corners*sizeof(double));
      tgt_grid_cell2[i]->coordinates_y   = (double*) malloc(tgt_num_cell_corners*sizeof(double));
      tgt_grid_cell2[i]->coordinates_xyz = (double*) malloc(3*tgt_num_cell_corners*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
    }

  struct grid_cell* src_grid_cells;
  struct grid_cell* overlap_buffer;
  struct grid_cell* src_grid_cells2[ompNumThreads];
  struct grid_cell* overlap_buffer2[ompNumThreads];
  for ( i = 0; i < ompNumThreads; ++i )
    {
      src_grid_cells2[i] = NULL;
      overlap_buffer2[i] = NULL;
    }

  double* partial_areas;
  double* partial_weights;
  double* partial_areas2[ompNumThreads];
  double* partial_weights2[ompNumThreads];
  for ( i = 0; i < ompNumThreads; ++i )
    {
      partial_areas2[i]   = NULL;
      partial_weights2[i] = NULL;
    }

  long max_srch_cells2[ompNumThreads];
  for ( i = 0; i < ompNumThreads; ++i )
    max_srch_cells2[i] = 0;

  int* srch_add2[ompNumThreads];
  for ( i = 0; i < ompNumThreads; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677
    srch_add2[i] = (int*) malloc(src_grid_size*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697

  srch_corners = src_num_cell_corners;

  if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
    {
      nx = src_grid->dims[0];
      ny = src_grid->dims[1];
     
      src_grid_bound_box[0] = src_grid->reg2d_corner_lat[0];
      src_grid_bound_box[1] = src_grid->reg2d_corner_lat[ny];
      if ( src_grid_bound_box[0] > src_grid_bound_box[1] )
	{
	  src_grid_bound_box[0] = src_grid->reg2d_corner_lat[ny];
	  src_grid_bound_box[1] = src_grid->reg2d_corner_lat[0];
	}
      src_grid_bound_box[2] = src_grid->reg2d_corner_lon[0];
      src_grid_bound_box[3] = src_grid->reg2d_corner_lon[nx];
      //printf("src_grid   lon: %g %g lat: %g %g\n", RAD2DEG*src_grid_bound_box[2],RAD2DEG*src_grid_bound_box[3],RAD2DEG*src_grid_bound_box[0],RAD2DEG*src_grid_bound_box[1] );
    }

698
  weightlinks_t *weightlinks = (weightlinks_t *) malloc(tgt_grid_size*sizeof(weightlinks_t));
699
  
700
  double findex = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
701
702
703
704

  int sum_srch_cells = 0;
  int sum_srch_cells2 = 0;

705
706
  /* Loop over destination grid */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
707
#if defined(_OPENMP)
708
#pragma omp parallel for default(shared) \
709
  shared(ompNumThreads, lyac, nbins, num_wts, src_remap_grid_type, tgt_remap_grid_type, src_grid_bound_box,	\
Uwe Schulzweida's avatar
Uwe Schulzweida committed
710
	 src_edge_type, tgt_edge_type, partial_areas2, partial_weights2,  \
711
         rv, cdoVerbose, max_srch_cells2, tgt_num_cell_corners, target_cell_type, \
712
         weightlinks, \
713
         srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size,	\
Uwe Schulzweida's avatar
Uwe Schulzweida committed
714
	 overlap_buffer2, src_grid_cells2, srch_add2, tgt_grid_cell2, findex, sum_srch_cells, sum_srch_cells2) \
715
  private(srch_add, tgt_grid_cell, tgt_area, k, num_srch_cells, max_srch_cells,  \
716
	  partial_areas, partial_weights, overlap_buffer, src_grid_cells, src_cell_add, tgt_cell_add, ioffset)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
717
#endif
718
  for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
719
    {
720
721
      double partial_weight;
      long n, num_weights, num_weights_old;
722
      int ompthID = cdo_omp_get_thread_num();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
723
724
725
726
      int lprogress = 1;
      if ( ompthID != 0 ) lprogress = 0;

#if defined(_OPENMP)
727
#include "pragma_omp_atomic_update.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
729
730
731
#endif
      findex++;
      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);

732
733
      weightlinks[tgt_cell_add].nlinks = 0;	

Uwe Schulzweida's avatar
Uwe Schulzweida committed
734
735
736
737
738
739
740
741
      srch_add = srch_add2[ompthID];
      tgt_grid_cell = tgt_grid_cell2[ompthID];

      /* Get search cells */

      if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D && tgt_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  double tgt_cell_bound_box[4];
742
	  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
743
	  restrict_boundbox(src_grid_bound_box, tgt_cell_bound_box);
744
	  if ( 0 && cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
	    printf("bound_box %ld  lon: %g %g lat: %g %g\n",
746
		   tgt_cell_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747
748
	  num_srch_cells = get_srch_cells_reg2d(src_grid->dims, src_grid->reg2d_corner_lat, src_grid->reg2d_corner_lon,
						tgt_cell_bound_box, srch_add);
749
750
751
752

	  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
753
754
755
756
	}
      else if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  double tgt_cell_bound_box[4];
757
	  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
758
	  restrict_boundbox(src_grid_bound_box, tgt_cell_bound_box);
759
	  if ( 0 && cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
760
	    printf("bound_box %ld  lon: %g %g lat: %g %g\n",
761
		   tgt_cell_add, RAD2DEG*tgt_cell_bound_box[2],RAD2DEG*tgt_cell_bound_box[3],RAD2DEG*tgt_cell_bound_box[0],RAD2DEG*tgt_cell_bound_box[1] );
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
763
	  num_srch_cells = get_srch_cells_reg2d(src_grid->dims, src_grid->reg2d_corner_lat, src_grid->reg2d_corner_lon,
						tgt_cell_bound_box, srch_add);
764
765
766
767

	  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
768
769
770
771
	}
      else
	{
	  restr_t tgt_cell_bound_box_r[4];
772
	  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
773

774
	  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
775
776
777
					  tgt_cell_bound_box_r, src_grid->cell_bound_box, src_grid_size, srch_add);
	}

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

780
      if ( 0 && cdoVerbose )
781
	printf("tgt_cell_add %ld  num_srch_cells %ld\n", tgt_cell_add, num_srch_cells);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
782
783
784
785
786
787

      if ( num_srch_cells == 0 ) continue;

      if ( tgt_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	{
	  long nx = tgt_grid->dims[0];
788
789
	  long iy = tgt_cell_add/nx;
	  long ix = tgt_cell_add - iy*nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
790
791
792
793
794
795
796
797
798
799
800
801
802
803

	  tgt_grid_cell->coordinates_x[0] = tgt_grid->reg2d_corner_lon[ix  ];
	  tgt_grid_cell->coordinates_y[0] = tgt_grid->reg2d_corner_lat[iy  ];
	  tgt_grid_cell->coordinates_x[1] = tgt_grid->reg2d_corner_lon[ix+1];
	  tgt_grid_cell->coordinates_y[1] = tgt_grid->reg2d_corner_lat[iy  ];
	  tgt_grid_cell->coordinates_x[2] = tgt_grid->reg2d_corner_lon[ix+1];
	  tgt_grid_cell->coordinates_y[2] = tgt_grid->reg2d_corner_lat[iy+1];
	  tgt_grid_cell->coordinates_x[3] = tgt_grid->reg2d_corner_lon[ix  ];
	  tgt_grid_cell->coordinates_y[3] = tgt_grid->reg2d_corner_lat[iy+1];
	}
      else
	{
	  for ( int ic = 0; ic < tgt_num_cell_corners; ++ic )
	    {
804
805
	      tgt_grid_cell->coordinates_x[ic] = tgt_grid->cell_corner_lon[tgt_cell_add*tgt_num_cell_corners+ic];
	      tgt_grid_cell->coordinates_y[ic] = tgt_grid->cell_corner_lat[tgt_cell_add*tgt_num_cell_corners+ic];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
806
807
808
809
810
811
	    }
	}
      
      for ( int ic = 0; ic < tgt_num_cell_corners; ++ic )
	LLtoXYZ(tgt_grid_cell->coordinates_x[ic], tgt_grid_cell->coordinates_y[ic], tgt_grid_cell->coordinates_xyz+ic*3);

812
      //printf("target: %ld\n", tgt_cell_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
813
      if ( lyac )
814
        if ( tgt_cell_add == 174752 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
	  {
	    for ( int n = 0; n < tgt_num_cell_corners; ++n )
	      {
		printf("  TargetCell.coordinates_x[%d] = %g*rad;\n", n, tgt_grid_cell->coordinates_x[n]/DEG2RAD);
		printf("  TargetCell.coordinates_y[%d] = %g*rad;\n", n, tgt_grid_cell->coordinates_y[n]/DEG2RAD);
	      }
	    /*
	    printf("> -Z1\n");
	    for ( int n = 0; n < tgt_num_cell_corners; ++n )
		printf("  %g %g\n", tgt_grid_cell->coordinates_x[n]/DEG2RAD, tgt_grid_cell->coordinates_y[n]/DEG2RAD);
	      printf("  %g %g\n", tgt_grid_cell->coordinates_x[0]/DEG2RAD, tgt_grid_cell->coordinates_y[0]/DEG2RAD);
	    */
	  }
      
      /* Create search arrays */

      max_srch_cells  = max_srch_cells2[ompthID];
      partial_areas   = partial_areas2[ompthID];
      partial_weights = partial_weights2[ompthID];
      overlap_buffer  = overlap_buffer2[ompthID];
      src_grid_cells  = src_grid_cells2[ompthID];

      if ( num_srch_cells > max_srch_cells )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
839
840
	  partial_areas   = (double*) realloc(partial_areas,   num_srch_cells*sizeof(double));
	  partial_weights = (double*) realloc(partial_weights, num_srch_cells*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
841

Uwe Schulzweida's avatar
Uwe Schulzweida committed
842
843
	  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));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859

	  for ( n = max_srch_cells; n < num_srch_cells; ++n )
	    {
	      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;
	    }

	  for ( n = max_srch_cells; n < num_srch_cells; ++n )
	    {
	      src_grid_cells[n].array_size      = srch_corners;
	      src_grid_cells[n].num_corners     = srch_corners;
	      src_grid_cells[n].edge_type       = src_edge_type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
860
861
862
	      src_grid_cells[n].coordinates_x   = (double*) malloc(srch_corners*sizeof(double));
	      src_grid_cells[n].coordinates_y   = (double*) malloc(srch_corners*sizeof(double));
	      src_grid_cells[n].coordinates_xyz = (double*) malloc(3*srch_corners*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
863
864
865
866
867
868
869
870
871
872
873
874
875
876
	    }

	  max_srch_cells = num_srch_cells;

	  max_srch_cells2[ompthID]  = max_srch_cells;
	  partial_areas2[ompthID]   = partial_areas;
	  partial_weights2[ompthID] = partial_weights;
	  overlap_buffer2[ompthID]  = overlap_buffer;
	  src_grid_cells2[ompthID]  = src_grid_cells;
	}

      // printf("  int ii = 0;\n");
      for ( n = 0; n < num_srch_cells; ++n )
	{
877
878
	  long srch_corners_new = srch_corners;

879
	  src_cell_add = srch_add[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
880
881
882
883
884

	  if ( src_remap_grid_type == REMAP_GRID_TYPE_REG2D )
	    {
	      int ix, iy;

885
886
	      iy = src_cell_add/nx;
	      ix = src_cell_add - iy*nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904

	      src_grid_cells[n].coordinates_x[0] = src_grid->reg2d_corner_lon[ix  ];
	      src_grid_cells[n].coordinates_y[0] = src_grid->reg2d_corner_lat[iy  ];
	      src_grid_cells[n].coordinates_x[1] = src_grid->reg2d_corner_lon[ix+1];
	      src_grid_cells[n].coordinates_y[1] = src_grid->reg2d_corner_lat[iy  ];
	      src_grid_cells[n].coordinates_x[2] = src_grid->reg2d_corner_lon[ix+1];
	      src_grid_cells[n].coordinates_y[2] = src_grid->reg2d_corner_lat[iy+1];
	      src_grid_cells[n].coordinates_x[3] = src_grid->reg2d_corner_lon[ix  ];
	      src_grid_cells[n].coordinates_y[3] = src_grid->reg2d_corner_lat[iy+1];
	      /*
	      printf("source1: %ld %ld", num_srch_cells, n);
	      for ( k = 0; k < srch_corners; ++k )
		printf(" %g %g", src_grid_cells[n].coordinates_x[k]/DEG2RAD, src_grid_cells[n].coordinates_y[k]/DEG2RAD);
	      printf("\n");
	      */
	    }
	  else
	    {
905
	      ioffset = src_cell_add*srch_corners;
906
907
908
909
910
911
912
	      /*
	      for ( k = srch_corners-1; k > 0; --k )
		{
		  if ( IS_NOT_EQUAL(src_grid->cell_corner_lon[ioffset+k], src_grid->cell_corner_lon[ioffset+k-1]) ||
		       IS_NOT_EQUAL(src_grid->cell_corner_lat[ioffset+k], src_grid->cell_corner_lat[ioffset+k-1]) )
		    break;
		}
913
	      if ( k != srch_corners-1 ) printf("%ld %ld %ld %ld\n", tgt_cell_add, n, srch_corners, k+1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
914

915
916
917
918
919
920
921
	      if ( k != srch_corners-1 )
		{
		  srch_corners_new = k+1;
		  src_grid_cells[n].num_corners = srch_corners_new;
		}
	      */
	      for ( k = 0; k < srch_corners_new; ++k )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
922
923
924
925
926
		{
		  src_grid_cells[n].coordinates_x[k] = src_grid->cell_corner_lon[ioffset+k];
		  src_grid_cells[n].coordinates_y[k] = src_grid->cell_corner_lat[ioffset+k];
		}
	      /*
927
928
929
930
931
932
933
	      for ( k = 0; k < srch_corners_new; ++k )
		{
		  printf("  SourceCell[ii].coordinates_x[%ld] = %g*rad;\n", k, src_grid_cells[n].coordinates_x[k]/DEG2RAD);
		  printf("  SourceCell[ii].coordinates_y[%ld] = %g*rad;\n", k, src_grid_cells[n].coordinates_y[k]/DEG2RAD);
		}
	      */
	      /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
934
	      printf("source2: %ld %ld", num_srch_cells, n);
935
	      for ( k = 0; k < srch_corners_new; ++k )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
936
937
938
939
940
		printf(" %g %g", src_grid_cells[n].coordinates_x[k]/DEG2RAD, src_grid_cells[n].coordinates_y[k]/DEG2RAD);
	      printf("\n");
	      */
	    }

941
	  for ( int ic = 0; ic < srch_corners_new; ++ic )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
942
943
944
	    LLtoXYZ(src_grid_cells[n].coordinates_x[ic], src_grid_cells[n].coordinates_y[ic], src_grid_cells[n].coordinates_xyz+ic*3);

	  if ( lyac )
945
	    if ( tgt_cell_add == 174752 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
946
947
	    {
	      // printf("n %d\n", (int)n);
948
	      for ( k = 0; k < srch_corners_new; ++k )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
949
950
951
952
953
954
955
		{
		  printf("  SourceCell[ii].coordinates_x[%ld] = %g*rad;\n", k, src_grid_cells[n].coordinates_x[k]/DEG2RAD);
		  printf("  SourceCell[ii].coordinates_y[%ld] = %g*rad;\n", k, src_grid_cells[n].coordinates_y[k]/DEG2RAD);
		}
	      printf("  ii++;\n");
	      /*
	      printf("> -Z1\n");
956
	      for ( k = 0; k < srch_corners_new; ++k )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
957
958
959
960
961
962
		printf("  %g %g\n", src_grid_cells[n].coordinates_x[k]/DEG2RAD, src_grid_cells[n].coordinates_y[k]/DEG2RAD);
	      printf("  %g %g\n", src_grid_cells[n].coordinates_x[0]/DEG2RAD, src_grid_cells[n].coordinates_y[0]/DEG2RAD);
	      */
	    }
	}

963
964
965
966
967
968
      if ( tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL )
	{
	  cdo_compute_overlap_areas(num_srch_cells, overlap_buffer, src_grid_cells, *tgt_grid_cell, partial_areas);
	}
      else
	{
969
970
	  double cell_center_lon = tgt_grid->cell_center_lon[tgt_cell_add];
	  double cell_center_lat = tgt_grid->cell_center_lat[tgt_cell_add];
971
972
	  cdo_compute_concave_overlap_areas(num_srch_cells, overlap_buffer, src_grid_cells, *tgt_grid_cell, cell_center_lon, cell_center_lat, partial_areas);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
973

974
      tgt_area = gridcell_area(*tgt_grid_cell);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
975
976
977
978
979
980
      // tgt_area = cell_area(tgt_grid_cell);

      for ( num_weights = 0, n = 0; n < num_srch_cells; ++n )
	{
	  if ( partial_areas[n] > 0 )
	    {
981
	      //printf(">>>>   %d %d %g %g\n", (int)tgt_cell_add, srch_add[n], tgt_area, partial_areas[n]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
982
983
984
985
986
987
	      partial_areas[num_weights] = partial_areas[n];
	      srch_add[num_weights] = srch_add[n];
	      num_weights++;
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
988
      if ( 0 && cdoVerbose ) sum_srch_cells2 += num_weights;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
989
990
991
992

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

993
      if ( rv->norm_opt == NORM_OPT_FRACAREA )
994
	yac_correct_weights((unsigned)num_weights, partial_weights);
995
996
997

      for ( n = 0; n < num_weights; ++n )
	partial_weights[n] *= tgt_area;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
998
999
      //#endif

1000
      num_weights_old = num_weights;
For faster browsing, not all history is shown. View entire blame