advection.f90 1.02 KB
subroutine advection

use parameters

implicit none

        double precision, allocatable :: diffF(:),diffl(:),diffr(:),phi(:),theta(:)
        double precision, allocatable :: F(:),diff1(:)

        
        allocate(diffl(nbin))
        allocate(diffr(nbin))
        allocate(diffF(nbin))
        allocate(theta(nbin))
        allocate(phi(nbin))
        allocate(F(nbin))
        allocate(diff1(nbin-1))
      
  
        do i=2,nbin
       diffl(i)=E(n-1,i,ii)-E(n-1,i-1,ii)
        end do

        do i=1,nbin-1
        diffr(i)=E(n-1,i+1,ii)-E(n-1,i,ii)
        end do
        diffr(nbin)=-E(n-1,nbin,ii)
        diffl(1)=E(n-1,1,ii)


        do i=1,nbin
        theta(i)=diffl(i)/(diffr(i)+3e-14)
        phi(i)=max(0d0,min(theta(i)*2,1d0),min(theta(i),2d0))
        end do


        F=E(n-1,1:nbin,ii)+0.5*(1-CN(ii))*diffr*phi
        do i=2,nbin
        diffF(i)=F(i)-F(i-1)
        end do
        diffF(1)=F(1)
        E(n,1:nbin,ii)=E(n-1,1:nbin,ii)-CN(ii)*diffF
      
      

end subroutine advection