subroutine inidoublelayer(u,imx1)

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,j,k,imx1
integer u(25)
double precision z0(imx1),u0(imx1),zz
double precision nudzt(nzst),nudzb(nzsb)
double precision kt(imx1),kb(imx1),lt(imx1),lb(imx1)
double precision qab,qat

! pore water pressure at ground water table
u0(:)=0.d0

!initiating z and t
z0(:)=0.d0

!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/)

do 1000, i=1,imx1

if  (slo(i)<slomin)     cycle
if  (ct(zo(i))>1.e6)    cycle

! steady infiltration rate
qab=rikzero(i)*kst(zo(i))/ksb(zo(i))
qat=rikzero(i)

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

! 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
! The sequence of sublayer computation is downward (From surface inward)
deltazt(i,:)=nudzt(:)*ltstar(i)
deltadzt(i,:)=alphat(zo(i))*deltazt(i,:)*cos(slo(i))**2.
zz=0.
do j=1,nzst+1
    if ((j-1)>0) zz=zz+deltazt(i,j-1)
    zt(i,j)=z0(i)+ltstar(i)-zz
    dzt(i,j)=alphat(zo(i))*zt(i,j)*cos(slo(i))**2.
    
    kkt0(i,j)=qat-(qat-(qab-(qab-exp(alphab(zo(i))*u0(i)))*exp(-lb(i)))**(alphat(zo(i))/alphab(zo(i))))*exp(-dzt(i,j))
    pt0(i,j)=1./alphat(zo(i))*dlog(kkt0(i,j))
    inithzt(i,j)=thresit(zo(i))+(thsatt(zo(i))-thresit(zo(i)))*exp(alphat(zo(i))*pt0(i,j))
    inidesatt(i,j)=inithzt(i,j)/thsatt(zo(i))
end do
! lower layer
deltazb(i,:)=nudzb(:)*lbstar(i)
deltadzb(i,:)=alphab(zo(i))*deltazb(i,:)*cos(slo(i))**2.
zz=0.
do j=1,nzsb+1
    if ((j-1)>0) zz=zz+deltazb(i,j-1)
    zb(i,j)=z0(i)-zz
    dzb(i,j)=alphab(zo(i))*zb(i,j)*cos(slo(i))**2.
    
    kkb0(i,j)=qab-(qab-exp(alphab(zo(i))*u0(i)))*exp(-(lb(i)+dzb(i,j)))
    pb0(i,j)=1./alphab(zo(i))*dlog(kkb0(i,j))
end do

! do 1000, i=1,imx1
1000 continue

pb=pb0
pt=pt0
kkb=kkb0
kkt=kkt0

pb1=pb
pt1=pt
kkb1=kkb
kkt1=kkt
kkb2=kkb
kkt2=kkt
! make kkb1(1) and kkb2(1) equal to 1, so that pb(1) is 0.

end subroutine inidoublelayer