ocean_surface_restore_apply_tracers Subroutine

public 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.

Arguments

Type IntentOptional Attributes Name
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(kind=wp), intent(in) :: dt
logical, intent(in), optional :: active

Optional thermo-cadence gate. Present-and-false ⇒ early return; absent ⇒ kernel runs.
real(kind=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.

Calls

proc~~ocean_surface_restore_apply_tracers~~CallsGraph proc~ocean_surface_restore_apply_tracers ocean_surface_restore_apply_tracers proc~apply_surface_restore_2d_cover_impl apply_surface_restore_2d_cover_impl proc~ocean_surface_restore_apply_tracers->proc~apply_surface_restore_2d_cover_impl proc~apply_surface_restore_2d_impl apply_surface_restore_2d_impl proc~ocean_surface_restore_apply_tracers->proc~apply_surface_restore_2d_impl local local proc~apply_surface_restore_2d_cover_impl->local proc~apply_surface_restore_2d_impl->local

Called by

proc~~ocean_surface_restore_apply_tracers~~CalledByGraph proc~ocean_surface_restore_apply_tracers ocean_surface_restore_apply_tracers proc~apply_sw_and_restore apply_sw_and_restore proc~apply_sw_and_restore->proc~ocean_surface_restore_apply_tracers proc~run_stage run_stage proc~run_stage->proc~apply_sw_and_restore proc~run_stage_split run_stage_split proc~run_stage_split->proc~apply_sw_and_restore proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step

Variables

Type Visibility Attributes Name Initial
integer, private :: idx_S
integer, private :: idx_T
logical, private :: masked
integer, private :: nx
integer, private :: ny
integer, private :: nz

Source Code

   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