From 7bb8d3e4d97f9ad4152e53a36be3647a907ac2ad Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 19 Jun 2026 14:18:42 -0700 Subject: [PATCH 01/11] 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 | 72 +++++++ 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, 597 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 94aa1a3b3daf..75cde0a2d0da 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -55,6 +55,8 @@ Omega: SurfaceTracerRestoringEnable: false BottomDragTendencyEnable: false BottomDragCoeff: 0.0 + SfcThicknessForcingTendencyEnable: false + SfcTracerForcingTendencyEnable: false TracerHorzAdvTendencyEnable: true TracerDiffTendencyEnable: true EddyDiff2: 10.0 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 e6ee85a5045c..d4565674752e 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 @@ -55,6 +57,8 @@ the currently available tendency terms: | SfcStressForcingOnEdge | SfcStressForcingTendencyEnable | enable/disable term | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term | | BottomDragCoeff | bottom drag coefficient +| SfcThicknessForcingOnCell | SfcThicknessForcingTendencyEnable | enable/disable term +| SfcTracerForcingOnCell | SfcTracerForcingTendencyEnable | enable/disable term | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term ## Second Order Horizontal Advection Algorithm @@ -140,5 +144,6 @@ Tracer higer order convergence example of a cosine bell advected on a sphere sho ## See Also -Additional information on forcing (currently wind forcing and surface tracer -restoring) is detailed in [](omega-user-forcing). +Additional information on forcing, including surface stress forcing, +surface 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 a85e24d37dc6..e526adae0046 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -245,6 +245,17 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config Err += TendConfig.get("BottomDragCoeff", this->BottomDrag.Coeff); CHECK_ERROR_ABORT(Err, "Tendencies: BottomDragCoeff not found in TendConfig"); + 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); @@ -418,6 +429,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), BottomDrag(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), @@ -467,6 +480,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); @@ -510,6 +524,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, @@ -729,6 +769,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); @@ -850,6 +891,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 ac826f915a85..3997ea5583f1 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -69,6 +69,8 @@ class Tendencies { VelocityHyperDiffOnEdge VelocityHyperDiff; SfcStressForcingOnEdge SfcStressForcing; BottomDragOnEdge BottomDrag; + 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 389d6bb1aafb..f1552d7b09b0 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -364,6 +364,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 ec458273c351..6d166f5baa95 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->BottomDrag.Enabled = false; DefVAdv->ThickVertAdvEnabled = false; From d5168b7b83106c104fcae41da9f9cf5923406d9d Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 09:18:03 -0700 Subject: [PATCH 02/11] 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 a530549d493d..86f5e31b2b6d 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -756,6 +756,12 @@ class Eos { /// Convert potential temperature to Conservative Temperature Real calcCtFromPt(const Real &Sa, const Real &Pt) const; + /// Calculate freezing Conservative Temperature for TEOS-10. + /// Aborts if EOS is not TEOS-10: CT freezing is not yet implemented + /// for other equation-of-state choices. + Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) const; + /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index e526adae0046..0c5f2106056c 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -430,7 +430,8 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), BottomDrag(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), @@ -908,16 +909,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 f1552d7b09b0..1b03115115b3 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" @@ -404,18 +405,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)) { @@ -423,31 +426,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 8c5db79b1fe4e768408edbb3a05bb812dd4fef1e Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:19:13 -0700 Subject: [PATCH 03/11] 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 2f9696f28c5d1217eb4b38ebbaa5fe29b488305c Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 13:46:07 -0700 Subject: [PATCH 04/11] 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 0c5f2106056c..a664429f3f94 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -915,15 +915,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 1b03115115b3..cd70a456e6b9 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -408,17 +408,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)) { @@ -433,20 +432,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 d00bcd252ad012c57f6f9069ccdffa9fcbcc5686 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:27:46 -0700 Subject: [PATCH 05/11] 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 f9d77f59f1c15d0ca7cbd3c043052ba2d1a2af05 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 26 Jun 2026 14:49:55 -0700 Subject: [PATCH 06/11] 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 d4565674752e..392d537cd1e2 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 b99c94d83d734dde826f72570bd99ddc16e31662 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:07:40 -0700 Subject: [PATCH 07/11] 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 a664429f3f94..0c5f2106056c 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -915,17 +915,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 cd70a456e6b9..1b03115115b3 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -408,16 +408,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)) { @@ -432,29 +433,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 ff48c66ba59008e70971a2737803b311701a43e1 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 10:53:52 -0700 Subject: [PATCH 08/11] 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 96c15a70e907a050ef85773f47e28752fb8619cd Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 6 Jul 2026 11:33:03 -0700 Subject: [PATCH 09/11] 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 1b03115115b3..3948576eaf0d 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -426,12 +426,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 2e1305260e717b34a94034a1536389b918204597 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 13:00:31 -0700 Subject: [PATCH 10/11] 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 392d537cd1e2..6cd59e91167c 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -145,5 +145,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 69ca6c24d99b44b2e369b1b4929a75df32079393 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 13 Jul 2026 14:41:04 -0700 Subject: [PATCH 11/11] 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; }