This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990),
Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90
vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum
in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization,
but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient
nu is computed from kappa_GM via thermal wind balance, and the following relation:

nu = kappa_GM * f^2 / N^2.

The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free
boundary conditions at the top and bottom.

In the current implementation, kappa_GM is assumed either (a) constant or as (b) having an
EBT structure. A third possible formulation of nu is depth-independent:

nu = f^2 * alpha

The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth.
Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode,
in which case we have 1/N^2 = h/g'.

More specifically, this commit adds a new subroutine that computes the couping coefficient
associated with GL90 via

a_cpl_gl90 = nu / h = kappa_GM * f^2 / g'

or

a_cpl_gl90 = nu / h = f^2 * alpha / h.

Further, a_cpl_gl90 is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary
layer, whose depth is set by Hbbl_gl90, and 1 otherwise. This modification is necessary to avlid fluxing
momentum into vanished layers that ride over steep topography.

Finally, a_cpl_gl90 is added to a_cpl, where the latter is the coupling coefficient associated with the
remaining vertical stresses, used in the vertical viscosity solver.

More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1),
Appendix B.

See PR mom-ocean#263

Change g' --> g^prime because doxygen and style test does not like
the apostrophe.

New diagnostics:
* au_gl90_visc: zonal viscous coupling coefficient associated with GL90, is contained in au_visc
* av_gl90_visc: meridional viscous coupling coefficient associated with GL90, is contained in av_visc
* Kv_gl90_u: GL90 vertical viscosity at u-points, is contained in Kv_u
* Kv_gl90_v: GL90 vertical viscosity at v-points, is contained in Kv_v

New diagnostics:
* du_dt_visc_gl90: zonal acceleration due to GL90 vertical viscosity, included in du_dt_visc
* dv_dt_visc_gl90: meridional acceleration due to GL90 vertical viscosity, included in dv_dt_visc

New diagnostics:
* GLwork: Kinetic Energy Source from GL90 Vertical Viscosity

The energetics of the GL90 parameterization (named "GLwork") are intentionally computed in
MOM_vert_friction, rather than in MOM_diagnostics, where the reamining kinetic energy budget
terms are computed. We have to do the computation in MOM_vert_friction to ensure sign-
definiteness when GLwork is summed in the vertical. Indeed, MOM_diagnostics does not have
access to the velocities and thicknesses used in the vertical solver, but rather uses a time-mean
barotropic transport [uv]h to compute the energy budget diagnostics. A detailed discussion and
exploration of this issue can be found in ocean-eddy-cpt#25.

As a result of not computing the energetics in MOM_diagnostics, GLwork is not exactly contained
in KE_visc. KE_visc represents the energetics of all vertical viscosity contributions, including
the GL90 vertical viscosity. We could implement a term "KE_visc_gl90" that can be 1-to-1 compared
to KE_visc; that is, KE_visc - KE_visc_gl90 would represent exactly the energetics of all viscosity
contributions EXCEPT the GL90 viscosity. If we implemented KE_visc_gl90, this term would in practice
be very similar as GLwork, but sign-definiteness is not ensured, see above.