set_bathymetry_isomip_plus Subroutine

   subroutine set_bathymetry_isomip_plus(b, grid, max_depth, x_origin, m_per_grid)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !!
      !! Fill `b(:,:)` with the MISMIP+ / ISOMIP+ analytic bedrock,
      !! Asay-Davis et al. (2016) Eqs. (1)-(4) + Table 1:
      !!
      !!     z_b(x,y) = max( Bx(x) + By(y), z_b,deep )      [Eq. (1)]
      !!
      !! `z_b` is an ELEVATION (positive up, sea level at 0) and is
      !! negative throughout the ISOMIP+ box; Roundabout's `b` is a
      !! DEPTH (positive down), so the last line is `b = -z_b`, floored
      !! at 0 so that a bed which the formula puts ABOVE sea level (it
      !! does for `x < ~140 km`, outside the ISOMIP+ box but reachable if
      !! a caller sets a smaller `x_origin`) is reported as dry land
      !! rather than as a negative depth.  `b = 0` is below
      !! `LAND_DEPTH_THRESHOLD`, so `seed_wet_mask_impl` masks the column
      !! out through the ordinary land path — there is no ISOMIP+ branch
      !! anywhere downstream.
      !!
      !! `max_depth` is the deep clip, i.e. `-z_b,deep`; the protocol
      !! value is `ISOMIP_ZB_DEEP` ⇒ `&ocean_topo_nml max_depth = 720.0`.
      !!
      !! MINIMUM WATER COLUMN.  The protocol (their Sect. 3.1.5) asks for
      !! "the minimum ocean column as thin as can reasonably be achieved"
      !! and leaves the value to the modeller, with the choice being
      !! either to modify the topography or to mark the column land.
      !! Roundabout takes the second option and it is NOT this routine's
      !! job: `&ocean_cavity_dyn_nml h_min_cavity` is the threshold and
      !! `seed_wet_mask_impl(water, land_cutoff=h_min_cavity)` is where it
      !! bites, on `water = b - z_draft`.
      !!
      !! UNITS.  Every Table-1 constant is METRES, as printed.  Grid
      !! positions are in GRID coordinate units (metres on Cartesian,
      !! DEGREES on spherical/curvilinear), so `m_per_grid` converts them
      !! to metres before the formula sees them — the inverse of the
      !! `topo_length_to_grid_units` conversion the spoon/seamount
      !! dispatch applies to `slope_scale`, and the same trap.  On a
      !! Cartesian grid `m_per_grid = 1` exactly and this is the identity.
      !! (The protocol prescribes a Cartesian box; the conversion exists
      !! so a curvilinear caller degrades predictably rather than
      !! silently collapsing the basin flat.)
      !!
      !! Fills the FULL array INCLUDING ghost rows by evaluating the
      !! formula at the ghost index — the CLAUDE.md rule every formula
      !! bathymetry setter follows; a ghost row left at the alloc-time
      !! zero sends the EOS into its `rho_0` vanishing-layer fallback and
      !! puts a spurious density jump at every wall-adjacent face.
      !!
      !! MPI: positions come off the GLOBAL index offsets + extents, so
      !! each rank fills its window of ONE global bed.  Single rank ⇒
      !! offsets 0 ⇒ byte-identical to the undecomposed formula.
      real(wp), intent(inout) :: b(:, :)
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: max_depth
         !! Deep clip (m, positive down) = `-z_b,deep`.  Protocol: 720.
      real(wp), intent(in) :: x_origin
         !! Absolute MISMIP+ x (m) of the domain's west edge.  Protocol
         !! (ISOMIP+): 320e3.
      real(wp), intent(in) :: m_per_grid
         !! Metres per grid coordinate unit (1 on Cartesian).

      real(wp) :: y_len, x_m, y_m, by, zb
      integer :: i, j, ng, nxt, nyt, ioff, joff

      ng = grid%nghost
      ioff = grid%i_offset_global
      joff = grid%j_offset_global
      nxt = size(b, 1)
      nyt = size(b, 2)
      y_len = real(grid%ny_global, wp)*grid%dy*m_per_grid

      do j = 1, nyt
         y_m = (real(j - ng + joff, wp) - 0.5_wp)*grid%dy*m_per_grid
         ! `By` depends on y alone — hoisted out of the i loop.
         by = isomip_plus_by(y_m, y_len)
         do i = 1, nxt
            x_m = x_origin + (real(i - ng + ioff, wp) - 0.5_wp)*grid%dx*m_per_grid
            zb = isomip_plus_bx(x_m) + by
            if (zb < -max_depth) zb = -max_depth     ! Eq. (1) deep clip
            b(i, j) = max(-zb, 0.0_wp)               ! elevation -> depth
         end do
      end do
   end subroutine set_bathymetry_isomip_plus