diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index e2c669fcc7..79fd4f3cbf 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -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)) :: & @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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