MOM_set_diffusivity.F90

Go to the documentation of this file.

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, September 1997 - June 2007                     *
!*                                                                     *
!*    This file contains the subroutines that sets the diapycnal       *
!*  diffusivity, perhaps adding up pieces that are calculated in other *
!*  files and passed in via the vertvisc type argument.                *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, buoy, ustar, T, S, Kd, ea, eb, etc.   *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_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 : diag_ctrl, time_type
use mom_diag_mediator,       only : safe_alloc_ptr, post_data, register_diag_field
use mom_diag_to_z,           only : diag_to_z_cs, register_zint_diag, calc_zint_diags
use mom_debugging,           only : hchksum, uvchksum
use mom_eos,                 only : calculate_density, calculate_density_derivs
use mom_error_handler,       only : mom_error, is_root_pe, fatal, warning, note
use mom_error_handler,       only : calltree_showquery
use mom_error_handler,       only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,         only : get_param, log_param, log_version, param_file_type
use mom_forcing_type,        only : forcing, optics_type
use mom_grid,                only : ocean_grid_type
use mom_internal_tides,      only : int_tide_cs, get_lowmode_loss
use mom_intrinsic_functions, only : invcosh
use mom_io,                  only : slasher, vardesc, var_desc
use mom_kappa_shear,         only : calculate_kappa_shear, kappa_shear_init, kappa_shear_cs
use mom_cvmix_shear,         only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_cs
use mom_string_functions,    only : uppercase
use mom_thickness_diffuse,   only : vert_fill_ts
use mom_variables,           only : thermo_var_ptrs, vertvisc_type, p3d
use mom_verticalgrid, only : verticalgrid_type

use user_change_diffusivity, only : user_change_diff, user_change_diff_init
use user_change_diffusivity, only : user_change_diff_end, user_change_diff_cs

use fms_mod,                 only : read_data

implicit none ; private

#include <MOM_memory.h>

public set_diffusivity
public set_bbl_tke
public set_diffusivity_init
public set_diffusivity_end

type, public :: set_diffusivity_cs ; private
  logical :: debug           ! If true, write verbose checksums for debugging.

  logical :: bulkmixedlayer  ! If true, a refined bulk mixed layer is used with
                             ! GV%nk_rho_varies variable density mixed & buffer
                             ! layers.
  real    :: fluxri_max      ! The flux Richardson number where the stratification is
                             ! large enough that N2 > omega2.  The full expression for
                             ! the Flux Richardson number is usually
                             ! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2.
  logical :: henyey_igw_background  ! If true, use a simplified variant of the
                             ! Henyey et al, JGR (1986) latitudinal scaling for
                             ! the background diapycnal diffusivity, which gives
                             ! a marked decrease in the diffusivity near the
                             ! equator.  The simplification here is to assume
                             ! that the in-situ stratification is the same as
                             ! the reference stratificaiton.
  logical :: henyey_igw_background_new ! same as Henyey_IGW_background
                             ! but incorporate the effect of
                             ! stratification on TKE dissipation,
                             !
                             ! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0
                             !
                             ! where e is the TKE dissipation, and N_0 and f_0 are the
                             ! reference buoyancy frequency and inertial frequencies respectively.
                             ! e_0 is the reference dissipation at (N_0,f_0). In the
                             ! previous version, N=N_0.
                             !
                             ! Additionally, the squared inverse relationship between
                             ! diapycnal diffusivities and stratification is included
                             !
                             ! kd = e/N^2
                             !
                             ! where kd is the diapycnal diffusivity.
                             ! This approach assumes that work done
                             ! against gravity is uniformly distributed
                             ! throughout the column. Whereas, kd=kd_0*e,
                             ! as in the original version, concentrates buoyancy
                             ! work in regions of strong stratification.

  logical :: kd_tanh_lat_fn  ! If true, use the tanh dependence of Kd_sfc on
                             ! latitude, like GFDL CM2.1/CM2M.  There is no physical
                             ! justification for this form, and it can not be
                             ! used with Henyey_IGW_background.
  real :: kd_tanh_lat_scale  ! A nondimensional scaling for the range of
                             ! diffusivities with Kd_tanh_lat_fn. Valid values
                             ! are in the range of -2 to 2; 0.4 reproduces CM2M.

  logical :: bottomdraglaw   ! If true, the  bottom stress is calculated with a
                             ! drag law c_drag*|u|*u.
  logical :: bbl_mixing_as_max !  If true, take the maximum of the diffusivity
                             ! from the BBL mixing and the other diffusivities.
                             ! Otherwise, diffusivities from the BBL_mixing is
                             ! added.
  logical :: use_lotw_bbl_diffusivity ! If true, use simpler/less precise, BBL diffusivity.
  logical :: lotw_bbl_use_omega ! If true, use simpler/less precise, BBL diffusivity.
  real    :: bbl_effic       ! efficiency with which the energy extracted
                             ! by bottom drag drives BBL diffusion (nondim)
  real    :: cdrag           ! quadratic drag coefficient (nondim)
  real    :: imax_decay      ! inverse of a maximum decay scale for
                             ! bottom-drag driven turbulence, (1/m)

  real    :: kd              ! interior diapycnal diffusivity (m2/s)
  real    :: kd_min          ! minimum diapycnal diffusivity (m2/s)
  real    :: kd_max          ! maximum increment for diapycnal diffusivity (m2/s)
                             ! Set to a negative value to have no limit.
  real    :: kd_add          ! uniform diffusivity added everywhere without
                             ! filtering or scaling (m2/s)
  real    :: kv              ! interior vertical viscosity (m2/s)
  real    :: kdml            ! mixed layer diapycnal diffusivity (m2/s)
                             ! when bulkmixedlayer==.false.
  real    :: hmix            ! mixed layer thickness (meter) when
                             ! bulkmixedlayer==.false.

  logical :: bryan_lewis_diffusivity ! If true, background vertical diffusivity
                                     ! uses Bryan-Lewis (1979) like tanh profile.
  real    :: kd_bryan_lewis_deep     ! abyssal value of Bryan-Lewis profile (m2/s)
  real    :: kd_bryan_lewis_surface  ! surface value of Bryan-Lewis profile (m2/s)
  real    :: bryan_lewis_depth_cent  ! center of transition depth in Bryan-Lewis (meter)
  real    :: bryan_lewis_width_trans ! width of transition for Bryan-Lewis (meter)

  real    :: n0_2omega       ! ratio of the typical Buoyancy frequency to
                             ! twice the Earth's rotation period, used with the
                             ! Henyey scaling from the mixing
  real    :: n2_floor_iomega2 ! floor applied to N2(k) scaled by Omega^2
                              ! If =0., N2(k) is positive definite
                              ! If =1., N2(k) > Omega^2 everywhere

  type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing

  real :: int_tide_decay_scale ! decay scale for internal wave TKE (meter)
  real :: mu_itides          ! efficiency for conversion of dissipation
                             ! to potential energy (nondimensional)
  real :: gamma_itides       ! fraction of local dissipation (nondimensional)
  real :: gamma_lee          ! fraction of local dissipation for lee waves
                             ! (Nikurashin's energy input) (nondimensional)
  real :: decay_scale_factor_lee ! Scaling factor for the decay scale of lee
                             ! wave energy dissipation (nondimensional)
  real :: min_zbot_itides    ! minimum depth for internal tide conversion (meter)
  logical :: int_tide_dissipation  ! Internal tide conversion (from barotropic) with
                                   ! the schemes of St Laurent et al (2002)/
                                   ! Simmons et al (2004)
  logical :: lowmode_itidal_dissipation ! Internal tide conversion (from low modes) with
                                        ! the schemes of St Laurent et al (2002)/
                                        ! Simmons et al (2004) !BDM
  integer :: int_tide_profile ! A coded integer indicating the vertical profile
                              ! for dissipation of the internal waves.  Schemes that
                              ! are currently encoded are St Laurent et al (2002) and
                              ! Polzin (2009).
  real :: nu_polzin      ! The non-dimensional constant used in Polzin form of
                         ! the vertical scale of decay of tidal dissipation
  real :: nbotref_polzin ! Reference value for the buoyancy frequency at the
                         ! ocean bottom used in Polzin formulation of the
                         ! vertical scale of decay of tidal dissipation (1/s)
  real :: polzin_decay_scale_factor ! Scaling factor for the decay length scale
                                    ! of the tidal dissipation profile in Polzin
                                    ! (nondimensional)
  real :: polzin_decay_scale_max_factor  ! The decay length scale of tidal
                         ! dissipation profile in Polzin formulation should not
                         ! exceed Polzin_decay_scale_max_factor * depth of the
                         ! ocean (nondimensional).
  real :: polzin_min_decay_scale ! minimum decay scale of the tidal dissipation
                                 ! profile in Polzin formulation (meter)
  logical :: lee_wave_dissipation ! Enable lee-wave driven mixing, following
                                  ! Nikurashin (2010), with a vertical energy
                                  ! deposition profile specified by Lee_wave_profile.
                                  ! St Laurent et al (2002) or
                                  ! Simmons et al (2004) scheme
  integer :: lee_wave_profile ! A coded integer indicating the vertical profile
                              ! for dissipation of the lee waves.  Schemes that are
                              ! currently encoded are St Laurent et al (2002) and
                              ! Polzin (2009).
  logical :: limit_dissipation ! If enabled, dissipation is limited to be larger
                               ! than the following:
  real :: dissip_min    ! Minimum dissipation (W/m3)
  real :: dissip_n0     ! Coefficient a in minimum dissipation = a+b*N (W/m3)
  real :: dissip_n1     ! Coefficient b in minimum dissipation = a+b*N (J/m3)
  real :: dissip_n2     ! Coefficient c in minimum dissipation = c*N2 (W m-3 s2)
  real :: dissip_kd_min ! Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2

  real :: tke_itide_max       ! maximum internal tide conversion (W m-2)
                              ! available to mix above the BBL
  real :: omega               ! Earth's rotation frequency (s-1)
  real :: utide               ! constant tidal amplitude (m s-1) used if
                              ! tidal amplitude file is not present
  real :: kappa_itides        ! topographic wavenumber and non-dimensional scaling
  real :: kappa_h2_factor     ! factor for the product of wavenumber * rms sgs height
  logical :: ml_radiation     ! allow a fraction of TKE available from wind work
                              ! to penetrate below mixed layer base with a vertical
                              ! decay scale determined by the minimum of
                              ! (1) The depth of the mixed layer, or
                              ! (2) An Ekman length scale.
                              ! Energy availble to drive mixing below the mixed layer is
                              ! given by E = ML_RAD_COEFF*MSTAR*USTAR**3.  Optionally, if
                              ! ML_rad_TKE_decay is true, this is further reduced by a factor
                              ! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is
                              ! calculated the same way as in the mixed layer code.
                              ! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2),
                              ! where N2 is the squared buoyancy frequency (s-2) and OMEGA2
                              ! is the rotation rate of the earth squared.
  real :: ml_rad_kd_max       ! Maximum diapycnal diffusivity due to turbulence
                              ! radiated from the base of the mixed layer (m2/s)
  real :: ml_rad_efold_coeff  ! non-dim coefficient to scale penetration depth
  real :: ml_rad_coeff        ! coefficient, which scales MSTAR*USTAR^3 to
                              ! obtain energy available for mixing below
                              ! mixed layer base (nondimensional)
  logical :: ml_rad_tke_decay ! If true, apply same exponential decay
                              ! to ML_rad as applied to the other surface
                              ! sources of TKE in the mixed layer code.
  real    :: ustar_min        ! A minimum value of ustar to avoid numerical
                              ! problems (m/s).  If the value is small enough,
                              ! this parameter should not affect the solution.
  real    :: tke_decay        ! ratio of natural Ekman depth to TKE decay scale (nondim)
  real    :: mstar            ! ratio of friction velocity cubed to
                              ! TKE input to the mixed layer (nondim)
  logical :: ml_use_omega     ! If true, use absolute rotation rate instead
                              ! of the vertical component of rotation when
                              ! setting the decay scale for mixed layer turbulence.
  real    :: ml_omega_frac    !   When setting the decay scale for turbulence, use
                              ! this fraction of the absolute rotation rate blended
                              ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2.
  logical :: user_change_diff ! If true, call user-defined code to change diffusivity.
  logical :: usekappashear    ! If true, use the kappa_shear module to find the
                              ! shear-driven diapycnal diffusivity.

  logical :: usecvmix         ! If true, use one of the CVMix modules to find
                              ! shear-driven diapycnal diffusivity.

  logical :: double_diffusion           ! If true, enable double-diffusive mixing.
  logical :: simple_tke_to_kd ! If true, uses a simple estimate of Kd/TKE that
                              ! does not rely on a layer-formulation.
  real    :: max_rrho_salt_fingers      ! max density ratio for salt fingering
  real    :: max_salt_diff_salt_fingers ! max salt diffusivity for salt fingers (m2/s)
  real    :: kv_molecular               ! molecular visc for double diff convect (m2/s)

  real, pointer, dimension(:,:) :: tke_niku    => null()
  real, pointer, dimension(:,:) :: tke_itidal  => null()
  real, pointer, dimension(:,:) :: nb          => null()
  real, pointer, dimension(:,:) :: mask_itidal => null()
  real, pointer, dimension(:,:) :: h2          => null()
  real, pointer, dimension(:,:) :: tideamp     => null() ! RMS tidal amplitude (m/s)

  character(len=200)                 :: inputdir
  type(user_change_diff_cs), pointer :: user_change_diff_csp => null()
  type(diag_to_z_cs),        pointer :: diag_to_z_csp        => null()
  type(kappa_shear_cs),      pointer :: kappashear_csp       => null()
  type(cvmix_shear_cs),      pointer :: cvmix_shear_csp      => null()
  type(int_tide_cs),         pointer :: int_tide_csp         => null()

  integer :: id_tke_itidal  = -1
  integer :: id_tke_leewave = -1
  integer :: id_maxtke      = -1
  integer :: id_tke_to_kd   = -1

  integer :: id_kd_itidal      = -1
  integer :: id_kd_niku        = -1
  integer :: id_kd_lowmode     = -1
  integer :: id_kd_user        = -1
  integer :: id_kd_layer       = -1
  integer :: id_kd_bbl         = -1
  integer :: id_kd_bbl_z       = -1
  integer :: id_kd_itidal_z    = -1
  integer :: id_kd_niku_z      = -1
  integer :: id_kd_lowmode_z   = -1
  integer :: id_kd_user_z      = -1
  integer :: id_kd_work        = -1
  integer :: id_kd_itidal_work = -1
  integer :: id_kd_niku_work   = -1
  integer :: id_kd_lowmode_work= -1

  integer :: id_fl_itidal                 = -1
  integer :: id_fl_lowmode                = -1
  integer :: id_polzin_decay_scale        = -1
  integer :: id_polzin_decay_scale_scaled = -1

  integer :: id_nb       = -1
  integer :: id_n2       = -1
  integer :: id_n2_z     = -1
  integer :: id_n2_bot   = -1
  integer :: id_n2_meanz = -1

  integer :: id_kt_extra   = -1
  integer :: id_ks_extra   = -1
  integer :: id_kt_extra_z = -1
  integer :: id_ks_extra_z = -1

end type set_diffusivity_cs

type diffusivity_diags
  real, pointer, dimension(:,:,:) :: &
    n2_3d          => null(),& ! squared buoyancy frequency at interfaces (1/s2)
    kd_itidal      => null(),& ! internal tide diffusivity at interfaces (m2/s)
    fl_itidal      => null(),& ! vertical flux of tidal turbulent dissipation (m3/s3)
    kd_lowmode     => null(),& ! internal tide diffusivity at interfaces
                               ! due to propagating low modes (m2/s) (BDM)
    fl_lowmode     => null(),& ! vertical flux of tidal turbulent dissipation
                               ! due to propagating low modes (m3/s3) (BDM)
    kd_niku        => null(),& ! lee-wave diffusivity at interfaces (m2/s)
    kd_user        => null(),& ! user-added diffusivity at interfaces (m2/s)
    kd_bbl         => null(),& ! BBL diffusivity at interfaces (m2/s)
    kd_work        => null(),& ! layer integrated work by diapycnal mixing (W/m2)
    kd_niku_work   => null(),& ! layer integrated work by lee-wave driven mixing (W/m2)
    kd_itidal_work => null(),& ! layer integrated work by int tide driven mixing (W/m2)
    kd_lowmode_work=> null(),& ! layer integrated work by low mode driven mixing (W/m2) BDM
    maxtke         => null(),& ! energy required to entrain to h_max (m3/s3)
    tke_to_kd      => null(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay))
                               ! between TKE dissipated within a layer and Kd
                               ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1
    kt_extra       => null(),& ! double diffusion diffusivity for temp (m2/s)
    ks_extra       => null()   ! double diffusion diffusivity for saln (m2/s)

  real, pointer, dimension(:,:) :: &
    tke_itidal_used           => null(),& ! internal tide TKE input at ocean bottom (W/m2)
    n2_bot                    => null(),& ! bottom squared buoyancy frequency (1/s2)
    n2_meanz                  => null(),& ! vertically averaged buoyancy frequency (1/s2)
    polzin_decay_scale_scaled => null(),& ! vertical scale of decay for tidal dissipation
    polzin_decay_scale        => null()   ! vertical decay scale for tidal diss with Polzin (meter)

end type diffusivity_diags

character*(20), parameter :: stlaurent_profile_string = "STLAURENT_02"
character*(20), parameter :: polzin_profile_string = "POLZIN_09"
integer,        parameter :: stlaurent_02 = 1
integer,        parameter :: polzin_09    = 2

! Clocks
integer :: id_clock_kappashear

contains

subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, &
                           G, GV, CS, Kd, Kd_int)
  type(ocean_grid_type),     intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type),   intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                             intent(in)    :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             intent(in)    :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                             intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                             intent(in)    :: u_h
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                             intent(in)    :: v_h
  type(thermo_var_ptrs),     intent(inout) :: tv   !< Structure with pointers to thermodynamic
                                                   !! fields. Out is for tv%TempxPmE.
  type(forcing),             intent(in)    :: fluxes !< Structure of surface fluxes that may be
                                                   !! used.
  type(optics_type),         pointer       :: optics
  type(vertvisc_type),       intent(inout) :: visc !< Structure containing vertical viscosities,
                                                   !! bottom boundary layer properies, and related
                                                   !! fields.
  real,                      intent(in)    :: dt   !< Time increment (sec).
  type(set_diffusivity_cs),  pointer       :: CS   !< Module control structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             intent(out)   :: Kd   !< Diapycnal diffusivity of each layer (m2/sec).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                   optional, intent(out)   :: Kd_int !< Diapycnal diffusivity at each interface
                                                   !! (m2/sec).

! Arguments:
!  (in)      u      - zonal velocity (m/s)
!  (in)      v      - meridional velocity (m/s)
!  (in)      h      - Layer thickness (m or kg/m2)
!  (in)      tv     - structure with pointers to thermodynamic fields
!  (in)      fluxes - structure of surface fluxes that may be used
!  (in)      visc   - structure containing vertical viscosities, bottom boundary
!                     layer properies, and related fields
!  (in)      dt     - time increment (sec)
!  (in)      G      - ocean grid structure
!  (in)      GV     - The ocean's vertical grid structure.
!  (in)      CS     - module control structure
!  (in)      j      - meridional index upon which to work
!  (out)     Kd     - diapycnal diffusivity of each layer (m2/sec)
!  (out,opt) Kd_int - diapycnal diffusivity at each interface (m2/sec)

  real, dimension(SZI_(G)) :: &
    depth, &      ! distance from surface of an interface (meter)
    N2_bot        ! bottom squared buoyancy frequency (1/s2)
  real, dimension(SZI_(G), SZJ_(G)) :: &
    Kd_sfc        ! surface value of the diffusivity (m2/s)

  type(diffusivity_diags) :: dd ! structure w/ arrays of pointers to avail diags

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    T_f, S_f      ! temperature and salinity (deg C and ppt);
                  ! massless layers filled vertically by diffusion.

  real, dimension(SZI_(G),SZK_(G)) :: &
    N2_lay, &     ! squared buoyancy frequency associated with layers (1/s2)
    maxTKE, &     ! energy required to entrain to h_max (m3/s3)
    TKE_to_Kd     ! conversion rate (~1.0 / (G_Earth + dRho_lay)) between
                  ! TKE dissipated within a layer and Kd in that layer, in
                  ! m2 s-1 / m3 s-3 = s2 m-1.

  real, dimension(SZI_(G),SZK_(G)+1) :: &
    N2_int,   &   ! squared buoyancy frequency associated at interfaces (1/s2)
    dRho_int, &   ! locally ref potential density difference across interfaces (in s-2) smg: or kg/m3?
    KT_extra, &   ! double difusion diffusivity on temperature (m2/sec)
    KS_extra      ! double difusion diffusivity on salinity (m2/sec)

  real :: I_trans       ! inverse of the transitional for Bryan-Lewis (1/m)
  real :: depth_c       ! depth of the center of a layer (meter)
  real :: I_Hmix        ! inverse of fixed mixed layer thickness (1/m)
  real :: I_Rho0        ! inverse of Boussinesq density (m3/kg)
  real :: I_x30         ! 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg)))
  real :: abs_sin       ! absolute value of sine of latitude (nondim)
  real :: atan_fn_sfc   ! surface value of Bryan-Lewis profile (nondim)
  real :: atan_fn_lay   ! value of Bryan-Lewis profile in layer middle (nondim)
  real :: I_atan_fn     ! inverse of change in Bryan-Lewis profile from surface to infinite depth (nondim)
  real :: deg_to_rad    ! factor converting degrees to radians, pi/180.
  real :: dissip        ! local variable for dissipation calculations (W/m3)
  real :: Omega2        ! squared absolute rotation rate (1/s2)
  real :: I_2Omega      ! 1/(2 Omega) (sec)
  real :: N_2Omega
  real :: N02_N2
  real :: epsilon

  logical   :: use_EOS      ! If true, compute density from T/S using equation of state.
  type(p3d) :: z_ptrs(6)    ! pointers to diagns to be interpolated into depth space
  integer   :: kb(szi_(g))  ! The index of the lightest layer denser than the
                            ! buffer layer.
  integer   :: num_z_diags  ! number of diagns to be interpolated to depth space
  integer   :: z_ids(6)     ! id numbers of diagns to be interpolated to depth space
  logical   :: showCallTree ! If true, show the call tree.

  integer :: i, j, k, is, ie, js, je, nz
  integer :: isd, ied, jsd, jed

  real      :: kappa_fill   ! diffusivity used to fill massless layers
  real      :: dt_fill      ! timestep used to fill massless layers

  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("set_diffusivity(), MOM_set_diffusivity.F90")

  if (.not.associated(cs)) call mom_error(fatal,"set_diffusivity: "//&
         "Module must be initialized before it is used.")

  i_rho0     = 1.0/gv%Rho0
  kappa_fill = 1.e-3 ! m2 s-1
  dt_fill    = 7200.
  deg_to_rad = atan(1.0)/45.0 ! = PI/180
  omega2     = cs%Omega*cs%Omega
  i_2omega   = 0.5/cs%Omega
  epsilon    = 1.e-10

  use_eos = associated(tv%eqn_of_state)

  if ((cs%double_diffusion) .and. &
      .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) &
    call mom_error(fatal, "set_diffusivity: visc%Kd_extra_T and "//&
         "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.")

  ! Set up arrays for diagnostics.

  if ((cs%id_N2 > 0) .or. (cs%id_N2_z > 0)) then
    allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1)) ; dd%N2_3d(:,:,:) = 0.0
  endif
  if ((cs%id_Kd_itidal > 0) .or. (cs%id_Kd_itidal_z > 0) .or. &
      (cs%id_Kd_Itidal_work > 0)) then
    allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Kd_itidal(:,:,:) = 0.0
  endif
  if ((cs%id_Kd_lowmode > 0) .or. (cs%id_Kd_lowmode_z > 0) .or. &
      (cs%id_Kd_lowmode_work > 0)) then
    allocate(dd%Kd_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Kd_lowmode(:,:,:) = 0.0
  endif
  if ( (cs%id_Fl_itidal > 0) ) then
    allocate(dd%Fl_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Fl_itidal(:,:,:) = 0.0
  endif
  if ( (cs%id_Fl_lowmode > 0) ) then
    allocate(dd%Fl_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Fl_lowmode(:,:,:) = 0.0
  endif
  if ( (cs%id_Polzin_decay_scale > 0) ) then
    allocate(dd%Polzin_decay_scale(isd:ied,jsd:jed))
    dd%Polzin_decay_scale(:,:) = 0.0
  endif
  if ( (cs%id_Polzin_decay_scale_scaled > 0) ) then
    allocate(dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed))
    dd%Polzin_decay_scale_scaled(:,:) = 0.0
  endif
  if ( (cs%id_N2_bot > 0) ) then
    allocate(dd%N2_bot(isd:ied,jsd:jed)) ; dd%N2_bot(:,:) = 0.0
  endif
  if ( (cs%id_N2_meanz > 0) ) then
    allocate(dd%N2_meanz(isd:ied,jsd:jed)) ; dd%N2_meanz(:,:) = 0.0
  endif
  if ((cs%id_Kd_Niku > 0) .or. (cs%id_Kd_Niku_z > 0) .or. &
      (cs%id_Kd_Niku_work > 0)) then
    allocate(dd%Kd_Niku(isd:ied,jsd:jed,nz+1)) ; dd%Kd_Niku(:,:,:) = 0.0
  endif
  if ((cs%id_Kd_user > 0) .or. (cs%id_Kd_user_z > 0)) then
    allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1)) ; dd%Kd_user(:,:,:) = 0.0
  endif
  if (cs%id_Kd_work > 0) then
    allocate(dd%Kd_work(isd:ied,jsd:jed,nz)) ; dd%Kd_work(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Niku_work > 0) then
    allocate(dd%Kd_Niku_work(isd:ied,jsd:jed,nz)) ; dd%Kd_Niku_work(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Itidal_work > 0) then
    allocate(dd%Kd_Itidal_work(isd:ied,jsd:jed,nz))
    dd%Kd_Itidal_work(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Lowmode_Work > 0) then
    allocate(dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz))
    dd%Kd_Lowmode_Work(:,:,:) = 0.0
  endif
  if (cs%id_TKE_itidal > 0) then
    allocate(dd%TKE_Itidal_used(isd:ied,jsd:jed)) ; dd%TKE_Itidal_used(:,:) = 0.
  endif
  if (cs%id_maxTKE > 0) then
    allocate(dd%maxTKE(isd:ied,jsd:jed,nz)) ; dd%maxTKE(:,:,:) = 0.0
  endif
  if (cs%id_TKE_to_Kd > 0) then
    allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0
  endif
  if ((cs%id_KT_extra > 0) .or. (cs%id_KT_extra_z > 0)) then
    allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0
  endif
  if ((cs%id_KS_extra > 0) .or. (cs%id_KS_extra_z > 0)) then
    allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0
  endif
  if ((cs%id_Kd_BBL > 0) .or. (cs%id_Kd_BBL_z > 0)) then
    allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0
  endif

  ! Smooth the properties through massless layers.
  if (use_eos) then
    if (cs%debug) then
      call hchksum(tv%T, "before vert_fill_TS tv%T",g%HI)
      call hchksum(tv%S, "before vert_fill_TS tv%S",g%HI)
      call hchksum(h, "before vert_fill_TS h",g%HI, scale=gv%H_to_m)
    endif
    call vert_fill_ts(h, tv%T, tv%S, kappa_fill, dt_fill, t_f, s_f, g, gv)
    if (cs%debug) then
      call hchksum(tv%T, "after vert_fill_TS tv%T",g%HI)
      call hchksum(tv%S, "after vert_fill_TS tv%S",g%HI)
      call hchksum(h, "after vert_fill_TS h",g%HI, scale=gv%H_to_m)
    endif
  endif

  if (cs%useKappaShear) then
    if (cs%debug) then
      call hchksum(u_h, "before calc_KS u_h",g%HI)
      call hchksum(v_h, "before calc_KS v_h",g%HI)
    endif
    call cpu_clock_begin(id_clock_kappashear)
    ! Changes: visc%Kd_turb, visc%TKE_turb (not clear that TKE_turb is used as input ????)
    ! Sets visc%Kv_turb
    call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_turb, visc%TKE_turb, &
                               visc%Kv_turb, dt, g, gv, cs%kappaShear_CSp)
    call cpu_clock_end(id_clock_kappashear)
    if (cs%debug) then
      call hchksum(visc%Kd_turb, "after calc_KS visc%Kd_turb",g%HI)
      call hchksum(visc%Kv_turb, "after calc_KS visc%Kv_turb",g%HI)
      call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb",g%HI)
    endif
    if (showcalltree) call calltree_waypoint("done with calculate_kappa_shear (set_diffusivity)")
 elseif (cs%useCVMix) then
    !NOTE{BGR}: this needs cleaned up.  Works in 1D case, not tested outside.
    call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_turb, visc%Kv_turb,g,gv,cs%CVMix_shear_CSp)
  elseif (associated(visc%Kv_turb)) then
    visc%Kv_turb(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled
  endif

!   Calculate the diffusivity, Kd, for each layer.  This would be
! the appropriate place to add a depth-dependent parameterization or
! another explicit parameterization of Kd.

  if (cs%Bryan_Lewis_diffusivity) then
!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc)
    do j=js,je ; do i=is,ie
      kd_sfc(i,j) = cs%Kd_Bryan_Lewis_surface
    enddo ; enddo
  else
!$OMP parallel do default(none) shared(is,ie,js,je,CS,Kd_sfc)
    do j=js,je ; do i=is,ie
      kd_sfc(i,j) = cs%Kd
    enddo ; enddo
  endif
  if (cs%Henyey_IGW_background) then
    i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G,deg_to_rad,epsilon,I_x30) &
!$OMP                          private(abs_sin)
    do j=js,je ; do i=is,ie
      abs_sin = abs(sin(g%geoLatT(i,j)*deg_to_rad))
      kd_sfc(i,j) = max(cs%Kd_min, kd_sfc(i,j) * &
           ((abs_sin * invcosh(cs%N0_2Omega/max(epsilon,abs_sin))) * i_x30) )
    enddo ; enddo
  elseif (cs%Kd_tanh_lat_fn) then
!$OMP parallel do default(none) shared(is,ie,js,je,Kd_sfc,CS,G)
    do j=js,je ; do i=is,ie
      !   The transition latitude and latitude range are hard-scaled here, since
      ! this is not really intended for wide-spread use, but rather for
      ! comparison with CM2M / CM2.1 settings.
      kd_sfc(i,j) = max(cs%Kd_min, kd_sfc(i,j) * (1.0 + &
          cs%Kd_tanh_lat_scale * 0.5*tanh((abs(g%geoLatT(i,j)) - 35.0)/5.0) ))
    enddo ; enddo
  endif

  if (cs%debug) call hchksum(kd_sfc,"Kd_sfc",g%HI,haloshift=0)
!$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, &
!$OMP                                  Kd,Kd_sfc,epsilon,deg_to_rad,I_2Omega,visc,    &
!$OMP                                  Kd_int,dt,u,v,Omega2)   &
!$OMP                          private(dRho_int,I_trans,atan_fn_sfc,I_atan_fn,atan_fn_lay, &
!$OMP                                  I_Hmix,depth_c,depth,N2_lay, N2_int, N2_bot,        &
!$OMP                                  I_x30,abs_sin,N_2Omega,N02_N2,KT_extra, KS_extra,   &
!$OMP                                  TKE_to_Kd,maxTKE,dissip,kb)
  do j=js,je
    ! Set up variables related to the stratification.
    call find_n2(h, tv, t_f, s_f, fluxes, j, g, gv, cs, drho_int, n2_lay, n2_int, n2_bot)
    if (associated(dd%N2_3d)) then
      do k=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,k) = n2_int(i,k) ; enddo ; enddo
    endif

    ! Set up the background diffusivity.
    if (cs%Bryan_Lewis_diffusivity) then
      i_trans = 1.0 / cs%Bryan_Lewis_width_trans
      atan_fn_sfc = atan(cs%Bryan_Lewis_depth_cent*i_trans)
      i_atan_fn = 1.0 / (2.0*atan(1.0) + atan_fn_sfc)
      do i=is,ie ; depth(i) = 0.0 ; enddo
      do k=1,nz ; do i=is,ie
        atan_fn_lay = atan((cs%Bryan_Lewis_depth_cent - &
                            (depth(i)+0.5*gv%H_to_m*h(i,j,k)))*i_trans)
        kd(i,j,k) = kd_sfc(i,j) + (cs%Kd_Bryan_Lewis_deep - kd_sfc(i,j)) * &
                                  (atan_fn_sfc - atan_fn_lay) * i_atan_fn
        depth(i) = depth(i) + gv%H_to_m*h(i,j,k)
      enddo ; enddo
    elseif ((.not.cs%bulkmixedlayer) .and. (cs%Kd /= cs%Kdml)) then
      i_hmix = 1.0 / cs%Hmix
      do i=is,ie ; depth(i) = 0.0 ; enddo
      do k=1,nz ; do i=is,ie
        depth_c = depth(i) + 0.5*gv%H_to_m*h(i,j,k)

        if (depth_c <= cs%Hmix) then ; kd(i,j,k) = cs%Kdml
        elseif (depth_c >= 2.0*cs%Hmix) then ; kd(i,j,k) = kd_sfc(i,j)
        else
          kd(i,j,k) = ((kd_sfc(i,j) - cs%Kdml) * i_hmix) * depth_c + &
                      (2.0*cs%Kdml - kd_sfc(i,j))
        endif

        depth(i) = depth(i) + gv%H_to_m*h(i,j,k)
      enddo ; enddo
    elseif (cs%Henyey_IGW_background_new) then
      i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
      do k=1,nz ; do i=is,ie
        abs_sin = max(epsilon,abs(sin(g%geoLatT(i,j)*deg_to_rad)))
        n_2omega = max(abs_sin,sqrt(n2_lay(i,k))*i_2omega)
        n02_n2 = (cs%N0_2Omega/n_2omega)**2
        kd(i,j,k) = max(cs%Kd_min, kd_sfc(i,j) * &
             ((abs_sin * invcosh(n_2omega/abs_sin)) * i_x30)*n02_n2)
      enddo ; enddo
    else
      do k=1,nz ; do i=is,ie
        kd(i,j,k) = kd_sfc(i,j)
      enddo ; enddo
    endif

    if (cs%double_diffusion) then
      call double_diffusion(tv, h, t_f, s_f, j, g, gv, cs, kt_extra, ks_extra)
      do k=2,nz ; do i=is,ie
        if (ks_extra(i,k) > kt_extra(i,k)) then ! salt fingering
          kd(i,j,k-1) = kd(i,j,k-1) + 0.5*kt_extra(i,k)
          kd(i,j,k)   = kd(i,j,k)   + 0.5*kt_extra(i,k)
          visc%Kd_extra_S(i,j,k) = ks_extra(i,k)-kt_extra(i,k)
          visc%Kd_extra_T(i,j,k) = 0.0
        elseif (kt_extra(i,k) > 0.0) then ! double-diffusive convection
          kd(i,j,k-1) = kd(i,j,k-1) + 0.5*ks_extra(i,k)
          kd(i,j,k)   = kd(i,j,k)   + 0.5*ks_extra(i,k)
          visc%Kd_extra_T(i,j,k) = kt_extra(i,k)-ks_extra(i,k)
          visc%Kd_extra_S(i,j,k) = 0.0
        else ! There is no double diffusion at this interface.
          visc%Kd_extra_T(i,j,k) = 0.0
          visc%Kd_extra_S(i,j,k) = 0.0
        endif
      enddo; enddo
      if (associated(dd%KT_extra)) then ; do k=1,nz+1 ; do i=is,ie
        dd%KT_extra(i,j,k) = kt_extra(i,k)
      enddo ; enddo ; endif

      if (associated(dd%KS_extra)) then ; do k=1,nz+1 ; do i=is,ie
        dd%KS_extra(i,j,k) = ks_extra(i,k)
      enddo ; enddo ; endif
    endif

  ! Add the input turbulent diffusivity.
    if (cs%useKappaShear .or. cs%useCVMix) then
      if (present(kd_int)) then
        do k=2,nz ; do i=is,ie
          kd_int(i,j,k) = visc%Kd_turb(i,j,k) + 0.5*(kd(i,j,k-1) + kd(i,j,k))
        enddo ; enddo
        do i=is,ie
          kd_int(i,j,1) = visc%Kd_turb(i,j,1) ! This isn't actually used. It could be 0.
          kd_int(i,j,nz+1) = 0.0
        enddo
      endif
      do k=1,nz ; do i=is,ie
        kd(i,j,k) = kd(i,j,k) + 0.5*(visc%Kd_turb(i,j,k) + visc%Kd_turb(i,j,k+1))
      enddo ; enddo
    else
      if (present(kd_int)) then
        do i=is,ie
          kd_int(i,j,1) = kd(i,j,1) ; kd_int(i,j,nz+1) = 0.0
        enddo
        do k=2,nz ; do i=is,ie
          kd_int(i,j,k) = 0.5*(kd(i,j,k-1) + kd(i,j,k))
        enddo ; enddo
      endif
    endif

    call find_tke_to_kd(h, tv, drho_int, n2_lay, j, dt, g, gv, cs, tke_to_kd, &
                        maxtke, kb)
    if (associated(dd%maxTKE)) then ; do k=1,nz ; do i=is,ie
      dd%maxTKE(i,j,k) = maxtke(i,k)
    enddo ; enddo ; endif
    if (associated(dd%TKE_to_Kd)) then ; do k=1,nz ; do i=is,ie
      dd%TKE_to_Kd(i,j,k) = tke_to_kd(i,k)
    enddo ; enddo ; endif

    ! Add the ML_Rad diffusivity.
    if (cs%ML_radiation) &
      call add_mlrad_diffusivity(h, fluxes, j, g, gv, cs, kd, tke_to_kd, kd_int)

    ! Add the Nikurashin and / or tidal bottom-driven mixing
    if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) &
      call add_int_tide_diffusivity(h, n2_bot, j, tke_to_kd, maxtke, g, gv, cs, &
                                    dd, n2_lay, kd, kd_int)

    ! This adds the diffusion sustained by the energy extracted from the flow
    ! by the bottom drag.
    if (cs%bottomdraglaw .and. (cs%BBL_effic>0.0)) then
      if (cs%use_LOTW_BBL_diffusivity) then
        call add_lotw_bbl_diffusivity(h, u, v, tv, fluxes, visc, j, n2_int, g, gv, cs,  &
                                      kd, kd_int, dd%Kd_BBL)
      else
        call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, tke_to_kd, &
                                  maxtke, kb, g, gv, cs, kd, kd_int, dd%Kd_BBL)
      endif
    endif

    if (cs%limit_dissipation) then
      do k=2,nz-1 ; do i=is,ie
      ! This calculates the dissipation ONLY from Kd calculated in this routine
      ! dissip has units of W/m3 (kg/m3 * m2/s * 1/s2 = J/s/m3)
      !   1) a global constant,
      !   2) a dissipation proportional to N (aka Gargett) and
      !   3) dissipation corresponding to a (nearly) constant diffusivity.
        dissip = max( cs%dissip_min, &   ! Const. floor on dissip.
                      cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_lay(i,k)), & ! Floor aka Gargett
                      cs%dissip_N2 * n2_lay(i,k) ) ! Floor of Kd_min*rho0/F_Ri
        kd(i,j,k) = max( kd(i,j,k) , &  ! Apply floor to Kd
           dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_lay(i,k) + omega2))) )
      enddo ; enddo

      if (present(kd_int)) then ; do k=2,nz ; do i=is,ie
      ! This calculates the dissipation ONLY from Kd calculated in this routine
      ! dissip has units of W/m3 (kg/m3 * m2/s * 1/s2 = J/s/m3)
      !   1) a global constant,
      !   2) a dissipation proportional to N (aka Gargett) and
      !   3) dissipation corresponding to a (nearly) constant diffusivity.
        dissip = max( cs%dissip_min, &   ! Const. floor on dissip.
                      cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_int(i,k)), & ! Floor aka Gargett
                      cs%dissip_N2 * n2_int(i,k) ) ! Floor of Kd_min*rho0/F_Ri
        kd_int(i,j,k) = max( kd_int(i,j,k) , &  ! Apply floor to Kd
           dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_int(i,k) + omega2))) )
      enddo ; enddo ; endif
    endif

    if (associated(dd%Kd_work)) then
      do k=1,nz ; do i=is,ie
        dd%Kd_Work(i,j,k)  = gv%Rho0 * kd(i,j,k) * n2_lay(i,k) * &
                             gv%H_to_m*h(i,j,k)  ! Watt m-2 s or kg s-3
      enddo ; enddo
    endif
  enddo ! j-loop

  if (cs%debug) then
    call hchksum(kd,"BBL Kd",g%HI,haloshift=0)
    if (cs%useKappaShear) call hchksum(visc%Kd_turb,"Turbulent Kd",g%HI,haloshift=0)
    if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then
      call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, &
                    g%HI, 0, symmetric=.true.)
    endif
    if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then
      call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, &
                    visc%bbl_thick_v, g%HI, 0, symmetric=.true.)
    endif
    if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then
      call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, g%HI, 0, symmetric=.true.)
    endif
  endif

  if (cs%Kd_add > 0.0) then
    if (present(kd_int)) then
!$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_int,CS,Kd)
      do k=1,nz ; do j=js,je ; do i=is,ie
        kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add
        kd(i,j,k) = kd(i,j,k) + cs%Kd_add
      enddo ; enddo ; enddo
    else
!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,Kd)
      do k=1,nz ; do j=js,je ; do i=is,ie
        kd(i,j,k) = kd(i,j,k) + cs%Kd_add
      enddo ; enddo ; enddo
    endif
  endif

  if (cs%user_change_diff) then
    call user_change_diff(h, tv, g, cs%user_change_diff_CSp, kd, kd_int, &
                          t_f, s_f, dd%Kd_user)
  endif

  if (cs%id_Kd_layer > 0) call post_data(cs%id_Kd_layer, kd, cs%diag)

  num_z_diags = 0
  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) then
    if (cs%id_TKE_itidal  > 0) call post_data(cs%id_TKE_itidal,  dd%TKE_itidal_used, cs%diag)
    if (cs%id_TKE_leewave > 0) call post_data(cs%id_TKE_leewave, cs%TKE_Niku,        cs%diag)
    if (cs%id_Nb          > 0) call post_data(cs%id_Nb,      cs%Nb,      cs%diag)
    if (cs%id_N2          > 0) call post_data(cs%id_N2,      dd%N2_3d,   cs%diag)
    if (cs%id_N2_bot      > 0) call post_data(cs%id_N2_bot,  dd%N2_bot,  cs%diag)
    if (cs%id_N2_meanz    > 0) call post_data(cs%id_N2_meanz,dd%N2_meanz,cs%diag)

    if (cs%id_Fl_itidal > 0) call post_data(cs%id_Fl_itidal, dd%Fl_itidal, cs%diag)
    if (cs%id_Kd_itidal > 0) call post_data(cs%id_Kd_itidal, dd%Kd_itidal, cs%diag)
    if (cs%id_Kd_Niku   > 0) call post_data(cs%id_Kd_Niku,   dd%Kd_Niku,   cs%diag)
    if (cs%id_Kd_lowmode> 0) call post_data(cs%id_Kd_lowmode, dd%Kd_lowmode, cs%diag)
    if (cs%id_Fl_lowmode> 0) call post_data(cs%id_Fl_lowmode, dd%Fl_lowmode, cs%diag)
    if (cs%id_Kd_user   > 0) call post_data(cs%id_Kd_user,   dd%Kd_user,   cs%diag)
    if (cs%id_Kd_Work   > 0) call post_data(cs%id_Kd_Work,   dd%Kd_Work,   cs%diag)
    if (cs%id_Kd_Itidal_Work > 0) &
      call post_data(cs%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, cs%diag)
    if (cs%id_Kd_Niku_Work > 0) call post_data(cs%id_Kd_Niku_Work, dd%Kd_Niku_Work, cs%diag)
    if (cs%id_Kd_Lowmode_Work > 0) &
      call post_data(cs%id_Kd_Lowmode_Work, dd%Kd_Lowmode_Work, cs%diag)
    if (cs%id_maxTKE > 0) call post_data(cs%id_maxTKE, dd%maxTKE, cs%diag)
    if (cs%id_TKE_to_Kd > 0) call post_data(cs%id_TKE_to_Kd, dd%TKE_to_Kd, cs%diag)

    if (cs%id_Polzin_decay_scale > 0 ) &
      call post_data(cs%id_Polzin_decay_scale, dd%Polzin_decay_scale, cs%diag)
    if (cs%id_Polzin_decay_scale_scaled > 0 ) &
      call post_data(cs%id_Polzin_decay_scale_scaled, dd%Polzin_decay_scale_scaled, cs%diag)

    if (cs%id_Kd_itidal_z > 0) then
      num_z_diags        = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_Kd_itidal_z
      z_ptrs(num_z_diags)%p => dd%Kd_itidal
    endif

    if (cs%id_Kd_Niku_z > 0) then
      num_z_diags        = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_Kd_Niku_z
      z_ptrs(num_z_diags)%p => dd%Kd_Niku
    endif

    if (cs%id_Kd_lowmode_z > 0) then
      num_z_diags        = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_Kd_lowmode_z
      z_ptrs(num_z_diags)%p => dd%Kd_lowmode
    endif

    if (cs%id_N2_z > 0) then
      num_z_diags = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_N2_z
      z_ptrs(num_z_diags)%p => dd%N2_3d
    endif

    if (cs%id_Kd_user_z > 0) then
      num_z_diags        = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_Kd_user_z
      z_ptrs(num_z_diags)%p => dd%Kd_user
    endif

  endif

  if (cs%id_KT_extra > 0) call post_data(cs%id_KT_extra, dd%KT_extra, cs%diag)
  if (cs%id_KS_extra > 0) call post_data(cs%id_KS_extra, dd%KS_extra, cs%diag)
  if (cs%id_Kd_BBL > 0)   call post_data(cs%id_Kd_BBL, dd%Kd_BBL, cs%diag)

  if (cs%id_KT_extra_z > 0) then
      num_z_diags = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_KT_extra_z
      z_ptrs(num_z_diags)%p => dd%KT_extra
  endif

  if (cs%id_KS_extra_z > 0) then
      num_z_diags = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_KS_extra_z
      z_ptrs(num_z_diags)%p => dd%KS_extra
  endif

  if (cs%id_Kd_BBL_z > 0) then
      num_z_diags = num_z_diags + 1
      z_ids(num_z_diags) = cs%id_Kd_BBL_z
      z_ptrs(num_z_diags)%p => dd%KS_extra
  endif

  if (num_z_diags > 0) &
    call calc_zint_diags(h, z_ptrs, z_ids, num_z_diags, g, gv, cs%diag_to_Z_CSp)

  if (associated(dd%N2_3d)) deallocate(dd%N2_3d)
  if (associated(dd%Kd_itidal)) deallocate(dd%Kd_itidal)
  if (associated(dd%Kd_lowmode)) deallocate(dd%Kd_lowmode)
  if (associated(dd%Fl_itidal)) deallocate(dd%Fl_itidal)
  if (associated(dd%Fl_lowmode)) deallocate(dd%Fl_lowmode)
  if (associated(dd%Polzin_decay_scale)) deallocate(dd%Polzin_decay_scale)
  if (associated(dd%Polzin_decay_scale_scaled)) deallocate(dd%Polzin_decay_scale_scaled)
  if (associated(dd%N2_bot)) deallocate(dd%N2_bot)
  if (associated(dd%N2_meanz)) deallocate(dd%N2_meanz)
  if (associated(dd%Kd_work)) deallocate(dd%Kd_work)
  if (associated(dd%Kd_user)) deallocate(dd%Kd_user)
  if (associated(dd%Kd_Niku)) deallocate(dd%Kd_Niku)
  if (associated(dd%Kd_Niku_work)) deallocate(dd%Kd_Niku_work)
  if (associated(dd%Kd_Itidal_Work))  deallocate(dd%Kd_Itidal_Work)
  if (associated(dd%Kd_Lowmode_Work)) deallocate(dd%Kd_Lowmode_Work)
  if (associated(dd%TKE_itidal_used)) deallocate(dd%TKE_itidal_used)
  if (associated(dd%maxTKE)) deallocate(dd%maxTKE)
  if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd)
  if (associated(dd%KT_extra)) deallocate(dd%KT_extra)
  if (associated(dd%KS_extra)) deallocate(dd%KS_extra)
  if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL)

  if (showcalltree) call calltree_leave("set_diffusivity()")

end subroutine set_diffusivity

subroutine find_tke_to_kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, CS, &
                          TKE_to_Kd, maxTKE, kb)
  type(ocean_grid_type),                   intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                 intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)   :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                   intent(in)    :: tv
  real, dimension(SZI_(G),SZK_(G)+1),      intent(in)    :: dRho_int
  real, dimension(SZI_(G),SZK_(G)),        intent(in)    :: N2_lay
  integer,                                 intent(in)    :: j
  real,                                    intent(in)    :: dt
  type(set_diffusivity_cs),                pointer       :: CS
  real, dimension(SZI_(G),SZK_(G)),        intent(out)   :: TKE_to_Kd, maxTKE
  integer, dimension(SZI_(G)),             intent(out)   :: kb

  real, dimension(SZI_(G),SZK_(G)) :: &
    ds_dsp1, &    ! coordinate variable (sigma-2) difference across an
                  ! interface divided by the difference across the interface
                  ! below it (nondimensional)
    dsp1_ds, &    ! inverse coordinate variable (sigma-2) difference
                  ! across an interface times the difference across the
                  ! interface above it (nondimensional)
    rho_0,   &    ! Layer potential densities relative to surface pressure (kg/m3)
    maxent        ! maxEnt is the maximum value of entrainment from below (with
                  ! compensating entrainment from above to keep the layer
                  ! density from changing) that will not deplete all of the
                  ! layers above or below a layer within a timestep (meter)
  real, dimension(SZI_(G)) :: &
    htot,    &    ! total thickness above or below a layer, or the
                  ! integrated thickness in the BBL (meter)
    mfkb,    &    ! total thickness in the mixed and buffer layers
                  ! times ds_dsp1 (meter)
    p_ref,   &    ! array of tv%P_Ref pressures
    rcv_kmb, &    ! coordinate density in the lowest buffer layer
    p_0           ! An array of 0 pressures

  real :: dh_max      ! maximum amount of entrainment a layer could
                      ! undergo before entraining all fluid in the layers
                      ! above or below (meter)
  real :: dRho_lay    ! density change across a layer (kg/m3)
  real :: Omega2      ! rotation rate squared (1/s2)
  real :: G_Rho0      ! gravitation accel divided by Bouss ref density (m4 s-2 kg-1)
  real :: I_Rho0      ! inverse of Boussinesq reference density (m3/kg)
  real :: I_dt        ! 1/dt (1/sec)
  real :: H_neglect   ! negligibly small thickness (units as h)
  real :: hN2pO2      ! h * (N^2 + Omega^2), in m s-2.
  logical :: do_i(szi_(g))

  integer :: i, k, is, ie, nz, i_rem, kmb, kb_min
  is = g%isc ; ie = g%iec ; nz = g%ke

  i_dt      = 1.0/dt
  omega2    = cs%Omega**2
  g_rho0    = gv%g_Earth / gv%Rho0
  h_neglect = gv%H_subroundoff
  i_rho0    = 1.0/gv%Rho0

  ! Simple but coordinate-independent estimate of Kd/TKE
  if (cs%simple_TKE_to_Kd) then
    do k=1,nz ; do i=is,ie
      hn2po2 = ( gv%H_to_m * h(i,j,k) ) * ( n2_lay(i,k) + omega2 ) ! Units of m s-2.
      if (hn2po2>0.) then
        tke_to_kd(i,k) = 1./ hn2po2 ! Units of s2 m-1.
      else; tke_to_kd(i,k) = 0.; endif
      ! The maximum TKE conversion we allow is really a statement
      ! about the upper diffusivity we allow. Kd_max must be set.
      maxtke(i,k) = hn2po2 * cs%Kd_max ! Units of m3 s-3.
    enddo ; enddo
    kb(is:ie) = -1 ! kb should not be used by any code in non-layered mode -AJA
    return
  endif

  ! Determine kb - the index of the shallowest active interior layer.
  if (cs%bulkmixedlayer) then
    kmb = gv%nk_rho_varies
    do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ; enddo
    do k=1,nz
      call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p_0,rho_0(:,k),&
                             is,ie-is+1,tv%eqn_of_state)
    enddo
    call calculate_density(tv%T(:,j,kmb),tv%S(:,j,kmb),p_ref,rcv_kmb,&
                           is,ie-is+1,tv%eqn_of_state)

    kb_min = kmb+1
    do i=is,ie
      !   Determine the next denser layer than the buffer layer in the
      ! coordinate density (sigma-2).
      do k=kmb+1,nz-1 ; if (rcv_kmb(i) <= gv%Rlay(k)) exit ; enddo
      kb(i) = k

    !   Backtrack, in case there are massive layers above that are stable
    ! in sigma-0.
      do k=kb(i)-1,kmb+1,-1
        if (rho_0(i,kmb) > rho_0(i,k)) exit
        if (h(i,j,k)>2.0*gv%Angstrom) kb(i) = k
      enddo
    enddo

    call set_density_ratios(h, tv, kb, g, gv, cs, j, ds_dsp1, rho_0)
  else ! not bulkmixedlayer
    kb_min = 2 ; kmb = 0
    do i=is,ie ; kb(i) = 1 ; enddo
    call set_density_ratios(h, tv, kb, g, gv, cs, j, ds_dsp1)
  endif

  ! Determine maxEnt - the maximum permitted entrainment from below by each
  ! interior layer.
  do k=2,nz-1 ; do i=is,ie
    dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
  enddo ; enddo
  do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo

  if (cs%bulkmixedlayer) then
    kmb = gv%nk_rho_varies
    do i=is,ie
      htot(i) = gv%H_to_m*h(i,j,kmb)
      mfkb(i) = 0.0
      if (kb(i) < nz) &
        mfkb(i) = ds_dsp1(i,kb(i)) * (gv%H_to_m*(h(i,j,kmb) - gv%Angstrom))
    enddo
    do k=1,kmb-1 ; do i=is,ie
      htot(i) = htot(i) + gv%H_to_m*h(i,j,k)
      mfkb(i) = mfkb(i) + ds_dsp1(i,k+1)*(gv%H_to_m*(h(i,j,k) - gv%Angstrom))
    enddo ; enddo
  else
    do i=is,i
      maxent(i,1) = 0.0 ; htot(i) = gv%H_to_m*(h(i,j,1) - gv%Angstrom)
    enddo
  endif
  do k=kb_min,nz-1 ; do i=is,ie
    if (k == kb(i)) then
      maxent(i,kb(i))= mfkb(i)
    elseif (k > kb(i)) then
      maxent(i,k) = (1.0/dsp1_ds(i,k))*(maxent(i,k-1) + htot(i))
!        maxEnt(i,k) = ds_dsp1(i,k)*(maxEnt(i,k-1) + htot(i)) ! BITWISE CHG
      htot(i) = htot(i) + gv%H_to_m*(h(i,j,k) - gv%Angstrom)
    endif
  enddo ; enddo

  do i=is,ie
    htot(i) = gv%H_to_m*(h(i,j,nz) - gv%Angstrom) ; maxent(i,nz) = 0.0
    do_i(i) = (g%mask2dT(i,j) > 0.5)
  enddo
  do k=nz-1,kb_min,-1
    i_rem = 0
    do i=is,ie ; if (do_i(i)) then
      if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
      i_rem = i_rem + 1  ! Count the i-rows that are still being worked on.
      maxent(i,k) = min(maxent(i,k),dsp1_ds(i,k+1)*maxent(i,k+1) + htot(i))
      htot(i) = htot(i) + gv%H_to_m*(h(i,j,k) - gv%Angstrom)
    endif ; enddo
    if (i_rem == 0) exit
  enddo ! k-loop

  ! Now set maxTKE and TKE_to_Kd.
  do i=is,ie
    maxtke(i,1) = 0.0 ; tke_to_kd(i,1) = 0.0
    maxtke(i,nz) = 0.0 ; tke_to_kd(i,nz) = 0.0
  enddo
  do k=2,kmb ; do i=is,ie
    maxtke(i,k) = 0.0
    tke_to_kd(i,k) = 1.0 / ((n2_lay(i,k) + omega2) * &
                            (gv%H_to_m*(h(i,j,k) + h_neglect)))
  enddo ; enddo
  do k=kmb+1,kb_min-1 ; do i=is,ie
    !   These are the properties in the deeper mixed and buffer layers, and
    ! should perhaps be revisited.
    maxtke(i,k) = 0.0 ; tke_to_kd(i,k) = 0.0
  enddo ; enddo
  do k=kb_min,nz-1 ; do i=is,ie
    if (k<kb(i)) then
      maxtke(i,k) = 0.0
      tke_to_kd(i,k) = 0.0
    else
      ! maxTKE is found by determining the kappa that gives maxEnt.
      ! ### This should be 1 / G_Earth * (delta rho_InSitu)
      !  kappa_max = I_dt * dRho_int(i,K+1) * maxEnt(i,k) * &
      !             (GV%H_to_m*h(i,j,k) + dh_max) / dRho_lay
      !  maxTKE(i,k) = GV%g_Earth * dRho_lay * kappa_max
      ! dRho_int should already be non-negative, so the max is redundant?
      dh_max = maxent(i,k) * (1.0 + dsp1_ds(i,k))
      drho_lay = 0.5 * max(drho_int(i,k) + drho_int(i,k+1), 0.0)
      maxtke(i,k) = i_dt * ((gv%g_Earth * i_rho0) * &
          (0.5*max(drho_int(i,k+1) + dsp1_ds(i,k)*drho_int(i,k),0.0))) * &
                   ((gv%H_to_m*h(i,j,k) + dh_max) * maxent(i,k))
      tke_to_kd(i,k) = 1.0 / (g_rho0 * drho_lay + &
                              cs%Omega**2 * gv%H_to_m*(h(i,j,k) + h_neglect))
    endif
  enddo ; enddo

end subroutine find_tke_to_kd

subroutine find_n2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, &
                   N2_lay, N2_int, N2_bot)
  type(ocean_grid_type),                    intent(in)   :: G    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)   :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)   :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                    intent(in)   :: tv
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)   :: T_f, S_f
  type(forcing),                            intent(in)   :: fluxes
  integer,                                  intent(in)   :: j
  type(set_diffusivity_cs),                 pointer      :: CS
  real, dimension(SZI_(G),SZK_(G)+1),       intent(out)  :: dRho_int, N2_int
  real, dimension(SZI_(G),SZK_(G)),         intent(out)  :: N2_lay
  real, dimension(SZI_(G)),                 intent(out)  :: N2_bot

  real, dimension(SZI_(G),SZK_(G)+1) :: &
    dRho_int_unfilt, & ! unfiltered density differences across interfaces
    dRho_dT,         & ! partial derivative of density wrt temp (kg m-3 degC-1)
    dRho_dS            ! partial derivative of density wrt saln (kg m-3 PPT-1)

  real, dimension(SZI_(G)) :: &
    pres,      &  ! pressure at each interface (Pa)
    Temp_int,  &  ! temperature at each interface (degC)
    Salin_int, &  ! salinity at each interface (PPT)
    drho_bot,  &
    h_amp,     &
    hb,        &
    z_from_bot

  real :: Rml_base  ! density of the deepest variable density layer
  real :: dz_int    ! thickness associated with an interface (meter)
  real :: G_Rho0    ! gravitation acceleration divided by Bouss reference density (m4 s-2 kg-1)
  real :: H_neglect ! negligibly small thickness, in the same units as h.

  logical :: do_i(szi_(g)), do_any
  integer :: i, k, is, ie, nz

  is = g%isc ; ie = g%iec ; nz = g%ke
  g_rho0    = gv%g_Earth / gv%Rho0
  h_neglect = gv%H_subroundoff

  ! Find the (limited) density jump across each interface.
  do i=is,ie
    drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
    drho_int_unfilt(i,1) = 0.0 ; drho_int_unfilt(i,nz+1) = 0.0
  enddo
  if (associated(tv%eqn_of_state)) then
    if (associated(fluxes%p_surf)) then
      do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
    else
      do i=is,ie ; pres(i) = 0.0 ; enddo
    endif
    do k=2,nz
      do i=is,ie
        pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
        temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
        salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
      enddo
      call calculate_density_derivs(temp_int, salin_int, pres, &
               drho_dt(:,k), drho_ds(:,k), is, ie-is+1, tv%eqn_of_state)
      do i=is,ie
        drho_int(i,k) = max(drho_dt(i,k)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
                            drho_ds(i,k)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
        drho_int_unfilt(i,k) = max(drho_dt(i,k)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + &
                            drho_ds(i,k)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0)
      enddo
    enddo
  else
    do k=2,nz ; do i=is,ie
      drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
    enddo ; enddo
  endif

  ! Set the buoyancy frequencies.
  do k=1,nz ; do i=is,ie
    n2_lay(i,k) = g_rho0 * 0.5*(drho_int(i,k) + drho_int(i,k+1)) / &
                  (gv%H_to_m*(h(i,j,k) + h_neglect))
  enddo ; enddo
  do i=is,ie ; n2_int(i,1) = 0.0 ; n2_int(i,nz+1) = 0.0 ; enddo
  do k=2,nz ; do i=is,ie
    n2_int(i,k) = g_rho0 * drho_int(i,k) / &
                  (0.5*gv%H_to_m*(h(i,j,k-1) + h(i,j,k) + h_neglect))
  enddo ; enddo

  ! Find the bottom boundary layer stratification, and use this in the deepest layers.
  do i=is,ie
    hb(i) = 0.0 ; drho_bot(i) = 0.0
    z_from_bot(i) = 0.5*gv%H_to_m*h(i,j,nz)
    do_i(i) = (g%mask2dT(i,j) > 0.5)

    if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) then
      h_amp(i) = sqrt(cs%h2(i,j)) ! for computing Nb
    else
      h_amp(i) = 0.0
    endif
  enddo

  do k=nz,2,-1
    do_any = .false.
    do i=is,ie ; if (do_i(i)) then
      dz_int = 0.5*gv%H_to_m*(h(i,j,k) + h(i,j,k-1))
      z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above

      hb(i) = hb(i) + dz_int
      drho_bot(i) = drho_bot(i) + drho_int(i,k)

      if (z_from_bot(i) > h_amp(i)) then
        if (k>2) then
          ! Always include at least one full layer.
          hb(i) = hb(i) + 0.5*gv%H_to_m*(h(i,j,k-1) + h(i,j,k-2))
          drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
        endif
        do_i(i) = .false.
      else
        do_any = .true.
      endif
    endif ; enddo
    if (.not.do_any) exit
  enddo

  do i=is,ie
    if (hb(i) > 0.0) then
      n2_bot(i) = (g_rho0 * drho_bot(i)) / hb(i)
    else ;  n2_bot(i) = 0.0 ; endif
    z_from_bot(i) = 0.5*gv%H_to_m*h(i,j,nz)
    do_i(i) = (g%mask2dT(i,j) > 0.5)
  enddo

  do k=nz,2,-1
    do_any = .false.
    do i=is,ie ; if (do_i(i)) then
      dz_int = 0.5*gv%H_to_m*(h(i,j,k) + h(i,j,k-1))
      z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above

      n2_int(i,k) = n2_bot(i)
      if (k>2) n2_lay(i,k-1) = n2_bot(i)

      if (z_from_bot(i) > h_amp(i)) then
        if (k>2) n2_int(i,k-1) = n2_bot(i)
        do_i(i) = .false.
      else
        do_any = .true.
      endif
    endif ; enddo
    if (.not.do_any) exit
  enddo

  if (associated(tv%eqn_of_state)) then
    do k=1,nz+1 ; do i=is,ie
      drho_int(i,k) = drho_int_unfilt(i,k)
    enddo ; enddo
  endif

end subroutine find_n2

!> This subroutine sets the additional diffusivities of temperature and
!! salinity due to double diffusion, using the same functional form as is
!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates
!! what was in Large et al. (1994).  All the coefficients here should probably
!! be made run-time variables rather than hard-coded constants.
!!
!! \todo Find reference for NCAR tech note above.
subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd)
  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(thermo_var_ptrs),    intent(in)  :: tv  !< Structure containing pointers to any available
                                               !! thermodynamic fields; absent fields have NULL
                                               !! ptrs.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)  :: h   !< Layer thicknesses, in H (usually m or kg m-2).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)  :: T_f !< layer temp in C with the values in massless layers
                                               !! filled vertically by diffusion.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)  :: S_f !< Layer salinities in PPT with values in massless
                                               !! layers filled vertically by diffusion.
  integer,                  intent(in)  :: j   !< Meridional index upon which to work.
  type(set_diffusivity_cs), pointer     :: CS  !< Module control structure.
  real, dimension(SZI_(G),SZK_(G)+1),       &
                            intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal
                                               !! diffusivity for temp (m2/sec).
  real, dimension(SZI_(G),SZK_(G)+1),       &
                            intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal
                                               !! diffusivity for saln (m2/sec).

! Arguments:
!  (in)      tv      - structure containing pointers to any available
!                      thermodynamic fields; absent fields have NULL ptrs
!  (in)      h       - layer thickness (m or kg m-2)
!  (in)      T_f     - layer temp in C with the values in massless layers
!                      filled vertically by diffusion
!  (in)      S_f     - layer salinities in PPT with values in massless layers
!                      filled vertically by diffusion
!  (in)      G       - ocean grid structure
!  (in)      GV      - The ocean's vertical grid structure.
!  (in)      CS      - module control structure
!  (in)      j       - meridional index upon which to work
!  (out)     Kd_T_dd - interface double diffusion diapycnal diffusivity for temp (m2/sec)
!  (out)     Kd_S_dd - interface double diffusion diapycnal diffusivity for saln (m2/sec)

! This subroutine sets the additional diffusivities of temperature and
! salinity due to double diffusion, using the same functional form as is
! used in MOM4.1, and taken from an NCAR technical note (###REF?) that updates
! what was in Large et al. (1994).  All the coefficients here should probably
! be made run-time variables rather than hard-coded constants.

  real, dimension(SZI_(G)) :: &
    dRho_dT,  &    ! partial derivatives of density wrt temp (kg m-3 degC-1)
    dRho_dS,  &    ! partial derivatives of density wrt saln (kg m-3 PPT-1)
    pres,     &    ! pressure at each interface (Pa)
    Temp_int, &    ! temp and saln at interfaces
    Salin_int

  real ::  alpha_dT ! density difference between layers due to temp diffs (kg/m3)
  real ::  beta_dS  ! density difference between layers due to saln diffs (kg/m3)

  real  :: Rrho    ! vertical density ratio
  real  :: diff_dd ! factor for double-diffusion
  real  :: prandtl ! flux ratio for diffusive convection regime

  real, parameter :: Rrho0  = 1.9          ! limit for double-diffusive density ratio
  real, parameter :: dsfmax = 1.e-4        ! max diffusivity in case of salt fingering
  real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity  (m2/sec)

  integer :: i, k, is, ie, nz
  is = g%isc ; ie = g%iec ; nz = g%ke

  if (associated(tv%eqn_of_state)) then
    do i=is,ie
      pres(i) = 0.0 ; kd_t_dd(i,1) = 0.0 ; kd_s_dd(i,1) = 0.0
      kd_t_dd(i,nz+1) = 0.0 ; kd_s_dd(i,nz+1) = 0.0
    enddo
    do k=2,nz
      do i=is,ie
        pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
        temp_int(i) = 0.5 * (t_f(i,j,k-1) + t_f(i,j,k))
        salin_int(i) = 0.5 * (s_f(i,j,k-1) + s_f(i,j,k))
      enddo
      call calculate_density_derivs(temp_int, salin_int, pres, &
             drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state)

      do i=is,ie
        alpha_dt = -1.0*drho_dt(i) * (t_f(i,j,k-1) - t_f(i,j,k))
        beta_ds  = drho_ds(i) * (s_f(i,j,k-1) - s_f(i,j,k))

        if ((alpha_dt > beta_ds) .and. (beta_ds > 0.0)) then  ! salt finger case
          rrho = min(alpha_dt/beta_ds,rrho0)
          diff_dd = 1.0 - ((rrho-1.0)/(rrho0-1.0))
          diff_dd = dsfmax*diff_dd*diff_dd*diff_dd
          kd_t_dd(i,k) = 0.7*diff_dd
          kd_s_dd(i,k) = diff_dd
        elseif ((alpha_dt < 0.) .and. (beta_ds < 0.) .and. (alpha_dt > beta_ds)) then ! diffusive convection
          rrho = alpha_dt/beta_ds
          diff_dd = kv_molecular*0.909*exp(4.6*exp(-0.54*(1/rrho-1)))
          prandtl = 0.15*rrho
          if (rrho > 0.5) prandtl = (1.85-0.85/rrho)*rrho
          kd_t_dd(i,k) = diff_dd
          kd_s_dd(i,k) = prandtl*diff_dd
        else
          kd_t_dd(i,k) = 0.0 ; kd_s_dd(i,k) = 0.0
        endif
      enddo
    enddo
  endif

end subroutine double_diffusion
!> This routine adds diffusion sustained by flow energy extracted by bottom drag.
subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, &
                                maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL)
  type(ocean_grid_type),                     intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                   intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                     intent(in)    :: tv
  type(forcing),                             intent(in)    :: fluxes
  type(vertvisc_type),                       intent(in)    :: visc
  integer,                                   intent(in)    :: j
  real, dimension(SZI_(G),SZK_(G)),          intent(in)    :: TKE_to_Kd, maxTKE
  integer, dimension(SZI_(G)),               intent(in)    :: kb
  type(set_diffusivity_cs),                  pointer       :: CS
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: Kd
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kd_int
  real, dimension(:,:,:),                    pointer       :: Kd_BBL

! This routine adds diffusion sustained by flow energy extracted by bottom drag.

  real, dimension(SZK_(G)+1) :: &
    Rint          ! coordinate density of an interface (kg/m3)
  real, dimension(SZI_(G)) :: &
    htot, &       ! total thickness above or below a layer, or the
                  ! integrated thickness in the BBL (meter)
    rho_htot, &   ! running integral with depth of density (kg/m2)
    gh_sum_top, & ! BBL value of g'h that can be supported by
                  ! the local ustar, times R0_g (kg/m2)
    rho_top, &    ! density at top of the BBL (kg/m3)
    tke, &        ! turbulent kinetic energy available to drive
                  ! bottom-boundary layer mixing in a layer (m3/s3)
    i2decay       ! inverse of twice the TKE decay scale (1/m)

  real    :: TKE_to_layer   ! TKE used to drive mixing in a layer (m3/s3)
  real    :: TKE_Ray        ! TKE from layer Rayleigh drag used to drive mixing in layer (m3/s3)
  real    :: TKE_here       ! TKE that goes into mixing in this layer (m3/s3)
  real    :: dRl, dRbot     ! temporaries holding density differences (kg/m3)
  real    :: cdrag_sqrt     ! square root of the drag coefficient (nondimensional)
  real    :: ustar_h        ! value of ustar at a thickness point (m/s)
  real    :: absf           ! average absolute Coriolis parameter around a thickness point (1/s)
  real    :: R0_g           ! Rho0 / G_Earth (kg s2 m-2)
  real    :: I_rho0         ! 1 / RHO0
  real    :: delta_Kd       ! increment to Kd from the bottom boundary layer mixing (m2/s)
  logical :: Rayleigh_drag  ! Set to true if Rayleigh drag velocities
                            ! defined in visc, on the assumption that this
                            ! extracted energy also drives diapycnal mixing.

  logical :: domore, do_i(szi_(g))
  logical :: do_diag_Kd_BBL

  integer :: i, k, is, ie, nz, i_rem, kb_min
  is = g%isc ; ie = g%iec ; nz = g%ke

  do_diag_kd_bbl = associated(kd_bbl)

  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return

  cdrag_sqrt = sqrt(cs%cdrag)
  tke_ray = 0.0 ; rayleigh_drag = .false.
  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.

  i_rho0 = 1.0/gv%Rho0
  r0_g = gv%Rho0/gv%g_Earth

  do k=2,nz ; rint(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo

  kb_min = max(gv%nk_rho_varies+1,2)

  ! The turbulence decay scale is 0.5*ustar/f from K&E & MOM_vertvisc.F90
  ! Any turbulence that makes it into the mixed layers is assumed
  ! to be relatively small and is discarded.
  do i=is,ie
    ustar_h = visc%ustar_BBL(i,j)
    if (ASSOCIATED(fluxes%ustar_tidal)) &
      ustar_h = ustar_h + fluxes%ustar_tidal(i,j)
    absf = 0.25*((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
                 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))))
    if ((ustar_h > 0.0) .and. (absf > 0.5*cs%IMax_decay*ustar_h))  then
      i2decay(i) = absf / ustar_h
    else
      ! The maximum decay scale should be something of order 200 m.
      ! If ustar_h = 0, this is land so this value doesn't matter.
      i2decay(i) = 0.5*cs%IMax_decay
    endif
    tke(i) = ((cs%BBL_effic * cdrag_sqrt) * &
              exp(-i2decay(i)*(gv%H_to_m*h(i,j,nz))) ) * &
             visc%TKE_BBL(i,j)

    if (ASSOCIATED(fluxes%TKE_tidal)) &
      tke(i) = tke(i) + fluxes%TKE_tidal(i,j) * i_rho0 * &
           (cs%BBL_effic * exp(-i2decay(i)*(gv%H_to_m*h(i,j,nz))))

    ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following
    ! Killworth & Edwards (1999) and Zilitikevich & Mironov (1996).
    ! Rho_top is determined by finding the density where
    ! integral(bottom, Z) (rho(z') - rho(Z)) dz' = rho_0 400 ustar^2 / g

    gh_sum_top(i) = r0_g * 400.0 * ustar_h**2

    do_i(i) = (g%mask2dT(i,j) > 0.5)
    htot(i) = gv%H_to_m*h(i,j,nz)
    rho_htot(i) = gv%Rlay(nz)*(gv%H_to_m*h(i,j,nz))
    rho_top(i) = gv%Rlay(1)
    if (cs%bulkmixedlayer .and. do_i(i)) rho_top(i) = gv%Rlay(kb(i)-1)
  enddo

  do k=nz-1,2,-1 ; domore = .false.
    do i=is,ie ; if (do_i(i)) then
      htot(i) = htot(i) + gv%H_to_m*h(i,j,k)
      rho_htot(i) = rho_htot(i) + gv%Rlay(k)*(gv%H_to_m*h(i,j,k))
      if (htot(i)*gv%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i))) then
        ! The top of the mixing is in the interface atop the current layer.
        rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i)
        do_i(i) = .false.
      elseif (k <= kb(i)) then ; do_i(i) = .false.
      else ; domore = .true. ; endif
    endif ; enddo
    if (.not.domore) exit
  enddo ! k-loop

  do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
  do k=nz-1,kb_min,-1
    i_rem = 0
    do i=is,ie ; if (do_i(i)) then
      if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
      i_rem = i_rem + 1  ! Count the i-rows that are still being worked on.
      !   Apply vertical decay of the turbulent energy.  This energy is
      ! simply lost.
      tke(i) = tke(i) * exp(-i2decay(i) * (gv%H_to_m*(h(i,j,k) + h(i,j,k+1))))

!      if (maxEnt(i,k) <= 0.0) cycle
      if (maxtke(i,k) <= 0.0) cycle

  ! This is an analytic integral where diffusity is a quadratic function of
  ! rho that goes asymptotically to 0 at Rho_top (vaguely following KPP?).
      if (tke(i) > 0.0) then
        if (rint(k) <= rho_top(i)) then
          tke_to_layer = tke(i)
        else
          drl = rint(k+1) - rint(k) ; drbot = rint(k+1) - rho_top(i)
          tke_to_layer = tke(i) * drl * (3.0*drbot*(rint(k) - rho_top(i)) +&
                                         drl**2) / drbot**3
        endif
      else ; tke_to_layer = 0.0 ; endif

      ! TKE_Ray has been initialized to 0 above.
      if (rayleigh_drag) tke_ray = 0.5*cs%BBL_effic * g%IareaT(i,j) * &
            ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
              g%areaCu(i,j)   * visc%Ray_u(i,j,k)   * u(i,j,k)**2) + &
             (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
              g%areaCv(i,j)   * visc%Ray_v(i,j,k)   * v(i,j,k)**2))

      if (tke_to_layer + tke_ray > 0.0) then
        if (cs%BBL_mixing_as_max) then
          if (tke_to_layer + tke_ray > maxtke(i,k)) &
              tke_to_layer = maxtke(i,k) - tke_ray

          tke(i) = tke(i) - tke_to_layer

          if (kd(i,j,k) < (tke_to_layer+tke_ray)*tke_to_kd(i,k)) then
            delta_kd = (tke_to_layer+tke_ray)*tke_to_kd(i,k) - kd(i,j,k)
            if ((cs%Kd_max >= 0.0) .and. (delta_kd > cs%Kd_max)) then
              delta_kd = cs%Kd_Max
              kd(i,j,k) = kd(i,j,k) + delta_kd
            else
              kd(i,j,k) = (tke_to_layer+tke_ray)*tke_to_kd(i,k)
            endif
            kd_int(i,j,k) = kd_int(i,j,k) + 0.5*delta_kd
            kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5*delta_kd
            if (do_diag_kd_bbl) then
              kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5*delta_kd
              kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5*delta_kd
            endif
          endif
        else
          if (kd(i,j,k) >= maxtke(i,k)*tke_to_kd(i,k)) then
            tke_here = 0.0
            tke(i) = tke(i) + tke_ray
          elseif (kd(i,j,k) + (tke_to_layer+tke_ray)*tke_to_kd(i,k) > &
                  maxtke(i,k)*tke_to_kd(i,k)) then
            tke_here = ((tke_to_layer+tke_ray) + kd(i,j,k)/tke_to_kd(i,k)) - &
                       maxtke(i,k)
            tke(i) = tke(i) - tke_here + tke_ray
          else
            tke_here = tke_to_layer + tke_ray
            tke(i) = tke(i) - tke_to_layer
          endif
          if (tke(i) < 0.0) tke(i) = 0.0 ! This should be unnecessary?

          if (tke_here > 0.0) then
            delta_kd = tke_here*tke_to_kd(i,k)
            if (cs%Kd_max >= 0.0) delta_kd = min(delta_kd, cs%Kd_max)
            kd(i,j,k) = kd(i,j,k) + delta_kd
            kd_int(i,j,k) = kd_int(i,j,k) + 0.5*delta_kd
            kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5*delta_kd
            if (do_diag_kd_bbl) then
              kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5*delta_kd
              kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5*delta_kd
            endif
          endif
        endif
      endif

      ! This may be risky - in the case that there are exactly zero
      ! velocities at 4 neighboring points, but nonzero velocities
      ! above the iterations would stop too soon. I don't see how this
      ! could happen in practice. RWH
      if ((tke(i)<= 0.0) .and. (tke_ray == 0.0)) then
        do_i(i) = .false. ; i_rem = i_rem - 1
      endif

    endif ; enddo
    if (i_rem == 0) exit
  enddo ! k-loop

end subroutine add_drag_diffusivity

!> Calculates a BBL diffusivity use a Prandtl number 1 diffusivitiy with a law of the
!! wall turbulent viscosity, up to a BBL height where the energy used for mixing has
!! consumed the mechanical TKE input.
subroutine add_lotw_bbl_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, &
                                    G, GV, CS, Kd, Kd_int, Kd_BBL)
  type(ocean_grid_type),                        intent(in)    :: G !< Grid structure
  type(verticalgrid_type),                      intent(in)    :: GV !< Vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)),    intent(in)    :: u !< u component of flow (m s-1)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)),    intent(in)    :: v !< v component of flow (m s-1)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),     intent(in)    :: h !< Layer thickness (m or kg m-2)
  type(thermo_var_ptrs),                        intent(in)    :: tv !< Thermodynamic variables structure
  type(forcing),                                intent(in)    :: fluxes !< Surface fluxes structure
  type(vertvisc_type),                          intent(in)    :: visc !< Vertical viscosity structure
  integer,                                      intent(in)    :: j !< j-index of row to work on
  real, dimension(SZI_(G),SZK_(G)+1),           intent(in)    :: N2_int !< Square of Brunt-Vaisala at interfaces (s-2)
  type(set_diffusivity_cs),                     pointer       :: CS !< Diffusivity control structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),     intent(inout) :: Kd !< Layer net diffusivity (m2 s-1)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1),   intent(inout) :: Kd_int !< Interface net diffusivity (m2 s-1)
  real, dimension(:,:,:),                       pointer       :: Kd_BBL !< Interface BBL diffusivity (m2 s-1)

  ! Local variables
  real :: TKE_column       ! net TKE input into the column (m3 s-3)
  real :: TKE_to_layer     ! TKE used to drive mixing in a layer (m3 s-3)
  real :: TKE_Ray          ! TKE from a layer Rayleigh drag used to drive mixing in that layer (m3 s-3)
  real :: TKE_remaining    ! remaining TKE available for mixing in this layer and above (m3 s-3)
  real :: TKE_consumed     ! TKE used for mixing in this layer (m3 s-3)
  real :: TKE_Kd_wall      ! TKE associated with unlimited law of the wall mixing (m3 s-3)
  real :: cdrag_sqrt       ! square root of the drag coefficient (nondimensional)
  real :: ustar            ! value of ustar at a thickness point (m/s)
  real :: ustar2           ! square of ustar, for convenience (m2/s2)
  real :: absf             ! average absolute value of Coriolis parameter around a thickness point (1/sec)
  real :: dh, dhm1         ! thickness of layers k and k-1, respecitvely (meter)
  real :: z                ! distance to interface k from bottom (meter)
  real :: D_minus_z        ! distance to interface k from surface (meter)
  real :: total_thickness  ! total thickness of water column (meter)
  real :: Idecay           ! inverse of decay scale used for "Joule heating" loss of TKE with height (1/m)
  real :: Kd_wall          ! Law of the wall diffusivity (m2/s)
  real :: Kd_lower         ! diffusivity for lower interface (m2/sec)
  real :: ustar_D          ! u* x D  (m2/s)
  real :: I_Rho0           ! 1 / rho0
  real :: N2_min           ! Minimum value of N2 to use in calculation of TKE_Kd_wall (1/s2)
  logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on
                           ! the assumption that this extracted energy also drives diapycnal mixing.
  integer :: i, k, km1
  real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant)
  logical :: do_diag_Kd_BBL

  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return
  do_diag_kd_bbl = associated(kd_bbl)

  n2_min = 0.
  if (cs%LOTW_BBL_use_omega) n2_min = (cs%omega**2)

  ! Determine whether to add Rayleigh drag contribution to TKE
  rayleigh_drag = .false.
  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.
  i_rho0 = 1.0/gv%Rho0
  cdrag_sqrt = sqrt(cs%cdrag)

  tke_ray = 0. ! In case Rayleigh_drag is not used.
  do i=g%isc,g%iec ! Developed in single-column mode

    ! Column-wise parameters.
    absf = 0.25*((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
                 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1)))) ! Non-zero on equator!

    ! u* at the bottom, in m s-1.
    ustar = visc%ustar_BBL(i,j)
    ustar2 = ustar**2
    ! In add_drag_diffusivity(), fluxes%ustar_tidal is added in. This might be double counting
    ! since ustar_BBL should already include all contributions to u*? -AJA
    if (ASSOCIATED(fluxes%ustar_tidal)) ustar = ustar + fluxes%ustar_tidal(i,j)

    ! The maximum decay scale should be something of order 200 m. We use the smaller of u*/f and
    ! (IMax_decay)^-1 as the decay scale. If ustar = 0, this is land so this value doesn't matter.
    idecay = cs%IMax_decay
    if ((ustar > 0.0) .and. (absf > cs%IMax_decay*ustar)) idecay = absf / ustar

    ! Energy input at the bottom, in m3 s-3.
    ! (Note that visc%TKE_BBL is in m3 s-3, set in set_BBL_TKE().)
    ! I am still unsure about sqrt(cdrag) in this expressions - AJA
    tke_column = cdrag_sqrt * visc%TKE_BBL(i,j)
    ! Add in tidal dissipation energy at the bottom, in m3 s-3.
    ! Note that TKE_tidal is in W m-2.
    if (ASSOCIATED(fluxes%TKE_tidal)) tke_column = tke_column + fluxes%TKE_tidal(i,j) * i_rho0
    tke_column = cs%BBL_effic * tke_column ! Only use a fraction of the mechanical dissipation for mixing.

    tke_remaining = tke_column
    total_thickness = ( sum(h(i,j,:)) + gv%H_subroundoff )* gv%H_to_m ! Total column thickness, in m.
    ustar_d = ustar * total_thickness
    z = 0.
    kd_lower = 0. ! Diffusivity on bottom boundary.

    ! Work upwards from the bottom, accumulating work used until it exceeds the available TKE input
    ! at the bottom.
    do k=g%ke,2,-1
      dh = gv%H_to_m * h(i,j,k) ! Thickness of this level in m.
      km1 = max(k-1, 1)
      dhm1 = gv%H_to_m * h(i,j,km1) ! Thickness of level above in m.

      ! Add in additional energy input from bottom-drag against slopes (sides)
      if (rayleigh_drag) tke_remaining = tke_remaining + 0.5*cs%BBL_effic * g%IareaT(i,j) * &
            ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
              g%areaCu(i,j)   * visc%Ray_u(i,j,k)   * u(i,j,k)**2) + &
             (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
              g%areaCv(i,j)   * visc%Ray_v(i,j,k)   * v(i,j,k)**2))

      ! Exponentially decay TKE across the thickness of the layer.
      ! This is energy loss in addition to work done as mixing, apparently to Joule heating.
      tke_remaining = exp(-idecay*dh) * tke_remaining

      z = z + h(i,j,k)*gv%H_to_m ! Distance between upper interface of layer and the bottom, in m.
      d_minus_z = max(total_thickness - z, 0.) ! Thickness above layer, m.

      ! Diffusivity using law of the wall, limited by rotation, at height z, in m2/s.
      ! This calculation is at the upper interface of the layer
      if ( ustar_d + absf * ( z * d_minus_z ) == 0.) then
        kd_wall = 0.
      else
        kd_wall = ( ( von_karm * ustar2 ) * ( z * d_minus_z ) )/( ustar_d + absf * ( z * d_minus_z ) )
      endif

      ! TKE associated with Kd_wall, in m3 s-2.
      ! This calculation if for the volume spanning the interface.
      tke_kd_wall = kd_wall * 0.5 * (dh + dhm1) * max(n2_int(i,k), n2_min)

      ! Now bound Kd such that the associated TKE is no greater than available TKE for mixing.
      if (tke_kd_wall>0.) then
        tke_consumed = min(tke_kd_wall, tke_remaining)
        kd_wall = (tke_consumed/tke_kd_wall) * kd_wall ! Scale Kd so that only TKE_consumed is used.
      else
        ! Either N2=0 or dh = 0.
        if (tke_remaining>0.) then
          kd_wall = cs%Kd_max
        else
          kd_wall = 0.
        endif
        tke_consumed = 0.
      endif

      ! Now use up the appropriate about of TKE associated with the diffusivity chosen
      tke_remaining = tke_remaining - tke_consumed ! Note this will be non-negative

      ! Add this BBL diffusivity to the model net diffusivity.
      kd_int(i,j,k) = kd_int(i,j,k) + kd_wall
      kd(i,j,k) = kd(i,j,k) + 0.5*(kd_wall + kd_lower)
      kd_lower = kd_wall ! Store for next level up.
      if (do_diag_kd_bbl) kd_bbl(i,j,k) = kd_wall
    enddo ! k
  enddo ! i

end subroutine add_lotw_bbl_diffusivity
!> This routine adds effects of mixed layer radiation to the layer diffusivities.
subroutine add_mlrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int)
  type(ocean_grid_type),                    intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(forcing),                            intent(in)    :: fluxes
  integer,                                  intent(in)    :: j
  type(set_diffusivity_cs),                 pointer       :: CS
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd
  real, dimension(SZI_(G),SZK_(G)),         intent(in)    :: TKE_to_Kd
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int

! This routine adds effects of mixed layer radiation to the layer diffusivities.

  real, dimension(SZI_(G)) ::         &
                         h_ml,        &
                         TKE_ml_flux, &
                         I_decay,     &
                         Kd_mlr_ml

  real :: f_sq, h_ml_sq, ustar_sq, Kd_mlr, C1_6
  real :: Omega2            ! rotation rate squared (1/s2)
  real :: z1                ! layer thickness times I_decay (nondim)
  real :: dzL               ! thickness converted to meter
  real :: I_decay_len2_TKE  ! squared inverse decay lengthscale for
                            ! TKE, as used in the mixed layer code (1/m2)
  real :: h_neglect         ! negligibly small thickness (meter)

  logical :: do_any, do_i(szi_(g))
  integer :: i, k, is, ie, nz, kml
  is = g%isc ; ie = g%iec ; nz = g%ke

  omega2    = cs%Omega**2
  c1_6      = 1.0 / 6.0
  kml       = gv%nkml
  h_neglect = gv%H_subroundoff*gv%H_to_m

  if (.not.cs%ML_radiation) return

  do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
  do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + gv%H_to_m*h(i,j,k) ; enddo ; enddo

  do i=is,ie ; if (do_i(i)) then
    if (cs%ML_omega_frac >= 1.0) then
      f_sq = 4.0*omega2
    else
      f_sq = 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))
      if (cs%ML_omega_frac > 0.0) &
        f_sq = cs%ML_omega_frac*4.0*omega2 + (1.0-cs%ML_omega_frac)*f_sq
    endif

    ustar_sq = max(fluxes%ustar(i,j), cs%ustar_min)**2

    tke_ml_flux(i) = (cs%mstar*cs%ML_rad_coeff)*(ustar_sq*fluxes%ustar(i,j))
    i_decay_len2_tke = cs%TKE_decay**2 * (f_sq / ustar_sq)

    if (cs%ML_rad_TKE_decay) &
      tke_ml_flux(i) = tke_ml_flux(i) * exp(-h_ml(i) * sqrt(i_decay_len2_tke))

    ! Calculate the inverse decay scale
    h_ml_sq = (cs%ML_rad_efold_coeff * (h_ml(i)+h_neglect))**2
    i_decay(i) = sqrt((i_decay_len2_tke * h_ml_sq + 1.0) / h_ml_sq)

    ! Average the dissipation layer kml+1, using
    ! a more accurate Taylor series approximations for very thin layers.
    z1 = (gv%H_to_m*h(i,j,kml+1)) * i_decay(i)
    if (z1 > 1e-5) then
      kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,kml+1)) * &
               (1.0 - exp(-z1))
    else
      kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,kml+1)) * &
               (z1 * (1.0 - z1 * (0.5 - c1_6*z1)))
    endif
    kd_mlr_ml(i) = min(kd_mlr,cs%ML_rad_kd_max)
    tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
  endif ; enddo

  do k=1,kml+1 ; do i=is,ie ; if (do_i(i)) then
    kd(i,j,k) = kd(i,j,k) + kd_mlr_ml(i)
  endif ; enddo ; enddo
  if (present(kd_int)) then
    do k=2,kml+1 ; do i=is,ie ; if (do_i(i)) then
      kd_int(i,j,k) = kd_int(i,j,k) + kd_mlr_ml(i)
    endif ; enddo ; enddo
    if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then
      kd_int(i,j,kml+2) = kd_int(i,j,kml+2) + 0.5*kd_mlr_ml(i)
    endif ; enddo ; endif
  endif

  do k=kml+2,nz-1
    do_any = .false.
    do i=is,ie ; if (do_i(i)) then
      dzl = gv%H_to_m*h(i,j,k) ;  z1 = dzl*i_decay(i)
      if (z1 > 1e-5) then
        kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * &
                 ((1.0 - exp(-z1)) / dzl)
      else
        kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * &
                 (i_decay(i) * (1.0 - z1 * (0.5 - c1_6*z1)))
      endif
      kd_mlr = min(kd_mlr,cs%ML_rad_kd_max)
      kd(i,j,k) =  kd(i,j,k) + kd_mlr
      if (present(kd_int)) then
        kd_int(i,j,k) = kd_int(i,j,k) + 0.5*kd_mlr
        kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5*kd_mlr
      endif

      tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
      if (tke_ml_flux(i) * i_decay(i) < 0.1*cs%Kd_min*omega2) then
        do_i(i) = .false.
      else ; do_any = .true. ; endif
    endif ; enddo
    if (.not.do_any) exit
  enddo

end subroutine add_mlrad_diffusivity

  !> This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities.
  !! The mechanisms considered are (1) local dissipation of internal waves generated by the
  !! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating
  !! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves.
  !! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction,
  !! Froude-number-depending breaking, PSI, etc.).
subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, &
                                    dd, N2_lay, Kd, Kd_int )
  type(ocean_grid_type),                    intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real, dimension(SZI_(G)),                 intent(in)    :: N2_bot
  real, dimension(SZI_(G),SZK_(G)),         intent(in)    :: N2_lay
  integer,                                  intent(in)    :: j
  real, dimension(SZI_(G),SZK_(G)),         intent(in)    :: TKE_to_Kd, max_TKE
  type(set_diffusivity_cs),                 pointer       :: CS
  type(diffusivity_diags),                  intent(inout) :: dd
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int

  ! This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities.
  ! The mechanisms considered are (1) local dissipation of internal waves generated by the
  ! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating
  ! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves.
  ! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction,
  ! Froude-number-depending breaking, PSI, etc.).

  real, dimension(SZI_(G)) :: &
    htot,             & ! total thickness above or below a layer, or the
                        ! integrated thickness in the BBL (meter)
    htot_wkb,         & ! distance from top to bottom (meter) WKB scaled
    tke_itidal_bot,   & ! internal tide TKE at ocean bottom (m3/s3)
    tke_niku_bot,     & ! lee-wave TKE at ocean bottom (m3/s3)
    tke_lowmode_bot,  & ! internal tide TKE at ocean bottom lost from all remote low modes (m3/s3) (BDM)
    inv_int,          & ! inverse of TKE decay for int tide over the depth of the ocean (nondim)
    inv_int_lee,      & ! inverse of TKE decay for lee waves over the depth of the ocean (nondim)
    inv_int_low,      & ! inverse of TKE decay for low modes over the depth of the ocean (nondim) (BDM)
    z0_polzin,        & ! TKE decay scale in Polzin formulation (meter)
    z0_polzin_scaled, & ! TKE decay scale in Polzin formulation (meter)
                        ! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z
                        ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz)
                        ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz
    n2_meanz,         & ! vertically averaged squared buoyancy frequency (1/s2) for WKB scaling
    tke_itidal_rem,   & ! remaining internal tide TKE (from barotropic source)
    tke_niku_rem,     & ! remaining lee-wave TKE
    tke_lowmode_rem,  & ! remaining internal tide TKE (from propagating low mode source) (BDM)
    tke_frac_top,     & ! fraction of bottom TKE that should appear at top of a layer (nondim)
    tke_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer (nondim)
    tke_frac_top_lowmode, &
                        ! fraction of bottom TKE that should appear at top of a layer (nondim) (BDM)
    z_from_bot,       & ! distance from bottom (meter)
    z_from_bot_wkb      ! distance from bottom (meter), WKB scaled

  real :: I_rho0        ! 1 / RHO0, (m3/kg)
  real :: Kd_add        ! diffusivity to add in a layer (m2/sec)
  real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) (m3/s3)
  real :: TKE_Niku_lay  ! lee-wave TKE imparted to a layer (m3/s3)
  real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) (m3/s3) (BDM)
  real :: frac_used     ! fraction of TKE that can be used in a layer (nondim)
  real :: Izeta         ! inverse of TKE decay scale (1/meter)
  real :: Izeta_lee     ! inverse of TKE decay scale for lee waves (1/meter)
  real :: z0_psl        ! temporary variable with units of meter
  real :: TKE_lowmode_tot ! TKE from all low modes (W/m2) (BDM)


  logical :: use_Polzin, use_Simmons
  integer :: i, k, is, ie, nz
  integer :: a, fr, m
  is = g%isc ; ie = g%iec ; nz = g%ke

  if (.not.(cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation)) return

  do i=is,ie ; htot(i) = 0.0 ; inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0 ;enddo
  do k=1,nz ; do i=is,ie
    htot(i) = htot(i) + gv%H_to_m*h(i,j,k)
  enddo ; enddo

  i_rho0 = 1.0/gv%Rho0

  use_polzin = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09)) .or. &
                (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09)) .or. &
                (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09)))
  use_simmons = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == stlaurent_02)) .or. &
                 (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == stlaurent_02)) .or. &
                 (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == stlaurent_02)))

  ! Calculate parameters for vertical structure of dissipation
  ! Simmons:
  if ( use_simmons ) then
    izeta = 1.0 / max(cs%Int_tide_decay_scale, gv%H_subroundoff*gv%H_to_m)
    izeta_lee = 1.0 / max(cs%Int_tide_decay_scale*cs%Decay_scale_factor_lee, &
                          gv%H_subroundoff*gv%H_to_m)
    do i=is,ie
      cs%Nb(i,j) = sqrt(n2_bot(i))
      if (associated(dd%N2_bot)) dd%N2_bot(i,j) = n2_bot(i)
      if ( cs%Int_tide_dissipation ) then
        if (izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
        endif
      endif
      if ( cs%Lee_wave_dissipation ) then
        if (izeta_lee*htot(i) > 1.0e-14) then  ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int_lee(i) = 1.0 / (1.0 - exp(-izeta_lee*htot(i)))
        endif
      endif
      if ( cs%Lowmode_itidal_dissipation) then
        if (izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int_low(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
        endif
      endif
      z_from_bot(i) = gv%H_to_m*h(i,j,nz)
    enddo
  endif ! Simmons

  ! Polzin:
  if ( use_polzin ) then
    ! WKB scaling of the vertical coordinate
    do i=is,ie ; n2_meanz(i)=0.0 ; enddo
    do k=1,nz ; do i=is,ie
      n2_meanz(i) = n2_meanz(i) + n2_lay(i,k)*gv%H_to_m*h(i,j,k)
    enddo ; enddo
    do i=is,ie
      n2_meanz(i) = n2_meanz(i) / (htot(i) + gv%H_subroundoff*gv%H_to_m)
      if (associated(dd%N2_meanz))  dd%N2_meanz(i,j) = n2_meanz(i)
    enddo

    ! WKB scaled z*(z=H) z* at the surface using the modified Polzin WKB scaling
    do i=is,ie ; htot_wkb(i) = htot(i) ; enddo
!    do i=is,ie ; htot_WKB(i) = 0.0 ; enddo
!    do k=1,nz ; do i=is,ie
!      htot_WKB(i) = htot_WKB(i) + GV%H_to_m*h(i,j,k)*N2_lay(i,k) / N2_meanz(i)
!    enddo ; enddo
    ! htot_WKB(i) = htot(i) ! Nearly equivalent and simpler

    do i=is,ie
      cs%Nb(i,j) = sqrt(n2_bot(i))
      if ((cs%tideamp(i,j) > 0.0) .and. &
          (cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 > 1.0e-14) ) then
        z0_polzin(i) = cs%Polzin_decay_scale_factor * cs%Nu_Polzin * &
                       cs%Nbotref_Polzin**2 * cs%tideamp(i,j) / &
                     ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 )
        if (z0_polzin(i) < cs%Polzin_min_decay_scale) &
          z0_polzin(i) = cs%Polzin_min_decay_scale
        if (n2_meanz(i) > 1.0e-14  ) then
          z0_polzin_scaled(i) = z0_polzin(i)*cs%Nb(i,j)**2 / n2_meanz(i)
        else
          z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
        endif
        if (z0_polzin_scaled(i) > (cs%Polzin_decay_scale_max_factor * htot(i)) ) &
          z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
      else
        z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
        z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
      endif

      if (associated(dd%Polzin_decay_scale)) &
        dd%Polzin_decay_scale(i,j) = z0_polzin(i)
      if (associated(dd%Polzin_decay_scale_scaled)) &
        dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i)
      if (associated(dd%N2_bot)) dd%N2_bot(i,j) = cs%Nb(i,j)*cs%Nb(i,j)

      if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
        ! For the Polzin formulation, this if loop prevents the vertical
        ! flux of energy dissipation from having NaN values
        if (htot_wkb(i) > 1.0e-14) then
          inv_int(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1
        endif
      endif
      if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09) ) then
        ! For the Polzin formulation, this if loop prevents the vertical
        ! flux of energy dissipation from having NaN values
        if (htot_wkb(i) > 1.0e-14) then
          inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / htot_wkb(i) ) + 1
        endif
      endif
      if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
        ! For the Polzin formulation, this if loop prevents the vertical
        ! flux of energy dissipation from having NaN values
        if (htot_wkb(i) > 1.0e-14) then
          inv_int_low(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1
        endif
      endif

      z_from_bot(i) = gv%H_to_m*h(i,j,nz)
      ! Use the new formulation for WKB scaling.  N2 is referenced to its
      ! vertical mean.
      if (n2_meanz(i) > 1.0e-14 ) then
        z_from_bot_wkb(i) = gv%H_to_m*h(i,j,nz)*n2_lay(i,nz) / n2_meanz(i)
      else ; z_from_bot_wkb(i) = 0 ; endif
    enddo
  endif  ! Polzin

  ! Calculate/get dissipation values at bottom
  ! Both Polzin and Simmons:
  do i=is,ie
    ! Dissipation of locally trapped internal tide (non-propagating high modes)
    tke_itidal_bot(i) = min(cs%TKE_itidal(i,j)*cs%Nb(i,j),cs%TKE_itide_max)
    if (associated(dd%TKE_itidal_used)) &
      dd%TKE_itidal_used(i,j) = tke_itidal_bot(i)
    tke_itidal_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_itides) * tke_itidal_bot(i)
    ! Dissipation of locally trapped lee waves
    tke_niku_bot(i) = 0.0
    if (cs%Lee_wave_dissipation) then
      tke_niku_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_lee) * cs%TKE_Niku(i,j)
    endif
    ! Dissipation of propagating internal tide (baroclinic low modes; rays) (BDM)
    tke_lowmode_tot    = 0.0
    tke_lowmode_bot(i) = 0.0
    if (cs%Lowmode_itidal_dissipation) then
      ! get loss rate due to wave drag on low modes (already multiplied by q)
      call get_lowmode_loss(i,j,g,cs%int_tide_CSp,"WaveDrag",tke_lowmode_tot)
      tke_lowmode_bot(i) = cs%Mu_itides * i_rho0 * tke_lowmode_tot
    endif
    ! Vertical energy flux at bottom
    tke_itidal_rem(i)  = inv_int(i)     * tke_itidal_bot(i)
    tke_niku_rem(i)    = inv_int_lee(i) * tke_niku_bot(i)
    tke_lowmode_rem(i) = inv_int_low(i) * tke_lowmode_bot(i)

    if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,nz) = tke_itidal_rem(i) !why is this here? BDM
  enddo

  ! Estimate the work that would be done by mixing in each layer.
  ! Simmons:
  if ( use_simmons ) then
    do k=nz-1,2,-1 ; do i=is,ie
      if (max_tke(i,k) <= 0.0) cycle
      z_from_bot(i) = z_from_bot(i) + gv%H_to_m*h(i,j,k)

      ! Fraction of bottom flux predicted to reach top of this layer
      tke_frac_top(i)         = inv_int(i)     * exp(-izeta * z_from_bot(i))
      tke_frac_top_lee(i)     = inv_int_lee(i) * exp(-izeta_lee * z_from_bot(i))
      tke_frac_top_lowmode(i) = inv_int_low(i) * exp(-izeta * z_from_bot(i))

      ! Actual influx at bottom of layer minus predicted outflux at top of layer to give
      ! predicted power expended
      tke_itide_lay   = tke_itidal_rem(i)  - tke_itidal_bot(i) * tke_frac_top(i)
      tke_niku_lay    = tke_niku_rem(i)    - tke_niku_bot(i)   * tke_frac_top_lee(i)
      tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)* tke_frac_top_lowmode(i)

      ! Actual power expended may be less than predicted if stratification is weak; adjust
      if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > max_tke(i,k)) then
         frac_used = max_tke(i,k) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
         tke_itide_lay   = frac_used * tke_itide_lay
         tke_niku_lay    = frac_used * tke_niku_lay
         tke_lowmode_lay = frac_used * tke_lowmode_lay
      endif

      ! Calculate vertical flux available to bottom of layer above
      tke_itidal_rem(i)  = tke_itidal_rem(i)  - tke_itide_lay
      tke_niku_rem(i)    = tke_niku_rem(i)    - tke_niku_lay
      tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay

      ! Convert power to diffusivity
      kd_add  = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)

      if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
      kd(i,j,k) = kd(i,j,k) + kd_add

      if (present(kd_int)) then
        kd_int(i,j,k)   = kd_int(i,j,k)   + 0.5*kd_add
        kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5*kd_add
      endif

      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        ! If at layers, dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_itide_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1)  dd%Kd_itidal(i,j,k)   = dd%Kd_itidal(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_Itidal_work)) &
        dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
      if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)

      if (associated(dd%Kd_Niku)) then
        ! If at layers, dd%Kd_Niku(i,j,K) is just TKE_to_Kd(i,k) * TKE_Niku_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_niku_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1) dd%Kd_Niku(i,j,k)    = dd%Kd_Niku(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
      endif
!     if (associated(dd%Kd_Niku)) dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
      if (associated(dd%Kd_Niku_work)) &
        dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay

      if (associated(dd%Kd_lowmode)) then
        ! If at layers, dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_lowmode_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1)  dd%Kd_lowmode(i,j,k)   = dd%Kd_lowmode(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_lowmode_work)) &
        dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
      if (associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)

    enddo ; enddo ;
  endif ! Simmons

  ! Polzin:
  if ( use_polzin ) then
    do k=nz-1,2,-1 ; do i=is,ie
      if (max_tke(i,k) <= 0.0) cycle
      z_from_bot(i) = z_from_bot(i) + gv%H_to_m*h(i,j,k)
      if (n2_meanz(i) > 1.0e-14 ) then
        z_from_bot_wkb(i) = z_from_bot_wkb(i) + gv%H_to_m*h(i,j,k)*n2_lay(i,k)/n2_meanz(i)
      else ; z_from_bot_wkb(i) = 0 ; endif

      ! Fraction of bottom flux predicted to reach top of this layer
      tke_frac_top(i)     = ( inv_int(i) * z0_polzin_scaled(i) ) / &
                            ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
      z0_psl = z0_polzin_scaled(i)*cs%Decay_scale_factor_lee
      tke_frac_top_lee(i) = (inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_wkb(i))
      tke_frac_top_lowmode(i) = ( inv_int_low(i) * z0_polzin_scaled(i) ) / &
                            ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )

      ! Actual influx at bottom of layer minus predicted outflux at top of layer to give
      ! predicted power expended
      tke_itide_lay   = tke_itidal_rem(i)  - tke_itidal_bot(i) *tke_frac_top(i)
      tke_niku_lay    = tke_niku_rem(i)    - tke_niku_bot(i)   * tke_frac_top_lee(i)
      tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)*tke_frac_top_lowmode(i)

      ! Actual power expended may be less than predicted if stratification is weak; adjust
      if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > max_tke(i,k)) then
        frac_used = max_tke(i,k) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
        tke_itide_lay   = frac_used * tke_itide_lay
        tke_niku_lay    = frac_used * tke_niku_lay
        tke_lowmode_lay = frac_used * tke_lowmode_lay
      endif

      ! Calculate vertical flux available to bottom of layer above
      tke_itidal_rem(i)  = tke_itidal_rem(i)  - tke_itide_lay
      tke_niku_rem(i)    = tke_niku_rem(i)    - tke_niku_lay
      tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay

      ! Convert power to diffusivity
      kd_add  = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)

      if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
      kd(i,j,k) = kd(i,j,k) + kd_add

      if (present(kd_int)) then
        kd_int(i,j,k)   = kd_int(i,j,k)   + 0.5*kd_add
        kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5*kd_add
      endif

      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        ! If at layers, this is just dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_itide_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1)  dd%Kd_itidal(i,j,k)   = dd%Kd_itidal(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_Itidal_work)) &
        dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
      if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)

      if (associated(dd%Kd_Niku)) then
        ! If at layers, this is just dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_niku_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1) dd%Kd_Niku(i,j,k)    = dd%Kd_Niku(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
      endif
   !  if (associated(dd%Kd_Niku)) dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
      if (associated(dd%Kd_Niku_work)) dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay

      if (associated(dd%Kd_lowmode)) then
        ! If at layers, dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_lowmode_lay
        if (cs%Kd_max >= 0.0) kd_add = min(kd_add, cs%Kd_max)
        if (k>1)  dd%Kd_lowmode(i,j,k)   = dd%Kd_lowmode(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_lowmode_work)) &
        dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
      if (associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)

    enddo ; enddo;
  endif ! Polzin

end subroutine add_int_tide_diffusivity
!> This subroutine calculates several properties related to bottom
!! boundary layer turbulence.
subroutine set_bbl_tke(u, v, h, fluxes, visc, G, GV, CS)
  type(ocean_grid_type),                     intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                   intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(forcing),                             intent(in)    :: fluxes
  type(vertvisc_type),                       intent(inout) :: visc
  type(set_diffusivity_cs),                  pointer       :: CS

  ! This subroutine calculates several properties related to bottom
  ! boundary layer turbulence.

  real, dimension(SZI_(G)) :: &
    htot          ! total thickness above or below a layer, or the
                  ! integrated thickness in the BBL (meter)

  real, dimension(SZIB_(G)) :: &
    uhtot, &      ! running integral of u in the BBL (m2/s)
    ustar, &      ! bottom boundary layer turbulence speed (m/s)
    u2_bbl        ! square of the mean zonal velocity in the BBL (m2/s2)

  real :: vhtot(szi_(g)) ! running integral of v in the BBL (m2/sec)

  real, dimension(SZI_(G),SZJB_(G)) :: &
    vstar, & ! ustar at at v-points in 2 j-rows (m/s)
    v2_bbl   ! square of average meridional velocity in BBL (m2/s2)

  real :: cdrag_sqrt  ! square root of the drag coefficient (nondim)
  real :: hvel        ! thickness at velocity points (meter)

  logical :: domore, do_i(szi_(g))
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (.not.associated(cs)) call mom_error(fatal,"set_BBL_TKE: "//&
         "Module must be initialized before it is used.")

  if (.not.cs%bottomdraglaw .or. (cs%BBL_effic<=0.0)) then
    if (associated(visc%ustar_BBL)) then
      do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo
    endif
    if (associated(visc%TKE_BBL)) then
      do j=js,je ; do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ; enddo ; enddo
    endif
    return
  endif

  cdrag_sqrt = sqrt(cs%cdrag)

!$OMP parallel default(none) shared(cdrag_sqrt,is,ie,js,je,nz,visc,CS,G,GV,vstar,h,v, &
!$OMP                               v2_bbl,u) &
!$OMP                       private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl)
!$OMP do
  do j=js-1,je
    ! Determine ustar and the square magnitude of the velocity in the
    ! bottom boundary layer. Together these give the TKE source and
    ! vertical decay scale.
    do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_v(i,j) > 0.0)) then
      do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0
      vstar(i,j) = visc%kv_bbl_v(i,j)/(cdrag_sqrt*visc%bbl_thick_v(i,j))
    else
      do_i(i) = .false. ; vstar(i,j) = 0.0 ; htot(i) = 0.0
    endif ; enddo
    do k=nz,1,-1
      domore = .false.
      do i=is,ie ; if (do_i(i)) then
        hvel = 0.5*gv%H_to_m*(h(i,j,k) + h(i,j+1,k))
        if ((htot(i) + hvel) >= visc%bbl_thick_v(i,j)) then
          vhtot(i) = vhtot(i) + (visc%bbl_thick_v(i,j) - htot(i))*v(i,j,k)
          htot(i) = visc%bbl_thick_v(i,j)
          do_i(i) = .false.
        else
          vhtot(i) = vhtot(i) + hvel*v(i,j,k)
          htot(i) = htot(i) + hvel
          domore = .true.
        endif
      endif ; enddo
      if (.not.domore) exit
    enddo
    do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (htot(i) > 0.0)) then
      v2_bbl(i,j) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i))
    else
      v2_bbl(i,j) = 0.0
    endif ; enddo
  enddo
!$OMP do
  do j=js,je
    do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_u(i,j) > 0.0))  then
      do_i(i) = .true. ; uhtot(i) = 0.0 ; htot(i) = 0.0
      ustar(i) = visc%kv_bbl_u(i,j)/(cdrag_sqrt*visc%bbl_thick_u(i,j))
    else
      do_i(i) = .false. ; ustar(i) = 0.0 ; htot(i) = 0.0
    endif ; enddo
    do k=nz,1,-1 ; domore = .false.
      do i=is-1,ie ; if (do_i(i)) then
        hvel = 0.5*gv%H_to_m*(h(i,j,k) + h(i+1,j,k))
        if ((htot(i) + hvel) >= visc%bbl_thick_u(i,j)) then
          uhtot(i) = uhtot(i) + (visc%bbl_thick_u(i,j) - htot(i))*u(i,j,k)
          htot(i) = visc%bbl_thick_u(i,j)
          do_i(i) = .false.
        else
          uhtot(i) = uhtot(i) + hvel*u(i,j,k)
          htot(i) = htot(i) + hvel
          domore = .true.
        endif
      endif ; enddo
      if (.not.domore) exit
    enddo
    do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (htot(i) > 0.0)) then
      u2_bbl(i) = (uhtot(i)*uhtot(i))/(htot(i)*htot(i))
    else
      u2_bbl(i) = 0.0
    endif ; enddo

    do i=is,ie
      visc%ustar_BBL(i,j) = sqrt(0.5*g%IareaT(i,j) * &
                ((g%areaCu(i-1,j)*(ustar(i-1)*ustar(i-1)) + &
                  g%areaCu(i,j)*(ustar(i)*ustar(i))) + &
                 (g%areaCv(i,j-1)*(vstar(i,j-1)*vstar(i,j-1)) + &
                  g%areaCv(i,j)*(vstar(i,j)*vstar(i,j))) ) )
      visc%TKE_BBL(i,j) = (((g%areaCu(i-1,j)*(ustar(i-1)*u2_bbl(i-1)) + &
                    g%areaCu(i,j) * (ustar(i)*u2_bbl(i))) + &
                   (g%areaCv(i,j-1)*(vstar(i,j-1)*v2_bbl(i,j-1)) + &
                    g%areaCv(i,j) * (vstar(i,j)*v2_bbl(i,j))))*g%IareaT(i,j))
    enddo
  enddo
!$OMP end parallel

end subroutine set_bbl_tke

subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0)
  type(ocean_grid_type),            intent(in)   :: G  !< The ocean's grid structure.
  type(verticalgrid_type),          intent(in)   :: GV !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(in)   :: h  !< Layer thicknesses, in H (usually m
                                                       !! or kg m-2).
  type(thermo_var_ptrs),            intent(in)   :: tv !< Structure containing pointers to any
                                                       !! available thermodynamic fields; absent
                                                       !! fields have NULL ptrs.
  integer, dimension(SZI_(G)),      intent(in)   :: kb !< Index of lightest layer denser than the
                                                       !! buffer layer.
  type(set_diffusivity_cs),         pointer      :: CS !< Control structure returned by previous
                                                       !! call to diabatic_entrain_init.
  integer,                          intent(in)   :: j  !< Meridional index upon which to work.
  real, dimension(SZI_(G),SZK_(G)), intent(out)  :: ds_dsp1 !< Coordinate variable (sigma-2)
                                                       !! difference across an interface divided by
                                                       !! the difference across the interface below
                                                       !! it (nondimensional)
  real, dimension(SZI_(G),SZK_(G)), &
                          optional, intent(in)   :: rho_0 !< Layer potential densities relative to
                                                       !! surface press (kg/m3).

! Arguments:
!  (in)      h       - layer thickness (meter)
!  (in)      tv      - structure containing pointers to any available
!                      thermodynamic fields; absent fields have NULL ptrs
!  (in)      kb      - index of lightest layer denser than the buffer layer
!  (in)      G       - ocean grid structure
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS      - control structure returned by previous call to diabatic_entrain_init
!  (in)      j       - meridional index upon which to work
!  (in)      ds_dsp1 - coordinate variable (sigma-2) difference across an
!                      interface divided by the difference across the interface
!                      below it (nondimensional)
!  (in)      rho_0   - layer potential densities relative to surface press (kg/m3)

  real :: g_R0                     ! g_R0 is g/Rho (m4 kg-1 s-2)
  real :: eps, tmp                 ! nondimensional temproray variables
  real :: a(szk_(g)), a_0(szk_(g)) ! nondimensional temporary variables
  real :: p_ref(szi_(g))           ! an array of tv%P_Ref pressures
  real :: Rcv(szi_(g),szk_(g))     ! coordinate density in the mixed and buffer layers (kg/m3)
  real :: I_Drho                   ! temporary variable (m3/kg)

  integer :: i, k, k3, is, ie, nz, kmb
  is = g%isc ; ie = g%iec ; nz = g%ke

  do k=2,nz-1
    if (gv%g_prime(k+1)/=0.) then
      do i=is,ie
        ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
      enddo
    else
      do i=is,ie
        ds_dsp1(i,k) = 1.
      enddo
    endif
  enddo

  if (cs%bulkmixedlayer) then
    g_r0 = gv%g_Earth/gv%Rho0
    kmb = gv%nk_rho_varies
    eps = 0.1
    do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
    do k=1,kmb
      call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p_ref,rcv(:,k),&
                             is,ie-is+1,tv%eqn_of_state)
    enddo
    do i=is,ie
      if (kb(i) <= nz-1) then
!   Set up appropriately limited ratios of the reduced gravities of the
! interfaces above and below the buffer layer and the next denser layer.
        k = kb(i)

        i_drho = g_r0 / gv%g_prime(k+1)
        ! The indexing convention for a is appropriate for the interfaces.
        do k3=1,kmb
          a(k3+1) = (gv%Rlay(k) - rcv(i,k3)) * i_drho
        enddo
        if ((present(rho_0)) .and. (a(kmb+1) < 2.0*eps*ds_dsp1(i,k))) then
!   If the buffer layer nearly matches the density of the layer below in the
! coordinate variable (sigma-2), use the sigma-0-based density ratio if it is
! greater (and stable).
          if ((rho_0(i,k) > rho_0(i,kmb)) .and. &
              (rho_0(i,k+1) > rho_0(i,k))) then
            i_drho = 1.0 / (rho_0(i,k+1)-rho_0(i,k))
            a_0(kmb+1) = min((rho_0(i,k)-rho_0(i,kmb)) * i_drho, ds_dsp1(i,k))
            if (a_0(kmb+1) > a(kmb+1)) then
              do k3=2,kmb
                a_0(k3) = a_0(kmb+1) + (rho_0(i,kmb)-rho_0(i,k3-1)) * i_drho
              enddo
              if (a(kmb+1) <= eps*ds_dsp1(i,k)) then
                do k3=2,kmb+1 ; a(k3) = a_0(k3) ; enddo
              else
! Alternative...  tmp = 0.5*(1.0 - cos(PI*(a(K2+1)/(eps*ds_dsp1(i,k)) - 1.0)) )
                tmp = a(kmb+1)/(eps*ds_dsp1(i,k)) - 1.0
                do k3=2,kmb+1 ; a(k3) = tmp*a(k3) + (1.0-tmp)*a_0(k3) ; enddo
              endif
            endif
          endif
        endif

        ds_dsp1(i,k) = max(a(kmb+1),1e-5)

        do k3=2,kmb
!           ds_dsp1(i,k3) = MAX(a(k3),1e-5)
          ! Deliberately treat convective instabilies of the upper mixed
          ! and buffer layers with respect to the deepest buffer layer as
          ! though they don't exist.  They will be eliminated by the upcoming
          ! call to the mixedlayer code anyway.
          ! The indexing convention is appropriate for the interfaces.
          ds_dsp1(i,k3) = max(a(k3),ds_dsp1(i,k))
        enddo
      endif ! (kb(i) <= nz-1)
    enddo ! I-loop.
  endif ! bulkmixedlayer

end subroutine set_density_ratios

subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp)
  type(time_type),          intent(in)    :: Time
  type(ocean_grid_type),    intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),  intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(param_file_type),    intent(in)    :: param_file !< A structure to parse for run-time
                                                  !! parameters.
  type(diag_ctrl), target,  intent(inout) :: diag !< structure used to regulate diagnostic output.
  type(set_diffusivity_cs), pointer       :: CS   !< pointer set to point to the module control
                                                  !! structure.
  type(diag_to_z_cs),       pointer       :: diag_to_Z_CSp !< pointer to the Z-diagnostics control
                                                  !! structure.
  type(int_tide_cs),        pointer       :: int_tide_CSp  !< pointer to the internal tides control
                                                  !! structure (BDM)

! Arguments:
!  (in)      Time          - current model time
!  (in)      G             - ocean grid structure
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file    - structure indicating open file to parse for params
!  (in)      diag          - structure used to regulate diagnostic output
!  (in/out)  CS            - pointer set to point to the module control structure
!  (in)      diag_to_Z_CSp - pointer to the Z-diagnostics control structure
!  (in)      int_tide_CSp  - pointer to the internal tides control structure (BDM)

  real :: decay_length, utide, zbot, hamp
  type(vardesc) :: vd
  logical :: read_tideamp, ML_use_omega
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_set_diffusivity"  ! This module's name.
  character(len=20)  :: tmpstr
  character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file
  real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data
  real :: omega_frac_dflt
  integer :: i, j, is, ie, js, je
  integer :: isd, ied, jsd, jed

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

  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  cs%diag => diag
  if (associated(int_tide_csp))  cs%int_tide_CSp  => int_tide_csp
  if (associated(diag_to_z_csp)) cs%diag_to_Z_CSp => diag_to_z_csp

  ! These default values always need to be set.
  cs%BBL_mixing_as_max = .true.
  cs%Kdml = 0.0 ; cs%cdrag = 0.003 ; cs%BBL_effic = 0.0 ;
  cs%bulkmixedlayer = (gv%nkml > 0)


  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")

  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
  cs%inputdir = slasher(cs%inputdir)
  call get_param(param_file, mdl, "FLUX_RI_MAX", cs%FluxRi_max, &
                 "The flux Richardson number where the stratification is \n"//&
                 "large enough that N2 > omega2.  The full expression for \n"//&
                 "the Flux Richardson number is usually \n"//&
                 "FLUX_RI_MAX*N2/(N2+OMEGA2).", default=0.2)
  call get_param(param_file, mdl, "OMEGA", cs%omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)

  call get_param(param_file, mdl, "ML_RADIATION", cs%ML_radiation, &
                 "If true, allow a fraction of TKE available from wind \n"//&
                 "work to penetrate below the base of the mixed layer \n"//&
                 "with a vertical decay scale determined by the minimum \n"//&
                 "of: (1) The depth of the mixed layer, (2) an Ekman \n"//&
                 "length scale.", default=.false.)
  if (cs%ML_radiation) then
    ! This give a minimum decay scale that is typically much less than Angstrom.
    cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom + gv%H_subroundoff)

    call get_param(param_file, mdl, "ML_RAD_EFOLD_COEFF", cs%ML_rad_efold_coeff, &
                 "A coefficient that is used to scale the penetration \n"//&
                 "depth for turbulence below the base of the mixed layer. \n"//&
                 "This is only used if ML_RADIATION is true.", units="nondim", &
                 default=0.2)
    call get_param(param_file, mdl, "ML_RAD_KD_MAX", cs%ML_rad_kd_max, &
                 "The maximum diapycnal diffusivity due to turbulence \n"//&
                 "radiated from the base of the mixed layer. \n"//&
                 "This is only used if ML_RADIATION is true.", units="m2 s-1", &
                 default=1.0e-3)
    call get_param(param_file, mdl, "ML_RAD_COEFF", cs%ML_rad_coeff, &
                 "The coefficient which scales MSTAR*USTAR^3 to obtain \n"//&
                 "the energy available for mixing below the base of the \n"//&
                 "mixed layer. This is only used if ML_RADIATION is true.", &
                 units="nondim", default=0.2)
    call get_param(param_file, mdl, "ML_RAD_APPLY_TKE_DECAY", cs%ML_rad_TKE_decay, &
                 "If true, apply the same exponential decay to ML_rad as \n"//&
                 "is applied to the other surface sources of TKE in the \n"//&
                 "mixed layer code. This is only used if ML_RADIATION is true.",&
                 default=.true.)
    call get_param(param_file, mdl, "MSTAR", cs%mstar, &
                 "The ratio of the friction velocity cubed to the TKE \n"//&
                 "input to the mixed layer.", "units=nondim", default=1.2)
    call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
                 "The ratio of the natural Ekman depth to the TKE decay scale.", &
                 units="nondim", default=2.5)
    call get_param(param_file, mdl, "ML_USE_OMEGA", ml_use_omega, &
                 "If true, use the absolute rotation rate instead of the \n"//&
                 "vertical component of rotation when setting the decay \n"//&
                 "scale for turbulence.", default=.false., do_not_log=.true.)
    omega_frac_dflt = 0.0
    if (ml_use_omega) then
      call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
      omega_frac_dflt = 1.0
    endif
    call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%ML_omega_frac, &
                   "When setting the decay scale for turbulence, use this \n"//&
                   "fraction of the absolute rotation rate blended with the \n"//&
                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
                   units="nondim", default=omega_frac_dflt)
  endif

  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
                 "If true, the bottom stress is calculated with a drag \n"//&
                 "law of the form c_drag*|u|*u. The velocity magnitude \n"//&
                 "may be an assumed value or it may be based on the \n"//&
                 "actual velocity in the bottommost HBBL, depending on \n"//&
                 "LINEAR_DRAG.", default=.true.)
  if  (cs%bottomdraglaw) then
    call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
                 "The drag coefficient relating the magnitude of the \n"//&
                 "velocity field to the bottom stress. CDRAG is only used \n"//&
                 "if BOTTOMDRAGLAW is true.", units="nondim", default=0.003)
    call get_param(param_file, mdl, "BBL_EFFIC", cs%BBL_effic, &
                 "The efficiency with which the energy extracted by \n"//&
                 "bottom drag drives BBL diffusion.  This is only \n"//&
                 "used if BOTTOMDRAGLAW is true.", units="nondim", default=0.20)
    call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, &
                 "The maximum decay scale for the BBL diffusion, or 0 \n"//&
                 "to allow the mixing to penetrate as far as \n"//&
                 "stratification and rotation permit.  The default is 0. \n"//&
                 "This is only used if BOTTOMDRAGLAW is true.", units="m", &
                 default=0.0)

    cs%IMax_decay = 1.0/200.0
    if (decay_length > 0.0) cs%IMax_decay = 1.0/decay_length
    call get_param(param_file, mdl, "BBL_MIXING_AS_MAX", cs%BBL_mixing_as_max, &
                 "If true, take the maximum of the diffusivity from the \n"//&
                 "BBL mixing and the other diffusivities. Otherwise, \n"//&
                 "diffusiviy from the BBL_mixing is simply added.", &
                 default=.true.)
    call get_param(param_file, mdl, "USE_LOTW_BBL_DIFFUSIVITY", cs%use_LOTW_BBL_diffusivity, &
                 "If true, uses a simple, imprecise but non-coordinate dependent, model\n"//&
                 "of BBL mixing diffusivity based on Law of the Wall. Otherwise, uses\n"//&
                 "the original BBL scheme.", default=.false.)
    if (cs%use_LOTW_BBL_diffusivity) then
      call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", cs%LOTW_BBL_use_omega, &
                 "If true, use the maximum of Omega and N for the TKE to diffusion\n"//&
                 "calculation. Otherwise, N is N.", default=.true.)
    endif
  else
    cs%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL
  endif
  cs%id_Kd_BBL = register_diag_field('ocean_model','Kd_BBL',diag%axesTi,time, &
       'Bottom Boundary Layer Diffusivity', 'meter2 sec-1')
  call get_param(param_file, mdl, "SIMPLE_TKE_TO_KD", cs%simple_TKE_to_Kd, &
                 "If true, uses a simple estimate of Kd/TKE that will\n"//&
                 "work for arbitrary vertical coordinates. If false,\n"//&
                 "calculates Kd/TKE and bounds based on exact energetics/n"//&
                 "for an isopycnal layer-formulation.", &
                 default=.false.)

  call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", &
                                cs%Bryan_Lewis_diffusivity, &
                 "If true, use a Bryan & Lewis (JGR 1979) like tanh \n"//&
                 "profile of background diapycnal diffusivity with depth.", &
                 default=.false.)
  if (cs%Bryan_Lewis_diffusivity) then
    call get_param(param_file, mdl, "KD_BRYAN_LEWIS_DEEP", &
                                  cs%Kd_Bryan_Lewis_deep, &
                 "The abyssal value of a Bryan-Lewis diffusivity profile. \n"//&
                 "KD_BRYAN_LEWIS_DEEP is only used if \n"//&
                 "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "KD_BRYAN_LEWIS_SURFACE", &
                                  cs%Kd_Bryan_Lewis_surface, &
                 "The surface value of a Bryan-Lewis diffusivity profile. \n"//&
                 "KD_BRYAN_LEWIS_SURFACE is only used if \n"//&
                 "BRYAN_LEWIS_DIFFUSIVITY is true.", units="m2 s-1", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "BRYAN_LEWIS_DEPTH_CENT", &
                                  cs%Bryan_Lewis_depth_cent, &
                 "The depth about which the transition in the Bryan-Lewis \n"//&
                 "profile is centered. BRYAN_LEWIS_DEPTH_CENT is only \n"//&
                 "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "BRYAN_LEWIS_WIDTH_TRANS", &
                                  cs%Bryan_Lewis_width_trans, &
                 "The width of the transition in the Bryan-Lewis \n"//&
                 "profile. BRYAN_LEWIS_WIDTH_TRANS is only \n"//&
                 "used if BRYAN_LEWIS_DIFFUSIVITY is true.", units="m", &
                 fail_if_missing=.true.)
  endif

  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", &
                                cs%Henyey_IGW_background, &
                 "If true, use a latitude-dependent scaling for the near \n"//&
                 "surface background diffusivity, as described in \n"//&
                 "Harrison & Hallberg, JPO 2008.", default=.false.)
  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", &
                                cs%Henyey_IGW_background_new, &
                 "If true, use a better latitude-dependent scaling for the\n"//&
                 "background diffusivity, as described in \n"//&
                 "Harrison & Hallberg, JPO 2008.", default=.false.)
  if (cs%Henyey_IGW_background .and. cs%Henyey_IGW_background_new) call mom_error(fatal, &
                 "set_diffusivity_init: HENYEY_IGW_BACKGROUND and HENYEY_IGW_BACKGROUND_NEW "// &
                 "are mutually exclusive. Set only one or none.")
  if (cs%Henyey_IGW_background) &
    call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", cs%N0_2Omega, &
                  "The ratio of the typical Buoyancy frequency to twice \n"//&
                  "the Earth's rotation period, used with the Henyey \n"//&
                  "scaling from the mixing.", units="nondim", default=20.0)
  call get_param(param_file, mdl, "N2_FLOOR_IOMEGA2", cs%N2_FLOOR_IOMEGA2, &
                  "The floor applied to N2(k) scaled by Omega^2:\n"//&
                  "\tIf =0., N2(k) is simply positive definite.\n"//&
                  "\tIf =1., N2(k) > Omega^2 everywhere.", units="nondim", &
                  default=1.0)

  call get_param(param_file, mdl, "KD_TANH_LAT_FN", &
                                  cs%Kd_tanh_lat_fn, &
                 "If true, use a tanh dependence of Kd_sfc on latitude, \n"//&
                 "like CM2.1/CM2M.  There is no physical justification \n"//&
                 "for this form, and it can not be used with \n"//&
                 "HENYEY_IGW_BACKGROUND.", default=.false.)
  if (cs%Kd_tanh_lat_fn) &
    call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", &
                                  cs%Kd_tanh_lat_scale, &
                 "A nondimensional scaling for the range ofdiffusivities \n"//&
                 "with KD_TANH_LAT_FN. Valid values are in the range of \n"//&
                 "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0)

  call get_param(param_file, mdl, "KV", cs%Kv, &
                 "The background kinematic viscosity in the interior. \n"//&
                 "The molecular value, ~1e-6 m2 s-1, may be used.", &
                 units="m2 s-1", fail_if_missing=.true.)

  call get_param(param_file, mdl, "KD", cs%Kd, &
                 "The background diapycnal diffusivity of density in the \n"//&
                 "interior. Zero or the molecular value, ~1e-7 m2 s-1, \n"//&
                 "may be used.", units="m2 s-1", fail_if_missing=.true.)
  call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
                 "The minimum diapycnal diffusivity.", &
                 units="m2 s-1", default=0.01*cs%Kd)
  call get_param(param_file, mdl, "KD_MAX", cs%Kd_max, &
                 "The maximum permitted increment for the diapycnal \n"//&
                 "diffusivity from TKE-based parameterizations, or a \n"//&
                 "negative value for no limit.", units="m2 s-1", default=-1.0)
  if (cs%simple_TKE_to_Kd .and. cs%Kd_max<=0.) call mom_error(fatal, &
         "set_diffusivity_init: To use SIMPLE_TKE_TO_KD, KD_MAX must be set to >0.")
  call get_param(param_file, mdl, "KD_ADD", cs%Kd_add, &
                 "A uniform diapycnal diffusivity that is added \n"//&
                 "everywhere without any filtering or scaling.", &
                 units="m2 s-1", default=0.0)
  if (cs%use_LOTW_BBL_diffusivity .and. cs%Kd_max<=0.) call mom_error(fatal, &
                 "set_diffusivity_init: KD_MAX must be set (positive) when "// &
                 "USE_LOTW_BBL_DIFFUSIVITY=True.")

  if (cs%bulkmixedlayer) then
    ! Check that Kdml is not set when using bulk mixed layer
    call get_param(param_file, mdl, "KDML", cs%Kdml, default=-1.)
    if (cs%Kdml>0.) call mom_error(fatal, &
                 "set_diffusivity_init: KDML cannot be set when using"// &
                 "bulk mixed layer.")
    cs%Kdml = cs%Kd ! This is not used with a bulk mixed layer, but also
                    ! cannot be a NaN.
  else
    call get_param(param_file, mdl, "KDML", cs%Kdml, &
                 "If BULKMIXEDLAYER is false, KDML is the elevated \n"//&
                 "diapycnal diffusivity in the topmost HMIX of fluid. \n"//&
                 "KDML is only used if BULKMIXEDLAYER is false.", &
                 units="m2 s-1", default=cs%Kd)
    call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
                 "The prescribed depth over which the near-surface \n"//&
                 "viscosity and diffusivity are elevated when the bulk \n"//&
                 "mixed layer is not used.", units="m", fail_if_missing=.true.)
  endif
  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", default=.false.)

  call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", cs%Int_tide_dissipation, &
                 "If true, use an internal tidal dissipation scheme to \n"//&
                 "drive diapycnal mixing, along the lines of St. Laurent \n"//&
                 "et al. (2002) and Simmons et al. (2004).", default=.false.)
  if (cs%Int_tide_dissipation) then
    call get_param(param_file, mdl, "INT_TIDE_PROFILE", tmpstr, &
                 "INT_TIDE_PROFILE selects the vertical profile of energy \n"//&
                 "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//&
                 "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
                 "\t                decay profile.\n"//&
                 "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//&
                 "\t                decay profile.", &
                 default=stlaurent_profile_string)
    tmpstr = uppercase(tmpstr)
    select case (tmpstr)
      case (stlaurent_profile_string) ; cs%int_tide_profile = stlaurent_02
      case (polzin_profile_string) ; cs%int_tide_profile = polzin_09
      case default
        call mom_error(fatal, "set_diffusivity_init: Unrecognized setting "// &
            "#define INT_TIDE_PROFILE "//trim(tmpstr)//" found in input file.")
    end select
  endif

  call get_param(param_file, mdl, "LEE_WAVE_DISSIPATION", cs%Lee_wave_dissipation, &
                 "If true, use an lee wave driven dissipation scheme to \n"//&
                 "drive diapycnal mixing, along the lines of Nikurashin \n"//&
                 "(2010) and using the St. Laurent et al. (2002) \n"//&
                 "and Simmons et al. (2004) vertical profile", default=.false.)
  if (cs%lee_wave_dissipation) then
    call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, &
                 "LEE_WAVE_PROFILE selects the vertical profile of energy \n"//&
                 "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//&
                 "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
                 "\t                decay profile.\n"//&
                 "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//&
                 "\t                decay profile.", &
                 default=stlaurent_profile_string)
    tmpstr = uppercase(tmpstr)
    select case (tmpstr)
      case (stlaurent_profile_string) ; cs%lee_wave_profile = stlaurent_02
      case (polzin_profile_string) ; cs%lee_wave_profile = polzin_09
      case default
        call mom_error(fatal, "set_diffusivity_init: Unrecognized setting "// &
            "#define LEE_WAVE_PROFILE "//trim(tmpstr)//" found in input file.")
    end select
  endif

  call get_param(param_file, mdl, "INT_TIDE_LOWMODE_DISSIPATION", cs%Lowmode_itidal_dissipation, &
                 "If true, consider mixing due to breaking low modes that \n"//&
                 "have been remotely generated; as with itidal drag on the \n"//&
                 "barotropic tide, use an internal tidal dissipation scheme to \n"//&
                 "drive diapycnal mixing, along the lines of St. Laurent \n"//&
                 "et al. (2002) and Simmons et al. (2004).", default=.false.)

  if ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09)) .or. &
      (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09))) then
    call get_param(param_file, mdl, "NU_POLZIN", cs%Nu_Polzin, &
                 "When the Polzin decay profile is used, this is a \n"//&
                 "non-dimensional constant in the expression for the \n"//&
                 "vertical scale of decay for the tidal energy dissipation.", &
                 units="nondim", default=0.0697)
    call get_param(param_file, mdl, "NBOTREF_POLZIN", cs%Nbotref_Polzin, &
                 "When the Polzin decay profile is used, this is the \n"//&
                 "Rreference value of the buoyancy frequency at the ocean \n"//&
                 "bottom in the Polzin formulation for the vertical \n"//&
                 "scale of decay for the tidal energy dissipation.", &
                 units="s-1", default=9.61e-4)
    call get_param(param_file, mdl, "POLZIN_DECAY_SCALE_FACTOR", &
                 cs%Polzin_decay_scale_factor, &
                 "When the Polzin decay profile is used, this is a \n"//&
                 "scale factor for the vertical scale of decay of the tidal \n"//&
                 "energy dissipation.", default=1.0, units="nondim")
    call get_param(param_file, mdl, "POLZIN_SCALE_MAX_FACTOR", &
                 cs%Polzin_decay_scale_max_factor, &
                 "When the Polzin decay profile is used, this is a factor \n"//&
                 "to limit the vertical scale of decay of the tidal \n"//&
                 "energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR \n"//&
                 "times the depth of the ocean.", units="nondim", default=1.0)
    call get_param(param_file, mdl, "POLZIN_MIN_DECAY_SCALE", cs%Polzin_min_decay_scale, &
                 "When the Polzin decay profile is used, this is the \n"//&
                 "minimum vertical decay scale for the vertical profile\n"//&
                 "of internal tide dissipation with the Polzin (2009) formulation", &
                 units="m", default=0.0)
  endif
  call get_param(param_file, mdl, "USER_CHANGE_DIFFUSIVITY", cs%user_change_diff, &
                 "If true, call user-defined code to change the diffusivity.", &
                 default=.false.)

  call get_param(param_file, mdl, "DISSIPATION_MIN", cs%dissip_min, &
                 "The minimum dissipation by which to determine a lower \n"//&
                 "bound of Kd (a floor).", units="W m-3", default=0.0)
  call get_param(param_file, mdl, "DISSIPATION_N0", cs%dissip_N0, &
                 "The intercept when N=0 of the N-dependent expression \n"//&
                 "used to set a minimum dissipation by which to determine \n"//&
                 "a lower bound of Kd (a floor): A in eps_min = A + B*N.", &
                 units="W m-3", default=0.0)
  call get_param(param_file, mdl, "DISSIPATION_N1", cs%dissip_N1, &
                 "The coefficient multiplying N, following Gargett, used to \n"//&
                 "set a minimum dissipation by which to determine a lower \n"//&
                 "bound of Kd (a floor): B in eps_min = A + B*N", &
                 units="J m-3", default=0.0)
  call get_param(param_file, mdl, "DISSIPATION_KD_MIN", cs%dissip_Kd_min, &
                 "The minimum vertical diffusivity applied as a floor.", &
                 units="m2 s-1", default=0.0)

  cs%limit_dissipation = (cs%dissip_min>0.) .or. (cs%dissip_N1>0.) .or. &
                         (cs%dissip_N0>0.) .or. (cs%dissip_Kd_min>0.)
  cs%dissip_N2 = 0.0
  if (cs%FluxRi_max > 0.0) &
    cs%dissip_N2 = cs%dissip_Kd_min * gv%Rho0 / cs%FluxRi_max

  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) then
    call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", cs%Int_tide_decay_scale, &
                 "The decay scale away from the bottom for tidal TKE with \n"//&
                 "the new coding when INT_TIDE_DISSIPATION is used.", &
                 units="m", default=0.0)
    call get_param(param_file, mdl, "MU_ITIDES", cs%Mu_itides, &
                 "A dimensionless turbulent mixing efficiency used with \n"//&
                 "INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2)
    call get_param(param_file, mdl, "GAMMA_ITIDES", cs%Gamma_itides, &
                 "The fraction of the internal tidal energy that is \n"//&
                 "dissipated locally with INT_TIDE_DISSIPATION.  \n"//&
                 "THIS NAME COULD BE BETTER.", &
                 units="nondim", default=0.3333)
    call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", cs%min_zbot_itides, &
                 "Turn off internal tidal dissipation when the total \n"//&
                 "ocean depth is less than this value.", units="m", default=0.0)

    call safe_alloc_ptr(cs%Nb,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%h2,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%TKE_itidal,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%mask_itidal,isd,ied,jsd,jed) ; cs%mask_itidal(:,:) = 1.0

    call get_param(param_file, mdl, "KAPPA_ITIDES", cs%kappa_itides, &
                 "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//&
                 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
                 units="m-1", default=8.e-4*atan(1.0))

    call get_param(param_file, mdl, "UTIDE", cs%utide, &
                 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
                 units="m s-1", default=0.0)
    call safe_alloc_ptr(cs%tideamp,is,ie,js,je) ; cs%tideamp(:,:) = cs%utide

    call get_param(param_file, mdl, "KAPPA_H2_FACTOR", cs%kappa_h2_factor, &
                 "A scaling factor for the roughness amplitude with n"//&
                 "INT_TIDE_DISSIPATION.",  units="nondim", default=1.0)
    call get_param(param_file, mdl, "TKE_ITIDE_MAX", cs%TKE_itide_max, &
                 "The maximum internal tide energy source availble to mix \n"//&
                 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
                 units="W m-2",  default=1.0e3)

    call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
                 "If true, read a file (given by TIDEAMP_FILE) containing \n"//&
                 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
    if (read_tideamp) then
      call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
                 "The path to the file containing the spatially varying \n"//&
                 "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
      filename = trim(cs%inputdir) // trim(tideamp_file)
      call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
      call read_data(filename, 'tideamp', cs%tideamp, &
                     domain=g%domain%mpp_domain, timelevel=1)
    endif

    call get_param(param_file, mdl, "H2_FILE", h2_file, &
                 "The path to the file containing the sub-grid-scale \n"//&
                 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
                 fail_if_missing=.true.)
    filename = trim(cs%inputdir) // trim(h2_file)
    call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
    call read_data(filename, 'h2', cs%h2, domain=g%domain%mpp_domain, &
                   timelevel=1)

    do j=js,je ; do i=is,ie
      if (g%bathyT(i,j) < cs%min_zbot_itides) cs%mask_itidal(i,j) = 0.0
      cs%tideamp(i,j) = cs%tideamp(i,j) * cs%mask_itidal(i,j) * g%mask2dT(i,j)

      ! Restrict rms topo to 10 percent of column depth.
      zbot = g%bathyT(i,j)
      hamp = sqrt(cs%h2(i,j))
      hamp = min(0.1*zbot,hamp)
      cs%h2(i,j) = hamp*hamp

      utide = cs%tideamp(i,j)
      ! Compute the fixed part of internal tidal forcing; units are [kg s-2] here.
      cs%TKE_itidal(i,j) = 0.5*cs%kappa_h2_factor*gv%Rho0*&
           cs%kappa_itides*cs%h2(i,j)*utide*utide
    enddo; enddo

  endif

  cs%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, time, &
      'Diapycnal diffusivity of layers (as set)', 'meter2 second-1')

  if (cs%Lee_wave_dissipation) then

    call get_param(param_file, mdl, "NIKURASHIN_TKE_INPUT_FILE",niku_tke_input_file, &
                 "The path to the file containing the TKE input from lee \n"//&
                 "wave driven mixing. Used with LEE_WAVE_DISSIPATION.", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "NIKURASHIN_SCALE",niku_scale, &
                 "A non-dimensional factor by which to scale the lee-wave \n"//&
                 "driven TKE input. Used with LEE_WAVE_DISSIPATION.", &
                 units="nondim", default=1.0)

    filename = trim(cs%inputdir) // trim(niku_tke_input_file)
    call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", &
                   filename)
    call safe_alloc_ptr(cs%TKE_Niku,is,ie,js,je); cs%TKE_Niku(:,:) = 0.0
    call read_data(filename, 'TKE_input', cs%TKE_Niku, &
                   domain=g%domain%mpp_domain, timelevel=1 ) ! ??? timelevel -aja
    cs%TKE_Niku(:,:) = niku_scale * cs%TKE_Niku(:,:)

    call get_param(param_file, mdl, "GAMMA_NIKURASHIN",cs%Gamma_lee, &
                 "The fraction of the lee wave energy that is dissipated \n"//&
                 "locally with LEE_WAVE_DISSIPATION.", units="nondim", &
                 default=0.3333)
    call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",cs%Decay_scale_factor_lee, &
                 "Scaling for the vertical decay scaleof the local \n"//&
                 "dissipation of lee waves dissipation.", units="nondim", &
                 default=1.0)

    cs%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,time, &
        'Lee wave Driven Turbulent Kinetic Energy', 'Watt meter-2')
    cs%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,time, &
         'Lee Wave Driven Diffusivity', 'meter2 sec-1')
  else
    cs%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False
  endif

  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. &
      cs%Lowmode_itidal_dissipation) then

    cs%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,time, &
        'Internal Tide Driven Turbulent Kinetic Energy', 'Watt meter-2')
    cs%id_maxTKE = register_diag_field('ocean_model','maxTKE',diag%axesTL,time, &
           'Maximum layer TKE', 'meter3 second-3')
    cs%id_TKE_to_Kd = register_diag_field('ocean_model','TKE_to_Kd',diag%axesTL,time, &
           'Convert TKE to Kd', 'second2 meter')

    cs%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,time, &
         'Bottom Buoyancy Frequency', 'sec-1')

    cs%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,time, &
         'Internal Tide Driven Diffusivity', 'meter2 sec-1')

    cs%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,time, &
         'Internal Tide Driven Diffusivity (from propagating low modes)', 'meter2 sec-1')

    cs%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,time, &
        'Vertical flux of tidal turbulent dissipation', 'meter3 sec-3')

    cs%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,time, &
         'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'meter3 sec-3')

    cs%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,time, &
         'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'meter')

    cs%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model','Polzin_decay_scale_scaled',diag%axesT1,time, &
         'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, scaled by N2_bot/N2_meanz', 'meter')

    cs%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,time, &
         'Bottom Buoyancy frequency squared', 's-2')

    cs%id_N2_meanz = register_diag_field('ocean_model','N2_meanz',diag%axesT1,time, &
         'Buoyancy frequency squared averaged over the water column', 's-2')

    cs%id_Kd_Work = register_diag_field('ocean_model','Kd_Work',diag%axesTL,time, &
         'Work done by Diapycnal Mixing', 'Watts m-2')

    cs%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,time, &
         'Work done by Internal Tide Diapycnal Mixing', 'Watts m-2')

    cs%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,time, &
         'Work done by Nikurashin Lee Wave Drag Scheme', 'Watts m-2')

    cs%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,time, &
         'Work done by Internal Tide Diapycnal Mixing (low modes)', 'Watts m-2')

    cs%id_N2 = register_diag_field('ocean_model','N2',diag%axesTi,time,            &
         'Buoyancy frequency squared', 'sec-2', cmor_field_name='obvfsq',          &
          cmor_units='s-2', cmor_long_name='Square of seawater buoyancy frequency',&
          cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water')

    if (cs%user_change_diff) &
      cs%id_Kd_user = register_diag_field('ocean_model','Kd_user',diag%axesTi,time, &
           'User-specified Extra Diffusivity', 'meter2 sec-1')

    if (associated(diag_to_z_csp)) then
      vd = var_desc("N2", "second-2",&
                    "Buoyancy frequency, interpolated to z", z_grid='z')
      cs%id_N2_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      vd = var_desc("Kd_itides","meter2 second-1", &
                    "Internal Tide Driven Diffusivity, interpolated to z", z_grid='z')
      cs%id_Kd_itidal_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      if (cs%Lee_wave_dissipation) then
         vd = var_desc("Kd_Nikurashin", "meter2 second-1", &
                       "Lee Wave Driven Diffusivity, interpolated to z", z_grid='z')
         cs%id_Kd_Niku_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      endif
      if (cs%Lowmode_itidal_dissipation) then
        vd = var_desc("Kd_lowmode","meter2 second-1", &
                  "Internal Tide Driven Diffusivity (from low modes), interpolated to z",&
                  z_grid='z')
        cs%id_Kd_lowmode_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      endif
      if (cs%user_change_diff) &
        cs%id_Kd_user_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    endif
  endif

  call get_param(param_file, mdl, "DOUBLE_DIFFUSION", cs%double_diffusion, &
                 "If true, increase diffusivitives for temperature or salt \n"//&
                 "based on double-diffusive paramaterization from MOM4/KPP.", &
                 default=.false.)
  if (cs%double_diffusion) then
    call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", cs%Max_Rrho_salt_fingers, &
                 "Maximum density ratio for salt fingering regime.", &
                 default=2.55, units="nondim")
    call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", cs%Max_salt_diff_salt_fingers, &
                 "Maximum salt diffusivity for salt fingering regime.", &
                 default=1.e-4, units="m2 s-1")
    call get_param(param_file, mdl, "KV_MOLECULAR", cs%Kv_molecular, &
                 "Molecular viscosity for calculation of fluxes under \n"//&
                 "double-diffusive convection.", default=1.5e-6, units="m2 s-1")
    ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults.

    cs%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,time, &
         'Double-diffusive diffusivity for temperature', 'meter2 sec-1')

    cs%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,time, &
         'Double-diffusive diffusivity for salinity', 'meter2 sec-1')

    if (associated(diag_to_z_csp)) then
      vd = var_desc("KT_extra", "meter2 second-1", &
                    "Double-Diffusive Temperature Diffusivity, interpolated to z", &
                    z_grid='z')
      cs%id_KT_extra_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      vd = var_desc("KS_extra", "meter2 second-1", &
                    "Double-Diffusive Salinity Diffusivity, interpolated to z",&
                    z_grid='z')
      cs%id_KS_extra_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
      vd = var_desc("Kd_BBL", "meter2 second-1", &
                    "Bottom Boundary Layer Diffusivity", z_grid='z')
      cs%id_Kd_BBL_z = register_zint_diag(vd, cs%diag_to_Z_CSp, time)
    endif
  endif

  if (cs%Int_tide_dissipation .and. cs%Bryan_Lewis_diffusivity) &
    call mom_error(fatal,"MOM_Set_Diffusivity: "// &
         "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.")
  if (cs%Henyey_IGW_background .and. cs%Kd_tanh_lat_fn) call mom_error(fatal, &
    "Set_diffusivity: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.")

  if (cs%user_change_diff) then
    call user_change_diff_init(time, g, param_file, diag, cs%user_change_diff_CSp)
  endif

  cs%useKappaShear = kappa_shear_init(time, g, gv, param_file, cs%diag, cs%kappaShear_CSp)
  if (cs%useKappaShear) &
    id_clock_kappashear = cpu_clock_id('(Ocean kappa_shear)', grain=clock_module)
  cs%useCVMix = cvmix_shear_init(time, g, gv, param_file, cs%diag, cs%CVMix_shear_CSp)



end subroutine set_diffusivity_init

subroutine set_diffusivity_end(CS)
  type(set_diffusivity_cs), pointer :: CS

  if (cs%user_change_diff) &
    call user_change_diff_end(cs%user_change_diff_CSp)

  if (associated(cs)) deallocate(cs)

end subroutine set_diffusivity_end

end module mom_set_diffusivity