From f0ed7f8760b36b578321facdcfed5b438dec0562 Mon Sep 17 00:00:00 2001 From: matthew harrison Date: Tue, 20 Aug 2019 13:36:02 -0400 Subject: [PATCH] changes needed for reproducibility across restarts using OBCs --- src/tracer/MOM_tracer_advect.F90 | 278 +++++++++++++++++++++---------- 1 file changed, 187 insertions(+), 91 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 201f8aeb6f..28312f1dca 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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,j,m) = segment%tr_Reg%Tr(m)%tres(i,j,k) + else + T_tmp(i,j+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k) + endif + else + if (segment%direction == OBC_DIRECTION_S) then + T_tmp(i,j,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc + else + T_tmp(i,j+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 j=segment%HI%JsdB-1,segment%HI%JsdB+1 + Tp = T_tmp(i,j+1,m) ; Tc = T_tmp(i,j,m) ; Tm = T_tmp(i,j-1,m) + 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 @@ -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 @@ -893,51 +989,51 @@ 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 + ! 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. do i=is,ie ; vhh(i,J) = 0.0 ; enddo do m=1,ntr ; do i=is,ie ; flux_y(i,m,J) = 0.0 ; enddo ; enddo