Skip to content

Commit

Permalink
Avoid copying the tracer registry in advect_tracer
Browse files Browse the repository at this point in the history 
  This change avoids an unnecessary copy of the tracer registry in
advect_tracer().  This should be more efficient, and it should facilitate the
development of improved tracer budget diagnostic capabilities.  It also adds
comments describing a few internal variables and their units.  All answers and
output are bitwise identical.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Mar 16, 2022
1 parent 5a8590b commit 14a28f4
Showing 1 changed file with 32 additions and 36 deletions.
68 changes: 32 additions & 36 deletions src/tracer/MOM_tracer_advect.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,6 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
optional, intent(out) :: vhr_out !< Remaining accumulated volume or mass fluxes
!! through the meridional faces [H L2 ~> m3 or kg]

type(tracer_type) :: Tr(MAX_FIELDS_) ! The array of registered tracers
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
hprev ! cell volume at the end of previous tracer change [H L2 ~> m3 or kg]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
Expand Down Expand Up @@ -127,7 +126,6 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3

ntr = Reg%ntr
do m=1,ntr ; Tr(m) = Reg%Tr(m) ; enddo
Idt = 1.0 / dt

max_iter = 2*INT(CEILING(dt/CS%dt)) + 1
Expand All @@ -138,7 +136,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
call create_group_pass(CS%pass_uhr_vhr_t_hprev, uhr, vhr, G%Domain)
call create_group_pass(CS%pass_uhr_vhr_t_hprev, hprev, G%Domain)
do m=1,ntr
call create_group_pass(CS%pass_uhr_vhr_t_hprev, Tr(m)%t, G%Domain)
call create_group_pass(CS%pass_uhr_vhr_t_hprev, Reg%Tr(m)%t, G%Domain)
enddo
call cpu_clock_end(id_clock_pass)

Expand Down Expand Up @@ -188,27 +186,11 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
! initialize diagnostic fluxes and tendencies
!$OMP do
do m=1,ntr
if (associated(Tr(m)%ad_x)) then
do k=1,nz ; do j=jsd,jed ; do i=isd,ied
Tr(m)%ad_x(I,j,k) = 0.0
enddo ; enddo ; enddo
endif
if (associated(Tr(m)%ad_y)) then
do k=1,nz ; do J=jsd,jed ; do i=isd,ied
Tr(m)%ad_y(i,J,k) = 0.0
enddo ; enddo ; enddo
endif
if (associated(Tr(m)%advection_xy)) then
do k=1,nz ; do j=jsd,jed ; do i=isd,ied
Tr(m)%advection_xy(i,j,k) = 0.0
enddo ; enddo ; enddo
endif
if (associated(Tr(m)%ad2d_x)) then
do j=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_x(I,j) = 0.0 ; enddo ; enddo
endif
if (associated(Tr(m)%ad2d_y)) then
do J=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_y(i,J) = 0.0 ; enddo ; enddo
endif
if (associated(Reg%Tr(m)%ad_x)) Reg%Tr(m)%ad_x(:,:,:) = 0.0
if (associated(Reg%Tr(m)%ad_y)) Reg%Tr(m)%ad_y(:,:,:) = 0.0
if (associated(Reg%Tr(m)%advection_xy)) Reg%Tr(m)%advection_xy(:,:,:) = 0.0
if (associated(Reg%Tr(m)%ad2d_x)) Reg%Tr(m)%ad2d_x(:,:) = 0.0
if (associated(Reg%Tr(m)%ad2d_y)) Reg%Tr(m)%ad2d_y(:,:) = 0.0
enddo
!$OMP end parallel

Expand Down Expand Up @@ -265,14 +247,14 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
!$OMP do ordered
do k=1,nz ; if (domore_k(k) > 0) then
! First, advect zonally.
call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
isv, iev, jsv-stencil, jev+stencil, k, G, GV, US, CS%usePPM, CS%useHuynh)
endif ; enddo

!$OMP do ordered
do k=1,nz ; if (domore_k(k) > 0) then
! Next, advect meridionally.
call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh)

! Update domore_k(k) for the next iteration
Expand All @@ -287,14 +269,14 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
!$OMP do ordered
do k=1,nz ; if (domore_k(k) > 0) then
! First, advect meridionally.
call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
isv-stencil, iev+stencil, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh)
endif ; enddo

!$OMP do ordered
do k=1,nz ; if (domore_k(k) > 0) then
! Next, advect zonally.
call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh)

! Update domore_k(k) for the next iteration
Expand Down Expand Up @@ -390,13 +372,20 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2].
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc]
real :: dMx ! Difference between the maximum of the surrounding cell concentrations and
! the value in the cell whose reconstruction is being found [conc]
real :: dMn ! Difference between the tracer concentration in the cell whose reconstruction
! is being found and the minimum of the surrounding values [conc]
real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc]
real :: dA ! Difference between the reconstruction tracer edge values [conc]
real :: mA ! Average of the reconstruction tracer edge values [conc]
real :: a6 ! Curvature of the reconstruction tracer values [conc]
logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points.
logical :: do_any_i
logical :: usePLMslope
integer :: i, j, m, n, i_up, stencil
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
real :: fac1,u_L_in,u_L_out ! terms used for time-stepping OBC reservoirs
type(OBC_segment_type), pointer :: segment=>NULL()
logical :: usePLMslope
logical, dimension(SZJ_(G),SZK_(GV)) :: domore_u_initial

! keep a local copy of the initial values of domore_u, which is to be used when computing ad2d_x
Expand Down Expand Up @@ -547,7 +536,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &

dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
Expand Down Expand Up @@ -754,14 +743,21 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2].
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc]
real :: dMx ! Difference between the maximum of the surrounding cell concentrations and
! the value in the cell whose reconstruction is being found [conc]
real :: dMn ! Difference between the tracer average in the cell whose reconstruction
! is being found and the minimum of the surrounding values [conc]
real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc]
real :: dA ! Difference between the reconstruction tracer edge values [conc]
real :: mA ! Average of the reconstruction tracer edge values [conc]
real :: a6 ! Curvature of the reconstruction tracer values [conc]
logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
logical :: do_i(SZIB_(G), SZJ_(G)) ! If true, work on given points.
logical :: do_any_i
logical :: usePLMslope
integer :: i, j, j2, m, n, j_up, stencil
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
real :: fac1,v_L_in,v_L_out ! terms used for time-stepping OBC reservoirs
type(OBC_segment_type), pointer :: segment=>NULL()
logical :: usePLMslope

usePLMslope = .not. (usePPM .and. useHuynh)
! stencil for calculating slope values
Expand Down Expand Up @@ -919,7 +915,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &

dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCv(i,J_up)*G%mask2dCv(i,J_up-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
Expand Down

0 comments on commit 14a28f4

Please sign in to comment.