Adds the thermo coupling

components/omega/configs/Default.yml

@@ -58,6 +58,8 @@ Omega:
       Mode: Implicit
       Type: Constant
       BottomDragCoeff: 1.0e-3
+    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 as a field originating from the coupler)
 
 ## Surface stress forcing design
 
@@ -37,6 +38,70 @@ 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:
+   - `LatentHeatFluxEvap`, `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{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{sw}}$,
+    the enthalpy flux from sea-ice interactions (conduction, phase-change and enthalpy of mass flux)  Q_{\text{ice}}, the enthalpy of added freshwater mass $(\text{RainFlux} + \text{RiverRunoffFlux}) c^0_{p,sw} max(C_T(0,0,0),C_T^{\text{top}})$ (freshwater $T >= 0$), the evaporation terms ($Q_{\text{latentEvap}} + EvapFlux c^0_{p,sw} C_T^{\text{top}})
+, the enthalpy change for frozen mass $(\text{SnowFlux} + \text{IceRunoffFlux}) PotEnthalpyIce$, where PotEnthalpyIce is the potential enthalpy of solid ice (which includes phase change and enthalpy of the melted mass). To first order, $PotEnthalpyIce$ is approximated by $-LatIce$, the engineering handbook value, neglecting 0.1%.
+    These potential enthalpy terms are then scaled by $H_{\text{FluxFac}}$ to provide the layer tracer tendency.
+  - For salinity: applies the mass salt flux with the appropriate 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
+
+## Notes
+
+- Currently all forcing is applied to the surface layer only. In the future, vertical spreading of river runoff contributions will be needed.
+- `SeaIceFreshWaterFlux` is the pure freshwater mass from sea ice. The full mass flux from sea ice is `SeaIceFreshWaterFlux + SeaIceSaltFlux`
+
 ## 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 mass and tracer 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/Enthalpy fluxes (W m⁻²):**
+- `LatentHeatFluxEvap`: latent heat flux from evaporation phase change
+- `SensibleHeatFlux`: sensible heat flux
+- `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. In the future, vertical spreading of contributions from river runoff will be needed.
+- 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
@@ -57,6 +59,8 @@ the currently available tendency terms:
 | | BottomDragTendency:Mode | bottom drag mode; `Implicit` or `Explicit`
 | | BottomDragTendency:Type | bottom drag type; `Constant`
 | | BottomDragTendency:BottomDragCoeff | bottom drag coefficient
+| SfcThicknessForcingOnCell | SfcThicknessForcingTendencyEnable | enable/disable term
+| SfcTracerForcingOnCell | SfcTracerForcingTendencyEnable | enable/disable term
 | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term
 
 ## Second Order Horizontal Advection Algorithm
@@ -142,5 +146,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

@@ -330,22 +330,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp,
    }
 }
 
-Real Eos::calcPtFromCt(const Real &Sa, const Real &Ct) const {
-   if (EosChoice == EosType::Teos10Eos) {
-      return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct);
-   }
-
-   return Ct;
-}
-
-Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const {
-   if (EosChoice == EosType::Teos10Eos) {
-      return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt);
-   }
-
-   return Pt;
-}
-
 /// Define IO fields and metadata for output
 void Eos::defineFields() {
 

components/omega/src/ocn/Eos.h

@@ -663,8 +663,8 @@ class Teos10Eos {
    /// (polynomial error in [-5e-4, 6e-4] K, from GSW package).
    /// P is relative pressure (gauge pressure in Pa, i.e., absolute pressure
    /// minus the standard atmosphere).
-   KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P,
-                                       const Real SaturationFract) const {
+   static KOKKOS_FUNCTION Real calcCtFreezingTeos10(
+       const Real Sa, const Real P, const Real SaturationFract) {
       constexpr Real Sso = 35.16504;
       constexpr Real C0  = 0.017947064327968736;
       constexpr Real C1  = -6.076099099929818;
@@ -1022,10 +1022,47 @@ class Eos {
                                   const Array2DReal &SpecVol);
 
    /// Convert Conservative Temperature to potential temperature
-   Real calcPtFromCt(const Real &Sa, const Real &Ct) const;
+   /// For TEOS-10, uses the TEOS-10 polynomial
+   /// For other EOS choices, conservative temperature is equal to potential
+   /// temperature
+   KOKKOS_FUNCTION Real calcPtFromCt(const Real &Sa, const Real &Ct) const {
+      if (EosChoice == EosType::Teos10Eos) {
+         return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct);
+      }
+      return Ct;
+   }
+
+   /// Convert potential temperature to Conservative Temperature.
+   /// For TEOS-10, uses the TEOS-10 polynomial
+   /// For other EOS choices, potential temperature equals conservative
+   /// temperature
+   KOKKOS_FUNCTION Real calcCtFromPt(const Real &Sa, const Real &Pt) const {
+      if (EosChoice == EosType::Teos10Eos) {
+         return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt);
+      }
+      return Pt;
+   }
 
-   /// Convert potential temperature to Conservative Temperature
-   Real calcCtFromPt(const Real &Sa, const Real &Pt) const;
+   /// Calculate freezing temperature of seawater.
+   /// For TEOS-10, uses the Roquet et al. 75-term polynomial.
+   /// For LinearEos, uses a simple linear salinity-dependent approximation.
+   /// For ConstantEos, returns a constant approximate ocean freezing point.
+   static KOKKOS_FUNCTION Real calcCtFreezing(EosType Choice, const Real Sa,
+                                              const Real P,
+                                              const Real SaturationFract) {
+      if (Choice == EosType::Teos10Eos) {
+         return Teos10Eos::calcCtFreezingTeos10(Sa, P, SaturationFract);
+      }
+      if (Choice == EosType::LinearEos) {
+         // Linear salinity-dependent freezing point; coefficient -0.054
+         // degC/PSU with absolute-to-practical salinity conversion (g/kg ->
+         // PSU).
+         constexpr Real Coeff = -0.054_Real;
+         return Coeff * Sa / Psu2Gpkg;
+      }
+      // ConstantEos: constant approximate ocean freezing point (degC)
+      return -1.9_Real;
+   }
 
    /// Initialize EOS from config and mesh
    static void init();