Detailed Description

Calculates horizontal viscosity and viscous stresses.

This module contains the subroutine horizontal_viscosity() that calculates the effects of horizontal viscosity, including parameterizations of the value of the viscosity itself. horizontal_viscosity() calculates the acceleration due to some combination of a biharmonic viscosity and a Laplacian viscosity. Either or both may use a coefficient that depends on the shear and strain of the flow. All metric terms are retained. The Laplacian is calculated as the divergence of a stress tensor, using the form suggested by Smagorinsky (1993). The biharmonic is calculated by twice applying the divergence of the stress tensor that is used to calculate the Laplacian, but without the dependence on thickness in the first pass. This form permits a variable viscosity, and indicates no acceleration for either resting fluid or solid body rotation.

The form of the viscous accelerations is discussed extensively in Griffies and Hallberg (2000), and the implementation here follows that discussion closely. We use the notation of Smith and McWilliams (2003) with the exception that the isotropic viscosity is $\kappa_h$.

Horizontal viscosity in MOM

In general, the horizontal stress tensor can be written as

\[ {\bf \sigma} = \begin{pmatrix} \frac{1}{2} \left( \sigma_D + \sigma_T \right) & \frac{1}{2} \sigma_S \\\\ \frac{1}{2} \sigma_S & \frac{1}{2} \left( \sigma_D - \sigma_T \right) \end{pmatrix} \]

where $\sigma_D$, $\sigma_T$ and $\sigma_S$ are stresses associated with invariant factors in the strain-rate tensor. For a Newtonian fluid, the stress tensor is usually linearly related to the strain-rate tensor. The horizontal strain-rate tensor is

\[ \dot{\bf e} = \begin{pmatrix} \frac{1}{2} \left( \dot{e}_D + \dot{e}_T \right) & \frac{1}{2} \dot{e}_S \\\\ \frac{1}{2} \dot{e}_S & \frac{1}{2} \left( \dot{e}_D - \dot{e}_T \right) \end{pmatrix} \]

where $\dot{e}_D = \partial_x u + \partial_y v$ is the horizontal divergence, $\dot{e}_T = \partial_x u - \partial_y v$ is the horizontal tension, and $\dot{e}_S = \partial_y u + \partial_x v$ is the horizontal shear strain.

The trace of the stress tensor, $tr(\bf \sigma) = \sigma_D$, is usually absorbed into the pressure and only the deviatoric stress tensor considered. From here on, we drop $\sigma_D$. The trace of the strain tensor, $tr(\bf e) = \dot{e}_D$ is non-zero for horizontally divergent flow but only enters the stress tensor through $\sigma_D$ and so we will drop $\sigma_D$ from calculations of the strain tensor in the code. Therefore the horizontal stress tensor can be considered to be

\[ {\bf \sigma} = \begin{pmatrix} \frac{1}{2} \sigma_T & \frac{1}{2} \sigma_S \\\\ \frac{1}{2} \sigma_S & - \frac{1}{2} \sigma_T \end{pmatrix} .\]

The stresses above are linearly related to the strain through a viscosity coefficient, $\kappa_h$:

\begin{eqnarray*} \sigma_T & = & 2 \kappa_h \dot{e}_T \\\\ \sigma_S & = & 2 \kappa_h \dot{e}_S . \end{eqnarray*}

The viscosity $\kappa_h$ may either be a constant or variable. For example, $\kappa_h$ may vary with the shear, as proposed by Smagorinsky (1993).

The accelerations resulting form the divergence of the stress tensor are

\begin{eqnarray*} \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) & = & \partial_x \left( \frac{1}{2} \sigma_T \right) + \partial_y \left( \frac{1}{2} \sigma_S \right) \\\\ & = & \partial_x \left( \kappa_h \dot{e}_T \right) + \partial_y \left( \kappa_h \dot{e}_S \right) \\\\ \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) & = & \partial_x \left( \frac{1}{2} \sigma_S \right) + \partial_y \left( \frac{1}{2} \sigma_T \right) \\\\ & = & \partial_x \left( \kappa_h \dot{e}_S \right) + \partial_y \left( - \kappa_h \dot{e}_T \right) . \end{eqnarray*}

The form of the Laplacian viscosity in general coordinates is:

\begin{eqnarray*} \hat{\bf x} \cdot \left( \nabla \cdot \sigma \right) & = & \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_T \right) + \partial_y \left( \kappa_h h \dot{e}_S \right) \right] \\\\ \hat{\bf y} \cdot \left( \nabla \cdot \sigma \right) & = & \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_S \right) - \partial_y \left( \kappa_h h \dot{e}_T \right) \right] . \end{eqnarray*}

Laplacian viscosity coefficient

The horizontal viscosity coefficient, $\kappa_h$, can have multiple components. The isotropic components are:

A uniform background component, $\kappa_{bg}$.
A constant but spatially variable 2D map, $\kappa_{2d}(x,y)$.
A ''MICOM'' viscosity, $U_\nu \Delta(x,y)$, which uses a constant velocity scale, $U_\nu$ and a measure of the grid-spacing $\Delta(x,y)^2 = \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}$.
A function of latitude, $\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n$.
A dynamic Smagorinsky viscosity, $\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}$.
A dynamic Leith viscosity, $\kappa_{Lth}(x,y,t) = C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}$.

A maximum stable viscosity, $\kappa_{max}(x,y)$ is calculated based on the grid-spacing and time-step and used to clip calculated viscosities.

The static components of $\kappa_h$ are first combined as follows:

\[ \kappa_{static} = \min \left[ \max\left( \kappa_{bg}, U_\nu \Delta(x,y), \kappa_{2d}(x,y), \kappa_\phi(x,y) \right) , \kappa_{max}(x,y) \right] \]

and stored in the module control structure as variables Kh_bg_xx and Kh_bg_xy for the tension (h-points) and shear (q-points) components respectively.

The full viscosity includes the dynamic components as follows:

\[ \kappa_h(x,y,t) = r(\Delta,L_d) \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right) \]

where $r(\Delta,L_d)$ is a resolution function.

The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each other.

Viscous boundary conditions

Free slip boundary conditions have been coded, although no slip boundary conditions can be used with the Laplacian viscosity based on the 2D land-sea mask. For a western boundary, for example, the boundary conditions with the biharmonic operator would be written as:

\[ \partial_x v = 0 ; \partial_x^3 v = 0 ; u = 0 ; \partial_x^2 u = 0 , \]

while for a Laplacian operator, they are simply

\[ \partial_x v = 0 ; u = 0 . \]

These boundary conditions are largely dictated by the use of an Arakawa C-grid and by the varying layer thickness.

Anisotropic viscosity

Large et al., 2001, proposed enhancing viscosity in a particular direction and the approach was generalized in Smith and McWilliams, 2003. We use the second form of their two coefficient anisotropic viscosity (section 4.3). We also replace their $A^\prime$ and $D$ such that $2A^\prime = 2 \kappa_h + D$ and $\kappa_a = D$ so that $\kappa_h$ can be considered the isotropic viscosity and $\kappa_a=D$ can be consider the anisotropic viscosity. The direction of anisotropy is defined by a unit vector $\hat{\bf n}=(n_1,n_2)$.

The contributions to the stress tensor are

\[ \begin{pmatrix} \sigma_T \\\\ \sigma_S \end{pmatrix} = \left[ \begin{pmatrix} 2 \kappa_h + \kappa_a & 0 \\\\ 0 & 2 \kappa_h \end{pmatrix} + 2 \kappa_a n_1 n_2 \begin{pmatrix} - 2 n_1 n_2 & n_1^2 - n_2^2 \\\\ n_1^2 - n_2^2 & 2 n_1 n_2 \end{pmatrix} \right] \begin{pmatrix} \dot{e}_T \\\\ \dot{e}_S \end{pmatrix} \]

Dissipation of kinetic energy requires $\kappa_h \geq 0$ and $2 \kappa_h + \kappa_a \geq 0$. Note that when anisotropy is aligned with the x-direction, $n_1 = \pm 1$, then $n_2 = 0$ and the cross terms vanish. The accelerations in this aligned limit with constant coefficients become

\begin{eqnarray*} \hat{\bf x} \cdot \nabla \cdot {\bf \sigma} & = & \partial_x \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right) + \partial_y \left( \kappa_h \dot{e}_S \right) \\\\ & = & \left( \kappa_h + \kappa_a \right) \partial_{xx} u + \kappa_h \partial_{yy} u - \frac{1}{2} \kappa_a \partial_x \left( \partial_x u + \partial_y v \right) \\\\ \hat{\bf y} \cdot \nabla \cdot {\bf \sigma} & = & \partial_x \left( \kappa_h \dot{e}_S \right) - \partial_y \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right) \\\\ & = & \kappa_h \partial_{xx} v + \left( \kappa_h + \kappa_a \right) \partial_{yy} v - \frac{1}{2} \kappa_a \partial_y \left( \partial_x u + \partial_y v \right) \end{eqnarray*}

which has contributions akin to a negative divergence damping (a divergence enhancement?) but which is weaker than the enhanced tension terms by half.

Discretization

The horizontal tension, $\dot{e}_T$, is stored in variable sh_xx and discretized as

\[ \dot{e}_T = \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} u \right) - \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} v \right) . \]

The horizontal divergent strain, $\dot{e}_D$, is stored in variable div_xx and discretized as

\[ \dot{e}_D = \frac{1}{h A} \left( \delta_i \left( \overline{h}^i \Delta y \, u \right) + \delta_j \left( \overline{h}^j\Delta x \, v \right) \right) . \]

Note that for expediency this is the exact discretization used in the continuity equation.

The horizontal shear strain, $\dot{e}_S$, is stored in variable sh_xy and discretized as

\[ \dot{e}_S = v_x + u_y \]

where

\begin{align*} v_x &= \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} v \right) \\\\ u_y &= \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} u \right) \end{align*}

which are calculated separately so that no-slip or free-slip boundary conditions can be applied to $v_x$ and $u_y$ where appropriate.

The tendency for the x-component of the divergence of stress is stored in variable diffu and discretized as

\[ \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) = \frac{1}{A \overline{h}^i} \left( \frac{1}{\Delta y} \delta_i \left( h \Delta y^2 \kappa_h \dot{e}_T \right) + \frac{1}{\Delta x} \delta_j \left( \tilde{h}^{ij} \Delta x^2 \kappa_h \dot{e}_S \right) \right) . \]

The tendency for the y-component of the divergence of stress is stored in variable diffv and discretized as

\[ \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) = \frac{1}{A \overline{h}^j} \left( \frac{1}{\Delta y} \delta_i \left( \tilde{h}^{ij} \Delta y^2 A_M \dot{e}_S \right) - \frac{1}{\Delta x} \delta_j \left( h \Delta x^2 A_M \dot{e}_T \right) \right) . \]

References

Griffies, S.M., and Hallberg, R.W., 2000: Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Monthly Weather Review, 128(8), 2935-2946. https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2

Large, W.G., Danabasoglu, G., McWilliams, J.C., Gent, P.R. and Bryan, F.O., 2001: Equatorial circulation of a global ocean climate model with anisotropic horizontal viscosity. Journal of Physical Oceanography, 31(2), pp.518-536. https://doi.org/10.1175/1520-0485(2001)031%3C0518:ECOAGO%3E2.0.CO;2

Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear viscosities. Large eddy simulation of complex engineering and geophysical flows, 1, 69-106.

Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for ocean models. Ocean Modelling, 5(2), 129-156. https://doi.org/10.1016/S1463-5003(02)00016-1

Data Types
type	hor_visc_cs
	Control structure for horizontal viscosity. More...

Functions/Subroutines

subroutine, public	horizontal_viscosity (u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS, OBC, BT)
	Calculates the acceleration due to the horizontal viscosity. More...

subroutine, public	hor_visc_init (Time, G, US, param_file, diag, CS, MEKE)
	Allocates space for and calculates static variables used by horizontal_viscosity(). hor_visc_init calculates and stores the values of a number of metric functions that are used in horizontal_viscosity(). More...

subroutine	align_aniso_tensor_to_grid (CS, n1, n2)
	Calculates factors in the anisotropic orientation tensor to be align with the grid. With n1=1 and n2=0, this recovers the approach of Large et al, 2001. More...

subroutine	smooth_gme (CS, G, GME_flux_h, GME_flux_q)
	Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any horizontal two-grid-point noise. More...

subroutine, public	hor_visc_end (CS)
	Deallocates any variables allocated in hor_visc_init. More...

Function/Subroutine Documentation

◆ align_aniso_tensor_to_grid()

subroutine mom_hor_visc::align_aniso_tensor_to_grid	(	type(hor_visc_cs), pointer	CS,
		real, intent(in)	n1,
		real, intent(in)	n2
	)

private

Calculates factors in the anisotropic orientation tensor to be align with the grid. With n1=1 and n2=0, this recovers the approach of Large et al, 2001.

Parameters

	cs	Control structure for horizontal viscosity
[in]	n1	i-component of direction vector [nondim]
[in]	n2	j-component of direction vector [nondim]

Definition at line 2123 of file MOM_hor_visc.F90.

   type(hor_visc_CS), pointer :: CS !< Control structure for horizontal viscosity
   real,              intent(in) :: n1 !< i-component of direction vector [nondim]
   real,              intent(in) :: n2 !< j-component of direction vector [nondim]
   ! Local variables
   real :: recip_n2_norm
  
   ! For normalizing n=(n1,n2) in case arguments are not a unit vector
   recip_n2_norm = n1**2 + n2**2
   if (recip_n2_norm > 0.) recip_n2_norm = 1./recip_n2_norm
  
   cs%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
   cs%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
   cs%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
   cs%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
  

Referenced by hor_visc_init().

Here is the caller graph for this function:

◆ hor_visc_end()

subroutine, public mom_hor_visc::hor_visc_end ( type(hor_visc_cs), pointer CS )

Deallocates any variables allocated in hor_visc_init.

Parameters

cs	The control structure returned by a previous call to hor_visc_init.

Definition at line 2215 of file MOM_hor_visc.F90.

   type(hor_visc_CS), pointer :: CS !< The control structure returned by a
                                    !! previous call to hor_visc_init.
  
   if (cs%Laplacian .or. cs%biharmonic) then
     dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)
     dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)
     dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)
   endif
  
   if (cs%Laplacian) then
     dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)
     if (cs%bound_Kh) then
       dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
     endif
     if (cs%Smagorinsky_Kh) then
       dealloc_(cs%Laplac2_const_xx) ; dealloc_(cs%Laplac2_const_xy)
     endif
     if (cs%Leith_Kh) then
       dealloc_(cs%Laplac3_const_xx) ; dealloc_(cs%Laplac3_const_xy)
     endif
   endif
  
   if (cs%biharmonic) then
     dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)
     dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)
     dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)
     if (cs%bound_Ah) then
       dealloc_(cs%Ah_Max_xx) ; dealloc_(cs%Ah_Max_xy)
     endif
     if (cs%Smagorinsky_Ah) then
       dealloc_(cs%Biharm5_const_xx) ; dealloc_(cs%Biharm5_const_xy)
       ! if (CS%bound_Coriolis) then
       !   DEALLOC_(CS%Biharm5_const2_xx) ; DEALLOC_(CS%Biharm5_const2_xy)
       ! endif
     endif
     if (cs%Leith_Ah) then
       dealloc_(cs%Biharm_const_xx) ; dealloc_(cs%Biharm_const_xy)
     endif
   endif
   if (cs%anisotropic) then
     dealloc_(cs%n1n2_h)
     dealloc_(cs%n1n2_q)
     dealloc_(cs%n1n1_m_n2n2_h)
     dealloc_(cs%n1n1_m_n2n2_q)
   endif
   deallocate(cs)
  

◆ hor_visc_init()

subroutine, public mom_hor_visc::hor_visc_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(inout)	G,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(hor_visc_cs), pointer	CS,
		type(meke_type), pointer	MEKE
	)

Allocates space for and calculates static variables used by horizontal_viscosity(). hor_visc_init calculates and stores the values of a number of metric functions that are used in horizontal_viscosity().

Parameters

[in]	time	Current model time.
[in,out]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
[in]	param_file	A structure to parse for run-time parameters.
[in,out]	diag	Structure to regulate diagnostic output.
	cs	Pointer to the control structure for this module
	meke	MEKE data

Definition at line 1349 of file MOM_hor_visc.F90.

   type(time_type),         intent(in)    :: Time !< Current model time.
   type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
   type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
   type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                  !! parameters.
   type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.
   type(hor_visc_CS), pointer             :: CS   !< Pointer to the control structure for this module
   type(MEKE_type), pointer               :: MEKE !< MEKE data
   ! Local variables
   real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v
   real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v
                 ! u0v is the Laplacian sensitivities to the v velocities
                 ! at u points [L-2 ~> m-2], with u0u, v0u, and v0v defined similarly.
   real :: grid_sp_h2       ! Harmonic mean of the squares of the grid [L2 ~> m2]
   real :: grid_sp_h3       ! Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
   real :: grid_sp_q2       ! spacings at h and q points [L2 ~> m2]
   real :: grid_sp_q3       ! spacings at h and q points^(3/2) [L3 ~> m3]
   real :: Kh_Limit         ! A coefficient [T-1 ~> s-1] used, along with the
                            ! grid spacing, to limit Laplacian viscosity.
   real :: fmax             ! maximum absolute value of f at the four
                            ! vorticity points around a thickness point [T-1 ~> s-1]
   real :: BoundCorConst    ! A constant used when using viscosity to bound the Coriolis accelerations
                            ! [T2 L-2 ~> s2 m-2]
   real :: Ah_Limit         ! coefficient [T-1 ~> s-1] used, along with the
                            ! grid spacing, to limit biharmonic viscosity
   real :: Kh               ! Lapacian horizontal viscosity [L2 T-1 ~> m2 s-1]
   real :: Ah               ! biharmonic horizontal viscosity [L4 T-1 ~> m4 s-1]
   real :: Kh_vel_scale     ! this speed [L T-1 ~> m s-1] times grid spacing gives Lap visc
   real :: Ah_vel_scale     ! this speed [L T-1 ~> m s-1] times grid spacing cubed gives bih visc
   real :: Ah_time_scale    ! damping time-scale for biharmonic visc [T ~> s]
   real :: Smag_Lap_const   ! nondimensional Laplacian Smagorinsky constant
   real :: Smag_bi_const    ! nondimensional biharmonic Smagorinsky constant
   real :: Leith_Lap_const  ! nondimensional Laplacian Leith constant
   real :: Leith_bi_const   ! nondimensional biharmonic Leith constant
   real :: dt               ! The dynamics time step [T ~> s]
   real :: Idt              ! The inverse of dt [T-1 ~> s-1]
   real :: denom            ! work variable; the denominator of a fraction
   real :: maxvel           ! largest permitted velocity components [m s-1]
   real :: bound_Cor_vel    ! grid-scale velocity variations at which value
                            ! the quadratically varying biharmonic viscosity
                            ! balances Coriolis acceleration [L T-1 ~> m s-1]
   real :: Kh_sin_lat       ! Amplitude of latitudinally dependent viscosity [L2 T-1 ~> m2 s-1]
   real :: Kh_pwr_of_sine   ! Power used to raise sin(lat) when using Kh_sin_lat
   logical :: bound_Cor_def ! parameter setting of BOUND_CORIOLIS
   logical :: get_all       ! If true, read and log all parameters, regardless of
                            ! whether they are used, to enable spell-checking of
                            ! valid parameters.
   logical :: split         ! If true, use the split time stepping scheme.
                            ! If false and USE_GME = True, issue a FATAL error.
   logical :: use_MEKE      ! If true, use the MEKE module for calculating eddy kinetic energy.
                            ! If false and USE_GME = True, issue a FATAL error.
   logical :: default_2018_answers
  
   character(len=64) :: inputdir, filename
   real    :: deg2rad       ! Converts degrees to radians
   real    :: slat_fn       ! sin(lat)**Kh_pwr_of_sine
   real    :: aniso_grid_dir(2) ! Vector (n1,n2) for anisotropic direction
   integer :: aniso_mode    ! Selects the mode for setting the anisotropic direction
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
   integer :: i, j
  
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mdl = "MOM_hor_visc"  ! module name
  
   is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec ; nz = g%ke
   isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
   isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed
   isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  
   if (associated(cs)) then
     call mom_error(warning, "hor_visc_init called with an associated "// &
                             "control structure.")
     return
   endif
   allocate(cs)
  
   cs%diag => diag
  
   ! Read parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
  
   !   It is not clear whether these initialization lines are needed for the
   ! cases where the corresponding parameters are not read.
   cs%bound_Kh = .false. ; cs%better_bound_Kh = .false. ; cs%Smagorinsky_Kh = .false. ; cs%Leith_Kh = .false.
   cs%bound_Ah = .false. ; cs%better_bound_Ah = .false. ; cs%Smagorinsky_Ah = .false. ; cs%Leith_Ah = .false.
   cs%use_QG_Leith_visc = .false.
   cs%bound_Coriolis = .false.
   cs%Modified_Leith = .false.
   cs%anisotropic = .false.
   cs%dynamic_aniso = .false.
  
   kh = 0.0 ; ah = 0.0
  
   !   If GET_ALL_PARAMS is true, all parameters are read in all cases to enable
   ! parameter spelling checks.
   call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.)
  
   call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                  "This sets the default value for the various _2018_ANSWERS parameters.", &
                  default=.true.)
   call get_param(param_file, mdl, "HOR_VISC_2018_ANSWERS", cs%answers_2018, &
                  "If true, use the order of arithmetic and expressions that recover the "//&
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
  
   call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
  
   call get_param(param_file, mdl, "LAPLACIAN", cs%Laplacian, &
                  "If true, use a Laplacian horizontal viscosity.", &
                  default=.false.)
   if (cs%Laplacian .or. get_all) then
     call get_param(param_file, mdl, "KH", kh,                      &
                  "The background Laplacian horizontal viscosity.", &
                  units = "m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
     call get_param(param_file, mdl, "KH_BG_MIN", cs%Kh_bg_min, &
                  "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
                  units = "m2 s-1",  default=0.0, scale=us%m_to_L**2*us%T_to_s)
     call get_param(param_file, mdl, "KH_VEL_SCALE", kh_vel_scale, &
                  "The velocity scale which is multiplied by the grid "//&
                  "spacing to calculate the Laplacian viscosity. "//&
                  "The final viscosity is the largest of this scaled "//&
                  "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
                  units="m s-1", default=0.0, scale=us%m_s_to_L_T)
     call get_param(param_file, mdl, "KH_SIN_LAT", kh_sin_lat, &
                  "The amplitude of a latitudinally-dependent background "//&
                  "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &
                  units = "m2 s-1",  default=0.0, scale=us%m_to_L**2*us%T_to_s)
     if (kh_sin_lat>0. .or. get_all) &
       call get_param(param_file, mdl, "KH_PWR_OF_SINE", kh_pwr_of_sine, &
                  "The power used to raise SIN(LAT) when using a latitudinally "//&
                  "dependent background viscosity.", &
                  units = "nondim",  default=4.0)
  
     call get_param(param_file, mdl, "SMAGORINSKY_KH", cs%Smagorinsky_Kh, &
                  "If true, use a Smagorinsky nonlinear eddy viscosity.", &
                  default=.false.)
     if (cs%Smagorinsky_Kh .or. get_all) &
       call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_lap_const, &
                  "The nondimensional Laplacian Smagorinsky constant, "//&
                  "often 0.15.", units="nondim", default=0.0, &
                   fail_if_missing = cs%Smagorinsky_Kh)
  
     call get_param(param_file, mdl, "LEITH_KH", cs%Leith_Kh, &
                  "If true, use a Leith nonlinear eddy viscosity.", &
                  default=.false.)
  
     call get_param(param_file, mdl, "MODIFIED_LEITH", cs%Modified_Leith, &
                  "If true, add a term to Leith viscosity which is "//&
                  "proportional to the gradient of divergence.", &
                  default=.false.)
     call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", cs%res_scale_MEKE, &
                  "If true, the viscosity contribution from MEKE is scaled by "//&
                  "the resolution function.", default=.false.)
  
     if (cs%Leith_Kh .or. get_all) then
       call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
                  "The nondimensional Laplacian Leith constant, "//&
                  "often set to 1.0", units="nondim", default=0.0, &
                   fail_if_missing = cs%Leith_Kh)
       call get_param(param_file, mdl, "USE_QG_LEITH_VISC", cs%use_QG_Leith_visc, &
                  "If true, use QG Leith nonlinear eddy viscosity.", &
                  default=.false.)
       if (cs%use_QG_Leith_visc .and. .not. cs%Leith_Kh) call mom_error(fatal, &
                  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
                  "LEITH_KH must be True when USE_QG_LEITH_VISC=True.")
     endif
     if (cs%Leith_Kh .or. cs%Leith_Ah .or. get_all) then
       call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_Leith, &
                  "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
                  default=cs%Leith_Kh)
       call get_param(param_file, mdl, "MODIFIED_LEITH", cs%modified_Leith, &
                  "If true, add a term to Leith viscosity which is \n"//&
                  "proportional to the gradient of divergence.", &
                  default=.false.)
     endif
     call get_param(param_file, mdl, "BOUND_KH", cs%bound_Kh, &
                  "If true, the Laplacian coefficient is locally limited "//&
                  "to be stable.", default=.true.)
     call get_param(param_file, mdl, "BETTER_BOUND_KH", cs%better_bound_Kh, &
                  "If true, the Laplacian coefficient is locally limited "//&
                  "to be stable with a better bounding than just BOUND_KH.", &
                  default=cs%bound_Kh)
     call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", cs%anisotropic, &
                  "If true, allow anistropic viscosity in the Laplacian "//&
                  "horizontal viscosity.", default=.false.)
     call get_param(param_file, mdl, "ADD_LES_VISCOSITY", cs%add_LES_viscosity, &
                  "If true, adds the viscosity from Smagorinsky and Leith to the "//&
                  "background viscosity instead of taking the maximum.", default=.false.)
   endif
   if (cs%anisotropic .or. get_all) then
     call get_param(param_file, mdl, "KH_ANISO", cs%Kh_aniso, &
                  "The background Laplacian anisotropic horizontal viscosity.", &
                  units = "m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
     call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, &
                  "Selects the mode for setting the direction of anistropy.\n"//&
                  "\t 0 - Points along the grid i-direction.\n"//&
                  "\t 1 - Points towards East.\n"//&
                  "\t 2 - Points along the flow direction, U/|U|.", &
                  default=0)
     select case (aniso_mode)
       case (0)
         call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
                  "The vector pointing in the direction of anistropy for "//&
                  "horizont viscosity. n1,n2 are the i,j components relative "//&
                  "to the grid.", units = "nondim", fail_if_missing=.true.)
       case (1)
         call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
                  "The vector pointing in the direction of anistropy for "//&
                  "horizont viscosity. n1,n2 are the i,j components relative "//&
                  "to the spherical coordinates.", units = "nondim", fail_if_missing=.true.)
     end select
   endif
  
   call get_param(param_file, mdl, "BIHARMONIC", cs%biharmonic, &
                  "If true, use a biharmonic horizontal viscosity. "//&
                  "BIHARMONIC may be used with LAPLACIAN.", &
                  default=.true.)
   if (cs%biharmonic .or. get_all) then
     call get_param(param_file, mdl, "AH", ah, &
                  "The background biharmonic horizontal viscosity.", &
                  units = "m4 s-1", default=0.0, scale=us%m_to_L**4*us%T_to_s)
     call get_param(param_file, mdl, "AH_VEL_SCALE", ah_vel_scale, &
                  "The velocity scale which is multiplied by the cube of "//&
                  "the grid spacing to calculate the biharmonic viscosity. "//&
                  "The final viscosity is the largest of this scaled "//&
                  "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
                  units="m s-1", default=0.0, scale=us%m_s_to_L_T)
     call get_param(param_file, mdl, "AH_TIME_SCALE", ah_time_scale, &
                  "A time scale whose inverse is multiplied by the fourth "//&
                  "power of the grid spacing to calculate biharmonic viscosity. "//&
                  "The final viscosity is the largest of all viscosity "//&
                  "formulations in use. 0.0 means that it's not used.", &
                  units="s", default=0.0, scale=us%s_to_T)
     call get_param(param_file, mdl, "SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
                  "If true, use a biharmonic Smagorinsky nonlinear eddy "//&
                  "viscosity.", default=.false.)
     call get_param(param_file, mdl, "LEITH_AH", cs%Leith_Ah, &
                  "If true, use a biharmonic Leith nonlinear eddy "//&
                  "viscosity.", default=.false.)
  
     call get_param(param_file, mdl, "BOUND_AH", cs%bound_Ah, &
                  "If true, the biharmonic coefficient is locally limited "//&
                  "to be stable.", default=.true.)
     call get_param(param_file, mdl, "BETTER_BOUND_AH", cs%better_bound_Ah, &
                  "If true, the biharmonic coefficient is locally limited "//&
                  "to be stable with a better bounding than just BOUND_AH.", &
                  default=cs%bound_Ah)
  
     if (cs%Smagorinsky_Ah .or. get_all) then
       call get_param(param_file, mdl, "SMAG_BI_CONST",smag_bi_const, &
                  "The nondimensional biharmonic Smagorinsky constant, "//&
                  "typically 0.015 - 0.06.", units="nondim", default=0.0, &
                  fail_if_missing = cs%Smagorinsky_Ah)
  
       call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_cor_def, default=.false.)
       call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &
                  "If true use a viscosity that increases with the square "//&
                  "of the velocity shears, so that the resulting viscous "//&
                  "drag is of comparable magnitude to the Coriolis terms "//&
                  "when the velocity differences between adjacent grid "//&
                  "points is 0.5*BOUND_CORIOLIS_VEL.  The default is the "//&
                  "value of BOUND_CORIOLIS (or false).", default=bound_cor_def)
       if (cs%bound_Coriolis .or. get_all) then
         call get_param(param_file, mdl, "MAXVEL", maxvel, default=3.0e8)
         bound_cor_vel = maxvel
         call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_cor_vel, &
                  "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&
                  "the biharmonic drag to have comparable magnitude to the "//&
                  "Coriolis acceleration.  The default is set by MAXVEL.", &
                  units="m s-1", default=maxvel, scale=us%m_s_to_L_T)
       endif
     endif
  
     if (cs%Leith_Ah .or. get_all) &
       call get_param(param_file, mdl, "LEITH_BI_CONST", leith_bi_const, &
                  "The nondimensional biharmonic Leith constant, "//&
                  "typical values are thus far undetermined.", units="nondim", default=0.0, &
                  fail_if_missing = cs%Leith_Ah)
  
   endif
  
   call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", cs%use_land_mask, &
                  "If true, use Use the land mask for the computation of thicknesses "//&
                  "at velocity locations. This eliminates the dependence on arbitrary "//&
                  "values over land or outside of the domain. Default is False in order to "//&
                  "maintain answers with legacy experiments but should be changed to True "//&
                  "for new experiments.", default=.false.)
  
   if (cs%better_bound_Ah .or. cs%better_bound_Kh .or. get_all) &
     call get_param(param_file, mdl, "HORVISC_BOUND_COEF", cs%bound_coef, &
                  "The nondimensional coefficient of the ratio of the "//&
                  "viscosity bounds to the theoretical maximum for "//&
                  "stability without considering other terms.", units="nondim", &
                  default=0.8)
  
   call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
                  "If true, no slip boundary conditions are used; otherwise "//&
                  "free slip boundary conditions are assumed. The "//&
                  "implementation of the free slip BCs on a C-grid is much "//&
                  "cleaner than the no slip BCs. The use of free slip BCs "//&
                  "is strongly encouraged, and no slip BCs are not used with "//&
                  "the biharmonic viscosity.", default=.false.)
  
   call get_param(param_file, mdl, "USE_KH_BG_2D", cs%use_Kh_bg_2d, &
                  "If true, read a file containing 2-d background harmonic "//&
                  "viscosities. The final viscosity is the maximum of the other "//&
                  "terms and this background value.", default=.false.)
  
   call get_param(param_file, mdl, "USE_GME", cs%use_GME, &
                  "If true, use the GM+E backscatter scheme in association \n"//&
                  "with the Gent and McWilliams parameterization.", default=.false.)
  
   if (cs%use_GME) then
     call get_param(param_file, mdl, "SPLIT", split, &
                  "Use the split time stepping if true.", default=.true., &
                   do_not_log=.true.)
     if (.not. split) call mom_error(fatal,"ERROR: Currently, USE_GME = True "// &
                                            "cannot be used with SPLIT=False.")
  
     call get_param(param_file, mdl, "USE_MEKE", use_meke, &
                  "If true, turns on the MEKE scheme which calculates\n"// &
                  "a sub-grid mesoscale eddy kinetic energy budget.", &
                  default=.false.)
  
     if (.not. use_meke) call mom_error(fatal,"ERROR: Currently, USE_GME = True "// &
                                            "cannot be used with USE_MEKE=False.")
  
     call get_param(param_file, mdl, "GME_H0", cs%GME_h0, &
                  "The strength of GME tapers quadratically to zero when the bathymetric "//&
                  "depth is shallower than GME_H0.", units="m", scale=us%m_to_Z, &
                  default=1000.0)
  
     call get_param(param_file, mdl, "GME_EFFICIENCY", cs%GME_efficiency, &
                  "The nondimensional prefactor multiplying the GME coefficient.", &
                  units="nondim", default=1.0)
  
     call get_param(param_file, mdl, "GME_LIMITER", cs%GME_limiter, &
                  "The absolute maximum value the GME coefficient is allowed to take.", &
                  units="m2 s-1", scale=us%m_to_L**2*us%T_to_s, default=1.0e7)
  
   endif
  
   if (cs%bound_Kh .or. cs%bound_Ah .or. cs%better_bound_Kh .or. cs%better_bound_Ah) &
     call get_param(param_file, mdl, "DT", dt, &
                  "The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T, &
                  fail_if_missing=.true.)
  
   if (cs%no_slip .and. cs%biharmonic) &
     call mom_error(fatal,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
                          "at the same time in MOM.")
  
   if (.not.(cs%Laplacian .or. cs%biharmonic)) then
     ! Only issue inviscid warning if not in single column mode (usually 2x2 domain)
     if ( max(g%domain%niglobal, g%domain%njglobal)>2 ) call mom_error(warning, &
       "hor_visc_init:  It is usually a very bad idea not to use either "//&
       "LAPLACIAN or BIHARMONIC viscosity.")
     return ! We are not using either Laplacian or Bi-harmonic lateral viscosity
   endif
  
   deg2rad = atan(1.0) / 45.
  
   alloc_(cs%dx2h(isd:ied,jsd:jed))        ; cs%dx2h(:,:)    = 0.0
   alloc_(cs%dy2h(isd:ied,jsd:jed))        ; cs%dy2h(:,:)    = 0.0
   alloc_(cs%dx2q(isdb:iedb,jsdb:jedb))    ; cs%dx2q(:,:)    = 0.0
   alloc_(cs%dy2q(isdb:iedb,jsdb:jedb))    ; cs%dy2q(:,:)    = 0.0
   alloc_(cs%dx_dyT(isd:ied,jsd:jed))      ; cs%dx_dyT(:,:)  = 0.0
   alloc_(cs%dy_dxT(isd:ied,jsd:jed))      ; cs%dy_dxT(:,:)  = 0.0
   alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0
   alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0
  
   if (cs%Laplacian) then
     alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed))     ; cs%Kh_bg_xx(:,:) = 0.0
     alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0
     if (cs%bound_Kh .or. cs%better_bound_Kh) then
       alloc_(cs%Kh_Max_xx(isd:ied,jsd:jed)) ; cs%Kh_Max_xx(:,:) = 0.0
       alloc_(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xy(:,:) = 0.0
     endif
     if (cs%Smagorinsky_Kh) then
       alloc_(cs%Laplac2_const_xx(isd:ied,jsd:jed))     ; cs%Laplac2_const_xx(:,:) = 0.0
       alloc_(cs%Laplac2_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac2_const_xy(:,:) = 0.0
     endif
     if (cs%Leith_Kh) then
       alloc_(cs%Laplac3_const_xx(isd:ied,jsd:jed)) ; cs%Laplac3_const_xx(:,:) = 0.0
       alloc_(cs%Laplac3_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac3_const_xy(:,:) = 0.0
     endif
   endif
   alloc_(cs%reduction_xx(isd:ied,jsd:jed))     ; cs%reduction_xx(:,:) = 0.0
   alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0
  
   if (cs%anisotropic) then
     alloc_(cs%n1n2_h(isd:ied,jsd:jed)) ; cs%n1n2_h(:,:) = 0.0
     alloc_(cs%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; cs%n1n1_m_n2n2_h(:,:) = 0.0
     alloc_(cs%n1n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n2_q(:,:) = 0.0
     alloc_(cs%n1n1_m_n2n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n1_m_n2n2_q(:,:) = 0.0
     select case (aniso_mode)
       case (0)
         call align_aniso_tensor_to_grid(cs, aniso_grid_dir(1), aniso_grid_dir(2))
       case (1)
       ! call align_aniso_tensor_to_grid(CS, aniso_grid_dir(1), aniso_grid_dir(2))
       case (2)
         cs%dynamic_aniso = .true.
       case default
         call mom_error(fatal, "MOM_hor_visc: "//&
              "Runtime parameter ANISOTROPIC_MODE is out of range.")
     end select
   endif
  
   if (cs%use_Kh_bg_2d) then
     alloc_(cs%Kh_bg_2d(isd:ied,jsd:jed))     ; cs%Kh_bg_2d(:,:) = 0.0
     call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, &
                  'The filename containing a 2d map of "Kh".', &
                  default='KH_background_2d.nc')
     call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
     inputdir = slasher(inputdir)
     call mom_read_data(trim(inputdir)//trim(filename), 'Kh', cs%Kh_bg_2d, &
                        g%domain, timelevel=1, scale=us%m_to_L**2*us%T_to_s)
     call pass_var(cs%Kh_bg_2d, g%domain)
   endif
  
   if (cs%biharmonic) then
     alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0
     alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0
     alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed))  ; cs%Idxdy2u(:,:)  = 0.0
     alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb))  ; cs%Idxdy2v(:,:)  = 0.0
  
     alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed))     ; cs%Ah_bg_xx(:,:) = 0.0
     alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0
     if (cs%bound_Ah .or. cs%better_bound_Ah) then
       alloc_(cs%Ah_Max_xx(isd:ied,jsd:jed))     ; cs%Ah_Max_xx(:,:) = 0.0
       alloc_(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_Max_xy(:,:) = 0.0
     endif
     if (cs%Smagorinsky_Ah) then
       alloc_(cs%Biharm_const_xx(isd:ied,jsd:jed))     ; cs%Biharm_const_xx(:,:) = 0.0
       alloc_(cs%Biharm_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const_xy(:,:) = 0.0
       if (cs%bound_Coriolis) then
         alloc_(cs%Biharm_const2_xx(isd:ied,jsd:jed))     ; cs%Biharm_const2_xx(:,:) = 0.0
         alloc_(cs%Biharm_const2_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const2_xy(:,:) = 0.0
       endif
     endif
     if (cs%Leith_Ah) then
         alloc_(cs%biharm5_const_xx(isd:ied,jsd:jed)) ; cs%biharm5_const_xx(:,:) = 0.0
         alloc_(cs%biharm5_const_xy(isdb:iedb,jsdb:jedb)) ; cs%biharm5_const_xy(:,:) = 0.0
     endif
   endif
  
   do j=js-2,jeq+1 ; do i=is-2,ieq+1
     cs%dx2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%dy2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)
     cs%DX_dyBu(i,j) = g%dxBu(i,j)*g%IdyBu(i,j) ; cs%DY_dxBu(i,j) = g%dyBu(i,j)*g%IdxBu(i,j)
   enddo ; enddo
   do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
     cs%dx2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%dy2h(i,j) = g%dyT(i,j)*g%dyT(i,j)
     cs%DX_dyT(i,j) = g%dxT(i,j)*g%IdyT(i,j) ; cs%DY_dxT(i,j) = g%dyT(i,j)*g%IdxT(i,j)
   enddo ; enddo
  
   do j=jsq,jeq+1 ; do i=isq,ieq+1
     cs%reduction_xx(i,j) = 1.0
     if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
         (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &
       cs%reduction_xx(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
     if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &
         (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &
       cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / (g%dyCu(i-1,j))
     if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
         (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &
       cs%reduction_xx(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
     if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &
         (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &
       cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / (g%dxCv(i,j-1))
   enddo ; enddo
  
   do j=js-1,jeq ; do i=is-1,ieq
     cs%reduction_xy(i,j) = 1.0
     if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
         (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &
       cs%reduction_xy(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
     if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &
         (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &
       cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / (g%dyCu(i,j+1))
     if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
         (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &
       cs%reduction_xy(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
     if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &
         (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &
       cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / (g%dxCv(i+1,j))
   enddo ; enddo
  
   if (cs%Laplacian) then
    ! The 0.3 below was 0.4 in MOM1.10.  The change in hq requires
    ! this to be less than 1/3, rather than 1/2 as before.
     if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)
  
     ! Calculate and store the background viscosity at h-points
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       ! Static factors in the Smagorinsky and Leith schemes
       grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j) + cs%dy2h(i,j))
       grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
       if (cs%Smagorinsky_Kh) cs%Laplac2_const_xx(i,j) = smag_lap_const * grid_sp_h2
       if (cs%Leith_Kh)       cs%Laplac3_const_xx(i,j) = leith_lap_const * grid_sp_h3
       ! Maximum of constant background and MICOM viscosity
       cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))
  
       ! Use the larger of the above and values read from a file
       if (cs%use_Kh_bg_2d) cs%Kh_bg_xx(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xx(i,j))
  
       ! Use the larger of the above and a function of sin(latitude)
       if (kh_sin_lat>0.) then
         slat_fn = abs( sin( deg2rad * g%geoLatT(i,j) ) ) ** kh_pwr_of_sine
         cs%Kh_bg_xx(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xx(i,j))
       endif
  
       if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
         ! Limit the background viscosity to be numerically stable
         cs%Kh_Max_xx(i,j) = kh_limit * grid_sp_h2
         cs%Kh_bg_xx(i,j) = min(cs%Kh_bg_xx(i,j), cs%Kh_Max_xx(i,j))
       endif
     enddo ; enddo
  
     ! Calculate and store the background viscosity at q-points
     do j=js-1,jeq ; do i=is-1,ieq
       ! Static factors in the Smagorinsky and Leith schemes
       grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j) + cs%dy2q(i,j))
       grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
       if (cs%Smagorinsky_Kh) cs%Laplac2_const_xy(i,j) = smag_lap_const * grid_sp_q2
       if (cs%Leith_Kh)       cs%Laplac3_const_xy(i,j) = leith_lap_const * grid_sp_q3
       ! Maximum of constant background and MICOM viscosity
       cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))
  
       ! Use the larger of the above and values read from a file
       !### This expression uses inconsistent staggering
       if (cs%use_Kh_bg_2d) cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xy(i,j))
  
       ! Use the larger of the above and a function of sin(latitude)
       if (kh_sin_lat>0.) then
         slat_fn = abs( sin( deg2rad * g%geoLatBu(i,j) ) ) ** kh_pwr_of_sine
         cs%Kh_bg_xy(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xy(i,j))
       endif
  
       if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
         ! Limit the background viscosity to be numerically stable
         cs%Kh_Max_xy(i,j) = kh_limit * grid_sp_q2
         cs%Kh_bg_xy(i,j) = min(cs%Kh_bg_xy(i,j), cs%Kh_Max_xy(i,j))
       endif
     enddo ; enddo
   endif
  
   if (cs%biharmonic) then
  
     do j=js-1,jeq+1 ; do i=isq-1,ieq+1
       cs%Idx2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)
       cs%Idxdy2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))
     enddo ; enddo
     do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
       cs%Idx2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)
       cs%Idxdy2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))
     enddo ; enddo
  
     cs%Ah_bg_xy(:,:) = 0.0
    ! The 0.3 below was 0.4 in MOM1.10.  The change in hq requires
    ! this to be less than 1/3, rather than 1/2 as before.
     if (cs%better_bound_Ah .or. cs%bound_Ah) ah_limit = 0.3 / (dt*64.0)
     if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &
       boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j)+cs%dy2h(i,j))
       grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
  
       if (cs%Smagorinsky_Ah) then
         cs%Biharm_const_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)
         if (cs%bound_Coriolis) then
           fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &
                      abs(g%CoriolisBu(i-1,j)),   abs(g%CoriolisBu(i,j)))
           cs%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
                                      (fmax * boundcorconst)
         endif
       endif
       if (cs%Leith_Ah) then
          cs%biharm5_const_xx(i,j) = leith_bi_const * (grid_sp_h3 * grid_sp_h2)
       endif
       cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
       if (ah_time_scale > 0.) cs%Ah_bg_xx(i,j) = &
             max(cs%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / ah_time_scale)
       if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
         cs%Ah_Max_xx(i,j) = ah_limit * (grid_sp_h2 * grid_sp_h2)
         cs%Ah_bg_xx(i,j) = min(cs%Ah_bg_xx(i,j), cs%Ah_Max_xx(i,j))
       endif
     enddo ; enddo
     do j=js-1,jeq ; do i=is-1,ieq
       grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j)+cs%dy2q(i,j))
       grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
  
       if (cs%Smagorinsky_Ah) then
         cs%Biharm_const_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)
         if (cs%bound_Coriolis) then
           cs%Biharm_const2_xy(i,j) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
                                      (abs(g%CoriolisBu(i,j)) * boundcorconst)
         endif
       endif
       if (cs%Leith_Ah) then
          cs%biharm5_const_xy(i,j) = leith_bi_const * (grid_sp_q3 * grid_sp_q2)
       endif
  
       cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
       if (ah_time_scale > 0.) cs%Ah_bg_xy(i,j) = &
            max(cs%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / ah_time_scale)
       if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
         cs%Ah_Max_xy(i,j) = ah_limit * (grid_sp_q2 * grid_sp_q2)
         cs%Ah_bg_xy(i,j) = min(cs%Ah_bg_xy(i,j), cs%Ah_Max_xy(i,j))
       endif
     enddo ; enddo
   endif
  
   ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.
   if (cs%Laplacian .and. cs%better_bound_Kh) then
     idt = 1.0 / dt
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       denom = max( &
          (cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &
           max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
          (cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &
           max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
       cs%Kh_Max_xx(i,j) = 0.0
       if (denom > 0.0) &
         cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom
     enddo ; enddo
     do j=js-1,jeq ; do i=is-1,ieq
       denom = max( &
          (cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &
           max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
          (cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &
           max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
       cs%Kh_Max_xy(i,j) = 0.0
       if (denom > 0.0) &
         cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom
     enddo ; enddo
     if (cs%debug) then
       call hchksum(cs%Kh_Max_xx, "Kh_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
       call bchksum(cs%Kh_Max_xy, "Kh_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
     endif
   endif
  
   ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but
   ! empirically work for CS%bound_coef <~ 1.0
   if (cs%biharmonic .and. cs%better_bound_Ah) then
     idt = 1.0 / dt
     do j=js-1,jeq+1 ; do i=isq-1,ieq+1
       u0u(i,j) = (cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*cs%DY_dxT(i+1,j)*(g%IdyCu(i+1,j) + g%IdyCu(i,j))   + &
                                    cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) + &
                  cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
                                    cs%dx2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1)) ) )
  
       u0v(i,j) = (cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*cs%DX_dyT(i+1,j)*(g%IdxCv(i+1,j) + g%IdxCv(i+1,j-1)) + &
                                    cs%dy2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) )   + &
                  cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j))   + &
                                    cs%dx2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1)) ) )
     enddo ; enddo
     do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
       v0u(i,j) = (cs%Idxdy2v(i,j)*(cs%dy2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j))       + &
                                    cs%dy2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j)) ) + &
                  cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*cs%DY_dxT(i,j+1)*(g%IdyCu(i,j+1) + g%IdyCu(i-1,j+1))   + &
                                    cs%dx2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) )
  
       v0v(i,j) = (cs%Idxdy2v(i,j)*(cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j))   + &
                                    cs%dy2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j)) ) + &
                  cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j))   + &
                                    cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) )
     enddo ; enddo
  
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       denom = max( &
          (cs%dy2h(i,j) * &
           (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0u(i,j) + g%IdyCu(i-1,j)*u0u(i-1,j))  + &
            cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0u(i,j) + g%IdxCv(i,j-1)*v0u(i,j-1))) * &
           max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ),   &
          (cs%dx2h(i,j) * &
           (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0v(i,j) + g%IdyCu(i-1,j)*u0v(i-1,j))  + &
            cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0v(i,j) + g%IdxCv(i,j-1)*v0v(i,j-1))) * &
           max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
       cs%Ah_Max_xx(i,j) = 0.0
       if (denom > 0.0) &
         cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom
     enddo ; enddo
  
     do j=js-1,jeq ; do i=is-1,ieq
       denom = max( &
          (cs%dx2q(i,j) * &
           (cs%DX_dyBu(i,j)*(u0u(i,j+1)*g%IdxCu(i,j+1) + u0u(i,j)*g%IdxCu(i,j))  + &
            cs%DY_dxBu(i,j)*(v0u(i+1,j)*g%IdyCv(i+1,j) + v0u(i,j)*g%IdyCv(i,j))) * &
           max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ),    &
          (cs%dy2q(i,j) * &
           (cs%DX_dyBu(i,j)*(u0v(i,j+1)*g%IdxCu(i,j+1) + u0v(i,j)*g%IdxCu(i,j))  + &
            cs%DY_dxBu(i,j)*(v0v(i+1,j)*g%IdyCv(i+1,j) + v0v(i,j)*g%IdyCv(i,j))) * &
           max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
       cs%Ah_Max_xy(i,j) = 0.0
       if (denom > 0.0) &
         cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom
     enddo ; enddo
     if (cs%debug) then
       call hchksum(cs%Ah_Max_xx, "Ah_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
       call bchksum(cs%Ah_Max_xy, "Ah_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
     endif
   endif
  
   ! Register fields for output from this module.
  
   cs%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, time, &
       'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
  
   cs%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, time, &
       'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
  
   if (cs%biharmonic) then
     cs%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, time,    &
         'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T, &
         cmor_field_name='difmxybo',                                             &
         cmor_long_name='Ocean lateral biharmonic viscosity',                     &
         cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')
  
     cs%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, time, &
         'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)
   endif
  
   if (cs%Laplacian) then
     cs%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, time,   &
         'Laplacian Horizontal Viscosity at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
         cmor_field_name='difmxylo',                                             &
         cmor_long_name='Ocean lateral Laplacian viscosity',                     &
         cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')
  
     cs%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, time, &
         'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  
     if (cs%Leith_Kh) then
       cs%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, time, &
         'Vertical vorticity at q Points', 's-1', conversion=us%s_to_T)
  
       cs%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, time, &
         'Horizontal divergence at h Points', 's-1', conversion=us%s_to_T)
     endif
  
   endif
  
   if (cs%use_GME) then
       cs%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, time, &
         'GME coefficient at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  
       cs%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, time, &
         'GME coefficient at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  
       cs%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,time,&
       'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', &
       'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2)
   endif
  
   cs%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,time,&
       'Integral work done by lateral friction terms', &
       'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2)
  
   cs%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,time,      &
       'Depth integrated work done by lateral friction', &
       'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2, &
       cmor_field_name='dispkexyfo',                                                              &
       cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',&
       cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
  
   if (cs%Laplacian .or. get_all) then
   endif
  

References align_aniso_tensor_to_grid(), and mom_error_handler::mom_error().

Referenced by mom_dynamics_split_rk2::initialize_dyn_split_rk2(), mom_dynamics_unsplit::initialize_dyn_unsplit(), and mom_dynamics_unsplit_rk2::initialize_dyn_unsplit_rk2().

Here is the call graph for this function:

Here is the caller graph for this function: