mom_tracer_hor_diff Module Reference

Detailed Description

Main routine for lateral (along surface or neutral) diffusion of tracers.

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.

Data Types
type	p2d
	A type that can be used to create arrays of pointers to 2D arrays. More...
type	p2di
	A type that can be used to create arrays of pointers to 2D integer arrays. More...
type	tracer_hor_diff_cs
	The ocntrol structure for along-layer and epineutral tracer diffusion. More...
integer	id_clock_diffuse
	CPU time clocks.
integer	id_clock_epimix
	CPU time clocks.
integer	id_clock_pass
	CPU time clocks.
integer	id_clock_sync
	CPU time clocks.
subroutine, public	tracer_hordiff (h, dt, MEKE, VarMix, G, GV, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y)
	Compute along-coordinate diffusion of all tracers using the diffusivity in CSKhTr, or using space-dependent diffusivity. Multiple iterations are used (if necessary) so that there is no limit on the acceptable time increment. More...
subroutine	tracer_epipycnal_ml_diff (h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, GV, CS, tv, num_itts)
	This subroutine does epipycnal diffusion of all tracers between the mixed and buffer layers and the interior, using the diffusivity in CSKhTr. Multiple iterations are used (if necessary) so that there is no limit on the acceptable time increment. More...
subroutine, public	tracer_hor_diff_init (Time, G, param_file, diag, EOS, CS)
	Initialize lateral tracer diffusion module. More...
subroutine, public	tracer_hor_diff_end (CS)

Function/Subroutine Documentation

tracer_epipycnal_ml_diff()

subroutine mom_tracer_hor_diff::tracer_epipycnal_ml_diff ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h, real, intent(in)  dt, type(tracer_type), dimension(:), intent(inout)  Tr, integer, intent(in)  ntr, real, dimension(szib_(g),szj_(g)), intent(in)  khdt_epi_x, real, dimension(szi_(g),szjb_(g)), intent(in)  khdt_epi_y, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(tracer_hor_diff_cs), intent(inout)  CS, type(thermo_var_ptrs), intent(in)  tv, integer, intent(in)  num_itts  )

private

This subroutine does epipycnal diffusion of all tracers between the mixed and buffer layers and the interior, using the diffusivity in CSKhTr. Multiple iterations are used (if necessary) so that there is no limit on the acceptable time increment.

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h	layer thickness [H ~> m or kg m-2]
[in]	dt	time step
[in,out]	tr	tracer array
[in]	ntr	number of tracers
[in]	khdt_epi_x	needs a comment
[in]	khdt_epi_y	needs a comment
[in,out]	cs	module control structure
[in]	tv	thermodynamic structure
[in]	num_itts	number of iterations (usually=1)

Definition at line 541 of file MOM_tracer_hor_diff.F90.

  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 [H ~> 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 [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 [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 [nondim].
    hP_Lu, hP_Ru       ! The total thickness on each side for each pair [H ~> 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 [nondim].
    hP_Lv, hP_Rv       ! The total thickness on each side for each pair [H ~> 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 [conc H m2 ~> conc m3 or conc 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 [kg m-3].
    h_srt      ! The thickness of each layer of the sorted columns [H ~> 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 [H ~> m or kg m-2].
    h_used_L, &   ! The summed thickness from the left or right columns that
    h_used_R, &   ! have actually been used [H ~> m or kg m-2].
    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 [H ~> m or kg m-2].
  real :: Idt        ! The inverse of the time step [s-1].
  real :: I_maxitt   ! The inverse of the maximum number of iterations.
  real :: rho_pair, rho_a, rho_b  ! Temporary densities [kg m-3].
  real :: Tr_min_face  ! The minimum and maximum tracer concentrations
  real :: Tr_max_face  ! associated with a pairing [Conc]
  real :: Tr_La, Tr_Lb ! The 4 tracer concentrations that might be
  real :: Tr_Ra, Tr_Rb ! associated with a pairing [Conc]
  real :: Tr_av_L    ! The average tracer concentrations on the left and right
  real :: Tr_av_R    ! sides of a pairing [Conc].
  real :: Tr_flux    ! The tracer flux from left to right in a pair [conc H m2 ~> conc m3 or conc kg].
  real :: Tr_adj_vert  ! A downward vertical adjustment to Tr_flux between the
                     ! two cells that make up one side of the pairing [conc H m2 ~> conc m3 or conc kg].
  real :: h_L, h_R   ! Thicknesses to the left and right [H ~> m or kg m-2].
  real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim].
  real :: vol        ! A cell volume or mass [H m2 ~> m3 or kg].
  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 do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
  ! Determine which layers the mixed- and buffer-layers map into...
  !$OMP parallel do default(shared)
  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_H + 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 do_group_pass(cs%pass_t, g%Domain, clock=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
 

tracer_hor_diff_end()

subroutine, public mom_tracer_hor_diff::tracer_hor_diff_end

(

type(tracer_hor_diff_cs), pointer

CS

)

Parameters

cs	module control structure

Definition at line 1488 of file MOM_tracer_hor_diff.F90.

  type(tracer_hor_diff_CS), pointer :: CS !< module control structure
 
  call neutral_diffusion_end(cs%neutral_diffusion_CSp)
  if (associated(cs)) deallocate(cs)
 

tracer_hor_diff_init()

subroutine, public mom_tracer_hor_diff::tracer_hor_diff_init	(	type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(eos_type), intent(in), target	EOS,
		type(tracer_hor_diff_cs), pointer	CS
	)

Initialize lateral tracer diffusion module.

Parameters

[in]	time	current model time
[in]	g	ocean grid structure
[in,out]	diag	diagnostic control
[in]	eos	Equation of state CS
[in]	param_file	parameter file
	cs	horz diffusion control structure

Definition at line 1379 of file MOM_tracer_hor_diff.F90.

  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(EOS_type),  target,    intent(in)    :: EOS        !< Equation of state CS
  type(param_file_type),      intent(in)    :: param_file !< parameter file
  type(tracer_hor_diff_CS),   pointer       :: CS         !< horz diffusion control structure
 
! 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 "//&
                 "diffusivity using a shear-based (Visbeck-like) "//&
                 "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 "//&
                 "grid-spacing in passivity, where passivity is the ratio "//&
                 "between along isopycnal mixing of tracers to thickness mixing. "//&
                 "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 "//&
                 "along isopycnal mixing of tracers to thickness mixing.", &
                 units="nondim", default=0.5)
  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
                 "If true, enable epipycnal mixing between the surface "//&
                 "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 "//&
                 "that the diffusive equivalent of the CFL limit is not "//&
                 "violated.  If false, always use the greater of 1 or "//&
                 "MAX_TR_DIFFUSION_CFL iteration.", default=.false.)
  call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", cs%max_diff_CFL, &
                 "If positive, locally limit the along-isopycnal tracer "//&
                 "diffusivity to keep the diffusive CFL locally at or "//&
                 "below this value.  The number of diffusive iterations "//&
                 "is often this value or the next greater integer.", &
                 units="nondim", default=-1.0)
  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 "//&
                 "horizontal diffusivity in the mixed layer to the "//&
                 "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, eos, 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', 'm2 s-1')
  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
     'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1')
  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time,&
     'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1',   &
     cmor_field_name='diftrelo',                                                &
     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', 'm2')
  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
     'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2')
  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', 'nondim')
  endif
 
 

tracer_hordiff()

subroutine, public mom_tracer_hor_diff::tracer_hordiff	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		real, intent(in)	dt,
		type(meke_type), pointer	MEKE,
		type(varmix_cs), pointer	VarMix,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(tracer_hor_diff_cs), pointer	CS,
		type(tracer_registry_type), pointer	Reg,
		type(thermo_var_ptrs), intent(in)	tv,
		logical, intent(in), optional	do_online_flag,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional	read_khdt_x,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in), optional	read_khdt_y
	)

Compute along-coordinate diffusion of all tracers using the diffusivity in CSKhTr, or using space-dependent diffusivity. Multiple iterations are used (if necessary) so that there is no limit on the acceptable time increment.

Parameters

[in,out]	g	Grid type
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	dt	time step [s]
	meke	MEKE type
	varmix	Variable mixing type
[in]	gv	ocean vertical grid structure
	cs	module control structure
	reg	registered tracers
[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.
[in]	do_online_flag	If present and true, do online tracer transport with stored velcities.
[in]	read_khdt_x	If present, these are the zonal
[in]	read_khdt_y	If present, these are the meridional

Definition at line 98 of file MOM_tracer_hor_diff.F90.

  type(ocean_grid_type),      intent(inout) :: G       !< Grid type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                              intent(in)    :: h       !< Layer thickness [H ~> m or kg m-2]
  real,                       intent(in)    :: dt      !< time step [s]
  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), pointer       :: 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, intent(in)    :: do_online_flag !< If present and true, do online
                                                       !! tracer transport with stored velcities.
  real, dimension(SZIB_(G),SZJ_(G)), &
                    optional, intent(in)    :: read_khdt_x !< If present, these are the zonal
                                                       !! diffusivities from previous run.
  real, dimension(SZI_(G),SZJB_(G)), &
                    optional, intent(in)    :: read_khdt_y !< If present, these are the meridional
                                                       !! diffusivities from previous run.
 
 
  real, dimension(SZI_(G),SZJ_(G)) :: &
    Ihdxdy, &     ! The inverse of the volume or mass of fluid in a layer in a
                  ! grid cell [H-1 m-2 ~> m-3 or kg-1].
    kh_h, &       ! The tracer diffusivity averaged to tracer points [m2 s-1].
    cfl, &        ! A diffusive CFL number for each cell [nondim].
    dtr           ! The change in a tracer's concentration, in units of concentration [Conc].
 
  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 [m2].
    coef_x, &     ! The coefficients relating zonal tracer differences
                  ! to time-integrated fluxes [H m2 ~> m3 or kg].
    kh_u          ! Tracer mixing coefficient at u-points [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 [m2].
    coef_y, &     ! The coefficients relating meridional tracer differences
                  ! to time-integrated fluxes [H m2 ~> m3 or kg].
    kh_v          ! Tracer mixing coefficient at u-points [m2 s-1].
 
  real :: khdt_max ! The local limiting value of khdt_x or khdt_y [m2].
  real :: max_CFL ! The global maximum of the diffusive CFL number.
  logical :: use_VarMix, Resoln_scaled, do_online, use_Eady
  integer :: S_idx, T_idx ! Indices for temperature and salinity if needed
  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 [s-1].
  real :: h_neglect  ! A thickness that is so small it is usually lost
                     ! in roundoff and can be neglected [H ~> m or kg m-2].
  real :: Kh_loc     ! The local value of Kh [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 do default(shared) private(Kh_loc,Rd_dx)
      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 parallel do default(shared) private(Kh_loc,Rd_dx)
      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 parallel do default(shared)
      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 parallel do default(shared)
      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
    elseif (resoln_scaled) then
      !$OMP parallel do default(shared) private(Res_fn)
      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 parallel do default(shared) private(Res_fn)
      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
    else  ! Use a simple constant diffusivity.
      if (cs%id_KhTr_u > 0) then
        !$OMP parallel do default(shared)
        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 parallel do default(shared)
        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 parallel do default(shared)
        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 parallel do default(shared)
        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
    endif ! VarMix
 
    if (cs%max_diff_CFL > 0.0) then
      if ((cs%id_KhTr_u > 0) .or. (cs%id_KhTr_h > 0)) then
        !$OMP parallel do default(shared) private(khdt_max)
        do j=js,je ; do i=is-1,ie
          khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
          if (khdt_x(i,j) > khdt_max) then
            khdt_x(i,j) = khdt_max
            if (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)) > 0.0) &
              kh_u(i,j) = khdt_x(i,j) / (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
          endif
        enddo ; enddo
      else
        !$OMP parallel do default(shared) private(khdt_max)
        do j=js,je ; do i=is-1,ie
          khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
          khdt_x(i,j) = min(khdt_x(i,j), khdt_max)
        enddo ; enddo
      endif
      if ((cs%id_KhTr_v > 0) .or. (cs%id_KhTr_h > 0)) then
        !$OMP parallel do default(shared) private(khdt_max)
        do j=js-1,je ; do i=is,ie
          khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
          if (khdt_y(i,j) > khdt_max) then
            khdt_y(i,j) = khdt_max
            if (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)) > 0.0) &
              kh_v(i,j) = khdt_y(i,j) / (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
          endif
        enddo ; enddo
      else
        !$OMP parallel do default(shared) private(khdt_max)
        do j=js-1,je ; do i=is,ie
          khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
          khdt_y(i,j) = min(khdt_y(i,j), khdt_max)
        enddo ; enddo
      endif
    endif
 
  else ! .not. do_online
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is-1,ie
      khdt_x(i,j) = read_khdt_x(i,j)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js-1,je ;  do i=is,ie
      khdt_y(i,j) = read_khdt_y(i,j)
    enddo ; enddo
    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 - 4.0*epsilon(max_cfl)))
    i_numitts = 1.0 / (real(num_itts))
    if (cs%id_CFL > 0) call post_data(cs%id_CFL, cfl, cs%diag, mask=g%mask2dT)
  elseif (cs%max_diff_CFL > 0.0) then
    num_itts = max(1, ceiling(cs%max_diff_CFL - 4.0*epsilon(cs%max_diff_CFL)))
    i_numitts = 1.0 / (real(num_itts))
  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 do_group_pass(cs%pass_t, g%Domain, clock=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, 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 do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
      endif
      call neutral_diffusion(g, gv,  h, coef_x, coef_y, i_numitts*dt, reg, cs%neutral_diffusion_CSp)
    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 do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
      !$OMP parallel do default(shared) 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.)
    if (cs%use_neutral_diffusion) then
      call uvchksum("After tracer diffusion Coef_[xy]", coef_x, coef_y, &
                    g%HI, haloshift=0, symmetric=.true.)
    endif
  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()")