Adds the thermo coupling

components/omega/configs/Default.yml

@@ -55,6 +55,8 @@ Omega:
     SurfaceTracerRestoringEnable: false
     BottomDragTendencyEnable: false
     BottomDragCoeff: 0.0
+    SfcThicknessForcingTendencyEnable: false
+    SfcTracerForcingTendencyEnable: false
     TracerHorzAdvTendencyEnable: true
     TracerDiffTendencyEnable: true
     EddyDiff2: 10.0

components/omega/doc/devGuide/Forcing.md

@@ -6,7 +6,8 @@ This page describes design and implementation details for forcing-related
 pathways in Omega, currently this includes:
 
 - Surface stress forcing (e.g. wind stress)
-- Surface tracer restoring
+- Surface thickness and tracer flux forcing (actively coupled or data-forced)
+- Surface tracer restoring (soon to be ported)
 
 ## Surface stress forcing design
 
@@ -37,6 +38,65 @@ pathways in Omega, currently this includes:
 - `Omega.Tendencies.SfcStressForcingTendencyEnable`
   - gates execution of surface stress forcing tendency kernel
 
+## Surface thickness and tracer flux forcing design
+
+### Surface thickness and tracer flux forcing data flow
+
+**Thickness equation pathway:**
+
+1. External fields provide freshwater and salt flux components:
+   - `SnowFlux`, `RainFlux`, `EvaporationFlux`
+   - `SeaIceFreshWaterFlux`, `IceRunoffFlux`, `RiverRunoffFlux`
+   - `SeaIceSaltFlux`
+2. `Forcing` stores the flux fields in `TracerForcingVars`
+3. The tendency term `SfcThicknessForcingOnCell` sums the freshwater and salt mass fluxes and applies them to
+the surface layer pseudo-thickness.
+
+**Tracer equation pathway:**
+
+1. External fields provide heat and salt flux components:
+   - `LatentHeatFlux`, `SensibleHeatFlux`
+   - `LongWaveHeatFluxUp`, `LongWaveHeatFluxDown`
+   - `SeaIceHeatFlux`, `ShortWaveHeatFlux`
+   - mass fluxes which add energy changes (`SnowFlux`, `RainFlux`, `IceRunoffFlux`, `RiverRunoffFlux`)
+   - `SeaIceSaltFlux`
+2. `Forcing` stores the flux fields in `TracerForcingVars`
+3. The tendency term `SfcTracerForcingOnCell` converts the summed external heat fluxes to a conservative-temperature tendency,
+  and applies the external sea-ice salt flux to the top layer salt content thus impacting salinity.
+
+### Surface thickness and tracer flux forcing key classes/components
+
+- `TracerForcingVars`
+  - Stores 13 coupled flux cell-centered fields: 6 freshwater fluxes, 6 heat
+    fluxes, and 1 salt flux component
+  - Fields initialized to zero and registered in `Forcing` field group
+- `SfcThicknessForcingOnCell` tendency term
+  - Computes the layer mass contribution (converted to pseudo-thickness): $\sum (\text{SnowFlux} + \text{RainFlux} + \text{EvaporationFlux} + \text{SeaIceFreshWaterFlux} + \text{IceRunoffFlux} + \text{RiverRunoffFlux} + \text{SeaIceSaltFlux}) / \rho_{sw}$
+  - Applied only at surface layer (top active layer) using `MinLayerCell`
+- `SfcTracerForcingOnCell` tendency term
+  - For temperature: adds the direct heat fluxes
+    $Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$
+, the phase change and enthalpy of added mass $(\text{RainFlux} + \text{RiverRunoffFlux}) c^0_{p,sw} C_T^{\text{top}} + (\text{SnowFlux} + \text{IceRunoffFlux})(c^0_{p,sw} C_T^{\text{frz}} - L_{\text{ice}})$,
+    (where $C_T^{\text{frz}}$ is from EOS at top-layer salinity and pressure),
+    and scales by $H_{\text{FluxFac}}$.
+  - For salinity: applies salt flux with unit conversion: $\text{SeaIceSaltFlux} \times S_{\text{FluxFac}}$
+  - Applied only at surface layer using `MinLayerCell`
+  - Uses tracer index validation to apply to specific tracers only
+- `Forcing`
+  - Manages `TracerForcingVars` instance
+- `Tendencies`
+  - Calls `SfcThicknessForcingOnCell` in `computePseudoThicknessTendenciesOnly`
+  - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring
+
+### Surface thickness and tracer flux forcing config coupling
+
+- `Omega.Tendencies.SfcThicknessForcingTendencyEnable`
+  - gates execution of coupled flux thickness kernel
+  - controls freshwater and salt flux forcing on sea surface height
+- `Omega.Tendencies.SfcTracerForcingTendencyEnable`
+  - gates execution of coupled flux tracer kernel
+  - controls direct heat flux forcing on temperature and salt flux forcing on salinity
+
 ## Surface tracer restoring design
 
 ### Surface tracer restoring data flow

components/omega/doc/devGuide/TendencyTerms.md

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

components/omega/doc/userGuide/Forcing.md

@@ -5,6 +5,7 @@
 This page documents the user-facing configuration and behavior for current forcing in Omega:
 
 - Surface stress forcing (e.g. wind stress)
+- Coupled tracer flux forcing (mass, energy and salt)
 - Surface tracer restoring
 
 ## Surface stress forcing
@@ -39,6 +40,74 @@ Surface stress forcing uses surface stress input fields:
 These are stored in forcing variables and used to form edge-normal stress
 (`NormalStressEdge`) that enters momentum tendencies.
 
+## Surface thickness and tracer flux forcing
+
+Surface thickness and tracer flux forcing applies ocean-atmosphere and ocean-sea ice fluxes from the other model
+components (atmosphere, sea ice) to the thickness and tracer equations. This enables
+the ocean to respond to heat, freshwater, and salt exchanges at the surface. These fluxes can be from data or (active) coupled components.
+
+### Surface thickness and tracer flux forcing configuration
+
+Surface thickness and tracer flux forcing is controlled by two configuration flags:
+
+```yaml
+Omega:
+  Tendencies:
+    SfcThicknessForcingTendencyEnable: false
+    SfcTracerForcingTendencyEnable: false
+```
+
+- `Tendencies.SfcThicknessForcingTendencyEnable`: enables coupled freshwater and salt flux forcing on thickness
+- `Tendencies.SfcTracerForcingTendencyEnable`: enables coupled heat and salt flux forcing on tracers
+
+
+### Required input fields
+
+Coupled tracer flux forcing uses 13 auxiliary fields organized by type:
+
+**Freshwater mass fluxes (kg m⁻² s⁻¹):**
+- `SnowFlux`: precipitation from snow
+- `RainFlux`: precipitation from rain
+- `EvaporationFlux`: evaporative water loss
+- `SeaIceFreshWaterFlux`: freshwater mass flux from sea-ice melt or formation
+- `IceRunoffFlux`: runoff from land ice
+- `RiverRunoffFlux`: runoff from rivers
+
+**Heat fluxes (W m⁻²):**
+- `LatentHeatFlux`: latent heat transfer
+- `SensibleHeatFlux`: sensible heat transfer
+- `LongWaveHeatFluxUp`: upward longwave radiation
+- `LongWaveHeatFluxDown`: downward longwave radiation
+- `SeaIceHeatFlux`: heat/energy from sea-ice interaction (incl. enthalpy of meltwater)
+- `ShortWaveHeatFlux`: shortwave (solar) radiation
+
+**Salt mass flux (kg m⁻² s⁻¹):**
+- `SeaIceSaltFlux`: salt flux from sea-ice formation/melt processes
+
+These fields are populated by external coupling components (typically atmosphere
+and ice models). Omega assumes the incoming values match the documented units.
+For now, there are assumed to come from a `forcing.nc` file, but later will be provided
+by the equivalent `ocn_comp_mct.F`.
+
+### Notes
+
+- Coupled fluxes are applied only at the surface layer (top active layer) for each cell.
+- Pseudo-thickness tendency is computed from the (six) freshwater mass fluxes and the salt mass flux
+  `SeaIceSaltFlux`, converted to a pseudo-thickness change.
+- Temperature tendency is computed from direct heat flux plus
+  mass-flux enthalpy terms, converted to conservative-temperature tendency via
+  $H_{\text{FluxFac}} = 1.0 / (\rho_{sw} c^0_{p,sw})$ where $c^0_{p,sw}$ is the reference
+  specific heat of seawater defined by TEOS-10. The enthalpy associated with mass fluxes is currently hard-coded to SST for liquid fluxes and the freezing temperature for solid fluxes (which are melted using a constant latent heat of fusion). Note that the enthalpy of liquid meltwater from sea ice is already included in `SeaIceHeatFlux`.
+- Salinity tendency from `SeaIceSaltFlux` is scaled by
+  $S_{\text{FluxFac}} = 1.0e3 / \rho_{sw}$ to account for unit conversion from
+  kg/(m²·s) to salinity units (g/kg).
+- Fluxes are assumed to be in the documented units (i.e. net mass fluxes);
+  any unit conversion should be performed by the coupling component before providing flux
+  values to Omega.
+- The reference density used here ($\rho_{sw}$) is not a Boussinesq density, it is the
+  conversion factor from mass to pseudo-thickness.
+- No iceberg fluxes are included for now.
+
 ## Surface tracer restoring
 
 Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean

components/omega/doc/userGuide/TendencyTerms.md

@@ -20,6 +20,8 @@ tendency terms are currently implemented:
 | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers
 | SfcStressForcingOnEdge | forcing by surface stress (e.g. wind), defined on edges
 | BottomDragOnEdge | bottom drag, defined on edges
+| SfcThicknessForcingOnCell | surface pseudo-thickness forcing from coupled freshwater and salt fluxes, defined on cells
+| SfcTracerForcingOnCell | surface tracer forcing from coupled heat and salt fluxes, with direct heat always and mass-flux enthalpy terms gated by thickness forcing, defined on cells
 | SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells
 
 Among the internal data stored by each functor is a `bool` which can enable or
@@ -55,6 +57,8 @@ the currently available tendency terms:
 | SfcStressForcingOnEdge | SfcStressForcingTendencyEnable | enable/disable term
 | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term
 | | BottomDragCoeff | bottom drag coefficient
+| SfcThicknessForcingOnCell | SfcThicknessForcingTendencyEnable | enable/disable term
+| SfcTracerForcingOnCell | SfcTracerForcingTendencyEnable | enable/disable term
 | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term
 
 ## Second Order Horizontal Advection Algorithm
@@ -140,5 +144,6 @@ Tracer higer order convergence example of a cosine bell advected on a sphere sho
 
 ## See Also
 
-Additional information on forcing (currently wind forcing and surface tracer
-restoring) is detailed in [](omega-user-forcing).
+Additional information on forcing, including surface stress forcing,
+surface thickness and tracer flux forcing, and surface tracer restoring, is detailed in
+[](omega-user-forcing).

components/omega/src/ocn/Eos.cpp

@@ -343,6 +343,20 @@ Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const {
    return Pt;
 }
 
+Real Eos::calcCtFreezing(const Real Sa, const Real P,
+                         const Real SaturationFract) const {
+   if (EosChoice == EosType::Teos10Eos) {
+      return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract);
+   }
+
+   ABORT_ERROR("Eos::calcCtFreezing: CT freezing temperature is only "
+               "implemented for TEOS-10. Support for the current EOS "
+               "choice has not yet been developed.");
+   // most likely I'd implement a polynomial here for non-teos10 e.g.
+   // return 0.0 - 0.0575 * Sa + 1.710523e-3 * sqrt(Sa^3) - 2.154996e-4 * Sa^2
+   return 0.0;
+}
+
 /// Define IO fields and metadata for output
 void Eos::defineFields() {
 

components/omega/src/ocn/Eos.h

@@ -756,6 +756,12 @@ class Eos {
    /// Convert potential temperature to Conservative Temperature
    Real calcCtFromPt(const Real &Sa, const Real &Pt) const;
 
+   /// Calculate freezing Conservative Temperature for TEOS-10.
+   /// Aborts if EOS is not TEOS-10: CT freezing is not yet implemented
+   /// for other equation-of-state choices.
+   Real calcCtFreezing(const Real Sa, const Real P,
+                       const Real SaturationFract) const;
+
    /// Initialize EOS from config and mesh
    static void init();
 

components/omega/src/ocn/Forcing.cpp

@@ -29,18 +29,22 @@ static std::string stripDefault(const std::string &Name) {
 // mesh/halo.
 Forcing::Forcing(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo)
     : Name(stripDefault(Name)), SfcStressForcing(stripDefault(Name), Mesh),
-      Mesh(Mesh), MeshHalo(MeshHalo) {}
+      TracerForcing(stripDefault(Name), Mesh), Mesh(Mesh), MeshHalo(MeshHalo) {}
 
 // Destructor. Unregisters fields from IO streams.
 Forcing::~Forcing() { unregisterFields(); }
 
 // Register surface stress fields with IO streams for a given mesh.
 void Forcing::registerFields(const std::string &MeshName) const {
    SfcStressForcing.registerFields(MeshName);
+   TracerForcing.registerFields(MeshName);
 }
 
 // Unregister surface stress fields from IO streams.
-void Forcing::unregisterFields() const { SfcStressForcing.unregisterFields(); }
+void Forcing::unregisterFields() const {
+   SfcStressForcing.unregisterFields();
+   TracerForcing.unregisterFields();
+}
 
 // Create and register a non-default forcing instance.
 Forcing *Forcing::create(const std::string &Name, const HorzMesh *Mesh,
@@ -153,7 +157,8 @@ void Forcing::computeSfcStressForcingOnEdge() const {
    Pacer::stop("Forcing:edge1", 2);
 }
 
-// Exchange halo for surface stress cell fields.
+// Exchange halo for surface stress cell fields. Only needed for variables that
+// need information beyond cell-centered values.
 I4 Forcing::exchangeHalo() const {
    I4 Err = 0;
 

components/omega/src/ocn/Forcing.h

@@ -17,6 +17,7 @@
 #include "Halo.h"
 #include "HorzMesh.h"
 #include "forcingVars/SfcStressForcingVars.h"
+#include "forcingVars/TracerForcingVars.h"
 
 #include <map>
 #include <memory>
@@ -32,6 +33,7 @@ class Forcing {
    std::string Name; ///< Name identifier for this forcing instance
 
    SfcStressForcingVars SfcStressForcing; ///< Surface stress forcing variables
+   TracerForcingVars TracerForcing; ///< Tracer forcing vars (thickness and T,S)
 
    ~Forcing();
 

components/omega/src/ocn/GlobalConstants.h

@@ -115,11 +115,12 @@ constexpr Real Pa2Db      = 1.0e-4;    // Pascal to Decibar
 constexpr Real Cm2M       = 1.0e-2;    // Centimeters to meters
 constexpr Real M2Cm       = 1.0e2;     // Meters to centimeters
 constexpr Real HFluxFac =
-    1.0 / (RhoSw * CpSw);         // Heat flux (W/m^2) to temp flux (C*m/s)
+    1.0 / (RhoSw * Cp0Sw); // Heat flux (W/m^2) to Conserv Temp flux (C*m/s)
 constexpr Real FwFluxFac = 1.e-6; // Fw flux (kg/m^2/s) to salt((msu/psu)*m/s)
 constexpr Real SaltFac =
-    -OcnRefSal * FwFluxFac;    // Fw flux (kg/m^2/s) to salt flux (msu*m/s)
-constexpr Real SFluxFac = 1.0; // Salt flux (kg/m^2/s) to salt flux (msu*m/s)
+    -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s)
+constexpr Real SFluxFac =
+    1.e3 / RhoSw; // Salt flux (kg/m^2/s) to salinity flux (m*(g/kg)/s)
 
 } // namespace OMEGA
 #endif