seed_h_layer_uniform_z_impl Subroutine

   pure subroutine seed_h_layer_uniform_z_impl(h_layer, b, nz_ml, max_depth, angstrom_h)
      !! `thickness_config = "uniform_z"`: MOM6 `initialize_thickness_uniform`
      !! port.  Lays uniform **z** interfaces
      !! over the GLOBAL `max_depth`, clips them bottom-up against the local
      !! bathymetry, and collapses whatever will not fit to a minimum-thickness
      !! floor:
      !!
      !!     z_top_target(k) = -max_depth * (nz_ml - k) / nz_ml
      !!     h(k) = max(z_top_target(k) - z_bot(k), h_floor)
      !!
      !! walking `k = 1 -> nz_ml` from the bed up (Roundabout's bottom-up
      !! convention; MOM6 walks `k = nz -> 1` from its surface-first index, the
      !! same sweep).  The surface layer's target is `z = 0`, so it absorbs the
      !! remainder and the telescoping sum is EXACTLY `b(i,j)` — the floor
      !! injects nothing, unlike the runtime `angstrom_h` clamp.
      !!
      !! Why this exists.  Under a horizontally-uniform density stack (the
      !! `rho_lightest`/`rho_range` linear-coordinate IC, MOM6
      !! `COORD_CONFIG="linear"`) the layer index IS the density, so the layer
      !! interfaces ARE the isopycnals.  Seeding `b/nz_ml` (`"sigma"`) therefore
      !! makes every isopycnal follow the bathymetry, which stands the whole
      !! `rho_range` contrast up across each shelf break at t=0: on the MOM6
      !! double-gyre spoon that is ~1.9 kg/m3 between the 100 m rim and the
      !! 2000 m interior, i.e. g' ~ 0.018 m/s2 driving a ~1.3 m/s gravity
      !! current over layers only `edge_depth/nz_ml` thick.  `"uniform_z"`
      !! instead gives flat resting isopycnals with the sub-bathymetry layers
      !! collapsed against the bed, which is the layered-model resting state.
      !!
      !! Deep columns are unaffected where `b >= max_depth`: every layer lands
      !! on its z target at `max_depth/nz_ml`, and any excess `b - max_depth`
      !! is absorbed by the bed-most layer (`k = 1`).
      !!
      !! `h_floor = max(angstrom_h, 2*H_VANISHED)`.  The `angstrom_h` term ties
      !! the seed to the same floor the Lagrangian continuity update uses, so a
      !! layer seeded as collapsed reads as vanished to
      !! `isopycnal_vanish_tol()` (`max(angstrom_h, H_VANISHED)`) and is picked
      !! up by `reset_vanished_u` / `cfl_ignore_vanished` from step 1.  The
      !! `2*H_VANISHED` term keeps the floor strictly positive when
      !! `angstrom_h = 0` (knob off) so no seeded layer is ever exactly zero.
      !!
      !! Degenerate columns (`b <= nz_ml * h_floor` — including dry/negative
      !! land bed) cannot hold `nz_ml` floored layers, so they fall back to the
      !! floored even split, matching `seed_h_layer_uniform_impl`'s wet/dry
      !! branch.  Land columns are masked out downstream either way.
      ! assumed-shape-ok: init routine called once at startup; size(b,1/2) used
      ! to derive loop bounds (flat-impl over registry-dereferenced allocatables).
      real(wp), intent(inout) :: h_layer(:, :, :)
      real(wp), intent(in)    :: b(:, :)  ! assumed-shape-ok: init routine; size(b,1/2) derives loop bounds
      integer, intent(in)     :: nz_ml
      real(wp), intent(in)    :: max_depth
         !! Global basin depth the z interfaces are laid over (`ocean_max_depth`).
      real(wp), intent(in)    :: angstrom_h
         !! Isopycnal minimum-thickness floor (`&ocean_isopycnal_nml angstrom_h`).
      integer  :: i, j, k, nx, ny
      real(wp) :: inv_nz, h_floor, depth, z_bot, z_top, z_top_target

      nx = size(b, 1)
      ny = size(b, 2)
      inv_nz = 1.0_wp/real(nz_ml, wp)
      ! Honour `angstrom_h` all the way down.  The previous
      ! `max(angstrom_h, 2*H_VANISHED)` pinned the seeded floor at 3e-4 m,
      ! which made angstrom_h = 1e-6 and 1e-10 indistinguishable.  That matters:
      ! the vanished layers stack along the bed, so each of their interfaces
      ! carries the full topographic slope and contributes g'*db/dx of spurious
      ! PGF at REST.  The acceleration is independent of h, but the transport it
      ! drives scales WITH h -- which is why MOM6 can tolerate the same
      ! acceleration at ANGSTROM = 1e-10.  Fall back to 2*H_VANISHED only when
      ! the knob is off, so no layer is ever seeded at exactly zero.
      if (angstrom_h > 0.0_wp) then
         h_floor = angstrom_h
      else
         h_floor = 2.0_wp*H_VANISHED
      end if

      ! Plain host loop ON PURPOSE — same reasoning as
      ! `seed_h_layer_uniform_impl`: this runs BEFORE enter_data, so a
      ! `do concurrent` here would make -stdpar=gpu round-trip the
      ! (unmapped) arrays through the device once per loop.
      do j = 1, ny
         do i = 1, nx
            depth = b(i, j)
            if (depth <= real(nz_ml, wp)*h_floor) then
               ! Too shallow (or dry) to hold nz_ml floored layers.
               do k = 1, nz_ml
                  h_layer(i, j, k) = max(depth*inv_nz, h_floor)
               end do
            else
               z_bot = -depth
               do k = 1, nz_ml
                  z_top_target = -max_depth*real(nz_ml - k, wp)*inv_nz
                  ! Never thinner than the floor ...
                  z_top = max(z_top_target, z_bot + h_floor)
                  ! ... and never so thick that the (nz_ml - k) layers still
                  ! above it cannot each clear the floor below z = 0.  Without
                  ! this the sweep can run past the surface whenever
                  ! `max_depth/nz_ml` is not comfortably above `h_floor`, and
                  ! the telescoping sum stops equalling `b`.  The column guard
                  ! (`depth > nz_ml*h_floor`) makes the two bounds consistent
                  ! at every k, and pins `z_top = 0` exactly at k = nz_ml.
                  z_top = min(z_top, -real(nz_ml - k, wp)*h_floor)
                  h_layer(i, j, k) = z_top - z_bot
                  z_bot = z_top
               end do
            end if
         end do
      end do
   end subroutine seed_h_layer_uniform_z_impl