cavity_mass_apply_impl Subroutine

   pure subroutine cavity_mass_apply_impl(nx, ny, nz, dt_over_rho0, inv_rho0_dt, s_ice, &
                                          active, melt, s_far, t_b, &
                                          h_layer, hTr_S, hTr_T, &
                                          salt_budget, heat_budget, k_top, n_thin)
      !! The real-mass top-layer source: add the meltwater VOLUME to
      !! `h_layer` at the first LIVE layer `k_top(i,j)`, replace the
      !! virtual salt flux the assembler
      !! stamped with the real advective salt `w*s_ice`, and add the
      !! enthalpy `dh*T_b` — mirroring both tracer increments into the
      !! existing surface budget contributors.
      !!
      !! See the module docstring for the derivation.  In one line:
      !! `dh = m*dt/rho_0`, `d(hS) = dh*s_ice`, `d(hT) = dh*T_b`, on top
      !! of the `-q_ocean` the assembler already delivers through
      !! `Q_heat`.
      !!
      !! `dt_over_rho0` and `inv_rho0_dt` are the SAME scalar, passed
      !! twice on purpose: the first multiplies `melt`, the second is the
      !! `dt/rho_0` the surface-flux apply used, and the salt undo is
      !! only the EXACT negation of that stamp if the two are the same
      !! bit pattern.  Passing one argument and reusing it makes that
      !! structural rather than a comment.
      !!
      !! FREEZING (`m < 0`) withdraws.  A column whose top layer cannot
      !! give up `|dh|` without falling through `H_CAVITY_FLOOR` is
      !! CLAMPED to what is there above that floor and COUNTED; the
      !! caller treats a non-zero count as fatal
      !! (`cavity_mass_thin_is_fatal`), because a clamped withdrawal no
      !! longer matches the tracked mass source.  The clamp is a plain
      !! `max` on a quantity that has already been range-tested, not a
      !! NaN-laundering `if/else` chain: a non-finite `melt` is caught
      !! upstream by the solver's fatal status.
      !!
      !! The floor is `H_CAVITY_FLOOR`, not `H_VANISHED`, and that gap is
      !! load-bearing: the abort fires AFTER this kernel has written the
      !! state, so whatever it leaves behind gets read at least once more.
      !! A layer left exactly ON the marker reads as VANISHED to every
      !! strict-`>` gate downstream, and `ocean_remap_tracer_column` then
      !! returns `hTr = 0` for it — the heat and salt the clamp was
      !! protecting are deleted by the very next remap.  See
      !! `H_CAVITY_FLOOR`.
      integer, intent(in) :: nx
         !! First dimension (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nz
         !! Layer count; `k = nz` is the surface / ice-base layer.
      real(wp), intent(in) :: dt_over_rho0
         !! `dt/rho_0` (m^3 s / kg) — multiplies `melt` to give `dh`.
      real(wp), intent(in) :: inv_rho0_dt
         !! The same `dt/rho_0`, used to undo the surface-flux stamp.
      real(wp), intent(in) :: s_ice
         !! Ice salinity (g/kg).
      real(wp), intent(in) :: active(nx, ny)
         !! Composed solve mask; 1 only on covered, wet, sampled columns.
      real(wp), intent(in) :: melt(nx, ny)
         !! Melt mass flux (kg/m^2/s), > 0 melting.
      real(wp), intent(in) :: s_far(nx, ny)
         !! Far-field salinity (g/kg) the solve used.  The VIRTUAL salt
         !! component `-melt*(s_far - s_ice)` is REBUILT from it here,
         !! with the same expression `cavity_flux_fill_impl` uses, rather
         !! than read back off `sf`: on a covered column the assembler's
         !! open-water factor is exactly zero, so `Q_salt` IS that
         !! component bit for bit, and rebuilding it keeps this kernel
         !! free of the surface-flux slot.
      real(wp), intent(in) :: t_b(nx, ny)
         !! Interface temperature (degC): the temperature the meltwater
         !! joins the column at.
      real(wp), intent(inout) :: h_layer(nx, ny, nz)
         !! Layer thickness (m); only `k = k_top(i,j)` is touched.
      real(wp), intent(inout) :: hTr_S(nx, ny, nz)
         !! `h*S` ((g/kg) m).
      real(wp), intent(inout) :: hTr_T(nx, ny, nz)
         !! `h*T` (degC m).
      real(wp), intent(inout) :: salt_budget(nx, ny, nz)
         !! `ms%salt_budget_surface`.
      real(wp), intent(inout) :: heat_budget(nx, ny, nz)
         !! `ms%heat_budget_surface`.
      integer, intent(in) :: k_top(nx, ny)
         !! `ms%k_top` — the first LIVE layer counting down from the
         !! top, `nz` wherever nothing vanishes against the top.  Under
         !! a quasi-geopotential coordinate beneath the shelf `k = nz` is
         !! an inert filler on every covered column: growing IT would put
         !! the meltwater where the remap drain deletes it, and shrinking
         !! it would trip the fatal thin-withdrawal clamp on every column
         !! at once (a filler is already AT the marker).
      integer, intent(out) :: n_thin
         !! Columns whose withdrawal had to be clamped.
      integer :: i, j, kt, c_thin
      real(wp) :: dh, h0, hn, cell_s, cell_t, virt

      c_thin = 0
      do concurrent(j=1:ny, i=1:nx) local(kt, dh, h0, hn, cell_s, cell_t, virt) reduce(+:c_thin)
         if (active(i, j) > 0.5_wp) then
            kt = k_top(i, j)
            h0 = h_layer(i, j, kt)
            dh = melt(i, j)*dt_over_rho0
            hn = h0 + dh
            if (dh < 0.0_wp .and. hn < H_CAVITY_FLOOR) then
               ! Pin the RESULT at the floor and derive the applied `dh`
               ! from it, rather than pinning `dh` and adding: `h0 +
               ! (H_CAVITY_FLOOR - h0)` rounds and can land a ulp BELOW
               ! the floor, which is the one place a thin-layer gate must
               ! never be.
               !
               ! `min(h0, ...)` takes only what sits ABOVE the floor: a
               ! column already at or below it has its withdrawal REFUSED
               ! (`hn = h0`, `dh = 0`) rather than turned into a deposit.
               ! Clamping `hn` UP to the floor there would invent mass the
               ! tracked source does not name — a worse failure than the
               ! one being prevented, and it is what the unconditional
               ! `hn = marker` form used to do.
               !
               ! Gated on `dh < 0` for the same reason: a MELTING column
               ! (`dh >= 0`) is adding mass, and a thin result then means
               ! the layer was already thin on arrival — not this kernel's
               ! doing, and not this kernel's to fabricate away.
               hn = min(h0, H_CAVITY_FLOOR)
               dh = hn - h0
               c_thin = c_thin + 1
            end if
            virt = -melt(i, j)*(s_far(i, j) - s_ice)
            cell_s = -inv_rho0_dt*virt + dh*s_ice
            cell_t = dh*t_b(i, j)
            h_layer(i, j, kt) = hn
            hTr_S(i, j, kt) = hTr_S(i, j, kt) + cell_s
            salt_budget(i, j, kt) = salt_budget(i, j, kt) + cell_s
            hTr_T(i, j, kt) = hTr_T(i, j, kt) + cell_t
            heat_budget(i, j, kt) = heat_budget(i, j, kt) + cell_t
         end if
      end do
      n_thin = c_thin
   end subroutine cavity_mass_apply_impl