MOM_set_viscosity.F90

Go to the documentation of this file.

module mom_set_visc
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - October 2006                      *
!*  Quadratic Bottom Drag by James Stephens and R. Hallberg.           *
!*                                                                     *
!*    This file contains the subroutine that calculates various values *
!*  related to the bottom boundary layer, such as the viscosity and    *
!*  thickness of the BBL (set_viscous_BBL).  This would also be the    *
!*  module in which other viscous quantities that are flow-independent *
!*  might be set.  This information is transmitted to other modules    *
!*  via a vertvisc type structure.                                     *
!*                                                                     *
!*    The same code is used for the two velocity components, by        *
!*  indirectly referencing the velocities and defining a handful of    *
!*  direction-specific defined variables.                              *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, frhatv, tauy                          *
!*    j    x ^ x ^ x   At >:  u, frhatu, taux                          *
!*    j    > o > o >   At o:  h                                        *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_debugging, only : uvchksum
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl, time_type
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type
use mom_kappa_shear, only : kappa_shear_is_used
use mom_cvmix_shear, only : cvmix_shear_is_used
use mom_io, only : vardesc, var_desc
use mom_restart, only : register_restart_field, mom_restart_cs
use mom_variables, only : thermo_var_ptrs
use mom_variables, only : vertvisc_type
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs
use mom_open_boundary, only : ocean_obc_type, obc_none, obc_direction_e
use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s
use mom_open_boundary, only : obc_segment_type
implicit none ; private

#include <MOM_memory.h>

public set_viscous_bbl, set_viscous_ml, set_visc_init, set_visc_end
public set_visc_register_restarts

type, public :: set_visc_cs ; private
  real    :: hbbl           ! The static bottom boundary layer thickness, in
                            ! the same units as thickness (m or kg m-2).
  real    :: cdrag          ! The quadratic drag coefficient.
  real    :: c_smag         ! The Laplacian Smagorinsky coefficient for
                            ! calculating the drag in channels.
  real    :: drag_bg_vel    ! An assumed unresolved background velocity for
                            ! calculating the bottom drag, in m s-1.
  real    :: bbl_thick_min  ! The minimum bottom boundary layer thickness in
                            ! the same units as thickness (m or kg m-2).
                            ! This might be Kv / (cdrag * drag_bg_vel) to give
                            ! Kv as the minimum near-bottom viscosity.
  real    :: htbl_shelf     ! A nominal thickness of the surface boundary layer
                            ! for use in calculating the near-surface velocity,
                            ! in units of m.
  real    :: htbl_shelf_min ! The minimum surface boundary layer thickness in m.
  real    :: kv_bbl_min     ! The minimum viscosities in the bottom and top
  real    :: kv_tbl_min     ! boundary layers, both in m2 s-1.

  logical :: bottomdraglaw  ! If true, the  bottom stress is calculated with a
                            ! drag law c_drag*|u|*u. The velocity magnitude
                            ! may be an assumed value or it may be based on the
                            ! actual velocity in the bottommost HBBL, depending
                            ! on whether linear_drag is true.
  logical :: bbl_use_eos    ! If true, use the equation of state in determining
                            ! the properties of the bottom boundary layer.
  logical :: linear_drag    ! If true, the drag law is cdrag*DRAG_BG_VEL*u.
  logical :: channel_drag   ! If true, the drag is exerted directly on each
                            ! layer according to what fraction of the bottom
                            ! they overlie.
  logical :: rino_mix       ! If true, use Richardson number dependent mixing.
  logical :: dynamic_viscous_ml  ! If true, use a bulk Richardson number criterion to
                            ! determine the mixed layer thickness for viscosity.
  real    :: bulk_ri_ml     ! The bulk mixed layer used to determine the
                            ! thickness of the viscous mixed layer.  Nondim.
  real    :: omega          !   The Earth's rotation rate, in s-1.
  real    :: ustar_min      ! A minimum value of ustar to avoid numerical
                            ! problems, in m s-1.  If the value is small enough,
                            ! this should not affect the solution.
  real    :: tke_decay      ! The ratio of the natural Ekman depth to the TKE
                            ! decay scale, nondimensional.
  real    :: omega_frac     !   When setting the decay scale for turbulence, use
                            ! this fraction of the absolute rotation rate blended
                            ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
  logical :: debug          ! If true, write verbose checksums for debugging purposes.
  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                            ! timing of diagnostic output.
  integer :: id_bbl_thick_u = -1, id_kv_bbl_u = -1
  integer :: id_bbl_thick_v = -1, id_kv_bbl_v = -1
  integer :: id_ray_u = -1, id_ray_v = -1
  integer :: id_nkml_visc_u = -1, id_nkml_visc_v = -1
  type(ocean_obc_type), pointer :: obc => null()
end type set_visc_cs

contains

!>   The following subroutine 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.
subroutine set_viscous_bbl(u, v, h, tv, visc, G, GV, CS)
  type(ocean_grid_type),    intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),  intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            intent(in)    :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                            intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  type(thermo_var_ptrs),    intent(in)    :: tv   !< A structure containing pointers to any
                                                  !! available thermodynamic fields. Absent fields
                                                  !! have NULL ptrs..
  type(vertvisc_type),      intent(inout) :: visc !< A structure containing vertical viscosities and
                                                  !! related fields.
  type(set_visc_cs),        pointer       :: CS   !< The control structure returned by a previous
                                                  !! call to vertvisc_init.
!   The following subroutine 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.
!
! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.  In the comments below,
!                the units of h are denoted as H.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (out)     visc - A structure containing vertical viscosities and related
!                   fields.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 vertvisc_init.
  real, dimension(SZIB_(G)) :: &
    ustar, &    !   The bottom friction velocity, in m s-1.
    T_EOS, &    !   The temperature used to calculate the partial derivatives
                ! of density with T and S, in deg C.
    s_eos, &    !   The salinity used to calculate the partial derivatives
                ! of density with T and S, in PSU.
    dr_dt, &    !   Partial derivative of the density in the bottom boundary
                ! layer with temperature, in units of kg m-3 K-1.
    dr_ds, &    !   Partial derivative of the density in the bottom boundary
                ! layer with salinity, in units of kg m-3 psu-1.
    press       !   The pressure at which dR_dT and dR_dS are evaluated, in Pa.
  real :: htot             ! Sum of the layer thicknesses up to some
                           ! point, in H (i.e., m or kg m-2).
  real :: htot_vel         ! Sum of the layer thicknesses up to some
                           ! point, in H (i.e., m or kg m-2).

  real :: Rhtot            ! Running sum of thicknesses times the
                           ! layer potential densities in H kg m-3.
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    D_u, &      ! Bottom depth interpolated to u points, in m.
    mask_u      ! A mask that disables any contributions from u points that
                ! are land or past open boundary conditions, nondim., 0 or 1.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    D_v, &      ! Bottom depth interpolated to v points, in m.
    mask_v      ! A mask that disables any contributions from v points that
                ! are land or past open boundary conditions, nondim., 0 or 1.
  real, dimension(SZIB_(G),SZK_(G)) :: &
    h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
                ! second order accurate estimate based on the previous velocity
                ! direction, in H.
    h_vel, &    ! Arithmetic mean of the layer thicknesses adjacent to a
                ! velocity point, in H.
    t_vel, &    ! Arithmetic mean of the layer temperatures adjacent to a
                ! velocity point, in deg C.
    s_vel, &    ! Arithmetic mean of the layer salinities adjacent to a
                ! velocity point, in PSU.
    rml_vel     ! Arithmetic mean of the layer coordinate densities adjacent
                ! to a velocity point, in kg m-3.

  real :: h_vel_pos        ! The arithmetic mean thickness at a velocity point
                           ! plus H_neglect to avoid 0 values, in H.
  real :: ustarsq          ! 400 times the square of ustar, times
                           ! Rho0 divided by G_Earth and the conversion
                           ! from m to thickness units, in kg m-2 or kg2 m-5.
  real :: cdrag_sqrt       ! Square root of the drag coefficient, nd.
  real :: oldfn            ! The integrated energy required to
                           ! entrain up to the bottom of the layer,
                           ! divided by G_Earth, in H kg m-3.
  real :: Dfn              ! The increment in oldfn for entraining
                           ! the layer, in H kg m-3.
  real :: Dh               ! The increment in layer thickness from
                           ! the present layer, in H.
  real :: bbl_thick        ! The thickness of the bottom boundary layer in m.
  real :: C2f              ! C2f = 2*f at velocity points.

  real :: U_bg_sq          ! The square of an assumed background
                           ! velocity, for calculating the mean
                           ! magnitude near the bottom for use in the
                           ! quadratic bottom drag, in m2.
  real :: hwtot            ! Sum of the thicknesses used to calculate
                           ! the near-bottom velocity magnitude, in H.
  real :: hutot            ! Running sum of thicknesses times the
                           ! velocity magnitudes, in H m s-1.
  real :: Thtot            ! Running sum of thickness times temperature, in H C.
  real :: Shtot            ! Running sum of thickness times salinity, in H psu.
  real :: hweight          ! The thickness of a layer that is within Hbbl
                           ! of the bottom, in H.
  real :: v_at_u, u_at_v   ! v at a u point or vice versa, m s-1.
  real :: Rho0x400_G       ! 400*Rho0/G_Earth, in kg s2 m-4.  The 400 is a
                           ! constant proposed by Killworth and Edwards, 1999.
  real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
    Rml                    ! The mixed layer coordinate density, in kg m-3.
  real :: p_ref(szi_(g))   !   The pressure used to calculate the coordinate
                           ! density, in Pa (usually set to 2e7 Pa = 2000 dbar).

  ! The units H in the following are thickness units - typically m or kg m-2.
  real :: D_vel            ! The bottom depth at a velocity point, in H.
  real :: Dp, Dm           ! The depths at the edges of a velocity cell, in H.
  real :: a                ! a is the curvature of the bottom depth across a
                           ! cell, times the cell width squared, in H.
  real :: a_3, a_12, C24_a ! a/3, a/12, and 24/a, in H, H, and H-1.
  real :: slope            ! The absolute value of the bottom depth slope across
                           ! a cell times the cell width, in H.
  real :: apb_4a, ax2_3apb ! Various nondimensional ratios of a and slope.
  real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of a and slope with units of H-1.
  ! All of the following "volumes" have units of meters as they are normalized
  ! by the full horizontal area of a velocity cell.
  real :: Vol_open         ! The cell volume above which it is open, in H.
  real :: Vol_direct       ! With less than Vol_direct (in H), there is a direct
                           ! solution of a cubic equation for L.
  real :: Vol_2_reg        ! The cell volume above which there are two separate
                           ! open areas that must be integrated, in H.
  real :: vol              ! The volume below the interface whose normalized
                           ! width is being sought, in H.
  real :: vol_below        ! The volume below the interface below the one that
                           ! is currently under consideration, in H.
  real :: Vol_err          ! The error in the volume with the latest estimate of
                           ! L, or the error for the interface below, in H.
  real :: Vol_quit         ! The volume error below which to quit iterating, in H.
  real :: Vol_tol          ! A volume error tolerance, in H.
  real :: L(szk_(g)+1)     ! The fraction of the full cell width that is open at
                           ! the depth of each interface, nondimensional.
  real :: L_direct         ! The value of L above volume Vol_direct, nondim.
  real :: L_max, L_min     ! Upper and lower bounds on the correct value for L.
  real :: Vol_err_max      ! The volume errors for the upper and lower bounds on
  real :: Vol_err_min      ! the correct value for L, in H.
  real :: Vol_0            ! A deeper volume with known width L0, in H.
  real :: L0               ! The value of L above volume Vol_0, nondim.
  real :: dVol             ! vol - Vol_0, in H.
  real :: dV_dL2           ! The partial derivative of volume with L squared
                           ! evaluated at L=L0, in H.
  real :: h_neglect        ! A thickness that is so small it is usually lost
                           ! in roundoff and can be neglected, in H.
  real :: ustH             ! ustar converted to units of H s-1.
  real :: root             ! A temporary variable with units of H s-1.
  real :: H_to_m, m_to_H   ! Local copies of unit conversion factors.

  real :: Cell_width       ! The transverse width of the velocity cell, in m.
  real :: Rayleigh         ! A nondimensional value that is multiplied by the
                           ! layer's velocity magnitude to give the Rayleigh
                           ! drag velocity.
  real :: gam              ! The ratio of the change in the open interface width
                           ! to the open interface width atop a cell, nondim.
  real :: BBL_frac         ! The fraction of a layer's drag that goes into the
                           ! viscous bottom boundary layer, nondim.
  real :: BBL_visc_frac    ! The fraction of all the drag that is expressed as
                           ! a viscous bottom boundary layer, nondim.
  real, parameter :: C1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0
  real :: C2pi_3           ! An irrational constant, 2/3 pi.
  real :: tmp              ! A temporary variable.
  real :: tmp_val_m1_to_p1
  logical :: use_BBL_EOS, do_i(szib_(g))
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml
  integer :: itt, maxitt=20
  type(ocean_obc_type), pointer :: OBC => null()

  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
  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
  h_neglect = gv%H_subroundoff
  rho0x400_g = 400.0*(gv%Rho0/gv%g_Earth)*gv%m_to_H
  vol_quit = 0.9*gv%Angstrom + h_neglect
  h_to_m = gv%H_to_m ; m_to_h = gv%m_to_H
  c2pi_3 = 8.0*atan(1.0)/3.0

  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(BBL): "//&
         "Module must be initialized before it is used.")
  if (.not.cs%bottomdraglaw) return

  use_bbl_eos = associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
  obc => cs%OBC

  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
  cdrag_sqrt=sqrt(cs%cdrag)
  k2 = max(nkmb+1, 2)

!  With a linear drag law, the friction velocity is already known.
!  if (CS%linear_drag) ustar(:) = cdrag_sqrt*CS%drag_bg_vel

  if ((nkml>0) .and. .not.use_bbl_eos) then
    do i=g%IscB,g%IecB+1 ; p_ref(i) = tv%P_ref ; enddo
    !$OMP parallel do default(shared)
    do j=jsq,jeq+1 ; do k=1,nkmb
      call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, &
                      rml(:,j,k), isq, ieq-isq+2, tv%eqn_of_state)
    enddo ; enddo
  endif

  !$OMP parallel do default(shared)
  do j=js-1,je ; do i=is-1,ie+1
    d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
    mask_v(i,j) = g%mask2dCv(i,j)
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je+1 ; do i=is-1,ie
    d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
    mask_u(i,j) = g%mask2dCu(i,j)
  enddo ; enddo

  if (associated(obc)) then ; do n=1,obc%number_of_segments
    if (.not. obc%segment(n)%on_pe) cycle
    ! Use a one-sided projection of bottom depths at OBC points.
    i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
    if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
      do i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
        if (obc%segment(n)%direction == obc_direction_n) d_v(i,j) = g%bathyT(i,j)
        if (obc%segment(n)%direction == obc_direction_s) d_v(i,j) = g%bathyT(i,j+1)
      enddo
    elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
      do j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
        if (obc%segment(n)%direction == obc_direction_e) d_u(i,j) = g%bathyT(i,j)
        if (obc%segment(n)%direction == obc_direction_w) d_u(i,j) = g%bathyT(i+1,j)
      enddo
    endif
  enddo ; endif
  if (associated(obc)) then ; do n=1,obc%number_of_segments
    ! Now project bottom depths across cell-corner points in the OBCs.  The two
    ! projections have to occur in sequence and can not be combined easily.
    if (.not. obc%segment(n)%on_pe) cycle
    ! Use a one-sided projection of bottom depths at OBC points.
    i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
    if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
      do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
        if (obc%segment(n)%direction == obc_direction_n) then
          d_u(i,j+1) = d_u(i,j) ;  mask_u(i,j+1) = 0.0
        elseif (obc%segment(n)%direction == obc_direction_s) then
          d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
        endif
      enddo
    elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
      do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
        if (obc%segment(n)%direction == obc_direction_e) then
          d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
        elseif (obc%segment(n)%direction == obc_direction_w) then
          d_v(i,j) = d_v(i+1,j) ;  mask_v(i,j) = 0.0
        endif
      enddo
    endif
  enddo ; endif

  if (.not.use_bbl_eos) rml_vel(:,:) = 0.0

!$OMP parallel do default(none) shared(u, v, h, tv, visc, G, GV, CS, Rml, is, ie, js, je,  &
!$OMP                                  nz, Isq, Ieq, Jsq, Jeq, nkmb, h_neglect, Rho0x400_G,&
!$OMP                                  C2pi_3, U_bg_sq, cdrag_sqrt,K2,use_BBL_EOS,OBC,     &
!$OMP                                  maxitt,nkml,m_to_H,H_to_m,Vol_quit,D_u,D_v,mask_u,mask_v) &
!$OMP                          private(do_i,h_at_vel,htot_vel,hwtot,hutot,Thtot,Shtot,     &
!$OMP                                  hweight,v_at_u,u_at_v,ustar,T_EOS,S_EOS,press,      &
!$OMP                                  dR_dT, dR_dS,ustarsq,htot,T_vel,S_vel,Rml_vel,      &
!$OMP                                  oldfn,Dfn,Dh,Rhtot,C2f,ustH,root,bbl_thick,         &
!$OMP                                  D_vel,tmp,Dp,Dm,a_3,a,a_12,slope,Vol_open,Vol_2_reg,&
!$OMP                                  C24_a,apb_4a,Iapb,a2x48_apb3,ax2_3apb,Vol_direct,   &
!$OMP                                  L_direct,Ibma_2,L,vol,vol_below,Vol_err,h_vel_pos,  &
!$OMP                                  BBL_visc_frac,h_vel,L0,Vol_0,dV_dL2,dVol,L_max,     &
!$OMP                                  L_min,Vol_err_min,Vol_err_max,BBL_frac,Cell_width,  &
!$OMP                                  gam,Rayleigh, Vol_tol, tmp_val_m1_to_p1)
  do j=g%JscB,g%JecB ; do m=1,2

    if (m==1) then
      if (j<g%Jsc) cycle
      is = g%IscB ; ie = g%IecB
      do i=is,ie
        do_i(i) = .false.
        if (g%mask2dCu(i,j) > 0) do_i(i) = .true.
      enddo
    else
      is = g%isc ; ie = g%iec
      do i=is,ie
        do_i(i) = .false.
        if (g%mask2dCv(i,j) > 0) do_i(i) = .true.
      enddo
    endif

    if (m==1) then
      do k=1,nz ; do i=is,ie
        if (do_i(i)) then
          if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
            h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
                            (h(i,j,k) + h(i+1,j,k) + h_neglect)
          else
            h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
          endif
        endif
        h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
      enddo ; enddo
      if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
        ! Perhaps these should be thickness weighted.
        t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
        s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
      enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
        rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
      enddo ; enddo ; endif
    else
      do k=1,nz ; do i=is,ie
        if (do_i(i)) then
          if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
            h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
                            (h(i,j,k) + h(i,j+1,k) + h_neglect)
          else
            h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
          endif
        endif
        h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
      enddo ; enddo
      if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
        t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
        s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
      enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
        rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
      enddo ; enddo ; endif
    endif

    if (associated(obc)) then ; if (obc%number_of_segments > 0) then
      ! Apply a zero gradient projection of thickness across OBC points.
      if (m==1) then
        do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
          if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
            do k=1,nz
              h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
            enddo
            if (use_bbl_eos) then
              do k=1,nz
                t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
              enddo
            else
              do k=1,nkmb
                rml_vel(i,k) = rml(i,j,k)
              enddo
            endif
          elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
            do k=1,nz
              h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
            enddo
            if (use_bbl_eos) then
              do k=1,nz
                t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
              enddo
            else
              do k=1,nkmb
                rml_vel(i,k) = rml(i+1,j,k)
              enddo
            endif
          endif
        endif ; enddo
      else
        do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
          if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
            do k=1,nz
              h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
            enddo
            if (use_bbl_eos) then
              do k=1,nz
                t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
              enddo
            else
              do k=1,nkmb
                rml_vel(i,k) = rml(i,j,k)
              enddo
            endif
          elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
            do k=1,nz
              h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
            enddo
            if (use_bbl_eos) then
              do k=1,nz
                t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
              enddo
            else
              do k=1,nkmb
                rml_vel(i,k) = rml(i,j+1,k)
              enddo
            endif
          endif
        endif ; enddo
      endif
    endif ; endif

    if (use_bbl_eos .or. .not.cs%linear_drag) then
      do i=is,ie ; if (do_i(i)) then
!   This block of code calculates the mean velocity magnitude over
! the bottommost CS%Hbbl of the water column for determining
! the quadratic bottom drag.
        htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
        thtot = 0.0 ; shtot = 0.0
        do k=nz,1,-1

          if (htot_vel>=cs%Hbbl) exit ! terminate the k loop

          hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
          if (hweight < 1.5*gv%Angstrom + h_neglect) cycle

          htot_vel  = htot_vel + h_at_vel(i,k)
          hwtot = hwtot + hweight

          if ((.not.cs%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
            v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v)
            hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
                                           v_at_u*v_at_u + u_bg_sq)
          else
            u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u)
            hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
                                           u_at_v*u_at_v + u_bg_sq)
          endif ; endif

          if (use_bbl_eos .and. (hweight >= 0.0)) then
            thtot = thtot + hweight * t_vel(i,k)
            shtot = shtot + hweight * s_vel(i,k)
          endif
        enddo ! end of k loop

        if (.not.cs%linear_drag .and. (hwtot > 0.0)) then
          ustar(i) = cdrag_sqrt*hutot/hwtot
        else
          ustar(i) = cdrag_sqrt*cs%drag_bg_vel
        endif

        if (use_bbl_eos) then ; if (hwtot > 0.0) then
          t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
        else
          t_eos(i) = 0.0 ; s_eos(i) = 0.0
        endif ; endif
      endif ; enddo
    else
      do i=is,ie ; ustar(i) = cdrag_sqrt*cs%drag_bg_vel ; enddo
    endif ! Not linear_drag

    if (use_bbl_eos) then
      do i=is,ie
        press(i) = 0.0 ! or = fluxes%p_surf(i,j)
        if (.not.do_i(i)) then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ; endif
      enddo
      do k=1,nz ; do i=is,ie
        press(i) = press(i) + gv%H_to_Pa * h_vel(i,k)
      enddo ; enddo
      call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
                                    is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
    endif

    do i=is,ie ; if (do_i(i)) then
!  The 400.0 in this expression is the square of a constant proposed
!  by Killworth and Edwards, 1999, in equation (2.20).
      ustarsq = rho0x400_g * ustar(i)**2
      htot = 0.0

!   This block of code calculates the thickness of a stratification
! limited bottom boundary layer, using the prescription from
! Killworth and Edwards, 1999, as described in Stephens and Hallberg
! 2000 (unpublished and lost manuscript).
      if (use_bbl_eos) then
        thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
        do k=nz,2,-1
          if (h_at_vel(i,k) <= 0.0) cycle

          oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
                  dr_ds(i)*(shtot - s_vel(i,k)*htot)
          if (oldfn >= ustarsq) exit

          dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
                 dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
                (h_at_vel(i,k) + htot)

          if ((oldfn + dfn) <= ustarsq) then
            dh = h_at_vel(i,k)
          else
            dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
          endif

          htot = htot + dh
          thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
        enddo
        if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0) then
          ! Layer 1 might be part of the BBL.
          if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
              dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) &
            htot = htot + h_at_vel(i,1)
        endif ! Examination of layer 1.
      else  ! Use Rlay and/or the coordinate density as density variables.
        rhtot = 0.0
        do k=nz,k2,-1
          oldfn = rhtot - gv%Rlay(k)*htot
          dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)

          if (oldfn >= ustarsq) then
            cycle
          else if ((oldfn + dfn) <= ustarsq) then
            dh = h_at_vel(i,k)
          else
            dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
          endif

          htot = htot + dh
          rhtot = rhtot + gv%Rlay(k)*dh
        enddo
        if (nkml>0) then
          do k=nkmb,2,-1
            oldfn = rhtot - rml_vel(i,k)*htot
            dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)

            if (oldfn >= ustarsq) then
              cycle
            else if ((oldfn + dfn) <= ustarsq) then
              dh = h_at_vel(i,k)
            else
              dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
            endif

            htot = htot + dh
            rhtot = rhtot + rml_vel(i,k)*dh
          enddo
          if (rhtot - rml_vel(i,1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
        else
          if (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
        endif
      endif ! use_BBL_EOS

! The Coriolis limit is 0.5*ustar/f. The buoyancy limit here is htot.
! The  bottom boundary layer thickness is found by solving the same
! equation as in Killworth and Edwards:    (h/h_f)^2 + h/h_N = 1.

      if (m==1) then ; c2f = (g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j))
      else ; c2f = (g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) ; endif

      if (cs%cdrag * u_bg_sq <= 0.0) then
        ! This avoids NaNs and overflows, and could be used in all cases,
        ! but is not bitwise identical to the current code.
        usth = ustar(i)*m_to_h ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
        if (htot*usth <= (cs%BBL_thick_min+h_neglect) * (0.5*usth + root)) then
          bbl_thick = cs%BBL_thick_min
        else
          bbl_thick = (htot * usth) / (0.5*usth + root)
        endif
      else
        bbl_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
          ((ustar(i)*ustar(i)) * (m_to_h**2) )))

        if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
      endif
! If there is Richardson number dependent mixing, that determines
! the vertical extent of the bottom boundary layer, and there is no
! need to set that scale here.  In fact, viscously reducing the
! shears over an excessively large region reduces the efficacy of
! the Richardson number dependent mixing.
      if ((bbl_thick > 0.5*cs%Hbbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%Hbbl

      if (cs%Channel_drag) then
        ! The drag within the bottommost bbl_thick is applied as a part of
        ! an enhanced bottom viscosity, while above this the drag is applied
        ! directly to the layers in question as a Rayleigh drag term.
        if (m==1) then
          d_vel = d_u(i,j)
          tmp = g%mask2dCu(i,j+1) * d_u(i,j+1)
          dp = 2.0 * d_vel * tmp / (d_vel + tmp)
          tmp = g%mask2dCu(i,j-1) * d_u(i,j-1)
          dm = 2.0 * d_vel * tmp / (d_vel + tmp)
        else
          d_vel = d_v(i,j)
          tmp = g%mask2dCv(i+1,j) * d_v(i+1,j)
          dp = 2.0 * d_vel * tmp / (d_vel + tmp)
          tmp = g%mask2dCv(i-1,j) * d_v(i-1,j)
          dm = 2.0 * d_vel * tmp / (d_vel + tmp)
        endif
        if (dm > dp) then ; tmp = dp ; dp = dm ; dm = tmp ; endif

        ! Convert the D's to the units of thickness.
        dp = m_to_h*dp ; dm = m_to_h*dm ; d_vel = m_to_h*d_vel

        a_3 = (dp + dm - 2.0*d_vel) ; a = 3.0*a_3 ; a_12 = 0.25*a_3
        slope = dp - dm
        ! If the curvature is small enough, there is no reason not to assume
        ! a uniformly sloping or flat bottom.
        if (abs(a) < 1e-2*(slope + cs%BBL_thick_min)) a = 0.0
        ! Each cell extends from x=-1/2 to 1/2, and has a topography
        ! given by D(x) = a*x^2 + b*x + D - a/12.

        ! Calculate the volume above which the entire cell is open and the
        ! other volumes at which the equation that is solved for L changes.
        if (a > 0.0) then
          if (slope >= a) then
            vol_open = d_vel - dm ; vol_2_reg = vol_open
          else
            tmp = slope/a
            vol_open = 0.25*slope*tmp + c1_12*a
            vol_2_reg = 0.5*tmp**2 * (a - c1_3*slope)
          endif
          ! Define some combinations of a & b for later use.
          c24_a = 24.0/a ; iapb = 1.0/(a+slope)
          apb_4a = (slope+a)/(4.0*a) ; a2x48_apb3 = (48.0*(a*a))*(iapb**3)
          ax2_3apb = 2.0*c1_3*a*iapb
        elseif (a == 0.0) then
          vol_open = 0.5*slope
          if (slope > 0) iapb = 1.0/slope
        else ! a < 0.0
          vol_open = d_vel - dm
          if (slope >= -a) then
            iapb = 1.0e30 ; if (slope+a /= 0.0) iapb = 1.0/(a+slope)
            vol_direct = 0.0 ; l_direct = 0.0 ; c24_a = 0.0
          else
            c24_a = 24.0/a ; iapb = 1.0/(a+slope)
            l_direct = 1.0 + slope/a ! L_direct < 1 because a < 0
            vol_direct = -c1_6*a*l_direct**3
          endif
          ibma_2 = 2.0 / (slope - a)
        endif

        l(nz+1) = 0.0 ; vol = 0.0 ; vol_err = 0.0 ; bbl_visc_frac = 0.0
        ! Determine the normalized open length at each interface.
        do k=nz,1,-1
          vol_below = vol

          vol = vol + h_vel(i,k)
          h_vel_pos = h_vel(i,k) + h_neglect

          if (vol >= vol_open) then ; l(k) = 1.0
          elseif (a == 0) then ! The bottom has no curvature.
            l(k) = sqrt(2.0*vol*iapb)
          elseif (a > 0) then
            ! There may be a minimum depth, and there are
            ! analytic expressions for L for all cases.
            if (vol < vol_2_reg) then
              ! In this case, there is a contiguous open region and
              !   vol = 0.5*L^2*(slope + a/3*(3-4L)).
              if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7?
                ! There is a very good approximation here for massless layers.
                l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2_3apb*l0)
              else
                l(k) = apb_4a * (1.0 - &
                         2.0 * cos(c1_3*acos(a2x48_apb3*vol - 1.0) - c2pi_3))
              endif
              ! To check the answers.
              ! Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
            else ! There are two separate open regions.
              !   vol = slope^2/4a + a/12 - (a/12)*(1-L)^2*(1+2L)
              ! At the deepest volume, L = slope/a, at the top L = 1.
              !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_a*(Vol_open - vol)) - C2pi_3)
              tmp_val_m1_to_p1 = 1.0 - c24_a*(vol_open - vol)
              tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
              l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
              ! To check the answers.
              ! Vol_err = Vol_open - a_12*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol
            endif
          else ! a < 0.
            if (vol <= vol_direct) then
              ! Both edges of the cell are bounded by walls.
              l(k) = (-0.25*c24_a*vol)**c1_3
            else
              ! x_R is at 1/2 but x_L is in the interior & L is found by solving
              !   vol = 0.5*L^2*(slope + a/3*(3-4L))

              !  Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + a_3*(3.0-4.0*L(K+1))) - vol_below
              ! Change to ...
              !   if (min(Vol_below + Vol_err, vol) <= Vol_direct) then ?
              if (vol_below + vol_err <= vol_direct) then
                l0 = l_direct ; vol_0 = vol_direct
              else
                l0 = l(k+1) ; vol_0 = vol_below + vol_err
                ! Change to   Vol_0 = min(Vol_below + Vol_err, vol) ?
              endif

              !   Try a relatively simple solution that usually works well
              ! for massless layers.
              dv_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol_0)
           !  dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol_0, 0.0)

           !### The following code is more robust when GV%Angstrom=0, but it
           !### changes answers.
           !   Vol_tol = max(0.5*GV%Angstrom + GV%H_subroundoff, 1e-14*vol)
           !   Vol_quit = max(0.9*GV%Angstrom + GV%H_subroundoff, 1e-14*vol)

           !   if (dVol <= 0.0) then
           !     L(K) = L0
           !     Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
           !   elseif (a*a*dVol**3 < Vol_tol*dV_dL2**2 * &
           !                     (dV_dL2*Vol_tol - 2.0*a*L0*dVol)) then
              if (a*a*dvol**3 < gv%Angstrom*dv_dl2**2 * &
                                (0.25*dv_dl2*gv%Angstrom - a*l0*dvol)) then
                ! One iteration of Newton's method should give an estimate
                ! that is accurate to within Vol_tol.
                l(k) = sqrt(l0*l0 + dvol / dv_dl2)
                vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
              else
                if (dv_dl2*(1.0-l0*l0) < dvol + &
                    dv_dl2 * (vol_open - vol)*ibma_2) then
                  l_max = sqrt(1.0 - (vol_open - vol)*ibma_2)
                else
                  l_max = sqrt(l0*l0 + dvol / dv_dl2)
                endif
                l_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l_max))

                vol_err_min = 0.5*(l_min**2)*(slope + a_3*(3.0-4.0*l_min)) - vol
                vol_err_max = 0.5*(l_max**2)*(slope + a_3*(3.0-4.0*l_max)) - vol
           !    if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
                if (abs(vol_err_min) <= vol_quit) then
                  l(k) = l_min ; vol_err = vol_err_min
                else
                  l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
                              (vol_err_max - vol_err_min))
                  do itt=1,maxitt
                    vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
                    if (abs(vol_err) <= vol_quit) exit
                    ! Take a Newton's method iteration. This equation has proven
                    ! robust enough not to need bracketing.
                    l(k) = l(k) - vol_err / (l(k)* (slope + a - 2.0*a*l(k)))
                    ! This would be a Newton's method iteration for L^2:
                    !   L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+a) - a*L(K)))
                  enddo
                endif ! end of iterative solver
              endif ! end of 1-boundary alternatives.
            endif ! end of a<0 cases.
          endif

          ! Determine the drag contributing to the bottom boundary layer
          ! and the Raleigh drag that acting 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
              bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
            else
              bbl_frac = 0.0
            endif

            if (m==1) then ; cell_width = g%dy_Cu(i,j)
            else ; cell_width = g%dx_Cv(i,j) ; endif
            gam = 1.0 - l(k+1)/l(k)
            rayleigh = cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
                (12.0*cs%c_Smag*h_vel_pos) /  (12.0*cs%c_Smag*h_vel_pos + m_to_h * &
                 cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
          else ! This layer feels no drag.
            rayleigh = 0.0
          endif

          if (m==1) then
            if (rayleigh > 0.0) then
              v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v)
              visc%Ray_u(i,j,k) = rayleigh*sqrt(u(i,j,k)*u(i,j,k) + &
                                                v_at_u*v_at_u + u_bg_sq)
            else ; visc%Ray_u(i,j,k) = 0.0 ; endif
          else
            if (rayleigh > 0.0) then
              u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u)
              visc%Ray_v(i,j,k) = rayleigh*sqrt(v(i,j,k)*v(i,j,k) + &
                                                u_at_v*u_at_v + u_bg_sq)
            else ; visc%Ray_v(i,j,k) = 0.0 ; endif
          endif

        enddo ! k loop to determine L(K).

        bbl_thick = bbl_thick * h_to_m
        if (m==1) then
          visc%kv_bbl_u(i,j) = max(cs%KV_BBL_min, &
                                   cdrag_sqrt*ustar(i)*bbl_thick*bbl_visc_frac)
          visc%bbl_thick_u(i,j) = bbl_thick
        else
          visc%kv_bbl_v(i,j) = max(cs%KV_BBL_min, &
                                   cdrag_sqrt*ustar(i)*bbl_thick*bbl_visc_frac)
          visc%bbl_thick_v(i,j) = bbl_thick
        endif

      else ! Not Channel_drag.
!   Here the near-bottom viscosity is set to a value which will give
! the correct stress when the shear occurs over bbl_thick.
        bbl_thick = bbl_thick * h_to_m
        if (m==1) then
          visc%kv_bbl_u(i,j) = max(cs%KV_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick)
          visc%bbl_thick_u(i,j) = bbl_thick
        else
          visc%kv_bbl_v(i,j) = max(cs%KV_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick)
          visc%bbl_thick_v(i,j) = bbl_thick
        endif
      endif
    endif ; enddo ! end of i loop
  enddo ; enddo ! end of m & j loops

! Offer diagnostics for averaging
  if (cs%id_bbl_thick_u > 0) &
    call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
  if (cs%id_kv_bbl_u > 0) &
    call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
  if (cs%id_bbl_thick_v > 0) &
    call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
  if (cs%id_kv_bbl_v > 0) &
    call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
  if (cs%id_Ray_u > 0) &
    call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
  if (cs%id_Ray_v > 0) &
    call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)

  if (cs%debug) then
    if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) &
        call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI,haloshift=0)
    if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) &
        call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI,haloshift=0)
    if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) &
        call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, &
                      visc%bbl_thick_v, g%HI,haloshift=0)
  endif

end subroutine set_viscous_bbl

!> This subroutine finds a thickness-weighted value of v at the u-points.
function set_v_at_u(v, h, G, i, j, k, mask2dCv)
  type(ocean_grid_type),                     intent(in) :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  integer,                                   intent(in) :: i, j, k
  real, dimension(SZI_(G),SZJB_(G)),         intent(in) :: mask2dCv
  real                                                  :: set_v_at_u

  ! This subroutine finds a thickness-weighted value of v at the u-points.
  real :: hwt(4)           ! Masked weights used to average v onto u, in H.
  real :: hwt_tot          ! The sum of the masked thicknesses, in H.

  hwt(1) = (h(i,j-1,k) + h(i,j,k)) * mask2dcv(i,j-1)
  hwt(2) = (h(i+1,j-1,k) + h(i+1,j,k)) * mask2dcv(i+1,j-1)
  hwt(3) = (h(i,j,k) + h(i,j+1,k)) * mask2dcv(i,j)
  hwt(4) = (h(i+1,j,k) + h(i+1,j+1,k)) * mask2dcv(i+1,j)

  hwt_tot = (hwt(1) + hwt(4)) + (hwt(2) + hwt(3))
  set_v_at_u = 0.0
  if (hwt_tot > 0.0) set_v_at_u = &
          ((hwt(3) * v(i,j,k) + hwt(2) * v(i+1,j-1,k)) + &
           (hwt(4) * v(i+1,j,k) + hwt(1) * v(i,j-1,k))) / hwt_tot
end function set_v_at_u

!> This subroutine finds a thickness-weighted value of u at the v-points.
function set_u_at_v(u, h, G, i, j, k, mask2dCu)
  type(ocean_grid_type),                  intent(in) :: G    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real, dimension(SZIB_(G),SZJ_(G)),  intent(in) :: mask2dCu
  integer,                                intent(in) :: i, j, k
  real                                               :: set_u_at_v

  ! This subroutine finds a thickness-weighted value of u at the v-points.
  real :: hwt(4)           ! Masked weights used to average u onto v, in H.
  real :: hwt_tot          ! The sum of the masked thicknesses, in H.

  hwt(1) = (h(i-1,j,k) + h(i,j,k)) * mask2dcu(i-1,j)
  hwt(2) = (h(i,j,k) + h(i+1,j,k)) * mask2dcu(i,j)
  hwt(3) = (h(i-1,j+1,k) + h(i,j+1,k)) * mask2dcu(i-1,j+1)
  hwt(4) = (h(i,j+1,k) + h(i+1,j+1,k)) * mask2dcu(i,j+1)

  hwt_tot = (hwt(1) + hwt(4)) + (hwt(2) + hwt(3))
  set_u_at_v = 0.0
  if (hwt_tot > 0.0) set_u_at_v = &
          ((hwt(2) * u(i,j,k) + hwt(3) * u(i-1,j+1,k)) + &
           (hwt(1) * u(i-1,j,k) + hwt(4) * u(i,j+1,k))) / hwt_tot
end function set_u_at_v

!>   The following subroutine calculates the thickness of the surface boundary
!! layer for applying an elevated viscosity.  A bulk Richardson criterion or
!! the thickness of the topmost NKML layers (with a bulk mixed layer) are
!! currently used.  The thicknesses are given in terms of fractional layers, so
!! that this thickness will move as the thickness of the topmost layers change.
subroutine set_viscous_ml(u, v, h, tv, fluxes, visc, dt, G, GV, CS)
  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure containing pointers to any available
                                                 !! thermodynamic fields. Absent fields have
                                                 !! NULL ptrs.
  type(forcing),           intent(in)    :: fluxes !< A structure containing pointers to any
                                                 !! possible forcing fields. Unused fields have
                                                 !! NULL ptrs.
  type(vertvisc_type),     intent(inout) :: visc !< A structure containing vertical viscosities and
                                                 !! related fields.
  real,                    intent(in)    :: dt   !< Time increment in s.
  type(set_visc_cs),       pointer       :: CS   !< The control structure returned by a previous
                                                 !! call to vertvisc_init.

!   The following subroutine calculates the thickness of the surface boundary
! layer for applying an elevated viscosity.  A bulk Richardson criterion or
! the thickness of the topmost NKML layers (with a bulk mixed layer) are
! currently used.  The thicknesses are given in terms of fractional layers, so
! that this thickness will move as the thickness of the topmost layers change.
!
! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.  In the comments below,
!                the units of h are denoted as H.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (out)     visc - A structure containing vertical viscosities and related
!                   fields.
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 vertvisc_init.

  real, dimension(SZIB_(G)) :: &
    htot, &     !   The total depth of the layers being that are within the
                ! surface mixed layer, in H.
    thtot, &    !   The integrated temperature of layers that are within the
                ! surface mixed layer, in H degC.
    shtot, &    !   The integrated salt of layers that are within the
                ! surface mixed layer H PSU.
    rhtot, &    !   The integrated density of layers that are within the
                ! surface mixed layer,  in H kg m-3.  Rhtot is only used if no
                ! equation of state is used.
    uhtot, &    !   The depth integrated zonal and meridional velocities within
    vhtot, &    ! the surface mixed layer, in H m s-1.
    idecay_len_tke, &  ! The inverse of a turbulence decay length scale, in H-1.
    dr_dt, &    !   Partial derivative of the density at the base of layer nkml
                ! (roughly the base of the mixed layer) with temperature, in
                ! units of kg m-3 K-1.
    dr_ds, &    !   Partial derivative of the density at the base of layer nkml
                ! (roughly the base of the mixed layer) with salinity, in units
                ! of kg m-3 psu-1.
    ustar, &    !   The surface friction velocity under ice shelves, in m s-1.
    press, &    ! The pressure at which dR_dT and dR_dS are evaluated, in Pa.
    t_eos, &    ! T_EOS and S_EOS are the potential temperature and salnity at which dR_dT and dR_dS
    s_eos       ! which dR_dT and dR_dS are evaluated, in degC and PSU.
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    mask_u      ! A mask that disables any contributions from u points that
                ! are land or past open boundary conditions, nondim., 0 or 1.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    mask_v      ! A mask that disables any contributions from v points that
                ! are land or past open boundary conditions, nondim., 0 or 1.
  real :: h_at_vel(szib_(g),szk_(g))! Layer thickness at velocity points,
                ! using an upwind-biased second order accurate estimate based
                ! on the previous velocity direction, in H.
  integer :: k_massive(szib_(g)) ! The k-index of the deepest layer yet found
                ! that has more than h_tiny thickness and will be in the
                ! viscous mixed layer.
  real :: Uh2   ! The squared magnitude of the difference between the velocity
                ! integrated through the mixed layer and the velocity of the
                ! interior layer layer times the depth of the the mixed layer,
                ! in H2 m2 s-2.
  real :: htot_vel  ! Sum of the layer thicknesses up to some
                    ! point, in H (i.e., m or kg m-2).
  real :: hwtot     ! Sum of the thicknesses used to calculate
                    ! the near-bottom velocity magnitude, in H.
  real :: hutot     ! Running sum of thicknesses times the
                    ! velocity magnitudes, in H m s-1.
  real :: hweight   ! The thickness of a layer that is within Hbbl
                    ! of the bottom, in H.

  real :: hlay      ! The layer thickness at velocity points, in H.
  real :: I_2hlay   ! 1 / 2*hlay, in H-1.
  real :: T_lay     ! The layer temperature at velocity points, in deg C.
  real :: S_lay     ! The layer salinity at velocity points, in PSU.
  real :: Rlay      ! The layer potential density at velocity points, in kg m-3.
  real :: Rlb       ! The potential density of the layer below, in kg m-3.
  real :: v_at_u    ! The meridonal velocity at a zonal velocity point in m s-1.
  real :: u_at_v    ! The zonal velocity at a meridonal velocity point in m s-1.
  real :: gHprime   ! The mixed-layer internal gravity wave speed squared, based
                    ! on the mixed layer thickness and density difference across
                    ! the base of the mixed layer, in m2 s-2.
  real :: RiBulk    ! The bulk Richardson number below which water is in the
                    ! viscous mixed layer, including reduction for turbulent
                    ! decay. Nondimensional.
  real :: dt_Rho0   ! The time step divided by the conversion from the layer
                    ! thickness to layer mass, in s H m2 kg-1.
  real :: g_H_Rho0  !   The gravitational acceleration times the conversion from
                    ! H to m divided by the mean density, in m5 s-2 H-1 kg-1.
  real :: ustarsq     ! 400 times the square of ustar, times
                      ! Rho0 divided by G_Earth and the conversion
                      ! from m to thickness units, in kg m-2 or kg2 m-5.
  real :: cdrag_sqrt  ! Square root of the drag coefficient, nd.
  real :: oldfn       ! The integrated energy required to
                      ! entrain up to the bottom of the layer,
                      ! divided by G_Earth, in H kg m-3.
  real :: Dfn         ! The increment in oldfn for entraining
                      ! the layer, in H kg m-3.
  real :: Dh          ! The increment in layer thickness from
                      ! the present layer, in H.
  real :: U_bg_sq   ! The square of an assumed background velocity, for
                    ! calculating the mean magnitude near the top for use in
                    ! the quadratic surface drag, in m2.
  real :: h_tiny    ! A very small thickness, in H. Layers that are less than
                    ! h_tiny can not be the deepest in the viscous mixed layer.
  real :: absf      ! The absolute value of f averaged to velocity points, s-1.
  real :: U_star    ! The friction velocity at velocity points, in m s-1.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in H.
  real :: Rho0x400_G  ! 400*Rho0/G_Earth, in kg s2 m-4.  The 400 is a
                      ! constant proposed by Killworth and Edwards, 1999.
  real :: H_to_m, m_to_H   ! Local copies of unit conversion factors.
  real :: ustar1    ! ustar in units of H/s
  real :: h2f2      ! (h*2*f)^2
  logical :: use_EOS, do_any, do_any_shelf, do_i(szib_(g))
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, K2, nkmb, nkml, n
  type(ocean_obc_type), pointer :: OBC => null()

  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
  nkmb = gv%nk_rho_varies ; nkml = gv%nkml

  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc_ML): "//&
         "Module must be initialized before it is used.")
  if (.not.(cs%dynamic_viscous_ML .or. associated(fluxes%frac_shelf_h))) return

  rho0x400_g = 400.0*(gv%Rho0/gv%g_Earth)*gv%m_to_H
  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
  cdrag_sqrt=sqrt(cs%cdrag)

  obc => cs%OBC
  use_eos = associated(tv%eqn_of_state)
  dt_rho0 = dt/gv%H_to_kg_m2
  h_neglect = gv%H_subroundoff
  h_tiny = 2.0*gv%Angstrom + h_neglect
  g_h_rho0 = (gv%g_Earth * gv%H_to_m) / gv%Rho0
  h_to_m = gv%H_to_m ; m_to_h = gv%m_to_H

  if (associated(fluxes%frac_shelf_h)) then
    ! This configuration has ice shelves, and the appropriate variables need to
    ! be allocated.
    if (.not.associated(fluxes%frac_shelf_u)) call mom_error(fatal, &
      "set_viscous_ML: fluxes%frac_shelf_h is associated, but " // &
      "fluxes%frac_shelf_u is not.")
    if (.not.associated(fluxes%frac_shelf_v)) call mom_error(fatal, &
      "set_viscous_ML: fluxes%frac_shelf_h is associated, but " // &
      "fluxes%frac_shelf_v is not.")
    if (.not.associated(visc%taux_shelf)) then
      allocate(visc%taux_shelf(g%IsdB:g%IedB,g%jsd:g%jed))
      visc%taux_shelf(:,:) = 0.0
    endif
    if (.not.associated(visc%tauy_shelf)) then
      allocate(visc%tauy_shelf(g%isd:g%ied,g%JsdB:g%JedB))
      visc%tauy_shelf(:,:) = 0.0
    endif
    if (.not.associated(visc%tbl_thick_shelf_u)) then
      allocate(visc%tbl_thick_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed))
      visc%tbl_thick_shelf_u(:,:) = 0.0
    endif
    if (.not.associated(visc%tbl_thick_shelf_v)) then
      allocate(visc%tbl_thick_shelf_v(g%isd:g%ied,g%JsdB:g%JedB))
      visc%tbl_thick_shelf_v(:,:) = 0.0
    endif
    if (.not.associated(visc%kv_tbl_shelf_u)) then
      allocate(visc%kv_tbl_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed))
      visc%kv_tbl_shelf_u(:,:) = 0.0
    endif
    if (.not.associated(visc%kv_tbl_shelf_v)) then
      allocate(visc%kv_tbl_shelf_v(g%isd:g%ied,g%JsdB:g%JedB))
      visc%kv_tbl_shelf_v(:,:) = 0.0
    endif

    !  With a linear drag law, the friction velocity is already known.
!    if (CS%linear_drag) ustar(:) = cdrag_sqrt*CS%drag_bg_vel
  endif

  !$OMP parallel do default(shared)
  do j=js-1,je ; do i=is-1,ie+1
    mask_v(i,j) = g%mask2dCv(i,j)
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je+1 ; do i=is-1,ie
    mask_u(i,j) = g%mask2dCu(i,j)
  enddo ; enddo

  if (associated(obc)) then ; do n=1,obc%number_of_segments
    ! Now project bottom depths across cell-corner points in the OBCs.  The two
    ! projections have to occur in sequence and can not be combined easily.
    if (.not. obc%segment(n)%on_pe) cycle
    ! Use a one-sided projection of bottom depths at OBC points.
    i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
    if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
      do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
        if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
        if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
      enddo
    elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= je)) then
      do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
        if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
        if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
      enddo
    endif
  enddo ; endif

!$OMP parallel do default(none) shared(u, v, h, tv, fluxes, visc, dt, G, GV, CS, use_EOS, &
!$OMP                                  dt_Rho0, h_neglect, h_tiny, g_H_Rho0,js,je,        &
!$OMP                                  H_to_m, m_to_H, Isq, Ieq, nz, U_bg_sq,mask_v,      &
!$OMP                                  cdrag_sqrt,Rho0x400_G,nkml) &
!$OMP                          private(do_any,htot,do_i,k_massive,Thtot,uhtot,vhtot,U_Star, &
!$OMP                                  Idecay_len_TKE,press,k2,I_2hlay,T_EOS,S_EOS,dR_dT,   &
!$OMP                                  dR_dS,hlay,v_at_u,Uh2,T_lay,S_lay,gHprime,           &
!$OMP                                  RiBulk,Shtot,Rhtot,absf,do_any_shelf,                &
!$OMP                                  h_at_vel,ustar,htot_vel,hwtot,hutot,hweight,ustarsq, &
!$OMP                                  oldfn,Dfn,Dh,Rlay,Rlb,h2f2,ustar1)
  do j=js,je  ! u-point loop
    if (cs%dynamic_viscous_ML) then
      do_any = .false.
      do i=isq,ieq
        htot(i) = 0.0
        if (g%mask2dCu(i,j) < 0.5) then
          do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
        else
          do_i(i) = .true. ; do_any = .true.
          k_massive(i) = nkml
          thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
          uhtot(i) = dt_rho0 * fluxes%taux(i,j)
          vhtot(i) = 0.25 * dt_rho0 * ((fluxes%tauy(i,j) + fluxes%tauy(i+1,j-1)) + &
                                       (fluxes%tauy(i,j-1) + fluxes%tauy(i+1,j)))

          if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
            absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
            if (cs%omega_frac > 0.0) &
              absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
          endif
          u_star = max(cs%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i+1,j)))
          idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * h_to_m
        endif
      enddo

      if (do_any) then ; do k=1,nz
        if (k > nkml) then
          do_any = .false.
          if (use_eos .and. (k==nkml+1)) then
            ! Find dRho/dT and dRho_dS.
            do i=isq,ieq
              press(i) = gv%H_to_Pa * htot(i)
              k2 = max(1,nkml)
              i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
              t_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i+1,j,k2)*tv%T(i+1,j,k2)) * i_2hlay
              s_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i+1,j,k2)*tv%S(i+1,j,k2)) * i_2hlay
            enddo
            call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
                                          isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
          endif

          do i=isq,ieq ; if (do_i(i)) then

            hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
            if (hlay > h_tiny) then ! Only consider non-vanished layers.
              i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
              v_at_u = 0.5 * (h(i,j,k)   * (v(i,j,k) + v(i,j-1,k)) + &
                              h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i_2hlay
              uh2 = (uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2

              if (use_eos) then
                t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i_2hlay
                s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i_2hlay
                ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
                                      dr_ds(i) * (s_lay*htot(i) - shtot(i)))
              else
                ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
              endif

              if (ghprime > 0.0) then
                ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
                if (ribulk * uh2 <= (htot(i)**2) * ghprime) then
                  visc%nkml_visc_u(i,j) = real(k_massive(i))
                  do_i(i) = .false.
                elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
                  visc%nkml_visc_u(i,j) = real(k-1) + &
                    ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
                  do_i(i) = .false.
                endif
              endif
              k_massive(i) = k
            endif ! hlay > h_tiny

            if (do_i(i)) do_any = .true.
          endif ; enddo

          if (.not.do_any) exit ! All columns are done.
        endif

        do i=isq,ieq ; if (do_i(i)) then
          htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
          uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
          vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
                                        h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
          if (use_eos) then
            thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k))
            shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k))
          else
            rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
          endif
        endif ; enddo
      enddo ; endif

      if (do_any) then ; do i=isq,ieq ; if (do_i(i)) then
        visc%nkml_visc_u(i,j) = k_massive(i)
      endif ; enddo ; endif
    endif ! dynamic_viscous_ML

    do_any_shelf = .false.
    if (associated(fluxes%frac_shelf_h)) then
      do i=isq,ieq
        if (fluxes%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0) then
          do_i(i) = .false.
          visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
        else
          do_i(i) = .true. ; do_any_shelf = .true.
        endif
      enddo
    endif

    if (do_any_shelf) then
      do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
        if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
          h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
                          (h(i,j,k) + h(i+1,j,k) + h_neglect)
        else
          h_at_vel(i,k) =  0.5 * (h(i,j,k) + h(i+1,j,k))
        endif
      else
        h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
      endif ; enddo ; enddo

      do i=isq,ieq ; if (do_i(i)) then
        htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
        thtot(i) = 0.0 ; shtot(i) = 0.0
        if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
          if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
          hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
          if (hweight <= 1.5*gv%Angstrom + h_neglect) cycle

          htot_vel  = htot_vel + h_at_vel(i,k)
          hwtot = hwtot + hweight

          if (.not.cs%linear_drag) then
            v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v)
            hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
                                           v_at_u**2 + u_bg_sq)
          endif
          if (use_eos) then
            thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
            shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
          endif
        enddo ; endif

        if ((.not.cs%linear_drag) .and. (hwtot > 0.0)) then
          ustar(i) = cdrag_sqrt*hutot/hwtot
        else
          ustar(i) = cdrag_sqrt*cs%drag_bg_vel
        endif

        if (use_eos) then ; if (hwtot > 0.0) then
          t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
        else
          t_eos(i) = 0.0 ; s_eos(i) = 0.0
        endif ; endif
      endif ; enddo ! I-loop

      if (use_eos) then
        call calculate_density_derivs(t_eos, s_eos, fluxes%p_surf(:,j), &
                     dr_dt, dr_ds, isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
      endif

      do i=isq,ieq ; if (do_i(i)) then
  !  The 400.0 in this expression is the square of a constant proposed
  !  by Killworth and Edwards, 1999, in equation (2.20).
        ustarsq = rho0x400_g * ustar(i)**2
        htot(i) = 0.0
        if (use_eos) then
          thtot(i) = 0.0 ; shtot(i) = 0.0
          do k=1,nz-1
            if (h_at_vel(i,k) <= 0.0) cycle
            t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
            s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
            oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
            if (oldfn >= ustarsq) exit

            dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
                   dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
                  (h_at_vel(i,k)+htot(i))
            if ((oldfn + dfn) <= ustarsq) then
              dh = h_at_vel(i,k)
            else
              dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
            endif

            htot(i) = htot(i) + dh
            thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
          enddo
          if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
            t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
            s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
            if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
                dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
              htot(i) = htot(i) + h_at_vel(i,nz)
          endif ! Examination of layer nz.
        else  ! Use Rlay as the density variable.
          rhtot = 0.0
          do k=1,nz-1
            rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)

            oldfn = rlay*htot(i) - rhtot(i)
            if (oldfn >= ustarsq) exit

            dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
            if ((oldfn + dfn) <= ustarsq) then
              dh = h_at_vel(i,k)
            else
              dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
            endif

            htot(i) = htot(i) + dh
            rhtot(i) = rhtot(i) + rlay*dh
          enddo
          if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
            htot(i) = htot(i) + h_at_vel(i,nz)
        endif ! use_EOS

       !visc%tbl_thick_shelf_u(I,j) = max(CS%Htbl_shelf_min, &
       !    htot(I) / (0.5 + sqrt(0.25 + &
       !                 (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &
       !                 (ustar(i)*m_to_H)**2 )) )
        ustar1 = ustar(i)*m_to_h
        h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%Omega)**2
        visc%tbl_thick_shelf_u(i,j) = max(cs%Htbl_shelf_min, &
            ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
        visc%kv_tbl_shelf_u(i,j) = max(cs%KV_TBL_min, &
                       cdrag_sqrt*ustar(i)*visc%tbl_thick_shelf_u(i,j))
      endif ; enddo ! I-loop
    endif ! do_any_shelf

  enddo ! j-loop at u-points

!$OMP parallel do default(none) shared(u, v, h, tv, fluxes, visc, dt, G, GV, CS, use_EOS,&
!$OMP                                  dt_Rho0, h_neglect, h_tiny, g_H_Rho0,is,ie,       &
!$OMP                                  Jsq,Jeq,nz,U_bg_sq,cdrag_sqrt,Rho0x400_G,nkml,    &
!$OMP                                  m_to_H,H_to_m,mask_u) &
!$OMP                          private(do_any,htot,do_i,k_massive,Thtot,vhtot,uhtot,absf,&
!$OMP                                  U_Star,Idecay_len_TKE,press,k2,I_2hlay,T_EOS,     &
!$OMP                                  S_EOS,dR_dT, dR_dS,hlay,u_at_v,Uh2,               &
!$OMP                                  T_lay,S_lay,gHprime,RiBulk,do_any_shelf,          &
!$OMP                                  Shtot,Rhtot,ustar,h_at_vel,htot_vel,hwtot,        &
!$OMP                                  hutot,hweight,ustarsq,oldfn,Dh,Rlay,Rlb,Dfn,      &
!$OMP                                  h2f2,ustar1)
  do j=jsq,jeq  ! v-point loop
    if (cs%dynamic_viscous_ML) then
      do_any = .false.
      do i=is,ie
        htot(i) = 0.0
        if (g%mask2dCv(i,j) < 0.5) then
          do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
        else
          do_i(i) = .true. ; do_any = .true.
          k_massive(i) = nkml
          thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
          vhtot(i) = dt_rho0 * fluxes%tauy(i,j)
          uhtot(i) = 0.25 * dt_rho0 * ((fluxes%taux(i,j) + fluxes%taux(i-1,j+1)) + &
                                       (fluxes%taux(i-1,j) + fluxes%taux(i,j+1)))

         if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
           absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
           if (cs%omega_frac > 0.0) &
             absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
         endif

         u_star = max(cs%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i,j+1)))
         idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * h_to_m

        endif
      enddo

      if (do_any) then ; do k=1,nz
        if (k > nkml) then
          do_any = .false.
          if (use_eos .and. (k==nkml+1)) then
            ! Find dRho/dT and dRho_dS.
            do i=is,ie
              press(i) = gv%H_to_Pa * htot(i)
              k2 = max(1,nkml)
              i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
              t_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i,j+1,k2)*tv%T(i,j+1,k2)) * i_2hlay
              s_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i,j+1,k2)*tv%S(i,j+1,k2)) * i_2hlay
            enddo
            call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
                                          is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
          endif

          do i=is,ie ; if (do_i(i)) then

            hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
            if (hlay > h_tiny) then ! Only consider non-vanished layers.
              i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
              u_at_v = 0.5 * (h(i,j,k)   * (u(i-1,j,k)   + u(i,j,k)) + &
                              h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i_2hlay
              uh2 = (uhtot(i) - htot(i)*u_at_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2

              if (use_eos) then
                t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i_2hlay
                s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i_2hlay
                ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
                                      dr_ds(i) * (s_lay*htot(i) - shtot(i)))
              else
                ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
              endif

              if (ghprime > 0.0) then
                ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
                if (ribulk * uh2 <= htot(i)**2 * ghprime) then
                  visc%nkml_visc_v(i,j) = real(k_massive(i))
                  do_i(i) = .false.
                elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
                  visc%nkml_visc_v(i,j) = real(k-1) + &
                    ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
                  do_i(i) = .false.
                endif
              endif
              k_massive(i) = k
            endif ! hlay > h_tiny

            if (do_i(i)) do_any = .true.
          endif ; enddo

          if (.not.do_any) exit ! All columns are done.
        endif

        do i=is,ie ; if (do_i(i)) then
          htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
          vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
          uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
                                        h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
          if (use_eos) then
            thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k))
            shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k))
          else
            rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
          endif
        endif ; enddo
      enddo ; endif

      if (do_any) then ; do i=is,ie ; if (do_i(i)) then
        visc%nkml_visc_v(i,j) = k_massive(i)
      endif ; enddo ; endif
    endif ! dynamic_viscous_ML

    do_any_shelf = .false.
    if (associated(fluxes%frac_shelf_h)) then
      do i=is,ie
        if (fluxes%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0) then
          do_i(i) = .false.
          visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
        else
          do_i(i) = .true. ; do_any_shelf = .true.
        endif
      enddo
    endif

    if (do_any_shelf) then
      do k=1,nz ; do i=is,ie ; if (do_i(i)) then
        if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
          h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
                          (h(i,j,k) + h(i,j+1,k) + h_neglect)
        else
          h_at_vel(i,k) =  0.5 * (h(i,j,k) + h(i,j+1,k))
        endif
      else
        h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
      endif ; enddo ; enddo

      do i=is,ie ; if (do_i(i)) then
        htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
        thtot(i) = 0.0 ; shtot(i) = 0.0
        if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
          if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
          hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
          if (hweight <= 1.5*gv%Angstrom + h_neglect) cycle

          htot_vel  = htot_vel + h_at_vel(i,k)
          hwtot = hwtot + hweight

          if (.not.cs%linear_drag) then
            u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u)
            hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
                                           u_at_v**2 + u_bg_sq)
          endif
          if (use_eos) then
            thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
            shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
          endif
        enddo ; endif

        if (.not.cs%linear_drag) then ; if (hwtot > 0.0) then
          ustar(i) = cdrag_sqrt*hutot/hwtot
        else
          ustar(i) = cdrag_sqrt*cs%drag_bg_vel
        endif ; endif

        if (use_eos) then ; if (hwtot > 0.0) then
          t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
        else
          t_eos(i) = 0.0 ; s_eos(i) = 0.0
        endif ; endif
      endif ; enddo ! I-loop

      if (use_eos) then
        call calculate_density_derivs(t_eos, s_eos, fluxes%p_surf(:,j), &
                     dr_dt, dr_ds, is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
      endif

      do i=is,ie ; if (do_i(i)) then
  !  The 400.0 in this expression is the square of a constant proposed
  !  by Killworth and Edwards, 1999, in equation (2.20).
        ustarsq = rho0x400_g * ustar(i)**2
        htot(i) = 0.0
        if (use_eos) then
          thtot(i) = 0.0 ; shtot(i) = 0.0
          do k=1,nz-1
            if (h_at_vel(i,k) <= 0.0) cycle
            t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
            s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
            oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
            if (oldfn >= ustarsq) exit

            dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
                   dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
                  (h_at_vel(i,k)+htot(i))
            if ((oldfn + dfn) <= ustarsq) then
              dh = h_at_vel(i,k)
            else
              dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
            endif

            htot(i) = htot(i) + dh
            thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
          enddo
          if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
            t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
            s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
            if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
                dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
              htot(i) = htot(i) + h_at_vel(i,nz)
          endif ! Examination of layer nz.
        else  ! Use Rlay as the density variable.
          rhtot = 0.0
          do k=1,nz-1
            rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)

            oldfn = rlay*htot(i) - rhtot(i)
            if (oldfn >= ustarsq) exit

            dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
            if ((oldfn + dfn) <= ustarsq) then
              dh = h_at_vel(i,k)
            else
              dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
            endif

            htot(i) = htot(i) + dh
            rhtot = rhtot + rlay*dh
          enddo
          if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
            htot(i) = htot(i) + h_at_vel(i,nz)
        endif ! use_EOS

       !visc%tbl_thick_shelf_v(i,J) = max(CS%Htbl_shelf_min, &
       !    htot(i) / (0.5 + sqrt(0.25 + &
       !        (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &
       !        (ustar(i)*m_to_H)**2 )) )
        ustar1 = ustar(i)*m_to_h
        h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%Omega)**2
        visc%tbl_thick_shelf_v(i,j) = max(cs%Htbl_shelf_min, &
            ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
        visc%kv_tbl_shelf_v(i,j) = max(cs%KV_TBL_min, &
                       cdrag_sqrt*ustar(i)*visc%tbl_thick_shelf_v(i,j))
      endif ; enddo ! i-loop
    endif ! do_any_shelf

  enddo ! J-loop at v-points

  if (cs%debug) then
    if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) &
      call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, &
                    visc%nkml_visc_v, g%HI,haloshift=0)
  endif
  if (cs%id_nkml_visc_u > 0) &
    call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
  if (cs%id_nkml_visc_v > 0) &
    call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)

end subroutine set_viscous_ml

!>   This subroutine is used to register any fields associated with the
!! vertvisc_type.
subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
  type(hor_index_type),    intent(in)    :: HI         !< A horizontal index type structure.
  type(verticalgrid_type), intent(in)    :: GV         !< The ocean's vertical grid structure.
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                       !! parameters.
  type(vertvisc_type),     intent(inout) :: visc       !< A structure containing vertical
                                                       !! viscosities and related fields.
                                                       !! Allocated here.
  type(mom_restart_cs),    pointer       :: restart_CS !< A pointer to the restart control
                                                       !! structure.
!   This subroutine is used to register any fields associated with the
! vertvisc_type.
! Arguments: HI - A horizontal index type structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (out)     visc - A structure containing vertical viscosities and related
!                   fields.  Allocated here.
!  (in)      restart_CS - A pointer to the restart control structure.
  type(vardesc) :: vd
  logical :: use_kappa_shear, adiabatic, useKPP, useEPBL
  logical :: use_CVMix, MLE_use_PBL_MLD
  integer :: isd, ied, jsd, jed, nz
  character(len=40)  :: mdl = "MOM_set_visc"  ! This module's name.
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
                 do_not_log=.true.)
  use_kappa_shear = .false. ; use_cvmix = .false. ;
  usekpp = .false. ; useepbl = .false.
  if (.not.adiabatic) then
    use_kappa_shear = kappa_shear_is_used(param_file)
    use_cvmix = cvmix_shear_is_used(param_file)
    call get_param(param_file, mdl, "USE_KPP", usekpp, &
                 "If true, turns on the [CVmix] KPP scheme of Large et al., 1984,\n"// &
                 "to calculate diffusivities and non-local transport in the OBL.", &
                 default=.false., do_not_log=.true.)
    call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, &
                 "If true, use an implied energetics planetary boundary \n"//&
                 "layer scheme to determine the diffusivity and viscosity \n"//&
                 "in the surface boundary layer.", default=.false., do_not_log=.true.)
  endif

  if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix) then
    allocate(visc%Kd_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kd_turb(:,:,:) = 0.0
    allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0
    allocate(visc%Kv_turb(isd:ied,jsd:jed,nz+1)) ; visc%Kv_turb(:,:,:) = 0.0

    vd = var_desc("Kd_turb","m2 s-1","Turbulent diffusivity at interfaces", &
                  hor_grid='h', z_grid='i')
    call register_restart_field(visc%Kd_turb, vd, .false., restart_cs)

    vd = var_desc("TKE_turb","m2 s-2","Turbulent kinetic energy per unit mass at interfaces", &
                  hor_grid='h', z_grid='i')
    call register_restart_field(visc%TKE_turb, vd, .false., restart_cs)
    vd = var_desc("Kv_turb","m2 s-1","Turbulent viscosity at interfaces", &
                  hor_grid='h', z_grid='i')
    call register_restart_field(visc%Kv_turb, vd, .false., restart_cs)
  endif

  ! visc%MLD is used to communicate the state of the (e)PBL to the rest of the model
  call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
                 default=.false., do_not_log=.true.)
  if (mle_use_pbl_mld) then
    allocate(visc%MLD(isd:ied,jsd:jed)) ; visc%MLD(:,:) = 0.0
    vd = var_desc("MLD","m","Instantaneous active mixing layer depth", &
                  hor_grid='h', z_grid='1')
    call register_restart_field(visc%MLD, vd, .false., restart_cs)
  endif

end subroutine set_visc_register_restarts

subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC)
  type(time_type), target, intent(in)    :: Time !< The current model time.
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(vertvisc_type),     intent(inout) :: visc !< A structure containing vertical viscosities and
                                                 !! related fields.  Allocated here.
  type(set_visc_cs),       pointer       :: CS   !< A pointer that is set to point to the control
                                                 !! structure for this module
  type(ocean_obc_type),    pointer       :: OBC
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (out)     visc - A structure containing vertical viscosities and related
!                   fields.  Allocated here.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
  real    :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt
  real    :: Kv_background
  real    :: omega_frac_dflt
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, i, j, n
  logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega
  type(obc_segment_type), pointer :: segment  ! pointer to OBC segment type
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_set_visc"  ! This module's name.

  if (associated(cs)) then
    call mom_error(warning, "set_visc_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  cs%OBC => obc

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB


  cs%diag => diag

! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")
  cs%RiNo_mix = .false.
  use_kappa_shear = .false. ; differential_diffusion = .false. !; adiabatic = .false.  ! Needed? -AJA
  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
                 "If true, the bottom stress is calculated with a drag \n"//&
                 "law of the form c_drag*|u|*u. The velocity magnitude \n"//&
                 "may be an assumed value or it may be based on the \n"//&
                 "actual velocity in the bottommost HBBL, depending on \n"//&
                 "LINEAR_DRAG.", default=.true.)
  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
                 "If true, the bottom drag is exerted directly on each \n"//&
                 "layer proportional to the fraction of the bottom it \n"//&
                 "overlies.", default=.false.)
  call get_param(param_file, mdl, "LINEAR_DRAG", cs%linear_drag, &
                 "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag \n"//&
                 "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
                 do_not_log=.true.)
  if (adiabatic) then
    call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
                 "There are no diapycnal mass fluxes if ADIABATIC is \n"//&
                 "true. This assumes that KD = KDML = 0.0 and that \n"//&
                 "there is no buoyancy forcing, but makes the model \n"//&
                 "faster by eliminating subroutine calls.", default=.false.)
  endif

  if (.not.adiabatic) then
    use_kappa_shear = kappa_shear_is_used(param_file)
    cs%RiNo_mix = use_kappa_shear
    call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, &
                 "If true, increase diffusivitives for temperature or salt \n"//&
                 "based on double-diffusive paramaterization from MOM4/KPP.", &
                 default=.false.)
  endif
  call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, &
                 "The turbulent Prandtl number applied to shear \n"//&
                 "instability.", units="nondim", default=1.0)
  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)

  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
                 "If true, use a bulk Richardson number criterion to \n"//&
                 "determine the mixed layer thickness for viscosity.", &
                 default=.false.)
  if (cs%dynamic_viscous_ML) then
    call get_param(param_file, mdl, "BULK_RI_ML", bulk_ri_ml_dflt, default=0.0)
    call get_param(param_file, mdl, "BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
                 "The efficiency with which mean kinetic energy released \n"//&
                 "by mechanically forced entrainment of the mixed layer \n"//&
                 "is converted to turbulent kinetic energy.  By default, \n"//&
                 "BULK_RI_ML_VISC = BULK_RI_ML or 0.", units="nondim", &
                 default=bulk_ri_ml_dflt)
    call get_param(param_file, mdl, "TKE_DECAY", tke_decay_dflt, default=0.0)
    call get_param(param_file, mdl, "TKE_DECAY_VISC", cs%TKE_decay, &
                 "TKE_DECAY_VISC relates the vertical rate of decay of \n"//&
                 "the TKE available for mechanical entrainment to the \n"//&
                 "natural Ekman depth for use in calculating the dynamic \n"//&
                 "mixed layer viscosity.  By default, \n"//&
                 "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
                 default=tke_decay_dflt)
    call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
                 "If true, use the absolute rotation rate instead of the \n"//&
                 "vertical component of rotation when setting the decay \n"//&
                   "scale for turbulence.", default=.false., do_not_log=.true.)
    omega_frac_dflt = 0.0
    if (use_omega) then
      call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
      omega_frac_dflt = 1.0
    endif
    call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
                   "When setting the decay scale for turbulence, use this \n"//&
                   "fraction of the absolute rotation rate blended with the \n"//&
                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
                   units="nondim", default=omega_frac_dflt)
    call get_param(param_file, mdl, "OMEGA", cs%omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)
    ! This give a minimum decay scale that is typically much less than Angstrom.
    cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_z + gv%H_to_m*gv%H_subroundoff)
  else
    call get_param(param_file, mdl, "OMEGA", cs%omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)
  endif

  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
                 "The thickness of a bottom boundary layer with a \n"//&
                 "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or \n"//&
                 "the thickness over which near-bottom velocities are \n"//&
                 "averaged for the drag law if BOTTOMDRAGLAW is defined \n"//&
                 "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.)
  if (cs%bottomdraglaw) then
    call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
                 "CDRAG is the drag coefficient relating the magnitude of \n"//&
                 "the velocity field to the bottom stress. CDRAG is only \n"//&
                 "used if BOTTOMDRAGLAW is defined.", units="nondim", &
                 default=0.003)
    call get_param(param_file, mdl, "DRAG_BG_VEL", cs%drag_bg_vel, &
                 "DRAG_BG_VEL is either the assumed bottom velocity (with \n"//&
                 "LINEAR_DRAG) or an unresolved  velocity that is \n"//&
                 "combined with the resolved velocity to estimate the \n"//&
                 "velocity magnitude.  DRAG_BG_VEL is only used when \n"//&
                 "BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0)
    call get_param(param_file, mdl, "BBL_USE_EOS", cs%BBL_use_EOS, &
                 "If true, use the equation of state in determining the \n"//&
                 "properties of the bottom boundary layer.  Otherwise use \n"//&
                 "the layer target potential densities.", default=.false.)
  endif
  call get_param(param_file, mdl, "BBL_THICK_MIN", cs%BBL_thick_min, &
                 "The minimum bottom boundary layer thickness that can be \n"//&
                 "used with BOTTOMDRAGLAW. This might be \n"//&
                 "Kv / (cdrag * drag_bg_vel) to give Kv as the minimum \n"//&
                 "near-bottom viscosity.", units="m", default=0.0)
  call get_param(param_file, mdl, "HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
                 "The minimum top boundary layer thickness that can be \n"//&
                 "used with BOTTOMDRAGLAW. This might be \n"//&
                 "Kv / (cdrag * drag_bg_vel) to give Kv as the minimum \n"//&
                 "near-top viscosity.", units="m", default=cs%BBL_thick_min)
  call get_param(param_file, mdl, "HTBL_SHELF", cs%Htbl_shelf, &
                 "The thickness over which near-surface velocities are \n"//&
                 "averaged for the drag law under an ice shelf.  By \n"//&
                 "default this is the same as HBBL", units="m", default=cs%Hbbl)

  call get_param(param_file, mdl, "KV", kv_background, &
                 "The background kinematic viscosity in the interior. \n"//&
                 "The molecular value, ~1e-6 m2 s-1, may be used.", &
                 units="m2 s-1", fail_if_missing=.true.)
  call get_param(param_file, mdl, "KV_BBL_MIN", cs%KV_BBL_min, &
                 "The minimum viscosities in the bottom boundary layer.", &
                 units="m2 s-1", default=kv_background)
  call get_param(param_file, mdl, "KV_TBL_MIN", cs%KV_TBL_min, &
                 "The minimum viscosities in the top boundary layer.", &
                 units="m2 s-1", default=kv_background)

  if (cs%Channel_drag) then
    call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, default=-1.0)

    csmag_chan_dflt = 0.15
    if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1

    call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", cs%c_Smag, &
                 "The nondimensional Laplacian Smagorinsky constant used \n"//&
                 "in calculating the channel drag if it is enabled.  The \n"//&
                 "default is to use the same value as SMAG_LAP_CONST if \n"//&
                 "it is defined, or 0.15 if it is not. The value used is \n"//&
                 "also 0.15 if the specified value is negative.", &
                 units="nondim", default=csmag_chan_dflt)
    if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
  endif

  if (cs%bottomdraglaw) then
    allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed)) ; visc%bbl_thick_u = 0.0
    allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed)) ; visc%kv_bbl_u = 0.0
    allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb)) ; visc%bbl_thick_v = 0.0
    allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb)) ; visc%kv_bbl_v = 0.0
    allocate(visc%ustar_bbl(isd:ied,jsd:jed)) ; visc%ustar_bbl = 0.0
    allocate(visc%TKE_bbl(isd:ied,jsd:jed)) ; visc%TKE_bbl = 0.0

    cs%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
       diag%axesCu1, time, 'BBL thickness at u points', 'meter')
    cs%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
       time, 'BBL viscosity at u points', 'meter2 second-1')
    cs%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
       diag%axesCv1, time, 'BBL thickness at v points', 'meter')
    cs%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
       time, 'BBL viscosity at v points', 'meter2 second-1')
  endif
  if (cs%Channel_drag) then
    allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray_u = 0.0
    allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray_v = 0.0
    cs%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
       time, 'Rayleigh drag velocity at u points', 'meter second-1')
    cs%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
       time, 'Rayleigh drag velocity at v points', 'meter second-1')
  endif

  if (differential_diffusion) then
    allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0
    allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0
  endif

  if (cs%dynamic_viscous_ML) then
    allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed)) ; visc%nkml_visc_u = 0.0
    allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb)) ; visc%nkml_visc_v = 0.0
    cs%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
       diag%axesCu1, time, 'Number of layers in viscous mixed layer at u points', 'meter')
    cs%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
       diag%axesCv1, time, 'Number of layers in viscous mixed layer at v points', 'meter')
  endif

  cs%Hbbl = cs%Hbbl * gv%m_to_H
  cs%BBL_thick_min = cs%BBL_thick_min * gv%m_to_H

end subroutine set_visc_init

subroutine set_visc_end(visc, CS)
  type(vertvisc_type), intent(inout) :: visc
  type(set_visc_cs),   pointer       :: CS
  if (cs%bottomdraglaw) then
    deallocate(visc%bbl_thick_u) ; deallocate(visc%bbl_thick_v)
    deallocate(visc%kv_bbl_u) ; deallocate(visc%kv_bbl_v)
  endif
  if (cs%Channel_drag) then
    deallocate(visc%Ray_u) ; deallocate(visc%Ray_v)
  endif
  if (cs%dynamic_viscous_ML) then
    deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v)
  endif
  if (associated(visc%Kd_turb)) deallocate(visc%Kd_turb)
  if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
  if (associated(visc%Kv_turb)) deallocate(visc%Kv_turb)
  if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
  if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl)
  if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf)
  if (associated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
  if (associated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
  if (associated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
  if (associated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
  if (associated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)

  deallocate(cs)
end subroutine set_visc_end

end module mom_set_visc