Add submesoscale eddy parametrization

components/omega/configs/Default.yml

@@ -98,6 +98,12 @@ Omega:
       RiCrit: 0.7
       Exponent: 3.0
       RiSmoothLoops: 2
+  Submeso:
+    Enable: false
+    Tau: 172800.0
+    Ce: 0.08
+    LfMin: 1.0e3
+    DsMax: 100.0e3
   IOStreams:
     HorzMeshIn:
       UsePointerFile: false

components/omega/doc/devGuide/SubmesoEddies.md

@@ -0,0 +1,73 @@
+(omega-dev-submeso-eddies)=
+
+# Submesoscale Eddy Parameterization
+
+Omega includes a `SubmesoEddies` class (`src/ocn/SubmesoEddies.{h,cpp}`) that implements the
+submesoscale mixed-layer instability closure used to produce an eddy-induced
+normal transport velocity. It provides methods to compute mixed-layer depth from
+a density threshold, buoyancy gradients, and eddy velocity.
+The current implementation follows the Fox-Kemper et al. (2011) (FK11) closure,
+as described in the
+{ref}`omega-design-submesoscale-eddies` design document.
+
+## Initialization
+
+The `SubmesoEddies` class is implemented as a singleton. Before creating it,
+[`HorzMesh`](#omega-dev-horz-mesh) and [`VertCoord`](#omega-dev-vert-coord)
+must be initialized. Create the instance with the static method
+```c++
+SubmesoEddies::init();
+```
+Retrieve the pointer at any time with:
+```c++
+SubmesoEddies* DefSubEddies = SubmesoEddies::getInstance();
+```
+
+## Data and algorithms
+
+Core public fields include:
+
+- `DenMixLayerIndex`, `DenMixLayerDepth`
+- `GradBuoyEdgeInterface`
+- `EddyVelocity`
+
+Key implementation details:
+
+- `computeDenMixLayerDepth` uses a density-threshold criterion referenced to a
+  fixed near-surface depth (`ReferenceDepth = 10 m`) and linear interpolation to
+  estimate the crossing depth.
+- `computeBuoyGrad` computes horizontal buoyancy gradients at edges and adds
+  the tilted-coordinate correction using `BruntVaisalaFreqSq`.
+- `computeEddyVelocity` forms mixed layer averaged buoyancy and stratification
+  terms, evaluates frontal-width limits (`LfMin`, `DsMax`), computes a
+  streamfunction with `shapeFunction()`, and takes its vertical divergence to
+  produce edge-normal eddy velocity.
+
+## Computation of mixed layer depth
+To compute the mixed layer depth `DenMixedLayerDepth` from the specific volume `SpecVol`, do
+```c++
+SubEddies.computeDenMixLayerDepth(SpecVol);
+```
+
+## Computation of buoyancy gradient
+To compute buoyancy gradient `GradBuoyEdgeInterface` from specific volume
+`SpecVol`, mean pseudo-thickness on edges `MeanPseudoThickEdge`, mid-layer
+`z` coordinate `GeomZMid`, and squared Brunt-Vaisala frequency
+`BruntVaisalaFreqSq`, use
+```c++
+SubEddies.computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid, BruntVaisalaFreqSq);
+```
+
+## Computation of eddy velocity
+To compute eddy velocity array `EddyVelocity` from squared Brunt-Vaisala
+frequency field `BruntVaisalaFreqSq` and mean pseudo-thickness on edges
+`MeanPseudoThickEdge`, use
+```c++
+SubEddies.computeEddyVelocity(BruntVaisalaFreqSq, MeanPseudoThickEdge);
+```
+
+## Finalization
+To clear the singleton instance, use the static method
+```c++
+SubmesoEddies::destroyInstance();
+```

components/omega/doc/index.md

@@ -53,6 +53,7 @@ userGuide/VertCoord
 userGuide/PGrad
 userGuide/Timing
 userGuide/VerticalMixingCoeff
+userGuide/SubmesoEddies
 userGuide/VertAdv
 userGuide/Forcing
 userGuide/SfcCoupling
@@ -102,6 +103,7 @@ devGuide/VertCoord
 devGuide/PGrad
 devGuide/Timing
 devGuide/VerticalMixingCoeff
+devGuide/SubmesoEddies
 devGuide/VertAdv
 devGuide/Forcing
 devGuide/SfcCoupling

components/omega/doc/userGuide/SubmesoEddies.md

@@ -0,0 +1,39 @@
+(omega-user-submeso-eddies)=
+
+# Submesoscale Eddy Parameterization
+
+Omega includes an optional submesoscale mixed layer instability (MLI)
+parameterization through the `SubmesoEddies` class. When enabled, the model
+computes an eddy-induced transport velocity on edges and adds it to the normal
+transport velocity used by thickness and tracer advection.
+
+The current implementation follows the Fox-Kemper et al. (2011) (FK11) closure
+described in the
+{ref}`omega-design-submesoscale-eddies` design document.
+
+## Configuration
+
+Configure the parameterization in the `Submeso` section of the YAML input:
+
+```yaml
+Submeso:
+  Enable: false
+  Tau: 172800.0
+  Ce: 0.08
+  LfMin: 1.0e3
+  DsMax: 100.0e3
+```
+
+- `Enable`: turns the parameterization on/off.
+- `Tau`: MLI timescale parameter (s).
+- `Ce`: nondimensional efficiency coefficient.
+- `LfMin`: minimum frontal width limiter (m).
+- `DsMax`: maximum grid-length limiter used in the closure (m).
+
+## Diagnostics
+
+When enabled, the following fields are available in the `Submeso` field group:
+
+- `DenMixLayerDepth` (m): density-threshold mixed-layer depth.
+- `GradBuoyEdgeInterface` (s^-2): buoyancy gradient on edge interfaces.
+- `EddyVelocity` (m/s): eddy-induced transport velocity.

components/omega/src/ocn/AuxiliaryState.cpp

@@ -4,6 +4,7 @@
 #include "Field.h"
 #include "Logging.h"
 #include "Pacer.h"
+#include "SubmesoEddies.h"
 #include "Tendencies.h"
 #include "TimeStepper.h"
 
@@ -30,7 +31,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
       VelocityDel2Aux(stripDefault(Name), Mesh, VCoord),
       SurfTracerRestAux(stripDefault(Name), Mesh, NTracers),
       TracerAux(stripDefault(Name), Mesh, VCoord, NTracers),
-      TimeStep(TimeStep) {
+      TransportAux(stripDefault(Name), Mesh, VCoord), TimeStep(TimeStep) {
 
    GroupName = "AuxiliaryState";
    if (Name != "Default") {
@@ -46,6 +47,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
    VelocityDel2Aux.registerFields(GroupName, AuxMeshName);
    SurfTracerRestAux.registerFields(GroupName, AuxMeshName);
    TracerAux.registerFields(GroupName, AuxMeshName);
+   TransportAux.registerFields(GroupName, AuxMeshName);
 }
 
 // Destructor. Unregisters the fields with IOStreams and destroys this auxiliary
@@ -57,6 +59,7 @@ AuxiliaryState::~AuxiliaryState() {
    VelocityDel2Aux.unregisterFields();
    SurfTracerRestAux.unregisterFields();
    TracerAux.unregisterFields();
+   TransportAux.unregisterFields();
 
    FieldGroup::destroy(GroupName);
 }
@@ -98,12 +101,104 @@ void AuxiliaryState::computeMomVertAux(const OceanState *State,
    // compute geometric height
    VCoord->computeGeomZHeight(PseudoThickCell, EosInstance->SpecVol);
 
+   // compute Brunt-Vaisala freqency squared
+   EosInstance->computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, PressureMid,
+                                          EosInstance->SpecVol);
+
    // compute target thickness
    VCoord->computeTargetThickness();
 
    Pacer::stop("AuxState:computeMomVertAux", 2);
 }
 
+// Compute transport velocity for pseudo-thickness and tracers
+void AuxiliaryState::computeTransportVelocity(const OceanState *State,
+                                              const Array3DReal &TracerArray,
+                                              int ThickTimeLevel,
+                                              int VelTimeLevel) const {
+   Pacer::start("AuxState:computeTransportVelocity", 2);
+
+   Array2DReal NormalVel = State->getNormalVelocity(VelTimeLevel);
+
+   const auto &NormalTransportVelocity = TransportAux.NormalTransportVelocity;
+
+   deepCopy(NormalTransportVelocity, NormalVel);
+
+   auto *SubEddies = SubmesoEddies::getInstance();
+
+   if (SubEddies && SubEddies->Enable) {
+
+      Eos *EosInstance = Eos::getInstance();
+
+      const auto &MeanPseudoThickEdge = PseudoThicknessAux.MeanPseudoThickEdge;
+      const auto &SpecVol             = EosInstance->SpecVol;
+      const auto &BVFreqSq            = EosInstance->BruntVaisalaFreqSq;
+      const auto &GeomZMid            = VCoord->GeomZMid;
+      const auto &MinLayerEdgeBot     = VCoord->MinLayerEdgeBot;
+      const auto &MaxLayerEdgeTop     = VCoord->MaxLayerEdgeTop;
+
+      SubEddies->computeDenMixLayerDepth(SpecVol);
+      SubEddies->computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid,
+                                 BVFreqSq);
+      SubEddies->computeEddyVelocity(BVFreqSq, MeanPseudoThickEdge);
+
+      const auto &EddyVelocity = SubEddies->EddyVelocity;
+
+      parallelForOuter(
+          {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+             const int MinLyrEdgeBot = MinLayerEdgeBot(IEdge);
+             const int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge);
+
+             parallelForInner(
+                 Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop},
+                 INNER_LAMBDA(int K) {
+                    NormalTransportVelocity(IEdge, K) += EddyVelocity(IEdge, K);
+                 });
+          });
+   }
+
+   Pacer::stop("AuxState:computeTransportVelocity", 2);
+}
+
+// Compute the auxiliary variables needed for pseudo-thickness equation
+void AuxiliaryState::computePseudoThicknessAux(const OceanState *State,
+                                               const Array3DReal &TracerArray,
+                                               int ThickTimeLevel,
+                                               int VelTimeLevel) const {
+
+   Array2DReal PseudoThick   = State->getPseudoThickness(ThickTimeLevel);
+   Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel);
+   OMEGA_SCOPE(LocPseudoThicknessAux, PseudoThicknessAux);
+   OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot);
+   OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop);
+
+   Pacer::start("Tend:computePseudoThickAux", 2);
+
+   parallelForOuter(
+       "computePseudoThickAux", {Mesh->NEdgesAll},
+       KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+          const int KMin   = MinLayerEdgeBot(IEdge);
+          const int KMax   = MaxLayerEdgeTop(IEdge);
+          const int KRange = vertRangeChunked(KMin, KMax);
+
+          parallelForInner(
+              Team, KRange, INNER_LAMBDA(int KChunk) {
+                 LocPseudoThicknessAux.computeVarsOnEdge(
+                     IEdge, KChunk, PseudoThick, NormalVelEdge);
+              });
+       });
+
+   auto *SubEddies = SubmesoEddies::getInstance();
+
+   if (SubEddies && SubEddies->Enable) {
+      computeMomVertAux(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+   }
+
+   computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+
+   Pacer::stop("Tend:computePseudoThickAux", 2);
+}
+
 // Compute the auxiliary variables needed for momentum equation
 void AuxiliaryState::computeMomAux(const OceanState *State,
                                    const Array3DReal &TracerArray,
@@ -196,6 +291,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
               });
        });
 
+   computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+
    parallelForOuter(
        "edgeAuxState2", {Mesh->NEdgesAll},
        KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
@@ -246,15 +343,57 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
 
    Pacer::start("AuxState:computeVerticalPseudoVelocity", 2);
 
-   const auto &FluxPseudoThickEdge = PseudoThicknessAux.FluxPseudoThickEdge;
-   VAdv->computeVerticalPseudoVelocity(NormalVelEdge, FluxPseudoThickEdge,
-                                       PseudoThickCell, ProjDtSeconds);
+   const auto &FluxPseudoThickEdge     = PseudoThicknessAux.FluxPseudoThickEdge;
+   const auto &NormalTransportVelocity = TransportAux.NormalTransportVelocity;
+   VAdv->computeVerticalPseudoVelocity(NormalTransportVelocity,
+                                       FluxPseudoThickEdge, PseudoThickCell,
+                                       ProjDtSeconds);
 
    Pacer::stop("AuxState:computeVerticalPseudoVelocity", 2);
 
    Pacer::stop("AuxState:computeMomAux", 1);
 }
 
+// Compute the auxiliary variables needed for tracer equation
+void AuxiliaryState::computeTracerAux(const OceanState *State,
+                                      const Array3DReal &TracerArray,
+                                      int ThickTimeLevel,
+                                      int VelTimeLevel) const {
+
+   OMEGA_SCOPE(LocTracerAux, TracerAux);
+   OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell);
+   OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell);
+
+   const auto &MeanPseudoThickEdge = PseudoThicknessAux.MeanPseudoThickEdge;
+
+   const int NTracers = Tracers::getNumTracers();
+
+   Pacer::start("Tend:computeTracerAuxCell", 2);
+
+   auto *SubEddies = SubmesoEddies::getInstance();
+
+   if (SubEddies && SubEddies->Enable) {
+      computeMomVertAux(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+   }
+
+   computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+
+   parallelForOuter(
+       "computeTracerAuxCell", {NTracers, Mesh->NCellsAll},
+       KOKKOS_LAMBDA(int LTracer, int ICell, const TeamMember &Team) {
+          const int KMin   = MinLayerCell(ICell);
+          const int KMax   = MaxLayerCell(ICell);
+          const int KRange = vertRangeChunked(KMin, KMax);
+
+          parallelForInner(
+              Team, KRange, INNER_LAMBDA(int KChunk) {
+                 LocTracerAux.computeVarsOnCells(
+                     LTracer, ICell, KChunk, MeanPseudoThickEdge, TracerArray);
+              });
+       });
+   Pacer::stop("Tend:computeTracerAuxCell", 2);
+}
+
 // Compute the auxiliary variables
 void AuxiliaryState::computeAll(const OceanState *State,
                                 const Array3DReal &TracerArray,

components/omega/src/ocn/AuxiliaryState.h

@@ -14,6 +14,7 @@
 #include "auxiliaryVars/PseudoThicknessAuxVars.h"
 #include "auxiliaryVars/SurfTracerRestAuxVars.h"
 #include "auxiliaryVars/TracerAuxVars.h"
+#include "auxiliaryVars/TransportAuxVars.h"
 #include "auxiliaryVars/VelocityDel2AuxVars.h"
 #include "auxiliaryVars/VorticityAuxVars.h"
 
@@ -42,6 +43,7 @@ class AuxiliaryState {
    VorticityAuxVars VorticityAux;
    VelocityDel2AuxVars VelocityDel2Aux;
    SurfTracerRestAuxVars SurfTracerRestAux;
+   TransportAuxVars TransportAux;
 
    ~AuxiliaryState();
 
@@ -74,16 +76,31 @@ class AuxiliaryState {
    /// Exchange halo
    I4 exchangeHalo();
 
+   // Compute all auxiliary variables needed for pseudo-thickness equation
+   void computePseudoThicknessAux(const OceanState *State,
+                                  const Array3DReal &TracerArray,
+                                  int ThickTimeLevel, int VelTimeLevel) const;
+
    // Compute auxiliary variables for vertical dynamics
    void computeMomVertAux(const OceanState *State,
                           const Array3DReal &TracerArray, int ThickTimeLevel,
                           int VelTimeLevel) const;
 
+   // Compute transport velocity for pseudo-thickness and tracers
+   void computeTransportVelocity(const OceanState *State,
+                                 const Array3DReal &TracerArray,
+                                 int ThickTimeLevel, int VelTimeLevel) const;
+
    // Compute all auxiliary variables needed for momentum equation
    void computeMomAux(const OceanState *State, const Array3DReal &TracerArray,
                       int ThickTimeLevel, int VelTimeLevel,
                       const TimeInterval ProjDt) const;
 
+   // Compute all auxiliary variables needed for tracer equation
+   void computeTracerAux(const OceanState *State,
+                         const Array3DReal &TracerArray, int ThickTimeLevel,
+                         int VelTimeLevel) const;
+
    /// Compute all auxiliary variables based on an ocean state at a given time
    /// level
    void computeAll(const OceanState *State, const Array3DReal &TracerArray,