MOM6
Data Types
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. More...
 
integer id_clock_epimix
 CPU time clocks. More...
 
integer id_clock_pass
 CPU time clocks. More...
 
integer id_clock_sync
 CPU time clocks. More...
 
subroutine, public tracer_hordiff (h, dt, MEKE, VarMix, G, GV, US, 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, US, 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, US, param_file, diag, EOS, diabatic_CSp, 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(unit_scale_type), intent(in)  US,
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]gocean grid structure
[in]gvocean vertical grid structure
[in]hlayer thickness [H ~> m or kg m-2]
[in]dttime step [T ~> s]
[in,out]trtracer array
[in]ntrnumber of tracers
[in]khdt_epi_xZonal epipycnal diffusivity times a time step and the ratio of the open face width over the distance between adjacent tracer points [L2 ~> m2]
[in]khdt_epi_yMeridional epipycnal diffusivity times a time step and the ratio of the open face width over the distance between adjacent tracer points [L2 ~> m2]
[in]usA dimensional unit scaling type
[in,out]csmodule control structure
[in]tvthermodynamic structure
[in]num_ittsnumber of iterations (usually=1)

Definition at line 579 of file MOM_tracer_hor_diff.F90.

579  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
580  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
581  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
582  real, intent(in) :: dt !< time step [T ~> s]
583  type(tracer_type), intent(inout) :: Tr(:) !< tracer array
584  integer, intent(in) :: ntr !< number of tracers
585  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: khdt_epi_x !< Zonal epipycnal diffusivity times
586  !! a time step and the ratio of the open face width over
587  !! the distance between adjacent tracer points [L2 ~> m2]
588  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times
589  !! a time step and the ratio of the open face width over
590  !! the distance between adjacent tracer points [L2 ~> m2]
591  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
592  type(tracer_hor_diff_CS), intent(inout) :: CS !< module control structure
593  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure
594  integer, intent(in) :: num_itts !< number of iterations (usually=1)
595 
596 
597  real, dimension(SZI_(G), SZJ_(G)) :: &
598  Rml_max ! The maximum coordinate density within the mixed layer [R ~> kg m-3].
599  real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: &
600  rho_coord ! The coordinate density that is used to mix along [R ~> kg m-3].
601 
602  ! The naming mnemnonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point.
603  ! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage.
604  type(p2d), dimension(SZJ_(G)) :: &
605  deep_wt_Lu, deep_wt_Ru, & ! The relative weighting of the deeper of a pair [nondim].
606  hP_Lu, hP_Ru ! The total thickness on each side for each pair [H ~> m or kg m-2].
607 
608  type(p2d), dimension(SZJB_(G)) :: &
609  deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair [nondim].
610  hP_Lv, hP_Rv ! The total thickness on each side for each pair [H ~> m or kg m-2].
611 
612  type(p2di), dimension(SZJ_(G)) :: &
613  k0b_Lu, k0a_Lu, & ! The original k-indices of the layers that participate
614  k0b_Ru, k0a_Ru ! in each pair of mixing at u-faces.
615  type(p2di), dimension(SZJB_(G)) :: &
616  k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate
617  k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces.
618 
619  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
620  tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
621  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R
622 
623  real, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
624  rho_srt, & ! The density of each layer of the sorted columns [R ~> kg m-3].
625  h_srt ! The thickness of each layer of the sorted columns [H ~> m or kg m-2].
626  integer, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
627  k0_srt ! The original k-index that each layer of the sorted column
628  ! corresponds to.
629 
630  real, dimension(SZK_(G)) :: &
631  h_demand_L, & ! The thickness in the left (_L) or right (_R) column that
632  h_demand_R, & ! is demanded to match the thickness in the counterpart [H ~> m or kg m-2].
633  h_used_L, & ! The summed thickness from the left or right columns that
634  h_used_R, & ! have actually been used [H ~> m or kg m-2].
635  h_supply_frac_L, & ! The fraction of the demanded thickness that can
636  h_supply_frac_R ! actually be supplied from a layer.
637  integer, dimension(SZK_(G)) :: &
638  kbs_Lp, & ! The sorted indicies of the Left and Right columns for
639  kbs_Rp ! each pairing.
640 
641  integer, dimension(SZI_(G), SZJ_(G)) :: &
642  num_srt, & ! The number of layers that are sorted in each column.
643  k_end_srt, & ! The maximum index in each column that might need to be
644  ! sorted, based on neighboring values of max_kRho
645  max_krho ! The index of the layer whose target density is just denser
646  ! than the densest part of the mixed layer.
647  integer, dimension(SZJ_(G)) :: &
648  max_srt ! The maximum value of num_srt in a k-row.
649  integer, dimension(SZIB_(G), SZJ_(G)) :: &
650  nPu ! The number of epipycnal pairings at each u-point.
651  integer, dimension(SZI_(G), SZJB_(G)) :: &
652  nPv ! The number of epipycnal pairings at each v-point.
653  real :: h_exclude ! A thickness that layers must attain to be considered
654  ! for inclusion in mixing [H ~> m or kg m-2].
655  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
656  real :: I_maxitt ! The inverse of the maximum number of iterations.
657  real :: rho_pair, rho_a, rho_b ! Temporary densities [R ~> kg m-3].
658  real :: Tr_min_face ! The minimum and maximum tracer concentrations
659  real :: Tr_max_face ! associated with a pairing [Conc]
660  real :: Tr_La, Tr_Lb ! The 4 tracer concentrations that might be
661  real :: Tr_Ra, Tr_Rb ! associated with a pairing [Conc]
662  real :: Tr_av_L ! The average tracer concentrations on the left and right
663  real :: Tr_av_R ! sides of a pairing [Conc].
664  real :: Tr_flux ! The tracer flux from left to right in a pair [conc H L2 ~> conc m3 or conc kg].
665  real :: Tr_adj_vert ! A downward vertical adjustment to Tr_flux between the
666  ! two cells that make up one side of the pairing [conc H L2 ~> conc m3 or conc kg].
667  real :: h_L, h_R ! Thicknesses to the left and right [H ~> m or kg m-2].
668  real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim].
669  real :: vol ! A cell volume or mass [H L2 ~> m3 or kg].
670  logical, dimension(SZK_(G)) :: &
671  left_set, & ! If true, the left or right point determines the density of
672  right_set ! of the trio. If densities are exactly equal, both are true.
673  real :: tmp
674  real :: p_ref_cv(SZI_(G))
675 
676  integer :: k_max, k_min, k_test, itmp
677  integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb
678  integer :: isd, ied, jsd, jed, IsdB, IedB, k_size
679  integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt
680  integer :: PEmax_kRho
681  integer :: isv, iev, jsv, jev ! The valid range of the indices.
682 
683  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
684  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
685  isdb = g%IsdB ; iedb = g%IedB
686  idt = 1.0 / dt
687  nkmb = gv%nk_rho_varies
688 
689  if (num_itts <= 1) then
690  max_itt = 1 ; i_maxitt = 1.0
691  else
692  max_itt = num_itts ; i_maxitt = 1.0 / (real(max_itt))
693  endif
694 
695  do i=is-2,ie+2 ; p_ref_cv(i) = tv%P_Ref ; enddo
696 
697  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
698  ! Determine which layers the mixed- and buffer-layers map into...
699  !$OMP parallel do default(shared)
700  do k=1,nkmb ; do j=js-2,je+2
701  call calculate_density(tv%T(:,j,k),tv%S(:,j,k), p_ref_cv, &
702  rho_coord(:,j,k), is-2, ie-is+5, tv%eqn_of_state, scale=us%kg_m3_to_R)
703  enddo ; enddo
704 
705  do j=js-2,je+2 ; do i=is-2,ie+2
706  rml_max(i,j) = rho_coord(i,j,1)
707  num_srt(i,j) = 0 ; max_krho(i,j) = 0
708  enddo ; enddo
709  do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2
710  if (rml_max(i,j) < rho_coord(i,j,k)) rml_max(i,j) = rho_coord(i,j,k)
711  enddo ; enddo ; enddo
712  ! Use bracketing and bisection to find the k-level that the densest of the
713  ! mixed and buffer layer corresponds to, such that:
714  ! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
715 !$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,G,GV,Rml_max,max_kRho) &
716 !$OMP private(k_min,k_max,k_test)
717  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (g%mask2dT(i,j) > 0.5) then
718  if (rml_max(i,j) > gv%Rlay(nz)) then ; max_krho(i,j) = nz+1
719  elseif (rml_max(i,j) <= gv%Rlay(nkmb+1)) then ; max_krho(i,j) = nkmb+1
720  else
721  k_min = nkmb+2 ; k_max = nz
722  do
723  k_test = (k_min + k_max) / 2
724  if (rml_max(i,j) <= gv%Rlay(k_test-1)) then ; k_max = k_test-1
725  elseif (gv%Rlay(k_test) < rml_max(i,j)) then ; k_min = k_test+1
726  else ; max_krho(i,j) = k_test ; exit ; endif
727 
728  if (k_min == k_max) then ; max_krho(i,j) = k_max ; exit ; endif
729  enddo
730  endif
731  endif ; enddo ; enddo
732 
733  pemax_krho = 0
734  do j=js-1,je+1 ; do i=is-1,ie+1
735  k_end_srt(i,j) = max(max_krho(i,j), max_krho(i-1,j), max_krho(i+1,j), &
736  max_krho(i,j-1), max_krho(i,j+1))
737  if (pemax_krho < k_end_srt(i,j)) pemax_krho = k_end_srt(i,j)
738  enddo ; enddo
739  if (pemax_krho > nz) pemax_krho = nz ! PEmax_kRho could have been nz+1.
740 
741  h_exclude = 10.0*(gv%Angstrom_H + gv%H_subroundoff)
742 !$OMP parallel default(none) shared(is,ie,js,je,nkmb,G,GV,h,h_exclude,num_srt,k0_srt, &
743 !$OMP rho_srt,h_srt,PEmax_kRho,k_end_srt,rho_coord,max_srt) &
744 !$OMP private(ns,tmp,itmp)
745 !$OMP do
746  do j=js-1,je+1
747  do k=1,nkmb ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
748  if (h(i,j,k) > h_exclude) then
749  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
750  k0_srt(i,ns,j) = k
751  rho_srt(i,ns,j) = rho_coord(i,j,k)
752  h_srt(i,ns,j) = h(i,j,k)
753  endif
754  endif ; enddo ; enddo
755  do k=nkmb+1,pemax_krho ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
756  if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
757  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
758  k0_srt(i,ns,j) = k
759  rho_srt(i,ns,j) = gv%Rlay(k)
760  h_srt(i,ns,j) = h(i,j,k)
761  endif
762  endif ; enddo ; enddo
763  enddo
764  ! Sort each column by increasing density. This should already be close,
765  ! and the size of the arrays are small, so straight insertion is used.
766 !$OMP do
767  do j=js-1,je+1; do i=is-1,ie+1
768  do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then
769  ! The last segment needs to be shuffled earlier in the list.
770  do k2 = k,2,-1 ; if (rho_srt(i,k2,j) >= rho_srt(i,k2-1,j)) exit
771  itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp
772  tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp
773  tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp
774  enddo
775  endif ; enddo
776  enddo ; enddo
777 !$OMP do
778  do j=js-1,je+1
779  max_srt(j) = 0
780  do i=is-1,ie+1 ; max_srt(j) = max(max_srt(j), num_srt(i,j)) ; enddo
781  enddo
782 !$OMP end parallel
783 
784  do j=js,je
785  k_size = max(2*max_srt(j),1)
786  allocate(deep_wt_lu(j)%p(isdb:iedb,k_size))
787  allocate(deep_wt_ru(j)%p(isdb:iedb,k_size))
788  allocate(hp_lu(j)%p(isdb:iedb,k_size))
789  allocate(hp_ru(j)%p(isdb:iedb,k_size))
790  allocate(k0a_lu(j)%p(isdb:iedb,k_size))
791  allocate(k0a_ru(j)%p(isdb:iedb,k_size))
792  allocate(k0b_lu(j)%p(isdb:iedb,k_size))
793  allocate(k0b_ru(j)%p(isdb:iedb,k_size))
794  enddo
795 
796 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lu,k0_srt, &
797 !$OMP k0b_Ru,k0a_Lu,k0a_Ru,deep_wt_Lu,deep_wt_Ru, &
798 !$OMP h_srt,nkmb,nPu,hP_Lu,hP_Ru) &
799 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
800 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
801 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
802 !$OMP h_supply_frac_L)
803  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
804  ! Set up the pairings for fluxes through the zonal faces.
805 
806  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
807  do k=1,num_srt(i+1,j) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
808 
809  ! First merge the left and right lists into a single, sorted list.
810 
811  ! Discard any layers that are lighter than the lightest in the other
812  ! column. They can only participate in mixing as the lighter part of a
813  ! pair of points.
814  if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then
815  kr = 1
816  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i+1,1,j)) exit ; enddo
817  elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then
818  kl = 1
819  do kr=2,num_srt(i+1,j) ; if (rho_srt(i+1,kr,j) >= rho_srt(i,1,j)) exit ; enddo
820  else
821  kl = 1 ; kr = 1
822  endif
823  np = 0
824  do ! Loop to accumulate pairs of columns.
825  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i+1,j))) exit
826 
827  if (rho_srt(i,kl,j) > rho_srt(i+1,kr,j)) then
828  ! The right point is lighter and defines the density for this trio.
829  np = np+1 ; k = np
830  rho_pair = rho_srt(i+1,kr,j)
831 
832  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
833  k0a_lu(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
834  kbs_lp(k) = kl ; kbs_rp(k) = kr
835 
836  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
837  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
838  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
839  deep_wt_lu(j)%p(i,k) = wt_b ; deep_wt_ru(j)%p(i,k) = 1.0
840 
841  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j) * wt_b
842  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i+1,kr,j) * (1.0-wt_b)
843 
844  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
845  elseif (rho_srt(i,kl,j) < rho_srt(i+1,kr,j)) then
846  ! The left point is lighter and defines the density for this trio.
847  np = np+1 ; k = np
848  rho_pair = rho_srt(i,kl,j)
849  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
850  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)
851 
852  kbs_lp(k) = kl ; kbs_rp(k) = kr
853 
854  rho_a = rho_srt(i+1,kr-1,j) ; rho_b = rho_srt(i+1,kr,j)
855  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
856  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
857  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = wt_b
858 
859  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
860  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
861 
862  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
863  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i+1,kr,j) <= nkmb)) then
864  ! The densities are exactly equal and one layer is above the interior.
865  np = np+1 ; k = np
866  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
867  k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
868  kbs_lp(k) = kl ; kbs_rp(k) = kr
869  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = 1.0
870 
871  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
872  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
873 
874  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
875  else ! The densities are exactly equal and in the interior.
876  ! Mixing in this case has already occurred, so accumulate the thickness
877  ! demanded for that mixing and skip onward.
878  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
879  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
880 
881  kl = kl+1 ; kr = kr+1
882  endif
883  enddo ! Loop to accumulate pairs of columns.
884  npu(i,j) = np ! This is the number of active pairings.
885 
886  ! Determine what fraction of the thickness "demand" can be supplied.
887  do k=1,num_srt(i+1,j)
888  h_supply_frac_r(k) = 1.0
889  if (h_demand_r(k) > 0.5*h_srt(i+1,k,j)) &
890  h_supply_frac_r(k) = 0.5*h_srt(i+1,k,j) / h_demand_r(k)
891  enddo
892  do k=1,num_srt(i,j)
893  h_supply_frac_l(k) = 1.0
894  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
895  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
896  enddo
897 
898  ! Distribute the "exported" thicknesses proportionately.
899  do k=1,npu(i,j)
900  kl = kbs_lp(k) ; kr = kbs_rp(k)
901  hp_lu(j)%p(i,k) = 0.0 ; hp_ru(j)%p(i,k) = 0.0
902  if (left_set(k)) then ! Add the contributing thicknesses on the right.
903  if (deep_wt_ru(j)%p(i,k) < 1.0) then
904  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))
905  wt_b = deep_wt_ru(j)%p(i,k)
906  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b)*hp_ru(j)%p(i,k)
907  h_used_r(kr) = h_used_r(kr) + wt_b*hp_ru(j)%p(i,k)
908  else
909  hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
910  h_used_r(kr) = h_used_r(kr) + hp_ru(j)%p(i,k)
911  endif
912  endif
913  if (right_set(k)) then ! Add the contributing thicknesses on the left.
914  if (deep_wt_lu(j)%p(i,k) < 1.0) then
915  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))
916  wt_b = deep_wt_lu(j)%p(i,k)
917  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b)*hp_lu(j)%p(i,k)
918  h_used_l(kl) = h_used_l(kl) + wt_b*hp_lu(j)%p(i,k)
919  else
920  hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * h_supply_frac_l(kl)
921  h_used_l(kl) = h_used_l(kl) + hp_lu(j)%p(i,k)
922  endif
923  endif
924  enddo
925 
926  ! The left-over thickness (at least half the layer thickness) is now
927  ! added to the thicknesses of the importing columns.
928  do k=1,npu(i,j)
929  if (left_set(k)) hp_lu(j)%p(i,k) = hp_lu(j)%p(i,k) + &
930  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
931  if (right_set(k)) hp_ru(j)%p(i,k) = hp_ru(j)%p(i,k) + &
932  (h_srt(i+1,kbs_rp(k),j) - h_used_r(kbs_rp(k)))
933  enddo
934 
935  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
936 
937  do j=js-1,je
938  k_size = max(max_srt(j)+max_srt(j+1),1)
939  allocate(deep_wt_lv(j)%p(isd:ied,k_size))
940  allocate(deep_wt_rv(j)%p(isd:ied,k_size))
941  allocate(hp_lv(j)%p(isd:ied,k_size))
942  allocate(hp_rv(j)%p(isd:ied,k_size))
943  allocate(k0a_lv(j)%p(isd:ied,k_size))
944  allocate(k0a_rv(j)%p(isd:ied,k_size))
945  allocate(k0b_lv(j)%p(isd:ied,k_size))
946  allocate(k0b_rv(j)%p(isd:ied,k_size))
947  enddo
948 
949 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lv,k0b_Rv, &
950 !$OMP k0_srt,k0a_Lv,k0a_Rv,deep_wt_Lv,deep_wt_Rv, &
951 !$OMP h_srt,nkmb,nPv,hP_Lv,hP_Rv) &
952 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
953 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
954 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
955 !$OMP h_supply_frac_L)
956  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
957  ! Set up the pairings for fluxes through the meridional faces.
958 
959  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
960  do k=1,num_srt(i,j+1) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
961 
962  ! First merge the left and right lists into a single, sorted list.
963 
964  ! Discard any layers that are lighter than the lightest in the other
965  ! column. They can only participate in mixing as the lighter part of a
966  ! pair of points.
967  if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then
968  kr = 1
969  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i,1,j+1)) exit ; enddo
970  elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then
971  kl = 1
972  do kr=2,num_srt(i,j+1) ; if (rho_srt(i,kr,j+1) >= rho_srt(i,1,j)) exit ; enddo
973  else
974  kl = 1 ; kr = 1
975  endif
976  np = 0
977  do ! Loop to accumulate pairs of columns.
978  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i,j+1))) exit
979 
980  if (rho_srt(i,kl,j) > rho_srt(i,kr,j+1)) then
981  ! The right point is lighter and defines the density for this trio.
982  np = np+1 ; k = np
983  rho_pair = rho_srt(i,kr,j+1)
984 
985  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
986  k0a_lv(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
987  kbs_lp(k) = kl ; kbs_rp(k) = kr
988 
989  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
990  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
991  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
992  deep_wt_lv(j)%p(i,k) = wt_b ; deep_wt_rv(j)%p(i,k) = 1.0
993 
994  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1) * wt_b
995  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i,kr,j+1) * (1.0-wt_b)
996 
997  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
998  elseif (rho_srt(i,kl,j) < rho_srt(i,kr,j+1)) then
999  ! The left point is lighter and defines the density for this trio.
1000  np = np+1 ; k = np
1001  rho_pair = rho_srt(i,kl,j)
1002  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1003  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)
1004 
1005  kbs_lp(k) = kl ; kbs_rp(k) = kr
1006 
1007  rho_a = rho_srt(i,kr-1,j+1) ; rho_b = rho_srt(i,kr,j+1)
1008  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1009  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1010  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = wt_b
1011 
1012  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
1013  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
1014 
1015  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
1016  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i,kr,j+1) <= nkmb)) then
1017  ! The densities are exactly equal and one layer is above the interior.
1018  np = np+1 ; k = np
1019  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1020  k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
1021  kbs_lp(k) = kl ; kbs_rp(k) = kr
1022  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = 1.0
1023 
1024  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1025  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1026 
1027  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
1028  else ! The densities are exactly equal and in the interior.
1029  ! Mixing in this case has already occurred, so accumulate the thickness
1030  ! demanded for that mixing and skip onward.
1031  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1032  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1033 
1034  kl = kl+1 ; kr = kr+1
1035  endif
1036  enddo ! Loop to accumulate pairs of columns.
1037  npv(i,j) = np ! This is the number of active pairings.
1038 
1039  ! Determine what fraction of the thickness "demand" can be supplied.
1040  do k=1,num_srt(i,j+1)
1041  h_supply_frac_r(k) = 1.0
1042  if (h_demand_r(k) > 0.5*h_srt(i,k,j+1)) &
1043  h_supply_frac_r(k) = 0.5*h_srt(i,k,j+1) / h_demand_r(k)
1044  enddo
1045  do k=1,num_srt(i,j)
1046  h_supply_frac_l(k) = 1.0
1047  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
1048  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
1049  enddo
1050 
1051  ! Distribute the "exported" thicknesses proportionately.
1052  do k=1,npv(i,j)
1053  kl = kbs_lp(k) ; kr = kbs_rp(k)
1054  hp_lv(j)%p(i,k) = 0.0 ; hp_rv(j)%p(i,k) = 0.0
1055  if (left_set(k)) then ! Add the contributing thicknesses on the right.
1056  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1057  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))
1058  wt_b = deep_wt_rv(j)%p(i,k)
1059  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b) * hp_rv(j)%p(i,k)
1060  h_used_r(kr) = h_used_r(kr) + wt_b * hp_rv(j)%p(i,k)
1061  else
1062  hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
1063  h_used_r(kr) = h_used_r(kr) + hp_rv(j)%p(i,k)
1064  endif
1065  endif
1066  if (right_set(k)) then ! Add the contributing thicknesses on the left.
1067  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1068  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))
1069  wt_b = deep_wt_lv(j)%p(i,k)
1070  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b) * hp_lv(j)%p(i,k)
1071  h_used_l(kl) = h_used_l(kl) + wt_b * hp_lv(j)%p(i,k)
1072  else
1073  hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * h_supply_frac_l(kl)
1074  h_used_l(kl) = h_used_l(kl) + hp_lv(j)%p(i,k)
1075  endif
1076  endif
1077  enddo
1078 
1079  ! The left-over thickness (at least half the layer thickness) is now
1080  ! added to the thicknesses of the importing columns.
1081  do k=1,npv(i,j)
1082  if (left_set(k)) hp_lv(j)%p(i,k) = hp_lv(j)%p(i,k) + &
1083  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
1084  if (right_set(k)) hp_rv(j)%p(i,k) = hp_rv(j)%p(i,k) + &
1085  (h_srt(i,kbs_rp(k),j+1) - h_used_r(kbs_rp(k)))
1086  enddo
1087 
1088 
1089  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1090 
1091 ! The tracer-specific calculations start here.
1092 
1093  ! Zero out tracer tendencies.
1094  do k=1,pemax_krho ; do j=js-1,je+1 ; do i=is-1,ie+1
1095  tr_flux_conv(i,j,k) = 0.0
1096  enddo ; enddo ; enddo
1097 
1098  do itt=1,max_itt
1099 
1100  if (itt > 1) then ! The halos have already been filled if itt==1.
1101  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
1102  endif
1103 
1104  do m=1,ntr
1105 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, &
1106 !$OMP k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru, &
1107 !$OMP hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) &
1108 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
1109 !$OMP Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
1110 !$OMP Tr_flux,Tr_adj_vert,wt_a,vol)
1111  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
1112  ! Determine the fluxes through the zonal faces.
1113 
1114  ! Find the acceptable range of tracer concentration around this face.
1115  if (npu(i,j) >= 1) then
1116  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1117  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1118  do k=2,nkmb
1119  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1120  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1121  enddo
1122 
1123  ! Include the next two layers denser than the densest buffer layer.
1124  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1125  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1126  kra = nkmb+1 ; if (max_krho(i+1,j) < nz+1) kra = max_krho(i+1,j)
1127  krb = kra ; if (max_krho(i+1,j) < nz) krb = max_krho(i+1,j)+1
1128  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1129  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1130  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1131  if (h(i+1,j,kra) > h_exclude) tr_ra = tr(m)%t(i+1,j,kra)
1132  if (h(i+1,j,krb) > h_exclude) tr_rb = tr(m)%t(i+1,j,krb)
1133  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1134  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1135 
1136  ! Include all points in diffusive pairings at this face.
1137  do k=1,npu(i,j)
1138  tr_lb = tr(m)%t(i,j,k0b_lu(j)%p(i,k))
1139  tr_rb = tr(m)%t(i+1,j,k0b_ru(j)%p(i,k))
1140  tr_la = tr_lb ; tr_ra = tr_rb
1141  if (deep_wt_lu(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lu(j)%p(i,k))
1142  if (deep_wt_ru(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i+1,j,k0a_ru(j)%p(i,k))
1143  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1144  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1145  enddo
1146  endif
1147 
1148  do k=1,npu(i,j)
1149  klb = k0b_lu(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1150  if (deep_wt_lu(j)%p(i,k) < 1.0) then
1151  kla = k0a_lu(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1152  wt_b = deep_wt_lu(j)%p(i,k)
1153  tr_av_l = wt_b*tr_lb + (1.0-wt_b)*tr_la
1154  endif
1155 
1156  krb = k0b_ru(j)%p(i,k) ; tr_rb = tr(m)%t(i+1,j,krb) ; tr_av_r = tr_rb
1157  if (deep_wt_ru(j)%p(i,k) < 1.0) then
1158  kra = k0a_ru(j)%p(i,k) ; tr_ra = tr(m)%t(i+1,j,kra)
1159  wt_b = deep_wt_ru(j)%p(i,k)
1160  tr_av_r = wt_b*tr_rb + (1.0-wt_b)*tr_ra
1161  endif
1162 
1163  h_l = hp_lu(j)%p(i,k) ; h_r = hp_ru(j)%p(i,k)
1164  tr_flux = i_maxitt * khdt_epi_x(i,j) * (tr_av_l - tr_av_r) * &
1165  ((2.0 * h_l * h_r) / (h_l + h_r))
1166 
1167 
1168  if (deep_wt_lu(j)%p(i,k) >= 1.0) then
1169  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux
1170  else
1171  tr_adj_vert = 0.0
1172  wt_b = deep_wt_lu(j)%p(i,k) ; wt_a = 1.0 - wt_b
1173  vol = hp_lu(j)%p(i,k) * g%areaT(i,j)
1174 
1175  ! Ensure that the tracer flux does not drive the tracer values
1176  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1177  ! does that the concentration in both contributing peices exceed
1178  ! this range equally. With downgradient fluxes and the initial tracer
1179  ! concentrations determining the valid range, the latter condition
1180  ! only enters for large values of the effective diffusive CFL number.
1181  if (tr_flux > 0.0) then
1182  if (tr_la < tr_lb) then ; if (vol*(tr_la-tr_min_face) < tr_flux) &
1183  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1184  (vol*wt_b) * (tr_lb - tr_la))
1185  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1186  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1187  (vol*wt_a) * (tr_la - tr_lb))
1188  endif
1189  elseif (tr_flux < 0.0) then
1190  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1191  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1192  (vol*wt_b) * (tr_la - tr_lb))
1193  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1194  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1195  (vol*wt_a)*(tr_lb - tr_la))
1196  endif
1197  endif
1198 
1199  tr_flux_conv(i,j,kla) = tr_flux_conv(i,j,kla) - (wt_a*tr_flux + tr_adj_vert)
1200  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - (wt_b*tr_flux - tr_adj_vert)
1201  endif
1202 
1203  if (deep_wt_ru(j)%p(i,k) >= 1.0) then
1204  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + tr_flux
1205  else
1206  tr_adj_vert = 0.0
1207  wt_b = deep_wt_ru(j)%p(i,k) ; wt_a = 1.0 - wt_b
1208  vol = hp_ru(j)%p(i,k) * g%areaT(i+1,j)
1209 
1210  ! Ensure that the tracer flux does not drive the tracer values
1211  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1212  ! does that the concentration in both contributing peices exceed
1213  ! this range equally. With downgradient fluxes and the initial tracer
1214  ! concentrations determining the valid range, the latter condition
1215  ! only enters for large values of the effective diffusive CFL number.
1216  if (tr_flux < 0.0) then
1217  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1218  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1219  (vol*wt_b) * (tr_rb - tr_ra))
1220  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1221  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1222  (vol*wt_a) * (tr_ra - tr_rb))
1223  endif
1224  elseif (tr_flux > 0.0) then
1225  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1226  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1227  (vol*wt_b) * (tr_ra - tr_rb))
1228  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1229  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1230  (vol*wt_a)*(tr_rb - tr_ra))
1231  endif
1232  endif
1233 
1234  tr_flux_conv(i+1,j,kra) = tr_flux_conv(i+1,j,kra) + &
1235  (wt_a*tr_flux - tr_adj_vert)
1236  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + &
1237  (wt_b*tr_flux + tr_adj_vert)
1238  endif
1239  if (associated(tr(m)%df2d_x)) &
1240  tr(m)%df2d_x(i,j) = tr(m)%df2d_x(i,j) + tr_flux * idt
1241  enddo ! Loop over pairings at faces.
1242  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
1243 
1244 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, &
1245 !$OMP k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv, &
1246 !$OMP hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d, &
1247 !$OMP Tr_adj_vert_L,Tr_adj_vert_R,Idt) &
1248 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, &
1249 !$OMP Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, &
1250 !$OMP h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
1251  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
1252  ! Determine the fluxes through the meridional faces.
1253 
1254  ! Find the acceptable range of tracer concentration around this face.
1255  if (npv(i,j) >= 1) then
1256  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1257  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1258  do k=2,nkmb
1259  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1260  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1261  enddo
1262 
1263  ! Include the next two layers denser than the densest buffer layer.
1264  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1265  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1266  kra = nkmb+1 ; if (max_krho(i,j+1) < nz+1) kra = max_krho(i,j+1)
1267  krb = kra ; if (max_krho(i,j+1) < nz) krb = max_krho(i,j+1)+1
1268  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1269  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1270  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1271  if (h(i,j+1,kra) > h_exclude) tr_ra = tr(m)%t(i,j+1,kra)
1272  if (h(i,j+1,krb) > h_exclude) tr_rb = tr(m)%t(i,j+1,krb)
1273  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1274  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1275 
1276  ! Include all points in diffusive pairings at this face.
1277  do k=1,npv(i,j)
1278  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))
1279  tr_la = tr_lb ; tr_ra = tr_rb
1280  if (deep_wt_lv(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lv(j)%p(i,k))
1281  if (deep_wt_rv(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i,j+1,k0a_rv(j)%p(i,k))
1282  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1283  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1284  enddo
1285  endif
1286 
1287  do k=1,npv(i,j)
1288  klb = k0b_lv(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1289  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1290  kla = k0a_lv(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1291  wt_b = deep_wt_lv(j)%p(i,k)
1292  tr_av_l = wt_b * tr_lb + (1.0-wt_b) * tr_la
1293  endif
1294 
1295  krb = k0b_rv(j)%p(i,k) ; tr_rb = tr(m)%t(i,j+1,krb) ; tr_av_r = tr_rb
1296  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1297  kra = k0a_rv(j)%p(i,k) ; tr_ra = tr(m)%t(i,j+1,kra)
1298  wt_b = deep_wt_rv(j)%p(i,k)
1299  tr_av_r = wt_b * tr_rb + (1.0-wt_b) * tr_ra
1300  endif
1301 
1302  h_l = hp_lv(j)%p(i,k) ; h_r = hp_rv(j)%p(i,k)
1303  tr_flux = i_maxitt * ((2.0 * h_l * h_r) / (h_l + h_r)) * &
1304  khdt_epi_y(i,j) * (tr_av_l - tr_av_r)
1305  tr_flux_3d(i,j,k) = tr_flux
1306 
1307  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1308  tr_adj_vert = 0.0
1309  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1310  vol = hp_lv(j)%p(i,k) * g%areaT(i,j)
1311 
1312  ! Ensure that the tracer flux does not drive the tracer values
1313  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1314  if (tr_flux > 0.0) then
1315  if (tr_la < tr_lb) then ; if (vol * (tr_la-tr_min_face) < tr_flux) &
1316  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1317  (vol*wt_b) * (tr_lb - tr_la))
1318  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1319  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1320  (vol*wt_a) * (tr_la - tr_lb))
1321  endif
1322  elseif (tr_flux < 0.0) then
1323  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1324  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1325  (vol*wt_b) * (tr_la - tr_lb))
1326  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1327  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1328  (vol*wt_a)*(tr_lb - tr_la))
1329  endif
1330  endif
1331  tr_adj_vert_l(i,j,k) = tr_adj_vert
1332  endif
1333 
1334  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1335  tr_adj_vert = 0.0
1336  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1337  vol = hp_rv(j)%p(i,k) * g%areaT(i,j+1)
1338 
1339  ! Ensure that the tracer flux does not drive the tracer values
1340  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1341  if (tr_flux < 0.0) then
1342  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1343  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1344  (vol*wt_b) * (tr_rb - tr_ra))
1345  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1346  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1347  (vol*wt_a) * (tr_ra - tr_rb))
1348  endif
1349  elseif (tr_flux > 0.0) then
1350  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1351  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1352  (vol*wt_b) * (tr_ra - tr_rb))
1353  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1354  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1355  (vol*wt_a)*(tr_rb - tr_ra))
1356  endif
1357  endif
1358  tr_adj_vert_r(i,j,k) = tr_adj_vert
1359  endif
1360  if (associated(tr(m)%df2d_y)) &
1361  tr(m)%df2d_y(i,j) = tr(m)%df2d_y(i,j) + tr_flux * idt
1362  enddo ! Loop over pairings at faces.
1363  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1364 !$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv, &
1365 !$OMP tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
1366 !$OMP deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
1367 !$OMP private(kLa,kLb,kRa,kRb,wt_b,wt_a)
1368  do i=is,ie ; do j=js-1,je ; if (g%mask2dCv(i,j) > 0.5) then
1369  do k=1,npv(i,j)
1370  klb = k0b_lv(j)%p(i,k); krb = k0b_rv(j)%p(i,k)
1371  if (deep_wt_lv(j)%p(i,k) >= 1.0) then
1372  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux_3d(i,j,k)
1373  else
1374  kla = k0a_lv(j)%p(i,k)
1375  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1376  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))
1377  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))
1378  endif
1379  if (deep_wt_rv(j)%p(i,k) >= 1.0) then
1380  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + tr_flux_3d(i,j,k)
1381  else
1382  kra = k0a_rv(j)%p(i,k)
1383  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1384  tr_flux_conv(i,j+1,kra) = tr_flux_conv(i,j+1,kra) + &
1385  (wt_a*tr_flux_3d(i,j,k) - tr_adj_vert_r(i,j,k))
1386  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + &
1387  (wt_b*tr_flux_3d(i,j,k) + tr_adj_vert_r(i,j,k))
1388  endif
1389  enddo
1390  endif ; enddo ; enddo
1391 !$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m)
1392  do k=1,pemax_krho ; do j=js,je ; do i=is,ie
1393  if ((g%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 0.0)) then
1394  tr(m)%t(i,j,k) = tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
1395  (h(i,j,k)*g%areaT(i,j))
1396  tr_flux_conv(i,j,k) = 0.0
1397  endif
1398  enddo ; enddo ; enddo
1399 
1400  enddo ! Loop over tracers
1401  enddo ! Loop over iterations
1402 
1403  do j=js,je
1404  deallocate(deep_wt_lu(j)%p) ; deallocate(deep_wt_ru(j)%p)
1405  deallocate(hp_lu(j)%p) ; deallocate(hp_ru(j)%p)
1406  deallocate(k0a_lu(j)%p) ; deallocate(k0a_ru(j)%p)
1407  deallocate(k0b_lu(j)%p) ; deallocate(k0b_ru(j)%p)
1408  enddo
1409 
1410  do j=js-1,je
1411  deallocate(deep_wt_lv(j)%p) ; deallocate(deep_wt_rv(j)%p)
1412  deallocate(hp_lv(j)%p) ; deallocate(hp_rv(j)%p)
1413  deallocate(k0a_lv(j)%p) ; deallocate(k0a_rv(j)%p)
1414  deallocate(k0b_lv(j)%p) ; deallocate(k0b_rv(j)%p)
1415  enddo
1416 

References id_clock_pass.

Referenced by tracer_hordiff().

Here is the caller graph for this function:

◆ tracer_hor_diff_end()

subroutine, public mom_tracer_hor_diff::tracer_hor_diff_end ( type(tracer_hor_diff_cs), pointer  CS)
Parameters
csmodule control structure

Definition at line 1542 of file MOM_tracer_hor_diff.F90.

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

Referenced by mom::mom_end().

Here is the caller graph for this function:

◆ 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(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(eos_type), intent(in), target  EOS,
type(diabatic_cs), intent(in), pointer  diabatic_CSp,
type(tracer_hor_diff_cs), pointer  CS 
)

Initialize lateral tracer diffusion module.

Parameters
[in]timecurrent model time
[in]gocean grid structure
[in]usA dimensional unit scaling type
[in,out]diagdiagnostic control
[in]eosEquation of state CS
[in]diabatic_cspEquation of state CS
[in]param_fileparameter file
cshorz diffusion control structure

Definition at line 1422 of file MOM_tracer_hor_diff.F90.

1422  type(time_type), target, intent(in) :: Time !< current model time
1423  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1424  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1425  type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control
1426  type(EOS_type), target, intent(in) :: EOS !< Equation of state CS
1427  type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS
1428  type(param_file_type), intent(in) :: param_file !< parameter file
1429  type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure
1430 
1431 ! This include declares and sets the variable "version".
1432 #include "version_variable.h"
1433  character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name.
1434  character(len=256) :: mesg ! Message for error messages.
1435 
1436  if (associated(cs)) then
1437  call mom_error(warning, "tracer_hor_diff_init called with associated control structure.")
1438  return
1439  endif
1440  allocate(cs)
1441 
1442  cs%diag => diag
1443  cs%show_call_tree = calltree_showquery()
1444 
1445  ! Read all relevant parameters and write them to the model log.
1446  call log_version(param_file, mdl, version, "")
1447  call get_param(param_file, mdl, "KHTR", cs%KhTr, &
1448  "The background along-isopycnal tracer diffusivity.", &
1449  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1450  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", cs%KhTr_Slope_Cff, &
1451  "The scaling coefficient for along-isopycnal tracer "//&
1452  "diffusivity using a shear-based (Visbeck-like) "//&
1453  "parameterization. A non-zero value enables this param.", &
1454  units="nondim", default=0.0)
1455  call get_param(param_file, mdl, "KHTR_MIN", cs%KhTr_Min, &
1456  "The minimum along-isopycnal tracer diffusivity.", &
1457  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1458  call get_param(param_file, mdl, "KHTR_MAX", cs%KhTr_Max, &
1459  "The maximum along-isopycnal tracer diffusivity.", &
1460  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1461  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", cs%KhTr_passivity_coeff, &
1462  "The coefficient that scales deformation radius over "//&
1463  "grid-spacing in passivity, where passivity is the ratio "//&
1464  "between along isopycnal mixing of tracers to thickness mixing. "//&
1465  "A non-zero value enables this parameterization.", &
1466  units="nondim", default=0.0)
1467  call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", cs%KhTr_passivity_min, &
1468  "The minimum passivity which is the ratio between "//&
1469  "along isopycnal mixing of tracers to thickness mixing.", &
1470  units="nondim", default=0.5)
1471  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
1472  "If true, enable epipycnal mixing between the surface "//&
1473  "boundary layer and the interior.", default=.false.)
1474  call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", cs%check_diffusive_CFL, &
1475  "If true, use enough iterations the diffusion to ensure "//&
1476  "that the diffusive equivalent of the CFL limit is not "//&
1477  "violated. If false, always use the greater of 1 or "//&
1478  "MAX_TR_DIFFUSION_CFL iteration.", default=.false.)
1479  call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", cs%max_diff_CFL, &
1480  "If positive, locally limit the along-isopycnal tracer "//&
1481  "diffusivity to keep the diffusive CFL locally at or "//&
1482  "below this value. The number of diffusive iterations "//&
1483  "is often this value or the next greater integer.", &
1484  units="nondim", default=-1.0)
1485  call get_param(param_file, mdl, "RECALC_NEUTRAL_SURF", cs%recalc_neutral_surf, &
1486  "If true, then recalculate the neutral surfaces if the \n"//&
1487  "diffusive CFL is exceeded. If false, assume that the \n"//&
1488  "positions of the surfaces do not change \n", default = .false.)
1489  cs%ML_KhTR_scale = 1.0
1490  if (cs%Diffuse_ML_interior) then
1491  call get_param(param_file, mdl, "ML_KHTR_SCALE", cs%ML_KhTR_scale, &
1492  "With Diffuse_ML_interior, the ratio of the truly "//&
1493  "horizontal diffusivity in the mixed layer to the "//&
1494  "epipycnal diffusivity. The valid range is 0 to 1.", &
1495  units="nondim", default=1.0)
1496  endif
1497 
1498  cs%use_neutral_diffusion = neutral_diffusion_init(time, g, param_file, diag, eos, diabatic_csp, &
1499  cs%neutral_diffusion_CSp )
1500  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1501  "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1502  cs%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(time, g, param_file, diag, diabatic_csp, &
1503  cs%lateral_boundary_diffusion_CSp)
1504  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1505  "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1506 
1507  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1508 
1509  id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=clock_module)
1510  id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=clock_module)
1511  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1512  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1513 
1514  cs%id_KhTr_u = -1
1515  cs%id_KhTr_v = -1
1516  cs%id_KhTr_h = -1
1517  cs%id_CFL = -1
1518 
1519  cs%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCu1, time, &
1520  'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1521  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
1522  'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1523  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time, &
1524  'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
1525  cmor_field_name='diftrelo', &
1526  cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', &
1527  cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')
1528 
1529  cs%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, time, &
1530  'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1531  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
1532  'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1533  if (cs%check_diffusive_CFL) then
1534  cs%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, time,&
1535  'Grid CFL number for lateral/neutral tracer diffusion', 'nondim')
1536  endif
1537 
1538 

References id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync, and mom_diag_mediator::register_diag_field().

Here is the call graph for this function:

◆ 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(unit_scale_type), intent(in)  US,
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]gGrid type
[in]hLayer thickness [H ~> m or kg m-2]
[in]dttime step [T ~> s]
mekeMEKE type
varmixVariable mixing type
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
csmodule control structure
regregistered tracers
[in]tvA 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_flagIf present and true, do online tracer transport with stored velcities.
[in]read_khdt_xIf present, these are the zonal
[in]read_khdt_yIf present, these are the meridional

Definition at line 107 of file MOM_tracer_hor_diff.F90.

107  type(ocean_grid_type), intent(inout) :: G !< Grid type
108  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
109  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
110  real, intent(in) :: dt !< time step [T ~> s]
111  type(MEKE_type), pointer :: MEKE !< MEKE type
112  type(VarMix_CS), pointer :: VarMix !< Variable mixing type
113  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
114  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
115  type(tracer_hor_diff_CS), pointer :: CS !< module control structure
116  type(tracer_registry_type), pointer :: Reg !< registered tracers
117  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
118  !! thermodynamic fields, including potential temp and
119  !! salinity or mixed layer density. Absent fields have
120  !! NULL ptrs, and these may (probably will) point to
121  !! some of the same arrays as Tr does. tv is required
122  !! for epipycnal mixing between mixed layer and the interior.
123  ! Optional inputs for offline tracer transport
124  logical, optional, intent(in) :: do_online_flag !< If present and true, do online
125  !! tracer transport with stored velcities.
126  real, dimension(SZIB_(G),SZJ_(G)), &
127  optional, intent(in) :: read_khdt_x !< If present, these are the zonal
128  !! diffusivities from previous run.
129  real, dimension(SZI_(G),SZJB_(G)), &
130  optional, intent(in) :: read_khdt_y !< If present, these are the meridional
131  !! diffusivities from previous run.
132 
133 
134  real, dimension(SZI_(G),SZJ_(G)) :: &
135  Ihdxdy, & ! The inverse of the volume or mass of fluid in a layer in a
136  ! grid cell [H-1 L-2 ~> m-3 or kg-1].
137  kh_h, & ! The tracer diffusivity averaged to tracer points [L2 T-1 ~> m2 s-1].
138  cfl, & ! A diffusive CFL number for each cell [nondim].
139  dtr ! The change in a tracer's concentration, in units of concentration [Conc].
140 
141  real, dimension(SZIB_(G),SZJ_(G)) :: &
142  khdt_x, & ! The value of Khtr*dt times the open face width divided by
143  ! the distance between adjacent tracer points [L2 ~> m2].
144  coef_x, & ! The coefficients relating zonal tracer differences
145  ! to time-integrated fluxes [H L2 ~> m3 or kg].
146  kh_u ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
147  real, dimension(SZI_(G),SZJB_(G)) :: &
148  khdt_y, & ! The value of Khtr*dt times the open face width divided by
149  ! the distance between adjacent tracer points [L2 ~> m2].
150  coef_y, & ! The coefficients relating meridional tracer differences
151  ! to time-integrated fluxes [H L2 ~> m3 or kg].
152  kh_v ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
153 
154  real :: khdt_max ! The local limiting value of khdt_x or khdt_y [L2 ~> m2].
155  real :: max_CFL ! The global maximum of the diffusive CFL number.
156  logical :: use_VarMix, Resoln_scaled, do_online, use_Eady
157  integer :: S_idx, T_idx ! Indices for temperature and salinity if needed
158  integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts
159  real :: I_numitts ! The inverse of the number of iterations, num_itts.
160  real :: scale ! The fraction of khdt_x or khdt_y that is applied in this
161  ! layer for this iteration [nondim].
162  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
163  real :: h_neglect ! A thickness that is so small it is usually lost
164  ! in roundoff and can be neglected [H ~> m or kg m-2].
165  real :: Kh_loc ! The local value of Kh [L2 T-1 ~> m2 s-1].
166  real :: Res_Fn ! The local value of the resolution function [nondim].
167  real :: Rd_dx ! The local value of deformation radius over grid-spacing [nondim].
168  real :: normalize ! normalization used for diagnostic Kh_h; diffusivity averaged to h-points.
169 
170  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
171 
172  do_online = .true.
173  if (present(do_online_flag)) do_online = do_online_flag
174 
175  if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
176  "register_tracer must be called before tracer_hordiff.")
177  if (loc(reg)==0) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
178  "register_tracer must be called before tracer_hordiff.")
179  if ((reg%ntr==0) .or. ((cs%KhTr <= 0.0) .and. .not.associated(varmix)) ) return
180 
181  if (cs%show_call_tree) call calltree_enter("tracer_hordiff(), MOM_tracer_hor_diff.F90")
182 
183  call cpu_clock_begin(id_clock_diffuse)
184 
185  ntr = reg%ntr
186  idt = 1.0 / dt
187  h_neglect = gv%H_subroundoff
188 
189  if (cs%Diffuse_ML_interior .and. cs%first_call) then ; if (is_root_pe()) then
190  do m=1,ntr ; if (associated(reg%Tr(m)%df_x) .or. associated(reg%Tr(m)%df_y)) &
191  call mom_error(warning, "tracer_hordiff: Tracer "//trim(reg%Tr(m)%name)// &
192  " has associated 3-d diffusive flux diagnostics. These are not "//&
193  "valid when DIFFUSE_ML_TO_INTERIOR is defined. Use 2-d tracer "//&
194  "diffusion diagnostics instead to get accurate total fluxes.")
195  enddo
196  endif ; endif
197  cs%first_call = .false.
198 
199  if (cs%debug) call mom_tracer_chksum("Before tracer diffusion ", reg%Tr, ntr, g)
200 
201  use_varmix = .false. ; resoln_scaled = .false. ; use_eady = .false.
202  if (Associated(varmix)) then
203  use_varmix = varmix%use_variable_mixing
204  resoln_scaled = varmix%Resoln_scaled_KhTr
205  use_eady = cs%KhTr_Slope_Cff > 0.
206  endif
207 
208  call cpu_clock_begin(id_clock_pass)
209  do m=1,ntr
210  call create_group_pass(cs%pass_t, reg%Tr(m)%t(:,:,:), g%Domain)
211  enddo
212  call cpu_clock_end(id_clock_pass)
213 
214  if (cs%show_call_tree) call calltree_waypoint("Calculating diffusivity (tracer_hordiff)")
215 
216  if (do_online) then
217  if (use_varmix) then
218  !$OMP parallel do default(shared) private(Kh_loc,Rd_dx)
219  do j=js,je ; do i=is-1,ie
220  kh_loc = cs%KhTr
221  if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2u(i,j)*varmix%SN_u(i,j)
222  if (associated(meke%Kh)) &
223  kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
224  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
225  if (resoln_scaled) &
226  kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
227  kh_u(i,j) = max(kh_loc, cs%KhTr_min)
228  if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
229  rd_dx=0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points
230  kh_loc = kh_u(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
231  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
232  kh_u(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
233  endif
234  enddo ; enddo
235  !$OMP parallel do default(shared) private(Kh_loc,Rd_dx)
236  do j=js-1,je ; do i=is,ie
237  kh_loc = cs%KhTr
238  if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
239  if (associated(meke%Kh)) &
240  kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
241  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
242  if (resoln_scaled) &
243  kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
244  kh_v(i,j) = max(kh_loc, cs%KhTr_min)
245  if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
246  rd_dx = 0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points
247  kh_loc = kh_v(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
248  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
249  kh_v(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
250  endif
251  enddo ; enddo
252 
253  !$OMP parallel do default(shared)
254  do j=js,je ; do i=is-1,ie
255  khdt_x(i,j) = dt*(kh_u(i,j)*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
256  enddo ; enddo
257  !$OMP parallel do default(shared)
258  do j=js-1,je ; do i=is,ie
259  khdt_y(i,j) = dt*(kh_v(i,j)*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
260  enddo ; enddo
261  elseif (resoln_scaled) then
262  !$OMP parallel do default(shared) private(Res_fn)
263  do j=js,je ; do i=is-1,ie
264  res_fn = 0.5 * (varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
265  kh_u(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
266  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j))) * res_fn
267  enddo ; enddo
268  !$OMP parallel do default(shared) private(Res_fn)
269  do j=js-1,je ; do i=is,ie
270  res_fn = 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
271  kh_v(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
272  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j))) * res_fn
273  enddo ; enddo
274  else ! Use a simple constant diffusivity.
275  if (cs%id_KhTr_u > 0) then
276  !$OMP parallel do default(shared)
277  do j=js,je ; do i=is-1,ie
278  kh_u(i,j) = cs%KhTr
279  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
280  enddo ; enddo
281  else
282  !$OMP parallel do default(shared)
283  do j=js,je ; do i=is-1,ie
284  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
285  enddo ; enddo
286  endif
287  if (cs%id_KhTr_v > 0) then
288  !$OMP parallel do default(shared)
289  do j=js-1,je ; do i=is,ie
290  kh_v(i,j) = cs%KhTr
291  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
292  enddo ; enddo
293  else
294  !$OMP parallel do default(shared)
295  do j=js-1,je ; do i=is,ie
296  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
297  enddo ; enddo
298  endif
299  endif ! VarMix
300 
301  if (cs%max_diff_CFL > 0.0) then
302  if ((cs%id_KhTr_u > 0) .or. (cs%id_KhTr_h > 0)) then
303  !$OMP parallel do default(shared) private(khdt_max)
304  do j=js,je ; do i=is-1,ie
305  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
306  if (khdt_x(i,j) > khdt_max) then
307  khdt_x(i,j) = khdt_max
308  if (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)) > 0.0) &
309  kh_u(i,j) = khdt_x(i,j) / (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
310  endif
311  enddo ; enddo
312  else
313  !$OMP parallel do default(shared) private(khdt_max)
314  do j=js,je ; do i=is-1,ie
315  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
316  khdt_x(i,j) = min(khdt_x(i,j), khdt_max)
317  enddo ; enddo
318  endif
319  if ((cs%id_KhTr_v > 0) .or. (cs%id_KhTr_h > 0)) then
320  !$OMP parallel do default(shared) private(khdt_max)
321  do j=js-1,je ; do i=is,ie
322  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
323  if (khdt_y(i,j) > khdt_max) then
324  khdt_y(i,j) = khdt_max
325  if (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)) > 0.0) &
326  kh_v(i,j) = khdt_y(i,j) / (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
327  endif
328  enddo ; enddo
329  else
330  !$OMP parallel do default(shared) private(khdt_max)
331  do j=js-1,je ; do i=is,ie
332  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
333  khdt_y(i,j) = min(khdt_y(i,j), khdt_max)
334  enddo ; enddo
335  endif
336  endif
337 
338  else ! .not. do_online
339  !$OMP parallel do default(shared)
340  do j=js,je ; do i=is-1,ie
341  khdt_x(i,j) = us%m_to_L**2*read_khdt_x(i,j)
342  enddo ; enddo
343  !$OMP parallel do default(shared)
344  do j=js-1,je ; do i=is,ie
345  khdt_y(i,j) = us%m_to_L**2*read_khdt_y(i,j)
346  enddo ; enddo
347  call pass_vector(khdt_x, khdt_y, g%Domain)
348  endif ! do_online
349 
350  if (cs%check_diffusive_CFL) then
351  if (cs%show_call_tree) call calltree_waypoint("Checking diffusive CFL (tracer_hordiff)")
352  max_cfl = 0.0
353  do j=js,je ; do i=is,ie
354  cfl(i,j) = 2.0*((khdt_x(i-1,j) + khdt_x(i,j)) + &
355  (khdt_y(i,j-1) + khdt_y(i,j))) * g%IareaT(i,j)
356  if (max_cfl < cfl(i,j)) max_cfl = cfl(i,j)
357  enddo ; enddo
358  call cpu_clock_begin(id_clock_sync)
359  call max_across_pes(max_cfl)
360  call cpu_clock_end(id_clock_sync)
361  num_itts = max(1, ceiling(max_cfl - 4.0*epsilon(max_cfl)))
362  i_numitts = 1.0 / (real(num_itts))
363  if (cs%id_CFL > 0) call post_data(cs%id_CFL, cfl, cs%diag, mask=g%mask2dT)
364  elseif (cs%max_diff_CFL > 0.0) then
365  num_itts = max(1, ceiling(cs%max_diff_CFL - 4.0*epsilon(cs%max_diff_CFL)))
366  i_numitts = 1.0 / (real(num_itts))
367  else
368  num_itts = 1 ; i_numitts = 1.0
369  endif
370 
371  do m=1,ntr
372  if (associated(reg%Tr(m)%df_x)) then
373  do k=1,nz ; do j=js,je ; do i=is-1,ie
374  reg%Tr(m)%df_x(i,j,k) = 0.0
375  enddo ; enddo ; enddo
376  endif
377  if (associated(reg%Tr(m)%df_y)) then
378  do k=1,nz ; do j=js-1,je ; do i=is,ie
379  reg%Tr(m)%df_y(i,j,k) = 0.0
380  enddo ; enddo ; enddo
381  endif
382  if (associated(reg%Tr(m)%df2d_x)) then
383  do j=js,je ; do i=is-1,ie ; reg%Tr(m)%df2d_x(i,j) = 0.0 ; enddo ; enddo
384  endif
385  if (associated(reg%Tr(m)%df2d_y)) then
386  do j=js-1,je ; do i=is,ie ; reg%Tr(m)%df2d_y(i,j) = 0.0 ; enddo ; enddo
387  endif
388  enddo
389 
390  if (cs%use_lateral_boundary_diffusion) then
391 
392  if (cs%show_call_tree) call calltree_waypoint("Calling lateral boundary mixing (tracer_hordiff)")
393 
394  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
395 
396  do j=js-1,je ; do i=is,ie
397  coef_y(i,j) = i_numitts * khdt_y(i,j)
398  enddo ; enddo
399  do j=js,je
400  do i=is-1,ie
401  coef_x(i,j) = i_numitts * khdt_x(i,j)
402  enddo
403  enddo
404 
405  do itt=1,num_itts
406  if (cs%show_call_tree) call calltree_waypoint("Calling lateral boundary diffusion (tracer_hordiff)",itt)
407  if (itt>1) then ! Update halos for subsequent iterations
408  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
409  endif
410  call lateral_boundary_diffusion(g, gv, us, h, coef_x, coef_y, i_numitts*dt, reg, &
411  cs%lateral_boundary_diffusion_CSp)
412  enddo ! itt
413  endif
414 
415  if (cs%use_neutral_diffusion) then
416 
417  if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)")
418 
419  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
420  ! We are assuming that neutral surfaces do not evolve (much) as a result of multiple
421  ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs()
422  ! would be inside the itt-loop. -AJA
423 
424  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
425  do j=js-1,je ; do i=is,ie
426  coef_y(i,j) = i_numitts * khdt_y(i,j)
427  enddo ; enddo
428  do j=js,je
429  do i=is-1,ie
430  coef_x(i,j) = i_numitts * khdt_x(i,j)
431  enddo
432  enddo
433 
434  do itt=1,num_itts
435  if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt)
436  if (itt>1) then ! Update halos for subsequent iterations
437  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
438  if (cs%recalc_neutral_surf) then
439  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
440  endif
441  endif
442  call neutral_diffusion(g, gv, h, coef_x, coef_y, i_numitts*dt, reg, us, cs%neutral_diffusion_CSp)
443  enddo ! itt
444 
445  else ! following if not using neutral diffusion, but instead along-surface diffusion
446 
447  if (cs%show_call_tree) call calltree_waypoint("Calculating horizontal diffusion (tracer_hordiff)")
448  do itt=1,num_itts
449  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
450  !$OMP parallel do default(shared) private(scale,Coef_y,Coef_x,Ihdxdy,dTr)
451  do k=1,nz
452  scale = i_numitts
453  if (cs%Diffuse_ML_interior) then
454  if (k<=gv%nkml) then
455  if (cs%ML_KhTr_scale <= 0.0) cycle
456  scale = i_numitts * cs%ML_KhTr_scale
457  endif
458  if ((k>gv%nkml) .and. (k<=gv%nk_rho_varies)) cycle
459  endif
460 
461  do j=js-1,je ; do i=is,ie
462  coef_y(i,j) = ((scale * khdt_y(i,j))*2.0*(h(i,j,k)*h(i,j+1,k))) / &
463  (h(i,j,k)+h(i,j+1,k)+h_neglect)
464  enddo ; enddo
465 
466  do j=js,je
467  do i=is-1,ie
468  coef_x(i,j) = ((scale * khdt_x(i,j))*2.0*(h(i,j,k)*h(i+1,j,k))) / &
469  (h(i,j,k)+h(i+1,j,k)+h_neglect)
470  enddo
471 
472  do i=is,ie
473  ihdxdy(i,j) = g%IareaT(i,j) / (h(i,j,k)+h_neglect)
474  enddo
475  enddo
476 
477  do m=1,ntr
478  do j=js,je ; do i=is,ie
479  dtr(i,j) = ihdxdy(i,j) * &
480  ((coef_x(i-1,j) * (reg%Tr(m)%t(i-1,j,k) - reg%Tr(m)%t(i,j,k)) - &
481  coef_x(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))) + &
482  (coef_y(i,j-1) * (reg%Tr(m)%t(i,j-1,k) - reg%Tr(m)%t(i,j,k)) - &
483  coef_y(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))))
484  enddo ; enddo
485  if (associated(reg%Tr(m)%df_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
486  reg%Tr(m)%df_x(i,j,k) = reg%Tr(m)%df_x(i,j,k) + coef_x(i,j) &
487  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k)) * idt
488  enddo ; enddo ; endif
489  if (associated(reg%Tr(m)%df_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
490  reg%Tr(m)%df_y(i,j,k) = reg%Tr(m)%df_y(i,j,k) + coef_y(i,j) &
491  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k)) * idt
492  enddo ; enddo ; endif
493  if (associated(reg%Tr(m)%df2d_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
494  reg%Tr(m)%df2d_x(i,j) = reg%Tr(m)%df2d_x(i,j) + coef_x(i,j) &
495  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k)) * idt
496  enddo ; enddo ; endif
497  if (associated(reg%Tr(m)%df2d_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
498  reg%Tr(m)%df2d_y(i,j) = reg%Tr(m)%df2d_y(i,j) + coef_y(i,j) &
499  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k)) * idt
500  enddo ; enddo ; endif
501  do j=js,je ; do i=is,ie
502  reg%Tr(m)%t(i,j,k) = reg%Tr(m)%t(i,j,k) + dtr(i,j)
503  enddo ; enddo
504  enddo
505 
506  enddo ! End of k loop.
507 
508  enddo ! End of "while" loop.
509 
510  endif ! endif for CS%use_neutral_diffusion
511  call cpu_clock_end(id_clock_diffuse)
512 
513 
514  if (cs%Diffuse_ML_interior) then
515  if (cs%show_call_tree) call calltree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)")
516  if (cs%debug) call mom_tracer_chksum("Before epipycnal diff ", reg%Tr, ntr, g)
517 
518  call cpu_clock_begin(id_clock_epimix)
519  call tracer_epipycnal_ml_diff(h, dt, reg%Tr, ntr, khdt_x, khdt_y, g, gv, us, &
520  cs, tv, num_itts)
521  call cpu_clock_end(id_clock_epimix)
522  endif
523 
524  if (cs%debug) call mom_tracer_chksum("After tracer diffusion ", reg%Tr, ntr, g)
525 
526  ! post diagnostics for 2d tracer diffusivity
527  if (cs%id_KhTr_u > 0) then
528  do j=js,je ; do i=is-1,ie
529  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
530  enddo ; enddo
531  call post_data(cs%id_KhTr_u, kh_u, cs%diag, mask=g%mask2dCu)
532  endif
533  if (cs%id_KhTr_v > 0) then
534  do j=js-1,je ; do i=is,ie
535  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
536  enddo ; enddo
537  call post_data(cs%id_KhTr_v, kh_v, cs%diag, mask=g%mask2dCv)
538  endif
539  if (cs%id_KhTr_h > 0) then
540  kh_h(:,:) = 0.0
541  do j=js,je ; do i=is-1,ie
542  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
543  enddo ; enddo
544  do j=js-1,je ; do i=is,ie
545  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
546  enddo ; enddo
547  do j=js,je ; do i=is,ie
548  normalize = 1.0 / ((g%mask2dCu(i-1,j)+g%mask2dCu(i,j)) + &
549  (g%mask2dCv(i,j-1)+g%mask2dCv(i,j)) + gv%H_subroundoff)
550  kh_h(i,j) = normalize*g%mask2dT(i,j)*((kh_u(i-1,j)+kh_u(i,j)) + &
551  (kh_v(i,j-1)+kh_v(i,j)))
552  enddo ; enddo
553  call post_data(cs%id_KhTr_h, kh_h, cs%diag, mask=g%mask2dT)
554  endif
555 
556 
557  if (cs%debug) then
558  call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, &
559  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2)
560  if (cs%use_neutral_diffusion) then
561  call uvchksum("After tracer diffusion Coef_[xy]", coef_x, coef_y, &
562  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2)
563  endif
564  endif
565 
566  if (cs%id_khdt_x > 0) call post_data(cs%id_khdt_x, khdt_x, cs%diag)
567  if (cs%id_khdt_y > 0) call post_data(cs%id_khdt_y, khdt_y, cs%diag)
568 
569  if (cs%show_call_tree) call calltree_leave("tracer_hordiff()")
570 

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync, mom_lateral_boundary_diffusion::lateral_boundary_diffusion(), mom_tracer_registry::mom_tracer_chksum(), mom_neutral_diffusion::neutral_diffusion(), mom_neutral_diffusion::neutral_diffusion_calc_coeffs(), and tracer_epipycnal_ml_diff().

Here is the call graph for this function:

Variable Documentation

◆ id_clock_diffuse

integer mom_tracer_hor_diff::id_clock_diffuse
private

CPU time clocks.

Definition at line 97 of file MOM_tracer_hor_diff.F90.

97 integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync

Referenced by tracer_hor_diff_init(), and tracer_hordiff().

◆ id_clock_epimix

integer mom_tracer_hor_diff::id_clock_epimix
private

CPU time clocks.

Definition at line 97 of file MOM_tracer_hor_diff.F90.

Referenced by tracer_hor_diff_init(), and tracer_hordiff().

◆ id_clock_pass

integer mom_tracer_hor_diff::id_clock_pass
private

CPU time clocks.

Definition at line 97 of file MOM_tracer_hor_diff.F90.

Referenced by tracer_epipycnal_ml_diff(), tracer_hor_diff_init(), and tracer_hordiff().

◆ id_clock_sync

integer mom_tracer_hor_diff::id_clock_sync
private

CPU time clocks.

Definition at line 97 of file MOM_tracer_hor_diff.F90.

Referenced by tracer_hor_diff_init(), and tracer_hordiff().