program EDDA
! EDDA:integrated simulation of debris flow erosion, deposition and property changes
! by H.X. CHEN, HKUST, 2015; P. SHEN, HKUST, 2016
use input_file_defs; use input_vars; use outflow_vars; use hydro_vars
use grids; use model_vars;use output_file_defs
implicit none
integer::grd
integer:: i,ii,j,k,l,m,n,imx1,mnd,maxzo,maxdirection
integer:: nodata,sctr
integer:: ncol,nrow,u(25)
character (len=255):: outfil 
character (len=14):: header(6)
character (len=8):: date
character (len=10):: time
character (len=4):: stp
character (len=31):: scratch,irfil
character (len=7):: vrsn
character (len=11):: bldate

! first executable statement ............	
call date_and_time(date,time)
u=(/11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,&
 &	28,29,30,31,32,33,34,35/)
test=-9999.D0; test1=-9999.
pi=3.141592653589793
eps=1.d-18 ! d: double precision  e: single precision
dg2rad=pi/180.D0
vrsn='1.5'; bldate='15 June 2016'
write (*,*) ''
write (*,*) 'EDDA:integrated simulation of debris flow erosion,'
write(*,*) 'deposition and property changes'
write (*,*) '       Version ', vrsn,', ',bldate
write (*,*) '  By Hongxin CHEN and Ping SHEN'
write (*,*) '       The Hong Kong University of Science and Technology'
write (*,*) '-----------------------------------------'
	
	
! open log file
outfil='EDDALog.txt'
open (u(19),file=trim(outfil),status='unknown',err=410)
write (u(19),*) ''
write (u(19),*) 'Starting EDDA ', vrsn,' ',bldate
write (u(19),*) 'Date ',date(5:6),'/',date(7:8),'/',date(1:4)
write (u(19),*) 'Time ',time(1:2),':',time(3:4),':',time(5:6)
!  read initialization file
call trini(u(19),u(1),dg2rad)

!allocate & initiate arrays for debris flow simulation
grd=row*col
imx1=imax  ! data nodes  imax <= grd
rhow=1000.
rhos=2650.
allocate (pf1(grd),pf2(grd))
allocate (manning(imax),ele(imax),ltstar(imax),lbstar(imax),inierodithick(imax),debdepothick(imax),fh(imax),frho(imax),zo(imax),sctrtrace(imax),slo(imax))
allocate (stemp(col),temp(col))
allocate (fp(imax,8),fv(imax,8),maxfv(imax),maxfh(imax))
allocate (ri(imax),outflow(imax))
allocate (tfg(imax))
allocate (indx(imax),nxt(imax),dir(imax))
allocate (dsctr(imax+1),dsc(nwf))
allocate (ele1(grd))
allocate (rizero(imax),ro(imax),ir(imax))
allocate (rik(imax*nper))
allocate (wf(nwf),rikzero(imax))
allocate (depth(imax),zmax(imax))
allocate (finalcell(imax),finalvolume(imax),finaldistance(imax),finaldeltah(imax))

pf1=0.
pf2=0 ! pf2 is integer
ele=0.
ele1=0.
manning=0.
ltstar=0.
lbstar=0.
inierodithick=0.
debdepothick=0.
fh=0.
frho=1000
zo=0
sctrtrace=0
slo=0.
maxfv=0.
maxfh=0.

stemp=0.
temp=0
fp=0
outflow=.false.
indx=0
nxt=0
dsctr=0
dsc=0

rizero=0.
ri=0.
rik=0.
ro=0.
wf=0.
rikzero=0.
ir=0.
depth=0.
zmax=0.



! fv(i,1) is the flow velocity cross the northern wall of cell i
! fv(i,2) is the flow velocity cross the northeastern wall of cell i
! fv(i,3) is the flow velocity cross the eastern wall of cell i
! fv(i,4) is the flow velocity cross the southeastern wall of cell i
! fv(i,5) is the flow velocity cross the southern wall of cell i
! fv(i,6) is the flow velocity cross the southwestern wall of cell i
! fv(i,7) is the flow velocity cross the western wall of cell i
! fv(i,8) is the flow velocity cross the northwestern wall of cell i
fv=0.

! *****************************************************************
! read gridded data from GIS
write (*,*) 'Reading input grids'
write(u(19),*) 'Input file name,            Cell count'

! read dem file
call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
      &ele,ele1,sctr,imax,stemp,u(8),demfil,param,header,u(19))
write(u(19),*) 'Dem file  grid'
write(u(19),*) trim(demfil),sctr,' data cells'
if(sctr/=imx1) then 
    write (u(19),*) 'Grid mismatch ',trim(demfil)
    write (*,*) 'Grid mismatch ',trim(demfil)
    pause 'Please check the file!'
    stop
end if
          
! read slope file
call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
      &slo,pf1,sctr,imax,stemp,u(8),slofil,param,header,u(19))
write(u(19),*) 'Slope file  grid'
write(u(19),*) trim(slofil),sctr,' data cells'
if(sctr/=imx1) then 
    write (u(19),*) 'Grid mismatch ',trim(slofil)
    write (*,*) 'Grid mismatch ',trim(slofil)
    pause 'Please check the file!'
    stop
end if 
slo=slo*dg2rad
     	  
! read upper layer soil depth file, ltstar
if (cltstar>0.) then
    ltstar=cltstar
else
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &ltstar,pf1,sctr,imax,stemp,u(8),ltstarfil,param,header,u(19))
    write(u(19),*) 'Upper layer soil thickness file  grid'
    write(u(19),*) trim(ltstarfil),sctr,' data cells'
    if(sctr/=imx1) then
        write (u(19),*) 'Grid mismatch ',trim(ltstarfil)
        write (*,*) 'Grid mismatch ',trim(ltstarfil)
        pause 'Please check the file!'
        stop
    end if 
end if  
inierodithick=ltstar

! read lower layer soil depth file, lbstar
if (clbstar>0.) then
    lbstar=clbstar
else
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &lbstar,pf1,sctr,imax,stemp,u(8),lbstarfil,param,header,u(19))
    write(u(19),*) 'Lower layer soil thickness file grid'
    write(u(19),*) trim(lbstarfil),sctr,' data cells'
    if(sctr/=imx1) then
        write (u(19),*) 'Grid mismatch ',trim(lbstarfil)
        write (*,*) 'Grid mismatch ',trim(lbstarfil)
        pause 'Please check the file!'
        stop
    end if 
end if  
     	  
! read manning coefficient
if (cmanning>=0) then
    manning=cmanning
else
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &manning,pf1,sctr,imax,stemp,u(8),manningfil,param,header,u(19))
	write(u(19),*) 'Manning file  grid'
 	write(u(19),*) trim(manningfil),sctr,' data cells'
    if(sctr/=imx1) then
        write (u(19),*) 'Grid mismatch ',trim(manningfil)
        write (*,*) 'Grid mismatch ',trim(manningfil)
        pause 'Please check the file!'
        stop
    end if 
end if

! read the zone file
if(nzon==1) then
    zo=1 ! if only one zone, all values of zone grid equal 1.
    write(*,*) 'One property zone, no grid required!'
    write(u(19),*) 'One property zone, no grid required!'
    parami=param ! added 7/29/2008 RLB
else
    call irdgrd(grd,col,ncol,nrow,celsiz,nodata,mnd,&
          &zo,pf2,sctrtrace,sctr,imax,temp,u(15),zonfil,parami,header,u(19))
    write(u(19),*) 'Property zone grid'
    write(u(19),*) trim(zonfil),sctr,' data cells'
    if(sctr/=imx1) then
        write (*,*) 'Grid mismatch ',trim(zonfil)
        write (*,*) 'Correct property-zone grid and/or initializtion file.' 
        write (u(19),*) 'Grid mismatch ',trim(zonfil)
        write (u(19),*) 'Correct property-zone grid and/or initializtion file.' 
        close(u(19))
        pause 'Press return/enter to quit'
        stop '-1'
    end if
    maxzo=maxval(zo)

    if (maxzo/=nzon) then
        write (*,*) 'Maximum zone number does not equal number of property zones!'
        write (*,*) 'Correct property-zone grid and/or initializtion file.' 
        write (u(19),*) 'Maximum zone number does not equal number of property zones!'
        write (u(19),*) 'Correct property-zone grid and/or initializtion file.' 
        close(u(19))
        pause 'Press return/enter to quit'
        stop '-1'
    end if
end if
	
!  read background infiltration rate, rizero (initial infiltration rate)
if (crizero.lt.0.) then 
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &rizero,pf1,sctr,imax,stemp,u(8),rizerofil,param,header,u(19))
    write(u(19),*) 'Background infiltration rate grid'
    write(u(19),*) trim(rizerofil),sctr,' data cells'
    if(sctr/=imx1) then 
        write (u(19),*) 'Grid mismatch ',trim(ltstarfil)
        write (*,*) 'Grid mismatch ',trim(ltstarfil)
        pause 'Please check the file!'
        stop
    end if
else
    rizero=crizero
end if

!  read initial depth to water table, depth
if (dep.lt.0.) then  
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &depth,pf1,sctr,imax,stemp,u(8),depfil,param,header,u(19))
    write(u(19),*) 'Initial water-table depth grid'
    write(u(19),*) trim(depfil),sctr,' data cells'
    if(sctr/=imx1) then
        write (u(19),*) 'Grid mismatch ',trim(depfil)
        write (*,*) 'Grid mismatch ',trim(depfil)
        pause 'Please check the file!'
        stop
    end if
else
    depth=dep
end if

!  read depth to base of potential slide, zmax
if (czmax.lt.0.) then 
    call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
          &zmax,pf1,sctr,imax,stemp,u(8),zfil,param,header,u(19))
    write(u(19),*) 'Maximum depth grid'
    write(u(19),*) trim(zfil),sctr,' data cells'
    if (sctr/=imx1) then
        write (u(19),*) 'Grid mismatch ',trim(zfil)
        write (*,*) 'Grid mismatch ',trim(zfil)
        pause 'Please check the file!'
        stop
    end if
else
    zmax=czmax
end if

! read rainfall information
! rideb is the rainfall intensity of each cell during each rainfall interval for debris flow simulation
allocate (rideb(imax,nper))
rideb=0.                 
do j=1,nper
    if (cri(j).lt.0) then
        call srdgrd(grd,col,ncol,nrow,celsiz,nodat,&
              &ri,pf1,sctr,imax,stemp,u(13),rifil(j),param,header,u(19))
        write(u(19),*) 'Precipitation intensity grid ',j
        write(u(19),*) trim(rifil(j)),sctr,' data cells'
        if(sctr/=imx1) then
            write (u(19),*) 'Grid mismatch ',trim(rifil(j))
            write (*,*) 'Grid mismatch ',trim(rifil(j))
     	    pause 'Please check the file!'
            stop
        end if 
    else
        ri=cri(j)
    end if
    rideb(:,j)=ri(:)
! end do j
end do

! for four flow direction system
! fp(i,1) is the northern nodal point of cell i
! fp(i,2) is the eastern nodal point of cell i
! fp(i,3) is the southern nodal point of cell i
! fp(i,4) is the western nodal point of cell i

! for eight flow direction system
! fp(i,1) is the northern nodal point of cell i
! fp(i,2) is the northeastern nodal point of cell i
! fp(i,3) is the eastern nodal point of cell i
! fp(i,4) is the southeastern nodal point of cell i
! fp(i,5) is the southern nodal point of cell i
! fp(i,6) is the southwestern nodal point of cell i
! fp(i,7) is the western nodal point of cell i
! fp(i,8) is the northwestern nodal point of cell i

call flodir(grd,col,ncol,nrow,celsiz,nodat,&
      &ele,pf1,sctr,imax,stemp,u(8),demfil,param,header,u(19),fp,outflow)
maxdirection=size(fp(1,:))  
    
! read information of inflow hydrographs
if (inflowsimul) call inflow_read(u(19),u(8))
! read outflow cells
if (outflowsimul) call outflow_read(u(19),u(8))
! read hydrograph cells
if (hydrosave) call hydro_read(u(19),u(8))

allocate (maxoutfq(nooutflow),maxoutft(nooutflow))
maxoutfq(:)=0.
maxoutft(:)=0.
allocate (maxhydrofq(nohydrocell),maxhydroft(nohydrocell))
maxhydrofq(:)=0.
maxhydroft(:)=0.

! *****************************************************************
write(u(19),*) '---------------******---------------'
write(*,*) '---------------******---------------'
! test and adjust (if necessary) steady background infiltration rates relvative to Kst, rikzero
call steady(u(19),imx1)

! parameters for slope stability analysis
! compute pore pressure distributions for either fully saturated or partially saturated conditions.
! Partially saturated zone overlies saturated zone
! Allocate and initialize new arrays
allocate (beta(imax))
allocate (fsmin(imax),pmin(imax),zfmin(imax),fdepth(imax),pfmin(imax))
allocate (zt(imx1,nzst+1),zb(imx1,nzsb+1),deltazt(imx1,nzst),deltazb(imx1,nzsb),dzt(imx1,nzst+1),dzb(imx1,nzsb+1),deltadzt(imx1,nzst),deltadzb(imx1,nzsb))
allocate (kkt(imx1,nzst+1),kkb(imx1,nzsb+1),kkt0(imx1,nzst+1),kkb0(imx1,nzsb+1),kkt1(imx1,nzst+1),kkb1(imx1,nzsb+1),kkt2(imx1,nzst+1),kkb2(imx1,nzsb+1))
allocate (pt(imax,nzst+1),pb(imax,nzsb+1),pt0(imax,nzst+1),pb0(imax,nzsb+1),pt1(imax,nzst+1),pb1(imax,nzsb+1),pt2(imax,nzst+1),pb2(imax,nzsb+1))
allocate (fc(imax,nzst+nzsb+1),fw(imax,nzst+nzsb+1),inithzt(imax,nzst+1),inithzb(imax,nzsb+1),inidesatt(imax,nzst+1),inidesatb(imax,nzsb+1))
allocate (desatt(imax,nzst+1),desatb(imax,nzsb+1),thzt(imax,nzst+1),thzb(imax,nzsb+1))
allocate (uwspt(nzst+1),uwspb(nzsb+1))
allocate (nvu(imax),nv(imax),gs(nzon),gst(nzon),gsb(nzon))
fsmin=0.; pmin=0.; zfmin=0.; fdepth=0.; pfmin=0.
zt=0.; zb=0.; deltazt=0.; deltazb=0.; dzt=0.; dzt=0.
kkb=0.; kkt=0.; kkb0=0.; kkt0=0.; kkb1=0.; kkt1=0.; kkb2=0.; kkt2=0.
pt=0.; pb=0.; pt0=0.; pb0=0.; pt1=0.; pb1=0.; pt2=0.; pb2=0.
fc=0.; fw=0.
inithzt=0.; inithzb=0.
inidesatt=0.; inidesatb=0.
desatt=0.; desatb=0.
thzb=0.; thzt=0.
uwspt=0.; uwspb=0.
fc=0.
fw=0.
nv=0
nvu=0

! Initialization of some parameters used in slope stability calculation
! Initial condition of pore water pressure distribution
call inidoublelayer(u,imx1)

!! determine the number of rainfall intervals, kper
!kper=nper
!if (t>capt(nper+1)) then
!    kper=nper+1 
!else
!    do k=1,nper ! find the period that contains t
!        if(t>=capt(k) .and. t<=capt(k+1)) kper=k
!    end do
!end if
!if (tx<1) tx=1
!nts=kper*tx ! number of time-steps from 0 to t
!tns=float(nts)
!tmin=0.
!tmax=t 
!tinc=(tmax-tmin)/tns
!! compute output pointers	
!allocate(jsav(nts+1))  
!jsav=0
!write (u(19),*) '******** Output times ********'
!write (u(19),*) 'number, timestep #,  time'
!write (*,*) '******** Output times ********'
!write (*,*) 'number, timestep #,  time'
!lwarn=.false.
!
!allocate (q(kper))
!
!tis=tiny(x1)
!write(*,*) 'Starting computations of pressure head and factor of safety'
!! we use finite depth unsaturated model in this analysis, other models can be used if the code is changed	
!if(unsat0) then
!call unsfin(imx1,u(19),u(2),ncc,nccs,profil)
!end if

t0=0. ! t0 is the starting time of simulation
dt=1.0 !dtmin ! dt is the time step of simulation, which is set to be dtmin at the beginning
ttout=tout ! tout is the output interval, ttout is the output time
tminimum=99. ! tminimum is used to record the minimum time step
tmaximum=-99 ! tmaximum is used to record the maximum time step
tnow=t0 !  tnow is the present time
ntsdeb=0 ! ntsdeb is used to record the number of time step of debris flow simulation
ntout=simul/tout ! star from 0 time
allocate (outflowhq(nooutflow,ntout+1),outflowht(ntout+1),outflowcv(nooutflow,ntout+1))
allocate (hydrohq(nohydrocell,ntout+1),hydroht(ntout+1),hydrocv(nohydrocell,ntout+1))
outflowhq=0.; outflowht=0.; outflowcv=0.
hydrohq=0.; hydroht=0.; hydrocv=0.


! compute the maximum possbile steps to finish the simulation
maxnts=2*simul/dtmin

if (debrissimul) then
call dfs(imx1,nrow,ncol,header,u,maxdirection)
else
call wfs(imx1,nrow,ncol,header,u,maxdirection)
end if

write (*,*) 'The total number of times is',' ',ntsdeb
avedt=tnow/float(ntsdeb)
write (*,*) 'The average time step interval is',' ',avedt
write (u(19),*) 'The total number of times is',' ',ntsdeb
write (u(19),*) 'The average time step interval is',' ',avedt
close (u(19))

pause 'The program EDDA runs successfully!'
stop '0'
! Error reporting 	
  410	continue
  	write (*,*) 'Error opening output file main program 410'
	write (*,*) '--> ',outfil
	write (*,*) 'Check file path and status'
  	write (u(19),*) 'Error opening output file main program 410'
	write (u(19),*) '--> ',outfil
	write (u(19),*) 'Check file path and status'
	pause 'Press RETURN to exit'
	stop '410'

end