!! Surface heat + salt fluxes for the ocean dynamical core.  Sister
!! module to `rdb_ocean_surface_stress` — same 2D-field shape, applies
!! a tracer flux at the surface layer (`k = nz` under the ROMS-style
!! k=1-bed convention).
!!
!! Heat: `dT/dt|_surface = Q_heat / (rho_0 * cp * h_top)` with
!! Q_heat in W/m^2 (positive downward into the ocean).  Salt:
!! `dS/dt|_surface = +Q_salt / (rho_0 * h_top)` with `Q_salt` the net
!! virtual salt flux (PSU.kg m-2 s-1, POSITIVE SALINIFIES — evaporation
!! excess, or sea-ice brine rejection via `rdb_ice_ocean_coupler`).
!!
!! The 2D fields `Q_heat(:,:)` / `Q_salt(:,:)` are seeded at configure
!! time from the `&ocean_thermo_nml q_heat / q_salt` scalars via
!! `set_surface_flux_const`.  They live on the device between
!! `enter_data` and `exit_data`; Area-A3 data-override and Area-A4
!! restoring write directly into these fields after seeding.
!!
!! `Q_heat_const` / `Q_salt_const` are kept as the fill-path scalars
!! (mirroring `set_wind_stress_const` on the stress side) so
!! `set_surface_flux_const` produces identical arithmetic to the old
!! scalar broadcast — existing nmls and tests are bit-identical.
!!
!! Phase 2 of the KPP / EPBL build-out also reads the heat/salt fluxes
!! to compute B_0.  Those shimmed paths now pick up `Q_heat(i,j)` /
!! `Q_salt(i,j)` per column instead of the broadcast scalar.
!!
!! Restart: the fields are configure-time-filled static forcing; they
!! are NOT registered in the restart registry (re-seeded from the
!! namelist scalar on every resume — see `rdb_driver.F90` for the
!! seeding call-site).
!!
!! **Component set (PR-12, `&ocean_forcing_nml enable_components`).**
!! `Q_heat` / `Q_salt` are the derived views every downstream kernel
!! keeps reading; when `use_components` is on, `ocean_surface_flux_assemble`
!! rebuilds them every thermo step from `Q_heat_const`/`Q_salt_const`
!! plus the per-component fields below.  Fill contract for any future
!! filler (a reader, the ice coupler, a rivers PR, ...):
!!   (a) write your OWN component(s) — NEVER `Q_heat` / `Q_salt` directly
!!       (those are assembler-owned; a second writer is a review-rejectable
!!       error, see `ocean_surface_flux_assemble`'s docstring);
!!   (b) `!$acc update device(...)` your write (or write on-device) before
!!       the assembler's next call — components are device-resident between
!!       `enter_data`/`exit_data`;
!!   (c) set `has_heat`/`has_salt` (existing contract) and, as needed,
!!       `has_mass_flux` (after any mass-flux write) / `has_q_sw` (after a
!!       `q_sw` write) host-side — never from a device reduction;
!!   (d) fill `heat_content_<flux>` for every mass flux you fill (the
!!       source owns the enthalpy of the mass it injects; the ocean can
!!       only compute the enthalpy of mass it loses, `heat_content_massout`,
!!       which the assembler derives from SST — never write it yourself);
!! **Ice-shelf cover (`&ocean_cavity_dyn_nml`).**  Under a shelf there
!! is no atmosphere, so every ATMOSPHERIC contribution must be zero in a
!! covered cell and unchanged everywhere else, while the cavity's own
!! `heat_cavity` / `salt_cavity` must NOT be masked.  The mask is an
!! OPTIONAL `cover_frac` argument (absent ⇒ the original kernel,
!! byte-identical) and it is applied in the **assembler**, not at apply
!! time.  That choice is load-bearing, for two reasons:
!!
!!   1. `Q_heat` / `Q_salt` are not only the apply kernel's input — they
!!      are what KPP and EPBL read to build `B_0`.  Masking at apply
!!      time would leave both boundary-layer schemes forced by an
!!      atmosphere that is not there, with the tracer deposit correct
!!      and the mixing wrong: a plausible, publishable, wrong answer.
!!   2. The assembler is the ONE place where the atmospheric bands and
!!      the cavity bands are still distinguishable.  After it, `Q_heat`
!!      is a sum, and any factor applied to the sum would also scale
!!      `heat_cavity` — i.e. mask away the melt flux the cover is
!!      supposed to admit.
!!
!! The two surface kernels that do NOT route through `Q_heat`/`Q_salt`
!! carry the same optional argument and mask themselves:
!! `ocean_surface_flux_apply_sw_penetration` (it reads `q_sw`, a
!! pristine INPUT component, on the `sw_source="q_sw"` branch — and even
!! on the `net_heat` branch it must not *move* heat that the masked
!! deposit never added) and `ocean_surface_restore_apply_tracers` (it
!! forms its flux in-kernel from the live SST/SSS).  With the component
!! set OFF there is no assembler, so the static `Q_heat`/`Q_salt` fill is
!! masked once at configure by `ocean_surface_flux_apply_cover_const`.
!!
!!   (e) register your component in the restart registry if it is
!!       time-varying (`registry_register_2d`, `optional=.true.`,
!!       `device_mapped=.true.`) — the assembler's own outputs
!!       (`Q_heat`/`Q_salt`, `heat_content_massin/massout`) are never
!!       registered (derived-field rule, `rdb_ocean_state.F90`).
module rdb_ocean_surface_flux
   use rdb_constants, only: wp
   use rdb_grid, only: hgrid_t
   use rdb_multilayer_state, only: multilayer_state_t
   use, intrinsic :: iso_fortran_env, only: int64
   use rdb_mem_report, only: arr_bytes
   implicit none
   private

   public :: ocean_surface_flux_t
   public :: ocean_surface_flux_apply_tracers
   public :: ocean_surface_flux_apply_sw_penetration
   public :: ocean_surface_restore_apply_tracers
   public :: ocean_surface_flux_assemble
   public :: ocean_surface_flux_apply_cover_const
   public :: sw_transmission
   public :: sw_pe_cost_shape
   public :: sw_source_is_implemented

   ! Seconds per day — converts the MOM6 `FLUXCONST` piston velocity
   ! (specified in m/day) to MKS m/s at the seed call.
   real(wp), parameter :: SECONDS_PER_DAY = 86400.0_wp

   ! Specific heat capacity of seawater (J/kg/K).  Wright (1997)
   ! Table A1 reference value; same as MOM6's `CP_SW` default.
   real(wp), parameter, public :: SEAWATER_CP = 3992.0_wp

   type :: ocean_surface_flux_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.
      logical :: has_heat = .false.
         !! True when `Q_heat` carries a non-zero fill (set by
         !! `set_surface_flux_const` when `q_heat_val /= 0`).
         !! Any future field-fill path (A3 data-override, A4 restoring)
         !! MUST set `has_heat = .true.` after writing into `Q_heat`
         !! so the apply-tracers kernel fires.  Host-side flag only
         !! (early-return guard in `ocean_surface_flux_apply_tracers`).
      logical :: has_salt = .false.
         !! True when `Q_salt` carries a non-zero fill (set by
         !! `set_surface_flux_const` when `q_salt_val /= 0`).
         !! Same contract as `has_heat` for any field-fill path.
      real(wp) :: rho0 = 1035.0_wp
         !! Boussinesq reference density (kg/m^3) — the `dt/(rho0*cp)` heat
         !! and `dt/rho0` salt divisors applied to EVERY surface tracer
         !! source, including whatever the sea-ice coupler writes into
         !! `Q_heat`/`Q_salt`.
         !!
         !! ASSIGNED FROM CONFIG by `configure_ocean_reference_density`,
         !! which copies the single rho0 of record (`&ocean_ic_nml rho_0`
         !! -> `eos%rho0`).  The literal here is only the pre-configure
         !! type default; do not read it as the value a run uses.  Host
         !! scalar: the divisor is folded into the `inv_scale` argument on
         !! the host, so the assignment owes no `!$acc update device`.
      real(wp) :: cp = SEAWATER_CP
         !! Specific heat capacity (J/kg/K).
      real(wp) :: h_min = 1.0e-3_wp
         !! Floor on the surface-layer thickness in the `1/h_top`
         !! division — keeps the kernel finite when the top layer
         !! pinches out.
      real(wp) :: Q_heat_const = 0.0_wp
         !! Scalar fill source for `Q_heat(:,:)`.  Seeded from
         !! `&ocean_thermo_nml q_heat` by `set_surface_flux_const`.
         !! Kept for diagnostic logging; kernels read `Q_heat` directly.
      real(wp) :: Q_salt_const = 0.0_wp
         !! Scalar fill source for `Q_salt(:,:)`.  Seeded from
         !! `&ocean_thermo_nml q_salt` by `set_surface_flux_const`.

      logical :: has_sw = .false.
         !! True when shortwave penetration is active (set by
         !! `set_sw_penetration` when `sw_pen_frac /= 0`).  Host-side
         !! gate only — `ocean_surface_flux_apply_sw_penetration`
         !! early-returns unless this is set, so the default-off path
         !! leaves the surface-flux deposition byte-for-byte unchanged.
      logical :: sw_from_qsw = .false.
         !! Selects the irradiance source for shortwave penetration and
         !! the boundary-layer SW coupling.  `.false.` (default): the
         !! source is the NET heat flux `Q_heat` (`sw_source="net_heat"`,
         !! legacy, bit-identical).  `.true.` (`sw_source="q_sw"`): the
         !! source is the dedicated `q_sw` component (>= 0), which
         !! removes the night-time negative-`I0` hazard where
         !! `sw_pen_frac*Q_heat < 0` drives unphysical negative
         !! irradiance down the two-band profile.  Set by
         !! `set_sw_penetration`; host-side gate only (never a device
         !! reduction) — the shim selects the source array on the host so
         !! the conditionally-allocated `q_sw` is never dereferenced in a
         !! device kernel.
      real(wp) :: sw_pen_frac = 0.0_wp
         !! Penetrating fraction of `Q_heat` carried below the surface
         !! layer as a two-band exponential (Paulson & Simpson 1977).
         !! 0 = off (all of `Q_heat` lands at `k = nz`, legacy path).
      real(wp) :: sw_band_ratio = 0.58_wp
         !! Band-1 weight `R` of the two-band irradiance decay.  Jerlov
         !! type I (clear open ocean) default.
      real(wp) :: sw_zeta1 = 0.35_wp
         !! Band-1 e-folding depth (m) — the rapidly-absorbed
         !! red/near-IR band.
      real(wp) :: sw_zeta2 = 23.0_wp
         !! Band-2 e-folding depth (m) — the slowly-absorbed
         !! blue/green band.
      logical :: has_restore_T = .false.
         !! True when SST restoring is active (`enable_restore_temp .and.
         !! restore_piston_T /= 0`).  Host-side gate only —
         !! `ocean_surface_restore_apply_tracers` early-returns unless
         !! this or `has_restore_S` is set, so the default-off path is
         !! byte-for-byte unchanged.
      logical :: has_restore_S = .false.
         !! True when SSS restoring is active (`enable_restore_salt .and.
         !! restore_piston_S /= 0`).
      real(wp) :: restore_piston_T = 0.0_wp
         !! SST piston velocity (m/s), seeded from `&ocean_restore_nml
         !! piston_t` (m/day) via `/86400`.  The surface relaxation rate
         !! for a top layer of thickness `h_top` is
         !! `lambda = restore_piston_T / h_top` [1/s].
      real(wp) :: restore_piston_S = 0.0_wp
         !! SSS piston velocity (m/s).
      real(wp) :: restore_T_target = 0.0_wp
         !! Scalar target SST (degC).  Read by-value into the device
         !! `_impl` kernel — no per-cell field (so no extra device
         !! array, the `enter_data` orchestrator is untouched).  A
         !! 2D-field target is the documented A4-v2 follow-up.
      real(wp) :: restore_S_target = 0.0_wp
         !! Scalar target SSS (PSU).
      real(wp), allocatable :: Q_heat(:, :)
         !! 2D net surface heat flux (W/m^2, positive downward),
         !! shape `(nx, ny)`.  Fill via `set_surface_flux_const` for
         !! spatially-uniform forcing (the default); Area-A3 override
         !! or Area-A4 restoring writes the field directly.
      real(wp), allocatable :: Q_salt(:, :)
         !! 2D net surface salt flux (kg salt/m^2/s, positive salinifies),
         !! shape `(nx, ny)`.

      ! ---- PR-12 component set — allocated iff `use_components` ----
      logical :: use_components = .false.
         !! Master gate (`&ocean_forcing_nml enable_components`).  `.false.`
         !! (default): none of the arrays below are allocated,
         !! `ocean_surface_flux_assemble` is a no-op, and `Q_heat`/`Q_salt`
         !! are filled exactly as today — bit-identical.  `.true.`:
         !! allocates the component set (`set_components`) and the
         !! assembler rebuilds `Q_heat`/`Q_salt` every thermo step.
      real(wp) :: q_assembled = 0.0_wp
         !! Host latch, 1 once `ocean_surface_flux_assemble` has derived
         !! `Q_heat`/`Q_salt` from the component set.  From then on the
         !! two arrays are CARRIED state: the assembler runs at the END of
         !! a thermo step and the steps up to the next one read what it
         !! left (SST-dependent terms included), so with `use_components`
         !! they are checkpointed, and this latch tells a warm restart
         !! that the checkpointed arrays are an assembly to resume from
         !! (0 => an older checkpoint or none yet: the configure-time
         !! seed stands).  A real, not a logical, because the restart
         !! registry carries real scalars.
      logical :: has_mass_flux = .false.
         !! Host-side latch — set by a filler after writing ANY of
         !! `evap`/`lprec`/`fprec`/`vprec`/`lrunoff`/`frunoff`/
         !! `seaice_melt`.  Cheap early-return gate for a future
         !! freshwater kernel (real-mass PR).  Same contract as
         !! `has_heat` (`:56-62`) — never set from a device reduction.
      logical :: has_q_sw = .false.
         !! Host-side latch — set by a filler after writing `q_sw`
         !! (e.g. an ice `sw_thru` coupler).  **Not** the same as
         !! `has_sw` below (that gates shortwave *penetration*, an
         !! unrelated pre-existing switch) — do not conflate the two.

      real(wp), allocatable :: q_sw(:, :)
         !! Shortwave into the ocean (W/m^2, **>= 0**, positive down).
      real(wp), allocatable :: q_lw(:, :)
         !! Net longwave (W/m^2, typically **< 0**, positive down).
      real(wp), allocatable :: q_lat(:, :)
         !! Latent heat flux (W/m^2, typically **< 0**, positive down).
      real(wp), allocatable :: q_sens(:, :)
         !! Sensible heat flux (W/m^2, typically **< 0**, positive down).
      real(wp), allocatable :: heat_added(:, :)
         !! Restoring / flux-adjustment / "other" net heat term not
         !! decomposed into the radiative/turbulent bands above (W/m^2,
         !! either sign; MOM6 `heat_added`).  The v1 ice coupler's
         !! `heat_flux_diag` lands here (§5.4 of the PR-12 plan) —
         !! it is already a net W/m^2, not further decomposable.
      real(wp), allocatable :: heat_cavity(:, :)
         !! **Ice-shelf cavity basal-melt heat component** (W/m^2, same
         !! positive-DOWN-into-the-ocean convention as every other heat
         !! band; `&ocean_cavity_melt_nml`).  OWNED by
         !! `rdb_ocean_cavity_flux`; written `heat_cavity = -q_ocean`,
         !! where `q_ocean = rho_w*c_w*gamma_t*(T_w - T_b) > 0` is the
         !! kernel's turbulent heat flux OCEAN -> INTERFACE, so warm
         !! water under a shelf COOLS the top of the column.  It is a
         !! SEPARATE field from `heat_added` precisely because the
         !! sea-ice coupler full-overwrites `heat_added` — two writers
         !! on one slot clobber silently (cavity x sea ice is refused
         !! today, but the ownership rule must not depend on that).
         !! Zero unless a cavity melt step ran.

      real(wp), allocatable :: evap(:, :)
         !! Evaporative mass flux (kg/m^2/s, **<= 0** — MOM6 convention,
         !! `(-1)*flux out of the ocean`).  v1: enthalpy + salt
         !! bookkeeping only — does NOT change column mass (real
         !! freshwater is a named follow-up, see the module docstring).
      real(wp), allocatable :: lprec(:, :)
         !! Liquid precipitation (kg/m^2/s, **>= 0** into the ocean).
      real(wp), allocatable :: fprec(:, :)
         !! Frozen precipitation / snowfall (kg/m^2/s, **>= 0**).
      real(wp), allocatable :: vprec(:, :)
         !! Virtual precipitation (kg/m^2/s, either sign — SSS-restoring
         !! convention; NOT wired to the restoring kernel in v1, see
         !! `set_restore` below and the module docstring's follow-up note).
      real(wp), allocatable :: lrunoff(:, :)
         !! Liquid river runoff (kg/m^2/s, **>= 0**).
      real(wp), allocatable :: frunoff(:, :)
         !! Frozen (calving/ice) runoff (kg/m^2/s, **>= 0**).
      real(wp), allocatable :: seaice_melt(:, :)
         !! Sea-ice melt-water mass flux (kg/m^2/s, **>0** = melt into
         !! the ocean, **<0** = formation / freezing withdraws mass).

      real(wp), allocatable :: heat_content_lprec(:, :)
         !! Enthalpy carried by `lprec` (W/m^2).  A filler that writes
         !! `lprec` MUST fill this — the v1 convenience is
         !! `SEAWATER_CP * T_source * lprec`.  The source (not the
         !! ocean) owns this enthalpy — see the module docstring §(d).
      real(wp), allocatable :: heat_content_fprec(:, :)
         !! Enthalpy carried by `fprec` (W/m^2).
      real(wp), allocatable :: heat_content_vprec(:, :)
         !! Enthalpy carried by `vprec` (W/m^2).
      real(wp), allocatable :: heat_content_lrunoff(:, :)
         !! Enthalpy carried by `lrunoff` (W/m^2).
      real(wp), allocatable :: heat_content_frunoff(:, :)
         !! Enthalpy carried by `frunoff` (W/m^2).
      real(wp), allocatable :: heat_content_seaice_melt(:, :)
         !! Enthalpy carried by `seaice_melt` (W/m^2).
      real(wp), allocatable :: heat_content_massin(:, :)
         !! **Assembler output — do NOT write.**  Sum of the six
         !! `heat_content_<flux>` companions above (W/m^2, >= 0 for warm
         !! inflow).  Filled by `ocean_surface_flux_assemble`.
      real(wp), allocatable :: heat_content_massout(:, :)
         !! **Assembler output — do NOT write.**  `SEAWATER_CP * T_sst *
         !! evap` (W/m^2, <= 0 since `evap <= 0`) — the enthalpy the ocean
         !! loses with evaporating mass, computed from the ocean's own
         !! surface temperature (the ocean, not a filler, owns this
         !! number).  There is deliberately NO `heat_content_evap` field
         !! — see the PR-12 plan §11.6.  Filled by
         !! `ocean_surface_flux_assemble`.

      real(wp), allocatable :: salt_flux(:, :)
         !! Net surface salt-flux COMPONENT (kg salt/m^2/s, **positive
         !! salinifies**) — a filler writes this (e.g. the ice brine
         !! coupler); the assembler adds `Q_salt_const` to produce
         !! `Q_salt`.  Virtual in v1 (no column-mass change).
      real(wp), allocatable :: salt_cavity(:, :)
         !! **Ice-shelf cavity basal-melt salt component**, same units
         !! and sign as `salt_flux` (positive salinifies;
         !! `&ocean_cavity_melt_nml`).  OWNED by `rdb_ocean_cavity_flux`
         !! and never written by the ice coupler, which full-overwrites
         !! `salt_flux`.
         !!
         !! **The fixed-mass dilution equivalent.**
         !!
         !!   `salt_cavity = -m_mass*(S_far - s_ice)`
         !!
         !! is the exact fixed-mass equivalent of adding mass `m_mass` at
         !! salinity `s_ice` — see the derivation in
         !! `rdb_ocean_cavity_flux`'s module docstring.  Melting
         !! (`m_mass > 0`, `S_far > s_ice`) therefore FRESHENS.
         !!
         !! Under `&ocean_cavity_melt_nml freshwater="virtual"` (the
         !! default) that IS the meltwater's whole effect: no mass moves.
         !!
         !! Under `freshwater="mass"` the meltwater is a REAL volume
         !! source on the top layer and this component is NOT the
         !! salinity tendency any more — but it is STILL assembled into
         !! `Q_salt` unchanged, because `Q_salt` is also what KPP and
         !! EPBL read to build `B_0`, and this term is the dominant
         !! (freshening) part of the surface buoyancy flux there.
         !! `ocean_cavity_mass_step` takes the increment back out of the
         !! SALINITY TRACER (and out of the pseudo-salt mirror) in the
         !! same stage, as the exact negation of what
         !! `apply_surface_src_2d_impl` stamped.  So: one field, two
         !! readers, and only the tracer reader is corrected.

      real(wp), allocatable :: p_surf_atm(:, :)
         !! **Input component.**  Atmospheric surface-pressure load
         !! (Pa, >= 0).  Filled by an external reader / configure-time
         !! scalar seed — the sea-ice path never writes this field.
         !! Ships zeroed with no consumer in this PR (the inverse-
         !! barometer PGF fold is a same-release-cycle follow-up).
      real(wp), allocatable :: p_surf(:, :)
         !! **Assembled total** (Pa, >= 0) — `p_surf_atm` plus any ice
         !! mass-loading term, **full overwrite, never `+=`** (a `+=`
         !! ratchets the load across outer steps with no bound).  No
         !! consumer in this PR; ships zeroed alongside `p_surf_atm` so
         !! the follow-up PGF fold needs no further plumbing.
   contains
      procedure, non_overridable :: init => ocean_surfflux_init
      procedure, non_overridable :: destroy => ocean_surfflux_destroy
      procedure, non_overridable :: enter_data => ocean_surfflux_enter_data
      procedure, non_overridable :: exit_data => ocean_surfflux_exit_data
      procedure, non_overridable :: set_surface_flux_const => ocean_surfflux_set_const
      procedure, non_overridable :: set_sw_penetration => ocean_surfflux_set_sw
      procedure, non_overridable :: set_restore => ocean_surfflux_set_restore
      procedure, non_overridable :: set_components => ocean_surfflux_set_components
      procedure, non_overridable :: set_p_surf_const => ocean_surfflux_set_p_surf_const
      procedure, non_overridable :: bytes => ocean_surface_flux_bytes
   end type ocean_surface_flux_t

contains

   subroutine ocean_surfflux_init(this, grid)
      class(ocean_surface_flux_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      integer :: nx, ny
      nx = grid%nx_total
      ny = grid%ny_total
      allocate (this%Q_heat(nx, ny), source=0.0_wp)
      allocate (this%Q_salt(nx, ny), source=0.0_wp)
      this%q_assembled = 0.0_wp
      this%is_init = .true.
   end subroutine ocean_surfflux_init

   subroutine ocean_surfflux_destroy(this)
      class(ocean_surface_flux_t), intent(inout) :: this
      this%is_init = .false.
      this%q_assembled = 0.0_wp
      if (allocated(this%Q_heat)) deallocate (this%Q_heat)
      if (allocated(this%Q_salt)) deallocate (this%Q_salt)
      call ocean_surfflux_dealloc_components(this)
   end subroutine ocean_surfflux_destroy

   subroutine ocean_surfflux_dealloc_components(this)
      !! Deallocate the component set (no-op on an already-unallocated
      !! slot — every `deallocate` is `if (allocated(...))`-guarded).
      !! Shared by `destroy` and by `set_components` re-entry.
      class(ocean_surface_flux_t), intent(inout) :: this
      if (allocated(this%q_sw)) deallocate (this%q_sw)
      if (allocated(this%q_lw)) deallocate (this%q_lw)
      if (allocated(this%q_lat)) deallocate (this%q_lat)
      if (allocated(this%q_sens)) deallocate (this%q_sens)
      if (allocated(this%heat_added)) deallocate (this%heat_added)
      if (allocated(this%heat_cavity)) deallocate (this%heat_cavity)
      if (allocated(this%evap)) deallocate (this%evap)
      if (allocated(this%lprec)) deallocate (this%lprec)
      if (allocated(this%fprec)) deallocate (this%fprec)
      if (allocated(this%vprec)) deallocate (this%vprec)
      if (allocated(this%lrunoff)) deallocate (this%lrunoff)
      if (allocated(this%frunoff)) deallocate (this%frunoff)
      if (allocated(this%seaice_melt)) deallocate (this%seaice_melt)
      if (allocated(this%heat_content_lprec)) deallocate (this%heat_content_lprec)
      if (allocated(this%heat_content_fprec)) deallocate (this%heat_content_fprec)
      if (allocated(this%heat_content_vprec)) deallocate (this%heat_content_vprec)
      if (allocated(this%heat_content_lrunoff)) deallocate (this%heat_content_lrunoff)
      if (allocated(this%heat_content_frunoff)) deallocate (this%heat_content_frunoff)
      if (allocated(this%heat_content_seaice_melt)) deallocate (this%heat_content_seaice_melt)
      if (allocated(this%heat_content_massin)) deallocate (this%heat_content_massin)
      if (allocated(this%heat_content_massout)) deallocate (this%heat_content_massout)
      if (allocated(this%salt_flux)) deallocate (this%salt_flux)
      if (allocated(this%salt_cavity)) deallocate (this%salt_cavity)
      if (allocated(this%p_surf_atm)) deallocate (this%p_surf_atm)
      if (allocated(this%p_surf)) deallocate (this%p_surf)
   end subroutine ocean_surfflux_dealloc_components

   subroutine ocean_surfflux_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 (AMD libomptarget cross-slot-overlap fix).
      class(ocean_surface_flux_t), intent(inout) :: this
      select type (this)
      type is (ocean_surface_flux_t)
         call ocean_surfflux_enter_data_impl(this)
      end select
   end subroutine ocean_surfflux_enter_data

   subroutine ocean_surfflux_enter_data_impl(this)
      type(ocean_surface_flux_t), intent(inout) :: this
      !$acc enter data copyin(this%Q_heat, this%Q_salt)
      !$acc update device(this%Q_heat, this%Q_salt)
      if (this%use_components) then
         !$acc enter data copyin(this%q_sw, this%q_lw, this%q_lat, this%q_sens, &
         !$acc&                  this%heat_added, this%heat_cavity, this%evap, &
         !$acc&                  this%lprec, this%fprec, &
         !$acc&                  this%vprec, this%lrunoff, this%frunoff, this%seaice_melt, &
         !$acc&                  this%heat_content_lprec, this%heat_content_fprec, &
         !$acc&                  this%heat_content_vprec, this%heat_content_lrunoff, &
         !$acc&                  this%heat_content_frunoff, this%heat_content_seaice_melt, &
         !$acc&                  this%heat_content_massin, this%heat_content_massout, &
         !$acc&                  this%salt_flux, this%salt_cavity, &
         !$acc&                  this%p_surf_atm, this%p_surf)
         !$acc update device(this%q_sw, this%q_lw, this%q_lat, this%q_sens, &
         !$acc&               this%heat_added, this%heat_cavity, this%evap, &
         !$acc&               this%lprec, this%fprec, &
         !$acc&               this%vprec, this%lrunoff, this%frunoff, this%seaice_melt, &
         !$acc&               this%heat_content_lprec, this%heat_content_fprec, &
         !$acc&               this%heat_content_vprec, this%heat_content_lrunoff, &
         !$acc&               this%heat_content_frunoff, this%heat_content_seaice_melt, &
         !$acc&               this%heat_content_massin, this%heat_content_massout, &
         !$acc&               this%salt_flux, this%salt_cavity, &
         !$acc&               this%p_surf_atm, this%p_surf)
      end if
   end subroutine ocean_surfflux_enter_data_impl

   subroutine ocean_surfflux_exit_data(this)
      class(ocean_surface_flux_t), intent(inout) :: this
      select type (this)
      type is (ocean_surface_flux_t)
         call ocean_surfflux_exit_data_impl(this)
      end select
   end subroutine ocean_surfflux_exit_data

   subroutine ocean_surfflux_exit_data_impl(this)
      type(ocean_surface_flux_t), intent(inout) :: this
      if (this%use_components) then
         !$acc exit data delete(this%q_sw, this%q_lw, this%q_lat, this%q_sens, &
         !$acc&                 this%heat_added, this%heat_cavity, this%evap, &
         !$acc&                 this%lprec, this%fprec, &
         !$acc&                 this%vprec, this%lrunoff, this%frunoff, this%seaice_melt, &
         !$acc&                 this%heat_content_lprec, this%heat_content_fprec, &
         !$acc&                 this%heat_content_vprec, this%heat_content_lrunoff, &
         !$acc&                 this%heat_content_frunoff, this%heat_content_seaice_melt, &
         !$acc&                 this%heat_content_massin, this%heat_content_massout, &
         !$acc&                 this%salt_flux, this%salt_cavity, &
         !$acc&                 this%p_surf_atm, this%p_surf)
      end if
      !$acc exit data delete(this%Q_heat, this%Q_salt)
   end subroutine ocean_surfflux_exit_data_impl

   subroutine ocean_surfflux_set_const(this, q_heat_val, q_salt_val)
      !! Fill `Q_heat` / `Q_salt` uniformly from scalar values and set
      !! the `has_heat` / `has_salt` flags so the apply-tracers kernel
      !! fires.  Mirrors `set_wind_stress_const` on the stress side.
      !! Host only — call `enter_data` afterwards (or `!$acc update
      !! device` if already mapped) to sync to the GPU.
      class(ocean_surface_flux_t), intent(inout) :: this
      real(wp), intent(in) :: q_heat_val, q_salt_val
      this%Q_heat_const = q_heat_val
      this%Q_salt_const = q_salt_val
      this%Q_heat = q_heat_val
      this%Q_salt = q_salt_val
      this%has_heat = (q_heat_val /= 0.0_wp)
      this%has_salt = (q_salt_val /= 0.0_wp)
   end subroutine ocean_surfflux_set_const

   subroutine ocean_surfflux_set_sw(this, sw_pen_frac, sw_band_ratio, &
                                    sw_zeta1, sw_zeta2, sw_source)
      !! Seed the shortwave-penetration band parameters and set the
      !! `has_sw` gate (`sw_pen_frac /= 0`).  Sibling to
      !! `set_surface_flux_const` — kept separate so existing callers of
      !! the heat/salt setter are unchanged.  Host only; the scalars are
      !! read host-side by the apply kernel (they parameterise the
      !! by-value arguments passed into the device `_impl`), so no extra
      !! device sync is needed beyond the existing `copyin(this)`.
      !!
      !! `sw_source` (optional; default `"net_heat"`) selects the
      !! irradiance source: `"net_heat"` ⇒ `I0 = sw_pen_frac*Q_heat`
      !! (legacy, bit-identical); `"q_sw"` ⇒ `I0 = sw_pen_frac*q_sw`
      !! (requires the PR-12 component set — the caller / `validate_config`
      !! guards allocation).  Sets the host-side `sw_from_qsw` gate.
      class(ocean_surface_flux_t), intent(inout) :: this
      real(wp), intent(in) :: sw_pen_frac, sw_band_ratio, sw_zeta1, sw_zeta2
      character(len=*), intent(in), optional :: sw_source
      this%sw_pen_frac = sw_pen_frac
      this%sw_band_ratio = sw_band_ratio
      this%sw_zeta1 = sw_zeta1
      this%sw_zeta2 = sw_zeta2
      this%has_sw = (sw_pen_frac /= 0.0_wp)
      this%sw_from_qsw = .false.
      if (present(sw_source)) this%sw_from_qsw = (trim(sw_source) == "q_sw")
   end subroutine ocean_surfflux_set_sw

   pure function sw_source_is_implemented(name) result(ok)
      !! Fail-loud predicate for the `&ocean_thermo_nml sw_source`
      !! selector — `validate_config` aborts on any string this rejects.
      !! The two recognised sources are the net-heat legacy path and the
      !! PR-12 `q_sw` component.
      character(len=*), intent(in) :: name
      logical :: ok
      ok = (trim(name) == "net_heat" .or. trim(name) == "q_sw")
   end function sw_source_is_implemented

   subroutine ocean_surfflux_set_restore(this, enable_T, enable_S, &
                                         piston_t_day, piston_s_day, &
                                         T_target, S_target)
      !! Seed the surface buoyancy restoring (MOM6 `RESTOREBUOY`)
      !! parameters and set the `has_restore_T` / `has_restore_S` gates.
      !! Sibling to `set_surface_flux_const` / `set_sw_penetration` —
      !! kept separate so existing callers are unchanged.  The piston
      !! velocities arrive in **m/day** (the MOM6 `FLUXCONST_*` unit) and
      !! are converted to MKS m/s here.  Effective-enable guard:
      !! `has_restore_* = enable_* .and. piston /= 0`, so an enabled
      !! switch with a zero piston is a silent no-op (rather than
      !! restoring everything toward the 0-degC / 0-PSU default target).
      !! Host only; all knobs are read host-side as by-value arguments to
      !! the device `_impl`, so no extra device sync beyond `copyin(this)`.
      class(ocean_surface_flux_t), intent(inout) :: this
      logical, intent(in) :: enable_T, enable_S
      real(wp), intent(in) :: piston_t_day, piston_s_day
      real(wp), intent(in) :: T_target, S_target
      this%restore_piston_T = piston_t_day/SECONDS_PER_DAY
      this%restore_piston_S = piston_s_day/SECONDS_PER_DAY
      this%restore_T_target = T_target
      this%restore_S_target = S_target
      this%has_restore_T = (enable_T .and. this%restore_piston_T /= 0.0_wp)
      this%has_restore_S = (enable_S .and. this%restore_piston_S /= 0.0_wp)
   end subroutine ocean_surfflux_set_restore

   subroutine ocean_surfflux_set_p_surf_const(this, p_surf_val)
      !! Seed the atmospheric surface-pressure INPUT component
      !! `p_surf_atm` (Pa) uniformly from a scalar namelist value (PR-17
      !! `&ocean_psurf_nml p_surf_const`).  Full overwrite of the pristine
      !! atmospheric base; the assembled total `p_surf` is built from it
      !! once per outer step in `p_surf_update_seam`.  No-op when the
      !! component set is not allocated (`use_components=.false.`) — the
      !! `&ocean_psurf_nml enable` guard in `validate_config` already
      !! requires `enable_components=.true.`, so a live consumer never hits
      !! the no-op.  Host only — call `enter_data` afterwards (or `!$acc
      !! update device` if already mapped) to sync to the GPU.
      !!
      !! Also seeds the assembled total `p_surf` to the same value — with no
      !! ice mass-loading `p_surf == p_surf_atm` and the atmospheric base is
      !! static, so this configure-time seed IS the assembly (PR-17 reads
      !! `p_surf` read-only in the dyn step).  PR-18's ice path overwrites
      !! `p_surf` per outer step via its own `inout` access.
      class(ocean_surface_flux_t), intent(inout) :: this
      real(wp), intent(in) :: p_surf_val
      if (.not. allocated(this%p_surf_atm)) return
      this%p_surf_atm = p_surf_val
      this%p_surf = p_surf_val
   end subroutine ocean_surfflux_set_p_surf_const

   subroutine ocean_surfflux_set_components(this, grid, enable)
      !! Configure-time gate for the PR-12 component set
      !! (`&ocean_forcing_nml enable_components`).  `init` runs before the
      !! namelist gate is known, so allocation happens HERE rather than in
      !! `init`: `enable = .false.` (default) leaves `use_components`
      !! false and allocates nothing — bit-identical, zero extra device
      !! memory.  `enable = .true.` allocates the full component set
      !! (`source=0.0_wp`) and flips the gate so
      !! `ocean_surface_flux_assemble` stops early-returning.  Must be
      !! called BEFORE `enter_data` (`rdb_ocean_state.F90`'s orchestrator)
      !! so the freshly-allocated arrays get mapped.  Re-entrant: calling
      !! again with a different `enable` deallocates first.
      class(ocean_surface_flux_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      logical, intent(in) :: enable
      call ocean_surfflux_dealloc_components(this)
      this%use_components = enable
      if (enable) call ocean_surfflux_alloc_components(this, grid)
   end subroutine ocean_surfflux_set_components

   subroutine ocean_surfflux_alloc_components(this, grid)
      !! Allocate the 22-field component set + the two `p_surf*` fields,
      !! all `source=0.0_wp`, shape `(nx_total, ny_total)`.  Private —
      !! called only from `set_components`.
      class(ocean_surface_flux_t), intent(inout) :: this
      type(hgrid_t), intent(in) :: grid
      integer :: nx, ny
      nx = grid%nx_total
      ny = grid%ny_total
      allocate (this%q_sw(nx, ny), source=0.0_wp)
      allocate (this%q_lw(nx, ny), source=0.0_wp)
      allocate (this%q_lat(nx, ny), source=0.0_wp)
      allocate (this%q_sens(nx, ny), source=0.0_wp)
      allocate (this%heat_added(nx, ny), source=0.0_wp)
      allocate (this%heat_cavity(nx, ny), source=0.0_wp)
      allocate (this%evap(nx, ny), source=0.0_wp)
      allocate (this%lprec(nx, ny), source=0.0_wp)
      allocate (this%fprec(nx, ny), source=0.0_wp)
      allocate (this%vprec(nx, ny), source=0.0_wp)
      allocate (this%lrunoff(nx, ny), source=0.0_wp)
      allocate (this%frunoff(nx, ny), source=0.0_wp)
      allocate (this%seaice_melt(nx, ny), source=0.0_wp)
      allocate (this%heat_content_lprec(nx, ny), source=0.0_wp)
      allocate (this%heat_content_fprec(nx, ny), source=0.0_wp)
      allocate (this%heat_content_vprec(nx, ny), source=0.0_wp)
      allocate (this%heat_content_lrunoff(nx, ny), source=0.0_wp)
      allocate (this%heat_content_frunoff(nx, ny), source=0.0_wp)
      allocate (this%heat_content_seaice_melt(nx, ny), source=0.0_wp)
      allocate (this%heat_content_massin(nx, ny), source=0.0_wp)
      allocate (this%heat_content_massout(nx, ny), source=0.0_wp)
      allocate (this%salt_flux(nx, ny), source=0.0_wp)
      allocate (this%salt_cavity(nx, ny), source=0.0_wp)
      allocate (this%p_surf_atm(nx, ny), source=0.0_wp)
      allocate (this%p_surf(nx, ny), source=0.0_wp)
   end subroutine ocean_surfflux_alloc_components

   subroutine ocean_surface_flux_apply_tracers(grid, sf, ms, dt, active, wet_dyn)
      !! Add the surface heat + salt fluxes directly to the top
      !! tracer layer.  Operates in `hTr` space (concentration·
      !! thickness): for temperature
      !!   d(hT_top)/dt = Q_heat(i,j) / (rho_0 · cp)
      !! For salinity
      !!   d(hS_top)/dt = Q_salt(i,j) / rho_0
      !! (Both expressed in units that match the `hTr` convention:
      !! `hTr = T·h` so the forcing has units of T·h/s = K·m/s.
      !! Q_heat / (rho_0·cp) has units (W/m^2)/(kg/m^3·J/kg/K) =
      !! K·m/s ✓.)
      !!
      !! Reads the 2D `Q_heat(:,:)` / `Q_salt(:,:)` fields per column
      !! (seeded uniformly from the scalar knobs by default; overwritten
      !! pointwise by Area-A3 data-override or Area-A4 restoring).
      !! No-op when `sf%has_heat` / `sf%has_salt` are false (set by
      !! `set_surface_flux_const` when the fill value is non-zero; any
      !! field-override path must set them before calling), or when no
      !! temperature / salinity tracer is registered.
      type(hgrid_t), intent(in) :: grid
      type(ocean_surface_flux_t), intent(in), optional :: sf
         !! Optional — when absent the kernel is a no-op (no surface
         !! forcing configured).
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      logical, intent(in), optional :: active
         !! Optional gate (thermo cadence).  Absent ⇒ kernel runs;
         !! present-and-false ⇒ early return.
      real(wp), intent(in), optional :: wet_dyn(:, :)
         !! Optional DYNAMIC cell wet mask (wet/dry,
         !! docs/ocean_wetdry_plan.md §4.4) composed multiplicatively
         !! with the static `ms%wet_mask` — surface fluxes must not
         !! enter a dynamically dry column (heating a mm-scale residual
         !! sliver blows its temperature up).  Absent ⇒ the original
         !! static-mask path, byte-identical.

      integer :: nx, ny, nz, idx_T, idx_S, idx_ps

      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. present(sf)) return
      ! Guard: has_heat / has_salt flags unset → skip (set by
      ! set_surface_flux_const or any field-override path).
      if (.not. sf%has_heat .and. .not. sf%has_salt) return
      if (.not. allocated(ms%tracers)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      idx_T = ms%idx_temperature
      idx_S = ms%idx_salinity
      idx_ps = ms%idx_pseudo_salt

      ! Cell-centred wet mask: heat / salt flux only enters ocean cells.
      ! Land cells (mask = 0) accumulate nothing — including in the
      ! budget contributor, so the conservation residual stays clean.
      ! The shim+_impl split keeps `tracers(idx)%hTr` deref on the host
      ! (array-of-DT registry indirection blocks NVHPC device codegen).
      ! The division Q(i,j)/(rho0*cp) [or Q(i,j)/rho0] is applied as
      ! the inv_scale multiplier inside the kernel.
      if (present(wet_dyn)) then
         if (idx_T > 0 .and. sf%has_heat) then
            call apply_surface_src_2d_dyn_impl(ms%tracers(idx_T)%hTr, &
                                               ms%heat_budget_surface, &
                                               ms%wet_mask, wet_dyn, sf%Q_heat, &
                                               dt/(sf%rho0*sf%cp), ms%k_top, nz, nx, ny)
         end if
         if (idx_S > 0 .and. sf%has_salt) then
            call apply_surface_src_2d_dyn_impl(ms%tracers(idx_S)%hTr, &
                                               ms%salt_budget_surface, &
                                               ms%wet_mask, wet_dyn, sf%Q_salt, &
                                               dt/sf%rho0, ms%k_top, nz, nx, ny)
         end if
         ! Pseudo-salt mirror: exactly salinity's surface salt flux,
         ! but through the NOBUDGET twin — budget_id = NONE so it must
         ! not add into salt_budget_surface (§5.6 / test
         ! pseudo_salt_no_budget_contribution).
         if (idx_ps > 0 .and. sf%has_salt) then
            call apply_surface_src_2d_dyn_nobudget_impl(ms%tracers(idx_ps)%hTr, &
                                                        ms%wet_mask, wet_dyn, sf%Q_salt, &
                                                        dt/sf%rho0, ms%k_top, nz, nx, ny)
         end if
         return
      end if
      if (idx_T > 0 .and. sf%has_heat) then
         call apply_surface_src_2d_impl(ms%tracers(idx_T)%hTr, &
                                        ms%heat_budget_surface, &
                                        ms%wet_mask, sf%Q_heat, &
                                        dt/(sf%rho0*sf%cp), ms%k_top, nz, nx, ny)
      end if
      if (idx_S > 0 .and. sf%has_salt) then
         call apply_surface_src_2d_impl(ms%tracers(idx_S)%hTr, &
                                        ms%salt_budget_surface, &
                                        ms%wet_mask, sf%Q_salt, &
                                        dt/sf%rho0, ms%k_top, nz, nx, ny)
      end if
      ! Pseudo-salt mirror (plain, static-mask path) — see note above.
      if (idx_ps > 0 .and. sf%has_salt) then
         call apply_surface_src_2d_nobudget_impl(ms%tracers(idx_ps)%hTr, &
                                                 ms%wet_mask, sf%Q_salt, &
                                                 dt/sf%rho0, ms%k_top, nz, nx, ny)
      end if
   end subroutine ocean_surface_flux_apply_tracers

   pure subroutine apply_surface_src_2d_impl(hTr, budget, wet_mask, Q_field, &
                                             inv_scale, k_top, nz, nx, ny)
      !! Stamp `inv_scale · Q_field(i,j) · wet_mask(i,j)` onto the first
      !! LIVE layer (`k_top(i,j)`) of a tracer's hTr array, mirror into
      !! the matching budget contributor.  Explicit-shape dummies so
      !! NVHPC stdpar can compile device kernels against static bounds.
      !!
      !! **Why `k_top` and not `nz`.** Under a quasi-geopotential
      !! coordinate beneath an ice shelf the layers inside the draft are
      !! inert fillers, so on a covered column `k = nz` carries
      !! `zstar_h_min` of water.  A flux stamped there is neither
      !! diffused down (the vdiff tracer matrix decouples a vanished row
      !! to the identity) nor kept (the next ALE remap drains it on
      !! `h_old <= H_FLOOR`) — but the budget mirror on the line below
      !! still counts it, so the column leaks exactly the deposit, every
      !! thermo step.  The atmospheric part of `Q_field` is zero under
      !! cover, but the cavity's own `heat_cavity`/`salt_cavity` pass
      !! through `ocean_surface_flux_assemble` UNMASKED, which is why
      !! this is the load-bearing site of the whole `k_top` slice.
      !!
      !! `inv_scale` = dt/(rho0·cp) for heat, dt/rho0 for salt — a
      !! column-invariant multiplier that the caller derives from `sf`.
      !! `Q_field` carries any (i,j) spatial variation; for the default
      !! constant-fill case it is uniform, giving arithmetic identical to
      !! the old scalar-broadcast path.
      !!
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: Q_field(nx, ny)
      real(wp), intent(in)    :: inv_scale
      integer, intent(in)    :: k_top(nx, ny)
         !! `ms%k_top` — the first LIVE layer counting down from the top.
         !! `nz` on every column that has no top-side filler (which is
         !! every column on every coordinate but `z_fixed` under a rigid
         !! top), so this reads the same memory as the literal `nz` it
         !! replaced and the arithmetic is bit-identical.
      integer :: i, j
      real(wp) :: cell
      do concurrent(j=1:ny, i=1:nx) local(cell)
         cell = inv_scale*Q_field(i, j)*wet_mask(i, j)
         hTr(i, j, k_top(i, j)) = hTr(i, j, k_top(i, j)) + cell
         budget(i, j, k_top(i, j)) = budget(i, j, k_top(i, j)) + cell
      end do
   end subroutine apply_surface_src_2d_impl

   pure subroutine apply_surface_src_2d_dyn_impl(hTr, budget, wet_mask, wet_dyn, &
                                                 Q_field, inv_scale, k_top, nz, nx, ny)
      !! Wet/dry variant of `apply_surface_src_2d_impl`: the DYNAMIC cell
      !! wet mask composes multiplicatively with the static one, so a
      !! dynamically dry column (total depth below `&ocean_wetdry_nml
      !! dry_depth`) receives NO surface flux — heating a mm-scale
      !! residual sliver would blow its temperature up
      !! (docs/ocean_wetdry_plan.md §4.4).  Separate _impl (not an
      !! in-loop optional test): the knob-off path keeps the original
      !! kernel untouched, byte-identical.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: wet_dyn(nx, ny)
      real(wp), intent(in)    :: Q_field(nx, ny)
      real(wp), intent(in)    :: inv_scale
      integer, intent(in)    :: k_top(nx, ny)
         !! See `apply_surface_src_2d_impl`.
      integer :: i, j
      real(wp) :: cell
      do concurrent(j=1:ny, i=1:nx) local(cell)
         cell = inv_scale*Q_field(i, j)*wet_mask(i, j)*wet_dyn(i, j)
         hTr(i, j, k_top(i, j)) = hTr(i, j, k_top(i, j)) + cell
         budget(i, j, k_top(i, j)) = budget(i, j, k_top(i, j)) + cell
      end do
   end subroutine apply_surface_src_2d_dyn_impl

   pure subroutine apply_surface_src_2d_nobudget_impl(hTr, wet_mask, Q_field, &
                                                      inv_scale, k_top, nz, nx, ny)
      !! Byte-for-byte copy of `apply_surface_src_2d_impl` with the
      !! `budget` dummy and its accumulation line removed — the
      !! pseudo-salt mirror of salinity's surface flux, which by
      !! contract (`budget_id = TRACER_BUDGET_NONE`) must not touch
      !! `salt_budget_surface`.  Separate `_impl`, not an in-loop
      !! `present(budget)` test (house idiom, see
      !! `apply_surface_src_2d_dyn_impl`'s docstring) — this keeps the
      !! production S/T impl untouched and the increment `hTr` receives
      !! bit-identical to salinity's.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: Q_field(nx, ny)
      real(wp), intent(in)    :: inv_scale
      integer, intent(in)    :: k_top(nx, ny)
         !! See `apply_surface_src_2d_impl`.  The SAME index salinity
         !! used, so the pseudo-salt increment stays bit-identical to
         !! salinity's and the deviation keeps measuring transport.
      integer :: i, j
      real(wp) :: cell
      do concurrent(j=1:ny, i=1:nx) local(cell)
         cell = inv_scale*Q_field(i, j)*wet_mask(i, j)
         hTr(i, j, k_top(i, j)) = hTr(i, j, k_top(i, j)) + cell
      end do
   end subroutine apply_surface_src_2d_nobudget_impl

   pure subroutine apply_surface_src_2d_dyn_nobudget_impl(hTr, wet_mask, wet_dyn, &
                                                          Q_field, inv_scale, k_top, nz, nx, ny)
      !! Wet/dry NOBUDGET twin — see `apply_surface_src_2d_nobudget_impl`
      !! and `apply_surface_src_2d_dyn_impl`.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: wet_dyn(nx, ny)
      real(wp), intent(in)    :: Q_field(nx, ny)
      real(wp), intent(in)    :: inv_scale
      integer, intent(in)    :: k_top(nx, ny)
         !! See `apply_surface_src_2d_impl`.
      integer :: i, j
      real(wp) :: cell
      do concurrent(j=1:ny, i=1:nx) local(cell)
         cell = inv_scale*Q_field(i, j)*wet_mask(i, j)*wet_dyn(i, j)
         hTr(i, j, k_top(i, j)) = hTr(i, j, k_top(i, j)) + cell
      end do
   end subroutine apply_surface_src_2d_dyn_nobudget_impl

   subroutine ocean_surface_flux_apply_sw_penetration(grid, sf, ms, dt, active, cover_frac)
      !! Additive correction that redistributes the penetrating
      !! shortwave fraction of `Q_heat` through the upper water column
      !! as a two-band exponential (Paulson & Simpson 1977; Jerlov
      !! types), instead of leaving all of it deposited at the surface
      !! layer by `ocean_surface_flux_apply_tracers`.
      !!
      !! Penetrating irradiance at downward depth `d`:
      !!   I(d) = I0 · [ R·exp(-d/zeta1) + (1-R)·exp(-d/zeta2) ]
      !! with `I0 = sw_pen_frac · Q_heat(i,j)`.  Per-layer absorbed SW =
      !! `I(d_top) - I(d_bot)`; the bed (`k = k_bot`, the first LIVE layer
      !! counting up — `1` off `z_fixed`) is treated as opaque
      !! (`I_bot := 0`) so the column absorbs all of `I0` and energy is
      !! conserved exactly (Σ_k absorbed_k = I0).
      !!
      !! Additive-correction structure: the surface kernel already
      !! deposited the full `I0` lump at `k = nz`; this kernel removes it
      !! there (`-inv_scale · I0`) and adds the distributed profile, so
      !! the net column heat change versus the legacy all-at-`nz`
      !! deposition is zero — shortwave only MOVES heat in depth.  Both
      !! the `hTr` and `heat_budget_surface` increments are mirrored, so
      !! the total surface heat budget is unchanged (just depth-spread).
      !!
      !! Bottom-up convention (`k = nz` surface, `k = 1` bed).  No-op
      !! unless `sf%has_sw .and. sf%has_heat`, a temperature tracer is
      !! registered, and (optionally) `active` is true.  Default-off
      !! (`sw_pen_frac = 0` ⇒ `has_sw = .false.`) leaves the path
      !! byte-for-byte unchanged.
      type(hgrid_t), intent(in) :: grid
      type(ocean_surface_flux_t), intent(in), optional :: sf
         !! Optional — absent ⇒ no-op (no forcing configured).
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      logical, intent(in), optional :: active
         !! Optional thermo-cadence gate.  Present-and-false ⇒ early
         !! return; absent ⇒ kernel runs.
      real(wp), intent(in), optional :: cover_frac(:, :)
         !! Optional ice-shelf cover fraction (`metrics%cover_frac`,
         !! v1 binary).  Present ⇒ the penetrating irradiance is scaled
         !! by `1 - cover_frac`, so a covered column absorbs NOTHING —
         !! no sunlight reaches the ocean through several hundred metres
         !! of ice.  Needed even on the `sw_source="net_heat"` branch,
         !! where `Q_heat` is already assembler-masked: this kernel's
         !! job is to MOVE a surface lump down the column, and on a
         !! masked column the lump it would remove was never deposited
         !! (the same argument `&ocean_wetdry_nml` makes for a dry
         !! column).  Absent ⇒ the original kernel, byte-identical.
      ! assumed-shape-ok: thermo-cadence shim, forwarded to an
      ! explicit-shape `_impl` before the device loop.

      integer :: nx, ny, nz, idx_T
      logical :: masked

      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. present(sf)) return
      if (.not. (sf%has_sw .and. sf%has_heat)) return
      if (.not. allocated(ms%tracers)) return

      idx_T = ms%idx_temperature
      if (idx_T <= 0) return

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

      ! Shim+_impl split: keep the `tracers(idx_T)%hTr` registry deref on
      ! the host (array-of-DT indirection blocks NVHPC device codegen).
      ! The irradiance source is selected HOST-SIDE (`sw_from_qsw`): the
      ! PR-12 `q_sw` component is allocated only under `use_components`,
      ! so passing it as an actual argument is only legal on the branch
      ! guarded by the host flag (validate_config forces
      ! enable_components when sw_source="q_sw", making this total).  The
      ! additive-correction identity `-I0 + Σ_k I0·(T_top - T_bot) = 0`
      ! holds for ANY I0, so the source swap cannot break conservation.
      masked = .false.
      if (present(cover_frac)) then
         masked = (size(cover_frac, 1) == nx .and. size(cover_frac, 2) == ny)
      end if
      if (masked) then
         if (sf%sw_from_qsw) then
            call apply_sw_penetration_cover_impl(ms%tracers(idx_T)%hTr, &
                                                 ms%heat_budget_surface, &
                                                 ms%h_layer, ms%k_bot, ms%wet_mask, cover_frac, sf%q_sw, &
                                                 dt/(sf%rho0*sf%cp), sf%sw_pen_frac, &
                                                 sf%sw_band_ratio, sf%sw_zeta1, sf%sw_zeta2, &
                                                 nz, nx, ny)
         else
            call apply_sw_penetration_cover_impl(ms%tracers(idx_T)%hTr, &
                                                 ms%heat_budget_surface, &
                                                 ms%h_layer, ms%k_bot, ms%wet_mask, cover_frac, sf%Q_heat, &
                                                 dt/(sf%rho0*sf%cp), sf%sw_pen_frac, &
                                                 sf%sw_band_ratio, sf%sw_zeta1, sf%sw_zeta2, &
                                                 nz, nx, ny)
         end if
         return
      end if
      if (sf%sw_from_qsw) then
         call apply_sw_penetration_impl(ms%tracers(idx_T)%hTr, &
                                        ms%heat_budget_surface, &
                                        ms%h_layer, ms%k_bot, ms%wet_mask, sf%q_sw, &
                                        dt/(sf%rho0*sf%cp), sf%sw_pen_frac, &
                                        sf%sw_band_ratio, sf%sw_zeta1, sf%sw_zeta2, &
                                        nz, nx, ny)
      else
         call apply_sw_penetration_impl(ms%tracers(idx_T)%hTr, &
                                        ms%heat_budget_surface, &
                                        ms%h_layer, ms%k_bot, ms%wet_mask, sf%Q_heat, &
                                        dt/(sf%rho0*sf%cp), sf%sw_pen_frac, &
                                        sf%sw_band_ratio, sf%sw_zeta1, sf%sw_zeta2, &
                                        nz, nx, ny)
      end if
   end subroutine ocean_surface_flux_apply_sw_penetration

   pure function sw_transmission(d, R, zeta1, zeta2) result(trans)
      !! Two-band (Paulson & Simpson 1977) normalised downward
      !! irradiance transmission at depth `d` below the free surface,
      !! `T(d) = R·exp(-d/zeta1) + (1-R)·exp(-d/zeta2)`, `T(0) = 1`.
      !! This is THE single shared definition — the SW deposition
      !! kernel, the KPP `MXL_SW`/`LV1_SW` boundary-layer correction, and
      !! the EPBL in-layer PE-cost ledger all consume it, so the three
      !! consumers cannot disagree about where the sunlight went.  Marked
      !! `!$acc routine seq` so it inlines into same-module device kernels
      !! and is callable from cross-module `do concurrent` kernels.
      !$acc routine seq
      real(wp), intent(in) :: d, R, zeta1, zeta2
      real(wp) :: trans
      trans = R*exp(-d/zeta1) + (1.0_wp - R)*exp(-d/zeta2)
   end function sw_transmission

   pure function sw_pe_cost_shape(tau) result(phi)
      !! In-layer potential-energy-cost shape function `Phi(tau)` for the
      !! EPBL TKE ledger, `tau = h/zeta` the in-layer optical depth of a
      !! single band.  It is the fraction of the pure-skin PE cost that
      !! homogenising an EXPONENTIALLY distributed in-layer heating
      !! actually incurs (Paulson & Simpson 1977 profile; the EPBL
      !! energetics of Reichl & Hallberg 2018):
      !!
      !!   Phi(tau) = [ tau·(1+e^-tau) - 2·(1-e^-tau) ] / [ tau·(1-e^-tau) ]
      !!
      !! Limits: `Phi(0) = 0` (heating already uniform through the layer
      !! ⇒ homogenising it costs nothing) and `Phi(∞) = 1` (all heating
      !! at the layer top ⇒ full skin cost, recovering Roundabout's existing
      !! `ctke_sfc` skin form exactly).  The closed form is `0/0` as
      !! `tau -> 0` and cancels catastrophically for small `tau`; a thin
      !! layer (`h << zeta2 = 23 m`) is the COMMON case, so the Taylor
      !! branch `Phi ≈ (tau/6)·(1 - tau²/60)` is mandatory below the
      !! `tau = 1e-2` seam.  `!$acc routine seq` for cross-module device
      !! calls (the EPBL prep sweep).
      !$acc routine seq
      real(wp), intent(in) :: tau
      real(wp) :: phi
      real(wp) :: em1
      real(wp), parameter :: TAU_TAYLOR = 1.0e-2_wp
      real(wp), parameter :: C1_6 = 1.0_wp/6.0_wp
      real(wp), parameter :: C1_60 = 1.0_wp/60.0_wp
      if (tau <= TAU_TAYLOR) then
         phi = C1_6*tau*(1.0_wp - C1_60*tau*tau)
      else
         em1 = 1.0_wp - exp(-tau)
         phi = (tau*(1.0_wp + exp(-tau)) - 2.0_wp*em1)/(tau*em1)
      end if
   end function sw_pe_cost_shape

   pure subroutine apply_sw_penetration_impl(hTr, budget, h_layer, k_bot, wet_mask, &
                                             sw_src, inv_scale, sw_pen_frac, R, &
                                             zeta1, zeta2, nz, nx, ny)
      !! Per-column two-band shortwave redistribution.  Explicit-shape
      !! dummies so NVHPC stdpar compiles device kernels against static
      !! bounds.  Difference form (`I(d_top) - I(d_bot)`) — no division
      !! by `h`, so a vanishing layer (`h → 0 ⇒ d_top == d_bot`) absorbs
      !! zero automatically with no guard.  The transmission `T(d)` is
      !! the shared `sw_transmission` (same-module ⇒ inlined by NVHPC),
      !! so the deposition and the boundary-layer coupling cannot diverge.
      !!
      !! `inv_scale` = dt/(rho0·cp); `R`/`zeta1`/`zeta2` are the two-band
      !! parameters; `sw_pen_frac` scales `sw_src` to the penetrating
      !! irradiance `I0`.  `sw_src` is the caller-selected source
      !! (`Q_heat` for the legacy net-heat path, `q_sw` for the PR-12
      !! component path); the additive-correction conservation identity
      !! `-I0 + Σ_k I0·(T_top - T_bot) = 0` holds for ANY `sw_src`, so the
      !! source swap cannot break conservation.  At `k = nz` the legacy
      !! `I0` lump is subtracted before adding the surface band, keeping
      !! the net column change against the all-at-`nz` baseline at zero.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      integer, intent(in)     :: k_bot(nx, ny)
         !! `ms%k_bot` — the first LIVE layer counting up from the bed.  The
         !! opaque-bed residual lands HERE, not on `k = 1`: under `z_fixed`
         !! the layers below are inert fillers (`1` elsewhere ⇒ unchanged).
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: sw_src(nx, ny)
      real(wp), intent(in)    :: inv_scale, sw_pen_frac, R, zeta1, zeta2
      integer :: i, j, k
      real(wp) :: i0col, d_top, d_bot, trans_top, trans_bot, absorbed, add

      do concurrent(j=1:ny, i=1:nx) local(i0col, d_top, d_bot, trans_top, &
                                          trans_bot, absorbed, add, k)
         i0col = sw_pen_frac*sw_src(i, j)*wet_mask(i, j)
         d_top = 0.0_wp
         do k = nz, k_bot(i, j), -1
            d_bot = d_top + h_layer(i, j, k)
            trans_top = sw_transmission(d_top, R, zeta1, zeta2)
            if (k > k_bot(i, j)) then
               trans_bot = sw_transmission(d_bot, R, zeta1, zeta2)
            else
               trans_bot = 0.0_wp   ! bed opaque: column absorbs all of I0
            end if
            absorbed = i0col*(trans_top - trans_bot)
            add = inv_scale*absorbed
            if (k == nz) add = add - inv_scale*i0col   ! remove the surface lump
            hTr(i, j, k) = hTr(i, j, k) + add
            budget(i, j, k) = budget(i, j, k) + add
            d_top = d_bot
         end do
      end do
   end subroutine apply_sw_penetration_impl

   pure subroutine apply_sw_penetration_cover_impl(hTr, budget, h_layer, k_bot, wet_mask, &
                                                   cover_frac, sw_src, inv_scale, &
                                                   sw_pen_frac, R, zeta1, zeta2, nz, nx, ny)
      !! Ice-shelf-cover twin of `apply_sw_penetration_impl`: the
      !! open-water factor `1 - cover_frac` composes multiplicatively
      !! with `wet_mask` into the column irradiance `I0`, so a fully
      !! covered column neither removes the surface lump nor deposits a
      !! profile — it is left EXACTLY untouched.  Separate `_impl`, not
      !! an in-loop `present()` test (house idiom, see
      !! `apply_surface_src_2d_dyn_impl`) — the cover-off path keeps the
      !! original kernel byte-identical.
      !!
      !! The additive-correction conservation identity
      !! `-I0 + Σ_k I0·(T_top - T_bot) = 0` holds for ANY `I0`, and
      !! `I0 = 0` is the degenerate case of it, so scaling the source
      !! cannot break column heat conservation.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      integer, intent(in)     :: k_bot(nx, ny)
         !! `ms%k_bot` — the first LIVE layer counting up from the bed.  The
         !! opaque-bed residual lands HERE, not on `k = 1`: under `z_fixed`
         !! the layers below are inert fillers (`1` elsewhere ⇒ unchanged).
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: cover_frac(nx, ny)
      real(wp), intent(in)    :: sw_src(nx, ny)
      real(wp), intent(in)    :: inv_scale, sw_pen_frac, R, zeta1, zeta2
      integer :: i, j, k
      real(wp) :: i0col, d_top, d_bot, trans_top, trans_bot, absorbed, add

      do concurrent(j=1:ny, i=1:nx) local(i0col, d_top, d_bot, trans_top, &
                                          trans_bot, absorbed, add, k)
         i0col = sw_pen_frac*sw_src(i, j)*wet_mask(i, j)*(1.0_wp - cover_frac(i, j))
         d_top = 0.0_wp
         do k = nz, k_bot(i, j), -1
            d_bot = d_top + h_layer(i, j, k)
            trans_top = sw_transmission(d_top, R, zeta1, zeta2)
            if (k > k_bot(i, j)) then
               trans_bot = sw_transmission(d_bot, R, zeta1, zeta2)
            else
               trans_bot = 0.0_wp   ! bed opaque: column absorbs all of I0
            end if
            absorbed = i0col*(trans_top - trans_bot)
            add = inv_scale*absorbed
            if (k == nz) add = add - inv_scale*i0col   ! remove the surface lump
            hTr(i, j, k) = hTr(i, j, k) + add
            budget(i, j, k) = budget(i, j, k) + add
            d_top = d_bot
         end do
      end do
   end subroutine apply_sw_penetration_cover_impl

   subroutine ocean_surface_restore_apply_tracers(grid, sf, ms, dt, active, cover_frac)
      !! Surface buoyancy restoring (MOM6 `RESTOREBUOY`): relax the
      !! top-layer (`k = nz`) temperature / salinity toward scalar
      !! targets with a piston velocity `p` [m/s].  Unlike
      !! `ocean_surface_flux_apply_tracers` (which reads a pre-filled
      !! static `Q_*` field), the restoring flux is DYNAMIC — it depends
      !! on the live SST / SSS each thermo step — so it is computed
      !! in-kernel from `(target - surface_concentration)` rather than a
      !! stored field.  This keeps the const-flux path byte-for-byte
      !! unchanged and adds no second device sync of `Q_*`.
      !!
      !! Bottom-up convention: surface = `k = nz`, bed = `k = 1`.  The
      !! surface concentration is `hTr(i,j,nz) / max(h_layer(i,j,nz),
      !! h_min)`.  Per thermo step, in hTr-space:
      !!   d(hT_top) = dt · p_T · (T_target - SST) · wet_mask   [K·m]
      !!   d(hS_top) = dt · p_S · (S_target - SSS) · wet_mask   [PSU·m]
      !! (The `rho0·cp` of the equivalent W/m^2 restoring heat flux
      !! cancels against the `dt/(rho0·cp)` apply scaling — see the spec.)
      !! The same increment is mirrored into `heat_budget_surface(:,:,nz)`
      !! / `salt_budget_surface(:,:,nz)` so the surface-budget diagnostics
      !! see restoring as an explicit (non-conservative) source.
      !!
      !! Restoring is a relaxation forcing, NOT a conservative process:
      !! it deliberately injects / removes heat + salt to nudge the
      !! surface.  The budget contributors account for it so the
      !! conservation diagnostics do not flag it as a leak.
      !!
      !! MOM6-fidelity note (salt path): the temperature branch is an
      !! exact analogue of MOM6 `RESTOREBUOY` `heat_added` (the `rho0*cp`
      !! cancellation above makes the tendency identical for a given
      !! `FLUXCONST_T`).  MOM6's SALT branch is instead a virtual
      !! freshwater flux `vprec = -rho0*p_S*(S*-SSS)/(0.5*(SSS+S*))` that
      !! changes the column mass and dilutes a CONSERVED salt content.
      !! We use a linearised salt-CONTENT injection (no thickness change),
      !! which reproduces MOM6's surface-salinity tendency TO FIRST ORDER
      !! — the dilution factor `SSS/(0.5*(SSS+S*)) ~ 1` for realistic
      !! anomalies — but is non-conservative and drops that second-order
      !! factor.  Exact `vprec` parity is a documented v2 follow-up.
      !!
      !! No-op unless `sf%has_restore_T .or. sf%has_restore_S` (set by
      !! `set_restore` when the matching switch is on AND the piston is
      !! non-zero), a matching tracer is registered, and (optionally)
      !! `active` is true.  Default-off leaves the surface-flux path
      !! byte-for-byte unchanged.
      type(hgrid_t), intent(in) :: grid
      type(ocean_surface_flux_t), intent(in), optional :: sf
         !! Optional — absent ⇒ no-op (no forcing configured).
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      logical, intent(in), optional :: active
         !! Optional thermo-cadence gate.  Present-and-false ⇒ early
         !! return; absent ⇒ kernel runs.
      real(wp), intent(in), optional :: cover_frac(:, :)
         !! Optional ice-shelf cover fraction (`metrics%cover_frac`,
         !! v1 binary).  Present ⇒ the restoring increment is scaled by
         !! `1 - cover_frac`, so a covered column is NOT relaxed toward
         !! an atmospheric target — under a shelf the surface is a
         !! melting ice interface, and restoring there would overwhelm
         !! the melt signal with a number the atmosphere never set.
         !! Absent ⇒ the original kernel, byte-identical.
      ! assumed-shape-ok: thermo-cadence shim, forwarded to an
      ! explicit-shape `_impl` before the device loop.

      integer :: nx, ny, nz, idx_T, idx_S
      logical :: masked

      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. present(sf)) return
      if (.not. sf%has_restore_T .and. .not. sf%has_restore_S) return
      if (.not. allocated(ms%tracers)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      idx_T = ms%idx_temperature
      idx_S = ms%idx_salinity

      masked = .false.
      if (present(cover_frac)) then
         masked = (size(cover_frac, 1) == nx .and. size(cover_frac, 2) == ny)
      end if
      if (masked) then
         if (idx_T > 0 .and. sf%has_restore_T) then
            call apply_surface_restore_2d_cover_impl(ms%tracers(idx_T)%hTr, &
                                                     ms%heat_budget_surface, &
                                                     ms%h_layer, ms%wet_mask, cover_frac, &
                                                     dt*sf%restore_piston_T, &
                                                     sf%restore_T_target, sf%h_min, &
                                                     nz, nx, ny)
         end if
         if (idx_S > 0 .and. sf%has_restore_S) then
            call apply_surface_restore_2d_cover_impl(ms%tracers(idx_S)%hTr, &
                                                     ms%salt_budget_surface, &
                                                     ms%h_layer, ms%wet_mask, cover_frac, &
                                                     dt*sf%restore_piston_S, &
                                                     sf%restore_S_target, sf%h_min, &
                                                     nz, nx, ny)
         end if
         return
      end if

      ! Shim+_impl split: keep the `tracers(idx)%hTr` registry deref on
      ! the host (array-of-DT indirection blocks NVHPC device codegen).
      ! Scalar targets / piston / h_min pass by value into the kernel —
      ! no per-cell field, so no extra device array.
      if (idx_T > 0 .and. sf%has_restore_T) then
         call apply_surface_restore_2d_impl(ms%tracers(idx_T)%hTr, &
                                            ms%heat_budget_surface, &
                                            ms%h_layer, ms%wet_mask, &
                                            dt*sf%restore_piston_T, &
                                            sf%restore_T_target, sf%h_min, &
                                            nz, nx, ny)
      end if
      if (idx_S > 0 .and. sf%has_restore_S) then
         call apply_surface_restore_2d_impl(ms%tracers(idx_S)%hTr, &
                                            ms%salt_budget_surface, &
                                            ms%h_layer, ms%wet_mask, &
                                            dt*sf%restore_piston_S, &
                                            sf%restore_S_target, sf%h_min, &
                                            nz, nx, ny)
      end if
   end subroutine ocean_surface_restore_apply_tracers

   pure subroutine apply_surface_restore_2d_impl(hTr, budget, h_layer, &
                                                 wet_mask, dt_piston, tgt, &
                                                 h_min, nz, nx, ny)
      !! Stamp the per-step restoring increment `dt·p·(target - surf)·
      !! wet_mask` onto the top layer (k = nz) of a tracer's hTr array
      !! and mirror it into the matching budget contributor.
      !! Explicit-shape dummies so NVHPC stdpar compiles device kernels
      !! against static bounds.
      !!
      !! `dt_piston` = dt·piston [m] is the column-invariant multiplier
      !! the caller derives from `sf` (the `rho0·cp` cancels — see the
      !! caller doc).  The surface concentration is recovered as
      !! `hTr(nz) / max(h_layer(nz), h_min)`; the `h_min` floor keeps the
      !! relaxation finite under a pinched top layer (rate `p/h`
      !! saturates rather than diverges).
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: dt_piston, tgt, h_min
      integer :: i, j
      real(wp) :: surf, inc
      do concurrent(j=1:ny, i=1:nx) local(surf, inc)
         surf = hTr(i, j, nz)/max(h_layer(i, j, nz), h_min)
         inc = dt_piston*(tgt - surf)*wet_mask(i, j)
         hTr(i, j, nz) = hTr(i, j, nz) + inc
         budget(i, j, nz) = budget(i, j, nz) + inc
      end do
   end subroutine apply_surface_restore_2d_impl

   pure subroutine apply_surface_restore_2d_cover_impl(hTr, budget, h_layer, &
                                                       wet_mask, cover_frac, dt_piston, &
                                                       tgt, h_min, nz, nx, ny)
      !! Ice-shelf-cover twin of `apply_surface_restore_2d_impl`: the
      !! open-water factor `1 - cover_frac` composes multiplicatively
      !! with `wet_mask`, so a covered column receives EXACTLY zero
      !! restoring — and, because the same factor multiplies the budget
      !! mirror, exactly zero restoring shows up in the heat/salt
      !! surface budget there too.  Separate `_impl`, not an in-loop
      !! `present()` test (house idiom).
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: budget(nx, ny, nz)
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      real(wp), intent(in)    :: wet_mask(nx, ny)
      real(wp), intent(in)    :: cover_frac(nx, ny)
      real(wp), intent(in)    :: dt_piston, tgt, h_min
      integer :: i, j
      real(wp) :: surf, inc
      do concurrent(j=1:ny, i=1:nx) local(surf, inc)
         ! vanished-ok: surface restoring measures the surface concentration against a
         ! CALLER-supplied floor `h_min` (a restoring depth), not the
         ! vanish marker; zero would restore toward a fictitious fresh sea.
         surf = hTr(i, j, nz)/max(h_layer(i, j, nz), h_min)
         inc = dt_piston*(tgt - surf)*wet_mask(i, j)*(1.0_wp - cover_frac(i, j))
         hTr(i, j, nz) = hTr(i, j, nz) + inc
         budget(i, j, nz) = budget(i, j, nz) + inc
      end do
   end subroutine apply_surface_restore_2d_cover_impl

   pure function ocean_surface_flux_bytes(this) result(nbytes)
      !! Counted allocatable footprint of the surface flux slot
      !! (0 when unallocated). One arr_bytes term per array — add a
      !! term here when a new allocatable joins the type.
      class(ocean_surface_flux_t), intent(in) :: this
      integer(int64) :: nbytes
      nbytes = arr_bytes(this%Q_heat) &
               + arr_bytes(this%Q_salt) &
               + arr_bytes(this%q_sw) + arr_bytes(this%q_lw) &
               + arr_bytes(this%q_lat) + arr_bytes(this%q_sens) &
               + arr_bytes(this%heat_added) &
               + arr_bytes(this%heat_cavity) &
               + arr_bytes(this%evap) + arr_bytes(this%lprec) &
               + arr_bytes(this%fprec) + arr_bytes(this%vprec) &
               + arr_bytes(this%lrunoff) + arr_bytes(this%frunoff) &
               + arr_bytes(this%seaice_melt) &
               + arr_bytes(this%heat_content_lprec) &
               + arr_bytes(this%heat_content_fprec) &
               + arr_bytes(this%heat_content_vprec) &
               + arr_bytes(this%heat_content_lrunoff) &
               + arr_bytes(this%heat_content_frunoff) &
               + arr_bytes(this%heat_content_seaice_melt) &
               + arr_bytes(this%heat_content_massin) &
               + arr_bytes(this%heat_content_massout) &
               + arr_bytes(this%salt_flux) &
               + arr_bytes(this%salt_cavity) &
               + arr_bytes(this%p_surf_atm) + arr_bytes(this%p_surf)
   end function ocean_surface_flux_bytes

   pure subroutine ocean_surface_flux_apply_cover_const(sf, cover_frac)
      !! Mask the STATIC scalar `Q_heat` / `Q_salt` fill with the
      !! ice-shelf cover, for the `use_components = .false.` path only.
      !!
      !! With the component set on, `Q_heat`/`Q_salt` are assembler
      !! outputs and the mask belongs there (`ocean_surface_flux_assemble`'s
      !! `cover_frac`); writing them here would be a second writer on an
      !! assembler-owned slot, which the fill contract forbids.  With the
      !! component set OFF there is no assembler, `Q_heat`/`Q_salt` are
      !! the configure-time `set_surface_flux_const` fill and nothing
      !! rewrites them per step — so masking them once, at configure
      !! after the cover is built, is the whole job.  Without this, a
      !! geometry-only cavity run (`&ocean_cavity_dyn_nml` with no melt)
      !! would still push a uniform `&ocean_thermo_nml q_heat` through
      !! the ice.
      !!
      !! No-op when the component set is on, when the fields are
      !! unallocated, or when `cover_frac` is the `(1,1)` placeholder.
      !! Idempotent (multiplies by 0 or 1).  Host-side at configure —
      !! call BEFORE `enter_data` (or follow with an `!$acc update
      !! device`).
      type(ocean_surface_flux_t), intent(inout) :: sf
      real(wp), intent(in) :: cover_frac(:, :)
         !! Ice-cover fraction at cell centres (`metrics%cover_frac`).
      ! assumed-shape-ok: configure-time, one call per run.
      integer :: nx, ny

      if (sf%use_components) return
      if (.not. allocated(sf%Q_heat) .or. .not. allocated(sf%Q_salt)) return
      nx = size(sf%Q_heat, 1)
      ny = size(sf%Q_heat, 2)
      if (size(cover_frac, 1) /= nx .or. size(cover_frac, 2) /= ny) return
      call ocean_surfflux_cover_const_impl(sf%Q_heat, sf%Q_salt, cover_frac, nx, ny)
   end subroutine ocean_surface_flux_apply_cover_const

   pure subroutine ocean_surfflux_cover_const_impl(Q_heat, Q_salt, cover_frac, nx, ny)
      !! Flat `do concurrent` kernel behind
      !! `ocean_surface_flux_apply_cover_const` — explicit-shape dummies,
      !! integer dims first (decl-order, ifx #8586).
      integer, intent(in)    :: nx, ny
      real(wp), intent(inout) :: Q_heat(nx, ny), Q_salt(nx, ny)
      real(wp), intent(in)    :: cover_frac(nx, ny)
      integer :: i, j
      real(wp) :: open_f
      do concurrent(j=1:ny, i=1:nx) local(open_f)
         open_f = 1.0_wp - cover_frac(i, j)
         Q_heat(i, j) = Q_heat(i, j)*open_f
         Q_salt(i, j) = Q_salt(i, j)*open_f
      end do
   end subroutine ocean_surfflux_cover_const_impl

   pure subroutine ocean_surface_flux_assemble(grid, sf, ms, active, cover_frac)
      !! **The single gate** that derives `Q_heat`/`Q_salt` from the
      !! component set (§3.1/§3.3 of the PR-12 plan) — the exact analogue
      !! of `vmix_assemble`: fillers contribute components, this routine
      !! alone derives the net fields every downstream kernel reads.  A
      !! no-op unless `sf%use_components` — with components off, `Q_heat`
      !! / `Q_salt` are exactly what `set_surface_flux_const` (or a
      !! field-override path) left them, byte-for-byte.
      !!
      !! Net surface heat into the ocean:
      !!   Q_heat = Q_heat_const + q_sw + q_lw + q_lat + q_sens + heat_added
      !!          + heat_cavity + heat_content_massin + heat_content_massout
      !!   heat_content_massin  = Σ heat_content_{lprec,fprec,vprec,
      !!                            lrunoff,frunoff,seaice_melt}
      !!   heat_content_massout = SEAWATER_CP * T_sst * evap
      !! Net surface salt flux:
      !!   Q_salt = Q_salt_const + salt_flux + salt_cavity
      !! All four outputs multiplied by `ms%wet_mask` (land carries
      !! exactly zero; interior loop bounds are NOT restricted — see
      !! CLAUDE.md's "nghost and the assembler loop bounds" gotcha).
      !!
      !! `has_heat`/`has_salt` are deliberately NOT touched here — a
      !! components-on run with no live filler must reproduce EXACTLY
      !! what `set_surface_flux_const` left them (the "+0.0 bit-identity
      !! trap": forcing them true would change `apply_surface_src_2d_impl`
      !! from an early-return to a `+0.0` stamp, altering a budget sum's
      !! operand COUNT even though the value is unchanged).  Setting
      !! `has_heat`/`has_salt` is the FILLER's job (see the module
      !! docstring's fill contract).
      !!
      !! Outer-shim + flat-impl (`ocean_surfflux_assemble_impl`): the SST
      !! read needs `ms%tracers(idx_temperature)%hTr`, an array-of-DT
      !! registry deref that must happen on the host before the `do
      !! concurrent` (NVHPC device codegen constraint — see
      !! `ocean_surface_flux_apply_tracers` for the identical pattern).
      !! No-op when no temperature tracer is registered (SST is
      !! undefined without one).
      type(hgrid_t), intent(in) :: grid
      type(ocean_surface_flux_t), intent(inout) :: sf
      type(multilayer_state_t), intent(in) :: ms
      logical, intent(in), optional :: active
         !! Optional thermo-cadence gate.  Present-and-false ⇒ early
         !! return; absent ⇒ kernel runs (matches
         !! `ocean_surface_flux_apply_tracers`'s convention).
      real(wp), intent(in), optional :: cover_frac(:, :)
         !! Optional ice-shelf cover fraction (`metrics%cover_frac`,
         !! v1 binary).  Present ⇒ every ATMOSPHERIC contribution
         !! (`Q_heat_const`, `q_sw`, `q_lw`, `q_lat`, `q_sens`,
         !! `heat_added`, both mass-enthalpy terms, `Q_salt_const`,
         !! `salt_flux`) is scaled by `1 - cover_frac`, while the
         !! cavity's OWN `heat_cavity` / `salt_cavity` pass through
         !! unmasked.  This is the single place those two groups are
         !! still distinguishable — see the module docstring for why the
         !! mask lives here and not at apply time.  Absent ⇒ the
         !! original kernel, byte-identical.
      ! assumed-shape-ok: thermo-cadence shim, forwarded to an
      ! explicit-shape `_impl` before the device loop.

      integer :: nx, ny, nz, idx_T
      logical :: masked

      if (.not. sf%use_components) return
      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. allocated(ms%tracers)) return
      idx_T = ms%idx_temperature
      if (idx_T <= 0) return

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

      masked = .false.
      if (present(cover_frac)) then
         masked = (size(cover_frac, 1) == nx .and. size(cover_frac, 2) == ny)
      end if
      if (masked) then
         call ocean_surfflux_assemble_cover_impl( &
            sf%heat_content_massin, sf%heat_content_massout, sf%Q_heat, sf%Q_salt, &
            sf%q_sw, sf%q_lw, sf%q_lat, sf%q_sens, sf%heat_added, sf%heat_cavity, &
            sf%heat_content_lprec, sf%heat_content_fprec, sf%heat_content_vprec, &
            sf%heat_content_lrunoff, sf%heat_content_frunoff, sf%heat_content_seaice_melt, &
            sf%evap, sf%salt_flux, sf%salt_cavity, &
            ms%tracers(idx_T)%hTr, ms%h_layer, ms%wet_mask, cover_frac, &
            sf%Q_heat_const, sf%Q_salt_const, sf%cp, sf%h_min, nz, nx, ny)
         sf%q_assembled = 1.0_wp
         return
      end if

      call ocean_surfflux_assemble_impl( &
         sf%heat_content_massin, sf%heat_content_massout, sf%Q_heat, sf%Q_salt, &
         sf%q_sw, sf%q_lw, sf%q_lat, sf%q_sens, sf%heat_added, sf%heat_cavity, &
         sf%heat_content_lprec, sf%heat_content_fprec, sf%heat_content_vprec, &
         sf%heat_content_lrunoff, sf%heat_content_frunoff, sf%heat_content_seaice_melt, &
         sf%evap, sf%salt_flux, sf%salt_cavity, &
         ms%tracers(idx_T)%hTr, ms%h_layer, ms%wet_mask, &
         sf%Q_heat_const, sf%Q_salt_const, sf%cp, sf%h_min, nz, nx, ny)
      sf%q_assembled = 1.0_wp
   end subroutine ocean_surface_flux_assemble

   pure subroutine ocean_surfflux_assemble_impl(heat_content_massin, heat_content_massout, &
                                                Q_heat, Q_salt, &
                                                q_sw, q_lw, q_lat, q_sens, heat_added, &
                                                heat_cavity, &
                                                heat_content_lprec, heat_content_fprec, &
                                                heat_content_vprec, heat_content_lrunoff, &
                                                heat_content_frunoff, heat_content_seaice_melt, &
                                                evap, salt_flux, salt_cavity, &
                                                hTr_T, h_layer, wet_mask, &
                                                Q_heat_const, Q_salt_const, cp, h_min, &
                                                nz, nx, ny)
      !! Flat `do concurrent` kernel — explicit-shape dummies, integer
      !! dims declared first (decl-order, ifx #8586).  See
      !! `ocean_surface_flux_assemble` for the physics; this is the
      !! arithmetic verbatim.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: heat_content_massin(nx, ny), heat_content_massout(nx, ny)
      real(wp), intent(inout) :: Q_heat(nx, ny), Q_salt(nx, ny)
      real(wp), intent(in)    :: q_sw(nx, ny), q_lw(nx, ny), q_lat(nx, ny), q_sens(nx, ny)
      real(wp), intent(in)    :: heat_added(nx, ny), heat_cavity(nx, ny)
      real(wp), intent(in)    :: heat_content_lprec(nx, ny), heat_content_fprec(nx, ny)
      real(wp), intent(in)    :: heat_content_vprec(nx, ny), heat_content_lrunoff(nx, ny)
      real(wp), intent(in)    :: heat_content_frunoff(nx, ny), heat_content_seaice_melt(nx, ny)
      real(wp), intent(in)    :: evap(nx, ny), salt_flux(nx, ny), salt_cavity(nx, ny)
      real(wp), intent(in)    :: hTr_T(nx, ny, nz), h_layer(nx, ny, nz), wet_mask(nx, ny)
      real(wp), intent(in)    :: Q_heat_const, Q_salt_const, cp, h_min
      integer :: i, j
      real(wp) :: sst, massin, massout

      do concurrent(j=1:ny, i=1:nx) local(sst, massin, massout)
         massin = heat_content_lprec(i, j) + heat_content_fprec(i, j) &
                  + heat_content_vprec(i, j) + heat_content_lrunoff(i, j) &
                  + heat_content_frunoff(i, j) + heat_content_seaice_melt(i, j)
         sst = hTr_T(i, j, nz)/max(h_layer(i, j, nz), h_min)
         massout = cp*sst*evap(i, j)
         heat_content_massin(i, j) = wet_mask(i, j)*massin
         heat_content_massout(i, j) = wet_mask(i, j)*massout
         Q_heat(i, j) = wet_mask(i, j)* &
                        (Q_heat_const + q_sw(i, j) + q_lw(i, j) + q_lat(i, j) + &
                         q_sens(i, j) + heat_added(i, j) + heat_cavity(i, j) + &
                         massin + massout)
         Q_salt(i, j) = wet_mask(i, j)*(Q_salt_const + salt_flux(i, j) + &
                                        salt_cavity(i, j))
      end do
   end subroutine ocean_surfflux_assemble_impl

   pure subroutine ocean_surfflux_assemble_cover_impl(heat_content_massin, &
                                                      heat_content_massout, &
                                                      Q_heat, Q_salt, &
                                                      q_sw, q_lw, q_lat, q_sens, heat_added, &
                                                      heat_cavity, &
                                                      heat_content_lprec, heat_content_fprec, &
                                                      heat_content_vprec, heat_content_lrunoff, &
                                                      heat_content_frunoff, heat_content_seaice_melt, &
                                                      evap, salt_flux, salt_cavity, &
                                                      hTr_T, h_layer, wet_mask, cover_frac, &
                                                      Q_heat_const, Q_salt_const, cp, h_min, &
                                                      nz, nx, ny)
      !! Ice-shelf-cover twin of `ocean_surfflux_assemble_impl`: the
      !! open-water factor `open_f = 1 - cover_frac` multiplies the
      !! ATMOSPHERIC group and NOT the cavity group.  Grouping, spelt
      !! out because it is the whole point of this kernel:
      !!
      !!   masked   — `Q_heat_const`, `q_sw`, `q_lw`, `q_lat`, `q_sens`,
      !!              `heat_added`, `heat_content_massin` (the six
      !!              mass-enthalpy companions), `heat_content_massout`
      !!              (`cp·SST·evap`), `Q_salt_const`, `salt_flux`;
      !!   UNmasked — `heat_cavity`, `salt_cavity`.
      !!
      !! The two `heat_content_mass*` OUTPUTS carry the factor too: they
      !! are diagnostics of atmospheric mass exchange, which under a
      !! shelf is zero, and reporting the unmasked value next to a
      !! masked `Q_heat` would make the ledger not add up.
      !!
      !! Separate `_impl`, not an in-loop `present()` test (house
      !! idiom) — the cover-off path keeps the production assembler
      !! byte-identical.
      integer, intent(in)    :: nz, nx, ny
      real(wp), intent(inout) :: heat_content_massin(nx, ny), heat_content_massout(nx, ny)
      real(wp), intent(inout) :: Q_heat(nx, ny), Q_salt(nx, ny)
      real(wp), intent(in)    :: q_sw(nx, ny), q_lw(nx, ny), q_lat(nx, ny), q_sens(nx, ny)
      real(wp), intent(in)    :: heat_added(nx, ny), heat_cavity(nx, ny)
      real(wp), intent(in)    :: heat_content_lprec(nx, ny), heat_content_fprec(nx, ny)
      real(wp), intent(in)    :: heat_content_vprec(nx, ny), heat_content_lrunoff(nx, ny)
      real(wp), intent(in)    :: heat_content_frunoff(nx, ny), heat_content_seaice_melt(nx, ny)
      real(wp), intent(in)    :: evap(nx, ny), salt_flux(nx, ny), salt_cavity(nx, ny)
      real(wp), intent(in)    :: hTr_T(nx, ny, nz), h_layer(nx, ny, nz), wet_mask(nx, ny)
      real(wp), intent(in)    :: cover_frac(nx, ny)
      real(wp), intent(in)    :: Q_heat_const, Q_salt_const, cp, h_min
      integer :: i, j
      real(wp) :: sst, massin, massout, open_f

      do concurrent(j=1:ny, i=1:nx) local(sst, massin, massout, open_f)
         open_f = 1.0_wp - cover_frac(i, j)
         massin = open_f*(heat_content_lprec(i, j) + heat_content_fprec(i, j) &
                          + heat_content_vprec(i, j) + heat_content_lrunoff(i, j) &
                          + heat_content_frunoff(i, j) + heat_content_seaice_melt(i, j))
         ! vanished-ok: same caller-supplied `h_min` floor — the SST carrying the
         ! evaporative heat content must stay a temperature.
         sst = hTr_T(i, j, nz)/max(h_layer(i, j, nz), h_min)
         massout = open_f*cp*sst*evap(i, j)
         heat_content_massin(i, j) = wet_mask(i, j)*massin
         heat_content_massout(i, j) = wet_mask(i, j)*massout
         Q_heat(i, j) = wet_mask(i, j)* &
                        (open_f*(Q_heat_const + q_sw(i, j) + q_lw(i, j) + q_lat(i, j) + &
                                 q_sens(i, j) + heat_added(i, j)) + &
                         heat_cavity(i, j) + massin + massout)
         Q_salt(i, j) = wet_mask(i, j)*(open_f*(Q_salt_const + salt_flux(i, j)) + &
                                        salt_cavity(i, j))
      end do
   end subroutine ocean_surfflux_assemble_cover_impl

end module rdb_ocean_surface_flux