MOM6
Data Types | Functions/Subroutines | Variables
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
 
type  p2di
 
type  tracer_hor_diff_cs
 

Functions/Subroutines

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, CS, CSnd)
 Initialize lateral tracer diffusion module. More...
 
subroutine, public tracer_hor_diff_end (CS)
 

Variables

integer id_clock_diffuse
 
integer id_clock_epimix
 
integer id_clock_pass
 
integer id_clock_sync
 

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]gocean grid structure
[in]gvocean vertical grid structure
[in]hlayer thickness (m or kg m-2)
[in]dttime step
[in,out]trtracer array
[in]ntrnumber of tracers
[in]khdt_epi_xneeds a comment
[in]khdt_epi_yneeds a comment
[in,out]csmodule control structure
[in]tvthermodynamic structure
[in]num_ittsnumber of iterations (usually=1)

Definition at line 497 of file MOM_tracer_hor_diff.F90.

References id_clock_pass.

Referenced by tracer_hordiff().

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

Definition at line 1448 of file MOM_tracer_hor_diff.F90.

1448  type(tracer_hor_diff_cs), pointer :: cs
1449 
1450  call neutral_diffusion_end(cs%neutral_diffusion_CSp)
1451  if (associated(cs)) deallocate(cs)
1452 

◆ 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(tracer_hor_diff_cs), pointer  CS,
type(neutral_diffusion_cs), pointer  CSnd 
)

Initialize lateral tracer diffusion module.

Parameters
[in]timecurrent model time
[in]gocean grid structure
[in,out]diagdiagnostic control
[in]param_fileparameter file
cshorz diffusion control structure
csndpointer to neutral diffusion CS

Definition at line 1341 of file MOM_tracer_hor_diff.F90.

References id_clock_diffuse, id_clock_epimix, id_clock_pass, and id_clock_sync.

1341  type(time_type), target, intent(in) :: time !< current model time
1342  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1343  type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control
1344  type(param_file_type), intent(in) :: param_file !< parameter file
1345  type(tracer_hor_diff_cs), pointer :: cs !< horz diffusion control structure
1346  type(neutral_diffusion_cs), pointer :: csnd !< pointer to neutral diffusion CS
1347 
1348 ! This include declares and sets the variable "version".
1349 #include "version_variable.h"
1350  character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name.
1351  character(len=256) :: mesg ! Message for error messages.
1352 
1353  if (associated(cs)) then
1354  call mom_error(warning, "tracer_hor_diff_init called with associated control structure.")
1355  return
1356  endif
1357  allocate(cs)
1358 
1359  cs%diag => diag
1360  cs%show_call_tree = calltree_showquery()
1361 
1362  ! Read all relevant parameters and write them to the model log.
1363  call log_version(param_file, mdl, version, "")
1364  call get_param(param_file, mdl, "KHTR", cs%KhTr, &
1365  "The background along-isopycnal tracer diffusivity.", &
1366  units="m2 s-1", default=0.0)
1367  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", cs%KhTr_Slope_Cff, &
1368  "The scaling coefficient for along-isopycnal tracer \n"//&
1369  "diffusivity using a shear-based (Visbeck-like) \n"//&
1370  "parameterization. A non-zero value enables this param.", &
1371  units="nondim", default=0.0)
1372  call get_param(param_file, mdl, "KHTR_MIN", cs%KhTr_Min, &
1373  "The minimum along-isopycnal tracer diffusivity.", &
1374  units="m2 s-1", default=0.0)
1375  call get_param(param_file, mdl, "KHTR_MAX", cs%KhTr_Max, &
1376  "The maximum along-isopycnal tracer diffusivity.", &
1377  units="m2 s-1", default=0.0)
1378  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", cs%KhTr_passivity_coeff, &
1379  "The coefficient that scales deformation radius over \n"//&
1380  "grid-spacing in passivity, where passiviity is the ratio \n"//&
1381  "between along isopycnal mxiing of tracers to thickness mixing. \n"//&
1382  "A non-zero value enables this parameterization.", &
1383  units="nondim", default=0.0)
1384  call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", cs%KhTr_passivity_min, &
1385  "The minimum passivity which is the ratio between \n"//&
1386  "along isopycnal mxiing of tracers to thickness mixing. \n", &
1387  units="nondim", default=0.5)
1388  call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
1389  desc="The (baroclinic) dynamics time step.", units="s")
1390 
1391 
1392  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
1393  "If true, enable epipycnal mixing between the surface \n"//&
1394  "boundary layer and the interior.", default=.false.)
1395  call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", cs%check_diffusive_CFL, &
1396  "If true, use enough iterations the diffusion to ensure \n"//&
1397  "that the diffusive equivalent of the CFL limit is not \n"//&
1398  "violated. If false, always use 1 iteration.", default=.false.)
1399  cs%ML_KhTR_scale = 1.0
1400  if (cs%Diffuse_ML_interior) then
1401  call get_param(param_file, mdl, "ML_KHTR_SCALE", cs%ML_KhTR_scale, &
1402  "With Diffuse_ML_interior, the ratio of the truly \n"//&
1403  "horizontal diffusivity in the mixed layer to the \n"//&
1404  "epipycnal diffusivity. The valid range is 0 to 1.", &
1405  units="nondim", default=1.0)
1406  endif
1407 
1408  cs%use_neutral_diffusion = neutral_diffusion_init(time, g, param_file, diag, cs%neutral_diffusion_CSp)
1409  csnd => cs%neutral_diffusion_CSp
1410  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1411  "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1412 
1413  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1414 
1415  id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=clock_module)
1416  id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=clock_module)
1417  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1418  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1419 
1420  cs%id_KhTr_u = -1
1421  cs%id_KhTr_v = -1
1422  cs%id_KhTr_h = -1
1423  cs%id_CFL = -1
1424 
1425  cs%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCu1, time, &
1426  'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'meter2 second-1')
1427  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
1428  'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'meter2 second-1')
1429  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time,&
1430  'Epipycnal tracer diffusivity at tracer cell center', 'meter2 second-1', &
1431  cmor_field_name='diftrelo', cmor_units='m2 sec-1', &
1432  cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', &
1433  cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')
1434 
1435  cs%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, time, &
1436  'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'meter2')
1437  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
1438  'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'meter2')
1439  if (cs%check_diffusive_CFL) then
1440  cs%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, time,&
1441  'Grid CFL number for lateral/neutral tracer diffusion', 'dimensionless')
1442  endif
1443 
1444 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19

◆ 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), intent(inout)  Reg,
type(thermo_var_ptrs), intent(in)  tv,
logical, 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 (m or kg m-2)
[in]dttime step (seconds)
mekeMEKE type
varmixVariable mixing type
[in]gvocean vertical grid structure
csmodule control structure
[in,out]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.

Definition at line 88 of file MOM_tracer_hor_diff.F90.

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_tracer_registry::mom_tracer_chksum(), mom_neutral_diffusion::neutral_diffusion(), mom_neutral_diffusion::neutral_diffusion_calc_coeffs(), and tracer_epipycnal_ml_diff().

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

Variable Documentation

◆ id_clock_diffuse

integer mom_tracer_hor_diff::id_clock_diffuse
private

Definition at line 79 of file MOM_tracer_hor_diff.F90.

Referenced by tracer_hor_diff_init(), and tracer_hordiff().

79 integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync

◆ id_clock_epimix

integer mom_tracer_hor_diff::id_clock_epimix
private

Definition at line 79 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

Definition at line 79 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

Definition at line 79 of file MOM_tracer_hor_diff.F90.

Referenced by tracer_hor_diff_init(), and tracer_hordiff().