diff --git a/docs/images/background_varying.png b/docs/images/background_varying.png
new file mode 100644
index 0000000000..44a65175a0
Binary files /dev/null and b/docs/images/background_varying.png differ
diff --git a/docs/parameterizations_vertical.rst b/docs/parameterizations_vertical.rst
index 0d22787294..c9404c5088 100644
--- a/docs/parameterizations_vertical.rst
+++ b/docs/parameterizations_vertical.rst
@@ -21,9 +21,12 @@ Interior and bottom-driven mixing
---------------------------------
Kappa-shear
- MOM_kappa_shear implement the shear-driven mixing of :cite:`jackson2008`.
+ MOM_kappa_shear implements the shear-driven mixing of :cite:`jackson2008`.
+
+ :ref:`Internal_Shear_Mixing`
Internal-tide driven mixing
+
The schemes of :cite:`st_laurent2002`, :cite:`polzin2009`, and :cite:`melet2012`, are all implemented through MOM_set_diffusivity and MOM_diabatic_driver.
:ref:`Internal_Tidal_Mixing`
@@ -33,6 +36,8 @@ Vertical friction
Vertical viscosity is implemented in MOM_vert_frict and coefficient computed in MOM_set_viscosity, although contributions to viscosity from other parameterizations are calculated in those respective modules (e.g. MOM_kappa_shear, MOM_KPP, MOM_energetic_PBL).
+ :ref:`Vertical_Viscosity`
+
Vertical diffusion
------------------
diff --git a/docs/zotero.bib b/docs/zotero.bib
index 957097f217..a00fe569bd 100644
--- a/docs/zotero.bib
+++ b/docs/zotero.bib
@@ -655,6 +655,30 @@ @article{killworth1992
pages = {1379--1387}
}
+@article{killworth1999,
+ doi = {10.1175/1520-0485(1999)029<1221:atbblc>2.0.co;2},
+ year = 1999,
+ publisher = {American Meteorological Society},
+ volume = {29},
+ number = {6},
+ pages = {1221--1238},
+ author = {P. D. Killworth and N. R. Edwards},
+ title = {A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models},
+ journal = {J. Phys. Oceanography}
+}
+
+@article{zilitinkevich1996,
+ doi = {10.1007/bf02430334},
+ year = 1996,
+ publisher = {Springer Science and Business Media {LLC}},
+ volume = {81},
+ number = {3-4},
+ pages = {325--351},
+ author = {S. Zilitinkevich and D. V. Mironov},
+ title = {A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer},
+ journal = {Boundary-Layer Meteorology}
+}
+
@article{gent1995,
title = {Parameterizing {Eddy}-{Induced} {Tracer} {Transports} in {Ocean} {Circulation} {Models}},
volume = {25},
@@ -800,6 +824,18 @@ @article{jackson2008
pages = {1033--1053}
}
+@article{turner1986,
+ doi = {10.1017/s0022112086001222},
+ year = 1986,
+ publisher = {Cambridge University Press ({CUP})},
+ volume = {173},
+ pages = {431--471},
+ author = {J. S. Turner},
+ title = {Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows},
+ journal = {J. Fluid Mech.}
+}
+
+
@article{reichl2018,
title = {A simplified energetics based planetary boundary layer ({ePBL}) approach for ocean climate simulations.},
volume = {132},
@@ -1426,6 +1462,18 @@ @article{harrison2008
pages = {1894--1912}
}
+@article{danabasoglu2012,
+ doi = {10.1175/jcli-d-11-00091.1},
+ year = 2012,
+ publisher = {American Meteorological Society},
+ volume = {25},
+ number = {5},
+ pages = {1361--1389},
+ author = {G. Danabasoglu and S. C. Bates and B. P. Briegleb and S. R. Jayne and M. Jochum and W. G. Large and S. Peacock and S. G. Yeager},
+ title = {The {CCSM}4 Ocean Component},
+ journal = {J. Climate}
+}
+
@article{henyey1986,
title = {Energy and action flow through the internal wave field: {An} eikonal approach},
volume = {91},
@@ -1761,6 +1809,18 @@ @article{large1994
pages = {363--403}
}
+@article{pacanowski1981,
+ doi = {10.1175/1520-0485(1981)011<1443:povmin>2.0.co;2},
+ year = 1981,
+ publisher = {American Meteorological Society},
+ volume = {11},
+ number = {11},
+ pages = {1443--1451},
+ author = {R. C. Pacanowski and S. G. H. Philander},
+ title = {Parameterization of Vertical Mixing in Numerical Models of Tropical Oceans},
+ journal = {J. Phys. Oceanography}
+}
+
@article{van_roekel2018,
title = {The {KPP} {Boundary} {Layer} {Scheme} for the {Ocean}: {Revisiting} {Its} {Formulation} and {Benchmarking} {One}-{Dimensional} {Simulations} {Relative} to {LES}},
volume = {10},
@@ -2343,6 +2403,19 @@ @article{hallberg2000
pages = {1402--1419}
}
+@article{umlauf2005,
+ doi = {10.1016/j.csr.2004.08.004},
+ year = 2005,
+ publisher = {Elsevier {BV}},
+ volume = {25},
+ number = {7-8},
+ pages = {795--827},
+ author = {L. Umlauf and H. Burchard},
+ title = {Second-order turbulence closure models for geophysical boundary layers. A review of recent work},
+ journal = {Continental Shelf Res.}
+}
+
+
@article{easter1993,
title = {Two Modified Versions of Bott's Positive-Definite Numerical
Advection Scheme},
@@ -2545,11 +2618,60 @@ @article{hallberg2005
}
@article{bell1975,
- author = {T. H. Bell},
- year = {1975},
- title = {Lee wavews in stratified flows with simple harmonic time dependence"},
- journal = {J. Fluid Mech.},
+ doi = {10.1017/s0022112075000560},
+ year = 1975,
+ publisher = {Cambridge University Press ({CUP})},
volume = {67},
- pages = {705--722}
+ number = {4},
+ pages = {705--722},
+ author = {T. H. Bell},
+ title = {Lee waves in stratified flows with simple harmonic time dependence},
+ journal = {J. Fluid Mech.}
+}
+
+@article{nikurashin2010a,
+ doi = {10.1175/2009jpo4199.1},
+ year = 2010,
+ publisher = {American Meteorological Society},
+ volume = {40},
+ number = {5},
+ pages = {1055--1074},
+ author = {M. Nikurashin and R. Ferrari},
+ title = {Radiation and Dissipation of Internal Waves Generated by Geostrophic Motions Impinging on Small-Scale Topography: Theory},
+ journal = {J. Phys. Oceanography}
+}
+
+@article{nikurashin2010b,
+ doi = {10.1175/2010jpo4315.1},
+ year = 2010,
+ publisher = {American Meteorological Society},
+ volume = {40},
+ number = {9},
+ pages = {2025--2042},
+ author = {M. Nikurashin and R. Ferrari},
+ title = {Radiation and Dissipation of Internal Waves Generated by Geostrophic Motions Impinging on Small-Scale Topography: Application to the Southern Ocean},
+ journal = {J. Phys. Oceanography}
+}
+
+@article{miles1961,
+ title = {On the stability of heterogeneous shear flows},
+ author = {JW Miles},
+ year = {1961},
+ journal = {J. of Fluid Mech.},
+ volume = {10},
+ pages = {496--508},
+ doi = {10.1017/S0022112061000305}
+}
+
+@article{bryan1979,
+ doi = {10.1029/jc084ic05p02503},
+ year = 1979,
+ publisher = {American Geophysical Union ({AGU})},
+ volume = {84},
+ number = {C5},
+ pages = {2503},
+ author = {K. Bryan and L. J. Lewis},
+ title = {A water mass model of the World Ocean},
+ journal = {J. Geophys. Res.}
}
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index f4874252f4..0d07f0fea4 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -198,16 +198,6 @@ module MOM_set_diffusivity
contains
-!> Sets the interior vertical diffusion of scalars due to the following processes:
-!! 1. Shear-driven mixing: two options, Jackson et at. and KPP interior;
-!! 2. Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by
-!! Harrison & Hallberg, JPO 2008;
-!! 3. Double-diffusion, old method and new method via CVMix;
-!! 4. Tidal mixing: many options available, see MOM_tidal_mixing.F90;
-!! In addition, this subroutine has the option to set the interior vertical
-!! viscosity associated with processes 1,2 and 4 listed above, which is stored in
-!! visc%Kv_slow. Vertical viscosity due to shear-driven mixing is passed via
-!! visc%Kv_shear
subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, &
G, GV, US, CS, Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 99bd91d8f8..9a2680ecc1 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -115,80 +115,6 @@ module MOM_set_visc
contains
!> Calculates the thickness of the bottom boundary layer and the viscosity within that layer.
-!!
-!! A drag law is used, either linearized about an assumed bottom velocity or using the
-!! actual near-bottom velocities combined with an assumed unresolved velocity. The bottom
-!! boundary layer thickness is limited by a combination of stratification and rotation, as
-!! in the paper of Killworth and Edwards, JPO 1999. It is not necessary to calculate the
-!! thickness and viscosity every time step; instead previous values may be used.
-!!
-!! \section set_viscous_BBL Viscous Bottom Boundary Layer
-!!
-!! If set_visc_cs.bottomdraglaw is True then a bottom boundary layer viscosity and thickness
-!! are calculated so that the bottom stress is
-!! \f[
-!! \mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl}
-!! \f]
-!! If set_visc_cs.bottomdraglaw is True then the term \f$|U_{bbl}|\f$ is set equal to the
-!! value in set_visc_cs.drag_bg_vel so that \f$C_d |U_{bbl}|\f$ becomes a Rayleigh bottom drag.
-!! Otherwise \f$|U_{bbl}|\f$ is found by averaging the flow over the bottom set_visc_cs.hbbl
-!! of the model, adding the amplitude of tides set_visc_cs.tideamp and a constant
-!! set_visc_cs.drag_bg_vel. For these calculations the vertical grid at the velocity
-!! component locations is found by
-!! \f[
-!! \begin{array}{ll}
-!! \frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0
-!! \\
-!! \frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0
-!! \end{array}
-!! \f]
-!! which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward
-!! thin upwind cells helps increase the effect of viscosity and inhibits flow out of these
-!! thin cells.
-!!
-!! After diagnosing \f$|U_{bbl}|\f$ over a fixed depth an active viscous boundary layer
-!! thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter KW99).
-!! KW99 solve the equation
-!! \f[
-!! \left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1
-!! \f]
-!! for the boundary layer depth \f$h_{bbl}\f$. Here
-!! \f[
-!! h_f = \frac{C_n u_*}{f}
-!! \f]
-!! is the rotation controlled boundary layer depth in the absence of stratification.
-!! \f$u_*\f$ is the surface friction speed given by
-!! \f[
-!! u_*^2 = C_d |U_{bbl}|^2
-!! \f]
-!! and is a function of near bottom model flow.
-!! \f[
-!! h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime}
-!! \f]
-!! is the stratification controlled boundary layer depth. The non-dimensional parameters
-!! \f$C_n=0.5\f$ and \f$C_i=20\f$ are suggested by Zilitinkevich and Mironov, 1996.
-!!
-!! If a Richardson number dependent mixing scheme is being used, as indicated by
-!! set_visc_cs.rino_mix, then the boundary layer thickness is bounded to be no larger
-!! than a half of set_visc_cs.hbbl .
-!!
-!! \todo Channel drag needs to be explained
-!!
-!! A BBL viscosity is calculated so that the no-slip boundary condition in the vertical
-!! viscosity solver implies the stress \f$\mathbf{\tau}_b\f$.
-!!
-!! \subsection set_viscous_BBL_ref References
-!!
-!! \arg Killworth, P. D., and N. R. Edwards, 1999:
-!! A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models.
-!! J. Phys. Oceanogr., 29, 1221-1238,
-!! doi:10.1175/1520-0485(1999)029<1221:ATBBLC>2.0.CO;2
-!! \arg Zilitinkevich, S., Mironov, D.V., 1996:
-!! A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer.
-!! Boundary-Layer Meteorology 81, 325-351.
-!! doi:10.1007/BF02430334
-!!
subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox
index 882b73dd1b..a07663d4a1 100644
--- a/src/parameterizations/vertical/_Internal_tides.dox
+++ b/src/parameterizations/vertical/_Internal_tides.dox
@@ -4,7 +4,7 @@ Two parameterizations of vertical mixing due to internal tides are
available with the option INT_TIDE_DISSIPATION. The first is that of
\cite st_laurent2002 while the second is that of \cite polzin2009. Choose
between them with the INT_TIDE_PROFILE option. There are other relevant
-paramters which can be seen in MOM_parameter_doc.all once the main tidal
+parameters which can be seen in MOM_parameter_doc.all once the main tidal
dissipation switch is turned on.
\section section_st_laurent St Laurent et al.
@@ -69,7 +69,7 @@ case the maximum of all the contributions is used.
The vertical diffusion profile of \cite polzin2009 is a WKB-stretched
algebraic decay profile. It is based on a radiation balance equation,
-which links the dissipation profile associtated with internal breaking to
+which links the dissipation profile associated with internal breaking to
the finescale internal wave shear producing that dissipation. The vertical
profile of internal-tide driven energy dissipation can then vary in time
and space, and evolve in a changing climate (\cite melet2012). \cite melet2012
@@ -135,9 +135,9 @@ at the ocean floor, so that in both formulations:
\int_{0}^{H} \epsilon (z) dz = \frac{qE}{\rho} .
\f]
-Whereas \cite polzin2009 assumed tthat the total dissipation was locally in balance with the
+Whereas \cite polzin2009 assumed that the total dissipation was locally in balance with the
barotropic to baroclinic energy conversion rate \f$(q=1)\f$, here we use the \cite simmons2004 value
-of \f$q=1/3\f$ to retain as much consistency as passible between both parameterizations.
+of \f$q=1/3\f$ to retain as much consistency as possible between both parameterizations.
\subsection subsection_vertical_decay_scale Vertical decay-scale reformulation
@@ -212,5 +212,16 @@ of the Earth. This allows the buoyancy fluxes to tend to zero in regions
of very weak stratification, allowing a no-flux bottom boundary condition
to be satisfied.
+\section Nikurashin Lee Wave Mixing
+
+If one has the INT_TIDE_DISSIPATION flag on, there is an option to also turn on
+LEE_WAVE_DISSIPATION. The theory is presented in \cite nikurashin2010a
+while the application of it is presented in \cite nikurashin2010b. For
+the implementation in MOM6, it is required that you provide an estimate
+of the TKE loss due to the Lee waves which is then applied with either
+the St. Laurent or the Polzin vertical profile.
+
+\todo Is there a script to produce this somewhere or what???
+
*/
diff --git a/src/parameterizations/vertical/_V_diffusivity.dox b/src/parameterizations/vertical/_V_diffusivity.dox
new file mode 100644
index 0000000000..4d671fec88
--- /dev/null
+++ b/src/parameterizations/vertical/_V_diffusivity.dox
@@ -0,0 +1,284 @@
+/*! \page Internal_Shear_Mixing Internal Vertical Mixing
+
+Sets the interior vertical diffusion of scalars due to the following processes:
+
+-# Shear-driven mixing: two options, \cite jackson2008 and KPP interior;
+-# Background mixing via CVMix (Bryan-Lewis profile), the scheme described by
+ \cite harrison2008, or that in \cite danabasoglu2012.
+-# Double-diffusion, old method and new method via CVMix;
+-# Tidal mixing: many options available, see \ref Internal_Tidal_Mixing.
+
+In addition, the MOM_set_diffusivity has the option to set the interior vertical
+viscosity associated with processes 1,2 and 4 listed above, which is stored in
+visc\%Kv\_slow. Vertical viscosity due to shear-driven mixing is passed via
+visc\%Kv\_shear
+
+The resulting diffusivity, \f$K_d\f$, is the sum of all the contributions
+unless you set BBL_MIXING_AS_MAX to True, in which case the maximum of
+all the contributions is used.
+
+In addition, \f$K_d\f$ is multiplied by the term:
+
+\f[
+ \frac{N^2}{N^2 + \Omega^2}
+\f]
+
+where \f$N\f$ is the buoyancy frequency and \f$\Omega\f$ is the angular velocity
+of the Earth. This allows the buoyancy fluxes to tend to zero in regions
+of very weak stratification, allowing a no-flux bottom boundary condition
+to be satisfied.
+
+\section section_Shear Shear-driven Mixing
+
+Below the surface mixed layer, there are places in the world's oceans
+where shear mixing is known to take place. This shear-driven mixing can
+be represented in MOM6 through either CVMix or the parameterization of
+\cite jackson2008.
+
+\subsection subsection_CVMix_shear Shear-driven mixing in CVMix
+
+The community vertical mixing (CVMix) code contains options for shear
+mixing from either \cite large1994 or from \cite pacanowski1981. In MOM6,
+CVMix is included via a git submodule which loads the external CVMix
+package. The shear mixing routine in CVMix was developed to reproduce the
+observed mixing of the equatorial undercurrent in the Pacific.
+
+We first compute the gradient Richardson number \f$\mbox{Ri} = N^2 / S^2\f$,
+where \f$S\f$ is the vertical shear (\f$S = ||\bf{u}_z ||\f$) and \f$N\f$
+is the buoyancy frequency (\f$N^2 = -g \rho_z / \rho_0\f$). The
+parameterization of \cite large1994 is as follows, where the diffusivity \f$\kappa\f$
+is given by
+
+\f[
+ \kappa = \kappa_0 \left[ 1 - \min \left( 1, \frac{\mbox{Ri}}{\mbox{Ri}_c} \right) ^2 \right] ^3 ,\
+\f]
+
+with \f$\kappa_0 = 5 \times 10^{-3}\, \mbox{m}^2 \,\mbox{s}^{-1}\f$ and \f$\mbox{Ri}_c = 0.7\f$.
+
+One can instead select the \cite pacanowski1981 scheme within CVMix. Unlike
+the \cite large1994 scheme, they propose that the\ vertical shear
+viscosity \f$\nu_{\mbox{shear}}\f$ be different from the vertical shear
+diffusivity \f$\kappa_{\mbox{shear}}\f$. For gravitationally stable
+profiles (i.e., \f$N^2 > 0\f$), they chose
+
+\f[
+ \nu_{\mbox{shear}} = \frac{\nu_0}{(1 + a \mbox{Ri})^n}
+\f]
+
+\f[
+ \kappa_{\mbox{shear}} = \frac{\nu_0}{(1 + a \mbox{Ri})^{n+1}}
+\f]
+
+where \f$\nu_0\f$, \f$a\f$ and \f$n\f$ are adjustable parameters. Common settings are \f$a = 5\f$
+and \f$n = 2\f$.
+
+For both CVMix shear mixing schemes, the mixing coefficients are set to
+a large value for gravitationally unstable profiles.
+
+\subsection subsection_kappa_shear Shear-driven mixing in Jackson
+
+While the above parameterization works well enough in the equatorial
+Pacific, another place one can expect shear-mixing to matter is
+in overflows of dense water. \cite jackson2008 proposes a new shear
+parameterization with the goal of working in both the equatorial undercurrent
+and for overflows, also to have smooth transitions between unstable and
+stable regions. Their scheme looks like:
+
+\f{eqnarray}
+ \frac{\partial^2 \kappa}{\partial z^2} - \frac{\kappa}{L^2_d} &= - 2 SF(\mbox{Ri}) .
+ \label{eq:Jackson_10}
+\f}
+
+This is similar to the locally constant stratification limit of
+\cite turner1986, but with the addition of a decay length scale
+\f$L_d = \lambda L_b\f$. Here \f$L_b = Q^{1/2} / N\f$ is the buoyancy
+length scale where \f$Q\f$ is the turbulent kinetic energy (TKE) per
+unit mass, and \f$\lambda\f$ is a nondimensional constant. The function
+\f$F(\mbox{Ri})\f$ is a function of the Richardson number that remains
+to be determined. As in \cite turner1986, there must be a critical
+value of \f$\mbox{Ri}\f$ above which \f$F(\mbox{Ri}) = 0\f$.
+For better agreement with observations in a law-of-the-wall configuration,
+we modify \f$L_d\f$ to be \f$\min (\lambda L_b, L_z)\f$, where \f$L_z\f$
+is the distance to the nearest solid boundary. This can be understood by
+considering \f$L_d\f$ to be the size of the largest turbulent eddies,
+whether they are constrained by the stratification (through \f$L_b\f$)
+or through the geometry (through \f$L_z\f$).
+
+There are two length scales: the width of the low Richardson number region
+as in \cite turner1986, and the buoyancy length scale, which is the
+length scale over which the TKE is affected by the stratification (see
+\cite jackson2008 for more details). In particular, the inclusion of a
+decay length scale means that the diffusivity decays exponentially away
+from the mixing region with a length scale of \f$L_d\f$. This is important
+since turbulent eddies generated in the low \f$\mbox{Ri}\f$ layer can
+be vertically self-advected and mix nearby regions. This method yields
+a smoother diffusivity than that in \cite hallberg2000, especially in
+areas where the Richardson number is noisy.
+
+This parameterization predicts the turbulent eddy diffusivity in terms
+of the vertical profiles of velocity and density, providing that the
+TKE is known. To complete the parameterization we use a TKE \f$Q\f$
+budget such as that used in second-order turbulence closure models
+(\cite umlauf2005). We make a few additional assumptions, however,
+and use the simplified form
+
+\f{eqnarray}
+ \frac{\partial}{\partial z} \left[ (\kappa + \nu_0) \frac{\partial Q}
+ {\partial z} \right] + \kappa (S^2 - N^2) - Q(c_N N + c_S S) &= 0.
+ \label{eq:Jackson_11}
+\f}
+
+The system is therefore in balance between a vertical diffusion of
+TKE caused by both the eddy and molecular viscosity \f$(\nu_0)\f$,
+the production of TKE by shear, a sink due to stratification, and the
+dissipation. Note that we are assuming a Prandtl number of 1, although a
+parameterization for the Prandtl number could be added. We have assumed
+that the TKE reaches a quasi-steady state faster than the flow is evolving
+and faster than it can be affected by mean-flow advection so that \f$DQ/Dt =
+0\f$. Since this parameterization is meant to be used in climate models
+with low horizontal resolution and large time steps compared to the
+mixing time scales, this is a reasonable assumption. The most tenuous
+assumption is in the form of the dissipation \f$\epsilon = Q(C_N N +
+c_S S)\f$ (where \f$c_N\f$ and \f$c_S\f$ are to be determined),
+which is assumed to be dependent on the buoyancy frequency (through loss
+of energy to internal waves) and the velocity shear (through the energy
+cascade to smaller scales).
+
+We can rewrite \eqref{eq:Jackson_10} as the steady "transport" equation
+for the turbulent diffusivity (i.e., with \f$D\kappa/Dt = 0\f$),
+
+\f[
+ \frac{\partial}{\partial z} \left( \kappa \frac{\partial \kappa}{\partial z}
+ \right) + 2\kappa SF(\mbox{Ri}) - \left( \frac{\kappa}{L_d} \right)^2 -
+ \left( \frac{\partial \kappa}{\partial z} \right) ^2 = 0 .
+\f]
+
+The first term on the left can be regarded as a vertical transport of
+diffusivity, the second term as a source, and the final two as sinks.
+This equation with \eqref{eq:Jackson_11} are simple enough to solve quickly
+using an iterative technique.
+
+We also need boundary conditions for \eqref{eq:Jackson_10}
+and \eqref{eq:Jackson_11}. For the turbulent diffusivity we use
+\f$\kappa = 0\f$ since our diffusivity is numerically defined on
+layer interfaces. This ensures that there is no turbulent flux across
+boundaries. For the TKE we use boundary conditions of \f$Q = Q_0\f$ where
+\f$Q_0\f$ is a constant value of TKE, used to prevent a singularity
+in \eqref{eq:Jackson_10}, that is chosen to be small enough to not
+influence results. Note that the value of \f$\kappa\f$ calculated here
+reflects shear-driven turbulent mixing only; the total diffusivity would
+be this value plus any diffusivities due to other turbulent processes
+or a background value.
+
+Based on \cite turner1986, we choose \f$F(\mbox{Ri})\f$ of the form
+
+\f[
+ F(\mbox{Ri}) = F_0 \left( \frac{1 - \mbox{Ri} / \mbox{Ri}_c}
+ {1 + \alpha \mbox{Ri} / \mbox{Ri}_c} \right) ,
+\f]
+
+where \f$\alpha\f$ is the curvature parameter. This table shows the default
+values of the relevant parameters:
+
+
+Shear mixing parameters
+Parameter | Default value | MOM6 parameter
+ |
---|
\f$\mbox{Ri}_c\f$ | 0.25 | RINO_CRIT
+ |
\f$\nu_0\f$ | \f$1.5 \times 10^{-5}\f$ | KD_KAPPA_SHEAR_0
+ |
\f$F_0\f$ | 0.089 | SHEARMIX_RATE
+ |
\f$\alpha\f$ | -0.97 | FRI_CURVATURE
+ |
\f$\lambda\f$ | 0.82 | KAPPA_BUOY_SCALE_COEF
+ |
\f$c_N\f$ | 0.24 | TKE_N_DECAY_CONST
+ |
\f$c_S\f$ | 0.14 | TKE_SHEAR_DECAY_CONST
+ |
+
+These can all be adjusted at run time, plus some other parameters such as the maximum number of iterations
+to perform.
+
+\section section_Background Background Mixing
+
+There are three choices for the vertical background mixing: that in
+CVMix (\cite bryan1979), that in \cite harrison2008, and that in
+\cite danabasoglu2012.
+
+\subsection subsection_bryan_lewis CVMix background mixing
+
+The background vertical mixing in \cite bryan1979 is of the form:
+
+\f[
+ \kappa = C_1 + C_2 \mbox{atan} [ C_3 ( |z| - C_4 )]
+\f]
+
+where the constants are runtime parameters as shown here:
+
+
+Bryan Lewis parameters
+Parameter | Units | MOM6 parameter
+ |
---|
\f$C_1\f$ | m2 s-1 | BRYAN_LEWIS_C1
+ |
\f$C_2\f$ | m2 s-1 | BRYAN_LEWIS_C2
+ |
\f$C_3\f$ | m-1 | BRYAN_LEWIS_C3
+ |
\f$C_4\f$ | m | BRYAN_LEWIS_C4
+ |
+
+\subsection subsection_henyey Henyey IGW background mixing
+
+\cite harrison2008 choose a vertical background mixing with a latitudinal
+dependence based on \cite henyey1986. Specifically, theory predicts
+a minimum in mixing due to wave-wave interactions at the equator and
+observations support that theory. In this option, the surface background
+diffusivity is
+
+\f[
+ \kappa_s (\phi) = \max \left[ 10^{-7}, \kappa_0 \left| \frac{f}{f_{30}} \right|
+ \frac{ \cosh^{-1} (1/f) }{ \cosh^{-1} (1/f_{30})} \right] ,
+\f]
+
+where \f$f_{30}\f$ is the Coriolis frequency at \f$30^\circ\f$ latitude. The two-dimensional equation for
+the diffusivity is
+
+\f[
+ \kappa(\phi, z) = \kappa_s + \Gamma \mbox{atan} \left( \frac{H_t}{\delta_t} \right) +
+ \Gamma \mbox{atan} \left( \frac{z - H_t}{\delta_t} \right) ,
+\f]
+\f[
+ \Gamma = \frac{(\kappa_d - \kappa_s) }{\left[ 0.5 \pi + \mbox{atan} \left( \frac{H_t}{\delta_t} \right)
+ \right] },
+\f]
+
+where \f$H_t = 2500\, \mbox{m}\f$, \f$\delta_t = 222\, \mbox{m}\f$, and
+\f$\kappa_d\f$ is the deep ocean diffusivity of \f$10^{-4}\, \mbox{m}^2
+\, \mbox{s}^{-1}\f$. Note that this is the vertical structure described
+in \cite harrison2008, but that isn't what is in the code. Instead, the surface
+value is propagated down, with the assumption that the tidal mixing parameterization
+will provide the deep mixing: \ref Internal_Tidal_Mixing.
+
+There is also a "new" Henyey version, taking into account the effect of stratification on
+TKE dissipation,
+
+\f[
+ \epsilon = \epsilon_0 \frac{f}{f_0} \frac{\mbox{acosh} (N/f)}{\mbox{acosh} (N_0 / f_0)}
+\f]
+
+where \f$N_0\f$ and \f$f_0\f$ are the reference buoyancy frequency and inertial frequencies, respectively
+and \f$\epsilon_0\f$ is the reference dissipation at \f$(N_0, f_0)\f$. In the previous version, \f$N =
+N_0\f$. Additionally, the relationship between diapycnal diffusivities and stratification is included:
+
+\f[
+ \kappa = \frac{\epsilon}{N^2}
+\f]
+This approach assumes that work done against gravity is uniformly distributed throughout the water column.
+The original version concentrates buoyancy work in regions of strong stratification.
+
+\subsection subsection_danabasoglu_back Danabasoglu background mixing
+
+The shape of the \cite danabasoglu2012 background mixing has a uniform background value, with a dip
+at the equator and a bump at \f$\pm 30^{\circ}\f$ degrees latitude. The form is shown in this figure
+
+\image html background_varying.png "Form of the vertically uniform background mixing in \cite danabasoglu2012. The values are symmetric about the equator."
+\imagelatex{background_varying.png,Form of the vertically uniform background mixing in \cite danabasoglu2012. The values are symmetric about the equator.,\includegraphics[width=\textwidth\,height=\textheight/2\,keepaspectratio=true]}
+
+Some parameters of this curve are set in the input file, some are hard-coded in calculate_bkgnd_mixing.
+
+\section section_Double_Diff Double Diffusion
+
+*/
diff --git a/src/parameterizations/vertical/_V_viscosity.dox b/src/parameterizations/vertical/_V_viscosity.dox
new file mode 100644
index 0000000000..cc59e83457
--- /dev/null
+++ b/src/parameterizations/vertical/_V_viscosity.dox
@@ -0,0 +1,64 @@
+/*! \page Vertical_Viscosity Viscous Bottom Boundary Layer
+
+A drag law is used, either linearized about an assumed bottom velocity or using the
+actual near-bottom velocities combined with an assumed unresolved velocity. The bottom
+boundary layer thickness is limited by a combination of stratification and rotation, as
+in the paper of \cite killworth1999. It is not necessary to calculate the
+thickness and viscosity every time step; instead previous values may be used.
+
+\section set_viscous_BBL Viscous Bottom Boundary Layer
+
+If set_visc_CS\%bottomdraglaw is True then a bottom boundary layer viscosity and thickness
+are calculated so that the bottom stress is
+\f[
+\mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl}
+\f]
+If set_visc_CS\%bottomdraglaw is True then the term \f$|U_{bbl}|\f$ is set equal to the
+value in set_visc_CS.drag_bg_vel so that \f$C_d |U_{bbl}|\f$ becomes a Rayleigh bottom drag.
+Otherwise \f$|U_{bbl}|\f$ is found by averaging the flow over the bottom set_visc_CS\%hbbl
+of the model, adding the amplitude of tides set_visc_CS\%tideamp and a constant
+set_visc_CS\%drag_bg_vel. For these calculations the vertical grid at the velocity
+component locations is found by
+\f[
+\begin{array}{ll}
+\frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0
+\\
+\frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0
+\end{array}
+\f]
+which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward
+thin upwind cells helps increase the effect of viscosity and inhibits flow out of these
+thin cells.
+
+After diagnosing \f$|U_{bbl}|\f$ over a fixed depth an active viscous boundary layer
+thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter KW99).
+KW99 solve the equation
+\f[
+\left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1
+\f]
+for the boundary layer depth \f$h_{bbl}\f$. Here
+\f[
+h_f = \frac{C_n u_*}{f}
+\f]
+is the rotation controlled boundary layer depth in the absence of stratification.
+\f$u_*\f$ is the surface friction speed given by
+\f[
+u_*^2 = C_d |U_{bbl}|^2
+\f]
+and is a function of near bottom model flow.
+\f[
+h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime}
+\f]
+is the stratification controlled boundary layer depth. The non-dimensional parameters
+\f$C_n=0.5\f$ and \f$C_i=20\f$ are suggested by \cite zilitinkevich1996.
+
+If a Richardson number dependent mixing scheme is being used, as indicated by
+set_visc_CS\%rino_mix, then the boundary layer thickness is bounded to be no larger
+than a half of set_visc_CS\%hbbl .
+
+\todo Channel drag needs to be explained
+
+A BBL viscosity is calculated so that the no-slip boundary condition in the vertical
+viscosity solver implies the stress \f$\mathbf{\tau}_b\f$.
+
+*/