Commit 3c16d90e authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remaplib: Changed remap weights from 2D to 1D array

parent 2adea761
2011-01-07 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* remaplib: Changed remap weights from 2D to 1D array
2011-01-06 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* using CDI library version 1.4.6
......
......@@ -331,7 +331,7 @@ void *Remap(void *argument)
{
int ival;
ival = atoi(envstr);
if ( ival > 0 )
if ( ival >= 0 )
{
remap_non_global = ival;
if ( cdoVerbose )
......@@ -344,7 +344,7 @@ void *Remap(void *argument)
{
int ival;
ival = atoi(envstr);
if ( ival > 0 )
if ( ival >= 0 )
{
remap_store_link_fast = ival;
if ( cdoVerbose )
......@@ -845,10 +845,10 @@ void *Remap(void *argument)
if ( operfunc == REMAPLAF )
remap_laf(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.grid2_add, remaps[r].vars.grid1_add, array1);
remaps[r].vars.num_wts, remaps[r].vars.grid2_add, remaps[r].vars.grid1_add, array1);
else if ( operfunc == REMAPSUM )
remap_sum(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.grid2_add, remaps[r].vars.grid1_add, array1);
remaps[r].vars.num_wts, remaps[r].vars.grid2_add, remaps[r].vars.grid1_add, array1);
else
remap(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.grid2_add, remaps[r].vars.grid1_add,
......
......@@ -127,7 +127,7 @@ typedef struct {
int *grid1_add; /* grid1 address for each link */
int *grid2_add; /* grid2 address for each link */
double *wts[4]; /* map weights for each link [max_links][num_wts] */
double *wts; /* map weights for each link [max_links*num_wts] */
remaplink_t links;
}
......@@ -150,16 +150,16 @@ void remapVarsInit(int map_type, remapgrid_t *rg, remapvars_t *rv);
void remapVarsFree(remapvars_t *rv);
void remapGridFree(remapgrid_t *rg);
void remap(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
void remap(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array,
const double *restrict src_grad1, const double *restrict src_grad2, const double *restrict src_grad3,
remaplink_t links);
void remap_laf(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_laf(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_sum(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_sum(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_bilin(remapgrid_t *rg, remapvars_t *rv);
void remap_bicub(remapgrid_t *rg, remapvars_t *rv);
......@@ -176,9 +176,8 @@ void remap_gradients(remapgrid_t rg, const double *restrict array, double *restr
void reorder_links(remapvars_t *rv);
void sort_add(long num_links, long num_wts, int *restrict add1, int *restrict add2, double *restrict *restrict weights);
void sort_add_test(long num_links, long num_wts, int *restrict add1, int *restrict add2, double *restrict *restrict weights);
void sort_iter(long num_links, long num_wts, int *restrict add1, int *restrict add2, double *restrict *restrict weights, int parent);
void sort_add(long num_links, long num_wts, int *restrict add1, int *restrict add2, double *restrict weights);
void sort_iter(long num_links, long num_wts, int *restrict add1, int *restrict add2, double *restrict weights, int parent);
void write_remap_scrip(const char *interp_file, int map_type, int submap_type,
int remap_order, remapgrid_t rg, remapvars_t rv);
......
......@@ -40,10 +40,10 @@
/*
2009-05-25 Uwe Schulzweida: Changed restict data type from double to int
2009-01-11 Uwe Schulzweida: OpenMP parallelization
2011-01-07 Uwe Schulzweida: Changed remap weights from 2D to 1D array
*/
#if defined (HAVE_CONFIG_H)
# include "config.h"
#endif
......@@ -116,9 +116,6 @@ static double north_thresh = 2.00; /* threshold for coord transf. */
static double south_thresh =-2.00; /* threshold for coord transf. */
double intlin(double x, double y1, double x1, double y2, double x2);
void sort_iter(long num_links, long num_wts, int *restrict add1, int *restrict add2,
double *restrict *restrict weights, int parent);
extern int timer_remap, timer_remap_sort, timer_remap_con, timer_remap_con2, timer_remap_con3;
......@@ -172,8 +169,7 @@ void remapVarsFree(remapvars_t *rv)
free(rv->grid1_add);
free(rv->grid2_add);
for ( i = 0; i < 4; i++ )
if ( rv->wts[i] ) free(rv->wts[i]);
free(rv->wts);
if ( rv->links.option == TRUE )
{
......@@ -1277,8 +1273,6 @@ void remapGridInit(int map_type, int lextrapolate, int gridID1, int gridID2, rem
*/
void remapVarsInit(int map_type, remapgrid_t *rg, remapvars_t *rv)
{
long i;
/* Initialize all pointer */
if ( rv->pinit == FALSE )
{
......@@ -1286,7 +1280,7 @@ void remapVarsInit(int map_type, remapgrid_t *rg, remapvars_t *rv)
rv->grid1_add = NULL;
rv->grid2_add = NULL;
for ( i = 0; i < 4; i++ ) rv->wts[i] = NULL;
rv->wts = NULL;
}
/* Determine the number of weights */
......@@ -1314,8 +1308,7 @@ void remapVarsInit(int map_type, remapgrid_t *rg, remapvars_t *rv)
rv->grid1_add = (int *) realloc(rv->grid1_add, rv->max_links*sizeof(int));
rv->grid2_add = (int *) realloc(rv->grid2_add, rv->max_links*sizeof(int));
for ( i = 0; i < rv->num_wts; i++ )
rv->wts[i] = (double *) realloc(rv->wts[i], rv->max_links*sizeof(double));
rv->wts = (double *) realloc(rv->wts, rv->num_wts*rv->max_links*sizeof(double));
rv->links.option = FALSE;
rv->links.max_links = 0;
......@@ -1339,7 +1332,6 @@ void resize_remap_vars(remapvars_t *rv, int increment)
Input variables:
int increment ! the number of links to add(subtract) to arrays
*/
long i;
/* Reallocate arrays at new size */
......@@ -1350,8 +1342,7 @@ void resize_remap_vars(remapvars_t *rv, int increment)
rv->grid1_add = (int *) realloc(rv->grid1_add, rv->max_links*sizeof(int));
rv->grid2_add = (int *) realloc(rv->grid2_add, rv->max_links*sizeof(int));
for ( i = 0; i < rv->num_wts; i++ )
rv->wts[i] = (double *) realloc(rv->wts[i], rv->max_links*sizeof(double));
rv->wts = (double *) realloc(rv->wts, rv->num_wts*rv->max_links*sizeof(double));
}
} /* resize_remap_vars */
......@@ -1364,7 +1355,7 @@ void resize_remap_vars(remapvars_t *rv, int increment)
-----------------------------------------------------------------------
*/
void remap(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
void remap(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array,
const double *restrict src_grad1, const double *restrict src_grad2, const double *restrict src_grad3,
remaplink_t links)
......@@ -1377,7 +1368,7 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
int num_wts ! num of weights used in remapping
double **map_wts ! remapping weights for each link
double *map_wts ! remapping weights for each link
double *src_array ! array with source field to be remapped
......@@ -1424,7 +1415,7 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
#endif
for ( n = 0; n < links.num_links[j]; n++ )
{
dst_array[links.dst_add[j][n]] += src_array[links.src_add[j][n]]*map_wts[0][links.w_index[j][n]];
dst_array[links.dst_add[j][n]] += src_array[links.src_add[j][n]]*map_wts[num_wts*links.w_index[j][n]];
}
}
}
......@@ -1433,9 +1424,9 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
for ( n = 0; n < num_links; n++ )
{
/*
printf("%5d %5d %5d %g # dst_add src_add n\n", dst_add[n], src_add[n], n, map_wts[0][n]);
printf("%5d %5d %5d %g # dst_add src_add n\n", dst_add[n], src_add[n], n, map_wts[num_wts*n]);
*/
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[0][n];
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[num_wts*n];
}
}
}
......@@ -1445,19 +1436,19 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
{
for ( n = 0; n < num_links; n++ )
{
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[0][n] +
src_grad1[src_add[n]]*map_wts[1][n] +
src_grad2[src_add[n]]*map_wts[2][n];
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[num_wts*n] +
src_grad1[src_add[n]]*map_wts[num_wts*n+1] +
src_grad2[src_add[n]]*map_wts[num_wts*n+2];
}
}
else if ( num_wts == 4 )
{
for ( n = 0; n < num_links; n++ )
{
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[0][n] +
src_grad1[src_add[n]]*map_wts[1][n] +
src_grad2[src_add[n]]*map_wts[2][n] +
src_grad3[src_add[n]]*map_wts[3][n];
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[num_wts*n] +
src_grad1[src_add[n]]*map_wts[num_wts*n+1] +
src_grad2[src_add[n]]*map_wts[num_wts*n+2] +
src_grad3[src_add[n]]*map_wts[num_wts*n+3];
}
}
}
......@@ -1465,7 +1456,7 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
if ( cdoTimer ) timer_stop(timer_remap);
}
static
long get_max_add(long num_links, long size, const int *restrict add)
{
long n, i;
......@@ -1491,8 +1482,8 @@ long get_max_add(long num_links, long size, const int *restrict add)
-----------------------------------------------------------------------
*/
void remap_laf(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array)
void remap_laf(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array)
{
/*
Input arrays:
......@@ -1500,7 +1491,9 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
int *dst_add ! destination address for each link
int *src_add ! source address for each link
double **map_wts ! remapping weights for each link
int num_wts ! num of weights used in remapping
double *map_wts ! remapping weights for each link
double *src_array ! array with source field to be remapped
......@@ -1544,7 +1537,7 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
#if defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(dst_size, src_cls2, src_wts2, num_links, dst_add, src_add, src_array, map_wts, \
dst_array, max_cls) \
num_wts, dst_array, max_cls) \
private(i, n, k, ompthID, src_cls, src_wts, ncls, imax, wts) \
schedule(dynamic,1)
#endif
......@@ -1572,7 +1565,7 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
ncls++;
}
src_wts[k] += map_wts[0][n];
src_wts[k] += map_wts[num_wts*n];
}
}
*/
......@@ -1605,7 +1598,7 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
ncls++;
}
src_wts[k] += map_wts[0][n];
src_wts[k] += map_wts[num_wts*n];
}
}
......@@ -1649,8 +1642,8 @@ void remap_laf(double *restrict dst_array, double missval, long dst_size, long n
-----------------------------------------------------------------------
*/
void remap_sum(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict *restrict map_wts,
const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array)
void remap_sum(double *restrict dst_array, double missval, long dst_size, long num_links, double *restrict map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array)
{
/*
Input arrays:
......@@ -1658,7 +1651,9 @@ void remap_sum(double *restrict dst_array, double missval, long dst_size, long n
int *dst_add ! destination address for each link
int *src_add ! source address for each link
double **map_wts ! remapping weights for each link
int num_wts ! num of weights used in remapping
double *map_wts ! remapping weights for each link
double *src_array ! array with source field to be remapped
......@@ -1680,12 +1675,12 @@ void remap_sum(double *restrict dst_array, double missval, long dst_size, long n
for ( n = 0; n < num_links; n++ )
{
/*
printf("%5d %5d %5d %g # dst_add src_add n\n", dst_add[n], src_add[n], n, map_wts[0][n]);
printf("%5d %5d %5d %g # dst_add src_add n\n", dst_add[n], src_add[n], n, map_wts[num_wts*n]);
*/
//dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[0][n];
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[0][n];
//dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[num_wts*n];
dst_array[dst_add[n]] += src_array[src_add[n]]*map_wts[num_wts*n];
printf("%ld %d %d %g %g %g\n", n, dst_add[n], src_add[n],
src_array[src_add[n]], map_wts[0][n], dst_array[dst_add[n]]);
src_array[src_add[n]], map_wts[num_wts*n], dst_array[dst_add[n]]);
}
}
......@@ -1991,7 +1986,7 @@ void store_link_bilin(remapvars_t *rv, int dst_add, const int *restrict src_add,
{
rv->grid1_add[num_links_old+n] = src_add[n]-1;
rv->grid2_add[num_links_old+n] = dst_add;
rv->wts [0][num_links_old+n] = weights[n];
rv->wts [num_links_old+n] = weights[n];
}
} /* store_link_bilin */
......@@ -2171,10 +2166,7 @@ void remap_bilin(remapgrid_t *rg, remapvars_t *rv)
/* Renormalize weights */
sum_wgts = 0.0;
for ( n = 0; n < 4; n++ ) sum_wgts += src_lats[n];
wgts[0] = src_lats[0]/sum_wgts;
wgts[1] = src_lats[1]/sum_wgts;
wgts[2] = src_lats[2]/sum_wgts;
wgts[3] = src_lats[3]/sum_wgts;
for ( n = 0; n < 4; n++ ) wgts[n] = src_lats[n]/sum_wgts;
rg->grid2_frac[dst_add] = ONE;
......@@ -2227,7 +2219,7 @@ void store_link_bicub(remapvars_t *rv, const int dst_add, const int *restrict sr
rv->grid1_add[num_links_old+n] = src_add[n]-1;
rv->grid2_add[num_links_old+n] = dst_add;
for ( k = 0; k < 4; k++ )
rv->wts[k][num_links_old+n] = weights[k][n];
rv->wts[4*(num_links_old+n)+k] = weights[k][n];
}
} /* store_link_bicub */
......@@ -2610,7 +2602,7 @@ void store_link_nbr(remapvars_t *rv, int add1, int add2, double weights)
Input variables:
int add1 ! address on grid1
int add2 ! address on grid2
double weights[] ! array of remapping weights for this link
double weights ! remapping weight for this link
*/
long nlink;
......@@ -2627,7 +2619,7 @@ void store_link_nbr(remapvars_t *rv, int add1, int add2, double weights)
rv->grid1_add[nlink] = add1;
rv->grid2_add[nlink] = add2;
rv->wts[0][nlink] = weights;
rv->wts[nlink] = weights;
} /* store_link_nbr */
......@@ -2904,7 +2896,7 @@ void store_link_nbr1(remapvars_t *rv, int add1, int add2, double weights)
Input variables:
int add1 ! address on grid1
int add2 ! address on grid2
double weights[] ! array of remapping weights for this link
double weights ! remapping weight for this link
*/
long nlink;
......@@ -2921,7 +2913,7 @@ void store_link_nbr1(remapvars_t *rv, int add1, int add2, double weights)
rv->grid1_add[nlink] = add1;
rv->grid2_add[nlink] = add2;
rv->wts[0][nlink] = weights;
rv->wts[nlink] = weights;
} /* store_link_nbr1 */
......@@ -3519,12 +3511,11 @@ void pole_intersection(long *location, double *intrsct_lat, double *intrsct_lon,
/*
This routine finds the next intersection of a destination grid
line with the line segment given by beglon, endlon, etc.
A coincidence flag is returned if the segment is entirely
coincident with an ocean grid line. The cells in which to search
for an intersection must have already been restricted in the
calling routine.
This routine finds the next intersection of a destination grid line with
the line segment given by beglon, endlon, etc.
A coincidence flag is returned if the segment is entirely coincident with
an ocean grid line. The cells in which to search for an intersection must
have already been restricted in the calling routine.
*/
static
void intersection(long *location, double *intrsct_lat, double *intrsct_lon, int *lcoinc,
......@@ -3920,19 +3911,18 @@ void intersection(long *location, double *intrsct_lat, double *intrsct_lon, int
by the input lat/lon of the endpoints.
*/
static
void line_integral(double *weights, int num_wts, double in_phi1, double in_phi2,
void line_integral(double *weights, double in_phi1, double in_phi2,
double theta1, double theta2, double grid1_lon, double grid2_lon)
{
/*
Intent(in):
int num_wts ! Number of weights to compute
double in_phi1, in_phi2, ! Longitude endpoints for the segment
double theta1, theta2, ! Latitude endpoints for the segment
double grid1_lon, ! Reference coordinates for each
double grid2_lon ! Grid (to ensure correct 0,2pi interv.)
Intent(out):
double weights[2*num_wts] ! Line integral contribution to weights
double weights[6] ! Line integral contribution to weights
*/
/* Local variables */
......@@ -3959,8 +3949,8 @@ void line_integral(double *weights, int num_wts, double in_phi1, double in_phi2,
*/
weights[0] = dphi*(sinth1 + sinth2);
weights[1] = dphi*(costh1 + costh2 + (theta1*sinth1 + theta2*sinth2));
weights[num_wts ] = weights[0];
weights[num_wts+1] = weights[1];
weights[3] = weights[0];
weights[4] = weights[1];
/*
The third and fifth weights are for the second-order phi gradient
......@@ -4001,21 +3991,20 @@ void line_integral(double *weights, int num_wts, double in_phi1, double in_phi2,
else if ( phi2 < -PI ) phi2 += PI2;
if ( (phi2-phi1) < PI && (phi2-phi1) > -PI )
weights[num_wts+2] = dphi*(phi1*f1 + phi2*f2);
weights[5] = dphi*(phi1*f1 + phi2*f2);
else
{
if ( phi1 > ZERO ) fac = PI;
else fac = -PI;
fint = f1 + (f2-f1)*(fac-phi1)/fabs(dphi);
weights[num_wts+2] = HALF*phi1*(phi1-fac)*f1 -
HALF*phi2*(phi2+fac)*f2 +
HALF*fac*(phi1+phi2)*fint;
weights[5] = HALF*phi1*(phi1-fac)*f1 -
HALF*phi2*(phi2+fac)*f2 +
HALF*fac*(phi1+phi2)*fint;
}
} /* line_integral */
static
void grid_store_init(grid_store_t *grid_store, long gridsize)
{
......@@ -4056,6 +4045,45 @@ void grid_store_init(grid_store_t *grid_store, long gridsize)
grid_store->blksize[grid_store->nblocks-1] = grid_store->max_size%grid_store->blk_size;
}
static
void grid_store_delete(grid_store_t *grid_store)
{
grid_layer_t *grid_layer, *grid_layer_f;
long ilayer;
long i, j;
long iblk;
for ( iblk = 0; iblk < grid_store->nblocks; ++iblk )
{
j = 0;
grid_layer = grid_store->layers[iblk];
for ( ilayer = 0; ilayer < grid_store->nlayers[iblk]; ++ilayer )
{
if ( cdoVerbose )
{
for ( i = 0; i < grid_store->blksize[iblk]; ++i )
if ( grid_layer->grid2_link[i] != -1 ) j++;
}
grid_layer_f = grid_layer;
free(grid_layer->grid2_link);
grid_layer = grid_layer->next;
free(grid_layer_f);
}
if ( cdoVerbose )
{
fprintf(stderr, "block = %ld nlayers = %d allocated = %d used = %ld\n",
iblk+1, grid_store->nlayers[iblk],
grid_store->nlayers[iblk]*grid_store->blksize[iblk], j);
}
}
free(grid_store->blksize);
free(grid_store->layers);
free(grid_store->nlayers);
}
/*
This routine stores the address and weight for this link in the appropriate
address and weight arrays and resizes those arrays if necessary.
......@@ -4067,7 +4095,7 @@ void store_link_cnsrv_fast(remapvars_t *rv, long add1, long add2, double *weight
Input variables:
int add1 ! address on grid1
int add2 ! address on grid2
double weights[] ! array of remapping weights for this link
double weights[6] ! array of remapping weights for this link
*/
/* Local variables */
long nlink; /* link index */
......@@ -4107,9 +4135,9 @@ void store_link_cnsrv_fast(remapvars_t *rv, long add1, long add2, double *weight
if ( lstore_link )
{
rv->wts[0][nlink] += weights[0];
rv->wts[1][nlink] += weights[1];
rv->wts[2][nlink] += weights[2];
rv->wts[3*nlink+0] += weights[0];
rv->wts[3*nlink+1] += weights[1];
rv->wts[3*nlink+2] += weights[2];
return;
}
......@@ -4147,9 +4175,9 @@ void store_link_cnsrv_fast(remapvars_t *rv, long add1, long add2, double *weight
rv->grid1_add[nlink] = add1;
rv->grid2_add[nlink] = add2;
rv->wts[0][nlink] = weights[0];
rv->wts[1][nlink] = weights[1];
rv->wts[2][nlink] = weights[2];
rv->wts[3*nlink+0] = weights[0];
rv->wts[3*nlink+1] = weights[1];
rv->wts[3*nlink+2] = weights[2];
} /* store_link_cnsrv_fast */
......@@ -4166,7 +4194,7 @@ void store_link_cnsrv(remapvars_t *rv, long add1, long add2, double *restrict we
Input variables:
int add1 ! address on grid1
int add2 ! address on grid2
double weights[] ! array of remapping weights for this link
double weights[6] ! array of remapping weights for this link
*/
/* Local variables */
long nlink, min_link, max_link; /* link index */
......@@ -4231,9 +4259,9 @@ void store_link_cnsrv(remapvars_t *rv, long add1, long add2, double *restrict we
if ( nlink <= max_link )
{
rv->wts[0][nlink] += weights[0];
rv->wts[1][nlink] += weights[1];
rv->wts[2][nlink] += weights[2];
rv->wts[3*nlink+0] += weights[0];
rv->wts[3*nlink+1] += weights[1];
rv->wts[3*nlink+2] += weights[2];
return;
}
......@@ -4252,9 +4280,9 @@ void store_link_cnsrv(remapvars_t *rv, long add1, long add2, double *restrict we
rv->grid1_add[nlink] = add1;
rv->grid2_add[nlink] = add2;
rv->wts[0][nlink] = weights[0];
rv->wts[1][nlink] = weights[1];
rv->wts[2][nlink] = weights[2];
rv->wts[3*nlink+0] = weights[0];
rv->wts[3*nlink+1] = weights[1];
rv->wts[3*nlink+2] = weights[2];
if ( link_add1[0][add1] == -1 ) link_add1[0][add1] = (int)nlink;
if ( link_add2[0][add2] == -1 ) link_add2[0][add2] = (int)nlink;
......@@ -4291,7 +4319,7 @@ void remap_conserv(remapgrid_t *rg, remapvars_t *rv)
long grid2_add; /* current linear address for grid2 cell */
long min_add; /* addresses for restricting search of */
long max_add; /* destination grid */
long n, n2, k; /* generic counters */
long n, n2, n3, k; /* generic counters */
long corner; /* corner of cell that segment starts from */
long next_corn; /* corner of cell that segment ends on */
long nbins, num_links;
......@@ -4615,10 +4643,10 @@ void remap_conserv(remapgrid_t *rg, remapvars_t *rv)
/* Compute line integral for this subsegment. */
if ( grid2_add != -1 )
line_integral(weights, rv->num_wts, beglon, intrsct_lon, beglat, intrsct_lat,
line_integral(weights, beglon, intrsct_lon, beglat, intrsct_lat,
rg->grid1_center_lon[grid1_add], rg->grid2_center_lon[grid2_add]);
else
line_integral(weights, rv->num_wts, beglon, intrsct_lon, beglat, intrsct_lat,
line_integral(weights, beglon, intrsct_lon, beglat, intrsct_lat,
rg->grid1_center_lon[grid1_add], rg->grid1_center_lon[grid1_add]);
/* If integrating in reverse order, change sign of weights */
......@@ -4641,7 +4669,7 @@ void remap_conserv(remapgrid_t *rg, remapvars_t *rv)
else
store_link_cnsrv(rv, grid1_add, grid2_add, weights, link_add1, link_add2);
rg->grid2_frac[grid2_add] += weights[rv->num_wts];
rg->grid2_frac[grid2_add] += weights[3];
}
rg->grid1_frac[grid1_add] += weights[0];
}
......@@ -4886,10 +4914,10 @@ void remap_conserv(remapgrid_t *rg, remapvars_t *rv)
/* Compute line integral for this subsegment. */
if ( grid1_add != -1 )
line_integral(weights, rv->num_wts, beglon, intrsct_lon, beglat, intrsct_lat,
line_integral(weights, beglon, intrsct_lon, beglat, intrsct_lat,
rg->grid1_center_lon[grid1_add], rg->grid2_center_lon[grid2_add]);
else
line_integral(weights, rv->num_wts, beglon, intrsct_lon, beglat, intrsct_lat,
line_integral(weights, beglon, intrsct_lon, beglat, intrsct_lat,
rg->grid2_center_lon[grid2_add], rg->grid2_center_lon[grid2_add]);
/* If integrating in reverse order, change sign of weights */
......@@ -4917,12 +4945,12 @@ void remap_conserv(remapgrid_t *rg, remapvars_t *rv)
rg->grid1_frac[grid1_add] += weights[0];
}
rg->grid2_frac[grid2_add] += weights[rv->num_wts];
rg->grid2_frac[grid2_add] += weights[3];
}
rg->grid2_area[grid2_add] += weights[rv->num_wts ];