Skip to content

Commit

Permalink
Merge pull request #68 from ESMG/tracer_advection_OBC
Browse files Browse the repository at this point in the history 
*Tracer advection obc
  • Loading branch information
@kshedstrom
kshedstrom authored Aug 21, 2019
2 parents 6b070fd + 6860cec commit 4af643c
Showing 1 changed file with 144 additions and 91 deletions.
235 changes: 144 additions & 91 deletions src/tracer/MOM_tracer_advect.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -353,6 +353,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
slope_x ! The concentration slope per grid point [conc].
real, dimension(SZIB_(G),ntr) :: &
flux_x ! The tracer flux across a boundary [H m2 conc ~> m3 conc or kg conc].
real, dimension(SZI_(G),ntr) :: &
T_tmp ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].
real :: maxslope ! The maximum concentration slope per grid point
! consistent with monotonicity [conc].
real :: hup, hlos ! hup is the upwind volume, hlos is the
Expand Down Expand Up @@ -420,6 +422,72 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
enddo ; enddo
endif ! usePLMslope

! make a copy of the tracers in case values need to be overridden for OBCs
do m = 1,ntr
do i=is-stencil,ie+stencil
T_tmp(i,m) = Tr(m)%t(i,j,k)
enddo
enddo
! loop through open boundaries and recalculate flux terms
if (associated(OBC)) then ; if (OBC%OBC_pe) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. associated(segment%tr_Reg)) cycle
if (segment%is_E_or_W) then
if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
I = segment%HI%IsdB

ishift=0 ! ishift+I corresponds to the nearest interior tracer cell index
idir=1 ! idir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_W) then
ishift=1
idir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
uhh(I)=uhr(I,j,k)
u_L_in=max(idir*uhh(I)*segment%Tr_InvLscale3_in,0.)
u_L_out=min(idir*uhh(I)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(u_L_in-u_L_out)
segment%tr_Reg%Tr(m)%tres(I,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
dt*(u_L_in*Tr(m)%t(I+ishift,j,k) - &
u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k)))
endif
enddo

do m = 1,ntr ! replace tracers with OBC values
if (associated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_W) then
T_tmp(i,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
else
T_tmp(I+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
endif
else
if (segment%direction == OBC_DIRECTION_W) then
T_tmp(i,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
else
T_tmp(I+1,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
endif
endif
enddo
do m = 1,ntr ! Apply update tracer values for slope calculation
do i=segment%HI%IsdB-1,segment%HI%IsdB+1
Tp = T_tmp(i+1,m) ; Tc = T_tmp(i,m) ; Tm = T_tmp(i-1,m)
dMx = max( Tp, Tc, Tm ) - Tc
dMn= Tc - min( Tp, Tc, Tm )
slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * &
sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm )
enddo
enddo

endif
endif
enddo
endif; endif


! Calculate the i-direction fluxes of each tracer, using as much
! the minimum of the remaining mass flux (uhr) and the half the mass
! in the cell plus whatever part of its half of the mass flux that
Expand Down Expand Up @@ -466,7 +534,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
endif

! Implementation of PPM-H3
Tp = Tr(m)%t(i_up+1,j,k) ; Tc = Tr(m)%t(i_up,j,k) ; Tm = Tr(m)%t(i_up-1,j,k)
Tp = T_tmp(i_up+1,m) ; Tc = T_tmp(i_up,m) ; Tm = T_tmp(i_up-1,m)

if (useHuynh) then
aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate
Expand Down Expand Up @@ -508,7 +576,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
!aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
!flux_x(I,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) )
! Alternative implementation of PLM
Tc = Tr(m)%t(i,j,k)
Tc = T_tmp(i,m)
flux_x(I,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) )
! Original implementation of PLM
!flux_x(I,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I))
Expand All @@ -521,7 +589,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
!aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
!flux_x(I,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) )
! Alternative implementation of PLM
Tc = Tr(m)%t(i+1,j,k)
Tc = T_tmp(i+1,m)
flux_x(I,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) )
! Original implementation of PLM
!flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I))
Expand All @@ -531,10 +599,9 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
endif ! usePPM

if (associated(OBC)) then ; if (OBC%OBC_pe) then
if (OBC%specified_u_BCs_exist_globally) then
if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. segment%specified) cycle
if (.not. associated(segment%tr_Reg)) cycle
if (segment%is_E_or_W) then
if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
Expand All @@ -554,48 +621,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
endif
endif
enddo
endif
endif

if (OBC%open_u_BCs_exist_globally) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
I = segment%HI%IsdB
if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
if (segment%specified) cycle
if (.not. associated(segment%tr_Reg)) cycle
ishift=0 ! ishift+I corresponds to the nearest interior tracer cell index
idir=1 ! idir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_W) then
ishift=1
idir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
uhh(I)=uhr(I,j,k)
u_L_in=max(idir*uhh(I)*segment%Tr_InvLscale3_in,0.)
u_L_out=min(idir*uhh(I)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(u_L_in-u_L_out)
segment%tr_Reg%Tr(m)%tres(I,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
dt*(u_L_in*Tr(m)%t(I+ishift,j,k) - &
u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k)))
endif
enddo

! Tracer fluxes are set to prescribed values only for inflows from masked areas.
if ((uhr(I,j,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. &
(uhr(I,j,k) < 0.0) .and. (G%mask2dT(i+1,j) < 0.5)) then
uhh(I) = uhr(I,j,k)
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k)
else; flux_x(I,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif
enddo
endif
endif
enddo
endif
endif ; endif

! Calculate new tracer concentration in each cell after accounting
Expand Down Expand Up @@ -680,7 +707,10 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
slope_y ! The concentration slope per grid point [conc].
real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
flux_y ! The tracer flux across a boundary [H m2 conc ~> m3 conc or kg conc].
flux_y ! The tracer flux across a boundary [H m2 conc ~> m3 conc or kg conc].
real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
T_tmp ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].

real :: maxslope ! The maximum concentration slope per grid point
! consistent with monotonicity [conc].
real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
Expand Down Expand Up @@ -757,6 +787,72 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
enddo ; enddo ; endif ; enddo ! End of i-, m-, & j- loops.
endif ! usePLMslope


! make a copy of the tracers in case values need to be overridden for OBCs

do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr ; do i=is,ie
T_tmp(i,m,j) = Tr(m)%t(i,j,k)
enddo ; enddo ; endif ; enddo

! loop through open boundaries and recalculate flux terms
if (associated(OBC)) then ; if (OBC%OBC_pe) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. associated(segment%tr_Reg)) cycle
do i=is,ie
if (segment%is_N_or_S) then
if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
J = segment%HI%JsdB
jshift=0 ! jshift+J corresponds to the nearest interior tracer cell index
jdir=1 ! jdir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_S) then
jshift=1
jdir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
vhh(i,J)=vhr(i,J,k)
v_L_in=max(jdir*vhh(i,J)*segment%Tr_InvLscale3_in,0.)
v_L_out=min(jdir*vhh(i,J)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(v_L_in-v_L_out)
segment%tr_Reg%Tr(m)%tres(i,J,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
dt*(v_L_in*Tr(m)%t(i,J+jshift,k) - &
v_L_out*segment%tr_Reg%Tr(m)%t(i,J,k)))
endif
enddo

do m = 1,ntr ! replace tracers with OBC values
if (associated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_S) then
T_tmp(i,m,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
else
T_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
endif
else
if (segment%direction == OBC_DIRECTION_S) then
T_tmp(i,m,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
else
T_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
endif
endif
enddo
do m = 1,ntr ! Apply update tracer values for slope calculation
do j=segment%HI%JsdB-1,segment%HI%JsdB+1
Tp = T_tmp(i,m,j+1) ; Tc = T_tmp(i,m,j) ; Tm = T_tmp(i,m,j-1)
dMx = max( Tp, Tc, Tm ) - Tc
dMn= Tc - min( Tp, Tc, Tm )
slope_y(i,m,j) = G%mask2dCv(i,J)*G%mask2dCv(i,J-1) * &
sign( min(0.5*abs(Tp-Tm), 2.0*dMx, 2.0*dMn), Tp-Tm )
enddo
enddo
endif
endif ! is_N_S
enddo ! i-loop
enddo ! segment loop
endif; endif

! Calculate the j-direction fluxes of each tracer, using as much
! the minimum of the remaining mass flux (vhr) and the half the mass
! in the cell plus whatever part of its half of the mass flux that
Expand Down Expand Up @@ -869,7 +965,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
endif ! usePPM

if (associated(OBC)) then ; if (OBC%OBC_pe) then
if (OBC%specified_v_BCs_exist_globally) then
if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. segment%specified) cycle
Expand All @@ -893,49 +989,6 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
endif
enddo
endif


if (OBC%open_v_BCs_exist_globally) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (segment%specified) cycle
if (.not. associated(segment%tr_Reg)) cycle
if (segment%is_N_or_S .and. &
(J >= segment%HI%JsdB .and. J<= segment%HI%JedB)) then
jshift=0
jdir=1
if (segment%direction == OBC_DIRECTION_S) then
jshift=1
jdir=-1
endif
do i=segment%HI%isd,segment%HI%ied
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
vhh(i,J)=vhr(i,J,k)
v_L_in=max(jdir*vhh(i,J)*segment%Tr_InvLscale3_in,0.)
v_L_out=min(jdir*vhh(i,J)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(v_L_in-v_L_out)
segment%tr_Reg%Tr(m)%tres(i,J,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
dt*v_L_in*Tr(m)%t(i,j+jshift,k) - &
dt*v_L_out*segment%tr_Reg%Tr(m)%t(i,j,k))
endif
enddo
! Tracer fluxes are set to prescribed values only for inflows from masked areas.
if ((vhr(i,J,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. &
(vhr(i,J,k) < 0.0) .and. (G%mask2dT(i,j+1) < 0.5)) then
vhh(i,J) = vhr(i,J,k)
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%t)) then
flux_y(i,m,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%tres(i,J,k)
else ; flux_y(i,m,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
enddo
endif
enddo
endif
enddo
endif
endif; endif

else ! not domore_v.
Expand Down

0 comments on commit 4af643c

Please sign in to comment.