Commit 5d814be4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remap_distwgt_scrip.c: cleanup

parent 77398174
......@@ -28,7 +28,7 @@ case "${HOSTNAME}" in
--with-udunits2=$HOME/local/udunits-2.1.24 \
--with-proj=/opt/local \
--with-curl=/opt/local \
CC=icc CFLAGS="-g -Wall -O3 -march=native -openmp -fp-model precise" \
CC=icc CFLAGS="-g -Wall -O3 -vec-report -march=native -openmp -fp-model precise" \
LIBS="-L/opt/local/lib -lopenjpeg"
;;
# i386-apple-darwin10
......
......@@ -206,6 +206,7 @@ void print_remap_info(int operfunc, remapgrid_t *src_grid, remapgrid_t *tgt_grid
double remap_threshhold = 2;
//double remap_search_radius = acos(0.01);
int remap_test = 0;
int remap_order = 1;
int remap_non_global = FALSE;
......@@ -314,6 +315,25 @@ void get_remap_env(void)
}
remap_set_threshhold(remap_threshhold);
/*
envstr = getenv("CDO_REMAP_SERACH_RADIUS");
if ( envstr )
{
double fval;
fval = atof(envstr);
if ( fval < 0 || fval > M_PI )
{
cdoAbort("CDO_REMAP_SERACH_RADIUS=%g out of bounds (0-%g)", fval, M_PI);
}
else
{
remap_search_radius = fval;
if ( cdoVerbose )
cdoPrint("Set CDO_REMAP_SEARCH_RADIUS to %g", remap_search_radius);
}
}
cdoPrint("Set CDO_REMAP_SEARCH_RADIUS to %g", remap_search_radius);
*/
envstr = getenv("REMAP_AREA_MIN");
if ( envstr )
......
......@@ -10,11 +10,11 @@
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
static
void get_restrict_add(remapgrid_t *src_grid, double plat, const int *restrict src_bin_add, long *minadd, long *maxadd)
void get_restrict_add(remapgrid_t *src_grid, double plat, const int *restrict src_bin_add, int *minadd, int *maxadd)
{
long n, n2;
long min_add = 0, max_add = 0, nm1, np1;
long nbins;
int n, n2;
int min_add = 0, max_add = 0, nm1, np1;
int nbins;
restr_t rlat;
restr_t *bin_lats = src_grid->bin_lats;
......@@ -49,6 +49,55 @@ void get_restrict_add(remapgrid_t *src_grid, double plat, const int *restrict sr
*/
}
static
void nbr_store_distance(int nadd, double distance, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist)
{
if ( num_neighbors == 1 )
{
if ( distance < nbr_dist[0] || (distance <= nbr_dist[0] && nadd < nbr_add[0]) )
{
nbr_add[0] = nadd;
nbr_dist[0] = distance;
}
}
else
{
int n, nchk;
for ( nchk = 0; nchk < num_neighbors; ++nchk )
{
if ( distance < nbr_dist[nchk] )
{
for ( n = num_neighbors-1; n > nchk; --n )
{
nbr_add[n] = nbr_add[n-1];
nbr_dist[n] = nbr_dist[n-1];
}
nbr_add[nchk] = nadd;
nbr_dist[nchk] = distance;
break;
}
}
}
}
static
void nbr_check_distance(int num_neighbors, const int *restrict nbr_add, double *restrict nbr_dist)
{
int nchk;
double distance;
/* Uwe Schulzweida: if distance is zero, set to small number */
for ( nchk = 0; nchk < num_neighbors; ++nchk )
{
if ( nbr_add[nchk] >= 0 )
{
distance = nbr_dist[nchk];
if ( IS_EQUAL(distance, 0.) ) distance = TINY;
nbr_dist[nchk] = distance;
}
}
}
/*
This routine finds the closest num_neighbor points to a search
point and computes a distance to each of the neighbors.
......@@ -74,12 +123,12 @@ void grid_search_nbr_reg2d(int num_neighbors, remapgrid_t *src_grid, int *restri
*/
/* Local variables */
int lfound;
long n, nadd, nchk;
long nx, nxm, ny;
int n, nadd;
int nx, nxm, ny;
long ii, jj;
long i, j, ix;
int i, j, ix;
int src_add[25];
long num_add = 0;
int num_add = 0;
double distance; // Angular distance
/*
double coslat_dst = cos(plat); // cos(lat) of the search point
......@@ -139,9 +188,9 @@ void grid_search_nbr_reg2d(int num_neighbors, remapgrid_t *src_grid, int *restri
if ( lfound )
{
long ix, iy;
int ix, iy;
for ( long na = 0; na < num_add; ++na )
for ( int na = 0; na < num_add; ++na )
{
nadd = src_add[na];
......@@ -161,30 +210,11 @@ void grid_search_nbr_reg2d(int num_neighbors, remapgrid_t *src_grid, int *restri
if ( distance < -1. ) distance = -1.;
distance = acos(distance);
/* Uwe Schulzweida: if distance is zero, set to small number */
if ( IS_EQUAL(distance, 0.) ) distance = TINY;
/* Store the address and distance if this is one of the smallest four so far */
for ( nchk = 0; nchk < num_neighbors; ++nchk )
{
if ( distance < nbr_dist[nchk] )
{
for ( n = num_neighbors-1; n > nchk; --n )
{
nbr_add[n] = nbr_add[n-1];
nbr_dist[n] = nbr_dist[n-1];
}
nbr_add[nchk] = nadd;
nbr_dist[nchk] = distance;
break;
}
else if ( num_neighbors == 1 && distance <= nbr_dist[0] && nadd < nbr_add[0] )
{
nbr_add[0] = nadd;
nbr_dist[0] = distance;
}
}
nbr_store_distance(nadd, distance, num_neighbors, nbr_add, nbr_dist);
}
nbr_check_distance(num_neighbors, nbr_add, nbr_dist);
}
else if ( src_grid->lextrapolate )
{
......@@ -217,8 +247,8 @@ void grid_search_nbr(int num_neighbors, remapgrid_t *src_grid, int *restrict nbr
double plon, ! longitude of the search point
*/
/* Local variables */
long n, nadd, nchk;
long min_add, max_add;
int n, nadd;
int min_add, max_add;
double distance; /* Angular distance */
/* result changed a little on a few points with high resolution grid
double xcoslat_dst = cos(plat); // cos(lat) of the search point
......@@ -245,30 +275,17 @@ void grid_search_nbr(int num_neighbors, remapgrid_t *src_grid, int *restrict nbr
(coslon_dst*coslon[nadd] + sinlon_dst*sinlon[nadd]);
/* 2008-07-30 Uwe Schulzweida: check that distance is inside the range of -1 to 1,
otherwise the result of acos(distance) is NaN */
if ( distance > 1 ) distance = 1;
if ( distance < -1 ) distance = -1;
//if ( distance < 0.99 ) continue;
//if ( distance < 0. ) continue;
if ( distance < -1. ) distance = -1.;
if ( distance > 1. ) distance = 1.;
distance = acos(distance);
/* Uwe Schulzweida: if distance is zero, set to small number */
if ( IS_EQUAL(distance, 0) ) distance = TINY;
/* Store the address and distance if this is one of the smallest four so far */
for ( nchk = 0; nchk < num_neighbors; ++nchk )
{
if ( distance < nbr_dist[nchk] )
{
for ( n = num_neighbors-1; n > nchk; --n )
{
nbr_add[n] = nbr_add[n-1];
nbr_dist[n] = nbr_dist[n-1];
}
nbr_add[nchk] = nadd;
nbr_dist[nchk] = distance;
break;
}
}
nbr_store_distance(nadd, distance, num_neighbors, nbr_add, nbr_dist);
}
nbr_check_distance(num_neighbors, nbr_add, nbr_dist);
} /* grid_search_nbr */
/*
......
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