Esmg docs (NOAA-EMC#57)

* additions for stochastic physics and ePBL perts

* cleanup of code and enhancement of ePBL perts

* Update MOM_diabatic_driver.F90

remove conflict with dev/emc

* Update MOM_diabatic_driver.F90

further resolve conflict

* Update MOM_diabatic_driver.F90

put id_sppt_wts, etc back.

* add stochy_restart writing to mom_cap

* additions for stochy restarts

* clean up debug statements

* clean up code

* fix non stochastic ePBL calculation

* re-write of stochastic code to remove CPP directives

* remove blank link in MOM_diagnostics

* clean up MOM_domains

* make stochastics optional

* correct coupled_driver/ocean_model_MOM.F90 and other cleanup

* clean up of code for MOM6 coding standards

* remove stochastics container

* revert MOM_domains.F90

* clean up of mom_ocean_model_nuopc.F90

* remove PE_here from mom_ocean_model_nuopc.F90

* remove debug statements

* stochastic physics re-write

* move stochastics to external directory

* doxygen cleanup

* add write_stoch_restart_ocn to MOM_stochastics

* add logic to remove incrments from restart if outside IAU window

* revert logic wrt increments

* add comments

* update to gfdl 20210806 (NOAA-EMC#74)

* remove white space and fix comment

* Update MOM_oda_incupd.F90

remove unused index bounds, and fix sum_h2 loop.

Co-authored-by: pjpegion <Philip.Pegion@noaa.gov>
Co-authored-by: Marshall Ward <marshall.ward@noaa.gov>

* Fussing with zotero.bib.

Getting a warning about a repeated bibliography entry for adcroft2004.
Rob thinks this is a hash failure.

* Still fussing with zotero.bib

- it was complaining about the (unused) Kasahara reference.

* Several little things, one is making sponge less verbose.

- Pointing to OBC wiki file from the lateral parameterizations doc.

- Using the MOM6 verbosity to control the time_interp verbosity.

- Making the check for negative water depths more informative.

* return a more accurate error message in MOM_stochasics

* Working on boundary layer docs.

* Done with EPBL docs?

* Undoing some patches from others

* Cleaning up too-new commits

* Adding in that SAL commit again.

* correction on type in directory name

* Added some to vertical viscisity doc.

* Cleaned up whitespace leftover from porous topomerge.

- Spacing within expressions was uneven and made multiplation look like
  POW functions. Leftover from merging NOAA-GFDL#3.
- No answer changes.

* Fix out-of-bounds k index in PPM flux

- An errant use of the porous face area led to an out-of-bounds k-index
  reported in NOAA-GFDL#19.
- Closes NOAA-EMC#19

* Adding Channel drag figure

* Take cite out of figure caption.

* Copyright year 2022

docs/conf.py

@@ -159,7 +159,7 @@ def latexPassthru(name, rawtext, text, lineno, inliner, options={}, content=[]):
 
 # General information about the project.
 project = u'MOM6'
-copyright = u'2017-2021, MOM6 developers'
+copyright = u'2017-2022, MOM6 developers'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the

src/parameterizations/vertical/MOM_CVMix_conv.F90

@@ -101,7 +101,7 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS)
 
   call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, &
                  "Diffusivity used in convective regime. Corresponding viscosity "//&
-                 "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", &
+                 "(KV_CONV) will be set to KD_CONV * PRANDTL_CONV.", &
                  units='m2/s', default=1.00)
 
   call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, &

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -911,7 +911,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
           else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K); endif
 
           ! Determine the drag contributing to the bottom boundary layer
-          ! and the Raleigh drag that acts on each layer.
+          ! and the Rayleigh drag that acts on each layer.
           if (L(K) > L(K+1)) then
             if (vol_below < bbl_thick) then
               BBL_frac = (1.0-vol_below/bbl_thick)**2

src/parameterizations/vertical/_V_diffusivity.dox

@@ -278,7 +278,7 @@ The original version concentrates buoyancy work in regions of strong stratificat
 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."
+\image html background_varying.png "Form of the vertically uniform background mixing in Danabasoglu [2012]. 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.

src/parameterizations/vertical/_V_viscosity.dox

@@ -1,13 +1,26 @@
-/*! \page Vertical_Viscosity Viscous Bottom Boundary Layer
+/*! \page Vertical_Viscosity Vertical Viscosity
+
+The vertical viscosity is composed of several components.
+
+-# The vertical diffusivity computations for the background and shear
+mixing all save contributions to the viscosity with an assumed turbulent
+Prandtl number of 1.0, though this can be changed with the PRANDTL_BKGND and
+PRANDTL_TURB parameters, respectively.
+-# If the ePBL scheme is used, it contributes to the vertical viscosity
+with a Prandtl number of PRANDTL_EPBL.
+-# If the CVMix scheme is used, it contributes to the vertical viscosity
+with a Prandtl number of PRANDTL_CONV.
+-# If the tidal mixing scheme is used, it contributes to the vertical
+viscosity with a Prandtl number of PRANDTL_TIDAL.
+
+\section set_viscous_BBL 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[
@@ -31,7 +44,7 @@ thin upwind cells helps increase the effect of viscosity and inhibits flow out o
 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).
+thickness is found using the ideas of \cite killworth1999 (hereafter KW99).
 KW99 solve the equation
 \f[
 \left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1
@@ -56,9 +69,54 @@ 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$.
+viscosity solver implies the stress \f$\mathbf{\tau}_b\f$:
+
+\f[
+    K_{bbl} = \frac{1}{2} h_{bbl} \sqrt{C_{drag}} \, u^\ast
+\f]
+
+\section section_Channel_drag Channel Drag
+
+The channel drag is an extra Rayleigh drag applied to those layers
+within the bottom boundary layer. It is called channel drag because it
+accounts for curvature of the bottom, applying the drag proportionally
+to how much of each cell is within the bottom boundary layer.
+The bottom shape is approximated as locally parabolic. The
+bottom drag is applied to each layer with a factor \f$R_k\f$, the sum
+of which is 1 over all the layers.
+
+\image html channel_drag.png "Example of layers intersecting a sloping bottom, with the blue showing the fraction of the cell over which bottom drag is applied."
+\imagelatex{channel_drag.png,Example of layers intersecting a sloping bottom\, with the blue showing the fraction of the cell over which bottom drag is applied.,\includegraphics[width=\textwidth\,height=\textheight/2\,keepaspectratio=true]}
+
+The velocity that is actually subject to the bottom drag may be
+substantially lower than the mean layer velocity, especially if only
+a small fraction of the layer's width is subject to the bottom drag.
+
+The code begins by finding the arithmetic mean of the water depths to
+find the depth at the velocity points. It then uses these to construct
+a parabolic bottom shape, valid for \f$I - \frac{1}{2}\f$ to \f$I +
+\frac{1}{2}\f$. The parabola is:
+
+\f[
+    D(x) = a x^2 + b x + D - \frac{a}{12}
+\f]
+
+For sufficiently small curvature \f$a\f$, one can drop the quadratic
+term and assume a linear function instead. We want a form that matches
+the traditional bottom drag when the bottom is flat.
+
+We defined the open fraction of each cell as \f$l(k) \equiv L(k)/L_{Tot}\f$,
+where terms of order \f$l^2\f$ will be dropped.
+
+Hallberg (personal communication) shows how they came up with the form used in the code, in which the
+\f$R_k\f$ above are set to:
+
+\f[
+   R_k = \gamma_k l_{k-1/2} \left[ \frac{12 c_{Smag} h_k}{12 c_{Smag} k_k + c_d \gamma_k (1 - \gamma_k)
+   (1 - \frac{3}{2} \gamma_k) l^2_{k-1/2} L_{Tot}} \right]
+\f]
+with the definition \f$\gamma_k \equiv (l_{k-1/2} - l_{k+1/2})/l_{k-1/2}\f$. This ensures that \f$\sum^N_{k=1}
+\gamma_k l_{k-1/2} = 1\f$ since \f$l_{1/2} = 1\f$ and \f$l_{N+1/2} = 0\f$.
 
 */