ocean_surface_flux_apply_tracers Subroutine

   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