rdb_ocean_cavity_flux.F90

Ice-shelf basal melt, COUPLED: the slot, the far-field sampler and the once-per-thermo-step driver that turns rdb_ocean_cavity_melt’s scalar three-equation kernel into two surface-flux components the ocean integrates.

&ocean_cavity_melt_nml enable (default .false.) is the gate; the geometry it stands on (z_draft, cover_frac, p_ice_ref) is &ocean_cavity_dyn_nml’s and is required. Knob off ⇒ the slot’s arrays stay at their (1,1) placeholders, no kernel is launched, and every path is byte-identical.

What this module owns

ocean_cavity_flux_t — the 2-D device-resident state of the interface: the sampled far field, u*, the interface (T_b, S_b), the canonical melt mass flux, the ocean → interface heat flux, and the per-column solver status.
cavity_far_field_impl — the far-field sample, in METRES below the ice base, never “layer nz”.
ocean_cavity_flux_step — the driver: sample, solve (through cavity_melt_columns_2d, which must own the do concurrent; see its docstring for the nvlink reason), deliver, count, fail loud.

Far-field sampling (why metres)

The melt rate a three-equation law returns is roughly linear in the thermal driving it is handed, and the thermal driving depends on HOW FAR FROM THE ICE the model sampled it. That is the dominant resolution artefact in the subject, not a refinement: Gwyther et al. (2020), Burchard et al. (2022) Table 2 p. 15 (the all-bulk melt-rate error GROWS under refinement — -8 % at 6.66 m, -30 % at 0.015 m) and Yung et al. (2026) p. 2074 all identify it. Sampling “the top layer” would therefore make the melt rate a function of the vertical coordinate’s cell thickness — which is precisely the quantity a coordinate study must hold fixed. So far_field_depth is in metres and the sampler is a thickness-weighted mean over the layers that span it, with a PARTIAL last layer.

Vanished layers (h <= H_VANISHED, the D4 skip/merge marker) carry no mass and are skipped rather than clamped. A column whose sample finds no mass at all is dropped from the cover mask — zero melt, OK status — not solved on a fabricated state.

The three pressures, and which one this is

The liquidus is evaluated at multilayer_state_t%p_top — THE interface pressure, the same field the FV_MOM6 pressure-stack boundary condition and the in-situ EOS read (src/core/ocean/README.md, the p_top seam contract). It is NOT eos%p_ref (a scalar potential -density reference, deliberately horizontally uniform) and NOT the eta_forcing seam (a gradient). One pressure, three consumers, and this module is a CONSUMER — it writes nothing to p_top. The sole producer is configure_ocean_cavity, which assembles p_top = metrics%p_ice_ref + sf%p_surf (re-assembled per outer step in ocean_dyn_step_split when the psurf seam makes sf%p_surf live); configure_ocean_cavity_melt asserts that it ran rather than seeding a second copy, because a second writer running after it would silently drop sf%p_surf. A consequence worth stating: an atmospheric load under the shelf depresses the freezing point exactly as the ice load does, with no extra wiring.

Sign and units — the table that decides whether the answer is right

quantity	symbol	unit	sign
melt	`melt` (`m_mass`)	kg/m^2/s	> 0 melting
kernel heat flux	`q_ocean`	W/m^2	> 0 warms the INTERFACE (the ocean cools)
delivered heat	`sf%heat_cavity`	W/m^2	positive DOWN into the ocean
delivered salt	`sf%salt_cavity`	same as `sf%salt_flux`	positive SALINIFIES

Heat. q_ocean = rho_w*c_w*gamma_t*(T_w - T_b) > 0 is the turbulent flux the interface takes FROM the ocean, so the component the ocean must be given is its negative:

heat_cavity = -q_ocean

The assembler adds it to Q_heat with a plain +, and apply_surface_src_2d_impl deposits dt/(rho0*cp) * Q_heat at k = nz — positive warms. Warm water under a shelf therefore COOLS the top of the column, which is the whole point.

Salt — and why it is not the same construction. Melting adds freshwater MASS the Boussinesq column does not yet carry (that is Phase 3). Emulate the dilution with a salt flux and the fixed-mass equivalent is exact, not approximate. For a column of mass M per unit area receiving mass m at salinity S_i:

d(M*S)/dt = m*S_i ,   dM/dt = m

=> MdS/dt = mS_i - Sm = -m(S - S_i)

so the virtual flux that reproduces dS/dt at FIXED M is

salt_cavity = -m_mass*(S_far - s_ice)

referenced to the SAME far-field salinity the solve was handed. Melting (m_mass > 0, S_far > s_ice) freshens.

The heat twin of that -S*dM/dt dilution term is -m*c_w*(T_w - T_b) and it is DELIBERATELY DROPPED. The asymmetry is physical, not sloppiness: S_far - s_ice ~ 34 g/kg is O(1) of the salinity itself, while T_w - T_b is a few hundredths of a degree, so the dropped heat term is m/(rho_w*gamma_t) ~ 1e-3 of q_ocean — a 0.1 % correction that belongs with the real-mass PR, whereas the salt term is the entire meltwater buoyancy signal and cannot wait for it.

NOTE the prototype’s column_demo.py removes f_turb = rho_w*gamma_s* (S_w - S_b) = m_mass*(S_b - s_ice) instead. That is the salt flux ACROSS the interface, not the dilution equivalent — it omits the meltwater that in reality returns to the column carrying S_b, and so under-freshens by (S_b - s_ice)/(S_far - s_ice). That demo is explicitly not part of the oracle (its README: “for reading”); the identity above is what the ocean’s fixed-mass salinity equation requires, and it is what the budget test asserts.

Phase 3 — REAL freshwater MASS (`&ocean_cavity_melt_nml freshwater`)

freshwater = "virtual" is the default and everything above is what runs. freshwater = "mass" moves the meltwater as a real Boussinesq VOLUME instead. The derivation, in full, because every sign and every omitted term below is load-bearing.

Which density converts kg to m

m_mass is a MASS flux (kg m^-2 s^-1). The model is BOUSSINESQ: it conserves VOLUME, h_layer is a thickness, and every extensive total the console prints is a volume integral scaled by a constant density. The one density that is allowed to convert a mass flux into a thickness tendency is therefore the Boussinesq reference rho_0 — the SAME rho_0 the surface-flux apply already divides by (dS/dt = Q_salt/(rho_0 h)), and the same one the EOS carries. Using the meltwater’s own freshwater density (~1000) instead would make a kilogram of meltwater occupy more volume than a kilogram of the sea water it joins, which in a Boussinesq model is not a refinement but an inconsistency: the column’s mass total is rho_0 * volume by construction, so a second density would put the tracked source and the tracked total on different scales. So

w  ==  m_mass / rho_0          [m/s, > 0 melting]
dh ==  w * dt                  added to `h_layer(:,:,nz)`

(The console’s total_mass = sum(h*areaT)*RHO_WATER uses the MODULE constant RHO_WATER, which a run may configure rho_0 away from. The tracked mass source is therefore RHO_WATER * sum(dh*areaT), i.e. (RHO_WATER/rho_0) * sum(m*area*dt) — exactly sum(m*area*dt) at the default rho_0 = RHO_WATER. The accumulator must measure the same mass the total does; it is not free to pick its own density.)

The three top-layer tendencies

Take the ocean column as the control volume. Its top boundary is the ice base, and the only things that cross it are the meltwater (mass m, salinity s_ice, temperature T_b) and the turbulent fluxes.

Mass. d h/dt = w.

Salt. Across the ocean’s top boundary, salt enters ADVECTIVELY with the water leaving the interface (m*S_b) and leaves DIFFUSIVELY into the interface (rho_w*gamma_S*(S_w - S_b)). The interface itself stores nothing, so its own balance — the third of the three equations — is

m*s_ice + rho_w*gamma_S*(S_w - S_b) = m*S_b

and the NET salt the ocean gains is

m*S_b - rho_w*gamma_S*(S_w - S_b) = m*S_b - m*(S_b - s_ice)
                                  = m*s_ice

— the ice is the only thing on the other side of the boundary, and it carries s_ice per kilogram. The advective and diffusive halves cancel EXACTLY; there is nothing to approximate. So

d(h*S)/dt = w*s_ice

and for the ISOMIP+ s_ice = 0 the salt content of the column does not change at all — the salinity falls purely because h grows.

Heat. The same control volume loses the turbulent flux q_ocean and gains the meltwater’s enthalpy. The model’s heat content is rho_0*c_p*integral(T dz) with an implicit 0 degC reference (total_heat = sum(h*T*areaT)*RHO_WATER carries no offset), so in hTr = h*T space the enthalpy of the added water is simply dh*T_b: NO heat capacity enters, and the model’s c_p and the melt law’s own c_w never have to agree.

d(h*T)/dt = -q_ocean/(rho_0*c_p) + w*T_b

Why the virtual salt flux must NOT also be applied — and why it

is nevertheless still assembled

Applying both would double-count: the dilution would happen once through the growing h and once through the removed salt. But the virtual flux is ALSO the entire meltwater buoyancy signal that KPP and EPBL read out of Q_salt to build B_0 (beta*Q_salt/rho_0 is ~6x the alpha*Q_heat/(rho_0*c_p) term at these temperatures), and the real volume source does not reach B_0 at all. Dropping it from Q_salt would leave both boundary-layer schemes blind to the dominant, destabilising half of the surface buoyancy flux.

So under "mass" the assembled Q_salt is UNCHANGED — salt_cavity keeps carrying -m*(S_far - s_ice) and B_0 is the same under both forms, which is the point — and ocean_cavity_mass_step removes that increment again from the TRACER in the same stage, adding the real w*s_ice in its place. The undo is the exact negation of what apply_surface_src_2d_impl stamped (same dt/rho_0, same field, and cover_frac = 1 makes the atmospheric group exactly zero on a covered column), so x + (-x) is zero in floating point and the salinity tracer load is left EXACTLY unchanged when s_ice = 0. Both increments are mirrored into salt_budget_surface in the same order, so the budget sees the identical pair and closes to round-off. The pseudo-salt mirror receives the same correction, so it keeps receiving “exactly salinity’s surface salt flux” — its documented purpose — under either form.

Passive tracers need NOTHING: h grows, h*C does not, so a tracer with zero concentration in the meltwater is diluted automatically. The ideal-age tracer is diluted the same way, which is the right answer (meltwater is new water).

Agreement with the virtual form, and the size of the difference

With eps = w*dt/h0 the two forms differ EXACTLY by

dS_mass - dS_virtual = +(S0 - s_ice) * eps^2/(1 + eps)
dT_mass - dT_virtual = -eps*(A + T0 - T_b)/(1 + eps),
   A = -q_ocean*dt/(rho_0*c_p*h0) = dT_virtual

so the salinity difference is second order in eps and the temperature difference is eps times the (tiny) dropped dilution term plus eps times the virtual increment itself. test_ocean_ cavity_freshwater asserts those two identities rather than asserting “the two agree”, because agreeing to zero is what a double count would also do.

Sea level (`volume_compensation`)

Real mass into a CLOSED domain raises it. ISOMIP+ Sect. 3.1.3 states the protocol’s position: the freshwater flux “is not compensated by a corresponding removal of water elsewhere in the domain” for the restored Ocean0-2 configurations — whose only open seam is the northern restoring band, which relaxes T and S and moves no volume — whereas the closed Ocean3/4 boxes may compensate. Scale, for Ocean0: ~30 m/yr of melt over ~1e10 m^2 of shelf into ~4e10 m^2 of surface is metres per year of sea level. Hence the knob.

volume_compensation = "uniform_open_ocean" removes the domain-integrated melt volume again each thermo step, spread uniformly PER UNIT AREA over the wet cells the ice does not cover, each parcel carrying that cell’s own T and S so no concentration there changes. It is tracked as a sink in all three budgets (a parcel that carries S does change the domain SALT total, so the salt and heat budgets must name it too).

Budgets

Mass gains a tracked SOURCE (ms%mass_src), accumulated with the same per-stage weight ms%mass_out uses, so the console residual (M - M0) + mass_out - mass_src stays at round-off while the total legitimately grows. Salt and heat need no new accumulator: every increment this path makes is mirrored into the EXISTING salt_budget_surface / heat_budget_surface contributors, by the same in-stage kernel and at the same cadence as the surface fluxes, so they carry the right ocean_budget_stage_weight by construction.

Budgets

Both components ride Q_heat/Q_salt, so they are integrated by ocean_surface_flux_apply_tracers and land in the EXISTING ms%heat_budget_surface / ms%salt_budget_surface contributors — which the console already folds in with the correct ocean_budget_stage_weight. No separate frazil-style accumulator is needed, and no full-weight/half-weight decision is taken here: the source is applied by the same per-stage kernel as every other surface flux, so it carries the same weight by construction.

`mem:separate`

Every array this module touches is device-resident and mapped by the slot’s enter_data. Nothing is copied host↔device per step except the four status COUNTS, which come back as do concurrent ... reduce results (host scalars), exactly like continuity_t%n_limited_step.

Citations (papers, never another model’s source)

Holland, D. M. and Jenkins, A. (1999): J. Phys. Oceanogr. 29, 1787-1800.
Jenkins, A., Nicholls, K. W. and Corr, H. F. J. (2010): J. Phys. Oceanogr. 40, 2298-2312.
Asay-Davis, X. S. et al. (2016): Geosci. Model Dev. 9, 2471-2497 (ISOMIP+).
Burchard, H. et al. (2022): Ocean Modelling 179, 102119.
Yung, C. K. et al. (2025): The Cryosphere 19, 5827-5861.
Yung, C. K. et al. (2026): The Cryosphere 20, 2053-2088.

Files dependent on this one

Help

Source Code

!! Ice-shelf basal melt, COUPLED: the slot, the far-field sampler and the
!! once-per-thermo-step driver that turns `rdb_ocean_cavity_melt`'s scalar
!! three-equation kernel into two surface-flux components the ocean
!! integrates.
!!
!! `&ocean_cavity_melt_nml enable` (default `.false.`) is the gate; the
!! geometry it stands on (`z_draft`, `cover_frac`, `p_ice_ref`) is
!! `&ocean_cavity_dyn_nml`'s and is required.  Knob off ⇒ the slot's
!! arrays stay at their `(1,1)` placeholders, no kernel is launched, and
!! every path is byte-identical.
!!
!! ## What this module owns
!!
!! 1. `ocean_cavity_flux_t` — the 2-D device-resident state of the
!!    interface: the sampled far field, `u*`, the interface `(T_b, S_b)`,
!!    the canonical melt mass flux, the ocean → interface heat flux, and
!!    the per-column solver status.
!! 2. `cavity_far_field_impl` — the far-field sample, in METRES below the
!!    ice base, never "layer `nz`".
!! 3. `ocean_cavity_flux_step` — the driver: sample, solve (through
!!    `cavity_melt_columns_2d`, which must own the `do concurrent`; see
!!    its docstring for the nvlink reason), deliver, count, fail loud.
!!
!! ## Far-field sampling (why metres)
!!
!! The melt rate a three-equation law returns is roughly linear in the
!! thermal driving it is handed, and the thermal driving depends on HOW
!! FAR FROM THE ICE the model sampled it.  That is the dominant
!! resolution artefact in the subject, not a refinement: Gwyther et al.
!! (2020), Burchard et al. (2022) Table 2 p. 15 (the all-bulk melt-rate
!! error GROWS under refinement — -8 % at 6.66 m, -30 % at 0.015 m) and
!! Yung et al. (2026) p. 2074 all identify it.  Sampling "the top layer"
!! would therefore make the melt rate a function of the vertical
!! coordinate's cell thickness — which is precisely the quantity a
!! coordinate study must hold fixed.  So `far_field_depth` is in metres
!! and the sampler is a thickness-weighted mean over the layers that
!! span it, with a PARTIAL last layer.
!!
!! Vanished layers (`h <= H_VANISHED`, the D4 skip/merge marker) carry no
!! mass and are skipped rather than clamped.  A column whose sample finds
!! no mass at all is dropped from the cover mask — zero melt, `OK`
!! status — not solved on a fabricated state.
!!
!! ## The three pressures, and which one this is
!!
!! The liquidus is evaluated at `multilayer_state_t%p_top` — THE
!! interface pressure, the same field the FV_MOM6 pressure-stack
!! boundary condition and the in-situ EOS read (`src/core/ocean/README.md`,
!! the `p_top` seam contract).  It is NOT `eos%p_ref` (a scalar potential
!! -density reference, deliberately horizontally uniform) and NOT the
!! `eta_forcing` seam (a gradient).  One pressure, three consumers, and
!! this module is a CONSUMER — it writes nothing to `p_top`.  The sole
!! producer is `configure_ocean_cavity`, which assembles
!! `p_top = metrics%p_ice_ref + sf%p_surf` (re-assembled per outer step
!! in `ocean_dyn_step_split` when the psurf seam makes `sf%p_surf`
!! live); `configure_ocean_cavity_melt` asserts that it ran rather than
!! seeding a second copy, because a second writer running after it would
!! silently drop `sf%p_surf`.  A consequence worth stating: an
!! atmospheric load under the shelf depresses the freezing point exactly
!! as the ice load does, with no extra wiring.
!!
!! ## Sign and units — the table that decides whether the answer is right
!!
!! | quantity | symbol | unit | sign |
!! |---|---|---|---|
!! | melt | `melt` (`m_mass`) | kg/m^2/s | **> 0 melting** |
!! | kernel heat flux | `q_ocean` | W/m^2 | **> 0 warms the INTERFACE** (the ocean cools) |
!! | delivered heat | `sf%heat_cavity` | W/m^2 | positive DOWN into the ocean |
!! | delivered salt | `sf%salt_cavity` | same as `sf%salt_flux` | positive SALINIFIES |
!!
!! **Heat.**  `q_ocean = rho_w*c_w*gamma_t*(T_w - T_b) > 0` is the
!! turbulent flux the interface takes FROM the ocean, so the component
!! the ocean must be given is its negative:
!!
!!     heat_cavity = -q_ocean
!!
!! The assembler adds it to `Q_heat` with a plain `+`, and
!! `apply_surface_src_2d_impl` deposits `dt/(rho0*cp) * Q_heat` at
!! `k = nz` — positive warms.  Warm water under a shelf therefore COOLS
!! the top of the column, which is the whole point.
!!
!! **Salt — and why it is not the same construction.**  Melting adds
!! freshwater MASS the Boussinesq column does not yet carry (that is
!! Phase 3).  Emulate the dilution with a salt flux and the fixed-mass
!! equivalent is exact, not approximate.  For a column of mass `M` per
!! unit area receiving mass `m` at salinity `S_i`:
!!
!!     d(M*S)/dt = m*S_i ,   dM/dt = m
!!  => M*dS/dt = m*S_i - S*m = -m*(S - S_i)
!!
!! so the virtual flux that reproduces `dS/dt` at FIXED `M` is
!!
!!     salt_cavity = -m_mass*(S_far - s_ice)
!!
!! referenced to the SAME far-field salinity the solve was handed.
!! Melting (`m_mass > 0`, `S_far > s_ice`) freshens.
!!
!! The heat twin of that `-S*dM/dt` dilution term is `-m*c_w*(T_w - T_b)`
!! and it is DELIBERATELY DROPPED.  The asymmetry is physical, not
!! sloppiness: `S_far - s_ice ~ 34 g/kg` is O(1) of the salinity itself,
!! while `T_w - T_b` is a few hundredths of a degree, so the dropped heat
!! term is `m/(rho_w*gamma_t) ~ 1e-3` of `q_ocean` — a 0.1 % correction
!! that belongs with the real-mass PR, whereas the salt term is the
!! entire meltwater buoyancy signal and cannot wait for it.
!!
!! NOTE the prototype's `column_demo.py` removes `f_turb = rho_w*gamma_s*
!! (S_w - S_b) = m_mass*(S_b - s_ice)` instead.  That is the salt flux
!! ACROSS the interface, not the dilution equivalent — it omits the
!! meltwater that in reality returns to the column carrying `S_b`, and so
!! under-freshens by `(S_b - s_ice)/(S_far - s_ice)`.  That demo is
!! explicitly not part of the oracle (its README: "for reading"); the
!! identity above is what the ocean's fixed-mass salinity equation
!! requires, and it is what the budget test asserts.
!!
!! ## Phase 3 — REAL freshwater MASS (`&ocean_cavity_melt_nml freshwater`)
!!
!! `freshwater = "virtual"` is the default and everything above is what
!! runs.  `freshwater = "mass"` moves the meltwater as a real Boussinesq
!! VOLUME instead.  The derivation, in full, because every sign and
!! every omitted term below is load-bearing.
!!
!! ### Which density converts kg to m
!!
!! `m_mass` is a MASS flux (kg m^-2 s^-1).  The model is BOUSSINESQ: it
!! conserves VOLUME, `h_layer` is a thickness, and every extensive total
!! the console prints is a volume integral scaled by a constant density.
!! The one density that is allowed to convert a mass flux into a
!! thickness tendency is therefore the Boussinesq reference `rho_0` —
!! the SAME `rho_0` the surface-flux apply already divides by
!! (`dS/dt = Q_salt/(rho_0 h)`), and the same one the EOS carries.
!! Using the meltwater's own freshwater density (~1000) instead would
!! make a kilogram of meltwater occupy more volume than a kilogram of
!! the sea water it joins, which in a Boussinesq model is not a
!! refinement but an inconsistency: the column's mass total is
!! `rho_0 * volume` by construction, so a second density would put the
!! tracked source and the tracked total on different scales.  So
!!
!!     w  ==  m_mass / rho_0          [m/s, > 0 melting]
!!     dh ==  w * dt                  added to `h_layer(:,:,nz)`
!!
!! (The console's `total_mass = sum(h*areaT)*RHO_WATER` uses the MODULE
!! constant `RHO_WATER`, which a run may configure `rho_0` away from.
!! The tracked mass source is therefore `RHO_WATER * sum(dh*areaT)`, i.e.
!! `(RHO_WATER/rho_0) * sum(m*area*dt)` — exactly `sum(m*area*dt)` at the
!! default `rho_0 = RHO_WATER`.  The accumulator must measure the same
!! mass the total does; it is not free to pick its own density.)
!!
!! ### The three top-layer tendencies
!!
!! Take the ocean column as the control volume.  Its top boundary is the
!! ice base, and the only things that cross it are the meltwater (mass
!! `m`, salinity `s_ice`, temperature `T_b`) and the turbulent fluxes.
!!
!! **Mass.**  `d h/dt = w`.
!!
!! **Salt.**  Across the ocean's top boundary, salt enters ADVECTIVELY
!! with the water leaving the interface (`m*S_b`) and leaves
!! DIFFUSIVELY into the interface (`rho_w*gamma_S*(S_w - S_b)`).  The
!! interface itself stores nothing, so its own balance — the third of
!! the three equations — is
!!
!!     m*s_ice + rho_w*gamma_S*(S_w - S_b) = m*S_b
!!
!! and the NET salt the ocean gains is
!!
!!     m*S_b - rho_w*gamma_S*(S_w - S_b) = m*S_b - m*(S_b - s_ice)
!!                                       = m*s_ice
!!
!! — the ice is the only thing on the other side of the boundary, and it
!! carries `s_ice` per kilogram.  The advective and diffusive halves
!! cancel EXACTLY; there is nothing to approximate.  So
!!
!!     d(h*S)/dt = w*s_ice
!!
!! and for the ISOMIP+ `s_ice = 0` the salt content of the column does
!! not change at all — the salinity falls purely because `h` grows.
!!
!! **Heat.**  The same control volume loses the turbulent flux
!! `q_ocean` and gains the meltwater's enthalpy.  The model's heat
!! content is `rho_0*c_p*integral(T dz)` with an implicit 0 degC
!! reference (`total_heat = sum(h*T*areaT)*RHO_WATER` carries no
!! offset), so in `hTr = h*T` space the enthalpy of the added water is
!! simply `dh*T_b`: NO heat capacity enters, and the model's `c_p` and
!! the melt law's own `c_w` never have to agree.
!!
!!     d(h*T)/dt = -q_ocean/(rho_0*c_p) + w*T_b
!!
!! ### Why the virtual salt flux must NOT also be applied — and why it
!! ### is nevertheless still assembled
!!
!! Applying both would double-count: the dilution would happen once
!! through the growing `h` and once through the removed salt.  But the
!! virtual flux is ALSO the entire meltwater buoyancy signal that KPP
!! and EPBL read out of `Q_salt` to build `B_0` (`beta*Q_salt/rho_0` is
!! ~6x the `alpha*Q_heat/(rho_0*c_p)` term at these temperatures), and
!! the real volume source does not reach `B_0` at all.  Dropping it from
!! `Q_salt` would leave both boundary-layer schemes blind to the
!! dominant, destabilising half of the surface buoyancy flux.
!!
!! So under `"mass"` the assembled `Q_salt` is UNCHANGED — `salt_cavity`
!! keeps carrying `-m*(S_far - s_ice)` and `B_0` is the same under both
!! forms, which is the point — and `ocean_cavity_mass_step` removes that
!! increment again from the TRACER in the same stage, adding the real
!! `w*s_ice` in its place.  The undo is the exact negation of what
!! `apply_surface_src_2d_impl` stamped (same `dt/rho_0`, same field, and
!! `cover_frac = 1` makes the atmospheric group exactly zero on a
!! covered column), so `x + (-x)` is zero in floating point and the
!! salinity tracer load is left EXACTLY unchanged when `s_ice = 0`.
!! Both increments are mirrored into `salt_budget_surface` in the same
!! order, so the budget sees the identical pair and closes to round-off.
!! The pseudo-salt mirror receives the same correction, so it keeps
!! receiving "exactly salinity's surface salt flux" — its documented
!! purpose — under either form.
!!
!! Passive tracers need NOTHING: `h` grows, `h*C` does not, so a tracer
!! with zero concentration in the meltwater is diluted automatically.
!! The ideal-age tracer is diluted the same way, which is the right
!! answer (meltwater is new water).
!!
!! ### Agreement with the virtual form, and the size of the difference
!!
!! With `eps = w*dt/h0` the two forms differ EXACTLY by
!!
!!     dS_mass - dS_virtual = +(S0 - s_ice) * eps^2/(1 + eps)
!!     dT_mass - dT_virtual = -eps*(A + T0 - T_b)/(1 + eps),
!!        A = -q_ocean*dt/(rho_0*c_p*h0) = dT_virtual
!!
!! so the salinity difference is second order in `eps` and the
!! temperature difference is `eps` times the (tiny) dropped dilution
!! term plus `eps` times the virtual increment itself.  `test_ocean_
!! cavity_freshwater` asserts those two identities rather than asserting
!! "the two agree", because agreeing to zero is what a double count
!! would also do.
!!
!! ### Sea level (`volume_compensation`)
!!
!! Real mass into a CLOSED domain raises it.  ISOMIP+ Sect. 3.1.3 states
!! the protocol's position: the freshwater flux "is not compensated by a
!! corresponding removal of water elsewhere in the domain" for the
!! restored Ocean0-2 configurations — whose only open seam is the
!! northern restoring band, which relaxes T and S and moves no volume —
!! whereas the closed Ocean3/4 boxes may compensate.  Scale, for Ocean0:
!! ~30 m/yr of melt over ~1e10 m^2 of shelf into ~4e10 m^2 of surface is
!! metres per year of sea level.  Hence the knob.
!!
!! `volume_compensation = "uniform_open_ocean"` removes the
!! domain-integrated melt volume again each thermo step, spread
!! uniformly PER UNIT AREA over the wet cells the ice does not cover,
!! each parcel carrying that cell's own `T` and `S` so no concentration
!! there changes.  It is tracked as a sink in all three budgets (a
!! parcel that carries `S` does change the domain SALT total, so the
!! salt and heat budgets must name it too).
!!
!! ### Budgets
!!
!! Mass gains a tracked SOURCE (`ms%mass_src`), accumulated with the same
!! per-stage weight `ms%mass_out` uses, so the console residual
!! `(M - M0) + mass_out - mass_src` stays at round-off while the total
!! legitimately grows.  Salt and heat need no new accumulator: every
!! increment this path makes is mirrored into the EXISTING
!! `salt_budget_surface` / `heat_budget_surface` contributors, by the
!! same in-stage kernel and at the same cadence as the surface fluxes,
!! so they carry the right `ocean_budget_stage_weight` by construction.
!!
!! ## Budgets
!!
!! Both components ride `Q_heat`/`Q_salt`, so they are integrated by
!! `ocean_surface_flux_apply_tracers` and land in the EXISTING
!! `ms%heat_budget_surface` / `ms%salt_budget_surface` contributors —
!! which the console already folds in with the correct
!! `ocean_budget_stage_weight`.  No separate frazil-style accumulator is
!! needed, and no full-weight/half-weight decision is taken here: the
!! source is applied by the same per-stage kernel as every other surface
!! flux, so it carries the same weight by construction.
!!
!! ## `mem:separate`
!!
!! Every array this module touches is device-resident and mapped by the
!! slot's `enter_data`.  Nothing is copied host↔device per step except
!! the four status COUNTS, which come back as `do concurrent ... reduce`
!! results (host scalars), exactly like `continuity_t%n_limited_step`.
!!
!! ## Citations (papers, never another model's source)
!!
!!   * Holland, D. M. and Jenkins, A. (1999): J. Phys. Oceanogr. 29,
!!     1787-1800.
!!   * Jenkins, A., Nicholls, K. W. and Corr, H. F. J. (2010): J. Phys.
!!     Oceanogr. 40, 2298-2312.
!!   * Asay-Davis, X. S. et al. (2016): Geosci. Model Dev. 9, 2471-2497
!!     (ISOMIP+).
!!   * Burchard, H. et al. (2022): Ocean Modelling 179, 102119.
!!   * Yung, C. K. et al. (2025): The Cryosphere 19, 5827-5861.
!!   * Yung, C. K. et al. (2026): The Cryosphere 20, 2053-2088.
module rdb_ocean_cavity_flux
   use, intrinsic :: iso_fortran_env, only: int64
   use rdb_constants, only: wp, H_VANISHED, H_DIV_EPS, RHO_WATER
   use rdb_grid, only: hgrid_t
   use rdb_eos, only: eos_t
   use rdb_mem_report, only: arr_bytes
   use rdb_multilayer_state, only: multilayer_state_t
   use rdb_ocean_metrics, only: ocean_metrics_t
   use rdb_ocean_surface_flux, only: ocean_surface_flux_t
   use rdb_ocean_cavity_melt, only: ocean_cavity_const_t, ocean_cavity_exchange_t, &
                                    ocean_cavity_ice_t, cavity_melt_columns_2d, &
                                    CAVITY_MELT_OK, CAVITY_MELT_NONFINITE_INPUT, &
                                    CAVITY_MELT_NONFINITE_STATE, &
                                    CAVITY_MELT_NOT_CONVERGED, &
                                    CAVITY_MELT_NO_PHYSICAL_ROOT, &
                                    CAVITY_FW_VIRTUAL, CAVITY_FW_MASS, &
                                    CAVITY_VC_NONE, CAVITY_VC_UNIFORM_OPEN
   implicit none
   private

   public :: ocean_cavity_flux_t
   public :: ocean_cavity_flux_step
   public :: ocean_cavity_mass_step
   public :: cavity_far_field_impl
   public :: cavity_flux_fill_impl
   public :: cavity_status_counts_impl
   public :: cavity_melt_status_is_fatal
   public :: cavity_mass_apply_impl
   public :: cavity_mass_salt_mirror_impl
   public :: cavity_mass_totals_impl
   public :: cavity_comp_apply_impl
   public :: cavity_comp_scale_tracer_impl
   public :: cavity_mass_thin_is_fatal
   public :: cavity_comp_withdrawal
   public :: H_CAVITY_FLOOR

   real(wp), parameter :: H_CAVITY_FLOOR = 2.0_wp*H_VANISHED
      !! Thickness (m) a clamped top-layer WITHDRAWAL is allowed to leave
      !! behind — one full marker ABOVE `H_VANISHED`, not on it.
      !!
      !! `H_VANISHED` is the D4 skip/merge MARKER, and every gate in the
      !! tree tests it with a strict `>`: a layer sitting exactly ON the
      !! marker reads as VANISHED.  So clamping a melt or compensation
      !! withdrawal to `H_VANISHED` exactly — which is what these kernels
      !! used to do, with an in-line comment arguing for landing on the
      !! marker — hands the next ALE remap a layer it treats as empty:
      !! `ocean_remap_tracer_column`'s `c = hTr/h` guard returns 0 there,
      !! so the layer's heat and salt content is DELETED, silently and at
      !! the next thermo step.  One ulp above the marker it is preserved.
      !! A clamped withdrawal is already a fail-loud condition
      !! (`cavity_mass_thin_is_fatal`), but the abort happens AFTER the
      !! state is written, so the value written has to survive being read.
      !!
      !! `2*H_VANISHED` is the convention the rest of the tree already
      !! uses where a producer must land ABOVE the marker on purpose —
      !! `compute_target_h_rho_impl`'s `h_floor_eff = max(zstar_h_min,
      !! 2*H_VANISHED)` and `seed_land_h_floor`'s `max(angstrom_h,
      !! 2*H_VANISHED)`.  Written as a multiple of the constant of record
      !! rather than as a literal, so it cannot drift from it.

   type :: ocean_cavity_flux_t
      !! Ice-shelf basal-melt slot: the 2-D interface state plus the
      !! parameter bundles the kernel is called with.
      !!
      !! Every array is `(nx_total, ny_total)` when `enable`, `(1,1)`
      !! otherwise — the `z_draft` gating convention, latched in
      !! `ocean_state_init_from_config` BEFORE `init` so the allocation
      !! gate can read it.
      logical :: is_init = .false.
         !! Set by `init` after every allocation succeeds; cleared by
         !! `destroy` first.  Always test this, never `allocated(...)`.
      logical :: enable = .false.
         !! `&ocean_cavity_melt_nml enable`, latched before `init`.

      ! ---- Knobs the kernel is called with (host scalars) ----
      real(wp) :: far_field_depth = 10.0_wp
         !! Thickness (m) below the ice base the far field is averaged
         !! over.  See the module docstring: metres, not layers.
      real(wp) :: cdrag_top = 2.5e-3_wp
         !! Top drag coefficient for the MELT friction velocity only.
      real(wp) :: u_tide = 1.0e-2_wp
         !! RMS tidal velocity (m/s) in the melt `u*`.
      real(wp) :: ustar_min = 1.0e-4_wp
         !! Friction-velocity floor (m/s).
      real(wp) :: s_ice = 0.0_wp
         !! Ice salinity (g/kg).
      integer :: freshwater = CAVITY_FW_VIRTUAL
         !! `&ocean_cavity_melt_nml freshwater`, parsed.
         !! `CAVITY_FW_VIRTUAL` (default) ⇒ `ocean_cavity_mass_step` is
         !! an immediate return and the path is byte-identical.
      integer :: volume_comp = CAVITY_VC_NONE
         !! `&ocean_cavity_melt_nml volume_compensation`, parsed.
         !! Requires `freshwater = CAVITY_FW_MASS` (validate_config).
      real(wp) :: rho0 = 1035.0_wp
         !! Boussinesq reference density (kg/m^3) — THE density that
         !! converts the melt MASS flux into a thickness tendency, and
         !! the same one `ocean_surface_flux_t%rho0` divides the surface
         !! tracer sources by.  Assigned at configure from the single
         !! `rho_0` of record; the literal is only the pre-configure type
         !! default.  The two MUST agree or the salt undo is not the
         !! exact negation of what the apply kernel stamped, which is
         !! asserted at configure.
      type(ocean_cavity_exchange_t) :: par
         !! Exchange-law selector + coefficients.  Flat POD; its
         !! `f_cor` member is unused here (the driver substitutes the
         !! per-column `f_cor` array).
      type(ocean_cavity_ice_t) :: ice
         !! Ice-conduction selector + `T_ice`.  Flat POD.
      type(ocean_cavity_const_t) :: const
         !! Thermodynamic + turbulence constants.  Flat POD, left at the
         !! ISOMIP+ protocol values (Asay-Davis et al. (2016) Table 4):
         !! `rho_w`, `c_w`, `alpha_T`, `beta_S` here are the MELT LAW's
         !! own calibrated constants, NOT copies of the model's
         !! Boussinesq `rho_0` — the same "out of scope on purpose"
         !! category as `RHO_WATER` and `&ocean_ice_nml rho_ocean` in
         !! `src/core/ocean/README.md`'s reference-density table.
         !! Overriding them would break ISOMIP+ comparability, which is
         !! the reason this path exists.

      ! ---- Static, filled once at configure ----
      real(wp), allocatable :: f_cor(:, :)
         !! Cell-centred Coriolis parameter (1/s), from the same
         !! `fill_coriolis_centre` every other slot uses.  Read by
         !! `hj99` only, as `|f|`; the law does not exist at `f = 0`.

      ! ---- Per-step interface state (all device-resident) ----
      real(wp), allocatable :: active(:, :)
         !! Composed solve mask (0/1): `cover_frac` AND wet AND "the
         !! far-field sample found mass".  This is what the kernel's
         !! `cover` argument is given.
      real(wp), allocatable :: t_far(:, :)
         !! Far-field temperature (degC), thickness-weighted over
         !! `far_field_depth`.
      real(wp), allocatable :: s_far(:, :)
         !! Far-field salinity (g/kg), same average.
      real(wp), allocatable :: u_far(:, :)
         !! Far-field x velocity at the cell CENTRE (m/s), same average.
      real(wp), allocatable :: v_far(:, :)
         !! Far-field y velocity at the cell CENTRE (m/s), same average.
      real(wp), allocatable :: ustar(:, :)
         !! Melt friction velocity (m/s).
      real(wp), allocatable :: t_b(:, :)
         !! Interface temperature (degC), on the liquidus.
      real(wp), allocatable :: s_b(:, :)
         !! Interface salinity (g/kg).
      real(wp), allocatable :: melt(:, :)
         !! **Canonical** melt mass flux (kg/m^2/s), > 0 melting.
      real(wp), allocatable :: q_ocean(:, :)
         !! Turbulent heat flux ocean → interface (W/m^2), > 0 cools the
         !! ocean.
      real(wp), allocatable :: gamma_t(:, :)
         !! Thermal exchange velocity (m/s) the column's solve converged
         !! on — NOT `par%gamma_t_coeff`, which is the dimensionless
         !! coefficient.  Under `hj99` / `yung25` it is an implicit
         !! function of the converged interface state, so it is stored
         !! rather than re-derived: a diagnostic that rebuilt it from
         !! `u*` alone would report the NEUTRAL value instead of the
         !! stratification-suppressed one.  Zero where inactive.
      real(wp), allocatable :: gamma_s(:, :)
         !! Haline exchange velocity (m/s), same convention.
      integer, allocatable :: status(:, :)
         !! `CAVITY_MELT_*` per column; `CAVITY_MELT_OK` where inactive.
      real(wp), allocatable :: comp_scale(:, :)
         !! Per-cell top-layer thickness RATIO `h_new/h_old` left behind
         !! by the `volume_compensation` sink, exactly `1` everywhere the
         !! sink did not act.  Every PASSIVE tracer's top-layer load is
         !! multiplied by it (`cavity_comp_scale_tracer_impl`) so the
         !! removed parcel carries that cell's own concentration and
         !! changes none of them.  Allocated with the rest of the slot
         !! (full size iff `enable`); left at `1` when
         !! `volume_comp = CAVITY_VC_NONE`.

      ! ---- Status counters (host scalars, reduced off the device) ----
      integer :: n_nonfinite_step = 0
         !! Active columns returning `NONFINITE_INPUT`/`NONFINITE_STATE`
         !! this step.  **FATAL** — the driver fails loud.
      integer :: n_not_converged_step = 0
         !! Active columns returning `NOT_CONVERGED` this step (zero melt
         !! applied there — the kernel's documented safe state).
      integer :: n_no_root_step = 0
         !! Active columns returning `NO_PHYSICAL_ROOT` this step.
      integer :: n_other_step = 0
         !! Active columns returning any other non-OK status
         !! (`BAD_INPUT`, `NO_CORIOLIS`, `LAW_DOMAIN`, ...).
      integer(int64) :: n_not_converged_total = 0_int64
         !! Running total, drained to the console like
         !! `continuity_t%n_limited_total`.  int64 because a
         !! long run with one stubborn column can exceed 2e9.
      integer(int64) :: n_no_root_total = 0_int64
         !! Running total of `NO_PHYSICAL_ROOT`.
      integer(int64) :: n_other_total = 0_int64
         !! Running total of every other non-OK, non-fatal status.
      integer :: n_thin_step = 0
         !! Columns this thermo step whose top layer could NOT give up
         !! the thickness the real-mass path asked of it without falling
         !! through `H_CAVITY_FLOOR` — a FREEZING column (`m < 0`) whose
         !! top layer is thinner than `|m|*dt/rho_0`, or a compensation
         !! sink deeper than the open-ocean top layer.  The withdrawal is
         !! CLAMPED to what is there ABOVE the floor (so `h` stays
         !! strictly above the vanish marker, never on or below it) and
         !! the count is FATAL: a clamped withdrawal is a
         !! withdrawal the tracked mass source no longer matches, so
         !! continuing would print a budget that silently stopped
         !! closing.  Same stance as the non-finite melt column.
      integer(int64) :: n_thin_total = 0_int64
         !! Running total of `n_thin_step`, for the record.
      real(wp) :: melt_volume_step = 0.0_wp
         !! Domain-integrated meltwater VOLUME (m^3) added this thermo
         !! step over the INTERIOR (ghosts excluded, globally combined).
         !! Diagnostic + the input the compensation sink spends.
      real(wp) :: open_area = 0.0_wp
         !! Interior wet area (m^2) the ice does NOT cover, globally
         !! combined — the denominator of the uniform sink.
      real(wp) :: comp_withdrawal_step = 0.0_wp
         !! Thickness (m) the compensation sink removed from each open
         !! -ocean column this thermo step; zero when compensation is off.
   contains
      procedure, non_overridable :: init => ocean_cavity_flux_init
      procedure, non_overridable :: destroy => ocean_cavity_flux_destroy
      procedure, non_overridable :: enter_data => ocean_cavity_flux_enter_data
      procedure, non_overridable :: exit_data => ocean_cavity_flux_exit_data
      procedure, non_overridable :: bytes => ocean_cavity_flux_bytes
   end type ocean_cavity_flux_t

contains

   ! ======================================================================
   ! Lifecycle
   ! ======================================================================

   subroutine ocean_cavity_flux_init(this, grid)
      !! Allocate the slot.  Gated on `enable` (latched by
      !! `ocean_state_init_from_config` before this runs), so a run
      !! without a cavity pays fourteen `(1,1)` placeholders.
      class(ocean_cavity_flux_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      integer :: nx, ny

      if (this%enable) then
         nx = grid%nx_total
         ny = grid%ny_total
      else
         nx = 1
         ny = 1
      end if

      allocate (this%f_cor(nx, ny), source=0.0_wp)
      allocate (this%active(nx, ny), source=0.0_wp)
      allocate (this%t_far(nx, ny), source=0.0_wp)
      allocate (this%s_far(nx, ny), source=0.0_wp)
      allocate (this%u_far(nx, ny), source=0.0_wp)
      allocate (this%v_far(nx, ny), source=0.0_wp)
      allocate (this%ustar(nx, ny), source=0.0_wp)
      allocate (this%t_b(nx, ny), source=0.0_wp)
      allocate (this%s_b(nx, ny), source=0.0_wp)
      allocate (this%melt(nx, ny), source=0.0_wp)
      allocate (this%q_ocean(nx, ny), source=0.0_wp)
      allocate (this%gamma_t(nx, ny), source=0.0_wp)
      allocate (this%gamma_s(nx, ny), source=0.0_wp)
      allocate (this%status(nx, ny), source=CAVITY_MELT_OK)
      ! `comp_scale` is a RATIO, so its inert value is 1, not 0 — a
      ! zero-initialised plane would annihilate every passive tracer's
      ! top layer the first time it was applied.
      allocate (this%comp_scale(nx, ny), source=1.0_wp)
      this%is_init = .true.
   end subroutine ocean_cavity_flux_init

   subroutine ocean_cavity_flux_destroy(this)
      !! Release the slot.  `is_init` is cleared FIRST.
      class(ocean_cavity_flux_t), intent(inout) :: this
      this%is_init = .false.
      if (allocated(this%f_cor)) deallocate (this%f_cor)
      if (allocated(this%active)) deallocate (this%active)
      if (allocated(this%t_far)) deallocate (this%t_far)
      if (allocated(this%s_far)) deallocate (this%s_far)
      if (allocated(this%u_far)) deallocate (this%u_far)
      if (allocated(this%v_far)) deallocate (this%v_far)
      if (allocated(this%ustar)) deallocate (this%ustar)
      if (allocated(this%t_b)) deallocate (this%t_b)
      if (allocated(this%s_b)) deallocate (this%s_b)
      if (allocated(this%melt)) deallocate (this%melt)
      if (allocated(this%q_ocean)) deallocate (this%q_ocean)
      if (allocated(this%gamma_t)) deallocate (this%gamma_t)
      if (allocated(this%gamma_s)) deallocate (this%gamma_s)
      if (allocated(this%status)) deallocate (this%status)
      if (allocated(this%comp_scale)) deallocate (this%comp_scale)
   end subroutine ocean_cavity_flux_destroy

   subroutine ocean_cavity_flux_enter_data(this)
      !! Type-bound wrapper — delegates to the non-polymorphic impl so
      !! the device-attach map base is the heap object, not a
      !! polymorphic stack box.
      class(ocean_cavity_flux_t), intent(inout) :: this
      select type (this)
      type is (ocean_cavity_flux_t)
         call ocean_cavity_flux_enter_data_impl(this)
      end select
   end subroutine ocean_cavity_flux_enter_data

   subroutine ocean_cavity_flux_enter_data_impl(this)
      !! `copyin` (not `create`) throughout: `f_cor` carries a
      !! configure-time host fill that MUST reach the device, and the
      !! rest carry the zero `init` promised on both toolchains.
      type(ocean_cavity_flux_t), intent(inout) :: this
      !$acc enter data copyin(this%f_cor, this%active, this%t_far, this%s_far, &
      !$acc&                  this%u_far, this%v_far, this%ustar, this%t_b, &
      !$acc&                  this%s_b, this%melt, this%q_ocean, this%gamma_t, &
      !$acc&                  this%gamma_s, this%status, this%comp_scale)
      !$acc update device(this%f_cor, this%active, this%t_far, this%s_far, &
      !$acc&              this%u_far, this%v_far, this%ustar, this%t_b, &
      !$acc&              this%s_b, this%melt, this%q_ocean, this%gamma_t, &
      !$acc&              this%gamma_s, this%status, this%comp_scale)
   end subroutine ocean_cavity_flux_enter_data_impl

   subroutine ocean_cavity_flux_exit_data(this)
      class(ocean_cavity_flux_t), intent(inout) :: this
      select type (this)
      type is (ocean_cavity_flux_t)
         call ocean_cavity_flux_exit_data_impl(this)
      end select
   end subroutine ocean_cavity_flux_exit_data

   subroutine ocean_cavity_flux_exit_data_impl(this)
      type(ocean_cavity_flux_t), intent(inout) :: this
      !$acc exit data delete(this%f_cor, this%active, this%t_far, this%s_far, &
      !$acc&                 this%u_far, this%v_far, this%ustar, this%t_b, &
      !$acc&                 this%s_b, this%melt, this%q_ocean, this%gamma_t, &
      !$acc&                 this%gamma_s, this%status, this%comp_scale)
   end subroutine ocean_cavity_flux_exit_data_impl

   pure function ocean_cavity_flux_bytes(this) result(nbytes)
      !! Counted allocatable footprint (0 when unallocated).  One
      !! `arr_bytes` term per array — add one here when an array joins
      !! the type.
      class(ocean_cavity_flux_t), intent(in) :: this
      integer(int64) :: nbytes
      nbytes = arr_bytes(this%f_cor) &
               + arr_bytes(this%active) &
               + arr_bytes(this%t_far) &
               + arr_bytes(this%s_far) &
               + arr_bytes(this%u_far) &
               + arr_bytes(this%v_far) &
               + arr_bytes(this%ustar) &
               + arr_bytes(this%t_b) &
               + arr_bytes(this%s_b) &
               + arr_bytes(this%melt) &
               + arr_bytes(this%q_ocean) &
               + arr_bytes(this%gamma_t) &
               + arr_bytes(this%gamma_s) &
               + arr_bytes(this%status) &
               + arr_bytes(this%comp_scale)
   end function ocean_cavity_flux_bytes

   ! ======================================================================
   ! Far-field sampling
   ! ======================================================================

   pure subroutine cavity_far_field_impl(nx, ny, nz, far_depth, cover, wet_mask, &
                                         h_layer, hTr_T, hTr_S, u_face_x, v_face_y, &
                                         active, t_far, s_far, u_far, v_far)
      !! Thickness-weighted mean of `(T, S, u, v)` over `far_depth`
      !! METRES below the ice base, with a PARTIAL last layer.
      !!
      !! Bottom-up stack (`k = nz` is the surface layer, i.e. the one
      !! against the ice base), so the walk runs `k = nz, 1, -1` and
      !! stops when the budget is spent.  `far_depth >= column
      !! thickness` ⇒ the whole column, which is the correct limit and
      !! not an error.
      !!
      !! Vanished layers (`h <= H_VANISHED`) are SKIPPED, not clamped —
      !! the D4 taxonomy's skip/merge marker.  They carry no mass, so
      !! including them would be dividing a zero tracer load by a
      !! near-zero thickness.
      !!
      !! Velocities are centred from the C-grid faces per layer BEFORE
      !! the vertical average (`u_c = 1/2 (u_{i} + u_{i+1})`), not after:
      !! the two commute only for a uniform column, and the order that
      !! matches "the speed the ice base feels" is the one that averages
      !! the centred, per-layer velocity.
      !!
      !! `active` is the composed solve mask handed to the melt kernel:
      !! covered AND wet AND the sample found mass.  A covered column
      !! whose entire sample is vanished gets `active = 0` and zeroed
      !! outputs — no melt, no status noise.
      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) :: far_depth
         !! Sampling thickness (m), > 0.
      real(wp), intent(in) :: cover(nx, ny)
         !! Ice-cover fraction (v1 binary).
      real(wp), intent(in) :: wet_mask(nx, ny)
         !! Static wet (1) / land (0) mask.
      real(wp), intent(in) :: h_layer(nx, ny, nz)
         !! Layer thickness (m).
      real(wp), intent(in) :: hTr_T(nx, ny, nz)
         !! `h*T` (degC m).
      real(wp), intent(in) :: hTr_S(nx, ny, nz)
         !! `h*S` ((g/kg) m).
      real(wp), intent(in) :: u_face_x(nx + 1, ny, nz)
         !! Zonal face velocity (m/s).
      real(wp), intent(in) :: v_face_y(nx, ny + 1, nz)
         !! Meridional face velocity (m/s).
      real(wp), intent(out) :: active(nx, ny)
         !! Composed solve mask, 0 or 1.
      real(wp), intent(out) :: t_far(nx, ny)
         !! Sampled temperature (degC); 0 where inactive.
      real(wp), intent(out) :: s_far(nx, ny)
         !! Sampled salinity (g/kg); 0 where inactive.
      real(wp), intent(out) :: u_far(nx, ny)
         !! Sampled centred x velocity (m/s); 0 where inactive.
      real(wp), intent(out) :: v_far(nx, ny)
         !! Sampled centred y velocity (m/s); 0 where inactive.
      integer :: i, j, k
      real(wp) :: acc_h, acc_t, acc_s, acc_u, acc_v, remain, h, w, inv_h

      do concurrent(j=1:ny, i=1:nx) &
         local(k, acc_h, acc_t, acc_s, acc_u, acc_v, remain, h, w, inv_h)
         acc_h = 0.0_wp
         acc_t = 0.0_wp
         acc_s = 0.0_wp
         acc_u = 0.0_wp
         acc_v = 0.0_wp
         remain = far_depth
         if (cover(i, j) > 0.5_wp .and. wet_mask(i, j) > 0.5_wp) then
            do k = nz, 1, -1
               if (remain <= 0.0_wp) exit
               h = h_layer(i, j, k)
               ! vanished-ok: the far-field sampler SKIPS a vanished layer entirely (it carries
               ! no water to melt against) rather than weighting a substituted
               ! concentration into the average.
               if (h <= H_VANISHED) cycle
               w = min(h, remain)
               inv_h = 1.0_wp/h
               acc_h = acc_h + w
               acc_t = acc_t + w*hTr_T(i, j, k)*inv_h
               acc_s = acc_s + w*hTr_S(i, j, k)*inv_h
               acc_u = acc_u + w*0.5_wp*(u_face_x(i, j, k) + u_face_x(i + 1, j, k))
               acc_v = acc_v + w*0.5_wp*(v_face_y(i, j, k) + v_face_y(i, j + 1, k))
               remain = remain - w
            end do
         end if
         if (acc_h > 0.0_wp) then
            active(i, j) = 1.0_wp
            t_far(i, j) = acc_t/acc_h
            s_far(i, j) = acc_s/acc_h
            u_far(i, j) = acc_u/acc_h
            v_far(i, j) = acc_v/acc_h
         else
            active(i, j) = 0.0_wp
            t_far(i, j) = 0.0_wp
            s_far(i, j) = 0.0_wp
            u_far(i, j) = 0.0_wp
            v_far(i, j) = 0.0_wp
         end if
      end do
   end subroutine cavity_far_field_impl

   ! ======================================================================
   ! Flux delivery
   ! ======================================================================

   pure subroutine cavity_flux_fill_impl(nx, ny, s_ice, active, melt, q_ocean, s_far, &
                                         heat_cavity, salt_cavity)
      !! Fill the two OWNED surface-flux components from the solved
      !! interface.  FULL OVERWRITE of the whole plane, never `+=` —
      !! the same rule the sea-ice coupler's fillers follow, and for the
      !! same reason (a `+=` ratchets across outer steps with no bound).
      !!
      !! See the module docstring for the sign derivation.  In one line:
      !! the ocean loses the heat the interface takes (`-q_ocean`), and
      !! the dilution by meltwater of salinity `s_ice` is emulated at
      !! fixed column mass by `-m_mass*(S_far - s_ice)`.
      !!
      !! No wet-mask factor here: `active` already carries it (the
      !! sampler composes it), and the assembler multiplies `Q_heat` /
      !! `Q_salt` by `ms%wet_mask` again.
      integer, intent(in) :: nx
         !! First dimension (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      real(wp), intent(in) :: s_ice
         !! Ice salinity (g/kg).
      real(wp), intent(in) :: active(nx, ny)
         !! Composed solve mask.
      real(wp), intent(in) :: melt(nx, ny)
         !! Melt mass flux (kg/m^2/s), > 0 melting.
      real(wp), intent(in) :: q_ocean(nx, ny)
         !! Ocean → interface heat flux (W/m^2).
      real(wp), intent(in) :: s_far(nx, ny)
         !! Far-field salinity the solve used (g/kg).
      real(wp), intent(out) :: heat_cavity(nx, ny)
         !! Heat component (W/m^2), positive down into the ocean.
      real(wp), intent(out) :: salt_cavity(nx, ny)
         !! Virtual salt component, positive salinifies.
      integer :: i, j

      do concurrent(j=1:ny, i=1:nx)
         if (active(i, j) > 0.5_wp) then
            heat_cavity(i, j) = -q_ocean(i, j)
            salt_cavity(i, j) = -melt(i, j)*(s_far(i, j) - s_ice)
         else
            heat_cavity(i, j) = 0.0_wp
            salt_cavity(i, j) = 0.0_wp
         end if
      end do
   end subroutine cavity_flux_fill_impl

   ! ======================================================================
   ! Status accounting
   ! ======================================================================

   pure subroutine cavity_status_counts_impl(nx, ny, active, status, n_nonfinite, &
                                             n_not_converged, n_no_root, n_other)
      !! Reduce the per-column status plane into four counts, ON DEVICE.
      !!
      !! `do concurrent ... reduce(+:)` and not `count(...)`: the array
      !! is device-resident under `mem:separate`, and a host intrinsic
      !! would silently read the stale host shadow — the same trap
      !! `continuity_t%n_limited_step` documents.
      !!
      !! Only ACTIVE columns are counted.  An inactive column always
      !! carries `CAVITY_MELT_OK`, but gating the count on the mask
      !! keeps that a property of the caller rather than of the kernel.
      integer, intent(in) :: nx
         !! First dimension.
      integer, intent(in) :: ny
         !! Second dimension.
      real(wp), intent(in) :: active(nx, ny)
         !! Composed solve mask.
      integer, intent(in) :: status(nx, ny)
         !! `CAVITY_MELT_*` per column.
      integer, intent(out) :: n_nonfinite
         !! Count of `NONFINITE_INPUT` + `NONFINITE_STATE`.
      integer, intent(out) :: n_not_converged
         !! Count of `NOT_CONVERGED`.
      integer, intent(out) :: n_no_root
         !! Count of `NO_PHYSICAL_ROOT`.
      integer, intent(out) :: n_other
         !! Count of every other non-OK status.
      integer :: i, j, c_nf, c_nc, c_nr, c_ot

      c_nf = 0
      c_nc = 0
      c_nr = 0
      c_ot = 0
      do concurrent(j=1:ny, i=1:nx) reduce(+:c_nf, c_nc, c_nr, c_ot)
         if (active(i, j) > 0.5_wp) then
            select case (status(i, j))
            case (CAVITY_MELT_OK)
               continue
            case (CAVITY_MELT_NONFINITE_INPUT, CAVITY_MELT_NONFINITE_STATE)
               c_nf = c_nf + 1
            case (CAVITY_MELT_NOT_CONVERGED)
               c_nc = c_nc + 1
            case (CAVITY_MELT_NO_PHYSICAL_ROOT)
               c_nr = c_nr + 1
            case default
               c_ot = c_ot + 1
            end select
         end if
      end do
      n_nonfinite = c_nf
      n_not_converged = c_nc
      n_no_root = c_nr
      n_other = c_ot
   end subroutine cavity_status_counts_impl

   pure function cavity_melt_status_is_fatal(n_nonfinite) result(is_fatal)
      !! Is this step's status tally a FAIL-LOUD condition?
      !!
      !! A non-finite input or intermediate on a covered column means the
      !! state feeding the interface is already corrupt; the kernel's
      !! safe state (zero melt) would hide it behind a plausible run.
      !! `NOT_CONVERGED` and `NO_PHYSICAL_ROOT` are NOT fatal — they are
      !! counted and warned, and those columns take the documented
      !! zero-melt safe state.
      !!
      !! A `pure` predicate so the test suite can assert the DECISION
      !! without provoking `error stop`, the way the repo tests every
      !! other fail-loud rule.
      integer, intent(in) :: n_nonfinite
         !! `n_nonfinite_step` from `cavity_status_counts_impl`.
      logical :: is_fatal
      is_fatal = (n_nonfinite > 0)
   end function cavity_melt_status_is_fatal

   ! ======================================================================
   ! Phase 3 — real freshwater MASS
   ! ======================================================================

   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

   pure subroutine cavity_mass_salt_mirror_impl(nx, ny, nz, dt_over_rho0, inv_rho0_dt, &
                                                s_ice, active, melt, s_far, hTr, k_top)
      !! The pseudo-salt mirror of `cavity_mass_apply_impl`'s SALT
      !! increment, with no budget accumulation — `budget_id = NONE`, so
      !! it must not touch `salt_budget_surface`.
      !!
      !! Pseudo-salt's documented purpose is to be a PASSIVE tracer given
      !! exactly salinity's surface salt flux, so the deviation measures
      !! the passive-vs-active transport-path error.  Under the real-mass
      !! form salinity's surface salt flux is `w*s_ice` (and the dilution
      !! by the added volume, which pseudo-salt gets for free because it
      !! shares `h`), so the mirror must receive the same correction — or
      !! it would keep the virtual flux salinity no longer has and the
      !! deviation would measure the bookkeeping instead of the transport.
      !!
      !! The arithmetic is the SAME EXPRESSION as the salt branch above,
      !! so the increment is bit-identical to salinity's.
      integer, intent(in) :: nx
         !! First dimension.
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nz
         !! Layer count.
      real(wp), intent(in) :: dt_over_rho0
         !! `dt/rho_0`.
      real(wp), intent(in) :: inv_rho0_dt
         !! The same `dt/rho_0` (see the sibling's docstring).
      real(wp), intent(in) :: s_ice
         !! Ice salinity (g/kg).
      real(wp), intent(in) :: active(nx, ny)
         !! Composed solve mask.
      real(wp), intent(in) :: melt(nx, ny)
         !! Melt mass flux (kg/m^2/s).
      real(wp), intent(in) :: s_far(nx, ny)
         !! Far-field salinity (g/kg); the virtual component being
         !! undone is rebuilt from it, exactly as in the sibling.
      real(wp), intent(inout) :: hTr(nx, ny, nz)
         !! Pseudo-salt `h*C`.
      integer, intent(in) :: k_top(nx, ny)
         !! The SAME first-live-layer index the salt branch used — the
         !! mirror's whole contract is that its increment is bit-identical
         !! to salinity's, which includes landing on the same row.
      integer :: i, j
      real(wp) :: cell_s, virt

      do concurrent(j=1:ny, i=1:nx) local(cell_s, virt)
         if (active(i, j) > 0.5_wp) then
            virt = -melt(i, j)*(s_far(i, j) - s_ice)
            cell_s = -inv_rho0_dt*virt + melt(i, j)*dt_over_rho0*s_ice
            hTr(i, j, k_top(i, j)) = hTr(i, j, k_top(i, j)) + cell_s
         end if
      end do
   end subroutine cavity_mass_salt_mirror_impl

   pure subroutine cavity_mass_totals_impl(nx, ny, nghost, dt_over_rho0, active, &
                                           wet_mask, cover_frac, melt, areaT, &
                                           vol_melt, area_open)
      !! The two INTERIOR integrals the mass budget and the compensation
      !! sink need, reduced on device in one pass:
      !!
      !!   `vol_melt`  — the meltwater volume (m^3) added this step,
      !!                 `sum(melt*dt/rho_0 * areaT)` over active columns;
      !!   `area_open` — the wet area (m^2) the ice does NOT cover.
      !!
      !! GHOSTS ARE EXCLUDED (`nghost+1 .. n-nghost`), exactly as
      !! `ocean_accumulate_mass_out` excludes them: the console's totals
      !! are interior integrals, so a source term that counted the halo
      !! would be measuring a different domain than the total it has to
      !! close.  The per-cell APPLY still runs over the whole array, like
      !! every other surface kernel — the halo exchange owns the ghosts.
      !!
      !! `melt*dt_over_rho0` is formed with the same scalar and in the
      !! same order as the apply kernel's `dh`, so the two agree to the
      !! last bit per cell and the residual is reduction round-off only.
      integer, intent(in) :: nx
         !! First dimension (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nghost
         !! Halo width to exclude on every side.
      real(wp), intent(in) :: dt_over_rho0
         !! `dt/rho_0`.
      real(wp), intent(in) :: active(nx, ny)
         !! Composed solve mask.
      real(wp), intent(in) :: wet_mask(nx, ny)
         !! Static wet (1) / land (0) mask.
      real(wp), intent(in) :: cover_frac(nx, ny)
         !! Ice-cover fraction (v1 binary).
      real(wp), intent(in) :: melt(nx, ny)
         !! Melt mass flux (kg/m^2/s).
      real(wp), intent(in) :: areaT(nx, ny)
         !! Cell area (m^2).
      real(wp), intent(out) :: vol_melt
         !! Meltwater volume this step (m^3).
      real(wp), intent(out) :: area_open
         !! Uncovered wet area (m^2).
      integer :: i, j, i_lo, i_hi, j_lo, j_hi
      real(wp) :: acc_v, acc_a

      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost
      acc_v = 0.0_wp
      acc_a = 0.0_wp
      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) reduce(+:acc_v, acc_a)
         if (active(i, j) > 0.5_wp) then
            acc_v = acc_v + melt(i, j)*dt_over_rho0*areaT(i, j)
         end if
         if (wet_mask(i, j) > 0.5_wp .and. cover_frac(i, j) < 0.5_wp) then
            acc_a = acc_a + areaT(i, j)
         end if
      end do
      vol_melt = acc_v
      area_open = acc_a
   end subroutine cavity_mass_totals_impl

   pure function cavity_comp_withdrawal(vol_melt, area_open) result(dw)
      !! The uniform per-unit-area thickness the `uniform_open_ocean`
      !! sink removes: `vol_melt/area_open`, and exactly zero when there
      !! is no open ocean to remove it from (a fully ice-covered domain —
      !! the sink then does nothing rather than dividing by zero, and the
      !! volume stays in, which the mass budget reports honestly as a
      !! growing total).
      !!
      !! `pure`, separate from the kernels, so the test suite can assert
      !! the arithmetic and the no-open-ocean fallback directly.
      real(wp), intent(in) :: vol_melt
         !! Meltwater volume this step (m^3); either sign.
      real(wp), intent(in) :: area_open
         !! Uncovered wet area (m^2), >= 0.
      real(wp) :: dw
      if (area_open > H_DIV_EPS) then
         dw = vol_melt/area_open
      else
         dw = 0.0_wp
      end if
   end function cavity_comp_withdrawal

   pure subroutine cavity_comp_apply_impl(nx, ny, nz, dw, wet_mask, cover_frac, &
                                          h_layer, hTr_S, hTr_T, &
                                          salt_budget, heat_budget, scale, n_thin)
      !! The `volume_compensation = "uniform_open_ocean"` sink: remove
      !! `dw` metres of the top layer from every wet cell the ice does
      !! NOT cover, the removed parcel carrying that cell's own `T` and
      !! `S` so no concentration there changes.
      !!
      !! "Carries its own concentration" is implemented as a RATIO:
      !! `hTr *= h_new/h_old`, which leaves `Tr = hTr/h` algebraically
      !! identical.  The ratio is published in `scale` so every PASSIVE
      !! tracer can be given exactly the same treatment
      !! (`cavity_comp_scale_tracer_impl`) — a tracer left alone here
      !! would be CONCENTRATED by the sink, which is a different physical
      !! statement from the one the knob makes.
      !!
      !! The salt and heat increments ARE mirrored into the surface
      !! budget contributors: a parcel that carries `S` out of the domain
      !! changes the domain salt total, so the budget has to name it.
      !!
      !! Same clamp (`H_CAVITY_FLOOR`, strictly above the vanish marker)
      !! + count + fatal policy as the mass source.
      integer, intent(in) :: nx
         !! First dimension (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nz
         !! Layer count; only `k = nz` is touched.
      real(wp), intent(in) :: dw
         !! Thickness (m) to remove from each open-ocean column.
      real(wp), intent(in) :: wet_mask(nx, ny)
         !! Static wet (1) / land (0) mask.
      real(wp), intent(in) :: cover_frac(nx, ny)
         !! Ice-cover fraction (v1 binary).
      real(wp), intent(inout) :: h_layer(nx, ny, nz)
         !! Layer thickness (m).
      real(wp), intent(inout) :: hTr_S(nx, ny, nz)
         !! `h*S`.
      real(wp), intent(inout) :: hTr_T(nx, ny, nz)
         !! `h*T`.
      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`.
      real(wp), intent(out) :: scale(nx, ny)
         !! `h_new/h_old` where the sink acted, exactly 1 elsewhere.
      integer, intent(out) :: n_thin
         !! Columns whose withdrawal had to be clamped.
      integer :: i, j, c_thin
      real(wp) :: h0, hn, d, f, cell_s, cell_t

      c_thin = 0
      do concurrent(j=1:ny, i=1:nx) local(h0, hn, d, f, cell_s, cell_t) reduce(+:c_thin)
         if (wet_mask(i, j) > 0.5_wp .and. cover_frac(i, j) < 0.5_wp) then
            h0 = h_layer(i, j, nz)
            d = dw
            hn = h0 - d
            if (hn < H_CAVITY_FLOOR) then
               ! Same floor and same refusal rule as the mass source: pin
               ! the RESULT (never the increment — `h0 - (h0 - floor)`
               ! rounds and can land a ulp below), take only what sits
               ! ABOVE the floor, and take nothing at all from a column
               ! already at or below it (`min(h0, ...)`, so `hn <= h0`
               ! always and the sink can never become a deposit).
               hn = min(h0, H_CAVITY_FLOOR)
               c_thin = c_thin + 1
            end if
            if (h0 > H_DIV_EPS) then
               f = hn/h0
            else
               f = 1.0_wp
            end if
            cell_s = hTr_S(i, j, nz)*(f - 1.0_wp)
            cell_t = hTr_T(i, j, nz)*(f - 1.0_wp)
            h_layer(i, j, nz) = hn
            hTr_S(i, j, nz) = hTr_S(i, j, nz) + cell_s
            salt_budget(i, j, nz) = salt_budget(i, j, nz) + cell_s
            hTr_T(i, j, nz) = hTr_T(i, j, nz) + cell_t
            heat_budget(i, j, nz) = heat_budget(i, j, nz) + cell_t
            scale(i, j) = f
         else
            scale(i, j) = 1.0_wp
         end if
      end do
      n_thin = c_thin
   end subroutine cavity_comp_apply_impl

   pure subroutine cavity_comp_scale_tracer_impl(nx, ny, nz, scale, hTr)
      !! Apply the compensation sink's top-layer ratio to one PASSIVE
      !! tracer's load.  Unconditional (`scale = 1` off the sink), so
      !! there is no mask branch inside the kernel.
      integer, intent(in) :: nx
         !! First dimension.
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nz
         !! Layer count.
      real(wp), intent(in) :: scale(nx, ny)
         !! `h_new/h_old` from `cavity_comp_apply_impl`.
      real(wp), intent(inout) :: hTr(nx, ny, nz)
         !! The tracer's `h*C`.
      integer :: i, j
      do concurrent(j=1:ny, i=1:nx)
         hTr(i, j, nz) = hTr(i, j, nz)*scale(i, j)
      end do
   end subroutine cavity_comp_scale_tracer_impl

   pure function cavity_mass_thin_is_fatal(n_thin) result(is_fatal)
      !! Is this step's clamped-withdrawal tally a FAIL-LOUD condition?
      !!
      !! Yes, on any count.  A clamped withdrawal is one the tracked mass
      !! source no longer matches, so from that step on the console's
      !! `Mass Error` would stop being a leak measurement while still
      !! being printed as one — the exact failure mode the budget gate
      !! exists to make impossible.  `h` is nevertheless clamped first,
      !! so the state handed to the abort path is finite and the last
      !! printed line is meaningful.
      !!
      !! A `pure` predicate, like `cavity_melt_status_is_fatal`, so the
      !! suite can assert the DECISION without provoking `error stop`.
      integer, intent(in) :: n_thin
         !! `n_thin_step`.
      logical :: is_fatal
      is_fatal = (n_thin > 0)
   end function cavity_mass_thin_is_fatal

   subroutine ocean_cavity_mass_step(grid, metrics, cav, ms, dt, weight, active)
      !! The real-freshwater MASS update — one call, at the THERMO
      !! cadence, from inside the RK2 stage immediately after
      !! `ocean_surface_flux_apply_tracers`.
      !!
      !! WHY THERE, and not in `engine_step_finalize` next to the melt
      !! solve.  Three things have to line up:
      !!
      !!   1. **The same melt value.**  The salt and heat halves of the
      !!      melt ride `Q_salt`/`Q_heat`, which the stage's surface-flux
      !!      apply spends.  Putting the volume anywhere else would spend
      !!      `melt(n)` for mass and `melt(n-1)` for salt.
      !!   2. **The same stage weight, by construction.**  Under
      !!      `ssp_rk2` both stages apply and `rk2_average` halves the
      !!      pair; under `pred_corr` `therm_active` is false in the
      !!      predictor and the corrector's single application IS the
      !!      step.  The caller passes the matching `weight` for the
      !!      budget accumulator, the same one `ocean_accumulate_mass_out`
      !!      takes.
      !!   3. **The barotropic mode sees it.**  `derive_bt_from_layers`
      !!      rebuilds `bt_eta = sum_k h - bt_H_ref` at the TOP of every
      !!      stage, so a thickness source applied here is in `bt_eta`
      !!      one stage later with no separate barotropic forcing term —
      !!      the free surface under the cavity datum simply rises.  The
      !!      BT substep's transport renormalisation (`bt_uhbt`) is a
      !!      constraint on the layer TRANSPORTS within a stage and never
      !!      reads a thickness source, so it is undisturbed.  This is
      !!      the documented `F_slow`-style operator split, and it is what
      !!      MOM6 does with its own surface mass fluxes.
      !!
      !! It runs BEFORE the vertical-mixing block in the same stage, so
      !! the implicit vdiff/drag solves see the thickened top layer —
      !! which is the right order: the added volume is part of the column
      !! before the column is mixed.  The ALE remap then redistributes it
      !! to the coordinate's target after `rk2_average`.
      !!
      !! NON-`pure`: it logs, it fails loud, it mutates `ms%mass_src`,
      !! and (with compensation on) it issues one collective.  Its six
      !! kernels are all `pure`.
      !!
      !! `mem:separate`: every array is mapped by its owning slot's
      !! `enter_data`; the tracer-registry derefs happen on the HOST
      !! before each kernel (the outer-shim rule).
      use pic_logger, only: global_logger
      use pic_strings, only: to_string
      use rdb_error_ring, only: fail
      use rdb_ocean_status, only: OCEAN_STATUS_ERR_SETUP
      use rdb_halo, only: halo_allreduce_sum
      type(hgrid_t), intent(in) :: grid
         !! Horizontal grid (`nx_total`/`ny_total`/`nghost`).
      type(ocean_metrics_t), intent(in) :: metrics
         !! Reads `cover_frac` and `areaT`.
      type(ocean_cavity_flux_t), intent(inout) :: cav
         !! The cavity-melt slot; reads `melt`/`active`/`t_b`, writes
         !! `comp_scale` and the thin counters.
      type(multilayer_state_t), intent(inout) :: ms
         !! Writes `h_layer(:,:,nz)`, the S/T/passive tracer loads, the
         !! two surface budget contributors and `mass_src`.
      real(wp), intent(in) :: dt
         !! Thermo timestep (s) — the same `therm_dt` the surface-flux
         !! apply was given.
      real(wp), intent(in) :: weight
         !! Per-stage weight for the `mass_src` accumulator (0.5 per
         !! SSP-RK2 stage; 0 / 1 for the pred_corr predictor /
         !! corrector), matching `ocean_accumulate_mass_out`.
      logical, intent(in), optional :: active
         !! Thermo-cadence gate.  Present-and-false ⇒ early return;
         !! absent ⇒ run.
      integer :: nx, ny, nz, idx_t, idx_s, idx_ps, it, n_thin_src, n_thin_comp
      real(wp) :: dt_over_rho0, vol_melt, area_open, dw, tmp

      if (.not. cav%is_init) return
      if (.not. cav%enable) return
      if (cav%freshwater /= CAVITY_FW_MASS) return
      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. allocated(ms%tracers)) return
      idx_t = ms%idx_temperature
      idx_s = ms%idx_salinity
      if (idx_t <= 0 .or. idx_s <= 0) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      idx_ps = ms%idx_pseudo_salt
      dt_over_rho0 = dt/cav%rho0

      ! (1) Interior integrals FIRST: the melt volume the budget tracks
      ! and the open-ocean area the sink would spend it over.  Combined
      ! across ranks so a decomposed domain removes the same total the
      ! whole domain gained.
      call cavity_mass_totals_impl(nx, ny, grid%nghost, dt_over_rho0, cav%active, &
                                   ms%wet_mask, metrics%cover_frac, cav%melt, &
                                   metrics%areaT, vol_melt, area_open)
      tmp = vol_melt
      call halo_allreduce_sum(tmp, vol_melt)
      tmp = area_open
      call halo_allreduce_sum(tmp, area_open)
      cav%melt_volume_step = vol_melt
      cav%open_area = area_open

      ! (2) The source.
      call cavity_mass_apply_impl(nx, ny, nz, dt_over_rho0, dt_over_rho0, cav%s_ice, &
                                  cav%active, cav%melt, cav%s_far, cav%t_b, &
                                  ms%h_layer, ms%tracers(idx_s)%hTr, ms%tracers(idx_t)%hTr, &
                                  ms%salt_budget_surface, ms%heat_budget_surface, &
                                  ms%k_top, n_thin_src)
      if (idx_ps > 0) then
         call cavity_mass_salt_mirror_impl(nx, ny, nz, dt_over_rho0, dt_over_rho0, &
                                           cav%s_ice, cav%active, cav%melt, &
                                           cav%s_far, ms%tracers(idx_ps)%hTr, ms%k_top)
      end if
      ms%mass_src = ms%mass_src + weight*RHO_WATER*vol_melt

      ! (3) The sink, when the knob asks for it.
      n_thin_comp = 0
      cav%comp_withdrawal_step = 0.0_wp
      if (cav%volume_comp == CAVITY_VC_UNIFORM_OPEN) then
         dw = cavity_comp_withdrawal(vol_melt, area_open)
         cav%comp_withdrawal_step = dw
         call cavity_comp_apply_impl(nx, ny, nz, dw, ms%wet_mask, metrics%cover_frac, &
                                     ms%h_layer, ms%tracers(idx_s)%hTr, &
                                     ms%tracers(idx_t)%hTr, ms%salt_budget_surface, &
                                     ms%heat_budget_surface, cav%comp_scale, n_thin_comp)
         ! NOTE the sink stays on `k = nz` and is NOT routed through
         ! `k_top`, deliberately: its own gate is `cover_frac < 0.5`, so
         ! it only ever acts on OPEN-OCEAN columns, and an open-ocean
         ! column has `z_top = 0` ⇒ no top-side filler ⇒ `k_top ≡ nz`.
         ! Routing it would be a provable no-op; leaving it spells out
         ! that "the top layer" and "the first live layer" are the same
         ! row wherever this kernel runs.  `cavity_comp_scale_tracer_impl`
         ! below rides the same argument (it is unconditional, but
         ! `comp_scale` is exactly 1 off the sink).
         do it = 1, size(ms%tracers)
            if (it == idx_s .or. it == idx_t) cycle
            call cavity_comp_scale_tracer_impl(nx, ny, nz, cav%comp_scale, &
                                               ms%tracers(it)%hTr)
         end do
         ! The removed interior volume is exactly `dw*area_open` — the
         ! same product the withdrawal was derived from — so with no
         ! clamping this cancels the source term to the last bit and the
         ! domain mass is constant.
         ms%mass_src = ms%mass_src - weight*RHO_WATER*dw*area_open
      end if

      ! (4) Accounting, then fail loud.
      cav%n_thin_step = n_thin_src + n_thin_comp
      cav%n_thin_total = cav%n_thin_total + int(cav%n_thin_step, int64)
      if (cavity_mass_thin_is_fatal(cav%n_thin_step)) then
         call global_logger%error("cavity real freshwater: "// &
                                  to_string(n_thin_src)//" melt column(s) and "// &
                                  to_string(n_thin_comp)//" compensation column(s) "// &
                                  "could not give up the requested thickness")
         call fail("&ocean_cavity_melt_nml freshwater='mass': "// &
                   to_string(cav%n_thin_step)//" column(s) would have been driven "// &
                   "below H_CAVITY_FLOOR (2*H_VANISHED) by the top-layer "// &
                   "withdrawal (a freezing "// &
                   "column thinner than |m|*dt/rho_0, or a compensation sink "// &
                   "deeper than the open-ocean top layer).  The withdrawal was "// &
                   "clamped so the state stays finite, but a clamped withdrawal "// &
                   "no longer matches the tracked mass source, so the console's "// &
                   "Mass Error would stop being a leak measurement while still "// &
                   "being printed as one.  Use a thicker top layer (a coordinate "// &
                   "with a larger surface target), a shorter dt, or "// &
                   "volume_compensation='none'.", code=OCEAN_STATUS_ERR_SETUP)
      end if
   end subroutine ocean_cavity_mass_step

   ! ======================================================================
   ! The driver
   ! ======================================================================

   subroutine ocean_cavity_flux_step(grid, cav, metrics, ms, eos, sf, active)
      !! One cavity basal-melt update: sample the far field, solve the
      !! three-equation interface on every covered column, deliver the
      !! two owned surface-flux components, and account for the solver
      !! status.
      !!
      !! CADENCE + PLACEMENT.  Called once per outer step at THERMO
      !! cadence from `engine_step_finalize`, immediately BEFORE
      !! `ocean_surface_flux_assemble` — the assembler must see the
      !! components this routine writes.  The tracers then integrate the
      !! assembled `Q_heat`/`Q_salt` on the NEXT outer step, which is the
      !! same one-step lag the sea-ice coupler documents.
      !!
      !! NON-`pure`, and it is the only non-pure procedure this module
      !! adds: it logs (`n_not_converged` warnings) and it FAILS LOUD on
      !! a non-finite column.  Its three kernels are all `pure`.
      !!
      !! `mem:separate`: every array read or written here is mapped by
      !! its owning slot's `enter_data`.  `ms%tracers(idx)%hTr` is
      !! dereferenced on the HOST before the kernels (the outer-shim
      !! rule for the array-of-derived-types registry).
      use pic_logger, only: global_logger
      use pic_strings, only: to_string
      use rdb_error_ring, only: fail
      use rdb_ocean_status, only: OCEAN_STATUS_ERR_SETUP
      type(hgrid_t), intent(in) :: grid
         !! Horizontal grid (for `nx_total`/`ny_total`).
      type(ocean_cavity_flux_t), intent(inout) :: cav
         !! The cavity-melt slot.
      type(ocean_metrics_t), intent(in) :: metrics
         !! Reads `cover_frac` only.
      type(multilayer_state_t), intent(in) :: ms
         !! Reads `h_layer`, the S/T tracer loads, the face velocities,
         !! `p_top` (THE interface pressure) and `wet_mask`.
      type(eos_t), intent(in) :: eos
         !! Shared EOS handle — the liquidus.  Flat POD, by value.
      type(ocean_surface_flux_t), intent(inout) :: sf
         !! Writes `heat_cavity`/`salt_cavity` and latches
         !! `has_heat`/`has_salt` host-side.
      logical, intent(in), optional :: active
         !! Thermo-cadence gate.  Present-and-false ⇒ early return;
         !! absent ⇒ run (the `ocean_surface_flux_assemble` convention).
      integer :: nx, ny, nz, idx_t, idx_s

      if (.not. cav%is_init) return
      if (.not. cav%enable) return
      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. allocated(ms%tracers)) return
      idx_t = ms%idx_temperature
      idx_s = ms%idx_salinity
      if (idx_t <= 0 .or. idx_s <= 0) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      call cavity_far_field_impl(nx, ny, nz, cav%far_field_depth, &
                                 metrics%cover_frac, ms%wet_mask, ms%h_layer, &
                                 ms%tracers(idx_t)%hTr, ms%tracers(idx_s)%hTr, &
                                 ms%u_face_x_layer, ms%v_face_y_layer, &
                                 cav%active, cav%t_far, cav%s_far, cav%u_far, cav%v_far)

      call cavity_melt_columns_2d(nx, ny, cav%active, cav%t_far, cav%s_far, ms%p_top, &
                                  cav%u_far, cav%v_far, cav%s_ice, cav%f_cor, &
                                  cav%cdrag_top, cav%u_tide, cav%ustar_min, &
                                  cav%par, cav%ice, eos, cav%const, &
                                  cav%ustar, cav%t_b, cav%s_b, cav%melt, cav%q_ocean, &
                                  cav%gamma_t, cav%gamma_s, cav%status)

      call cavity_flux_fill_impl(nx, ny, cav%s_ice, cav%active, cav%melt, cav%q_ocean, &
                                 cav%s_far, sf%heat_cavity, sf%salt_cavity)

      ! The filler's own obligations (PR-12 fill contract): latch the
      ! has_* flags HOST-side, never from a device reduction.  They stay
      ! latched for the rest of the run — the melt rate can legitimately
      ! pass through zero, and un-latching there would change the
      ! apply-path operand count mid-run.
      sf%has_heat = .true.
      sf%has_salt = .true.

      call cavity_status_counts_impl(nx, ny, cav%active, cav%status, &
                                     cav%n_nonfinite_step, cav%n_not_converged_step, &
                                     cav%n_no_root_step, cav%n_other_step)
      cav%n_not_converged_total = cav%n_not_converged_total &
                                  + int(cav%n_not_converged_step, int64)
      cav%n_no_root_total = cav%n_no_root_total + int(cav%n_no_root_step, int64)
      cav%n_other_total = cav%n_other_total + int(cav%n_other_step, int64)

      if (cavity_melt_status_is_fatal(cav%n_nonfinite_step)) then
         call fail("&ocean_cavity_melt_nml: "//to_string(cav%n_nonfinite_step)// &
                   " ice-covered column(s) fed the basal-melt solver a non-finite "// &
                   "far-field temperature, salinity, velocity or interface pressure "// &
                   "(CAVITY_MELT_NONFINITE_*).  The kernel's safe state is zero melt, "// &
                   "which would hide an already-corrupt column behind a plausible "// &
                   "run, so this is fatal.  NOT_CONVERGED / NO_PHYSICAL_ROOT are "// &
                   "counted and warned instead.", code=OCEAN_STATUS_ERR_SETUP)
      end if

      if (cav%n_not_converged_step > 0 .or. cav%n_no_root_step > 0 .or. &
          cav%n_other_step > 0) then
         call global_logger%warning("cavity melt: zero melt applied on "// &
                                    to_string(cav%n_not_converged_step)// &
                                    " non-converged, "//to_string(cav%n_no_root_step)// &
                                    " no-physical-root and "// &
                                    to_string(cav%n_other_step)// &
                                    " otherwise-refused column(s) this step")
      end if
   end subroutine ocean_cavity_flux_step

end module rdb_ocean_cavity_flux

rdb_ocean_cavity_flux.F90 Source File

Contents

Modules

Source Code

What this module owns

Far-field sampling (why metres)

The three pressures, and which one this is

Sign and units — the table that decides whether the answer is right

Phase 3 — REAL freshwater MASS (`&ocean_cavity_melt_nml freshwater`)

Which density converts kg to m

The three top-layer tendencies

Why the virtual salt flux must NOT also be applied — and why it

is nevertheless still assembled

Agreement with the virtual form, and the size of the difference

Sea level (`volume_compensation`)

Budgets

Budgets

`mem:separate`

Citations (papers, never another model’s source)

This file depends on

Files dependent on this one

Source Code

rdb_ocean_cavity_flux.F90 Source File

Contents

Modules

Source Code

What this module owns

Far-field sampling (why metres)

The three pressures, and which one this is

Sign and units — the table that decides whether the answer is right

Phase 3 — REAL freshwater MASS (&ocean_cavity_melt_nml freshwater)

Which density converts kg to m

The three top-layer tendencies

Why the virtual salt flux must NOT also be applied — and why it

is nevertheless still assembled

Agreement with the virtual form, and the size of the difference

Sea level (volume_compensation)

Budgets

Budgets

mem:separate

Citations (papers, never another model’s source)

This file depends on

Files dependent on this one

Source Code

Phase 3 — REAL freshwater MASS (`&ocean_cavity_melt_nml freshwater`)

Sea level (`volume_compensation`)

`mem:separate`