Change shear mixing to LMD94 formulation and add ri smoothing

components/omega/configs/Default.yml

@@ -69,9 +69,10 @@ Omega:
       TriggerBVF: 0.0
     Shear:
       Enable: true
-      NuZero: 0.005
-      Alpha: 5.0
-      Exponent: 2.0
+      BaseShearValue: 0.005
+      RiCrit: 0.7
+      Exponent: 3.0
+      RiSmoothLoops: 2
   IOStreams:
     InitialVertCoord:
       UsePointerFile: false

components/omega/doc/devGuide/VerticalMixingCoeff.md

@@ -2,12 +2,12 @@
 
 # Vertical Mixing Coefficients
 
-Omega includes a `VertMix` class that provides functions that compute `VertDiff` and `VertVisc`, the
-vertical diffusivity and viscosity, where both are defined at the center of the cell and top of the layer.
+Omega includes a `VertMix` class that provides functions that compute `VertDiff`, `VertVisc`, `GradRichNum`, and `GradRichNumSmoothed`, the
+vertical diffusivity and viscosity, the gradient Richardson number, a smoothed gradient Richardson number, where all are defined at the center of the cell and top of the layer.
 Currently the values of `VertDiff` and `VertVisc` are calculated using the linear combination of three options: (1) a
 constant background mixing value, (2) a convective instability mixing value, and (3) a Richardson
-number dependent shear mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear mixing values `BruntVaisalaFreqSq` is needed, which
-is calculated by the `EOS` class.
+number dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear instability driven mixing parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear instability driven mixing values `BruntVaisalaFreqSq` is needed, which
+is calculated by the `EOS` class. `GradRichNum` is smoothed using a 1-2-1 filter to produce `GradRichNumSmoothed` which is used by the LMD94 shear instability driven mixing formulation.
 
 ## Initialization and Usage
 
@@ -37,9 +37,10 @@ VertMix:
     TriggerBVF: 0.0
   Shear:
     Enable: true
-    NuZero: 0.005
-    Alpha: 5.0
-    Exponent: 2.0
+    BaseShearValue: 0.005
+    RiCrit: 0.7
+    Exponent: 3.0
+    RiSmoothLoops: 2
 ```
 
 ## Class Structure
@@ -48,6 +49,8 @@ VertMix:
 
 - `VertDiff`: 2D array storing vertical diffusivity coefficients (m²/s)
 - `VertVisc`: 2D array storing vertical viscosity coefficients (m²/s)
+- `GradRichNum`: 2D array storing the gradient Richardson number
+- `GradRichNumSmoothed`: 2D array storing the smoothed gradient Richardson number
 
 ### Mixing Parameters
 
@@ -60,11 +63,12 @@ VertMix:
    - `ConvDiff`: Convective mixing coefficient (m²/s; Default: 1.0)
    - `ConvTriggerBVF`: Trigger threshold for convective mixing (Default: 0.0)
 
-3. Shear Mixing:
-   - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True)
-   - `ShearNuZero`: Base coefficient for Pacanowski-Philander scheme (Default: 0.005)
-   - `ShearAlpha`: Alpha parameter for P-P scheme (Default: 5.0)
-   - `ShearExponent`: Exponent parameter for P-P scheme (Default: 2.0)
+3. Shear Instability Driven Mixing:
+   - `EnableShearMix`: Flag to enable/disable shear instability driven mixing (Default: True)
+   - `BaseShearValue`: Base values of maximum viscosity and diffusivity for the LMD94 interior shear instability driven mixing formulation (Default: 0.005)
+   - `RiCrit`: Critical Richardson number for the LMD94 formulation (Default: 0.7)
+   - `ShearExponent`: Exponent parameter number for the LMD94 formulation (Default: 3.0)
+   - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 2)
 
 ## Core Functionality (Vertical Mixing Coefficient Calculation)
 
@@ -75,8 +79,3 @@ void computeVertMix(const Array2DReal &NormalVelocity,
                    const Array2DReal &TangentialVelocity,
                    const Array2DReal &BruntVaisalaFreqSq);
 ```
-
-This method combines the effects of:
-- Background mixing (constant coefficients)
-- Convective mixing (triggered by static instability)
-- Shear instability driven mixing (Pacanowski-Philander scheme; to be changed to [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) shear mixing scheme in a later development)

components/omega/doc/userGuide/VerticalMixingCoeff.md

@@ -6,9 +6,8 @@ The vertical mixing module in Omega handles the parameterization of unresolved v
 processes in the ocean. It calculates vertical diffusivity and viscosity coefficients that
 determine how properties (like momentum, heat, salt, and biogeochemical tracers) mix vertically
 in the ocean model. Currently, Omega offers three different mixing processes within the water column: (1) a constant
-background mixing value, (2) a convective instability mixing value, and (3) a Richardson number
-dependent shear instability driven mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These are linearly additive and are describe a bit
-more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future.
+background mixing value, (2) a convective instability mixing value, and (3) a Richardson number-dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 shear instability driven mixing parameterization. These are linearly additive and are describe a bit
+more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation.
 
 The user-configurable options are the following parameters in the yaml configuration file:
 
@@ -23,17 +22,17 @@ VertMix:
     TriggerBVF: 0.0      # Squared Brunt-Vaisala frequency threshold
   Shear:
     Enable: true         # Enables the shear-instability driven mixing option
-    NuZero: 1.0e-2       # Base viscosity coefficient
-    Alpha: 5.0           # Stability parameter
-    Exponent: 2          # Richardson number exponent
+    BaseShearValue: 0.005 # Base viscosity/diffusivity value
+    RiCrit: 0.7          # Critical Richardson number
+    Exponent: 3.0        # Richardson number exponent
+    RiSmoothLoops: 2     # Number of Richardson number smoothing loops
 ```
 
 ## Vertical Mixing Processes/Types
 
 ### 1. Background Mixing
 
-A constant background mixing value that represents small-scale mixing processes not explicitly resolved by the model. Typically, this is chosen to represent low values of vertical mixing
-happening in the ocean's stratified interior.
+A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior and is assumed roughly equivalent to the globally averaged interior mixing from all sources (e.g., from internal wave breaking).
 
 ### 2. Convective Mixing
 
@@ -49,22 +48,23 @@ $$
 
 This is different than some current implementations (i.e. in MPAS-Ocean and the CVMix package), where convective adjustment occurs both with unstable and neutral conditions ($N^2 \leq N^2_{crit}$). $\kappa_{conv}$ and $N^2_{crit}$ are constant parameters set in the `VertMix` section of the yaml file (`Diffusivity` and `TriggerBVF` under the `Convective` header).
 
-### 3. Shear Mixing
+### 3. Shear-Instability-Driven Mixing
 
-Mixing induced by vertical velocity shear, implemented using the Pacanowski-Philander scheme, through the gradient Richardson number (ratio of buoyancy to shear).
+Mixing induced by vertical velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear).
 
 $$
-\nu = \frac{\nu_o}{(1+\alpha Ri)^n} + \nu_b\,,
-$$
-
-$$
-\kappa = \frac{\nu}{(1+\alpha Ri)} + \kappa_b\,.
+\nu_{shear} = = \kappa_{shear} = =
+\begin{cases}
+\nu_o \quad \text{ if } Ri_g < 0\\
+\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 \leq Ri_g < Ri_{crit}\\
+0.0 \quad \text{ if } Ri_{crit} \leq Ri_g
+\end{cases}
 $$
 
-where $Ri$ is defined as:
+where $\nu_o$, $Ri_{crit}$, and $p$ are constant parameters set in the `VertMix` section of the yaml file (`BaseShearValue`, `RiCrit`, and `Exponent` under the `Shear` header). $Ri$ is defined as:
 
 $$
 Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,,
 $$
 
-where $\nu_o$, $\alpha$, $n$, $\nu_b$, $\kappa_b$ are constant parameters set in the `VertMix` section of the yaml file (`NuZero`, `Alpha`, `Exponent` under the `Shear` header and `Viscosity`, `Diffusivity` under the `Background` header). $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`.  $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero for the surface layer. In a later development, the shear mixing option will be changed to the interior shear mixing formulation in [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872).
+where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 (vertical) filter before being used to calculate the shear-instability-driven mixing. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`. $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero at the surface.