diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c3673bb1..7359a19e9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -203,7 +203,7 @@ list(APPEND fms_c_src_files grid_utils/gradient_c2l.c grid_utils/grid_utils.c grid_utils/tree_utils.c - horiz_interp/include/horiz_interp_conserve_xgrid.c + horiz_interp/create_xgrid_utils.c mpp/mpp_memuse.c parser/yaml_parser_binding.c parser/yaml_output_functions.c @@ -302,7 +302,8 @@ foreach(kind ${kinds}) add_library(${libTgt}_c OBJECT ${fms_c_src_files}) target_include_directories(${libTgt}_c PRIVATE include - grid_utils) + grid_utils + horiz_interp/include) target_compile_definitions(${libTgt}_c PRIVATE "${fms_defs}") target_link_libraries(${libTgt}_c PRIVATE NetCDF::NetCDF_C MPI::MPI_C) @@ -442,7 +443,8 @@ if(NOT kinds) add_library(${libTgt}_c OBJECT ${fms_c_src_files}) target_include_directories(${libTgt}_c PRIVATE include - grid_utils) + grid_utils + horiz_interp/include) target_compile_definitions(${libTgt}_c PRIVATE "${fms_defs}") target_link_libraries(${libTgt}_c PRIVATE NetCDF::NetCDF_C MPI::MPI_C) diff --git a/grid_utils/grid_utils.c b/grid_utils/grid_utils.c index b19279b09..e3d3284c0 100644 --- a/grid_utils/grid_utils.c +++ b/grid_utils/grid_utils.c @@ -476,7 +476,6 @@ int invert_matrix_3x3(long double m[], long double m_inv[]) { return 1; } - int inside_a_polygon(double *lon1, double *lat1, int *npts, double *lon2, double *lat2) { @@ -521,17 +520,6 @@ int inside_a_polygon(double *lon1, double *lat1, int *npts, double *lon2, double } -int inside_a_polygon_(double *lon1, double *lat1, int *npts, double *lon2, double *lat2) -{ - - int isinside; - - isinside = inside_a_polygon(lon1, lat1, npts, lon2, lat2); - - return isinside; - -} - double get_global_area(void) { double garea; @@ -540,14 +528,6 @@ double get_global_area(void) return garea; } -double get_global_area_(void) -{ - double garea; - garea = 4*M_PI*RADIUS*RADIUS; - - return garea; -} - double poly_area(const double x[], const double y[], int n) { double area = 0.0; @@ -718,12 +698,7 @@ double spherical_excess_area(const double* p_ll, const double* p_ul, void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, const double *area) return the grid area. *******************************************************************************/ -void get_grid_area_(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area) -{ - get_grid_area(nlon, nlat, lon, lat, area); -} - -void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area) +void get_grid_area(const int* nlon, const int* nlat, const double *lon, const double *lat, double *area) { int nx, ny, nxp, i, j, n_in; double x_in[20], y_in[20]; @@ -736,7 +711,7 @@ void get_grid_area(const int *nlon, const int *nlat, const double *lon, const do x_in[0] = lon[j*nxp+i]; x_in[1] = lon[j*nxp+i+1]; x_in[2] = lon[(j+1)*nxp+i+1]; - x_in[3] = lon[(j+1)*nxp+i]; + x_in[3] = lon[(j+1)*nxp+i]; // fails here y_in[0] = lat[j*nxp+i]; y_in[1] = lat[j*nxp+i+1]; y_in[2] = lat[(j+1)*nxp+i+1]; @@ -752,12 +727,7 @@ void get_grid_area(const int *nlon, const int *nlat, const double *lon, const do void get_grid_area_ug(const int *npts, const double *lon, const double *lat, const double *area) return the grid area. *******************************************************************************/ -void get_grid_area_ug_(const int *npts, const double *lon, const double *lat, double *area) -{ - get_grid_area_ug(npts, lon, lat, area); -} - -void get_grid_area_ug(const int *npts, const double *lon, const double *lat, double *area) +void get_grid_area_ug(const int* npts, const double *lon, const double *lat, double *area) { int nl, l, n_in, nv; double x_in[20], y_in[20]; @@ -780,21 +750,13 @@ void get_grid_area_ug(const int *npts, const double *lon, const double *lat, dou } /* get_grid_area_ug */ - -void get_grid_great_circle_area_(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area) -{ - get_grid_great_circle_area(nlon, nlat, lon, lat, area); - -} - -void get_grid_great_circle_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) { int nx, ny, nxp, nyp, i, j; int n0, n1, n2, n3; struct Node *grid=NULL; double *x=NULL, *y=NULL, *z=NULL; - nx = *nlon; ny = *nlat; nxp = nx + 1; @@ -827,13 +789,7 @@ void get_grid_great_circle_area(const int *nlon, const int *nlat, const double * } /* get_grid_great_circle_area */ -void get_grid_great_circle_area_ug_(const int *npts, const double *lon, const double *lat, double *area) -{ - get_grid_great_circle_area_ug(npts, lon, lat, area); - -} - -void get_grid_great_circle_area_ug(const int *npts, const double *lon, const double *lat, double *area) +void get_grid_great_circle_area_ug(const int* npts, const double *lon, const double *lat, double *area) { int l, nl, nv; int n0, n1, n2, n3; @@ -870,7 +826,7 @@ void get_grid_great_circle_area_ug(const int *npts, const double *lon, const dou } /* get_grid_great_circle_area_ug */ -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) { int nx, ny, nxp, i, j, n_in; double x_in[20], y_in[20]; @@ -896,7 +852,7 @@ void get_grid_area_dimensionless(const int *nlon, const int *nlat, const double -void get_grid_area_no_adjust(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 nx, ny, nxp, i, j, n_in; double x_in[20], y_in[20]; @@ -1734,3 +1690,5 @@ int inside_edge(double x0, double y0, double x1, double y1, double x, double y) return (product<=SMALL) ? 1:0; } /* inside_edge */ + + diff --git a/grid_utils/grid_utils.h b/grid_utils/grid_utils.h index 563044798..9be6c0e33 100644 --- a/grid_utils/grid_utils.h +++ b/grid_utils/grid_utils.h @@ -35,6 +35,15 @@ #define max(a,b) (a>b ? a:b) #define SMALL_VALUE ( 1.e-10 ) +// append underscore so we can call from fortran +#define get_maxxgrid get_maxxgrid_ +#define get_global_area get_global_area_ +#define get_grid_area get_grid_area_ +#define get_grid_great_circle_area get_grid_great_circle_area_ +#define get_grid_area_no_adjust get_grid_area_no_adjust_ +#define get_grid_area_ug get_grid_area_ug_ +#define inside_a_polygon inside_a_polygon_ + void error_handler(const char *msg); int lon_fix(double *x, double *y, int n_in, double tlon); @@ -84,8 +93,6 @@ int inside_a_polygon( double *lon1, double *lat1, int *npts, double *lon2, doubl int samePoint(double x1, double y1, double z1, double x2, double y2, double z2); -int inside_a_polygon_(double *lon1, double *lat1, int *npts, double *lon2, double *lat2); - int inside_edge(double x0, double y0, double x1, double y1, double x, double y); int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double *q3, @@ -97,12 +104,8 @@ double poly_ctrlat(const double lon[], const double lat[], int n); int get_maxxgrid(void); -int get_maxxgrid_(void); - double get_global_area(void); -double get_global_area_(void); - double poly_area(const double lon[], const double lat[], int n); double poly_area_dimensionless(const double x[], const double y[], int n); @@ -110,11 +113,11 @@ double poly_area_dimensionless(const double x[], const double y[], int n); double spherical_excess_area(const double* p_ll, const double* p_ul, const double* p_lr, const double* p_ur, double radius); -void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area); +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_great_circle_area(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); +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[]); @@ -127,16 +130,8 @@ int clip_2dx2d_great_circle(const double x1_in[], const double y1_in[], const do const double x2_in[], const double y2_in[], const double z2_in [], int n2_in, double x_out[], double y_out[], double z_out[]); -void get_grid_area_ug(const int *npts, const double *lon, const double *lat, double *area); - -void get_grid_great_circle_area_ug(const int *npts, const double *lon, const double *lat, double *area); - -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_ug_(const int *npts, const double *lon, const double *lat, double *area); +void get_grid_area_ug(const int* npts, const double *lon, const double *lat, double *area); -void get_grid_great_circle_area_ug_(const int *npts, const double *lon, const double *lat, double *area); +void get_grid_great_circle_area_ug(const int* npts, const double *lon, const double *lat, double *area); #endif diff --git a/horiz_interp/Makefile.am b/horiz_interp/Makefile.am index 3c5289e62..64227f8cf 100644 --- a/horiz_interp/Makefile.am +++ b/horiz_interp/Makefile.am @@ -37,6 +37,7 @@ libhoriz_interp_la_SOURCES = \ horiz_interp.F90 \ horiz_interp_spherical.F90 \ horiz_interp_type.F90 \ + create_xgrid_utils.c \ include/horiz_interp_bicubic.inc \ include/horiz_interp_bilinear.inc \ include/horiz_interp_conserve.inc \ @@ -44,8 +45,7 @@ libhoriz_interp_la_SOURCES = \ include/horiz_interp_spherical.inc \ include/horiz_interp_type.inc \ include/horiz_interp_bicubic_r4.fh \ - include/horiz_interp_conserve_xgrid.h \ - include/horiz_interp_conserve_xgrid.c \ + include/create_xgrid.h \ include/horiz_interp_bilinear_r4.fh \ include/horiz_interp_conserve_r4.fh \ include/horiz_interp_r4.fh \ diff --git a/horiz_interp/include/horiz_interp_conserve_xgrid.c b/horiz_interp/create_xgrid_utils.c similarity index 72% rename from horiz_interp/include/horiz_interp_conserve_xgrid.c rename to horiz_interp/create_xgrid_utils.c index 51c484478..64d142bf7 100644 --- a/horiz_interp/include/horiz_interp_conserve_xgrid.c +++ b/horiz_interp/create_xgrid_utils.c @@ -20,7 +20,7 @@ #include #include "grid_utils.h" #include "tree_utils.h" -#include "horiz_interp_conserve_xgrid.h" +#include "create_xgrid.h" #include "constant.h" #if defined(_OPENMP) @@ -30,7 +30,17 @@ /** \file * \ingroup mosaic * \brief Grid creation and calculation functions for use in @ref mosaic_mod - * / + * This file provides routines to be called via interfaces in create_xgrid_mod. + **/ + +/******************************************************************************* + int get_maxxgrid + return constants MAXXGRID. +*******************************************************************************/ +int get_maxxgrid(void) +{ + return MAXXGRID; +} /******************************************************************************* void create_xgrid_1dx2d_order1 @@ -38,19 +48,7 @@ conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds. *******************************************************************************/ -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, - const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) -{ - int nxgrid; - - nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, - i_in, j_in, i_out, j_out, xgrid_area); - return nxgrid; - -} - -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, +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, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) @@ -142,19 +140,7 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds. *******************************************************************************/ -int create_xgrid_1dx2d_order1_ug_(const int *nlon_in, const int *nlat_in, const int *npts_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 *l_out, double *xgrid_area) -{ - int nxgrid; - - nxgrid = create_xgrid_1dx2d_order1_ug(nlon_in, nlat_in, npts_out, lon_in, lat_in, lon_out, lat_out, mask_in, - i_in, j_in, l_out, xgrid_area); - return nxgrid; - -} - -int create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const int *npts_out, const double *lon_in, +int create_xgrid_1dx2d_order1_ug(const int* nlon_in, const int* nlat_in, const int* npts_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 *l_out, double *xgrid_area) { @@ -242,18 +228,7 @@ int create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const i conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds. ********************************************************************************/ -int create_xgrid_1dx2d_order2_(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) -{ - int nxgrid; - nxgrid = create_xgrid_1dx2d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, - j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); - return nxgrid; - -} -int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, +int create_xgrid_1dx2d_order2(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) @@ -342,19 +317,7 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. *******************************************************************************/ -int create_xgrid_2dx1d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area) -{ - int nxgrid; - - nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, - i_in, j_in, i_out, j_out, xgrid_area); - return nxgrid; - -} -int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, +int create_xgrid_2dx1d_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, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) @@ -444,19 +407,7 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. ********************************************************************************/ -int create_xgrid_2dx1d_order2_(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) -{ - int nxgrid; - nxgrid = create_xgrid_2dx1d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, - j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); - return nxgrid; - -} - -int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, +int create_xgrid_2dx1d_order2(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) @@ -543,231 +494,110 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int } /* create_xgrid_2dx1d_order2 */ -/******************************************************************************* - void create_xgrid_2DX2D_order1 - This routine generate exchange grids between two grids for the first order - conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell - and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds. - mask is on grid lon_in/lat_in. -*******************************************************************************/ -int create_xgrid_2dx2d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area) -{ - int nxgrid; - nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, - i_in, j_in, i_out, j_out, xgrid_area); - return nxgrid; -} -int create_xgrid_2dx2d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, +/** + void create_xgrid_2DX2D_order1 + This routine generate exchange grids between two grids for the first order + conservative interpolation. + + It takes the arguments nlon_in, nlat_in, nlon_out, and nlat_out for amount of grid cells in each direction + and lon_in, lat_in, lon_out, lat_out for the geographic grid locations of each cell's bounds. + Values can also be masked out via mask_in/mask_out. + + From this grid information, it generates an exchange grid (finest refinement over the input/output grids), + in the form of i/j indices +**/ +int create_xgrid_2dx2d_order1(const int* nlon_input_cells, const int* nlat_input_cells, + const int* nlon_output_cells, const int* nlat_output_cells, + const double *input_grid_lon, const double *input_grid_lat, + const double *output_grid_lon, const double *output_grid_lat, + const double *skip_input_cells, const double *skip_output_cells, + int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) { - int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid; +#define MAX_V 8 + int nx_input_cells, nx_output_cells, ny_input_cells, ny_output_cells, nx_input_points, nx_output_points, nxgrid; double *area_in, *area_out; - int nblocks =1; + int nthreads =1; int *istart2=NULL, *iend2=NULL; - int npts_left, nblks_left, pos, m, npts_my, ij; + int npts_left, nblks_left, pos, curr_thread, 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; 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; -#if defined(_OPENMP) -#pragma omp parallel - nthreads = omp_get_num_threads(); -#endif - - 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 MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V"); - lon_out_min_list[n] = minval_double(n2_in, x2_in); - lon_out_max_list[n] = maxval_double(n2_in, x2_in); - lon_out_avg[n] = avgval_double(n2_in, x2_in); - n2_list[n] = n2_in; - for(l=0; l 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(4, y1_in); - lat_in_max = maxval_double(4, y1_in); - n1_in = fix_lon(x1_in, y1_in, 4, M_PI); - lon_in_min = minval_double(n1_in, x1_in); - lon_in_max = maxval_double(n1_in, x1_in); - lon_in_avg = avgval_double(n1_in, x1_in); - for(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 M_PI) { - lon_out_min -= TPI; - lon_out_max -= TPI; - for (l=0; l= lon_in_max || lon_out_max <= lon_in_min ) continue; - if ( (n_out = clip_2dx2d( 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 (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; - min_area = 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("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; - pi_in[nn] = i1; - pj_in[nn] = j1; - pi_out[nn] = i2; - pj_out[nn] = j2; - } - - } - - } - } + // clip the edge of each input cell against output cells in the bounding box and calculate areas for the exchange grid + #if defined(_OPENMP) + #pragma omp parallel for default(none) shared(nx_input_cells, ny_input_cells, nx_input_points, \ + input_grid_lon, input_grid_lat, skip_input_cells, \ + nx_output_cells, lat_out_min_list, lat_out_max_list, \ + n2_list, lon_out_list, lat_out_list, lon_out_min_list, \ + lon_out_max_list, lon_out_avg, area_in, area_out, \ + pxgrid_area, pnxgrid, pi_in, pj_in, pi_out, pj_out, pstart, \ + istart2, iend2, skip_output_cells, nthreads) + #endif + for(curr_thread=0; curr_thread 1 */ - if(nblocks == 1) { + // if using openmp, copy data over to the original arrays + if(nthreads == 1) { nxgrid = pnxgrid[0]; pi_in = NULL; pj_in = NULL; @@ -778,9 +608,9 @@ int create_xgrid_2dx2d_order1(const int *nlon_in, const int *nlat_in, const int else { int nn, i; nxgrid = 0; - for(m=0; m MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V"); + lon_out_min_list[n] = minval_double(n2_in, x2_in); + lon_out_max_list[n] = maxval_double(n2_in, x2_in); + lon_out_avg[n] = avgval_double(n2_in, x2_in); + n2_list[n] = n2_in; + for(l=0; l 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*nx_input_points+i1; n1 = j1*nx_input_points+i1+1; + n2 = (j1+1)*nx_input_points+i1+1; n3 = (j1+1)*nx_input_points+i1; + x1_in[0] = input_grid_lon[n0]; y1_in[0] = input_grid_lat[n0]; + x1_in[1] = input_grid_lon[n1]; y1_in[1] = input_grid_lat[n1]; + x1_in[2] = input_grid_lon[n2]; y1_in[2] = input_grid_lat[n2]; + x1_in[3] = input_grid_lon[n3]; y1_in[3] = input_grid_lat[n3]; + lat_in_min = minval_double(4, y1_in); + lat_in_max = maxval_double(4, y1_in); + n1_in = fix_lon(x1_in, y1_in, 4, M_PI); + lon_in_min = minval_double(n1_in, x1_in); + lon_in_max = maxval_double(n1_in, x1_in); + lon_in_avg = avgval_double(n1_in, x1_in); + + for(ij=istart2[curr_thread]; ij<=iend2[curr_thread]; ij++) { + + if(skip_output_cells[ij]>MASK_THRESH) { + int n_in, 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%nx_output_cells; + j2 = ij/nx_output_cells; + + 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 M_PI) { + lon_out_min -= TPI; + lon_out_max -= TPI; + for (l=0; l= lon_in_max || lon_out_max <= lon_in_min ) continue; + if ( (n_out = clip_2dx2d( 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 (x_out, y_out, n_out ) * skip_input_cells[j1*nx_input_cells+i1]; + min_area = min(area_in[j1*nx_input_cells+i1], area_out[j2*nx_output_cells+i2]); + + if( xarea/min_area > AREA_RATIO_THRESH ) { + pnxgrid[curr_thread]++; + if(pnxgrid[curr_thread]>= MAXXGRID/nthreads) + error_handler("The xgrid size is too large for resources.\n" + " nxgrid is greater than MAXXGRID/nthreads; increase MAXXGRID,\n" + " decrease nthreads, or increase number of MPI ranks."); + nn = pstart[curr_thread] + pnxgrid[curr_thread]-1; + + pxgrid_area[nn] = xarea; + pi_in[nn] = i1; + pj_in[nn] = j1; + pi_out[nn] = i2; + pj_out[nn] = j2; + + } + + } + } + } + } } + diff --git a/horiz_interp/include/create_xgrid.h b/horiz_interp/include/create_xgrid.h new file mode 100644 index 000000000..9a9f90d6b --- /dev/null +++ b/horiz_interp/include/create_xgrid.h @@ -0,0 +1,119 @@ +/*********************************************************************** + * GNU Lesser General Public License + * + * This file is part of the GFDL Flexible Modeling System (FMS). + * + * FMS 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. + * + * FMS 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 FMS. If not, see . + **********************************************************************/ +#ifndef HORIZ_INTERP_CREATE_XGRID_H_ +#define HORIZ_INTERP_CREATE_XGRID_H_ + +#ifndef MAXXGRID +#define MAXXGRID 1e6 +#endif + +#define AREA_RATIO_THRESH (1.e-6) +#define MASK_THRESH (0.5) +#define MAX_V 8 + +// redefine functions for fortran usage (ie. append underscore) +// TODO proper fortran bind(C) routines +#define get_maxxgrid get_maxxgrid_ +#define create_xgrid_1dx2d_order1 create_xgrid_1dx2d_order1_ +#define create_xgrid_2dx1d_order1 create_xgrid_2dx1d_order1_ +#define create_xgrid_2dx2d_order1 create_xgrid_2dx2d_order1_ +#define create_xgrid_2dx2d_order2 create_xgrid_2dx2d_order2_ +#define create_xgrid_great_circle create_xgrid_great_circle_ + +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, + const double *mask_in, int *i_in, int *j_in, int *i_out, + int *j_out, double *xgrid_area); + +int create_xgrid_1dx2d_order2(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); + +int create_xgrid_2dx1d_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, + const double *mask_in, int *i_in, int *j_in, int *i_out, + int *j_out, double *xgrid_area); + +int create_xgrid_2dx1d_order2(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); + +int create_xgrid_2dx2d_order1(const int* nlon_input_cells, const int* nlat_input_cells, + const int* nlon_output_cells, const int* nlat_output_cells, + const double *input_grid_lon, const double *input_grid_lat, + const double *output_grid_lon, const double *output_grid_lat, + const double *skip_input_cells, const double *skip_output_cells, + int *i_in, int *j_in, int *i_out, + int *j_out, double *xgrid_area); + +int create_xgrid_2dx2d_order2(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); + +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, + 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 create_xgrid_1dx2d_order1_ug(const int* nlon_in, const int* nlat_in, const int* npts_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 *l_out, double *xgrid_area); + +int create_xgrid_great_circle_ug(const int* nlon_in, const int* nlat_in, const int* npts_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 *l_out, + double *xgrid_area, double *xgrid_clon, double *xgrid_clat); + +int get_maxxgrid(void); + +int block_setup(int ** i_in, int ** j_in, int ** i_out, int ** j_out, + double** xgrid_area, + int **pi_in, int **pj_in, int **pi_out, int **pj_out, + double **pxgrid_area, + int **istart2, int **iend2, int **pstart, int **pnxgrid, + const int nx2, const int ny2); + +void compute_output_cell_bounds(const int nx_output_cells, const int ny_output_cells, + const int nx_output_points, + const double *output_grid_lon, const double *output_grid_lat, + double *lon_out_min_list, double *lon_out_max_list, + double *lat_out_min_list, double *lat_out_max_list, + double *lon_out_avg, int *n2_list, + double *lon_out_list, double *lat_out_list); + + +void clip_and_calc_2dx2d_order1(int curr_thread, + int nx_input_cells, int ny_input_cells, int nx_input_points, + const double *input_grid_lon, const double *input_grid_lat, + const double *skip_input_cells, + int nx_output_cells, + const double *lat_out_min_list, const double *lat_out_max_list, + const int *n2_list, + const double *lon_out_list, const double *lat_out_list, + const double *lon_out_min_list, const double *lon_out_max_list, + const double *lon_out_avg, + const double *area_in, const double *area_out, + double *pxgrid_area, int *pnxgrid, int *pi_in, int *pj_in, int *pi_out, int *pj_out, int *pstart, + const int *istart2, const int *iend2, + const double *skip_output_cells, int nthreads); + +#endif diff --git a/horiz_interp/include/horiz_interp_conserve_xgrid.h b/horiz_interp/include/horiz_interp_conserve_xgrid.h deleted file mode 100644 index 4aeab07e8..000000000 --- a/horiz_interp/include/horiz_interp_conserve_xgrid.h +++ /dev/null @@ -1,114 +0,0 @@ -/*********************************************************************** - * Apache License 2.0 - * - * This file is part of the GFDL Flexible Modeling System (FMS). - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * FMS is distributed in the hope that it will be useful, but WITHOUT - * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied; - * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A - * PARTICULAR PURPOSE. See the License for the specific language - * governing permissions and limitations under the License. - ***********************************************************************/ -#ifndef HORIZ_INTERP_CREATE_XGRID_H_ -#define HORIZ_INTERP_CREATE_XGRID_H_ - -#ifndef MAXXGRID -#define MAXXGRID 1e6 -#endif - -#define AREA_RATIO_THRESH (1.e-6) -#define MASK_THRESH (0.5) -#define MAX_V 8 - -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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area); - -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, - const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area); - -int create_xgrid_1dx2d_order2(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); - -int create_xgrid_1dx2d_order2_(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); - -int create_xgrid_2dx1d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area); - -int create_xgrid_2dx1d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area); - -int create_xgrid_2dx1d_order2(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); - -int create_xgrid_2dx1d_order2_(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); - -int create_xgrid_2dx2d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area); - -int create_xgrid_2dx2d_order2(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); - -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, - 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 create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const int *npts_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 *l_out, double *xgrid_area); - -int create_xgrid_great_circle_ug(const int *nlon_in, const int *nlat_in, const int *npts_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 *l_out, - double *xgrid_area, double *xgrid_clon, double *xgrid_clat); - -int create_xgrid_2dx2d_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, - const double *mask_in, int *i_in, int *j_in, int *i_out, - int *j_out, double *xgrid_area); - -int create_xgrid_2dx2d_order2_(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); - -int create_xgrid_1dx2d_order1_ug_(const int *nlon_in, const int *nlat_in, const int *npts_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 *l_out, double *xgrid_area); - -int create_xgrid_great_circle_ug_(const int *nlon_in, const int *nlat_in, const int *npts_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 *l_out, - double *xgrid_area, double *xgrid_clon, double *xgrid_clat); - -int get_maxxgrid(void); -int get_maxxgrid_(void); - -#endif