diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..971dc4f15ad3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -58,6 +58,8 @@ Omega: Mode: Implicit Type: Constant BottomDragCoeff: 1.0e-3 + SfcThicknessForcingTendencyEnable: false + SfcTracerForcingTendencyEnable: false TracerHorzAdvTendencyEnable: true TracerDiffTendencyEnable: true EddyDiff2: 10.0 @@ -69,6 +71,7 @@ Omega: VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true PressureGradTendencyEnable: true + FrazilTendencyEnable: false VelVertMixTendencyEnable: true TracerVertMixTendencyEnable: true ManufacturedSolution: @@ -84,6 +87,11 @@ Omega: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + Frazil: + FrazilType: basic + MassLimit: 0.1 + DepthLimit: -1.0 + ConservationCheck: false VertMix: Background: Diffusivity: 1.0e-5 diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 97d38d8ae42d..7e69dc302bc6 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 thickness and tracer flux forcing (actively coupled or data-forced) +- Surface tracer restoring (soon to be ported as a field originating from the coupler) ## Surface stress forcing design @@ -37,6 +38,70 @@ pathways in Omega, currently this includes: - `Omega.Tendencies.SfcStressForcingTendencyEnable` - gates execution of surface stress forcing tendency kernel +## Surface thickness and tracer flux forcing design + +### Surface thickness and tracer flux forcing data flow + +**Thickness equation pathway:** + +1. External fields provide freshwater and salt flux components: + - `SnowFlux`, `RainFlux`, `EvaporationFlux` + - `SeaIceFreshWaterFlux`, `IceRunoffFlux`, `RiverRunoffFlux` + - `SeaIceSaltFlux` +2. `Forcing` stores the flux fields in `TracerForcingVars` +3. The tendency term `SfcThicknessForcingOnCell` sums the freshwater and salt mass fluxes and applies them to +the surface layer pseudo-thickness. + +**Tracer equation pathway:** + +1. External fields provide heat and salt flux components: + - `LatentHeatFlux`, `SensibleHeatFlux` + - `LongWaveHeatFluxUp`, `LongWaveHeatFluxDown` + - `SeaIceHeatFlux`, `ShortWaveHeatFlux` + - mass fluxes which add energy changes (`SnowFlux`, `RainFlux`, `IceRunoffFlux`, `RiverRunoffFlux`) + - `SeaIceSaltFlux` +2. `Forcing` stores the flux fields in `TracerForcingVars` +3. The tendency term `SfcTracerForcingOnCell` converts the summed external heat fluxes to a conservative-temperature tendency, + and applies the external sea-ice salt flux to the top layer salt content thus impacting salinity. + +### Surface thickness and tracer flux forcing key classes/components + +- `TracerForcingVars` + - Stores 13 coupled flux cell-centered fields: 6 freshwater fluxes, 6 heat + fluxes, and 1 salt flux component + - Fields initialized to zero and registered in `Forcing` field group +- `SfcThicknessForcingOnCell` tendency term + - Computes the layer mass contribution (converted to pseudo-thickness): $\sum (\text{SnowFlux} + \text{RainFlux} + \text{EvaporationFlux} + \text{SeaIceFreshWaterFlux} + \text{IceRunoffFlux} + \text{RiverRunoffFlux} + \text{SeaIceSaltFlux}) / \rho_{sw}$ + - Applied only at surface layer (top active layer) using `MinLayerCell` +- `SfcTracerForcingOnCell` tendency term + - For temperature: adds the direct heat fluxes + $Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$ +, the phase change and enthalpy of added mass $(\text{RainFlux} + \text{RiverRunoffFlux}) c^0_{p,sw} C_T^{\text{top}} + (\text{SnowFlux} + \text{IceRunoffFlux})(c^0_{p,sw} C_T^{\text{frz}} - L_{\text{ice}})$, + (where $C_T^{\text{frz}}$ is from EOS at top-layer salinity and pressure), + and scales by $H_{\text{FluxFac}}$. + - For salinity: applies salt flux with unit conversion: $\text{SeaIceSaltFlux} \times S_{\text{FluxFac}}$ + - Applied only at surface layer using `MinLayerCell` + - Uses tracer index validation to apply to specific tracers only +- `Forcing` + - Manages `TracerForcingVars` instance +- `Tendencies` + - Calls `SfcThicknessForcingOnCell` in `computePseudoThicknessTendenciesOnly` + - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring + +### Surface thickness and tracer flux forcing config coupling + +- `Omega.Tendencies.SfcThicknessForcingTendencyEnable` + - gates execution of coupled flux thickness kernel + - controls freshwater and salt flux forcing on sea surface height +- `Omega.Tendencies.SfcTracerForcingTendencyEnable` + - gates execution of coupled flux tracer kernel + - controls direct heat flux forcing on temperature and salt flux forcing on salinity + +## Notes + +- Currently all forcing is applied to the surface layer only. In the future, vertical spreading of river runoff contributions will be needed. +- `SeaIceFreshWaterFlux` is the pure freshwater mass from sea ice. The full mass flux from sea ice is `SeaIceFreshWaterFlux + SeaIceSaltFlux` + ## Surface tracer restoring design ### Surface tracer restoring data flow diff --git a/components/omega/doc/devGuide/Frazil.md b/components/omega/doc/devGuide/Frazil.md new file mode 100644 index 000000000000..4186324be612 --- /dev/null +++ b/components/omega/doc/devGuide/Frazil.md @@ -0,0 +1,93 @@ +(omega-dev-frazil)= + +# Frazil + +This page describes frazil design and implementation details in Omega, +for the current basic pathway. A TEOS-10 option will be added soon. + +## Purpose and coupling points + +Frazil computes phase-change-related tendencies that modify: + +- pseudo-thickness tendency +- temperature tracer tendency +- salinity tracer tendency + +The tendency hook-up is implemented through `FrazilOnCell` in the tracer +tendency phase, where frazil contributions are added to accumulated +`PseudoThicknessTend` and `TracerTend` arrays. + +## Data flow and call sequence + +1. `Tendencies.computeTracerTendenciesOnly` checks + `Tendencies.FrazilTendencyEnable`. +2. If enabled, `FrazilOnCell.operator()` retrieves the default Frazil object + and zeros frazil tendency and accumulator arrays. +3. `FrazilOnCell` extracts `Temperature` and `Salinity` tracer subviews, then calls + `Frazil.computeFrazil(CT, SA, PressureMid, PseudoThickness)`. +4. `Frazil.computeFrazil` dispatches using FrazilType. The only current supported + option is basic. A TEOS-10 option will be added soon. +5. Returned frazil tendencies are added into `PseudoThicknessTend` and + `TracerTend` for active cell layers. + +## Configuration coupling + +Frazil behavior is configured with: + +- `Omega.Tendencies.FrazilTendencyEnable` + - Global switch for applying frazil tendency terms. +- `Omega.Frazil.FrazilType` + - Implementation choice. Only current supported value: `basic`. +- `Omega.Frazil.MassLimit` + - Per-layer mass and thickness limiter used by basic formation and melt. +- `Omega.Frazil.DepthLimit` + - Optional depth cutoff for frazil activity. Negative means no cutoff. +- `Omega.Frazil.ConservationCheck` + - Optional post-compute column conservation diagnostic logging. + +## Basic pathway summary + +- Freezing-point checks are based on conservative temperature and absolute + salinity with pressure-dependent freezing temperature. +- Vertical accumulation order is bottom-to-top within each active column. +- Frazil tendencies are not time-step scaled inside `Frazil`; they are + accumulated as tendency contributions. +- Surface salt redistribution is applied at the top active layer. +- Column accumulators are converted to coupler units at the end of the + per-column loop. + +### Basic pathway (`FrazilType: basic`) + +- Uses simplified energetics: the energy of the super-cooled water sets the +amount of solid ice formed (used constant latent heat of fusion of fresh ice). +Salt is added based on a constant bulk salinity `IceRefSal` (default) or a +manual toggle (for now) using the local salinity and the frazil porosity. +Melting of existing frazil is set by the amount of pure ice that can be melted +by the warm layer, and the enthalpy of melt water at the local freezing point +is added to the temperature tendency. +- Computes local layer tendencies (`HTend`, `TTend`, `STend`) and updates + accumulated frazil stores. +- Applies surface salt redistribution adjustment and converts accumulators to + coupler units at the end of the column loop. +- Warnings: 1) basic frazil formation/melt is does not conserve energy. + 2) Using porosity to set the salt content includes a redistribution of excess + salt in the surface layer, which can be very significant. + +## Existing ctest coverage + +The existing frazil test driver is in +`components/omega/test/ocn/FrazilTest.cpp` and covers: + +- basic frazil formation in cold and warm single-layer states +- mixed warm/cold column behavior with sign checks for branch switching +- depth-limit behavior ensuring excluded deep layers have zero frazil tendency + +## Extensibility + +The frazil choice dispatch and configuration plumbing are retained so future +frazil options can be added without changing Tendencies call sites. + +## Related pages + +- User-facing options: [User Frazil Guide](../userGuide/Frazil.md) +- Tendency container and hook-up: [Tendencies](Tendencies.md) diff --git a/components/omega/doc/devGuide/Tendencies.md b/components/omega/doc/devGuide/Tendencies.md index 443e8e6ddcd5..32ce44642451 100644 --- a/components/omega/doc/devGuide/Tendencies.md +++ b/components/omega/doc/devGuide/Tendencies.md @@ -64,6 +64,22 @@ Tendencies.computeVelocityTendencies(State, AuxState, ThickTimeLevel, VelTimeLev To call only the tracer tendency terms: Tendencies.computeTracerTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel); +## Frazil tendencies + +Frazil tendencies are calculated by the Frazil class, not the TendencyTerms. +The tendencies are integrated in the tracer tendency compute path and read in from +the frazil computation. When +`Omega.Tendencies.FrazilTendencyEnable` is true, +`computeTracerTendenciesOnly` invokes `FrazilOnCell`, which calls the default +`Frazil` object and adds frazil contributions to: + +- `PseudoThicknessTend` +- `TracerTend` for `Temperature` +- `TracerTend` for `Salinity` + +The frazil implementation choice and parameters are read from the +`Omega.Frazil` block. Detailed algorithm notes are in [Frazil](Frazil.md). + ## Removal of tendencies To erase a specific named tendencies instance use `erase` ```c++ diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 5fa72197132f..2fd762e28b3e 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 mass and tracer flux forcing, and +surface tracer restoring) is detailed in [](omega-dev-forcing). diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index dfc81a284670..30f0a777ba22 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -54,6 +54,7 @@ userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv userGuide/Forcing +userGuide/Frazil userGuide/SfcCoupling userGuide/Analysis ``` @@ -102,6 +103,7 @@ devGuide/Timing devGuide/VerticalMixingCoeff devGuide/VertAdv devGuide/Forcing +devGuide/Frazil devGuide/SfcCoupling devGuide/Analysis ``` diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index d8fac4383730..01eeff4b5e3d 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 tracer flux forcing (mass, energy and salt) - Surface tracer restoring ## Surface stress forcing @@ -39,6 +40,74 @@ Surface stress forcing uses surface stress input fields: These are stored in forcing variables and used to form edge-normal stress (`NormalStressEdge`) that enters momentum tendencies. +## Surface thickness and tracer flux forcing + +Surface thickness and tracer flux forcing applies ocean-atmosphere and ocean-sea ice fluxes from the other model +components (atmosphere, sea ice) to the thickness and tracer equations. This enables +the ocean to respond to heat, freshwater, and salt exchanges at the surface. These fluxes can be from data or (active) coupled components. + +### Surface thickness and tracer flux forcing configuration + +Surface thickness and tracer flux forcing is controlled by two configuration flags: + +```yaml +Omega: + Tendencies: + SfcThicknessForcingTendencyEnable: false + SfcTracerForcingTendencyEnable: false +``` + +- `Tendencies.SfcThicknessForcingTendencyEnable`: enables coupled freshwater and salt flux forcing on thickness +- `Tendencies.SfcTracerForcingTendencyEnable`: enables coupled heat and salt flux forcing on tracers + + +### Required input fields + +Coupled tracer flux forcing uses 13 auxiliary fields organized by type: + +**Freshwater mass fluxes (kg m⁻² s⁻¹):** +- `SnowFlux`: precipitation from snow +- `RainFlux`: precipitation from rain +- `EvaporationFlux`: evaporative water loss +- `SeaIceFreshWaterFlux`: freshwater mass flux from sea-ice melt or formation +- `IceRunoffFlux`: runoff from land ice +- `RiverRunoffFlux`: runoff from rivers + +**Heat fluxes (W m⁻²):** +- `LatentHeatFlux`: latent heat transfer +- `SensibleHeatFlux`: sensible heat transfer +- `LongWaveHeatFluxUp`: upward longwave radiation +- `LongWaveHeatFluxDown`: downward longwave radiation +- `SeaIceHeatFlux`: heat/energy from sea-ice interaction (incl. enthalpy of meltwater) +- `ShortWaveHeatFlux`: shortwave (solar) radiation + +**Salt mass flux (kg m⁻² s⁻¹):** +- `SeaIceSaltFlux`: salt flux from sea-ice formation/melt processes + +These fields are populated by external coupling components (typically atmosphere +and ice models). Omega assumes the incoming values match the documented units. +For now, there are assumed to come from a `forcing.nc` file, but later will be provided +by the equivalent `ocn_comp_mct.F`. + +### Notes + +- Coupled fluxes are applied only at the surface layer (top active layer) for each cell. In the future, vertical spreading of contributions from river runoff will be needed. +- Pseudo-thickness tendency is computed from the (six) freshwater mass fluxes and the salt mass flux + `SeaIceSaltFlux`, converted to a pseudo-thickness change. +- Temperature tendency is computed from direct heat flux plus + mass-flux enthalpy terms, converted to conservative-temperature tendency via + $H_{\text{FluxFac}} = 1.0 / (\rho_{sw} c^0_{p,sw})$ where $c^0_{p,sw}$ is the reference + specific heat of seawater defined by TEOS-10. The enthalpy associated with mass fluxes is currently hard-coded to SST for liquid fluxes and the freezing temperature for solid fluxes (which are melted using a constant latent heat of fusion). Note that the enthalpy of liquid meltwater from sea ice is already included in `SeaIceHeatFlux`. +- Salinity tendency from `SeaIceSaltFlux` is scaled by + $S_{\text{FluxFac}} = 1.0e3 / \rho_{sw}$ to account for unit conversion from + kg/(m²·s) to salinity units (g/kg). +- Fluxes are assumed to be in the documented units (i.e. net mass fluxes); + any unit conversion should be performed by the coupling component before providing flux + values to Omega. +- The reference density used here ($\rho_{sw}$) is not a Boussinesq density, it is the + conversion factor from mass to pseudo-thickness. +- No iceberg fluxes are included for now. + ## Surface tracer restoring Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean diff --git a/components/omega/doc/userGuide/Frazil.md b/components/omega/doc/userGuide/Frazil.md new file mode 100644 index 000000000000..e4078a8ca141 --- /dev/null +++ b/components/omega/doc/userGuide/Frazil.md @@ -0,0 +1,63 @@ +(omega-user-frazil)= + +# Frazil + +This page describes user-facing configuration for frazil tendencies in Omega. +Frazil physics represents the formation and melt of frazil ice within +the ocean water column. It impacts local layer pseudo-thickness (i.e. mass), +temperature, and salinity tendencies. The vertical sum of the frazil energy, +mass of water and mass of salt are passed to the coupler (if coupled) or +discarded (in ocean standalone mode). + +## Configuration overview + +Frazil behavior is controlled by one enable switch in Tendencies and one +Frazil configuration block: + +```yaml +Omega: + Tendencies: + FrazilTendencyEnable: true + + Frazil: + FrazilType: basic + MassLimit: 0.1 + DepthLimit: -1.0 + ConservationCheck: false +``` + +- `Tendencies.FrazilTendencyEnable` + - Enables/disables application of frazil tendency contributions. +- `Frazil.FrazilType` + - Selects frazil option. + - Current supported option: `basic`. +- `Frazil.MassLimit` + - limits per-layer frazil mass/thickness tendency magnitude (applied + to formation and melt) +- `Frazil.DepthLimit` + - Limits depth range where frazil is computed. + - Negative values mean no depth limit. +- `Frazil.ConservationCheck` + - Enables a column-level diagnostic conservation check with logging. + +## Available frazil option + +Omega currently supports one active frazil pathway: + +- `basic` + - freezing is based on the formation of fresh solid ice, to which + salt is added (similar the mpas-ocean implementation). + +This pathway contributes to: + +- pseudo-thickness tendency +- temperature tracer tendency +- salinity tracer tendency + +## Notes + +- Frazil tendencies are applied through the `Tendencies` tracer-step workflow. +- The frazil tendency hook assumes tracer names include `Temperature` and + `Salinity`. +- For implementation and algorithm details, see + [Developer Frazil Guide](../devGuide/Frazil.md). diff --git a/components/omega/doc/userGuide/Tendencies.md b/components/omega/doc/userGuide/Tendencies.md index 6fc1803b176d..00597e75245c 100644 --- a/components/omega/doc/userGuide/Tendencies.md +++ b/components/omega/doc/userGuide/Tendencies.md @@ -6,3 +6,9 @@ The `Tendencies` class provides a container for the [tendency terms](#omega-user Upon creation of an `Tendencies` instance, these functors are initialized and arrays for the accumulated tendencies are allocated. There are no user-configurable options beyond those for the tendency term functors. + +Frazil is one of the configurable tendency contributions. It is controlled by +`Omega.Tendencies.FrazilTendencyEnable` together with the `Omega.Frazil` +configuration block (`FrazilType`, `MassLimit`, `DepthLimit`, and +`ConservationCheck`). For operational guidance and options, see +[Frazil](Frazil.md). diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 3db2b4c10098..c2cd4b3b4327 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, with direct heat always and mass-flux enthalpy terms gated by thickness forcing, defined on cells | SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells Among the internal data stored by each functor is a `bool` which can enable or @@ -57,6 +59,8 @@ the currently available tendency terms: | | BottomDragTendency:Mode | bottom drag mode; `Implicit` or `Explicit` | | BottomDragTendency:Type | bottom drag type; `Constant` | | BottomDragTendency:BottomDragCoeff | bottom drag coefficient +| SfcThicknessForcingOnCell | SfcThicknessForcingTendencyEnable | enable/disable term +| SfcTracerForcingOnCell | SfcTracerForcingTendencyEnable | enable/disable term | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term ## Second Order Horizontal Advection Algorithm @@ -142,5 +146,6 @@ Tracer higer order convergence example of a cosine bell advected on a sphere sho ## See Also -Additional information on forcing (currently wind forcing and surface tracer -restoring) is detailed in [](omega-user-forcing). +Additional information on forcing, including surface stress forcing, +surface thickness and tracer flux forcing, and surface tracer restoring, is detailed in +[](omega-user-forcing). diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 16ea4da9f5a2..ba36a2bb61ac 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,22 +327,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } -Real Eos::calcPtFromCt(const Real &Sa, const Real &Ct) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct); - } - - return Ct; -} - -Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); - } - - return Pt; -} - /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 5e4fd89cae21..2494d6a2d5d8 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -358,8 +358,8 @@ class Teos10Eos { /// (polynomial error in [-5e-4, 6e-4] K, from GSW package). /// P is relative pressure (gauge pressure in Pa, i.e., absolute pressure /// minus the standard atmosphere). - KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { + static KOKKOS_FUNCTION Real calcCtFreezingTeos10( + const Real Sa, const Real P, const Real SaturationFract) { constexpr Real Sso = 35.16504; constexpr Real C0 = 0.017947064327968736; constexpr Real C1 = -6.076099099929818; @@ -756,10 +756,47 @@ class Eos { const Array2DReal &SpecVol); /// Convert Conservative Temperature to potential temperature - Real calcPtFromCt(const Real &Sa, const Real &Ct) const; + /// For TEOS-10, uses the TEOS-10 polynomial + /// For other EOS choices, conservative temperature is equal to potential + /// temperature + KOKKOS_FUNCTION Real calcPtFromCt(const Real &Sa, const Real &Ct) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct); + } + return Ct; + } + + /// Convert potential temperature to Conservative Temperature. + /// For TEOS-10, uses the TEOS-10 polynomial + /// For other EOS choices, potential temperature equals conservative + /// temperature + KOKKOS_FUNCTION Real calcCtFromPt(const Real &Sa, const Real &Pt) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); + } + return Pt; + } - /// Convert potential temperature to Conservative Temperature - Real calcCtFromPt(const Real &Sa, const Real &Pt) const; + /// Calculate freezing temperature of seawater. + /// For TEOS-10, uses the Roquet et al. 75-term polynomial. + /// For LinearEos, uses a simple linear salinity-dependent approximation. + /// For ConstantEos, returns a constant approximate ocean freezing point. + static KOKKOS_FUNCTION Real calcCtFreezing(EosType Choice, const Real Sa, + const Real P, + const Real SaturationFract) { + if (Choice == EosType::Teos10Eos) { + return Teos10Eos::calcCtFreezingTeos10(Sa, P, SaturationFract); + } + if (Choice == EosType::LinearEos) { + // Linear salinity-dependent freezing point; coefficient -0.054 + // degC/PSU with absolute-to-practical salinity conversion (g/kg -> + // PSU). + constexpr Real Coeff = -0.054_Real; + return Coeff * Sa / Psu2Gpkg; + } + // ConstantEos: constant approximate ocean freezing point (degC) + return -1.9_Real; + } /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 42643dba5385..46af26094113 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -29,18 +29,30 @@ static std::string stripDefault(const std::string &Name) { // mesh/halo. Forcing::Forcing(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo) : Name(stripDefault(Name)), SfcStressForcing(stripDefault(Name), Mesh), - Mesh(Mesh), MeshHalo(MeshHalo) {} + TracerForcing(stripDefault(Name), Mesh), Mesh(Mesh), MeshHalo(MeshHalo) {} // Destructor. Unregisters fields from IO streams. Forcing::~Forcing() { unregisterFields(); } // Register surface stress fields with IO streams for a given mesh. void Forcing::registerFields(const std::string &MeshName) const { - SfcStressForcing.registerFields(MeshName); + if (SfcStressFieldsEnabled) { + SfcStressForcing.registerFields(MeshName); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.registerFields(MeshName); + } } // Unregister surface stress fields from IO streams. -void Forcing::unregisterFields() const { SfcStressForcing.unregisterFields(); } +void Forcing::unregisterFields() const { + if (SfcStressFieldsEnabled) { + SfcStressForcing.unregisterFields(); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.unregisterFields(); + } +} // Create and register a non-default forcing instance. Forcing *Forcing::create(const std::string &Name, const HorzMesh *Mesh, @@ -76,10 +88,9 @@ void Forcing::init() { ABORT_ERROR("Forcing: failed to initialize default forcing state"); } - DefaultForcing->registerFields(DefMesh->MeshName); - Config *OmegaConfig = Config::getOmegaConfig(); DefaultForcing->readConfigOptions(OmegaConfig); + DefaultForcing->registerFields(DefMesh->MeshName); // for now, forcing fields are read at start-up only. // to be extended to include switch from standalone to coupled. // to be moved to a Forcing->prepareForStep(SimTime) method later. @@ -136,10 +147,62 @@ void Forcing::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("Forcing: Unknown InterpType requested"); } + + Config TendConfig("Tendencies"); + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, "Forcing: Tendencies group not found in Config"); + + Err += + TendConfig.get("SfcStressForcingTendencyEnable", SfcStressFieldsEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcStressForcingTendencyEnable not found " + "in Tendencies config"); + + bool SfcThicknessForcingEnabled = false; + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + SfcThicknessForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcThicknessForcingTendencyEnable not " + "found in Tendencies config"); + + bool SfcTracerForcingEnabled = false; + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + SfcTracerForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcTracerForcingTendencyEnable not found " + "in Tendencies config"); + + TracerForcingFieldsEnabled = + SfcThicknessForcingEnabled || SfcTracerForcingEnabled; } // Compute all forcing variables (dispatches to specific computations). -void Forcing::computeAll() const { computeSfcStressForcingOnEdge(); } +void Forcing::computeAll() const { + exchangeHalo(); + computeSfcStressForcingOnEdge(); +} + +// Reset forcing arrays so omitted optional fields remain zero after read. +void Forcing::resetArrays() { + if (SfcStressFieldsEnabled) { + deepCopy(SfcStressForcing.NormalStressEdge, 0.0_Real); + deepCopy(SfcStressForcing.ZonalStressCell, 0.0_Real); + deepCopy(SfcStressForcing.MeridStressCell, 0.0_Real); + } + + if (TracerForcingFieldsEnabled) { + deepCopy(TracerForcing.SnowFluxCell, 0.0_Real); + deepCopy(TracerForcing.RainFluxCell, 0.0_Real); + deepCopy(TracerForcing.EvaporationFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); + } +} // Compute edge-normal stress from cell-center zonal and meridional components. void Forcing::computeSfcStressForcingOnEdge() const { @@ -153,14 +216,17 @@ 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; - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, - OnCell); + if (SfcStressFieldsEnabled) { + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, + OnCell); + } return Err; } @@ -172,13 +238,20 @@ void Forcing::readStreamIntoArrays() { std::string StreamName = "Forcing"; + resetArrays(); + + // Nothing to read if neither stress nor tracer forcing tendencies are + // enabled. + if (!SfcStressFieldsEnabled && !TracerForcingFieldsEnabled) { + return; + } + // Attempt to read stream; if unavailable, log and fall back to zero forcing. Err = IOStream::read(StreamName); if (Err.isFail()) { LOG_INFO("Forcing: Error while reading {} stream, using zero forcing", StreamName); - deepCopy(SfcStressForcing.ZonalStressCell, 0._Real); - deepCopy(SfcStressForcing.MeridStressCell, 0._Real); + resetArrays(); } I4 HaloErr = exchangeHalo(); diff --git a/components/omega/src/ocn/Forcing.h b/components/omega/src/ocn/Forcing.h index 5fdae7e550b8..b061bdcf25d7 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(); @@ -69,6 +71,9 @@ class Forcing { /// Read forcing fields from input stream at startup void readStreamIntoArrays(); + /// Reset all forcing arrays to zero before reading optional fields + void resetArrays(); + /// Compute all forcing variables void computeAll() const; @@ -86,6 +91,8 @@ class Forcing { const HorzMesh *Mesh; Halo *MeshHalo; + bool SfcStressFieldsEnabled = false; + bool TracerForcingFieldsEnabled = false; static Forcing *DefaultForcing; static std::map> AllForcing; diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp new file mode 100644 index 000000000000..c0e968adb0c9 --- /dev/null +++ b/components/omega/src/ocn/Frazil.cpp @@ -0,0 +1,369 @@ +//===-- ocn/Frazil.cpp - Frazil Ice Formation -----------------*- C++ -*-===// +// +// The Frazil class manages frazil tendencies and accumulators. +// This initial implementation only has a basic configuration. +// +//===----------------------------------------------------------------------===// + +#include "Frazil.h" +#include "Config.h" +#include "Eos.h" +#include "Error.h" +#include "Logging.h" + +namespace OMEGA { + +namespace { + +KOKKOS_INLINE_FUNCTION bool isApprox(const Real X, const Real Y, + const Real RTol, const Real ATol = 0) { + if (Kokkos::isnan(X) || Kokkos::isnan(Y) || Kokkos::isinf(X) || + Kokkos::isinf(Y)) { + return false; + } + + return Kokkos::abs(X - Y) <= + Kokkos::max(ATol, RTol * Kokkos::max(Kokkos::abs(X), Kokkos::abs(Y))); +} + +} // namespace + +Frazil *Frazil::DefaultFrazil = nullptr; +std::map> Frazil::AllFrazil; + +/// Constructor for BasicFrazilFormation +BasicFrazilFormation::BasicFrazilFormation() {} + +/// Constructor for BasicFrazilMelt +BasicFrazilMelt::BasicFrazilMelt() {} + +void Frazil::init() { + + if (!HorzMesh::getDefault() || !VertCoord::getDefault()) { + ABORT_ERROR("Frazil::init: HorzMesh and VertCoord must be initialized"); + } + + // Frazil freezing-temperature calculations depend on EOS configuration. + Eos::init(); + + if (!DefaultFrazil) { + Error Err; + bool FrazilTendencyEnable = false; + Config *OmegaConfig = Config::getOmegaConfig(); + Config TendConfig("Tendencies"); + + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, + "Frazil::init: Tendencies group not found in Config"); + + Err += TendConfig.get("FrazilTendencyEnable", FrazilTendencyEnable); + CHECK_ERROR_ABORT( + Err, "Frazil::init: FrazilTendencyEnable not found in Tendencies"); + + if (!FrazilTendencyEnable) { + LOG_INFO("Frazil::init: Frazil tendency disabled; skipping default " + "frazil object creation"); + LOG_INFO("All frazil is off - frazil parameters will be ignored"); + return; + } + + DefaultFrazil = create("Default"); + } +} + +Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) + : frazilChoice(FrazilType::BasicFrazil), computeBasicFrazilFormation(), + computeBasicFrazilMelt(), NCellsAll(Mesh->NCellsAll), + NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh), + VCoordPtr(VCoord) { + + FrazilTTend = + Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilSTend = + Array2DReal("FrazilSTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilHTend = + Array2DReal("FrazilHTend", Mesh->NCellsSize, VCoord->NVertLayers); + + AccMIce = Array1DReal("AccMIce", Mesh->NCellsSize); + AccEIce = Array1DReal("AccEIce", Mesh->NCellsSize); + AccMLiq = Array1DReal("AccMLiq", Mesh->NCellsSize); + AccELiq = Array1DReal("AccELiq", Mesh->NCellsSize); + AccMSalt = Array1DReal("AccMSalt", Mesh->NCellsSize); + + deepCopy(FrazilTTend, 0.0_Real); + deepCopy(FrazilSTend, 0.0_Real); + deepCopy(FrazilHTend, 0.0_Real); + deepCopy(AccMIce, 0.0_Real); + deepCopy(AccEIce, 0.0_Real); + deepCopy(AccMLiq, 0.0_Real); + deepCopy(AccELiq, 0.0_Real); + deepCopy(AccMSalt, 0.0_Real); +} + +Frazil::~Frazil() {} + +Frazil *Frazil::create(const std::string &Name) { + if (AllFrazil.find(Name) != AllFrazil.end()) { + LOG_ERROR("Attempted to create Frazil {} but it already exists", Name); + return nullptr; + } + + auto *NewFrazil = + new Frazil(HorzMesh::getDefault(), VertCoord::getDefault()); + AllFrazil.emplace(Name, NewFrazil); + + Error Err; + Config *OmegaConfig = Config::getOmegaConfig(); + Config FrazilConfig("Frazil"); + Err += OmegaConfig->get(FrazilConfig); + CHECK_ERROR_ABORT(Err, "Frazil::create: Frazil group not found in Config"); + + std::string FrazilTypeStr; + Err += FrazilConfig.get("FrazilType", FrazilTypeStr); + CHECK_ERROR_ABORT(Err, + "Frazil::create: FrazilType not found in Frazil config"); + + if (((FrazilTypeStr == "Basic") || (FrazilTypeStr == "basic") || + (FrazilTypeStr == "BasicFrazil"))) { + NewFrazil->frazilChoice = FrazilType::BasicFrazil; + } else { + ABORT_ERROR( + "Frazil::create: Only basic frazil is supported at the moment"); + } + + Err += FrazilConfig.get("MassLimit", + NewFrazil->computeBasicFrazilFormation.massLimit); + Err += FrazilConfig.get("MassLimit", + NewFrazil->computeBasicFrazilMelt.massLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: MassLimit not found in Frazil config"); + + NewFrazil->computeBasicFrazilMelt.massLimit = + NewFrazil->computeBasicFrazilFormation.massLimit; + + Err += FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); + CHECK_ERROR_ABORT( + Err, "Frazil::create: ConservationCheck not found in Frazil config"); + + Err += FrazilConfig.get("DepthLimit", NewFrazil->depthLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: DepthLimit not found in Frazil config"); + + if (Name == "Default") { + DefaultFrazil = NewFrazil; + } + + return NewFrazil; +} + +Frazil *Frazil::getDefault() { return DefaultFrazil; } + +Frazil *Frazil::get(const std::string &Name) { + auto it = AllFrazil.find(Name); + if (it != AllFrazil.end()) { + return it->second.get(); + } + + LOG_ERROR("Frazil::get: Attempted to retrieve non-existent Frazil {}", Name); + return nullptr; +} + +void Frazil::erase(std::string InName) { + auto *ToErase = get(InName); + AllFrazil.erase(InName); + if (ToErase == DefaultFrazil) { + DefaultFrazil = nullptr; + } +} + +void Frazil::clear() { + AllFrazil.clear(); + DefaultFrazil = nullptr; +} + +void Frazil::checkColumnConservation() const { + auto MinLayerCellH = createHostMirrorCopy(VCoordPtr->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoordPtr->MaxLayerCell); + auto FrazilHTendH = createHostMirrorCopy(FrazilHTend); + auto FrazilTTendH = createHostMirrorCopy(FrazilTTend); + auto FrazilSTendH = createHostMirrorCopy(FrazilSTend); + auto AccMIceH = createHostMirrorCopy(AccMIce); + auto AccMLiqH = createHostMirrorCopy(AccMLiq); + auto AccMSaltH = createHostMirrorCopy(AccMSalt); + auto AccELiqH = createHostMirrorCopy(AccELiq); + auto AccEIceH = createHostMirrorCopy(AccEIce); + + constexpr Real RTol = 1.0e-10_Real; + + for (I4 ICell = 0; ICell < NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + Real MassTend = 0.0_Real; + Real EnergyTend = 0.0_Real; + Real SaltTend = 0.0_Real; + + for (I4 K = KMin; K <= KMax; ++K) { + MassTend += FrazilHTendH(ICell, K); + EnergyTend += FrazilTTendH(ICell, K); + SaltTend += FrazilSTendH(ICell, K); + } + + const Real MassTotal = AccMIceH(ICell) + AccMLiqH(ICell); + const Real EnergyTotal = AccELiqH(ICell) + AccEIceH(ICell); + const Real SaltTotal = AccMSaltH(ICell); + + if (ICell == 0) { + LOG_INFO("Frazil column conservation check: cell {} MassTend={} " + "MassTotal={} " + "EnergyTend={} EnergyTotal={} SaltTend={} SaltTotal={}", + ICell, MassTend * RhoSw, MassTotal, + EnergyTend * Cp0Sw * RhoSw, EnergyTotal, + SaltTend * RhoSw * PPt2Salt, SaltTotal); + LOG_INFO("Frazil column conservation check: cell {} EpsMass={} " + "EpsE={} EpsS={} ", + ICell, MassTend * RhoSw + MassTotal, + EnergyTend * Cp0Sw * RhoSw + EnergyTotal, + SaltTend * RhoSw * PPt2Salt + SaltTotal); + } + + if (!isApprox(-MassTend * RhoSw, MassTotal, RTol)) { + LOG_INFO( + "Frazil column mass check failed: cell {} tendency={} total={}", + ICell, -MassTend * RhoSw, MassTotal); + } + if (!isApprox(-EnergyTend * Cp0Sw * RhoSw, EnergyTotal, RTol)) { + LOG_INFO( + "Frazil column energy check failed: cell {} tendency={} total={}", + ICell, -EnergyTend * Cp0Sw * RhoSw, EnergyTotal); + } + if (!isApprox(-SaltTend * RhoSw * PPt2Salt, SaltTotal, RTol)) { + LOG_INFO( + "Frazil column salt check failed: cell {} tendency={} total={}", + ICell, -SaltTend * RhoSw * PPt2Salt, SaltTotal); + } + } +} + +void Frazil::computeFrazilBasicImpl(const Array2DReal &CT, + const Array2DReal &SA, const Array2DReal &P, + const Array2DReal &LayerH) { + const EosType LocEosChoice = Eos::getInstance()->EosChoice; + const Real LocDepthLimit = depthLimit; + + OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell); + OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid); + + OMEGA_SCOPE(LocComputeBasicFrazilFormation, computeBasicFrazilFormation); + OMEGA_SCOPE(LocComputeBasicFrazilMelt, computeBasicFrazilMelt); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocAccMIce, AccMIce); + OMEGA_SCOPE(LocAccEIce, AccEIce); + OMEGA_SCOPE(LocAccMLiq, AccMLiq); + OMEGA_SCOPE(LocAccELiq, AccELiq); + OMEGA_SCOPE(LocAccMSalt, AccMSalt); + OMEGA_SCOPE(LocIceRefSal, IceRefSal); + OMEGA_SCOPE(LocLatIce, LatIce); + + parallelFor( + {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + I4 Klim = KMax; + bool HasKlim = true; + const bool Limit = (LocDepthLimit >= 0.0_Real); + + if (Limit) { + HasKlim = false; + for (I4 K = KMax; K >= KMin; --K) { + if (Kokkos::abs(LocGeomZMid(ICell, K)) <= LocDepthLimit) { + Klim = K; + HasKlim = true; + break; + } + } + } + + // Explicit accumulation order: bottom layer to top layer. + for (I4 K = KMax; K >= KMin; --K) { + if (!HasKlim || K > Klim) { + LocFrazilHTend(ICell, K) = 0.0_Real; + LocFrazilTTend(ICell, K) = 0.0_Real; + LocFrazilSTend(ICell, K) = 0.0_Real; + continue; + } + + const Real SAIn = SA(ICell, K); + const Real CTIn = CT(ICell, K); + const Real PIn = P(ICell, K); + const Real PDb = PIn * Pa2Db; + const Real H = LayerH(ICell, K); + + const Real Tfrz = + Eos::calcCtFreezing(LocEosChoice, SAIn, PDb, 0.0_Real); + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + if (CTIn < Tfrz) { + LocComputeBasicFrazilFormation( + SAIn, CTIn, PDb, H, LocAccMIce(ICell), LocAccMSalt(ICell), + LocAccEIce(ICell), HTend, TTend, STend, Tfrz); + } else if (LocAccMIce(ICell) > 0.0_Real) { + LocComputeBasicFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMSalt(ICell), LocAccEIce(ICell), + HTend, TTend, STend, Tfrz); + } + + LocFrazilHTend(ICell, K) = HTend; + LocFrazilTTend(ICell, K) = TTend; + LocFrazilSTend(ICell, K) = STend; + } + + // Redistribute excess salt at the surface. No treatment of low + // salinity frazil for now. + LocFrazilSTend(ICell, KMin) += Kokkos::max( + 0.0_Real, LocAccMIce(ICell) * LocIceRefSal - LocAccMSalt(ICell)); + // hijack total terms before the coupling + LocAccMSalt(ICell) = LocAccMIce(ICell) * LocIceRefSal; + LocAccEIce(ICell) = -LocAccMIce(ICell) * LocLatIce; + + // Convert to coupler units + LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; + LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; + LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; + LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; + LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; + }); +} + +void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &LayerH) { + Eos *DefEos = Eos::getInstance(); + if (!DefEos) { + ABORT_ERROR("Frazil::computeFrazil: Eos must be initialized before " + "computeFrazil"); + } + + switch (frazilChoice) { + case FrazilType::BasicFrazil: + computeFrazilBasicImpl(CT, SA, P, LayerH); + break; + // TODO: add additional frazil choices here as new implementations are + // introduced + default: + ABORT_ERROR("Frazil::computeFrazil: Unsupported frazilChoice"); + break; + } + + if (conservationCheck) { + checkColumnConservation(); + } +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h new file mode 100644 index 000000000000..c5e7f11e1b40 --- /dev/null +++ b/components/omega/src/ocn/Frazil.h @@ -0,0 +1,207 @@ +#ifndef OMEGA_FRAZIL_H +#define OMEGA_FRAZIL_H +//===-- ocn/Frazil.h - Frazil Ice Formation -------------------*- C++ -*-===// +// +// The Frazil class manages frazil tendencies and accumulators. +// This initial implementation only has a basic configuration. +// but carries scaffolding for other implementations. +// +//===----------------------------------------------------------------------===// + +#include "GlobalConstants.h" +#include "HorzMesh.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include +#include +#include + +namespace OMEGA { + +enum class FrazilType { + BasicFrazil ///< MPAS-O style basic frazil option + // Additional options can be added here in the future. +}; + +class BasicFrazilFormation { + public: + BasicFrazilFormation(); + + Real massLimit = 0.1_Real; // to do: remove default + Real FrazilIceSalinity = IceRefSal; // Global constant + Real LatFrazil = LatIce; // Global constant + + Real FrazilPorosity = + 1.0_Real; // Internal for now; can move to config later. + + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb, + const Real H, Real &SumIceThickness, + Real &SumSalt, Real &SumEnergy, Real &HTend, + Real &TTend, Real &STend, + const Real Tfrz) const { + + const Real potential = H * Cp0Sw * RhoSw * (CT - Tfrz); + const Real freezingEnergy = Kokkos::max(0.0_Real, -potential); + + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + + Real newFrzThickness = + freezingEnergy / + (LatFrazil * RhoSw); // frazil (ice) mass in pseudo-thickness terms + newFrzThickness = Kokkos::min(newFrzThickness, H * massLimit); + const Real newFrzEnergy = + -newFrzThickness * LatFrazil; // (<0; enthalpy of ice) + + // MANUAL TOGGLE: uncomment line below to use porosity + // const Real FrazilIceSalinity = FrazilPorosity * SA; + + const Real frazilSalinity = Kokkos::min(FrazilIceSalinity, SA); + const Real newSaltContent = + newFrzThickness * frazilSalinity; // in m.(g/kg) + + HTend = -newFrzThickness - newSaltContent; + // TTend below should include the enthalpy associated with the mass flux + // this is intentionally not added here to match the mpas-o implementation + TTend = + -(newFrzEnergy) / Cp0Sw; // (E< 0 so TTend>0) // scaled to h.CT tend + STend = -newSaltContent; + + SumIceThickness += newFrzThickness; + // non-conservation between Sum and HTend by construction in the original + SumSalt += newSaltContent; + SumEnergy += newFrzEnergy; + } +}; + +// TO-DO: is there a new for salt reconciliation at the surface? +// frazilSalinityTendency(minLevelCell(iCell),iCell) = +// frazilSalinityTendency(minLevelCell(iCell),iCell) + & +// max(0.0_RKIND,(sumNewThicknessWeightedSaltContent - +// newThicknessWeightedSaltContent) ) / dt +// accumulatedFrazilIceSalinityNew(iCell) = +// accumulatedFrazilIceSalinityOld(iCell) + newThicknessWeightedSaltContent + +class BasicFrazilMelt { + public: + BasicFrazilMelt(); + + Real massLimit = 0.1_Real; // to do: remove default + Real FrazilIceSalinity = IceRefSal; // Global constant + Real LatFrazil = LatIce; // Global constant + + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb, + const Real H, Real &SumIceThickness, + Real &SumSalt, Real &SumEnergy, Real &HTend, + Real &TTend, Real &STend, + const Real Tfrz) const { + constexpr Real Eps = 1.0e-12_Real; + + if (SumIceThickness <= Eps) { // potential leak if we dont redistribute + SumIceThickness = 0.0_Real; + SumSalt = 0.0_Real; + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + return; + } + + const Real potential = H * Cp0Sw * RhoSw * (CT - Tfrz); + const Real availableEnergy = Kokkos::max(0.0_Real, potential); + + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + + Real meltThickness = + availableEnergy / + (LatFrazil * RhoSw); // mass in pseudo-thickness units + meltThickness = Kokkos::min(meltThickness, SumIceThickness); + meltThickness = Kokkos::min(meltThickness, + H * massLimit); // also 0.1h lim on added mass + const Real meltAverageSalinity = SumSalt / SumIceThickness; + const Real meltingEnergy = + meltThickness * LatFrazil; // (energy needed to melt) + + HTend = +meltThickness * (1 + meltAverageSalinity); // (>0 so HTend>0) + TTend = + -meltingEnergy / Cp0Sw + // (TTend < 0 when melting for phase change) + meltThickness * + Tfrz; // (added the enthalpy of the melt water at freezing temp) + STend = +meltAverageSalinity * meltThickness; // (STend > 0 when melting) + + SumIceThickness -= meltThickness; + SumSalt -= meltAverageSalinity * meltThickness; + SumEnergy -= + meltingEnergy; // necessarily non-conservative by construction + } +}; + +class Frazil { + public: + static void init(); + /// Creates a new frazil object and stores it in the AllFrazil map. + static Frazil *create(const std::string &Name); + + /// Retrieve frazil object by name. + static Frazil *get(const std::string &Name); + + /// Retrieve default frazil object. + static Frazil *getDefault(); + + /// Destructor + ~Frazil(); + + /// Deallocates arrays + static void clear(); + + /// Remove frazil object by name. + static void erase(std::string InName); ///< [in] name to remove + + Array2DReal FrazilTTend; + Array2DReal FrazilSTend; + Array2DReal FrazilHTend; + Array1DReal AccMIce; + Array1DReal AccEIce; + Array1DReal AccMLiq; + Array1DReal AccELiq; + Array1DReal AccMSalt; + + void computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &H); + void computeFrazilBasicImpl(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &LayerH); + + bool conservationCheck = false; + Real depthLimit = -1.0_Real; + + private: + static Frazil *DefaultFrazil; + static std::map> AllFrazil; + + Frazil(const HorzMesh *Mesh, const VertCoord *VCoord); + + // Forbid copy and move construction/assignment. + Frazil(const Frazil &) = delete; + Frazil &operator=(const Frazil &) = delete; + Frazil(Frazil &&) = delete; + Frazil &operator=(Frazil &&) = delete; + + FrazilType frazilChoice; + BasicFrazilFormation computeBasicFrazilFormation; + BasicFrazilMelt computeBasicFrazilMelt; + + I4 NCellsAll; + I4 NChunks; + + const HorzMesh *MeshPtr; + const VertCoord *VCoordPtr; + + void checkColumnConservation() const; +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/src/ocn/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/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index ea1e61043f06..4c6421a7bf6a 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -11,6 +11,7 @@ #include "Eos.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -50,6 +51,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time AuxiliaryState::clear(); Forcing::clear(); OceanState::clear(); + Frazil::clear(); VertAdv::clear(); VertCoord::clear(); Dimension::clear(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 7218941ee1c2..2fd751ea692e 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -16,6 +16,7 @@ #include "Error.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -208,6 +209,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { Forcing::init(); AuxiliaryState::init(); Eos::init(); + Frazil::init(); PressureGrad::init(); VertMix::init(); Tendencies::init(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index d8273a844de0..86403d407647 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -283,6 +283,18 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config } } + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + this->SfcThicknessForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcThicknessForcingTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + this->SfcTracerForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcTracerForcingTendencyEnable not found in TendConfig"); + if (this->TracerDiffusion.Enabled) { Err += TendConfig.get("EddyDiff2", this->TracerDiffusion.EddyDiff2); CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff2 not found in TendConfig"); @@ -303,6 +315,10 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); + Err += TendConfig.get("FrazilTendencyEnable", this->FrazilTerm.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: FrazilTendencyEnable not found in TendConfig"); + Err += TendConfig.get("SurfaceTracerRestoringEnable", this->SurfaceTracerRestoring.Enabled); CHECK_ERROR_ABORT( @@ -461,9 +477,12 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord), VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), ExplicitBottomDrag(Mesh, VCoord), + SfcThicknessForcing(Mesh, VCoord), + SfcTracerForcing(Mesh, VCoord, Tracers::IndxTemp, Tracers::IndxSalt, + EqState), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), - CustomThicknessTend(InCustomThicknessTend), + FrazilTerm(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad), VMix(VMix) { @@ -510,6 +529,7 @@ void Tendencies::computePseudoThicknessTendenciesOnly( OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); OMEGA_SCOPE(LocThicknessFluxDiv, PseudoThicknessFluxDiv); + OMEGA_SCOPE(LocSfcThicknessForcing, SfcThicknessForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -553,6 +573,32 @@ void Tendencies::computePseudoThicknessTendenciesOnly( VAdv->computePseudoThicknessVAdvTend(PseudoThicknessTend); Pacer::stop("Tend:computePseudoThicknessVAdvTend", 2); + if (LocSfcThicknessForcing.Enabled) { + Pacer::start("Tend:sfcThicknessForcing", 2); + const auto *ForcingState = Forcing::getDefault(); + + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &EvaporationFlux = + ForcingState->TracerForcing.EvaporationFluxCell; + const auto &SeaIceFreshWaterFlux = + ForcingState->TracerForcing.SeaIceFreshWaterFluxCell; + const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &RiverRunoffFlux = + ForcingState->TracerForcing.RiverRunoffFluxCell; + const auto &SeaIceSaltFlux = + ForcingState->TracerForcing.SeaIceSaltFluxCell; + + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + LocSfcThicknessForcing(LocPseudoThicknessTend, ICell, SnowFlux, + RainFlux, EvaporationFlux, + SeaIceFreshWaterFlux, IceRunoffFlux, + RiverRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcThicknessForcing", 2); + } + if (CustomThicknessTend) { Pacer::start("Tend:customThicknessTend", 2); CustomThicknessTend(LocPseudoThicknessTend, State, AuxState, @@ -772,6 +818,7 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); + OMEGA_SCOPE(LocSfcTracerForcing, SfcTracerForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -893,6 +940,51 @@ 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 &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, TracerArray, PressureMid, LatentHeatFlux, + SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, + IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcTracerForcing", 2); + } + + if (FrazilTerm.Enabled) { + Pacer::start("Tend:frazil", 2); + const auto &PressureMid = VCoord->PressureMid; + Array2DReal PseudoThickness = State->getPseudoThickness(ThickTimeLevel); + FrazilTerm(PseudoThicknessTend, TracerTend, TracerArray, PressureMid, + PseudoThickness); + Pacer::stop("Tend:frazil", 2); + } + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 877ba22faf2b..4c17069184e4 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -74,10 +74,13 @@ class Tendencies { VelocityHyperDiffOnEdge VelocityHyperDiff; SfcStressForcingOnEdge SfcStressForcing; BottomDragOnEdge ExplicitBottomDrag; + SfcThicknessForcingOnCell SfcThicknessForcing; + SfcTracerForcingOnCell SfcTracerForcing; TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; SurfaceTracerRestoringOnCell SurfaceTracerRestoring; + FrazilOnCell FrazilTerm; std::string Name; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 5b142888142c..09a2d6121740 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -11,6 +11,8 @@ #include "TendencyTerms.h" #include "AuxiliaryState.h" #include "DataTypes.h" +#include "Eos.h" +#include "Error.h" #include "HorzMesh.h" #include "HorzOperators.h" #include "OceanState.h" @@ -71,6 +73,19 @@ 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, + const Eos *EosInst) + : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + EosChoice(EosInst->EosChoice) {} + TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) : HorzontalMesh(Mesh), VerticalCoord(VCoord), @@ -113,6 +128,73 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( const HorzMesh *Mesh) {} +FrazilOnCell::FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) + : NCellsAll(Mesh->NCellsAll), TempTracerIndex(-1), SaltTracerIndex(-1), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) { + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + OMEGA_REQUIRE(TempTracerIndex >= 0, + "FrazilOnCell: Temperature tracer index is undefined"); + OMEGA_REQUIRE(SaltTracerIndex >= 0, + "FrazilOnCell: Salinity tracer index is undefined"); +} + +void FrazilOnCell::operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const { + auto *Frazil = Frazil::getDefault(); + if (!Enabled || !Frazil) { + return; + } + + deepCopy(Frazil->FrazilTTend, 0.0_Real); + deepCopy(Frazil->FrazilSTend, 0.0_Real); + deepCopy(Frazil->FrazilHTend, 0.0_Real); + deepCopy(Frazil->AccMIce, 0.0_Real); + deepCopy(Frazil->AccEIce, 0.0_Real); + deepCopy(Frazil->AccMLiq, 0.0_Real); + deepCopy(Frazil->AccELiq, 0.0_Real); + deepCopy(Frazil->AccMSalt, 0.0_Real); + + const auto ConservTemp = + Kokkos::subview(TracerArray, TempTracerIndex, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, SaltTracerIndex, Kokkos::ALL, Kokkos::ALL); + + Frazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid, + PseudoThickness); + + const auto FrazilHTend = Frazil->FrazilHTend; + const auto FrazilTTend = Frazil->FrazilTTend; + const auto FrazilSTend = Frazil->FrazilSTend; + const I4 TempIndex = TempTracerIndex; + const I4 SaltIndex = SaltTracerIndex; + + OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, TracerTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); + + parallelForOuter( + {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTend(ICell, K) += LocFrazilHTend(ICell, K); + LocTracerTend(TempIndex, ICell, K) += LocFrazilTTend(ICell, K); + LocTracerTend(SaltIndex, ICell, K) += LocFrazilSTend(ICell, K); + }); + }); +} + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const VertCoord *VCoord = this->VerticalCoord; diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 548832b290bf..ff99feb8c92a 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -11,6 +11,8 @@ //===----------------------------------------------------------------------===// #include "AuxiliaryState.h" +#include "Eos.h" +#include "Frazil.h" #include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" @@ -372,6 +374,106 @@ 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, + 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)) { + return; + } + + if (TempIndex >= 0) { + const Real PTopDb = PressureMid(ICell, KTop) * Pa2Db; + const Real SaTop = SaltIndex >= 0 + ? TracerCell(SaltIndex, ICell, KTop) + : 0.0_Real; // not sure we want zero here? + const Real CtFrz = + Eos::calcCtFreezing(EosChoice, SaTop, PTopDb, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); + + // Heat tendencies are due to direct heat fluxes + enthalpy fluxes + // The enthalpy of liquid water is assumed to be: + // - local SST for liquid mass fluxes (rain, rivers) + // - local freezing point for solid --> liq mass fluxes (snow, frozen + // 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) * SFluxFac; + } + } + + private: + I4 TempIndex; + I4 SaltIndex; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + EosType EosChoice; +}; + // Tracer horizontal advection term class TracerHorzAdvOnCell { public: @@ -607,5 +709,26 @@ class SurfaceTracerRestoringOnCell { } }; +/// Frazil tendency hook-up term +class FrazilOnCell { + public: + bool Enabled = false; + + FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + void operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const; + + private: + I4 NCellsAll; + I4 TempTracerIndex; + I4 SaltTracerIndex; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp new file mode 100644 index 000000000000..a6478bb612c4 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -0,0 +1,155 @@ +#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("SeaIceSaltFlux" + 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); +} + +void TracerForcingVars::registerFields(const std::string &MeshName) const { + 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(), 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(), 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(), 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(), NDims, DimNames); + auto LongWaveHeatFluxDownField = Field::create( + LongWaveHeatFluxDownCell.label(), "downward longwave heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + 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); + + 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); + 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()); +} +} // 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..e38d9948f672 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.h @@ -0,0 +1,42 @@ +#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; + + TracerForcingVars(const std::string &Suffix, const HorzMesh *Mesh); + + void registerFields(const std::string &MeshName) const; + void unregisterFields() const; +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 5ff4c95e5ea5..aff93f2fa86d 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -512,6 +512,13 @@ add_omega_test( "-n;2" ) +add_omega_test( + FRAZIL_TEST + testFrazil.exe + ocn/FrazilTest.cpp + "-n;2" +) + ################ # VertCoord test ################ diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 8143c9fcbfb3..7c7cfa155ee7 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -117,12 +117,13 @@ void initIOStreamTest(Clock *&ModelClock // Model clock // Initialize VertMix VertMix::init(); + // Initialize Tracers before Tendencies so tracer indices are available + // during FrazilOnCell construction inside Tendencies::init(). + Tracers::init(); + // Intialize Tendencies Tendencies::init(); - // Initialize Tracers - Tracers::init(); - // IOStream::validateAll() depends on Forcing::init() so Forcing fields // are registered before stream validation. Forcing::init(); diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 0e7eb2c5be76..f8ca084d4559 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -706,6 +706,42 @@ void testBruntVaisalaFreqSqTeos10() { return; } +/// Test all Eos::calcCtFreezing pathways (Teos10, Linear, Constant) +void testCalcCtFreezing() { + const Real RTol = 1e-10; + + constexpr Real SaturationFrac = 0.0; + constexpr Real PDb = 500.0; // pressure in dbar for GSW pathway + constexpr Real SaLocal = 32.0; + + const Real CtTeosExpected = + gsw_ct_freezing_poly(SaLocal, PDb, SaturationFrac); + const Real CtTeos = + Eos::calcCtFreezing(EosType::Teos10Eos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtTeos, CtTeosExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Teos10 FAIL, expected {}, got {}", + CtTeosExpected, CtTeos); + } + + const Real CtLinearExpected = -0.054_Real * SaLocal / Psu2Gpkg; + const Real CtLinear = + Eos::calcCtFreezing(EosType::LinearEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtLinear, CtLinearExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Linear FAIL, expected {}, got {}", + CtLinearExpected, CtLinear); + } + + const Real CtConstExpected = -1.9_Real; + const Real CtConst = + Eos::calcCtFreezing(EosType::ConstantEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtConst, CtConstExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Constant FAIL, expected {}, got {}", + CtConstExpected, CtConst); + } + + return; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { Eos::destroyInstance(); @@ -767,22 +803,22 @@ void checkValueGswcN2() { } /// Test that the calcCtFreezing function returns the expected value -void checkValueCtFreezing() { +void checkValueGswcCtFreezing() { const Real RTol = 1e-10; - Teos10Eos TestEos(VertCoord::getDefault()); constexpr Real SaturationFrac = 0.0; constexpr Real P = 500.0 * Db2Pa; // Convert dbar to Pa constexpr Real Sa = 32.0; /// Get freezing temperature from GSW-C library double CtFreezGswc = gsw_ct_freezing_poly(Sa, P * Pa2Db, SaturationFrac); - double CtFreez = TestEos.calcCtFreezing(Sa, P * Pa2Db, SaturationFrac); + double CtFreez = + Teos10Eos::calcCtFreezingTeos10(Sa, P * Pa2Db, SaturationFrac); /// Check the value against the GSW-C value bool Check = isApprox(CtFreezGswc, CtFreez, RTol); if (!Check) { - ABORT_ERROR("checkValueCtFreezing: CtFreez FAIL, expected {}, got {}", + ABORT_ERROR("checkValueGswcCtFreezing: CtFreez FAIL, expected {}, got {}", CtFreezGswc, CtFreez); } return; @@ -837,7 +873,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); - checkValueCtFreezing(); + checkValueGswcCtFreezing(); checkValueGswcCtFromPt(); checkValueGswcPtFromCt(); @@ -848,6 +884,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { testEosTeos10(); testEosTeos10Displaced(); testBruntVaisalaFreqSqTeos10(); + testCalcCtFreezing(); finalizeEosTest(); diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp new file mode 100644 index 000000000000..9c159073b420 --- /dev/null +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -0,0 +1,497 @@ +//===-- Test driver for OMEGA Frazil (basic option) --------------*- C++ +//-*-===// +// +/// \file +/// \brief Minimal test driver for OMEGA frazil formation functor +// +//===-----------------------------------------------------------------------===/ + +#include "Frazil.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Eos.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanTestCommon.h" +#include "OmegaKokkos.h" +#include "Pacer.h" +#include "TimeMgr.h" +#include "VertCoord.h" +#include "mpi.h" + +#include + +using namespace OMEGA; + +constexpr int NVertLayers = 60; + +void initFrazilTest(const std::string &mesh) { + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + LOG_INFO("------ Frazil Unit Tests ------"); + + Config("Omega"); + Config::readAll("omega.yml"); + + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + IO::init(DefComm); + Decomp::init(mesh); + Field::init(ModelClock); + IOStream::init(ModelClock); + Halo::init(); + HorzMesh::init(ModelClock); + VertCoord::init(false); + Eos::init(); + Frazil::init(); + + // Frazil::init only creates the default object when + // Tendencies.FrazilTendencyEnable is true in config. + // Unit tests need a frazil object regardless of that runtime switch. + if (Frazil::getDefault() == nullptr) { + Frazil *CreatedDefault = Frazil::create("Default"); + if (CreatedDefault == nullptr) { + ABORT_ERROR("FrazilTest: failed to create default Frazil object"); + } + } +} + +void finalizeFrazilTest() { + Frazil::clear(); + VertCoord::clear(); + HorzMesh::clear(); + Halo::clear(); + Decomp::clear(); + Field::clear(); + Dimension::clear(); + IOStream::finalize(); + MachEnv::removeAll(); +} + +// this test only excercises the frazil formation functor (no melt) +// in a warm case: the frazil FORMATION terms should all be zero +void testBasicFrazilFormationWarm() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = 10.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + BasicFrazilFormation ComputeFrazilFormation; + ComputeFrazilFormation.massLimit = 0.1_Real; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend, + TTend, STend, CTfrz); + + if (!isApprox(AccMIce, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccMIce, got {}", + AccMIce); + } + + if (!isApprox(AccMSalt, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccMSalt, got {}", + AccMSalt); + } + + if (!isApprox(AccEIce, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilBasicFormationTest warm: expected zero AccEIce, got {}", + AccEIce); + } + if (!isApprox(HTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero HTend, got {}", + HTend); + } + + if (!isApprox(TTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero TTend, got {}", + TTend); + } + + if (!isApprox(STend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilBasicFormationTest warm: expected zero STend, got {}", + STend); + } + LOG_INFO("FrazilBasicFormationTestWarm: AccMIce = {}, AccMLiq = {}, " + "AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +// this test only excercises the frazil formation functor (no melt) +// in a cold case: the frazil terms should be +// - strictly positive for ice, liquid, and salt mass +//- strictly negative for ice and liquid energy +// - positive for T tendency and negative for S, H tendencies +void testBasicFrazilFormationCold() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = -2.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + BasicFrazilFormation ComputeFrazilFormation; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend, + TTend, STend, CTfrz); + + if (AccMIce <= 0.0_Real) { + ABORT_ERROR("FrazilBasicFormationTestCold: accumulated ice mass is " + "non-positive: {}", + AccMIce); + } + + if (AccMSalt <= 0.0_Real) { + ABORT_ERROR("FrazilBasicFormationTestCold: accumulated salt mass is " + "non-positive: {}", + AccMSalt); + } + + if (AccEIce >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: accumulated ice energy is positive " + "(exp. negative): {}", + AccEIce); + } + if (HTend >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: HTend is positive (exp. negative): {}", + HTend); + } + if (TTend <= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: TTend is negative (exp. positive): {}", + TTend); + } + if (STend >= 0.0_Real) { + ABORT_ERROR( + "FrazilBasicFormationTestCold: STend is positive (exp. negative): {}", + STend); + } + LOG_INFO("FrazilBasicFormationTestCold: AccMIce = {}, AccMLiq = {}, " + "AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +// this test exercises the frazil formation and melt functors +// in a column of water with both cold and warm layers. +// It turns to frazil column conservation check. +// In dev, there is extra verbose logging in the frazil code (TBRemoved). +void testComputeFrazilColumn() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilTestColumn: test frazil object is null"); + } + + const Real RTol = 1e-10_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 4) { + ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - 0.5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; + deepCopy(CT, CTH); + + const bool SavedConservationCheck = TestFrazil->conservationCheck; + TestFrazil->conservationCheck = false; + TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + if (HTendH(ICell, KBottom0) >= 0.0_Real || + TTendH(ICell, KBottom0) <= 0.0_Real || + STendH(ICell, KBottom0) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: bottom cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KBottom1) >= 0.0_Real || + TTendH(ICell, KBottom1) <= 0.0_Real || + STendH(ICell, KBottom1) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: second cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KWarm) < 0.0_Real || TTendH(ICell, KWarm) > 0.0_Real || + STendH(ICell, KWarm) < 0.0_Real) { + ABORT_ERROR("FrazilTestColumn: warm layer sign check failed (HTend>=0, " + "TTend<=0, STend>=0 expected)"); + } + + if (HTendH(ICell, KTopCold) >= 0.0_Real || + TTendH(ICell, KTopCold) <= 0.0_Real || + STendH(ICell, KTopCold) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: top cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + LOG_INFO( + "FrazilTestColumn: XTend branch-switching checks passed for ICell={}", + ICell); +} + +// this test exercises the frazil formation and melt functors +// with a depth limit set. Layers deeper than the depth limit +// should have zero frazil tendencies. +void testComputeFrazilDepthLimit() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilDepthLimitTest: test frazil object is null"); + } + + const Real RTol = 1e-12_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 10) { + ABORT_ERROR("FrazilDepthLimitTest: cell {} has fewer than 10 layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - 0.5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; + deepCopy(CT, CTH); + + const bool SavedConservationCheck = TestFrazil->conservationCheck; + const Real SavedDepthLimit = TestFrazil->depthLimit; + const Real TestDepthLimit = 35.0_Real; // this needs to be positive + // if TestDepthLimit is negative, test will fail: + // - the code assume depthlimit < 0 mean no limit (i.e. full depth frazil) + // - the test below will exclude all layers and fail because Tend !=0. + + // Populate GeomZMid explicitly for the test column. + auto GeomZMidH = createHostMirrorCopy(VCoord->GeomZMid); + for (I4 K = KMin; K <= KMax; ++K) { + GeomZMidH(ICell, K) = -10.0_Real * (K - KMin + 1); + } + deepCopy(VCoord->GeomZMid, GeomZMidH); + + TestFrazil->conservationCheck = false; + TestFrazil->depthLimit = TestDepthLimit; + TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; + TestFrazil->depthLimit = SavedDepthLimit; + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + bool FoundExcludedLayer = false; + for (I4 K = KMin; K <= KMax; ++K) { + const Real Depth = GeomZMidH(ICell, K); + const Real AbsDepth = Depth < 0.0_Real ? -Depth : Depth; + + if (AbsDepth > TestDepthLimit) { + FoundExcludedLayer = true; + if (!isApprox(HTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(TTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(STendH(ICell, K), 0.0_Real, RTol)) { + ABORT_ERROR("FrazilDepthLimitTest: excluded layer K={} has " + "non-zero tendencies (HTend={}, TTend={}, STend={})", + K, HTendH(ICell, K), TTendH(ICell, K), + STendH(ICell, K)); + } + } + } + + if (!FoundExcludedLayer) { + ABORT_ERROR("FrazilDepthLimitTest: no layers excluded for ICell={} with " + "depthLimit={}", + ICell, TestDepthLimit); + } + + LOG_INFO("FrazilDepthLimitTest: depthLimit={} exclusion check passed for " + "ICell={}", + TestDepthLimit, ICell); +} + +void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { + initFrazilTest(MeshFile); + testBasicFrazilFormationCold(); + testBasicFrazilFormationWarm(); + testComputeFrazilColumn(); + testComputeFrazilDepthLimit(); + finalizeFrazilTest(); +} + +int main(int argc, char *argv[]) { + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + frazilTest(); + + LOG_INFO("------ Frazil Unit Tests Successful ------"); + + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + return 0; +} diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5268d0a436ea..ac11def2793d 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; @@ -139,6 +142,14 @@ int initTendenciesTest(const std::string &mesh) { VertCoord::init(); Tracers::init(); + Frazil::init(); + if (Frazil::getDefault() == nullptr) { + Frazil *CreatedDefault = Frazil::create("Default"); + if (CreatedDefault == nullptr) { + Err++; + LOG_ERROR("TendenciesTest: failed to create default Frazil object"); + } + } VertAdv::init(); PressureGrad::init(); Eos::init(); @@ -303,11 +314,162 @@ int testTendencies() { "NormalVelocityTend"); } + const Real NormVelTendSum = + sum(DefTendencies->NormalVelocityTend, Mesh->NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); + if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { + Err++; + LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); + } + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + // Verify frazil tendencies are plumbed into pseudo-thickness and tracer + // tendencies by comparing runs with FrazilTerm disabled/enabled. + I4 TempTracerIndex = -1; + I4 SaltTracerIndex = -1; + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + if (TempTracerIndex < 0 || SaltTracerIndex < 0) { + Err++; + LOG_ERROR("TendenciesTest: missing Temperature/Salinity tracer indices"); + } else { + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + auto TracerArrayH = createHostMirrorCopy(TracerArray); + + // Set Salinity in active layers to [31, 33], top-to-bottom ramp. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + for (I4 K = KMin; K <= KMax; ++K) { + Real frac = 0.0_Real; + if (KMax > KMin) { + frac = + static_cast(K - KMin) / static_cast(KMax - KMin); + } + TracerArrayH(SaltTracerIndex, ICell, K) = + 31.0_Real + 2.0_Real * frac; + } + } + + // Keep Temperature warm everywhere first, then add one cold point. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + for (I4 K = KMin; K <= KMax; ++K) { + TracerArrayH(TempTracerIndex, ICell, K) = 8.0_Real; + } + } + const I4 FrazilTestCell = 0; + const I4 FrazilTestLayer = MinLayerCellH(FrazilTestCell); + TracerArrayH(TempTracerIndex, FrazilTestCell, FrazilTestLayer) = + -2.0_Real; + deepCopy(TracerArray, TracerArrayH); + + Array2DReal BaselinePseudoThicknessTend( + "BaselinePseudoThicknessTend", Mesh->NCellsSize, VCoord->NVertLayers); + Array3DReal BaselineTracerTend("BaselineTracerTend", + Tracers::getNumTracers(), Mesh->NCellsSize, + VCoord->NVertLayers); + + const bool OrigFrazilEnabled = DefTendencies->FrazilTerm.Enabled; + + DefTendencies->FrazilTerm.Enabled = false; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + deepCopy(BaselinePseudoThicknessTend, DefTendencies->PseudoThicknessTend); + deepCopy(BaselineTracerTend, DefTendencies->TracerTend); + + deepCopy(DefTendencies->PseudoThicknessTend, 0.0_Real); + deepCopy(DefTendencies->NormalVelocityTend, 0.0_Real); + deepCopy(DefTendencies->TracerTend, 0.0_Real); + + DefTendencies->FrazilTerm.Enabled = true; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + Array2DReal PseudoThicknessTendDiff( + "PseudoThicknessTendDiff", Mesh->NCellsSize, VCoord->NVertLayers); + Array2DReal TempTracerTendDiff("TempTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + Array2DReal SaltTracerTendDiff("SaltTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + + OMEGA_SCOPE(LocPseudoThicknessTendDiff, PseudoThicknessTendDiff); + OMEGA_SCOPE(LocTempTracerTendDiff, TempTracerTendDiff); + OMEGA_SCOPE(LocSaltTracerTendDiff, SaltTracerTendDiff); + OMEGA_SCOPE(LocPseudoThicknessTend, DefTendencies->PseudoThicknessTend); + OMEGA_SCOPE(LocBaselinePseudoThicknessTend, BaselinePseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, DefTendencies->TracerTend); + OMEGA_SCOPE(LocBaselineTracerTend, BaselineTracerTend); + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + + parallelForOuter( + "TendenciesTest:FrazilTendDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTendDiff(ICell, K) = + Kokkos::abs(LocPseudoThicknessTend(ICell, K) - + LocBaselinePseudoThicknessTend(ICell, K)); + LocTempTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(TempTracerIndex, ICell, K) - + LocBaselineTracerTend(TempTracerIndex, ICell, K)); + LocSaltTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(SaltTracerIndex, ICell, K) - + LocBaselineTracerTend(SaltTracerIndex, ICell, K)); + }); + }); + + const Real PseudoDelta = sum(PseudoThicknessTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real TempDelta = sum(TempTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real SaltDelta = sum(SaltTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + + constexpr Real FrazilDeltaATol = 1e-12_Real; + if (!Kokkos::isfinite(PseudoDelta) || + isApprox(PseudoDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change " + "PseudoThicknessTend"); + } + if (!Kokkos::isfinite(TempDelta) || + isApprox(TempDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Temperature tracer " + "tendency"); + } + if (!Kokkos::isfinite(SaltDelta) || + isApprox(SaltDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Salinity tracer " + "tendency"); + } + + DefTendencies->FrazilTerm.Enabled = OrigFrazilEnabled; + } + + // Test surface tracer forcing with enthalpy terms (TEOS-10 CtFrz path) + const int TracerForcingErr = testSfcTracerForcing(); + Err += TracerForcingErr; + + // Test surface thickness forcing with freshwater terms + const int ThicknessForcingErr = testSfcThicknessForcing(); + Err += ThicknessForcingErr; + // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; - int NEdgesOwned = Mesh->NEdgesOwned; int NTracers = Tracers::getNumTracers(); const Real PseudoThickTendSum = @@ -318,23 +480,391 @@ int testTendencies() { LOG_ERROR("TendenciesTest: PseudoThickTend FAIL"); } - const Real NormVelTendSum = - sum(DefTendencies->NormalVelocityTend, NEdgesOwned, - VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); - if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { - Err++; - LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); - } - - const Real TraceTendSum = + const Real TracerTendSum = sum(DefTendencies->TracerTend, NTracers, NCellsOwned, VCoord->MinLayerCell, VCoord->MaxLayerCell); - if (!Kokkos::isfinite(TraceTendSum) || TraceTendSum == 0) { + if (!Kokkos::isfinite(TracerTendSum) || TracerTendSum == 0) { Err++; - LOG_ERROR("TendenciesTest: TraceTendSum FAIL"); + LOG_ERROR("TendenciesTest: TracerTendSum FAIL"); } Tendencies::clear(); + 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->MinLayerCellH(ICellTest); + + if (KTop > VCoord->MaxLayerCellH(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 a reference expectations for temperature tendency: + // using TEOS-10 freezing CT (expected to pass under strict tolerance). + + HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); + deepCopy(PressureMidH, VCoord->PressureMid); + const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; + const Real CtFrz = + Eos::calcCtFreezing(EosInst->EosChoice, SaTopValue, PTopDb, 0.0_Real); + const Real ExpectedTempTend = + (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + + TestSnow * (Cp0Sw * CtFrz - 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-pass check with TEOS freezing CT reference. + if (!isApprox(ComputedTempTend, ExpectedTempTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); + LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", + ExpectedTempTend, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTend)); + } 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_INFO(" 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->MinLayerCellH(ICellTest); + + if (KTop > VCoord->MaxLayerCellH(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_INFO(" 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; } @@ -342,6 +872,7 @@ int testTendencies() { void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); + Frazil::clear(); PressureGrad::clear(); VertMix::destroyInstance(); Eos::destroyInstance(); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 50b4c12b3eda..9cf3737da573 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -254,6 +254,8 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerDiffusion.Enabled = false; TestTendencies->TracerHyperDiff.Enabled = false; TestTendencies->SfcStressForcing.Enabled = false; + TestTendencies->SfcTracerForcing.Enabled = false; + TestTendencies->SfcThicknessForcing.Enabled = false; TestTendencies->SurfaceTracerRestoring.Enabled = false; TestTendencies->ExplicitBottomDrag.Enabled = false; DefVAdv->ThickVertAdvEnabled = false;