subroutine advection(N, Y, dt, h, hint, v, Bcup, Bcdown, & & Yup, Ydown, VFint, VF, Aint, A, method, mode, dY, cu, it) c----------------------------------------------------------------------------------------- c based on the advection scheme in the GOTM model, code from 2006-11-06 c fluxes are defined on the interfaces, in an upstream-biased way. c slope-delimeters are applied to obtain monotonic and positive schemes c also in the presence of large gradients. c there are 5 different slope delimeters; first-order upstream, c 3rd order upstream-biased polynomial, 3rd order with superbee delimiter, c 3rd order with MUSCL limiter c as described in Pietrzak 98 c Karline: made changes to make it work for negative ww... c added volume fraction, surface area; these two properties will generally = 1. cc----------------------------------------------------------------------------------------- IMPLICIT NONE c number of vertical layers, time step INTEGER :: N DOUBLE PRECISION :: dt c layer thickness (m), distance from mid to mid of each layer DOUBLE PRECISION :: h(N), hint(0:N) c advection speed in the direction of the axis DOUBLE PRECISION :: v(0:N), ww(0:N) c volume fraction and surface at interface and in middle of layers DOUBLE PRECISION :: VFint(0:N), Aint(0:N), VF(N), A(N) c type of upper and lower Boundary Condition (only 1 and 2 used in R) INTEGER :: Bcdown, Bcup integer, parameter :: Flux = 1 integer, parameter :: Value = 2 integer, parameter :: OneSided = 3 integer, parameter :: zeroDivergence = 4 c value of upper and lower bnd conc DOUBLE PRECISION :: Ydown, Yup c type of advection scheme, slope delimeters INTEGER :: method integer, parameter :: UPSTREAM =5 integer, parameter :: P2 =4 integer, parameter :: P2_PDM =3 integer, parameter :: Superbee =2 integer, parameter :: MUSCL =1 c advection mode: 0= non-conservative (e.g.water flow), 1= conservative, e.g. sinking INTEGER :: mode c concentration to be transported DOUBLE PRECISION :: Y(N) c rate of change due to advection DOUBLE PRECISION :: dY(N) DOUBLE PRECISION :: one6th=1.0d0/6.0d0 c maximal number of iterations INTEGER, parameter :: itmax=100 c LOCAL VARIABLES: integer :: i,k,it DOUBLE PRECISION :: x,r,Phi,limit=0.d0 DOUBLE PRECISION :: Yu,Yc,Yd DOUBLE PRECISION :: c,cmax DOUBLE PRECISION :: cu(0:N) c----------------------------------------------------------------------------------------- c initialize upstream interface fluxes with zero cu = 0.d0 c initialize maximum Courant number cmax = 0.d0 c copy of current value of state variables do k =1,N dy(k) = y(k) enddo c compute maximum Courant number; estimate nr of iterations do k=0,N ww(k) = v(k)*VFint(k) ! convert to per fraction (in case VFint !=1) c = dabs(ww(k))*dt/hint(k) if (c.gt.cmax) cmax=c enddo it=min(itmax,int(cmax)+1) c (time) splitting loop do i=1,it c spatial loop - karline : changed into 1:N do k=1,N c positive speed if (ww(k) .gt. 0.d0) then c=ww(k)/dble(it)*dt/hint(k-1) ! courant number if (k .gt. 1) then Yu=Y(k-1) ! upstream value else Yu=Y(k) end if Yc=Y(k ) ! central value if (k .lt. N) then Yd=Y(k+1) ! downstream value else Yd=Y(k) endif if (abs(Yd-Yc) .gt. 1e-10) then ! slope ratio r=(Yc-Yu)/(Yd-Yc) else r=(Yc-Yu)*1.e10 end if c negative speed else c=-ww(k)/dble(it)*dt/hint(k) ! courant number if (k .lt. N) then Yu=Y(k+1) ! upstream value else Yu=Y(k) end if Yc=Y(k) ! central value if (k .gt. 1) then Yd=Y(k-1) ! downstream value else Yd=Y(k) endif if (abs(Yc-Yd) .gt. 1e-10) then ! slope ratio r=(Yu-Yc)/(Yc-Yd) else r=(Yu-Yc)*1.e10 end if end if c limit the flux according to different suggestions select case (method) case (UPSTREAM) limit=0.d0 case ((P2),(P2_PDM)) c the flux-factor phi x = one6th*(1.-2.0*c) Phi = (0.5+x)+(0.5-x)*r if (method.eq.P2) then limit=Phi else limit=max(0.d0,min(Phi,2./(1.d0-c),2.*r/(c+1.e-10))) end if case (Superbee) limit=max(0.d0, min(1.d0, 2.0*r), min(r,2.*1.d0) ) case (MUSCL) limit=max(0.d0,min(2.*1.d0,2.0*r,0.5*(1.d0+r))) case default c call rerror( 'unkown advection method') ! should not happen end select c compute the limited flux KARLINE: changed for negative ww(k) ! if (ww(k) .gt. 0.) THEN cu(k) = ww(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc)) else cu(k-1)= ww(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc)) endif end do c downstream boundary conditions select case (Bcdown) case (flux) cu(N) = Ydown ! flux OUT of the domain is positive case (value) if (ww(N).lt. 0.d0) then cu(N) = ww(N)*Ydown else cu(N) = ww(N)*Y(N) end if case (oneSided) if (ww(N).lt. 0.d0) then cu(N) = ww(N)*Y(N) else cu(N) = 0.d0 end if case (zeroDivergence) cu(N) = cu(N-1) case default call rwarn('unkown downstream boundary condition type') end select c upstream boundary conditions select case (Bcup) case (flux) cu(0) = Yup ! flux into the domain is positive case (value) if(ww(0) .gt. 0.d0) then cu(0) = ww(0)*Yup else cu(0) = ww(0)*Y(1) end if case (oneSided) if(ww(0) .ge. 0.d0) then cu(0) = ww(0)*Y(1) else cu(0) = 0.d0 end if case (zeroDivergence) cu(0) = cu(1) case default call rwarn('unkown upstream boundary condition type') end select c the advection step if (mode.eq.0) then ! non-conservative - KARLINE CHECK... do k=1,N Y(k)=Y(k)-1.d0/dble(it)*dt*((Aint(k)*cu(k)- & & Aint(k-1)*cu(k-1))/ h(k)/A(k)/VF(k) & & -Y(k)*(ww(k)-ww(k-1))/h(k)/A(k)/VF(k)) enddo else ! conservative do k=1,N Y(k)=Y(k)-1.d0/dble(it)*dt*((Aint(k)*cu(k)- & & Aint(k-1)*cu(k-1))/h(k)/A(k)/VF(k)) enddo end if end do ! end of the iteration loop c rate of change due to advection do k =1,N dy(k) = (y(k)-dy(k))/dt enddo c Still to do: integrate fluxes in time (now cu = cu of last step ) c flux = 0. ! at start c do k=1,N c flux(k)=flux(k)+1.d0/dble(it)*dt*cu(k) return end subroutine advection !----------------------------------------------------------------------- ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org ! ... extensively modified by Karline Soetaert !----------------------------------------------------------------------- subroutine advectvol(N, Y, dt, V, Vint, flow, Bcup, Bcdown, & & Yup, Ydown, method,mode,dY, cu, it) c----------------------------------------------------------------------------------------- c Similar as above, but for volumetric transport: c use flow = v*A rather than v c volume = h*A rather than h c hint -> c(volume(1),volume) ... for now ... not so important c----------------------------------------------------------------------------------------- IMPLICIT NONE c number of vertical layers, time step INTEGER :: N DOUBLE PRECISION :: dt c layer thickness (m), distance from mid to mid DOUBLE PRECISION :: V(N), Vint(0:N) c vertical advection speed DOUBLE PRECISION :: flow(0:N) c type of upper and lower Boundary Condition (1 and 2 used) INTEGER :: Bcdown, Bcup integer, parameter :: Flux =1 integer, parameter :: Value =2 integer, parameter :: OneSided =3 ! not used integer, parameter :: zeroDivergence =4 ! not used c value of upper and lower bnd conc DOUBLE PRECISION :: Ydown, Yup c type of advection scheme, slope delimeters INTEGER :: method integer, parameter :: UPSTREAM =5 integer, parameter :: P2 =4 integer, parameter :: P2_PDM =3 integer, parameter :: Superbee =2 integer, parameter :: MUSCL =1 c advection mode 0: non-conservative (e.g. water flow), 1: conservative, eg.g. sinking INTEGER :: mode c concentration to be transported DOUBLE PRECISION :: Y(N) c rate of change due to advection DOUBLE PRECISION :: dY(N) DOUBLE PRECISION :: one6th=1.0d0/6.0d0 INTEGER, parameter :: itmax=100 c LOCAL VARIABLES: integer :: i,k,it DOUBLE PRECISION :: x,r,Phi,limit=0.d0 DOUBLE PRECISION :: Yu,Yc,Yd DOUBLE PRECISION :: c,cmax DOUBLE PRECISION :: cu(0:N) c----------------------------------------------------------------------------------------- c initialize interface fluxes with zero cu = 0.d0 c initialize maximum Courant number cmax = 0.d0 c copy of current value of state variables do k =1,N dy(k) = y(k) enddo c compute maximum Courant number; estimate nr of iterations do k=0,N c=dabs(flow(k))*dt/Vint(k) if (c.gt.cmax) cmax=c enddo it=min(itmax,int(cmax)+1) c (time) splitting loop do i=1,it c spatial loop - karline : changed into 1:N do k=1,N c positive speed if (flow(k) .gt. 0.d0) then c=flow(k)/dble(it)*dt/Vint(k-1) ! courant number if (k .gt. 1) then Yu=Y(k-1) ! upstream value else Yu=Y(k) end if Yc=Y(k ) ! central value if (k .lt. N) then Yd=Y(k+1) ! downstream value else Yd=Y(k) endif if (abs(Yd-Yc) .gt. 1e-10) then ! slope ratio r=(Yc-Yu)/(Yd-Yc) else r=(Yc-Yu)*1.e10 end if c negative speed else c=-flow(k)/dble(it)*dt/Vint(k) ! courant number if (k .lt. N) then Yu=Y(k+1) ! upstream value else Yu=Y(k) end if Yc=Y(k) ! central value if (k .gt. 1) then Yd=Y(k-1) ! downstream value else Yd=Y(k) endif if (dabs(Yc-Yd) .gt. 1e-10) then ! slope ratio r=(Yu-Yc)/(Yc-Yd) else r=(Yu-Yc)*1.e10 ! CHECK THIS: not -??? end if end if c limit the flux according to different suggestions, phi = flux-factor select case (method) case (UPSTREAM) limit=0.d0 case ((P2),(P2_PDM)) c - for quickest x = one6th*(1.d0-2.d0*c) Phi = (0.5d0+x)+(0.5d0-x)*r if (method.eq.P2) then limit=Phi else limit=max(0.d0,min(Phi,2./(1.-c),2.*r/(c+1.e-10))) end if case (Superbee) limit=max(0.d0, min(1.d0, 2.0*r), min(r,2.*1.d0) ) case (MUSCL) limit=max(0.d0,min(2.*1.d0,2.0*r,0.5*(1.0+r))) case default c call rerror( 'unkown advection method') end select c compute the limited flux KARLINE: changed for negative flow(k) ! if (flow(k) .gt. 0.) THEN cu(k) = flow(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc)) else cu(k-1)= flow(k) *(Yc+0.5d0*limit*(1.-c)*(Yd-Yc)) endif end do c downstream boundary conditions select case (Bcdown) case (flux) cu(N) = Ydown ! flux OUT of the domain is positive case (value) if (flow(N).lt. 0.d0) then cu(N) = flow(N)*Ydown else cu(N) = flow(N)*Y(N) end if case (oneSided) if (flow(N).lt. 0.d0) then cu(N) = flow(N)*Y(N) else cu(N) = 0.d0 end if case (zeroDivergence) cu(N) = cu(N-1) case default call rwarn('unkown downstream boundary condition type') end select c upstream boundary conditions select case (Bcup) case (flux) cu(0) = Yup ! flux into the domain is positive case (value) if(flow(0) .gt. 0.d0) then ! Karline: CHECK! cu(0) = flow(0)*Yup else cu(0) = flow(0)*Y(1) end if case (oneSided) if(flow(0) .ge. 0.d0) then cu(0) = flow(0)*Y(1) else cu(0) = 0.d0 end if case (zeroDivergence) cu(0) = cu(1) case default call rwarn('unkown upstream boundary condition type') end select c the advection step if (mode.eq.0) then ! non-conservative - KARLINE CHECK... do k=1,N Y(k)=Y(k)-1.d0/dble(it)*dt*((cu(k)- cu(k-1))/ V(k) & & -Y(k)*(flow(k)-flow(k-1))/V(k)) enddo else ! conservative - this is actually used do k=1,N Y(k)=Y(k)-1.d0/dble(it)*dt*((cu(k)- cu(k-1))/V(k)) enddo end if end do ! end of the iteration loop c rate of change due to advection do k =1,N dy(k) = (y(k)-dy(k))/dt enddo c Still to do: integrate fluxes in time (now cu = cu of last step ) c flux = 0. ! at start c do k=1,N c flux(k)=flux(k)+1.d0/real(it)*dt*cu(k) return end subroutine advectvol