From ebc9f7552e06ff0d25946ab957e1953d7c83f8a5 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 19 Jun 2026 14:18:42 -0700 Subject: [PATCH 01/23] 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 42643dba5385..2c51af33ae7d 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -29,7 +29,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(); } @@ -37,10 +37,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, @@ -162,6 +166,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 d8273a844de0..562badac7437 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -283,6 +283,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"); @@ -461,6 +473,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), @@ -510,6 +524,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); @@ -553,6 +568,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, @@ -772,6 +813,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); @@ -893,6 +935,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 877ba22faf2b..c60a2783ddb3 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 50b4c12b3eda..9cf3737da573 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -254,6 +254,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 a39f44d09f6be49b526324097c1fa462ae2fa5ff Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 09:18:03 -0700 Subject: [PATCH 02/23] 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 16ea4da9f5a2..516cb629e154 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -343,6 +343,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 5e4fd89cae21..2b3d6d78f462 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -761,6 +761,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 562badac7437..419c105292a6 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -474,7 +474,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), @@ -952,16 +953,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 52e66c42d9ca83aafdf7e62ff2f70c898c402436 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:19:13 -0700 Subject: [PATCH 03/23] 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 01b9a4e688e059fcbae81a6bc34a81a75eb105ba Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:46:07 -0700 Subject: [PATCH 04/23] 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 419c105292a6..3533cbc5aa0f 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -959,15 +959,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 69889f5792de175471d8abc311752dc4978ba561 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:27:46 -0700 Subject: [PATCH 05/23] 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 516cb629e154..49ea504fb84c 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -352,7 +352,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 8f637a493245d01f73c254736c3668933ec0cfa8 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:49:55 -0700 Subject: [PATCH 06/23] 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 3f5be9e475c149e7ced14c79160a0d5061e42c16 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:07:40 -0700 Subject: [PATCH 07/23] 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 3533cbc5aa0f..419c105292a6 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -959,17 +959,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 0ab11b631e45750bacd5b772e19c4f125c928f60 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:53:52 -0700 Subject: [PATCH 08/23] 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 83768fbd01e89fec563fe372350bb9b2aedfd858 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 11:33:03 -0700 Subject: [PATCH 09/23] 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 b62da63286fc21ee84392952bac599488f7332e3 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 13:00:31 -0700 Subject: [PATCH 10/23] 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 2c51af33ae7d..95ba94b3f305 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -157,7 +157,8 @@ void Forcing::computeSfcStressForcingOnEdge() const { Pacer::stop("Forcing:edge1", 2); } -// Exchange halo for surface stress cell fields. +// Exchange halo for surface stress cell fields. Only needed for variables that +// need information beyond cell-centered values. I4 Forcing::exchangeHalo() const { I4 Err = 0; @@ -166,33 +167,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 cb37ca54de078f1121ea4a06cd9931457e73a4a2 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 14:41:04 -0700 Subject: [PATCH 11/23] 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 0402dbeb27342aaf386908a11ecb60a7aac05caf Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Thu, 16 Jul 2026 13:08:27 -0700 Subject: [PATCH 12/23] 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 49ea504fb84c..9a3428f537f7 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,14 +327,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 2b3d6d78f462..0911910a79f4 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -756,7 +756,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 6ec6568a993c10d66ffa74521e16d935f6069607 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 20 Jul 2026 09:38:07 -0700 Subject: [PATCH 13/23] 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 9a3428f537f7..ba36a2bb61ac 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,28 +327,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 0911910a79f4..61d5171a2031 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -756,6 +756,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); @@ -763,14 +766,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 ae341ebb3e6e1a8fc908d3c58adfa95c8b718847 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 20 Jul 2026 09:48:41 -0700 Subject: [PATCH 14/23] 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 5b1b932f5dc5934fc43102971df66404bb3e0436 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 20 Jul 2026 21:26:40 -0700 Subject: [PATCH 15/23] 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 95ba94b3f305..e316456ddbd7 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -143,7 +143,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 { @@ -166,6 +191,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; } @@ -177,13 +226,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 162b0074f4f66c40e20fd0ddcd1df359274fbab9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 22 Jul 2026 17:57:57 -0600 Subject: [PATCH 16/23] 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 173f735027011a53374db8ef1a19f017b231a925 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 22 Jul 2026 18:03:52 -0600 Subject: [PATCH 17/23] 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 d840a9da4b6b5bae19a89cd2efd3e203e76640a3 Mon Sep 17 00:00:00 2001 From: Katherine Smith Date: Thu, 23 Jul 2026 19:26:08 -0400 Subject: [PATCH 18/23] 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 e316456ddbd7..5b7fef628433 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -167,7 +167,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 3e855db968d8349be3445706a2e992ab20ce1038 Mon Sep 17 00:00:00 2001 From: Katherine Smith Date: Fri, 24 Jul 2026 01:18:27 -0400 Subject: [PATCH 19/23] 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 ba36a2bb61ac..6dc9b444be31 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,6 +327,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 61d5171a2031..b77a9b60cb89 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -777,26 +777,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 409ff6e3b68c204243187fd345cdb1d747eda7c3 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Tue, 28 Jul 2026 17:47:10 -0400 Subject: [PATCH 20/23] 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 6dc9b444be31..ba36a2bb61ac 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,29 +327,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 b77a9b60cb89..61d5171a2031 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -777,8 +777,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 fd25d6b17c9416ea86870e2b35345624318c1716 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Thu, 30 Jul 2026 16:35:31 -0400 Subject: [PATCH 21/23] 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 61d5171a2031..2494d6a2d5d8 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -358,8 +358,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; @@ -777,17 +777,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 0e7eb2c5be76..f8ca084d4559 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -706,6 +706,42 @@ void testBruntVaisalaFreqSqTeos10() { 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(); @@ -767,22 +803,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; @@ -837,7 +873,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); - checkValueCtFreezing(); + checkValueGswcCtFreezing(); checkValueGswcCtFromPt(); checkValueGswcPtFromCt(); @@ -848,6 +884,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { testEosTeos10(); testEosTeos10Displaced(); testBruntVaisalaFreqSqTeos10(); + 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 c44f365d088299a5db0d86afcdff48183d176f66 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Fri, 31 Jul 2026 15:00:26 -0400 Subject: [PATCH 22/23] 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 5b7fef628433..46af26094113 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -36,14 +36,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. @@ -80,10 +88,9 @@ void Forcing::init() { ABORT_ERROR("Forcing: failed to initialize default forcing state"); } - DefaultForcing->registerFields(DefMesh->MeshName); - Config *OmegaConfig = Config::getOmegaConfig(); 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. @@ -140,6 +147,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). @@ -150,23 +181,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. @@ -186,34 +221,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; } @@ -227,6 +240,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 8f0c28e384b15c4dfa6642b6a6c69783a90ae5a0 Mon Sep 17 00:00:00 2001 From: katsmith133 Date: Mon, 10 Aug 2026 17:28:33 -0400 Subject: [PATCH 23/23] adds basic frazil option --- components/omega/configs/Default.yml | 6 + components/omega/doc/devGuide/Frazil.md | 93 ++++ components/omega/doc/devGuide/Tendencies.md | 16 + components/omega/doc/index.md | 2 + components/omega/doc/userGuide/Frazil.md | 63 +++ components/omega/doc/userGuide/Tendencies.md | 6 + components/omega/src/ocn/Frazil.cpp | 369 ++++++++++++++ components/omega/src/ocn/Frazil.h | 207 ++++++++ components/omega/src/ocn/OceanFinal.cpp | 2 + components/omega/src/ocn/OceanInit.cpp | 2 + components/omega/src/ocn/Tendencies.cpp | 15 +- components/omega/src/ocn/Tendencies.h | 1 + components/omega/src/ocn/TendencyTerms.cpp | 68 +++ components/omega/src/ocn/TendencyTerms.h | 22 + components/omega/test/CMakeLists.txt | 7 + components/omega/test/infra/IOStreamTest.cpp | 7 +- components/omega/test/ocn/FrazilTest.cpp | 497 +++++++++++++++++++ components/omega/test/ocn/TendenciesTest.cpp | 151 +++++- 18 files changed, 1527 insertions(+), 7 deletions(-) create mode 100644 components/omega/doc/devGuide/Frazil.md create mode 100644 components/omega/doc/userGuide/Frazil.md create mode 100644 components/omega/src/ocn/Frazil.cpp create mode 100644 components/omega/src/ocn/Frazil.h create mode 100644 components/omega/test/ocn/FrazilTest.cpp diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index e33403c7cc91..971dc4f15ad3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -71,6 +71,7 @@ Omega: VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true PressureGradTendencyEnable: true + FrazilTendencyEnable: false VelVertMixTendencyEnable: true TracerVertMixTendencyEnable: true ManufacturedSolution: @@ -86,6 +87,11 @@ Omega: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + Frazil: + FrazilType: basic + MassLimit: 0.1 + DepthLimit: -1.0 + ConservationCheck: false VertMix: Background: Diffusivity: 1.0e-5 diff --git a/components/omega/doc/devGuide/Frazil.md b/components/omega/doc/devGuide/Frazil.md new file mode 100644 index 000000000000..4186324be612 --- /dev/null +++ b/components/omega/doc/devGuide/Frazil.md @@ -0,0 +1,93 @@ +(omega-dev-frazil)= + +# Frazil + +This page describes frazil design and implementation details in Omega, +for the current basic pathway. A TEOS-10 option will be added soon. + +## Purpose and coupling points + +Frazil computes phase-change-related tendencies that modify: + +- pseudo-thickness tendency +- temperature tracer tendency +- salinity tracer tendency + +The tendency hook-up is implemented through `FrazilOnCell` in the tracer +tendency phase, where frazil contributions are added to accumulated +`PseudoThicknessTend` and `TracerTend` arrays. + +## Data flow and call sequence + +1. `Tendencies.computeTracerTendenciesOnly` checks + `Tendencies.FrazilTendencyEnable`. +2. If enabled, `FrazilOnCell.operator()` retrieves the default Frazil object + and zeros frazil tendency and accumulator arrays. +3. `FrazilOnCell` extracts `Temperature` and `Salinity` tracer subviews, then calls + `Frazil.computeFrazil(CT, SA, PressureMid, PseudoThickness)`. +4. `Frazil.computeFrazil` dispatches using FrazilType. The only current supported + option is basic. A TEOS-10 option will be added soon. +5. Returned frazil tendencies are added into `PseudoThicknessTend` and + `TracerTend` for active cell layers. + +## Configuration coupling + +Frazil behavior is configured with: + +- `Omega.Tendencies.FrazilTendencyEnable` + - Global switch for applying frazil tendency terms. +- `Omega.Frazil.FrazilType` + - Implementation choice. Only current supported value: `basic`. +- `Omega.Frazil.MassLimit` + - Per-layer mass and thickness limiter used by basic formation and melt. +- `Omega.Frazil.DepthLimit` + - Optional depth cutoff for frazil activity. Negative means no cutoff. +- `Omega.Frazil.ConservationCheck` + - Optional post-compute column conservation diagnostic logging. + +## Basic pathway summary + +- Freezing-point checks are based on conservative temperature and absolute + salinity with pressure-dependent freezing temperature. +- Vertical accumulation order is bottom-to-top within each active column. +- Frazil tendencies are not time-step scaled inside `Frazil`; they are + accumulated as tendency contributions. +- Surface salt redistribution is applied at the top active layer. +- Column accumulators are converted to coupler units at the end of the + per-column loop. + +### Basic pathway (`FrazilType: basic`) + +- Uses simplified energetics: the energy of the super-cooled water sets the +amount of solid ice formed (used constant latent heat of fusion of fresh ice). +Salt is added based on a constant bulk salinity `IceRefSal` (default) or a +manual toggle (for now) using the local salinity and the frazil porosity. +Melting of existing frazil is set by the amount of pure ice that can be melted +by the warm layer, and the enthalpy of melt water at the local freezing point +is added to the temperature tendency. +- Computes local layer tendencies (`HTend`, `TTend`, `STend`) and updates + accumulated frazil stores. +- Applies surface salt redistribution adjustment and converts accumulators to + coupler units at the end of the column loop. +- Warnings: 1) basic frazil formation/melt is does not conserve energy. + 2) Using porosity to set the salt content includes a redistribution of excess + salt in the surface layer, which can be very significant. + +## Existing ctest coverage + +The existing frazil test driver is in +`components/omega/test/ocn/FrazilTest.cpp` and covers: + +- basic frazil formation in cold and warm single-layer states +- mixed warm/cold column behavior with sign checks for branch switching +- depth-limit behavior ensuring excluded deep layers have zero frazil tendency + +## Extensibility + +The frazil choice dispatch and configuration plumbing are retained so future +frazil options can be added without changing Tendencies call sites. + +## Related pages + +- User-facing options: [User Frazil Guide](../userGuide/Frazil.md) +- Tendency container and hook-up: [Tendencies](Tendencies.md) diff --git a/components/omega/doc/devGuide/Tendencies.md b/components/omega/doc/devGuide/Tendencies.md index 443e8e6ddcd5..32ce44642451 100644 --- a/components/omega/doc/devGuide/Tendencies.md +++ b/components/omega/doc/devGuide/Tendencies.md @@ -64,6 +64,22 @@ Tendencies.computeVelocityTendencies(State, AuxState, ThickTimeLevel, VelTimeLev To call only the tracer tendency terms: Tendencies.computeTracerTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel); +## Frazil tendencies + +Frazil tendencies are calculated by the Frazil class, not the TendencyTerms. +The tendencies are integrated in the tracer tendency compute path and read in from +the frazil computation. When +`Omega.Tendencies.FrazilTendencyEnable` is true, +`computeTracerTendenciesOnly` invokes `FrazilOnCell`, which calls the default +`Frazil` object and adds frazil contributions to: + +- `PseudoThicknessTend` +- `TracerTend` for `Temperature` +- `TracerTend` for `Salinity` + +The frazil implementation choice and parameters are read from the +`Omega.Frazil` block. Detailed algorithm notes are in [Frazil](Frazil.md). + ## Removal of tendencies To erase a specific named tendencies instance use `erase` ```c++ diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index dfc81a284670..30f0a777ba22 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -54,6 +54,7 @@ userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv userGuide/Forcing +userGuide/Frazil userGuide/SfcCoupling userGuide/Analysis ``` @@ -102,6 +103,7 @@ devGuide/Timing devGuide/VerticalMixingCoeff devGuide/VertAdv devGuide/Forcing +devGuide/Frazil devGuide/SfcCoupling devGuide/Analysis ``` diff --git a/components/omega/doc/userGuide/Frazil.md b/components/omega/doc/userGuide/Frazil.md new file mode 100644 index 000000000000..e4078a8ca141 --- /dev/null +++ b/components/omega/doc/userGuide/Frazil.md @@ -0,0 +1,63 @@ +(omega-user-frazil)= + +# Frazil + +This page describes user-facing configuration for frazil tendencies in Omega. +Frazil physics represents the formation and melt of frazil ice within +the ocean water column. It impacts local layer pseudo-thickness (i.e. mass), +temperature, and salinity tendencies. The vertical sum of the frazil energy, +mass of water and mass of salt are passed to the coupler (if coupled) or +discarded (in ocean standalone mode). + +## Configuration overview + +Frazil behavior is controlled by one enable switch in Tendencies and one +Frazil configuration block: + +```yaml +Omega: + Tendencies: + FrazilTendencyEnable: true + + Frazil: + FrazilType: basic + MassLimit: 0.1 + DepthLimit: -1.0 + ConservationCheck: false +``` + +- `Tendencies.FrazilTendencyEnable` + - Enables/disables application of frazil tendency contributions. +- `Frazil.FrazilType` + - Selects frazil option. + - Current supported option: `basic`. +- `Frazil.MassLimit` + - limits per-layer frazil mass/thickness tendency magnitude (applied + to formation and melt) +- `Frazil.DepthLimit` + - Limits depth range where frazil is computed. + - Negative values mean no depth limit. +- `Frazil.ConservationCheck` + - Enables a column-level diagnostic conservation check with logging. + +## Available frazil option + +Omega currently supports one active frazil pathway: + +- `basic` + - freezing is based on the formation of fresh solid ice, to which + salt is added (similar the mpas-ocean implementation). + +This pathway contributes to: + +- pseudo-thickness tendency +- temperature tracer tendency +- salinity tracer tendency + +## Notes + +- Frazil tendencies are applied through the `Tendencies` tracer-step workflow. +- The frazil tendency hook assumes tracer names include `Temperature` and + `Salinity`. +- For implementation and algorithm details, see + [Developer Frazil Guide](../devGuide/Frazil.md). diff --git a/components/omega/doc/userGuide/Tendencies.md b/components/omega/doc/userGuide/Tendencies.md index 6fc1803b176d..00597e75245c 100644 --- a/components/omega/doc/userGuide/Tendencies.md +++ b/components/omega/doc/userGuide/Tendencies.md @@ -6,3 +6,9 @@ The `Tendencies` class provides a container for the [tendency terms](#omega-user Upon creation of an `Tendencies` instance, these functors are initialized and arrays for the accumulated tendencies are allocated. There are no user-configurable options beyond those for the tendency term functors. + +Frazil is one of the configurable tendency contributions. It is controlled by +`Omega.Tendencies.FrazilTendencyEnable` together with the `Omega.Frazil` +configuration block (`FrazilType`, `MassLimit`, `DepthLimit`, and +`ConservationCheck`). For operational guidance and options, see +[Frazil](Frazil.md). diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp new file mode 100644 index 000000000000..c0e968adb0c9 --- /dev/null +++ b/components/omega/src/ocn/Frazil.cpp @@ -0,0 +1,369 @@ +//===-- ocn/Frazil.cpp - Frazil Ice Formation -----------------*- C++ -*-===// +// +// The Frazil class manages frazil tendencies and accumulators. +// This initial implementation only has a basic configuration. +// +//===----------------------------------------------------------------------===// + +#include "Frazil.h" +#include "Config.h" +#include "Eos.h" +#include "Error.h" +#include "Logging.h" + +namespace OMEGA { + +namespace { + +KOKKOS_INLINE_FUNCTION bool isApprox(const Real X, const Real Y, + const Real RTol, const Real ATol = 0) { + if (Kokkos::isnan(X) || Kokkos::isnan(Y) || Kokkos::isinf(X) || + Kokkos::isinf(Y)) { + return false; + } + + return Kokkos::abs(X - Y) <= + Kokkos::max(ATol, RTol * Kokkos::max(Kokkos::abs(X), Kokkos::abs(Y))); +} + +} // namespace + +Frazil *Frazil::DefaultFrazil = nullptr; +std::map> Frazil::AllFrazil; + +/// Constructor for BasicFrazilFormation +BasicFrazilFormation::BasicFrazilFormation() {} + +/// Constructor for BasicFrazilMelt +BasicFrazilMelt::BasicFrazilMelt() {} + +void Frazil::init() { + + if (!HorzMesh::getDefault() || !VertCoord::getDefault()) { + ABORT_ERROR("Frazil::init: HorzMesh and VertCoord must be initialized"); + } + + // Frazil freezing-temperature calculations depend on EOS configuration. + Eos::init(); + + if (!DefaultFrazil) { + Error Err; + bool FrazilTendencyEnable = false; + Config *OmegaConfig = Config::getOmegaConfig(); + Config TendConfig("Tendencies"); + + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, + "Frazil::init: Tendencies group not found in Config"); + + Err += TendConfig.get("FrazilTendencyEnable", FrazilTendencyEnable); + CHECK_ERROR_ABORT( + Err, "Frazil::init: FrazilTendencyEnable not found in Tendencies"); + + if (!FrazilTendencyEnable) { + LOG_INFO("Frazil::init: Frazil tendency disabled; skipping default " + "frazil object creation"); + LOG_INFO("All frazil is off - frazil parameters will be ignored"); + return; + } + + DefaultFrazil = create("Default"); + } +} + +Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) + : frazilChoice(FrazilType::BasicFrazil), computeBasicFrazilFormation(), + computeBasicFrazilMelt(), NCellsAll(Mesh->NCellsAll), + NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh), + VCoordPtr(VCoord) { + + FrazilTTend = + Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilSTend = + Array2DReal("FrazilSTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilHTend = + Array2DReal("FrazilHTend", Mesh->NCellsSize, VCoord->NVertLayers); + + AccMIce = Array1DReal("AccMIce", Mesh->NCellsSize); + AccEIce = Array1DReal("AccEIce", Mesh->NCellsSize); + AccMLiq = Array1DReal("AccMLiq", Mesh->NCellsSize); + AccELiq = Array1DReal("AccELiq", Mesh->NCellsSize); + AccMSalt = Array1DReal("AccMSalt", Mesh->NCellsSize); + + deepCopy(FrazilTTend, 0.0_Real); + deepCopy(FrazilSTend, 0.0_Real); + deepCopy(FrazilHTend, 0.0_Real); + deepCopy(AccMIce, 0.0_Real); + deepCopy(AccEIce, 0.0_Real); + deepCopy(AccMLiq, 0.0_Real); + deepCopy(AccELiq, 0.0_Real); + deepCopy(AccMSalt, 0.0_Real); +} + +Frazil::~Frazil() {} + +Frazil *Frazil::create(const std::string &Name) { + if (AllFrazil.find(Name) != AllFrazil.end()) { + LOG_ERROR("Attempted to create Frazil {} but it already exists", Name); + return nullptr; + } + + auto *NewFrazil = + new Frazil(HorzMesh::getDefault(), VertCoord::getDefault()); + AllFrazil.emplace(Name, NewFrazil); + + Error Err; + Config *OmegaConfig = Config::getOmegaConfig(); + Config FrazilConfig("Frazil"); + Err += OmegaConfig->get(FrazilConfig); + CHECK_ERROR_ABORT(Err, "Frazil::create: Frazil group not found in Config"); + + std::string FrazilTypeStr; + Err += FrazilConfig.get("FrazilType", FrazilTypeStr); + CHECK_ERROR_ABORT(Err, + "Frazil::create: FrazilType not found in Frazil config"); + + if (((FrazilTypeStr == "Basic") || (FrazilTypeStr == "basic") || + (FrazilTypeStr == "BasicFrazil"))) { + NewFrazil->frazilChoice = FrazilType::BasicFrazil; + } else { + ABORT_ERROR( + "Frazil::create: Only basic frazil is supported at the moment"); + } + + Err += FrazilConfig.get("MassLimit", + NewFrazil->computeBasicFrazilFormation.massLimit); + Err += FrazilConfig.get("MassLimit", + NewFrazil->computeBasicFrazilMelt.massLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: MassLimit not found in Frazil config"); + + NewFrazil->computeBasicFrazilMelt.massLimit = + NewFrazil->computeBasicFrazilFormation.massLimit; + + Err += FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); + CHECK_ERROR_ABORT( + Err, "Frazil::create: ConservationCheck not found in Frazil config"); + + Err += FrazilConfig.get("DepthLimit", NewFrazil->depthLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: DepthLimit not found in Frazil config"); + + if (Name == "Default") { + DefaultFrazil = NewFrazil; + } + + return NewFrazil; +} + +Frazil *Frazil::getDefault() { return DefaultFrazil; } + +Frazil *Frazil::get(const std::string &Name) { + auto it = AllFrazil.find(Name); + if (it != AllFrazil.end()) { + return it->second.get(); + } + + LOG_ERROR("Frazil::get: Attempted to retrieve non-existent Frazil {}", Name); + return nullptr; +} + +void Frazil::erase(std::string InName) { + auto *ToErase = get(InName); + AllFrazil.erase(InName); + if (ToErase == DefaultFrazil) { + DefaultFrazil = nullptr; + } +} + +void Frazil::clear() { + AllFrazil.clear(); + DefaultFrazil = nullptr; +} + +void Frazil::checkColumnConservation() const { + auto MinLayerCellH = createHostMirrorCopy(VCoordPtr->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoordPtr->MaxLayerCell); + auto FrazilHTendH = createHostMirrorCopy(FrazilHTend); + auto FrazilTTendH = createHostMirrorCopy(FrazilTTend); + auto FrazilSTendH = createHostMirrorCopy(FrazilSTend); + auto AccMIceH = createHostMirrorCopy(AccMIce); + auto AccMLiqH = createHostMirrorCopy(AccMLiq); + auto AccMSaltH = createHostMirrorCopy(AccMSalt); + auto AccELiqH = createHostMirrorCopy(AccELiq); + auto AccEIceH = createHostMirrorCopy(AccEIce); + + constexpr Real RTol = 1.0e-10_Real; + + for (I4 ICell = 0; ICell < NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + Real MassTend = 0.0_Real; + Real EnergyTend = 0.0_Real; + Real SaltTend = 0.0_Real; + + for (I4 K = KMin; K <= KMax; ++K) { + MassTend += FrazilHTendH(ICell, K); + EnergyTend += FrazilTTendH(ICell, K); + SaltTend += FrazilSTendH(ICell, K); + } + + const Real MassTotal = AccMIceH(ICell) + AccMLiqH(ICell); + const Real EnergyTotal = AccELiqH(ICell) + AccEIceH(ICell); + const Real SaltTotal = AccMSaltH(ICell); + + if (ICell == 0) { + LOG_INFO("Frazil column conservation check: cell {} MassTend={} " + "MassTotal={} " + "EnergyTend={} EnergyTotal={} SaltTend={} SaltTotal={}", + ICell, MassTend * RhoSw, MassTotal, + EnergyTend * Cp0Sw * RhoSw, EnergyTotal, + SaltTend * RhoSw * PPt2Salt, SaltTotal); + LOG_INFO("Frazil column conservation check: cell {} EpsMass={} " + "EpsE={} EpsS={} ", + ICell, MassTend * RhoSw + MassTotal, + EnergyTend * Cp0Sw * RhoSw + EnergyTotal, + SaltTend * RhoSw * PPt2Salt + SaltTotal); + } + + if (!isApprox(-MassTend * RhoSw, MassTotal, RTol)) { + LOG_INFO( + "Frazil column mass check failed: cell {} tendency={} total={}", + ICell, -MassTend * RhoSw, MassTotal); + } + if (!isApprox(-EnergyTend * Cp0Sw * RhoSw, EnergyTotal, RTol)) { + LOG_INFO( + "Frazil column energy check failed: cell {} tendency={} total={}", + ICell, -EnergyTend * Cp0Sw * RhoSw, EnergyTotal); + } + if (!isApprox(-SaltTend * RhoSw * PPt2Salt, SaltTotal, RTol)) { + LOG_INFO( + "Frazil column salt check failed: cell {} tendency={} total={}", + ICell, -SaltTend * RhoSw * PPt2Salt, SaltTotal); + } + } +} + +void Frazil::computeFrazilBasicImpl(const Array2DReal &CT, + const Array2DReal &SA, const Array2DReal &P, + const Array2DReal &LayerH) { + const EosType LocEosChoice = Eos::getInstance()->EosChoice; + const Real LocDepthLimit = depthLimit; + + OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell); + OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid); + + OMEGA_SCOPE(LocComputeBasicFrazilFormation, computeBasicFrazilFormation); + OMEGA_SCOPE(LocComputeBasicFrazilMelt, computeBasicFrazilMelt); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocAccMIce, AccMIce); + OMEGA_SCOPE(LocAccEIce, AccEIce); + OMEGA_SCOPE(LocAccMLiq, AccMLiq); + OMEGA_SCOPE(LocAccELiq, AccELiq); + OMEGA_SCOPE(LocAccMSalt, AccMSalt); + OMEGA_SCOPE(LocIceRefSal, IceRefSal); + OMEGA_SCOPE(LocLatIce, LatIce); + + parallelFor( + {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + I4 Klim = KMax; + bool HasKlim = true; + const bool Limit = (LocDepthLimit >= 0.0_Real); + + if (Limit) { + HasKlim = false; + for (I4 K = KMax; K >= KMin; --K) { + if (Kokkos::abs(LocGeomZMid(ICell, K)) <= LocDepthLimit) { + Klim = K; + HasKlim = true; + break; + } + } + } + + // Explicit accumulation order: bottom layer to top layer. + for (I4 K = KMax; K >= KMin; --K) { + if (!HasKlim || K > Klim) { + LocFrazilHTend(ICell, K) = 0.0_Real; + LocFrazilTTend(ICell, K) = 0.0_Real; + LocFrazilSTend(ICell, K) = 0.0_Real; + continue; + } + + const Real SAIn = SA(ICell, K); + const Real CTIn = CT(ICell, K); + const Real PIn = P(ICell, K); + const Real PDb = PIn * Pa2Db; + const Real H = LayerH(ICell, K); + + const Real Tfrz = + Eos::calcCtFreezing(LocEosChoice, SAIn, PDb, 0.0_Real); + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + if (CTIn < Tfrz) { + LocComputeBasicFrazilFormation( + SAIn, CTIn, PDb, H, LocAccMIce(ICell), LocAccMSalt(ICell), + LocAccEIce(ICell), HTend, TTend, STend, Tfrz); + } else if (LocAccMIce(ICell) > 0.0_Real) { + LocComputeBasicFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMSalt(ICell), LocAccEIce(ICell), + HTend, TTend, STend, Tfrz); + } + + LocFrazilHTend(ICell, K) = HTend; + LocFrazilTTend(ICell, K) = TTend; + LocFrazilSTend(ICell, K) = STend; + } + + // Redistribute excess salt at the surface. No treatment of low + // salinity frazil for now. + LocFrazilSTend(ICell, KMin) += Kokkos::max( + 0.0_Real, LocAccMIce(ICell) * LocIceRefSal - LocAccMSalt(ICell)); + // hijack total terms before the coupling + LocAccMSalt(ICell) = LocAccMIce(ICell) * LocIceRefSal; + LocAccEIce(ICell) = -LocAccMIce(ICell) * LocLatIce; + + // Convert to coupler units + LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; + LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; + LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; + LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; + LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; + }); +} + +void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &LayerH) { + Eos *DefEos = Eos::getInstance(); + if (!DefEos) { + ABORT_ERROR("Frazil::computeFrazil: Eos must be initialized before " + "computeFrazil"); + } + + switch (frazilChoice) { + case FrazilType::BasicFrazil: + computeFrazilBasicImpl(CT, SA, P, LayerH); + break; + // TODO: add additional frazil choices here as new implementations are + // introduced + default: + ABORT_ERROR("Frazil::computeFrazil: Unsupported frazilChoice"); + break; + } + + if (conservationCheck) { + checkColumnConservation(); + } +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h new file mode 100644 index 000000000000..c5e7f11e1b40 --- /dev/null +++ b/components/omega/src/ocn/Frazil.h @@ -0,0 +1,207 @@ +#ifndef OMEGA_FRAZIL_H +#define OMEGA_FRAZIL_H +//===-- ocn/Frazil.h - Frazil Ice Formation -------------------*- C++ -*-===// +// +// The Frazil class manages frazil tendencies and accumulators. +// This initial implementation only has a basic configuration. +// but carries scaffolding for other implementations. +// +//===----------------------------------------------------------------------===// + +#include "GlobalConstants.h" +#include "HorzMesh.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include +#include +#include + +namespace OMEGA { + +enum class FrazilType { + BasicFrazil ///< MPAS-O style basic frazil option + // Additional options can be added here in the future. +}; + +class BasicFrazilFormation { + public: + BasicFrazilFormation(); + + Real massLimit = 0.1_Real; // to do: remove default + Real FrazilIceSalinity = IceRefSal; // Global constant + Real LatFrazil = LatIce; // Global constant + + Real FrazilPorosity = + 1.0_Real; // Internal for now; can move to config later. + + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb, + const Real H, Real &SumIceThickness, + Real &SumSalt, Real &SumEnergy, Real &HTend, + Real &TTend, Real &STend, + const Real Tfrz) const { + + const Real potential = H * Cp0Sw * RhoSw * (CT - Tfrz); + const Real freezingEnergy = Kokkos::max(0.0_Real, -potential); + + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + + Real newFrzThickness = + freezingEnergy / + (LatFrazil * RhoSw); // frazil (ice) mass in pseudo-thickness terms + newFrzThickness = Kokkos::min(newFrzThickness, H * massLimit); + const Real newFrzEnergy = + -newFrzThickness * LatFrazil; // (<0; enthalpy of ice) + + // MANUAL TOGGLE: uncomment line below to use porosity + // const Real FrazilIceSalinity = FrazilPorosity * SA; + + const Real frazilSalinity = Kokkos::min(FrazilIceSalinity, SA); + const Real newSaltContent = + newFrzThickness * frazilSalinity; // in m.(g/kg) + + HTend = -newFrzThickness - newSaltContent; + // TTend below should include the enthalpy associated with the mass flux + // this is intentionally not added here to match the mpas-o implementation + TTend = + -(newFrzEnergy) / Cp0Sw; // (E< 0 so TTend>0) // scaled to h.CT tend + STend = -newSaltContent; + + SumIceThickness += newFrzThickness; + // non-conservation between Sum and HTend by construction in the original + SumSalt += newSaltContent; + SumEnergy += newFrzEnergy; + } +}; + +// TO-DO: is there a new for salt reconciliation at the surface? +// frazilSalinityTendency(minLevelCell(iCell),iCell) = +// frazilSalinityTendency(minLevelCell(iCell),iCell) + & +// max(0.0_RKIND,(sumNewThicknessWeightedSaltContent - +// newThicknessWeightedSaltContent) ) / dt +// accumulatedFrazilIceSalinityNew(iCell) = +// accumulatedFrazilIceSalinityOld(iCell) + newThicknessWeightedSaltContent + +class BasicFrazilMelt { + public: + BasicFrazilMelt(); + + Real massLimit = 0.1_Real; // to do: remove default + Real FrazilIceSalinity = IceRefSal; // Global constant + Real LatFrazil = LatIce; // Global constant + + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb, + const Real H, Real &SumIceThickness, + Real &SumSalt, Real &SumEnergy, Real &HTend, + Real &TTend, Real &STend, + const Real Tfrz) const { + constexpr Real Eps = 1.0e-12_Real; + + if (SumIceThickness <= Eps) { // potential leak if we dont redistribute + SumIceThickness = 0.0_Real; + SumSalt = 0.0_Real; + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + return; + } + + const Real potential = H * Cp0Sw * RhoSw * (CT - Tfrz); + const Real availableEnergy = Kokkos::max(0.0_Real, potential); + + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + + Real meltThickness = + availableEnergy / + (LatFrazil * RhoSw); // mass in pseudo-thickness units + meltThickness = Kokkos::min(meltThickness, SumIceThickness); + meltThickness = Kokkos::min(meltThickness, + H * massLimit); // also 0.1h lim on added mass + const Real meltAverageSalinity = SumSalt / SumIceThickness; + const Real meltingEnergy = + meltThickness * LatFrazil; // (energy needed to melt) + + HTend = +meltThickness * (1 + meltAverageSalinity); // (>0 so HTend>0) + TTend = + -meltingEnergy / Cp0Sw + // (TTend < 0 when melting for phase change) + meltThickness * + Tfrz; // (added the enthalpy of the melt water at freezing temp) + STend = +meltAverageSalinity * meltThickness; // (STend > 0 when melting) + + SumIceThickness -= meltThickness; + SumSalt -= meltAverageSalinity * meltThickness; + SumEnergy -= + meltingEnergy; // necessarily non-conservative by construction + } +}; + +class Frazil { + public: + static void init(); + /// Creates a new frazil object and stores it in the AllFrazil map. + static Frazil *create(const std::string &Name); + + /// Retrieve frazil object by name. + static Frazil *get(const std::string &Name); + + /// Retrieve default frazil object. + static Frazil *getDefault(); + + /// Destructor + ~Frazil(); + + /// Deallocates arrays + static void clear(); + + /// Remove frazil object by name. + static void erase(std::string InName); ///< [in] name to remove + + Array2DReal FrazilTTend; + Array2DReal FrazilSTend; + Array2DReal FrazilHTend; + Array1DReal AccMIce; + Array1DReal AccEIce; + Array1DReal AccMLiq; + Array1DReal AccELiq; + Array1DReal AccMSalt; + + void computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &H); + void computeFrazilBasicImpl(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &LayerH); + + bool conservationCheck = false; + Real depthLimit = -1.0_Real; + + private: + static Frazil *DefaultFrazil; + static std::map> AllFrazil; + + Frazil(const HorzMesh *Mesh, const VertCoord *VCoord); + + // Forbid copy and move construction/assignment. + Frazil(const Frazil &) = delete; + Frazil &operator=(const Frazil &) = delete; + Frazil(Frazil &&) = delete; + Frazil &operator=(Frazil &&) = delete; + + FrazilType frazilChoice; + BasicFrazilFormation computeBasicFrazilFormation; + BasicFrazilMelt computeBasicFrazilMelt; + + I4 NCellsAll; + I4 NChunks; + + const HorzMesh *MeshPtr; + const VertCoord *VCoordPtr; + + void checkColumnConservation() const; +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index ea1e61043f06..4c6421a7bf6a 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -11,6 +11,7 @@ #include "Eos.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -50,6 +51,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time AuxiliaryState::clear(); Forcing::clear(); OceanState::clear(); + Frazil::clear(); VertAdv::clear(); VertCoord::clear(); Dimension::clear(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 7218941ee1c2..2fd751ea692e 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -16,6 +16,7 @@ #include "Error.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -208,6 +209,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { Forcing::init(); AuxiliaryState::init(); Eos::init(); + Frazil::init(); PressureGrad::init(); VertMix::init(); Tendencies::init(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 419c105292a6..86403d407647 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -315,6 +315,10 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); + Err += TendConfig.get("FrazilTendencyEnable", this->FrazilTerm.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: FrazilTendencyEnable not found in TendConfig"); + Err += TendConfig.get("SurfaceTracerRestoringEnable", this->SurfaceTracerRestoring.Enabled); CHECK_ERROR_ABORT( @@ -478,7 +482,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies EqState), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), - CustomThicknessTend(InCustomThicknessTend), + FrazilTerm(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad), VMix(VMix) { @@ -972,6 +976,15 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:sfcTracerForcing", 2); } + if (FrazilTerm.Enabled) { + Pacer::start("Tend:frazil", 2); + const auto &PressureMid = VCoord->PressureMid; + Array2DReal PseudoThickness = State->getPseudoThickness(ThickTimeLevel); + FrazilTerm(PseudoThicknessTend, TracerTend, TracerArray, PressureMid, + PseudoThickness); + Pacer::stop("Tend:frazil", 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 c60a2783ddb3..4c17069184e4 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -80,6 +80,7 @@ class Tendencies { TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; SurfaceTracerRestoringOnCell SurfaceTracerRestoring; + FrazilOnCell FrazilTerm; std::string Name; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 2353f38049bd..09a2d6121740 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -12,6 +12,7 @@ #include "AuxiliaryState.h" #include "DataTypes.h" #include "Eos.h" +#include "Error.h" #include "HorzMesh.h" #include "HorzOperators.h" #include "OceanState.h" @@ -127,6 +128,73 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( const HorzMesh *Mesh) {} +FrazilOnCell::FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) + : NCellsAll(Mesh->NCellsAll), TempTracerIndex(-1), SaltTracerIndex(-1), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) { + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + OMEGA_REQUIRE(TempTracerIndex >= 0, + "FrazilOnCell: Temperature tracer index is undefined"); + OMEGA_REQUIRE(SaltTracerIndex >= 0, + "FrazilOnCell: Salinity tracer index is undefined"); +} + +void FrazilOnCell::operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const { + auto *Frazil = Frazil::getDefault(); + if (!Enabled || !Frazil) { + return; + } + + deepCopy(Frazil->FrazilTTend, 0.0_Real); + deepCopy(Frazil->FrazilSTend, 0.0_Real); + deepCopy(Frazil->FrazilHTend, 0.0_Real); + deepCopy(Frazil->AccMIce, 0.0_Real); + deepCopy(Frazil->AccEIce, 0.0_Real); + deepCopy(Frazil->AccMLiq, 0.0_Real); + deepCopy(Frazil->AccELiq, 0.0_Real); + deepCopy(Frazil->AccMSalt, 0.0_Real); + + const auto ConservTemp = + Kokkos::subview(TracerArray, TempTracerIndex, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, SaltTracerIndex, Kokkos::ALL, Kokkos::ALL); + + Frazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid, + PseudoThickness); + + const auto FrazilHTend = Frazil->FrazilHTend; + const auto FrazilTTend = Frazil->FrazilTTend; + const auto FrazilSTend = Frazil->FrazilSTend; + const I4 TempIndex = TempTracerIndex; + const I4 SaltIndex = SaltTracerIndex; + + OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, TracerTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); + + parallelForOuter( + {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTend(ICell, K) += LocFrazilHTend(ICell, K); + LocTracerTend(TempIndex, ICell, K) += LocFrazilTTend(ICell, K); + LocTracerTend(SaltIndex, ICell, K) += LocFrazilSTend(ICell, K); + }); + }); +} + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const VertCoord *VCoord = this->VerticalCoord; diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index a2a8dc09492c..ff99feb8c92a 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -12,6 +12,7 @@ #include "AuxiliaryState.h" #include "Eos.h" +#include "Frazil.h" #include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" @@ -708,5 +709,26 @@ class SurfaceTracerRestoringOnCell { } }; +/// Frazil tendency hook-up term +class FrazilOnCell { + public: + bool Enabled = false; + + FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + void operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const; + + private: + I4 NCellsAll; + I4 TempTracerIndex; + I4 SaltTracerIndex; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + } // namespace OMEGA #endif diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 5ff4c95e5ea5..aff93f2fa86d 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -512,6 +512,13 @@ add_omega_test( "-n;2" ) +add_omega_test( + FRAZIL_TEST + testFrazil.exe + ocn/FrazilTest.cpp + "-n;2" +) + ################ # VertCoord test ################ diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 8143c9fcbfb3..7c7cfa155ee7 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -117,12 +117,13 @@ void initIOStreamTest(Clock *&ModelClock // Model clock // Initialize VertMix VertMix::init(); + // Initialize Tracers before Tendencies so tracer indices are available + // during FrazilOnCell construction inside Tendencies::init(). + Tracers::init(); + // Intialize Tendencies Tendencies::init(); - // Initialize Tracers - Tracers::init(); - // IOStream::validateAll() depends on Forcing::init() so Forcing fields // are registered before stream validation. Forcing::init(); diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp new file mode 100644 index 000000000000..9c159073b420 --- /dev/null +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -0,0 +1,497 @@ +//===-- Test driver for OMEGA Frazil (basic option) --------------*- C++ +//-*-===// +// +/// \file +/// \brief Minimal test driver for OMEGA frazil formation functor +// +//===-----------------------------------------------------------------------===/ + +#include "Frazil.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Eos.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanTestCommon.h" +#include "OmegaKokkos.h" +#include "Pacer.h" +#include "TimeMgr.h" +#include "VertCoord.h" +#include "mpi.h" + +#include + +using namespace OMEGA; + +constexpr int NVertLayers = 60; + +void initFrazilTest(const std::string &mesh) { + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + LOG_INFO("------ Frazil Unit Tests ------"); + + Config("Omega"); + Config::readAll("omega.yml"); + + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + IO::init(DefComm); + Decomp::init(mesh); + Field::init(ModelClock); + IOStream::init(ModelClock); + Halo::init(); + HorzMesh::init(ModelClock); + VertCoord::init(false); + Eos::init(); + Frazil::init(); + + // Frazil::init only creates the default object when + // Tendencies.FrazilTendencyEnable is true in config. + // Unit tests need a frazil object regardless of that runtime switch. + if (Frazil::getDefault() == nullptr) { + Frazil *CreatedDefault = Frazil::create("Default"); + if (CreatedDefault == nullptr) { + ABORT_ERROR("FrazilTest: failed to create default Frazil object"); + } + } +} + +void finalizeFrazilTest() { + Frazil::clear(); + VertCoord::clear(); + HorzMesh::clear(); + Halo::clear(); + Decomp::clear(); + Field::clear(); + Dimension::clear(); + IOStream::finalize(); + MachEnv::removeAll(); +} + +// this test only excercises the frazil formation functor (no melt) +// in a warm case: the frazil FORMATION terms should all be zero +void testBasicFrazilFormationWarm() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = 10.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + BasicFrazilFormation ComputeFrazilFormation; + ComputeFrazilFormation.massLimit = 0.1_Real; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend, + TTend, STend, CTfrz); + + if (!isApprox(AccMIce, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccMIce, got {}", + AccMIce); + } + + if (!isApprox(AccMSalt, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccMSalt, got {}", + AccMSalt); + } + + if (!isApprox(AccEIce, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccEIce, got {}", + AccEIce); + } + if (!isApprox(HTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero HTend, got {}", + HTend); + } + + if (!isApprox(TTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero TTend, got {}", + TTend); + } + + if (!isApprox(STend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero STend, got {}", + STend); + } + LOG_INFO("FrazilBasicFormationTestWarm: AccMIce = {}, AccMLiq = {}, " + "AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +// this test only excercises the frazil formation functor (no melt) +// in a cold case: the frazil terms should be +// - strictly positive for ice, liquid, and salt mass +//- strictly negative for ice and liquid energy +// - positive for T tendency and negative for S, H tendencies +void testBasicFrazilFormationCold() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = -2.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + BasicFrazilFormation ComputeFrazilFormation; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend, + TTend, STend, CTfrz); + + if (AccMIce <= 0.0_Real) { + ABORT_ERROR("FrazilBasicFormationTestCold: accumulated ice mass is " + "non-positive: {}", + AccMIce); + } + + if (AccMSalt <= 0.0_Real) { + ABORT_ERROR("FrazilBasicFormationTestCold: accumulated salt mass is " + "non-positive: {}", + AccMSalt); + } + + if (AccEIce >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: accumulated ice energy is positive " + "(exp. negative): {}", + AccEIce); + } + if (HTend >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: HTend is positive (exp. negative): {}", + HTend); + } + if (TTend <= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: TTend is negative (exp. positive): {}", + TTend); + } + if (STend >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: STend is positive (exp. negative): {}", + STend); + } + LOG_INFO("FrazilBasicFormationTestCold: AccMIce = {}, AccMLiq = {}, " + "AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +// this test exercises the frazil formation and melt functors +// in a column of water with both cold and warm layers. +// It turns to frazil column conservation check. +// In dev, there is extra verbose logging in the frazil code (TBRemoved). +void testComputeFrazilColumn() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilTestColumn: test frazil object is null"); + } + + const Real RTol = 1e-10_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 4) { + ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - 0.5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; + deepCopy(CT, CTH); + + const bool SavedConservationCheck = TestFrazil->conservationCheck; + TestFrazil->conservationCheck = false; + TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + if (HTendH(ICell, KBottom0) >= 0.0_Real || + TTendH(ICell, KBottom0) <= 0.0_Real || + STendH(ICell, KBottom0) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: bottom cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KBottom1) >= 0.0_Real || + TTendH(ICell, KBottom1) <= 0.0_Real || + STendH(ICell, KBottom1) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: second cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KWarm) < 0.0_Real || TTendH(ICell, KWarm) > 0.0_Real || + STendH(ICell, KWarm) < 0.0_Real) { + ABORT_ERROR("FrazilTestColumn: warm layer sign check failed (HTend>=0, " + "TTend<=0, STend>=0 expected)"); + } + + if (HTendH(ICell, KTopCold) >= 0.0_Real || + TTendH(ICell, KTopCold) <= 0.0_Real || + STendH(ICell, KTopCold) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: top cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + LOG_INFO( + "FrazilTestColumn: XTend branch-switching checks passed for ICell={}", + ICell); +} + +// this test exercises the frazil formation and melt functors +// with a depth limit set. Layers deeper than the depth limit +// should have zero frazil tendencies. +void testComputeFrazilDepthLimit() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilDepthLimitTest: test frazil object is null"); + } + + const Real RTol = 1e-12_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 10) { + ABORT_ERROR("FrazilDepthLimitTest: cell {} has fewer than 10 layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - 0.5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; + deepCopy(CT, CTH); + + const bool SavedConservationCheck = TestFrazil->conservationCheck; + const Real SavedDepthLimit = TestFrazil->depthLimit; + const Real TestDepthLimit = 35.0_Real; // this needs to be positive + // if TestDepthLimit is negative, test will fail: + // - the code assume depthlimit < 0 mean no limit (i.e. full depth frazil) + // - the test below will exclude all layers and fail because Tend !=0. + + // Populate GeomZMid explicitly for the test column. + auto GeomZMidH = createHostMirrorCopy(VCoord->GeomZMid); + for (I4 K = KMin; K <= KMax; ++K) { + GeomZMidH(ICell, K) = -10.0_Real * (K - KMin + 1); + } + deepCopy(VCoord->GeomZMid, GeomZMidH); + + TestFrazil->conservationCheck = false; + TestFrazil->depthLimit = TestDepthLimit; + TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; + TestFrazil->depthLimit = SavedDepthLimit; + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + bool FoundExcludedLayer = false; + for (I4 K = KMin; K <= KMax; ++K) { + const Real Depth = GeomZMidH(ICell, K); + const Real AbsDepth = Depth < 0.0_Real ? -Depth : Depth; + + if (AbsDepth > TestDepthLimit) { + FoundExcludedLayer = true; + if (!isApprox(HTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(TTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(STendH(ICell, K), 0.0_Real, RTol)) { + ABORT_ERROR("FrazilDepthLimitTest: excluded layer K={} has " + "non-zero tendencies (HTend={}, TTend={}, STend={})", + K, HTendH(ICell, K), TTendH(ICell, K), + STendH(ICell, K)); + } + } + } + + if (!FoundExcludedLayer) { + ABORT_ERROR("FrazilDepthLimitTest: no layers excluded for ICell={} with " + "depthLimit={}", + ICell, TestDepthLimit); + } + + LOG_INFO("FrazilDepthLimitTest: depthLimit={} exclusion check passed for " + "ICell={}", + TestDepthLimit, ICell); +} + +void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { + initFrazilTest(MeshFile); + testBasicFrazilFormationCold(); + testBasicFrazilFormationWarm(); + testComputeFrazilColumn(); + testComputeFrazilDepthLimit(); + finalizeFrazilTest(); +} + +int main(int argc, char *argv[]) { + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + frazilTest(); + + LOG_INFO("------ Frazil Unit Tests Successful ------"); + + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + return 0; +} diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 13f8ed2ba2a4..ac11def2793d 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -142,6 +142,14 @@ int initTendenciesTest(const std::string &mesh) { VertCoord::init(); Tracers::init(); + Frazil::init(); + if (Frazil::getDefault() == nullptr) { + Frazil *CreatedDefault = Frazil::create("Default"); + if (CreatedDefault == nullptr) { + Err++; + LOG_ERROR("TendenciesTest: failed to create default Frazil object"); + } + } VertAdv::init(); PressureGrad::init(); Eos::init(); @@ -316,6 +324,142 @@ int testTendencies() { DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + // Verify frazil tendencies are plumbed into pseudo-thickness and tracer + // tendencies by comparing runs with FrazilTerm disabled/enabled. + I4 TempTracerIndex = -1; + I4 SaltTracerIndex = -1; + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + if (TempTracerIndex < 0 || SaltTracerIndex < 0) { + Err++; + LOG_ERROR("TendenciesTest: missing Temperature/Salinity tracer indices"); + } else { + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + auto TracerArrayH = createHostMirrorCopy(TracerArray); + + // Set Salinity in active layers to [31, 33], top-to-bottom ramp. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + for (I4 K = KMin; K <= KMax; ++K) { + Real frac = 0.0_Real; + if (KMax > KMin) { + frac = + static_cast(K - KMin) / static_cast(KMax - KMin); + } + TracerArrayH(SaltTracerIndex, ICell, K) = + 31.0_Real + 2.0_Real * frac; + } + } + + // Keep Temperature warm everywhere first, then add one cold point. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + for (I4 K = KMin; K <= KMax; ++K) { + TracerArrayH(TempTracerIndex, ICell, K) = 8.0_Real; + } + } + const I4 FrazilTestCell = 0; + const I4 FrazilTestLayer = MinLayerCellH(FrazilTestCell); + TracerArrayH(TempTracerIndex, FrazilTestCell, FrazilTestLayer) = + -2.0_Real; + deepCopy(TracerArray, TracerArrayH); + + Array2DReal BaselinePseudoThicknessTend( + "BaselinePseudoThicknessTend", Mesh->NCellsSize, VCoord->NVertLayers); + Array3DReal BaselineTracerTend("BaselineTracerTend", + Tracers::getNumTracers(), Mesh->NCellsSize, + VCoord->NVertLayers); + + const bool OrigFrazilEnabled = DefTendencies->FrazilTerm.Enabled; + + DefTendencies->FrazilTerm.Enabled = false; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + deepCopy(BaselinePseudoThicknessTend, DefTendencies->PseudoThicknessTend); + deepCopy(BaselineTracerTend, DefTendencies->TracerTend); + + deepCopy(DefTendencies->PseudoThicknessTend, 0.0_Real); + deepCopy(DefTendencies->NormalVelocityTend, 0.0_Real); + deepCopy(DefTendencies->TracerTend, 0.0_Real); + + DefTendencies->FrazilTerm.Enabled = true; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + Array2DReal PseudoThicknessTendDiff( + "PseudoThicknessTendDiff", Mesh->NCellsSize, VCoord->NVertLayers); + Array2DReal TempTracerTendDiff("TempTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + Array2DReal SaltTracerTendDiff("SaltTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + + OMEGA_SCOPE(LocPseudoThicknessTendDiff, PseudoThicknessTendDiff); + OMEGA_SCOPE(LocTempTracerTendDiff, TempTracerTendDiff); + OMEGA_SCOPE(LocSaltTracerTendDiff, SaltTracerTendDiff); + OMEGA_SCOPE(LocPseudoThicknessTend, DefTendencies->PseudoThicknessTend); + OMEGA_SCOPE(LocBaselinePseudoThicknessTend, BaselinePseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, DefTendencies->TracerTend); + OMEGA_SCOPE(LocBaselineTracerTend, BaselineTracerTend); + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + + parallelForOuter( + "TendenciesTest:FrazilTendDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTendDiff(ICell, K) = + Kokkos::abs(LocPseudoThicknessTend(ICell, K) - + LocBaselinePseudoThicknessTend(ICell, K)); + LocTempTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(TempTracerIndex, ICell, K) - + LocBaselineTracerTend(TempTracerIndex, ICell, K)); + LocSaltTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(SaltTracerIndex, ICell, K) - + LocBaselineTracerTend(SaltTracerIndex, ICell, K)); + }); + }); + + const Real PseudoDelta = sum(PseudoThicknessTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real TempDelta = sum(TempTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real SaltDelta = sum(SaltTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + + constexpr Real FrazilDeltaATol = 1e-12_Real; + if (!Kokkos::isfinite(PseudoDelta) || + isApprox(PseudoDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change " + "PseudoThicknessTend"); + } + if (!Kokkos::isfinite(TempDelta) || + isApprox(TempDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Temperature tracer " + "tendency"); + } + if (!Kokkos::isfinite(SaltDelta) || + isApprox(SaltDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Salinity tracer " + "tendency"); + } + + DefTendencies->FrazilTerm.Enabled = OrigFrazilEnabled; + } + // Test surface tracer forcing with enthalpy terms (TEOS-10 CtFrz path) const int TracerForcingErr = testSfcTracerForcing(); Err += TracerForcingErr; @@ -336,12 +480,12 @@ int testTendencies() { LOG_ERROR("TendenciesTest: PseudoThickTend FAIL"); } - const Real TraceTendSum = + const Real TracerTendSum = sum(DefTendencies->TracerTend, NTracers, NCellsOwned, VCoord->MinLayerCell, VCoord->MaxLayerCell); - if (!Kokkos::isfinite(TraceTendSum) || TraceTendSum == 0) { + if (!Kokkos::isfinite(TracerTendSum) || TracerTendSum == 0) { Err++; - LOG_ERROR("TendenciesTest: TraceTendSum FAIL"); + LOG_ERROR("TendenciesTest: TracerTendSum FAIL"); } Tendencies::clear(); @@ -728,6 +872,7 @@ int testSfcThicknessForcing() { void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); + Frazil::clear(); PressureGrad::clear(); VertMix::destroyInstance(); Eos::destroyInstance();