set_bathymetry_double_drake Subroutine

   subroutine set_bathymetry_double_drake(b, grid, max_depth, half_frac)
      !! "Double Drake" idealised supercontinent (Ferreira, Marshall &
      !! Campin 2010, *J. Climate*): a flat-bottom global ocean at
      !! `max_depth` with TWO thin meridional wall-continents 90° of
      !! longitude apart, each running from the north pole down to a
      !! southern-channel latitude, leaving a reentrant circumpolar
      !! channel (Drake-Passage analogue) to the south.  Used as the
      !! static-land-mask integration showcase on a spherical periodic-x
      !! sector grid.
      !!
      !! Index-space geometry (mapped to the Ferreira config by the
      !! caller's grid lon/lat extent): physical-x is the 360° periodic
      !! longitude; the two walls sit at the western seam (`i_phys=1`)
      !! and a quarter of the way across (`i_phys≈nx_phys/4`, ≡ 90°
      !! apart on a 360° domain).  Each wall is one cell wide and spans
      !! the NORTHERN `(1-2*half_frac)` fraction of the physical-y axis
      !! (north = pole), leaving the southern `2*half_frac` band fully
      !! open as the reentrant channel.  `half_frac` is taken from
      !! `&ocean_topo_nml slope_scale` (no new knob), default 0.2
      !! ⇒ southern 40 % of the basin is the open channel.
      !!
      !! Fills the FULL array incl. ghosts (formula-bathy ghost-fill
      !! gotcha): the walls extend through the y-ghost rows so the
      !! periodic-x seam wall is land in the halo too; the southern
      !! channel + interior basin stay ocean.
      !!
      !! MPI: the global physical-index offsets and the global domain
      !! extents come off `grid` (`i_offset_global` / `j_offset_global`,
      !! `nx_global` / `ny_global`), so the two wall meridians and the
      !! channel latitude are located on the WHOLE domain and each rank
      !! carves only the part inside its own tile.  On a single rank the
      !! offsets are 0 and global == local, so the fill is byte-identical.
      real(wp), intent(inout) :: b(:, :)
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: max_depth, half_frac
      integer :: i, j, i_phys, j_phys, ng, iw1, iw2, j_chan_top
      integer :: ioff, joff, nxg, nyg
      real(wp) :: hf

      ioff = grid%i_offset_global
      joff = grid%j_offset_global
      nxg = grid%nx_global
      nyg = grid%ny_global

      ng = grid%nghost
      hf = half_frac
      if (hf <= 0.0_wp) hf = 0.2_wp
      if (hf > 0.49_wp) hf = 0.49_wp

      ! Two one-cell-wide meridional walls, 90° apart on the GLOBAL 360° axis.
      iw1 = 1                   ! western-seam global wall at global i_phys = 1
      iw2 = max(2, nxg/4)       ! a quarter across ≡ 90° of longitude globally
      ! Channel occupies the southern `2*hf` fraction of the GLOBAL domain;
      ! walls cover the rest (the northern band up to the pole).
      j_chan_top = nint(2.0_wp*hf*real(nyg, wp))

      b = max_depth
      do j = 1, size(b, 2)
         j_phys = j - ng
         do i = 1, size(b, 1)
            i_phys = i - ng
            ! Wall cells: on either meridional line AND north of the
            ! reentrant channel.  Compare GLOBAL physical indices.
            ! j-ghosts north of the channel are wall (walls reach the pole),
            ! south ghosts stay ocean (open channel).
            if (((i_phys + ioff) == iw1 .or. (i_phys + ioff) == iw2) .and. &
                (j_phys + joff) > j_chan_top) then
               b(i, j) = 0.0_wp   ! LAND (< LAND_DEPTH_THRESHOLD)
            end if
         end do
      end do
   end subroutine set_bathymetry_double_drake