From f329b393ab5527470810280b7ce0c7b1bfdebc83 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Sat, 4 Nov 2023 09:38:28 -0600 Subject: [PATCH] While determining masked blocks, take reentrancy and tripolar stitch into account --- config_src/infra/FMS1/MOM_domain_infra.F90 | 13 ++- src/framework/MOM_domains.F90 | 102 ++++++++++++++------- 2 files changed, 81 insertions(+), 34 deletions(-) diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90 index 2c97a0bb31..08555436f7 100644 --- a/config_src/infra/FMS1/MOM_domain_infra.F90 +++ b/config_src/infra/FMS1/MOM_domain_infra.F90 @@ -40,7 +40,7 @@ module MOM_domain_infra public :: domain2D, domain1D, group_pass_type ! These interfaces are actually implemented or have explicit interfaces in this file. public :: create_MOM_domain, clone_MOM_domain, get_domain_components, get_domain_extent -public :: deallocate_MOM_domain, get_global_shape, compute_block_extent +public :: deallocate_MOM_domain, get_global_shape, compute_block_extent, compute_extent public :: pass_var, pass_vector, fill_symmetric_edges, rescale_comp_data public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete public :: create_group_pass, do_group_pass, start_group_pass, complete_group_pass @@ -1945,6 +1945,17 @@ subroutine compute_block_extent(isg, ieg, ndivs, ibegin, iend) call mpp_compute_block_extent(isg, ieg, ndivs, ibegin, iend) end subroutine compute_block_extent +!> Get the array ranges in one dimension for the divisions of a global index space +subroutine compute_extent(isg, ieg, ndivs, ibegin, iend) + integer, intent(in) :: isg !< The starting index of the global index space + integer, intent(in) :: ieg !< The ending index of the global index space + integer, intent(in) :: ndivs !< The number of divisions + integer, dimension(:), intent(out) :: ibegin !< The starting index of each division + integer, dimension(:), intent(out) :: iend !< The ending index of each division + + call mpp_compute_extent(isg, ieg, ndivs, ibegin, iend) +end subroutine compute_extent + !> Broadcast a 2-d domain from the root PE to the other PEs subroutine broadcast_domain(domain) type(domain2d), intent(inout) :: domain !< The domain2d type that will be shared across PEs. diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index fe191e166c..fa46299520 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -299,7 +299,8 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & ! Auto-generate a mask file and determine the layout call cpu_clock_begin(id_clock_auto_mask) if (is_root_PE()) then - call gen_auto_mask_table(n_global, PEs_used, param_file, inputdir, auto_mask_table_fname, layout) + call gen_auto_mask_table(n_global, reentrant, tripolar_N, PEs_used, param_file, inputdir, & + auto_mask_table_fname, layout) endif call broadcast(layout, length=2) call cpu_clock_end(id_clock_auto_mask) @@ -434,38 +435,42 @@ subroutine MOM_define_layout(n_global, ndivs, layout) end subroutine MOM_define_layout !> Given a desired number of active npes, generate a layout and mask_table -subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, layout) +subroutine gen_auto_mask_table(n_global, reentrant, tripolar_N, npes, param_file, inputdir, filename, layout) integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions - integer, intent(in) :: npes !< The desired number of active PEs. + logical, dimension(2), intent(in) :: reentrant !< True if the x- and y- directions are periodic. + logical :: tripolar_N !< A flag indicating whether there is n. tripolar connectivity + integer, intent(in) :: npes !< The desired number of active PEs. type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters character(len=128), intent(in) :: inputdir !< INPUTDIR parameter> character(len=:), allocatable, intent(in) :: filename !< Mask table file path (to be auto-generated.) integer, dimension(2), intent(out) :: layout !< The generated layout of PEs (incl. masked blocks) !local - real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE) - integer, dimension(n_global(1), n_global(2)) :: mask ! Cell masks (based on D and MINIMUM_DEPTH) + real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE) + integer, dimension(:,:), allocatable :: mask ! Cell masks (based on D and MINIMUM_DEPTH) character(len=200) :: topo_filepath, topo_file ! Strings for file/path character(len=200) :: topo_varname ! Variable name in file character(len=200) :: topo_config character(len=40) :: mdl = "gen_auto_mask_table" ! This subroutine's name. integer :: i, j, p - real :: Dmask ! The depth for masking in the same units as D - real :: min_depth ! The minimum ocean depth in the same units as D - real :: mask_depth ! The depth shallower than which to mask a point as land. + real :: Dmask ! The depth for masking in the same units as D + real :: min_depth ! The minimum ocean depth in the same units as D + real :: mask_depth ! The depth shallower than which to mask a point as land. real :: glob_ocn_frac ! ratio of ocean points to total number of points - real :: r_p ! aspect ratio for division count p. + real :: r_p ! aspect ratio for division count p. + integer :: nx, ny ! global domain sizes + integer, parameter :: ibuf=2, jbuf=2 real, parameter :: r_extreme = 10.0 ! aspect ratio limit (>1) for a layout to be considered. integer :: num_masked_blocks integer, allocatable :: mask_table(:,:) ! Read in params necessary for auto-masking - call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, do_not_log=.true., units="m", default=0.0) + call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, do_not_log=.true., units="m", default=0.0) call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, do_not_log=.true., units="m", default=-9999.0) call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, do_not_log=.true., fail_if_missing=.true.) call get_param(param_file, mdl, "TOPO_FILE", topo_file, do_not_log=.true., default="topog.nc") call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, do_not_log=.true., default="depth") topo_filepath = trim(inputdir)//trim(topo_file) - + ! Sanity checks if (.not. is_root_pe()) then call MOM_error(FATAL, 'gen_auto_mask_table should only be called by the root PE.') @@ -477,36 +482,62 @@ subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, l call MOM_error(FATAL, " gen_auto_mask_table: Unable to open "//trim(topo_filepath)) endif + nx = n_global(1) + ny = n_global(2) + ! Read in bathymetric depth. D(:,:) = -9.0e30 ! Initializing to a very large negative depth (tall mountains) everywhere. call read_field(topo_filepath, trim(topo_varname), D, start=(/1, 1/), nread=n_global, no_domain=.true.) + allocate( mask(nx+2*ibuf, ny+2*jbuf), source=0) + ! Determine cell masks Dmask = mask_depth if (mask_depth == -9999.0) Dmask = min_depth - do i=1,n_global(1) ; do j=1,n_global(2) + do i=1,nx ; do j=1,ny if (D(i,j) <= Dmask) then - mask(i,j) = 0 + mask(i+ibuf,j+jbuf) = 0 else - mask(i,j) = 1 + mask(i+ibuf,j+jbuf) = 1 endif enddo ; enddo - glob_ocn_frac = real(sum(mask)) / size(mask) + ! fill in buffer cells + + if (reentrant(1)) then ! REENTRANT_X + mask(1:ibuf, :) = mask(nx+1:nx+ibuf, :) + mask(ibuf+nx+1:nx+2*ibuf, :) = mask(ibuf+1:2*ibuf, :) + endif + + if (reentrant(2)) then ! REENTRANT_Y + mask(:, 1:jbuf) = mask(:, ny+1:ny+jbuf) + mask(:, jbuf+ny+1:ny+2*jbuf) = mask(:, jbuf+1:2*jbuf) + endif + + if (tripolar_N) then ! TRIPOLAR_N + do i=1,nx+2*ibuf + do j=1,jbuf + mask(i, jbuf+ny+j) = mask(nx+2*ibuf+1-i, jbuf+ny+1-j) + enddo + enddo + !mask(:, jbuf+ny) = 1 ! TODO - REMOVE + endif + + glob_ocn_frac = real(sum(mask(1+ibuf:nx+ibuf, 1+jbuf:ny+jbuf))) / (nx * ny) ! Iteratively check for all possible division counts starting from the upper bound of npes/glob_ocn_frac, - ! which is over-optimistic for realistic domains, but may be satisfied with idealized domains. + ! which is over-optimistic for realistic domains, but may be satisfied with idealized domains. do p = ceiling(npes/glob_ocn_frac), npes, -1 - + ! compute the layout for the current division count, p call MOM_define_layout(n_global, p, layout) ! don't bother checking this p if the aspect ratio is extreme - r_p = (real(n_global(1))/layout(1)) / (real(n_global(2))/layout(2)) + r_p = (real(nx)/layout(1)) / (real(ny)/layout(2)) if ( r_p * r_extreme < 1 .or. r_extreme < r_p ) cycle ! Get the number of masked_blocks for this particular division count - call determine_land_blocks(mask, n_global, layout(1), layout(2), num_masked_blocks) + call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks) ! If we can eliminate enough blocks to reach the target npes, adopt ! this p (and the associated layout) and terminate the iteration. @@ -521,21 +552,27 @@ subroutine gen_auto_mask_table(n_global, npes, param_file, inputdir, filename, l "of MOM6 PEs or set AUTO_MASKTABLE to False.") endif - ! Call determine_land_blocks once again, this time to retrieve the mask_table. + ! Call determine_land_blocks once again, this time to retrieve and write out the mask_table. allocate(mask_table(num_masked_blocks,2)) - call determine_land_blocks(mask, n_global, layout(1), layout(2), num_masked_blocks, mask_table) + call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks, mask_table) call write_auto_mask_file(mask_table, layout, npes, filename) deallocate(mask_table) + deallocate(mask) end subroutine gen_auto_mask_table -! Given number of domain divisions, compute the max number of land blocks that can be eliminated, +! Given a number of domain divisions, compute the max number of land blocks that can be eliminated, ! and return the resulting mask table if requested. -subroutine determine_land_blocks(mask, n_global, idiv, jdiv, num_masked_blocks, mask_table) +subroutine determine_land_blocks(mask, nx, ny, idiv, jdiv, ibuf, jbuf, num_masked_blocks, mask_table) integer, dimension(:,:), intent(in) :: mask !< cell masks based on depth and MINIMUM_DEPTH - integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions - integer, intent(in) :: idiv !< number of divisions along x axis - integer, intent(in) :: jdiv !< number of divisions along y axis + integer, intent(in) :: nx !< Total number of gridpoints in x-dir (global) + integer, intent(in) :: ny !< Total number of gridpoints in y-dir (global) + integer, intent(in) :: idiv !< number of divisions along x-dir + integer, intent(in) :: jdiv !< number of divisions along y-dir + integer, intent(in) :: ibuf !< number of buffer cells in x-dir. + !! (not necessarily the same as NIHALO) + integer, intent(in) :: jbuf !< number of buffer cells in y-dir. + !! (not necessarily the same as NJHALO) integer, intent(out) :: num_masked_blocks integer, intent(out), optional :: mask_table(:,:) ! integer @@ -543,20 +580,19 @@ subroutine determine_land_blocks(mask, n_global, idiv, jdiv, num_masked_blocks, integer, dimension(idiv) :: iend !< The ending index of each division along x axis integer, dimension(jdiv) :: jbegin !< The starting index of each division along y axis integer, dimension(jdiv) :: jend !< The ending index of each division along y axis - integer, parameter :: ibuf=2, jbuf=2 integer :: i, j, ib, ie, jb,je - call compute_extent(1, n_global(1), idiv, ibegin, iend) - call compute_extent(1, n_global(2), jdiv, jbegin, jend) + call compute_extent(1, nx, idiv, ibegin, iend) + call compute_extent(1, ny, jdiv, jbegin, jend) num_masked_blocks = 0 do i=1,idiv - ib = max(ibegin(i)-ibuf, 1) - ie = min(iend(i)+ibuf, n_global(1)) + ib = ibegin(i) + ie = iend(i) + 2 * ibuf do j=1,jdiv - jb = max(jbegin(j)-jbuf, 1) - je = min(jend(j)+jbuf, n_global(2)) + jb = jbegin(j) + je = jend(j) + 2 * jbuf if (any(mask(ib:ie,jb:je)==1)) cycle