subroutine dfs(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; use hydro_vars
implicit none
integer:: i,ii,j,k,m,n,nsteps,ncol,nrow,trace,imx1,nq,inflowsize,vv,vvmax,maxdirection,nq1,mntout
integer(8):: nt
integer:: u(25),ncc,nccs,gindx(imx1),eroindx(imx1)
double precision:: tnext,tempdt
double precision:: tempri(imx1),tempinflowh(imx1),tempinflowrho(imx1),inflx(imx1),fhw(imx1),tempfsh(imx1),tempfsrho(imx1),fhpredi1(imx1),frhopredi1(imx1)
double precision:: tempele(imx1),eleori(imx1),fdepothick(imx1),ferothick(imx1),fdepovolume,ferovolume
double precision:: ci(imx1),tempci(imx1),tempir(imx1)
double precision:: fvpredi(imx1,maxdirection),fvpredi2(imx1,maxdirection),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,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),qqt(imx1,maxdirection),qtnet(imx1),hinflow,cellarea,qqmass(imx1,maxdirection),qmassnet(imx1),dttest,dfhtest,dpfhtest,fvlimit,dt0,localvdiff,dv
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,totalfsvolume,tempfsvolume,totaldepovolume,tempdepovolume
double precision:: volumeerror,volumerelaerror
double precision:: tempoutflowh(nooutflow),tempoutflowhq(nooutflow),tempoutflowcv(nooutflow),temphydrohq(nohydrocell),temphydrocv(nohydrocell)
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.
totalfsvolume=0.
totalinfilvolume=0.; totalinflowvolume=0.; totalrivolume=0.; tempoutflowvolume=0.
mntout=1

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

! initialize
tempdt=0.
tempri=0.
fhw=0.
ci=0; tempci=0.;
tempir=0.
gindx=0 ! record slope failure cells
eroindx=0 ! record erosion cells
tempfsh=0.
tempfsrho=0.
tempinflowh=0.
tempinflowrho=0.
tempinierodithick=inierodithick
tempdebdepothick=0.
tempele=ele
eleori=ele ! record original elevation data

! tol is the minimum flow depth for velocity computation
tol=0.03
cvtol=0.1
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


if (erosionsimul .or. sepdepositionsimul) then 
    vvmax=1
    do i=1,imx1
        cvlimit(i)=rhow*tan(slo(i))/((rhos-rhow)*(tan(phit(zo(i)))-tan(slo(i))))
        if (cvlimit(i)<0.15) cvlimit(i)=6.7*cvlimit(i)**2.
        if (cvlimit(i)<0. .or. cvlimit(i)>1.) cvlimit(i)=cvstar
        rholimit(i)=cvlimit(i)*(rhos-rhow)+rhow
    end do
else
    vvmax=1 
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.
tempfsh=0.
tempfsrho=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
!! 8/18/2009 RLB added optional offset of background flux to prevent excessive drying during periods of zero infiltration.	     
!    if(bkgrof) then
!        do i=1,imx1
!            tempri(i)=kst(zo(i))*rikzero(i)+tempri(i)
!        end do
!    end if
    
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
    fhw(i)=fh(i)*(1-cv(i)/cvstar)+tempri(i)*dt+tempinflowh(i)
    if (fhw(i)<tol) fhw(i)=0.
    inflx(i)=fhw(i)/dt
    ! 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))<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


! *****************************************************************
do vv=1,vvmax
!fvpredi2=0.5*(fv+fvpredi)
fvpredi2(:,:)=fv(:,:)
! compute the limit flow velocity for erosion and deposition processes
if (erosionsimul .or. sepdepositionsimul) then 
do i=1,imx1
    cv(i)=(frhopredi1(i)-rhow)/(rhos-rhow)
    if (cv(i)<0.) cv(i)=0
    !write (*,*) 'i',i,'cv<0',cv(i),'frhopredi1',frhopredi1(i),'fhpredi1',fhpredi1(i)
    !pause
    !end if

    if (isnan(cv(i))) then
    write (*,*) 'i',i,'cv is NAN',cv(i),'frhopredi1',frhopredi1(i),'fhpredi1',fhpredi1(i)
    pause
    end if

    ! computing the flow velocity for computation of erosion and deposition
    vx=(abs(fvpredi2(i,2))+abs(fvpredi2(i,3))+abs(fvpredi2(i,6))+abs(fvpredi2(i,7)))/4.
    vy=(abs(fvpredi2(i,4))+abs(fvpredi2(i,5))+abs(fvpredi2(i,8))+abs(fvpredi2(i,1)))/4.
    absubar(i)=(vx**2.+vy**2.)**0.5

    lambdainverse=(cvstar/cv(i))**0.333-1
    tanthetae=(cv(i)*(rhos-rhow)*tan(phit(zo(i))))/(cv(i)*(rhos-rhow)+rhow)
    sinthetae=sin(atan(tanthetae))
    fvdepo(i)=2./5./d50*(grav*sinthetae*frhopredi1(i)/0.02/rhos)**0.5*lambdainverse*fhpredi1(i)**1.5
end do
end if


! *****************************************************************
! determine the erosion rate of each cell
if (erosionsimul) then

do i=1,imx1
    ! cv must be smaller than cvlimit
    if (cv(i)<cvlimit(i)) then
        gammadeb=frhopredi1(i)*grav
        normfriccoe=cos(slo(i))**2.*tan(phit(zo(i)))
        if (cv(i)<=cvtol) then
            sfy=0.
        else
            if (slo(i)>0.175) then
            ! slope is larger than 10 degree
            ! use Mohr-Coulomb yield stress considering particle suspension in flume and channel
            sfy=(1-cs)*cv(i)*(rhos-rhow)/frhopredi1(i)*normfriccoe
            else
            ! use O'Brien yield stress on debris fan
            sfy=alpha1*exp(beta1*cv(i))/frhopredi1(i)/grav/fhpredi1(i)
            end if
        end if
!        if (cv(i)<=cvtol) then
!        miubar=0.
!        elseif (cvtol<cv(i) .and. cv(i)<=0.1) then
        if (cv(i)<=0.1) then
        miudebris=0.001+cv(i)/0.1*(alpha2*exp(beta2*0.1)-0.001)
        else
        miudebris=alpha2*exp(beta2*cv(i))
        end if
        coemiu=kresis*miudebris/8./gammadeb/fhpredi1(i)**2.
        sfmiu=coemiu*absubar(i)
        if (cv(i)>cvtol) then
        manningbar=manning(i)*manningb*exp(manningm*cv(i))
        else
        manningbar=manning(i)
        end if
        coemanning=manningbar**2./fhpredi1(i)**1.333
        sfmanning=coemanning*absubar(i)**2.
        tao=(sfmanning+sfy+sfmiu)*gammadeb*fhpredi1(i)
        taoc=8.7

        if (fhpredi1(i)>0.02 .and. tao>taoc) then
            erorate(i)=kero(zo(i))*(tao-taoc)
            eroindx(i)=1
        end if
        
        if (erorate(i)<0.) then
        write (*,*) 'i',i,'erorate',erorate(i),'cvlimit',cvlimit(i),'cv',cv(i),'fhpredi1',fhpredi1(i),'absubar',absubar(i),'cv(i)<cvlimit(i)'
        pause 
        end if

        !if (cv(i)<cvlimit) then
    else
        cycle
        !if (cv(i)<cvlimit) then
    end if

    ! to make sure the cv is not larger than cvlimit due to erosion in one time step

    if ((frhopredi1(i)*fhpredi1(i)+erorate(i)*dt*rhoero)>(rholimit(i)*(fhpredi1(i)+erorate(i)*dt))) then
    erorate(i)=(rholimit(i)-frhopredi1(i))*fhpredi1(i)/(rhoero-rholimit(i))/dt
    end if

    ! update the thickness of the initial erodible layer and the deposit layer
    if (erorate(i)*dt<=inierodithick(i)) then
    tempinierodithick(i)=inierodithick(i)-erorate(i)*dt
    else
    erorate(i)=(inierodithick(i))/dt
    tempinierodithick(i)=0.
    end if
    !write (*,*) 'i',i,'erorate', erorate(i)
end do 
!if (erosionsimul) then
end if


! *****************************************************************
! determine the deposition rate of each cell
if (sepdepositionsimul) then

do i=1,imx1
    if (cv(i)>cvlimit(i) .and. absubar(i)<2./3.*fvdepo(i)) then
    deporate(i)=coedepo*(1.-3./2.*absubar(i)/fvdepo(i))*(cvlimit(i)-cv(i))/cvstar*absubar(i)

    ! check the deporate rate to make sure it is not too large
    if (abs(deporate(i)*dt)>fhpredi1(i)) deporate(i)=-fhpredi1(i)/dt
    if (abs(deporate(i)*dt*rhodepo(i))>fhpredi1(i)*frhopredi1(i)) rhodepo(i)=-fhpredi1(i)*frhopredi1(i)/deporate(i)/dt

    if ((frhopredi1(i)*fhpredi1(i)+deporate(i)*dt*rhodepo(i))<(rhow*(fhpredi1(i)+deporate(i)*dt))) then
    deporate(i)=(rhow-frhopredi1(i))*fhpredi1(i)/(rhodepo(i)-rhow)/dt
    end if


    tempdebdepothick(i)=debdepothick(i)+abs(deporate(i)*dt)

    !if (tempdebdepothick(i)>0.) then
    !write (*,*) 'i',i,'tempdebdepothick',tempdebdepothick(i),'deporate',deporate(i)
    !pause
    !end if


    end if
end do
!if (sepdepositionsimul) then
end if

tfg=eleori-ele
do i = 1,imx1
    if (tfg(i)>10) then
    write (*,*) i,tfg(i)
    pause
    end if
end do

! *****************************************************************
! compute shallow rain-induced landslide using numerical scheme
! 8/18/2009 RLB added optional offset of background flux to prevent excessive drying during periods of zero infiltration.
tempir=ir    
if(bkgrof) then
    do i=1,imx1
        if (tempir(i)==0.) tempir(i)=kst(zo(i))*rikzero(i)
    end do
end if
if (fssimul) then
call doublelayer(imx1,tempir,tempfsh,tempfsrho,dt,gindx,eroindx,u)
end if

! *****************************************************************
! elevation of each cell
tempele(:)=ele(:)-erorate(:)*dt+abs(deporate(:))*dt-tempfsh(:)

! *****************************************************************
! predict the debris flow density and flow depth in the first predicting step
fhpredi(:)=fhpredi1(:)+(erorate(:)+deporate(:))*dt+tempfsh(:)
frhopredi(:)=(frhopredi1(:)*fhpredi1(:)+erorate(:)*dt*rhoero+deporate(:)*dt*rhodepo(:)+tempfsh(:)*tempfsrho(:))/fhpredi(:)

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
    if (i==7085 .or. i==8350) then
        continue
    end if
    if (outflow(i)==.true.) then
        continue
    end if
    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 (fhpredi(i)<=tol .and. fhpredi(nq)<=tol) then
        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/=0.) then 

            ! average value of cv
            parai=(frhopredi(i)-rhow)/(rhos-rhow)
            paran=(frhopredi(nq)-rhow)/(rhos-rhow)
            if (parai<0.) parai=0.
            if (paran<0.) paran=0.
            cvbar=0.5*(parai+paran)
            ! average value of dynamic viscosity
!            if (cvbar<=cvtol) then
!            miubar=0.
!            elseif (cvbar<0.1 .and. cvbar>cvtol) then
            if (cvbar<0.1) then
            miubar=0.001+cvbar/0.1*(alpha2*exp(beta2*0.1)-0.001)
            else
            miubar=alpha2*exp(beta2*cvbar)
            end if
            ! 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) 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
            ! average value  
!            manningbar=0.5*(abs(manningi)+abs(manningnq))
            manningbar=0.5*(abs(manning(i))+abs(manning(nq)))
            if (cvbar>cvtol) then
            manningbar=manningbar*manningb*exp(manningm*cvbar)
            end if
            ! average value of debris flow density
            frhobar=0.5*(frhopredi(i)+frhopredi(nq))
            if (frhobar<rhow) frhobar=rhow
            ! average value of normal friction coefficient
            cosslope=cos(atan(abs(grad)))
            parai=cosslope**2.*tan(phit(zo(i)))
            paran=cosslope**2.*tan(phit(zo(nq)))
            normfriccoebar=0.5*(parai+paran)


            if (cvbar<=cvtol) then
                sfy=0.
            else
                if (slo(i)>0.175) then
                ! slope is larger than 10 degree
                ! use Mohr-Coulomb yield stress considering particle suspension in flume and channel
                sfy=(1-cs)*cvbar*(rhos-rhow)/frhobar*normfriccoebar
                else
                ! use O'Brien yield stress on debris fan
                sfy=alpha1*exp(beta1*cvbar)/frhobar/grav/ybar
                end if
            end if


            if (sfy>=abs(grad) .and. abs(fv(i,ii))<=eps) then
                fvpredi(i,ii)=0.
            else
                sfmiu=kresis*miubar/(8.*frhobar*grav*ybar**2.)*abs(fv(i,ii))
                sfmanning=manningbar**2./ybar**1.333*abs(fv(i,ii))**2.
                sf=sfy+sfmiu+sfmanning
                !sf=manningbar**2./ybar**1.333*abs(fv(i,ii))**2.

                if (fv(i,ii)==0.) then
                sf=sign(sf,-grad)
                else
                sf=sign(sf,fv(i,ii))
                end if

                if (ii==1) localvdiff=0.5*(fv(nq,1)+fv(i,5))
                if (ii==2) localvdiff=0.5*(fv(nq,2)+fv(i,6))
                if (ii==3) localvdiff=0.5*(fv(nq,3)+fv(i,7))
                if (ii==4) localvdiff=0.5*(fv(nq,4)+fv(i,8))
                if (ii==5) localvdiff=0.5*(fv(nq,5)+fv(i,1))
                if (ii==6) localvdiff=0.5*(fv(nq,6)+fv(i,2))
                if (ii==7) localvdiff=0.5*(fv(nq,7)+fv(i,3))
                if (ii==8) localvdiff=0.5*(fv(nq,8)+fv(i,4))


                ! artificial viscosity
                if (ii==1) artivis=fv(nq,1)-2.*fv(i,1)-fv(i,5)
                if (ii==2) artivis=fv(nq,2)-2.*fv(i,2)-fv(i,6)
                if (ii==3) artivis=fv(nq,3)-2.*fv(i,3)-fv(i,7)
                if (ii==4) artivis=fv(nq,4)-2.*fv(i,4)-fv(i,8)
                if (ii==5) artivis=fv(nq,5)-2.*fv(i,5)-fv(i,1)
                if (ii==6) artivis=fv(nq,6)-2.*fv(i,6)-fv(i,2)
                if (ii==7) artivis=fv(nq,7)-2.*fv(i,7)-fv(i,3)
                if (ii==8) artivis=fv(nq,8)-2.*fv(i,8)-fv(i,4)

                if (ii==1 .or. ii==3 .or. ii==5 .or. ii==7) then
                dv=(-grad-sf-fv(i,ii)*localvdiff/celsiz/grav)*grav*dt+0.02*abs(fhpredi(i)-fhpredi(nq))/(fhpredi(i)+fhpredi(nq))*artivis-fv(i,ii)*(tempfsh(i)/dt+erorate(i)+deporate(i))*dt/ybar
                else
                dv=(-grad-sf-fv(i,ii)*localvdiff/celsiz/grav/2.**0.5)*grav*dt+0.02*abs(fhpredi(i)-fhpredi(nq))/(fhpredi(i)+fhpredi(nq))*artivis-fv(i,ii)*(tempfsh(i)/dt+erorate(i)+deporate(i))*dt/ybar
                end if
                
                fvpredi(i,ii)=dv+fv(i,ii)

                ! check whether the signs of fv and fvpredi are not the same
                if (fv(i,ii)*fvpredi(i,ii)<0) then
                
                    dt0=-fv(i,ii)/(dv/dt)

                    ! if the yield stress is larger than the total head gradient, fvpredi is zero
                    if (sfy>=abs(grad)) then
                    fvpredi(i,ii)=0.
                    else
                    ! acceleration in the other direction
                    sfy=sign(sfy,-grad)
                    fvpredi(i,ii)=(-grad-sfy)*grav*(dt-dt0)
                    end if
                
                end if

            end if

!            currentfr=fvpredi(i,ii)/sqrt(grav*ybar)
!            if (currentfr>limitfr .and. fhpredi(nq)>=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
!            else
!                manning(i)=manning(i)-0.0001
!                if (manning(i)<manningori(i)) manning(i)=manningori(i)
!            end if
            
            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
            continue
            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/(vel+sqrt(grav*ybar))
            else
            dttest=0.6*celsiz*2.**0.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) dt=dtmin
            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=atan(abs((tempele(nq)-tempele(i))/celsiz))
!                if (bedslope==0.) then
!                    dttest=dt+1
!                else
!                    dttest=wavemax*bedslope*celsiz**2./(vel*ybar)
!                end if
!            else
!                bedslope=atan(abs((tempele(nq)-tempele(i))/celsiz*2.**0.5))
!                if (bedslope==0.) then
!                    dttest=dt+1
!                else
!                    dttest=wavemax*bedslope*(celsiz*2.**0.5)**2./(vel*ybar)
!                end if
!            end if
!            
!            if (dt>dttest) then
!            dt=dt-dtd
!            if (dt<dtmin) dt=dtmin
!            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/=0.)
        else
            fvpredi(i,ii)=0.
        ! if (ybar/=0.)
        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
        qqt(i,ii)=fvpredi(i,ii)*ybar*width
        qq(i,ii)=qqt(i,ii)*(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
        qqt(nq,5)=-qqt(i,ii)
        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
        qqt(nq,6)=-qqt(i,ii)
        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
        qqt(nq,7)=-qqt(i,ii)
        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
        qqt(nq,8)=-qqt(i,ii)
        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
        qqt(nq,1)=-qqt(i,ii)
        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
        qqt(nq,2)=-qqt(i,ii)
        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
        qqt(nq,3)=-qqt(i,ii)
        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
        qqt(nq,4)=-qqt(i,ii)
        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
    qtnet(i)=-(qqt(i,1)+qqt(i,2)+qqt(i,3)+qqt(i,4)+qqt(i,5)+qqt(i,6)+qqt(i,7)+qqt(i,8))
    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 (frhopredi2(i)<995.0) then
    write (*,*) 'frhopredi2(i)', frhopredi2(i),'<995'
    write (*,*) 'frhopredi(i)', frhopredi(i)
    dt=dt-dtd
    if (dt<dtmin) dt=dtmin
    goto 1000
    end if
    
    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
    write (*,*) 'i',i
    write (*,*) 'fhpredi2<0'
    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 (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
    write (*,*) 'qq 1 2 3 4 are', qq(i,1),qq(i,2),qq(i,3),qq(i,4)
    !pause
    dt=dt-dtd
    if (dt<dtmin) dt=dtmin
    goto 1000
    end if

    ! determine whether the time step is too large using the change of flow depth
    if (outflow(i)==.true.) then
    dfhtest=abs(fhpredi2(i)-fhpredi(i))
    dpfhtest=0.1
    else
    dfhtest=abs(fhpredi2(i)-fhpredi(i))
    dpfhtest=abs((fhpredi2(i)-fhpredi(i))/fhpredi(i))
    end if
    if (dfhtest>toldh .and. dpfhtest>toldhp) 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 (*,*) 'fvpredi 1 2 3 4 are', fvpredi(i,1),fvpredi(i,2),fvpredi(i,3),fvpredi(i,4)
    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 do

! end do vv=1,vvmax
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
    tempoutflowcv(k)=(frhopredi2(outflowcell(k))-rhow)/(rhos-rhow)
    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


! **********************HYDROGRAPH ELEMENT******************************
temphydrohq(:)=0.
do k=1,nohydrocell
    do ii=1,maxdirection
        if (-qqt(hydrocell(k),ii)>0.) temphydrohq(k)=temphydrohq(k)+(-qqt(hydrocell(k),ii))
    end do
    temphydrocv(k)=(frhopredi2(hydrocell(k))-rhow)/(rhos-rhow)
    
    if (isnan(temphydrocv(k))) then
    write (*,*) 'k',k,'temphydrocv is NAN',temphydrocv(k),'frhopredi2',frhopredi2(hydrocell(k)),'fhpredi2',fhpredi2(hydrocell(k))
    pause
    end if
end do


!! 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+abs(deporate(outflowcell(k)))*dt*cellarea
end do
tempdepovolume=totaldepovolume+sum(abs(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
!------------------------------------------------------------------
! The landslide volume
tempvolume=0.
do k=1,nooutflow
    tempvolume=tempvolume+tempfsh(outflowcell(k))*cellarea
end do
tempfsvolume=totalfsvolume+sum(tempfsh*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+tempfsvolume-tempinfilvolume-tempoutflowvolume-tempflowvolume-tempdepositvolume
volumerelaerror=volumeerror/(temprivolume+tempinflowvolume+temperosionvolume+tempfsvolume)
if (abs(volumerelaerror)>0.00001) then 
write (*,*) 'VOLUME CONSERVATION FAILED, TIME STEP DECREASES'
!pause
dt=dt-dtd
if (dt<dtmin) dt=dtmin
goto 1000
end if

! Volume conserved
totaloutflowvolume=tempoutflowvolume
totalinfilvolume=tempinfilvolume
totalinflowvolume=tempinflowvolume
totalrivolume=temprivolume
totalerosionvolume=temperosionvolume
totalfsvolume=tempfsvolume
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)
kkt=kkt2
kkb=kkb2
pt=pt2
pb=pb2
!ci=tempci
where (fh<eps) fh=0.

! record the maximum flow velocity of each cell
do i=1,imx1
do ii=1,maxdirection
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

! record the maximum flux of each hydrograph element
do k=1,nohydrocell
if (temphydrohq(k)>maxhydrofq(k)) then 
    maxhydrofq(k)=temphydrohq(k)
    maxhydroft(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)
	
	! minimum factor of safety
    if (fsminsave) then
        tfg=0.
	    do i=1,imx1
	      tfg(i)=fsmin(i)
	    end do
	    outfil=trim(folder)//trim(fminfil)//trim(suffix)//'_'//trim(stp)//'.txt'
   	    call ssvgrd(tfg,imax,pf1,nrow,ncol,u(4),test1,param,u(19),&
            &outfil,ti,header)
	end if
	
	! 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)=cv(i)
	    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(:)
	outflowcv(:,mntout)=tempoutflowcv(:)
	outflowht(mntout)=ttout/3600
	if (mntout>ntout+1) then
	pause 'check the output times of outflow process'
	end if
	
	! save the hydrograph of specified cells
	hydrohq(:,mntout)=temphydrohq(:)
	hydrocv(:,mntout)=temphydrocv(:)
	hydroht(mntout)=ttout/3600
	if (mntout>ntout+1) then
	pause 'check the output times of htdrograph 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

! output the hydrograph process
if (hydrosave) then 
    outfil=trim(folder)//trim(hydrofil)//trim(suffix)//'.txt'
    call shydro(ntout,u(4),u(19),outfil)
end if

write (*,*) 'The time is larger than the simulation time, program ends'

! Final true erosion volume
ferothick=eleori-ele
where (ferothick<0.) ferothick=0.
ferovolume=sum(ferothick*cellarea)
write (*,*) 'The final true erosion volume is ',ferovolume
write (u(19),*) 'The final true erosion volume is ',ferovolume

! 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