From 2a602341e1b5d3cd2dc0e9cd674d76bae26d677e Mon Sep 17 00:00:00 2001 From: mlee03 Date: Tue, 10 Oct 2023 13:08:46 -0400 Subject: [PATCH 01/12] introduce structs for out lists --- tools/fregrid/conserve_interp.c | 36 +- tools/fregrid/globals.h | 17 +- tools/libfrencutils/Makefile.am | 2 +- tools/libfrencutils/create_xgrid.c | 558 ++++++++++++------------ tools/libfrencutils/create_xgrid_acc.c | 32 +- tools/libfrencutils/create_xgrid_acc.h | 7 +- tools/libfrencutils/create_xgrid_util.c | 77 ++-- tools/libfrencutils/create_xgrid_util.h | 10 +- 8 files changed, 361 insertions(+), 378 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index ec3cb737..7b8924f0 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -84,6 +84,8 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles //START NTILES_OUT for(n=0; n #include #include +#include "globals.h" #include "mosaic_util.h" #include "create_xgrid.h" #include "create_xgrid_util.h" @@ -38,21 +39,21 @@ 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -102,9 +103,9 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_out[(j2+1)*nx2p+i2+1]; y_in[3] = lat_out[(j2+1)*nx2p+i2]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; @@ -113,17 +114,17 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = Xarea; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = Xarea; + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -143,9 +144,9 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int 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) + 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, @@ -154,9 +155,9 @@ int create_xgrid_1dx2d_order2_(const int *nlon_in, const int *nlat_in, const int }; 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -200,9 +201,9 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_out[(j2+1)*nx2p+i2+1]; y_in[3] = lat_out[(j2+1)*nx2p+i2]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; @@ -212,19 +213,19 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; + 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + if(xarea/min_area > AREA_RATIO_THRESH ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -243,21 +244,21 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -301,9 +302,9 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_in[(j1+1)*nx1p+i1+1]; y_in[3] = lat_in[(j1+1)*nx1p+i1]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_in[j1*nx1p+i1]; x_in[1] = lon_in[j1*nx1p+i1+1]; @@ -313,17 +314,17 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = Xarea; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = Xarea; + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -344,9 +345,9 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int 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) + 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, @@ -356,9 +357,9 @@ int create_xgrid_2dx1d_order2_(const int *nlon_in, const int *nlat_in, const int }; 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -404,9 +405,9 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_in[(j1+1)*nx1p+i1+1]; y_in[3] = lat_in[(j1+1)*nx1p+i1]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_in[j1*nx1p+i1]; x_in[1] = lon_in[j1*nx1p+i1+1]; @@ -417,19 +418,19 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -450,22 +451,22 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *******************************************************************************/ #ifndef __AIX 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) + 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); + i_in, j_in, i_out, j_out, xgrid_area); return nxgrid; }; #endif 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) + 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) { #define MAX_V 8 @@ -611,57 +612,57 @@ nxgrid = 0; 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= 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) { + 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]++; + 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; - } + 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; + } - } + } } } @@ -681,13 +682,13 @@ nxgrid = 0; nxgrid = 0; for(m=0; m= 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]; + 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]; #pragma acc loop seq - for(l=0; l M_PI) { - lon_out_min -= TPI; - lon_out_max -= TPI; + lon_out_min -= TPI; + lon_out_max -= TPI; #pragma acc loop seq - for (l=0; l= lon_in_max || lon_out_max <= lon_in_min ) continue; - n_out = 1; - if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) { + for (l=0; l= lon_in_max || lon_out_max <= lon_in_min ) continue; + n_out = 1; + 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; - 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; + 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 ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; #pragma atomic update - nxgrid++; - } - } + nxgrid++; + } + } } } } @@ -903,9 +904,9 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int };/* get_xgrid_2Dx2D_order2 */ #else 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) + 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) { #define MAX_V 8 @@ -1053,57 +1054,57 @@ nxgrid = 0; 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_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%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= 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) { + 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]++; + 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; - pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out ); - pi_in[nn] = i1; - pj_in[nn] = j1; - pi_out[nn] = i2; - pj_out[nn] = j2; - } - } + 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; + pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out ); + pi_in[nn] = i1; + pj_in[nn] = j1; + pi_out[nn] = i2; + pj_out[nn] = j2; + } + } } } } @@ -1124,15 +1125,15 @@ nxgrid = 0; nxgrid = 0; for(m=0; m 0) { - xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; - min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); - if( xarea/min_area > AREA_RATIO_THRESH ) { + x_out, y_out, z_out)) > 0) { + xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; + min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); + if( xarea/min_area > AREA_RATIO_THRESH ) { #ifdef debug_test_create_xgrid - printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea); + printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea); #endif - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */ - xgrid_clat[nxgrid] = 0; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */ + xgrid_clat[nxgrid] = 0; + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -1274,4 +1275,3 @@ int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int return nxgrid; };/* create_xgrid_great_circle */ - diff --git a/tools/libfrencutils/create_xgrid_acc.c b/tools/libfrencutils/create_xgrid_acc.c index b1903369..b4a3bd1e 100644 --- a/tools/libfrencutils/create_xgrid_acc.c +++ b/tools/libfrencutils/create_xgrid_acc.c @@ -20,6 +20,7 @@ #include #include #include +#include "globals.h" #include "mosaic_util.h" #include "create_xgrid.h" #include "create_xgrid_util.h" @@ -37,10 +38,7 @@ *******************************************************************************/ int create_xgrid_2dx2d_order2_acc(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 *lon_out_min_list, const double *lon_out_max_list, const double *lon_out_avg, - 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 *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, + Minmaxavglists *out_minmaxavglists, 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) { @@ -68,16 +66,16 @@ int create_xgrid_2dx2d_order2_acc(const int *nlon_in, const int *nlat_in, const #pragma acc data present(lon_out[0:(nx2+1)*(ny2+1)], lat_out[0:(nx2+1)*(ny2+1)]) #pragma acc data present(lon_in[0:(nx1+1)*(ny1+1)], lat_in[0:(nx1+1)*(ny1+1)], mask_in[0:nx1*ny1]) -#pragma acc data present(lon_out_list[0:MAX_V*nx2*ny2], lat_out_list[0:MAX_V*nx2*ny2], \ - lat_out_min_list[0:nx2*ny2], lat_out_max_list[0:nx2*ny2], \ - lon_out_min_list[0:nx2*ny2], lon_out_max_list[0:nx2*ny2], \ - lon_out_avg[0:nx2*ny2], n2_list[0:nx2*ny2]) +#pragma acc data present(out_minmaxavglists->lon_list[0:MAX_V*nx2*ny2], out_minmaxavglists->lat_list[0:MAX_V*nx2*ny2]) +#pragma acc data present(out_minmaxavglists->n_list[0:nx2*ny2], out_minmaxavglists->lon_avg[0:nx2*ny2]) +#pragma acc data present(out_minmaxavglists->lat_min_list[0:nx2*ny2], out_minmaxavglists->lat_max_list[0:nx2*ny2]) +#pragma acc data present(out_minmaxavglists->lon_min_list[0:nx2*ny2], out_minmaxavglists->lon_max_list[0:nx2*ny2]) #pragma acc data present(xgrid_area[0:mxxgrid], xgrid_clon[0:mxxgrid], xgrid_clat[0:mxxgrid], \ i_in[0:mxxgrid], j_in[0:mxxgrid], i_out[0:mxxgrid],j_out[0:mxxgrid]) #pragma acc data copyin(area_in[0:nx1*ny1], area_out[0:nx2*ny2]) #pragma acc kernels copy(nxgrid) { -#pragma acc loop independent collapse(2) //reduction(+:nxgrid) +#pragma acc loop independent collapse(2) reduction(+:nxgrid) for(j1=0; j1 MASK_THRESH ) { int n0, n1, n2, n3, n1_in; double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg; @@ -94,7 +92,7 @@ int create_xgrid_2dx2d_order2_acc(const int *nlon_in, const int *nlat_in, const 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); -#pragma acc loop independent //reduction(+:nxgrid) +#pragma acc loop independent reduction(+:nxgrid) for(ij=0; ij= lat_in_max || lat_out_max_list[ij] <= lat_in_min ) continue; + if(out_minmaxavglists->lat_min_list[ij] >= lat_in_max || out_minmaxavglists->lat_max_list[ij] <= lat_in_min ) continue; /* adjust x2_in according to lon_in_avg*/ - n2_in = n2_list[ij]; + n2_in = out_minmaxavglists->n_list[ij]; #pragma acc loop seq for(l=0; llon_list[ij*MAX_V+l]; + y2_in[l] = out_minmaxavglists->lat_list[ij*MAX_V+l]; } - lon_out_min = lon_out_min_list[ij]; - lon_out_max = lon_out_max_list[ij]; - dx = lon_out_avg[ij] - lon_in_avg; + lon_out_min = out_minmaxavglists->lon_min_list[ij]; + lon_out_max = out_minmaxavglists->lon_max_list[ij]; + dx = out_minmaxavglists->lon_avg[ij] - lon_in_avg; if(dx < -M_PI ) { lon_out_min += TPI; lon_out_max += TPI; diff --git a/tools/libfrencutils/create_xgrid_acc.h b/tools/libfrencutils/create_xgrid_acc.h index 823faf18..ac679107 100644 --- a/tools/libfrencutils/create_xgrid_acc.h +++ b/tools/libfrencutils/create_xgrid_acc.h @@ -20,11 +20,10 @@ #ifndef CREATE_XGRID_VER_H_ #define CREATE_XGRID_VER_H_ +#include "globals.h" + int create_xgrid_2dx2d_order2_acc(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 *lon_out_min_list, const double *lon_out_max_list, const double *lon_out_avg, - 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 *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, + Minmaxavglists *out_minmaxavglists, 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); #endif diff --git a/tools/libfrencutils/create_xgrid_util.c b/tools/libfrencutils/create_xgrid_util.c index 03224195..11e283e1 100644 --- a/tools/libfrencutils/create_xgrid_util.c +++ b/tools/libfrencutils/create_xgrid_util.c @@ -20,6 +20,7 @@ #include #include #include +#include "globals.h" #include "mosaic_util.h" #include "create_xgrid_util.h" #include "constant.h" @@ -43,47 +44,22 @@ int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double void malloc_minmaxavg_lists allocates lists to hold min, max, avg values of lat/lon coordinates *******************************************************************************/ -void malloc_minmaxavg_lists(const int n, - double **lon_min_list, double **lon_max_list, double **lat_min_list, double **lat_max_list, - int **n_list, double **lon_avg, double **lon_list, double **lat_list) +Minmaxavglists malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) { - if(*lon_min_list!=NULL){ - free(*lon_min_list) ; *lon_min_list=NULL; - } - if(*lon_max_list!=NULL){ - free(*lon_max_list) ; *lon_max_list=NULL; - } - if(*lat_min_list!=NULL){ - free(*lat_min_list) ; *lat_min_list=NULL; - } - if(*lat_max_list!=NULL){ - free(*lat_max_list) ; *lat_max_list=NULL; - } - if(*n_list!=NULL){ - free(*n_list) ; *n_list=NULL; - } - if(*lon_avg!=NULL){ - free(*lon_avg) ; *lon_avg=NULL; - } - if(*lon_list!=NULL){ - free(*lon_list) ; *lon_list=NULL; - } - if(*lat_list!=NULL){ - free(*lat_list) ; *lat_list=NULL; - } - if(n>0){ - *lon_min_list=(double *)malloc(n*sizeof(double)); - *lon_max_list=(double *)malloc(n*sizeof(double)); - *lat_min_list=(double *)malloc(n*sizeof(double)); - *lat_max_list=(double *)malloc(n*sizeof(double)); - *n_list=(int *)malloc(n*sizeof(int)); - *lon_avg=(double *)malloc(n*sizeof(double)); - *lon_list=(double *)malloc(MAX_V*n*sizeof(double)); - *lat_list=(double *)malloc(MAX_V*n*sizeof(double)); + minmaxavglists->lon_min_list=(double *)malloc(n*sizeof(double)); + minmaxavglists->lon_max_list=(double *)malloc(n*sizeof(double)); + minmaxavglists->lat_min_list=(double *)malloc(n*sizeof(double)); + minmaxavglists->lat_max_list=(double *)malloc(n*sizeof(double)); + minmaxavglists->n_list=(int *)malloc(n*sizeof(int)); + minmaxavglists->lon_avg=(double *)malloc(n*sizeof(double)); + minmaxavglists->lon_list=(double *)malloc(MAX_V*n*sizeof(double)); + minmaxavglists->lat_list=(double *)malloc(MAX_V*n*sizeof(double)); } + return *minmaxavglists; + }//malloc_minmaxavg_lists /******************************************************************************* @@ -91,8 +67,7 @@ void malloc_minmaxavg_lists(const int n, computes lists to hold min, max, avg values of lat/lon coordinates *******************************************************************************/ void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, - double *lon_min_list, double *lon_max_list, double *lat_min_list, double *lat_max_list, - int *n_list, double *lon_avg, double *lon_list, double *lat_list) + Minmaxavglists *minmaxavglists) { int nxp, nyp; @@ -100,9 +75,12 @@ void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const do nxp = nx+1; nyp = ny+1; -#pragma acc data present(lon[0:nxp*nyp], lat[0:nxp*nyp], lon_min_list[0:nx*ny], lon_max_list[0:nx*ny],\ - lat_min_list[0:nx*ny], lat_max_list[0:nx*ny], n_list[0:nx*ny],\ - lon_avg[0:nx*ny], lon_list[0:MAX_V*nx*ny], lat_list[0:MAX_V*nx*ny]) +#pragma acc data present(lon[0:nxp*nyp], lat[0:nxp*nyp]) +#pragma acc data present(minmaxavglists->lon_list[0:MAX_V*nx*ny], minmaxavglists->lat_list[0:MAX_V*nx*ny]) +#pragma acc data present(minmaxavglists->n_list[0:nx*ny], minmaxavglists->lon_avg[0:nx*ny]) +#pragma acc data present(minmaxavglists->lat_min_list[0:nx*ny], minmaxavglists->lat_max_list[0:nx*ny]) +#pragma acc data present(minmaxavglists->lon_min_list[0:nx*ny], minmaxavglists->lon_max_list[0:nx*ny]) + #pragma acc parallel loop independent for(int ij=0; ijlat_min_list[n] = minval_double(4, y_in); + minmaxavglists->lat_max_list[n] = maxval_double(4, y_in); n_in = fix_lon(x_in, y_in, 4, M_PI); + //Commented out for now. OpenACC does not like error_handler //if(n2_in > MAX_V) error_handler("create_xgrid.c: n_in is greater than MAX_V"); - lon_min_list[n] = minval_double(n_in, x_in); - lon_max_list[n] = maxval_double(n_in, x_in); - lon_avg[n] = avgval_double(n_in, x_in); - n_list[n] = n_in; + minmaxavglists->lon_min_list[n] = minval_double(n_in, x_in); + minmaxavglists->lon_max_list[n] = maxval_double(n_in, x_in); + minmaxavglists->lon_avg[n] = avgval_double(n_in, x_in); + minmaxavglists->n_list[n] = n_in; #pragma acc loop independent for(l=0; llon_list[n*MAX_V+l] = x_in[l]; + minmaxavglists->lat_list[n*MAX_V+l] = y_in[l]; } } diff --git a/tools/libfrencutils/create_xgrid_util.h b/tools/libfrencutils/create_xgrid_util.h index c3209599..77abf19b 100644 --- a/tools/libfrencutils/create_xgrid_util.h +++ b/tools/libfrencutils/create_xgrid_util.h @@ -23,15 +23,13 @@ #define MAXXGRID 1e6 #endif +#include "globals.h" + #define MV 50 /* this value is small compare to earth area */ -void malloc_minmaxavg_lists(const int n, - double **lon_min_list, double **lon_max_list, double **lat_min_list, double **lat_max_list, - int **n_list, double **lon_avg, double **lon_list, double **lat_list); -void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, - double *lon_min_list, double *lon_max_list, double *lat_min_list, double *lat_max_list, - int *n_list, double *lon_avg, double *lon_list, double *lat_list); +Minmaxavglists malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists); +void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, Minmaxavglists *minmaxavglists); #pragma acc routine seq double poly_ctrlon(const double lon[], const double lat[], int n, double clon); #pragma acc routine seq From 0a82676c0648d1ffe36a7b61c62729ee65e56e06 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Tue, 10 Oct 2023 13:10:47 -0400 Subject: [PATCH 02/12] remove original lists --- tools/fregrid/conserve_interp.c | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index 7b8924f0..fdbba27b 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -86,11 +86,6 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles Minmaxavglists out_minmaxavglists; - double *lon_out_min_list=NULL, *lon_out_max_list=NULL, *lon_out_avg=NULL; - double *lat_out_min_list=NULL, *lat_out_max_list=NULL; - double *lon_out_list=NULL, *lat_out_list=NULL; - int *n2_list=NULL; - nx_out = grid_out[n].nxc; ny_out = grid_out[n].nyc; interp[n].nxgrid = 0; From 519fe502c54afbd8a9e7adc047380ad055c23e20 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Wed, 11 Oct 2023 14:44:35 -0400 Subject: [PATCH 03/12] comment out reduction because it breaks the code until further cleanup --- tools/fregrid/conserve_interp.c | 5 ++++- tools/libfrencutils/create_xgrid_acc.c | 7 ++++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index fdbba27b..36e24482 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -93,6 +93,7 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles #pragma acc enter data copyin(grid_out[n].lonc[0:(nx_out+1)*(ny_out+1)], \ grid_out[n].latc[0:(nx_out+1)*(ny_out+1)]) + //allocate memory for the lists out_minmaxavglists = malloc_minmaxavg_lists(nx_out*ny_out, &out_minmaxavglists); @@ -160,6 +161,7 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles mask, i_in, j_in, i_out, j_out, xgrid_area); for(i=0; i MASK_THRESH ) { int n0, n1, n2, n3, n1_in; double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg; @@ -92,7 +93,7 @@ int create_xgrid_2dx2d_order2_acc(const int *nlon_in, const int *nlat_in, const 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); -#pragma acc loop independent reduction(+:nxgrid) +#pragma acc loop independent //reduction(+:nxgrid) for(ij=0; ij Date: Wed, 11 Oct 2023 14:51:53 -0400 Subject: [PATCH 04/12] change malloc function to void again --- tools/fregrid/conserve_interp.c | 2 +- tools/libfrencutils/create_xgrid_util.c | 2 +- tools/libfrencutils/create_xgrid_util.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index 36e24482..1bb6ca00 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -95,7 +95,7 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles //allocate memory for the lists - out_minmaxavglists = malloc_minmaxavg_lists(nx_out*ny_out, &out_minmaxavglists); + malloc_minmaxavg_lists(nx_out*ny_out, &out_minmaxavglists); #define MAX_V 8 #pragma acc enter data create(out_minmaxavglists) diff --git a/tools/libfrencutils/create_xgrid_util.c b/tools/libfrencutils/create_xgrid_util.c index 11e283e1..366ca47e 100644 --- a/tools/libfrencutils/create_xgrid_util.c +++ b/tools/libfrencutils/create_xgrid_util.c @@ -44,7 +44,7 @@ int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double void malloc_minmaxavg_lists allocates lists to hold min, max, avg values of lat/lon coordinates *******************************************************************************/ -Minmaxavglists malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) +void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) { if(n>0){ diff --git a/tools/libfrencutils/create_xgrid_util.h b/tools/libfrencutils/create_xgrid_util.h index 77abf19b..6eace6f4 100644 --- a/tools/libfrencutils/create_xgrid_util.h +++ b/tools/libfrencutils/create_xgrid_util.h @@ -28,7 +28,7 @@ #define MV 50 /* this value is small compare to earth area */ -Minmaxavglists malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists); +void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists); void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, Minmaxavglists *minmaxavglists); #pragma acc routine seq double poly_ctrlon(const double lon[], const double lat[], int n, double clon); From ccede600b5780d8092bb9ef3375c90661f21d94e Mon Sep 17 00:00:00 2001 From: mlee03 Date: Wed, 11 Oct 2023 14:58:12 -0400 Subject: [PATCH 05/12] add free --- tools/libfrencutils/create_xgrid_util.c | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tools/libfrencutils/create_xgrid_util.c b/tools/libfrencutils/create_xgrid_util.c index 366ca47e..c9a2aa5f 100644 --- a/tools/libfrencutils/create_xgrid_util.c +++ b/tools/libfrencutils/create_xgrid_util.c @@ -47,6 +47,15 @@ int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) { + if( !minmaxavglists->lon_min_list ) free(minmaxavglists->lon_min_list); + if( !minmaxavglists->lon_max_list ) free(minmaxavglists->lon_max_list); + if( !minmaxavglists->lat_min_list ) free(minmaxavglists->lat_min_list); + if( !minmaxavglists->lat_max_list ) free(minmaxavglists->lat_max_list); + if( !minmaxavglists->n_list) free(minmaxavglists->n_list); + if( !minmaxavglists->lon_avg) free(minmaxavglists->lon_avg); + if( !minmaxavglists->lon_list) free(minmaxavglists->lon_list); + if( !minmaxavglists->lat_list) free(minmaxavglists->lat_list); + if(n>0){ minmaxavglists->lon_min_list=(double *)malloc(n*sizeof(double)); minmaxavglists->lon_max_list=(double *)malloc(n*sizeof(double)); From b3b01a562999c35f767704a84312ba212bb675d7 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Wed, 18 Oct 2023 08:09:36 -0400 Subject: [PATCH 06/12] change minmaxavlists to minmaxavg_lists --- tools/fregrid/conserve_interp.c | 30 ++++++------ tools/fregrid/globals.h | 2 +- tools/libfrencutils/create_xgrid_acc.c | 25 +++++----- tools/libfrencutils/create_xgrid_acc.h | 3 +- tools/libfrencutils/create_xgrid_util.c | 62 ++++++++++++------------- tools/libfrencutils/create_xgrid_util.h | 4 +- 6 files changed, 63 insertions(+), 63 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index 1bb6ca00..b692643d 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -84,7 +84,7 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles //START NTILES_OUT for(n=0; nlon_list[0:MAX_V*nx2*ny2], out_minmaxavglists->lat_list[0:MAX_V*nx2*ny2]) -#pragma acc data present(out_minmaxavglists->n_list[0:nx2*ny2], out_minmaxavglists->lon_avg[0:nx2*ny2]) -#pragma acc data present(out_minmaxavglists->lat_min_list[0:nx2*ny2], out_minmaxavglists->lat_max_list[0:nx2*ny2]) -#pragma acc data present(out_minmaxavglists->lon_min_list[0:nx2*ny2], out_minmaxavglists->lon_max_list[0:nx2*ny2]) +#pragma acc data present(out_minmaxavg_lists->lon_list[0:MAX_V*nx2*ny2], out_minmaxavg_lists->lat_list[0:MAX_V*nx2*ny2]) +#pragma acc data present(out_minmaxavg_lists->n_list[0:nx2*ny2], out_minmaxavg_lists->lon_avg[0:nx2*ny2]) +#pragma acc data present(out_minmaxavg_lists->lat_min_list[0:nx2*ny2], out_minmaxavg_lists->lat_max_list[0:nx2*ny2]) +#pragma acc data present(out_minmaxavg_lists->lon_min_list[0:nx2*ny2], out_minmaxavg_lists->lon_max_list[0:nx2*ny2]) #pragma acc data present(xgrid_area[0:mxxgrid], xgrid_clon[0:mxxgrid], xgrid_clat[0:mxxgrid], \ i_in[0:mxxgrid], j_in[0:mxxgrid], i_out[0:mxxgrid],j_out[0:mxxgrid]) #pragma acc data copyin(area_in[0:nx1*ny1], area_out[0:nx2*ny2]) @@ -102,17 +103,17 @@ int create_xgrid_2dx2d_order2_acc(const int *nlon_in, const int *nlat_in, const i2 = ij%nx2; j2 = ij/nx2; - if(out_minmaxavglists->lat_min_list[ij] >= lat_in_max || out_minmaxavglists->lat_max_list[ij] <= lat_in_min ) continue; + if(out_minmaxavg_lists->lat_min_list[ij] >= lat_in_max || out_minmaxavg_lists->lat_max_list[ij] <= lat_in_min ) continue; /* adjust x2_in according to lon_in_avg*/ - n2_in = out_minmaxavglists->n_list[ij]; + n2_in = out_minmaxavg_lists->n_list[ij]; #pragma acc loop seq for(l=0; llon_list[ij*MAX_V+l]; - y2_in[l] = out_minmaxavglists->lat_list[ij*MAX_V+l]; + x2_in[l] = out_minmaxavg_lists->lon_list[ij*MAX_V+l]; + y2_in[l] = out_minmaxavg_lists->lat_list[ij*MAX_V+l]; } - lon_out_min = out_minmaxavglists->lon_min_list[ij]; - lon_out_max = out_minmaxavglists->lon_max_list[ij]; - dx = out_minmaxavglists->lon_avg[ij] - lon_in_avg; + lon_out_min = out_minmaxavg_lists->lon_min_list[ij]; + lon_out_max = out_minmaxavg_lists->lon_max_list[ij]; + dx = out_minmaxavg_lists->lon_avg[ij] - lon_in_avg; if(dx < -M_PI ) { lon_out_min += TPI; lon_out_max += TPI; diff --git a/tools/libfrencutils/create_xgrid_acc.h b/tools/libfrencutils/create_xgrid_acc.h index ac679107..0913f0ff 100644 --- a/tools/libfrencutils/create_xgrid_acc.h +++ b/tools/libfrencutils/create_xgrid_acc.h @@ -24,6 +24,7 @@ int create_xgrid_2dx2d_order2_acc(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, - Minmaxavglists *out_minmaxavglists, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, + Minmaxavg_lists *out_minmaxavg_lists, 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); #endif diff --git a/tools/libfrencutils/create_xgrid_util.c b/tools/libfrencutils/create_xgrid_util.c index c9a2aa5f..3a73c26e 100644 --- a/tools/libfrencutils/create_xgrid_util.c +++ b/tools/libfrencutils/create_xgrid_util.c @@ -44,31 +44,29 @@ int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double void malloc_minmaxavg_lists allocates lists to hold min, max, avg values of lat/lon coordinates *******************************************************************************/ -void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) +void malloc_minmaxavg_lists(const int n, Minmaxavg_lists *minmaxavg_lists) { - if( !minmaxavglists->lon_min_list ) free(minmaxavglists->lon_min_list); - if( !minmaxavglists->lon_max_list ) free(minmaxavglists->lon_max_list); - if( !minmaxavglists->lat_min_list ) free(minmaxavglists->lat_min_list); - if( !minmaxavglists->lat_max_list ) free(minmaxavglists->lat_max_list); - if( !minmaxavglists->n_list) free(minmaxavglists->n_list); - if( !minmaxavglists->lon_avg) free(minmaxavglists->lon_avg); - if( !minmaxavglists->lon_list) free(minmaxavglists->lon_list); - if( !minmaxavglists->lat_list) free(minmaxavglists->lat_list); + if( !minmaxavg_lists->lon_min_list ) free(minmaxavg_lists->lon_min_list); + if( !minmaxavg_lists->lon_max_list ) free(minmaxavg_lists->lon_max_list); + if( !minmaxavg_lists->lat_min_list ) free(minmaxavg_lists->lat_min_list); + if( !minmaxavg_lists->lat_max_list ) free(minmaxavg_lists->lat_max_list); + if( !minmaxavg_lists->n_list) free(minmaxavg_lists->n_list); + if( !minmaxavg_lists->lon_avg) free(minmaxavg_lists->lon_avg); + if( !minmaxavg_lists->lon_list) free(minmaxavg_lists->lon_list); + if( !minmaxavg_lists->lat_list) free(minmaxavg_lists->lat_list); if(n>0){ - minmaxavglists->lon_min_list=(double *)malloc(n*sizeof(double)); - minmaxavglists->lon_max_list=(double *)malloc(n*sizeof(double)); - minmaxavglists->lat_min_list=(double *)malloc(n*sizeof(double)); - minmaxavglists->lat_max_list=(double *)malloc(n*sizeof(double)); - minmaxavglists->n_list=(int *)malloc(n*sizeof(int)); - minmaxavglists->lon_avg=(double *)malloc(n*sizeof(double)); - minmaxavglists->lon_list=(double *)malloc(MAX_V*n*sizeof(double)); - minmaxavglists->lat_list=(double *)malloc(MAX_V*n*sizeof(double)); + minmaxavg_lists->lon_min_list=(double *)malloc(n*sizeof(double)); + minmaxavg_lists->lon_max_list=(double *)malloc(n*sizeof(double)); + minmaxavg_lists->lat_min_list=(double *)malloc(n*sizeof(double)); + minmaxavg_lists->lat_max_list=(double *)malloc(n*sizeof(double)); + minmaxavg_lists->n_list=(int *)malloc(n*sizeof(int)); + minmaxavg_lists->lon_avg=(double *)malloc(n*sizeof(double)); + minmaxavg_lists->lon_list=(double *)malloc(MAX_V*n*sizeof(double)); + minmaxavg_lists->lat_list=(double *)malloc(MAX_V*n*sizeof(double)); } - return *minmaxavglists; - }//malloc_minmaxavg_lists /******************************************************************************* @@ -76,7 +74,7 @@ void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists) computes lists to hold min, max, avg values of lat/lon coordinates *******************************************************************************/ void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, - Minmaxavglists *minmaxavglists) + Minmaxavg_lists *minmaxavg_lists) { int nxp, nyp; @@ -85,10 +83,10 @@ void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const do nyp = ny+1; #pragma acc data present(lon[0:nxp*nyp], lat[0:nxp*nyp]) -#pragma acc data present(minmaxavglists->lon_list[0:MAX_V*nx*ny], minmaxavglists->lat_list[0:MAX_V*nx*ny]) -#pragma acc data present(minmaxavglists->n_list[0:nx*ny], minmaxavglists->lon_avg[0:nx*ny]) -#pragma acc data present(minmaxavglists->lat_min_list[0:nx*ny], minmaxavglists->lat_max_list[0:nx*ny]) -#pragma acc data present(minmaxavglists->lon_min_list[0:nx*ny], minmaxavglists->lon_max_list[0:nx*ny]) +#pragma acc data present(minmaxavg_lists->lon_list[0:MAX_V*nx*ny], minmaxavg_lists->lat_list[0:MAX_V*nx*ny]) +#pragma acc data present(minmaxavg_lists->n_list[0:nx*ny], minmaxavg_lists->lon_avg[0:nx*ny]) +#pragma acc data present(minmaxavg_lists->lat_min_list[0:nx*ny], minmaxavg_lists->lat_max_list[0:nx*ny]) +#pragma acc data present(minmaxavg_lists->lon_min_list[0:nx*ny], minmaxavg_lists->lon_max_list[0:nx*ny]) #pragma acc parallel loop independent for(int ij=0; ijlat_min_list[n] = minval_double(4, y_in); - minmaxavglists->lat_max_list[n] = maxval_double(4, y_in); + minmaxavg_lists->lat_min_list[n] = minval_double(4, y_in); + minmaxavg_lists->lat_max_list[n] = maxval_double(4, y_in); n_in = fix_lon(x_in, y_in, 4, M_PI); //Commented out for now. OpenACC does not like error_handler //if(n2_in > MAX_V) error_handler("create_xgrid.c: n_in is greater than MAX_V"); - minmaxavglists->lon_min_list[n] = minval_double(n_in, x_in); - minmaxavglists->lon_max_list[n] = maxval_double(n_in, x_in); - minmaxavglists->lon_avg[n] = avgval_double(n_in, x_in); - minmaxavglists->n_list[n] = n_in; + minmaxavg_lists->lon_min_list[n] = minval_double(n_in, x_in); + minmaxavg_lists->lon_max_list[n] = maxval_double(n_in, x_in); + minmaxavg_lists->lon_avg[n] = avgval_double(n_in, x_in); + minmaxavg_lists->n_list[n] = n_in; #pragma acc loop independent for(l=0; llon_list[n*MAX_V+l] = x_in[l]; - minmaxavglists->lat_list[n*MAX_V+l] = y_in[l]; + minmaxavg_lists->lon_list[n*MAX_V+l] = x_in[l]; + minmaxavg_lists->lat_list[n*MAX_V+l] = y_in[l]; } } diff --git a/tools/libfrencutils/create_xgrid_util.h b/tools/libfrencutils/create_xgrid_util.h index 6eace6f4..ff39b164 100644 --- a/tools/libfrencutils/create_xgrid_util.h +++ b/tools/libfrencutils/create_xgrid_util.h @@ -28,8 +28,8 @@ #define MV 50 /* this value is small compare to earth area */ -void malloc_minmaxavg_lists(const int n, Minmaxavglists *minmaxavglists); -void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, Minmaxavglists *minmaxavglists); +void malloc_minmaxavg_lists(const int n, Minmaxavg_lists *minmaxavg_lists); +void get_minmaxavg_lists(const int nx, const int ny, const double *lon, const double *lat, Minmaxavg_lists *minmaxavg_lists); #pragma acc routine seq double poly_ctrlon(const double lon[], const double lat[], int n, double clon); #pragma acc routine seq From ae674a77c32607a4affe1fa50ad11757a0ba6c58 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 09:41:05 -0400 Subject: [PATCH 07/12] test spacing --- tools/libfrencutils/create_xgrid.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index 50c9ade2..ae43860b 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -39,8 +39,8 @@ 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) + 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; From 598ab47f0e8a0d40027e5cf083058c26de0bc69b Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 10:07:31 -0400 Subject: [PATCH 08/12] test spacing --- tools/libfrencutils/create_xgrid.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index ae43860b..166d7ff0 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -39,7 +39,7 @@ 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 *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; @@ -51,9 +51,9 @@ int create_xgrid_1dx2d_order1_(const int *nlon_in, const int *nlat_in, const int }; 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; From 4a008b89bb6cf90157d1077b82cd177d713c6765 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 10:12:31 -0400 Subject: [PATCH 09/12] testing spacing --- tools/libfrencutils/create_xgrid.c | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index 166d7ff0..8ef79ff2 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -39,19 +39,19 @@ 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 *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); + 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, - const double *lat_in, const double *lon_out, const double *lat_out, + 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) { @@ -103,9 +103,9 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_out[(j2+1)*nx2p+i2+1]; y_in[3] = lat_out[(j2+1)*nx2p+i2]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; @@ -114,17 +114,17 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { + 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 ) { xgrid_area[nxgrid] = Xarea; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } From 70f6ebd8bf362fab49d1166a433d670b5efbf434 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 10:37:16 -0400 Subject: [PATCH 10/12] spacing omg --- tools/libfrencutils/create_xgrid.c | 142 ++++++++++++++--------------- 1 file changed, 71 insertions(+), 71 deletions(-) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index 8ef79ff2..4ffc2443 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -150,14 +150,14 @@ int create_xgrid_1dx2d_order2_(const int *nlon_in, const int *nlat_in, const int { 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -214,7 +214,7 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { 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]); + min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if(xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = xarea; xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); @@ -244,21 +244,21 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -317,7 +317,7 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int 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 ) { - xgrid_area[nxgrid] = Xarea; + xgrid_area[nxgrid] = Xarea; i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; @@ -345,21 +345,21 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -451,22 +451,22 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *******************************************************************************/ #ifndef __AIX 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) + 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); + i_in, j_in, i_out, j_out, xgrid_area); return nxgrid; }; #endif 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) + 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) { #define MAX_V 8 @@ -553,8 +553,8 @@ int create_xgrid_2dx2d_order1(const int *nlon_in, const int *nlat_in, const int lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double)); #if defined(_OPENMP) #pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \ - lat_out_max_list,lon_out_min_list,lon_out_max_list, \ - lon_out_avg,n2_list,lon_out_list,lat_out_list) + lat_out_max_list,lon_out_min_list,lon_out_max_list, \ + lon_out_avg,n2_list,lon_out_list,lat_out_list) #endif for(ij=0; ij M_PI) { + else if (dx > 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; + 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"); + 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; @@ -722,22 +722,22 @@ nxgrid = 0; ********************************************************************************/ #ifndef __AIX 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) + 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_2dx2d_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); + j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); return nxgrid; }; #endif #ifdef _OPENACC 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) + 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) { #define MAX_V 8 @@ -848,14 +848,14 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int } lon_out_min = lon_out_min_list[ij]; lon_out_max = lon_out_max_list[ij]; - dx = lon_out_avg[ij] - lon_in_avg; + dx = lon_out_avg[ij] - lon_in_avg; if(dx < -M_PI ) { lon_out_min += TPI; lon_out_max += TPI; #pragma acc loop seq for (l=0; l M_PI) { + else if (dx > M_PI) { lon_out_min -= TPI; lon_out_max -= TPI; #pragma acc loop seq @@ -868,7 +868,7 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue; n_out = 1; 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; + double min_area; 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 ) { @@ -904,9 +904,9 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int };/* get_xgrid_2Dx2D_order2 */ #else 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) + 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) { #define MAX_V 8 @@ -997,8 +997,8 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double)); lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double)); #pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \ - lat_out_max_list,lon_out_min_list,lon_out_max_list, \ - lon_out_avg,n2_list,lon_out_list,lat_out_list) + lat_out_max_list,lon_out_min_list,lon_out_max_list, \ + lon_out_avg,n2_list,lon_out_list,lat_out_list) for(ij=0; ij MASK_THRESH ) { @@ -1070,13 +1070,13 @@ nxgrid = 0; } lon_out_min = lon_out_min_list[ij]; lon_out_max = lon_out_max_list[ij]; - dx = lon_out_avg[ij] - lon_in_avg; + dx = lon_out_avg[ij] - lon_in_avg; if(dx < -M_PI ) { lon_out_min += TPI; lon_out_max += TPI; for (l=0; l M_PI) { + else if (dx > 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; + 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"); + 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; pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); @@ -1163,22 +1163,22 @@ nxgrid = 0; #ifndef __AIX 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) + 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_great_circle(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); + mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); return nxgrid; }; #endif 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) + 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 nx1, nx2, ny1, ny2, nx1p, nx2p, ny1p, ny2p, nxgrid, n1_in, n2_in; @@ -1240,7 +1240,7 @@ int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int x2_in[3] = x2[n3]; y2_in[3] = y2[n3]; z2_in[3] = z2[n3]; if ( (n_out = clip_2dx2d_great_circle( x1_in, y1_in, z1_in, n1_in, x2_in, y2_in, z2_in, n2_in, - x_out, y_out, z_out)) > 0) { + x_out, y_out, z_out)) > 0) { xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); if( xarea/min_area > AREA_RATIO_THRESH ) { From 9460734c122c76136f0e6181593b5ba6201c2006 Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 11:14:25 -0400 Subject: [PATCH 11/12] spacing --- tools/libfrencutils/create_xgrid.c | 566 ++++++++++++++--------------- 1 file changed, 282 insertions(+), 284 deletions(-) diff --git a/tools/libfrencutils/create_xgrid.c b/tools/libfrencutils/create_xgrid.c index 4ffc2443..1cf4c528 100644 --- a/tools/libfrencutils/create_xgrid.c +++ b/tools/libfrencutils/create_xgrid.c @@ -82,7 +82,7 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int /* This is just a temporary fix to solve the issue that there is one point in zonal direction */ // TODO: Finish this "temporary fix" if(nx1 > 1) - get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in); + get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in); else get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in); @@ -92,9 +92,9 @@ int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int for(j1=0; j1 MASK_THRESH ) { - ll_lon = lon_in[i1]; ll_lat = lat_in[j1]; - ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1]; - for(j2=0; j2=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; @@ -213,19 +213,19 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -244,21 +244,21 @@ int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -302,9 +302,9 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int y_in[2] = lat_in[(j1+1)*nx1p+i1+1]; y_in[3] = lat_in[(j1+1)*nx1p+i1]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) - && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; + && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) - && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; + && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_in[j1*nx1p+i1]; x_in[1] = lon_in[j1*nx1p+i1+1]; @@ -314,17 +314,17 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = Xarea; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = Xarea; + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -345,21 +345,21 @@ int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int 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) + 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); + 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, - 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) + 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 nx1, ny1, nx2, ny2, nx1p, nx2p; @@ -418,19 +418,19 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { - 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + 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 ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } @@ -451,22 +451,22 @@ int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *******************************************************************************/ #ifndef __AIX 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) + 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); + i_in, j_in, i_out, j_out, xgrid_area); return nxgrid; }; #endif 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) + 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) { #define MAX_V 8 @@ -587,10 +587,10 @@ nxgrid = 0; #if defined(_OPENMP) #pragma omp parallel for default(none) shared(nblocks,nx1,ny1,nx1p,mask_in,lon_in,lat_in, \ - istart2,iend2,nx2,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,nthreads) + istart2,iend2,nx2,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,nthreads) #endif for(m=0; m= 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= 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; - } - - } - + if(lon_out_min >= 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; + } + } } } } @@ -682,13 +680,13 @@ nxgrid = 0; nxgrid = 0; for(m=0; m= 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]; + 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]; #pragma acc loop seq - for(l=0; l M_PI) { - lon_out_min -= TPI; - lon_out_max -= TPI; + for (l=0; l M_PI) { + lon_out_min -= TPI; + lon_out_max -= TPI; #pragma acc loop seq - for (l=0; l= lon_in_max || lon_out_max <= lon_in_min ) continue; - n_out = 1; - 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; - 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 ) { - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; + /* x2_in should in the same range as x1_in after lon_fix, so no need to + consider cyclic condition + */ + if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue; + n_out = 1; + 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; + 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 ) { + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; #pragma atomic update - nxgrid++; - } - } + nxgrid++; + } + } } } } @@ -904,9 +902,9 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int };/* get_xgrid_2Dx2D_order2 */ #else 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) + 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) { #define MAX_V 8 @@ -997,8 +995,8 @@ int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double)); lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double)); #pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \ - lat_out_max_list,lon_out_min_list,lon_out_max_list, \ - lon_out_avg,n2_list,lon_out_list,lat_out_list) + lat_out_max_list,lon_out_min_list,lon_out_max_list, \ + lon_out_avg,n2_list,lon_out_list,lat_out_list) for(ij=0; ij MASK_THRESH ) { int n0, n1, n2, n3, l,n1_in; @@ -1054,57 +1052,57 @@ nxgrid = 0; 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_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%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= 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; - pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); - pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out ); - pi_in[nn] = i1; - pj_in[nn] = j1; - pi_out[nn] = i2; - pj_out[nn] = j2; - } - } + /* x2_in should in the same range as x1_in after lon_fix, so no need to + * consider cyclic condition + * */ + if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue; + if ( (n_out = clip_2dx2d( 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; + pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); + pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out ); + pi_in[nn] = i1; + pj_in[nn] = j1; + pi_out[nn] = i2; + pj_out[nn] = j2; + } + } } } } @@ -1125,15 +1123,15 @@ nxgrid = 0; nxgrid = 0; for(m=0; m 0) { - xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; - min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); - if( xarea/min_area > AREA_RATIO_THRESH ) { + x_out, y_out, z_out)) > 0) { + xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; + min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); + if( xarea/min_area > AREA_RATIO_THRESH ) { #ifdef debug_test_create_xgrid - printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea); + printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea); #endif - xgrid_area[nxgrid] = xarea; - xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */ - xgrid_clat[nxgrid] = 0; - i_in[nxgrid] = i1; - j_in[nxgrid] = j1; - i_out[nxgrid] = i2; - j_out[nxgrid] = j2; - ++nxgrid; - if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); - } + xgrid_area[nxgrid] = xarea; + xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */ + xgrid_clat[nxgrid] = 0; + i_in[nxgrid] = i1; + j_in[nxgrid] = j1; + i_out[nxgrid] = i2; + j_out[nxgrid] = j2; + ++nxgrid; + if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); + } } } } From bfd5ce0c6b23d3998f4ee91944773ac6c2e2139f Mon Sep 17 00:00:00 2001 From: mlee03 Date: Thu, 19 Oct 2023 11:40:29 -0400 Subject: [PATCH 12/12] nullify the lists --- tools/fregrid/conserve_interp.c | 8 +++++ tools/libfrencutils/create_xgrid_util.c | 40 ++++++++++++++++++++----- 2 files changed, 40 insertions(+), 8 deletions(-) diff --git a/tools/fregrid/conserve_interp.c b/tools/fregrid/conserve_interp.c index b692643d..ae90ed5f 100644 --- a/tools/fregrid/conserve_interp.c +++ b/tools/fregrid/conserve_interp.c @@ -85,6 +85,14 @@ void setup_conserve_interp(int ntiles_in, const Grid_config *grid_in, int ntiles for(n=0; nlon_min_list ) free(minmaxavg_lists->lon_min_list); - if( !minmaxavg_lists->lon_max_list ) free(minmaxavg_lists->lon_max_list); - if( !minmaxavg_lists->lat_min_list ) free(minmaxavg_lists->lat_min_list); - if( !minmaxavg_lists->lat_max_list ) free(minmaxavg_lists->lat_max_list); - if( !minmaxavg_lists->n_list) free(minmaxavg_lists->n_list); - if( !minmaxavg_lists->lon_avg) free(minmaxavg_lists->lon_avg); - if( !minmaxavg_lists->lon_list) free(minmaxavg_lists->lon_list); - if( !minmaxavg_lists->lat_list) free(minmaxavg_lists->lat_list); + if( minmaxavg_lists->lon_min_list != NULL ) { + free(minmaxavg_lists->lon_min_list); + minmaxavg_lists-> lon_min_list = NULL; + } + if( minmaxavg_lists->lon_max_list != NULL ) { + free(minmaxavg_lists->lon_max_list); + minmaxavg_lists->lon_max_list = NULL; + } + if( minmaxavg_lists->lat_min_list != NULL ) { + free(minmaxavg_lists->lat_min_list); + minmaxavg_lists->lat_min_list = NULL; + } + if( minmaxavg_lists->lat_max_list != NULL ) { + free(minmaxavg_lists->lat_max_list); + minmaxavg_lists->lat_max_list = NULL; + } + if( minmaxavg_lists->n_list != NULL ) { + free(minmaxavg_lists->n_list); + minmaxavg_lists->n_list = NULL; + } + if( minmaxavg_lists->lon_avg != NULL ) { + free(minmaxavg_lists->lon_avg); + minmaxavg_lists->lon_avg = NULL; + } + if( minmaxavg_lists->lon_list != NULL ) { + free(minmaxavg_lists->lon_list); + minmaxavg_lists->lon_list = NULL; + } + if( minmaxavg_lists->lat_list != NULL ) { + free(minmaxavg_lists->lat_list); + minmaxavg_lists->lat_list = NULL; + } if(n>0){ minmaxavg_lists->lon_min_list=(double *)malloc(n*sizeof(double));