Commit be66f43a authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added scrip_remap_bilinear()

parent a54a9678
......@@ -162,21 +162,21 @@ int maptype2operfunc(int map_type, int submap_type, int num_neighbors, int remap
}
static
void print_remap_info(int operfunc, int map_type, remapgrid_t *src_grid, remapgrid_t *tgt_grid, int nmiss)
void print_remap_info(int operfunc, remapgrid_t *src_grid, remapgrid_t *tgt_grid, int nmiss)
{
char line[256];
char tmpstr[256];
line[0] = 0;
if ( operfunc == REMAPBIL || operfunc == GENBIL ) strcpy(line, "Bilinear");
else if ( operfunc == REMAPBIC || operfunc == GENBIC ) strcpy(line, "Bicubic");
else if ( operfunc == REMAPNN || operfunc == GENNN ) strcpy(line, "Nearest neighbor");
else if ( operfunc == REMAPDIS || operfunc == GENDIS ) strcpy(line, "Distance-weighted average");
else if ( operfunc == REMAPCON || operfunc == GENCON ) strcpy(line, "First order conservative");
else if ( operfunc == REMAPCON2 || operfunc == GENCON2 ) strcpy(line, "Second order conservative");
else if ( operfunc == REMAPLAF || operfunc == GENLAF ) strcpy(line, "Largest area fraction");
else if ( operfunc == REMAPCONS || operfunc == GENCONS ) strcpy(line, "First order conservative (sphere)");
if ( operfunc == REMAPBIL || operfunc == GENBIL ) strcpy(line, "SCRIP bilinear");
else if ( operfunc == REMAPBIC || operfunc == GENBIC ) strcpy(line, "SCRIP bicubic");
else if ( operfunc == REMAPNN || operfunc == GENNN ) strcpy(line, "SCRIP nearest neighbor");
else if ( operfunc == REMAPDIS || operfunc == GENDIS ) strcpy(line, "SCRIP distance-weighted average");
else if ( operfunc == REMAPCON || operfunc == GENCON ) strcpy(line, "SCRIP first order conservative");
else if ( operfunc == REMAPCON2 || operfunc == GENCON2 ) strcpy(line, "SCRIP second order conservative");
else if ( operfunc == REMAPLAF || operfunc == GENLAF ) strcpy(line, "SCRIP largest area fraction");
else if ( operfunc == REMAPCONS || operfunc == GENCONS ) strcpy(line, "First order conservative");
else strcpy(line, "Unknown");
strcat(line, " remapping from ");
......@@ -213,7 +213,7 @@ int remap_num_srch_bins = 180;
int lremap_num_srch_bins = FALSE;
int remap_extrapolate = FALSE;
int lextrapolate = FALSE;
int max_remaps = 0;
int max_remaps = -1;
int sort_mode = HEAP_SORT;
double remap_area_min = 0;
......@@ -445,7 +445,7 @@ void *Remap(void *argument)
int taxisID1, taxisID2;
int gridID1 = -1, gridID2;
int gridtype;
int nmiss1, nmiss2, i, j, r;
int nmiss1, nmiss2, i, j, r = -1;
int *imask = NULL;
int nremaps = 0;
int norm_opt = NORM_OPT_NONE;
......@@ -460,9 +460,10 @@ void *Remap(void *argument)
double missval;
double *array1 = NULL, *array2 = NULL;
double *grad1_lat = NULL, *grad1_lon = NULL, *grad1_latlon = NULL;
remap_t *remaps;
remap_t *remaps = NULL;
char *envstr;
char *remap_file = NULL;
int remap_weights = 1;
int lwrite_remap;
if ( cdoTimer )
......@@ -580,7 +581,7 @@ void *Remap(void *argument)
gridID1 = vlistGrid(vlistID1, index);
if ( max_remaps == 0 )
if ( max_remaps == -1 )
{
nzaxis = vlistNzaxis(vlistID1);
for ( index = 0; index < nzaxis; index++ )
......@@ -599,12 +600,15 @@ void *Remap(void *argument)
cdoPrint("Set max_remaps to %d", max_remaps);
}
remaps = malloc(max_remaps*sizeof(remap_t));
for ( r = 0; r < max_remaps; r++ )
if ( max_remaps > 0 )
{
remaps[r].gridID = -1;
remaps[r].gridsize = 0;
remaps[r].nmiss = 0;
remaps = malloc(max_remaps*sizeof(remap_t));
for ( r = 0; r < max_remaps; r++ )
{
remaps[r].gridID = -1;
remaps[r].gridsize = 0;
remaps[r].nmiss = 0;
}
}
if ( operfunc == REMAPXXX )
......@@ -923,42 +927,45 @@ void *Remap(void *argument)
remap_vars_init(map_type, remaps[r].src_grid.size, remaps[r].tgt_grid.size, &remaps[r].vars);
if ( cdoTimer ) timer_stop(timer_remap_init);
print_remap_info(operfunc, map_type, &remaps[r].src_grid, &remaps[r].tgt_grid, nmiss1);
if ( map_type == MAP_TYPE_CONSERV ) remap_conserv(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BILINEAR ) remap_bilinear(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BICUBIC ) remap_bicubic(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_DISTWGT ) remap_distwgt(num_neighbors, &remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_CONSPHERE ) remap_consphere(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
if ( remap_weights )
{
print_remap_info(operfunc, &remaps[r].src_grid, &remaps[r].tgt_grid, nmiss1);
if ( remaps[r].vars.num_links != remaps[r].vars.max_links )
resize_remap_vars(&remaps[r].vars, remaps[r].vars.num_links-remaps[r].vars.max_links);
if ( map_type == MAP_TYPE_CONSERV ) scrip_remap_weights_conserv(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BILINEAR ) scrip_remap_weights_bilinear(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_BICUBIC ) scrip_remap_weights_bicubic(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_DISTWGT ) scrip_remap_weights_distwgt(num_neighbors, &remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
else if ( map_type == MAP_TYPE_CONSPHERE ) remap_weights_conserv(&remaps[r].src_grid, &remaps[r].tgt_grid, &remaps[r].vars);
if ( remaps[r].vars.sort_add )
{
if ( cdoTimer ) timer_start(timer_remap_sort);
if ( sort_mode == MERGE_SORT )
{ /*
** use a combination of the old sort_add and a split and merge approach.
** The chunk size is determined by MERGE_SORT_LIMIT_SIZE in remaplib.c.
** OpenMP parallelism is supported
*/
sort_iter(remaps[r].vars.num_links, remaps[r].vars.num_wts,
remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add,
remaps[r].vars.wts, ompNumThreads);
}
else
{ /* use a pure heap sort without any support of parallelism */
sort_add(remaps[r].vars.num_links, remaps[r].vars.num_wts,
remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add,
remaps[r].vars.wts);
if ( remaps[r].vars.num_links != remaps[r].vars.max_links )
resize_remap_vars(&remaps[r].vars, remaps[r].vars.num_links-remaps[r].vars.max_links);
if ( remaps[r].vars.sort_add )
{
if ( cdoTimer ) timer_start(timer_remap_sort);
if ( sort_mode == MERGE_SORT )
{ /*
** use a combination of the old sort_add and a split and merge approach.
** The chunk size is determined by MERGE_SORT_LIMIT_SIZE in remaplib.c.
** OpenMP parallelism is supported
*/
sort_iter(remaps[r].vars.num_links, remaps[r].vars.num_wts,
remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add,
remaps[r].vars.wts, ompNumThreads);
}
else
{ /* use a pure heap sort without any support of parallelism */
sort_add(remaps[r].vars.num_links, remaps[r].vars.num_wts,
remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add,
remaps[r].vars.wts);
}
if ( cdoTimer ) timer_stop(timer_remap_sort);
}
if ( cdoTimer ) timer_stop(timer_remap_sort);
}
if ( lwrite_remap ) goto WRITE_REMAP;
if ( lwrite_remap ) goto WRITE_REMAP;
if ( remap_test ) reorder_links(&remaps[r].vars);
if ( remap_test ) reorder_links(&remaps[r].vars);
}
}
if ( gridInqType(gridID1) == GRID_GME )
......@@ -976,16 +983,23 @@ void *Remap(void *argument)
remap_gradients(remaps[r].src_grid, array1, grad1_lat, grad1_lon, grad1_latlon);
}
if ( operfunc == REMAPLAF )
remap_laf(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add, array1);
else if ( operfunc == REMAPSUM )
remap_sum(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add, array1);
if ( remap_weights )
{
if ( operfunc == REMAPLAF )
remap_laf(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add, array1);
else if ( operfunc == REMAPSUM )
remap_sum(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_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.tgt_grid_add, remaps[r].vars.src_grid_add,
array1, grad1_lat, grad1_lon, grad1_latlon, remaps[r].vars.links);
}
else
remap(array2, missval, gridInqSize(gridID2), remaps[r].vars.num_links, remaps[r].vars.wts,
remaps[r].vars.num_wts, remaps[r].vars.tgt_grid_add, remaps[r].vars.src_grid_add,
array1, grad1_lat, grad1_lon, grad1_latlon, remaps[r].vars.links);
{
if ( map_type == MAP_TYPE_BILINEAR ) scrip_remap_bilinear(&remaps[r].src_grid, &remaps[r].tgt_grid, array1, array2, missval);
}
gridsize2 = gridInqSize(gridID2);
......@@ -1105,15 +1119,18 @@ void *Remap(void *argument)
if ( grad1_lon ) free(grad1_lon);
if ( grad1_lat ) free(grad1_lat);
for ( r = 0; r < nremaps; r++ )
if ( max_remaps > 0 )
{
remapVarsFree(&remaps[r].vars);
remapGridFree(&remaps[r].src_grid);
remapGridFree(&remaps[r].tgt_grid);
for ( r = 0; r < nremaps; r++ )
{
remapVarsFree(&remaps[r].vars);
remapGridFree(&remaps[r].src_grid);
remapGridFree(&remaps[r].tgt_grid);
}
if ( remaps ) free(remaps);
}
if ( remaps ) free(remaps);
cdoFinish();
return (0);
......
......@@ -183,11 +183,13 @@ 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 map_wts,
long num_wts, const int *restrict dst_add, const int *restrict src_add, const double *restrict src_array);
void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_consphere(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_distwgt(int num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void scrip_remap_weights_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void scrip_remap_weights_distwgt(int num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void scrip_remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void scrip_remap_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval);
void resize_remap_vars(remapvars_t *rv, int increment);
......
......@@ -76,7 +76,7 @@ void store_link_bicub(remapvars_t *rv, int dst_add, int *restrict src_add, doubl
-----------------------------------------------------------------------
*/
void remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
void scrip_remap_weights_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* Local variables */
int lwarn = TRUE;
......
......@@ -142,7 +142,7 @@ long find_ij_weights(double plon, double plat, double* restrict src_lats, double
-----------------------------------------------------------------------
*/
void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
void scrip_remap_weights_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* Local variables */
int lwarn = TRUE;
......@@ -183,7 +183,7 @@ void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *r
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, remap_max_iter, lwarn, findex) \
private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(dynamic,1)
schedule(static)
#endif
/* grid_loop1 */
for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
......@@ -226,7 +226,7 @@ void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *r
{
double iw, jw; /* current guess for bilinear coordinate */
tgt_grid->cell_frac[dst_add] = ONE;
tgt_grid->cell_frac[dst_add] = 1.;
iter = find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw);
......@@ -284,7 +284,7 @@ void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *r
if ( src_grid->mask[src_add[n]] )
icount++;
else
src_lats[n] = ZERO;
src_lats[n] = 0.;
}
if ( icount > 0 )
......@@ -295,7 +295,7 @@ void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *r
for ( n = 0; n < 4; ++n ) sum_wgts += fabs(src_lats[n]);
for ( n = 0; n < 4; ++n ) wgts[n] = fabs(src_lats[n])/sum_wgts;
tgt_grid->cell_frac[dst_add] = ONE;
tgt_grid->cell_frac[dst_add] = 1.;
#if defined(_OPENMP)
#pragma omp critical
......@@ -306,4 +306,169 @@ void remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *r
} /* grid_loop1 */
if ( cdoTimer ) timer_stop(timer_remap_bil);
} /* remap_bilinear */
} /* scrip_remap_weights_bilinear */
void scrip_remap_bilinear(remapgrid_t* src_grid, remapgrid_t* tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval)
{
/* Local variables */
int lwarn = TRUE;
int search_result;
long tgt_grid_size;
long dst_add; /* destination addresss */
long n, icount;
long iter; /* iteration counters */
int src_add[4]; /* address for the four source points */
double src_lats[4]; /* latitudes of four bilinear corners */
double src_lons[4]; /* longitudes of four bilinear corners */
double wgts[4]; /* bilinear weights for four corners */
double plat, plon; /* lat/lon coords of destination point */
double findex = 0;
extern int timer_remap_bil;
extern long remap_max_iter;
int remap_grid_type = src_grid->remap_grid_type;
if ( cdoVerbose ) cdoPrint("Called %s()", __func__);
if ( cdoTimer ) timer_start(timer_remap_bil);
progressInit();
tgt_grid_size = tgt_grid->size;
/* Compute mappings from source to target grid */
if ( src_grid->rank != 2 )
cdoAbort("Can not do bilinear interpolation when source grid rank != 2");
/* Loop over destination grid */
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, remap_max_iter, lwarn, findex) \
private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
schedule(static)
#endif
/* grid_loop1 */
for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
{
int lprogress = 1;
#if defined(_OPENMP)
if ( omp_get_thread_num() != 0 ) lprogress = 0;
#endif
#if defined(_OPENMP)
#pragma omp atomic
#endif
findex++;
if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
tgt_array[dst_add] = missval;
if ( ! tgt_grid->mask[dst_add] ) continue;
plat = tgt_grid->cell_center_lat[dst_add];
plon = tgt_grid->cell_center_lon[dst_add];
/* Find nearest square of grid points on source grid */
if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons,
plat, plon, src_grid->dims,
src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
else
search_result = grid_search(src_grid, src_add, src_lats, src_lons,
plat, plon, src_grid->dims,
src_grid->cell_center_lat, src_grid->cell_center_lon,
src_grid->cell_bound_box, src_grid->bin_addr);
/* Check to see if points are land points */
if ( search_result > 0 )
{
for ( n = 0; n < 4; ++n )
if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
}
/* If point found, find local iw,jw coordinates for weights */
if ( search_result > 0 )
{
double iw, jw; /* current guess for bilinear coordinate */
tgt_grid->cell_frac[dst_add] = 1.;
iter = find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw);
if ( iter < remap_max_iter )
{
/* Successfully found iw,jw - compute weights */
wgts[0] = (ONE-iw) * (ONE-jw);
wgts[1] = iw * (ONE-jw);
wgts[2] = iw * jw;
wgts[3] = (ONE-iw) * jw;
tgt_array[dst_add] = 0.;
for ( n = 0; n < 4; ++n ) tgt_array[dst_add] += src_array[src_add[n]]*wgts[n];
}
else
{
if ( cdoVerbose )
{
cdoPrint("Point coords: %g %g", plat, plon);
cdoPrint("Src grid lats: %g %g %g %g", src_lats[0], src_lats[1], src_lats[2], src_lats[3]);
cdoPrint("Src grid lons: %g %g %g %g", src_lons[0], src_lons[1], src_lons[2], src_lons[3]);
cdoPrint("Src grid addresses: %d %d %d %d", src_add[0], src_add[1], src_add[2], src_add[3]);
cdoPrint("Src grid lats: %g %g %g %g",
src_grid->cell_center_lat[src_add[0]], src_grid->cell_center_lat[src_add[1]],
src_grid->cell_center_lat[src_add[2]], src_grid->cell_center_lat[src_add[3]]);
cdoPrint("Src grid lons: %g %g %g %g",
src_grid->cell_center_lon[src_add[0]], src_grid->cell_center_lon[src_add[1]],
src_grid->cell_center_lon[src_add[2]], src_grid->cell_center_lon[src_add[3]]);
cdoPrint("Current iw,jw : %g %g", iw, jw);
}
if ( cdoVerbose || lwarn )
{
lwarn = FALSE;
// cdoWarning("Iteration for iw,jw exceed max iteration count of %d!", remap_max_iter);
cdoWarning("Bilinear interpolation failed for some grid points - used a distance-weighted average instead!");
}
search_result = -1;
}
}
/*
Search for bilinear failed - use a distance-weighted average instead (this is typically near the pole)
Distance was stored in src_lats!
*/
if ( search_result < 0 )
{
icount = 0;
for ( n = 0; n < 4; ++n )
{
if ( src_grid->mask[src_add[n]] )
icount++;
else
src_lats[n] = 0.;
}
if ( icount > 0 )
{
/* Renormalize weights */
double sum_wgts = 0.0; /* sum of weights for normalization */
/* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
for ( n = 0; n < 4; ++n ) sum_wgts += fabs(src_lats[n]);
for ( n = 0; n < 4; ++n ) wgts[n] = fabs(src_lats[n])/sum_wgts;
tgt_grid->cell_frac[dst_add] = 1.;
tgt_array[dst_add] = 0.;
for ( n = 0; n < 4; ++n ) tgt_array[dst_add] += src_array[src_add[n]]*wgts[n];
}
}
} /* grid_loop1 */
if ( cdoTimer ) timer_stop(timer_remap_bil);
} /* scrip_remap_bilinear */
......@@ -335,7 +335,7 @@ int get_lonlat_circle_index(remapgrid_t *remap_grid)
}
void remap_consphere(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
void remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* local variables */
......
......@@ -997,7 +997,7 @@ void line_integral(double *weights, double in_phi1, double in_phi2,
-----------------------------------------------------------------------
*/
void remap_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
void scrip_remap_weights_conserv(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* local variables */
......
......@@ -310,7 +310,7 @@ void store_link_nbr(remapvars_t *rv, int add1, int add2, double weights)
-----------------------------------------------------------------------
*/
void remap_distwgt(int num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
void scrip_remap_weights_distwgt(int num_neighbors, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
{
/* Local variables */
......
......@@ -1114,7 +1114,7 @@ void remap(double *restrict dst_array, double missval, long dst_size, long num_l
#ifdef SX
#pragma cdir nodep
#endif
for ( n = 0; n < num_links; ++n ) dst_array[dst_add[n]] = ZERO;
for ( n = 0; n < num_links; ++n ) dst_array[dst_add[n]] = 0.;
if ( iorder == 1 ) /* First order remapping */
{
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment