MOM_kappa_shear.F90

!> Shear-dependent mixing following Jackson et al. 2008.
module mom_kappa_shear
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_module_driver, clock_module, 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_debugging, only : hchksum, bchksum
use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs
 
implicit none ; private
 
#include <MOM_memory.h>
#ifdef use_netCDF
#include <netcdf.inc>
#endif
 
public calculate_kappa_shear, calc_kappa_shear_vertex, kappa_shear_init
public kappa_shear_is_used, kappa_shear_at_vertex
 
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
 
!> This control structure holds the parameters that regulate shear mixing
type, public :: kappa_shear_cs ; private
  real    :: rino_crit       !< The critical shear Richardson number for
                             !! shear-entrainment. The theoretical value is 0.25.
                             !! The values found by Jackson et al. are 0.25-0.35.
  real    :: shearmix_rate   !< A nondimensional rate scale for shear-driven
                             !! entrainment.  The value given by Jackson et al.
                             !! is 0.085-0.089.
  real    :: fri_curvature   !<   A constant giving the curvature of the function
                             !! of the Richardson number that relates shear to
                             !! sources in the kappa equation [nondim].
                             !! The values found by Jackson et al. are -0.97 - -0.89.
  real    :: c_n             !<   The coefficient for the decay of TKE due to
                             !! stratification (i.e. proportional to N*tke) [nondim].
                             !! The values found by Jackson et al. are 0.24-0.28.
  real    :: c_s             !<   The coefficient for the decay of TKE due to
                             !! shear (i.e. proportional to |S|*tke) [nondim].
                             !! The values found by Jackson et al. are 0.14-0.12.
  real    :: lambda          !<   The coefficient for the buoyancy length scale
                             !! in the kappa equation.  Nondimensional.
                             !! The values found by Jackson et al. are 0.82-0.81.
  real    :: lambda2_n_s     !<   The square of the ratio of the coefficients of
                             !! the buoyancy and shear scales in the diffusivity
                             !! equation, 0 to eliminate the shear scale. Nondim.
  real    :: tke_bg          !<   The background level of TKE [Z2 T-2 ~> m2 s-2].
  real    :: kappa_0         !<   The background diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
  real    :: kappa_tol_err   !<   The fractional error in kappa that is tolerated.
  real    :: prandtl_turb    !< Prandtl number used to convert Kd_shear into viscosity.
  integer :: nkml            !<   The number of layers in the mixed layer, as
                             !! treated in this routine.  If the pieces of the
                             !! mixed layer are not to be treated collectively,
                             !! nkml is set to 1.
  integer :: max_rino_it     !< The maximum number of iterations that may be used
                             !! to estimate the instantaneous shear-driven mixing.
  integer :: max_ks_it       !< The maximum number of iterations that may be used
                             !! to estimate the time-averaged diffusivity.
  logical :: ks_at_vertex    !< If true, do the calculations of the shear-driven mixing
                             !! at the cell vertices (i.e., the vorticity points).
  logical :: eliminate_massless !< If true, massless layers are merged with neighboring
                             !! massive layers in this calculation.
                             !  I can think of no good reason why this should be false. - RWH
  real    :: vel_underflow   !< Velocity components smaller than vel_underflow
                             !! are set to 0 [Z T-1 ~> m s-1].
!  logical :: layer_stagger = .false. ! If true, do the calculations centered at
                             !  layers, rather than the interfaces.
  logical :: debug = .false. !< If true, write verbose debugging messages.
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                             !! regulate the timing of diagnostic output.
  !>@{ Diagnostic IDs
  integer :: id_kd_shear = -1, id_tke = -1, id_ild2 = -1, id_dz_int = -1
  !!@}
end type kappa_shear_cs
 
! integer :: id_clock_project, id_clock_KQ, id_clock_avg, id_clock_setup
 
#undef  DEBUG
#undef  ADD_DIAGNOSTICS
 
contains
 
!> Subroutine for calculating shear-driven diffusivity and TKE in tracer columns
subroutine calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
                                 kv_io, dt, G, GV, US, CS, initialize_all)
  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(unit_scale_type),   intent(in)    :: us     !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: u_in   !< Initial zonal velocity [m s-1]. (Intent in)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: v_in   !< Initial meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: h      !< Layer thicknesses [H ~> 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.
  real, dimension(:,:),    pointer       :: p_surf !< The pressure at the ocean surface [Pa] (or NULL).
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(inout) :: kappa_io !< The diapycnal diffusivity at each interface
                                                   !! (not layer!) [Z2 T-1 ~> m2 s-1].  Initially this is the
                                                   !! value from the previous timestep, which may
                                                   !! accelerate the iteration toward convergence.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(out) :: tke_io   !< The turbulent kinetic energy per unit mass at
                                                   !! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(inout) :: kv_io  !< The vertical viscosity at each interface
                                                   !! (not layer!) [Z2 T-1 ~> m2 s-1]. This discards any
                                                   !! previous value (i.e. it is intent out) and
                                                   !! simply sets Kv = Prandtl * Kd_shear
  real,                    intent(in)    :: dt     !< Time increment [T ~> s].
  type(kappa_shear_cs),    pointer       :: cs     !< The control structure returned by a previous
                                                   !! call to kappa_shear_init.
  logical,       optional, intent(in)    :: initialize_all !< If present and false, the previous
                                                   !! value of kappa is used to start the iterations
 
  ! Local variables
  real, dimension(SZI_(G),SZK_(GV)) :: &
    h_2d, &             ! A 2-D version of h, but converted to [Z ~> m].
    u_2d, v_2d, &       ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
    t_2d, s_2d, rho_2d  ! 2-D versions of T, S, and rho.
  real, dimension(SZI_(G),SZK_(GV)+1) :: &
    kappa_2d, & ! 2-D version of kappa_io [Z2 T-1 ~> m2 s-1].
    tke_2d      ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
  real, dimension(SZK_(GV)) :: &
    idz, &      ! The inverse of the distance between TKE points [Z-1 ~> m-1].
    dz, &       ! The layer thickness [Z ~> m].
    u0xdz, &    ! The initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
    v0xdz, &    ! The initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
    t0xdz, &    ! The initial temperature times dz [degC Z ~> degC m].
    s0xdz       ! The initial salinity times dz [ppt Z ~> ppt m].
  real, dimension(SZK_(GV)+1) :: &
    kappa, &    ! The shear-driven diapycnal diffusivity at an interface [Z2 T-1 ~> m2 s-1].
    tke, &      ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
    kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
    tke_avg     ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
  real :: f2   ! The squared Coriolis parameter of each column [T-2 ~> s-2].
  real :: surface_pres  ! The top surface pressure [Pa].
 
  real :: dz_in_lay     !   The running sum of the thickness in a layer [Z ~> m].
  real :: k0dt          ! The background diffusivity times the timestep [Z2 ~> m2].
  real :: dz_massless   ! A layer thickness that is considered massless [Z ~> m].
  logical :: use_temperature  !  If true, temperature and salinity have been
                        ! allocated and are being used as state variables.
  logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
                        ! last call to this subroutine.
 
  integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
                        ! interfaces and the interfaces with massless layers
                        ! merged into nearby massive layers.
  real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
                        ! interpolating back to the original index space [nondim].
  integer :: is, ie, js, je, i, j, k, nz, nzc
 
  ! Diagnostics that should be deleted?
#ifdef ADD_DIAGNOSTICS
  real, dimension(SZK_(GV)+1) :: &  ! Additional diagnostics.
    i_ld2_1d, dz_int_1d
  real, dimension(SZI_(G),SZK_(GV)+1) :: & ! 2-D versions of diagnostics.
    i_ld2_2d, dz_int_2d
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & ! 3-D versions of diagnostics.
    i_ld2_3d, dz_int_3d
#endif
  is = g%isc ; ie = g%iec; js = g%jsc ; je = g%jec ; nz = gv%ke
 
  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
  new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all
 
  k0dt = dt*cs%kappa_0
  dz_massless = 0.1*sqrt(k0dt)
 
  !$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,new_kappa, &
#ifdef ADD_DIAGNOSTICS
  !$OMP                                I_Ld2_3d,dz_Int_3d, &
#endif
  !$OMP                                tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io)
  do j=js,je
    do k=1,nz ; do i=is,ie
      h_2d(i,k) = h(i,j,k)*gv%H_to_Z
      u_2d(i,k) = u_in(i,j,k)*us%m_s_to_L_T ; v_2d(i,k) = v_in(i,j,k)*us%m_s_to_L_T
    enddo ; enddo
    if (use_temperature) then ; do k=1,nz ; do i=is,ie
      t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
    enddo ; enddo ; else ; do k=1,nz ; do i=is,ie
      rho_2d(i,k) = gv%Rlay(k) ! Could be tv%Rho(i,j,k) ?
    enddo ; enddo ; endif
    if (.not.new_kappa) then ; do k=1,nz+1 ; do i=is,ie
      kappa_2d(i,k) = kappa_io(i,j,k)
    enddo ; enddo ; endif
 
!---------------------------------------
! Work on each column.
!---------------------------------------
    do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
    ! call cpu_clock_begin(id_clock_setup)
      ! Store a transposed version of the initial arrays.
      ! Any elimination of massless layers would occur here.
      if (cs%eliminate_massless) then
        nzc = 1
        do k=1,nz
          ! Zero out the thicknesses of all layers, even if they are unused.
          dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
          t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
 
          ! Add a new layer if this one has mass.
!          if ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless)) nzc = nzc+1
          if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
              (h_2d(i,k) > dz_massless)) nzc = nzc+1
 
          ! Only merge clusters of massless layers.
!         if ((dz(nzc) > dz_massless) .or. &
!             ((dz(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless))) nzc = nzc+1
 
          kc(k) = nzc
          dz(nzc) = dz(nzc) + h_2d(i,k)
          u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
          v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
          if (use_temperature) then
            t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
            s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
          else
            t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
            s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
          endif
        enddo
        kc(nz+1) = nzc+1
 
        ! Set up Idz as the inverse of layer thicknesses.
        do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
 
        !   Now determine kf, the fractional weight of interface kc when
        ! interpolating between interfaces kc and kc+1.
        kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
        do k=2,nz
          if (kc(k) > kc(k-1)) then
            kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
          else
            kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
          endif
        enddo
        kf(nz+1) = 0.0
      else
        do k=1,nz
          dz(k) = h_2d(i,k)
          u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
        enddo
        if (use_temperature) then
          do k=1,nz
            t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
          enddo
        else
          do k=1,nz
            t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
          enddo
        endif
        nzc = nz
        do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
      endif
      f2 = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
                   (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
      surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j)
 
    ! ----------------------------------------------------    I_Ld2_1d, dz_Int_1d
 
    ! Set the initial guess for kappa, here defined at interfaces.
    ! ----------------------------------------------------
      if (new_kappa) then
        do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ; enddo
      else
        do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k) ; enddo
      endif
 
#ifdef ADD_DIAGNOSTICS
      call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
                              dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
                              tke_avg, tv, cs, gv, us, i_ld2_1d, dz_int_1d)
#else
      call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
                              dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
                              tke_avg, tv, cs, gv, us)
#endif
 
    ! call cpu_clock_begin(id_clock_setup)
    ! Extrapolate from the vertically reduced grid back to the original layers.
      if (nz == nzc) then
        do k=1,nz+1
          kappa_2d(i,k) = kappa_avg(k)
          !### Should this be tke_avg?
          tke_2d(i,k) = tke(k)
        enddo
      else
        do k=1,nz+1
          if (kf(k) == 0.0) then
            kappa_2d(i,k) = kappa_avg(kc(k))
            tke_2d(i,k) = tke_avg(kc(k))
          else
            kappa_2d(i,k) = (1.0-kf(k)) * kappa_avg(kc(k)) + &
                             kf(k) * kappa_avg(kc(k)+1)
            tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + &
                           kf(k) * tke_avg(kc(k)+1)
          endif
        enddo
      endif
#ifdef ADD_DIAGNOSTICS
      do k=1,nz+1
        i_ld2_2d(i,k) = i_ld2_1d(k) ; dz_int_2d(i,k) = dz_int_1d(k)
      enddo
#endif
    ! call cpu_clock_end(id_clock_setup)
    else  ! Land points, still inside the i-loop.
      do k=1,nz+1
        kappa_2d(i,k) = 0.0 ; tke_2d(i,k) = 0.0
#ifdef ADD_DIAGNOSTICS
        i_ld2_2d(i,k) = 0.0  ; dz_int_2d(i,k) = 0.0
#endif
      enddo
    endif ; enddo ! i-loop
 
    do k=1,nz+1 ; do i=is,ie
      kappa_io(i,j,k) = g%mask2dT(i,j) * kappa_2d(i,k)
      tke_io(i,j,k) = g%mask2dT(i,j) * tke_2d(i,k)
      kv_io(i,j,k) = ( g%mask2dT(i,j) * kappa_2d(i,k) ) * cs%Prandtl_turb
#ifdef ADD_DIAGNOSTICS
      i_ld2_3d(i,j,k) = i_ld2_2d(i,k) ; dz_int_3d(i,j,k) = dz_int_2d(i,k)
#endif
    enddo ; enddo
 
  enddo ! end of j-loop
 
  if (cs%debug) then
    call hchksum(kappa_io, "kappa", g%HI, scale=us%Z2_T_to_m2_s)
    call hchksum(tke_io, "tke", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
  endif
 
  if (cs%id_Kd_shear > 0) call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
  if (cs%id_TKE > 0) call post_data(cs%id_TKE, tke_io, cs%diag)
#ifdef ADD_DIAGNOSTICS
  if (cs%id_ILd2 > 0) call post_data(cs%id_ILd2, i_ld2_3d, cs%diag)
  if (cs%id_dz_Int > 0) call post_data(cs%id_dz_Int, dz_int_3d, cs%diag)
#endif
 
end subroutine calculate_kappa_shear
 
 
!> Subroutine for calculating shear-driven diffusivity and TKE in corner columns
subroutine calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, &
                                   kv_io, dt, G, GV, US, CS, initialize_all)
  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(unit_scale_type),    intent(in)   :: us     !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: u_in   !< Initial zonal velocity [m s-1]. (Intent in)
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),   &
                           intent(in)    :: v_in   !< Initial meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: t_in   !< Layer potential temperatures [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   &
                           intent(in)    :: s_in   !< Layer salinities in ppt.
  type(thermo_var_ptrs),   intent(in)    :: tv     !< A structure containing pointers to any
                                                   !! available thermodynamic fields. Absent fields
                                                   !! have NULL ptrs.
  real, dimension(:,:),    pointer       :: p_surf !< The pressure at the ocean surface [Pa]
                                                   !! (or NULL).
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(out)   :: kappa_io !< The diapycnal diffusivity at each interface
                                                   !! (not layer!) [Z2 T-1 ~> m2 s-1].
  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
                           intent(out)   :: tke_io !< The turbulent kinetic energy per unit mass at
                                                   !! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
  real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
                           intent(inout) :: kv_io  !< The vertical viscosity at each interface [Z2 T-1 ~> m2 s-1].
                                                   !! The previous value is used to initialize kappa
                                                   !! in the vertex columes as Kappa = Kv/Prandtl
                                                   !! to accelerate the iteration toward covergence.
  real,                    intent(in)    :: dt     !< Time increment [T ~> s].
  type(kappa_shear_cs),    pointer       :: cs     !< The control structure returned by a previous
                                                   !! call to kappa_shear_init.
  logical,       optional, intent(in)    :: initialize_all !< If present and false, the previous
                                                   !! value of kappa is used to start the iterations
 
  ! Local variables
  real, dimension(SZIB_(G),SZK_(GV)) :: &
    h_2d, &             ! A 2-D version of h, but converted to [Z ~> m].
    u_2d, v_2d, &       ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
    t_2d, s_2d, rho_2d  ! 2-D versions of T, S, and rho.
  real, dimension(SZIB_(G),SZK_(GV)+1,2) :: &
    kappa_2d    ! Quasi 2-D versions of kappa_io [Z2 T-1 ~> m2 s-1].
  real, dimension(SZIB_(G),SZK_(GV)+1) :: &
    tke_2d      ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
  real, dimension(SZK_(GV)) :: &
    idz, &      ! The inverse of the distance between TKE points [Z-1 ~> m-1].
    dz, &       ! The layer thickness [Z ~> m].
    u0xdz, &    ! The initial zonal velocity times dz [L Z T-1 ~> m2 s-1].
    v0xdz, &    ! The initial meridional velocity times dz [L Z T-1 ~> m2 s-1].
    t0xdz, &    ! The initial temperature times dz [degC Z ~> degC m].
    s0xdz       ! The initial salinity times dz [ppt Z ~> ppt m].
  real, dimension(SZK_(GV)+1) :: &
    kappa, &    ! The shear-driven diapycnal diffusivity at an interface [Z2 T-1 ~> m2 s-1].
    tke, &      ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
    kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
    tke_avg     ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
  real :: f2   ! The squared Coriolis parameter of each column [T-2 ~> s-2].
  real :: surface_pres  ! The top surface pressure [Pa].
 
  real :: dz_in_lay     !   The running sum of the thickness in a layer [Z ~> m].
  real :: k0dt          ! The background diffusivity times the timestep [Z2 ~> m2].
  real :: dz_massless   ! A layer thickness that is considered massless [Z ~> m].
  real :: i_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1].
  real :: i_prandtl
  logical :: use_temperature  !  If true, temperature and salinity have been
                        ! allocated and are being used as state variables.
  logical :: new_kappa = .true. ! If true, ignore the value of kappa from the
                        ! last call to this subroutine.
  logical :: do_i       ! If true, work on this column.
 
  integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
                        ! interfaces and the interfaces with massless layers
                        ! merged into nearby massive layers.
  real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
                        ! interpolating back to the original index space [nondim].
  integer :: isb, ieb, jsb, jeb, i, j, k, nz, nzc, j2, j2m1
 
  ! Diagnostics that should be deleted?
#ifdef ADD_DIAGNOSTICS
  real, dimension(SZK_(GV)+1) :: &  ! Additional diagnostics.
    i_ld2_1d, dz_int_1d
  real, dimension(SZI_(G),SZK_(GV)+1) :: & ! 2-D versions of diagnostics.
    i_ld2_2d, dz_int_2d
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & ! 3-D versions of diagnostics.
    i_ld2_3d, dz_int_3d
#endif
  isb = g%isc-1 ; ieb = g%iecB ; jsb = g%jsc-1 ; jeb = g%jecB ; nz = gv%ke
 
  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
  new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all
 
  k0dt =  dt*cs%kappa_0
  dz_massless = 0.1*sqrt(k0dt)
  i_prandtl = 0.0 ; if (cs%Prandtl_turb > 0.0) i_prandtl = 1.0 / cs%Prandtl_turb
 
  !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,new_kappa, &
#ifdef ADD_DIAGNOSTICS
  !$OMP                                I_Ld2_3d,dz_Int_3d, &
#endif
  !$OMP                                tv,G,GV,US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io)
  do j=jsb,jeb
    j2 = mod(j,2)+1 ; j2m1 = 3-j2 ! = mod(J-1,2)+1
 
    ! Interpolate the various quantities to the corners, using masks.
    do k=1,nz ; do i=isb,ieb
      u_2d(i,k) = us%m_s_to_L_T * &
                  (u_in(i,j,k)   * (g%mask2dCu(i,j)   * (h(i,j,k)   + h(i+1,j,k))) + &
                   u_in(i,j+1,k) * (g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) ) / &
                  ((g%mask2dCu(i,j)   * (h(i,j,k)   + h(i+1,j,k)) + &
                    g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
      v_2d(i,k) = us%m_s_to_L_T * &
                  (v_in(i,j,k)   * (g%mask2dCv(i,j)   * (h(i,j,k)   + h(i,j+1,k))) + &
                   v_in(i+1,j,k) * (g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) ) / &
                  ((g%mask2dCv(i,j)   * (h(i,j,k)   + h(i,j+1,k)) + &
                    g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
      i_hwt = 1.0 / (((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
                      (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k))) + &
                     gv%H_subroundoff)
      if (use_temperature) then
        t_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * t_in(i,j,k) + &
                       (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * t_in(i+1,j+1,k)) + &
                      ((g%mask2dT(i+1,j) * h(i+1,j,k)) * t_in(i+1,j,k) + &
                       (g%mask2dT(i,j+1) * h(i,j+1,k)) * t_in(i,j+1,k)) ) * i_hwt
        s_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * s_in(i,j,k) + &
                       (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * s_in(i+1,j+1,k)) + &
                      ((g%mask2dT(i+1,j) * h(i+1,j,k)) * s_in(i+1,j,k) + &
                       (g%mask2dT(i,j+1) * h(i,j+1,k)) * s_in(i,j+1,k)) ) * i_hwt
      endif
      h_2d(i,k) = gv%H_to_Z * ((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
                               (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k)) ) / &
                              ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + &
                               (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) + 1.0e-36 )
!      h_2d(I,k) = 0.25*((h(i,j,k) + h(i+1,j+1,k)) + (h(i+1,j,k) + h(i,j+1,k)))*GV%H_to_Z
!      h_2d(I,k) = ((h(i,j,k)**2 + h(i+1,j+1,k)**2) + &
!                   (h(i+1,j,k)**2 + h(i,j+1,k)**2))*GV%H_to_Z * I_hwt
    enddo ; enddo
    if (.not.use_temperature) then ; do k=1,nz ; do i=isb,ieb
      rho_2d(i,k) = gv%Rlay(k)
    enddo ; enddo ; endif
    if (.not.new_kappa) then ; do k=1,nz+1 ; do i=isb,ieb
      kappa_2d(i,k,j2) = kv_io(i,j,k) * i_prandtl
    enddo ; enddo ; endif
 
!---------------------------------------
! Work on each column.
!---------------------------------------
    do i=isb,ieb ; if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
                       (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0) then
    ! call cpu_clock_begin(Id_clock_setup)
      ! Store a transposed version of the initial arrays.
      ! Any elimination of massless layers would occur here.
      if (cs%eliminate_massless) then
        nzc = 1
        do k=1,nz
          ! Zero out the thicknesses of all layers, even if they are unused.
          dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
          t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
 
          ! Add a new layer if this one has mass.
!          if ((dz(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1
          if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
              (h_2d(i,k) > dz_massless)) nzc = nzc+1
 
          ! Only merge clusters of massless layers.
!         if ((dz(nzc) > dz_massless) .or. &
!             ((dz(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1
 
          kc(k) = nzc
          dz(nzc) = dz(nzc) + h_2d(i,k)
          u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
          v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
          if (use_temperature) then
            t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
            s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
          else
            t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
            s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
          endif
        enddo
        kc(nz+1) = nzc+1
 
        ! Set up Idz as the inverse of layer thicknesses.
        do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
 
        !   Now determine kf, the fractional weight of interface kc when
        ! interpolating between interfaces kc and kc+1.
        kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
        do k=2,nz
          if (kc(k) > kc(k-1)) then
            kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
          else
            kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
          endif
        enddo
        kf(nz+1) = 0.0
      else
        do k=1,nz
          dz(k) = h_2d(i,k)
          u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
        enddo
        if (use_temperature) then
          do k=1,nz
            t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
          enddo
        else
          do k=1,nz
            t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
          enddo
        endif
        nzc = nz
        do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
      endif
      f2 = g%CoriolisBu(i,j)**2
      surface_pres = 0.0 ; if (associated(p_surf)) &
        surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + &
                               (p_surf(i+1,j) + p_surf(i,j+1)))
 
    ! ----------------------------------------------------
    ! Set the initial guess for kappa, here defined at interfaces.
    ! ----------------------------------------------------
      if (new_kappa) then
        do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ; enddo
      else
        do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k,j2) ; enddo
      endif
 
#ifdef ADD_DIAGNOSTICS
      call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
                              dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
                              tke_avg, tv, cs, gv, us, i_ld2_1d, dz_int_1d)
#else
      call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
                              dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
                              tke_avg, tv, cs, gv, us)
#endif
    ! call cpu_clock_begin(Id_clock_setup)
    ! Extrapolate from the vertically reduced grid back to the original layers.
      if (nz == nzc) then
        do k=1,nz+1
          kappa_2d(i,k,j2) = kappa_avg(k)
          !### Should this be tke_avg?
          tke_2d(i,k) = tke(k)
        enddo
      else
        do k=1,nz+1
          if (kf(k) == 0.0) then
            kappa_2d(i,k,j2) = kappa_avg(kc(k))
            tke_2d(i,k) = tke_avg(kc(k))
          else
            kappa_2d(i,k,j2) = (1.0-kf(k)) * kappa_avg(kc(k)) + kf(k) * kappa_avg(kc(k)+1)
            tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + kf(k) * tke_avg(kc(k)+1)
          endif
        enddo
      endif
#ifdef ADD_DIAGNOSTICS
      do k=1,nz+1
        i_ld2_2d(i,k) = i_ld2_1d(k) ; dz_int_2d(i,k) = dz_int_1d(k)
      enddo
#endif
    ! call cpu_clock_end(Id_clock_setup)
    else  ! Land points, still inside the i-loop.
      do k=1,nz+1
        kappa_2d(i,k,j2) = 0.0 ; tke_2d(i,k) = 0.0
#ifdef ADD_DIAGNOSTICS
        i_ld2_2d(i,k) = 0.0 ; dz_int_2d(i,k) = 0.0
#endif
      enddo
    endif ; enddo ! i-loop
 
    do k=1,nz+1 ; do i=isb,ieb
      tke_io(i,j,k) = g%mask2dBu(i,j) * tke_2d(i,k)
      kv_io(i,j,k) = ( g%mask2dBu(i,j) * kappa_2d(i,k,j2) ) * cs%Prandtl_turb
#ifdef ADD_DIAGNOSTICS
      i_ld2_3d(i,j,k) = i_ld2_2d(i,k) ; dz_int_3d(i,j,k) = dz_int_2d(i,k)
#endif
    enddo ; enddo
    if (j>=g%jsc) then ; do k=1,nz+1 ; do i=g%isc,g%iec
      ! Set the diffusivities in tracer columns from the values at vertices.
      kappa_io(i,j,k) = g%mask2dT(i,j) * 0.25 * &
                        ((kappa_2d(i-1,k,j2m1) + kappa_2d(i,k,j2)) + &
                         (kappa_2d(i-1,k,j2)   + kappa_2d(i,k,j2m1)))
    enddo ; enddo ; endif
 
  enddo ! end of J-loop
 
  if (cs%debug) then
    call hchksum(kappa_io, "kappa", g%HI, scale=us%Z2_T_to_m2_s)
    call bchksum(tke_io, "tke", g%HI)
  endif
 
  if (cs%id_Kd_shear > 0) call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
  if (cs%id_TKE > 0) call post_data(cs%id_TKE, tke_io, cs%diag)
#ifdef ADD_DIAGNOSTICS
  if (cs%id_ILd2 > 0) call post_data(cs%id_ILd2, i_ld2_3d, cs%diag)
  if (cs%id_dz_Int > 0) call post_data(cs%id_dz_Int, dz_int_3d, cs%diag)
#endif
 
end subroutine calc_kappa_shear_vertex
 
 
!> This subroutine calculates shear-driven diffusivity and TKE in a single column
subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
                              dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
                              tke_avg, tv, CS, GV, US, I_Ld2_1d, dz_Int_1d)
  type(verticalgrid_type), intent(in)    :: GV !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
  real, dimension(SZK_(GV)+1), &
                     intent(inout) :: kappa !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)+1), &
                     intent(out)   :: tke  !< The Turbulent Kinetic Energy per unit mass at
                                           !! an interface [Z2 T-2 ~> m2 s-2].
  integer,           intent(in)    :: nzc  !< The number of active layers in the column.
  real,              intent(in)    :: f2   !< The square of the Coriolis parameter [T-2 ~> s-2].
  real,              intent(in)    :: surface_pres  !< The surface pressure [Pa].
  real, dimension(SZK_(GV)), &
                     intent(in)    :: dz   !< The layer thickness [Z ~> m].
  real, dimension(SZK_(GV)), &
                     intent(in)    :: u0xdz !< The initial zonal velocity times dz [Z L T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)), &
                     intent(in)    :: v0xdz !< The initial meridional velocity times dz [Z L T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)), &
                     intent(in)    :: T0xdz !< The initial temperature times dz [degC Z ~> degC m].
  real, dimension(SZK_(GV)), &
                     intent(in)    :: S0xdz !< The initial salinity times dz [ppt Z ~> ppt m].
  real, dimension(SZK_(GV)+1), &
                     intent(out)   :: kappa_avg !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)+1), &
                     intent(out)   :: tke_avg  !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
  real,                    intent(in)    :: dt !< Time increment [T ~> s].
  type(thermo_var_ptrs),   intent(in)    :: tv !< A structure containing pointers to any
                                               !! available thermodynamic fields. Absent fields
                                               !! have NULL ptrs.
  type(kappa_shear_cs),    pointer       :: CS !< The control structure returned by a previous
                                               !! call to kappa_shear_init.
  real,  dimension(SZK_(GV)+1), &
           optional, intent(out)   :: I_Ld2_1d !< The inverse of the squared mixing length [Z-2 ~> m-2].
  real,  dimension(SZK_(GV)+1), &
           optional, intent(out)   :: dz_Int_1d !< The extent of a finite-volume space surrounding an interface,
                                               !! as used in calculating kappa and TKE [Z ~> m].
 
  real, dimension(nzc) :: &
    u, &        ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
    v, &        ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
    Idz, &      ! The inverse of the distance between TKE points [Z-1 ~> m-1].
    T, &        ! The potential temperature after a timestep of mixing [degC].
    Sal, &      ! The salinity after a timestep of mixing [ppt].
    u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
    T_test, S_test ! Temporary temperatures [degC] and salinities [ppt].
 
  real, dimension(nzc+1) :: &
    N2, &       ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
    dz_Int, &   ! The extent of a finite-volume space surrounding an interface,
                ! as used in calculating kappa and TKE [Z ~> m].
    i_dz_int, & ! The inverse of the distance between velocity & density points
                ! above and below an interface [Z-1 ~> m-1].  This is used to
                ! calculate N2, shear, and fluxes, and it might differ from
                ! 1/dz_Int, as they have different uses.
    s2, &       ! The squared shear at an interface [T-2 ~> s-2].
    a1, &       ! a1 is the coupling between adjacent interfaces in the TKE,
                ! velocity, and density equations [Z s-1 ~> m s-1] or [Z ~> m]
    c1, &       ! c1 is used in the tridiagonal (and similar) solvers.
    k_src, &    ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
    kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1].
    kappa_out, & ! The kappa that results from the kappa equation [Z2 T-1 ~> m2 s-1].
    kappa_mid, & ! The average of the initial and predictor estimates of kappa [Z2 T-1 ~> m2 s-1].
    tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
    kappa_pred, & ! The value of kappa from a predictor step [Z2 T-1 ~> m2 s-1].
    pressure, & ! The pressure at an interface [Pa].
    t_int, &    ! The temperature interpolated to an interface [degC].
    sal_int, &  ! The salinity interpolated to an interface [ppt].
    dbuoy_dt, & ! The partial derivatives of buoyancy with changes in temperature
    dbuoy_ds, & ! and salinity, [Z T-2 degC-1 ~> m s-2 degC-1] and [Z T-2 ppt-1 ~> m s-2 ppt-1].
    i_l2_bdry, &   ! The inverse of the square of twice the harmonic mean
                   ! distance to the top and bottom boundaries [Z-2 ~> m-2].
    k_q, &         ! Diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s].
    k_q_tmp, &     ! A temporary copy of diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s].
    local_src_avg, & ! The time-integral of the local source [nondim].
    tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
    tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
    tol_chg, & ! The tolerated change integrated in time [s T-nondim].
    dist_from_top, &  ! The distance from the top surface [Z ~> m].
    local_src     ! The sum of all sources of kappa, including kappa_src and
                  ! sources from the elliptic term [T-1 ~> s-1].
 
  real :: dist_from_bot ! The distance from the bottom surface [Z ~> m].
  real :: b1            ! The inverse of the pivot in the tridiagonal equations.
  real :: bd1           ! A term in the denominator of b1.
  real :: d1            ! 1 - c1 in the tridiagonal equations.
  real :: gR0           ! A conversion factor from Z to Pa equal to Rho_0 times g
                        ! [kg m-1 Z-1 s-2 ~> kg m-2 s-2].
  real :: g_R0          ! g_R0 is a rescaled version of g/Rho [Z m3 kg-1 T-2 ~> m4 kg-1 s-2].
  real :: Norm          ! A factor that normalizes two weights to 1 [Z-2 ~> m-2].
  real :: tol_dksrc, tol2  ! ### Tolerances that need to be set better later.
  real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
                        ! within an iteration.  0 < tol_dksrc_low < 1.
  real :: Ri_crit       !   The critical shear Richardson number for shear-
                        ! driven mixing. The theoretical value is 0.25.
  real :: dt_rem        !   The remaining time to advance the solution [T ~> s].
  real :: dt_now        !   The time step used in the current iteration [T ~> s].
  real :: dt_wt         !   The fractional weight of the current iteration [nondim].
  real :: dt_test       !   A time-step that is being tested for whether it
                        ! gives acceptably small changes in k_src [T ~> s].
  real :: Idtt          !   Idtt = 1 / dt_test [T-1 ~> s-1].
  real :: dt_inc        !   An increment to dt_test that is being tested [T ~> s].
 
  real :: k0dt          ! The background diffusivity times the timestep [Z2 ~> m2].
  logical :: valid_dt   ! If true, all levels so far exhibit acceptably small
                        ! changes in k_src.
  logical :: use_temperature  !  If true, temperature and salinity have been
                        ! allocated and are being used as state variables.
  integer :: ks_kappa, ke_kappa  ! The k-range with nonzero kappas.
  integer :: dt_halvings   ! The number of times that the time-step is halved
                           ! in seeking an acceptable timestep.  If none is
                           ! found, dt_rem*0.5^dt_halvings is used.
  integer :: dt_refinements ! The number of 2-fold refinements that will be used
                           ! to estimate the maximum permitted time step.  I.e.,
                           ! the resolution is 1/2^dt_refinements.
  integer :: k, itt, itt_dt
#ifdef DEBUG
  integer :: max_debug_itt ; parameter(max_debug_itt=20)
  real :: wt(SZK_(GV)+1), wt_tot, I_wt_tot, wt_itt
  real, dimension(SZK_(GV)+1) :: &
    Ri_k, tke_prev, dtke, dkappa, dtke_norm, &
    N2_debug, & ! A version of N2 for debugging [T-2 ~> s-2]
    ksrc_av     ! The average through the iterations of k_src [T-1 ~> s-1].
  real, dimension(SZK_(GV)+1,0:max_debug_itt) :: &
    tke_it1, N2_it1, Sh2_it1, ksrc_it1, kappa_it1, kprev_it1
  real, dimension(SZK_(GV)+1,1:max_debug_itt) :: &
    dkappa_it1, wt_it1, K_Q_it1, d_dkappa_it1, dkappa_norm
  real, dimension(SZK_(GV),0:max_debug_itt) :: &
    u_it1, v_it1, rho_it1, T_it1, S_it1
  real, dimension(0:max_debug_itt) :: &
    dk_wt_it1, dkpos_wt_it1, dkneg_wt_it1, k_mag
  real, dimension(max_debug_itt) ::  dt_it1
#endif
 
  ri_crit = cs%Rino_crit
  gr0 = gv%z_to_H*gv%H_to_Pa
  g_r0 = (us%L_to_Z**2 * gv%g_Earth) / gv%Rho0
  k0dt = dt*cs%kappa_0
  ! These are hard-coded for now.  Perhaps these could be made dynamic later?
  ! tol_dksrc = 0.5*tol_ksrc_chg ; tol_dksrc_low = 1.0 - 1.0/tol_ksrc_chg ?
  tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*cs%kappa_tol_err
  dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
  use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
 
 
  ! Set up Idz as the inverse of layer thicknesses.
  do k=1,nzc ; idz(k) = 1.0 / dz(k) ; enddo
  !   Set up I_dz_int as the inverse of the distance between
  ! adjacent layer centers.
  i_dz_int(1) = 2.0 / dz(1)
  dist_from_top(1) = 0.0
  do k=2,nzc
    i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
    dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
  enddo
  i_dz_int(nzc+1) = 2.0 / dz(nzc)
 
  !   Determine the velocities and thicknesses after eliminating massless
  ! layers and applying a time-step of background diffusion.
  if (nzc > 1) then
    a1(2) = k0dt*i_dz_int(2)
    b1 = 1.0 / (dz(1) + a1(2))
    u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
    t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
    c1(2) = a1(2) * b1 ; d1 = dz(1) * b1 ! = 1 - c1
    do k=2,nzc-1
      bd1 = dz(k) + d1*a1(k)
      a1(k+1) = k0dt*i_dz_int(k+1)
      b1 = 1.0 / (bd1 + a1(k+1))
      u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
      v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
      t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
      sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
      c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
    enddo
    ! rho or T and S have insulating boundary conditions, u & v use no-slip
    ! bottom boundary conditions (if kappa0 > 0).
    ! For no-slip bottom boundary conditions
    b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
    u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
    v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
    ! For insulating boundary conditions
    b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
    t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
    sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
    do k=nzc-1,1,-1
      u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
      t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
    enddo
  else
    ! This is correct, but probably unnecessary.
    b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
    u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
    b1 = 1.0 / dz(1)
    t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
  endif
 
  ! This uses half the harmonic mean of thicknesses to provide two estimates
  ! of the boundary between cells, and the inverse of the harmonic mean to
  ! weight the two estimates.  The net effect is that interfaces around thin
  ! layers have thin cells, and the total thickness adds up properly.
  ! The top- and bottom- interfaces have zero thickness, consistent with
  ! adding additional zero thickness layers.
  dz_int(1) = 0.0 ; dz_int(2) = dz(1)
  do k=2,nzc-1
    norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
    dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
    dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
  enddo
  dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
 
  dist_from_bot = 0.0
  do k=nzc,2,-1
    dist_from_bot = dist_from_bot + dz(k)
    i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
                   (dist_from_top(k) * dist_from_bot)**2
  enddo
 
  ! Calculate thermodynamic coefficients and an initial estimate of N2.
  if (use_temperature) then
    pressure(1) = surface_pres
    do k=2,nzc
      pressure(k) = pressure(k-1) + gr0*dz(k-1)
      t_int(k) = 0.5*(t(k-1) + t(k))
      sal_int(k) = 0.5*(sal(k-1) + sal(k))
    enddo
    call calculate_density_derivs(t_int, sal_int, pressure, dbuoy_dt, &
                                  dbuoy_ds, 2, nzc-1, tv%eqn_of_state)
    do k=2,nzc
      dbuoy_dt(k) = -g_r0*dbuoy_dt(k)
      dbuoy_ds(k) = -g_r0*dbuoy_ds(k)
    enddo
  else
    do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ; enddo
  endif
 
#ifdef DEBUG
  n2_debug(1) = 0.0 ; n2_debug(nzc+1) = 0.0
  do k=2,nzc
    n2_debug(k) = max((dbuoy_dt(k) * (t0xdz(k-1)*idz(k-1) - t0xdz(k)*idz(k)) + &
                       dbuoy_ds(k) * (s0xdz(k-1)*idz(k-1) - s0xdz(k)*idz(k))) * &
                       i_dz_int(k), 0.0)
  enddo
  do k=1,nzc
    u_it1(k,0) = u0xdz(k)*idz(k) ; v_it1(k,0) = v0xdz(k)*idz(k)
    t_it1(k,0) = t0xdz(k)*idz(k) ; s_it1(k,0) = s0xdz(k)*idz(k)
  enddo
  do k=1,nzc+1
    kprev_it1(k,0) = kappa(k) ; kappa_it1(k,0) = kappa(k)
    tke_it1(k,0) = 0.0
    n2_it1(k,0) = n2_debug(k) ; sh2_it1(k,0) = s2(k) ; ksrc_it1(k,0) = k_src(k)
  enddo
  do k=nzc+1,gv%ke
    u_it1(k,0) = 0.0 ; v_it1(k,0) = 0.0
    t_it1(k,0) = 0.0 ; s_it1(k,0) = 0.0
    kprev_it1(k+1,0) = 0.0 ; kappa_it1(k+1,0) = 0.0 ; tke_it1(k+1,0) = 0.0
    n2_it1(k+1,0) = 0.0 ; sh2_it1(k+1,0) = 0.0 ; ksrc_it1(k+1,0) = 0.0
  enddo
  do itt=1,max_debug_itt
    dt_it1(itt) = 0.0
    do k=1,gv%ke
      u_it1(k,itt) = 0.0 ; v_it1(k,itt) = 0.0
      t_it1(k,itt) = 0.0 ; s_it1(k,itt) = 0.0
      rho_it1(k,itt) = 0.0
    enddo
    do k=1,gv%ke+1
      kprev_it1(k,itt) = 0.0 ; kappa_it1(k,itt) = 0.0 ; tke_it1(k,itt) = 0.0
      n2_it1(k,itt) = 0.0 ; sh2_it1(k,itt) = 0.0
      ksrc_it1(k,itt) = 0.0
      dkappa_it1(k,itt) = 0.0 ; wt_it1(k,itt) = 0.0
      k_q_it1(k,itt) = 0.0 ; d_dkappa_it1(k,itt) = 0.0
    enddo
  enddo
  do k=1,gv%ke+1 ; ksrc_av(k) = 0.0 ; enddo
#endif
 
  ! This call just calculates N2 and S2.
  call calculate_projected_state(kappa, u, v, t, sal, 0.0, nzc, dz, i_dz_int, &
                                 dbuoy_dt, dbuoy_ds, u, v, t, sal, gv, us, &
                                 n2=n2, s2=s2, vel_underflow=cs%vel_underflow)
! ----------------------------------------------------
! Iterate
! ----------------------------------------------------
  dt_rem = dt
  do k=1,nzc+1
    k_q(k) = 0.0
    kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
    local_src_avg(k) = 0.0
    ! Use the grid spacings to scale errors in the source.
    if ( dz_int(k) > 0.0 ) &
      local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
  enddo
 
! call cpu_clock_end(id_clock_setup)
 
! do itt=1,CS%max_RiNo_it
  do itt=1,cs%max_KS_it
 
! ----------------------------------------------------
! Calculate new values of u, v, rho, N^2 and S.
! ----------------------------------------------------
#ifdef DEBUG
    do k=1,nzc+1
      ri_k(k) = 1e3 ; if (s2(k) > 1e-3*n2(k)) ri_k(k) = n2(k) / s2(k)
      if (itt > 1) then ; tke_prev(k) = tke(k) ; else ; tke_prev(k) = 0.0 ; endif
    enddo
#endif
 
  ! call cpu_clock_begin(id_clock_KQ)
    call find_kappa_tke(n2, s2, kappa, idz, dz_int, i_l2_bdry, f2, &
                        nzc, cs, gv, us, k_q, tke, kappa_out, kappa_src, local_src)
  ! call cpu_clock_end(id_clock_KQ)
 
  ! call cpu_clock_begin(id_clock_avg)
    ! Determine the range of non-zero values of kappa_out.
    ks_kappa = gv%ke+1 ; ke_kappa = 0
    do k=2,nzc ; if (kappa_out(k) > 0.0) then
      ks_kappa = k ; exit
    endif ; enddo
    do k=nzc,ks_kappa,-1 ; if (kappa_out(k) > 0.0) then
      ke_kappa = k ; exit
    endif ; enddo
    if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
  ! call cpu_clock_end(id_clock_avg)
 
    ! Determine how long to use this value of kappa (dt_now).
 
  ! call cpu_clock_begin(id_clock_project)
    if ((ke_kappa < ks_kappa) .or. (itt==cs%max_RiNo_it)) then
      dt_now = dt_rem
    else
      ! Limit dt_now so that |k_src(k)-kappa_src(k)| < tol * local_src(k)
      dt_test = dt_rem
      do k=2,nzc
        tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
        tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
        tol_chg(k) = tol2 * local_src_avg(k)
      enddo
 
      do itt_dt=1,(cs%max_KS_it+1-itt)/2
        !   The maximum number of times that the time-step is halved in
        ! seeking an acceptable timestep is reduced with each iteration,
        ! so that as the maximum number of iterations is approached, the
        ! whole remaining timestep is used.  Typically, an acceptable
        ! timestep is found long before the minimum is reached, so the
        ! value of max_KS_it may be unimportant, especially if it is large
        ! enough.
        call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*dt_test, nzc, dz, i_dz_int, &
                                       dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
                                       gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, &
                                       vel_underflow=cs%vel_underflow)
        valid_dt = .true.
        idtt = 1.0 / dt_test
        do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
          if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
            k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
                       ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
            if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
                (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
              valid_dt = .false. ; exit
            endif
          else
            if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
              valid_dt = .false. ; k_src(k) = 0.0 ; exit
            endif
          endif
        enddo
 
        if (valid_dt) exit
        dt_test = 0.5*dt_test
      enddo
      if ((dt_test < dt_rem) .and. valid_dt) then
        dt_inc = 0.5*dt_test
        do itt_dt=1,dt_refinements
          call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*(dt_test+dt_inc), &
                   nzc, dz, i_dz_int, dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
                   gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=cs%vel_underflow)
          valid_dt = .true.
          idtt = 1.0 / (dt_test+dt_inc)
          do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
            if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
              k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
                         ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
              if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
                  (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))) then
                valid_dt = .false. ; exit
              endif
            else
              if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))) then
                valid_dt = .false. ; k_src(k) = 0.0 ; exit
              endif
            endif
          enddo
 
          if (valid_dt) dt_test = dt_test + dt_inc
          dt_inc = 0.5*dt_inc
        enddo
      else
        dt_inc = 0.0
      endif
 
      dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc, dt_rem)
      do k=2,nzc
        local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
      enddo
    endif  ! Are all the values of kappa_out 0?
  ! call cpu_clock_end(id_clock_project)
 
    ! The state has already been projected forward. Now find new values of kappa.
 
    if (ke_kappa < ks_kappa) then
      ! There is no mixing now, and will not be again.
    ! call cpu_clock_begin(id_clock_avg)
      dt_wt = dt_rem / dt ; dt_rem = 0.0
      do k=1,nzc+1
        kappa_mid(k) = 0.0
        ! This would be here but does nothing.
        ! kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
        tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
#ifdef DEBUG
        tke_pred(k) = tke(k) ; kappa_pred(k) = 0.0 ; kappa(k) = 0.0
#endif
      enddo
    ! call cpu_clock_end(id_clock_avg)
    else
    ! call cpu_clock_begin(id_clock_project)
      call calculate_projected_state(kappa_out, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
                                     dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
                                     gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
                                     vel_underflow=cs%vel_underflow)
    ! call cpu_clock_end(id_clock_project)
 
    ! call cpu_clock_begin(id_clock_KQ)
      do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ; enddo
      call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
                          nzc, cs, gv, us, k_q_tmp, tke_pred, kappa_pred)
    ! call cpu_clock_end(id_clock_KQ)
 
      ks_kappa = gv%ke+1 ; ke_kappa = 0
      do k=1,nzc+1
        kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
        if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
        if (kappa_mid(k) > 0.0) ke_kappa = k
      enddo
 
    ! call cpu_clock_begin(id_clock_project)
      call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
                                     dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
                                     gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
                                     vel_underflow=cs%vel_underflow)
    ! call cpu_clock_end(id_clock_project)
 
    ! call cpu_clock_begin(id_clock_KQ)
      call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
                          nzc, cs, gv, us, k_q, tke_pred, kappa_pred)
    ! call cpu_clock_end(id_clock_KQ)
 
    ! call cpu_clock_begin(id_clock_avg)
      dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
      do k=1,nzc+1
        kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
        kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
        tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
        kappa(k) = kappa_pred(k) ! First guess for the next iteration.
      enddo
    ! call cpu_clock_end(id_clock_avg)
    endif
 
    if (dt_rem > 0.0) then
      ! Update the values of u, v, T, Sal, N2, and S2 for the next iteration.
    ! call cpu_clock_begin(id_clock_project)
      call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
                                     dz, i_dz_int, dbuoy_dt, dbuoy_ds, u, v, t, sal, &
                                     gv, us, n2, s2, vel_underflow=cs%vel_underflow)
    ! call cpu_clock_end(id_clock_project)
    endif
 
#ifdef DEBUG
    if (itt <= max_debug_itt) then
      dt_it1(itt) = dt_now
      dk_wt_it1(itt) = 0.0 ; dkpos_wt_it1(itt) = 0.0 ;  dkneg_wt_it1(itt) = 0.0
      k_mag(itt) = 0.0
      wt_itt = 1.0/real(itt) ; wt_tot = 0.0
      do k=1,nzc+1
        ksrc_av(k) = (1.0-wt_itt)*ksrc_av(k) + wt_itt*k_src(k)
        wt_tot = wt_tot + dz_int(k) * ksrc_av(k)
      enddo
      ! Use the 1/0=0 convention.
      i_wt_tot = 0.0 ; if (wt_tot > 0.0) i_wt_tot = 1.0/wt_tot
 
      do k=1,nzc+1
        wt(k) = (dz_int(k)*ksrc_av(k)) * i_wt_tot
        k_mag(itt) = k_mag(itt) + wt(k)*kappa_mid(k)
        dkappa_it1(k,itt) = kappa_pred(k) - kappa_out(k)
        dk_wt_it1(itt) = dk_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
        if (dkappa_it1(k,itt) > 0.0) then
          dkpos_wt_it1(itt) = dkpos_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
        else
          dkneg_wt_it1(itt) = dkneg_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
        endif
        wt_it1(k,itt) = wt(k)
      enddo
    endif
    do k=1,nzc+1
      ri_k(k) = 1e3 ; if (n2(k) < 1e3 * s2(k)) ri_k(k) = n2(k) / s2(k)
      dtke(k) = tke_pred(k) - tke(k)
      dtke_norm(k) = dtke(k) / (0.5*(tke(k) + tke_pred(k)))
      dkappa(k) = kappa_pred(k) - kappa_out(k)
    enddo
    if (itt <= max_debug_itt) then
      do k=1,nzc
        u_it1(k,itt) = u(k) ; v_it1(k,itt) = v(k)
        t_it1(k,itt) = t(k) ; s_it1(k,itt) = sal(k)
      enddo
      do k=1,nzc+1
        kprev_it1(k,itt) = kappa_out(k)
        kappa_it1(k,itt) = kappa_mid(k) ; tke_it1(k,itt) = 0.5*(tke(k)+tke_pred(k))
        n2_it1(k,itt)=n2(k) ; sh2_it1(k,itt)=s2(k)
        ksrc_it1(k,itt) = kappa_src(k)
        k_q_it1(k,itt) = kappa_out(k) / (tke(k))
        if (itt > 1) then
          if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
            d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
        endif
        dkappa_norm(k,itt) = dkappa(k) / max(0.5*(kappa_pred(k) + kappa_out(k)), us%m2_s_to_Z2_T*1e-100)
      enddo
    endif
#endif
 
    if (dt_rem <= 0.0) exit
 
  enddo ! end itt loop
 
#ifdef ADD_DIAGNOSTICS
  if (present(i_ld2_1d)) then
    do k=1,gv%ke+1 ; i_ld2_1d(k) = 0.0 ; enddo
    do k=2,nzc ; if (tke(k) > 0.0) &
      i_ld2_1d(k) = i_l2_bdry(k) + (n2(k) / cs%lambda**2 + f2) / tke(k)
    enddo
  endif
  if (present(dz_int_1d)) then
    do k=1,nzc+1 ; dz_int_1d(k) = dz_int(k) ; enddo
    do k=nzc+2,gv%ke ; dz_int_1d(k) = 0.0 ; enddo
  endif
#endif
 
end subroutine kappa_shear_column
 
!>   This subroutine calculates the velocities, temperature and salinity that
!! the water column will have after mixing for dt with diffusivities kappa.  It
!! may also calculate the projected buoyancy frequency and shear.
subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, &
                                     dz, I_dz_int, dbuoy_dT, dbuoy_dS, &
                                     u, v, T, Sal, GV, US, N2, S2, ks_int, ke_int, vel_underflow)
  integer,               intent(in)    :: nz  !< The number of layers (after eliminating massless
                                              !! layers?).
  real, dimension(nz+1), intent(in)    :: kappa !< The diapycnal diffusivity at interfaces,
                                              !! [Z2 T-1 ~> m2 s-1].
  real, dimension(nz),   intent(in)    :: u0  !< The initial zonal velocity [m s-1].
  real, dimension(nz),   intent(in)    :: v0  !< The initial meridional velocity [m s-1].
  real, dimension(nz),   intent(in)    :: T0  !< The initial temperature [degC].
  real, dimension(nz),   intent(in)    :: S0  !< The initial salinity [ppt].
  real, dimension(nz),   intent(in)    :: dz  !< The grid spacing of layers [Z ~> m].
  real, dimension(nz+1), intent(in)    :: I_dz_int !< The inverse of the layer's thicknesses
                                              !! [Z-1 ~> m-1].
  real, dimension(nz+1), intent(in)    :: dbuoy_dT !< The partial derivative of buoyancy with
                                              !! temperature [Z T-2 degC-1 ~> m s-2 degC-1].
  real, dimension(nz+1), intent(in)    :: dbuoy_dS !< The partial derivative of buoyancy with
                                              !! salinity [Z T-2 ppt-1 ~> m s-2 ppt-1].
  real,                  intent(in)    :: dt  !< The time step [T ~> s].
  real, dimension(nz),   intent(inout) :: u   !< The zonal velocity after dt [m s-1].
  real, dimension(nz),   intent(inout) :: v   !< The meridional velocity after dt [m s-1].
  real, dimension(nz),   intent(inout) :: T   !< The temperature after dt [degC].
  real, dimension(nz),   intent(inout) :: Sal !< The salinity after dt [ppt].
  type(verticalgrid_type), intent(in)  :: GV  !< The ocean's vertical grid structure.
  type(unit_scale_type), intent(in)    :: US  !< A dimensional unit scaling type
  real, dimension(nz+1), optional, &
                         intent(inout) :: N2  !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
  real, dimension(nz+1), optional, &
                         intent(inout) :: S2  !< The squared shear at interfaces [T-2 ~> s-2].
  integer, optional,     intent(in)    :: ks_int !< The topmost k-index with a non-zero diffusivity.
  integer, optional,     intent(in)    :: ke_int !< The bottommost k-index with a non-zero
                                              !! diffusivity.
  real,    optional,     intent(in)    :: vel_underflow !< If present and true, any velocities that
                                              !! are smaller in magnitude than this value are
                                              !! set to 0 [m s-1].
 
  ! Local variables
  real, dimension(nz+1) :: c1
  real :: L2_to_Z2       ! A conversion factor from horizontal length units to vertical depth
                         ! units squared [Z2 s2 T-2 m-2 ~> 1].
  real :: underflow_vel  ! Velocities smaller in magnitude than underflow_vel are set to 0 [m s-1].
  real :: a_a, a_b, b1, d1, bd1, b1nz_0
  integer :: k, ks, ke
 
  ks = 1 ; ke = nz
  if (present(ks_int)) ks = max(ks_int-1,1)
  if (present(ke_int)) ke = min(ke_int,nz)
  underflow_vel = 0.0 ; if (present(vel_underflow)) underflow_vel = vel_underflow
 
  if (ks > ke) return
 
  if (dt > 0.0) then
    a_b = dt*(kappa(ks+1)*i_dz_int(ks+1))
    b1 = 1.0 / (dz(ks) + a_b)
    c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
 
    u(ks) = (b1 * dz(ks))*u0(ks) ; v(ks) = (b1 * dz(ks))*v0(ks)
    t(ks) = (b1 * dz(ks))*t0(ks) ; sal(ks) = (b1 * dz(ks))*s0(ks)
    do k=ks+1,ke-1
      a_a = a_b
      a_b = dt*(kappa(k+1)*i_dz_int(k+1))
      bd1 = dz(k) + d1*a_a
      b1 = 1.0 / (bd1 + a_b)
      c1(k+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
 
      u(k) = b1 * (dz(k)*u0(k) + a_a*u(k-1))
      v(k) = b1 * (dz(k)*v0(k) + a_a*v(k-1))
      t(k) = b1 * (dz(k)*t0(k) + a_a*t(k-1))
      sal(k) = b1 * (dz(k)*s0(k) + a_a*sal(k-1))
    enddo
    !   T and S have insulating boundary conditions, u & v use no-slip
    ! bottom boundary conditions at the solid bottom.
 
    ! For insulating boundary conditions or mixing simply stopping, use...
    a_a = a_b
    b1 = 1.0 / (dz(ke) + d1*a_a)
    t(ke) = b1 * (dz(ke)*t0(ke) + a_a*t(ke-1))
    sal(ke) = b1 * (dz(ke)*s0(ke) + a_a*sal(ke-1))
 
    !   There is no distinction between the effective boundary conditions for
    ! tracers and velocities if the mixing is separated from the bottom, but if
    ! the mixing goes all the way to the bottom, use no-slip BCs for velocities.
    if (ke == nz) then
      a_b = dt*(kappa(nz+1)*i_dz_int(nz+1))
      b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
    else
      b1nz_0 = b1
    endif
    u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1))
    v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1))
    if (abs(u(ke)) < underflow_vel) u(ke) = 0.0
    if (abs(v(ke)) < underflow_vel) v(ke) = 0.0
 
    do k=ke-1,ks,-1
      u(k) = u(k) + c1(k+1)*u(k+1)
      v(k) = v(k) + c1(k+1)*v(k+1)
      if (abs(u(k)) < underflow_vel) u(k) = 0.0
      if (abs(v(k)) < underflow_vel) v(k) = 0.0
      t(k) = t(k) + c1(k+1)*t(k+1)
      sal(k) = sal(k) + c1(k+1)*sal(k+1)
    enddo
  else ! dt <= 0.0
    do k=1,nz
      u(k) = u0(k) ; v(k) = v0(k) ; t(k) = t0(k) ; sal(k) = s0(k)
      if (abs(u(k)) < underflow_vel) u(k) = 0.0
      if (abs(v(k)) < underflow_vel) v(k) = 0.0
    enddo
  endif
 
  if (present(s2)) then
    ! L2_to_Z2 = US%m_to_Z**2 * US%T_to_s**2
    l2_to_z2 = us%L_to_Z**2
    s2(1) = 0.0 ; s2(nz+1) = 0.0
    if (ks > 1) &
      s2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (l2_to_z2*i_dz_int(ks)**2)
    do k=ks+1,ke
      s2(k) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (l2_to_z2*i_dz_int(k)**2)
    enddo
    if (ke<nz) &
      s2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (l2_to_z2*i_dz_int(ke+1)**2)
  endif
 
  if (present(n2)) then
    n2(1) = 0.0 ; n2(nz+1) = 0.0
    if (ks > 1) &
      n2(ks) = max(0.0, i_dz_int(ks) * &
        (dbuoy_dt(ks) * (t0(ks-1)-t(ks)) + dbuoy_ds(ks) * (s0(ks-1)-sal(ks))))
    do k=ks+1,ke
      n2(k) = max(0.0, i_dz_int(k) * &
        (dbuoy_dt(k) * (t(k-1)-t(k)) + dbuoy_ds(k) * (sal(k-1)-sal(k))))
    enddo
    if (ke<nz) &
      n2(ke+1) = max(0.0, i_dz_int(ke+1) * &
        (dbuoy_dt(ke+1) * (t(ke)-t0(ke+1)) + dbuoy_ds(ke+1) * (sal(ke)-s0(ke+1))))
  endif
 
end subroutine calculate_projected_state
 
!> This subroutine calculates new, consistent estimates of TKE and kappa.
subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, &
                          nz, CS, GV, US, K_Q, tke, kappa, kappa_src, local_src)
  integer,               intent(in)    :: nz  !< The number of layers to work on.
  real, dimension(nz+1), intent(in)    :: N2  !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
  real, dimension(nz+1), intent(in)    :: S2  !< The squared shear at interfaces [T-2 ~> s-2].
  real, dimension(nz+1), intent(in)    :: kappa_in  !< The initial guess at the diffusivity
                                              !! [Z2 T-1 ~> m2 s-1].
  real, dimension(nz+1), intent(in)    :: dz_Int !< The thicknesses associated with interfaces
                                              !! [Z-1 ~> m-1].
  real, dimension(nz+1), intent(in)    :: I_L2_bdry !< The inverse of the squared distance to
                                              !! boundaries [Z-2 !> m-2].
  real, dimension(nz),   intent(in)    :: Idz !< The inverse grid spacing of layers [Z-1 ~> m-1].
  real,                  intent(in)    :: f2  !< The squared Coriolis parameter [T-2 ~> s-2].
  type(kappa_shear_cs),  pointer       :: CS  !< A pointer to this module's control structure.
  type(verticalgrid_type), intent(in)  :: GV  !< The ocean's vertical grid structure.
  type(unit_scale_type), intent(in)    :: US  !< A dimensional unit scaling type
  real, dimension(nz+1), intent(inout) :: K_Q !< The shear-driven diapycnal diffusivity divided by
                                              !! the turbulent kinetic energy per unit mass at
                                              !! interfaces [Z2 m-2 s2 T-1 ~> s].
  real, dimension(nz+1), intent(out)   :: tke !< The turbulent kinetic energy per unit mass at
                                              !! interfaces [Z2 T-2 ~> m2 s-2].
  real, dimension(nz+1), intent(out)   :: kappa  !< The diapycnal diffusivity at interfaces
                                              !! [Z2 T-1 ~> m2 s-1].
  real, dimension(nz+1), optional, &
                         intent(out)   :: kappa_src !< The source term for kappa [T-1].
  real, dimension(nz+1), optional, &
                         intent(out)   :: local_src !< The sum of all local sources for kappa,
                                              !! [T-1 ~> s-1].
!   This subroutine calculates new, consistent estimates of TKE and kappa.
 
  ! Local variables
  real, dimension(nz) :: &
    aQ, &       ! aQ is the coupling between adjacent interfaces in the TKE equations [Z T-1 ~> m s-1].
    dQdz        ! Half the partial derivative of TKE with depth [Z T-2 ~> m s-2].
  real, dimension(nz+1) :: &
    dK, &         ! The change in kappa [Z2 T-1 ~> m2 s-1].
    dQ, &         ! The change in TKE [Z2 T-2 ~> m2 s-2].
    cQ, cK, &     ! cQ and cK are the upward influences in the tridiagonal and
                  ! hexadiagonal solvers for the TKE and kappa equations [nondim].
    i_ld2, &      ! 1/Ld^2, where Ld is the effective decay length scale for kappa [Z-2 ~> m-2].
    tke_decay, &  ! The local TKE decay rate [T-1 ~> s-1].
    k_src, &      ! The source term in the kappa equation [T-1 ~> s-1].
    dqmdk, &      ! With Newton's method the change in dQ(k-1) due to dK(k) [T ~> s].
    dkdq, &       ! With Newton's method the change in dK(k) due to dQ(k) [T-1 ~> s-1].
    e1            ! The fractional change in a layer TKE due to a change in the
                  ! TKE of the layer above when all the kappas below are 0.
                  ! e1 is nondimensional, and 0 < e1 < 1.
  real :: tke_src       ! The net source of TKE due to mixing against the shear
                        ! and stratification [Z2 T-3 ~> m2 s-3].  (For convenience,
                        ! a term involving the non-dissipation of q0 is also
                        ! included here.)
  real :: bQ            ! The inverse of the pivot in the tridiagonal equations [T Z-1 ~> s m-1].
  real :: bK            ! The inverse of the pivot in the tridiagonal equations [Z-1 ~> m-1].
  real :: bQd1          ! A term in the denominator of bQ [Z T-1 ~> m s-1].
  real :: bKd1          ! A term in the denominator of bK [Z ~> m].
  real :: cQcomp, cKcomp ! 1 - cQ or 1 - cK in the tridiagonal equations.
  real :: c_s2          !   The coefficient for the decay of TKE due to
                        ! shear (i.e. proportional to |S|*tke), nondimensional.
  real :: c_n2          !   The coefficient for the decay of TKE due to
                        ! stratification (i.e. proportional to N*tke) [nondim].
  real :: Ri_crit       !   The critical shear Richardson number for shear-
                        ! driven mixing. The theoretical value is 0.25.
  real :: q0            !   The background level of TKE [Z2 T-2 ~> m2 s-2].
  real :: Ilambda2      ! 1.0 / CS%lambda**2 [nondim]
  real :: TKE_min       !   The minimum value of shear-driven TKE that can be
                        ! solved for [Z2 T-2 ~> m2 s-2].
  real :: kappa0        ! The background diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
  real :: kappa_trunc   ! Diffusivities smaller than this are rounded to 0 [Z2 T-1 ~> m2 s-1].
 
  real :: eden1, eden2, I_eden, ome  ! Variables used in calculating e1.
  real :: diffusive_src ! The diffusive source in the kappa equation [Z T-1 ~> m s-1].
  real :: chg_by_k0     ! The value of k_src that leads to an increase of
                        ! kappa_0 if only the diffusive term is a sink [T-1 ~> s-1].
 
  real :: kappa_mean    ! A mean value of kappa [Z2 T-1 ~> m2 s-1].
  real :: Newton_test   ! The value of relative error that will cause the next
                        ! iteration to use Newton's method.
  ! Temporary variables used in the Newton's method iterations.
  real :: decay_term_k  ! The decay term in the diffusivity equation
  real :: decay_term_Q  ! The decay term in the TKE equation - proportional to [T-1 ~> s-1]
  real :: I_Q           ! The inverse of TKE [s2 m-2]
  real :: kap_src
  real :: v1            ! A temporary variable proportional to [T-1 ~> s-1]
  real :: v2
  real :: tol_err        ! The tolerance for max_err that determines when to
                         ! stop iterating.
  real :: Newton_err     ! The tolerance for max_err that determines when to
                         ! start using Newton's method.  Empirically, an initial
                         ! value of about 0.2 seems to be most efficient.
  real, parameter :: roundoff = 1.0e-16 ! A negligible fractional change in TKE.
                         ! This could be larger but performance gains are small.
 
  logical :: tke_noflux_bottom_BC = .false. ! Specify the boundary conditions
  logical :: tke_noflux_top_BC = .false.    ! that are applied to the TKE eqns.
  logical :: do_Newton    ! If .true., use Newton's method for the next iteration.
  logical :: abort_Newton ! If .true., an Newton's method has encountered a 0
                          ! pivot, and should not have been used.
  logical :: was_Newton   ! The value of do_Newton before checking convergence.
  logical :: within_tolerance ! If .true., all points are within tolerance to
                          ! enable this subroutine to return.
  integer :: ks_src, ke_src ! The range indices that have nonzero k_src.
  integer :: ks_kappa, ke_kappa, ke_tke   ! The ranges of k-indices that are or
  integer :: ks_kappa_prev, ke_kappa_prev ! were being worked on.
  integer :: itt, k, k2
#ifdef DEBUG
  integer :: max_debug_itt ; parameter(max_debug_itt=20)
  real :: K_err_lin, Q_err_lin, TKE_src_norm
  real, dimension(nz+1) :: &
    I_Ld2_debug, & ! A separate version of I_Ld2 for debugging [Z-2 ~> m-2].
    kappa_prev, & ! The value of kappa at the start of the current iteration [Z2 T-1 ~> m2 s-1].
    TKE_prev   ! The value of TKE at the start of the current iteration [Z2 T-2 ~> m2 s-2].
  real, dimension(nz+1,1:max_debug_itt) :: &
    tke_it1, kappa_it1, kprev_it1, &  ! Various values from each iteration.
    dkappa_it1, K_Q_it1, d_dkappa_it1, dkappa_norm_it1
  integer :: it2
#endif
 
  c_n2 = cs%C_N**2 ; c_s2 = cs%C_S**2
  q0 = cs%TKE_bg ; kappa0 = cs%kappa_0
  tke_min = max(cs%TKE_bg, 1.0e-20*us%m_to_Z**2*us%T_to_s**2)
  ri_crit = cs%Rino_crit
  ilambda2 = 1.0 / cs%lambda**2
  kappa_trunc = 0.01*kappa0  ! ### CHANGE THIS HARD-WIRING LATER?
  do_newton = .false. ; abort_newton = .false.
  tol_err = cs%kappa_tol_err
  newton_err = 0.2     ! This initial value may be automatically reduced later.
 
  ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
 
  ke_src = 0 ; ks_src = nz+1
  do k=2,nz
    if (n2(k) < ri_crit * s2(k)) then ! Equivalent to Ri < Ri_crit.
!       Ri = N2(K) / S2(K)
!       k_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
!                  ((Ri_crit - Ri) / (Ri_crit + CS%FRi_curvature*Ri))
      k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
                 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
      ke_src = k
      if (ks_src > k) ks_src = k
    else
      k_src(k) = 0.0
    endif
  enddo
 
  ! If there is no source anywhere, return kappa(K) = 0.
  if (ks_src > ke_src) then
    do k=1,nz+1
      kappa(k) = 0.0 ; k_q(k) = 0.0 ; tke(k) = tke_min
    enddo
    if (present(kappa_src)) then ; do k=1,nz+1 ; kappa_src(k) = 0.0 ; enddo ; endif
    if (present(local_src)) then ; do k=1,nz+1 ; local_src(k) = 0.0 ; enddo ; endif
    return
  endif
 
  do k=1,nz+1
    kappa(k) = kappa_in(k)
!     TKE_decay(K) = c_n*sqrt(N2(K)) + c_s*sqrt(S2(K)) ! The expression in JHL.
    tke_decay(k) = sqrt(c_n2*n2(k) + c_s2*s2(k))
    if ((kappa(k) > 0.0) .and. (k_q(k) > 0.0)) then
      tke(k) = kappa(k) / k_q(k)
    else
      tke(k) = tke_min
    endif
  enddo
  ! Apply boundary conditions to kappa.
  kappa(1) = 0.0 ; kappa(nz+1) = 0.0
 
  ! Calculate the term (e1) that allows changes in TKE to be calculated quickly
  ! below the deepest nonzero value of kappa.  If kappa = 0, below interface
  ! k-1, the final changes in TKE are related by dQ(K+1) = e1(K+1)*dQ(K).
  eden2 = kappa0 * idz(nz)
  if (tke_noflux_bottom_bc) then
    eden1 = dz_int(nz+1)*tke_decay(nz+1)
    i_eden = 1.0 / (eden2 + eden1)
    e1(nz+1) = eden2 * i_eden ; ome = eden1 * i_eden
  else
    e1(nz+1) = 0.0 ; ome = 1.0
  endif
  do k=nz,2,-1
    eden1 = dz_int(k)*tke_decay(k) + ome * eden2
    eden2 = kappa0 * idz(k-1)
    i_eden = 1.0 / (eden2 + eden1)
    e1(k) = eden2 * i_eden ; ome = eden1 * i_eden ! = 1-e1
  enddo
  e1(1) = 0.0
 
 
  ! Iterate here to convergence to within some tolerance of order tol_err.
  do itt=1,cs%max_RiNo_it
 
  ! ----------------------------------------------------
  ! Calculate TKE
  ! ----------------------------------------------------
 
#ifdef DEBUG
    do k=1,nz+1 ; kappa_prev(k) = kappa(k) ; tke_prev(k) = tke(k) ; enddo
#endif
 
    if (.not.do_newton) then
      !   Use separate steps of the TKE and kappa equations, that are
      ! explicit in the nonlinear source terms, implicit in a linearized
      ! version of the nonlinear sink terms, and implicit in the linear
      ! terms.
 
      ke_tke = max(ke_kappa,ke_kappa_prev)+1
      ! aQ is the coupling between adjacent interfaces [Z T-1 ~> m s-1].
      do k=1,min(ke_tke,nz)
        aq(k) = (0.5*(kappa(k)+kappa(k+1)) + kappa0) * idz(k)
      enddo
      dq(1) = -tke(1)
      if (tke_noflux_top_bc) then
        tke_src = kappa0*s2(1) + q0 * tke_decay(1) ! Uses that kappa(1) = 0
        bqd1 = dz_int(1) * tke_decay(1)
        bq = 1.0 / (bqd1 +  aq(1))
        tke(1) = bq * (dz_int(1)*tke_src)
        cq(2) = aq(1) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
      else
        tke(1) = q0 ; cq(2) = 0.0 ; cqcomp = 1.0
      endif
      do k=2,ke_tke-1
        dq(k) = -tke(k)
        tke_src = (kappa(k) + kappa0)*s2(k) + q0*tke_decay(k)
        bqd1 = dz_int(k)*(tke_decay(k) + n2(k)*k_q(k)) + cqcomp*aq(k-1)
        bq = 1.0 / (bqd1 + aq(k))
        tke(k) = bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1))
        cq(k+1) = aq(k) * bq ; cqcomp = bqd1 * bq ! = 1 - cQ
      enddo
      if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_bc)) then
        tke(nz+1) = tke_min
        dq(nz+1) = 0.0
      else
        k = ke_tke
        tke_src = kappa0*s2(k) + q0*tke_decay(k) ! Uses that kappa(ke_tke) = 0
        if (k == nz+1) then
          dq(k) = -tke(k)
          bq = 1.0 / (dz_int(k)*tke_decay(k) + cqcomp*aq(k-1))
          tke(k) = max(tke_min, bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)))
          dq(k) = tke(k) + dq(k)
        else
          bq = 1.0 / ((dz_int(k)*tke_decay(k) + cqcomp*aq(k-1)) + aq(k))
          cq(k+1) = aq(k) * bq
          ! Account for all changes deeper in the water column.
          dq(k) = -tke(k)
          tke(k) = max((bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)) + &
                        cq(k+1)*(tke(k+1) - e1(k+1)*tke(k))) / (1.0 - cq(k+1)*e1(k+1)), tke_min)
          dq(k) = tke(k) + dq(k)
 
          ! Adjust TKE deeper in the water column in case ke_tke increases.
          ! This might not be strictly necessary?
          do k=ke_tke+1,nz+1
            dq(k) = e1(k)*dq(k-1)
            tke(k) = max(tke(k) + dq(k), tke_min)
            if (abs(dq(k)) < roundoff*tke(k)) exit
          enddo
          do k2=k+1,nz
            if (dq(k2) == 0.0) exit
            dq(k2) = 0.0
          enddo
        endif
      endif
      do k=ke_tke-1,1,-1
        tke(k) = max(tke(k) + cq(k+1)*tke(k+1), tke_min)
        dq(k) = tke(k) + dq(k)
      enddo
 
  ! ----------------------------------------------------
  ! Calculate kappa, here defined at interfaces.
  ! ----------------------------------------------------
 
      ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
 
      dk(1) = 0.0 ! kappa takes boundary values of 0.
      ck(2) = 0.0 ; ckcomp = 1.0
      if (itt == 1) then ; dO k=2,nz
        i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
      enddo ; endif
      do k=2,nz
        dk(k) = -kappa(k)
        if (itt>1) &
          i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
        bkd1 = dz_int(k)*i_ld2(k) + ckcomp*idz(k-1)
        bk = 1.0 / (bkd1 + idz(k))
 
        kappa(k) = bk * (idz(k-1)*kappa(k-1) + dz_int(k) * k_src(k))
        ck(k+1) = idz(k) * bk ; ckcomp = bkd1 * bk ! = 1 - cK(K+1)
 
        ! Neglect values that are smaller than kappa_trunc.
        if (kappa(k) < ckcomp*kappa_trunc) then
          kappa(k) = 0.0
          if (k > ke_src) then ; ke_kappa = k-1 ; k_q(k) = 0.0 ; exit ; endif
        elseif (kappa(k) < 2.0*ckcomp*kappa_trunc) then
          kappa(k) = 2.0 * (kappa(k) - ckcomp*kappa_trunc)
        endif
      enddo
      k_q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
      dk(ke_kappa) = dk(ke_kappa) + kappa(ke_kappa)
      do k=ke_kappa+2,ke_kappa_prev
        dk(k) = -kappa(k) ; kappa(k) = 0.0 ; k_q(k) = 0.0
      enddo
      do k=ke_kappa-1,2,-1
        kappa(k) = kappa(k) + ck(k+1)*kappa(k+1)
        ! Neglect values that are smaller than kappa_trunc.
        if (kappa(k) <= kappa_trunc) then
          kappa(k) = 0.0
          if (k < ks_src) then ; ks_kappa = k+1 ; k_q(k) = 0.0 ; exit ; endif
        elseif (kappa(k) < 2.0*kappa_trunc) then
          kappa(k) = 2.0 * (kappa(k) - kappa_trunc)
        endif
 
        dk(k) = dk(k) + kappa(k)
        k_q(k) = kappa(k) / tke(k)
      enddo
      do k=ks_kappa_prev,ks_kappa-2 ; kappa(k) = 0.0 ; k_q(k) = 0.0 ; enddo
 
    else ! do_Newton is .true.
!   Once the solutions are close enough, use a Newton's method solver of the
!  whole system to accelerate convergence.
      ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
      ks_kappa = 2
      dk(1) = 0.0 ; ck(2) = 0.0 ; ckcomp = 1.0 ; dkdq(1) = 0.0
      aq(1) = (0.5*(kappa(1)+kappa(2))+kappa0) * idz(1)
      dqdz(1) = 0.5*(tke(1) - tke(2))*idz(1)
      if (tke_noflux_top_bc) then
        tke_src = dz_int(1) * (kappa0*s2(1) - (tke(1) - q0)*tke_decay(1)) - &
                  aq(1) * (tke(1) - tke(2))
 
        bq = 1.0 / (aq(1) + dz_int(1)*tke_decay(1))
        cq(2) = aq(1) * bq
        cqcomp = (dz_int(1)*tke_decay(1)) * bq ! = 1 - cQ(2)
        dqmdk(2) = -dqdz(1) * bq
        dq(1) = bq * tke_src
      else
        dq(1) = 0.0 ; cq(2) = 0.0 ; cqcomp = 1.0 ; dqmdk(2) = 0.0
      endif
      do k=2,nz
        i_q = 1.0 / tke(k)
        i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
 
        kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa(k)) + &
                            idz(k-1)*(kappa(k-1)-kappa(k)) - idz(k)*(kappa(k)-kappa(k+1))
 
        ! Ensure that the pivot is always positive, and that 0 <= cK <= 1.
        ! Otherwise do not use Newton's method.
        decay_term_k = -idz(k-1)*dqmdk(k)*dkdq(k-1) + dz_int(k)*i_ld2(k)
        if (decay_term_k < 0.0) then ; abort_newton = .true. ; exit ; endif
        bk = 1.0 / (idz(k) + idz(k-1)*ckcomp + decay_term_k)
 
        ck(k+1) = bk * idz(k)
        ckcomp = bk * (idz(k-1)*ckcomp + decay_term_k) ! = 1-cK(K+1)
        !### The following expression appears to be dimensionally inconsistent in length. -RWH
        dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
                      us%m_to_Z*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
        !  I think that the second term needs to be multiplied by dz_Int(K):
        !             dz_Int(K)*(N2(K)*Ilambda2 + f2) * I_Q**2 * kappa(K) )
        dk(k) = bk * (kap_src + idz(k-1)*dk(k-1) + idz(k-1)*dkdq(k-1)*dq(k-1))
 
        ! Truncate away negligibly small values of kappa.
        if (dk(k) <= ckcomp*(kappa_trunc - kappa(k))) then
          dk(k) = -ckcomp*kappa(k)
!         if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
        elseif (dk(k) < ckcomp*(2.0*kappa_trunc - kappa(k))) then
          dk(k) = 2.0 * dk(k) - ckcomp*(2.0*kappa_trunc - kappa(k))
        endif
 
        ! Solve for dQ(K)...
        aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
        dqdz(k) = 0.5*(tke(k) - tke(k+1))*idz(k)
        tke_src = dz_int(k) * (((kappa(k) + kappa0)*s2(k) - kappa(k)*n2(k)) - &
                                  (tke(k) - q0)*tke_decay(k)) - &
                  (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k)))
        v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
        v2 = (v1*dqmdk(k) + dqdz(k-1)*ck(k)) + &
             ((dqdz(k-1) - dqdz(k)) + dz_int(k)*(s2(k) - n2(k)))
 
        ! Ensure that the pivot is always positive, and that 0 <= cQ <= 1.
        ! Otherwise do not use Newton's method.
        decay_term_q = dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k) - v2*dkdq(k)
        if (decay_term_q < 0.0) then ; abort_newton = .true. ; exit ; endif
        bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
 
        cq(k+1) = aq(k) * bq
        cqcomp = (cqcomp*aq(k-1) + decay_term_q) * bq
        dqmdk(k+1) = (v2 * ck(k+1) - dqdz(k)) * bq
 
        ! Ensure that TKE+dQ will not drop below 0.5*TKE.
        dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + &
                          (v2 * dk(k) + tke_src)), cqcomp*(-0.5*tke(k)))
 
        ! Check whether the next layer will be affected by any nonzero kappas.
        if ((itt > 1) .and. (k > ke_src) .and. (dk(k) == 0.0) .and. &
            ((kappa(k) + kappa(k+1)) == 0.0)) then
        ! Could also do  .and. (bQ*abs(tke_src) < roundoff*TKE(K)) then
          ke_kappa = k-1 ; exit
        endif
      enddo
      if ((ke_kappa == nz) .and. (.not. abort_newton)) then
        dk(nz+1) = 0.0 ; dkdq(nz+1) = 0.0
        if (tke_noflux_bottom_bc) then
          k = nz+1
          tke_src = dz_int(k) * (kappa0*s2(k) - (tke(k) - q0)*tke_decay(k)) + &
                    aq(k-1) * (tke(k-1) - tke(k))
 
          v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
          decay_term_q = max(0.0, dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k))
          if (decay_term_q < 0.0) then
            abort_newton = .true.
          else
            bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
          ! Ensure that TKE+dQ will not drop below 0.5*TKE.
            dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + tke_src), -0.5*tke(k))
            tke(k) = max(tke(k) + dq(k), tke_min)
          endif
        else
          dq(nz+1) = 0.0
        endif
      elseif (.not. abort_newton) then
        ! Alter the first-guess determination of dQ(K).
        dq(ke_kappa+1) = dq(ke_kappa+1) / (1.0 - cq(ke_kappa+2)*e1(ke_kappa+2))
        tke(ke_kappa+1) = max(tke(ke_kappa+1) + dq(ke_kappa+1), tke_min)
        do k=ke_kappa+2,nz+1
#ifdef DEBUG
          if (k < nz+1) then
          ! Ignore this source?
            aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
            tke_src_norm = (dz_int(k) * (kappa0*s2(k) - (tke(k)-q0)*tke_decay(k)) - &
                           (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k))) ) / &
                           (aq(k) + (aq(k-1) + dz_int(k)*tke_decay(k)))
          endif
#endif
          dk(k) = 0.0
        ! Ensure that TKE+dQ will not drop below 0.5*TKE.
          dq(k) = max(e1(k)*dq(k-1),-0.5*tke(k))
          tke(k) = max(tke(k) + dq(k), tke_min)
          if (abs(dq(k)) < roundoff*tke(k)) exit
        enddo
#ifdef DEBUG
        do k2=k+1,ke_kappa_prev+1 ; dq(k2) = 0.0 ; dk(k2) = 0.0 ; enddo
        do k=k2,nz+1 ; if (dq(k) == 0.0) exit ; dq(k) = 0.0 ; dk(k) = 0.0 ; enddo
#endif
      endif
      if (.not. abort_newton) then
        do k=ke_kappa,2,-1
          ! Ensure that TKE+dQ will not drop below 0.5*TKE.
          dq(k) = max(dq(k) + (cq(k+1)*dq(k+1) + dqmdk(k+1) * dk(k+1)), -0.5*tke(k))
          tke(k) = max(tke(k) + dq(k), tke_min)
          dk(k) = dk(k) + (ck(k+1)*dk(k+1) + dkdq(k) * dq(k))
          ! Truncate away negligibly small values of kappa.
          if (dk(k) <= kappa_trunc - kappa(k)) then
            dk(k) = -kappa(k)
            kappa(k) = 0.0
            if ((k < ks_src) .and. (k+1 > ks_kappa)) ks_kappa = k+1
          elseif (dk(k) < 2.0*kappa_trunc - kappa(k)) then
            dk(k) =  2.0*dk(k) - (2.0*kappa_trunc - kappa(k))
            kappa(k) = max(kappa(k) + dk(k), 0.0) ! The max is for paranoia.
            if (k<=ks_kappa) ks_kappa = 2
          else
            kappa(k) = kappa(k) + dk(k)
            if (k<=ks_kappa) ks_kappa = 2
          endif
        enddo
        dq(1) = max(dq(1) + cq(2)*dq(2) + dqmdk(2) * dk(2), tke_min - tke(1))
        tke(1) = max(tke(1) + dq(1), tke_min)
        dk(1) = 0.0
      endif
 
#ifdef DEBUG
      ! Check these solutions for consistency.
      !  The unit conversions here have not been carefully tested.
      do k=2,nz
        ! In these equations, K_err_lin and Q_err_lin should be at round-off levels
        ! compared with the dominant terms, perhaps, dz_Int*I_Ld2*kappa and
        ! dz_Int*TKE_decay*TKE.  The exception is where, either 1) the decay term has been
        ! been increased to ensure a positive pivot, or 2) negative TKEs have been
        ! truncated, or 3) small or negative kappas have been rounded toward 0.
        i_q = 1.0 / tke(k)
        i_ld2_debug(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
 
        kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa_prev(k)) + &
                            (idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
                             idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
        !### The last line of the following appears to be dimensionally inconsistent with the first two.
        !  I think that the term on the last line needs to be multiplied by dz_Int(K).
        k_err_lin = -idz(k-1)*(dk(k-1)-dk(k)) + idz(k)*(dk(k)-dk(k+1)) + &
                     dz_int(k)*i_ld2_debug(k)*dk(k) - kap_src - &
                     us%m_to_Z*(n2(k)*ilambda2 + f2)*i_q**2*kappa_prev(k) * dq(k)
 
        tke_src = dz_int(k) * ((kappa_prev(k) + kappa0)*s2(k) - &
                     kappa_prev(k)*n2(k) - (tke_prev(k) - q0)*tke_decay(k)) - &
                  (aq(k) * (tke_prev(k) - tke_prev(k+1)) -  aq(k-1) * (tke_prev(k-1) - tke_prev(k)))
        q_err_lin = tke_src + (aq(k-1) * (dq(k-1)-dq(k)) - aq(k) * (dq(k)-dq(k+1))) - &
                    0.5*(tke_prev(k)-tke_prev(k+1))*idz(k)  * (dk(k) + dk(k+1)) - &
                    0.5*(tke_prev(k)-tke_prev(k-1))*idz(k-1)* (dk(k-1) + dk(k)) + &
                    dz_int(k) * (dk(k) * (s2(k) - n2(k)) - dq(k)*tke_decay(k))
      enddo
#endif
    endif  ! End of the Newton's method solver.
 
    ! Test kappa for convergence...
    if ((tol_err < newton_err) .and. (.not.abort_newton)) then
      !  A lower tolerance is used to switch to Newton's method than to switch back.
      newton_test = newton_err ; if (do_newton) newton_test = 2.0*newton_err
      was_newton = do_newton
      within_tolerance = .true. ; do_newton = .true.
      do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
        kappa_mean = kappa0 + (kappa(k) - 0.5*dk(k))
        if (abs(dk(k)) > newton_test * kappa_mean) then
          if (do_newton) abort_newton = .true.
          within_tolerance = .false. ; do_newton = .false. ; exit
        elseif (abs(dk(k)) > tol_err * kappa_mean) then
          within_tolerance = .false. ; if (.not.do_newton) exit
        endif
        if (abs(dq(k)) > newton_test*(tke(k) - 0.5*dq(k))) then
          if (do_newton) abort_newton = .true.
          do_newton = .false. ; if (.not.within_tolerance) exit
        endif
      enddo
 
    else  ! Newton's method will not be used again, so no need to check.
      within_tolerance = .true.
      do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
        if (abs(dk(k)) > tol_err * (kappa0 + (kappa(k) - 0.5*dk(k)))) then
          within_tolerance = .false. ;  exit
        endif
      enddo
    endif
 
    if (abort_newton) then
      do_newton = .false. ; abort_newton = .false.
      ! We went to Newton too quickly last time, so restrict the tolerance.
      newton_err = 0.5*newton_err
      ke_kappa_prev = nz
      do k=2,nz ; k_q(k) = kappa(k) / max(tke(k), tke_min) ; enddo
    endif
 
#ifdef DEBUG
    if (itt <= max_debug_itt) then
      do k=1,nz+1
        kprev_it1(k,itt) = kappa_prev(k)
        kappa_it1(k,itt) = kappa(k) ; tke_it1(k,itt) = tke(k)
        dkappa_it1(k,itt) = kappa(k) - kappa_prev(k)
        dkappa_norm_it1(k,itt) = (kappa(k) - kappa_prev(k)) / &
            (kappa0 + 0.5*(kappa(k) + kappa_prev(k)))
        k_q_it1(k,itt) = kappa(k) / max(tke(k),tke_min)
        d_dkappa_it1(k,itt) = 0.0
        if (itt > 1) then ; if (abs(dkappa_it1(k,itt-1)) > 1e-20*us%m2_s_to_Z2_T) &
            d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
        endif
      enddo
    endif
#endif
 
    if (within_tolerance) exit
 
  enddo
 
#ifdef DEBUG
  do it2=itt+1,max_debug_itt ; do k=1,nz+1
    kprev_it1(k,it2) = 0.0 ; kappa_it1(k,it2) = 0.0 ; tke_it1(k,it2) = 0.0
    dkappa_it1(k,it2) = 0.0 ; k_q_it1(k,it2) = 0.0 ; d_dkappa_it1(k,it2) = 0.0
  enddo ; enddo
#endif
 
  if (do_newton) then  ! K_Q needs to be calculated.
    do k=1,ks_kappa-1 ;  k_q(k) = 0.0 ; enddo
    do k=ks_kappa,ke_kappa ; k_q(k) = kappa(k) / tke(k) ; enddo
    do k=ke_kappa+1,nz+1 ; k_q(k) = 0.0 ; enddo
  endif
 
  if (present(local_src)) then
    local_src(1) = 0.0 ; local_src(nz+1) = 0.0
    do k=2,nz
      diffusive_src = idz(k-1)*(kappa(k-1)-kappa(k)) + idz(k)*(kappa(k+1)-kappa(k))
      chg_by_k0 = kappa0 * ((idz(k-1)+idz(k)) / dz_int(k) + i_ld2(k))
      if (diffusive_src <= 0.0) then
        local_src(k) = k_src(k) + chg_by_k0
      else
        local_src(k) = (k_src(k) + chg_by_k0) + diffusive_src / dz_int(k)
      endif
    enddo
  endif
  if (present(kappa_src)) then
    kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
    do k=2,nz
      kappa_src(k) = k_src(k)
    enddo
  endif
 
end subroutine find_kappa_tke
 
!> This subroutineinitializesthe parameters that regulate shear-driven mixing
function kappa_shear_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type),         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(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(kappa_shear_cs),    pointer       :: cs   !< A pointer that is set to point to the control
                                                 !! structure for this module
  logical :: kappa_shear_init !< True if module is to be used, False otherwise
 
  ! Local variables
  logical :: merge_mixedlayer
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_kappa_shear"  ! This module's name.
  real :: kd_normal ! The KD of the main model, read here only as a parameter
                    ! for setting the default of KD_SMOOTH
 
  if (associated(cs)) then
    call mom_error(warning, "kappa_shear_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)
 
  !   The Jackson-Hallberg-Legg shear mixing parameterization uses the following
  ! 6 nondimensional coefficients.  That paper gives 3 best fit parameter sets.
  !    Ri_Crit  Rate    FRi_Curv  K_buoy  TKE_N  TKE_Shear
  ! p1: 0.25    0.089    -0.97     0.82    0.24    0.14
  ! p2: 0.30    0.085    -0.94     0.86    0.26    0.13
  ! p3: 0.35    0.088    -0.89     0.81    0.28    0.12
  !   Future research will reveal how these should be modified to take
  ! subgridscale inhomogeneity into account.
 
! Set default, read and log parameters
  call log_version(param_file, mdl, version, &
    "Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008")
  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, &
                 "If true, use the Jackson-Hallberg-Legg (JPO 2008) "//&
                 "shear mixing parameterization.", default=.false.)
  call get_param(param_file, mdl, "VERTEX_SHEAR", cs%KS_at_vertex, &
                 "If true, do the calculations of the shear-driven mixing "//&
                 "at the cell vertices (i.e., the vorticity points).", &
                 default=.false.)
  call get_param(param_file, mdl, "RINO_CRIT", cs%RiNo_crit, &
                 "The critical Richardson number for shear mixing.", &
                 units="nondim", default=0.25)
  call get_param(param_file, mdl, "SHEARMIX_RATE", cs%Shearmix_rate, &
                 "A nondimensional rate scale for shear-driven entrainment. "//&
                 "Jackson et al find values in the range of 0.085-0.089.", &
                 units="nondim", default=0.089)
  call get_param(param_file, mdl, "MAX_RINO_IT", cs%max_RiNo_it, &
                 "The maximum number of iterations that may be used to "//&
                 "estimate the Richardson number driven mixing.", &
                 units="nondim", default=50)
  call get_param(param_file, mdl, "KD", kd_normal, default=1.0e-7, do_not_log=.true.)
  call get_param(param_file, mdl, "KD_KAPPA_SHEAR_0", cs%kappa_0, &
                 "The background diffusivity that is used to smooth the "//&
                 "density and shear profiles before solving for the "//&
                 "diffusivities. Defaults to value of KD.", &
                 units="m2 s-1", default=kd_normal, scale=us%m2_s_to_Z2_T)
  call get_param(param_file, mdl, "FRI_CURVATURE", cs%FRi_curvature, &
                 "The nondimensional curvature of the function of the "//&
                 "Richardson number in the kappa source term in the "//&
                 "Jackson et al. scheme.", units="nondim", default=-0.97)
  call get_param(param_file, mdl, "TKE_N_DECAY_CONST", cs%C_N, &
                 "The coefficient for the decay of TKE due to "//&
                 "stratification (i.e. proportional to N*tke). "//&
                 "The values found by Jackson et al. are 0.24-0.28.", &
                 units="nondim", default=0.24)
!  call get_param(param_file, mdl, "LAYER_KAPPA_STAGGER", CS%layer_stagger, &
!                 default=.false.)
  call get_param(param_file, mdl, "TKE_SHEAR_DECAY_CONST", cs%C_S, &
                 "The coefficient for the decay of TKE due to shear (i.e. "//&
                 "proportional to |S|*tke). The values found by Jackson "//&
                 "et al. are 0.14-0.12.", units="nondim", default=0.14)
  call get_param(param_file, mdl, "KAPPA_BUOY_SCALE_COEF", cs%lambda, &
                 "The coefficient for the buoyancy length scale in the "//&
                 "kappa equation.  The values found by Jackson et al. are "//&
                 "in the range of 0.81-0.86.", units="nondim", default=0.82)
  call get_param(param_file, mdl, "KAPPA_N_OVER_S_SCALE_COEF2", cs%lambda2_N_S, &
                 "The square of the ratio of the coefficients of the "//&
                 "buoyancy and shear scales in the diffusivity equation, "//&
                 "Set this to 0 (the default) to eliminate the shear scale. "//&
                 "This is only used if USE_JACKSON_PARAM is true.", &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "KAPPA_SHEAR_TOL_ERR", cs%kappa_tol_err, &
                 "The fractional error in kappa that is tolerated. "//&
                 "Iteration stops when changes between subsequent "//&
                 "iterations are smaller than this everywhere in a "//&
                 "column.  The peak diffusivities usually converge most "//&
                 "rapidly, and have much smaller errors than this.", &
                 units="nondim", default=0.1)
  call get_param(param_file, mdl, "TKE_BACKGROUND", cs%TKE_bg, &
                 "A background level of TKE used in the first iteration "//&
                 "of the kappa equation.  TKE_BACKGROUND could be 0.", &
                 units="m2 s-2", default=0.0, scale=us%m_to_Z**2*us%T_to_s**2)
  call get_param(param_file, mdl, "KAPPA_SHEAR_ELIM_MASSLESS", cs%eliminate_massless, &
                 "If true, massless layers are merged with neighboring "//&
                 "massive layers in this calculation.  The default is "//&
                 "true and I can think of no good reason why it should "//&
                 "be false. This is only used if USE_JACKSON_PARAM is true.", &
                 default=.true.)
  call get_param(param_file, mdl, "MAX_KAPPA_SHEAR_IT", cs%max_KS_it, &
                 "The maximum number of iterations that may be used to "//&
                 "estimate the time-averaged diffusivity.", units="nondim", &
                 default=13)
  call get_param(param_file, mdl, "PRANDTL_TURB", cs%Prandtl_turb, &
                 "The turbulent Prandtl number applied to shear "//&
                 "instability.", units="nondim", default=1.0, do_not_log=.true.)
  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
                 "A negligibly small velocity magnitude below which velocity "//&
                 "components are set to 0.  A reasonable value might be "//&
                 "1e-30 m/s, which is less than an Angstrom divided by "//&
                 "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
  call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", cs%debug, &
                 "If true, write debugging data for the kappa-shear code. \n"//&
                 "Caution: this option is _very_ verbose and should only "//&
                 "be used in single-column mode!", &
                 default=.false., debuggingparam=.true.)
 
!    id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear', grain=CLOCK_ROUTINE)
!    id_clock_avg = cpu_clock_id('Ocean KS avg', grain=CLOCK_ROUTINE)
!    id_clock_project = cpu_clock_id('Ocean KS project', grain=CLOCK_ROUTINE)
!    id_clock_setup = cpu_clock_id('Ocean KS setup', grain=CLOCK_ROUTINE)
 
  cs%nkml = 1
  if (gv%nkml>0) then
    call get_param(param_file, mdl, "KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
                 "If true, combine the mixed layers together before "//&
                 "solving the kappa-shear equations.", default=.true.)
    if (merge_mixedlayer) cs%nkml = gv%nkml
  endif
 
! Forego remainder of initialization if not using this scheme
  if (.not. kappa_shear_init) return
 
  cs%diag => diag
 
  cs%id_Kd_shear = register_diag_field('ocean_model','Kd_shear',diag%axesTi,time, &
      'Shear-driven Diapycnal Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
  cs%id_TKE = register_diag_field('ocean_model','TKE_shear',diag%axesTi,time, &
      'Shear-driven Turbulent Kinetic Energy', 'm2 s-2', conversion=us%Z_to_m**2*us%s_to_T**2)
#ifdef ADD_DIAGNOSTICS
  cs%id_ILd2 = register_diag_field('ocean_model','ILd2_shear',diag%axesTi,time, &
      'Inverse kappa decay scale at interfaces', 'm-2', conversion=us%m_to_Z**2)
  cs%id_dz_Int = register_diag_field('ocean_model','dz_Int_shear',diag%axesTi,time, &
      'Finite volume thickness of interfaces', 'm', conversion=us%Z_to_m)
#endif
 
end function kappa_shear_init
 
!> This function indicates to other modules whether the Jackson et al shear mixing
!! parameterization will be used without needing to duplicate the log entry.
logical function kappa_shear_is_used(param_file)
  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
! Reads the parameter "USE_JACKSON_PARAM" and returns state.
  character(len=40)  :: mdl = "MOM_kappa_shear"  ! This module's name.
 
  call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_is_used, &
                 default=.false., do_not_log=.true.)
end function kappa_shear_is_used
 
!> This function indicates to other modules whether the Jackson et al shear mixing
!! parameterization will be used without needing to duplicate the log entry.
logical function kappa_shear_at_vertex(param_file)
  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
! Reads the parameter "USE_JACKSON_PARAM" and returns state.
  character(len=40)  :: mdl = "MOM_kappa_shear"  ! This module's name.
 
  logical :: do_kappa_shear
 
  call get_param(param_file, mdl, "USE_JACKSON_PARAM", do_kappa_shear, &
                 default=.false., do_not_log=.true.)
  kappa_shear_at_vertex = .false.
  if (do_kappa_shear) &
    call get_param(param_file, mdl, "VERTEX_SHEAR", kappa_shear_at_vertex, &
                 "If true, do the calculations of the shear-driven mixing "//&
                 "at the cell vertices (i.e., the vorticity points).", &
                 default=.false., do_not_log=.true.)
 
end function kappa_shear_at_vertex
 
!> \namespace mom_kappa_shear
!!
!! By Laura Jackson and Robert Hallberg, 2006-2008
!!
!!   This file contains the subroutines that determine the diapycnal
!! diffusivity driven by resolved shears, as specified by the
!! parameterizations described in Jackson and Hallberg (JPO, 2008).
!!
!!   The technique by which the 6 equations (for kappa, TKE, u, v, T,
!! and S) are solved simultaneously has been dramatically revised
!! from the previous version. The previous version was not converging
!! in some cases, especially near the surface mixed layer, while the
!! revised version does.  The revised version solves for kappa and
!! TKE with shear and stratification fixed, then marches the density
!! and velocities forward with an adaptive (and aggressive) time step
!! in a predictor-corrector-corrector emulation of a trapezoidal
!! scheme.  Run-time-settable parameters determine the tolerence to
!! which the kappa and TKE equations are solved and the minimum time
!! step that can be taken.
 
end module mom_kappa_shear