While determining masked blocks, take reentrancy and tripolar stitch …

…into account

config_src/infra/FMS1/MOM_domain_infra.F90

@@ -16,7 +16,7 @@ module MOM_domain_infra
 use mpp_domains_mod, only : mpp_create_group_update, mpp_do_group_update
 use mpp_domains_mod, only : mpp_reset_group_update_field, mpp_group_update_initialized
 use mpp_domains_mod, only : mpp_start_group_update, mpp_complete_group_update
-use mpp_domains_mod, only : mpp_compute_block_extent
+use mpp_domains_mod, only : mpp_compute_block_extent, mpp_compute_extent
 use mpp_domains_mod, only : mpp_broadcast_domain, mpp_redistribute, mpp_global_field
 use mpp_domains_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM
 use mpp_domains_mod, only : CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE
@@ -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.

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,42 +552,47 @@ 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
   integer, dimension(idiv) :: ibegin   !< The starting index of each division along x axis
   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