subroutine doublelayer(u1,tout,finf,fmn)

use input_vars
use model_vars
implicit none


integer nra,nrb,nrc
integer i,j,k,m,u1,nt
double precision tnext,slope
double precision z0,dtp,finf,fmn,tout(nout),zz
double precision,allocatable:: zt(:),zb(:),deltazt(:),deltazb(:),dzt(:),dzb(:),deltadzt(:),deltadzb(:)
double precision kt,kb,u0,lt,lb,rka,rkb,rkst,rksb,vt,vb,rt,rb
double precision,allocatable:: kkb(:),kkt(:),kkb0(:),kkt0(:),kkb1(:),kkt1(:),kkb2(:),kkt2(:)
double precision qab,qat,qbb,qbt,rqa,rqaold,deltarqa,rqb,rqbold,deltarqb,rqc,rqcold,deltarqc,fs,uwsum,uwt1,uwt2,fft,ffb,fsw,fsc
double precision sdfs,partfs,uwsfin,xindx,pf
double precision nudzt(nzst),nudzb(nzsb)

!constant parameter
pi=3.141592653587973d0
eps=1.0d-18

!parameters of double layer
slope=slo
u0=0.d0

qab=qa/ksb
qat=qa/kst

qbb=qb/ksb
qbt=qb/kst

!fft=tan(phit)/tan(slope)
!ffb=tan(phib)/tan(slope)

!initiating z and t
z0=0.d0
allocate (kkt(nzst+1),kkb(nzsb+1),kkt0(nzst+1),kkb0(nzsb+1),kkt1(nzst+1),kkb1(nzsb+1),kkt2(nzst+1),kkb2(nzsb+1))
allocate (zt(nzst+1),zb(nzsb+1),deltazt(nzst),deltazb(nzsb),dzt(nzst+1),dzb(nzsb+1),deltadzt(nzst),deltadzb(nzsb))
kkb=0.; kkt=0.
kkb0=0.; kkt0=0.
kkb1=0.; kkt1=0.
kkb2=0.; kkt2=0.

!interval of t and z
nudzt=(/0.01,0.01,0.01,0.015,0.015,0.02,0.03,0.04,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07, &
 &	0.04,0.03,0.02,0.015,0.015,0.01,0.01,0.01/)
nudzb=(/0.01,0.01,0.01,0.015,0.015,0.02,0.03,0.04,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07, &
 &	0.04,0.03,0.02,0.015,0.015,0.01,0.01,0.01/)
!nudzt=(/0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05/)
!nudzb=(/0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05/)


!parameters of double layer
beta=alphab*ksb*(thsatt-thresit)/alphat/kst/(thsatb-thresib) 
rkst=kst
rksb=ksb
lt=alphat*ltstar*cos(slo)**2.
lb=alphab*lbstar*cos(slo)**2.

dtp=alphab*ksb*dt/(thsatb-thresib)*cos(slope)**2.

! initiating unit weight
uwsum=0.
uwspt=uwst
uwspb=uwsb
uwt1=0.


! initial condition
! The foward direction of z is upward
! The datum (i.e., z0) is assumed to be at the interface of the two layers.
! upper layer
zz=z0
zt(1)=z0
dzt(1)=alphat*zt(1)*cos(slope)**2.
kkt0(1)=qat-(qat-(qab-(qab-exp(alphab*u0))*exp(-lb))**(alphat/alphab))*exp(-dzt(1))
pt0(1)=1./alphat*dlog(kkt0(1))
do i=1,nzst 
    deltazt(i)=nudzt(i)*ltstar 
    zz=zz+deltazt(i)
    zt(i+1)=zz
    
    deltadzt(i)=alphat*deltazt(i)*cos(slope)**2.
    dzt(i+1)=alphat*zt(i+1)*cos(slope)**2.
    
    kkt0(i+1)=qat-(qat-(qab-(qab-exp(alphab*u0))*exp(-lb))**(alphat/alphab))*exp(-dzt(i+1))
    pt0(i+1)=1./alphat*dlog(kkt0(i+1))
end do
! lower layer
zz=-lbstar+z0
zb(1)=-lbstar+z0
dzb(1)=alphab*zb(1)*cos(slope)**2.
kkb0(1)=qab-(qab-exp(alphab*u0))*exp(-(lb+dzb(1)))
pb0(1)=1./alphab*dlog(kkb0(1))
do i=1,nzsb 
    deltazb(i)=nudzb(i)*lbstar
    zz=zz+deltazb(i)
    zb(i+1)=zz
    
    deltadzb(i)=alphab*deltazb(i)*cos(slope)**2.
    dzb(i+1)=alphab*zb(i+1)*cos(slope)**2.
    
    kkb0(i+1)=qab-(qab-exp(alphab*u0))*exp(-(lb+dzb(i+1)))
    pb0(i+1)=1./alphab*dlog(kkb0(i+1))
end do

write (u1,*) 'tsave',tnow
do k=1,nzst+1
    write (u1,*) pt0(nzst+2-k)
end do
do k=1,nzsb+1
    write (u1,*) pb0(nzsb+2-k)
end do

pt1=pt0
pb1=pb0
pt2=pt0
pb2=pb0
kkt1=kkt0
kkb1=kkb0
kkt2=kkt0
kkb2=kkb0

! time to save
j=1
ttout=tout(j)

if (dtp>deltadzt(1)**2./2. .or. dtp>deltadzb(1)**2./2.) then
    pause 'Decrease time step!'
    stop '0'
end if
    
! main loop
do 1000, nt=1,maxnts

tnext=tnow+dt

if (tnext<=simul) then

    ! bottom layer
    do i=2,nzsb
        kkb2(i)=kkb1(i)+dtp/(deltadzb(i)+deltadzb(i-1))*(kkb1(i+1)-kkb1(i-1))+&
             &  dtp/(deltadzt(i-1)*deltadzb(i)**2.)*(deltadzb(i-1)*kkb1(i+1)-(deltadzb(i)+deltadzb(i-1))*kkb1(i)+deltadzb(i)*kkb1(i-1))
    end do

    ! top layer, with beta
    do i=2,nzst
        kkt2(i)=kkt1(i)+dtp/beta/(deltadzt(i)+deltadzt(i-1))*(kkt1(i+1)-kkt1(i-1))+&
             &  dtp/beta/(deltadzt(i-1)*deltadzt(i)**2.)*(deltadzt(i-1)*kkt1(i+1)-(deltadzt(i)+deltadzt(i-1))*kkt1(i)+deltadzt(i)*kkt1(i-1))
    end do
    
    kkt2(nzst+1)=(qbt*deltadzt(nzst)+kkt2(nzst))/(deltadzt(nzst)+1)
    
    ! interface of the two layers
    kkt2(1)=(rksb/deltadzb(nzsb)*kkb2(nzsb)+rkst/deltadzt(1)*kkt2(2))/(rksb/deltadzb(nzsb)*(1+deltadzb(nzsb))+rkst/deltadzt(1)*(1-deltadzt(1)))
    kkb2(nzsb+1)=kkt2(1)

    pt2=1./alphat*dlog(kkt2)
    pb2=1./alphab*dlog(kkb2)

    !thzt(i)=thresit+(thsatt-thresit)*exp(alphat*pt(i))
    !desatt(i)=thzt(i)/thsatt
    !thzb(i)=thresib+(thsatb-thresib)*exp(alphab*pb(i))
    !desatb(i)=thzb(i)/thsatb

    if (tnext>=ttout) then
        write (u1,*) 'tsave',ttout
        do k=1,nzst+1
            write (u1,*) pt2(nzst+2-k)
        end do
        do k=1,nzsb+1
            write (u1,*) pb2(nzsb+2-k)
        end do
        j=j+1
        if (j<=nout) then
            ttout=tout(j)
        else
            ttout=simul
        end if
    end if

    kkt1=kkt2
    kkb1=kkb2
    pt1=pt2
    pb1=pb2
    tnow=tnext

else

    exit
   
! if (tnext<=simul) then
end if

! end do i=1,maxnts
1000 continue

end subroutine doublelayer