MOM_tracer_hor_diff.F90

Go to the documentation of this file.

!> Main routine for lateral (along surface or neutral) diffusion of tracers
module mom_tracer_hor_diff

! This file is part of MOM6. See LICENSE.md for the license.

use mom_cpu_clock,             only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,             only : clock_module, clock_routine
use mom_diag_mediator,         only : post_data, diag_ctrl
use mom_diag_mediator,         only : register_diag_field, safe_alloc_ptr, time_type
use mom_domains,               only : sum_across_pes, max_across_pes
use mom_domains,               only : create_group_pass, do_group_pass, group_pass_type
use mom_domains,               only : pass_vector
use mom_debugging,             only : hchksum, uvchksum
use mom_eos,                   only : calculate_density
use mom_error_handler,         only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_error_handler,         only : mom_set_verbosity, calltree_showquery
use mom_error_handler,         only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,           only : get_param, log_version, param_file_type
use mom_grid,                  only : ocean_grid_type
use mom_lateral_mixing_coeffs, only : varmix_cs
use mom_meke_types,            only : meke_type
use mom_neutral_diffusion,     only : neutral_diffusion_init, neutral_diffusion_end
use mom_neutral_diffusion,     only : neutral_diffusion_cs
use mom_neutral_diffusion,     only : neutral_diffusion_calc_coeffs, neutral_diffusion
use mom_tracer_registry,       only : tracer_registry_type, tracer_type, mom_tracer_chksum
use mom_variables,             only : thermo_var_ptrs
use mom_verticalgrid,          only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public tracer_hordiff, tracer_hor_diff_init, tracer_hor_diff_end

type, public :: tracer_hor_diff_cs ; private
  real    :: dt             ! The baroclinic dynamics time step, in s.
  real    :: khtr           ! The along-isopycnal tracer diffusivity in m2/s.
  real    :: khtr_slope_cff ! The non-dimensional coefficient in KhTr formula
  real    :: khtr_min       ! Minimum along-isopycnal tracer diffusivity in m2/s.
  real    :: khtr_max       ! Maximum along-isopycnal tracer diffusivity in m2/s.
  real    :: khtr_passivity_coeff ! Passivity coefficient that scales Rd/dx (default = 0)
                                  ! where passivity is the ratio between along-isopycnal
                                  ! tracer mixing and thickness mixing
  real    :: khtr_passivity_min   ! Passivity minimum (default = 1/2)
  real    :: ml_khtr_scale  ! With Diffuse_ML_interior, the ratio of the truly
                            ! horizontal diffusivity in the mixed layer to the
                            ! epipycnal diffusivity.  Nondim.
  logical :: diffuse_ml_interior  ! If true, diffuse along isopycnals between
                                  ! the mixed layer and the interior.
  logical :: check_diffusive_cfl  ! If true, automatically iterate the diffusion
                                  ! to ensure that the diffusive equivalent of the CFL
                                  ! limit is not violated.
  logical :: use_neutral_diffusion ! If true, use the neutral_diffusion module from within
                                   ! tracer_hor_diff.

  type(neutral_diffusion_cs), pointer :: neutral_diffusion_csp => null() ! Control structure for neutral diffusion.
  type(diag_ctrl), pointer :: diag ! structure to regulate timing of diagnostic output.
  logical :: debug                 ! If true, write verbose checksums for debugging purposes.
  logical :: show_call_tree        ! Display the call tree while running. Set by VERBOSITY level.
  logical :: first_call = .true.
  integer :: id_khtr_u  = -1
  integer :: id_khtr_v  = -1
  integer :: id_khtr_h  = -1
  integer :: id_cfl     = -1
  integer :: id_khdt_x  = -1
  integer :: id_khdt_y  = -1

  type(group_pass_type) :: pass_t !For group halo pass, used in both
                                  !tracer_hordiff and tracer_epipycnal_ML_diff
end type tracer_hor_diff_cs

type p2d
  real, dimension(:,:), pointer :: p => null()
end type p2d
type p2di
  integer, dimension(:,:), pointer :: p => null()
end type p2di

integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync

contains

!> Compute along-coordinate diffusion of all tracers
!! using the diffusivity in CS%KhTr, or using space-dependent diffusivity.
!! Multiple iterations are used (if necessary) so that there is no limit
!! on the acceptable time increment.
subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y)
  type(ocean_grid_type),                 intent(inout) :: G       !< Grid type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h       !< Layer thickness (m or kg m-2)
  real,                                  intent(in)    :: dt      !< time step (seconds)
  type(meke_type),                       pointer       :: MEKE    !< MEKE type
  type(varmix_cs),                       pointer       :: VarMix  !< Variable mixing type
  type(verticalgrid_type),               intent(in)    :: GV      !< ocean vertical grid structure
  type(tracer_hor_diff_cs),              pointer       :: CS      !< module control structure
  type(tracer_registry_type),            intent(inout) :: Reg     !< registered tracers
  type(thermo_var_ptrs),                 intent(in)    :: tv      !< A structure containing pointers to any available
                                                                  !! thermodynamic fields, including potential temp and
                                                                  !! salinity or mixed layer density. Absent fields have
                                                                  !! NULL ptrs, and these may (probably will) point to
                                                                  !! some of the same arrays as Tr does.  tv is required
                                                                  !! for epipycnal mixing between mixed layer and the interior.
  ! Optional inputs for offline tracer transport
  logical,                           optional             :: do_online_flag
  real, dimension(SZIB_(G),SZJ_(G)), optional, intent(in) :: read_khdt_x
  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: read_khdt_y


  real, dimension(SZI_(G),SZJ_(G)) :: &
    Ihdxdy, &     ! The inverse of the volume or mass of fluid in a layer in a
                  ! grid cell, in m-3 or kg-1.
    kh_h, &       ! The tracer diffusivity averaged to tracer points, in m2 s-1.
    cfl, &        ! A diffusive CFL number for each cell, nondim.
    dtr           ! The change in a tracer's concentration, in units of
                  ! concentration.

  real, dimension(SZIB_(G),SZJ_(G)) :: &
    khdt_x, &     ! The value of Khtr*dt times the open face width divided by
                  ! the distance between adjacent tracer points, in m2.
    coef_x, &     ! The coefficients relating zonal tracer differences
                  ! to time-integrated fluxes, in m3 or kg.
    kh_u          ! Tracer mixing coefficient at u-points, in m2 s-1.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    khdt_y, &     ! The value of Khtr*dt times the open face width divided by
                  ! the distance between adjacent tracer points, in m2.
    coef_y, &     ! The coefficients relating meridional tracer differences
                  ! to time-integrated fluxes, in m3 or kg.
    kh_v          ! Tracer mixing coefficient at u-points, in m2 s-1.

  real :: max_CFL ! The global maximum of the diffusive CFL number.
  logical :: use_VarMix, Resoln_scaled, do_online, use_Eady
  integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts
  real :: I_numitts  ! The inverse of the number of iterations, num_itts.
  real :: scale      ! The fraction of khdt_x or khdt_y that is applied in this
                     ! layer for this iteration, nondim.
  real :: Idt        ! The inverse of the time step, in s-1.
  real :: h_neglect  ! A thickness that is so small it is usually lost
                     ! in roundoff and can be neglected, in m.
  real :: Kh_loc     ! The local value of Kh, in m2 s-1.
  real :: Res_Fn     ! The local value of the resolution function, nondim.
  real :: Rd_dx      ! The local value of deformation radius over grid-spacing, nondim.
  real :: normalize  ! normalization used for diagnostic Kh_h; diffusivity averaged to h-points.

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  do_online = .true.
  if (present(do_online_flag)) do_online = do_online_flag

  if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
       "register_tracer must be called before tracer_hordiff.")
  if (loc(reg)==0) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
       "register_tracer must be called before tracer_hordiff.")
  if ((reg%ntr==0) .or. ((cs%KhTr <= 0.0) .and. .not.associated(varmix)) ) return

  if (cs%show_call_tree) call calltree_enter("tracer_hordiff(), MOM_tracer_hor_diff.F90")

  call cpu_clock_begin(id_clock_diffuse)

  ntr = reg%ntr
  idt = 1.0/dt
  h_neglect = gv%H_subroundoff

  if (cs%Diffuse_ML_interior .and. cs%first_call) then ; if (is_root_pe()) then
    do m=1,ntr ; if (associated(reg%Tr(m)%df_x) .or. associated(reg%Tr(m)%df_y)) &
      call mom_error(warning, "tracer_hordiff: Tracer "//trim(reg%Tr(m)%name)// &
          " has associated 3-d diffusive flux diagnostics.  These are not "//&
          "valid when DIFFUSE_ML_TO_INTERIOR is defined. Use 2-d tracer "//&
          "diffusion diagnostics instead to get accurate total fluxes.")
    enddo
  endif ; endif
  cs%first_call = .false.

  if (cs%debug) call mom_tracer_chksum("Before tracer diffusion ", reg%Tr, ntr, g)

  use_varmix = .false. ; resoln_scaled = .false. ; use_eady = .false.
  if (Associated(varmix)) then
    use_varmix = varmix%use_variable_mixing
    resoln_scaled = varmix%Resoln_scaled_KhTr
    use_eady = cs%KhTr_Slope_Cff > 0.
  endif

  call cpu_clock_begin(id_clock_pass)
  do m=1,ntr
    call create_group_pass(cs%pass_t, reg%Tr(m)%t(:,:,:), g%Domain)
  enddo
  call cpu_clock_end(id_clock_pass)

  if (cs%show_call_tree) call calltree_waypoint("Calculating diffusivity (tracer_hordiff)")

  if (do_online) then
      if (use_varmix) then
    !$OMP parallel default(none) shared(is,ie,js,je,CS,VarMix,MEKE,Resoln_scaled, &
    !$OMP                               Kh_u,Kh_v,khdt_x,dt,G,khdt_y,use_Eady)    &
    !$OMP                       private(Kh_loc,Rd_dx)
    !$OMP do
        do j=js,je ; do i=is-1,ie
          kh_loc = cs%KhTr
          if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2u(i,j)*varmix%SN_u(i,j)
          if (associated(meke%Kh)) &
            kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
          if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
          if (resoln_scaled) &
            kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
          kh_u(i,j) = max(kh_loc, cs%KhTr_min)
          if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
            rd_dx=0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points
            kh_loc=kh_u(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
            if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
            kh_u(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
          endif
        enddo ; enddo
    !$OMP do
        do j=js-1,je ;  do i=is,ie
          kh_loc = cs%KhTr
          if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
          if (associated(meke%Kh)) &
            kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
          if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
          if (resoln_scaled) &
            kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
          kh_v(i,j) = max(kh_loc, cs%KhTr_min)
          if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
            rd_dx=0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points
            kh_loc=kh_v(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
            if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
            kh_v(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
          endif
        enddo ; enddo

    !$OMP do
        do j=js,je ; do i=is-1,ie
          khdt_x(i,j) = dt*(kh_u(i,j)*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
        enddo ; enddo
    !$OMP do
        do j=js-1,je ; do i=is,ie
          khdt_y(i,j) = dt*(kh_v(i,j)*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
        enddo ; enddo
    !$OMP end parallel
      elseif (resoln_scaled) then
    !$OMP parallel default(none) shared(is,ie,js,je,VarMix,Kh_u,Kh_v,khdt_x,khdt_y,CS,dt,G) &
    !$OMP                       private(Res_fn)
    !$OMP do
        do j=js,je ; do i=is-1,ie
          res_fn = 0.5 * (varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
          kh_u(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
          khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j))) * res_fn
        enddo ; enddo
    !$OMP do
        do j=js-1,je ;  do i=is,ie
          res_fn = 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
          kh_v(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
          khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j))) * res_fn
        enddo ; enddo
    !$OMP end parallel
      else
    !$OMP parallel default(none) shared(is,ie,js,je,Kh_u,Kh_v,khdt_x,khdt_y,CS,G,dt)
        if (cs%id_KhTr_u > 0) then
    !$OMP do
          do j=js,je ; do i=is-1,ie
            kh_u(i,j) = cs%KhTr
            khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
          enddo ; enddo
        else
    !$OMP do
          do j=js,je ; do i=is-1,ie
            khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
          enddo ; enddo
        endif
        if (cs%id_KhTr_v > 0) then
    !$OMP do
          do j=js-1,je ;  do i=is,ie
            kh_v(i,j) = cs%KhTr
            khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
          enddo ; enddo
        else
    !$OMP do
          do j=js-1,je ;  do i=is,ie
            khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
          enddo ; enddo
        endif
    !$OMP end parallel
      endif ! VarMix
  else ! .not. do_online
    khdt_x = read_khdt_x
    khdt_y = read_khdt_y
    call pass_vector(khdt_x,khdt_y,g%Domain)
  endif ! do_online



  if (cs%check_diffusive_CFL) then
    if (cs%show_call_tree) call calltree_waypoint("Checking diffusive CFL (tracer_hordiff)")
    max_cfl = 0.0
    do j=js,je ; do i=is,ie
      cfl(i,j) = 2.0*((khdt_x(i-1,j) + khdt_x(i,j)) + &
                      (khdt_y(i,j-1) + khdt_y(i,j))) * g%IareaT(i,j)
      if (max_cfl < cfl(i,j)) max_cfl = cfl(i,j)
    enddo ; enddo
    call cpu_clock_begin(id_clock_sync)
    call max_across_pes(max_cfl)
    call cpu_clock_end(id_clock_sync)
    num_itts = max(1,ceiling(max_cfl))
    i_numitts = 1.0 ; if (num_itts > 1) i_numitts = 1.0 / (real(num_itts))
    if(cs%id_CFL > 0) call post_data(cs%id_CFL, cfl, cs%diag, mask=g%mask2dT)
  else
    num_itts = 1 ; i_numitts = 1.0
  endif

  do m=1,ntr
    if (associated(reg%Tr(m)%df_x)) then
      do k=1,nz ; do j=js,je ; do i=is-1,ie
        reg%Tr(m)%df_x(i,j,k) = 0.0
      enddo ; enddo ; enddo
    endif
    if (associated(reg%Tr(m)%df_y)) then
      do k=1,nz ; do j=js-1,je ; do i=is,ie
        reg%Tr(m)%df_y(i,j,k) = 0.0
      enddo ; enddo ; enddo
    endif
    if (associated(reg%Tr(m)%df2d_x)) then
      do j=js,je ; do i=is-1,ie ; reg%Tr(m)%df2d_x(i,j) = 0.0 ; enddo ; enddo
    endif
    if (associated(reg%Tr(m)%df2d_y)) then
      do j=js-1,je ; do i=is,ie ; reg%Tr(m)%df2d_y(i,j) = 0.0 ; enddo ; enddo
    endif
  enddo

  if (cs%use_neutral_diffusion) then

    if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)")
    call cpu_clock_begin(id_clock_pass)
    call do_group_pass(cs%pass_t, g%Domain)
    call cpu_clock_end(id_clock_pass)
    ! We are assuming that neutral surfaces do not evolve (much) as a result of multiple
    ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs()
    ! would be inside the itt-loop. -AJA

    call neutral_diffusion_calc_coeffs(g, gv, h, tv%T, tv%S, tv%eqn_of_state, &
                                       cs%neutral_diffusion_CSp)
    do j=js-1,je ; do i=is,ie
      coef_y(i,j) = i_numitts * khdt_y(i,j)
    enddo ; enddo
    do j=js,je
      do i=is-1,ie
        coef_x(i,j) = i_numitts * khdt_x(i,j)
      enddo
    enddo

    do itt=1,num_itts
      if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt)
      if (itt>1) then ! Update halos for subsequent iterations
        call cpu_clock_begin(id_clock_pass)
        call do_group_pass(cs%pass_t, g%Domain)
        call cpu_clock_end(id_clock_pass)
      endif
      do m=1,ntr ! for each tracer
        call neutral_diffusion(g, gv,  h, coef_x, coef_y, reg%Tr(m)%t, m, i_numitts*dt, &
                               reg%Tr(m)%name, cs%neutral_diffusion_CSp)
      enddo ! m
    enddo ! itt

  else    ! following if not using neutral diffusion, but instead along-surface diffusion

    if (cs%show_call_tree) call calltree_waypoint("Calculating horizontal diffusion (tracer_hordiff)")
    do itt=1,num_itts
      call cpu_clock_begin(id_clock_pass)
      call do_group_pass(cs%pass_t, g%Domain)
      call cpu_clock_end(id_clock_pass)
!$OMP parallel do default(none) shared(is,ie,js,je,nz,I_numitts,CS,G,GV,khdt_y,h, &
!$OMP                                    h_neglect,khdt_x,ntr,Idt,Reg)           &
!$OMP                            private(scale,Coef_y,Coef_x,Ihdxdy,dTr)
      do k=1,nz
        scale = i_numitts
        if (cs%Diffuse_ML_interior) then
          if (k<=gv%nkml) then
            if (cs%ML_KhTr_scale <= 0.0) cycle
            scale = i_numitts * cs%ML_KhTr_scale
          endif
          if ((k>gv%nkml) .and. (k<=gv%nk_rho_varies)) cycle
        endif

        do j=js-1,je ; do i=is,ie
          coef_y(i,j) = ((scale * khdt_y(i,j))*2.0*(h(i,j,k)*h(i,j+1,k))) / &
                                                   (h(i,j,k)+h(i,j+1,k)+h_neglect)
        enddo ; enddo

        do j=js,je
          do i=is-1,ie
            coef_x(i,j) = ((scale * khdt_x(i,j))*2.0*(h(i,j,k)*h(i+1,j,k))) / &
                                                     (h(i,j,k)+h(i+1,j,k)+h_neglect)
          enddo

          do i=is,ie
            ihdxdy(i,j) = g%IareaT(i,j) / (h(i,j,k)+h_neglect)
          enddo
        enddo

        do m=1,ntr
          do j=js,je ; do i=is,ie
            dtr(i,j) = ihdxdy(i,j) * &
              ((coef_x(i-1,j) * (reg%Tr(m)%t(i-1,j,k) - reg%Tr(m)%t(i,j,k)) - &
                coef_x(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))) + &
               (coef_y(i,j-1) * (reg%Tr(m)%t(i,j-1,k) - reg%Tr(m)%t(i,j,k)) - &
                coef_y(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))))
          enddo ; enddo
          if (associated(reg%Tr(m)%df_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
            reg%Tr(m)%df_x(i,j,k) = reg%Tr(m)%df_x(i,j,k) + coef_x(i,j) * &
                                (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))*idt
          enddo ; enddo ; endif
          if (associated(reg%Tr(m)%df_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
            reg%Tr(m)%df_y(i,j,k) = reg%Tr(m)%df_y(i,j,k) + coef_y(i,j) * &
                                (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))*idt
          enddo ; enddo ; endif
          if (associated(reg%Tr(m)%df2d_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
            reg%Tr(m)%df2d_x(i,j) = reg%Tr(m)%df2d_x(i,j) + coef_x(i,j) * &
                                (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))*idt
          enddo ; enddo ; endif
          if (associated(reg%Tr(m)%df2d_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
            reg%Tr(m)%df2d_y(i,j) = reg%Tr(m)%df2d_y(i,j) + coef_y(i,j) * &
                                (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))*idt
          enddo ; enddo ; endif
          do j=js,je ; do i=is,ie
            reg%Tr(m)%t(i,j,k) = reg%Tr(m)%t(i,j,k) + dtr(i,j)
          enddo ; enddo
        enddo

      enddo ! End of k loop.

    enddo ! End of "while" loop.

  endif   ! endif for CS%use_neutral_diffusion
  call cpu_clock_end(id_clock_diffuse)


  if (cs%Diffuse_ML_interior) then
    if (cs%show_call_tree) call calltree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)")
    if (cs%debug) call mom_tracer_chksum("Before epipycnal diff ", reg%Tr, ntr, g)

    call cpu_clock_begin(id_clock_epimix)
    call tracer_epipycnal_ml_diff(h, dt, reg%Tr, ntr, khdt_x, khdt_y, g, gv, &
                                  cs, tv, num_itts)
    call cpu_clock_end(id_clock_epimix)
  endif

  if (cs%debug) call mom_tracer_chksum("After tracer diffusion ", reg%Tr, ntr, g)

  ! post diagnostics for 2d tracer diffusivity
  if (cs%id_KhTr_u > 0) then
    do j=js,je ; do i=is-1,ie
      kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
    enddo ; enddo
    call post_data(cs%id_KhTr_u, kh_u, cs%diag, mask=g%mask2dCu)
  endif
  if (cs%id_KhTr_v > 0) then
    do j=js-1,je ; do i=is,ie
      kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
    enddo ; enddo
    call post_data(cs%id_KhTr_v, kh_v, cs%diag, mask=g%mask2dCv)
  endif
  if (cs%id_KhTr_h > 0) then
    kh_h(:,:) = 0.0
    do j=js,je ; do i=is-1,ie
      kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
    enddo ; enddo
    do j=js,je ; do i=is,ie
      normalize = 1.0 / ((g%mask2dCu(i-1,j)+g%mask2dCu(i,j)) + &
                  (g%mask2dCv(i,j-1)+g%mask2dCv(i,j)) + gv%H_subroundoff)
      kh_h(i,j) = normalize*g%mask2dT(i,j)*((kh_u(i-1,j)+kh_u(i,j)) + &
                                            (kh_v(i,j-1)+kh_v(i,j)))
    enddo ; enddo
    call post_data(cs%id_KhTr_h, kh_h, cs%diag, mask=g%mask2dT)
  endif


  if (cs%debug) then
    call uvchksum("After tracer diffusion khdt_[xy]", &
                  khdt_x, khdt_y, g%HI, haloshift=0, symmetric=.true.)
    call uvchksum("After tracer diffusion Coef_[xy]", &
                  coef_x, coef_y, g%HI, haloshift=0, symmetric=.true.)
  endif

  if (cs%id_khdt_x > 0) call post_data(cs%id_khdt_x, khdt_x, cs%diag)
  if (cs%id_khdt_y > 0) call post_data(cs%id_khdt_y, khdt_y, cs%diag)

  if (cs%show_call_tree) call calltree_leave("tracer_hordiff()")

end subroutine tracer_hordiff

!> This subroutine does epipycnal diffusion of all tracers between the mixed
!! and buffer layers and the interior, using the diffusivity in CS%KhTr.
!! Multiple iterations are used (if necessary) so that there is no limit on the
!! acceptable time increment.
subroutine tracer_epipycnal_ml_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
                                    GV, CS, tv, num_itts)
  type(ocean_grid_type),                    intent(inout) :: G          !< ocean grid structure
  type(verticalgrid_type),                  intent(in)    :: GV         !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h          !< layer thickness (m or kg m-2)
  real,                                     intent(in)    :: dt         !< time step
  type(tracer_type),                        intent(inout) :: Tr(:)      !< tracer array
  integer,                                  intent(in)    :: ntr        !< number of tracers
  real, dimension(SZIB_(G),SZJ_(G)),        intent(in)    :: khdt_epi_x !< needs a comment
  real, dimension(SZI_(G),SZJB_(G)),        intent(in)    :: khdt_epi_y !< needs a comment
  type(tracer_hor_diff_cs),                 intent(inout) :: CS         !< module control structure
  type(thermo_var_ptrs),                    intent(in)    :: tv         !< thermodynamic structure
  integer,                                  intent(in)    :: num_itts   !< number of iterations (usually=1)


  real, dimension(SZI_(G), SZJ_(G)) :: &
    Rml_max  ! The maximum coordinate density within the mixed layer, in kg m-3.
  real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: &
    rho_coord ! The coordinate density that is used to mix along, in kg m-3.

  ! The naming mnemnonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point.
  ! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage.
  type(p2d), dimension(SZJ_(G)) :: &
    deep_wt_Lu, deep_wt_Ru, &  ! The relative weighting of the deeper of a pair, ND.
    hP_Lu, hP_Ru       ! The total thickness on each side for each pair, in m or kg m-2.

  type(p2d), dimension(SZJB_(G)) :: &
    deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair, ND.
    hP_Lv, hP_Rv       ! The total thickness on each side for each pair, in m or kg m-2.

  type(p2di), dimension(SZJ_(G)) :: &
    k0b_Lu, k0a_Lu, &  ! The original k-indices of the layers that participate
    k0b_Ru, k0a_Ru     ! in each pair of mixing at u-faces.
  type(p2di), dimension(SZJB_(G)) :: &
    k0b_Lv, k0a_Lv, &  ! The original k-indices of the layers that participate
    k0b_Rv, k0a_Rv     ! in each pair of mixing at v-faces.

  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
    tr_flux_conv  ! The flux convergence of tracers, in TR m3 or TR kg.
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R

  real, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
    rho_srt, & ! The density of each layer of the sorted columns, in kg m-3.
    h_srt      ! The thickness of each layer of the sorted columns, in m or kg m-2.
  integer, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
    k0_srt     ! The original k-index that each layer of the sorted column
               ! corresponds to.

  real, dimension(SZK_(G)) :: &
    h_demand_L, & ! The thickness in the left (_L) or right (_R) column that
    h_demand_R, & ! is demanded to match the thickness in the counterpart, in H.
    h_used_L, &   ! The summed thickness from the left or right columns that
    h_used_R, &   ! have actually been used, in m or kg m-2 (H).
    h_supply_frac_L, &  ! The fraction of the demanded thickness that can
    h_supply_frac_R     ! actually be supplied from a layer.
  integer, dimension(SZK_(G)) :: &
    kbs_Lp, &   ! The sorted indicies of the Left and Right columns for
    kbs_Rp      ! each pairing.

  integer, dimension(SZI_(G), SZJ_(G))  :: &
    num_srt, &   ! The number of layers that are sorted in each column.
    k_end_srt, & ! The maximum index in each column that might need to be
                 ! sorted, based on neighboring values of max_kRho
    max_krho     ! The index of the layer whose target density is just denser
                 ! than the densest part of the mixed layer.
  integer, dimension(SZJ_(G))           :: &
    max_srt      ! The maximum value of num_srt in a k-row.
  integer, dimension(SZIB_(G), SZJ_(G)) :: &
    nPu          ! The number of epipycnal pairings at each u-point.
  integer, dimension(SZI_(G), SZJB_(G)) :: &
    nPv          ! The number of epipycnal pairings at each v-point.
  real :: h_exclude    ! A thickness that layers must attain to be considered
                       ! for inclusion in mixing, in m.
  real :: Idt        ! The inverse of the time step, in s-1.
  real :: I_maxitt   ! The inverse of the maximum number of iterations.
  real :: rho_pair, rho_a, rho_b  ! Temporary densities, in kg m-3.
  real :: Tr_min_face  ! The minimum and maximum tracer concentrations
  real :: Tr_max_face  ! associated with a pairing, in conc.
  real :: Tr_La, Tr_Lb ! The 4 tracer concentrations that might be
  real :: Tr_Ra, Tr_Rb ! associated with a pairing, in conc.
  real :: Tr_av_L    ! The average tracer concentrations on the left and right
  real :: Tr_av_R    ! sides of a pairing, in conc.
  real :: Tr_flux    ! The tracer flux from left to right in a pair, in conc m3.
  real :: Tr_adj_vert  ! A downward vertical adjustment to Tr_flux between the
                     ! two cells that make up one side of the pairing, in conc m3.
  real :: h_L, h_R   ! Thicknesses to the left and right, in m or kg m-2 (H).
  real :: wt_a, wt_b ! Fractional weights of layers above and below, ND.
  real :: vol        ! A cell volume or mass, in m3 or kg (H m2).
  logical, dimension(SZK_(G)) :: &
    left_set, &  ! If true, the left or right point determines the density of
    right_set    ! of the trio.  If densities are exactly equal, both are true.
  real :: tmp
  real :: p_ref_cv(szi_(g))

  integer :: k_max, k_min, k_test, itmp
  integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb
  integer :: isd, ied, jsd, jed, IsdB, IedB, k_size
  integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt
  integer :: PEmax_kRho
  integer :: isv, iev, jsv, jev ! The valid range of the indices.

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB
  idt = 1.0/dt
  nkmb = gv%nk_rho_varies

  if (num_itts <= 1) then
    max_itt = 1 ; i_maxitt = 1.0
  else
    max_itt = num_itts ; i_maxitt = 1.0 / (real(max_itt))
  endif

  do i=is-2,ie+2 ; p_ref_cv(i) = tv%P_Ref ; enddo

  call cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_t, g%Domain)
  call cpu_clock_end(id_clock_pass)
  ! Determine which layers the mixed- and buffer-layers map into...
!$OMP parallel do default(none) shared(nkmb,is,ie,js,je,tv,p_ref_cv,rho_coord)
  do k=1,nkmb
    do j=js-2,je+2
       call calculate_density(tv%T(:,j,k),tv%S(:,j,k), p_ref_cv, &
                           rho_coord(:,j,k), is-2, ie-is+5, tv%eqn_of_state)
    enddo
 enddo

  do j=js-2,je+2 ; do i=is-2,ie+2
    rml_max(i,j) = rho_coord(i,j,1)
    num_srt(i,j) = 0 ; max_krho(i,j) = 0
  enddo ; enddo
  do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2
    if (rml_max(i,j) < rho_coord(i,j,k)) rml_max(i,j) = rho_coord(i,j,k)
  enddo ; enddo ; enddo
  !   Use bracketing and bisection to find the k-level that the densest of the
  ! mixed and buffer layer corresponds to, such that:
  !     GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
!$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,G,GV,Rml_max,max_kRho) &
!$OMP                          private(k_min,k_max,k_test)
  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (g%mask2dT(i,j) > 0.5) then
    if (rml_max(i,j) > gv%Rlay(nz)) then ; max_krho(i,j) = nz+1
    elseif (rml_max(i,j) <= gv%Rlay(nkmb+1)) then ; max_krho(i,j) = nkmb+1
    else
      k_min = nkmb+2 ; k_max = nz
      do
        k_test = (k_min + k_max) / 2
        if (rml_max(i,j) <= gv%Rlay(k_test-1)) then ; k_max = k_test-1
        elseif (gv%Rlay(k_test) < rml_max(i,j)) then ; k_min = k_test+1
        else ; max_krho(i,j) = k_test ; exit ; endif

        if (k_min == k_max) then ; max_krho(i,j) = k_max ; exit ; endif
      enddo
    endif
  endif ; enddo ; enddo

  pemax_krho = 0
  do j=js-1,je+1 ; do i=is-1,ie+1
    k_end_srt(i,j) = max(max_krho(i,j), max_krho(i-1,j), max_krho(i+1,j), &
                         max_krho(i,j-1), max_krho(i,j+1))
    if (pemax_krho < k_end_srt(i,j)) pemax_krho = k_end_srt(i,j)
  enddo ; enddo
  if (pemax_krho > nz) pemax_krho = nz ! PEmax_kRho could have been nz+1.

  h_exclude = 10.0*(gv%Angstrom + gv%H_subroundoff)
!$OMP parallel default(none) shared(is,ie,js,je,nkmb,G,GV,h,h_exclude,num_srt,k0_srt, &
!$OMP                               rho_srt,h_srt,PEmax_kRho,k_end_srt,rho_coord,max_srt) &
!$OMP                       private(ns,tmp,itmp)
!$OMP do
  do j=js-1,je+1
    do k=1,nkmb ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
      if (h(i,j,k) > h_exclude) then
        num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
        k0_srt(i,ns,j) = k
        rho_srt(i,ns,j) = rho_coord(i,j,k)
        h_srt(i,ns,j) = h(i,j,k)
      endif
    endif ; enddo ; enddo
    do k=nkmb+1,pemax_krho ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
      if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
        num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
        k0_srt(i,ns,j) = k
        rho_srt(i,ns,j) = gv%Rlay(k)
        h_srt(i,ns,j) = h(i,j,k)
      endif
    endif ; enddo ; enddo
  enddo
  ! Sort each column by increasing density.  This should already be close,
  ! and the size of the arrays are small, so straight insertion is used.
!$OMP do
   do j=js-1,je+1; do i=is-1,ie+1
    do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then
      ! The last segment needs to be shuffled earlier in the list.
      do k2 = k,2,-1 ; if (rho_srt(i,k2,j) >= rho_srt(i,k2-1,j)) exit
        itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp
        tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp
        tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp
      enddo
    endif ; enddo
  enddo; enddo
!$OMP do
  do j=js-1,je+1
    max_srt(j) = 0
    do i=is-1,ie+1 ; max_srt(j) = max(max_srt(j), num_srt(i,j)) ; enddo
  enddo
!$OMP end parallel

  do j=js,je
    k_size = max(2*max_srt(j),1)
    allocate(deep_wt_lu(j)%p(isdb:iedb,k_size))
    allocate(deep_wt_ru(j)%p(isdb:iedb,k_size))
    allocate(hp_lu(j)%p(isdb:iedb,k_size))
    allocate(hp_ru(j)%p(isdb:iedb,k_size))
    allocate(k0a_lu(j)%p(isdb:iedb,k_size))
    allocate(k0a_ru(j)%p(isdb:iedb,k_size))
    allocate(k0b_lu(j)%p(isdb:iedb,k_size))
    allocate(k0b_ru(j)%p(isdb:iedb,k_size))
  enddo

!$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lu,k0_srt, &
!$OMP                                  k0b_Ru,k0a_Lu,k0a_Ru,deep_wt_Lu,deep_wt_Ru,  &
!$OMP                                  h_srt,nkmb,nPu,hP_Lu,hP_Ru)                  &
!$OMP                          private(h_demand_L,h_used_L,h_demand_R,h_used_R,     &
!$OMP                                  kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
!$OMP                                  wt_b,left_set,right_set,h_supply_frac_R,     &
!$OMP                                  h_supply_frac_L)
  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
    ! Set up the pairings for fluxes through the zonal faces.

    do k=1,num_srt(i,j)   ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
    do k=1,num_srt(i+1,j) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo

    ! First merge the left and right lists into a single, sorted list.

    !   Discard any layers that are lighter than the lightest in the other
    ! column.  They can only participate in mixing as the lighter part of a
    ! pair of points.
    if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then
      kr = 1
      do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i+1,1,j)) exit ; enddo
    elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then
      kl = 1
      do kr=2,num_srt(i+1,j) ; if (rho_srt(i+1,kr,j) >= rho_srt(i,1,j)) exit ; enddo
    else
      kl = 1 ; kr = 1
    endif
    np = 0
    do ! Loop to accumulate pairs of columns.
      if ((kl > num_srt(i,j)) .or. (kr > num_srt(i+1,j))) exit

      if (rho_srt(i,kl,j) > rho_srt(i+1,kr,j)) then
      ! The right point is lighter and defines the density for this trio.
        np = np+1 ; k = np
        rho_pair = rho_srt(i+1,kr,j)

        k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
        k0a_lu(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
        kbs_lp(k) = kl ; kbs_rp(k) = kr

        rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
        wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
          wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
        deep_wt_lu(j)%p(i,k) = wt_b ; deep_wt_ru(j)%p(i,k) = 1.0

        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j) * wt_b
        h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i+1,kr,j) * (1.0-wt_b)

        kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
      elseif (rho_srt(i,kl,j) < rho_srt(i+1,kr,j)) then
      ! The left point is lighter and defines the density for this trio.
        np = np+1 ; k = np
        rho_pair = rho_srt(i,kl,j)
        k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
        k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0_srt(i+1,kr-1,j)

        kbs_lp(k) = kl ; kbs_rp(k) = kr

        rho_a = rho_srt(i+1,kr-1,j) ; rho_b = rho_srt(i+1,kr,j)
        wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
          wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
        deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = wt_b

        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
        h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)

        kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
      elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i+1,kr,j) <= nkmb)) then
        ! The densities are exactly equal and one layer is above the interior.
        np = np+1 ; k = np
        k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
        k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
        kbs_lp(k) = kl ; kbs_rp(k) = kr
        deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = 1.0

        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)

        kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
      else ! The densities are exactly equal and in the interior.
        ! Mixing in this case has already occurred, so accumulate the thickness
        ! demanded for that mixing and skip onward.
        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)

        kl = kl+1 ; kr = kr+1
      endif
    enddo ! Loop to accumulate pairs of columns.
    npu(i,j) = np ! This is the number of active pairings.

    ! Determine what fraction of the thickness "demand" can be supplied.
    do k=1,num_srt(i+1,j)
      h_supply_frac_r(k) = 1.0
      if (h_demand_r(k) > 0.5*h_srt(i+1,k,j)) &
        h_supply_frac_r(k) = 0.5*h_srt(i+1,k,j) / h_demand_r(k)
    enddo
    do k=1,num_srt(i,j)
      h_supply_frac_l(k) = 1.0
      if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
        h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
    enddo

    !  Distribute the "exported" thicknesses proportionately.
    do k=1,npu(i,j)
      kl = kbs_lp(k) ; kr = kbs_rp(k)
      hp_lu(j)%p(i,k) = 0.0 ; hp_ru(j)%p(i,k) = 0.0
      if (left_set(k)) then ! Add the contributing thicknesses on the right.
        if (deep_wt_ru(j)%p(i,k) < 1.0) then
          hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * min(h_supply_frac_r(kr), h_supply_frac_r(kr-1))
          wt_b = deep_wt_ru(j)%p(i,k)
          h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b)*hp_ru(j)%p(i,k)
          h_used_r(kr) = h_used_r(kr) + wt_b*hp_ru(j)%p(i,k)
        else
          hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
          h_used_r(kr) = h_used_r(kr) + hp_ru(j)%p(i,k)
        endif
      endif
      if (right_set(k)) then ! Add the contributing thicknesses on the left.
        if (deep_wt_lu(j)%p(i,k) < 1.0) then
          hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * min(h_supply_frac_l(kl), h_supply_frac_l(kl-1))
          wt_b = deep_wt_lu(j)%p(i,k)
          h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b)*hp_lu(j)%p(i,k)
          h_used_l(kl) = h_used_l(kl) + wt_b*hp_lu(j)%p(i,k)
        else
          hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * h_supply_frac_l(kl)
          h_used_l(kl) = h_used_l(kl) + hp_lu(j)%p(i,k)
        endif
      endif
    enddo

    !   The left-over thickness (at least half the layer thickness) is now
    ! added to the thicknesses of the importing columns.
    do k=1,npu(i,j)
      if (left_set(k)) hp_lu(j)%p(i,k) = hp_lu(j)%p(i,k) + &
                           (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
      if (right_set(k)) hp_ru(j)%p(i,k) = hp_ru(j)%p(i,k) + &
                            (h_srt(i+1,kbs_rp(k),j) - h_used_r(kbs_rp(k)))
    enddo

  endif ; enddo ; enddo ! i- & j- loops over zonal faces.

  do j=js-1,je
    k_size = max(max_srt(j)+max_srt(j+1),1)
    allocate(deep_wt_lv(j)%p(isd:ied,k_size))
    allocate(deep_wt_rv(j)%p(isd:ied,k_size))
    allocate(hp_lv(j)%p(isd:ied,k_size))
    allocate(hp_rv(j)%p(isd:ied,k_size))
    allocate(k0a_lv(j)%p(isd:ied,k_size))
    allocate(k0a_rv(j)%p(isd:ied,k_size))
    allocate(k0b_lv(j)%p(isd:ied,k_size))
    allocate(k0b_rv(j)%p(isd:ied,k_size))
  enddo

!$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lv,k0b_Rv, &
!$OMP                                  k0_srt,k0a_Lv,k0a_Rv,deep_wt_Lv,deep_wt_Rv,  &
!$OMP                                  h_srt,nkmb,nPv,hP_Lv,hP_Rv)                  &
!$OMP                          private(h_demand_L,h_used_L,h_demand_R,h_used_R,     &
!$OMP                                  kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
!$OMP                                  wt_b,left_set,right_set,h_supply_frac_R,     &
!$OMP                                  h_supply_frac_L)
  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
    ! Set up the pairings for fluxes through the meridional faces.

    do k=1,num_srt(i,j)   ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
    do k=1,num_srt(i,j+1) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo

    ! First merge the left and right lists into a single, sorted list.

    !   Discard any layers that are lighter than the lightest in the other
    ! column.  They can only participate in mixing as the lighter part of a
    ! pair of points.
    if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then
      kr = 1
      do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i,1,j+1)) exit ; enddo
    elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then
      kl = 1
      do kr=2,num_srt(i,j+1) ; if (rho_srt(i,kr,j+1) >= rho_srt(i,1,j)) exit ; enddo
    else
      kl = 1 ; kr = 1
    endif
    np = 0
    do ! Loop to accumulate pairs of columns.
      if ((kl > num_srt(i,j)) .or. (kr > num_srt(i,j+1))) exit

      if (rho_srt(i,kl,j) > rho_srt(i,kr,j+1)) then
      ! The right point is lighter and defines the density for this trio.
        np = np+1 ; k = np
        rho_pair = rho_srt(i,kr,j+1)

        k0b_lv(j)%p(i,k) = k0_srt(i,kl,j)   ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
        k0a_lv(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
        kbs_lp(k) = kl ; kbs_rp(k) = kr

        rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
        wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
          wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
        deep_wt_lv(j)%p(i,k) = wt_b ; deep_wt_rv(j)%p(i,k) = 1.0

        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1) * wt_b
        h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i,kr,j+1) * (1.0-wt_b)

        kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
      elseif (rho_srt(i,kl,j) < rho_srt(i,kr,j+1)) then
      ! The left point is lighter and defines the density for this trio.
        np = np+1 ; k = np
        rho_pair = rho_srt(i,kl,j)
        k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
        k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k) ; k0a_rv(j)%p(i,k) = k0_srt(i,kr-1,j+1)

        kbs_lp(k) = kl ; kbs_rp(k) = kr

        rho_a = rho_srt(i,kr-1,j+1) ; rho_b = rho_srt(i,kr,j+1)
        wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
          wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
        deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = wt_b

        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
        h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)

        kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
      elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i,kr,j+1) <= nkmb)) then
        ! The densities are exactly equal and one layer is above the interior.
        np = np+1 ; k = np
        k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
        k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k)  ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
        kbs_lp(k) = kl ; kbs_rp(k) = kr
        deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = 1.0

        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)

        kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
      else ! The densities are exactly equal and in the interior.
        ! Mixing in this case has already occurred, so accumulate the thickness
        ! demanded for that mixing and skip onward.
        h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
        h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)

        kl = kl+1 ; kr = kr+1
      endif
    enddo ! Loop to accumulate pairs of columns.
    npv(i,j) = np ! This is the number of active pairings.

    ! Determine what fraction of the thickness "demand" can be supplied.
    do k=1,num_srt(i,j+1)
      h_supply_frac_r(k) = 1.0
      if (h_demand_r(k) > 0.5*h_srt(i,k,j+1)) &
        h_supply_frac_r(k) = 0.5*h_srt(i,k,j+1) / h_demand_r(k)
    enddo
    do k=1,num_srt(i,j)
      h_supply_frac_l(k) = 1.0
      if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
        h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
    enddo

    !  Distribute the "exported" thicknesses proportionately.
    do k=1,npv(i,j)
      kl = kbs_lp(k) ; kr = kbs_rp(k)
      hp_lv(j)%p(i,k) = 0.0 ; hp_rv(j)%p(i,k) = 0.0
      if (left_set(k)) then ! Add the contributing thicknesses on the right.
        if (deep_wt_rv(j)%p(i,k) < 1.0) then
          hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * min(h_supply_frac_r(kr), h_supply_frac_r(kr-1))
          wt_b = deep_wt_rv(j)%p(i,k)
          h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b) * hp_rv(j)%p(i,k)
          h_used_r(kr) = h_used_r(kr) + wt_b * hp_rv(j)%p(i,k)
        else
          hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
          h_used_r(kr) = h_used_r(kr) + hp_rv(j)%p(i,k)
        endif
      endif
      if (right_set(k)) then ! Add the contributing thicknesses on the left.
        if (deep_wt_lv(j)%p(i,k) < 1.0) then
          hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * min(h_supply_frac_l(kl), h_supply_frac_l(kl-1))
          wt_b = deep_wt_lv(j)%p(i,k)
          h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b) * hp_lv(j)%p(i,k)
          h_used_l(kl) = h_used_l(kl) + wt_b * hp_lv(j)%p(i,k)
        else
          hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * h_supply_frac_l(kl)
          h_used_l(kl) = h_used_l(kl) + hp_lv(j)%p(i,k)
        endif
      endif
    enddo

    !   The left-over thickness (at least half the layer thickness) is now
    ! added to the thicknesses of the importing columns.
    do k=1,npv(i,j)
      if (left_set(k)) hp_lv(j)%p(i,k) = hp_lv(j)%p(i,k) + &
                            (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
      if (right_set(k)) hp_rv(j)%p(i,k) = hp_rv(j)%p(i,k) + &
                             (h_srt(i,kbs_rp(k),j+1) - h_used_r(kbs_rp(k)))
    enddo


  endif ; enddo ; enddo ! i- & j- loops over meridional faces.

! The tracer-specific calculations start here.

  ! Zero out tracer tendencies.
  do k=1,pemax_krho ; do j=js-1,je+1 ; do i=is-1,ie+1
    tr_flux_conv(i,j,k) = 0.0
  enddo ; enddo ; enddo

  do itt=1,max_itt

    if (itt > 1) then ! The halos have already been filled if itt==1.
      call cpu_clock_begin(id_clock_pass)
      call do_group_pass(cs%pass_t, g%Domain)
      call cpu_clock_end(id_clock_pass)
    endif

    do m=1,ntr
!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, &
!$OMP                                  k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru,   &
!$OMP                                  hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) &
!$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
!$OMP                                     Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
!$OMP                                     Tr_flux,Tr_adj_vert,wt_a,vol)
      do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
        ! Determine the fluxes through the zonal faces.

        ! Find the acceptable range of tracer concentration around this face.
        if (npu(i,j) >= 1) then
          tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
          tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
          do k=2,nkmb
            tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
            tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
          enddo

          ! Include the next two layers denser than the densest buffer layer.
          kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
          klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
          kra = nkmb+1 ; if (max_krho(i+1,j) < nz+1) kra = max_krho(i+1,j)
          krb = kra ; if (max_krho(i+1,j) < nz) krb = max_krho(i+1,j)+1
          tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
          if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
          if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
          if (h(i+1,j,kra) > h_exclude) tr_ra = tr(m)%t(i+1,j,kra)
          if (h(i+1,j,krb) > h_exclude) tr_rb = tr(m)%t(i+1,j,krb)
          tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
          tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)

          ! Include all points in diffusive pairings at this face.
          do k=1,npu(i,j)
            tr_lb = tr(m)%t(i,j,k0b_lu(j)%p(i,k))
            tr_rb = tr(m)%t(i+1,j,k0b_ru(j)%p(i,k))
            tr_la = tr_lb ; tr_ra = tr_rb
            if (deep_wt_lu(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lu(j)%p(i,k))
            if (deep_wt_ru(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i+1,j,k0a_ru(j)%p(i,k))
            tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
            tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
          enddo
        endif

        do k=1,npu(i,j)
          klb = k0b_lu(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
          if (deep_wt_lu(j)%p(i,k) < 1.0) then
            kla = k0a_lu(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
            wt_b = deep_wt_lu(j)%p(i,k)
            tr_av_l = wt_b*tr_lb + (1.0-wt_b)*tr_la
          endif

          krb = k0b_ru(j)%p(i,k) ; tr_rb = tr(m)%t(i+1,j,krb) ; tr_av_r = tr_rb
          if (deep_wt_ru(j)%p(i,k) < 1.0) then
            kra = k0a_ru(j)%p(i,k) ; tr_ra = tr(m)%t(i+1,j,kra)
            wt_b = deep_wt_ru(j)%p(i,k)
            tr_av_r = wt_b*tr_rb + (1.0-wt_b)*tr_ra
          endif

          h_l = hp_lu(j)%p(i,k) ; h_r = hp_ru(j)%p(i,k)
          tr_flux = i_maxitt * khdt_epi_x(i,j) * (tr_av_l - tr_av_r) * &
            ((2.0 * h_l * h_r) / (h_l + h_r))


          if (deep_wt_lu(j)%p(i,k) >= 1.0) then
            tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux
          else
            tr_adj_vert = 0.0
            wt_b = deep_wt_lu(j)%p(i,k) ; wt_a = 1.0 - wt_b
            vol = hp_lu(j)%p(i,k) * g%areaT(i,j)

            !   Ensure that the tracer flux does not drive the tracer values
            ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
            ! does that the concentration in both contributing peices exceed
            ! this range equally. With downgradient fluxes and the initial tracer
            ! concentrations determining the valid range, the latter condition
            ! only enters for large values of the effective diffusive CFL number.
            if (tr_flux > 0.0) then
              if (tr_la < tr_lb) then ; if (vol*(tr_la-tr_min_face) < tr_flux) &
                tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
                                          (vol*wt_b) * (tr_lb - tr_la))
              else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
                tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
                                         (vol*wt_a) * (tr_la - tr_lb))
              endif
            elseif (tr_flux < 0.0) then
              if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
                tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
                                         (vol*wt_b) * (tr_la - tr_lb))
              else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
                tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
                                          (vol*wt_a)*(tr_lb - tr_la))
              endif
            endif

            tr_flux_conv(i,j,kla) = tr_flux_conv(i,j,kla) - (wt_a*tr_flux + tr_adj_vert)
            tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - (wt_b*tr_flux - tr_adj_vert)
          endif

          if (deep_wt_ru(j)%p(i,k) >= 1.0) then
            tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + tr_flux
          else
            tr_adj_vert = 0.0
            wt_b = deep_wt_ru(j)%p(i,k) ; wt_a = 1.0 - wt_b
            vol = hp_ru(j)%p(i,k) * g%areaT(i+1,j)

            !   Ensure that the tracer flux does not drive the tracer values
            ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
            ! does that the concentration in both contributing peices exceed
            ! this range equally. With downgradient fluxes and the initial tracer
            ! concentrations determining the valid range, the latter condition
            ! only enters for large values of the effective diffusive CFL number.
            if (tr_flux < 0.0) then
              if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
                tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
                                          (vol*wt_b) * (tr_rb - tr_ra))
              else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
                tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
                                         (vol*wt_a) * (tr_ra - tr_rb))
              endif
            elseif (tr_flux > 0.0) then
              if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
                tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
                                         (vol*wt_b) * (tr_ra - tr_rb))
              else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
                tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
                                          (vol*wt_a)*(tr_rb - tr_ra))
              endif
            endif

            tr_flux_conv(i+1,j,kra) = tr_flux_conv(i+1,j,kra) + &
                                            (wt_a*tr_flux - tr_adj_vert)
            tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + &
                                            (wt_b*tr_flux + tr_adj_vert)
          endif
          if (associated(tr(m)%df2d_x)) &
            tr(m)%df2d_x(i,j) = tr(m)%df2d_x(i,j) + tr_flux * idt
        enddo ! Loop over pairings at faces.
      endif ; enddo ; enddo ! i- & j- loops over zonal faces.

!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, &
!$OMP                                  k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv,   &
!$OMP                                  hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d,          &
!$OMP                                  Tr_adj_vert_L,Tr_adj_vert_R,Idt)                     &
!$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,             &
!$OMP                                  Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,        &
!$OMP                                  h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
      do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
        ! Determine the fluxes through the meridional faces.

        ! Find the acceptable range of tracer concentration around this face.
        if (npv(i,j) >= 1) then
          tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
          tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
          do k=2,nkmb
            tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
            tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
          enddo

          ! Include the next two layers denser than the densest buffer layer.
          kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
          klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
          kra = nkmb+1 ; if (max_krho(i,j+1) < nz+1) kra = max_krho(i,j+1)
          krb = kra ; if (max_krho(i,j+1) < nz) krb = max_krho(i,j+1)+1
          tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
          if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
          if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
          if (h(i,j+1,kra) > h_exclude) tr_ra = tr(m)%t(i,j+1,kra)
          if (h(i,j+1,krb) > h_exclude) tr_rb = tr(m)%t(i,j+1,krb)
          tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
          tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)

          ! Include all points in diffusive pairings at this face.
          do k=1,npv(i,j)
            tr_lb = tr(m)%t(i,j,k0b_lv(j)%p(i,k)) ; tr_rb = tr(m)%t(i,j+1,k0b_rv(j)%p(i,k))
            tr_la = tr_lb ; tr_ra = tr_rb
            if (deep_wt_lv(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lv(j)%p(i,k))
            if (deep_wt_rv(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i,j+1,k0a_rv(j)%p(i,k))
            tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
            tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
          enddo
        endif

        do k=1,npv(i,j)
          klb = k0b_lv(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
          if (deep_wt_lv(j)%p(i,k) < 1.0) then
            kla = k0a_lv(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
            wt_b = deep_wt_lv(j)%p(i,k)
            tr_av_l = wt_b * tr_lb + (1.0-wt_b) * tr_la
          endif

          krb = k0b_rv(j)%p(i,k) ; tr_rb = tr(m)%t(i,j+1,krb) ; tr_av_r = tr_rb
          if (deep_wt_rv(j)%p(i,k) < 1.0) then
            kra = k0a_rv(j)%p(i,k) ; tr_ra = tr(m)%t(i,j+1,kra)
            wt_b = deep_wt_rv(j)%p(i,k)
            tr_av_r = wt_b * tr_rb + (1.0-wt_b) * tr_ra
          endif

          h_l = hp_lv(j)%p(i,k) ; h_r = hp_rv(j)%p(i,k)
          tr_flux = i_maxitt * ((2.0 * h_l * h_r) / (h_l + h_r)) * &
                    khdt_epi_y(i,j) * (tr_av_l - tr_av_r)
          tr_flux_3d(i,j,k) = tr_flux

          if (deep_wt_lv(j)%p(i,k) < 1.0) then
            tr_adj_vert = 0.0
            wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
            vol = hp_lv(j)%p(i,k) * g%areaT(i,j)

            !   Ensure that the tracer flux does not drive the tracer values
            ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
            if (tr_flux > 0.0) then
              if (tr_la < tr_lb) then ; if (vol * (tr_la-tr_min_face) < tr_flux) &
                tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
                                          (vol*wt_b) * (tr_lb - tr_la))
              else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
                tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
                                         (vol*wt_a) * (tr_la - tr_lb))
              endif
            elseif (tr_flux < 0.0) then
              if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
                tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
                                         (vol*wt_b) * (tr_la - tr_lb))
              else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
                tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
                                          (vol*wt_a)*(tr_lb - tr_la))
              endif
            endif
            tr_adj_vert_l(i,j,k) = tr_adj_vert
          endif

          if (deep_wt_rv(j)%p(i,k) < 1.0) then
            tr_adj_vert = 0.0
            wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
            vol = hp_rv(j)%p(i,k) * g%areaT(i,j+1)

            !   Ensure that the tracer flux does not drive the tracer values
            ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
            if (tr_flux < 0.0) then
              if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
                tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
                                          (vol*wt_b) * (tr_rb - tr_ra))
              else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
                tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
                                         (vol*wt_a) * (tr_ra - tr_rb))
              endif
            elseif (tr_flux > 0.0) then
              if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
                tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
                                         (vol*wt_b) * (tr_ra - tr_rb))
              else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
                tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
                                          (vol*wt_a)*(tr_rb - tr_ra))
              endif
            endif
            tr_adj_vert_r(i,j,k) = tr_adj_vert
          endif
          if (associated(tr(m)%df2d_y)) &
            tr(m)%df2d_y(i,j) = tr(m)%df2d_y(i,j) + tr_flux * idt
        enddo ! Loop over pairings at faces.
      endif ; enddo ; enddo ! i- & j- loops over meridional faces.
!$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv,  &
!$OMP                                  tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
!$OMP                                  deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
!$OMP                          private(kLa,kLb,kRa,kRb,wt_b,wt_a)
      do i=is,ie ; do j=js-1,je ; if (g%mask2dCv(i,j) > 0.5) then
        do k=1,npv(i,j)
          klb = k0b_lv(j)%p(i,k); krb = k0b_rv(j)%p(i,k)
          if (deep_wt_lv(j)%p(i,k) >= 1.0) then
            tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux_3d(i,j,k)
          else
            kla = k0a_lv(j)%p(i,k)
            wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
            tr_flux_conv(i,j,kla) = tr_flux_conv(i,j,kla) - (wt_a*tr_flux_3d(i,j,k) + tr_adj_vert_l(i,j,k))
            tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - (wt_b*tr_flux_3d(i,j,k) - tr_adj_vert_l(i,j,k))
          endif
          if (deep_wt_rv(j)%p(i,k) >= 1.0) then
            tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + tr_flux_3d(i,j,k)
          else
            kra = k0a_rv(j)%p(i,k)
            wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
            tr_flux_conv(i,j+1,kra) = tr_flux_conv(i,j+1,kra) + &
                                            (wt_a*tr_flux_3d(i,j,k) - tr_adj_vert_r(i,j,k))
            tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + &
                                            (wt_b*tr_flux_3d(i,j,k) + tr_adj_vert_r(i,j,k))
          endif
        enddo
      endif ; enddo ; enddo
!$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m)
      do k=1,pemax_krho ; do j=js,je ; do i=is,ie
        if ((g%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 0.0)) then
          tr(m)%t(i,j,k) = tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
                                            (h(i,j,k)*g%areaT(i,j))
          tr_flux_conv(i,j,k) = 0.0
        endif
      enddo ; enddo ; enddo

    enddo ! Loop over tracers
  enddo ! Loop over iterations

  do j=js,je
    deallocate(deep_wt_lu(j)%p) ; deallocate(deep_wt_ru(j)%p)
    deallocate(hp_lu(j)%p)  ; deallocate(hp_ru(j)%p)
    deallocate(k0a_lu(j)%p) ; deallocate(k0a_ru(j)%p)
    deallocate(k0b_lu(j)%p) ; deallocate(k0b_ru(j)%p)
  enddo

  do j=js-1,je
    deallocate(deep_wt_lv(j)%p) ; deallocate(deep_wt_rv(j)%p)
    deallocate(hp_lv(j)%p)  ; deallocate(hp_rv(j)%p)
    deallocate(k0a_lv(j)%p) ; deallocate(k0a_rv(j)%p)
    deallocate(k0b_lv(j)%p) ; deallocate(k0b_rv(j)%p)
  enddo

end subroutine tracer_epipycnal_ml_diff


!> Initialize lateral tracer diffusion module
subroutine tracer_hor_diff_init(Time, G, param_file, diag, CS, CSnd)
  type(time_type), target,    intent(in)    :: Time       !< current model time
  type(ocean_grid_type),      intent(in)    :: G          !< ocean grid structure
  type(diag_ctrl), target,    intent(inout) :: diag       !< diagnostic control
  type(param_file_type),      intent(in)    :: param_file !< parameter file
  type(tracer_hor_diff_cs),   pointer       :: CS         !< horz diffusion control structure
  type(neutral_diffusion_cs), pointer       :: CSnd       !< pointer to neutral diffusion CS

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_tracer_hor_diff" ! This module's name.
  character(len=256) :: mesg    ! Message for error messages.

  if (associated(cs)) then
    call mom_error(warning, "tracer_hor_diff_init called with associated control structure.")
    return
  endif
  allocate(cs)

  cs%diag => diag
  cs%show_call_tree = calltree_showquery()

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "KHTR", cs%KhTr, &
                 "The background along-isopycnal tracer diffusivity.", &
                 units="m2 s-1", default=0.0)
  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", cs%KhTr_Slope_Cff, &
                 "The scaling coefficient for along-isopycnal tracer \n"//&
                 "diffusivity using a shear-based (Visbeck-like) \n"//&
                 "parameterization.  A non-zero value enables this param.", &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "KHTR_MIN", cs%KhTr_Min, &
                 "The minimum along-isopycnal tracer diffusivity.", &
                 units="m2 s-1", default=0.0)
  call get_param(param_file, mdl, "KHTR_MAX", cs%KhTr_Max, &
                 "The maximum along-isopycnal tracer diffusivity.", &
                 units="m2 s-1", default=0.0)
  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", cs%KhTr_passivity_coeff, &
               "The coefficient that scales deformation radius over \n"//&
               "grid-spacing in passivity, where passiviity is the ratio \n"//&
               "between along isopycnal mxiing of tracers to thickness mixing. \n"//&
               "A non-zero value enables this parameterization.", &
               units="nondim", default=0.0)
  call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", cs%KhTr_passivity_min, &
               "The minimum passivity which is the ratio between \n"//&
               "along isopycnal mxiing of tracers to thickness mixing. \n", &
               units="nondim", default=0.5)
  call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
          desc="The (baroclinic) dynamics time step.", units="s")


  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
                 "If true, enable epipycnal mixing between the surface \n"//&
                 "boundary layer and the interior.", default=.false.)
  call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", cs%check_diffusive_CFL, &
                 "If true, use enough iterations the diffusion to ensure \n"//&
                 "that the diffusive equivalent of the CFL limit is not \n"//&
                 "violated.  If false, always use 1 iteration.", default=.false.)
  cs%ML_KhTR_scale = 1.0
  if (cs%Diffuse_ML_interior) then
    call get_param(param_file, mdl, "ML_KHTR_SCALE", cs%ML_KhTR_scale, &
                 "With Diffuse_ML_interior, the ratio of the truly \n"//&
                 "horizontal diffusivity in the mixed layer to the \n"//&
                 "epipycnal diffusivity.  The valid range is 0 to 1.", &
                 units="nondim", default=1.0)
  endif

  cs%use_neutral_diffusion = neutral_diffusion_init(time, g, param_file, diag, cs%neutral_diffusion_CSp)
  csnd => cs%neutral_diffusion_CSp
  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
       "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")

  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)

  id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)',          grain=clock_module)
  id_clock_epimix  = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=clock_module)
  id_clock_pass    = cpu_clock_id('(Ocean tracer halo updates)',     grain=clock_routine)
  id_clock_sync    = cpu_clock_id('(Ocean tracer global synch)',     grain=clock_routine)

  cs%id_KhTr_u = -1
  cs%id_KhTr_v = -1
  cs%id_KhTr_h = -1
  cs%id_CFL    = -1

  cs%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCu1, time, &
     'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'meter2 second-1')
  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
     'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'meter2 second-1')
  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time,&
     'Epipycnal tracer diffusivity at tracer cell center', 'meter2 second-1',   &
     cmor_field_name='diftrelo', cmor_units='m2 sec-1',                         &
     cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity',       &
     cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')

  cs%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, time, &
     'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'meter2')
  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
     'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'meter2')
  if (cs%check_diffusive_CFL) then
    cs%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, time,&
       'Grid CFL number for lateral/neutral tracer diffusion', 'dimensionless')
  endif


end subroutine tracer_hor_diff_init

subroutine tracer_hor_diff_end(CS)
  type(tracer_hor_diff_cs), pointer :: CS

  call neutral_diffusion_end(cs%neutral_diffusion_CSp)
  if (associated(cs)) deallocate(cs)

end subroutine tracer_hor_diff_end


!> \namespace mom_tracer_hor_diff
!!
!! \section section_intro Introduction to the module
!!
!!    This module contains subroutines that handle horizontal
!!  diffusion (i.e., isoneutral or along layer) of tracers.
!!
!!    Each of the tracers are subject to Fickian along-coordinate
!!  diffusion if Khtr is defined and positive.  The tracer diffusion
!!  can use a suitable number of iterations to guarantee stability
!!  with an arbitrarily large time step.

end module mom_tracer_hor_diff