From a3735ae82cb4c7843a53f2fd2bc28d8a3a836071 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Fri, 21 Nov 2025 17:04:33 -0500 Subject: [PATCH 1/7] add in refactored create_xgrid_2dx2d_order1 and create new module for bindings --- grid_utils/grid_utils.c | 42 +- grid_utils/grid_utils.h | 18 +- horiz_interp/Makefile.am | 6 +- horiz_interp/create_xgrid.F90 | 91 +++ ..._conserve_xgrid.c => create_xgrid_utils.c} | 711 +++++++++--------- horiz_interp/include/create_xgrid.h | 118 +++ .../include/horiz_interp_conserve_xgrid.h | 115 --- 7 files changed, 598 insertions(+), 503 deletions(-) create mode 100644 horiz_interp/create_xgrid.F90 rename horiz_interp/{include/horiz_interp_conserve_xgrid.c => create_xgrid_utils.c} (71%) create mode 100644 horiz_interp/include/create_xgrid.h delete mode 100644 horiz_interp/include/horiz_interp_conserve_xgrid.h diff --git a/grid_utils/grid_utils.c b/grid_utils/grid_utils.c index b2693054e5..2f489a550f 100644 --- a/grid_utils/grid_utils.c +++ b/grid_utils/grid_utils.c @@ -719,18 +719,18 @@ 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) +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]; - nx = *nlon; - ny = *nlat; + nx = nlon; + ny = nlat; nxp = nx + 1; for(j=0; j #include "grid_utils.h" #include "tree_utils.h" -#include "horiz_interp_conserve_xgrid.h" +#include "create_xgrid.h" #include "constant.h" #if defined(_OPENMP) @@ -31,7 +31,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 @@ -39,19 +49,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) @@ -63,10 +61,10 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = *nlon_in; - ny1 = *nlat_in; - nx2 = *nlon_out; - ny2 = *nlat_out; + nx1 = nlon_in; + ny1 = nlat_in; + nx2 = nlon_out; + ny2 = nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -143,19 +141,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) { @@ -166,10 +152,10 @@ int create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const i double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = *nlon_in; - ny1 = *nlat_in; + nx1 = nlon_in; + ny1 = nlat_in; nv = 4; - npts2 = *npts_out; + npts2 = npts_out; nxgrid = 0; nx1p = nx1 + 1; @@ -243,18 +229,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) @@ -266,10 +241,10 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = *nlon_in; - ny1 = *nlat_in; - nx2 = *nlon_out; - ny2 = *nlat_out; + nx1 = nlon_in; + ny1 = nlat_in; + nx2 = nlon_out; + ny2 = nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -343,19 +318,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) @@ -370,10 +333,10 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int double Xarea; - nx1 = *nlon_in; - ny1 = *nlat_in; - nx2 = *nlon_out; - ny2 = *nlat_out; + nx1 = nlon_in; + ny1 = nlat_in; + nx2 = nlon_out; + ny2 = nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -445,19 +408,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) @@ -473,10 +424,10 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int double xarea; - nx1 = *nlon_in; - ny1 = *nlat_in; - nx2 = *nlon_out; - ny2 = *nlat_out; + nx1 = nlon_in; + ny1 = nlat_in; + nx2 = nlon_out; + ny2 = nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -544,231 +495,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; @@ -779,9 +609,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 0000000000..d0e379b3f3 --- /dev/null +++ b/horiz_interp/include/create_xgrid.h @@ -0,0 +1,118 @@ +/*********************************************************************** + * 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_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 4711723357..0000000000 --- a/horiz_interp/include/horiz_interp_conserve_xgrid.h +++ /dev/null @@ -1,115 +0,0 @@ -/*********************************************************************** - * 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 - -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 From 38b2dee81107d97f80ab49a08e6aa8055b4ef8b7 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Fri, 21 Nov 2025 17:30:05 -0500 Subject: [PATCH 2/7] remove added module for now --- horiz_interp/Makefile.am | 1 - horiz_interp/create_xgrid.F90 | 91 ----------------------------------- 2 files changed, 92 deletions(-) delete mode 100644 horiz_interp/create_xgrid.F90 diff --git a/horiz_interp/Makefile.am b/horiz_interp/Makefile.am index 8cbf194ff0..9f92685e25 100644 --- a/horiz_interp/Makefile.am +++ b/horiz_interp/Makefile.am @@ -38,7 +38,6 @@ libhoriz_interp_la_SOURCES = \ horiz_interp_spherical.F90 \ horiz_interp_type.F90 \ create_xgrid_utils.c \ - create_xgrid.F90 \ include/horiz_interp_bicubic.inc \ include/horiz_interp_bilinear.inc \ include/horiz_interp_conserve.inc \ diff --git a/horiz_interp/create_xgrid.F90 b/horiz_interp/create_xgrid.F90 deleted file mode 100644 index 5b3cc81c08..0000000000 --- a/horiz_interp/create_xgrid.F90 +++ /dev/null @@ -1,91 +0,0 @@ -module create_xgrid_mod - - use, intrinsic :: iso_c_binding - - implicit none - - interface - - function 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) BIND(C, NAME="create_xgrid_1dx2d_order1") - import :: C_INT, C_DOUBLE - integer(C_INT), value :: nlon_in - integer(C_INT), value :: nlat_in - integer(C_INT), value :: nlon_out - integer(C_INT), value :: nlat_out - real(C_DOUBLE), dimension(*) :: lon_in - real(C_DOUBLE), dimension(*) :: lat_in - real(C_DOUBLE), dimension(*) :: lon_out - real(C_DOUBLE), dimension(*) :: lat_out - real(C_DOUBLE), dimension(*) :: mask_in - integer(C_INT), dimension(*) :: i_in - integer(C_INT), dimension(*) :: j_in - integer(C_INT), dimension(*) :: i_out - integer(C_INT), dimension(*) :: j_out - real(C_DOUBLE), dimension(*) :: xgrid_area - end function - - function 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) BIND(C, NAME="create_xgrid_2dx1d_order1") - import :: C_INT, C_DOUBLE - integer(C_INT), value :: nlon_in - integer(C_INT), value :: nlat_in - integer(C_INT), value :: nlon_out - integer(C_INT), value :: nlat_out - real(C_DOUBLE), dimension(*) :: lon_in - real(C_DOUBLE), dimension(*) :: lat_in - real(C_DOUBLE), dimension(*) :: lon_out - real(C_DOUBLE), dimension(*) :: lat_out - real(C_DOUBLE), dimension(*) :: mask_in - integer(C_INT), dimension(*) :: i_in - integer(C_INT), dimension(*) :: j_in - integer(C_INT), dimension(*) :: i_out - integer(C_INT), dimension(*) :: j_out - real(C_DOUBLE), dimension(*) :: xgrid_area - end function - - function 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) BIND(C, NAME="create_xgrid_2dx2d_order1") - import :: C_INT, C_DOUBLE - integer(C_INT), value :: nlon_in - integer(C_INT), value :: nlat_in - integer(C_INT), value :: nlon_out - integer(C_INT), value :: nlat_out - real(C_DOUBLE), dimension(*) :: lon_in - real(C_DOUBLE), dimension(*) :: lat_in - real(C_DOUBLE), dimension(*) :: lon_out - real(C_DOUBLE), dimension(*) :: lat_out - real(C_DOUBLE), dimension(*) :: mask_in - integer(C_INT), dimension(*) :: i_in - integer(C_INT), dimension(*) :: j_in - integer(C_INT), dimension(*) :: i_out - integer(C_INT), dimension(*) :: j_out - real(C_DOUBLE), dimension(*) :: xgrid_area - end function - - - function get_grid_area(nlon, nlat, lon, lat, area) BIND(C, NAME="get_grid_area") - import :: C_INT, C_DOUBLE - integer(C_INT), value :: nlon - integer(C_INT), value :: nlat - real(C_DOUBLE), dimension(*) :: lon - real(C_DOUBLE), dimension(*) :: lat - real(C_DOUBLE), dimension(*) :: area - end function - - end interface - -! the routines used in fms: -!../../libFMS/.libs/libFMS.so: undefined reference to `get_grid_area_' -!../../libFMS/.libs/libFMS.so: undefined reference to `create_xgrid_great_circle_' -!../../libFMS/.libs/libFMS.so: undefined reference to `create_xgrid_2dx2d_order1_' -!../../libFMS/.libs/libFMS.so: undefined reference to `create_xgrid_2dx1d_order1_' -!../../libFMS/.libs/libFMS.so: undefined reference to `get_grid_great_circle_area_' -!../../libFMS/.libs/libFMS.so: undefined reference to `get_maxxgrid_' -!../../libFMS/.libs/libFMS.so: undefined reference to `create_xgrid_1dx2d_order1_' - - -end module From 29bde6cbd5a496441b08226d56eb2d464c130329 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Wed, 3 Dec 2025 15:44:04 -0500 Subject: [PATCH 3/7] fix build issues from ci --- horiz_interp/Makefile.am | 1 - horiz_interp/include/create_xgrid.h | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/horiz_interp/Makefile.am b/horiz_interp/Makefile.am index 9f92685e25..64227f8cfa 100644 --- a/horiz_interp/Makefile.am +++ b/horiz_interp/Makefile.am @@ -74,7 +74,6 @@ MODFILES = \ horiz_interp_bilinear_mod.$(FC_MODEXT) \ horiz_interp_conserve_mod.$(FC_MODEXT) \ horiz_interp_spherical_mod.$(FC_MODEXT) \ - create_xgrid.$(FC_MODEXT) \ horiz_interp_mod.$(FC_MODEXT) nodist_include_HEADERS = $(MODFILES) BUILT_SOURCES = $(MODFILES) diff --git a/horiz_interp/include/create_xgrid.h b/horiz_interp/include/create_xgrid.h index d0e379b3f3..3cc182a4f5 100644 --- a/horiz_interp/include/create_xgrid.h +++ b/horiz_interp/include/create_xgrid.h @@ -33,6 +33,7 @@ #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, From 5c9460857a3766f3b8fc23d911efc3d490e99827 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Wed, 3 Dec 2025 16:06:15 -0500 Subject: [PATCH 4/7] fix cmake build --- CMakeLists.txt | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c3673bb10..7359a19e93 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) From be5f7b250697cccea69cb11bfb8c584e2cadfdf3 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Thu, 4 Dec 2025 14:15:33 -0500 Subject: [PATCH 5/7] change int args back to pointers for now --- grid_utils/grid_utils.c | 69 ++++++---------------- grid_utils/grid_utils.h | 30 +++++----- horiz_interp/create_xgrid_utils.c | 88 ++++++++++++++--------------- horiz_interp/include/create_xgrid.h | 20 +++---- 4 files changed, 84 insertions(+), 123 deletions(-) diff --git a/grid_utils/grid_utils.c b/grid_utils/grid_utils.c index 67efa62200..e5e7bd31ff 100644 --- a/grid_utils/grid_utils.c +++ b/grid_utils/grid_utils.c @@ -521,17 +521,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; @@ -718,25 +707,20 @@ 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]; - nx = nlon; - ny = nlat; + nx = *nlon; + ny = *nlat; nxp = nx + 1; for(j=0; jb ? 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_ + void error_handler(const char *msg); int lon_fix(double *x, double *y, int n_in, double tlon); @@ -97,12 +105,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 +114,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 +131,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/create_xgrid_utils.c b/horiz_interp/create_xgrid_utils.c index db6dff66eb..64d142bf7c 100644 --- a/horiz_interp/create_xgrid_utils.c +++ b/horiz_interp/create_xgrid_utils.c @@ -48,7 +48,7 @@ int get_maxxgrid(void) 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, +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) @@ -60,10 +60,10 @@ int create_xgrid_1dx2d_order1(const int nlon_in, const int nlat_in, const int nl double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -140,7 +140,7 @@ int create_xgrid_1dx2d_order1(const int nlon_in, const int nlat_in, const int nl 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, +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) { @@ -151,10 +151,10 @@ int create_xgrid_1dx2d_order1_ug(const int nlon_in, const int nlat_in, const int double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = nlon_in; - ny1 = nlat_in; + nx1 = *nlon_in; + ny1 = *nlat_in; nv = 4; - npts2 = npts_out; + npts2 = *npts_out; nxgrid = 0; nx1p = nx1 + 1; @@ -228,7 +228,7 @@ int create_xgrid_1dx2d_order1_ug(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_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) @@ -240,10 +240,10 @@ int create_xgrid_1dx2d_order2(const int nlon_in, const int nlat_in, const int nl double *area_in, *area_out, min_area; double *tmpx, *tmpy; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -317,7 +317,7 @@ int create_xgrid_1dx2d_order2(const int nlon_in, const int nlat_in, const int nl 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, +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) @@ -332,10 +332,10 @@ int create_xgrid_2dx1d_order1(const int nlon_in, const int nlat_in, const int nl double Xarea; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -407,7 +407,7 @@ int create_xgrid_2dx1d_order1(const int nlon_in, const int nlat_in, const int nl 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, +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) @@ -423,10 +423,10 @@ int create_xgrid_2dx1d_order2(const int nlon_in, const int nlat_in, const int nl double xarea; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; @@ -508,8 +508,8 @@ int create_xgrid_2dx1d_order2(const int nlon_in, const int nlat_in, const int nl 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, +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, @@ -531,10 +531,10 @@ int create_xgrid_2dx2d_order1(const int nlon_input_cells, const int nlat_input_c int *n2_list; int nxgrid_block_max; - nx_input_cells = nlon_input_cells; - ny_input_cells = nlat_input_cells; - nx_output_cells = nlon_output_cells; - ny_output_cells = nlat_output_cells; + nx_input_cells = *nlon_input_cells; + ny_input_cells = *nlat_input_cells; + nx_output_cells = *nlon_output_cells; + ny_output_cells = *nlat_output_cells; nx_input_points = nx_input_cells + 1; nx_output_points = nx_output_cells + 1; @@ -648,7 +648,7 @@ int create_xgrid_2dx2d_order1(const int nlon_input_cells, const int nlat_input_c 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_order2(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out, +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) @@ -667,10 +667,10 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl int *n2_list; int nthreads, nxgrid_block_max; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nx1p = nx1 + 1; nx2p = nx2 + 1; @@ -910,7 +910,7 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl /** * Variant of create_xgrid that uses the great circle algorithm to take account for curvature */ -int create_xgrid_great_circle(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_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, 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) @@ -926,10 +926,10 @@ int create_xgrid_great_circle(const int nlon_in, const int nlat_in, const int nl double *area1, *area2, min_area; - nx1 = nlon_in; - ny1 = nlat_in; - nx2 = nlon_out; - ny2 = nlat_out; + nx1 = *nlon_in; + ny1 = *nlat_in; + nx2 = *nlon_out; + ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; @@ -1012,7 +1012,7 @@ int create_xgrid_great_circle(const int nlon_in, const int nlat_in, const int nl * Variant of create_xgrid that uses the great circle algorithm, specifically for a unstructured grid * TODO: currently unused in FMS, should consider removal if not going to be tested */ -int create_xgrid_great_circle_ug(const int nlon_in, const int nlat_in, const int npts_out, +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) @@ -1028,10 +1028,10 @@ int create_xgrid_great_circle_ug(const int nlon_in, const int nlat_in, const int double *area1, *area2, min_area; - nx1 = nlon_in; - ny1 = nlat_in; + nx1 = *nlon_in; + ny1 = *nlat_in; nv = 4; - npts2 = npts_out; + npts2 = *npts_out; nxgrid = 0; nx1p = nx1 + 1; ny1p = ny1 + 1; diff --git a/horiz_interp/include/create_xgrid.h b/horiz_interp/include/create_xgrid.h index 3cc182a4f5..b990c595ee 100644 --- a/horiz_interp/include/create_xgrid.h +++ b/horiz_interp/include/create_xgrid.h @@ -36,49 +36,49 @@ #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, +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, +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, +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, +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, +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, +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, +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, +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, +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); From 029b18379a612e3dd08b388ade7d3caa09af8dda Mon Sep 17 00:00:00 2001 From: rem1776 Date: Thu, 4 Dec 2025 14:16:36 -0500 Subject: [PATCH 6/7] remove some unneeded wrapper functions --- grid_utils/grid_utils.c | 8 -------- 1 file changed, 8 deletions(-) diff --git a/grid_utils/grid_utils.c b/grid_utils/grid_utils.c index e5e7bd31ff..5740dfe42c 100644 --- a/grid_utils/grid_utils.c +++ b/grid_utils/grid_utils.c @@ -529,14 +529,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; From 40208a35780f57b7a41ec07234d25062d0231813 Mon Sep 17 00:00:00 2001 From: rem1776 Date: Thu, 4 Dec 2025 14:23:35 -0500 Subject: [PATCH 7/7] clean up spaces and add macro for inside_a_polygon --- grid_utils/grid_utils.c | 1 - grid_utils/grid_utils.h | 5 ++--- horiz_interp/include/create_xgrid.h | 2 +- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/grid_utils/grid_utils.c b/grid_utils/grid_utils.c index 5740dfe42c..e3d3284c01 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) { diff --git a/grid_utils/grid_utils.h b/grid_utils/grid_utils.h index 837b83424a..9be6c0e33f 100644 --- a/grid_utils/grid_utils.h +++ b/grid_utils/grid_utils.h @@ -41,7 +41,8 @@ #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 get_grid_area_ug get_grid_area_ug_ +#define inside_a_polygon inside_a_polygon_ void error_handler(const char *msg); @@ -92,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, diff --git a/horiz_interp/include/create_xgrid.h b/horiz_interp/include/create_xgrid.h index b990c595ee..9a9f90d6b9 100644 --- a/horiz_interp/include/create_xgrid.h +++ b/horiz_interp/include/create_xgrid.h @@ -31,7 +31,7 @@ // 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_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_