ocean_vcoord_rho_target_column Subroutine

   pure subroutine ocean_vcoord_rho_target_column(i, j, nx, ny, nz, target_h, remap_h_old, &
                                                  total_h, eta, t_conc, s_conc, dsig, &
                                                  floor_dz, rho_target, eos, p_ref, h_min, &
                                                  h_floor_eff, hybrid, floor_mode, h_nominal)
      !! One column of the RHO / HYCOM regrid (steps 0-5 of
      !! `ocean_vcoord_compute_target_h_rho_impl`) — the per-thread body of
      !! `ocean_vcoord_rho_target`.  Same module as its caller (the
      !! project rule for `!$acc routine seq` callees of a `do concurrent`).
      !$acc routine seq
      integer, intent(in), value :: i
         !! Column i-index.
      integer, intent(in), value :: j
         !! Column j-index.
      integer, intent(in), value :: nx
         !! i-extent of every horizontal array (total, incl. halos).
      integer, intent(in), value :: ny
         !! j-extent of every horizontal array (total, incl. halos).
      integer, intent(in), value :: nz
         !! Number of layers; `k = 1` is the bed, `k = nz` the surface.
      real(wp), intent(inout) :: target_h(nx, ny, nz)
         !! Target layer thickness (m), bottom-up.
      real(wp), intent(in) :: remap_h_old(nx, ny, nz)
         !! Pre-remap layer thickness snapshot (m), bottom-up.
      real(wp), intent(in) :: total_h(nx, ny)
         !! Column reference depth H (m).
      real(wp), intent(in) :: eta(nx, ny)
         !! Free-surface anomaly η (m).
      real(wp), intent(in) :: t_conc(nx, ny, nz)
         !! Layer-mean potential temperature (°C), bottom-up.
      real(wp), intent(in) :: s_conc(nx, ny, nz)
         !! Layer-mean salinity (PSU), bottom-up.
      real(wp), intent(in) :: dsig(nz)
         !! Nominal layer fractions, bottom-up (`dsig(nz)` = surface) —
         !! the HYCOM floor only under `HYCOM_FLOOR_SIGMA`.
      real(wp), intent(in) :: floor_dz(nz)
         !! HYCOM z* nominal layer thicknesses (m), bottom-up
         !! (`floor_dz(nz)` = surface) — read only under `HYCOM_FLOOR_PROFILE`.
      real(wp), intent(in) :: rho_target(0:nz)
         !! Target potential densities (kg/m³), `0` = lightest = surface.
      type(eos_t), intent(in) :: eos
         !! Shared EOS handle (flat POD).
      real(wp), intent(in), value :: p_ref
         !! Coordinate reference pressure (Pa) — `rho_ref_pressure`.
      real(wp), intent(in), value :: h_min
         !! Pre-compaction strip threshold (m) — `zstar_h_min`.
      real(wp), intent(in), value :: h_floor_eff
         !! Min-thickness inflation floor (m), `> H_VANISHED`.
      logical, intent(in), value :: hybrid
         !! `.true.` = HYCOM deltas; `.false.` = pure RHO.
      integer, intent(in), value :: floor_mode
         !! HYCOM z* floor source: `HYCOM_FLOOR_PROFILE` (`floor_dz`),
         !! `HYCOM_FLOOR_UNIFORM` (`h_nominal`) or `HYCOM_FLOOR_SIGMA`
         !! (`dsig·H`, the unconfigured-slot fallback).
      real(wp), intent(in), value :: h_nominal
         !! Uniform z* nominal thickness (m) for `HYCOM_FLOOR_UNIFORM`.

      integer :: k, kk, nk, ns, idx_thick, src, ii
      integer :: mapping(NZ_STACK_MAX)
      real(wp) :: h_col(NZ_STACK_MAX), t_col(NZ_STACK_MAX), s_col(NZ_STACK_MAX)
      real(wp) :: hc(NZ_STACK_MAX), rhoc(NZ_STACK_MAX), rtgt(NZ_STACK_MAX)
      real(wp) :: z_new(NZ_STACK_MAX + 1), h_new(NZ_STACK_MAX)
      real(wp) :: col_extent, donate
      real(wp) :: total_need, thick_max, excess, frac
      real(wp) :: nominal_z, stretching, h_ref_col

      ! One column per (j,i).  NZ_STACK_MAX fixed-size locals; no name
      ! shadows a Fortran intrinsic; cross-module pure EOS helper carries
      ! its own `!$acc routine seq`.

      ! --- gather TOP-DOWN: working index 1 = surface = state k=nz ---
      ! Source thicknesses come from the remap snapshot the
      ! orchestrator placed in `remap_h_old` (the live, pre-remap
      ! `h_layer`); T/S are the layer-mean concentrations.  Flip the
      ! bottom-up state index (k=nz surface) into the top-down work
      ! frame (work index 1 = surface).
      col_extent = max(total_h(i, j) + eta(i, j), 0.0_wp)
      do k = 1, nz
         ii = nz - k + 1                 ! state (bottom-up) index
         h_col(k) = remap_h_old(i, j, ii)
         t_col(k) = t_conc(i, j, ii)
         s_col(k) = s_conc(i, j, ii)
      end do

      ! --- too-thin column: cannot carry the coordinate, keep h_old ---
      ! A column thinner than `nz*h_floor_eff` (land / dry columns hold
      ! `nz*H_VANISHED`; a wet/dry or ice-cavity sliver can be any size)
      ! cannot have every layer at the inflation floor.  Step 5 then either
      ! minted mass (`ns == 0`: every layer to the floor) or debited the one
      ! surviving layer below zero (`ns > 0`: a 50 x 1.5e-4 m land column
      ! collapsed into one layer gave 0.0075 - 49*3e-4 = -0.0072 m).  Leave
      ! such a column exactly where it is — the same no-motion answer as the
      ! `nk <= 1` fast path below — so the remap is the identity on it: no
      ! negative thickness, no created mass, `sum(h)` preserved bit-for-bit.
      ! Every column the guard catches was mis-handled by step 5, so it is
      ! inert on every column that step handled correctly.
      if (col_extent < real(nz, wp)*h_floor_eff) then
         do k = 1, nz
            target_h(i, j, k) = remap_h_old(i, j, k)
         end do
         return
      end if

      ! --- step 0: pre-compaction (strip h <= h_min, donate) ---
      nk = 0
      do k = 1, nz
         if (h_col(k) > h_min) then
            nk = nk + 1
            mapping(nk) = k
            hc(nk) = h_col(k)
         end if
      end do
      if (nk <= 1) then
         ! Fast path: <= 1 finite layer.  nz == nk_state here, so
         ! keep the source thicknesses unchanged (h_new = h_old),
         ! flipped back into the bottom-up state.
         do k = 1, nz
            target_h(i, j, k) = remap_h_old(i, j, k)
         end do
         return
      end if

      ! Donate the stripped volume to the thickest survivor.
      donate = col_extent
      do kk = 1, nk
         donate = donate - hc(kk)
      end do
      if (donate > 0.0_wp) then
         idx_thick = 1
         do kk = 2, nk
            if (hc(kk) > hc(idx_thick)) idx_thick = kk
         end do
         hc(idx_thick) = hc(idx_thick) + donate
      end if

      ! --- step 1: layer potential densities on the compacted column ---
      do kk = 1, nk
         src = mapping(kk)
         rhoc(kk) = eos_density_point(eos, t_col(src), s_col(src), p_ref)
      end do

      ! --- step 1b (HYCOM only): bottom-up density monotonize ---
      ! Work frame is top-down (index 1 = surface, nk = bed): cap each
      ! cell by the one below it sweeping bed-up so density is
      ! non-decreasing downward.  Pure RHO (hybrid=.false.) skips this
      ! and stays bit-identical to the merged RHO regrid.
      if (hybrid) then
         do kk = nk - 1, 1, -1
            rhoc(kk) = min(rhoc(kk), rhoc(kk + 1))
         end do
      end if

      ! --- steps 2-4: PPM reconstruct + invert each interior target
      !     density to an interface depth + monotone interfaces.  Shared
      !     density-space inversion (also used by the DENSITY diagnostic
      !     remap, `rdb_ocean_diag_fills`) — single source of truth for
      !     the bracket + fixed-iter-Newton solve.  Copy the interior
      !     targets into a stack array so the device call passes a whole
      !     fixed-size local (no derived-type section descriptor in the
      !     hot per-column kernel).  For HYCOM the rhoc fed in was
      !     monotonized above; the z* floor below then lifts the result.
      do kk = 1, nz - 1
         rtgt(kk) = rho_target(kk)
      end do
      call invert_density_targets(nk, hc, rhoc, nz - 1, rtgt, z_new)

      ! --- step 4b (HYCOM only): z* nominal-floor sweep ---
      ! Surface-side minimum-depth floor on the isopycnal interfaces
      ! (MOM6 `build_hycom1_column`).  Walk interfaces from the surface
      ! down, accumulating the nominal z* depth and pushing any
      ! too-shallow interface DOWN to it (clamped to the column bottom);
      ! deep interfaces already below the floor are untouched.
      ! stretching = col_extent/total_h = (H+η)/H, the z* stretch.  The
      ! nominal increment is a THICKNESS IN METRES — the z* coordinate
      ! resolution (`floor_dz`, or uniform `h_nominal`) — so the band is
      ! the same depth range in every column and a shallow column clamps
      ! its deeper interfaces onto the bed, exactly like a z* grid.  Only
      ! the `HYCOM_FLOOR_SIGMA` fallback accumulates the column FRACTION
      ! `dsig·H` (a sigma floor: k/nz of every column's depth).  Both
      ! tables are bottom-up (index nz = surface layer); the work layer
      ! above interface kk maps to bottom-up index nz-kk+2.  The floor can
      ! break monotonicity, so re-monotonize after it.  Pure RHO
      ! (hybrid=.false.) skips this and is bit-identical.
      if (hybrid) then
         h_ref_col = total_h(i, j)
         if (h_ref_col > 0.0_wp) then
            stretching = col_extent/h_ref_col
         else
            stretching = 1.0_wp
         end if
         nominal_z = 0.0_wp
         do kk = 2, nz + 1
            select case (floor_mode)
            case (HYCOM_FLOOR_PROFILE)
               nominal_z = nominal_z + floor_dz(nz - kk + 2)*stretching
            case (HYCOM_FLOOR_UNIFORM)
               nominal_z = nominal_z + h_nominal*stretching
            case default
               nominal_z = nominal_z + dsig(nz - kk + 2)*h_ref_col*stretching
            end select
            if (z_new(kk) < nominal_z) z_new(kk) = nominal_z
            if (z_new(kk) > col_extent) z_new(kk) = col_extent
         end do
         do kk = 2, nz + 1
            if (z_new(kk) < z_new(kk - 1)) z_new(kk) = z_new(kk - 1)
         end do
      end if
      do kk = 1, nz
         h_new(kk) = z_new(kk + 1) - z_new(kk)
      end do

      ! --- step 5: MOM6 min-thickness inflation (floor h_floor_eff) ---
      ns = 0
      do kk = 1, nz
         if (h_new(kk) > h_floor_eff) ns = ns + 1
      end do
      if (ns == nz) then
         ! all OK
      else if (ns == 0) then
         do kk = 1, nz
            h_new(kk) = h_floor_eff
         end do
      else
         total_need = 0.0_wp
         do kk = 1, nz
            if (h_new(kk) <= h_floor_eff) then
               total_need = total_need + (h_floor_eff - h_new(kk))
               h_new(kk) = h_floor_eff
            end if
         end do
         ! debit the single thickest layer once
         idx_thick = 1
         thick_max = h_new(1)
         do kk = 2, nz
            if (h_new(kk) > thick_max) then
               thick_max = h_new(kk)
               idx_thick = kk
            end if
         end do
         if (thick_max - total_need >= h_floor_eff) then
            h_new(idx_thick) = h_new(idx_thick) - total_need
         else
            ! The thickest layer alone cannot pay (several comparably thin
            ! survivors): debit EVERY above-floor layer in proportion to its
            ! excess over the floor.  `col_extent >= nz*h_floor_eff` (guard
            ! above) makes the total excess >= total_need, so every layer
            ! stays >= h_floor_eff and the sum is unchanged to round-off.
            ! The single-layer debit used to drive the thickest one below
            ! the floor, or negative.
            excess = 0.0_wp
            do kk = 1, nz
               if (h_new(kk) > h_floor_eff) excess = excess + (h_new(kk) - h_floor_eff)
            end do
            if (excess > 0.0_wp) then
               frac = total_need/excess
               do kk = 1, nz
                  if (h_new(kk) > h_floor_eff) then
                     h_new(kk) = h_new(kk) - frac*(h_new(kk) - h_floor_eff)
                  end if
               end do
            end if
         end if
      end if

      ! --- assignment: FLIP top-down working -> bottom-up state ---
      do k = 1, nz
         target_h(i, j, k) = h_new(nz - k + 1)
      end do
   end subroutine ocean_vcoord_rho_target_column