ocean_vcoord_z_fixed_target Subroutine

   pure subroutine ocean_vcoord_z_fixed_target(target_h, total_h, eta, z_top, &
                                               nx, ny, nz, h_nominal, use_profile, &
                                               zi, dz_nom, h_min)
      !! `VCOORD_Z_FIXED` target grid — quasi-geopotential interfaces
      !! under a rigid top, with inert fillers and a partial cell at BOTH
      !! ends (Yung, Hallberg, Adcroft & Morrison 2026, JAMES, Fig. 1b:
      !! quasi-z layers are geopotential and VANISH where they outcrop
      !! into the ice base; Asay-Davis et al. 2016 §3.1.5: z-level models
      !! use both partial top and bottom cells).
      !!
      !! ### The rule
      !!
      !! Nominal interface depths are GEOPOTENTIAL and unchanged by the
      !! rigid top: the interface above layer `k` sits at depth
      !! `(nz − k)·h_nominal − η` below `z = 0`, i.e. at
      !! `(nz − k)·h_nominal − z_top` below the column TOP, which is
      !! itself at depth `z_top − η`.  The walk is bed-up (`k = 1` is the
      !! bed) in "depth below the column top", `z_below_loc` tracking the
      !! bottom interface of the layer being laid:
      !!
      !!   * **bed side.** A layer whose nominal top interface is deeper
      !!     than the remaining column collapses to the inert filler
      !!     `h_min` and hands its water UP; the lowest live layer is the
      !!     partial BOTTOM cell and absorbs `η` plus the bed fillers'
      !!     `h_min` debt.  A cut that would leave the partial bottom
      !!     cell at or below `Z_FIXED_BED_PARTIAL_MIN` collapses it too
      !!     — see that constant: no LIVE thickness may land in the
      !!     `(h_min, H_VANISHED]` band.
      !!   * **top side, new.** A layer whose nominal range lies entirely
      !!     above the column top — i.e. inside the ice — collapses to
      !!     `h_min` and stacks immediately under the ice base.  The
      !!     layer whose nominal range STRADDLES the column top is the
      !!     partial TOP cell: its thickness is the part of its nominal
      !!     range below the top, less the top fillers' `h_min` debt.
      !!     If that cut would leave less than
      !!     `max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*h_nominal)`, the sliver
      !!     is merged into the layer BELOW and the index becomes another
      !!     filler (`k_live_top` moves down one) — the mirror of the
      !!     bed's own sliver rule.
      !!
      !! ### Stretched nominal profile
      !!
      !! With `use_profile` the nominal interface above layer `k` sits at
      !! `zi(k)` (bottom-up table, `zi(nz) = 0`) instead of
      !! `(nz − k)·h_nominal`, its bottom interface at `zi(k − 1)`, and the
      !! partial-top threshold is the straddling layer's OWN nominal
      !! thickness, `max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*dz_nom(k))`.
      !! Everything else — the bed rule, the top filler debt, the closing
      !! rule — is shared.  Each branch keeps its whole expression
      !! (`a*b − c` stays one expression in the uniform branch) so an
      !! FMA-contracting build contracts exactly what it did before.
      !!
      !! `Σ_k target_h = total_h + eta` is closed by construction: every
      !! branch assigns exactly what it subtracts from `z_below_loc`, so
      !! each END pays for its own fillers and there is ONE closing rule
      !! that works when both ends vanish.  (A degenerate column thinner
      !! than `nz·h_min` overshoots to `nz·h_min`, exactly as the
      !! pre-cavity code did.)
      !!
      !! ### Bit-identity
      !!
      !! `z_top ≡ 0` ⇒ `k_live_top = nz`, so the top branch is taken ONLY
      !! at `k = nz`, where it evaluates `real(0, wp)*h_min = +0.0` — the
      !! same value as the pre-cavity `real(nz − nz, wp)*h_nominal`.  Every
      !! other layer evaluates `real(nz − k, wp)*h_nominal − 0.0_wp`, which
      !! is bit-for-bit `real(nz − k, wp)*h_nominal` in IEEE round-to-
      !! nearest (and under FMA contraction, since `fma(a, b, −0.0)`
      !! rounds `a·b` once).  Pinned by
      !! `tests/test_ocean_vcoord_zfixed_cavity.F90`.
      integer, intent(in) :: nx
         !! i-extent of every array (total, incl. halos).
      integer, intent(in) :: ny
         !! j-extent of every array (total, incl. halos).
      integer, intent(in) :: nz
         !! Number of layers; `k = 1` is the bed, `k = nz` the top.
      real(wp), intent(out) :: target_h(nx, ny, nz)
         !! Target layer thickness (m).
      real(wp), intent(in) :: total_h(nx, ny)
         !! Column reference thickness `H` (m) — `Σ h_layer − η`.
      real(wp), intent(in) :: eta(nx, ny)
         !! Free-surface anomaly `η` (m).
      real(wp), intent(in) :: z_top(nx, ny)
         !! Geopotential depth of the column top (m, positive down,
         !! `>= 0`).  `0` ⇒ the pre-cavity arithmetic, bit-for-bit.
      real(wp), intent(in) :: h_nominal
         !! Nominal layer spacing `z_fixed_h_ref/nz` (m), `> 0`.  Unused
         !! when `use_profile`.
      logical, intent(in) :: use_profile
         !! Take the nominal interfaces from `zi` / `dz_nom`.
      real(wp), intent(in) :: zi(0:nz)
         !! Nominal interface depths (m), bottom-up, `zi(k)` = top of
         !! layer `k`, `zi(nz) = 0`.  Read only when `use_profile`.
      real(wp), intent(in) :: dz_nom(nz)
         !! Nominal layer thicknesses (m), bottom-up.  Read only when
         !! `use_profile`.
      real(wp), intent(in) :: h_min
         !! Inert-filler thickness (`zstar_h_min`, `<= H_VANISHED`).
      integer :: i, j, k, k_live_top
      real(wp) :: z_below_loc, z_above_nominal_loc, z_top_loc, partial_min
      real(wp) :: bed_partial_min
      logical :: vanish_loc

      partial_min = max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*h_nominal)
      bed_partial_min = max(h_min, Z_FIXED_BED_PARTIAL_MIN)
      do concurrent(j=1:ny, i=1:nx) &
         local(k, k_live_top, z_below_loc, z_above_nominal_loc, z_top_loc, vanish_loc)
         z_top_loc = z_top(i, j)
         ! Index of the shallowest layer the rigid top leaves live: the
         ! largest k whose nominal BOTTOM interface, at depth
         ! `(nz − k + 1)*h_nominal` below z = 0, clears the top by more
         ! than the minimum partial-cell thickness.  Monotone in k, so
         ! the last k that passes is the largest.  Skipped entirely when
         ! there is no rigid top (the bit-identity gate).
         k_live_top = nz
         if (z_top_loc > 0.0_wp) then
            k_live_top = 1
            do k = 1, nz
               if (use_profile) then
                  if (zi(k - 1) - z_top_loc > &
                      max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*dz_nom(k))) then
                     k_live_top = k
                  end if
               else
                  if (real(nz - k + 1, wp)*h_nominal - z_top_loc > partial_min) then
                     k_live_top = k
                  end if
               end if
            end do
         end if
         z_below_loc = total_h(i, j) + eta(i, j)
         do k = 1, nz
            if (k >= k_live_top) then
               ! At and above the first live layer the stack is measured
               ! from the column TOP: `(nz − k)*h_min` is the debt the
               ! top fillers above this layer still owe, so the partial
               ! top cell (k = k_live_top) is cut at the ice base and
               ! pays for them, and every filler above lands on h_min.
               z_above_nominal_loc = real(nz - k, wp)*h_min
               vanish_loc = z_above_nominal_loc > z_below_loc - h_min
            else
               if (use_profile) then
                  z_above_nominal_loc = zi(k) - z_top_loc
               else
                  z_above_nominal_loc = real(nz - k, wp)*h_nominal - z_top_loc
               end if
               ! Bed side: a LIVE partial bottom cell must clear the
               ! REMAP's vanish marker, not merely `h_min`.  `>=` (not
               ! `>`) because the marker itself reads as vanished — the
               ! strict-`>` convention of `H_VANISHED`.
               vanish_loc = z_above_nominal_loc >= z_below_loc - bed_partial_min
            end if
            if (vanish_loc) then
               ! Below the bed / would be a sub-marker sliver — vanish.
               target_h(i, j, k) = h_min
               z_below_loc = z_below_loc - h_min
            else
               ! Layer fits — nominal thickness, or the end residual.
               target_h(i, j, k) = z_below_loc - z_above_nominal_loc
               z_below_loc = z_above_nominal_loc
            end if
         end do
      end do
   end subroutine ocean_vcoord_z_fixed_target