diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index 35f1a2c3..00000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "tools/ocean_model_grid_generator"] - path = tools/ocean_model_grid_generator - url = https://github.com/NOAA-GFDL/ocean_model_grid_generator.git diff --git a/Makefile.am b/Makefile.am index 8af0bab0..696bcdd9 100644 --- a/Makefile.am +++ b/Makefile.am @@ -52,10 +52,6 @@ SUBDIRS = postprocessing/combine_restarts \ tools/simple_hydrog/rmvpr \ tools/nc_null_check -if ENABLE_OCEAN_MODEL_GRID_GENERATOR -SUBDIRS += tools/ocean_model_grid_generator -endif - if WITH_CHECK_PROGS SUBDIRS += t endif diff --git a/README.md b/README.md index 319165ec..cf26a607 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ These tools were largely written by members of the GFDL primarily for use in the [Flexible Modeling System](https://www.gfdl.noaa.gov/fms) (FMS) [Runtime Environment](https://www.gfdl.noaa.gov/fre) (FRE) supporting the -work of the +work of the [Geophysical Fluid Dynamics Laboratory](https://www.gfdl.noaa.gov) (GFDL). @@ -58,16 +58,8 @@ The tools available in FRE-NCtools are: ### Other Tools -The [Ocean Model Grid Generator](https://github.com/NOAA-GFDL/ocean_model_grid_generator) -is a collection of tools for creating finite element spherical tripolar grids for -GFDL's MOM based ocean models. Unlike the other tools, NCTools includes it as a submodule, -and also it is a Python3 project. Because of the former attribute, recursive -cloning (see below) is recommended. Because of the latter attribute, -the users python environment may need modification and/or options to the -autotools configure command may need to be specified. This includes the -``--disable-ocean-model-grid-generator`` option (default is enable) -and the ```--enable-venv``` option to build a Python venv containing the -Ocean Model Grid Generator script and all python dependencies. +The [Ocean Model Grid Generator](https://github.com/NOAA-GFDL/ocean_model_grid_generator) can be copied or cloned from its GFDL homepage. + ### User Documentation Documentation on using individual tools may be obtained by running @@ -98,24 +90,13 @@ Contribution requirements include : * The passing of existing unit tests when they are run in the CI system. * The (potential) addition of unit tests when adding new functionality. * For external projects, unit tests may be required and the unit tests - should be run on the external CI system. + should be run on the external CI system. Additionally, since NCTools is distributable via the Spack package manager, the NCTools team will need to be able to compile and distribute via Spack any submodules and their dependencies. Contributors are encouraged to provide Spack recipes for their projects. -## Cloning and submodules -The NCTools github repository contains Ocean Model Grid Generator's repository -as a submodule. After cloning NCTools, it must be initialized and updated for -its submodules: - -``` -git clone --recursive https://github.com/NOAA-GFDL/FRE-NCtools -cd FRE-NCtools -git submodule update --init --recursive -``` - ## Building and Installation - General Information FRE-NCtools has a collection of C and Fortran sources. Within GFDL, FRE-NCtools is built using a recent version of the GNU and Intel C and Fortran compilers. @@ -164,26 +145,13 @@ configure --help=recursive It is common to compile into a build directory (e.g. named `build`) and install into an installation directory (e.g. with full path ``). -If the ocean_model_grid_generator is desired, it may be convenient to allow -the build system to set up a Python venv. These three choices can be -done with these steps: -``` -cd FRE-NCtools -autoreconf -i -mkdir build && cd build -../configure --prefix= --enable-venv -make -make install -``` - -If the ocean_model_grid_generator is not desired, a similar configuration would -be achieved by : +For example: ``` cd FRE-NCtools autoreconf -i mkdir build && cd build -../configure --prefix= --disable-ocean-model-grid-generator +../configure --prefix= make make install ``` diff --git a/configure.ac b/configure.ac index 46c5a111..5931508f 100644 --- a/configure.ac +++ b/configure.ac @@ -22,7 +22,7 @@ AC_PREREQ([2.63]) AC_INIT( [FRE NC Tools], - [2022.02], + [2023.01], [oar.gfdl.help@noaa.gov], [fre-nctools], [https://github.com/NOAA-GFDL/FRE-NCtools]) @@ -53,15 +53,6 @@ AC_ARG_ENABLE([quad-precision], [], [enable_quad_precision=no]) -AC_ARG_ENABLE([ocean-model-grid-generator], - [AS_HELP_STRING([--enable-ocean-model-grid-generator], - [build and install the ocean model grid generator tool])]) -AS_IF([test ${enable_ocean_model_grid_generator:-yes} = yes], - [enable_ocean_model_grid_generator=yes], - [enable_ocean_model_grid_generator=no]) -AM_CONDITIONAL([ENABLE_OCEAN_MODEL_GRID_GENERATOR], - [test $enable_ocean_model_grid_generator = yes]) - AC_PROG_CC([icc gcc]) AM_PROG_CC_C_O ## When autoconf v2.70 is more available, this can be replaced with: @@ -139,6 +130,9 @@ if test "$with_netcdf_fortran" = "no"; then fi AC_LANG_POP([Fortran]) +# Check if Linux sched_getaffinity is available +AC_CHECK_FUNCS([sched_getaffinity], [], []) + # Check for typedefs, structures, and compiler characteristics # Check if the C compiler supports a working `long double` with more range # or precision than `double`. @@ -160,9 +154,6 @@ AC_DEFINE([GIT_REVISION], "git_revision", AC_DEFINE([GIT_HEADHASH], "git_hashval", [Holds the 'git rev-parse HEAD' information if configure ran within a git working directory]) -AS_IF([test $enable_ocean_model_grid_generator = yes], - AC_CONFIG_SUBDIRS([tools/ocean_model_grid_generator]) -) #Code for setting rpath based ncview's configure.in code. echo "ac_computer_gnu: $ac_compiler_gnu" diff --git a/cpp/libfrencutils/create_xgrid.C b/cpp/libfrencutils/create_xgrid.C index 3e30aa40..8d1326b3 100644 --- a/cpp/libfrencutils/create_xgrid.C +++ b/cpp/libfrencutils/create_xgrid.C @@ -37,6 +37,7 @@ #include "DITree.h" #include "mosaic_util.h" #include "create_xgrid.h" +#include "create_xgrid_aux.h" #include "constant.h" #define AREA_RATIO_THRESH (1.e-6) @@ -45,6 +46,15 @@ #define EPSLN30 (1.0e-30) #define EPSLN10 (1.0e-10) +using namespace std; +using namespace nct; +using BBox_t = nct::BBox3D; +using BPair_t = nct::BoxAndId; +using Poly_t = nct::MeshPolygon; +using Point_t = nct::Point3D; +using std::vector; +using std::string; + double grid_box_radius(const double *x, const double *y, const double *z, int n); double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1, const double *x2, const double *y2, const double *z2, int n2); @@ -2341,13 +2351,6 @@ int inside_edge(double x0, double y0, double x1, double y1, double x, double y) // Code below for using search algorithms in create_xgrid_... -using namespace std; -using namespace nct; -using BBox_t = nct::BBox3D; -using BPair_t = nct::BoxAndId; -using Poly_t = nct::MeshPolygon; -using Point_t = nct::Point3D; - inline size_t pt_idx(const size_t i, const size_t j, const size_t nx) { return ( j * nx + i); @@ -2436,8 +2439,9 @@ size_t latlons_outside_ccd_domain(const unsigned int NV4, const double *yv, dou * @param debugf * @return */ + BBox_t getBoxForSphericalPolygon(const double lat_m[], const double lon_m[], - const array &is, bool debugf) { + const array &is, bool debugf) { constexpr unsigned int NV4{4}; // xlons are the longitudes that define the X=0 and Y=0 planes(see ll2xyz function) // where its possible to have an extrema in X or Y when an edge crosses them. @@ -2521,6 +2525,7 @@ BBox_t getBoxForSphericalPolygon(const double lat_m[], const double lon_m[], } + /** * Generate the exchange grid between two grids for the 2nd order * conservative interpolation. This version explicitly uses a search @@ -2634,6 +2639,7 @@ void create_xgrid_2dx2d_order2_ws(const int nlon_in, const int nlat_in, const i for (auto l = 0; l < n2_in; l++) x2_in[l] -= TPI; } + //Call the 2D_by_2D clipping algorithm auto n_out = clip_2dx2d(x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out); if (n_out > 0) { @@ -2687,7 +2693,12 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl }else{ - create_xgrid_2dx2d_order2_ws(nlon_in, nlat_in, nlon_out, nlat_out, + /*create_xgrid_2dx2d_order2_ws(nlon_in, nlat_in, nlon_out, nlat_out, + lon_in, lat_in, lon_out, lat_out, mask_in, + i_in_r, j_in_r, i_out_r, j_out_r, + xgrid_area_r, xgrid_clon_r, xgrid_clat_r);*/ + + create_xgrid_2dx2d_order2_bfwbb(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in_r, j_in_r, i_out_r, j_out_r, xgrid_area_r, xgrid_clon_r, xgrid_clat_r); @@ -2730,6 +2741,7 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl } } +//TODO: this is to be removed: void create_xgrid_2dx2d_order2_check(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, diff --git a/cpp/libfrencutils/create_xgrid.h b/cpp/libfrencutils/create_xgrid.h index fcfcbaf8..9a8dcbc6 100644 --- a/cpp/libfrencutils/create_xgrid.h +++ b/cpp/libfrencutils/create_xgrid.h @@ -23,14 +23,22 @@ #define MAXXGRID 1e7 #endif -#include -#include "BBox3D.h" -#define MV 50 #include +#include #include +#define MV 50 + +#include "BBox3D.h" +using std::vector; + + + /* this value is small compare to earth area */ + + + double poly_ctrlon(const double lon[], const double lat[], int n, double clon); double poly_ctrlat(const double lon[], const double lat[], int n); double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon); @@ -79,15 +87,20 @@ int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int double *xgrid_area, double *xgrid_clon, double *xgrid_clat); void latlon2xyz(const double lat, const double lon, std::array & v); +std::array +get_cell_idxs_ccw_4(const size_t i, const size_t j, const size_t nx); nct::BBox3D getBoxForSphericalPolygon(const double lat_m[], const double lon_m[], - const std::array &is, bool debugf = false); + const std::array &is, bool debugf = false); int search_grids(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, std::vector> & results1) ; -int create_xgrid_2dx2d_order2_ws(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, - const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, - const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, - double *xgrid_area, double *xgrid_clon, double *xgrid_clat); + +void create_xgrid_2dx2d_order2_ws(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, + const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, + const double *mask_in, std::vector& i_in, std::vector& j_in, + std::vector& i_out, std::vector& j_out, std::vector& xgrid_area, + std::vector& xgrid_clon, std::vector& xgrid_clat); + void create_xgrid_2dx2d_order2_ws_check(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, int nxgrid1); @@ -103,5 +116,11 @@ void create_xgrid_2dx2d_order2_check(const int nlon_in, const int nlat_in, const template void printPolygon(std::ostream &os, std::span lonv, std::span latv) ; +void create_xgrid_2dx2d_order2_bfwbb(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, + const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, + const double *mask_in, std::vector& i_in, std::vector& j_in, + std::vector& i_out, std::vector& j_out, std::vector& xgrid_area, + std::vector& xgrid_clon, std::vector& xgrid_clat); + #endif diff --git a/cpp/libfrencutils/create_xgrid_aux.h b/cpp/libfrencutils/create_xgrid_aux.h index 37aeaa45..b12b9d09 100644 --- a/cpp/libfrencutils/create_xgrid_aux.h +++ b/cpp/libfrencutils/create_xgrid_aux.h @@ -8,10 +8,17 @@ #include #include #include + +#include "constant.h" +#include "mpp.h" +#include "create_xgrid.h" + + #include "BBox3D.h" +#include "BoxedObj.h" #include "Polygon.h" #include "mosaic_util.h" -#include "create_xgrid.h" + bool checkBBoxViaPolySamples(std::span yv, std::span xv, nct::BBox3D & box, unsigned int npoints1D = 5 ) { @@ -35,11 +42,11 @@ bool checkBBoxViaPolySamples(std::span yv, std::span xv, auto contains_point = nct::BBox3D::contains(box, pt); if (!contains_point) { passed = false; - //TODO: the polygon and the box can be saved to strstream - printPolygon(std::cout, xv, yv); - std::cout << box << std::endl; - auto str = std::format("< {:16.10e}, {:16.10e}, {:16.10e}>", pt[0],pt[1],pt[2]); - std::cout << str <(std::cout, xv, yv); + // std::cout << box << std::endl; + // auto str = std::format("< {:16.10e}, {:16.10e}, {:16.10e}>", pt[0],pt[1],pt[2]); + // std::cout << str < #include #include +#include +#include +#include "mpp.h" +#include "BBox3D.h" +#include "Polygon.h" +#include "BoxedObj.h" #include "create_xgrid.h" +//#include "cartesian_product.hpp" + +/** + * README IMPORTANT + * Most of the functions in this file are direct copies of those from the + * create_xgrid.C file. The reason is that the nvc++ compiler using C++ std + * parallelism requires many functions to be defined in the same file as the + * function targeted for GPU execution (e.g. create_xgrid_2dx2d_order2_bfwbb + * near the bottom of this file). The amount of source copying within this file + * will be improved ASAP. + */ + using std::sin; +using std::vector; +using std::string; +using std::array; + +using BBox_t = nct::BBox3D; +using BPair_t = nct::BoxAndId; +using Point_t = nct::Point3D; + #include "constant.h" #ifndef MAXXGRID @@ -25,16 +51,50 @@ using std::sin; #define AREA_RATIO_THRESH (1.e-6) #define MASK_THRESH (0.5) -// #define EPSLN8 (1.e-8) #define EPSLN30 (1.0e-30) -// #define EPSLN10 (1.0e-10) constexpr double MU_TOLORENCE{1.e-6}; constexpr double SMALL_VALUE_PA{1.0E-10}; -// These functions are defined here as they were moved from the .C file to the the h -// file in order to use the std parallelization for GPUs with nvc++ 23.x compiers: +//TODO: place elsewhere +//NOTE: some of these are are negative values: +constexpr double FILL_VALUE_DOUBLE{ -std::numeric_limits::max()}; +constexpr int FILL_VALUE_INT{ std::numeric_limits::max() }; +constexpr size_t FILL_VALUE_SIZE_T{ std::numeric_limits::max() }; + + +void gpu_error(const string& str) +{ + gpu_error(str.c_str()); +} + +void gpu_error(const char * str) +{ + fprintf(stderr, "Error from GPU: %s\n", str ); + //exit(1); //TODO: +} + +void latlon2xyz_gpu(const double lat, const double lon, std::array & v){ + v[0] = RADIUS * cos(lat) * cos(lon ); + v[1] = RADIUS * cos(lat) * sin(lon); + v[2] = RADIUS * sin(lat); +} + + +size_t latlons_outside_ccd_domain_gpu(const unsigned int NV4, const double *yv, double *xv) { + size_t count {0}; + for (unsigned int i = 0; i < NV4; i++) { + if (xv[i] == 2 * M_PI) xv[i] = 0; + if (xv[i] >= 2 * M_PI || xv[i] < 0.) { + count++; + } + if (yv[i] < -M_PI_2 || yv[i] > M_PI_2) { + count++; + } + } + return count; +} extern double maxval_double_gpu(int size, const double *data) { int n; double maxval; @@ -365,6 +425,7 @@ int clip_2dx2d_gpu(const double lon1_in[], const double lat1_in[], int n1_in, int gttwopi = 0; /* clip polygon with each boundary of the polygon */ /* We treat lon1_in/lat1_in as clip polygon and lon2_in/lat2_in as subject polygon */ + //Unexpected branch type (/home/mzuniga/nct_search/FRE-NCtools/cpp/libfrencutils/create_xgrid_gpu.C: 428) n_out = n1_in; for (i1 = 0; i1 < n1_in; i1++) { lon_tmp[i1] = lon1_in[i1]; @@ -432,263 +493,388 @@ int clip_2dx2d_gpu(const double lon1_in[], const double lat1_in[], int n1_in, return (n_out); } + +inline +size_t pt_idx_gpu(const size_t i, const size_t j, const size_t nx) { + return ( j * nx + i); +} +std::array +get_cell_idxs_ccw_4_gpu(const size_t i, const size_t j, const size_t nx) { + std::array idxs; + idxs[0] = pt_idx_gpu(i, j, nx); //ll + idxs[1] = pt_idx_gpu(i + 1, j , nx); //lr + idxs[2] = pt_idx_gpu(i + 1, j + 1, nx); //ur + idxs[3] = pt_idx_gpu(i, j + 1, nx);//ul + return idxs; +} + +std::tuple + getPolygonMinMax_gpu(const double lat_m[], const double lon_m[], + const array &is, bool debugf=false) { + constexpr unsigned int NV4{4}; + auto tempy = lat_m[is[0]]; + double maxl = tempy; + double minl = tempy; + for (auto i = 1; i < NV4; i++) { + tempy = lat_m[is[i]]; + if (tempy > maxl) maxl = tempy; + if (tempy < minl) minl = tempy; + } + return {minl, maxl}; +} + + +BBox_t getBoxForSphericalPolygon_gpu(const double lat_m[], const double lon_m[], + const array &is, bool debugf=false) { + constexpr unsigned int NV4{4}; + // xlons are the longitudes that define the X=0 and Y=0 planes(see ll2xyz function) + // where its possible to have an extrema in X or Y when an edge crosses them. + constexpr array xlons {0, M_PI_2, M_PI, 3. * M_PI_2, 2. * M_PI}; + std::array pt{}; // a 3D point. + BBox_t box; //the thing we are calculating. + bool crosses_equator = false; + double yv[ NV4 ], xv[ NV4 ]; //the latitudes (yv) and longitudes (xv) of one polygon/cell. + for (auto i = 0; i< NV4 ; i++) { + xv[i] = lon_m[is[i]]; + yv[i] = lat_m[is[i]]; + } + + auto occd = latlons_outside_ccd_domain_gpu(4, yv, xv); + if(occd>0){ + gpu_error("There are lats or lons outside CCD"); + } + + //Expand the bounding box in consideration of the max and mins of the lats and lons: + const auto [miny_it, maxy_it] = std::minmax_element(std::begin(yv), std::end(yv)); + const double miny = *miny_it; + const double maxy = *maxy_it; + const auto [minx_it, maxx_it] = std::minmax_element(std::begin(xv), std::end(xv)); + const double minx = *minx_it; + const double maxx = *maxx_it; + latlon2xyz_gpu( miny, minx, pt); + box.expand( pt ); + latlon2xyz_gpu( maxy, minx, pt); + box.expand( pt ); + latlon2xyz_gpu( miny, maxx, pt); + box.expand( pt ); + latlon2xyz_gpu( maxy, maxx, pt); + box.expand( pt ); + + //For case where a pole is inside a polygon. Note that + // any Z=k plane of the bbox calculated from above + // (or even from the actual vertices) + // would be sufficient to determine this: + latlon2xyz_gpu(M_PI_2, 0.0, pt); //the north pole + if (BBox_t::contains_zk(box, pt) && (yv[0] >= 0)){ + box.expand(pt); + } + latlon2xyz_gpu(-M_PI_2, 0.0, pt); //the south pole + if (BBox_t::contains_zk(box, pt) && (yv[0] < 0)){ + box.expand(pt); + } + + //Expand the box for coordinate value extrema when edge crosses the equator: + for (auto i = 0; i < NV4; i++) { + auto equator_lat = 0.; + auto t1 = yv[i]; + auto t2 = yv[(i + 1) % NV4]; + if (std::signbit(t1) != std::signbit(t2 )) { + crosses_equator = true; + latlon2xyz_gpu(equator_lat, xv[i], pt); + box.expand(pt); + latlon2xyz_gpu(equator_lat, xv[(i + 1) % NV4], pt); + box.expand(pt); + } + } + + //Expand the box for coordinate value extrema when edge crosses XZ or YZ plane + for (auto i = 0; i < NV4; i++) { + auto t1 = xv[i]; + auto t2 = xv[(i + 1) % NV4]; + for (const auto &xang: xlons) { + if (std::signbit(t1 - xang) != std::signbit(t2 - xang)) { + latlon2xyz_gpu(yv[i], xang, pt); + box.expand(pt); + latlon2xyz_gpu(yv[(i + 1) % NV4], xang, pt); + box.expand(pt); + if(crosses_equator) { + latlon2xyz_gpu(0.0, xang, pt); + box.expand(pt); + } + } + } + } + box.expand_for_doubles(); + return box; +} + +// grid suffix synonyms: '2'/out/target +// grid suffix synonyms: '1'/in/source + +//TODO: This formula for the maximum number of cells in one grid to any other in another grid is +// merely a heuristic. A better one might involve max and min cell areas. Consider how many +//cells from one grid can overlap the largest cell from another. The factor of 1000 below +//should take care of really wierd cells. +size_t get_max_cell_nns(const size_t nx1, const size_t nx2, const size_t ny1, const size_t ny2) +{ return 1000 * std::max(nx1 * ny1, nx2 * ny2) / std::min(nx1 * ny1, nx2 * ny2); } + +//Return the max number of near neighbors for a grid. +size_t +get_max_grid_nns(const size_t nx1, const size_t nx2, const size_t ny1, const size_t ny2) +{ return 10 * std::max(nx1 * ny1, nx2 * ny2); } + +/* + * create_xgrid_2dx2d_order2 - brute force with bounding box usage + */ + +void print(std::tuple t) +{ + const auto& [a, b] = t; + std::cout << '(' << a << ' ' << b << ')' << std::endl; +} + +std::tuple +index_pair_from_combo( const int idx_pair, const int nxy2){ + int ij1 = idx_pair / nxy2; //First index + int ij2 = idx_pair % nxy2; + return {ij1, ij2}; +} + int create_xgrid_2dx2d_order2_legacy_gpu(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { + return (0); +} + +//V3 + +/* + * function create_xgrid_2dx2d_order2 (bfwbb) + * This function uses the std::copy_if function to collect the [ij1 , ij2] pairs of indexes of + * those polygons that pass the "is near neighbors" filters. The filters arr bounding-box intersections + * (and (possibly TBD) also lat/lon overlap tests). The sole purpose of using copy_if + * is to parallelize the brute-force (O(n1 x n2) ) computations. Note that copy_if is merely (with the + * cartesian_product class or the iota class) generating the possible index pairs and copying (or saving) + * only those pairs that pass the neighbors tests. + * + * ISSUES: There have been some trouble with copy_if using cartesian_product class (which at this + * time is not yet available in the C++23 lib). Additionally, there was a compiler bug in creating + * iotas templated to size_t; and this has forced all the following: + * A) A single iota is used (instead of the cartesian_product). A trick was used to combine + * two integers into one (e.g. ij = i * MAXJ + j). When later needed the 'i' and 'j' are recouped from + * from 'ij' (using modulo and division operators). + * B) The iota is templated on int. Because of the trick in A) above, it would have been better to + * use size_t, but unfortunately that does not compile with nvc++. + * C) To make sure the combine ints fit into one int, an outer loop was introduced. So now a + * combine int is made up of ij1 and i2. + * TODO: Change to the C++23 cartesian_product when avaialble. + * NOTE: Recall C++ is row-major order, so column index (J) should be outer index. + */ + +void create_xgrid_2dx2d_order2_bfwbb(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, + const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, + const double *mask_in, vector& i_in, vector& j_in, + vector& i_out, vector& j_out, vector& xgrid_area, + vector& xgrid_clon, vector& xgrid_clat) { #define MAX_V 8 - int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid; - double *area_in, *area_out; - int nblocks = 1; - int *istart2 = NULL, *iend2 = NULL; - int npts_left, nblks_left, pos, m, npts_my, ij; - double *lon_out_min_list, *lon_out_max_list, *lon_out_avg, *lat_out_min_list, *lat_out_max_list; - double *lon_out_list, *lat_out_list; - int *pnxgrid = NULL, *pstart; - int *pi_in = NULL, *pj_in = NULL, *pi_out = NULL, *pj_out = NULL; - double *pxgrid_area = NULL, *pxgrid_clon = NULL, *pxgrid_clat = NULL; - int *n2_list; - int nthreads, nxgrid_block_max; - - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; - nx1p = nx1 + 1; - nx2p = nx2 + 1; - - area_in = (double *)malloc(nx1 * ny1 * sizeof(double)); - area_out = (double *)malloc(nx2 * ny2 * sizeof(double)); - get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in); - get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out); - - nthreads = 1; - - nblocks = nthreads; - - istart2 = (int *)malloc(nblocks * sizeof(int)); - iend2 = (int *)malloc(nblocks * sizeof(int)); - - pstart = (int *)malloc(nblocks * sizeof(int)); - pnxgrid = (int *)malloc(nblocks * sizeof(int)); - - nxgrid_block_max = MAXXGRID / nblocks; - - for (m = 0; m < nblocks; m++) { - pnxgrid[m] = 0; - pstart[m] = m * nxgrid_block_max; + const size_t nx1{(size_t) nlon_in}, nx2{(size_t) nlon_out}, ny1{(size_t) nlat_in}, ny2{(size_t) nlat_out}; + const size_t nx1p{nx1 + 1}; + const size_t nx2p{nx2 + 1}; + + //These two are sized int because their use in copy_if + const int nxy1{nx1 * ny1}; + const int nxy2{nx2 * ny2}; + + if ((size_t) nxy1 * ny2 >= std::numeric_limits::max()) { + /* Grids are too big - but only for the current 23.7 nvc++ limitations on using copy_if + * with iotas templated with ints (i.e. cant use size_t). + */ + std::cerr << "Grids too big nxy1 * ny2 ) >= std::numeric_limits::max()) " << std::endl; + //exit here is possible } + const int max_grid_nns = get_max_grid_nns(nx1, nx2, ny1, ny2); + std::cout << "max_grid_nns estimated at :" << max_grid_nns << std::endl; + + std::vector area_in(nx1 * ny1); + std::vector area_out(nx2 * ny2); + //TODO: note get_grid_area can be parallelized. + get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in.data()); + get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out.data()); + + // grid suffix synonyms: '2'/out/target + // grid suffix synonyms: '1'/in/source + //The out or target pairs + vector boxes_2(nx2 * ny2); + auto ids2 = std::views::iota(0); + std::for_each_n(std::execution::par, ids2.begin(), nx2 * ny2, + [=, boxes_2 = boxes_2.data()](int ij) { + auto i2 = ij % nx2; + auto j2 = ij / nx2; + auto ip = get_cell_idxs_ccw_4_gpu(i2, j2, nx2p); + boxes_2[ij] = getBoxForSphericalPolygon_gpu(lat_out, lon_out, ip); + }); + //The "in" or source pairs + vector boxes_1(nx1 * ny1); + auto ids1 = std::views::iota(0); + std::for_each_n(std::execution::par, ids1.begin(), nx1 * ny1, + [=, boxes_1 = boxes_1.data()](int ij) { + auto i1 = ij % nx1; + auto j1 = ij / nx1; + auto ip = get_cell_idxs_ccw_4_gpu(i1, j1, nx1p); + boxes_1[ij] = getBoxForSphericalPolygon_gpu(lat_in, lon_in, ip); + }); - if (nblocks == 1) { - pi_in = i_in; - pj_in = j_in; - pi_out = i_out; - pj_out = j_out; - pxgrid_area = xgrid_area; - pxgrid_clon = xgrid_clon; - pxgrid_clat = xgrid_clat; - } else { - pi_in = (int *)malloc(MAXXGRID * sizeof(int)); - pj_in = (int *)malloc(MAXXGRID * sizeof(int)); - pi_out = (int *)malloc(MAXXGRID * sizeof(int)); - pj_out = (int *)malloc(MAXXGRID * sizeof(int)); - pxgrid_area = (double *)malloc(MAXXGRID * sizeof(double)); - pxgrid_clon = (double *)malloc(MAXXGRID * sizeof(double)); - pxgrid_clat = (double *)malloc(MAXXGRID * sizeof(double)); + std::cout << "BBox array sizes: " << boxes_1.size() << " ; " << boxes_2.size() << std::endl; + + vector> nn_pairs; + for (int i2 = 0; i2 < nx2; i2++) { + vector nn_pairs1(max_grid_nns, FILL_VALUE_INT); + //An iota that generates all the continuous integers in [0, nxy1 * nx2): + auto g12_ids = std::views::iota((int) 0, (int) (nxy1 * ny2)); + std::copy_if(std::execution::par, //TODO: experiment with par_unseq + g12_ids.begin(), g12_ids.end(), + nn_pairs1.begin(), + [=, boxes_1 = boxes_1.data(), boxes_2 = boxes_2.data()](auto ij12) { + auto [ij1, j2] = index_pair_from_combo(ij12, ny2); + auto ij2 = j2 * nx2 + i2; + return (nct::BBox3D::intersect(boxes_1[ij1], boxes_2[ij2])); + }); + + //TRIM nn_pairs to only hold the real pairs. + for (auto &ij12: nn_pairs1) { + if (ij12 == FILL_VALUE_INT) break; + auto [ij1, j2] = index_pair_from_combo(ij12, ny2); + int ij2 = j2 * nx2 + i2; + nn_pairs.push_back({(int)ij1, ij2}); + } } - npts_left = nx2 * ny2; - nblks_left = nblocks; - pos = 0; - for (m = 0; m < nblocks; m++) { - istart2[m] = pos; - npts_my = npts_left / nblks_left; - iend2[m] = istart2[m] + npts_my - 1; - pos = iend2[m] + 1; - npts_left -= npts_my; - nblks_left--; - } + //NOTE For reproducibility with baseline, results are ordered by increasing : i1; j1; ij2 + auto iidx_cmp = [](const std::tuple& a, const std::tuple& b) -> bool + { + if(std::get<0>(a) == std::get<0>(b)){ + return std::get<1>(a) < std::get<1>(b); + }else{ + return std::get<0>(a) < std::get<0>(b); + } + }; + std::sort(nn_pairs.begin(), nn_pairs.end(), iidx_cmp); + + + //std::cout << "nn_pairs1[0]= " << nn_pairs1[0] << " nn_pairs1[size()-1]"<< nn_pairs1[nn_pairs1.size()-1] << std::endl; + //std::cout << "nn_pairs[0]= " << nn_pairs[0] << " nn_pairs[size()-1]"<< nn_pairs[nn_pairs.size()-1] << std::endl; + //std::cout << " nn_pairs.size() : " << nn_pairs.size() << std::endl; + + //Given the initial neighbor pairs, calculate the final pair intersections. + + std::vector xarea_v(max_grid_nns, FILL_VALUE_DOUBLE); + std::vector clon_v(max_grid_nns, FILL_VALUE_DOUBLE); + std::vector clat_v(max_grid_nns, FILL_VALUE_DOUBLE); + vectori_in_v(max_grid_nns, FILL_VALUE_INT); + vectorj_in_v(max_grid_nns, FILL_VALUE_INT); + vectori_out_v(max_grid_nns, FILL_VALUE_INT); + vectorj_out_v(max_grid_nns, FILL_VALUE_INT); + + //NOTE: Parallelizing the loop below may not buy much since size of nn_pairs should + // already be linear in the size of the larger grid. + for(size_t ij = 0; ij < nn_pairs.size(); ij++) { + auto [ij1, ij2] = nn_pairs[ij]; + // auto [ij1, ij2] = index_pair_from_combo(nn_pairs[ij], ny2); + auto i1 = ij1 % nx1; + auto j1 = ij1 / nx1; + auto i2 = ij2 % nx2; + auto j2 = ij2 / nx2; + + //TODO: refactor into a function ? + double x1_in[MAX_V], y1_in[MAX_V], x2_in[MAX_V], y2_in[MAX_V], x_out[MAX_V], y_out[MAX_V]; + auto idxs_1 = get_cell_idxs_ccw_4_gpu(i1, j1, nx1p); + for (int i = 0; i < 4; i++) { + x1_in[i] = lon_in[idxs_1[i]]; + y1_in[i] = lat_in[idxs_1[i]]; + } - lon_out_min_list = (double *)malloc(nx2 * ny2 * sizeof(double)); - lon_out_max_list = (double *)malloc(nx2 * ny2 * sizeof(double)); - lat_out_min_list = (double *)malloc(nx2 * ny2 * sizeof(double)); - lat_out_max_list = (double *)malloc(nx2 * ny2 * sizeof(double)); - lon_out_avg = (double *)malloc(nx2 * ny2 * sizeof(double)); - n2_list = (int *)malloc(nx2 * ny2 * sizeof(int)); - lon_out_list = (double *)malloc(MAX_V * nx2 * ny2 * sizeof(double)); - lat_out_list = (double *)malloc(MAX_V * nx2 * ny2 * sizeof(double)); - - auto ids = std::views::iota(0, (nx2 * ny2)); - std::for_each_n(std::execution::par, - ids.begin(), (nx2 * ny2), - [=](int ij) { - // for (ij = 0; ij < nx2 * ny2; ij++) { - int i2, j2, n, n0, n1, n2, n3, n2_in, l; - double x2_in[MV], y2_in[MV]; - i2 = ij % nx2; - j2 = ij / nx2; - n = j2 * nx2 + i2; - n0 = j2 * nx2p + i2; - n1 = j2 * nx2p + i2 + 1; - n2 = (j2 + 1) * nx2p + i2 + 1; - n3 = (j2 + 1) * nx2p + i2; - x2_in[0] = lon_out[n0]; - y2_in[0] = lat_out[n0]; - x2_in[1] = lon_out[n1]; - y2_in[1] = lat_out[n1]; - x2_in[2] = lon_out[n2]; - y2_in[2] = lat_out[n2]; - x2_in[3] = lon_out[n3]; - y2_in[3] = lat_out[n3]; - - lat_out_min_list[n] = minval_double_gpu(4, y2_in); - lat_out_max_list[n] = maxval_double_gpu(4, y2_in); - n2_in = fix_lon_gpu(x2_in, y2_in, 4, M_PI); - if (n2_in > MAX_V) error_handler_gpu("create_xgrid.c: n2_in is greater than MAX_V"); - lon_out_min_list[n] = minval_double_gpu(n2_in, x2_in); - lon_out_max_list[n] = maxval_double_gpu(n2_in, x2_in); - lon_out_avg[n] = avgval_double_gpu(n2_in, x2_in); - n2_list[n] = n2_in; - for (l = 0; l < n2_in; l++) { - lon_out_list[n * MAX_V + l] = x2_in[l]; - lat_out_list[n * MAX_V + l] = y2_in[l]; - } - }); + auto idxs_2 = get_cell_idxs_ccw_4_gpu(i2, j2, nx2p); + //The legacy polygon lat-lon representation: + for (int i = 0; i < 4; i++) { + x2_in[i] = lon_out[idxs_2[i]]; + y2_in[i] = lat_out[idxs_2[i]]; + } - nxgrid = 0; - - for (m = 0; m < nblocks; m++) { - auto ids = std::views::iota(0, ny1); - // for(j1=0; j1 MASK_THRESH) { - int n0, n1, n2, n3, l, n1_in; - double lat_in_min, lat_in_max, lon_in_min, lon_in_max, lon_in_avg; - double x1_in[MV], y1_in[MV], x_out[MV], y_out[MV]; - - n0 = j1 * nx1p + i1; - n1 = j1 * nx1p + i1 + 1; - n2 = (j1 + 1) * nx1p + i1 + 1; - n3 = (j1 + 1) * nx1p + i1; - x1_in[0] = lon_in[n0]; - y1_in[0] = lat_in[n0]; - x1_in[1] = lon_in[n1]; - y1_in[1] = lat_in[n1]; - x1_in[2] = lon_in[n2]; - y1_in[2] = lat_in[n2]; - x1_in[3] = lon_in[n3]; - y1_in[3] = lat_in[n3]; - lat_in_min = minval_double_gpu(4, y1_in); - lat_in_max = maxval_double_gpu(4, y1_in); - n1_in = fix_lon_gpu(x1_in, y1_in, 4, M_PI); - lon_in_min = minval_double_gpu(n1_in, x1_in); - lon_in_max = maxval_double_gpu(n1_in, x1_in); - lon_in_avg = avgval_double_gpu(n1_in, x1_in); - for (int ij = istart2[m]; ij <= iend2[m]; ij++) { - int n_out, i2, j2, n2_in; - double xarea, dx, lon_out_min, lon_out_max; - double x2_in[MAX_V], y2_in[MAX_V]; - - i2 = ij % nx2; - j2 = ij / nx2; - - if (lat_out_min_list[ij] >= lat_in_max || lat_out_max_list[ij] <= lat_in_min) continue; - /* adjust x2_in according to lon_in_avg*/ - n2_in = n2_list[ij]; - for (l = 0; l < n2_in; l++) { - x2_in[l] = lon_out_list[ij * MAX_V + l]; - y2_in[l] = lat_out_list[ij * MAX_V + l]; - } - lon_out_min = lon_out_min_list[ij]; - lon_out_max = lon_out_max_list[ij]; - dx = lon_out_avg[ij] - lon_in_avg; - if (dx < -M_PI) { - lon_out_min += TPI; - lon_out_max += TPI; - for (l = 0; l < n2_in; l++) x2_in[l] += TPI; - } else if (dx > M_PI) { - lon_out_min -= TPI; - lon_out_max -= TPI; - for (l = 0; l < n2_in; l++) x2_in[l] -= TPI; - } - - /* x2_in should in the same range as x1_in after lon_fix, so no need to - consider cyclic condition - */ - if (lon_out_min >= lon_in_max || lon_out_max <= lon_in_min) continue; - if ((n_out = clip_2dx2d_gpu(x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out)) > 0) { - double min_area; - int nn; - xarea = poly_area_gpu(x_out, y_out, n_out) * mask_in[j1 * nx1 + i1]; - min_area = std::min(area_in[j1 * nx1 + i1], area_out[j2 * nx2 + i2]); - if (xarea / min_area > AREA_RATIO_THRESH) { - pnxgrid[m]++; - if (pnxgrid[m] >= MAXXGRID / nthreads) - error_handler_gpu( - "nxgrid is greater than MAXXGRID/nthreads, increase MAXXGRID, decrease nthreads, or increase number of MPI ranks"); - nn = pstart[m] + pnxgrid[m] - 1; - pxgrid_area[nn] = xarea; - pxgrid_clon[nn] = poly_ctrlon_gpu(x_out, y_out, n_out, lon_in_avg); - pxgrid_clat[nn] = poly_ctrlat_gpu(x_out, y_out, n_out); - pi_in[nn] = i1; - pj_in[nn] = j1; - pi_out[nn] = i2; - pj_out[nn] = j2; - } - } - } - } - } - }); + //Some polys require "fixing" before calling area (and clipping?) function. + auto n1_in = fix_lon_gpu(x1_in, y1_in, 4, M_PI); + auto n2_in = fix_lon_gpu(x2_in, y2_in, 4, M_PI); //TODO: does this one get fixed + auto lon_in_avg = avgval_double_gpu(n1_in, x1_in); + + //Call the 2D_by_2D clipping algorithm + auto n_out = clip_2dx2d_gpu(x1_in, y1_in, n1_in, x2_in, + y2_in, n2_in, x_out, y_out); + if (n_out > 0) { + auto xarea = poly_area_gpu(x_out, y_out, n_out) * mask_in[j1 * nx1 + i1]; + auto min_area = std::min(area_in[j1 * nx1 + i1], area_out[j2 * nx2 + i2]); + if (xarea / min_area > AREA_RATIO_THRESH) { + xarea_v[ij] = xarea; + clon_v[ij] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + clat_v[ij] = poly_ctrlat(x_out, y_out, n_out); + i_in_v[ij] = i1; //stores the lower corner of poly + j_in_v[ij] = j1; + i_out_v[ij] = i2; + j_out_v[ij] = j2; + } //else leave in the FILL_VALUE + } } - /*copy data if nblocks > 1 */ - if (nblocks == 1) { - nxgrid = pnxgrid[0]; - pi_in = NULL; - pj_in = NULL; - pi_out = NULL; - pj_out = NULL; - pxgrid_area = NULL; - pxgrid_clon = NULL; - pxgrid_clat = NULL; - } else { - int nn, i; - nxgrid = 0; - for (m = 0; m < nblocks; m++) { - for (i = 0; i < pnxgrid[m]; i++) { - nn = pstart[m] + i; - i_in[nxgrid] = pi_in[nn]; - j_in[nxgrid] = pj_in[nn]; - i_out[nxgrid] = pi_out[nn]; - j_out[nxgrid] = pj_out[nn]; - xgrid_area[nxgrid] = pxgrid_area[nn]; - xgrid_clon[nxgrid] = pxgrid_clon[nn]; - xgrid_clat[nxgrid] = pxgrid_clat[nn]; - nxgrid++; - } + //Save final answers (though copy_if can be used again); but don't use more memory than needed. + auto nxgrid = count_if(xarea_v.begin(), xarea_v.end(), [](double x) { return x > 0; }); + + xgrid_area.reserve(nxgrid); //TODO: is assign(capacity(, FILL_VALUE) needed? + i_in.reserve(nxgrid); + j_in.reserve(nxgrid); + i_out.reserve(nxgrid); + j_out.reserve(nxgrid); + xgrid_clon.reserve(nxgrid); + xgrid_clat.reserve(nxgrid); + + //auto ir = 0; + for(int i = 0; i < nn_pairs.size(); i ++) { + if (xarea_v[i] > 0.0) { + xgrid_area.emplace_back(xarea_v[ i ]); + xgrid_clon.emplace_back(clon_v[i]); + xgrid_clat.emplace_back(clat_v[i]); + i_in.emplace_back(i_in_v[i]); + j_in.emplace_back(j_in_v[i]); + i_out .emplace_back(i_out_v[i]); + j_out.emplace_back(j_out_v[i]); } - free(pi_in); - free(pj_in); - free(pi_out); - free(pj_out); - free(pxgrid_area); - free(pxgrid_clon); - free(pxgrid_clat); } - free(area_in); - free(area_out); - free(lon_out_min_list); - free(lon_out_max_list); - free(lat_out_min_list); - free(lat_out_max_list); - free(lon_out_avg); - free(n2_list); - free(lon_out_list); - free(lat_out_list); - - return nxgrid; + std::cout <<"create_xgrid_2dx2d_order2_bff2 end; xgrid_area.size= " << xgrid_area.size() < lats_min_1(nx1 * ny1); +vector lats_max_1(nx1 * ny1); +ids1 = std::views::iota(0); +std::for_each_n(std::execution::par, ids1.begin(), nx1 * ny1, +[=, lats_min_1 = lats_min_1.data(), lats_max_1 = lats_max_1.data()](int ij) { +auto i1 = ij % nx1; +auto j1 = ij / nx1; +auto ip = get_cell_idxs_ccw_4_gpu(i1, j1, nx1p); +auto [minl, maxl] = getPolygonMinMax_gpu(lat_out, lon_out, ip); +lats_min_1[ij] = minl; +lats_max_1[ij] = maxl; +}); + + vector lats_min_2(nx2 * ny2); + vector lats_max_2(nx2 * ny2); + for (int ij = 0; ij < nxy2; ij++) { + auto i2 = ij % nx2; + auto j2 = ij / nx2; + auto ip = get_cell_idxs_ccw_4_gpu(i2, j2, nx2p); + auto [minl, maxl] = getPolygonMinMax_gpu(lat_out, lon_out, ip); + lats_min_2[ij] = minl; + lats_max_2[ij] = maxl; + } +*/ \ No newline at end of file diff --git a/cpp/search/BBox3D.h b/cpp/search/BBox3D.h index bb8a23ba..f28a9a88 100644 --- a/cpp/search/BBox3D.h +++ b/cpp/search/BBox3D.h @@ -30,6 +30,15 @@ namespace nct { inline float getLo(int dim) { return lo[dim]; } inline float getHi(int dim) { return hi[dim]; } + void operator=(const BBox3D& b) + { + for (auto i = 0; i < 3; i++) { + lo[i] = b.lo[i]; + hi[i] = b.hi[i]; + } + } + + //Constructor from a mesh of polygons, which in turn is really // just a vector of pointers to points. // This is the constructor for use with legacy nctools diff --git a/site-configs/gfdl-ws/config.site b/site-configs/gfdl-ws/config.site index 2d3cc342..413ee7b0 100644 --- a/site-configs/gfdl-ws/config.site +++ b/site-configs/gfdl-ws/config.site @@ -21,7 +21,6 @@ # # configure options test -z "$with_mpi" && with_mpi=yes -test -z "$enable_venv" && enable_venv=yes # Standard prefix location for FRE distribution test "$prefix" = NONE && prefix=/home/fms/local/opt/fre-nctools/${PACKAGE_VERSION}/gfdl-ws diff --git a/site-configs/gfdl/config.site b/site-configs/gfdl/config.site index 0c59604c..e4d33ac8 100644 --- a/site-configs/gfdl/config.site +++ b/site-configs/gfdl/config.site @@ -21,7 +21,6 @@ # # configure options test -z "$with_mpi" && with_mpi=yes -test -z "$enable_venv" && enable_venv=yes # Standard prefix location for FRE distribution test "$prefix" = NONE && prefix=/home/fms/local/opt/fre-nctools/${PACKAGE_VERSION}/gfdl diff --git a/site-configs/ncrc/config.site b/site-configs/ncrc/config.site index d4f8722b..ae3e45da 100644 --- a/site-configs/ncrc/config.site +++ b/site-configs/ncrc/config.site @@ -21,7 +21,6 @@ # # configure options test -z "$with_mpi" && with_mpi=yes -test -z "$enable_venv" && enable_venv=yes # Standard prefix location for FRE distribution test "$prefix" = NONE && prefix=/ncrc/home2/fms/local/opt/fre-nctools/${PACKAGE_VERSION}/ncrc diff --git a/site-configs/ncrc5/config.site b/site-configs/ncrc5/config.site index d4f8722b..ae3e45da 100644 --- a/site-configs/ncrc5/config.site +++ b/site-configs/ncrc5/config.site @@ -21,7 +21,6 @@ # # configure options test -z "$with_mpi" && with_mpi=yes -test -z "$enable_venv" && enable_venv=yes # Standard prefix location for FRE distribution test "$prefix" = NONE && prefix=/ncrc/home2/fms/local/opt/fre-nctools/${PACKAGE_VERSION}/ncrc diff --git a/t/Test01-check_programs_exist.sh b/t/Test01-check_programs_exist.sh index 3d8d9478..480d5a3d 100755 --- a/t/Test01-check_programs_exist.sh +++ b/t/Test01-check_programs_exist.sh @@ -22,42 +22,37 @@ load test_utils @test "Check make_hgrid exists and is executable" { - run command -v make_hgrid - [ "$status" -eq 0 ] + run_and_check -v make_hgrid run make_hgrid -h [ "$status" -eq 2 ] } @test "Check make_vgrid exists and is executable" { - run command -v make_vgrid - [ "$status" -eq 0 ] + run_and_check -v make_vgrid run make_vgrid -h [ "$status" -eq 1 ] } @test "Check make_solo_mosaic exists and is executable" { - run command -v make_solo_mosaic - [ "$status" -eq 0 ] + run_and_check -v make_solo_mosaic run make_solo_mosaic -h [ "$status" -eq 2 ] } @test "Check make_topog exists and is executable" { - run command -v make_topog - [ "$status" -eq 0 ] + run_and_check -v make_topog run make_topog -h [ "$status" -eq 1 ] } # only run with mpi @test "Check make_topog_parallel exists and is executable" { [ ! -z $skip_mpi ] && skip "not built with MPI" - run command -v make_topog_parallel - [ "$status" -eq 0 ] + run_and_check -v make_topog_parallel + } @test "Check coupler_mosaic exists and is executable" { - run command -v make_coupler_mosaic - [ "$status" -eq 0 ] + run_and_check -v make_coupler_mosaic run make_coupler_mosaic -h [ "$status" -eq 2 ] } @@ -65,13 +60,12 @@ load test_utils # only run with mpi @test "Check coupler_mosaic_parallel exists and is executable" { [ ! -z $skip_mpi ] && skip "not built with MPI" - run command -v make_coupler_mosaic_parallel - [ "$status" -eq 0 ] + run_and_check -v make_coupler_mosaic_parallel + } @test "Check fregrid exists and is executable" { - run command -v fregrid - [ "$status" -eq 0 ] + run_and_check -v fregrid run fregrid -h [ "$status" -eq 2 ] } @@ -79,34 +73,30 @@ load test_utils # only run with mpi @test "Check fregrid_parallel exists and is executable" { [ ! -z $skip_mpi ] && skip "not built with MPI" - run command -v fregrid_parallel - [ "$status" -eq 0 ] + run_and_check -v fregrid_parallel + } @test "Check runoff_regrid exists and is executable" { - run command -v runoff_regrid - [ "$status" -eq 0 ] + run_and_check -v runoff_regrid run runoff_regrid -h [ "$status" -eq 2 ] } @test "Check river_regrid exists and is executable" { - run command -v river_regrid - [ "$status" -eq 0 ] + run_and_check -v river_regrid run river_regrid -h [ "$status" -eq 2 ] } @test "Check check_mask exists and is executable" { - run command -v check_mask - [ "$status" -eq 0 ] + run_and_check -v check_mask run check_mask -h [ "$status" -eq 2 ] } @test "Check remap_land exists and is executable" { - run command -v remap_land - [ "$status" -eq 0 ] + run_and_check -v remap_land run remap_land -h [ "$status" -eq 1 ] } @@ -114,51 +104,44 @@ load test_utils # only run with mpi @test "Check remap_land_parallel exists and is executable" { [ ! -z $skip_mpi ] && skip "not built with MPI" - run command -v remap_land_parallel - [ "$status" -eq 0 ] + run_and_check -v remap_land_parallel + } @test "Check make_regional_mosaic exists and is executable" { - run command -v make_regional_mosaic - [ "$status" -eq 0 ] + run_and_check -v make_regional_mosaic run make_regional_mosaic -h [ "$status" -eq 2 ] } @test "Check mppncscatter exists and is executable" { - run command -v mppncscatter - [ "$status" -eq 0 ] + run_and_check -v mppncscatter } @test "Check mppnccombine exists and is executable" { - run command -v mppnccombine - [ "$status" -eq 0 ] + run_and_check -v mppnccombine run mppnccombine -h [ "$status" -eq 1 ] } @test "Check combine-ncc exists and is executable" { - run command -v combine-ncc - [ "$status" -eq 0 ] + run_and_check -v combine-ncc } @test "Check decompress-ncc exists and is executable" { - run command -v decompress-ncc - [ "$status" -eq 0 ] + run_and_check -v decompress-ncc } @test "Check cr_lake_files exists and is executable" { - run command -v cr_lake_files - [ "$status" -eq 0 ] + run_and_check -v cr_lake_files + } @test "Check cp_river_vars exists and is executable" { - run command -v cp_river_vars - [ "$status" -eq 0 ] + run_and_check -v cp_river_vars } @test "Check rmv_parallel_rivers exists and is executable" { - run command -v rmv_parallel_rivers - [ "$status" -eq 0 ] + run_and_check -v rmv_parallel_rivers } diff --git a/t/Test02-mppnccombine.sh b/t/Test02-mppnccombine.sh index 5550d841..6785e76a 100755 --- a/t/Test02-mppnccombine.sh +++ b/t/Test02-mppnccombine.sh @@ -30,6 +30,6 @@ load test_utils mppnccombine \ mppnccombine_output.nc \ mppnccombine.nc.???? - [ -e mppnccombine_output.nc ] + run_and_check [ -e mppnccombine_output.nc ] ncdump -h mppnccombine_output.nc } diff --git a/t/Test03-grid_coupled_model.sh b/t/Test03-grid_coupled_model.sh index 0af3eb55..90e9a98e 100755 --- a/t/Test03-grid_coupled_model.sh +++ b/t/Test03-grid_coupled_model.sh @@ -76,7 +76,7 @@ load test_utils --vgrid ocean_vgrid.nc \ --output topog_parallel.nc - nccmp -md topog.nc topog_parallel.nc + run_and_check nccmp -md topog.nc topog_parallel.nc fi @@ -113,8 +113,9 @@ load test_utils --ocean_topog ../topog.nc \ --mosaic_name grid_spec \ --area_ratio_thresh 1.e-10 + # compare any created files to non-parallel (exclude directory differences) - nccmp -md --exclude=atm_mosaic_dir --exclude=lnd_mosaic_dir --exclude=ocn_mosaic_dir \ + run_and_check nccmp -md --exclude=atm_mosaic_dir --exclude=lnd_mosaic_dir --exclude=ocn_mosaic_dir \ --exclude=ocn_topog_dir grid_spec.nc ../grid_spec.nc fi } diff --git a/t/Test04-grid_coupled_nest.sh b/t/Test04-grid_coupled_nest.sh index b7a04587..86e2680c 100755 --- a/t/Test04-grid_coupled_nest.sh +++ b/t/Test04-grid_coupled_nest.sh @@ -128,7 +128,7 @@ ncgen -o OCCAM_p5degree.nc $BATS_TEST_DIRNAME/Test03-input/OCCAM_p5degree.ncl --land_mosaic ../land_mosaic.nc --ocean_mosaic ../ocean_mosaic.nc \ --ocean_topog ../topog.nc --interp_order 1 --mosaic_name grid_spec # directory paths should differ - nccmp -md --exclude=atm_mosaic_dir --exclude=lnd_mosaic_dir --exclude=ocn_mosaic_dir \ + run_and_check nccmp -md --exclude=atm_mosaic_dir --exclude=lnd_mosaic_dir --exclude=ocn_mosaic_dir \ --exclude=ocn_topog_dir grid_spec.nc ../grid_spec.nc fi } diff --git a/t/Test12-mppnscatter.sh b/t/Test12-mppnscatter.sh index 30dc5bbd..96ef35da 100755 --- a/t/Test12-mppnscatter.sh +++ b/t/Test12-mppnscatter.sh @@ -48,6 +48,6 @@ prepare_input_data () mppnccombine -64 output/$test_file ${test_file}.???? # Compare output - nccmp -w format -md output/$test_file input/$test_file + run_and_check nccmp -w format -md output/$test_file input/$test_file } diff --git a/t/Test15-regrid_land.sh b/t/Test15-regrid_land.sh index 97ad80cc..cb29c96b 100755 --- a/t/Test15-regrid_land.sh +++ b/t/Test15-regrid_land.sh @@ -64,7 +64,7 @@ load test_utils --output_file out_parallel.nc \ --remap_file remap_file.nc - nccmp -md out.nc out_parallel.nc + run_and_check nccmp -md out.nc out_parallel.nc fi # remap other fields diff --git a/t/Test17-combine-ncc.sh b/t/Test17-combine-ncc.sh index 8db34afd..8a562d57 100755 --- a/t/Test17-combine-ncc.sh +++ b/t/Test17-combine-ncc.sh @@ -29,7 +29,7 @@ load test_utils combine-ncc \ combine-ncc.atmos_daily.nc.copy \ combine-ncc_output.nc - [ -e combine-ncc_output.nc ] - ncdump -h combine-ncc_output.nc + run_and_check [ -e combine-ncc_output.nc ] + run_and_check ncdump -h combine-ncc_output.nc } diff --git a/t/Test18-decompress-ncc.sh b/t/Test18-decompress-ncc.sh index efb5856b..596ce56e 100755 --- a/t/Test18-decompress-ncc.sh +++ b/t/Test18-decompress-ncc.sh @@ -29,6 +29,6 @@ load test_utils decompress-ncc \ decompress-ncc.atmos_daily.nc.copy \ decompress-ncc_output.nc - [ -e decompress-ncc_output.nc ] - ncdump -h decompress-ncc_output.nc + run_and_check [ -e decompress-ncc_output.nc ] + run_and_check ncdump -h decompress-ncc_output.nc } diff --git a/t/Test21-reference_mppnccombine.sh b/t/Test21-reference_mppnccombine.sh index 72818613..4697b82b 100755 --- a/t/Test21-reference_mppnccombine.sh +++ b/t/Test21-reference_mppnccombine.sh @@ -33,12 +33,12 @@ load test_utils mppnccombine \ mppnccombine_output.nc \ mppnccombine.nc.???? - [ -e mppnccombine_output.nc ] + run_and_check [ -e mppnccombine_output.nc ] ncdump -h mppnccombine_output.nc - [ -e $top_srcdir/t/Test02-reference/mppnccombine_output.nc ] + run_and_check [ -e $top_srcdir/t/Test02-reference/mppnccombine_output.nc ] nccmp -V - nccmp -d mppnccombine_output.nc $top_srcdir/t/Test02-reference/mppnccombine_output.nc + run_and_check nccmp -d mppnccombine_output.nc $top_srcdir/t/Test02-reference/mppnccombine_output.nc } diff --git a/t/Test22-reference_combine-ncc.sh b/t/Test22-reference_combine-ncc.sh index 4f9e7e22..d4f05d3f 100755 --- a/t/Test22-reference_combine-ncc.sh +++ b/t/Test22-reference_combine-ncc.sh @@ -30,12 +30,12 @@ load test_utils combine-ncc \ combine-ncc.atmos_daily.nc.copy \ combine-ncc_output.nc - [ -e combine-ncc_output.nc ] + run_and_check [ -e combine-ncc_output.nc ] ncdump -h combine-ncc_output.nc - [ -e $top_srcdir/t/Test17-reference/combine-ncc_output.nc ] + run_and_check [ -e $top_srcdir/t/Test17-reference/combine-ncc_output.nc ] - nccmp -V + run_and_check nccmp -V - nccmp -d combine-ncc_output.nc $top_srcdir/t/Test17-reference/combine-ncc_output.nc + run_and_check nccmp -d combine-ncc_output.nc $top_srcdir/t/Test17-reference/combine-ncc_output.nc } diff --git a/t/Test23-reference_decompress-ncc.sh b/t/Test23-reference_decompress-ncc.sh index 66adfbbf..587e39ee 100755 --- a/t/Test23-reference_decompress-ncc.sh +++ b/t/Test23-reference_decompress-ncc.sh @@ -29,12 +29,12 @@ load test_utils decompress-ncc \ decompress-ncc.atmos_daily.nc.copy \ decompress-ncc_output.nc - [ -e decompress-ncc_output.nc ] + run_and_check [ -e decompress-ncc_output.nc ] ncdump -h decompress-ncc_output.nc - [ -e $top_srcdir/t/Test18-reference/decompress-ncc_output.nc ] + run_and_check [ -e $top_srcdir/t/Test18-reference/decompress-ncc_output.nc ] - nccmp -V + run_and_check nccmp -V - nccmp -d decompress-ncc_output.nc $top_srcdir/t/Test18-reference/decompress-ncc_output.nc + run_and_check nccmp -d decompress-ncc_output.nc $top_srcdir/t/Test18-reference/decompress-ncc_output.nc } diff --git a/t/Test24-reference_fregrid.sh b/t/Test24-reference_fregrid.sh index ecdd13af..887a2b07 100755 --- a/t/Test24-reference_fregrid.sh +++ b/t/Test24-reference_fregrid.sh @@ -56,15 +56,15 @@ load test_utils --output_mosaic latlon_mosaic.nc \ --check_conserve - [ -e ocean_temp_salt.res.latlon.nc ] - [ -e $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc ] + run_and_check [ -e ocean_temp_salt.res.latlon.nc ] + run_and_check [ -e $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc ] - nccmp -V + run_and_check nccmp -V - nccmp -d ocean_temp_salt.res.latlon.nc $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc + run_and_check nccmp -d ocean_temp_salt.res.latlon.nc $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc - [ -e latlon_mosaic.nc ] - [ -e $top_srcdir/t/Test20-reference/latlon_mosaic.nc ] + run_and_check [ -e latlon_mosaic.nc ] + run_and_check [ -e $top_srcdir/t/Test20-reference/latlon_mosaic.nc ] - nccmp -d latlon_mosaic.nc $top_srcdir/t/Test20-reference/latlon_mosaic.nc + run_and_check nccmp -d latlon_mosaic.nc $top_srcdir/t/Test20-reference/latlon_mosaic.nc } diff --git a/t/Test25-hydrology.sh b/t/Test25-hydrology.sh index d1945e3e..ffd955ec 100755 --- a/t/Test25-hydrology.sh +++ b/t/Test25-hydrology.sh @@ -47,6 +47,6 @@ load test_utils #Run wrapper hydrology script run $top_srcdir/tools/simple_hydrog/share/make_simple_hydrog.csh -f 0. -t 1.e-5 -m $top_srcdir/t/Test25-input/grid_spec.nc - [ "$status" -eq 0 ] + run_and_check [ "$status" -eq 0 ] } diff --git a/t/Test26-reference_fregrid_gcc.sh b/t/Test26-reference_fregrid_gcc.sh index 4cecfc93..0eaea53b 100755 --- a/t/Test26-reference_fregrid_gcc.sh +++ b/t/Test26-reference_fregrid_gcc.sh @@ -56,16 +56,16 @@ load test_utils --output_mosaic latlon_mosaic.nc \ --check_conserve - [ -e ocean_temp_salt.res.latlon.nc ] - [ -e $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc ] + run_and_check [ -e ocean_temp_salt.res.latlon.nc ] + run_and_check [ -e $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc ] - nccmp -V + run_and_check nccmp -V # checks values with float tolerance to bypass any result differences from compilers - nccmp -d -t 0.000001 ocean_temp_salt.res.latlon.nc $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc + run_and_check nccmp -d -t 0.000001 ocean_temp_salt.res.latlon.nc $top_srcdir/t/Test20-reference/ocean_temp_salt.res.latlon.nc - [ -e latlon_mosaic.nc ] - [ -e $top_srcdir/t/Test20-reference/latlon_mosaic.nc ] + run_and_check [ -e latlon_mosaic.nc ] + run_and_check [ -e $top_srcdir/t/Test20-reference/latlon_mosaic.nc ] - nccmp -d latlon_mosaic.nc $top_srcdir/t/Test20-reference/latlon_mosaic.nc + run_and_check nccmp -d latlon_mosaic.nc $top_srcdir/t/Test20-reference/latlon_mosaic.nc } diff --git a/t/Test31-fregrid_stretched.sh b/t/Test31-fregrid_stretched.sh index 157d8f0b..4c2aa1a6 100755 --- a/t/Test31-fregrid_stretched.sh +++ b/t/Test31-fregrid_stretched.sh @@ -1,65 +1,84 @@ #!/usr/bin/env bats + +#*********************************************************************** +# GNU Lesser General Public License +# +# This file is part of the GFDL FRE NetCDF tools package (FRE-NCTools). +# +# FRE-NCTools is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or (at +# your option) any later version. +# +# FRE-NCTools is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with FRE-NCTools. If not, see +# . +#*********************************************************************** + +load test_utils + # test stretched grid @test "Test stretched grid data lats 32.0 34.0 35.4" { - if [ ! -d "Test31" ] + if [ ! -d "Test31" ] then mkdir Test31 fi cd Test31 - cp $top_srcdir/t/Test31-input/ocean_hgrid.nc . + cp $top_srcdir/t/Test31-input/ocean_hgrid.nc . cp $top_srcdir/t/Test31-input/ocean_mosaic.nc . cp $top_srcdir/t/Test31-input/ocean_topog.nc . -#Make streetched grid - run command make_hgrid \ +#Make streetched grid + run_and_check make_hgrid \ --grid_type gnomonic_ed --do_schmidt \ --stretch_factor 2.5 \ --target_lon 262.4 \ --target_lat 32.0 \ --nlon 512 \ --grid_name C256_grid_32.0 - [ "$status" -eq 0 ] - run command make_hgrid \ + run_and_check make_hgrid \ --grid_type gnomonic_ed --do_schmidt \ --stretch_factor 2.5 \ --target_lon 262.4 \ --target_lat 32.0 \ --nlon 512 \ --grid_name C256_grid_34.0 - [ "$status" -eq 0 ] - run command make_hgrid \ + + run_and_check make_hgrid \ --grid_type gnomonic_ed --do_schmidt \ --stretch_factor 2.5 \ --target_lon 262.4 \ --target_lat 35.4 \ --nlon 512 \ --grid_name C256_grid_35.4 - [ "$status" -eq 0 ] - #Create stretched grid mosaic - run command make_solo_mosaic \ + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_32.0 --tile_file C256_grid_32.0.tile1.nc,C256_grid_32.0.tile2.nc,C256_grid_32.0.tile3.nc,C256_grid_32.0.tile4.nc,C256_grid_32.0.tile5.nc,C256_grid_32.0.tile6.nc - [ "$status" -eq 0 ] - run command make_solo_mosaic \ + + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_34.0 --tile_file C256_grid_34.0.tile1.nc,C256_grid_34.0.tile2.nc,C256_grid_34.0.tile3.nc,C256_grid_34.0.tile4.nc,C256_grid_34.0.tile5.nc,C256_grid_34.0.tile6.nc - [ "$status" -eq 0 ] - run command make_solo_mosaic \ + + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_35.4 --tile_file C256_grid_35.4.tile1.nc,C256_grid_35.4.tile2.nc,C256_grid_35.4.tile3.nc,C256_grid_35.4.tile4.nc,C256_grid_35.4.tile5.nc,C256_grid_35.4.tile6.nc - [ "$status" -eq 0 ] # stretched grid lats 32.0 34.0 35.4 run bash -c 'fregrid \ @@ -75,8 +94,7 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_32_0=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_32_0} \< 0.00001 - [ "$status" -eq 0 ] + run_and_check expr ${var_32_0} \< 0.00001 run bash -c 'fregrid \ --input_mosaic C256_mosaic_34.0.nc \ @@ -91,8 +109,8 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_34_0=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_34_0} \< 0.00001 - [ "$status" -eq 0 ] + run_and_check expr ${var_34_0} \< 0.00001 + run bash -c 'fregrid \ --input_mosaic C256_mosaic_35.4.nc \ @@ -107,9 +125,8 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_35_4=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_35_4} \< 0.00001 - [ "$status" -eq 0 ] - + run_and_check expr ${var_35_4} \< 0.00001 + cd .. # rm -rf Test31 } diff --git a/t/Test32-fregrid_no_stretched.sh b/t/Test32-fregrid_no_stretched.sh index 232831ac..28e9b76f 100755 --- a/t/Test32-fregrid_no_stretched.sh +++ b/t/Test32-fregrid_no_stretched.sh @@ -1,55 +1,72 @@ #!/usr/bin/env bats -# test no stretched grid + +#*********************************************************************** +# GNU Lesser General Public License +# +# This file is part of the GFDL FRE NetCDF tools package (FRE-NCTools). +# +# FRE-NCTools is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or (at +# your option) any later version. +# +# FRE-NCTools is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with FRE-NCTools. If not, see +# . +#*********************************************************************** + +load test_utils + +# test no stretched grid @test "Test no stretched grid data lats 32.0 34.0 35.4" { - if [ ! -d "Test32" ] + if [ ! -d "Test32" ] then mkdir Test32 fi cd Test32 - cp $top_srcdir/t/Test32-input/ocean_hgrid.nc . + cp $top_srcdir/t/Test32-input/ocean_hgrid.nc . cp $top_srcdir/t/Test32-input/ocean_mosaic.nc . cp $top_srcdir/t/Test32-input/topog.nc . -#Make no stretched grid - run command make_hgrid \ +#Make no stretched grid + run_and_check make_hgrid \ --grid_type gnomonic_ed \ --nlon 512 \ --grid_name C256_grid_32.0 - [ "$status" -eq 0 ] - run command make_hgrid \ + run_and_check make_hgrid \ --grid_type gnomonic_ed \ --nlon 512 \ --grid_name C256_grid_34.0 - [ "$status" -eq 0 ] - run command make_hgrid \ + run_and_check make_hgrid \ --grid_type gnomonic_ed \ --nlon 512 \ --grid_name C256_grid_35.4 - [ "$status" -eq 0 ] #Create no stretched grid mosaic - run command make_solo_mosaic \ + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_32.0 --tile_file C256_grid_32.0.tile1.nc,C256_grid_32.0.tile2.nc,C256_grid_32.0.tile3.nc,C256_grid_32.0.tile4.nc,C256_grid_32.0.tile5.nc,C256_grid_32.0.tile6.nc - [ "$status" -eq 0 ] - run command make_solo_mosaic \ + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_34.0 --tile_file C256_grid_34.0.tile1.nc,C256_grid_34.0.tile2.nc,C256_grid_34.0.tile3.nc,C256_grid_34.0.tile4.nc,C256_grid_34.0.tile5.nc,C256_grid_34.0.tile6.nc - [ "$status" -eq 0 ] - run command make_solo_mosaic \ + run_and_check make_solo_mosaic \ --num_tiles 6 \ --dir ./ \ --mosaic_name C256_mosaic_35.4 --tile_file C256_grid_35.4.tile1.nc,C256_grid_35.4.tile2.nc,C256_grid_35.4.tile3.nc,C256_grid_35.4.tile4.nc,C256_grid_35.4.tile5.nc,C256_grid_35.4.tile6.nc - [ "$status" -eq 0 ] # no stretched grid lats 32.0 34.0 35.4 run bash -c 'fregrid \ @@ -65,8 +82,7 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_32_0=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_32_0} \< 0.00001 - [ "$status" -eq 0 ] + run_and_check expr ${var_32_0} \< 0.00001 run bash -c 'fregrid \ --input_mosaic C256_mosaic_34.0.nc \ @@ -81,8 +97,8 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_34_0=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_34_0} \< 0.00001 - [ "$status" -eq 0 ] + run_and_check expr ${var_34_0} \< 0.00001 + run bash -c 'fregrid \ --input_mosaic C256_mosaic_35.4.nc \ @@ -97,8 +113,7 @@ | awk 'NR==4' | rev | cut -c39-49 | rev' var_35_4=$(echo $output | awk '{ print sprintf("%.9f", $1); }') echo $output - expr ${var_35_4} \< 0.00001 - [ "$status" -eq 0 ] + run_and_check expr ${var_35_4} \< 0.00001 cd .. # rm -rf Test32 diff --git a/t/Test33-reference_make_hgrid.sh b/t/Test33-reference_make_hgrid.sh index 1ef229a1..25c17def 100755 --- a/t/Test33-reference_make_hgrid.sh +++ b/t/Test33-reference_make_hgrid.sh @@ -40,7 +40,7 @@ fv3_grids_filelist="" for i in 1 2 3 4 5 6 do fv3_file=$fv3_grid_name".tile"$i".nc" - [ -e $top_srcdir/t/Test33-reference/$fv3_file ] + run_and_check [ -e $top_srcdir/t/Test33-reference/$fv3_file ] done cp $top_srcdir/t/Test33-reference/*.nc . diff --git a/t/test_utils.bash b/t/test_utils.bash index 504c0295..51ecfd4c 100755 --- a/t/test_utils.bash +++ b/t/test_utils.bash @@ -83,8 +83,8 @@ generate_all_from_ncl(){ done } -# run the command (the argument list; 1st in list is the command) using the bats "run command" -# its status function, and theck the status is 0 +# run the command (the argument list; 1st in list is the command) using the bats "run command", +# theck the status is 0, and echo the output if not. function run_and_check() { local cmd="$*" #Expands the list into a single string; spearating parms with space. diff --git a/tools/fregrid/fregrid.c b/tools/fregrid/fregrid.c index 6cf1c54f..66e067ed 100644 --- a/tools/fregrid/fregrid.c +++ b/tools/fregrid/fregrid.c @@ -23,20 +23,20 @@ AUTHOR: Zhi Liang (Zhi.Liang@noaa.gov) NOAA Geophysical Fluid Dynamics Lab, Princeton, NJ - + This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. - + This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + For the full text of the GNU General Public License, write to: Free Software Foundation, Inc., - 675 Mass Ave, Cambridge, MA 02139, USA. + 675 Mass Ave, Cambridge, MA 02139, USA. ----------------------------------------------------------------------- */ #include @@ -73,20 +73,20 @@ char *usage[] = { " [--associated_file_dir dir] [--format format] ", " [--deflation #] [--shuffle 1|0] ", " ", - "fregrid remaps data (scalar or vector) from input_mosaic onto ", - "output_mosaic. Note that the target grid also could be specified ", - "through lonBegin, lonEnd, latBegin, latEnd, nlon and nlat. Currently ", - "only T-cell scalar regridding and AGRID vector regridding is ", - "available. Bilinear interpolation is implemented only for cubic grid ", - "vector interpolation. The interpolation algorithm used is controlled ", - "by --interp_method with default 'conserve_order1'. Currently the ", - "'conserve_order1', 'conserve_order2' and 'bilinear' remapping schemes ", - "are implemented. 'bilinear' is only used to remap data from cubic grid ", - "to latlon grid. Alternative schemes can be added if needed. fregrid ", - "expects NetCDF format input. scalar_field and/or u_field/v_field must ", - "be specified. u_fld and v_fld must be paired together. ", - "If using NetCDF4, output will have the same deflation and shuffle settings ", - "unless specified. ", + "fregrid remaps data (scalar or vector) from an input_mosaic onto an output_mosaic, ", + "but the output (or target) grid can also be specified through the inputs ", + "lonBegin, lonEnd, latBegin, latEnd, nlon and nlat. Currently only T-cell scalar ", + "regridding and AGRID vector regridding (where the vector field is expressed in ", + "spherical coordinates) is available. Vector fields can be mapped ", + "as vectors only from the cubic grid to the lat-lon grid, and while using bilinear ", + "interpolation; but users can instead specify the vector components as (independent) ", + "scalar components and thereby choose a conservative interpolation. The interpolation ", + "algorithm used is controled by --interp_method with default 'conserve_order1'. ", + " Currently the implemented remapping schemes are 'conserve_order1', 'conserve_order2' ", + "and 'bilinear'. Alternative schemes can be added if needed. fregrid expects NetCDF ", + "format input. scalar_field and/or u_field/v_field must be specified. u_fld and v_fld ", + "must be paired together. If using NetCDF4, output will have the same deflation and ", + "shuffle settings unless specified. ", " ", "fregrid takes the following flags: ", " ", @@ -97,7 +97,7 @@ char *usage[] = { " information for each tile. ", " ", "OPTIONAL FLAGS ", - " ", + " ", "--input_file input_file specify the input file name. The suffix '.nc' can be ", " omitted. The suffix 'tile#' should not present for ", " multiple-tile files. The number of files must be 1 for ", @@ -106,15 +106,19 @@ char *usage[] = { " ", "--scalar_field scalar_fld specify the scalar field name to be regridded. The ", " multiple entry field names are seperated by comma. ", - " ", + " ", "--u_field u_fld specify the vector field u-componentname to be ", " regridded. The multiple entry field names are seperated ", " by comma. u_field must be paired together with v_field. ", - " ", + " u_field corresponds to the zonal component of the ", + " vector field ", + " ", "--v_field v_fld specify the vector field v-componentname to be ", " regridded. The multiple entry field names are seperated ", " by comma. v_field must be paired together with u_field. ", - " ", + " v_field corresponds to the meridional component of the ", + " vector field. ", + " ", "--output_mosaic output_mosaic specify the output mosaic information. This file ", " contains list of tile files which specify the grid ", " information for each tile. If output_mosaic is not ", @@ -134,32 +138,32 @@ char *usage[] = { " ", "--latEnd #decimal specify the ending latitude(in degree) of the ", " geographical region of the target grid on which the ", - " output is desired. The default value is 90. ", + " output is desired. The default value is 90. ", " ", "--nlon #integer specify number of grid box cells in x-direction for a ", " regular lat-lon grid. ", " ", "--nlat #integer specify number of grid box cells in y-direction for a ", " regular lat-lon grid. ", - " ", + " ", "--KlevelBegin #integer specify begin index of the k-level (depth axis) that ", " to be regridded. ", " ", "--KlevelEnd #integer specify end index of the k-level (depth axis) that ", - " to be regridded. ", + " to be regridded. ", " ", "--LstepBegin #integer specify the begin index of L-step (time axis) that ", " to be regridded. ", " ", "--LstepEnd #integer specify the end index of L-step (time axis) that ", - " to be regridded. ", - " ", + " to be regridded. ", + " ", "--output_file output_file specify the output file name. If not presented, ", " output_file will take the value of input_file. The ", " suffix '.nc' can be omitted. The suffix 'tile#' should ", " not present for multiple-tile files. The number of ", " files must be 1 for scalar regridding and can be 1 or 2 ", - " for vector regridding. File path should not be includes.", + " for vector regridding. File path should not be includes.", " ", "--input_dir input_dir specify the path that stores input_file. If not ", " presented, the input file is assumed to be stored in ", @@ -262,7 +266,7 @@ char *usage[] = { "--shuffle # If using NetCDF4 , use shuffle if 1 and don't use if 0 ", " Defaults to input file settings. ", " ", - " Example 1: Remap C48 data onto N45 grid. ", + " Example 1: Remap C48 data onto N45 grid. ", " (use GFDL-CM3 data as example) ", " fregrid --input_mosaic C48_mosaic.nc --input_dir input_dir --input_file input_file ", " --scalar_field temp,salt --nlon 144 --nlat 90 ", @@ -303,9 +307,9 @@ int main(int argc, char* argv[]) char *associated_file_dir = NULL; int check_conserve = 0; /* 0 means no check */ double lonbegin = 0, lonend = 360; - double latbegin = -90, latend = 90; + double latbegin = -90, latend = 90; int nlon = 0, nlat = 0; - int kbegin = 0, kend = -1; + int kbegin = 0, kend = -1; int lbegin = 0, lend = -1; char *remap_file = NULL; char interp_method[STRING] = "conserve_order1"; @@ -328,7 +332,7 @@ int main(int argc, char* argv[]) int deflation = -1; int shuffle = -1; char *format=NULL; - + char wt_file_obj[512]; char *weight_file=NULL; char *weight_field = NULL; @@ -351,14 +355,14 @@ int main(int argc, char* argv[]) Interp_config *interp = NULL; /* store remapping information */ int save_weight_only = 0; int nthreads = 1; - + double time_get_in_grid=0, time_get_out_grid=0, time_get_input=0; double time_setup_interp=0, time_do_interp=0, time_write=0; clock_t time_start, time_end; - + int errflg = (argc == 1); int fid; - + static struct option long_options[] = { {"input_mosaic", required_argument, NULL, 'a'}, {"output_mosaic", required_argument, NULL, 'b'}, @@ -404,12 +408,12 @@ int main(int argc, char* argv[]) {"format", required_argument, NULL, 'U'}, {"help", no_argument, NULL, 'h'}, {0, 0, 0, 0}, - }; - + }; + /* start parallel */ mpp_init(&argc, &argv); mpp_domain_init(); - + while ((c = getopt_long(argc, argv, "", long_options, &option_index)) != -1) { switch (c) { case 'a': @@ -430,7 +434,7 @@ int main(int argc, char* argv[]) tokenize(entry, ",", STRING, NFILE, (char *)input_file, &nfiles); break; case 'f': - if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -f"); + if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -f"); strcpy(entry, optarg); tokenize(entry, ",", STRING, NFILE, (char *)output_file, &nfiles_out); break; @@ -438,20 +442,20 @@ int main(int argc, char* argv[]) remap_file = optarg; break; case 's': - if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -s"); + if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -s"); strcpy(entry, optarg); tokenize(entry, ",", STRING, NVAR, (char *)scalar_name, &nscalar); break; case 'u': - if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -u"); + if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -u"); strcpy(entry, optarg); tokenize(entry, ",", STRING, NVAR, (char *)u_name, &nvector); - break; + break; case 'v': - if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -v"); + if(strlen(optarg) >= MAXSTRING) mpp_error("fregrid: the entry is not long for option -v"); strcpy(entry, optarg); tokenize(entry, ",", STRING, NVAR, (char *)v_name, &nvector2); - break; + break; case 'j': strcpy(interp_method, optarg); break; @@ -460,7 +464,7 @@ int main(int argc, char* argv[]) break; case 'k': test_param = atof(optarg); - break; + break; case 'l': opcode |= SYMMETRY; break; @@ -538,7 +542,7 @@ int main(int argc, char* argv[]) break; case 'P': debug = 1; - break; + break; case 'Q': nthreads = atoi(optarg); break; @@ -564,7 +568,7 @@ int main(int argc, char* argv[]) char **u = usage; while (*u) { fprintf(stderr, "%s\n", *u); u++; } exit(2); - } + } /* check the arguments */ if( !mosaic_in ) mpp_error("fregrid: input_mosaic is not specified"); if( !mosaic_out ) { @@ -590,7 +594,7 @@ int main(int argc, char* argv[]) opcode |= MONOTONIC; } else if(!strcmp(interp_method, "bilinear") ) { - if(mpp_pe() == mpp_root_pe())printf("****fregrid: bilinear remapping scheme will be used for regridding.\n"); + if(mpp_pe() == mpp_root_pe())printf("****fregrid: bilinear remapping scheme will be used for regridding.\n"); opcode |= BILINEAR; } else @@ -632,7 +636,7 @@ int main(int argc, char* argv[]) if(weight_field) { if(nvector >0) mpp_error("fregrid: weight_field should not be specified for vector interpolation, contact developer"); if(!weight_file) { - + if(nfiles==0) mpp_error("fregrid: weight_field is specified, but both weight_file and input_file are not specified"); if(dir_in) sprintf(wt_file_obj, "%s/%s", dir_in, input_file[0]); @@ -654,7 +658,7 @@ int main(int argc, char* argv[]) mpp_error("fregrid: shuffle must be 0 (off) or 1 (on)"); if (deflation < -1 || deflation > 9) mpp_error("fregrid: deflation must be between 0 (off) and 9"); - + /* define history to be the history in the grid file */ strcpy(history,argv[0]); @@ -669,7 +673,7 @@ int main(int argc, char* argv[]) else strcat(history, argv[i]); } - + { int base_cpu; @@ -688,7 +692,7 @@ int main(int argc, char* argv[]) mpp_close(fid); /* second order conservative interpolation is only avail for the cubic sphere input grid */ - if( ntiles_in != 6 && (opcode & CONSERVE_ORDER2) ) + if( ntiles_in != 6 && (opcode & CONSERVE_ORDER2) ) mpp_error("fregrid: when the input grid is not cubic sphere grid, interp_method can not be conserve_order2"); if(mosaic_out) { @@ -707,7 +711,7 @@ int main(int argc, char* argv[]) if(check_conserve) opcode |= CHECK_CONSERVE; if( opcode & STANDARD_DIMENSION ) printf("fregrid: --standard_dimension is set\n"); - + if( opcode & BILINEAR ) { int ncontact; ncontact = read_mosaic_ncontacts(mosaic_in); @@ -716,11 +720,11 @@ int main(int argc, char* argv[]) if(ncontact !=12) mpp_error("fregrid: when interp_method is bilinear, the input mosaic should be 12 contact cubic grid"); if(mpp_npes() > 1) mpp_error("fregrid: parallel is not implemented for bilinear remapping"); } - else + else y_at_center = 1; if(extrapolate) opcode |= EXTRAPOLATE; - + /* memory allocation for data structure */ grid_in = (Grid_config *)malloc(ntiles_in *sizeof(Grid_config)); grid_out = (Grid_config *)malloc(ntiles_out*sizeof(Grid_config)); @@ -738,7 +742,7 @@ int main(int argc, char* argv[]) print_mem_usage("After calling get_input_grid"); time_start = clock(); } - if(mosaic_out) + if(mosaic_out) get_output_grid_from_mosaic( ntiles_out, grid_out, mosaic_out, opcode, &great_circle_algorithm_out ); else { great_circle_algorithm_out = 0; @@ -751,7 +755,7 @@ int main(int argc, char* argv[]) print_mem_usage("After calling get_output_grid"); } /* find out if great_circle algorithm is used in the input grid or output grid */ - + if( great_circle_algorithm_in == 0 && great_circle_algorithm_out == 0 ) opcode |= LEGACY_CLIP; else { @@ -763,7 +767,7 @@ int main(int argc, char* argv[]) /* get the grid cell_area */ get_input_output_cell_area(ntiles_in, grid_in, ntiles_out, grid_out, opcode); - if(debug) print_mem_usage("After get_input_output_cell_area"); + if(debug) print_mem_usage("After get_input_output_cell_area"); /* currently extrapolate are limited to ntiles = 1. extrapolate are limited to lat-lon input grid */ if( extrapolate ) { int i, j, ind0, ind1, ind2; @@ -777,7 +781,7 @@ int main(int argc, char* argv[]) if(fabs( grid_in[0].lont[ind0]-grid_in[0].lont[ind2] ) > EPSLN10 || fabs( grid_in[0].latt[ind0]-grid_in[0].latt[ind1] ) > EPSLN10 ) mpp_error("fregrid: extrapolate is limited to lat-lon grid"); - + } } @@ -788,13 +792,13 @@ int main(int argc, char* argv[]) if(vertical_interp) mpp_error("fregrid: vertical_interp is not supported for vector fields"); if(extrapolate) mpp_error("fregrid: extrapolate is not supported for vector fields"); } - - if(remap_file) set_remap_file(ntiles_out, mosaic_out, remap_file, interp, &opcode, save_weight_only); + + if(remap_file) set_remap_file(ntiles_out, mosaic_out, remap_file, interp, &opcode, save_weight_only); if(!save_weight_only) { file_in = (File_config *)malloc(ntiles_in *sizeof(File_config)); file_out = (File_config *)malloc(ntiles_out*sizeof(File_config)); - + if(nfiles == 2) { file2_in = (File_config *)malloc(ntiles_in *sizeof(File_config)); file2_out = (File_config *)malloc(ntiles_out*sizeof(File_config)); @@ -810,10 +814,10 @@ int main(int argc, char* argv[]) get_input_vgrid(&vgrid_in, file_in[0].name, scalar_name[0]); setup_vertical_interp(&vgrid_in, &vgrid_out); } - + if(nfiles == 2) { set_mosaic_data_file(ntiles_in, mosaic_in, dir_in, file2_in, input_file[1]); - set_mosaic_data_file(ntiles_out, mosaic_out, dir_out, file2_out, output_file[1]); + set_mosaic_data_file(ntiles_out, mosaic_out, dir_out, file2_out, output_file[1]); } //Open the input files. Save the nc format of the first one. @@ -848,21 +852,29 @@ int main(int argc, char* argv[]) if(nscalar == 0 && nvector == 0) { if(mpp_pe() == mpp_root_pe()) printf("NOTE from fregrid: no scalar and vector field need to be regridded.\n"); mpp_end(); - return 0; + return 0; } - + if(nscalar > 0) { scalar_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); scalar_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); } if(nvector > 0) { - mpp_error("fregrid: currently does not support vertical interpolation, contact developer"); - u_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); - u_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); - v_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); - v_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); + if ((opcode & CONSERVE_ORDER1) || (opcode & CONSERVE_ORDER2 )){ + mpp_error("fregrid: conservative interpolation of vector fields is not supported. \n" + "Use bilinear interpolation or regrid the vector components independently as scalars.\n"); + } + else if (!(opcode & BILINEAR)){ + mpp_error("fregrid: For vector fields, the interpolation method must be bilinear interpolation.\n"); + } + else{ + u_in = (Field_config *)malloc(ntiles_in * sizeof(Field_config)); + u_out = (Field_config *)malloc(ntiles_out * sizeof(Field_config)); + v_in = (Field_config *)malloc(ntiles_in * sizeof(Field_config)); + v_out = (Field_config *)malloc(ntiles_out * sizeof(Field_config)); + } } - + set_field_struct ( ntiles_in, scalar_in, nscalar, scalar_name_remap[0], file_in); set_field_struct ( ntiles_out, scalar_out, nscalar, scalar_name_remap[0], file_out); set_field_struct ( ntiles_in, u_in, nvector, u_name[0], file_in); @@ -890,7 +902,7 @@ int main(int argc, char* argv[]) }else{ printf("WARNING: fregrid could not set in_format"); } - + set_output_metadata(ntiles_in, nfiles, file_in, file2_in, scalar_in, u_in, v_in, ntiles_out, file_out, file2_out, scalar_out, u_out, v_out, grid_out, &vgrid_out, history, tagname, opcode, deflation, shuffle); @@ -921,7 +933,7 @@ int main(int argc, char* argv[]) break; } } - } + } } /* preparing for the interpolation, if remapping information exist, read it from remap_file, @@ -949,7 +961,7 @@ int main(int argc, char* argv[]) latbegin_in = -0.5*M_PI; else latbegin_in = latbegin*D2R; - + setup_bilinear_interp(ntiles_in, grid_in, ntiles_out, grid_out, interp, opcode, dlon_in, dlat_in, lonbegin_in, latbegin_in ); } else @@ -972,14 +984,14 @@ int main(int argc, char* argv[]) } } mpp_end(); - return 0; + return 0; } - + if(nscalar > 0) { get_field_attribute(ntiles_in, scalar_in); copy_field_attribute(ntiles_out, scalar_in, scalar_out); } - + if(nvector > 0) { get_field_attribute(ntiles_in, u_in); get_field_attribute(ntiles_in, v_in); @@ -988,7 +1000,7 @@ int main(int argc, char* argv[]) } - + /* set time step to 1, only test scalar field now, nz need to be 1 */ if(test_case) { if(nscalar != 1 || nvector != 0) mpp_error("fregrid: when test_case is specified, nscalar must be 1 and nvector must be 0"); @@ -996,14 +1008,14 @@ int main(int argc, char* argv[]) file_in->nt = 1; file_out->nt = 1; } - + /* Then doing the regridding */ for(m=0; mnt; m++) { int memsize, level_z, level_n, level_t; write_output_time(ntiles_out, file_out, m); if(nfiles > 1) write_output_time(ntiles_out, file2_out, m); - + /* first interp scalar variable */ for(l=0; lvar[l].has_taxis && m>0) continue; @@ -1014,7 +1026,7 @@ int main(int argc, char* argv[]) if(extrapolate) { get_input_data(ntiles_in, scalar_in, grid_in, bound_T, l, -1, level_n, level_t, extrapolate, stop_crit); allocate_field_data(ntiles_out, scalar_out, grid_out, scalar_in->var[l].nz); - if( opcode & BILINEAR ) + if( opcode & BILINEAR ) do_scalar_bilinear_interp(interp, l, ntiles_in, grid_in, grid_out, scalar_in, scalar_out, finer_step, fill_missing); else do_scalar_conserve_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, scalar_in, scalar_out, opcode, scalar_in->var[l].nz); @@ -1031,7 +1043,7 @@ int main(int argc, char* argv[]) } else { for(level_z=scalar_in->var[l].kstart; level_z <= scalar_in->var[l].kend; level_z++) - { + { if(debug) time_start = clock(); if(test_case) get_test_input_data(test_case, test_param, ntiles_in, scalar_in, grid_in, bound_T, opcode); @@ -1044,7 +1056,7 @@ int main(int argc, char* argv[]) allocate_field_data(ntiles_out, scalar_out, grid_out, 1); if(debug) time_start = clock(); - if( opcode & BILINEAR ) + if( opcode & BILINEAR ) do_scalar_bilinear_interp(interp, l, ntiles_in, grid_in, grid_out, scalar_in, scalar_out, finer_step, fill_missing); else do_scalar_conserve_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, scalar_in, scalar_out, opcode,1); @@ -1085,7 +1097,7 @@ int main(int argc, char* argv[]) do_vector_bilinear_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, u_in, v_in, u_out, v_out, finer_step, fill_missing); else do_vector_conserve_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, u_in, v_in, u_out, v_out, opcode); - + write_field_data(ntiles_out, u_out, grid_out, l, level_z, level_n, m); write_field_data(ntiles_out, v_out, grid_out, l, level_z, level_n, m); for(n=0; n. + * License along with FMS. If not, see . **********************************************************************/ +#ifndef _GNU_SOURCE #define _GNU_SOURCE +#endif #include #include #include +#include #include #include #include @@ -30,23 +34,35 @@ #include #include "config.h" +#ifdef __APPLE__ +#include +#endif + +/** + * gettid function for systems that do not have this function (e.g. on Mac OS.) + */ #ifndef HAVE_GETTID -#ifndef __APPLE__ static pid_t gettid(void) { - return syscall(__NR_gettid); -} +#if defined(__APPLE__) + uint64_t tid64; + pthread_threadid_np(NULL, &tid64); + pid_t tid = (pid_t)tid64; +#else + pid_t tid = syscall(__NR_gettid); #endif + return tid; +} #endif -/* +/** * Returns this thread's CPU affinity, if bound to a single core, * or else -1. */ int get_cpu_affinity(void) { -#if defined(use_libMPI) && !defined(__APPLE__) +#ifdef HAVE_SCHED_GETAFFINITY cpu_set_t coremask; /* core affinity mask */ CPU_ZERO(&coremask); @@ -64,17 +80,14 @@ int get_cpu_affinity(void) return -1; } -int get_cpu_affinity_(void) { return get_cpu_affinity(); } /* Fortran interface */ - - -/* +/** * Returns this groups CPUSET * and also the CPUSET size or -1 (in case of a storage error) */ int get_cpuset(int fsz, int *output, int pe, _Bool debug) { -#ifndef __APPLE__ - cpu_set_t coremask; /* core affinity mask */ +#ifdef HAVE_SCHED_GETAFFINITY + cpu_set_t coremask; /* core affinity mask */ CPU_ZERO(&coremask); if (sched_getaffinity(gettid(),sizeof(cpu_set_t),&coremask) != 0) { @@ -108,16 +121,13 @@ int get_cpuset(int fsz, int *output, int pe, _Bool debug) #endif } -int get_cpuset_(int *fsz, int *output, int *pe, _Bool *debug) { return get_cpuset(*fsz, output, *pe, *debug); } /* Fortran interface */ - - -/* +/** * Set CPU affinity to one core. */ int set_cpu_affinity(int cpu) { -#ifndef __APPLE__ - cpu_set_t coremask; /* core affinity mask */ +#ifdef HAVE_SCHED_GETAFFINITY + cpu_set_t coremask; /* core affinity mask */ CPU_ZERO(&coremask); CPU_SET(cpu,&coremask); @@ -127,5 +137,3 @@ int set_cpu_affinity(int cpu) #endif return 0; } - -int set_cpu_affinity_(int *cpu) { return set_cpu_affinity(*cpu); } /* Fortran interface */ diff --git a/tools/libfrencutils/constant.h b/tools/libfrencutils/constant.h index f6eacf94..82cc19e6 100644 --- a/tools/libfrencutils/constant.h +++ b/tools/libfrencutils/constant.h @@ -24,8 +24,18 @@ #define RADIUS (6371000.) #define STRING 255 +#include + #ifndef M_PI -#define M_PI (3.14159265358979323846) +#define M_PI (3.14159265358979323846) +#endif + +#ifndef M_PI_2 +#define M_PI_2 (1.57079632679489661923) +#endif + +#ifndef M_SQRT2 +#define M_SQRT2 (1.41421356237309504880) #endif #define R2D (180/M_PI) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index c9780669..69b6acfe 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -235,6 +235,7 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int tmpy[j1*nx1p+i1] = lat_in[j1]; } /* This is just a temporary fix to solve the issue that there is one point in zonal direction */ + // TODO: Finish this "temporary fix" if(nx1 > 1) get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in); else @@ -1284,7 +1285,7 @@ int clip_2dx2d(const double lon1_in[], const double lat1_in[], int n1_in, } //Some grid boxes near North Pole are clipped wrong (issue #42 ) //The following heuristic fix seems to work. Why? - if(gttwopi){pimod(lon_tmp,n1_in);pimod(lon2_tmp,n2_in);} + if(gttwopi){pimod(lon_tmp,n1_in);pimod(lon2_tmp,n2_in);} x2_0 = lon2_tmp[n2_in-1]; y2_0 = lat2_tmp[n2_in-1]; @@ -2785,7 +2786,7 @@ int main(int argc, char* argv[]) case 14: /**************************************************************** test clip_2dx2d_great_cirle case 14: Cubic sphere grid at tile = 3, point (i=24,j=1) - identical grid boxes + identical grid boxes ****************************************************************/ /* nlon1 = 1; @@ -2835,7 +2836,7 @@ int main(int argc, char* argv[]) case 16: /*Must give [[-57.748, -30, -30, -97.6, -97.6], - [89.876, 89.891, 90, 90, 89.9183]]*/ + [89.876, 89.891, 90, 90, 89.9183]]*/ n1_in = 6; n2_in = 5; double lon1_16[] = {82.400, 82.400, 262.400, 262.400, 326.498, 379.641}; double lat1_16[] = {89.835, 90.000, 90.000, 89.847, 89.648, 89.642}; @@ -2877,11 +2878,11 @@ int main(int argc, char* argv[]) /*Must give nothing*/ case 19: /**************************************************************** - test clip_2dx2d 2: two boxes that include the North Pole + test clip_2dx2d 2: two boxes that include the North Pole one has vertices on the tripolar fold the other is totally outside the first This actually happens for some stretched grid - configurations mosaic_c256r25tlat32.0_om4p25 + configurations mosaic_c256r25tlat32.0_om4p25 The test gives wrong answers! ****************************************************************/ n1_in = 6; n2_in = 5; @@ -2898,7 +2899,7 @@ int main(int argc, char* argv[]) case 20: /*Must give - n_out= 5 + n_out= 5 122.176, 150, 150, 82.4, 82.4, 89.761, 89.789, 90, 90, 89.8429, */ n1_in = 6; n2_in = 5; @@ -2915,10 +2916,10 @@ int main(int argc, char* argv[]) case 21: /*Must give - n_out= 5 + n_out= 5 60.000, 82.400, 82.400, 60.000], 89.889, 89.843, 90.000, 90.000] - */ + */ n1_in = 6; n2_in = 5; double lon1_21[] = {82.400, 82.400, 262.400, 262.400, 326.498, 379.641}; double lat1_21[] = {89.835, 90.000, 90.000, 89.847, 89.648, 89.642}; @@ -2934,7 +2935,7 @@ int main(int argc, char* argv[]) case 26: /*Side crosses SP (Right cell). Must give same box - */ + */ n1_in = 4; n2_in = 4; double lon1_22[] = {209.68793552504,158.60256162113,82.40000000000,262.40000000000}; double lat1_22[] = {-89.11514201451,-89.26896927380,-89.82370183256, -89.46584623220}; @@ -2950,8 +2951,8 @@ int main(int argc, char* argv[]) case 23: /*Side does not cross SP (Right cell). Must give same box - */ - + */ + n1_in = 4; n2_in = 4; double lon1_23[] = {158.60256162113,121.19651597620,82.40000000000,82.40000000000}; double lat1_23[] = {-89.26896927380,-88.85737639760,-89.10746816044,-89.82370183256}; @@ -2967,7 +2968,7 @@ int main(int argc, char* argv[]) case 24: /*Side crosses SP (Left cell). Added twin poles. Must give the same box - */ + */ n1_in = 6; n2_in = 6; double lon1_24[] = {262.40000000000,262.40000000000,82.4,82.4,6.19743837887,-44.88793552504}; double lat1_24[] = {-89.46584623220,-90.0, -90.0,-89.82370183256, -89.26896927380, -89.11514201451}; @@ -2980,9 +2981,9 @@ int main(int argc, char* argv[]) memcpy(lat2_in,lat2_24,sizeof(lat2_in)); break; case 25: - /*Side crosses SP (Left cell). + /*Side crosses SP (Left cell). Must givethe same box - */ + */ n1_in = 4; n2_in = 4; double lon1_25[] = {262.40000000000,82.4,6.19743837887,-44.88793552504}; double lat1_25[] = {-89.46584623220, -89.82370183256, -89.26896927380, -89.11514201451}; @@ -2997,7 +2998,7 @@ int main(int argc, char* argv[]) case 22: /*Side does not cross SP (Left cell). Must give same box - */ + */ n1_in = 4; n2_in = 4; double lon1_26[] = {82.4,82.4,43.60348402380,6.19743837887}; double lat1_26[] = {-89.82370183256, -89.10746816044, -88.85737639760, -89.26896927380}; @@ -3110,7 +3111,7 @@ int main(int argc, char* argv[]) printf("\n"); for(i=0; i1.0e14 || area2>1.0e14 || area_out>1.0e14) printf("Error in calculating area !\n"); + if(area1>1.0e14 || area2>1.0e14 || area_out>1.0e14) printf("Error in calculating area !\n"); if(n==16 || n==20) printf("Must result n_out=5!\n"); if(n==21) printf("Must result n_out=4!\n"); if(n==15 || n==17) printf("Must result the second box!\n"); diff --git a/tools/libfrencutils/create_xgrid.h b/tools/libfrencutils/create_xgrid.h index 7d13a005..88fdb048 100644 --- a/tools/libfrencutils/create_xgrid.h +++ b/tools/libfrencutils/create_xgrid.h @@ -33,13 +33,13 @@ double box_ctrlat(double ll_lon, double ll_lat, double ur_lon, double ur_lat); int get_maxxgrid(void); void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); void get_grid_great_circle_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); -void get_grid_area_dimensionless(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); +//void get_grid_area_dimensionless(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); void get_grid_area_no_adjust(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); int clip(const double lon_in[], const double lat_in[], int n_in, double ll_lon, double ll_lat, double ur_lon, double ur_lat, double lon_out[], double lat_out[]); void pimod(double x[],int nn); -int clip_2dx2d(const double lon1_in[], const double lat1_in[], int n1_in, - const double lon2_in[], const double lat2_in[], int n2_in, +int clip_2dx2d(const double lon1_in[], const double lat1_in[], int n1_in, + const double lon2_in[], const double lat2_in[], int n2_in, double lon_out[], double lat_out[]); int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, @@ -65,8 +65,8 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat); -int clip_2dx2d_great_circle(const double x1_in[], const double y1_in[], const double z1_in[], int n1_in, - const double x2_in[], const double y2_in[], const double z2_in [], int n2_in, +int clip_2dx2d_great_circle(const double x1_in[], const double y1_in[], const double z1_in[], int n1_in, + const double x2_in[], const double y2_in[], const double z2_in [], int n2_in, double x_out[], double y_out[], double z_out[]); int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, diff --git a/tools/libfrencutils/mosaic_util.c b/tools/libfrencutils/mosaic_util.c index 0bea93a1..edb56fd5 100644 --- a/tools/libfrencutils/mosaic_util.c +++ b/tools/libfrencutils/mosaic_util.c @@ -37,17 +37,28 @@ #define EPSLN10 (1.e-10) #define EPSLN15 (1.e-15) #define EPSLN30 (1.e-30) -/*********************************************************** - void error_handler(char *str) - error handler: will print out error message and then abort -***********************************************************/ + +const double from_pole_threshold_rad = 0.0174533; // 1.0 deg + int reproduce_siena = 0; +int rotate_poly_flag = 0; +double the_rotation_matrix[3][3] = { 0 }; -void set_reproduce_siena_true(void) -{ +void set_reproduce_siena_true(void){ reproduce_siena = 1; } +void set_rotate_poly_true(void){ + rotate_poly_flag = 1; + set_the_rotation_matrix(); +} + + + +/*********************************************************** + void error_handler(char *str) + error handler: will print out error message and then abort +***********************************************************/ void error_handler(const char *msg) { fprintf(stderr, "FATAL Error: %s\n", msg ); @@ -56,7 +67,7 @@ void error_handler(const char *msg) #else exit(1); #endif -}; /* error_handler */ +} /* error_handler */ /********************************************************************* @@ -100,15 +111,15 @@ int nearest_index(double value, const double *array, int ia) } return index; -}; +} /******************************************************************/ void tokenize(const char * const string, const char *tokens, unsigned int varlen, unsigned int maxvar, char * pstring, unsigned int * const nstr) { - size_t i, j, nvar, len, ntoken; - int found, n; + size_t i, j, nvar, len, ntoken, n; + int found; nvar = 0; j = 0; len = strlen(string); @@ -160,7 +171,7 @@ double maxval_double(int size, const double *data) return maxval; -}; /* maxval_double */ +} /* maxval_double */ /******************************************************************************* @@ -179,7 +190,7 @@ double minval_double(int size, const double *data) return minval; -}; /* minval_double */ +} /* minval_double */ /******************************************************************************* double avgval_double(int size, double *data) @@ -196,7 +207,7 @@ double avgval_double(int size, const double *data) return avgval; -}; /* avgval_double */ +} /* avgval_double */ /******************************************************************************* @@ -223,7 +234,7 @@ void xyz2latlon( int np, const double *x, const double *y, const double *z, doub { double xx, yy, zz; - double dist, sinp; + double dist; int i; for(i=0; i M_PI) dx = dx - 2.0*M_PI; if(dx < -M_PI) dx = dx + 2.0*M_PI; return (dx*(sin(ur_lat)-sin(ll_lat))*RADIUS*RADIUS ) ; -}; /* box_area */ +} /* box_area */ +//TODO: Determine if poly_area_dimensionless can be deleted. +// Possibly functions that call it (e.g. get_grid_area_dimensionless) +// can also be deleted. /*------------------------------------------------------------------------------ double poly_area(const x[], const y[], int n) obtains area of input polygon by line integrating -sin(lat)d(lon) @@ -287,7 +300,7 @@ double poly_area_dimensionless(const double x[], const double y[], int n) if(dx < -M_PI) dx = dx + 2.0*M_PI; if (dx==0.0) continue; - if ( fabs(lat1-lat2) < SMALL_VALUE) /* cheap area calculation along latitude */ + if ( fabs(lat1-lat2) < SMALL_VALUE) // cheap area calculation along latitude area -= dx*sin(0.5*(lat1+lat2)); else { if(reproduce_siena) { @@ -301,14 +314,13 @@ double poly_area_dimensionless(const double x[], const double y[], int n) } } if(fabs(area) > HPI) { - printf("Error in poly_area_dimensionless: Large values for poly_area_dimensionless: %19.15f\n", area); - } + printf("Error in poly_area_dimensionless: Large values for poly_area_dimensionless: %19.15f\n", area); + } if(area < 0) return (-area/(4*M_PI)); else return (area/(4*M_PI)); - -}; /* poly_area */ +} /*------------------------------------------------------------------------------ double poly_area(const x[], const y[], int n) @@ -321,25 +333,25 @@ double poly_area_dimensionless(const double x[], const double y[], int n) Spherical Coordinates, P. Jones, Monthly Weather Review, 1998, vol127, p2204 The following is an implementation of equation (12) in the above paper: \int dA = \int_c [-sin(lat)] dlon - + An alternative derivation of line integrating -sin(lat)d(lon) formula: - Consider a vector function in spherical coordinates (r,lon,lat) + Consider a vector function in spherical coordinates (r,lon,lat) with only a lon component : - A=(0, (1-sin(lat))/cos(lat)/r , 0) - Then - Curl(A)=(1/r**2 , 0, 0) . - Now consider any loop C on the suface of the sphere enclosing an area S - and apply the Stokes theorem: - \integral_surface_S Curl(A).da = \integral_loop_C A.dl + A=(0, (1-sin(lat))/cos(lat)/r , 0) + Then + Curl(A)=(1/r**2 , 0, 0) . + Now consider any loop C on the suface of the sphere enclosing an area S + and apply the Stokes theorem: + \integral_surface_S Curl(A).da = \integral_loop_C A.dl where da and dl are the vectorial suface and line elements on sphere: da=(da, 0, 0) and dl=(0, R*cos(lat)*d(lon), R*d(lat)). We get \integral_surface_S Curl(A) = \integral_surface_S (da/r**2) = S/R**2 and - \integral_loop_C A.dl = \integral_loop_C (1-sin(lat))*d(lon) + \integral_loop_C A.dl = \integral_loop_C (1-sin(lat))*d(lon) Hence per Stokes formula: - S/R**2 = \integral_loop_C d(lon) - \integral_loop_C sin(lat)*d(lon). + S/R**2 = \integral_loop_C d(lon) - \integral_loop_C sin(lat)*d(lon). = I1 - I2 Now the approximation used for the second loop integral @@ -350,7 +362,7 @@ double poly_area_dimensionless(const double x[], const double y[], int n) If d(lon)/d(lat) is assumed constant over a loop segment, given that the segments used here are the sides of a grid cell, then we can take it outside the integral sum_over_loop_segments \integral_loop_C sin(lat)*d(lat) *d(lon)/d(lat) - =sum_over_loop_segments delta(lon)/delta(lat) \int_segment sin(lat)*d(lat) + =sum_over_loop_segments delta(lon)/delta(lat) \int_segment sin(lat)*d(lat) =sum_over_grid_cell_sides (lon2-lon1)/(lat2-lat1) * (-(cos(lat2)-cos(lat1))) =sum_over_grid_cell_sides (lon2-lon1)/(lat2-lat1) * (2*sin(0.5*(lat2+lat1))*sin(0.5*(lat2-lat1))) @@ -376,7 +388,7 @@ double poly_area_dimensionless(const double x[], const double y[], int n) as the next vertex in the cell (which is 90 degrees or close to 90 degrees appart) the contribution to I1 shifts into -I2 and we can again assume I1=0. This scheme is implemented in subroutine fix_lon() which is always called before poly_area(). That is why in the implementation below I1 is totally absent. - + In some pathological grid cells that do not have a pole vertex, I1 may not come out zero due to ambiguities in choosing d(lon) mod 2pi and we have to be careful in the implementation to deal with those grid cells. As we saw the approximation for I2 is based on the assumtion that d(lon)/d(lat) is @@ -385,12 +397,12 @@ double poly_area_dimensionless(const double x[], const double y[], int n) i.e., lat = lat1 + (lon-lon1)*(lat2-lat1)/(lon2-lon1) This assumption breaks down if the segment passes through a pole, e.g. along a longitude great circle where lon abruptly changes by pi as it goes through the pole like in the following two grid cells: - x----x----x + x----x----x | | | | Pole | | | | | | | - x----x----x + x----x----x x denotes the grid points. We can diagnose such situations either via I1=sum(dlons) coming out as 2*pi instead of 0, in that case we can correct the area by 2*pi @@ -412,66 +424,125 @@ double poly_area_dimensionless(const double x[], const double y[], int n) So, I2=\integral dx sin(arctan(tan(lat0)*cos(x-lon0))) can be shown to give an accurate estimate of grid cell areas surrounding the pole without a bump. ----------------------------------------------------------------------------*/ -double poly_area(const double x[], const double y[], int n) -{ +double poly_area_main(const double x[], const double y[], int n) { double area = 0.0; - int i; + int i; - for (i=0;i M_PI) dx = dx - 2.0*M_PI; - if(dx < -M_PI) dx = dx + 2.0*M_PI; + if (dx > M_PI) + dx = dx - 2.0 * M_PI; + if (dx < -M_PI) + dx = dx + 2.0 * M_PI; /*Sides that go through a pole contribute PI regardless of their direction and extent.*/ - if(fabs(dx+M_PI)< SMALL_VALUE || fabs(dx-M_PI)< SMALL_VALUE){ + if (fabs(dx + M_PI) < SMALL_VALUE || fabs(dx - M_PI) < SMALL_VALUE) { area += M_PI; - continue; //next i + continue; // next i } /* We have to be careful in implementing sin(0.5*(lat2-lat1))/(0.5*(lat2-lat1)) in the limit of lat1=lat2 where it becomes 1. Note that in that limit we get the well known formula for the area of a section of sphere bounded by two longitudes and two latitude circles. */ - if ( fabs(lat1-lat2) < SMALL_VALUE) /* cheap area calculation along latitude */ - area -= dx*sin(0.5*(lat1+lat2)); + if (fabs(lat1 - lat2) < SMALL_VALUE) /* cheap area calculation along latitude */ + area -= dx * sin(0.5 * (lat1 + lat2)); else { - if(reproduce_siena) { - area += dx*(cos(lat1)-cos(lat2))/(lat1-lat2); - } - else { - //This expression is a trig identity with the above reproduce_siena case - dy = 0.5*(lat1-lat2); - dat = sin(dy)/dy; - area -= dx*sin(0.5*(lat1+lat2))*dat; + if (reproduce_siena) { + area += dx * (cos(lat1) - cos(lat2)) / (lat1 - lat2); + } else { + // This expression is a trig identity with the above reproduce_siena case + dy = 0.5 * (lat1 - lat2); + dat = sin(dy) / dy; + area -= dx * sin(0.5 * (lat1 + lat2)) * dat; } } } - if(fabs(area) > HPI) printf("WARNING poly_area: Large values for area: %19.15f\n", area); - if(area < 0) - return -area*RADIUS*RADIUS; + + + if (fabs(area) > HPI) + printf("WARNING poly_area: Large values for area: %19.15f\n", area); + if (area < 0) + return -area * RADIUS * RADIUS; else - return area*RADIUS*RADIUS; + return area * RADIUS * RADIUS; +} /* poly_area_main */ + +/* + Calculate the area of a polygon with the original poly_area function, + exept that for those polygons that are polar, if the global poly_are_flag + is set, the area is calculated by first rotating a copy of the polygon away + from the pole and calculating the area of the rotated polygon instead. + Polar means having any of its verices cloase to a pole, or an edge crossing + a pole. + + This routine is here merely for improving the accuracy of the are calculation + in the given conditions. + TODO: The tiling error reported by make_coupler mosaic may be non-zero when + using this feature. This may be developped in the future. +*/ +double poly_area(const double xo[], const double yo[], int n) { + double area_pa = 0.0; + double area_par = 0.0; + int pole = 0; + int crosses = 0; + + double xr[8]; // rotated lon + double yr[8]; // rotated lat + + if (rotate_poly_flag == 0) { + area_pa = poly_area_main(xo, yo, n); + return area_pa; + } else { + // anything near enough to the pole gets rotated + pole = is_near_pole(yo, n); + crosses = crosses_pole(xo, n); + if (crosses == 1 && pole == 0) { + error_handler("crosses == 1 && pole == 0"); + } + + if (pole == 1) { + if (n > 8) { + error_handler("poly_area: n > 8. n=%d,n"); + } + rotate_poly(xo, yo, n, xr, yr); -}; /* poly_area */ + int pole2 = is_near_pole(yr, n); + if (pole2 == 1) { + error_handler("poly_area: pole2 == 1"); + } + area_par = poly_area_main(xr, yr, n); + } else { + area_pa = poly_area_main(xo, yo, n); + } + + if (pole == 1) { + return area_par; + } else { + return area_pa; + } + } +} /*An alternate implementation of poly_area for future developments. Under construction.*/ +//TODO: double poly_area2(const double x[], const double y[], int n) { double area = 0.0; double dx,dy,dat,lat1,lat2,avg_y,hdy,da,dxs= 0.0; - int i,j,ip; - int hasPole=0, hasBadxm=0, hasBadxp=0; + int i, ip; + int hasBadxm=0, hasBadxp=0; for (i=0;i=HPI-TOLORENCE) { int im=(i+nn-1)%nn, ip=(i+1)%nn; @@ -646,7 +717,8 @@ int fix_lon(double x[], double y[], int n, double tlon) /*If a polygon side passes through a Pole insert twin vertices at the Pole*/ /*A fix is also directly applied to poly_area to handle this case.*/ for (i=0;i M_PI_2 ) { + pole = 1; + break; + } + } + return pole; +} + +/* + crosses_pole() returns 1 iff one line segment of the polygon has its end at opposit + sides of a pole. i.e. if the longitudes are seperated by about Pi. Note, for realistic + data (not huge polygons), if crosses_pole() reutrns 1, so should is_near_pole(). +*/ +int crosses_pole(const double x[] , int n) { + int has_cl = 0; + for (int i = 0; i < n; i++) { + int im = (i + n - 1) % n; + //int ip = (i + 1) % n; + double dx = x[i] - x[im]; + if (fabs(dx + M_PI) < SMALL_VALUE || fabs(dx - M_PI) < SMALL_VALUE) { + has_cl = 1; + break; + } + } + return has_cl; +} + +/* + Set the_rotation_matrix global variable. + The rotation is 45 degrees and about the vector with orign at + earths center and the direction <0,1,1>/SQRT(2). I.e. a big + rotation away from the pole if what is being rotaed is near a pole. + For rotation matricies formulas and examples, see F.S.Hill, Computer + Graphics Using OpenGL, @nd ed., Chapter 5.3. +*/ +void set_the_rotation_matrix() { + static const double is2 = 1.0 /M_SQRT2; + + static const double m00 = 0; + static const double m01 = - is2; + static const double m02 = is2; + static const double m11 = 1.0/2; + static const double m12 = 0.5; + + static const double m[3][3] = { {m00, m01, m02}, {m02, m11, m12},{m01, m12, m11} }; + + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + the_rotation_matrix[i][j] = m[i][j]; + } + } +} + +/* Rotate point given the passed in rotation matrix */ +void rotate_point(double rv[], double rmat [][3]) { + double v[3]; + + for (int i = 0; i < 3; i++) { + v[i] = 0.0; + for (int j = 0; j < 3; j++) { + v[i] += rmat[i][j] * rv[j]; + } + } + for (int i = 0; i < 3; i++) { + rv[i] = v[i]; + } +} + +/* + Rotate polygon defined by x[], y[] points and store in xr[], yr[]*/ +void rotate_poly(const double x[], const double y[], const int n, double xr[], double yr[]) { + double sv[2]; //a rotated lat/lon + double rv[3]; //rotated xyz point + for (int i = 0; i < n; i++) { + latlon2xyz(1, &x[i], &y[i], &rv[0], &rv[1], &rv[2]); + rotate_point(rv, the_rotation_matrix); + xyz2latlon(1, &rv[0], &rv[1], &rv[2], &sv[0], &sv[1]); + xr[i] = sv[0]; + yr[i] = sv[1]; + } +} diff --git a/tools/libfrencutils/mosaic_util.h b/tools/libfrencutils/mosaic_util.h index fbdb733d..a559f0f9 100644 --- a/tools/libfrencutils/mosaic_util.h +++ b/tools/libfrencutils/mosaic_util.h @@ -20,15 +20,12 @@ #ifndef _MOSAIC_UTIL_H #define _MOSAIC_UTIL_H 1 +#endif #ifndef RANGE_CHECK_CRITERIA #define RANGE_CHECK_CRITERIA 0.05 #endif -#ifndef _MATH_H -#include -#endif - /* override the `fabs` function based on the type */ #define fabs(X) _Generic((X), \ long double: fabsl, \ @@ -113,5 +110,12 @@ void setInbound(struct Node *interList, struct Node *list); int isInside(struct Node *node); void set_reproduce_siena_true(void); +void set_rotate_poly_true(void); +int is_near_pole(const double y[], int n); +int crosses_pole(const double x[], int n); +void rotate_point( double rv[], double rmat [][3]); +void rotate_poly(const double x[], const double y[], const int n, + double xr[], double yr[]); +void set_the_rotation_matrix(); + -#endif diff --git a/tools/make_coupler_mosaic/make_coupler_mosaic.c b/tools/make_coupler_mosaic/make_coupler_mosaic.c index c56f1d64..68e22e4d 100644 --- a/tools/make_coupler_mosaic/make_coupler_mosaic.c +++ b/tools/make_coupler_mosaic/make_coupler_mosaic.c @@ -125,11 +125,14 @@ char *usage[] = { " The default value is 1.e-6 ", " ", "--check check the tiling error", - "", + " ", "--print_memory debug memory usage when it is set ", - " ", + " ", "--reproduce_siena Set to reproduce siena shared codes results ", " ", + "--rotate_poly Set to calculate polar polygon areas by caculating the area of a copy ", + " of the polygon, with the copy being rotated far away from the pole. ", + " ", "--verbose Set --verbose to print out messages during running. ", "", @@ -405,6 +408,7 @@ int main (int argc, char *argv[]) int tile_parent, is_parent, ie_parent, js_parent, je_parent; int print_memory=0; int reproduce_siena=0; + int rotate_poly=0; static struct option long_options[] = { {"atmos_mosaic", required_argument, NULL, 'a'}, @@ -420,6 +424,7 @@ int main (int argc, char *argv[]) {"verbose", no_argument, NULL, 'v'}, {"print_memory", no_argument, NULL, 'p'}, {"reproduce_siena", no_argument, NULL, 'q'}, + {"rotate_poly", no_argument, NULL, 'u'}, {NULL, 0, NULL, 0} }; @@ -470,6 +475,9 @@ int main (int argc, char *argv[]) case 'q': reproduce_siena = 1; break; + case 'u': + rotate_poly = 1; + break; case '?': errflg++; } @@ -494,6 +502,8 @@ int main (int argc, char *argv[]) if(reproduce_siena) set_reproduce_siena_true(); + if(rotate_poly) set_rotate_poly_true(); + /*mosaic_file can not have the same name as amosaic, lmosaic or omosaic, also the file name of amosaic, lmosaic, omosaic can not be "mosaic.nc" */ diff --git a/tools/make_hgrid/make_hgrid.c b/tools/make_hgrid/make_hgrid.c index aa55c914..df1d35b7 100644 --- a/tools/make_hgrid/make_hgrid.c +++ b/tools/make_hgrid/make_hgrid.c @@ -307,6 +307,10 @@ char *usage[] = { " --non_length_angle When specified, will not output length(dx,dy) and ", " angle (angle_dx, angle_dy) ", " ", + " --rotate_poly Set to calculate polar polygon areas by caculating ", + " the area of a copy of the polygon, with the copy ", + " being rotated far away from the pole. ", + " ", " --verbose Will print out running time message when this ", " option is set. Otherwise the run will be silent ", " when there is no error. ", @@ -566,6 +570,7 @@ int main(int argc, char* argv[]) char entry[MAXBOUNDS*STRINGLEN]; int n, errflg, c, i; int option_index; + int rotate_poly=0; static struct option long_options[] = { {"grid_type", required_argument, NULL, 'a'}, @@ -605,6 +610,7 @@ int main(int argc, char* argv[]) {"no_length_angle", no_argument, NULL, 'M'}, {"angular_midpoint", no_argument, NULL, 'N'}, {"help", no_argument, NULL, 'h'}, + {"rotate_poly", no_argument, NULL, 'n'}, {"verbose", no_argument, NULL, 'v'}, {0, 0, 0, 0}, @@ -744,6 +750,9 @@ int main(int argc, char* argv[]) case 'N': use_angular_midpoint = 1; break; + case 'n': + rotate_poly = 1; + break; case 'v': verbose = 1; break; @@ -770,6 +779,8 @@ int main(int argc, char* argv[]) strcat(history, argv[i]); } + if(rotate_poly) set_rotate_poly_true(); + if(mpp_pe() == mpp_root_pe() && verbose) printf("==>NOTE: the grid type is %s\n",grid_type); if(strcmp(grid_type,"regular_lonlat_grid") ==0 ) diff --git a/tools/make_topog/make_topog.c b/tools/make_topog/make_topog.c index d209cce5..d6c361e4 100644 --- a/tools/make_topog/make_topog.c +++ b/tools/make_topog/make_topog.c @@ -237,6 +237,10 @@ char *usage[] = { " --min_thickness # inimum vertical thickness allowed. with default value ", " 0.1. Increase or decrease this number as needed. ", " ", + " --rotate_poly Set to calculate polar polygon areas by caculating the ", + " area of a copy of the polygon, with the copy being ", + " rotated far away from the pole. ", + " ", " --help Print out this message and then exit. ", " ", " --verbose Will print out running time message when this ", @@ -314,6 +318,7 @@ int main(int argc, char* argv[]) int cyclic_x, cyclic_y, tripolar_grid; int errflg = (argc == 1); int option_index, i, c; + int rotate_poly=0; unsigned int verbose = 0; /* @@ -366,6 +371,7 @@ int main(int argc, char* argv[]) {"fraction_full_cell", required_argument, NULL, 'R'}, {"dont_open_very_this_cell", no_argument, NULL, 'S'}, {"min_thickness", required_argument, NULL, 'T'}, + {"rotate_poly", no_argument, NULL, 'z'}, {"help", no_argument, NULL, 'h'}, {"verbose", no_argument, NULL, 'V'}, {0, 0, 0, 0}, @@ -514,6 +520,9 @@ int main(int argc, char* argv[]) case 'T': min_thickness = atof(optarg); break; + case 'z': + rotate_poly = 1; + break; case 'V': verbose = 1; break; @@ -530,6 +539,7 @@ int main(int argc, char* argv[]) } } + if(rotate_poly) set_rotate_poly_true(); /* Write out arguments value */ if(mpp_pe() == mpp_root_pe() && verbose) printf("NOTE from make_topog ==> the topog_type is: %s\n",topog_type); diff --git a/tools/ocean_model_grid_generator b/tools/ocean_model_grid_generator deleted file mode 160000 index 2b076607..00000000 --- a/tools/ocean_model_grid_generator +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 2b0766078502e2e87721dca79ee9dc5e376ed436