Add wind forcing and bottom drag tendency terms

components/omega/configs/Default.yml

@@ -26,6 +26,9 @@ Omega:
     VelHyperDiffTendencyEnable: true
     ViscDel4: 1.2e11
     DivFactor: 1.0
+    WindForcingTendencyEnable: false
+    BottomDragTendencyEnable: false
+    BottomDragCoeff: 0.0
     TracerHorzAdvTendencyEnable: true
     TracerDiffTendencyEnable: true
     EddyDiff2: 10.0
@@ -91,7 +94,7 @@ Omega:
       Contents:
         - Tracers
         - State
-        - SshCellDefault
+        - SshCell
     Highfreq:
       UsePointerFile: false
       Filename: ocn.hifreq.$Y-$M

components/omega/doc/devGuide/AuxiliaryVariables.md

@@ -101,3 +101,6 @@ The following auxiliary variable groups are currently implemented:
 ||  Del2RelVortVertex ||
 | TracerAuxVars | HTracersEdge | Center or Upwind|
 || Del2TracersCell ||
+| WindForcingAuxVars |  ZonalStressEdge ||
+||  MeridStressCell ||
+||  NormalStressCell ||

components/omega/doc/devGuide/TendencyTerms.md

@@ -38,3 +38,5 @@ implemented:
 - `TracerHorzAdvOnCell`
 - `TracerDiffOnCell`
 - `TracerHyperDiffOnCell`
+- `WindForcingOnEdge`
+- `BottomDragOnEdge`

components/omega/doc/userGuide/AuxiliaryVariables.md

@@ -29,3 +29,6 @@ The following auxiliary variables are currently available:
 | VelDel2RelVortVertex | curl of laplacian of horizontal velocity on cells
 | HTracersEdge | thickness-weighted tracers used for fluxes through edges. May be centered, upwinded or a combination of the two
 | Del2TracersCell | laplacian of tracers on cells
+| ZonalStressEdge | zonal component of wind stress on cells
+| MeridStressEdge | meridional component of wind stress on cells
+| NormalStressEdge | normal component of wind stress on edge

components/omega/doc/userGuide/TendencyTerms.md

@@ -17,6 +17,8 @@ tendency terms are currently implemented:
 | TracerHorzAdvOnCell | horizontal advection of thickness-weighted tracers
 | TracerDiffOnCell | horizontal diffusion of thickness-weighted tracers
 | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers
+| WindForcingOnEdge | forcing by wind stress, defined on edges
+| BottomDragOnEdge | bottom drag, defined on edges
 
 Among the internal data stored by each functor is a `bool` which can enable or
 disable the contribution of that particular term to the tendency. These flags
@@ -45,3 +47,6 @@ the currently available tendency terms:
 | | EddyDiff2 | horizontal diffusion coefficient
 | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term
 | | EddyDiff4 | biharmonic horizontal mixing coeffienct for tracers
+| WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term
+| BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term
+| | BottomDragCoeff | bottom drag coefficient

components/omega/src/ocn/AuxiliaryState.cpp

@@ -10,15 +10,21 @@ AuxiliaryState *AuxiliaryState::DefaultAuxState = nullptr;
 std::map<std::string, std::unique_ptr<AuxiliaryState>>
     AuxiliaryState::AllAuxStates;
 
+static std::string stripDefault(const std::string &Name) {
+   return Name != "Default" ? Name : "";
+}
+
 // Constructor. Constructs the member auxiliary variables and registers their
 // fields with IOStreams
 AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
-                               int NVertLevels, int NTracers)
-    : Mesh(Mesh), Name(Name), KineticAux(Name, Mesh, NVertLevels),
-      LayerThicknessAux(Name, Mesh, NVertLevels),
-      VorticityAux(Name, Mesh, NVertLevels),
-      VelocityDel2Aux(Name, Mesh, NVertLevels),
-      TracerAux(Name, Mesh, NVertLevels, NTracers) {
+                               Halo *MeshHalo, int NVertLevels, int NTracers)
+    : Mesh(Mesh), MeshHalo(MeshHalo), Name(stripDefault(Name)),
+      KineticAux(stripDefault(Name), Mesh, NVertLevels),
+      LayerThicknessAux(stripDefault(Name), Mesh, NVertLevels),
+      VorticityAux(stripDefault(Name), Mesh, NVertLevels),
+      VelocityDel2Aux(stripDefault(Name), Mesh, NVertLevels),
+      WindForcingAux(stripDefault(Name), Mesh, NVertLevels),
+      TracerAux(stripDefault(Name), Mesh, NVertLevels, NTracers) {
 
    GroupName = "AuxiliaryState";
    if (Name != "Default") {
@@ -32,6 +38,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
    LayerThicknessAux.registerFields(GroupName, AuxMeshName);
    VorticityAux.registerFields(GroupName, AuxMeshName);
    VelocityDel2Aux.registerFields(GroupName, AuxMeshName);
+   WindForcingAux.registerFields(GroupName, AuxMeshName);
    TracerAux.registerFields(GroupName, AuxMeshName);
 }
 
@@ -42,6 +49,7 @@ AuxiliaryState::~AuxiliaryState() {
    LayerThicknessAux.unregisterFields();
    VorticityAux.unregisterFields();
    VelocityDel2Aux.unregisterFields();
+   WindForcingAux.unregisterFields();
    TracerAux.unregisterFields();
 
    int Err = FieldGroup::destroy(GroupName);
@@ -64,6 +72,7 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel,
    OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux);
    OMEGA_SCOPE(LocVorticityAux, VorticityAux);
    OMEGA_SCOPE(LocVelocityDel2Aux, VelocityDel2Aux);
+   OMEGA_SCOPE(LocWindForcingAux, WindForcingAux);
 
    parallelFor(
        "vertexAuxState1", {Mesh->NVerticesAll, NChunks},
@@ -81,6 +90,10 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel,
    const auto &VelocityDivCell = KineticAux.VelocityDivCell;
    const auto &RelVortVertex   = VorticityAux.RelVortVertex;
 
+   auto MeshHalo = Halo::getDefault();
+   MeshHalo->exchangeFullArrayHalo(LocWindForcingAux.ZonalStressCell, OnCell);
+   MeshHalo->exchangeFullArrayHalo(LocWindForcingAux.MeridStressCell, OnCell);
+
    parallelFor(
        "edgeAuxState1", {Mesh->NEdgesAll, NChunks},
        KOKKOS_LAMBDA(int IEdge, int KChunk) {
@@ -89,6 +102,7 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel,
                                                  NormalVelEdge);
           LocVelocityDel2Aux.computeVarsOnEdge(IEdge, KChunk, VelocityDivCell,
                                                RelVortVertex);
+          LocWindForcingAux.computeVarsOnEdge(IEdge, KChunk);
        });
 
    parallelFor(
@@ -153,16 +167,17 @@ void AuxiliaryState::computeAll(const OceanState *State,
 
 // Create a non-default auxiliary state
 AuxiliaryState *AuxiliaryState::create(const std::string &Name,
-                                       const HorzMesh *Mesh, int NVertLevels,
-                                       const int NTracers) {
+                                       const HorzMesh *Mesh, Halo *MeshHalo,
+                                       int NVertLevels, const int NTracers) {
    if (AllAuxStates.find(Name) != AllAuxStates.end()) {
       LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it "
                 "already exists",
                 Name);
       return nullptr;
    }
 
-   auto *NewAuxState = new AuxiliaryState(Name, Mesh, NVertLevels, NTracers);
+   auto *NewAuxState =
+       new AuxiliaryState(Name, Mesh, MeshHalo, NVertLevels, NTracers);
    AllAuxStates.emplace(Name, NewAuxState);
 
    return NewAuxState;
@@ -173,12 +188,13 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name,
 int AuxiliaryState::init() {
    int Err                 = 0;
    const HorzMesh *DefMesh = HorzMesh::getDefault();
+   Halo *DefHalo           = Halo::getDefault();
 
    int NVertLevels = DefMesh->NVertLevels;
    int NTracers    = Tracers::getNumTracers();
 
-   AuxiliaryState::DefaultAuxState =
-       AuxiliaryState::create("Default", DefMesh, NVertLevels, NTracers);
+   AuxiliaryState::DefaultAuxState = AuxiliaryState::create(
+       "Default", DefMesh, DefHalo, NVertLevels, NTracers);
 
    Config *OmegaConfig = Config::getOmegaConfig();
    Err                 = DefaultAuxState->readConfigOptions(OmegaConfig);
@@ -267,4 +283,19 @@ int AuxiliaryState::readConfigOptions(Config *OmegaConfig) {
    return Err;
 }
 
+//------------------------------------------------------------------------------
+// Perform auxiliary state halo exchange
+// Note that only non-computed auxiliary variables needs to be exchanged
+I4 AuxiliaryState::exchangeHalo() {
+   I4 Err = 0;
+
+   Err +=
+       MeshHalo->exchangeFullArrayHalo(WindForcingAux.ZonalStressCell, OnCell);
+   Err +=
+       MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell);
+
+   return Err;
+
+} // end exchangeHalo
+
 } // namespace OMEGA

components/omega/src/ocn/AuxiliaryState.h

@@ -3,6 +3,7 @@
 
 #include "Config.h"
 #include "DataTypes.h"
+#include "Halo.h"
 #include "HorzMesh.h"
 #include "OceanState.h"
 #include "Tracers.h"
@@ -11,6 +12,7 @@
 #include "auxiliaryVars/TracerAuxVars.h"
 #include "auxiliaryVars/VelocityDel2AuxVars.h"
 #include "auxiliaryVars/VorticityAuxVars.h"
+#include "auxiliaryVars/WindForcingAuxVars.h"
 
 #include <memory>
 #include <string>
@@ -36,6 +38,7 @@ class AuxiliaryState {
    TracerAuxVars TracerAux;
    VorticityAuxVars VorticityAux;
    VelocityDel2AuxVars VelocityDel2Aux;
+   WindForcingAuxVars WindForcingAux;
 
    ~AuxiliaryState();
 
@@ -46,7 +49,7 @@ class AuxiliaryState {
 
    // Create a non-default auxiliary state
    static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh,
-                                 int NVertLevels, int NTracers);
+                                 Halo *MeshHalo, int NVertLevels, int NTracers);
 
    /// Get the default auxiliary state
    static AuxiliaryState *getDefault();
@@ -63,6 +66,9 @@ class AuxiliaryState {
    /// Read and set config options
    int readConfigOptions(Config *OmegaConfig);
 
+   /// Exchange halo
+   I4 exchangeHalo();
+
    // Compute all auxiliary variables needed for momentum equation
    void computeMomAux(const OceanState *State, int ThickTimeLevel,
                       int VelTimeLevel) const;
@@ -75,13 +81,14 @@ class AuxiliaryState {
                    int TimeLevel) const;
 
  private:
-   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
+   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo,
                   int NVertLevels, int NTracers);
 
    AuxiliaryState(const AuxiliaryState &) = delete;
    AuxiliaryState(AuxiliaryState &&)      = delete;
 
    const HorzMesh *Mesh;
+   Halo *MeshHalo;
    static AuxiliaryState *DefaultAuxState;
    static std::map<std::string, std::unique_ptr<AuxiliaryState>> AllAuxStates;
 };

components/omega/src/ocn/OceanInit.cpp

@@ -203,13 +203,17 @@ int initOmegaModules(MPI_Comm Comm) {
       }
    }
 
-   // Update Halos and Host arrays with new state and tracer fields
+   // Update Halos and Host arrays with new state, auxiliary state, and tracer
+   // fields
 
    OceanState *DefState = OceanState::getDefault();
    I4 CurTimeLevel      = 0;
    DefState->exchangeHalo(CurTimeLevel);
    DefState->copyToHost(CurTimeLevel);
 
+   AuxiliaryState *DefAuxState = AuxiliaryState::getDefault();
+   DefAuxState->exchangeHalo();
+
    // Now update tracers - assume using same time level index
    Err = Tracers::exchangeHalo(CurTimeLevel);
    if (Err != 0) {

components/omega/src/ocn/Tendencies.cpp

@@ -203,6 +203,29 @@ int Tendencies::readTendConfig(Config *TendConfig ///< [in] Tendencies subconfig
       return DivFactorErr;
    }
 
+   I4 WindForcingErr =
+       TendConfig->get("WindForcingTendencyEnable", this->WindForcing.Enabled);
+   if (WindForcingErr != 0) {
+      LOG_CRITICAL("Tendencies: WindForcingTendencyEnable not found in "
+                   "TendConfig");
+      return WindForcingErr;
+   }
+
+   I4 BottomDragErr =
+       TendConfig->get("BottomDragTendencyEnable", this->BottomDrag.Enabled);
+   if (BottomDragErr != 0) {
+      LOG_CRITICAL("Tendencies: BottomDragTendencyEnable not found in "
+                   "TendConfig");
+      return BottomDragErr;
+   }
+
+   I4 BottomDragCoeffErr =
+       TendConfig->get("BottomDragCoeff", this->BottomDrag.Coeff);
+   if (BottomDragCoeffErr != 0 && this->BottomDrag.Enabled) {
+      LOG_CRITICAL("Tendencies: BottomDragCoeff not found in TendConfig");
+      return BottomDragCoeffErr;
+   }
+
    I4 TrHAdvErr = TendConfig->get("TracerHorzAdvTendencyEnable",
                                   this->TracerHorzAdv.Enabled);
    if (TrHAdvErr != 0) {
@@ -258,7 +281,8 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies
                        CustomTendencyType InCustomVelocityTend)
     : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh), KEGrad(Mesh),
       SSHGrad(Mesh), VelocityDiffusion(Mesh), VelocityHyperDiff(Mesh),
-      TracerHorzAdv(Mesh), TracerDiffusion(Mesh), TracerHyperDiff(Mesh),
+      WindForcing(Mesh), BottomDrag(Mesh), TracerHorzAdv(Mesh),
+      TracerDiffusion(Mesh), TracerHyperDiff(Mesh),
       CustomThicknessTend(InCustomThicknessTend),
       CustomVelocityTend(InCustomVelocityTend) {
 
@@ -338,9 +362,13 @@ void Tendencies::computeVelocityTendenciesOnly(
    OMEGA_SCOPE(LocSSHGrad, SSHGrad);
    OMEGA_SCOPE(LocVelocityDiffusion, VelocityDiffusion);
    OMEGA_SCOPE(LocVelocityHyperDiff, VelocityHyperDiff);
+   OMEGA_SCOPE(LocWindForcing, WindForcing);
+   OMEGA_SCOPE(LocBottomDrag, BottomDrag);
 
    deepCopy(LocNormalVelocityTend, 0);
 
+   const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel];
+
    // Compute potential vorticity horizontal advection
    const Array2DReal &FluxLayerThickEdge =
        AuxState->LayerThicknessAux.FluxLayerThickEdge;
@@ -398,6 +426,27 @@ void Tendencies::computeVelocityTendenciesOnly(
           });
    }
 
+   // Compute wind forcing
+   const auto &NormalStressEdge = AuxState->WindForcingAux.NormalStressEdge;
+   const auto &MeanLayerThickEdge =
+       AuxState->LayerThicknessAux.MeanLayerThickEdge;
+   if (LocWindForcing.Enabled) {
+      parallelFor(
+          {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) {
+             LocWindForcing(LocNormalVelocityTend, IEdge, KChunk,
+                            NormalStressEdge, MeanLayerThickEdge);
+          });
+   }
+
+   // Compute bottom drag
+   if (LocBottomDrag.Enabled) {
+      parallelFor(
+          {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) {
+             LocBottomDrag(LocNormalVelocityTend, IEdge, KChunk, NormalVelEdge,
+                           KECell, MeanLayerThickEdge);
+          });
+   }
+
    if (CustomVelocityTend) {
       CustomVelocityTend(LocNormalVelocityTend, State, AuxState, ThickTimeLevel,
                          VelTimeLevel, Time);

components/omega/src/ocn/Tendencies.h

@@ -41,6 +41,8 @@ class Tendencies {
    SSHGradOnEdge SSHGrad;
    VelocityDiffusionOnEdge VelocityDiffusion;
    VelocityHyperDiffOnEdge VelocityHyperDiff;
+   WindForcingOnEdge WindForcing;
+   BottomDragOnEdge BottomDrag;
    TracerHorzAdvOnCell TracerHorzAdv;
    TracerDiffOnCell TracerDiffusion;
    TracerHyperDiffOnCell TracerHyperDiff;

components/omega/src/ocn/TendencyTerms.cpp

@@ -43,6 +43,13 @@
       DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge),
       MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
 
+WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh)
+    : SaltWaterDensity(1.026e3) {}
+
+BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh)
+    : Coeff(0), CellsOnEdge(Mesh->CellsOnEdge), NVertLevels(Mesh->NVertLevels) {
+}
+
 TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),