Add surface tracer restoring

components/omega/configs/Default.yml

@@ -31,6 +31,9 @@ Omega:
     VerticalTracerFluxOrder: 3
   WindStress:
     InterpType: Isotropic
+  SurfaceRestoring:
+    TracersToRestore: [Temperature, Salinity]
+    PistonVelocity: 1.585e-5
   VertCoord:
     MovementWeightType: Uniform
   PressureGrad:
@@ -46,6 +49,7 @@ Omega:
     ViscDel4: 1.2e11
     DivFactor: 1.0
     WindForcingTendencyEnable: false
+    SurfaceTracerRestoringEnable: false
     BottomDragTendencyEnable: false
     BottomDragCoeff: 0.0
     TracerHorzAdvTendencyEnable: true

components/omega/doc/devGuide/AuxiliaryVariables.md

@@ -104,3 +104,9 @@ The following auxiliary variable groups are currently implemented:
 | WindForcingAuxVars |  ZonalStressCell ||
 ||  MeridStressCell ||
 ||  NormalStressEdge ||
+| SurfTracerRestAuxVars | TracersMonthlySurfClimoCell ||
+
+## See Also
+
+Additional information on forcing (currently wind forcing and surface tracer
+restoring) is detailed in [](omega-dev-forcing).

components/omega/doc/devGuide/Forcing.md

@@ -0,0 +1,79 @@
+(omega-dev-forcing)=
+
+# Forcing
+
+This page describes design and implementation details for forcing-related
+pathways in Omega, currently this includes:
+
+- Wind forcing
+- Surface tracer restoring
+
+## Wind forcing design
+
+### Wind forcing data flow
+
+1. External fields provide:
+   - `WindStressZonal`
+   - `WindStressMeridional`
+2. Auxiliary-state compute builds `NormalStressEdge` from those fields.
+3. Tendency term applies wind-stress forcing to edge-normal velocity tendency.
+
+### Wind forcing key classes/components
+
+- `WindForcingAuxVars`
+  - Stores wind-stress cell fields and computed `NormalStressEdge`
+  - Applies configured interpolation choice (`InterpType`)
+- `AuxiliaryState::computeMomAux`
+  - Calls `WindForcingAuxVars::computeVarsOnEdge`
+- `WindForcingOnEdge` tendency term
+  - Adds contribution proportional to normal stress and inverse layer
+    thickness in the surface layer
+
+### Wind forcing config coupling
+
+- `Omega.WindStress.InterpType`
+  - mapped to `InterpCellToEdgeOption`
+- `Omega.Tendencies.WindForcingTendencyEnable`
+  - gates execution of wind forcing tendency kernel
+
+## Surface tracer restoring design
+
+### Surface tracer restoring data flow
+
+1. External fields provide target values: `TracersMonthlySurfClimoCell` (values and units should match the state variables)
+2. Auxiliary-state compute forms restoring differences: `SurfTracerRestoringDiffsCell = target - tracer_surface`
+3. Tendency term applies restoring only at surface layer and only for tracers selected from `SurfaceRestoring.TracersToRestore`.
+
+### Surface tracer restoring key classes/components
+
+- `SurfTracerRestAuxVars`
+  - Inputs: `TracersMonthlySurfClimoCell`, tracer state array
+  - Output: `SurfTracerRestoringDiffsCell`
+  - Uses `MinLayerCell` to select surface layer index
+- `SurfaceTracerRestoringOnCell` tendency term
+  - Applies `PistonVelocity * SurfTracerRestoringDiffsCell` at surface
+- `Tendencies`
+  - Parses `SurfaceRestoring.TracersToRestore` and resolves tracer indices
+  - Builds `TracerIdsToRestore` and `NTracersToRestore`
+  - Applies tracer-selection logic at call site in
+    `computeTracerTendenciesOnly`
+  - Aborts if restoring is enabled but no tracer IDs are available
+
+### Surface tracer restoring config coupling
+
+- `Omega.SurfaceRestoring.PistonVelocity`
+  - tendency scaling
+- `Omega.SurfaceRestoring.TracersToRestore`
+  - tracer-level enable list used to build `TracerIdsToRestore`
+- `Omega.Tendencies.SurfaceTracerRestoringEnable`
+  - gates restoring tendency execution
+
+## Notes
+
+- If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer.
+- If restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error.
+- It is assumed that the incoming `TracersMonthlySurfClimoCell` fields (values and units) match the Omega state variables (i.e. conservative temperature and absolute salinity for Teos-10). If not, a pre-processing conversion should be implemented.
+- Surface tracer restoring is active everywhere if enabled. A flag to turn it off under sea ice will need to be added in later development if this feature is desired.
+- Unlike MPAS-Ocean, a `MaxDiff` clamping is not applied here. This check should instead be implemented in Ocean Validate when that is available.
+- A global scaling to ensure zero-sum has not been implemented for the surface tracer restoring, but should be added in later development.
+- At this stage, there is no temporal interpolation applied to the restoring targets, the raw `TracersMonthlySurfClimoCell` snapshot is used.

components/omega/doc/devGuide/TendencyTerms.md

@@ -41,3 +41,9 @@ implemented:
 - `TracerHighOrderHorzAdvOnCell`
 - `TracerDiffOnCell`
 - `TracerHyperDiffOnCell`
+- `SurfaceTracerRestoringOnCell`
+
+## See Also
+
+Additional information on forcing (currently wind forcing and surface tracer
+restoring) is detailed in [](omega-dev-forcing).

components/omega/doc/index.md

@@ -53,6 +53,7 @@ userGuide/PGrad
 userGuide/Timing
 userGuide/VerticalMixingCoeff
 userGuide/VertAdv
+userGuide/Forcing
 ```
 
 ```{toctree}
@@ -98,6 +99,7 @@ devGuide/PGrad
 devGuide/Timing
 devGuide/VerticalMixingCoeff
 devGuide/VertAdv
+devGuide/Forcing
 ```
 
 ```{toctree}

components/omega/doc/userGuide/AuxiliaryVariables.md

@@ -32,3 +32,10 @@ The following auxiliary variables are currently available:
 | ZonalStressCell | zonal component of wind stress on cells
 | MeridStressCell | meridional component of wind stress on cells
 | NormalStressEdge | normal component of wind stress on edge
+| SurfTracerRestoringDiffsCell | surface tracer restoring differences on cells
+| TracersMonthlySurfClimoCell | monthly climatology values to restore to for surface tracer on cells
+
+## See Also
+
+Additional information on forcing (currently wind forcing and surface tracer
+restoring) is detailed in [](omega-user-forcing).

components/omega/doc/userGuide/Forcing.md

@@ -0,0 +1,86 @@
+(omega-user-forcing)=
+
+# Forcing
+
+This page documents the user-facing configuration and behavior for current forcing in Omega:
+
+- Wind forcing
+- Surface tracer restoring
+
+## Wind forcing
+
+Wind forcing adds momentum tendency from surface wind stress.
+
+### Wind forcing configuration
+
+Wind forcing behavior is controlled by two configuration blocks:
+
+```yaml
+Omega:
+  WindStress:
+    InterpType: Isotropic
+
+  Tendencies:
+    WindForcingTendencyEnable: true
+```
+
+- `WindStress.InterpType`
+  - `Isotropic`: isotropic cell-to-edge interpolation for wind stress
+  - `Anisotropic`: anisotropic interpolation option
+- `Tendencies.WindForcingTendencyEnable`: switch to enable wind forcing tendency
+
+### Required input fields
+
+Wind forcing uses auxiliary wind-stress fields:
+
+- `WindStressZonal`
+- `WindStressMeridional`
+
+These are used to form edge-normal stress (`NormalStressEdge`) that enters
+momentum tendencies.
+
+## Surface tracer restoring
+
+Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean
+surface for selected tracers. This is implemented to mitigate drifts in chosen tracers
+(most often salinity) by nudging the model's simulated tracer values towards observed climatological values.
+This process prevents oceanic regimes from shifting away from reality due to errors in surface freshwater
+forcing (in the case of salinity restoring). Currently, it is applied everywhere when enabled.
+
+### Surface tracer restoring configuration
+
+Surface tracer restoring is controlled by two configuration blocks:
+
+```yaml
+Omega:
+  SurfaceRestoring:
+    TracersToRestore: [Temperature, Salinity]
+    PistonVelocity: 1.585e-5
+
+  Tendencies:
+    SurfaceTracerRestoringEnable: true
+```
+
+- `TracersToRestore`: list of tracer names that restoring is applied to
+- `PistonVelocity`: restoring rate coefficient
+- `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring
+
+When restoring is enabled, Omega resolves `TracersToRestore` into an internal
+list of tracer IDs and applies restoring only to tracers in that list.
+
+### Restoring target fields
+
+Surface restoring uses auxiliary fields:
+
+- `TracersMonthlySurfClimoCell`: restoring target climatological values
+- `SurfTracerRestoringDiffsCell`: computed target-minus-state differences
+
+The restoring tendency is computed at the surface layer only and is limited by
+the configured `PistonVelocity` and target-minus-state difference.
+
+## Notes
+
+- If a tracer is not listed in `TracersToRestore`, no restoring tendency is
+  applied to that tracer.
+- If surface restoring is enabled but no tracer IDs are available at tendency
+  compute-time, Omega aborts with an error.

components/omega/doc/userGuide/TendencyTerms.md

@@ -20,6 +20,7 @@ tendency terms are currently implemented:
 | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers
 | WindForcingOnEdge | forcing by wind stress, defined on edges
 | BottomDragOnEdge | bottom drag, defined on edges
+| SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells
 
 Among the internal data stored by each functor is a `bool` which can enable or
 disable the contribution of that particular term to the tendency. These flags
@@ -54,7 +55,7 @@ the currently available tendency terms:
 | WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term
 | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term
 | | BottomDragCoeff | bottom drag coefficient
-
+| SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term
 
 ## Second Order Horizontal Advection Algorithm
 
@@ -136,3 +137,8 @@ is about order 1.7 as shown in  {numref}`tracer-higher-order-convergence`:
 :width: 600 px
 Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate
 ```
+
+## See Also
+
+Additional information on forcing (currently wind forcing and surface tracer
+restoring) is detailed in [](omega-user-forcing).

components/omega/src/ocn/AuxiliaryState.cpp

@@ -3,6 +3,7 @@
 #include "Field.h"
 #include "Logging.h"
 #include "Pacer.h"
+#include "Tendencies.h"
 #include "TimeStepper.h"
 
 namespace OMEGA {
@@ -28,6 +29,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
       VorticityAux(stripDefault(Name), Mesh, VCoord),
       VelocityDel2Aux(stripDefault(Name), Mesh, VCoord),
       WindForcingAux(stripDefault(Name), Mesh),
+      SurfTracerRestAux(stripDefault(Name), Mesh, NTracers),
       TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) {
 
    GroupName = "AuxiliaryState";
@@ -43,6 +45,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
    VorticityAux.registerFields(GroupName, AuxMeshName);
    VelocityDel2Aux.registerFields(GroupName, AuxMeshName);
    WindForcingAux.registerFields(GroupName, AuxMeshName);
+   SurfTracerRestAux.registerFields(GroupName, AuxMeshName);
    TracerAux.registerFields(GroupName, AuxMeshName);
 }
 
@@ -54,6 +57,7 @@ AuxiliaryState::~AuxiliaryState() {
    VorticityAux.unregisterFields();
    VelocityDel2Aux.unregisterFields();
    WindForcingAux.unregisterFields();
+   SurfTracerRestAux.unregisterFields();
    TracerAux.unregisterFields();
 
    FieldGroup::destroy(GroupName);
@@ -406,6 +410,16 @@ I4 AuxiliaryState::exchangeHalo() {
    Err +=
        MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell);
 
+   // Performing halo exchange on individual tracers because full halo exchange
+   // on a 2D array assumes the first dimension is the vertical
+   const I4 NTracers =
+       SurfTracerRestAux.TracersMonthlySurfClimoCell.extent_int(0);
+   for (I4 LTracer = 0; LTracer < NTracers; ++LTracer) {
+      auto TracerSurfClimoCell = Kokkos::subview(
+          SurfTracerRestAux.TracersMonthlySurfClimoCell, LTracer, Kokkos::ALL);
+      Err += MeshHalo->exchangeFullArrayHalo(TracerSurfClimoCell, OnCell);
+   }
+
    return Err;
 
 } // end exchangeHalo

components/omega/src/ocn/AuxiliaryState.h

@@ -12,6 +12,7 @@
 #include "VertCoord.h"
 #include "auxiliaryVars/KineticAuxVars.h"
 #include "auxiliaryVars/LayerThicknessAuxVars.h"
+#include "auxiliaryVars/SurfTracerRestAuxVars.h"
 #include "auxiliaryVars/TracerAuxVars.h"
 #include "auxiliaryVars/VelocityDel2AuxVars.h"
 #include "auxiliaryVars/VorticityAuxVars.h"
@@ -42,6 +43,7 @@ class AuxiliaryState {
    VorticityAuxVars VorticityAux;
    VelocityDel2AuxVars VelocityDel2Aux;
    WindForcingAuxVars WindForcingAux;
+   SurfTracerRestAuxVars SurfTracerRestAux;
 
    ~AuxiliaryState();
 

components/omega/src/ocn/OceanInit.cpp

@@ -139,6 +139,15 @@ int initOmegaModules(MPI_Comm Comm) {
    PressureGrad::init();
    Eos::init();
    Tendencies::init();
+
+   // Validate SurfaceTracerRestoring configuration
+   Tendencies *DefTend = Tendencies::getDefault();
+   if (DefTend->SurfaceTracerRestoring.Enabled &&
+       DefTend->SurfaceTracerRestoring.NTracersToRestore == 0) {
+      ABORT_ERROR("OceanInit: SurfaceTracerRestoring is enabled but "
+                  "TracersToRestore is empty");
+   }
+
    TimeStepper::init2();
 
    Err = OceanState::init();