From ec5655910d3d3a2e9753ef3819fd64b83d418c42 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 19 Jun 2026 14:18:42 -0700 Subject: [PATCH 01/27] port forcing tendencies (thickness and tracers) + update doc, yml --- components/omega/configs/Default.yml | 2 + components/omega/doc/devGuide/Forcing.md | 57 ++++- components/omega/doc/userGuide/Forcing.md | 68 ++++++ .../omega/doc/userGuide/TendencyTerms.md | 9 +- components/omega/src/ocn/Forcing.cpp | 35 ++- components/omega/src/ocn/Forcing.h | 2 + components/omega/src/ocn/GlobalConstants.h | 7 +- components/omega/src/ocn/Tendencies.cpp | 73 +++++++ components/omega/src/ocn/Tendencies.h | 2 + components/omega/src/ocn/TendencyTerms.cpp | 11 + components/omega/src/ocn/TendencyTerms.h | 86 ++++++++ .../src/ocn/forcingVars/TracerForcingVars.cpp | 203 ++++++++++++++++++ .../src/ocn/forcingVars/TracerForcingVars.h | 49 +++++ .../test/timeStepping/TimeStepperTest.cpp | 2 + 14 files changed, 598 insertions(+), 8 deletions(-) create mode 100644 components/omega/src/ocn/forcingVars/TracerForcingVars.cpp create mode 100644 components/omega/src/ocn/forcingVars/TracerForcingVars.h diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..e33403c7cc91 100644 --- a/components/omega/configs/Default.yml +++ b/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 diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 97d38d8ae42d..18e9c2072190 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/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 flux forcing (actively coupled or data-forced) +- Surface tracer restoring (soon to be ported) ## Surface stress forcing design @@ -37,6 +38,60 @@ pathways in Omega, currently this includes: - `Omega.Tendencies.SfcStressForcingTendencyEnable` - gates execution of surface stress forcing tendency kernel +## Surface flux forcing design + +### Surface 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` + - `SeaIceSaltFlux`, `SnowFlux`, `IceRunoffFlux` +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 salinity (g/kg) in the surface layer. + +### Surface 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 freshwater flux contribution: $\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: computes the sum of the six heat-flux fields and scales it 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 `computeThicknessTendenciesOnly` + - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring + +### Surface 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 heat flux forcing on temperature and salt flux forcing on salinity + ## Surface tracer restoring design ### Surface tracer restoring data flow diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index d8fac4383730..b49eabd0586f 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/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 flux forcing - Surface tracer restoring ## Surface stress forcing @@ -39,6 +40,73 @@ 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 flux forcing + +Surface 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 flux forcing configuration + +Surface 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 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 input 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 from sea-ice interaction +- `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 the sum of the six heat-flux fields, + 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. [soon to be updated with latent heat and enthalpy of liquid water] +- 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 diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 3db2b4c10098..1259d3387054 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/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, 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 flux forcing, and surface tracer restoring, is detailed in +[](omega-user-forcing). diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 39cbe9311b18..3a3a54064dff 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -30,7 +30,7 @@ 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(); } @@ -38,10 +38,14 @@ 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, @@ -171,6 +175,33 @@ I4 Forcing::exchangeHalo() const { Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SnowFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RainFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.EvaporationFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo( + TracerForcing.SeaIceFreshWaterFluxCell, OnCell); + Err += + MeshHalo->exchangeFullArrayHalo(TracerForcing.IceRunoffFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RiverRunoffFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LatentHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SensibleHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LongWaveHeatFluxUpCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo( + TracerForcing.LongWaveHeatFluxDownCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.ShortWaveHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceSaltFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SurfInsituTemperature, + OnCell); + return Err; } diff --git a/components/omega/src/ocn/Forcing.h b/components/omega/src/ocn/Forcing.h index 5fdae7e550b8..fda7b91d414e 100644 --- a/components/omega/src/ocn/Forcing.h +++ b/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 #include @@ -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(); diff --git a/components/omega/src/ocn/GlobalConstants.h b/components/omega/src/ocn/GlobalConstants.h index d31bae11f489..0fe2adb1f54e 100644 --- a/components/omega/src/ocn/GlobalConstants.h +++ b/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 diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 5bb131f32ab6..34e91a59053c 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -297,6 +297,18 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config } } + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + this->SfcThicknessForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcThicknessForcingTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + this->SfcTracerForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcTracerForcingTendencyEnable not found in TendConfig"); + if (this->TracerDiffusion.Enabled) { Err += TendConfig.get("EddyDiff2", this->TracerDiffusion.EddyDiff2); CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff2 not found in TendConfig"); @@ -475,6 +487,8 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord), VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), ExplicitBottomDrag(Mesh, VCoord), + SfcThicknessForcing(Mesh, VCoord), + SfcTracerForcing(Mesh, VCoord, Tracers::IndxTemp, Tracers::IndxSalt), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), @@ -563,6 +577,7 @@ void Tendencies::computePseudoThicknessTendenciesOnly( OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); OMEGA_SCOPE(LocThicknessFluxDiv, PseudoThicknessFluxDiv); + OMEGA_SCOPE(LocSfcThicknessForcing, SfcThicknessForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -606,6 +621,32 @@ void Tendencies::computePseudoThicknessTendenciesOnly( VAdv->computePseudoThicknessVAdvTend(PseudoThicknessTend); Pacer::stop("Tend:computePseudoThicknessVAdvTend", 2); + if (LocSfcThicknessForcing.Enabled) { + Pacer::start("Tend:sfcThicknessForcing", 2); + const auto *ForcingState = Forcing::getDefault(); + + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &EvaporationFlux = + ForcingState->TracerForcing.EvaporationFluxCell; + const auto &SeaIceFreshWaterFlux = + ForcingState->TracerForcing.SeaIceFreshWaterFluxCell; + const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &RiverRunoffFlux = + ForcingState->TracerForcing.RiverRunoffFluxCell; + const auto &SeaIceSaltFlux = + ForcingState->TracerForcing.SeaIceSaltFluxCell; + + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + LocSfcThicknessForcing(LocPseudoThicknessTend, ICell, SnowFlux, + RainFlux, EvaporationFlux, + SeaIceFreshWaterFlux, IceRunoffFlux, + RiverRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcThicknessForcing", 2); + } + if (CustomThicknessTend) { Pacer::start("Tend:customThicknessTend", 2); CustomThicknessTend(LocPseudoThicknessTend, State, AuxState, @@ -825,6 +866,7 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); + OMEGA_SCOPE(LocSfcTracerForcing, SfcTracerForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -946,6 +988,37 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:surfaceTracerRestoring", 2); } + // compute tracer forcing tendency + if (LocSfcTracerForcing.Enabled) { + Pacer::start("Tend:sfcTracerForcing", 2); + const auto *ForcingState = Forcing::getDefault(); + const auto &LatentHeatFlux = + ForcingState->TracerForcing.LatentHeatFluxCell; + const auto &SensibleHeatFlux = + ForcingState->TracerForcing.SensibleHeatFluxCell; + const auto &LongWaveHeatFluxUp = + ForcingState->TracerForcing.LongWaveHeatFluxUpCell; + const auto &LongWaveHeatFluxDown = + ForcingState->TracerForcing.LongWaveHeatFluxDownCell; + const auto &SeaIceHeatFlux = + ForcingState->TracerForcing.SeaIceHeatFluxCell; + const auto &ShortWaveHeatFlux = + ForcingState->TracerForcing.ShortWaveHeatFluxCell; + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &SeaIceSaltFlux = + ForcingState->TracerForcing.SeaIceSaltFluxCell; + + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + LocSfcTracerForcing( + LocTracerTend, ICell, LatentHeatFlux, SensibleHeatFlux, + LongWaveHeatFluxUp, LongWaveHeatFluxDown, SeaIceHeatFlux, + ShortWaveHeatFlux, SnowFlux, IceRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcTracerForcing", 2); + } + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index c277547c9e62..13be243aa8b7 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -74,6 +74,8 @@ class Tendencies { VelocityHyperDiffOnEdge VelocityHyperDiff; SfcStressForcingOnEdge SfcStressForcing; BottomDragOnEdge ExplicitBottomDrag; + SfcThicknessForcingOnCell SfcThicknessForcing; + SfcTracerForcingOnCell SfcTracerForcing; TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 5b142888142c..f2c6bdca73ad 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -71,6 +71,17 @@ BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh, NVertLayers(VCoord->NVertLayers), EdgeMask(VCoord->EdgeMask), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +SfcThicknessForcingOnCell::SfcThicknessForcingOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + +SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord, + I4 TempTracerIndex, + I4 SaltTracerIndex) + : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) : HorzontalMesh(Mesh), VerticalCoord(VCoord), diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 548832b290bf..1745ee713734 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -372,6 +372,92 @@ class BottomDragOnEdge { Array1DI4 MaxLayerEdgeTop; }; +/// Coupled freshwater flux forcing for thickness equation. +class SfcThicknessForcingOnCell { + public: + bool Enabled = false; + + SfcThicknessForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 ICell, + const Array1DReal &SnowFlux, + const Array1DReal &RainFlux, + const Array1DReal &EvaporationFlux, + const Array1DReal &SeaIceFreshWaterFlux, + const Array1DReal &IceRunoffFlux, + const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { + + const I4 KTop = MinLayerCell(ICell); + if (KTop > MaxLayerCell(ICell)) { + return; + } + + const Real FreshWaterFlux = SnowFlux(ICell) + RainFlux(ICell) + + EvaporationFlux(ICell) + + SeaIceFreshWaterFlux(ICell) + + IceRunoffFlux(ICell) + RiverRunoffFlux(ICell); + + Tend(ICell, KTop) += (FreshWaterFlux + SeaIceSaltFlux(ICell)) / RhoSw; + } + + private: + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + +/// Coupled surface flux forcing for active tracers. +class SfcTracerForcingOnCell { + public: + bool Enabled = false; + + SfcTracerForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord, + I4 TempTracerIndex, I4 SaltTracerIndex); + + KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 ICell, + const Array1DReal &LatentHeatFlux, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, + const Array1DReal &ShortWaveHeatFlux, + const Array1DReal &SnowFlux, + const Array1DReal &IceRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { + + const I4 KTop = MinLayerCell(ICell); + if (KTop > MaxLayerCell(ICell)) { + return; + } + + if (TempIndex >= 0) { + const Real HeatFlux = LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + + LongWaveHeatFluxUp(ICell) + + LongWaveHeatFluxDown(ICell) + + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell); + // + + // (RainFlux(ICell) + RiverRunoffFlux(ICell)) * + // Cp0Sw * TracerCell(TempIndex, ICell, KTop) + + // (SnowFlux(ICell) + IceRunoffFlux(ICell)) * + // (Cp0Sw * Eos.Ctfreez - LatIce; + + Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFactor; + } + + if (SaltIndex >= 0) { + Tend(SaltIndex, ICell, KTop) += SeaIceSaltFlux(ICell) * SFluxFactor; + } + } + + private: + I4 TempIndex; + I4 SaltIndex; + Real HFluxFactor; + Real SFluxFactor; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + // Tracer horizontal advection term class TracerHorzAdvOnCell { public: diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp new file mode 100644 index 000000000000..38016216a410 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -0,0 +1,203 @@ +#include "TracerForcingVars.h" +#include "Eos.h" +#include "Field.h" +#include "Tracers.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +TracerForcingVars::TracerForcingVars(const std::string &Suffix, + const HorzMesh *Mesh) + : SnowFluxCell("snowFlux" + Suffix, Mesh->NCellsSize), + RainFluxCell("rainFlux" + Suffix, Mesh->NCellsSize), + EvaporationFluxCell("evaporationFlux" + Suffix, Mesh->NCellsSize), + SeaIceFreshWaterFluxCell("seaIceFreshWaterFlux" + Suffix, + Mesh->NCellsSize), + IceRunoffFluxCell("iceRunoffFlux" + Suffix, Mesh->NCellsSize), + RiverRunoffFluxCell("riverRunoffFlux" + Suffix, Mesh->NCellsSize), + LatentHeatFluxCell("latentHeatFlux" + Suffix, Mesh->NCellsSize), + SensibleHeatFluxCell("sensibleHeatFlux" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxUpCell("longWaveHeatFluxUp" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxDownCell("longWaveHeatFluxDown" + Suffix, + Mesh->NCellsSize), + SeaIceHeatFluxCell("seaIceHeatFlux" + Suffix, Mesh->NCellsSize), + ShortWaveHeatFluxCell("shortWaveHeatFlux" + Suffix, Mesh->NCellsSize), + SeaIceSaltFluxCell("seaIceSalinityFlux" + Suffix, Mesh->NCellsSize), + SurfInsituTemperature("surfInsituTemperature" + Suffix, + Mesh->NCellsSize) { + deepCopy(SnowFluxCell, 0.0_Real); + deepCopy(RainFluxCell, 0.0_Real); + deepCopy(EvaporationFluxCell, 0.0_Real); + deepCopy(SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(IceRunoffFluxCell, 0.0_Real); + deepCopy(RiverRunoffFluxCell, 0.0_Real); + deepCopy(LatentHeatFluxCell, 0.0_Real); + deepCopy(SensibleHeatFluxCell, 0.0_Real); + deepCopy(LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(SeaIceHeatFluxCell, 0.0_Real); + deepCopy(ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(SeaIceSaltFluxCell, 0.0_Real); + deepCopy(SurfInsituTemperature, 0.0_Real); +} + +void TracerForcingVars::registerFields(const std::string &MeshName) const { + const Real FillValue = -9.99e30; + const int NDims = 1; + std::vector DimNames(NDims); + + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + DimNames[0] = "NCells" + DimSuffix; + + auto SnowFluxField = Field::create( + SnowFluxCell.label(), "snow freshwater flux", "kg m^-2 s^-1", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + auto RainFluxField = Field::create( + RainFluxCell.label(), "rain freshwater flux", "kg m^-2 s^-1", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + auto EvaporationFluxField = Field::create( + EvaporationFluxCell.label(), "evaporation freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + auto SeaIceFreshWaterFluxField = Field::create( + SeaIceFreshWaterFluxCell.label(), "sea-ice freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + auto IceRunoffFluxField = Field::create( + IceRunoffFluxCell.label(), "ice runoff freshwater flux", "kg m^-2 s^-1", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + auto RiverRunoffFluxField = Field::create( + RiverRunoffFluxCell.label(), "river runoff freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + + auto LatentHeatFluxField = Field::create( + LatentHeatFluxCell.label(), "latent heat flux", "W m^-2", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + auto SensibleHeatFluxField = Field::create( + SensibleHeatFluxCell.label(), "sensible heat flux", "W m^-2", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + auto LongWaveHeatFluxUpField = Field::create( + LongWaveHeatFluxUpCell.label(), "upward longwave heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + auto LongWaveHeatFluxDownField = Field::create( + LongWaveHeatFluxDownCell.label(), "downward longwave heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), FillValue, NDims, DimNames); + auto SeaIceHeatFluxField = Field::create( + SeaIceHeatFluxCell.label(), "sea-ice heat flux", "W m^-2", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + auto ShortWaveHeatFluxField = Field::create( + ShortWaveHeatFluxCell.label(), "shortwave heat flux", "W m^-2", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + + auto SeaIceSaltFluxField = Field::create( + SeaIceSaltFluxCell.label(), "sea-ice salt flux", "kg m^-2 s^-1", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + + auto SurfInsituTemperatureField = Field::create( + SurfInsituTemperature.label(), + "insitu (potential) temperature at surface layer", "degrees Celsius", "", + std::numeric_limits::lowest(), std::numeric_limits::max(), + FillValue, NDims, DimNames); + + FieldGroup::addFieldToGroup(SnowFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(RainFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(EvaporationFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceFreshWaterFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(IceRunoffFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(RiverRunoffFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LatentHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SensibleHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LongWaveHeatFluxUpCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LongWaveHeatFluxDownCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(ShortWaveHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceSaltFluxCell.label(), "Forcing"); + + SnowFluxField->attachData(SnowFluxCell); + RainFluxField->attachData(RainFluxCell); + EvaporationFluxField->attachData(EvaporationFluxCell); + SeaIceFreshWaterFluxField->attachData(SeaIceFreshWaterFluxCell); + IceRunoffFluxField->attachData(IceRunoffFluxCell); + RiverRunoffFluxField->attachData(RiverRunoffFluxCell); + LatentHeatFluxField->attachData(LatentHeatFluxCell); + SensibleHeatFluxField->attachData(SensibleHeatFluxCell); + LongWaveHeatFluxUpField->attachData(LongWaveHeatFluxUpCell); + LongWaveHeatFluxDownField->attachData(LongWaveHeatFluxDownCell); + SeaIceHeatFluxField->attachData(SeaIceHeatFluxCell); + ShortWaveHeatFluxField->attachData(ShortWaveHeatFluxCell); + SurfInsituTemperatureField->attachData(SurfInsituTemperature); + SeaIceSaltFluxField->attachData(SeaIceSaltFluxCell); +} + +void TracerForcingVars::unregisterFields() const { + Field::destroy(SnowFluxCell.label()); + Field::destroy(RainFluxCell.label()); + Field::destroy(EvaporationFluxCell.label()); + Field::destroy(SeaIceFreshWaterFluxCell.label()); + Field::destroy(IceRunoffFluxCell.label()); + Field::destroy(RiverRunoffFluxCell.label()); + Field::destroy(LatentHeatFluxCell.label()); + Field::destroy(SensibleHeatFluxCell.label()); + Field::destroy(LongWaveHeatFluxUpCell.label()); + Field::destroy(LongWaveHeatFluxDownCell.label()); + Field::destroy(SeaIceHeatFluxCell.label()); + Field::destroy(ShortWaveHeatFluxCell.label()); + Field::destroy(SeaIceSaltFluxCell.label()); + Field::destroy(SurfInsituTemperature.label()); +} + +void TracerForcingVars::computeSurfInsituTemp(const Array3DReal &TracerArray, + const VertCoord *VCoord, + const Eos *EosInst) const { + const int IndxTemp = Tracers::IndxTemp; + const int IndxSalt = Tracers::IndxSalt; + + // Skip computation if temperature or salinity tracers are not defined + if (IndxTemp < 0 || IndxSalt < 0) { + return; + } + + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocSurfInsituTemp, SurfInsituTemperature); + + int NCellsOwned = SurfInsituTemperature.extent_int(0); + + parallelFor( + "TracerForcing:computeSurfInsituTemp", {NCellsOwned}, + KOKKOS_LAMBDA(int ICell) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + // Only compute for valid ocean cells + if (KMin <= KMax) { + const Real ConservTemp = TracerArray(IndxTemp, ICell, KMin); + const Real AbsSalinity = TracerArray(IndxSalt, ICell, KMin); + + // Call EOS function to compute potential temperature from + // conservative temperature at surface (reference pressure = 0) + LocSurfInsituTemp(ICell) = + EosInst->calcPtFromCt(AbsSalinity, ConservTemp); + } + }); +} +} // namespace OMEGA diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.h b/components/omega/src/ocn/forcingVars/TracerForcingVars.h new file mode 100644 index 000000000000..1a0747121ea2 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.h @@ -0,0 +1,49 @@ +#ifndef OMEGA_TRACER_FORCING_H +#define OMEGA_TRACER_FORCING_H + +#include "DataTypes.h" +#include "HorzMesh.h" + +#include + +namespace OMEGA { + +// Forward declarations. Full definitions not needed in this header since only +// pointers are used. +class VertCoord; +class Eos; + +class TracerForcingVars { + public: + Array1DReal SnowFluxCell; + Array1DReal RainFluxCell; + Array1DReal EvaporationFluxCell; + Array1DReal SeaIceFreshWaterFluxCell; + Array1DReal IceRunoffFluxCell; + Array1DReal RiverRunoffFluxCell; + + Array1DReal LatentHeatFluxCell; + Array1DReal SensibleHeatFluxCell; + Array1DReal LongWaveHeatFluxUpCell; + Array1DReal LongWaveHeatFluxDownCell; + Array1DReal SeaIceHeatFluxCell; + Array1DReal ShortWaveHeatFluxCell; + + Array1DReal SeaIceSaltFluxCell; + + Array1DReal SurfInsituTemperature; + + TracerForcingVars(const std::string &Suffix, const HorzMesh *Mesh); + + void registerFields(const std::string &MeshName) const; + void unregisterFields() const; + + /// Compute surface insitu temperature from conservative temperature + void computeSurfInsituTemp(const Array3DReal &TracerArray, + const VertCoord *VCoord, + const Eos *EosInst) const; +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index af84424217ff..d64b461896a5 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -255,6 +255,8 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerDiffusion.Enabled = false; TestTendencies->TracerHyperDiff.Enabled = false; TestTendencies->SfcStressForcing.Enabled = false; + TestTendencies->SfcTracerForcing.Enabled = false; + TestTendencies->SfcThicknessForcing.Enabled = false; TestTendencies->SurfaceTracerRestoring.Enabled = false; TestTendencies->ExplicitBottomDrag.Enabled = false; DefVAdv->ThickVertAdvEnabled = false; From 4956cf880de8c379fb078c2a2a1ea677967b4270 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 09:18:03 -0700 Subject: [PATCH 02/27] added the enthalpy of mass fluxes; CtFrz has public interface --- components/omega/src/ocn/Eos.cpp | 12 ++++ components/omega/src/ocn/Eos.h | 6 ++ components/omega/src/ocn/Tendencies.cpp | 14 +++-- components/omega/src/ocn/TendencyTerms.cpp | 7 ++- components/omega/src/ocn/TendencyTerms.h | 66 +++++++++++++--------- 5 files changed, 73 insertions(+), 32 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 0a98eede897d..542a0aa20686 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -346,6 +346,18 @@ 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."); + return 0; +} + /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 462103c2c3f2..ca2ef47b352b 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -1027,6 +1027,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(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 34e91a59053c..c5d079d8e59e 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -488,7 +488,8 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), ExplicitBottomDrag(Mesh, VCoord), SfcThicknessForcing(Mesh, VCoord), - SfcTracerForcing(Mesh, VCoord, Tracers::IndxTemp, Tracers::IndxSalt), + SfcTracerForcing(Mesh, VCoord, Tracers::IndxTemp, Tracers::IndxSalt, + EqState), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), @@ -1005,16 +1006,21 @@ void Tendencies::computeTracerTendenciesOnly( const auto &ShortWaveHeatFlux = ForcingState->TracerForcing.ShortWaveHeatFluxCell; const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &RiverRunoffFlux = + ForcingState->TracerForcing.RiverRunoffFluxCell; const auto &SeaIceSaltFlux = ForcingState->TracerForcing.SeaIceSaltFluxCell; + const auto &PressureMid = VCoord->PressureMid; parallelFor( {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { LocSfcTracerForcing( - LocTracerTend, ICell, LatentHeatFlux, SensibleHeatFlux, - LongWaveHeatFluxUp, LongWaveHeatFluxDown, SeaIceHeatFlux, - ShortWaveHeatFlux, SnowFlux, IceRunoffFlux, SeaIceSaltFlux); + LocTracerTend, ICell, TracerArray, PressureMid, LatentHeatFlux, + SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, + IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); }); Pacer::stop("Tend:sfcTracerForcing", 2); } diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index f2c6bdca73ad..37bfe6ee0500 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -11,6 +11,7 @@ #include "TendencyTerms.h" #include "AuxiliaryState.h" #include "DataTypes.h" +#include "Eos.h" #include "HorzMesh.h" #include "HorzOperators.h" #include "OceanState.h" @@ -78,9 +79,11 @@ SfcThicknessForcingOnCell::SfcThicknessForcingOnCell(const HorzMesh *Mesh, SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord, I4 TempTracerIndex, - I4 SaltTracerIndex) + I4 SaltTracerIndex, + const Eos *EosInst) : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), - MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + EosImpl(VCoord) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 1745ee713734..77b8ce93df04 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -11,6 +11,7 @@ //===----------------------------------------------------------------------===// #include "AuxiliaryState.h" +#include "Eos.h" #include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" @@ -412,18 +413,20 @@ class SfcTracerForcingOnCell { bool Enabled = false; SfcTracerForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord, - I4 TempTracerIndex, I4 SaltTracerIndex); - - KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 ICell, - const Array1DReal &LatentHeatFlux, - const Array1DReal &SensibleHeatFlux, - const Array1DReal &LongWaveHeatFluxUp, - const Array1DReal &LongWaveHeatFluxDown, - const Array1DReal &SeaIceHeatFlux, - const Array1DReal &ShortWaveHeatFlux, - const Array1DReal &SnowFlux, - const Array1DReal &IceRunoffFlux, - const Array1DReal &SeaIceSaltFlux) const { + I4 TempTracerIndex, I4 SaltTracerIndex, + const Eos *EosInst); + + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, + const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, + const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, + const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, + const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { const I4 KTop = MinLayerCell(ICell); if (KTop > MaxLayerCell(ICell)) { @@ -431,31 +434,42 @@ class SfcTracerForcingOnCell { } if (TempIndex >= 0) { - const Real HeatFlux = LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + - LongWaveHeatFluxUp(ICell) + - LongWaveHeatFluxDown(ICell) + - SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell); - // + - // (RainFlux(ICell) + RiverRunoffFlux(ICell)) * - // Cp0Sw * TracerCell(TempIndex, ICell, KTop) + - // (SnowFlux(ICell) + IceRunoffFlux(ICell)) * - // (Cp0Sw * Eos.Ctfreez - LatIce; - - Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFactor; + const Real PTop = PressureMid(ICell, KTop); + const Real SaTop = SaltIndex >= 0 + ? TracerCell(SaltIndex, ICell, KTop) + : 0.0_Real; // not sure we want zero here? + const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTop, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); + + // Heat tendencies are due to direct heat fluxes + enthalpy fluxes + // The enthalpy of liquid water is assumed to be: + // - local SST for liquid mass fluxes (rain, rivers) + // - local freezing point for solid --> liq mass fluxes (snow, frozen + // runoff) + // - solid mass fluxes are locally melted by the ocean (constant Lat + // heat of fusion) + const Real HeatFlux = + LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * + (Cp0Sw * CtFrz - LatIce); + + Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; } if (SaltIndex >= 0) { - Tend(SaltIndex, ICell, KTop) += SeaIceSaltFlux(ICell) * SFluxFactor; + Tend(SaltIndex, ICell, KTop) += SeaIceSaltFlux(ICell) * SFluxFac; } } private: I4 TempIndex; I4 SaltIndex; - Real HFluxFactor; - Real SFluxFactor; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; + Teos10Eos EosImpl; }; // Tracer horizontal advection term From 72f00b3561cbadf865285b75dd4fcbe9cf210cec Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:19:13 -0700 Subject: [PATCH 03/27] added a test for thermo forcing tendencies --- components/omega/test/ocn/TendenciesTest.cpp | 405 +++++++++++++++++++ 1 file changed, 405 insertions(+) diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5268d0a436ea..1045701339dd 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -54,6 +54,9 @@ struct TestSetup { constexpr Geometry Geom = Geometry::Spherical; constexpr int NVertLayers = 60; +int testSfcTracerForcing(); +int testSfcThicknessForcing(); + int initState() { int Err = 0; @@ -305,6 +308,12 @@ int testTendencies() { DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + // Test surface tracer forcing with enthalpy terms + Err += testSfcTracerForcing(); + + // Test surface thickness forcing with freshwater terms + Err += testSfcThicknessForcing(); + // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; int NEdgesOwned = Mesh->NEdgesOwned; @@ -339,6 +348,402 @@ int testTendencies() { return Err; } +int testSfcTracerForcing() { + int Err = 0; + + auto *VCoord = VertCoord::getDefault(); + auto *DefTendencies = Tendencies::getDefault(); + auto *State = OceanState::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + auto *DefForcing = Forcing::getDefault(); + auto *EosInst = Eos::getInstance(); + + Array3DReal TracerArray = Tracers::getAll(0); + + const I4 TempIndex = Tracers::IndxTemp; + const I4 SaltIndex = Tracers::IndxSalt; + + if (TempIndex < 0 || SaltIndex < 0) { + LOG_ERROR("TendenciesTest: Invalid tracer indices for SfcTracerForcing"); + return -1; + } + + deepCopy(DefTendencies->TracerTend, 0._Real); + + // Set up single test cell at top layer + const I4 ICellTest = 0; + const I4 KTop = VCoord->MinLayerCell(ICellTest); + + if (KTop > VCoord->MaxLayerCell(ICellTest)) { + LOG_ERROR("TendenciesTest: Test cell has no layers"); + return -1; + } + + // Known tracer values for testing + const Real CtTopValue = 15.0_Real; // °C (conservative temperature) + const Real SaTopValue = 35.0_Real; // g/kg (salinity) + + // Set tracer values at test cell + OMEGA_SCOPE(LocTracerArray, TracerArray); + Kokkos::parallel_for( + "SetTestTracersForcing", 1, KOKKOS_LAMBDA(int i) { + LocTracerArray(TempIndex, ICellTest, KTop) = CtTopValue; + LocTracerArray(SaltIndex, ICellTest, KTop) = SaTopValue; + }); + + // Retrieve forcing field views + auto &SensibleHeatFlux = DefForcing->TracerForcing.SensibleHeatFluxCell; + auto &LatentHeatFlux = DefForcing->TracerForcing.LatentHeatFluxCell; + auto &LongWaveHeatFluxUp = DefForcing->TracerForcing.LongWaveHeatFluxUpCell; + auto &LongWaveHeatFluxDown = + DefForcing->TracerForcing.LongWaveHeatFluxDownCell; + auto &SeaIceHeatFlux = DefForcing->TracerForcing.SeaIceHeatFluxCell; + auto &ShortWaveHeatFlux = DefForcing->TracerForcing.ShortWaveHeatFluxCell; + auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; + auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; + auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; + auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; + auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; + + // Initialize all fluxes to zero + deepCopy(SensibleHeatFlux, 0._Real); + deepCopy(LatentHeatFlux, 0._Real); + deepCopy(LongWaveHeatFluxUp, 0._Real); + deepCopy(LongWaveHeatFluxDown, 0._Real); + deepCopy(SeaIceHeatFlux, 0._Real); + deepCopy(ShortWaveHeatFlux, 0._Real); + deepCopy(RainFlux, 0._Real); + deepCopy(RiverRunoffFlux, 0._Real); + deepCopy(SnowFlux, 0._Real); + deepCopy(IceRunoffFlux, 0._Real); + deepCopy(SeaIceSaltFlux, 0._Real); + + // Set test forcing values + // Non-zero sensible heat: 100 W/m² + const Real TestSensibleHeat = 100.0_Real; + // Non-zero rain: 1e-8 kg/m²/s + const Real TestRain = 1.0e-8_Real; + // Non-zero snow: 5e-9 kg/m²/s + const Real TestSnow = 5.0e-9_Real; + // Sea ice salt flux: 1e-4 kg/m²/s + const Real TestSeaIceSaltFlux = 1.0e-4_Real; + + OMEGA_SCOPE(LocSensibleHeatFlux, SensibleHeatFlux); + OMEGA_SCOPE(LocRainFlux, RainFlux); + OMEGA_SCOPE(LocSnowFlux, SnowFlux); + OMEGA_SCOPE(LocSeaIceSaltFlux, SeaIceSaltFlux); + Kokkos::parallel_for( + "SetTestForcingTracer", 1, KOKKOS_LAMBDA(int i) { + LocSensibleHeatFlux(ICellTest) = TestSensibleHeat; + LocRainFlux(ICellTest) = TestRain; + LocSnowFlux(ICellTest) = TestSnow; + LocSeaIceSaltFlux(ICellTest) = TestSeaIceSaltFlux; + }); + + DefForcing->computeAll(); + + // Disable all tendencies except SfcTracerForcing + const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; + const bool OrigSfcThicknessEnabled = + DefTendencies->SfcThicknessForcing.Enabled; + const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; + const bool OrigPseudoThicknessDiv = + DefTendencies->PseudoThicknessFluxDiv.Enabled; + const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; + const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; + const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; + const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; + const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; + const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; + const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; + const bool OrigSurfaceTracerRestoring = + DefTendencies->SurfaceTracerRestoring.Enabled; + + DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->SfcThicknessForcing.Enabled = false; + DefTendencies->SfcTracerForcing.Enabled = false; + DefTendencies->PseudoThicknessFluxDiv.Enabled = false; + DefTendencies->PotentialVortHAdv.Enabled = false; + DefTendencies->KEGrad.Enabled = false; + DefTendencies->VelocityDiffusion.Enabled = false; + DefTendencies->VelocityHyperDiff.Enabled = false; + DefTendencies->TracerHorzAdv.Enabled = false; + DefTendencies->TracerDiffusion.Enabled = false; + DefTendencies->TracerHyperDiff.Enabled = false; + DefTendencies->SurfaceTracerRestoring.Enabled = false; + + // Compute tendencies + int ThickTimeLevel = 0; + int VelTimeLevel = 0; + int TracerTimeLevel = 0; + TimeInstant Time; + TimeInterval Interval(1., TimeUnits::Seconds); + + // because vertical advection tendencies are always on, we need to compute a + // baseline first. the actual test is whether the total tendencies change + // with the flag toggling. + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + HostArray3DReal TracerTendBaseH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendBaseH, DefTendencies->TracerTend); + const Real BaselineTempTend = TracerTendBaseH(TempIndex, ICellTest, KTop); + const Real BaselineSaltTend = TracerTendBaseH(SaltIndex, ICellTest, KTop); + // Now enable SfcTracerForcing and compute again + DefTendencies->SfcTracerForcing.Enabled = true; + + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + // Build two reference expectations for temperature tendency: + // 1) fixed estimate (expected to fail under strict tolerance), + // 2) TEOS-10 freezing CT (expected to pass under strict tolerance). + const Real CtFrzEstimate = -2.0_Real; + const Real ExpectedTempTendEstimate = + (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + + TestSnow * (Cp0Sw * CtFrzEstimate - LatIce)) * + HFluxFac; + + HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); + deepCopy(PressureMidH, VCoord->PressureMid); + const Real PTop = PressureMidH(ICellTest, KTop); + const Real CtFrzTeos = EosInst->calcCtFreezing(SaTopValue, PTop, 0.0_Real); + const Real ExpectedTempTendTeos = + (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + + TestSnow * (Cp0Sw * CtFrzTeos - LatIce)) * + HFluxFac; + + // SaltTend = SeaIceSaltFlux * SFluxFac + const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; + + HostArray3DReal TracerTendH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendH, DefTendencies->TracerTend); + const Real ComputedTempTend = + TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; + const Real ComputedSaltTend = + TracerTendH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 + + // Expected-fail check with fixed CtFrz estimate. + if (!isApprox(ComputedTempTend, ExpectedTempTendEstimate, RelTol, AbsTol)) { + LOG_INFO( + "TendenciesTest: expected tempTend fail because CtFrzEstimate != EOS " + "CtFrz - PASS"); + LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", + ExpectedTempTendEstimate, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTendEstimate)); + } else { + Err++; + LOG_ERROR("TendenciesTest: CtFrz estimate unexpectedly matched strict " + "reference - FAIL"); + } + + // Expected-pass check with TEOS freezing CT reference. + if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); + LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", + ExpectedTempTendTeos, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTendTeos)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); + } + + // Check salinity tendency + if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); + LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, + ComputedSaltTend, + Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); + } + + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; + DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; + DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; + DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; + DefTendencies->KEGrad.Enabled = OrigKEGrad; + DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; + DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; + DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; + DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; + DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; + DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; + + return Err; +} + +int testSfcThicknessForcing() { + int Err = 0; + + auto *VCoord = VertCoord::getDefault(); + auto *DefTendencies = Tendencies::getDefault(); + auto *State = OceanState::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + auto *DefForcing = Forcing::getDefault(); + + Array3DReal TracerArray = Tracers::getAll(0); + + deepCopy(DefTendencies->PseudoThicknessTend, 0._Real); + + // Set up single test cell at top layer + const I4 ICellTest = 0; + const I4 KTop = VCoord->MinLayerCell(ICellTest); + + if (KTop > VCoord->MaxLayerCell(ICellTest)) { + LOG_ERROR("TendenciesTest: Test cell has no layers for thickness test"); + return -1; + } + + // Retrieve forcing field views for thickness + auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; + auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; + auto &EvaporationFlux = DefForcing->TracerForcing.EvaporationFluxCell; + auto &SeaIceFreshWater = DefForcing->TracerForcing.SeaIceFreshWaterFluxCell; + auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; + auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; + auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; + + // Initialize all fluxes to zero + deepCopy(SnowFlux, 0._Real); + deepCopy(RainFlux, 0._Real); + deepCopy(EvaporationFlux, 0._Real); + deepCopy(SeaIceFreshWater, 0._Real); + deepCopy(IceRunoffFlux, 0._Real); + deepCopy(RiverRunoffFlux, 0._Real); + deepCopy(SeaIceSaltFlux, 0._Real); + + // Set test freshwater flux values + // Rain: 1e-8 kg/m²/s + const Real TestRain = 1.0e-8_Real; + // Snow: 5e-9 kg/m²/s + const Real TestSnow = 5.0e-9_Real; + // Ice runoff: 2e-9 kg/m²/s + const Real TestIceRunoff = 2.0e-9_Real; + // River runoff: 3e-9 kg/m²/s + const Real TestRiverRunoff = 3.0e-9_Real; + // Sea ice freshwater: 1e-9 kg/m²/s + const Real TestSeaIceFreshWater = 1.0e-9_Real; + // Sea ice salt flux: 1e-4 kg/m²/s (affects thickness via salt) + const Real TestSeaIceSaltFlux = 1.0e-4_Real; + + OMEGA_SCOPE(LocSnowFlux, SnowFlux); + OMEGA_SCOPE(LocRainFlux, RainFlux); + OMEGA_SCOPE(LocIceRunoffFlux, IceRunoffFlux); + OMEGA_SCOPE(LocRiverRunoffFlux, RiverRunoffFlux); + OMEGA_SCOPE(LocSeaIceFreshWater, SeaIceFreshWater); + OMEGA_SCOPE(LocSeaIceSaltFlux, SeaIceSaltFlux); + Kokkos::parallel_for( + "SetTestForcingThickness", 1, KOKKOS_LAMBDA(int i) { + LocRainFlux(ICellTest) = TestRain; + LocSnowFlux(ICellTest) = TestSnow; + LocIceRunoffFlux(ICellTest) = TestIceRunoff; + LocRiverRunoffFlux(ICellTest) = TestRiverRunoff; + LocSeaIceFreshWater(ICellTest) = TestSeaIceFreshWater; + LocSeaIceSaltFlux(ICellTest) = TestSeaIceSaltFlux; + }); + + DefForcing->computeAll(); + + const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; + const bool OrigSfcThicknessEnabled = + DefTendencies->SfcThicknessForcing.Enabled; + const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; + const bool OrigPseudoThicknessDiv = + DefTendencies->PseudoThicknessFluxDiv.Enabled; + const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; + const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; + const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; + const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; + const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; + const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; + const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; + const bool OrigSurfaceTracerRestoring = + DefTendencies->SurfaceTracerRestoring.Enabled; + + DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->SfcThicknessForcing.Enabled = false; + DefTendencies->SfcTracerForcing.Enabled = false; + DefTendencies->PseudoThicknessFluxDiv.Enabled = false; + DefTendencies->PotentialVortHAdv.Enabled = false; + DefTendencies->KEGrad.Enabled = false; + DefTendencies->VelocityDiffusion.Enabled = false; + DefTendencies->VelocityHyperDiff.Enabled = false; + DefTendencies->TracerHorzAdv.Enabled = false; + DefTendencies->TracerDiffusion.Enabled = false; + DefTendencies->TracerHyperDiff.Enabled = false; + DefTendencies->SurfaceTracerRestoring.Enabled = false; + + // Compute baseline tendencies (vertical advection is always on) + int ThickTimeLevel = 0; + int VelTimeLevel = 0; + TimeInstant Time; + DefTendencies->computePseudoThicknessTendenciesOnly( + State, AuxState, ThickTimeLevel, VelTimeLevel, Time); + + HostArray2DReal PseudoThicknessTendBaseH = + createHostMirrorCopy(DefTendencies->PseudoThicknessTend); + deepCopy(PseudoThicknessTendBaseH, DefTendencies->PseudoThicknessTend); + const Real BaselineThickTend = PseudoThicknessTendBaseH(ICellTest, KTop); + + // Now enable SfcThicknessForcing and compute again + DefTendencies->SfcThicknessForcing.Enabled = true; + DefTendencies->computePseudoThicknessTendenciesOnly( + State, AuxState, ThickTimeLevel, VelTimeLevel, Time); + + // Calculate expected thickness tendency + // ThickTend = (Rain + Snow + IceRunoff + RiverRunoff + SeaIceFreshWater + + // SeaIceSaltFlux) / RhoSw + const Real ExpectedThickTend = + (TestRain + TestSnow + TestIceRunoff + TestRiverRunoff + + TestSeaIceFreshWater + TestSeaIceSaltFlux) / + RhoSw; + + HostArray2DReal PseudoThicknessTendH = + createHostMirrorCopy(DefTendencies->PseudoThicknessTend); + deepCopy(PseudoThicknessTendH, DefTendencies->PseudoThicknessTend); + const Real ComputedThickTend = + PseudoThicknessTendH(ICellTest, KTop) - BaselineThickTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; + + // Check thickness tendency + if (!isApprox(ComputedThickTend, ExpectedThickTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcThicknessForcing thickness tendency FAIL"); + LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedThickTend, + ComputedThickTend, + Kokkos::abs(ComputedThickTend - ExpectedThickTend)); + } else { + LOG_INFO("TendenciesTest: SfcThicknessForcing thickness tendency PASS"); + } + + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; + DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; + DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; + DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; + DefTendencies->KEGrad.Enabled = OrigKEGrad; + DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; + DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; + DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; + DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; + DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; + DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; + + return Err; +} + void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); From 45c792f2571a97976703bf0b742476eccc1f2b1a Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:46:07 -0700 Subject: [PATCH 04/27] made mass enthalpy flux dependent on thickness flag - under discussion --- components/omega/src/ocn/Tendencies.cpp | 14 +- components/omega/src/ocn/TendencyTerms.h | 46 ++++--- components/omega/test/ocn/TendenciesTest.cpp | 133 +++++++++++++++---- 3 files changed, 141 insertions(+), 52 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index c5d079d8e59e..a24334c2323a 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -1012,15 +1012,17 @@ void Tendencies::computeTracerTendenciesOnly( ForcingState->TracerForcing.RiverRunoffFluxCell; const auto &SeaIceSaltFlux = ForcingState->TracerForcing.SeaIceSaltFluxCell; - const auto &PressureMid = VCoord->PressureMid; + const auto &PressureMid = VCoord->PressureMid; + const bool UseMassFluxHeat = SfcThicknessForcing.Enabled; parallelFor( {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { - LocSfcTracerForcing( - LocTracerTend, ICell, TracerArray, PressureMid, LatentHeatFlux, - SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, - SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, - IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); + LocSfcTracerForcing(LocTracerTend, ICell, TracerArray, PressureMid, + LatentHeatFlux, SensibleHeatFlux, + LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, + RainFlux, IceRunoffFlux, RiverRunoffFlux, + SeaIceSaltFlux, UseMassFluxHeat); }); Pacer::stop("Tend:sfcTracerForcing", 2); } diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 77b8ce93df04..30ceafde45d0 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -416,17 +416,16 @@ class SfcTracerForcingOnCell { I4 TempTracerIndex, I4 SaltTracerIndex, const Eos *EosInst); - KOKKOS_FUNCTION void - operator()(const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, - const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, - const Array1DReal &SensibleHeatFlux, - const Array1DReal &LongWaveHeatFluxUp, - const Array1DReal &LongWaveHeatFluxDown, - const Array1DReal &SeaIceHeatFlux, - const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, - const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, - const Array1DReal &RiverRunoffFlux, - const Array1DReal &SeaIceSaltFlux) const { + KOKKOS_FUNCTION void operator()( + const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, + const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, const Array1DReal &ShortWaveHeatFlux, + const Array1DReal &SnowFlux, const Array1DReal &RainFlux, + const Array1DReal &IceRunoffFlux, const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux, const bool UseMassFluxHeat) const { const I4 KTop = MinLayerCell(ICell); if (KTop > MaxLayerCell(ICell)) { @@ -441,20 +440,29 @@ class SfcTracerForcingOnCell { const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTop, 0.0_Real); const Real CtTop = TracerCell(TempIndex, ICell, KTop); - // Heat tendencies are due to direct heat fluxes + enthalpy fluxes - // The enthalpy of liquid water is assumed to be: + // Always include direct surface heat fluxes. + const Real DirectHeatFlux = + LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell); + + // Apply enthalpy of mass fluxes only when thickness forcing is + // enabled. + const Real MassFluxHeat = + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * + (Cp0Sw * CtFrz - LatIce); + // Note: the enthalpy of liquid water above is assumed to be: // - local SST for liquid mass fluxes (rain, rivers) // - local freezing point for solid --> liq mass fluxes (snow, frozen // runoff) // - solid mass fluxes are locally melted by the ocean (constant Lat // heat of fusion) + // - meltwater enthalpy from sea ice is already included in + // SeaIceHeatFlux + const Real HeatFlux = - LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + - LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + - SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + - (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + - (SnowFlux(ICell) + IceRunoffFlux(ICell)) * - (Cp0Sw * CtFrz - LatIce); + DirectHeatFlux + (UseMassFluxHeat ? MassFluxHeat : 0.0_Real); Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; } diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 1045701339dd..d30babd08797 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -494,11 +494,71 @@ int testSfcTracerForcing() { // Now enable SfcTracerForcing and compute again DefTendencies->SfcTracerForcing.Enabled = true; + // First pass: thickness forcing disabled, so only direct heat flux should + // contribute to temperature tendency. + DefTendencies->SfcThicknessForcing.Enabled = false; DefTendencies->computeAllTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, TracerTimeLevel, Time, Interval); - // Build two reference expectations for temperature tendency: + HostArray3DReal TracerTendNoMassH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendNoMassH, DefTendencies->TracerTend); + const Real ComputedTempTendNoMass = + TracerTendNoMassH(TempIndex, ICellTest, KTop) - BaselineTempTend; + const Real ComputedSaltTendNoMass = + TracerTendNoMassH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; + + // With thickness forcing disabled, only direct heat flux terms are applied. + const Real ExpectedTempTendNoMass = TestSensibleHeat * HFluxFac; + + // SaltTend = SeaIceSaltFlux * SFluxFac + const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 + + if (!isApprox(ComputedTempTendNoMass, ExpectedTempTendNoMass, RelTol, + AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL with " + "SfcThicknessForcing disabled"); + LOG_ERROR(" Expected (direct only): {}, Computed: {}, Diff: {}", + ExpectedTempTendNoMass, ComputedTempTendNoMass, + Kokkos::abs(ComputedTempTendNoMass - ExpectedTempTendNoMass)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS with " + "SfcThicknessForcing disabled"); + } + + if (!isApprox(ComputedSaltTendNoMass, ExpectedSaltTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL with " + "SfcThicknessForcing disabled"); + LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, + ComputedSaltTendNoMass, + Kokkos::abs(ComputedSaltTendNoMass - ExpectedSaltTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS with " + "SfcThicknessForcing disabled"); + } + + // Second pass: thickness forcing enabled, so mass-flux enthalpy terms are + // also included in temperature tendency. + DefTendencies->SfcThicknessForcing.Enabled = true; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + HostArray3DReal TracerTendMassH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendMassH, DefTendencies->TracerTend); + const Real ComputedTempTendMass = + TracerTendMassH(TempIndex, ICellTest, KTop) - BaselineTempTend; + const Real ComputedSaltTendMass = + TracerTendMassH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; + + // Build two reference expectations for the mass-on case: // 1) fixed estimate (expected to fail under strict tolerance), // 2) TEOS-10 freezing CT (expected to pass under strict tolerance). const Real CtFrzEstimate = -2.0_Real; @@ -516,28 +576,31 @@ int testSfcTracerForcing() { TestSnow * (Cp0Sw * CtFrzTeos - LatIce)) * HFluxFac; - // SaltTend = SeaIceSaltFlux * SFluxFac - const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; - - HostArray3DReal TracerTendH = - createHostMirrorCopy(DefTendencies->TracerTend); - deepCopy(TracerTendH, DefTendencies->TracerTend); - const Real ComputedTempTend = - TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; - const Real ComputedSaltTend = - TracerTendH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; - - constexpr Real RelTol = 1.0e-10_Real; - constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 + // Expected-fail check: no-mass expectation should fail when mass-flux + // terms are enabled. + if (!isApprox(ComputedTempTendMass, ExpectedTempTendNoMass, RelTol, + AbsTol)) { + LOG_INFO( + "TendenciesTest: expected tempTend fail because mass-flux heat is " + "enabled but compared against direct-only reference - PASS"); + LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", + ExpectedTempTendNoMass, ComputedTempTendMass, + Kokkos::abs(ComputedTempTendMass - ExpectedTempTendNoMass)); + } else { + Err++; + LOG_ERROR("TendenciesTest: mass-flux-enabled run unexpectedly matched " + "direct-only reference - FAIL"); + } // Expected-fail check with fixed CtFrz estimate. - if (!isApprox(ComputedTempTend, ExpectedTempTendEstimate, RelTol, AbsTol)) { + if (!isApprox(ComputedTempTendMass, ExpectedTempTendEstimate, RelTol, + AbsTol)) { LOG_INFO( "TendenciesTest: expected tempTend fail because CtFrzEstimate != EOS " "CtFrz - PASS"); LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendEstimate, ComputedTempTend, - Kokkos::abs(ComputedTempTend - ExpectedTempTendEstimate)); + ExpectedTempTendEstimate, ComputedTempTendMass, + Kokkos::abs(ComputedTempTendMass - ExpectedTempTendEstimate)); } else { Err++; LOG_ERROR("TendenciesTest: CtFrz estimate unexpectedly matched strict " @@ -545,25 +608,41 @@ int testSfcTracerForcing() { } // Expected-pass check with TEOS freezing CT reference. - if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { + if (!isApprox(ComputedTempTendMass, ExpectedTempTendTeos, RelTol, AbsTol)) { Err++; LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendTeos, ComputedTempTend, - Kokkos::abs(ComputedTempTend - ExpectedTempTendTeos)); + ExpectedTempTendTeos, ComputedTempTendMass, + Kokkos::abs(ComputedTempTendMass - ExpectedTempTendTeos)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASSwith " + "SfcThicknessForcing enabled"); } - // Check salinity tendency - if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { + // Check salinity tendency for mass-on pass + if (!isApprox(ComputedSaltTendMass, ExpectedSaltTend, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL with " + "SfcThicknessForcing enabled"); LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, - ComputedSaltTend, - Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); + ComputedSaltTendMass, + Kokkos::abs(ComputedSaltTendMass - ExpectedSaltTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS with " + "SfcThicknessForcing enabled"); + } + + if (!isApprox(ComputedSaltTendNoMass, ComputedSaltTendMass, RelTol, + AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency changed with " + "SfcThicknessForcing toggle - FAIL"); + LOG_ERROR(" Off: {}, On: {}, Diff: {}", ComputedSaltTendNoMass, + ComputedSaltTendMass, + Kokkos::abs(ComputedSaltTendNoMass - ComputedSaltTendMass)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency invariant under " + "SfcThicknessForcing toggle PASS"); } DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; From 36c9a37d1a41dc973b516f96ee68366a7c9d8035 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:27:46 -0700 Subject: [PATCH 05/27] draft a non-teos10 CtFrz in comments - WIP --- components/omega/src/ocn/Eos.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 542a0aa20686..b2d7905e6789 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -355,7 +355,9 @@ Real Eos::calcCtFreezing(const Real Sa, const Real P, ABORT_ERROR("Eos::calcCtFreezing: CT freezing temperature is only " "implemented for TEOS-10. Support for the current EOS " "choice has not yet been developed."); - return 0; + // 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 From 9b8776978bfea62584a165a3037e0e4621bca425 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:49:55 -0700 Subject: [PATCH 06/27] updated the documentation --- components/omega/doc/devGuide/Forcing.md | 19 +++++++++++++----- .../omega/doc/devGuide/TendencyTerms.md | 6 ++++-- components/omega/doc/userGuide/Forcing.md | 20 ++++++++++++++++--- .../omega/doc/userGuide/TendencyTerms.md | 2 +- 4 files changed, 36 insertions(+), 11 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 18e9c2072190..700711e86761 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -58,10 +58,11 @@ the surface layer pseudo-thickness. - `LatentHeatFlux`, `SensibleHeatFlux` - `LongWaveHeatFluxUp`, `LongWaveHeatFluxDown` - `SeaIceHeatFlux`, `ShortWaveHeatFlux` - - `SeaIceSaltFlux`, `SnowFlux`, `IceRunoffFlux` + - 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 salinity (g/kg) in the surface layer. + and applies the external sea-ice salt flux to salinity (g/kg) in the surface layer. [under discussion: in the latest implementation, if the thickness tendencies are turned off, the temperature tendency does not include the enthalpy associated with explicit mass fluxes] ### Surface flux forcing key classes/components @@ -73,14 +74,20 @@ the surface layer pseudo-thickness. - Computes freshwater flux contribution: $\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: computes the sum of the six heat-flux fields and scales it by $H_{\text{FluxFac}}$ + - For temperature: computes + $Q_{\text{direct}} = Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$ + and scales by $H_{\text{FluxFac}}$. + - For temperature: when `SfcThicknessForcing` is enabled, also adds + mass-flux enthalpy + $(\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. - 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 `computeThicknessTendenciesOnly` + - Calls `SfcThicknessForcingOnCell` in `computePseudoThicknessTendenciesOnly` - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring ### Surface flux forcing config coupling @@ -88,9 +95,11 @@ the surface layer pseudo-thickness. - `Omega.Tendencies.SfcThicknessForcingTendencyEnable` - gates execution of coupled flux thickness kernel - controls freshwater and salt flux forcing on sea surface height + - also gates whether mass-flux enthalpy terms are added in tracer + temperature forcing - `Omega.Tendencies.SfcTracerForcingTendencyEnable` - gates execution of coupled flux tracer kernel - - controls heat flux forcing on temperature and salt flux forcing on salinity + - controls direct heat flux forcing on temperature and salt flux forcing on salinity ## Surface tracer restoring design diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 5fa72197132f..028ba9c02032 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/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). diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index b49eabd0586f..34a964a61510 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -60,6 +60,11 @@ Omega: - `Tendencies.SfcThicknessForcingTendencyEnable`: enables coupled freshwater and salt flux forcing on thickness - `Tendencies.SfcTracerForcingTendencyEnable`: enables coupled heat and salt flux forcing on tracers +When `Tendencies.SfcTracerForcingTendencyEnable` is enabled, direct surface heat +flux terms are always applied to temperature. Additional mass-flux enthalpy +terms (rain/river and snow/ice runoff) are applied only when +`Tendencies.SfcThicknessForcingTendencyEnable` is also enabled. + ### Required input fields Coupled flux forcing uses 13 auxiliary fields organized by type: @@ -93,10 +98,19 @@ by the equivalent `ocn_comp_mct.F`. - 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 the sum of the six heat-flux fields, - converted to conservative-temperature tendency via +- Temperature tendency is computed from direct heat flux plus optional + 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. [soon to be updated with latent heat and enthalpy of liquid water] + specific heat of seawater defined by TEOS-10. + The direct heat part is + $Q_{\text{direct}} = Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$. + The mass-flux enthalpy part is + $Q_{\text{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 computed from EOS at top-layer salinity and pressure. + The applied heat flux is + $Q_{\text{direct}} + Q_{\text{mass}}$ when + `Tendencies.SfcThicknessForcingTendencyEnable` is true, and + $Q_{\text{direct}}$ otherwise. - 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). diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 1259d3387054..7847f0cfb357 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -21,7 +21,7 @@ tendency terms are currently implemented: | 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, 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 From fae5a4a8fe6ff1313dc71c749ce7ad529b9b5573 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:07:40 -0700 Subject: [PATCH 07/27] Revert "made mass enthalpy flux dependent on thickness flag - under discussion" This reverts commit 70b0ca28a7939588d536496e9084c38b0663c5e1. --- components/omega/src/ocn/Tendencies.cpp | 14 +- components/omega/src/ocn/TendencyTerms.h | 46 +++---- components/omega/test/ocn/TendenciesTest.cpp | 133 ++++--------------- 3 files changed, 52 insertions(+), 141 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index a24334c2323a..c5d079d8e59e 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -1012,17 +1012,15 @@ void Tendencies::computeTracerTendenciesOnly( ForcingState->TracerForcing.RiverRunoffFluxCell; const auto &SeaIceSaltFlux = ForcingState->TracerForcing.SeaIceSaltFluxCell; - const auto &PressureMid = VCoord->PressureMid; - const bool UseMassFluxHeat = SfcThicknessForcing.Enabled; + const auto &PressureMid = VCoord->PressureMid; parallelFor( {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { - LocSfcTracerForcing(LocTracerTend, ICell, TracerArray, PressureMid, - LatentHeatFlux, SensibleHeatFlux, - LongWaveHeatFluxUp, LongWaveHeatFluxDown, - SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, - RainFlux, IceRunoffFlux, RiverRunoffFlux, - SeaIceSaltFlux, UseMassFluxHeat); + LocSfcTracerForcing( + LocTracerTend, ICell, TracerArray, PressureMid, LatentHeatFlux, + SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, + IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); }); Pacer::stop("Tend:sfcTracerForcing", 2); } diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 30ceafde45d0..77b8ce93df04 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -416,16 +416,17 @@ class SfcTracerForcingOnCell { I4 TempTracerIndex, I4 SaltTracerIndex, const Eos *EosInst); - KOKKOS_FUNCTION void operator()( - const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, - const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, - const Array1DReal &SensibleHeatFlux, - const Array1DReal &LongWaveHeatFluxUp, - const Array1DReal &LongWaveHeatFluxDown, - const Array1DReal &SeaIceHeatFlux, const Array1DReal &ShortWaveHeatFlux, - const Array1DReal &SnowFlux, const Array1DReal &RainFlux, - const Array1DReal &IceRunoffFlux, const Array1DReal &RiverRunoffFlux, - const Array1DReal &SeaIceSaltFlux, const bool UseMassFluxHeat) const { + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, + const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, + const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, + const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, + const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { const I4 KTop = MinLayerCell(ICell); if (KTop > MaxLayerCell(ICell)) { @@ -440,29 +441,20 @@ class SfcTracerForcingOnCell { const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTop, 0.0_Real); const Real CtTop = TracerCell(TempIndex, ICell, KTop); - // Always include direct surface heat fluxes. - const Real DirectHeatFlux = - LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + - LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + - SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell); - - // Apply enthalpy of mass fluxes only when thickness forcing is - // enabled. - const Real MassFluxHeat = - (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + - (SnowFlux(ICell) + IceRunoffFlux(ICell)) * - (Cp0Sw * CtFrz - LatIce); - // Note: the enthalpy of liquid water above is assumed to be: + // Heat tendencies are due to direct heat fluxes + enthalpy fluxes + // The enthalpy of liquid water is assumed to be: // - local SST for liquid mass fluxes (rain, rivers) // - local freezing point for solid --> liq mass fluxes (snow, frozen // runoff) // - solid mass fluxes are locally melted by the ocean (constant Lat // heat of fusion) - // - meltwater enthalpy from sea ice is already included in - // SeaIceHeatFlux - const Real HeatFlux = - DirectHeatFlux + (UseMassFluxHeat ? MassFluxHeat : 0.0_Real); + LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * + (Cp0Sw * CtFrz - LatIce); Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; } diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index d30babd08797..1045701339dd 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -494,71 +494,11 @@ int testSfcTracerForcing() { // Now enable SfcTracerForcing and compute again DefTendencies->SfcTracerForcing.Enabled = true; - // First pass: thickness forcing disabled, so only direct heat flux should - // contribute to temperature tendency. - DefTendencies->SfcThicknessForcing.Enabled = false; DefTendencies->computeAllTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, TracerTimeLevel, Time, Interval); - HostArray3DReal TracerTendNoMassH = - createHostMirrorCopy(DefTendencies->TracerTend); - deepCopy(TracerTendNoMassH, DefTendencies->TracerTend); - const Real ComputedTempTendNoMass = - TracerTendNoMassH(TempIndex, ICellTest, KTop) - BaselineTempTend; - const Real ComputedSaltTendNoMass = - TracerTendNoMassH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; - - // With thickness forcing disabled, only direct heat flux terms are applied. - const Real ExpectedTempTendNoMass = TestSensibleHeat * HFluxFac; - - // SaltTend = SeaIceSaltFlux * SFluxFac - const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; - - constexpr Real RelTol = 1.0e-10_Real; - constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 - - if (!isApprox(ComputedTempTendNoMass, ExpectedTempTendNoMass, RelTol, - AbsTol)) { - Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL with " - "SfcThicknessForcing disabled"); - LOG_ERROR(" Expected (direct only): {}, Computed: {}, Diff: {}", - ExpectedTempTendNoMass, ComputedTempTendNoMass, - Kokkos::abs(ComputedTempTendNoMass - ExpectedTempTendNoMass)); - } else { - LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS with " - "SfcThicknessForcing disabled"); - } - - if (!isApprox(ComputedSaltTendNoMass, ExpectedSaltTend, RelTol, AbsTol)) { - Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL with " - "SfcThicknessForcing disabled"); - LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, - ComputedSaltTendNoMass, - Kokkos::abs(ComputedSaltTendNoMass - ExpectedSaltTend)); - } else { - LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS with " - "SfcThicknessForcing disabled"); - } - - // Second pass: thickness forcing enabled, so mass-flux enthalpy terms are - // also included in temperature tendency. - DefTendencies->SfcThicknessForcing.Enabled = true; - DefTendencies->computeAllTendencies(State, AuxState, TracerArray, - ThickTimeLevel, VelTimeLevel, - TracerTimeLevel, Time, Interval); - - HostArray3DReal TracerTendMassH = - createHostMirrorCopy(DefTendencies->TracerTend); - deepCopy(TracerTendMassH, DefTendencies->TracerTend); - const Real ComputedTempTendMass = - TracerTendMassH(TempIndex, ICellTest, KTop) - BaselineTempTend; - const Real ComputedSaltTendMass = - TracerTendMassH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; - - // Build two reference expectations for the mass-on case: + // Build two reference expectations for temperature tendency: // 1) fixed estimate (expected to fail under strict tolerance), // 2) TEOS-10 freezing CT (expected to pass under strict tolerance). const Real CtFrzEstimate = -2.0_Real; @@ -576,31 +516,28 @@ int testSfcTracerForcing() { TestSnow * (Cp0Sw * CtFrzTeos - LatIce)) * HFluxFac; - // Expected-fail check: no-mass expectation should fail when mass-flux - // terms are enabled. - if (!isApprox(ComputedTempTendMass, ExpectedTempTendNoMass, RelTol, - AbsTol)) { - LOG_INFO( - "TendenciesTest: expected tempTend fail because mass-flux heat is " - "enabled but compared against direct-only reference - PASS"); - LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendNoMass, ComputedTempTendMass, - Kokkos::abs(ComputedTempTendMass - ExpectedTempTendNoMass)); - } else { - Err++; - LOG_ERROR("TendenciesTest: mass-flux-enabled run unexpectedly matched " - "direct-only reference - FAIL"); - } + // SaltTend = SeaIceSaltFlux * SFluxFac + const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; + + HostArray3DReal TracerTendH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendH, DefTendencies->TracerTend); + const Real ComputedTempTend = + TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; + const Real ComputedSaltTend = + TracerTendH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 // Expected-fail check with fixed CtFrz estimate. - if (!isApprox(ComputedTempTendMass, ExpectedTempTendEstimate, RelTol, - AbsTol)) { + if (!isApprox(ComputedTempTend, ExpectedTempTendEstimate, RelTol, AbsTol)) { LOG_INFO( "TendenciesTest: expected tempTend fail because CtFrzEstimate != EOS " "CtFrz - PASS"); LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendEstimate, ComputedTempTendMass, - Kokkos::abs(ComputedTempTendMass - ExpectedTempTendEstimate)); + ExpectedTempTendEstimate, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTendEstimate)); } else { Err++; LOG_ERROR("TendenciesTest: CtFrz estimate unexpectedly matched strict " @@ -608,41 +545,25 @@ int testSfcTracerForcing() { } // Expected-pass check with TEOS freezing CT reference. - if (!isApprox(ComputedTempTendMass, ExpectedTempTendTeos, RelTol, AbsTol)) { + if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { Err++; LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendTeos, ComputedTempTendMass, - Kokkos::abs(ComputedTempTendMass - ExpectedTempTendTeos)); + ExpectedTempTendTeos, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTendTeos)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASSwith " - "SfcThicknessForcing enabled"); + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); } - // Check salinity tendency for mass-on pass - if (!isApprox(ComputedSaltTendMass, ExpectedSaltTend, RelTol, AbsTol)) { + // Check salinity tendency + if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL with " - "SfcThicknessForcing enabled"); + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, - ComputedSaltTendMass, - Kokkos::abs(ComputedSaltTendMass - ExpectedSaltTend)); - } else { - LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS with " - "SfcThicknessForcing enabled"); - } - - if (!isApprox(ComputedSaltTendNoMass, ComputedSaltTendMass, RelTol, - AbsTol)) { - Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency changed with " - "SfcThicknessForcing toggle - FAIL"); - LOG_ERROR(" Off: {}, On: {}, Diff: {}", ComputedSaltTendNoMass, - ComputedSaltTendMass, - Kokkos::abs(ComputedSaltTendNoMass - ComputedSaltTendMass)); + ComputedSaltTend, + Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency invariant under " - "SfcThicknessForcing toggle PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); } DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; From 90f045ecb0169b562c42f918b648fc5596f67e32 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:53:52 -0700 Subject: [PATCH 08/27] updated the documentation --- components/omega/doc/devGuide/Forcing.md | 14 +++++--------- components/omega/doc/userGuide/Forcing.md | 17 ++--------------- 2 files changed, 7 insertions(+), 24 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 700711e86761..869019e2469e 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -62,7 +62,7 @@ the surface layer pseudo-thickness. - `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 salinity (g/kg) in the surface layer. [under discussion: in the latest implementation, if the thickness tendencies are turned off, the temperature tendency does not include the enthalpy associated with explicit mass fluxes] + and applies the external sea-ice salt flux to salinity (g/kg) in the surface layer. ### Surface flux forcing key classes/components @@ -74,13 +74,11 @@ the surface layer pseudo-thickness. - Computes freshwater flux contribution: $\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: computes - $Q_{\text{direct}} = Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$ + - 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 temperature: when `SfcThicknessForcing` is enabled, also adds - mass-flux enthalpy - $(\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. - 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 @@ -95,8 +93,6 @@ the surface layer pseudo-thickness. - `Omega.Tendencies.SfcThicknessForcingTendencyEnable` - gates execution of coupled flux thickness kernel - controls freshwater and salt flux forcing on sea surface height - - also gates whether mass-flux enthalpy terms are added in tracer - temperature forcing - `Omega.Tendencies.SfcTracerForcingTendencyEnable` - gates execution of coupled flux tracer kernel - controls direct heat flux forcing on temperature and salt flux forcing on salinity diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index 34a964a61510..e08754ee2651 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -60,10 +60,6 @@ Omega: - `Tendencies.SfcThicknessForcingTendencyEnable`: enables coupled freshwater and salt flux forcing on thickness - `Tendencies.SfcTracerForcingTendencyEnable`: enables coupled heat and salt flux forcing on tracers -When `Tendencies.SfcTracerForcingTendencyEnable` is enabled, direct surface heat -flux terms are always applied to temperature. Additional mass-flux enthalpy -terms (rain/river and snow/ice runoff) are applied only when -`Tendencies.SfcThicknessForcingTendencyEnable` is also enabled. ### Required input fields @@ -98,19 +94,10 @@ by the equivalent `ocn_comp_mct.F`. - 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 optional +- 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 direct heat part is - $Q_{\text{direct}} = Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$. - The mass-flux enthalpy part is - $Q_{\text{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 computed from EOS at top-layer salinity and pressure. - The applied heat flux is - $Q_{\text{direct}} + Q_{\text{mass}}$ when - `Tendencies.SfcThicknessForcingTendencyEnable` is true, and - $Q_{\text{direct}}$ otherwise. + 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). From 1226b96273b128d202f5b8ad8a531539c004454a Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 11:33:03 -0700 Subject: [PATCH 09/27] correction to pressure units and ctest --- components/omega/src/ocn/TendencyTerms.h | 12 ++--- components/omega/test/ocn/TendenciesTest.cpp | 49 ++++++-------------- 2 files changed, 21 insertions(+), 40 deletions(-) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 77b8ce93df04..b0a6591add83 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -434,12 +434,12 @@ class SfcTracerForcingOnCell { } if (TempIndex >= 0) { - const Real PTop = PressureMid(ICell, KTop); - const Real SaTop = SaltIndex >= 0 - ? TracerCell(SaltIndex, ICell, KTop) - : 0.0_Real; // not sure we want zero here? - const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTop, 0.0_Real); - const Real CtTop = TracerCell(TempIndex, ICell, KTop); + const Real PTopDb = PressureMid(ICell, KTop) * Pa2Db; + const Real SaTop = SaltIndex >= 0 + ? TracerCell(SaltIndex, ICell, KTop) + : 0.0_Real; // not sure we want zero here? + const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTopDb, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes // The enthalpy of liquid water is assumed to be: diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 1045701339dd..9114d4c7a116 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -309,10 +309,12 @@ int testTendencies() { DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; // Test surface tracer forcing with enthalpy terms - Err += testSfcTracerForcing(); + const int TracerForcingErr = testSfcTracerForcing(); + Err += TracerForcingErr; // Test surface thickness forcing with freshwater terms - Err += testSfcThicknessForcing(); + const int ThicknessForcingErr = testSfcThicknessForcing(); + Err += ThicknessForcingErr; // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; @@ -344,7 +346,6 @@ int testTendencies() { } Tendencies::clear(); - return Err; } @@ -491,6 +492,7 @@ int testSfcTracerForcing() { deepCopy(TracerTendBaseH, DefTendencies->TracerTend); const Real BaselineTempTend = TracerTendBaseH(TempIndex, ICellTest, KTop); const Real BaselineSaltTend = TracerTendBaseH(SaltIndex, ICellTest, KTop); + // Now enable SfcTracerForcing and compute again DefTendencies->SfcTracerForcing.Enabled = true; @@ -498,19 +500,13 @@ int testSfcTracerForcing() { ThickTimeLevel, VelTimeLevel, TracerTimeLevel, Time, Interval); - // Build two reference expectations for temperature tendency: - // 1) fixed estimate (expected to fail under strict tolerance), - // 2) TEOS-10 freezing CT (expected to pass under strict tolerance). - const Real CtFrzEstimate = -2.0_Real; - const Real ExpectedTempTendEstimate = - (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + - TestSnow * (Cp0Sw * CtFrzEstimate - LatIce)) * - HFluxFac; + // Build a reference expectations for temperature tendency: + // using TEOS-10 freezing CT (expected to pass under strict tolerance). HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); deepCopy(PressureMidH, VCoord->PressureMid); - const Real PTop = PressureMidH(ICellTest, KTop); - const Real CtFrzTeos = EosInst->calcCtFreezing(SaTopValue, PTop, 0.0_Real); + const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; + const Real CtFrzTeos = EosInst->calcCtFreezing(SaTopValue, PTopDb, 0.0_Real); const Real ExpectedTempTendTeos = (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + TestSnow * (Cp0Sw * CtFrzTeos - LatIce)) * @@ -530,20 +526,6 @@ int testSfcTracerForcing() { constexpr Real RelTol = 1.0e-10_Real; constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 - // Expected-fail check with fixed CtFrz estimate. - if (!isApprox(ComputedTempTend, ExpectedTempTendEstimate, RelTol, AbsTol)) { - LOG_INFO( - "TendenciesTest: expected tempTend fail because CtFrzEstimate != EOS " - "CtFrz - PASS"); - LOG_INFO("tempTend Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendEstimate, ComputedTempTend, - Kokkos::abs(ComputedTempTend - ExpectedTempTendEstimate)); - } else { - Err++; - LOG_ERROR("TendenciesTest: CtFrz estimate unexpectedly matched strict " - "reference - FAIL"); - } - // Expected-pass check with TEOS freezing CT reference. if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { Err++; @@ -559,9 +541,9 @@ int testSfcTracerForcing() { if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { Err++; LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); - LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, - ComputedSaltTend, - Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); + LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, + ComputedSaltTend, + Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); } else { LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); } @@ -651,7 +633,6 @@ int testSfcThicknessForcing() { LocSeaIceFreshWater(ICellTest) = TestSeaIceFreshWater; LocSeaIceSaltFlux(ICellTest) = TestSeaIceSaltFlux; }); - DefForcing->computeAll(); const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; @@ -721,9 +702,9 @@ int testSfcThicknessForcing() { if (!isApprox(ComputedThickTend, ExpectedThickTend, RelTol, AbsTol)) { Err++; LOG_ERROR("TendenciesTest: SfcThicknessForcing thickness tendency FAIL"); - LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedThickTend, - ComputedThickTend, - Kokkos::abs(ComputedThickTend - ExpectedThickTend)); + LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedThickTend, + ComputedThickTend, + Kokkos::abs(ComputedThickTend - ExpectedThickTend)); } else { LOG_INFO("TendenciesTest: SfcThicknessForcing thickness tendency PASS"); } From db218acfe3016a434292b8a0f10e138ee87e00f3 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 13:00:31 -0700 Subject: [PATCH 10/27] update due to fill values and review comments --- components/omega/doc/devGuide/Forcing.md | 14 +-- components/omega/doc/userGuide/Forcing.md | 16 +-- .../omega/doc/userGuide/TendencyTerms.md | 2 +- components/omega/src/ocn/Forcing.cpp | 30 +----- .../src/ocn/forcingVars/TracerForcingVars.cpp | 101 +++++++++--------- 5 files changed, 68 insertions(+), 95 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 869019e2469e..2a05a96b4c99 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -6,7 +6,7 @@ This page describes design and implementation details for forcing-related pathways in Omega, currently this includes: - Surface stress forcing (e.g. wind stress) -- Surface flux forcing (actively coupled or data-forced) +- Surface thickness and tracer flux forcing (actively coupled or data-forced) - Surface tracer restoring (soon to be ported) ## Surface stress forcing design @@ -38,9 +38,9 @@ pathways in Omega, currently this includes: - `Omega.Tendencies.SfcStressForcingTendencyEnable` - gates execution of surface stress forcing tendency kernel -## Surface flux forcing design +## Surface thickness and tracer flux forcing design -### Surface flux forcing data flow +### Surface thickness and tracer flux forcing data flow **Thickness equation pathway:** @@ -62,16 +62,16 @@ the surface layer pseudo-thickness. - `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 salinity (g/kg) in the surface layer. + and applies the external sea-ice salt flux to the top layer salt content thus impacting salinity. -### Surface flux forcing key classes/components +### 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 freshwater flux contribution: $\sum (\text{SnowFlux} + \text{RainFlux} + \text{EvaporationFlux} + \text{SeaIceFreshWaterFlux} + \text{IceRunoffFlux} + \text{RiverRunoffFlux} + \text{SeaIceSaltFlux}) / \rho_{sw}$ + - 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 @@ -88,7 +88,7 @@ the surface layer pseudo-thickness. - Calls `SfcThicknessForcingOnCell` in `computePseudoThicknessTendenciesOnly` - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring -### Surface flux forcing config coupling +### Surface thickness and tracer flux forcing config coupling - `Omega.Tendencies.SfcThicknessForcingTendencyEnable` - gates execution of coupled flux thickness kernel diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index e08754ee2651..85ecdfe1e1d8 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -5,7 +5,7 @@ This page documents the user-facing configuration and behavior for current forcing in Omega: - Surface stress forcing (e.g. wind stress) -- Coupled flux forcing +- Coupled tracer flux forcing (mass, energy and salt) - Surface tracer restoring ## Surface stress forcing @@ -40,15 +40,15 @@ 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 flux forcing +## Surface thickness and tracer flux forcing -Surface flux forcing applies ocean-atmosphere and ocean-sea ice fluxes from the other model +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 flux forcing configuration +### Surface thickness and tracer flux forcing configuration -Surface flux forcing is controlled by two configuration flags: +Surface thickness and tracer flux forcing is controlled by two configuration flags: ```yaml Omega: @@ -63,13 +63,13 @@ Omega: ### Required input fields -Coupled flux forcing uses 13 auxiliary fields organized by type: +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 input from sea-ice melt or formation +- `SeaIceFreshWaterFlux`: freshwater mass flux from sea-ice melt or formation - `IceRunoffFlux`: runoff from land ice - `RiverRunoffFlux`: runoff from rivers @@ -78,7 +78,7 @@ Coupled flux forcing uses 13 auxiliary fields organized by type: - `SensibleHeatFlux`: sensible heat transfer - `LongWaveHeatFluxUp`: upward longwave radiation - `LongWaveHeatFluxDown`: downward longwave radiation -- `SeaIceHeatFlux`: heat from sea-ice interaction +- `SeaIceHeatFlux`: heat/energy from sea-ice interaction (incl. enthalpy of meltwater) - `ShortWaveHeatFlux`: shortwave (solar) radiation **Salt mass flux (kg m⁻² s⁻¹):** diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 7847f0cfb357..c2cd4b3b4327 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -147,5 +147,5 @@ Tracer higer order convergence example of a cosine bell advected on a sphere sho ## See Also Additional information on forcing, including surface stress forcing, -surface flux forcing, and surface tracer restoring, is detailed in +surface thickness and tracer flux forcing, and surface tracer restoring, is detailed in [](omega-user-forcing). diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 3a3a54064dff..184ab5c30c13 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -166,7 +166,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; @@ -175,33 +176,6 @@ I4 Forcing::exchangeHalo() const { Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SnowFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RainFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.EvaporationFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo( - TracerForcing.SeaIceFreshWaterFluxCell, OnCell); - Err += - MeshHalo->exchangeFullArrayHalo(TracerForcing.IceRunoffFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RiverRunoffFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LatentHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SensibleHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LongWaveHeatFluxUpCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo( - TracerForcing.LongWaveHeatFluxDownCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.ShortWaveHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceSaltFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SurfInsituTemperature, - OnCell); - return Err; } diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp index 38016216a410..a00b4e9ee95c 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -24,7 +24,7 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, Mesh->NCellsSize), SeaIceHeatFluxCell("seaIceHeatFlux" + Suffix, Mesh->NCellsSize), ShortWaveHeatFluxCell("shortWaveHeatFlux" + Suffix, Mesh->NCellsSize), - SeaIceSaltFluxCell("seaIceSalinityFlux" + Suffix, Mesh->NCellsSize), + SeaIceSaltFluxCell("seaIceSaltFlux" + Suffix, Mesh->NCellsSize), SurfInsituTemperature("surfInsituTemperature" + Suffix, Mesh->NCellsSize) { deepCopy(SnowFluxCell, 0.0_Real); @@ -44,8 +44,7 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, } void TracerForcingVars::registerFields(const std::string &MeshName) const { - const Real FillValue = -9.99e30; - const int NDims = 1; + const int NDims = 1; std::vector DimNames(NDims); std::string DimSuffix; @@ -57,66 +56,66 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { DimNames[0] = "NCells" + DimSuffix; - auto SnowFluxField = Field::create( - SnowFluxCell.label(), "snow freshwater flux", "kg m^-2 s^-1", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - auto RainFluxField = Field::create( - RainFluxCell.label(), "rain freshwater flux", "kg m^-2 s^-1", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - auto EvaporationFluxField = Field::create( - EvaporationFluxCell.label(), "evaporation freshwater flux", - "kg m^-2 s^-1", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); + auto SnowFluxField = + Field::create(SnowFluxCell.label(), "snow freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto RainFluxField = + Field::create(RainFluxCell.label(), "rain freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto EvaporationFluxField = + Field::create(EvaporationFluxCell.label(), "evaporation freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); auto SeaIceFreshWaterFluxField = Field::create( SeaIceFreshWaterFluxCell.label(), "sea-ice freshwater flux", "kg m^-2 s^-1", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); - auto IceRunoffFluxField = Field::create( - IceRunoffFluxCell.label(), "ice runoff freshwater flux", "kg m^-2 s^-1", - "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); + std::numeric_limits::max(), NDims, DimNames); + auto IceRunoffFluxField = + Field::create(IceRunoffFluxCell.label(), "ice runoff freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); auto RiverRunoffFluxField = Field::create( RiverRunoffFluxCell.label(), "river runoff freshwater flux", "kg m^-2 s^-1", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); - - auto LatentHeatFluxField = Field::create( - LatentHeatFluxCell.label(), "latent heat flux", "W m^-2", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - auto SensibleHeatFluxField = Field::create( - SensibleHeatFluxCell.label(), "sensible heat flux", "W m^-2", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); + std::numeric_limits::max(), NDims, DimNames); + + auto LatentHeatFluxField = + Field::create(LatentHeatFluxCell.label(), "latent heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto SensibleHeatFluxField = + Field::create(SensibleHeatFluxCell.label(), "sensible heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); auto LongWaveHeatFluxUpField = Field::create( LongWaveHeatFluxUpCell.label(), "upward longwave heat flux", "W m^-2", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); + std::numeric_limits::max(), NDims, DimNames); auto LongWaveHeatFluxDownField = Field::create( LongWaveHeatFluxDownCell.label(), "downward longwave heat flux", "W m^-2", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), FillValue, NDims, DimNames); - auto SeaIceHeatFluxField = Field::create( - SeaIceHeatFluxCell.label(), "sea-ice heat flux", "W m^-2", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - auto ShortWaveHeatFluxField = Field::create( - ShortWaveHeatFluxCell.label(), "shortwave heat flux", "W m^-2", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - - auto SeaIceSaltFluxField = Field::create( - SeaIceSaltFluxCell.label(), "sea-ice salt flux", "kg m^-2 s^-1", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); - - auto SurfInsituTemperatureField = Field::create( - SurfInsituTemperature.label(), - "insitu (potential) temperature at surface layer", "degrees Celsius", "", - std::numeric_limits::lowest(), std::numeric_limits::max(), - FillValue, NDims, DimNames); + std::numeric_limits::max(), NDims, DimNames); + auto SeaIceHeatFluxField = + Field::create(SeaIceHeatFluxCell.label(), "sea-ice heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto ShortWaveHeatFluxField = + Field::create(ShortWaveHeatFluxCell.label(), "shortwave heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + + auto SeaIceSaltFluxField = + Field::create(SeaIceSaltFluxCell.label(), "sea-ice salt flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + + auto SurfInsituTemperatureField = + Field::create(SurfInsituTemperature.label(), + "insitu (potential) temperature at surface layer", + "degrees Celsius", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); FieldGroup::addFieldToGroup(SnowFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(RainFluxCell.label(), "Forcing"); From 965b5490812d1908c0beab69d9e500a23587c27f Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 14:41:04 -0700 Subject: [PATCH 11/27] resolve memory issue on GPUs --- components/omega/test/ocn/TendenciesTest.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 9114d4c7a116..fde3863c3dbb 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -373,9 +373,9 @@ int testSfcTracerForcing() { // Set up single test cell at top layer const I4 ICellTest = 0; - const I4 KTop = VCoord->MinLayerCell(ICellTest); + const I4 KTop = VCoord->MinLayerCellH(ICellTest); - if (KTop > VCoord->MaxLayerCell(ICellTest)) { + if (KTop > VCoord->MaxLayerCellH(ICellTest)) { LOG_ERROR("TendenciesTest: Test cell has no layers"); return -1; } @@ -579,9 +579,9 @@ int testSfcThicknessForcing() { // Set up single test cell at top layer const I4 ICellTest = 0; - const I4 KTop = VCoord->MinLayerCell(ICellTest); + const I4 KTop = VCoord->MinLayerCellH(ICellTest); - if (KTop > VCoord->MaxLayerCell(ICellTest)) { + if (KTop > VCoord->MaxLayerCellH(ICellTest)) { LOG_ERROR("TendenciesTest: Test cell has no layers for thickness test"); return -1; } From 554a9fff50b30202093395a30e5d68655ecb8f20 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Thu, 16 Jul 2026 13:08:27 -0700 Subject: [PATCH 12/27] inlines eos::calcPtFromCt in header with kokkos_function to fix gpu warning --- components/omega/src/ocn/Eos.cpp | 8 -------- components/omega/src/ocn/Eos.h | 7 ++++++- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index b2d7905e6789..848b95abe499 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -330,14 +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); diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index ca2ef47b352b..dd53799a2f56 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -1022,7 +1022,12 @@ class Eos { const Array2DReal &SpecVol); /// Convert Conservative Temperature to potential temperature - Real calcPtFromCt(const Real &Sa, const Real &Ct) const; + 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 Real calcCtFromPt(const Real &Sa, const Real &Pt) const; From 0fd8e353d7c38f2215210c7e90d0e7afb7bbed19 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 20 Jul 2026 09:38:07 -0700 Subject: [PATCH 13/27] adds linear and constant eos options to thermal forcing --- components/omega/src/ocn/Eos.cpp | 22 -- components/omega/src/ocn/Eos.h | 40 +++- components/omega/src/ocn/TendencyTerms.cpp | 2 +- components/omega/src/ocn/TendencyTerms.h | 4 +- components/omega/test/ocn/TendenciesTest.cpp | 208 ++++++++++++++++++- 5 files changed, 234 insertions(+), 42 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 848b95abe499..9620fe41f3f4 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -330,28 +330,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } -Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); - } - - 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() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index dd53799a2f56..4ac59871f7c2 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -1022,6 +1022,9 @@ class Eos { const Array2DReal &SpecVol); /// Convert Conservative Temperature to potential temperature + /// 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); @@ -1029,14 +1032,37 @@ class Eos { return Ct; } - /// Convert potential temperature to Conservative Temperature - Real calcCtFromPt(const Real &Sa, const Real &Pt) const; + /// 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; + } - /// 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; + /// Calculate freezing Conservative Temperature. + /// For TEOS-10, uses the Roquet et al. 75-term polynomial. + /// For LinearEos, uses a simple linear salinity-dependent approximation + /// consistent with the linear EOS philosophy (Sa in g/kg converted to PSU). + /// For ConstantEos, returns a constant approximate ocean freezing point. + KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); + } + if (EosChoice == 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(); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 37bfe6ee0500..65df82c4eb00 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -83,7 +83,7 @@ SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, const Eos *EosInst) : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - EosImpl(VCoord) {} + EosImpl(EosInst) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index b0a6591add83..efd1e43a9698 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -438,7 +438,7 @@ class SfcTracerForcingOnCell { const Real SaTop = SaltIndex >= 0 ? TracerCell(SaltIndex, ICell, KTop) : 0.0_Real; // not sure we want zero here? - const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTopDb, 0.0_Real); + const Real CtFrz = EosImpl->calcCtFreezing(SaTop, PTopDb, 0.0_Real); const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes @@ -469,7 +469,7 @@ class SfcTracerForcingOnCell { I4 SaltIndex; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; - Teos10Eos EosImpl; + const Eos *EosImpl; }; // Tracer horizontal advection term diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index fde3863c3dbb..01b0bff4cbdf 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -54,7 +54,8 @@ struct TestSetup { constexpr Geometry Geom = Geometry::Spherical; constexpr int NVertLayers = 60; -int testSfcTracerForcing(); +int testSfcTracerForcingTeos10(); +int testSfcTracerForcingLinear(); int testSfcThicknessForcing(); int initState() { @@ -308,9 +309,13 @@ int testTendencies() { DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; - // Test surface tracer forcing with enthalpy terms - const int TracerForcingErr = testSfcTracerForcing(); - Err += TracerForcingErr; + // Test surface tracer forcing with enthalpy terms (TEOS-10 CtFrz path) + const int TracerForcingTeos10Err = testSfcTracerForcingTeos10(); + Err += TracerForcingTeos10Err; + + // Test surface tracer forcing with LinearEos (linear CtFrz path) + const int TracerForcingLinearErr = testSfcTracerForcingLinear(); + Err += TracerForcingLinearErr; // Test surface thickness forcing with freshwater terms const int ThicknessForcingErr = testSfcThicknessForcing(); @@ -349,7 +354,7 @@ int testTendencies() { return Err; } -int testSfcTracerForcing() { +int testSfcTracerForcingTeos10() { int Err = 0; auto *VCoord = VertCoord::getDefault(); @@ -365,7 +370,8 @@ int testSfcTracerForcing() { const I4 SaltIndex = Tracers::IndxSalt; if (TempIndex < 0 || SaltIndex < 0) { - LOG_ERROR("TendenciesTest: Invalid tracer indices for SfcTracerForcing"); + LOG_ERROR( + "TendenciesTest: Invalid tracer indices for SfcTracerForcingTeos10"); return -1; } @@ -529,25 +535,207 @@ int testSfcTracerForcing() { // Expected-pass check with TEOS freezing CT reference. if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); + LOG_ERROR("TendenciesTest: SfcTracerForcingTeos10 temp tendency FAIL"); LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", ExpectedTempTendTeos, ComputedTempTend, Kokkos::abs(ComputedTempTend - ExpectedTempTendTeos)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcingTeos10 temp tendency PASS"); } // Check salinity tendency if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); + LOG_ERROR("TendenciesTest: SfcTracerForcingTeos10 salt tendency FAIL"); LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, ComputedSaltTend, Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcingTeos10 salt tendency PASS"); + } + + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; + DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; + DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; + DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; + DefTendencies->KEGrad.Enabled = OrigKEGrad; + DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; + DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; + DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; + DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; + DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; + DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; + + return Err; +} + +// Tests the SfcTracerForcing path using LinearEos. The EosChoice is +// temporarily set to LinearEos so that calcCtFreezing uses the linear +// salinity-dependent approximation instead of the TEOS-10 polynomial. +// Snow flux is applied so the CtFrz term is exercised. +int testSfcTracerForcingLinear() { + int Err = 0; + + auto *VCoord = VertCoord::getDefault(); + auto *DefTendencies = Tendencies::getDefault(); + auto *State = OceanState::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + auto *DefForcing = Forcing::getDefault(); + auto *EosInst = Eos::getInstance(); + + Array3DReal TracerArray = Tracers::getAll(0); + + const I4 TempIndex = Tracers::IndxTemp; + const I4 SaltIndex = Tracers::IndxSalt; + + if (TempIndex < 0 || SaltIndex < 0) { + LOG_ERROR("TendenciesTest: Invalid tracer indices for " + "SfcTracerForcingLinear"); + return -1; + } + + deepCopy(DefTendencies->TracerTend, 0._Real); + + const I4 ICellTest = 0; + const I4 KTop = VCoord->MinLayerCellH(ICellTest); + + if (KTop > VCoord->MaxLayerCellH(ICellTest)) { + LOG_ERROR("TendenciesTest: Test cell has no layers"); + return -1; + } + + const Real CtTopValue = 10.0_Real; // conservative temperature (degC) + const Real SaTopValue = 34.0_Real; // absolute salinity (g/kg) + + OMEGA_SCOPE(LocTracerArray, TracerArray); + Kokkos::parallel_for( + "SetTestTracersForcingNonTeos10", 1, KOKKOS_LAMBDA(int i) { + LocTracerArray(TempIndex, ICellTest, KTop) = CtTopValue; + LocTracerArray(SaltIndex, ICellTest, KTop) = SaTopValue; + }); + + auto &SensibleHeatFlux = DefForcing->TracerForcing.SensibleHeatFluxCell; + auto &LatentHeatFlux = DefForcing->TracerForcing.LatentHeatFluxCell; + auto &LongWaveHeatFluxUp = DefForcing->TracerForcing.LongWaveHeatFluxUpCell; + auto &LongWaveHeatFluxDown = + DefForcing->TracerForcing.LongWaveHeatFluxDownCell; + auto &SeaIceHeatFlux = DefForcing->TracerForcing.SeaIceHeatFluxCell; + auto &ShortWaveHeatFlux = DefForcing->TracerForcing.ShortWaveHeatFluxCell; + auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; + auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; + auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; + auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; + auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; + + deepCopy(SensibleHeatFlux, 0._Real); + deepCopy(LatentHeatFlux, 0._Real); + deepCopy(LongWaveHeatFluxUp, 0._Real); + deepCopy(LongWaveHeatFluxDown, 0._Real); + deepCopy(SeaIceHeatFlux, 0._Real); + deepCopy(ShortWaveHeatFlux, 0._Real); + deepCopy(RainFlux, 0._Real); + deepCopy(RiverRunoffFlux, 0._Real); + deepCopy(SnowFlux, 0._Real); + deepCopy(IceRunoffFlux, 0._Real); + deepCopy(SeaIceSaltFlux, 0._Real); + + // Only snow flux so the expected value depends solely on CtFrz. + const Real TestSnow = 5.0e-9_Real; // kg/m2/s + + OMEGA_SCOPE(LocSnowFlux, SnowFlux); + Kokkos::parallel_for( + "SetTestForcingNonTeos10", 1, + KOKKOS_LAMBDA(int i) { LocSnowFlux(ICellTest) = TestSnow; }); + + DefForcing->computeAll(); + + // Switch EOS to LinearEos so calcCtFreezing uses the linear approximation. + const EosType OrigEosChoice = EosInst->EosChoice; + EosInst->EosChoice = EosType::LinearEos; + + const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; + const bool OrigSfcThicknessEnabled = + DefTendencies->SfcThicknessForcing.Enabled; + const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; + const bool OrigPseudoThicknessDiv = + DefTendencies->PseudoThicknessFluxDiv.Enabled; + const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; + const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; + const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; + const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; + const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; + const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; + const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; + const bool OrigSurfaceTracerRestoring = + DefTendencies->SurfaceTracerRestoring.Enabled; + + DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->SfcThicknessForcing.Enabled = false; + DefTendencies->SfcTracerForcing.Enabled = false; + DefTendencies->PseudoThicknessFluxDiv.Enabled = false; + DefTendencies->PotentialVortHAdv.Enabled = false; + DefTendencies->KEGrad.Enabled = false; + DefTendencies->VelocityDiffusion.Enabled = false; + DefTendencies->VelocityHyperDiff.Enabled = false; + DefTendencies->TracerHorzAdv.Enabled = false; + DefTendencies->TracerDiffusion.Enabled = false; + DefTendencies->TracerHyperDiff.Enabled = false; + DefTendencies->SurfaceTracerRestoring.Enabled = false; + + int ThickTimeLevel = 0; + int VelTimeLevel = 0; + int TracerTimeLevel = 0; + TimeInstant Time; + TimeInterval Interval(1., TimeUnits::Seconds); + + // Compute baseline (vertical advection always on) + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + HostArray3DReal TracerTendBaseH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendBaseH, DefTendencies->TracerTend); + const Real BaselineTempTend = TracerTendBaseH(TempIndex, ICellTest, KTop); + + // Enable SfcTracerForcing and compute again + DefTendencies->SfcTracerForcing.Enabled = true; + + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + // Expected CtFrz from LinearEos path in Eos::calcCtFreezing: + // Tf = -0.054 * Sa * (35.0/35.16504) (no pressure dependence) + const Real CtFrzNonTeos = + -0.054_Real * SaTopValue * (35.0_Real / 35.16504_Real); + + // HeatFlux = Snow * (Cp0Sw * CtFrz - LatIce) + const Real ExpectedTempTend = + TestSnow * (Cp0Sw * CtFrzNonTeos - LatIce) * HFluxFac; + + HostArray3DReal TracerTendH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendH, DefTendencies->TracerTend); + const Real ComputedTempTend = + TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; + + if (!isApprox(ComputedTempTend, ExpectedTempTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcingLinear temp tendency FAIL"); + LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedTempTend, + ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcingLinear temp tendency PASS"); } + // Restore EOS choice and tendency flags + EosInst->EosChoice = OrigEosChoice; DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; From b5c586f784f65128f9c5b3cdd4c9429c556c5c68 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 20 Jul 2026 09:48:41 -0700 Subject: [PATCH 14/27] adds notes to docs and adds suggestions from review --- components/omega/doc/devGuide/Forcing.md | 4 +++ components/omega/doc/userGuide/Forcing.md | 2 +- .../src/ocn/forcingVars/TracerForcingVars.cpp | 28 +++++++++---------- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 2a05a96b4c99..92aff1fd9c8e 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -97,6 +97,10 @@ the surface layer pseudo-thickness. - 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. + ## Surface tracer restoring design ### Surface tracer restoring data flow diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index 85ecdfe1e1d8..01eeff4b5e3d 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -91,7 +91,7 @@ by the equivalent `ocn_comp_mct.F`. ### Notes -- Coupled fluxes are applied only at the surface layer (top active layer) for each cell. +- 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 diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp index a00b4e9ee95c..3deb7dfcb5da 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -10,22 +10,22 @@ namespace OMEGA { TracerForcingVars::TracerForcingVars(const std::string &Suffix, const HorzMesh *Mesh) - : SnowFluxCell("snowFlux" + Suffix, Mesh->NCellsSize), - RainFluxCell("rainFlux" + Suffix, Mesh->NCellsSize), - EvaporationFluxCell("evaporationFlux" + Suffix, Mesh->NCellsSize), - SeaIceFreshWaterFluxCell("seaIceFreshWaterFlux" + Suffix, + : SnowFluxCell("SnowFlux" + Suffix, Mesh->NCellsSize), + RainFluxCell("RainFlux" + Suffix, Mesh->NCellsSize), + EvaporationFluxCell("EvaporationFlux" + Suffix, Mesh->NCellsSize), + SeaIceFreshWaterFluxCell("SeaIceFreshWaterFlux" + Suffix, Mesh->NCellsSize), - IceRunoffFluxCell("iceRunoffFlux" + Suffix, Mesh->NCellsSize), - RiverRunoffFluxCell("riverRunoffFlux" + Suffix, Mesh->NCellsSize), - LatentHeatFluxCell("latentHeatFlux" + Suffix, Mesh->NCellsSize), - SensibleHeatFluxCell("sensibleHeatFlux" + Suffix, Mesh->NCellsSize), - LongWaveHeatFluxUpCell("longWaveHeatFluxUp" + Suffix, Mesh->NCellsSize), - LongWaveHeatFluxDownCell("longWaveHeatFluxDown" + Suffix, + IceRunoffFluxCell("IceRunoffFlux" + Suffix, Mesh->NCellsSize), + RiverRunoffFluxCell("RiverRunoffFlux" + Suffix, Mesh->NCellsSize), + LatentHeatFluxCell("LatentHeatFlux" + Suffix, Mesh->NCellsSize), + SensibleHeatFluxCell("SensibleHeatFlux" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxUpCell("LongWaveHeatFluxUp" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxDownCell("LongWaveHeatFluxDown" + Suffix, Mesh->NCellsSize), - SeaIceHeatFluxCell("seaIceHeatFlux" + Suffix, Mesh->NCellsSize), - ShortWaveHeatFluxCell("shortWaveHeatFlux" + Suffix, Mesh->NCellsSize), - SeaIceSaltFluxCell("seaIceSaltFlux" + Suffix, Mesh->NCellsSize), - SurfInsituTemperature("surfInsituTemperature" + Suffix, + SeaIceHeatFluxCell("SeaIceHeatFlux" + Suffix, Mesh->NCellsSize), + ShortWaveHeatFluxCell("ShortWaveHeatFlux" + Suffix, Mesh->NCellsSize), + SeaIceSaltFluxCell("SeaIceSaltFlux" + Suffix, Mesh->NCellsSize), + SurfInsituTemperature("SurfInsituTemperature" + Suffix, Mesh->NCellsSize) { deepCopy(SnowFluxCell, 0.0_Real); deepCopy(RainFluxCell, 0.0_Real); From a739dfc9073248cdbc46204e5b1f89605005b598 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 20 Jul 2026 21:26:40 -0700 Subject: [PATCH 15/27] Adds reset of forcing fields if not in stream --- components/omega/src/ocn/Forcing.cpp | 56 ++++++++++++++++++++++++++-- components/omega/src/ocn/Forcing.h | 3 ++ 2 files changed, 56 insertions(+), 3 deletions(-) diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 184ab5c30c13..5637f98191e9 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -152,7 +152,32 @@ void Forcing::readConfigOptions(Config *OmegaConfig) { } // Compute all forcing variables (dispatches to specific computations). -void Forcing::computeAll() const { computeSfcStressForcingOnEdge(); } +void Forcing::computeAll() const { + exchangeHalo(); + computeSfcStressForcingOnEdge(); +} + +// Reset forcing arrays so omitted optional fields remain zero after read. +void Forcing::resetArrays() { + deepCopy(SfcStressForcing.NormalStressEdge, 0.0_Real); + deepCopy(SfcStressForcing.ZonalStressCell, 0.0_Real); + deepCopy(SfcStressForcing.MeridStressCell, 0.0_Real); + + deepCopy(TracerForcing.SnowFluxCell, 0.0_Real); + deepCopy(TracerForcing.RainFluxCell, 0.0_Real); + deepCopy(TracerForcing.EvaporationFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); + deepCopy(TracerForcing.SurfInsituTemperature, 0.0_Real); +} // Compute edge-normal stress from cell-center zonal and meridional components. void Forcing::computeSfcStressForcingOnEdge() const { @@ -175,6 +200,30 @@ I4 Forcing::exchangeHalo() const { OnCell); Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SnowFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RainFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.EvaporationFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo( + TracerForcing.SeaIceFreshWaterFluxCell, OnCell); + Err += + MeshHalo->exchangeFullArrayHalo(TracerForcing.IceRunoffFluxCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RiverRunoffFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LatentHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SensibleHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LongWaveHeatFluxUpCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo( + TracerForcing.LongWaveHeatFluxDownCell, OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.ShortWaveHeatFluxCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceSaltFluxCell, + OnCell); return Err; } @@ -186,13 +235,14 @@ void Forcing::readStreamIntoArrays() { std::string StreamName = "Forcing"; + resetArrays(); + // Attempt to read stream; if unavailable, log and fall back to zero forcing. Err = IOStream::read(StreamName); if (Err.isFail()) { LOG_INFO("Forcing: Error while reading {} stream, using zero forcing", StreamName); - deepCopy(SfcStressForcing.ZonalStressCell, 0._Real); - deepCopy(SfcStressForcing.MeridStressCell, 0._Real); + resetArrays(); } I4 HaloErr = exchangeHalo(); diff --git a/components/omega/src/ocn/Forcing.h b/components/omega/src/ocn/Forcing.h index fda7b91d414e..de9749fecf28 100644 --- a/components/omega/src/ocn/Forcing.h +++ b/components/omega/src/ocn/Forcing.h @@ -71,6 +71,9 @@ class Forcing { /// Read forcing fields from input stream at startup void readStreamIntoArrays(); + /// Reset all forcing arrays to zero before reading optional fields + void resetArrays(); + /// Compute all forcing variables void computeAll() const; From 14065fcc07c4e7a5389d2a0772146251d70f491b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 22 Jul 2026 17:57:57 -0600 Subject: [PATCH 16/27] Update components/omega/doc/devGuide/Forcing.md --- components/omega/doc/devGuide/Forcing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 92aff1fd9c8e..8f9c97ead56d 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -7,7 +7,7 @@ pathways in Omega, currently this includes: - Surface stress forcing (e.g. wind stress) - Surface thickness and tracer flux forcing (actively coupled or data-forced) -- Surface tracer restoring (soon to be ported) +- Surface tracer restoring (soon to be ported as a field originating from the coupler) ## Surface stress forcing design From fded0f4b81cc4c43e6b034e01b2cb18e2c91bc4b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 22 Jul 2026 18:03:52 -0600 Subject: [PATCH 17/27] Fixup documentation --- components/omega/doc/devGuide/TendencyTerms.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 028ba9c02032..2fd762e28b3e 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -47,5 +47,5 @@ implemented: ## See Also -Additional information on forcing (surface stress, surface flux forcing, and +Additional information on forcing (surface stress, surface mass and tracer flux forcing, and surface tracer restoring) is detailed in [](omega-dev-forcing). From 9fa07ce44c3ac171ae58d5f6a23295d2bd2046cc Mon Sep 17 00:00:00 2001 From: Katherine Smith Date: Thu, 23 Jul 2026 19:26:08 -0400 Subject: [PATCH 18/27] remove SurfInsituTemp calcs --- components/omega/doc/devGuide/Forcing.md | 1 + components/omega/src/ocn/Forcing.cpp | 1 - .../src/ocn/forcingVars/TracerForcingVars.cpp | 49 +------------------ .../src/ocn/forcingVars/TracerForcingVars.h | 7 --- 4 files changed, 2 insertions(+), 56 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 8f9c97ead56d..7e69dc302bc6 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -100,6 +100,7 @@ the surface layer pseudo-thickness. ## 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 diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 5637f98191e9..bb005ad3565d 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -176,7 +176,6 @@ void Forcing::resetArrays() { deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); - deepCopy(TracerForcing.SurfInsituTemperature, 0.0_Real); } // Compute edge-normal stress from cell-center zonal and meridional components. diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp index 3deb7dfcb5da..a6478bb612c4 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -24,9 +24,7 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, Mesh->NCellsSize), SeaIceHeatFluxCell("SeaIceHeatFlux" + Suffix, Mesh->NCellsSize), ShortWaveHeatFluxCell("ShortWaveHeatFlux" + Suffix, Mesh->NCellsSize), - SeaIceSaltFluxCell("SeaIceSaltFlux" + Suffix, Mesh->NCellsSize), - SurfInsituTemperature("SurfInsituTemperature" + Suffix, - Mesh->NCellsSize) { + SeaIceSaltFluxCell("SeaIceSaltFlux" + Suffix, Mesh->NCellsSize) { deepCopy(SnowFluxCell, 0.0_Real); deepCopy(RainFluxCell, 0.0_Real); deepCopy(EvaporationFluxCell, 0.0_Real); @@ -40,7 +38,6 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, deepCopy(SeaIceHeatFluxCell, 0.0_Real); deepCopy(ShortWaveHeatFluxCell, 0.0_Real); deepCopy(SeaIceSaltFluxCell, 0.0_Real); - deepCopy(SurfInsituTemperature, 0.0_Real); } void TracerForcingVars::registerFields(const std::string &MeshName) const { @@ -111,12 +108,6 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { "kg m^-2 s^-1", "", std::numeric_limits::lowest(), std::numeric_limits::max(), NDims, DimNames); - auto SurfInsituTemperatureField = - Field::create(SurfInsituTemperature.label(), - "insitu (potential) temperature at surface layer", - "degrees Celsius", "", std::numeric_limits::lowest(), - std::numeric_limits::max(), NDims, DimNames); - FieldGroup::addFieldToGroup(SnowFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(RainFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(EvaporationFluxCell.label(), "Forcing"); @@ -143,7 +134,6 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { LongWaveHeatFluxDownField->attachData(LongWaveHeatFluxDownCell); SeaIceHeatFluxField->attachData(SeaIceHeatFluxCell); ShortWaveHeatFluxField->attachData(ShortWaveHeatFluxCell); - SurfInsituTemperatureField->attachData(SurfInsituTemperature); SeaIceSaltFluxField->attachData(SeaIceSaltFluxCell); } @@ -161,42 +151,5 @@ void TracerForcingVars::unregisterFields() const { Field::destroy(SeaIceHeatFluxCell.label()); Field::destroy(ShortWaveHeatFluxCell.label()); Field::destroy(SeaIceSaltFluxCell.label()); - Field::destroy(SurfInsituTemperature.label()); -} - -void TracerForcingVars::computeSurfInsituTemp(const Array3DReal &TracerArray, - const VertCoord *VCoord, - const Eos *EosInst) const { - const int IndxTemp = Tracers::IndxTemp; - const int IndxSalt = Tracers::IndxSalt; - - // Skip computation if temperature or salinity tracers are not defined - if (IndxTemp < 0 || IndxSalt < 0) { - return; - } - - OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); - OMEGA_SCOPE(LocSurfInsituTemp, SurfInsituTemperature); - - int NCellsOwned = SurfInsituTemperature.extent_int(0); - - parallelFor( - "TracerForcing:computeSurfInsituTemp", {NCellsOwned}, - KOKKOS_LAMBDA(int ICell) { - const int KMin = LocMinLayerCell(ICell); - const int KMax = LocMaxLayerCell(ICell); - - // Only compute for valid ocean cells - if (KMin <= KMax) { - const Real ConservTemp = TracerArray(IndxTemp, ICell, KMin); - const Real AbsSalinity = TracerArray(IndxSalt, ICell, KMin); - - // Call EOS function to compute potential temperature from - // conservative temperature at surface (reference pressure = 0) - LocSurfInsituTemp(ICell) = - EosInst->calcPtFromCt(AbsSalinity, ConservTemp); - } - }); } } // namespace OMEGA diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.h b/components/omega/src/ocn/forcingVars/TracerForcingVars.h index 1a0747121ea2..e38d9948f672 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.h +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.h @@ -31,17 +31,10 @@ class TracerForcingVars { Array1DReal SeaIceSaltFluxCell; - Array1DReal SurfInsituTemperature; - TracerForcingVars(const std::string &Suffix, const HorzMesh *Mesh); void registerFields(const std::string &MeshName) const; void unregisterFields() const; - - /// Compute surface insitu temperature from conservative temperature - void computeSurfInsituTemp(const Array3DReal &TracerArray, - const VertCoord *VCoord, - const Eos *EosInst) const; }; } // namespace OMEGA From 5bed0595bf448be8b77c89e46fd93b93699c1e0b Mon Sep 17 00:00:00 2001 From: Katherine Smith Date: Fri, 24 Jul 2026 01:18:27 -0400 Subject: [PATCH 19/27] fixes GPU failures on frontier --- components/omega/src/ocn/Eos.cpp | 23 ++++++++++++++++++++++ components/omega/src/ocn/Eos.h | 22 ++------------------- components/omega/src/ocn/TendencyTerms.cpp | 2 +- components/omega/src/ocn/TendencyTerms.h | 4 ++-- 4 files changed, 28 insertions(+), 23 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 9620fe41f3f4..6db728db57a0 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -330,6 +330,29 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } +Real Eos::calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); + } + if (EosChoice == 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; + } + if (EosChoice == EosType::ConstantEos) { + // Constant approximate ocean freezing point (degC) + return -1.9_Real; + } + ABORT_ERROR( + "Eos::calcCtFreezing: CT freezing temperature is only " + "implemented for TEOS-10, Linear, and Constant EOS types. " + "Support for the current EOS choice has not yet been developed."); + return 0; +} + /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 4ac59871f7c2..e0d3b5140129 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -1043,26 +1043,8 @@ class Eos { return Pt; } - /// Calculate freezing Conservative Temperature. - /// For TEOS-10, uses the Roquet et al. 75-term polynomial. - /// For LinearEos, uses a simple linear salinity-dependent approximation - /// consistent with the linear EOS philosophy (Sa in g/kg converted to PSU). - /// For ConstantEos, returns a constant approximate ocean freezing point. - KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); - } - if (EosChoice == 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; - } + Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const; /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 65df82c4eb00..37bfe6ee0500 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -83,7 +83,7 @@ SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, const Eos *EosInst) : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - EosImpl(EosInst) {} + EosImpl(VCoord) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index efd1e43a9698..b0a6591add83 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -438,7 +438,7 @@ class SfcTracerForcingOnCell { const Real SaTop = SaltIndex >= 0 ? TracerCell(SaltIndex, ICell, KTop) : 0.0_Real; // not sure we want zero here? - const Real CtFrz = EosImpl->calcCtFreezing(SaTop, PTopDb, 0.0_Real); + const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTopDb, 0.0_Real); const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes @@ -469,7 +469,7 @@ class SfcTracerForcingOnCell { I4 SaltIndex; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; - const Eos *EosImpl; + Teos10Eos EosImpl; }; // Tracer horizontal advection term From 08640506976e90b2176b7def18fde791db87a79e Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Tue, 28 Jul 2026 17:47:10 -0400 Subject: [PATCH 20/27] Revert "fixes GPU failures on frontier" This reverts commit d32ee4d7352fc7fed80766ac2c7e0943161d99b3. --- components/omega/src/ocn/Eos.cpp | 23 ---------------------- components/omega/src/ocn/Eos.h | 22 +++++++++++++++++++-- components/omega/src/ocn/TendencyTerms.cpp | 2 +- components/omega/src/ocn/TendencyTerms.h | 4 ++-- 4 files changed, 23 insertions(+), 28 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 6db728db57a0..9620fe41f3f4 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -330,29 +330,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } -Real Eos::calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); - } - if (EosChoice == 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; - } - if (EosChoice == EosType::ConstantEos) { - // Constant approximate ocean freezing point (degC) - return -1.9_Real; - } - ABORT_ERROR( - "Eos::calcCtFreezing: CT freezing temperature is only " - "implemented for TEOS-10, Linear, and Constant EOS types. " - "Support for the current EOS choice has not yet been developed."); - return 0; -} - /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index e0d3b5140129..4ac59871f7c2 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -1043,8 +1043,26 @@ class Eos { return Pt; } - Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const; + /// Calculate freezing Conservative Temperature. + /// For TEOS-10, uses the Roquet et al. 75-term polynomial. + /// For LinearEos, uses a simple linear salinity-dependent approximation + /// consistent with the linear EOS philosophy (Sa in g/kg converted to PSU). + /// For ConstantEos, returns a constant approximate ocean freezing point. + KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); + } + if (EosChoice == 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(); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 37bfe6ee0500..65df82c4eb00 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -83,7 +83,7 @@ SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, const Eos *EosInst) : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - EosImpl(VCoord) {} + EosImpl(EosInst) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index b0a6591add83..efd1e43a9698 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -438,7 +438,7 @@ class SfcTracerForcingOnCell { const Real SaTop = SaltIndex >= 0 ? TracerCell(SaltIndex, ICell, KTop) : 0.0_Real; // not sure we want zero here? - const Real CtFrz = EosImpl.calcCtFreezing(SaTop, PTopDb, 0.0_Real); + const Real CtFrz = EosImpl->calcCtFreezing(SaTop, PTopDb, 0.0_Real); const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes @@ -469,7 +469,7 @@ class SfcTracerForcingOnCell { I4 SaltIndex; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; - Teos10Eos EosImpl; + const Eos *EosImpl; }; // Tracer horizontal advection term From b35ff14c32b9774ddfcd20f91dd2ea2121dfef25 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Thu, 30 Jul 2026 16:35:31 -0400 Subject: [PATCH 21/27] fixed GPU isssues on Frontier --- components/omega/src/ocn/Eos.h | 20 +- components/omega/src/ocn/TendencyTerms.cpp | 2 +- components/omega/src/ocn/TendencyTerms.h | 7 +- components/omega/test/ocn/EosTest.cpp | 47 +++- components/omega/test/ocn/TendenciesTest.cpp | 238 ++----------------- 5 files changed, 82 insertions(+), 232 deletions(-) diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 4ac59871f7c2..dceb0d1a9cc0 100644 --- a/components/omega/src/ocn/Eos.h +++ b/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; @@ -1043,17 +1043,17 @@ class Eos { return Pt; } - /// Calculate freezing Conservative Temperature. + /// 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 - /// consistent with the linear EOS philosophy (Sa in g/kg converted to PSU). + /// For LinearEos, uses a simple linear salinity-dependent approximation. /// For ConstantEos, returns a constant approximate ocean freezing point. - KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFreezing(Sa, P, SaturationFract); + 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 (EosChoice == EosType::LinearEos) { + if (Choice == EosType::LinearEos) { // Linear salinity-dependent freezing point; coefficient -0.054 // degC/PSU with absolute-to-practical salinity conversion (g/kg -> // PSU). diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 65df82c4eb00..2353f38049bd 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -83,7 +83,7 @@ SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, const Eos *EosInst) : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - EosImpl(EosInst) {} + EosChoice(EosInst->EosChoice) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index efd1e43a9698..a2a8dc09492c 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -438,8 +438,9 @@ class SfcTracerForcingOnCell { const Real SaTop = SaltIndex >= 0 ? TracerCell(SaltIndex, ICell, KTop) : 0.0_Real; // not sure we want zero here? - const Real CtFrz = EosImpl->calcCtFreezing(SaTop, PTopDb, 0.0_Real); - const Real CtTop = TracerCell(TempIndex, ICell, KTop); + const Real CtFrz = + Eos::calcCtFreezing(EosChoice, SaTop, PTopDb, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes // The enthalpy of liquid water is assumed to be: @@ -469,7 +470,7 @@ class SfcTracerForcingOnCell { I4 SaltIndex; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; - const Eos *EosImpl; + EosType EosChoice; }; // Tracer horizontal advection term diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index ed12201a6750..48e474491986 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -1089,6 +1089,42 @@ void testEosConstantDerivs() { return; } +/// Test all Eos::calcCtFreezing pathways (Teos10, Linear, Constant) +void testCalcCtFreezing() { + const Real RTol = 1e-10; + + constexpr Real SaturationFrac = 0.0; + constexpr Real PDb = 500.0; // pressure in dbar for GSW pathway + constexpr Real SaLocal = 32.0; + + const Real CtTeosExpected = + gsw_ct_freezing_poly(SaLocal, PDb, SaturationFrac); + const Real CtTeos = + Eos::calcCtFreezing(EosType::Teos10Eos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtTeos, CtTeosExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Teos10 FAIL, expected {}, got {}", + CtTeosExpected, CtTeos); + } + + const Real CtLinearExpected = -0.054_Real * SaLocal / Psu2Gpkg; + const Real CtLinear = + Eos::calcCtFreezing(EosType::LinearEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtLinear, CtLinearExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Linear FAIL, expected {}, got {}", + CtLinearExpected, CtLinear); + } + + const Real CtConstExpected = -1.9_Real; + const Real CtConst = + Eos::calcCtFreezing(EosType::ConstantEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtConst, CtConstExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Constant FAIL, expected {}, got {}", + CtConstExpected, CtConst); + } + + return; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { Eos::destroyInstance(); @@ -1150,22 +1186,22 @@ void checkValueGswcN2() { } /// Test that the calcCtFreezing function returns the expected value -void checkValueCtFreezing() { +void checkValueGswcCtFreezing() { const Real RTol = 1e-10; - Teos10Eos TestEos(VertCoord::getDefault()); constexpr Real SaturationFrac = 0.0; constexpr Real P = 500.0 * Db2Pa; // Convert dbar to Pa constexpr Real Sa = 32.0; /// Get freezing temperature from GSW-C library double CtFreezGswc = gsw_ct_freezing_poly(Sa, P * Pa2Db, SaturationFrac); - double CtFreez = TestEos.calcCtFreezing(Sa, P * Pa2Db, SaturationFrac); + double CtFreez = + Teos10Eos::calcCtFreezingTeos10(Sa, P * Pa2Db, SaturationFrac); /// Check the value against the GSW-C value bool Check = isApprox(CtFreezGswc, CtFreez, RTol); if (!Check) { - ABORT_ERROR("checkValueCtFreezing: CtFreez FAIL, expected {}, got {}", + ABORT_ERROR("checkValueGswcCtFreezing: CtFreez FAIL, expected {}, got {}", CtFreezGswc, CtFreez); } return; @@ -1483,7 +1519,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); - checkValueCtFreezing(); + checkValueGswcCtFreezing(); checkValueGswcCtFromPt(); checkValueGswcPtFromCt(); checkValueGswcSpecVolDerivs(); @@ -1500,6 +1536,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { testEosTeos10Derivs(); testEosLinearDerivs(); testEosConstantDerivs(); + testCalcCtFreezing(); finalizeEosTest(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 01b0bff4cbdf..13f8ed2ba2a4 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -54,8 +54,7 @@ struct TestSetup { constexpr Geometry Geom = Geometry::Spherical; constexpr int NVertLayers = 60; -int testSfcTracerForcingTeos10(); -int testSfcTracerForcingLinear(); +int testSfcTracerForcing(); int testSfcThicknessForcing(); int initState() { @@ -307,15 +306,19 @@ int testTendencies() { "NormalVelocityTend"); } + const Real NormVelTendSum = + sum(DefTendencies->NormalVelocityTend, Mesh->NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); + if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { + Err++; + LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); + } + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; // Test surface tracer forcing with enthalpy terms (TEOS-10 CtFrz path) - const int TracerForcingTeos10Err = testSfcTracerForcingTeos10(); - Err += TracerForcingTeos10Err; - - // Test surface tracer forcing with LinearEos (linear CtFrz path) - const int TracerForcingLinearErr = testSfcTracerForcingLinear(); - Err += TracerForcingLinearErr; + const int TracerForcingErr = testSfcTracerForcing(); + Err += TracerForcingErr; // Test surface thickness forcing with freshwater terms const int ThicknessForcingErr = testSfcThicknessForcing(); @@ -323,7 +326,6 @@ int testTendencies() { // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; - int NEdgesOwned = Mesh->NEdgesOwned; int NTracers = Tracers::getNumTracers(); const Real PseudoThickTendSum = @@ -334,14 +336,6 @@ int testTendencies() { LOG_ERROR("TendenciesTest: PseudoThickTend FAIL"); } - const Real NormVelTendSum = - sum(DefTendencies->NormalVelocityTend, NEdgesOwned, - VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); - if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { - Err++; - LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); - } - const Real TraceTendSum = sum(DefTendencies->TracerTend, NTracers, NCellsOwned, VCoord->MinLayerCell, VCoord->MaxLayerCell); @@ -354,7 +348,7 @@ int testTendencies() { return Err; } -int testSfcTracerForcingTeos10() { +int testSfcTracerForcing() { int Err = 0; auto *VCoord = VertCoord::getDefault(); @@ -370,8 +364,7 @@ int testSfcTracerForcingTeos10() { const I4 SaltIndex = Tracers::IndxSalt; if (TempIndex < 0 || SaltIndex < 0) { - LOG_ERROR( - "TendenciesTest: Invalid tracer indices for SfcTracerForcingTeos10"); + LOG_ERROR("TendenciesTest: Invalid tracer indices for SfcTracerForcing"); return -1; } @@ -511,11 +504,12 @@ int testSfcTracerForcingTeos10() { HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); deepCopy(PressureMidH, VCoord->PressureMid); - const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; - const Real CtFrzTeos = EosInst->calcCtFreezing(SaTopValue, PTopDb, 0.0_Real); - const Real ExpectedTempTendTeos = + const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; + const Real CtFrz = + Eos::calcCtFreezing(EosInst->EosChoice, SaTopValue, PTopDb, 0.0_Real); + const Real ExpectedTempTend = (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + - TestSnow * (Cp0Sw * CtFrzTeos - LatIce)) * + TestSnow * (Cp0Sw * CtFrz - LatIce)) * HFluxFac; // SaltTend = SeaIceSaltFlux * SFluxFac @@ -533,209 +527,27 @@ int testSfcTracerForcingTeos10() { constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 // Expected-pass check with TEOS freezing CT reference. - if (!isApprox(ComputedTempTend, ExpectedTempTendTeos, RelTol, AbsTol)) { + if (!isApprox(ComputedTempTend, ExpectedTempTend, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcingTeos10 temp tendency FAIL"); + LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", - ExpectedTempTendTeos, ComputedTempTend, - Kokkos::abs(ComputedTempTend - ExpectedTempTendTeos)); + ExpectedTempTend, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTend)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcingTeos10 temp tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); } // Check salinity tendency if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcingTeos10 salt tendency FAIL"); + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, ComputedSaltTend, Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); } else { - LOG_INFO("TendenciesTest: SfcTracerForcingTeos10 salt tendency PASS"); - } - - DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; - DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; - DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; - DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; - DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; - DefTendencies->KEGrad.Enabled = OrigKEGrad; - DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; - DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; - DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; - DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; - DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; - DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; - - return Err; -} - -// Tests the SfcTracerForcing path using LinearEos. The EosChoice is -// temporarily set to LinearEos so that calcCtFreezing uses the linear -// salinity-dependent approximation instead of the TEOS-10 polynomial. -// Snow flux is applied so the CtFrz term is exercised. -int testSfcTracerForcingLinear() { - int Err = 0; - - auto *VCoord = VertCoord::getDefault(); - auto *DefTendencies = Tendencies::getDefault(); - auto *State = OceanState::getDefault(); - auto *AuxState = AuxiliaryState::getDefault(); - auto *DefForcing = Forcing::getDefault(); - auto *EosInst = Eos::getInstance(); - - Array3DReal TracerArray = Tracers::getAll(0); - - const I4 TempIndex = Tracers::IndxTemp; - const I4 SaltIndex = Tracers::IndxSalt; - - if (TempIndex < 0 || SaltIndex < 0) { - LOG_ERROR("TendenciesTest: Invalid tracer indices for " - "SfcTracerForcingLinear"); - return -1; - } - - deepCopy(DefTendencies->TracerTend, 0._Real); - - const I4 ICellTest = 0; - const I4 KTop = VCoord->MinLayerCellH(ICellTest); - - if (KTop > VCoord->MaxLayerCellH(ICellTest)) { - LOG_ERROR("TendenciesTest: Test cell has no layers"); - return -1; - } - - const Real CtTopValue = 10.0_Real; // conservative temperature (degC) - const Real SaTopValue = 34.0_Real; // absolute salinity (g/kg) - - OMEGA_SCOPE(LocTracerArray, TracerArray); - Kokkos::parallel_for( - "SetTestTracersForcingNonTeos10", 1, KOKKOS_LAMBDA(int i) { - LocTracerArray(TempIndex, ICellTest, KTop) = CtTopValue; - LocTracerArray(SaltIndex, ICellTest, KTop) = SaTopValue; - }); - - auto &SensibleHeatFlux = DefForcing->TracerForcing.SensibleHeatFluxCell; - auto &LatentHeatFlux = DefForcing->TracerForcing.LatentHeatFluxCell; - auto &LongWaveHeatFluxUp = DefForcing->TracerForcing.LongWaveHeatFluxUpCell; - auto &LongWaveHeatFluxDown = - DefForcing->TracerForcing.LongWaveHeatFluxDownCell; - auto &SeaIceHeatFlux = DefForcing->TracerForcing.SeaIceHeatFluxCell; - auto &ShortWaveHeatFlux = DefForcing->TracerForcing.ShortWaveHeatFluxCell; - auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; - auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; - auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; - auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; - auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; - - deepCopy(SensibleHeatFlux, 0._Real); - deepCopy(LatentHeatFlux, 0._Real); - deepCopy(LongWaveHeatFluxUp, 0._Real); - deepCopy(LongWaveHeatFluxDown, 0._Real); - deepCopy(SeaIceHeatFlux, 0._Real); - deepCopy(ShortWaveHeatFlux, 0._Real); - deepCopy(RainFlux, 0._Real); - deepCopy(RiverRunoffFlux, 0._Real); - deepCopy(SnowFlux, 0._Real); - deepCopy(IceRunoffFlux, 0._Real); - deepCopy(SeaIceSaltFlux, 0._Real); - - // Only snow flux so the expected value depends solely on CtFrz. - const Real TestSnow = 5.0e-9_Real; // kg/m2/s - - OMEGA_SCOPE(LocSnowFlux, SnowFlux); - Kokkos::parallel_for( - "SetTestForcingNonTeos10", 1, - KOKKOS_LAMBDA(int i) { LocSnowFlux(ICellTest) = TestSnow; }); - - DefForcing->computeAll(); - - // Switch EOS to LinearEos so calcCtFreezing uses the linear approximation. - const EosType OrigEosChoice = EosInst->EosChoice; - EosInst->EosChoice = EosType::LinearEos; - - const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; - const bool OrigSfcThicknessEnabled = - DefTendencies->SfcThicknessForcing.Enabled; - const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; - const bool OrigPseudoThicknessDiv = - DefTendencies->PseudoThicknessFluxDiv.Enabled; - const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; - const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; - const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; - const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; - const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; - const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; - const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; - const bool OrigSurfaceTracerRestoring = - DefTendencies->SurfaceTracerRestoring.Enabled; - - DefTendencies->SfcStressForcing.Enabled = false; - DefTendencies->SfcThicknessForcing.Enabled = false; - DefTendencies->SfcTracerForcing.Enabled = false; - DefTendencies->PseudoThicknessFluxDiv.Enabled = false; - DefTendencies->PotentialVortHAdv.Enabled = false; - DefTendencies->KEGrad.Enabled = false; - DefTendencies->VelocityDiffusion.Enabled = false; - DefTendencies->VelocityHyperDiff.Enabled = false; - DefTendencies->TracerHorzAdv.Enabled = false; - DefTendencies->TracerDiffusion.Enabled = false; - DefTendencies->TracerHyperDiff.Enabled = false; - DefTendencies->SurfaceTracerRestoring.Enabled = false; - - int ThickTimeLevel = 0; - int VelTimeLevel = 0; - int TracerTimeLevel = 0; - TimeInstant Time; - TimeInterval Interval(1., TimeUnits::Seconds); - - // Compute baseline (vertical advection always on) - DefTendencies->computeAllTendencies(State, AuxState, TracerArray, - ThickTimeLevel, VelTimeLevel, - TracerTimeLevel, Time, Interval); - - HostArray3DReal TracerTendBaseH = - createHostMirrorCopy(DefTendencies->TracerTend); - deepCopy(TracerTendBaseH, DefTendencies->TracerTend); - const Real BaselineTempTend = TracerTendBaseH(TempIndex, ICellTest, KTop); - - // Enable SfcTracerForcing and compute again - DefTendencies->SfcTracerForcing.Enabled = true; - - DefTendencies->computeAllTendencies(State, AuxState, TracerArray, - ThickTimeLevel, VelTimeLevel, - TracerTimeLevel, Time, Interval); - - // Expected CtFrz from LinearEos path in Eos::calcCtFreezing: - // Tf = -0.054 * Sa * (35.0/35.16504) (no pressure dependence) - const Real CtFrzNonTeos = - -0.054_Real * SaTopValue * (35.0_Real / 35.16504_Real); - - // HeatFlux = Snow * (Cp0Sw * CtFrz - LatIce) - const Real ExpectedTempTend = - TestSnow * (Cp0Sw * CtFrzNonTeos - LatIce) * HFluxFac; - - HostArray3DReal TracerTendH = - createHostMirrorCopy(DefTendencies->TracerTend); - deepCopy(TracerTendH, DefTendencies->TracerTend); - const Real ComputedTempTend = - TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; - - constexpr Real RelTol = 1.0e-10_Real; - constexpr Real AbsTol = 1.0e-12_Real; - - if (!isApprox(ComputedTempTend, ExpectedTempTend, RelTol, AbsTol)) { - Err++; - LOG_ERROR("TendenciesTest: SfcTracerForcingLinear temp tendency FAIL"); - LOG_ERROR(" Expected: {}, Computed: {}, Diff: {}", ExpectedTempTend, - ComputedTempTend, - Kokkos::abs(ComputedTempTend - ExpectedTempTend)); - } else { - LOG_INFO("TendenciesTest: SfcTracerForcingLinear temp tendency PASS"); + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); } - // Restore EOS choice and tendency flags - EosInst->EosChoice = OrigEosChoice; DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; From c5611014f967a4ad0ed9c68bfcbba095faf356ad Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Fri, 31 Jul 2026 15:00:26 -0400 Subject: [PATCH 22/27] fixes omega_pr errors on Frontier --- components/omega/src/ocn/Forcing.cpp | 121 ++++++++++++++++----------- components/omega/src/ocn/Forcing.h | 2 + 2 files changed, 72 insertions(+), 51 deletions(-) diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index bb005ad3565d..69e7c502ff64 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -37,14 +37,22 @@ 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); + if (SfcStressFieldsEnabled) { + SfcStressForcing.registerFields(MeshName); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.registerFields(MeshName); + } } // Unregister surface stress fields from IO streams. void Forcing::unregisterFields() const { - SfcStressForcing.unregisterFields(); - TracerForcing.unregisterFields(); + if (SfcStressFieldsEnabled) { + SfcStressForcing.unregisterFields(); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.unregisterFields(); + } } // Create and register a non-default forcing instance. @@ -88,11 +96,10 @@ void Forcing::init() { ABORT_ERROR("Forcing: failed to initialize default forcing state"); } - DefaultForcing->registerFields(DefMesh->MeshName); - Config *OmegaConfig = Config::getOmegaConfig(); OMEGA_REQUIRE(OmegaConfig, "Null OmegaConfig pointer in Forcing::init"); DefaultForcing->readConfigOptions(OmegaConfig); + DefaultForcing->registerFields(DefMesh->MeshName); // for now, forcing fields are read at start-up only. // to be extended to include switch from standalone to coupled. // to be moved to a Forcing->prepareForStep(SimTime) method later. @@ -149,6 +156,30 @@ void Forcing::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("Forcing: Unknown InterpType requested"); } + + Config TendConfig("Tendencies"); + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, "Forcing: Tendencies group not found in Config"); + + Err += + TendConfig.get("SfcStressForcingTendencyEnable", SfcStressFieldsEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcStressForcingTendencyEnable not found " + "in Tendencies config"); + + bool SfcThicknessForcingEnabled = false; + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + SfcThicknessForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcThicknessForcingTendencyEnable not " + "found in Tendencies config"); + + bool SfcTracerForcingEnabled = false; + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + SfcTracerForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcTracerForcingTendencyEnable not found " + "in Tendencies config"); + + TracerForcingFieldsEnabled = + SfcThicknessForcingEnabled || SfcTracerForcingEnabled; } // Compute all forcing variables (dispatches to specific computations). @@ -159,23 +190,27 @@ void Forcing::computeAll() const { // Reset forcing arrays so omitted optional fields remain zero after read. void Forcing::resetArrays() { - deepCopy(SfcStressForcing.NormalStressEdge, 0.0_Real); - deepCopy(SfcStressForcing.ZonalStressCell, 0.0_Real); - deepCopy(SfcStressForcing.MeridStressCell, 0.0_Real); - - deepCopy(TracerForcing.SnowFluxCell, 0.0_Real); - deepCopy(TracerForcing.RainFluxCell, 0.0_Real); - deepCopy(TracerForcing.EvaporationFluxCell, 0.0_Real); - deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); - deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); - deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); - deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); - deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); - deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); - deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); - deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); - deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); - deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); + if (SfcStressFieldsEnabled) { + deepCopy(SfcStressForcing.NormalStressEdge, 0.0_Real); + deepCopy(SfcStressForcing.ZonalStressCell, 0.0_Real); + deepCopy(SfcStressForcing.MeridStressCell, 0.0_Real); + } + + if (TracerForcingFieldsEnabled) { + deepCopy(TracerForcing.SnowFluxCell, 0.0_Real); + deepCopy(TracerForcing.RainFluxCell, 0.0_Real); + deepCopy(TracerForcing.EvaporationFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); + } } // Compute edge-normal stress from cell-center zonal and meridional components. @@ -195,34 +230,12 @@ void Forcing::computeSfcStressForcingOnEdge() const { I4 Forcing::exchangeHalo() const { I4 Err = 0; - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SnowFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RainFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.EvaporationFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo( - TracerForcing.SeaIceFreshWaterFluxCell, OnCell); - Err += - MeshHalo->exchangeFullArrayHalo(TracerForcing.IceRunoffFluxCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.RiverRunoffFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LatentHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SensibleHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.LongWaveHeatFluxUpCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo( - TracerForcing.LongWaveHeatFluxDownCell, OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.ShortWaveHeatFluxCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(TracerForcing.SeaIceSaltFluxCell, - OnCell); + if (SfcStressFieldsEnabled) { + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, + OnCell); + } return Err; } @@ -236,6 +249,12 @@ void Forcing::readStreamIntoArrays() { resetArrays(); + // Nothing to read if neither stress nor tracer forcing tendencies are + // enabled. + if (!SfcStressFieldsEnabled && !TracerForcingFieldsEnabled) { + return; + } + // Attempt to read stream; if unavailable, log and fall back to zero forcing. Err = IOStream::read(StreamName); if (Err.isFail()) { diff --git a/components/omega/src/ocn/Forcing.h b/components/omega/src/ocn/Forcing.h index de9749fecf28..b061bdcf25d7 100644 --- a/components/omega/src/ocn/Forcing.h +++ b/components/omega/src/ocn/Forcing.h @@ -91,6 +91,8 @@ class Forcing { const HorzMesh *Mesh; Halo *MeshHalo; + bool SfcStressFieldsEnabled = false; + bool TracerForcingFieldsEnabled = false; static Forcing *DefaultForcing; static std::map> AllForcing; From e7999077a41d07eddbaede052725a1400548ad45 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Wed, 12 Aug 2026 18:35:07 -0400 Subject: [PATCH 23/27] changes assumptions about enthalpy calculation, per PR review comments --- components/omega/src/ocn/TendencyTerms.h | 22 +++++++++++--------- components/omega/test/ocn/TendenciesTest.cpp | 12 ++--------- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index a2a8dc09492c..8c45b7bee392 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -434,28 +434,30 @@ class SfcTracerForcingOnCell { } if (TempIndex >= 0) { - const Real PTopDb = PressureMid(ICell, KTop) * Pa2Db; - const Real SaTop = SaltIndex >= 0 - ? TracerCell(SaltIndex, ICell, KTop) - : 0.0_Real; // not sure we want zero here? - const Real CtFrz = - Eos::calcCtFreezing(EosChoice, SaTop, PTopDb, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); // Heat tendencies are due to direct heat fluxes + enthalpy fluxes // The enthalpy of liquid water is assumed to be: // - local SST for liquid mass fluxes (rain, rivers) - // - local freezing point for solid --> liq mass fluxes (snow, frozen + // - zero degrees for solid --> liq mass fluxes (snow, frozen // runoff) // - solid mass fluxes are locally melted by the ocean (constant Lat // heat of fusion) + // Technically there is an enthalpy flux associated with solid water + // associated with the energy needed to bring the solid water to the + // freezing point, captured by a Cp0Sw * CtFrz term. We make the + // assumption that the incoming solid fresh water from snow and ice + // is at 0C and that this term is 0. Thus, the enthalpy flux associated + // with solid water is only captured by the latent heat of fusion term, + // LatIce. const Real HeatFlux = LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + - (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + - (SnowFlux(ICell) + IceRunoffFlux(ICell)) * - (Cp0Sw * CtFrz - LatIce); + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * + Kokkos::max(0.0_Real, CtTop) - + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * LatIce; Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; } diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 13f8ed2ba2a4..eb68daa03751 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -499,17 +499,9 @@ int testSfcTracerForcing() { ThickTimeLevel, VelTimeLevel, TracerTimeLevel, Time, Interval); - // Build a reference expectations for temperature tendency: - // using TEOS-10 freezing CT (expected to pass under strict tolerance). - - HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); - deepCopy(PressureMidH, VCoord->PressureMid); - const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; - const Real CtFrz = - Eos::calcCtFreezing(EosInst->EosChoice, SaTopValue, PTopDb, 0.0_Real); + // Build a reference expectations for temperature tendency const Real ExpectedTempTend = - (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + - TestSnow * (Cp0Sw * CtFrz - LatIce)) * + (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue - TestSnow * LatIce) * HFluxFac; // SaltTend = SeaIceSaltFlux * SFluxFac From 5772ba34715d47d0b93a95b39d9dcdf1787dd4c1 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 17 Aug 2026 08:07:01 -0700 Subject: [PATCH 24/27] added limiter and Evap Enthalpy to TendencyTerms --- components/omega/src/ocn/GlobalConstants.h | 2 ++ components/omega/src/ocn/TendencyTerms.h | 42 ++++++++++++---------- 2 files changed, 25 insertions(+), 19 deletions(-) diff --git a/components/omega/src/ocn/GlobalConstants.h b/components/omega/src/ocn/GlobalConstants.h index 0fe2adb1f54e..51a9054d6630 100644 --- a/components/omega/src/ocn/GlobalConstants.h +++ b/components/omega/src/ocn/GlobalConstants.h @@ -52,6 +52,8 @@ constexpr Real TkTrip = constexpr Real TkFrz = pcd::pure_water_freezing_temperature_reference; // Freezing point of fresh water ~ K (from Physical Constants Dictionary) constexpr Real TkFrzSw = TkFrz - 1.8; // Freezing point of seawater ~ K +constexpr Real Ct0Fw = 0.0153; +// Conservative Temperature of freshwater at 0 C (ref pressure) ~ C constexpr Real RhoAir = pcd::dry_air_density_at_standard_temperature_and_pressure; // Density of air ~ kg/m^3 (from Physical Constants Dictionary) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 8c45b7bee392..29088b032d65 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -437,27 +437,31 @@ class SfcTracerForcingOnCell { const Real CtTop = TracerCell(TempIndex, ICell, KTop); - // Heat tendencies are due to direct heat fluxes + enthalpy fluxes - // The enthalpy of liquid water is assumed to be: - // - local SST for liquid mass fluxes (rain, rivers) - // - zero degrees for solid --> liq mass fluxes (snow, frozen - // runoff) - // - solid mass fluxes are locally melted by the ocean (constant Lat - // heat of fusion) - // Technically there is an enthalpy flux associated with solid water - // associated with the energy needed to bring the solid water to the - // freezing point, captured by a Cp0Sw * CtFrz term. We make the - // assumption that the incoming solid fresh water from snow and ice - // is at 0C and that this term is 0. Thus, the enthalpy flux associated - // with solid water is only captured by the latent heat of fusion term, - // LatIce. + // CT tendencies are due to direct heat fluxes + pot enthalpy fluxes + // Each mass flux has an associated potential enthalpy flux. + // Levels of simplification can be done here. For now: + // - We approximate PotEnthalpyIce(Tinsitu, P=0) ~ -LatIce (constant); + // Altrntively, we could use cnst PotEnthalpyIce(0, 0) ​= −333360 + // J/kg a 0.1% / 340 J/kg difference with the LatIce value from pcd. + // The full expression is gsw_pot_enthalpy_ice(T, P). + const Real PotEnthalpyIce = -LatIce; + // - We assume dry snow and use the same PotEnthalpyIce for snow. + // - We assume liquid water comes in at the specific enthalpy as the + // top ocean layer, (i.e. mass flux is CT-neutral). The enthalpy is + // capped by a lower bound of pot enthalpy of freshwater at 0.0 C. + // Another choice would be to add it at the same in situ temp as ocean + // i.e. CT(Sa=0, max(0, T)). + //- We assume evaporation removes the same specific enthalpy as the + // top ocean layer; not capped, to keep the mass flux CT-neutral. + const Real PotEnthalpyFwIn = Cp0Sw * Kokkos::max(Ct0Fw, CtTop); + const Real PotEnthalpyFwout = Cp0Sw * CtTop; const Real HeatFlux = - LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + - SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + - (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * - Kokkos::max(0.0_Real, CtTop) - - (SnowFlux(ICell) + IceRunoffFlux(ICell)) * LatIce; + ShortWaveHeatFlux(ICell) + SensibleHeatFlux(ICell) + + SeaIceHeatFlux(ICell) + // includes enthalpy of meltwater already + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * PotEnthalpyFwIn + + LatentHeatFlux(ICell) + EvaporationFlux(ICell) * PotEnthalpyFwOut + + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * PotEnthalpyIce; Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; } From d1212bcc6f97fb5f7394966176a76dd19a1f0317 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 17 Aug 2026 09:32:18 -0700 Subject: [PATCH 25/27] renamed latent to specify evap and updated the docs --- components/omega/doc/devGuide/Forcing.md | 12 +++++----- components/omega/doc/userGuide/Forcing.md | 6 ++--- components/omega/src/ocn/Forcing.cpp | 2 +- components/omega/src/ocn/Tendencies.cpp | 15 ++++++------ components/omega/src/ocn/TendencyTerms.h | 24 +++++++++---------- .../src/ocn/forcingVars/TracerForcingVars.cpp | 16 ++++++------- .../src/ocn/forcingVars/TracerForcingVars.h | 2 +- components/omega/test/ocn/TendenciesTest.cpp | 4 ++-- 8 files changed, 41 insertions(+), 40 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 7e69dc302bc6..0a6cae92d1b3 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -55,7 +55,7 @@ the surface layer pseudo-thickness. **Tracer equation pathway:** 1. External fields provide heat and salt flux components: - - `LatentHeatFlux`, `SensibleHeatFlux` + - `LatentHeatFluxEvap`, `SensibleHeatFlux` - `LongWaveHeatFluxUp`, `LongWaveHeatFluxDown` - `SeaIceHeatFlux`, `ShortWaveHeatFlux` - mass fluxes which add energy changes (`SnowFlux`, `RainFlux`, `IceRunoffFlux`, `RiverRunoffFlux`) @@ -75,11 +75,11 @@ the surface layer pseudo-thickness. - 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}}$ + 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` diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index 01eeff4b5e3d..bd397e13bbeb 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -73,9 +73,9 @@ Coupled tracer flux forcing uses 13 auxiliary fields organized by type: - `IceRunoffFlux`: runoff from land ice - `RiverRunoffFlux`: runoff from rivers -**Heat fluxes (W m⁻²):** -- `LatentHeatFlux`: latent heat transfer -- `SensibleHeatFlux`: sensible heat transfer +**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) diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 69e7c502ff64..6bf0c7b9b7fd 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -203,7 +203,7 @@ void Forcing::resetArrays() { deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); - deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.LatentHeatFluxEvapCell, 0.0_Real); deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index c5d079d8e59e..9840b18d5688 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -993,8 +993,8 @@ void Tendencies::computeTracerTendenciesOnly( if (LocSfcTracerForcing.Enabled) { Pacer::start("Tend:sfcTracerForcing", 2); const auto *ForcingState = Forcing::getDefault(); - const auto &LatentHeatFlux = - ForcingState->TracerForcing.LatentHeatFluxCell; + const auto &LatentHeatFluxEvap = + ForcingState->TracerForcing.LatentHeatFluxEvapCell; const auto &SensibleHeatFlux = ForcingState->TracerForcing.SensibleHeatFluxCell; const auto &LongWaveHeatFluxUp = @@ -1016,11 +1016,12 @@ void Tendencies::computeTracerTendenciesOnly( parallelFor( {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { - LocSfcTracerForcing( - LocTracerTend, ICell, TracerArray, PressureMid, LatentHeatFlux, - SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, - SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, - IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); + LocSfcTracerForcing(LocTracerTend, ICell, TracerArray, PressureMid, + LatentHeatFluxEvap, SensibleHeatFlux, + LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, + RainFlux, IceRunoffFlux, RiverRunoffFlux, + SeaIceSaltFlux); }); Pacer::stop("Tend:sfcTracerForcing", 2); } diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 29088b032d65..9d2b76c4bb0a 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -416,17 +416,16 @@ class SfcTracerForcingOnCell { I4 TempTracerIndex, I4 SaltTracerIndex, const Eos *EosInst); - KOKKOS_FUNCTION void - operator()(const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, - const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, - const Array1DReal &SensibleHeatFlux, - const Array1DReal &LongWaveHeatFluxUp, - const Array1DReal &LongWaveHeatFluxDown, - const Array1DReal &SeaIceHeatFlux, - const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, - const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, - const Array1DReal &RiverRunoffFlux, - const Array1DReal &SeaIceSaltFlux) const { + KOKKOS_FUNCTION void operator()( + const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, + const Array2DReal &PressureMid, const Array1DReal &LatentHeatFluxEvap, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, const Array1DReal &ShortWaveHeatFlux, + const Array1DReal &SnowFlux, const Array1DReal &RainFlux, + const Array1DReal &IceRunoffFlux, const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { const I4 KTop = MinLayerCell(ICell); if (KTop > MaxLayerCell(ICell)) { @@ -460,7 +459,8 @@ class SfcTracerForcingOnCell { ShortWaveHeatFlux(ICell) + SensibleHeatFlux(ICell) + SeaIceHeatFlux(ICell) + // includes enthalpy of meltwater already (RainFlux(ICell) + RiverRunoffFlux(ICell)) * PotEnthalpyFwIn + - LatentHeatFlux(ICell) + EvaporationFlux(ICell) * PotEnthalpyFwOut + + LatentHeatFluxEvap(ICell) + + EvaporationFlux(ICell) * PotEnthalpyFwOut + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * PotEnthalpyIce; Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp index a6478bb612c4..9e8278b1e7d1 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -17,7 +17,7 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, Mesh->NCellsSize), IceRunoffFluxCell("IceRunoffFlux" + Suffix, Mesh->NCellsSize), RiverRunoffFluxCell("RiverRunoffFlux" + Suffix, Mesh->NCellsSize), - LatentHeatFluxCell("LatentHeatFlux" + Suffix, Mesh->NCellsSize), + LatentHeatFluxEvapCell("LatentHeatFluxEvap" + Suffix, Mesh->NCellsSize), SensibleHeatFluxCell("SensibleHeatFlux" + Suffix, Mesh->NCellsSize), LongWaveHeatFluxUpCell("LongWaveHeatFluxUp" + Suffix, Mesh->NCellsSize), LongWaveHeatFluxDownCell("LongWaveHeatFluxDown" + Suffix, @@ -31,7 +31,7 @@ TracerForcingVars::TracerForcingVars(const std::string &Suffix, deepCopy(SeaIceFreshWaterFluxCell, 0.0_Real); deepCopy(IceRunoffFluxCell, 0.0_Real); deepCopy(RiverRunoffFluxCell, 0.0_Real); - deepCopy(LatentHeatFluxCell, 0.0_Real); + deepCopy(LatentHeatFluxEvapCell, 0.0_Real); deepCopy(SensibleHeatFluxCell, 0.0_Real); deepCopy(LongWaveHeatFluxUpCell, 0.0_Real); deepCopy(LongWaveHeatFluxDownCell, 0.0_Real); @@ -78,9 +78,9 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { "kg m^-2 s^-1", "", std::numeric_limits::lowest(), std::numeric_limits::max(), NDims, DimNames); - auto LatentHeatFluxField = - Field::create(LatentHeatFluxCell.label(), "latent heat flux", "W m^-2", - "", std::numeric_limits::lowest(), + auto LatentHeatFluxEvapField = + Field::create(LatentHeatFluxEvapCell.label(), "latent heat flux", + "W m^-2", "", std::numeric_limits::lowest(), std::numeric_limits::max(), NDims, DimNames); auto SensibleHeatFluxField = Field::create(SensibleHeatFluxCell.label(), "sensible heat flux", @@ -114,7 +114,7 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { FieldGroup::addFieldToGroup(SeaIceFreshWaterFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(IceRunoffFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(RiverRunoffFluxCell.label(), "Forcing"); - FieldGroup::addFieldToGroup(LatentHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LatentHeatFluxEvapCell.label(), "Forcing"); FieldGroup::addFieldToGroup(SensibleHeatFluxCell.label(), "Forcing"); FieldGroup::addFieldToGroup(LongWaveHeatFluxUpCell.label(), "Forcing"); FieldGroup::addFieldToGroup(LongWaveHeatFluxDownCell.label(), "Forcing"); @@ -128,7 +128,7 @@ void TracerForcingVars::registerFields(const std::string &MeshName) const { SeaIceFreshWaterFluxField->attachData(SeaIceFreshWaterFluxCell); IceRunoffFluxField->attachData(IceRunoffFluxCell); RiverRunoffFluxField->attachData(RiverRunoffFluxCell); - LatentHeatFluxField->attachData(LatentHeatFluxCell); + LatentHeatFluxEvapField->attachData(LatentHeatFluxEvapCell); SensibleHeatFluxField->attachData(SensibleHeatFluxCell); LongWaveHeatFluxUpField->attachData(LongWaveHeatFluxUpCell); LongWaveHeatFluxDownField->attachData(LongWaveHeatFluxDownCell); @@ -144,7 +144,7 @@ void TracerForcingVars::unregisterFields() const { Field::destroy(SeaIceFreshWaterFluxCell.label()); Field::destroy(IceRunoffFluxCell.label()); Field::destroy(RiverRunoffFluxCell.label()); - Field::destroy(LatentHeatFluxCell.label()); + Field::destroy(LatentHeatFluxEvapCell.label()); Field::destroy(SensibleHeatFluxCell.label()); Field::destroy(LongWaveHeatFluxUpCell.label()); Field::destroy(LongWaveHeatFluxDownCell.label()); diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.h b/components/omega/src/ocn/forcingVars/TracerForcingVars.h index e38d9948f672..8e6f15bf3f9f 100644 --- a/components/omega/src/ocn/forcingVars/TracerForcingVars.h +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.h @@ -22,7 +22,7 @@ class TracerForcingVars { Array1DReal IceRunoffFluxCell; Array1DReal RiverRunoffFluxCell; - Array1DReal LatentHeatFluxCell; + Array1DReal LatentHeatFluxEvapCell; Array1DReal SensibleHeatFluxCell; Array1DReal LongWaveHeatFluxUpCell; Array1DReal LongWaveHeatFluxDownCell; diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index eb68daa03751..1b61c0266687 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -393,7 +393,7 @@ int testSfcTracerForcing() { // Retrieve forcing field views auto &SensibleHeatFlux = DefForcing->TracerForcing.SensibleHeatFluxCell; - auto &LatentHeatFlux = DefForcing->TracerForcing.LatentHeatFluxCell; + auto &LatentHeatFluxEvap = DefForcing->TracerForcing.LatentHeatFluxEvapCell; auto &LongWaveHeatFluxUp = DefForcing->TracerForcing.LongWaveHeatFluxUpCell; auto &LongWaveHeatFluxDown = DefForcing->TracerForcing.LongWaveHeatFluxDownCell; @@ -407,7 +407,7 @@ int testSfcTracerForcing() { // Initialize all fluxes to zero deepCopy(SensibleHeatFlux, 0._Real); - deepCopy(LatentHeatFlux, 0._Real); + deepCopy(LatentHeatFluxEvap, 0._Real); deepCopy(LongWaveHeatFluxUp, 0._Real); deepCopy(LongWaveHeatFluxDown, 0._Real); deepCopy(SeaIceHeatFlux, 0._Real); From 85cff877ff77552240904ef5eb7fd4fb712f3539 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 17 Aug 2026 09:52:49 -0700 Subject: [PATCH 26/27] added missing arguments --- components/omega/src/ocn/Tendencies.cpp | 8 +++++--- components/omega/src/ocn/TendencyTerms.h | 3 ++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 9840b18d5688..cf5ee9e1823e 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -1005,8 +1005,10 @@ void Tendencies::computeTracerTendenciesOnly( ForcingState->TracerForcing.SeaIceHeatFluxCell; const auto &ShortWaveHeatFlux = ForcingState->TracerForcing.ShortWaveHeatFluxCell; - const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; - const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &EvaporationFlux = + ForcingState->TracerForcing.EvaporationFluxCell; const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; const auto &RiverRunoffFlux = ForcingState->TracerForcing.RiverRunoffFluxCell; @@ -1021,7 +1023,7 @@ void Tendencies::computeTracerTendenciesOnly( LongWaveHeatFluxUp, LongWaveHeatFluxDown, SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, IceRunoffFlux, RiverRunoffFlux, - SeaIceSaltFlux); + EvaporationFlux, SeaIceSaltFlux); }); Pacer::stop("Tend:sfcTracerForcing", 2); } diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 9d2b76c4bb0a..30d3e609bd1e 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -425,6 +425,7 @@ class SfcTracerForcingOnCell { const Array1DReal &SeaIceHeatFlux, const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, const Array1DReal &RiverRunoffFlux, + const Array1DReal &EvaporationFlux, const Array1DReal &SeaIceSaltFlux) const { const I4 KTop = MinLayerCell(ICell); @@ -453,7 +454,7 @@ class SfcTracerForcingOnCell { //- We assume evaporation removes the same specific enthalpy as the // top ocean layer; not capped, to keep the mass flux CT-neutral. const Real PotEnthalpyFwIn = Cp0Sw * Kokkos::max(Ct0Fw, CtTop); - const Real PotEnthalpyFwout = Cp0Sw * CtTop; + const Real PotEnthalpyFwOut = Cp0Sw * CtTop; const Real HeatFlux = LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + ShortWaveHeatFlux(ICell) + SensibleHeatFlux(ICell) + From 808185a4c96677e40f264fe46a7800ab6c4a1ed4 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Mon, 17 Aug 2026 18:41:08 -0400 Subject: [PATCH 27/27] adds non-teos option to limiter --- components/omega/src/ocn/TendencyTerms.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 30d3e609bd1e..9253026d75f2 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -453,7 +453,9 @@ class SfcTracerForcingOnCell { // i.e. CT(Sa=0, max(0, T)). //- We assume evaporation removes the same specific enthalpy as the // top ocean layer; not capped, to keep the mass flux CT-neutral. - const Real PotEnthalpyFwIn = Cp0Sw * Kokkos::max(Ct0Fw, CtTop); + const Real CtLim = + (EosChoice == EosType::Teos10Eos) ? Ct0Fw : 0.0_Real; + const Real PotEnthalpyFwIn = Cp0Sw * Kokkos::max(CtLim, CtTop); const Real PotEnthalpyFwOut = Cp0Sw * CtTop; const Real HeatFlux = LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +