set_bathymetry_spoon Subroutine

   subroutine set_bathymetry_spoon(b, grid, max_depth, edge_depth, slope_scale)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Fill `b(:,:)` with the MOM6 spoon bathymetry.  In Cartesian
      !! terms (we collapse MOM6's lat/lon factors of `R_earth · π / 180`
      !! into a direct meters scale), the local depth is
      !!
      !!     D(i,j) = D_edge
      !!            + D_0 · sin(π · x_phys / x_len)
      !!                  · (1 − exp((y_phys − y_len) / slope_scale))
      !!
      !! with `D_0 = (max_depth − D_edge) / (1 − exp(−0.5 · y_len /
      !! slope_scale))²`.  Zero at east/west walls, full depth at the
      !! south wall, exponentially decaying to `D_edge` at the north
      !! wall.  Physical interior only; ghost cells retain whatever
      !! value `b` had on entry.
      !!
      !! 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 `x_len` / `y_len` and the
      !! index→position map describe the WHOLE domain on every rank.  On a
      !! single rank the offsets are 0 and global == local, so the fill is
      !! byte-identical to the undecomposed formula.
      real(wp), intent(inout) :: b(:, :)
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: max_depth, edge_depth, slope_scale
      real(wp), parameter :: PI = 4.0_wp*atan(1.0_wp)
      real(wp) :: x_len, y_len, x_phys, y_phys, D_0, denom, D_local
      integer :: i, j, i_phys, j_phys, ng
      integer :: ioff, joff

      ioff = grid%i_offset_global
      joff = grid%j_offset_global

      ng = grid%nghost
      x_len = real(grid%nx_global, wp)*grid%dx
      y_len = real(grid%ny_global, wp)*grid%dy
      denom = 1.0_wp - exp(-0.5_wp*y_len/slope_scale)
      D_0 = (max_depth - edge_depth)/(denom*denom)

      ! Fill the full array including ghost rows.  Ghost cells left at
      ! the alloc-time zero cause `h_layer = 0` there, which sends the
      ! EOS into its vanishing-layer fallback (`rho_layer = rho_0`).
      ! That puts a non-physical density jump at every wall-adjacent
      ! face whenever the IC has `T_init /= T_ref` or `S_init /= S_ref`,
      ! and the BPG operator picks it up as a spurious horizontal
      ! pressure gradient.  Confirmed via the seamount test.
      do j = 1, size(b, 2)
         j_phys = j - ng
         ! Global physical y position: local j_phys shifted by joff.
         y_phys = (real(j_phys + joff, wp) - 0.5_wp)*grid%dy
         do i = 1, size(b, 1)
            i_phys = i - ng
            ! Global physical x position: local i_phys shifted by ioff.
            x_phys = (real(i_phys + ioff, wp) - 0.5_wp)*grid%dx
            D_local = edge_depth &
                      + D_0*sin(PI*x_phys/x_len) &
                      *(1.0_wp - exp((y_phys - y_len)/slope_scale))
            ! Safety: the sin term goes negative outside [0, x_len],
            ! and the exp term saturates the depth at the north wall.
            ! Clamp to [edge_depth, max_depth].
            if (D_local < edge_depth) D_local = edge_depth
            if (D_local > max_depth) D_local = max_depth
            b(i, j) = D_local
         end do
      end do
   end subroutine set_bathymetry_spoon