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

added remap_search_latbins.c

parent e848eeb0
......@@ -587,6 +587,7 @@ src/readline.c -text
src/realtime.c -text
src/remap.h -text
src/remap_scrip_io.c -text
src/remap_search_latbins.c -text
src/remaplib.c -text
src/remapsort.c -text
src/results_template_parser.c -text
......
......@@ -36,6 +36,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
#include "statistic.h"
enum T_EIGEN_MODE {JACOBI, DANIELSON_LANCZOS};
......
......@@ -36,6 +36,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
#include "statistic.h"
......
......@@ -26,6 +26,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
// NO MISSING VALUE SUPPORT ADDED SO FAR
......
......@@ -26,6 +26,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
// NO MISSING VALUE SUPPORT ADDED SO FAR
......
......@@ -25,6 +25,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
void *Fldrms(void *argument)
......
......@@ -35,6 +35,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
static
void print_location_LL(int operfunc, int vlistID, int varID, int levelID, int gridID, double sglval, double *fieldptr,
......
......@@ -26,6 +26,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
/* routine corr copied from PINGO */
......
......@@ -25,6 +25,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
#include "gradsdeslib.h"
......
......@@ -270,6 +270,7 @@ cdo_SOURCES += Adisit.c \
remaplib.c \
remapsort.c \
remap_scrip_io.c \
remap_search_latbins.c \
stdnametable.c \
stdnametable.h \
specspace.c \
......
......@@ -134,17 +134,17 @@ am__cdo_SOURCES_DIST = cdo.c Adisit.c Arith.c Arithc.c Arithdays.c \
percentiles.h pipe.c pipe.h printinfo.h process.c process.h \
pstream.c pstream.h pstream_int.h pthread_debug.c \
pthread_debug.h readline.c realtime.c remap.h remaplib.c \
remapsort.c remap_scrip_io.c stdnametable.c stdnametable.h \
specspace.c specspace.h statistic.c statistic.h table.c \
timebase.h timer.c userlog.c util.c util.h vinterp.c vinterp.h \
zaxis.c clipping/clipping.c clipping/clipping.h \
clipping/area.c clipping/area.h clipping/ensure_array_size.c \
clipping/ensure_array_size.h clipping/geometry_tools.c \
clipping/geometry.h clipping/grid.h clipping/points.h \
clipping/dep_list.h clipping/grid_cell.c clipping/grid_cell.h \
clipping/intersection.c clipping/utils.c clipping/utils.h \
Magplot.c Magvector.c Maggraph.c template_parser.h \
template_parser.c results_template_parser.h \
remapsort.c remap_scrip_io.c remap_search_latbins.c \
stdnametable.c stdnametable.h specspace.c specspace.h \
statistic.c statistic.h table.c timebase.h timer.c userlog.c \
util.c util.h vinterp.c vinterp.h zaxis.c clipping/clipping.c \
clipping/clipping.h clipping/area.c clipping/area.h \
clipping/ensure_array_size.c clipping/ensure_array_size.h \
clipping/geometry_tools.c clipping/geometry.h clipping/grid.h \
clipping/points.h clipping/dep_list.h clipping/grid_cell.c \
clipping/grid_cell.h clipping/intersection.c clipping/utils.c \
clipping/utils.h Magplot.c Magvector.c Maggraph.c \
template_parser.h template_parser.c results_template_parser.h \
results_template_parser.c magics_template_parser.h \
magics_template_parser.c StringUtilities.h StringUtilities.c \
CdoMagicsMapper.h CdoMagicsMapper.c
......@@ -263,10 +263,10 @@ am_cdo_OBJECTS = cdo-cdo.$(OBJEXT) cdo-Adisit.$(OBJEXT) \
cdo-pthread_debug.$(OBJEXT) cdo-readline.$(OBJEXT) \
cdo-realtime.$(OBJEXT) cdo-remaplib.$(OBJEXT) \
cdo-remapsort.$(OBJEXT) cdo-remap_scrip_io.$(OBJEXT) \
cdo-stdnametable.$(OBJEXT) cdo-specspace.$(OBJEXT) \
cdo-statistic.$(OBJEXT) cdo-table.$(OBJEXT) \
cdo-timer.$(OBJEXT) cdo-userlog.$(OBJEXT) cdo-util.$(OBJEXT) \
cdo-vinterp.$(OBJEXT) cdo-zaxis.$(OBJEXT) \
cdo-remap_search_latbins.$(OBJEXT) cdo-stdnametable.$(OBJEXT) \
cdo-specspace.$(OBJEXT) cdo-statistic.$(OBJEXT) \
cdo-table.$(OBJEXT) cdo-timer.$(OBJEXT) cdo-userlog.$(OBJEXT) \
cdo-util.$(OBJEXT) cdo-vinterp.$(OBJEXT) cdo-zaxis.$(OBJEXT) \
cdo-clipping.$(OBJEXT) cdo-area.$(OBJEXT) \
cdo-ensure_array_size.$(OBJEXT) cdo-geometry_tools.$(OBJEXT) \
cdo-grid_cell.$(OBJEXT) cdo-intersection.$(OBJEXT) \
......@@ -517,11 +517,11 @@ cdo_SOURCES = cdo.c Adisit.c Arith.c Arithc.c Arithdays.c Arithlat.c \
pipe.c pipe.h printinfo.h process.c process.h pstream.c \
pstream.h pstream_int.h pthread_debug.c pthread_debug.h \
readline.c realtime.c remap.h remaplib.c remapsort.c \
remap_scrip_io.c stdnametable.c stdnametable.h specspace.c \
specspace.h statistic.c statistic.h table.c timebase.h timer.c \
userlog.c util.c util.h vinterp.c vinterp.h zaxis.c \
clipping/clipping.c clipping/clipping.h clipping/area.c \
clipping/area.h clipping/ensure_array_size.c \
remap_scrip_io.c remap_search_latbins.c stdnametable.c \
stdnametable.h specspace.c specspace.h statistic.c statistic.h \
table.c timebase.h timer.c userlog.c util.c util.h vinterp.c \
vinterp.h zaxis.c clipping/clipping.c clipping/clipping.h \
clipping/area.c clipping/area.h clipping/ensure_array_size.c \
clipping/ensure_array_size.h clipping/geometry_tools.c \
clipping/geometry.h clipping/grid.h clipping/points.h \
clipping/dep_list.h clipping/grid_cell.c clipping/grid_cell.h \
......@@ -896,6 +896,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-readline.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-realtime.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_scrip_io.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remap_search_latbins.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remaplib.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-remapsort.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cdo-results_template_parser.Po@am__quote@
......@@ -4153,6 +4154,20 @@ cdo-remap_scrip_io.obj: remap_scrip_io.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_scrip_io.obj `if test -f 'remap_scrip_io.c'; then $(CYGPATH_W) 'remap_scrip_io.c'; else $(CYGPATH_W) '$(srcdir)/remap_scrip_io.c'; fi`
cdo-remap_search_latbins.o: remap_search_latbins.c
@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_search_latbins.o -MD -MP -MF $(DEPDIR)/cdo-remap_search_latbins.Tpo -c -o cdo-remap_search_latbins.o `test -f 'remap_search_latbins.c' || echo '$(srcdir)/'`remap_search_latbins.c
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/cdo-remap_search_latbins.Tpo $(DEPDIR)/cdo-remap_search_latbins.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='remap_search_latbins.c' object='cdo-remap_search_latbins.o' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_search_latbins.o `test -f 'remap_search_latbins.c' || echo '$(srcdir)/'`remap_search_latbins.c
cdo-remap_search_latbins.obj: remap_search_latbins.c
@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-remap_search_latbins.obj -MD -MP -MF $(DEPDIR)/cdo-remap_search_latbins.Tpo -c -o cdo-remap_search_latbins.obj `if test -f 'remap_search_latbins.c'; then $(CYGPATH_W) 'remap_search_latbins.c'; else $(CYGPATH_W) '$(srcdir)/remap_search_latbins.c'; fi`
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/cdo-remap_search_latbins.Tpo $(DEPDIR)/cdo-remap_search_latbins.Po
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='remap_search_latbins.c' object='cdo-remap_search_latbins.obj' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o cdo-remap_search_latbins.obj `if test -f 'remap_search_latbins.c'; then $(CYGPATH_W) 'remap_search_latbins.c'; else $(CYGPATH_W) '$(srcdir)/remap_search_latbins.c'; fi`
cdo-stdnametable.o: stdnametable.c
@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(cdo_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT cdo-stdnametable.o -MD -MP -MF $(DEPDIR)/cdo-stdnametable.Tpo -c -o cdo-stdnametable.o `test -f 'stdnametable.c' || echo '$(srcdir)/'`stdnametable.c
@am__fastdepCC_TRUE@ $(am__mv) $(DEPDIR)/cdo-stdnametable.Tpo $(DEPDIR)/cdo-stdnametable.Po
......
......@@ -770,8 +770,12 @@ void *Outputgmt(void *argument)
{
if ( lgrid_gen_bounds )
{
if ( ! lzon ) genXbounds(nlon, nlat, grid_center_lon, grid_corner_lon, 0);
if ( ! lmer ) genYbounds(nlon, nlat, grid_center_lat, grid_corner_lat);
char xunitstr[CDI_MAX_NAME];
char yunitstr[CDI_MAX_NAME];
gridInqXunits(gridID, xunitstr);
gridInqYunits(gridID, yunitstr);
if ( ! lzon ) grid_cell_center_to_bounds_X2D(xunitstr, nlon, nlat, grid_center_lon, grid_corner_lon, 0);
if ( ! lmer ) grid_cell_center_to_bounds_Y2D(yunitstr, nlon, nlat, grid_center_lat, grid_corner_lat);
}
else
cdoAbort("Grid corner missing!");
......
......@@ -25,6 +25,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
void trms(field_t field1, field_t field2, double *dp, field_t *field3)
......
......@@ -25,6 +25,7 @@
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "grid.h"
void *Varrms(void *argument)
......
......@@ -117,18 +117,9 @@ int vlistIsSzipped(int vlistID);
void cdoGenFileSuffix(char *filesuffix, size_t maxlen, int filetype, int vlistID, const char *refname);
int gridWeights(int gridID, double *weights);
int gridGenArea(int gridID, double *area);
void gaussaw(double pa[], double pw[], int nlat);
void genXbounds(long xsize, long ysize, const double * restrict grid_center_lon,
double * restrict grid_corner_lon, double dlon);
void genYbounds(long xsize, long ysize, const double * restrict grid_center_lat,
double * restrict grid_corner_lat);
void writeNCgrid(const char *gridfile, int gridID, int *imask);
void defineZaxis(const char *zaxisarg);
void cdiDefTableID(int tableID);
void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *xvals);
void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *yvals);
int gridFromName(const char *gridname);
int zaxisFromName(const char *zaxisname);
......
......@@ -209,6 +209,127 @@ void grid_check_lat_borders(int n, double *ybounds)
}
void grid_cell_center_to_bounds_X2D(const char* xunitstr, long xsize, long ysize, const double* restrict grid_center_lon,
double* restrict grid_corner_lon, double dlon)
{
long i, j, index;
double minlon, maxlon;
if ( ! (dlon > 0) ) dlon = 360./xsize;
/*
if ( xsize == 1 || (grid_center_lon[xsize-1]-grid_center_lon[0]+dlon) < 359 )
cdoAbort("Cannot calculate Xbounds for %d vals with dlon = %g", xsize, dlon);
*/
for ( i = 0; i < xsize; ++i )
{
minlon = grid_center_lon[i] - 0.5*dlon;
maxlon = grid_center_lon[i] + 0.5*dlon;
for ( j = 0; j < ysize; ++j )
{
index = (j<<2)*xsize + (i<<2);
grid_corner_lon[index ] = minlon;
grid_corner_lon[index+1] = maxlon;
grid_corner_lon[index+2] = maxlon;
grid_corner_lon[index+3] = minlon;
}
}
}
static
double genYmin(double y1, double y2)
{
double ymin, dy;
dy = y2 - y1;
ymin = y1 - dy/2;
if ( y1 < -85 && ymin < -87.5 ) ymin = -90;
if ( cdoVerbose )
cdoPrint("genYmin: y1 = %g y2 = %g dy = %g ymin = %g", y1, y2, dy, ymin);
return (ymin);
}
static
double genYmax(double y1, double y2)
{
double ymax, dy;
dy = y1 - y2;
ymax = y1 + dy/2;
if ( y1 > 85 && ymax > 87.5 ) ymax = 90;
if ( cdoVerbose )
cdoPrint("genYmax: y1 = %g y2 = %g dy = %g ymax = %g", y1, y2, dy, ymax);
return (ymax);
}
/*****************************************************************************/
void grid_cell_center_to_bounds_Y2D(const char* yunitstr, long xsize, long ysize, const double* restrict grid_center_lat, double* restrict grid_corner_lat)
{
long i, j, index;
double minlat, maxlat;
double firstlat, lastlat;
firstlat = grid_center_lat[0];
lastlat = grid_center_lat[xsize*ysize-1];
// if ( ysize == 1 ) cdoAbort("Cannot calculate Ybounds for 1 value!");
for ( j = 0; j < ysize; ++j )
{
if ( ysize == 1 )
{
minlat = grid_center_lat[0] - 360./ysize;
maxlat = grid_center_lat[0] + 360./ysize;
}
else
{
index = j*xsize;
if ( firstlat > lastlat )
{
if ( j == 0 )
maxlat = genYmax(grid_center_lat[index], grid_center_lat[index+xsize]);
else
maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);
if ( j == (ysize-1) )
minlat = genYmin(grid_center_lat[index], grid_center_lat[index-xsize]);
else
minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
}
else
{
if ( j == 0 )
minlat = genYmin(grid_center_lat[index], grid_center_lat[index+xsize]);
else
minlat = 0.5*(grid_center_lat[index]+grid_center_lat[index-xsize]);
if ( j == (ysize-1) )
maxlat = genYmax(grid_center_lat[index], grid_center_lat[index-xsize]);
else
maxlat = 0.5*(grid_center_lat[index]+grid_center_lat[index+xsize]);
}
}
for ( i = 0; i < xsize; ++i )
{
index = (j<<2)*xsize + (i<<2);
grid_corner_lat[index ] = minlat;
grid_corner_lat[index+1] = minlat;
grid_corner_lat[index+2] = maxlat;
grid_corner_lat[index+3] = maxlat;
}
}
}
void gridGenRotBounds(int gridID, int nx, int ny,
double *xbounds, double *ybounds, double *xbounds2D, double *ybounds2D)
{
......@@ -254,7 +375,7 @@ void gridGenRotBounds(int gridID, int nx, int ny,
}
static
void gridGenXbounds2D(int nx, int ny, const double * restrict xbounds, double * restrict xbounds2D)
void gridGenXbounds2D(long nx, long ny, const double* restrict xbounds, double* restrict xbounds2D)
{
long i, j, index;
double minlon, maxlon;
......@@ -281,7 +402,7 @@ void gridGenXbounds2D(int nx, int ny, const double * restrict xbounds, double *
}
static
void gridGenYbounds2D(int nx, int ny, const double * restrict ybounds, double * restrict ybounds2D)
void gridGenYbounds2D(long nx, long ny, const double* restrict ybounds, double* restrict ybounds2D)
{
long i, j, index;
double minlat, maxlat;
......
......@@ -19,6 +19,17 @@ void grid_to_radian(const char *units, long nvals, double *restrict values, cons
void grid_to_degree(const char *units, long nvals, double *restrict values, const char *description);
void grid_gen_corners(long n, const double* restrict vals, double* restrict corners);
void grid_cell_center_to_bounds_X2D(const char* xunitstr, long xsize, long ysize,
const double* restrict grid_center_lon, double* restrict grid_corner_lon, double dlon);
void grid_cell_center_to_bounds_Y2D(const char* yunitstr, long xsize, long ysize,
const double* restrict grid_center_lat, double* restrict grid_corner_lat);
void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *xvals);
void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *yvals);
int gridWeights(int gridID, double *weights);
int gridGenArea(int gridID, double *area);
void gaussaw(double pa[], double pw[], int nlat);
int referenceToGrid(int gridID);
int gridToZonal(int gridID);
......
......@@ -223,7 +223,7 @@ double huiliers_area(int num_corners, double *cell_corner_lon, double *cell_corn
}
int gridGenArea(int gridID, double *area)
int gridGenArea(int gridID, double* area)
{
int status = 0;
int gridtype;
......@@ -231,9 +231,9 @@ int gridGenArea(int gridID, double *area)
int lgriddestroy = FALSE;
long i;
long nv, gridsize;
double *grid_corner_lon = NULL;
double *grid_corner_lat = NULL;
int *grid_mask = NULL;
int* grid_mask = NULL;
double* grid_corner_lon = NULL;
double* grid_corner_lat = NULL;
gridsize = gridInqSize(gridID);
gridtype = gridInqType(gridID);
......@@ -314,13 +314,12 @@ int gridGenArea(int gridID, double *area)
{
if ( lgrid_gen_bounds )
{
char xunitstr[CDI_MAX_NAME];
char yunitstr[CDI_MAX_NAME];
int nlon = gridInqXsize(gridID);
int nlat = gridInqYsize(gridID);
double dlon = 0;
if ( nlon == 1 )
{
dlon = 1;
}
if ( nlon == 1 ) dlon = 1;
double *grid_center_lon = NULL;
double *grid_center_lat = NULL;
......@@ -331,8 +330,11 @@ int gridGenArea(int gridID, double *area)
gridInqXvals(gridID, grid_center_lon);
gridInqYvals(gridID, grid_center_lat);
genXbounds(nlon, nlat, grid_center_lon, grid_corner_lon, dlon);
genYbounds(nlon, nlat, grid_center_lat, grid_corner_lat);
gridInqXunits(gridID, xunitstr);
gridInqYunits(gridID, yunitstr);
grid_cell_center_to_bounds_X2D(xunitstr, nlon, nlat, grid_center_lon, grid_corner_lon, dlon);
grid_cell_center_to_bounds_Y2D(yunitstr, nlon, nlat, grid_center_lat, grid_corner_lat);
free(grid_center_lon);
free(grid_center_lat);
......
......@@ -2,6 +2,12 @@
#include <omp.h>
#endif
#include <math.h>
#define PI M_PI
#define PI2 2.0*PI
#define PIH 0.5*PI
#define REMAP_GRID_BASIS_SRC 1
#define REMAP_GRID_BASIS_TGT 2
......@@ -184,3 +190,7 @@ void read_remap_scrip(const char *interp_file, int gridID1, int gridID2, int *ma
int *remap_order, remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv);
void store_link_bilin(remapvars_t *rv, int dst_add, int *restrict src_add, double *restrict weights);
void calc_bin_addr(long gridsize, long nbins, const restr_t* restrict bin_lats, const restr_t* restrict cell_bound_box, int* restrict bin_addr);
void calc_lat_bins(remapgrid_t* src_grid, remapgrid_t* tgt_grid, int map_type);
#include "cdo.h"
#include "remap.h"
void calc_bin_addr(long gridsize, long nbins, const restr_t* restrict bin_lats, const restr_t* restrict cell_bound_box, int* restrict bin_addr)
{
long n, n2, nele, nele4;
restr_t cell_bound_box_lat1, cell_bound_box_lat2;
for ( n = 0; n < nbins; ++n )
{
n2 = n<<1;
bin_addr[n2 ] = gridsize;
bin_addr[n2+1] = 0;
}
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
private(n, n2, nele4, cell_bound_box_lat1, cell_bound_box_lat2) \
shared(gridsize, nbins, bin_lats, cell_bound_box, bin_addr)
#endif
for ( nele = 0; nele < gridsize; ++nele )
{
nele4 = nele<<2;
cell_bound_box_lat1 = cell_bound_box[nele4 ];
cell_bound_box_lat2 = cell_bound_box[nele4+1];
for ( n = 0; n < nbins; ++n )
{
n2 = n<<1;
if ( cell_bound_box_lat1 <= bin_lats[n2+1] &&
cell_bound_box_lat2 >= bin_lats[n2 ] )
{
/*
#if defined(_OPENMP)
if ( nele < bin_addr[n2 ] || nele > bin_addr[n2+1] )
#pragma omp critical
#endif
*/
{
bin_addr[n2 ] = MIN(nele, bin_addr[n2 ]);
bin_addr[n2+1] = MAX(nele, bin_addr[n2+1]);
}
}
}
}
}
/*
static
void calc_bin_addr(long gridsize, long nbins, const restr_t* restrict bin_lats, const restr_t* restrict cell_bound_box, int* restrict bin_addr)
{
long n, n2, nele, nele4;
for ( n = 0; n < nbins; ++n )
{
n2 = n<<1;
bin_addr[n2 ] = gridsize;
bin_addr[n2+1] = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
private(nele4) \
shared(n2, gridsize, bin_lats, cell_bound_box, bin_addr)
#endif
for ( nele = 0; nele < gridsize; ++nele )
{
nele4 = nele<<2;
if ( cell_bound_box[nele4 ] <= bin_lats[n2+1] &&
cell_bound_box[nele4+1] >= bin_lats[n2 ] )
{
bin_addr[n2 ] = MIN(nele, bin_addr[n2 ]);
bin_addr[n2+1] = MAX(nele, bin_addr[n2+1]);
}
}
}
}
*/
void calc_lat_bins(remapgrid_t* src_grid, remapgrid_t* tgt_grid, int map_type)
{
long nbins;
long n; /* Loop counter */
long n2;
double dlat; /* lat/lon intervals for search bins */
restr_t *bin_lats = NULL;
nbins = src_grid->num_srch_bins;
dlat = PI/nbins;
if ( cdoVerbose ) cdoPrint("Using %d latitude bins to restrict search.", nbins);
if ( nbins > 0 )
{
bin_lats = src_grid->bin_lats = realloc(src_grid->bin_lats, 2*nbins*sizeof(restr_t));
for ( n = 0; n < nbins; ++n )
{
n2 = n<<1;
bin_lats[n2 ] = RESTR_SCALE((n )*dlat - PIH);
bin_lats[n2+1] = RESTR_SCALE((n+1)*dlat - PIH);
}
src_grid->bin_addr = realloc(src_grid->bin_addr, 2*nbins*sizeof(int));
calc_bin_addr(src_grid->size, nbins, bin_lats, src_grid->cell_bound_box, src_grid->bin_addr);
if ( map_type == MAP_TYPE_CONSERV || map_type == MAP_TYPE_CONSPHERE )
{
tgt_grid->bin_addr = realloc(tgt_grid->bin_addr, 2*nbins*sizeof(int));
calc_bin_addr(tgt_grid->size, nbins, bin_lats, tgt_grid->cell_bound_box, tgt_grid->bin_addr);
free(src_grid->bin_lats); src_grid->bin_lats = NULL;
}
}
if ( map_type == MAP_TYPE_CONSPHERE )
{
free(tgt_grid->cell_bound_box); tgt_grid->cell_bound_box = NULL;
}
if ( map_type == MAP_TYPE_DISTWGT )
{
free(src_grid->cell_bound_box); src_grid->cell_bound_box = NULL;
}
}
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