subroutine wfs(imx1,nrow,ncol,header,u,maxdirection)
! subroutine for debris flow simulation
! When calculating, only the cells containing data are considered. Nodata cells are added in the output process.
use model_vars;use input_vars; use grids;use input_file_defs;use output_file_defs
use inflow_vars; use outflow_vars
implicit none
integer:: i,ii,j,k,m,n,nsteps,ncol,nrow,trace,imx1,nq,mq,inflowsize,vv,vvmax,maxdirection,nq1,mntout
integer(8):: nt
integer:: u(25),ncc,nccs
double precision:: tnext,tempdt
double precision:: tempri(imx1),tempinflowh(imx1),tempinflowrho(imx1),inflx(imx1),fhw(imx1),fhpredi1(imx1),frhopredi1(imx1)
double precision:: tempele(imx1),eleori(imx1),fdepothick(imx1),fdepovolume
double precision:: ci(imx1),tempci(imx1)
double precision:: fvpredi(imx1,maxdirection),fvpredi2(imx1,maxdirection),fvesti,fvold,fvnew,ffv,ffvprime,convv,localvdiff,fvlimit,dt0,dv,maxfvv
double precision:: fhpredi(imx1),frhopredi(imx1),cv(imx1),erorate(imx1),deporate(imx1),fhpredi2(imx1),frhopredi2(imx1),fybar(imx1,maxdirection),fq(imx1,maxdirection)
double precision:: sf,sfy,sfmiu,sfmanning,coemiu,coemanning,tao,taoc,grav,normfric,absubar(imx1),normfriccoe,miudebris,gammadeb,rhoero,rhodepo(imx1),manningb,manningm
double precision:: tol,cvtol,manningi,manningnq,fhmax,currentfr,manningori(imx1),dmanning
double precision:: hi,hn,si,gd,grad,bedslope,ybar,agrad,cvbar,normfriccoebar,miubar,hbar,manningbar,frhobar,frhoflux
double precision:: parai,paran,aa,bb,cc,dd,vel,qq(imx1,maxdirection),qnet(imx1),hinflow,cellarea,qqmass(imx1,maxdirection),qmassnet(imx1),dttest,dfhtest,dpfhtest,fsign
double precision:: hour,cosslope
character (len=25):: stp
character (len=255):: outfil
character (len=14)::  header(6)
double precision:: totalvolume,totalflowvolume,totaldepositvolume,totaloutflowvolume,totalinfilvolume,totalinflowvolume,totalrivolume,tempvolume
double precision:: tempinfilvolume,tempinflowvolume,temprivolume,tempoutflowvolume,tempflowvolume,tempdepositvolume,totalerosionvolume,temperosionvolume,totaldepovolume,tempdepovolume
double precision:: volumeerror,volumerelaerror
double precision:: tempoutflowh(nooutflow),tempoutflowhq(nooutflow)
double precision:: lambdainverse,tanthetae,sinthetae,cvlimit(imx1),fvdepo(imx1)
double precision:: tempinierodithick(imx1),tempdebdepothick(imx1)
double precision:: artivis
double precision:: width,rholimit(imx1),vx,vy,fhtemp,fvtemp,cvtemp

cv=0.
grav=9.81
manningb=0.0538
manningm=6.0896
totalvolume=0.
totalflowvolume=0.
totaldepositvolume=0.
totaloutflowvolume=0.
totalerosionvolume=0.
totalinfilvolume=0.; totalinflowvolume=0.; totalrivolume=0.; tempoutflowvolume=0.
mntout=1
tempoutflowh=0.; tempoutflowhq=0.
outflowhq=0.; outflowht=0.


! rhoero is the density of the bed
rhoero=cvstar*(rhos-rhow)+rhow
rhodepo=rhoero

! initialize
tempri=0.
ci=0;
tempinflowh=0.
tempinflowrho=0.
tempinierodithick=inierodithick
tempdebdepothick=0.
tempele=ele
eleori=ele ! record original elevation data
fvpredi=0.

! tol is the minimum flow depth for velocity computation
tol=0.03
fhmax=1.0
manningori=manning

cellarea=celsiz**2.

outflow=.false.

if (noinflow>0) inflowsize=size(inflowht(noinflow,:))

if (outflowsimul) then
    do k=1,nooutflow
    do i=1,imx1
        if (i==outflowcell(k)) then
            outflow(i)=.true. 
        end if 
    end do
    end do
end if


! main loop
do 1000, nt=1,maxnts

! erorate is the erosion rate
erorate=0.
deporate=0.
fvdepo=0.
rhodepo=rhoero
tnext=tnow+dt
qq=0.
qqmass=0.
fybar=0.


if (tnow<simul) then

if (tnext>ttout) then
tempdt=dt
tnext=ttout
dt=ttout-tnow
end if
if (tnext>simul) then 
tnext=simul
dt=tnext-tnow
end if

! *****************************************************************
! determine the excessive rainfall intensity of each cell
if (rainsimul) then
    do j=1,nper
        if (capt(j)<=tnow .and. tnext<=capt(j+1)) then
            tempri(:)=rideb(:,j)
            exit
        elseif (tnow<=capt(j+1) .and. capt(j+1)<=tnext) then
            if (j<nper) then
                tempri(:)=((capt(j+1)-tnow)*rideb(:,j)+(tnext-capt(j+1))*rideb(:,j+1))/dt
                exit
            else
                tempri(:)=(capt(j+1)-tnow)*rideb(:,j)/dt
                exit
            end if
        elseif (j==nper .and. capt(j+1)<tnow) then
            tempri(:)=0
            exit
        end if
    end do
end if

! *****************************************************************
! determine the information of inflow hydrograph for the needed cell
if (inflowsimul) then
if (noinflow>0) then

    do 1010, k=1,noinflow
    do i=1,imx1
    
        if (i==inflowcell(k)) then
        if (tnow<inflowht(k,inflowsize)) then

            do j=1,inflowsize-1

                if (inflowht(k,j)<=tnow .and. tnext<=inflowht(k,j+1)) then 
                    tempinflowh(i)=inflowhq(k,j+1)*dt/cellarea
                    tempinflowrho(i)=(rhos-rhow)*inflowhcv(k,j+1)+rhow

                    goto 1010
                elseif (tnow<=inflowht(k,j+1) .and. inflowht(k,j+1)<=tnext) then

                    if (j<=inflowsize-2) then
                    tempinflowh(i)=((inflowht(k,j+1)-tnow)*inflowhq(k,j+1)+(tnext-inflowht(k,j+1))*inflowhq(k,j+2))/cellarea
                    tempinflowrho(i)=(rhos-rhow)*inflowhcv(k,j+1)+rhow

                    goto 1010
                    else
                    tempinflowh(i)=(inflowht(k,j+1)-tnow)*inflowhq(k,j+1)/cellarea
                    tempinflowrho(i)=(rhos-rhow)*inflowhcv(k,j+1)+rhow

                    goto 1010
                    end if
                    
                    goto 1010
                end if
                
            end do

        else 

            tempinflowh(i)=0.
            tempinflowrho(i)=0.

        end if  
        end if
        
    end do
    1010 continue
! if (noinflow>0) then
end if
! if (inflowsimul) then
end if

! *****************************************************************
! compute normalized transient infiltration rate ir & rik=Itransient/Ks (rik=ir/Kst)
! compute infiltration at each cell for each time increment dt
if (infilsimul) then
!do i=1,imx1
!    inflx(i)=tempri(i)+(tempinflowh(i)+fh(i)*(1-cv(i)/cvstar))/dt
!    if (inflx(i)<0.) inflx(i)=0.
!    if (cv(i)>0.1) inflx(i)=0.
!    ! case 1. exfiltration at cells where water table is initially at
!    ! at the ground surface; does not track outflow from cells where
!    ! the water table was initially below the surface and later filled up.
!    if (depth(i)==0.0 .and. rizero(i)<0.0) then
!        ir(i)=0.
!    ! case 2. infiltration, Green-Ampt model  
!    else if (kst(zo(i)).lt.inflx(i)) then
!        ir(i)=kst(zo(i))
!    ! case 3. available water < Ks
!    else
!        ir(i)=inflx(i)
!    end if
!end do !the water table was initially below the surface and later filled up.
fhw=fh*(1-cv/cvstar)+tempinflowh+tempri*dt
inflx=(fhw-tol)/dt
where (inflx<0.) inflx=0.
where (cv>0.1) inflx=0.
call infr(imx1,ir,ci,tempci,inflx,dt)
end if

! *****************************************************************
! update debris flow density and flow depth before erosion and deposition calculation
fhpredi1(:)=fh(:)+(tempri(:)-ir(:))*dt+tempinflowh(:)
frhopredi1(:)=(frho(:)*fh(:)+(tempri(:)-ir(:))*dt*rhow+tempinflowh(:)*tempinflowrho(:))/fhpredi1(:)
where (fhpredi1<=0.) fhpredi1=0.
where (fhpredi1<=eps) frhopredi1=rhow

where (outflow==.true.) fhpredi1=0.
where (outflow==.true.) frhopredi1=rhow

! *****************************************************************
! elevation of each cell
tempele=ele

! *****************************************************************
! predict the debris flow density and flow depth in the first predicting step
fhpredi(:)=fhpredi1(:)
frhopredi(:)=frhopredi1(:)

where (fhpredi<=0.) fhpredi=0.
where (fhpredi<=eps) frhopredi=rhow

! out flow
where (outflow==.true.) fhpredi=0.
where (outflow==.true.) frhopredi=rhow


! calculate flow velocities and flow rates
! using the total friction slope
do i=1,imx1
    do ii=1,maxdirection
        nq=fp(i,ii)
        dt0=0.
        if (nq==0) cycle
        if (qq(i,ii)/=0.) cycle
        if (nq<i) cycle
        
        ! hi is the water surface elevation of cell i
        ! hn is the water surface elevation of cell nq
!        if (erosionsimul .or. sepdepositionsimul) then
!        hi=fhpredi(i)+ele(i)+tempinierodithick(i)+tempdebdepothick(i)
!        hn=fhpredi(nq)+ele(nq)+tempinierodithick(i)+tempdebdepothick(i)
!        else
        hi=fhpredi(i)+tempele(i)
        hn=fhpredi(nq)+tempele(nq)
!        end if

        if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
            si=sin(atan((tempele(nq)-tempele(i))/celsiz))
            gd=(fhpredi(nq)-fhpredi(i))/celsiz
        else
            si=sin(atan((tempele(nq)-tempele(i))/celsiz/2.**0.5))
            gd=(fhpredi(nq)-fhpredi(i))/celsiz/2.**0.5
        end if
        
        if ((fhpredi(i)<=tol .and. hi>=hn) .or. (fhpredi(nq)<=tol .and. hn>=hi)) then
        fvpredi(i,ii)=0.
        qq(i,ii)=0.
        qqmass(i,ii)=0.
        cycle
        end if
      
        if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
        grad=(hn-hi)/celsiz
        else
        grad=(hn-hi)/celsiz/2.**0.5
        end if
        ! hbar is the average flow depth of the two cells
        hbar=0.5*(fhpredi(i)+fhpredi(nq))
        ybar=hbar
        fybar(i,ii)=ybar

        if (ybar>tol) then 
            ! average value of equivelant manning
            ! manning coefficient of cell i
!            if (fhpredi(i)<0.06) then
!                manningi=shallown
!            elseif (0.06<=fhpredi(i) .and. fhpredi(i)<0.15) then
!                manningi=shallown/2.
!            elseif (0.15<=fhpredi(i) .and. fhpredi(i)<1.0) then
!                manningi=manning(i)*1.5*exp(-0.4*fhpredi(i)/fhmax)
!            else
!                manningi=manning(i)
!            end if
!            ! manning coefficient of cell nq
!            if (fhpredi(nq)<0.06) then
!                manningnq=shallown
!            elseif (0.06<=fhpredi(nq) .and. fhpredi(nq)<0.15) then
!                manningnq=shallown/2.
!            elseif (0.15<=fhpredi(nq) .and. fhpredi(nq)<1.0) then
!                manningnq=manning(nq)*1.5*exp(-0.4*fhpredi(nq)/fhmax)
!            else
!                manningnq=manning(nq)
!            end if
            if (fhpredi(i)<1.0) then
                manningi=manning(i)*1.5*exp(-0.4*fhpredi(i)/fhmax)
            else
                manningi=manning(i)
            end if
            ! manning coefficient of cell nq
            if (fhpredi(nq)<1.0) then
                manningnq=manning(nq)*1.5*exp(-0.4*fhpredi(nq)/fhmax)
            else
                manningnq=manning(nq)
            end if
            ! average value  
            manningbar=0.5*(abs(manningi)+abs(manningnq))
!            manningbar=0.5*(manning(i)+manning(nq))
           
            ! diffusive wave solution as a first estimate
!            fvesti=sqrt(ybar**(4./3.)/manningbar**2.*abs(grad))
!            fvesti=sign(fvesti,-grad)
            fvesti=sqrt(ybar**(4./3.)/manningbar**2.*abs(si+gd))
            fvesti=sign(fvesti,-(si+gd))
            
            ! If the Newton-Raphson solution fails to converge after 3 iterations, use the diffusive wave solution.
            fvpredi(i,ii)=fvesti
            
            if (abs(fvesti)>3.5) then
            continue
            end if
            
            ! convective acceleration term
            ! outflow velocity is defined as positive, inflow is defined as negative
!            if (ii==1) convv=(fv(i,1)+fv(i,5))*fv(i,1)
!            if (ii==2) convv=(fv(i,2)+fv(i,6))*fv(i,2)
!            if (ii==3) convv=(fv(i,3)+fv(i,7))*fv(i,3)
!            if (ii==4) convv=(fv(i,4)+fv(i,8))*fv(i,4)
!            if (ii==5) convv=(fv(i,5)+fv(i,1))*fv(i,5)
!            if (ii==6) convv=(fv(i,6)+fv(i,2))*fv(i,6)
!            if (ii==7) convv=(fv(i,7)+fv(i,3))*fv(i,7)
!            if (ii==8) convv=(fv(i,8)+fv(i,4))*fv(i,8)
            
            ! Newton-Raphson method to solve for fvpredi(i,ii)
            fvold=fvesti
            do k=1,6
                if (ii==1) convv=(fvold+fvpredi(i,5))*fvold
                if (ii==2) convv=(fvold+fvpredi(i,6))*fvold
                if (ii==3) convv=(fvold+fvpredi(i,7))*fvold
                if (ii==4) convv=(fvold+fvpredi(i,8))*fvold
                if (ii==5) convv=(fvold+fvpredi(i,1))*fvold
                if (ii==6) convv=(fvold+fvpredi(i,2))*fvold
                if (ii==7) convv=(fvold+fvpredi(i,3))*fvold
                if (ii==8) convv=(fvold+fvpredi(i,4))*fvold
                
                if (ii==1) localvdiff=(2*fvold+fvpredi(i,5))
                if (ii==2) localvdiff=(2*fvold+fvpredi(i,6))
                if (ii==3) localvdiff=(2*fvold+fvpredi(i,7))
                if (ii==4) localvdiff=(2*fvold+fvpredi(i,8))
                if (ii==5) localvdiff=(2*fvold+fvpredi(i,1))
                if (ii==6) localvdiff=(2*fvold+fvpredi(i,2))
                if (ii==7) localvdiff=(2*fvold+fvpredi(i,3))
                if (ii==8) localvdiff=(2*fvold+fvpredi(i,4))
                
                if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
                ffv=fvold*abs(fvold)*manningbar**2./ybar**(4./3.)+(fvold-fv(i,ii))/grav/dt+convv/grav/celsiz+si+gd
                ffvprime=abs(fvold)*manningbar**2./ybar**(4./3.)+1.0/grav/dt+localvdiff/grav/celsiz
                else
                ffv=fvold*abs(fvold)*manningbar**2./ybar**(4./3.)+(fvold-fv(i,ii))/grav/dt+convv/grav/celsiz/2.**0.5+si+gd
                ffvprime=abs(fvold)*manningbar**2./ybar**(4./3.)+1.0/grav/dt+localvdiff/grav/celsiz/2.**0.5
                end if
                
                if (abs(ffvprime)<eps) exit
                
                fvnew=fvold-ffv/ffvprime
               
                if (abs(fvnew-fvold)<=10.e-4*abs(fvold)) then
                    fvpredi(i,ii)=fvnew
                    exit
                end if
                fvold=fvnew
            end do
            
!            fvlimit=limitfr*sqrt(grav*ybar)
!            if (abs(fvpredi(i,ii))>fvlimit) fvpredi(i,ii)=sign(fvlimit,fvpredi(i,ii))
            vel=abs(fvpredi(i,ii))
            if (vel>3.5) then
            maxfvv=maxval(maxfv(:))
            continue
            end if
            
            currentfr=vel/sqrt(grav*ybar)
            if (currentfr>limitfr) then
!                if (fhpredi(i)>=1.0) then
                    dmanning=(manning(i)-manningori(i))/manningori(i)
                    if (dmanning<0.002) then
                        manning(i)=manning(i)+0.0002
                    elseif (0.002<=dmanning .and. dmanning<0.005) then
                        manning(i)=manning(i)+0.0001
                    elseif (0.005<=dmanning .and. dmanning<0.01) then
                        manning(i)=manning(i)+0.00002
                    else
                        manning(i)=manning(i)+0.000002
                    end if
                    
                    dmanning=(manning(nq)-manningori(nq))/manningori(nq)
                    if (dmanning<0.002) then
                        manning(nq)=manning(nq)+0.0002
                    elseif (0.002<=dmanning .and. dmanning<0.005) then
                        manning(nq)=manning(nq)+0.0001
                    elseif (0.005<=dmanning .and. dmanning<0.01) then
                        manning(nq)=manning(nq)+0.00002
                    else
                        manning(nq)=manning(nq)+0.000002
                    end if
            else
                    manning(i)=manning(i)-0.0001
                    if (manning(i)<manningori(i)) manning(i)=manningori(i)
                    
                    manning(nq)=manning(nq)-0.0001
                    if (manning(nq)<manningori(nq)) manning(nq)=manningori(nq)

            end if
            
            ! jude whether the time-step is too large using CFL criterion
            if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
            dttest=0.6*celsiz/(1.5*vel+sqrt(grav*ybar))
            else
            dttest=0.6*celsiz*2.**0.5/(1.5*vel+sqrt(grav*ybar))
            end if

            if (isnan(vel)) then
            write (*,*) 'i',i,'ii',ii
            write (*,*) 'aa',aa,'bb',bb,'cc',cc
            write (*,*) 'ybar',ybar,'fhpredi(i)',fhpredi(i),'fhpredi(nq)',fhpredi(nq),'frhobar',frhobar
            write (*,*) 'vel',vel
            pause 'please check the vel'
            end if

            if (dt>dttest) then
            dt=dt-dtd
            if (dt<dtmin) then
            write (*,*) 'dt<dtmin, CFL criterion violated'
            write (*,*) 'i',i
            write (*,*) 'dttest',dttest
            write (*,*) 'fhpredi',fhpredi(i)
            write (*,*) 'fybar 1 2 3 4 are', fybar(i,1),fybar(i,2),fybar(i,3),fybar(i,4)
            write (*,*) 'fybar 5 6 7 8 are', fybar(i,5),fybar(i,6),fybar(i,7),fybar(i,8)
            write (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
            write (*,*) 'fvpredi 5 6 7 8 are', fvpredi(i,5),fvpredi(i,6),fvpredi(i,7),fvpredi(i,8)
            pause
            dt=dtmin
            end if
            goto 1000
            end if
            
            ! jude whether the time-step is too large using the full dynamic wave equation numerical stability criteria
            if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
                !bedslope=sin(atan((tempele(nq)-tempele(i))/celsiz))
                if (si==0.) then
                    dttest=dt+1
                else
                    dttest=wavemax*abs(si)*celsiz**2./(vel*ybar)
                end if
            else
                !bedslope=atan(abs((hn-hi)/celsiz/2.**0.5))
                if (si==0.) then
                    dttest=dt+1
                else
                    dttest=wavemax*abs(si)*(celsiz*2.**0.5)**2./(vel*ybar)
                end if
            end if
            
            if (dt>dttest) then
            dt=dt-dtd
            if (dt<dtmin) then
            write (*,*) 'dt<dtmin, full dynamic wave stability violated'
            write (*,*) 'i',i
            write (*,*) 'dttest',dttest
            write (*,*) 'fhpredi',fhpredi(i)
            write (*,*) 'fybar 1 2 3 4 are', fybar(i,1),fybar(i,2),fybar(i,3),fybar(i,4)
            write (*,*) 'fybar 5 6 7 8 are', fybar(i,5),fybar(i,6),fybar(i,7),fybar(i,8)
            write (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
            write (*,*) 'fvpredi 5 6 7 8 are', fvpredi(i,5),fvpredi(i,6),fvpredi(i,7),fvpredi(i,8)
            pause
            dt=dtmin
            end if
            goto 1000
            end if
            
            ! estimate the debris flow density of the flux
            if (fvpredi(i,ii)>=0.) then
            ybar=min(fhpredi(i),hbar)
            frhoflux=frhopredi(i)
            else
            ybar=min(fhpredi(nq),hbar)
            frhoflux=frhopredi(nq)
            end if

        ! if (ybar>tol)
        else
            fvpredi(i,ii)=0.
        ! if (ybar>tol)
        end if
        
        if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
        width=celsiz/2.
        else
        width=celsiz/(2.**0.5)
        end if
        qq(i,ii)=fvpredi(i,ii)*ybar*width*(dt-dt0)
        qqmass(i,ii)=frhoflux*qq(i,ii)

        ! computing the flux and mass flux of the cell nq to same time
        if (ii==1) then
        qq(nq,5)=-qq(i,ii)
        qqmass(nq,5)=-qqmass(i,ii)
        fvpredi(nq,5)=-fvpredi(i,ii)
        fybar(nq,5)=fybar(i,ii)
        elseif (ii==2) then
        qq(nq,6)=-qq(i,ii)
        qqmass(nq,6)=-qqmass(i,ii)
        fvpredi(nq,6)=-fvpredi(i,ii)
        fybar(nq,6)=fybar(i,ii)
        elseif (ii==3) then
        qq(nq,7)=-qq(i,ii)
        qqmass(nq,7)=-qqmass(i,ii)
        fvpredi(nq,7)=-fvpredi(i,ii)
        fybar(nq,7)=fybar(i,ii)
        elseif (ii==4) then
        qq(nq,8)=-qq(i,ii)
        qqmass(nq,8)=-qqmass(i,ii)
        fvpredi(nq,8)=-fvpredi(i,ii)
        fybar(nq,8)=fybar(i,ii)
        elseif (ii==5) then
        qq(nq,1)=-qq(i,ii)
        qqmass(nq,1)=-qqmass(i,ii)
        fvpredi(nq,1)=-fvpredi(i,ii)
        fybar(nq,1)=fybar(i,ii)
        elseif (ii==6) then
        qq(nq,2)=-qq(i,ii)
        qqmass(nq,2)=-qqmass(i,ii)
        fvpredi(nq,2)=-fvpredi(i,ii)
        fybar(nq,2)=fybar(i,ii)
        elseif (ii==7) then
        qq(nq,3)=-qq(i,ii)
        qqmass(nq,3)=-qqmass(i,ii)
        fvpredi(nq,3)=-fvpredi(i,ii)
        fybar(nq,3)=fybar(i,ii)
        elseif (ii==8) then
        qq(nq,4)=-qq(i,ii)
        qqmass(nq,4)=-qqmass(i,ii)
        fvpredi(nq,4)=-fvpredi(i,ii)
        fybar(nq,4)=fybar(i,ii)
        end if
    ! end do ii
    end do

! end do i
end do

    ! estimate accumulation of inflow or outflow
do i=1,imx1
    if ((i==8808 .or. i==8841) .and. tnow>=34630.) then
        continue
    end if
    qnet(i)=-(qq(i,1)+qq(i,2)+qq(i,3)+qq(i,4)+qq(i,5)+qq(i,6)+qq(i,7)+qq(i,8))
    qmassnet(i)=-(qqmass(i,1)+qqmass(i,2)+qqmass(i,3)+qqmass(i,4)+qqmass(i,5)+qqmass(i,6)+qqmass(i,7)+qqmass(i,8))
    ! estimate the change of flow depth because of inflow or outflow
    hinflow=qnet(i)/cellarea
    fhpredi2(i)=fhpredi(i)+hinflow
    frhopredi2(i)=(frhopredi(i)*fhpredi(i)*cellarea+qmassnet(i))/fhpredi2(i)/cellarea

    if (isnan(fhpredi2(i))) then
    write (*,*) 'fhpredi2',i,'is not a number'
    write (*,*) 'fhpredi',fhpredi(i)
    write (*,*) 'frhopredi',frhopredi(i)
    write (*,*) 'hinflow',hinflow
    write (*,*) 'fhpredi',fhpredi(i)
    pause 'please check the fhpredi2'
    end if

    ! if the flow depth is smaller than 0, decrease the time step
    if (fhpredi2(i)<0.) then
        dt=dt-dtd
        if (dt<dtmin) then
        write (*,*) 'dt<dtmin'
        write (*,*) 'i',i
        write (*,*) 'fhpredi2<0'
        write (*,*) 'fhpredi',fhpredi(i),'fhpredi2',fhpredi2(i)
        write (*,*) 'fybar 1 2 3 4 are', fybar(i,1),fybar(i,2),fybar(i,3),fybar(i,4)
        write (*,*) 'fybar 5 6 7 8 are', fybar(i,5),fybar(i,6),fybar(i,7),fybar(i,8)
        write (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
        write (*,*) 'fvpredi 5 6 7 8 are', fvpredi(i,5),fvpredi(i,6),fvpredi(i,7),fvpredi(i,8)
        write (*,*) 'qq 1 2 3 4 are', qq(i,1),qq(i,2),qq(i,3),qq(i,4)
        pause
        dt=dtmin
        end if
        goto 1000
    end if

    ! determine whether the time step is too large using the change of flow depth
    if (outflow(i)==.true. .or. fhpredi(i)==0.) then
    dfhtest=abs(fhpredi2(i)-fhpredi(i))
    dpfhtest=1.0
    else
    dfhtest=abs(fhpredi2(i)-fhpredi(i))
    dpfhtest=abs((fhpredi2(i)-fhpredi(i))/fhpredi(i))
    end if
    
    if (outflow(i)==.true. .or. fhpredi(i)==0.) then
        if (dfhtest>tol*toldhp) then
        write (*,*) 'dfhtest>toldhp*tol'
        dt=dt-dtd
        if (dt<dtmin) then
        write (*,*) 'dt<dtmin'
        write (*,*) 'i',i
        write (*,*) 'dfhtest',dfhtest,'dpfhtest',dpfhtest
        write (*,*) 'fhpredi',fhpredi(i),'fhpredi2',fhpredi2(i)
        write (*,*) 'fybar 1 2 3 4 are', fybar(i,1),fybar(i,2),fybar(i,3),fybar(i,4)
        write (*,*) 'fybar 5 6 7 8 are', fybar(i,5),fybar(i,6),fybar(i,7),fybar(i,8)
        write (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
        write (*,*) 'fvpredi 5 6 7 8 are', fvpredi(i,5),fvpredi(i,6),fvpredi(i,7),fvpredi(i,8)
        write (*,*) 'qq 1 2 3 4 are', qq(i,1),qq(i,2),qq(i,3),qq(i,4)
        pause
        dt=dtmin
        end if
        goto 1000
        end if
    
    elseif (fhpredi(i)/=0.) then
        if (dpfhtest>toldhp) then
        write (*,*) 'dpfhtest>toldhp'
        dt=dt-dtd
        if (dt<dtmin) then
        write (*,*) 'dt<dtmin'
        write (*,*) 'i',i
        write (*,*) 'dfhtest',dfhtest,'dpfhtest',dpfhtest
        write (*,*) 'fhpredi',fhpredi(i),'fhpredi2',fhpredi2(i)
        write (*,*) 'fybar 1 2 3 4 are', fybar(i,1),fybar(i,2),fybar(i,3),fybar(i,4)
        write (*,*) 'fybar 5 6 7 8 are', fybar(i,5),fybar(i,6),fybar(i,7),fybar(i,8)
        write (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
        write (*,*) 'fvpredi 5 6 7 8 are', fvpredi(i,5),fvpredi(i,6),fvpredi(i,7),fvpredi(i,8)
        write (*,*) 'qq 1 2 3 4 are', qq(i,1),qq(i,2),qq(i,3),qq(i,4)
        pause
        dt=dtmin
        end if
        goto 1000
        end if
    end if
end do

! **********************OUTFLOW VOLUME******************************
! The outflow volume 1
tempvolume=0.
do k=1,nooutflow
    tempoutflowh(k)=fhpredi2(outflowcell(k))
    tempoutflowhq(k)=tempoutflowh(k)*cellarea/dt
    tempvolume=tempvolume+tempoutflowh(k)*cellarea
end do
tempoutflowvolume=totaloutflowvolume+tempvolume

! outflow
where (outflow==.true.) fhpredi2=0.
where (outflow==.true.) frhopredi2=rhow
where (fhpredi2<eps) frhopredi2=rhow

!! The outflow volume 2
!do k=1,nooutflow
!do ii=1,8
!    nq1=fp(outflowcell(k),ii)
!    if (nq1==0) cycle
!    if (outflow(nq1)==.true.) cycle
!    totaloutflowvolume=totaloutflowvolume-qq(outflowcell(k),ii)
!end do
!end do
!------------------------------------------------------------------
! The infiltration volume
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+ir(outflowcell(k))*dt*cellarea
end do
tempinfilvolume=totalinfilvolume+sum(ir*dt*cellarea)-tempvolume
!------------------------------------------------------------------
! The deposition volume
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+deporate(outflowcell(k))*dt*cellarea
end do
tempdepovolume=totaldepovolume+sum(deporate*dt*cellarea)-tempvolume

! **********************INFLOW VOLUME******************************
! The inflow volue
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+tempinflowh(outflowcell(k))*cellarea
end do
tempinflowvolume=totalinflowvolume+sum(tempinflowh*cellarea)-tempvolume
!------------------------------------------------------------------
! The rainfall volume
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+tempri(outflowcell(k))*dt*cellarea
end do
temprivolume=totalrivolume+sum(tempri*dt*cellarea)-tempvolume
!------------------------------------------------------------------
! The erosion volume
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+erorate(outflowcell(k))*dt*cellarea
end do
temperosionvolume=totalerosionvolume+sum(erorate*dt*cellarea)-tempvolume

! **********************VOLUME CONSERVATION************************
! compute the volume conservation of each time step
tempflowvolume=0.
tempdepositvolume=0.
do i=1,imx1
tempflowvolume=tempflowvolume+fhpredi2(i)*cellarea
tempdepositvolume=tempdepositvolume+tempdebdepothick(i)*cellarea
end do
volumeerror=temprivolume+tempinflowvolume+temperosionvolume-tempinfilvolume-tempoutflowvolume-tempflowvolume-tempdepositvolume
volumerelaerror=volumeerror/(temprivolume+tempinflowvolume+temperosionvolume)
if (abs(volumerelaerror)>0.00001) then 
write (*,*) 'VOLUME CONSERVATION FAILED, TIME STEP DECREASES'
dt=dt-dtd
if (dt<dtmin) then
write (*,*) 'dt<dtmin, Volume conservation failed'
write (*,*) 'RAINFALL',temprivolume,' INFLOW',tempinflowvolume
write (*,*) 'INFILTRATION',tempinfilvolume,' OUTFLOW',tempoutflowvolume,' STORAGE',tempflowvolume
write (*,*) 'VOLUME ERROR',volumeerror,' RELATIVE ERROR',volumerelaerror
pause
dt=dtmin
end if
goto 1000
end if

! Volume conserved
totaloutflowvolume=tempoutflowvolume
totalinfilvolume=tempinfilvolume
totalinflowvolume=tempinflowvolume
totalrivolume=temprivolume
totalerosionvolume=temperosionvolume
totaldepovolume=tempdepovolume

! *****************************************************************
! the time step goes forward successfully, the time step can be increased to the next time step
tnow=tnext
dt=dt+dti
if (dt>dtmax) dt=dtmax
ntsdeb=ntsdeb+1
write (*,*) 'nt',nt
write (*,*) 'tnow',tnow,'dt',dt

fh=fhpredi2
frho=frhopredi2
fv=fvpredi
!inierodithick=tempinierodithick
inierodithick=tempinierodithick+abs(deporate)*dt ! depostion becomes erodible material
debdepothick=tempdebdepothick
ele=tempele
cv=(frho-rhow)/(rhos-rhow)
ci=tempci
where (fh<eps) fh=0.

! record the maximum flow velocity of each cell
do i=1,imx1
!do ii=1,maxdirection
do ii=3,6
if (abs(fv(i,ii))>maxfv(i))  maxfv(i)=abs(fv(i,ii))
if (fh(i)>maxfh(i))  maxfh(i)=fh(i)
end do
end do
! record the maximum flux of each outflow element
do k=1,nooutflow
if (tempoutflowhq(k)>maxoutfq(k)) then 
    maxoutfq(k)=tempoutflowhq(k)
    maxoutft(k)=tnow/3600
end if
end do


! output resutls based on different option 
if (tnext>=ttout) then
mntout=mntout+1

ti=tiny(param(m))
	write(stp,'(F10.1)') ttout ! F(10.1) is necessary
	stp=adjustl(stp)
	! save the grid of flow depth
    if (flowdepthsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=fh(i)
	    end do
	    outfil=trim(folder)//trim(fhfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! save the grid of maxximum flow depth
	if (maxflowdepthsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=maxfh(i)
	    end do
	    outfil=trim(folder)//trim(maxfhfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! save the grid of flow velocity
	if (fvsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=0.5*(abs(fv(i,1))+abs(fv(i,2))+abs(fv(i,3))+abs(fv(i,4)))
	    end do
	    outfil=trim(folder)//trim(fvfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! save the grid of maximum flow velocity
	if (maxfvsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=maxfv(i)
	    end do
	    outfil=trim(folder)//trim(maxfvfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	   ! save the grid of erosion depth
	if (erosionsimul .and. erodepthsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=eleori(i)-ele(i)
	      if (tfg(i)<0.) tfg(i)=0.
	    end do
	    outfil=trim(folder)//trim(erodepthfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	! save the grid of deposit depth when simulating water and soil deposition seperately
	if (sepdepositionsimul .and. debdepodepthsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=ele(i)-eleori(i) 
	      if (tfg(i)<0.) tfg(i)=0.
	    end do
	    outfil=trim(folder)//trim(debdepodepthfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
   ! save the grid of deposit depth when simulating water and soil deposition seperately
	if (sepdepositionsimul .and. totaldepthsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=fh(i)+debdepothick(i)
	    end do
	    outfil=trim(folder)//trim(totaldepthfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! save the grid of volumetric sediment concentration
	if (cvsave) then 
	    tfg=0.
	    do i=1,imx1
	      tfg(i)=(frho(i)-rhow)/(rhos-rhow)
	    end do
	    outfil=trim(folder)//trim(cvfil)//trim(suffix)//'_'//trim(stp)//'.txt'
	    write (*,*) 'ttout',ttout
	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! save the outflow process
	outflowhq(:,mntout)=tempoutflowhq
	outflowht(mntout)=ttout/3600
	if (mntout>ntout+1) then
	pause 'check the output times of outflow process'
	end if
	
! update the output time
ttout=ttout+tout
if (tempdt>dt) dt=tempdt ! reuse the previous larger dt

! if (tnext>=ttout) then
end if


! if (tnext<=simul) then
else

! output the outflow process
if (outflowsave) then 
    outfil=trim(folder)//trim(outffil)//trim(suffix)//'.txt'
    call soutf(ntout,u(4),u(19),outfil)
end if

write (*,*) 'The time is larger than the simulation time, program ends'
! Final true deposition volume
fdepothick=ele-eleori
where (fdepothick<0.) fdepothick=0.
fdepovolume=sum(fdepothick*cellarea)
write (*,*) 'The final true deposition volume is ',fdepovolume
write (u(19),*) 'The final true deposition volume is ',fdepovolume

! compute the total volume of debris flow mixture in the study area
totalflowvolume=0.
totaldepositvolume=0.
do i=1,imx1
totalflowvolume=totalflowvolume+fh(i)*cellarea
totaldepositvolume=totaldepositvolume+debdepothick(i)*cellarea
end do
totalvolume=totalflowvolume+totaldepositvolume
write (*,*) '=================================================================='
volumeerror=totalrivolume+totalinflowvolume+totalerosionvolume-totalinfilvolume-totaloutflowvolume-totalflowvolume-totaldepositvolume
volumerelaerror=volumeerror/(totalrivolume+totalinflowvolume+totalerosionvolume)

write (*,*) 'tempoutflowvolume ',tempoutflowvolume,' totaloutflowvolume ',totaloutflowvolume
write (*,*) '=================================================================='
write (*,*) 'Total flow volume is ', totalflowvolume
write (*,*) 'Total deposit volume is ', totaldepositvolume
write (*,*) 'Total volume is ', totalvolume
write (*,*) '=================================================================='
write (*,*) 'Mass balance  inflow - outflow volume'
write (*,*) '*** Inflow (Cubic Meters) ***'
write (*,*) '                                                   Water'
write (*,*) 'Rainfall volume                                    ',totalrivolume
write (*,*) 'Inflow Hydrograph                                  ',totalinflowvolume
write (*,*) 'Erosion volume                                     ',totalerosionvolume
write (*,*) 'Inflow Hydrograph + Rainfall                       ',totalinflowvolume+totalrivolume+totalerosionvolume
write (*,*) '*** Outflow (Cubic Meters) ***'
write (*,*) '                                                   Water'
write (*,*) 'Water lost to infiltration & interception          ',totalinfilvolume
write (*,*) 'Floodplain storage                                 ',totalflowvolume
write (*,*) 'Floodplain deposite                                ',totaldepositvolume
write (*,*) 'Floodplain outflow hydrograph                      ',totaloutflowvolume
write (*,*) 'FLoodplain outflow, infiltration & storage         ',totaloutflowvolume+totalinfilvolume+totalflowvolume+totaldepositvolume
write (*,*) '=================================================================='
write (*,*) 'Volume error is ', volumeerror, ' Volume relative error is ', volumerelaerror

write (u(19),*) 'tempoutflowvolume ',tempoutflowvolume,' totaloutflowvolume ',totaloutflowvolume
write (U(19),*) '=================================================================='
write (u(19),*) 'Total flow volume is ', totalflowvolume
write (u(19),*) 'Total deposit volume is ', totaldepositvolume
write (u(19),*) 'Total volume is ', totalvolume
write (u(19),*) '=================================================================='
write (u(19),*) 'Mass balance  inflow - outflow volume'
write (u(19),*) '*** Inflow (Cubic Meters) ***'
write (u(19),*) '                                               Water'
write (u(19),*) 'Rainfall volume                                ',totalrivolume
write (u(19),*) 'Inflow Hydrograph                              ',totalinflowvolume
write (u(19),*) 'Erosion volume                                 ',totalerosionvolume
write (u(19),*) 'Inflow Hydrograph + Rainfall                   ',totalinflowvolume+totalrivolume+totalerosionvolume
write (u(19),*) '*** Outflow (Cubic Meters) ***'
write (u(19),*) '                                               Water'
write (u(19),*) 'Water lost to infiltration & interception      ',totalinfilvolume
write (u(19),*) 'Floodplain storage                             ',totalflowvolume
write (u(19),*) 'Floodplain deposite                            ',totaldepositvolume
write (u(19),*) 'Floodplain outflow hydrograph                  ',totaloutflowvolume
write (u(19),*) 'FLoodplain outflow, infiltration & storage     ',totaloutflowvolume+totalinfilvolume+totalflowvolume+totaldepositvolume
write (u(19),*) '=================================================================='
write (u(19),*) 'Volume error is ', volumeerror,' Volume relative error is ', volumerelaerror
exit

! if (tnext<=simul) then
end if

! end do i=1,maxnts
1000 continue
end