Add submesoscale eddy parametrization

components/omega/configs/Default.yml

@@ -89,6 +89,12 @@ Omega:
       NuZero: 0.005
       Alpha: 5.0
       Exponent: 2.0
+  Submeso:
+    Enable: false
+    Tau: 172800.0
+    Ce: 0.08
+    LfMin: 1.0e3
+    DsMax: 100.0e3
   IOStreams:
     HorzMeshIn:
       UsePointerFile: false

components/omega/src/ocn/AuxiliaryState.cpp

@@ -3,6 +3,7 @@
 #include "Field.h"
 #include "Logging.h"
 #include "Pacer.h"
+#include "SubmesoEddies.h"
 #include "Tendencies.h"
 #include "TimeStepper.h"
 
@@ -29,7 +30,8 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
       VelocityDel2Aux(stripDefault(Name), Mesh, VCoord),
       WindForcingAux(stripDefault(Name), Mesh),
       SurfTracerRestAux(stripDefault(Name), Mesh, NTracers),
-      TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) {
+      TracerAux(stripDefault(Name), Mesh, VCoord, NTracers),
+      TransportAux(stripDefault(Name), Mesh, VCoord) {
 
    GroupName = "AuxiliaryState";
    if (Name != "Default") {
@@ -46,6 +48,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
    WindForcingAux.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
@@ -58,6 +61,7 @@ AuxiliaryState::~AuxiliaryState() {
    WindForcingAux.unregisterFields();
    SurfTracerRestAux.unregisterFields();
    TracerAux.unregisterFields();
+   TransportAux.unregisterFields();
 
    FieldGroup::destroy(GroupName);
 }
@@ -99,12 +103,94 @@ 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->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;
+
+      SubEddies->computeDenMixLayerDepth(SpecVol);
+      SubEddies->computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid,
+                                 BVFreqSq);
+      SubEddies->computeEddyVelocity(BVFreqSq, MeanPseudoThickEdge);
+
+      const auto &EddyVelocity = SubEddies->EddyVelocity;
+
+      parallelFor(
+          {Mesh->NEdgesAll, VCoord->NVertLayers},
+          KOKKOS_LAMBDA(int IEdge, 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);
+              });
+       });
+
+   if (SubmesoEddies::getInstance()->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,
@@ -197,6 +283,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State,
               });
        });
 
+   computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel);
+
    parallelForOuter(
        "edgeAuxState2", {Mesh->NEdgesAll},
        KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
@@ -243,15 +331,55 @@ 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);
+
+   if (SubmesoEddies::getInstance()->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"
 #include "auxiliaryVars/WindForcingAuxVars.h"
@@ -44,6 +45,7 @@ class AuxiliaryState {
    VelocityDel2AuxVars VelocityDel2Aux;
    WindForcingAuxVars WindForcingAux;
    SurfTracerRestAuxVars SurfTracerRestAux;
+   TransportAuxVars TransportAux;
 
    ~AuxiliaryState();
 
@@ -76,16 +78,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,

components/omega/src/ocn/OceanFinal.cpp

@@ -17,6 +17,7 @@
 #include "OceanDriver.h"
 #include "OceanState.h"
 #include "PGrad.h"
+#include "SubmesoEddies.h"
 #include "Tendencies.h"
 #include "TimeMgr.h"
 #include "TimeStepper.h"
@@ -41,6 +42,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time
 
    Tracers::clear();
    TimeStepper::clear();
+   SubmesoEddies::destroyInstance();
    PressureGrad::clear();
    Eos::destroyInstance();
    Tendencies::clear();

components/omega/src/ocn/OceanInit.cpp

@@ -24,6 +24,7 @@
 #include "OceanState.h"
 #include "PGrad.h"
 #include "Pacer.h"
+#include "SubmesoEddies.h"
 #include "Tendencies.h"
 #include "TimeMgr.h"
 #include "TimeStepper.h"
@@ -195,6 +196,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) {
    AuxiliaryState::init();
    Eos::init();
    PressureGrad::init();
+   SubmesoEddies::init();
    Tendencies::init();
 
    // Validate SurfaceTracerRestoring configuration