set_bathymetry_seamount Subroutine

   subroutine set_bathymetry_seamount(b, grid, max_depth, peak_depth, half_width)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Fill `b(:,:)` with a centred Gaussian seamount bathymetry:
      !!
      !!     D(i,j) = max_depth − (max_depth − peak_depth)
      !!                          · exp(−((x − xc)² + (y − yc)²) / L²)
      !!
      !! with `(xc, yc)` at the centre of the physical domain and
      !! `L = half_width` (gaussian e-folding distance) — expressed in the
      !! GRID coordinate units (metres on Cartesian, degrees on spherical/
      !! curvilinear).  The production dispatch converts the metres
      !! `slope_scale` knob via `topo_length_to_grid_units`; callers that
      !! pass `half_width` directly must match the grid's units.
      !! Reaches `peak_depth` at the bump centre, asymptotes to
      !! `max_depth` far from the bump.
      !!
      !! Designed as a Tier-1.5 bridge test between flat-bottom
      !! analytical setups and full real-bathymetry regional runs.
      !! Tests σ-coord pressure-gradient and Coriolis-advection
      !! kernels under a varying h_layer without the confounds of
      !! a 50× depth ratio + IC stratification + closed-basin
      !! resonance that `tasman_2km.nml` introduces.
      !!
      !! Standard test protocol: uniform T,S (no APE), zero
      !! Coriolis or f-plane only, no wind, walls all sides.  The
      !! expected steady state is **identically zero motion** — any
      !! non-zero u, v, η is a σ-coord-related numerical artefact.
      !!
      !! 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 `(xc, yc)` sits at the centre of the
      !! WHOLE domain on every rank rather than the centre of each tile.
      !! 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, peak_depth, half_width
      real(wp) :: x_len, y_len, xc, yc, x_phys, y_phys, r2, depression
      integer :: i, j, i_phys, j_phys, ng, nx_total, ny_total
      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
      xc = 0.5_wp*x_len
      yc = 0.5_wp*y_len
      depression = max_depth - peak_depth
      nx_total = size(b, 1)
      ny_total = size(b, 2)

      ! Fill the full array including ghost rows so the EOS, BPG, and
      ! mass-flux operators see a smooth bathymetry on every face they
      ! touch.  Leaving ghosts at zero injects a spurious density jump
      ! at wall-adjacent faces (EOS falls back to `rho_0` where
      ! `h_layer = 0`) — the smoking-gun bug the seamount test was built
      ! to expose.
      do j = 1, ny_total
         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, nx_total
            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
            r2 = (x_phys - xc)**2 + (y_phys - yc)**2
            b(i, j) = max_depth - depression*exp(-r2/(half_width*half_width))
         end do
      end do
   end subroutine set_bathymetry_seamount