ocean_vcoord_compute_target_h_rho_impl Subroutine

   pure subroutine ocean_vcoord_compute_target_h_rho_impl(this, total_h, eta, T, S, eos, hybrid)
      !! Isopycnal regrid: place layer interfaces on the prescribed
      !! `rho_target(0:nz)` potential-density surfaces.  ONE
      !! `do concurrent(j,i)` over columns, each running the full
      !! density-space inversion on `NZ_STACK_MAX` fixed-size locals —
      !! no host loop, no per-call allocate.
      !!
      !! Algorithm (clean-room from Bleck 2002 / White & Adcroft 2008;
      !! MOM6 `coord_rho` is the behavioural oracle), per the spec:
      !!
      !!   0. Pre-compaction — strip source layers `h <= H_MIN`, donate
      !!      their volume to the thickest survivor (volume-conserving),
      !!      build the survivor count.  Fast path: `<= 1` survivor →
      !!      `h_new = h_old` (no inversion).
      !!   1. Layer potential densities via `eos_density_point` at
      !!      `rho_ref_pressure` on the compacted column.
      !!   2. PPM (Colella-Woodward monotone) reconstruction of the
      !!      density profile over the compacted thicknesses.
      !!   3. Per interior target, bracket-ordered inversion: light
      !!      boundary → surface; discontinuous-jump sweep; dense
      !!      boundary → bed; else fixed-8-iter Newton on `xi in [0,1]`
      !!      (convergence on `|delta| < NR_TOL` AFTER `xi += delta`;
      !!      zero-gradient `NR_OFFSET` escape at both ends; masked
      !!      fallback to the previous interface on no-bracket).
      !!   4. Monotone non-decreasing interfaces → `h_new`.
      !!   5. MOM6 min-thickness inflation, floor = `max(zstar_h_min,
      !!      H_VANISHED)` (MUST-HAVE #2: keeps RHO-collapsed layers
      !!      above the remap-drain `H_FLOOR` so the next regrid does
      !!      not zero their tracer mass); debit the single thickest
      !!      layer once — or, when that would take it below the floor,
      !!      every above-floor layer in proportion to its excess.
      !!   (guard) A column thinner than `nz·h_floor_eff` cannot hold
      !!      every layer at the floor (land columns hold `nz·H_VANISHED`):
      !!      it keeps `h_new = h_old` (the remap is the identity on it).
      !!      Without the guard step 5 wrote a negative thickness or
      !!      minted mass on every such column (audit finding H3).
      !!
      !! Internal working frame is TOP-DOWN (index 1 = surface, +down),
      !! matching the validated prototype and the `rho_target(0)` =
      !! lightest = surface convention.  The final assignment FLIPS to
      !! the bottom-up state (MUST-HAVE #3): `target_h(:,:,k) =
      !! h_new_td(nz - k + 1)`, so the lightest target lands at the
      !! surface (k=nz) and the densest at the bed (k=1).  Sum is
      !! conserved exactly so `sum_k target_h = H` (η is implicit in
      !! `total_h` here — the caller passes the live column total as
      !! `total_h`, see `ocean_apply_ale_remap_step`).
      !!
      !! HYCOM hybrid (`hybrid = .true.`, Bleck 2002 / MOM6
      !! `build_hycom1_column`): two deltas around the unchanged RHO
      !! inversion, both inside this same column kernel.
      !!   (1b) BOTTOM-UP density monotonize before the PPM reconstruction:
      !!        in the top-down work frame, cap each cell by the one below
      !!        it (toward the bed) — `do k=nk-1,1,-1: rhoc(k)=min(rhoc(k),
      !!        rhoc(k+1))`.  Pure RHO omits this (assumes a monotone
      !!        profile + leans on the PPM limiter); HYCOM enforces it so
      !!        the inversion always sees a monotone column.
      !!   (4b) z* NOMINAL-FLOOR sweep after the inversion, before the
      !!        monotone/inflation: walk interior+bottom interfaces from
      !!        the surface down, accumulating the z* NOMINAL THICKNESS IN
      !!        METRES times `stretching = (H+η)/H`, and push each interface
      !!        DOWN to at least that depth (clamped to the column bottom).
      !!        The nominal thicknesses are the z* coordinate resolution —
      !!        the same `&vcoord_nml z_fixed_profile` table `z_fixed` uses
      !!        (`z_fixed_dz` for "list"/"tanh", `max_depth/nz` for
      !!        "uniform"), as MOM6 HYCOM1 takes its `coordinateResolution`
      !!        from the ALE_COORDINATE_CONFIG that would define a z* grid.
      !!        So the band is a fixed depth range in every column (a 2 m
      !!        surface layer stays 2 m over the shelf and the abyss alike)
      !!        and a shallow column's deeper interfaces clamp onto its bed;
      !!        deep isopycnal interfaces already below the floor are
      !!        untouched.  Until 2026-10-02 the increment was the column
      !!        FRACTION `dsig·(H+η)` with `dsig ≡ 1/nz` — a sigma floor that
      !!        set 95 % of the 1-degree Southern Ocean's interfaces (audit
      !!        finding H1); it survives only as the fallback for a slot no
      !!        setup path configured (`z_fixed_h_ref <= 0`).
      !!        CRITICAL: `stretching` multiplies a nominal thickness that is
      !!        referenced to `H` (= `total_h`), not `H+η`; scaling by
      !!        `(H+η)` twice over-stretches by `(H+η)/H` — invisible at
      !!        η=0, wrong with a free surface.  No renormalize after the
      !!        floor sweep (MOM6 pins the bottom interface + uses the
      !!        debit-thickest inflation instead).
      !! `hybrid = .false.` (the `VCOORD_RHO` path) skips BOTH deltas and
      !! is bit-identical to the P2 kernel.
      !!
      !! Host dispatcher: resolves the scalar knobs and hands every
      !! component to the flat `ocean_vcoord_rho_target` kernel as an
      !! explicit-shape / scalar dummy, so no `this%` reference and no
      !! `associate`-name reaches the `do concurrent` (see that kernel's
      !! docstring for why that is load-bearing on the GPU build).
      type(ocean_vcoord_t), intent(inout) :: this
      ! assumed-shape-ok: cadence-bounded (once per outer ALE step); forwarded
      ! to the explicit-shape kernel below, which is where the loop runs.
      real(wp), intent(in) :: total_h(:, :)
         !! Column reference depth H(i, j) (m).  Caller passes the live
         !! column total (sum of h_layer) so the new grid spans it exactly.
      real(wp), intent(in) :: eta(:, :)  ! assumed-shape-ok: see total_h
         !! Free-surface anomaly η(i, j) (m).  Added to `total_h` to form
         !! the column extent the new interfaces span.
      real(wp), intent(in) :: T(:, :, :)  ! assumed-shape-ok: see total_h
         !! Layer-mean potential temperature concentration (°C), `(nx,ny,nz)`.
      real(wp), intent(in) :: S(:, :, :)  ! assumed-shape-ok: see total_h
         !! Layer-mean salinity concentration (PSU), `(nx,ny,nz)`.
      type(eos_t), intent(in) :: eos
         !! Shared device-resident EOS handle (flat POD, by value).
      logical, intent(in) :: hybrid
         !! `.true.` = HYCOM (apply the monotonize + z*-floor deltas);
         !! `.false.` = pure RHO (bit-identical with the P2 kernel).
      real(wp) :: h_floor_eff, h_nominal
      integer :: floor_mode

      if (.not. this%is_init) return
      ! Inflation floor must be STRICTLY above H_VANISHED: the remap drain
      ! (`ocean_remap_tracer_field`) gates on `h_old > H_FLOOR` (== H_VANISHED)
      ! with a strict `>`, so a layer sitting exactly at H_VANISHED has its
      ! tracer concentration zeroed when this column is fed back as `h_old`
      ! on the next regrid.  Floor at 2·H_VANISHED so inflated layers always
      ! survive the drain (closes the multi-regrid mass-loss footgun).
      h_floor_eff = max(this%zstar_h_min, 2.0_wp*H_VANISHED)
      ! HYCOM z* nominal floor: METRES from the z* coordinate resolution
      ! (the `z_fixed` nominal profile — MOM6 HYCOM1 reads its
      ! `coordinateResolution` from the same ALE_COORDINATE_CONFIG that
      ! sets a z* grid), stretched by (H+η)/H.  A stretched profile
      ! (`z_fixed_use_profile`) gives per-layer `z_fixed_dz`; otherwise the
      ! uniform `z_fixed_h_ref/nz` (setup always writes `max_depth` there).
      ! Only a slot nobody configured (`z_fixed_h_ref <= 0`: unit tests that
      ! build the vcoord by hand) falls back to the historical column-
      ! fraction floor `dsig·(H+η)`.
      if (this%z_fixed_use_profile) then
         floor_mode = HYCOM_FLOOR_PROFILE
         h_nominal = 0.0_wp
      else if (this%z_fixed_h_ref > 0.0_wp) then
         floor_mode = HYCOM_FLOOR_UNIFORM
         h_nominal = this%z_fixed_h_ref/real(this%nz_ml, wp)
      else
         floor_mode = HYCOM_FLOOR_SIGMA
         h_nominal = 0.0_wp
      end if
      call ocean_vcoord_rho_target(this%nx_total, this%ny_total, this%nz_ml, &
                                   this%target_h, this%remap_h_old, total_h, eta, &
                                   T, S, this%dsig, this%z_fixed_dz, this%rho_target, eos, &
                                   this%rho_ref_pressure, this%zstar_h_min, &
                                   h_floor_eff, hybrid, floor_mode, h_nominal)
   end subroutine ocean_vcoord_compute_target_h_rho_impl