[Fortran] 纯文本查看 复制代码
subroutine doublelayer(imx1,fave,tempfsh,tempfsrho,dt0,gindx,eroindx,u)
use model_vars; use input_vars;
use grids;
implicit none
integer i,j,k,m,nt,imx1,jmin
integer u(25),gindx(imx1),eroindx(imx1)
double precision z0(imx1),u0(imx1)
double precision finf,fmn,flm(nzst+1)
double precision fave(imx1),tempfsh(imx1),tempfsrho(imx1),inflowh,inflowrho,fmrho
double precision dt0,dtn
double precision dtp,vt,vb,rt,rb
double precision rkst,rksb
double precision qbb,qbt,uwsum,uwt1,uwt2,fs,fft,ffb,fsw,fsc
finf=10.
! pore water pressure at ground water table
u0=0.d0
!initiating z and t
z0=0.d0
nt=5
dtn=dt0/nt
pb1=pb
pt1=pt
kkb1=kkb
kkt1=kkt
!$omp parallel do
do i=1,imx1
if (gindx(i)==1) cycle ! no futher slope failure in cells with failure during past time steps
if (eroindx(i)==1) cycle ! no futher slope failure in eroded cells
if(slo(i)<slomin) then ! default values for gently sloping cells
fsmin(i)=finf+1
zfmin(i)=ltstar(i)
pmin(i)=0.
cycle
end if
if (ct(zo(i))>1.e6) then
fsmin(i)=finf+1
zfmin(i)=ltstar(i)
pmin(i)=0.
cycle
end if
inflowh=0.
inflowrho=0.
fmrho=0.
rksb=ksb(zo(i))
rkst=kst(zo(i))
qbb=fave(i)/rksb
qbt=fave(i)/rkst
fft=tan(phit(zo(i)))/tan(slo(i))
ffb=tan(phib(zo(i)))/tan(slo(i))
!parameters of time
dtp=alphab(zo(i))*rksb*dtn/(thsatb(zo(i))-thresib(zo(i)))*cos(slo(i))**2.
if (dtp>deltadzt(i,1)**2./2. .or. dtp>deltadzb(i,1)**2./2.) then
pause 'Decrease time step!'
stop '0'
end if
!vb=dtp/dzb(i,j)
!rb=dtp/dzb(i,j)**2.
!vt=dtp/dzt(i,j)
!rt=dtp/dzt(i,j)**2.
! initiating unit weight
uwsum=0.
uwspt=uwst(zo(i)) ! allocate (uwspt(nzst+1),uwspb(nzsb+1))
uwspb=uwsb(zo(i)) ! allocate (uwspt(nzst+1),uwspb(nzsb+1))
uwt1=0.
! numercial loop for nt times
do m=1,nt
! bottom layer
do j=2,nzsb
! kkb2(i,j)=(vb/2.+rb)*kkb1(i,j+1)+(1-2*rb)*kkb1(i,j)+(rb-vb/2.)*kkb1(i,j-1)
kkb2(i,j)=kkb1(i,j)+dtp/(deltadzb(i,j)+deltadzb(i,j-1))*(kkb1(i,j+1)-kkb1(i,j-1))+&
& dtp/(deltadzb(i,j-1)*deltadzb(i,j)**2.)*(deltadzb(i,j-1)*kkb1(i,j+1)-(deltadzb(i,j)+deltadzb(i,j-1))*kkb1(i,j)+deltadzb(i,j)*kkb1(i,j-1))
end do
! top layer, with beta
do j=2,nzst
! kkt2(i,j)=(vt/2.+rt)/beta(i)*kkt1(i,j+1)+(beta(i)-2*rt)/beta(i)*kkt1(i,j)+(rt-vt/2.)/beta(i)*kkt1(i,j-1)
kkt2(i,j)=kkt1(i,j)+dtp/beta(i)/(deltadzt(i,j)+deltadzt(i,j-1))*(kkt1(i,j+1)-kkt1(i,j-1))+&
& dtp/beta(i)/(deltadzt(i,j-1)*deltadzt(i,j)**2.)*(deltadzt(i,j-1)*kkt1(i,j+1)-(deltadzt(i,j)+deltadzt(i,j-1))*kkt1(i,j)+deltadzt(i,j)*kkt1(i,j-1))
end do
! Surface point
! kkt2(i,nzst+1)=(qbt*dzt(i,j)+kkt2(i,nzst))/(dzt(i,j)+1)
kkt2(i,1)=(qbt*deltadzt(i,1)+kkt2(i,2))/(deltadzt(i,1)+1.)
! interface of the two layers
! kkt2(i,1)=(rksb/dzb(i,j)*kkb2(i,nzsb)+rkst/dzt(i,j)*kkt2(i,2))/(rksb/dzb(i,j)*(1+dzb(i,j))+rkst/dzt(i,j)*(1-dzt(i,j)))
kkt2(i,nzst+1)=(rksb/deltadzb(i,1)*kkb2(i,2)+rkst/deltadzt(i,nzst)*kkt2(i,nzst))/(rksb/deltadzb(i,1)*(1.+deltadzb(i,1))+rkst/deltadzt(i,nzst)*(1.-deltadzt(i,nzst)))
kkb2(i,1)=kkt2(i,nzst+1)
! kkt1=kkt2
! kkb1=kkb2
end do
!1000 continue
fmn=10.
flm(:)=10.
! top layer, find the minimum factor of safety among all sublayers.
do 1100, j=1,nzst+1
pt2(i,j)=1./alphat(zo(i))*dlog(kkt2(i,j))
thzt(i,j)=thresit(zo(i))+(thsatt(zo(i))-thresit(zo(i)))*exp(alphat(zo(i))*pt2(i,j))
desatt(i,j)=thzt(i,j)/thsatt(zo(i))
! The unit weight of top soil at depth (i-1)*deltazt at time t0, for cell i.
if (pt2(i,j)<-1.d0/alphat(zo(i))) then
uwt1=(uwst(zo(i))/uww-thsatt(zo(i))+thzt(i,j))*uww
else
uwt1=uwst(zo(i))
end if
uwsum=uwsum+uwt1
uwspt(j)=uwsum/float(j) ! average unit weight of top i layers soil at time t0.
if (ltstar(i)-zt(i,j)>zmin) then
fsc=ct(zo(i))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
else
fsc=ct(zo(i))/uwspt(j)/(ltstar(i)-zt(i,j)+zmin)/sin(slo(i))/cos(slo(i))
end if
! compute factor of safety for top layer
if (ltstar(i)-zt(i,j)>zmin) then
if (pt2(i,j)<0.) then
fsw=-pt2(i,j)*uww*tan(phibt(zo(i)))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
else
fsw=-pt2(i,j)*uww*tan(phit(zo(i)))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
end if
fs=fft+fsc+fsw
else
fs=finf
end if
! frictional strength cannot be less then zero
if (fs<fsc) fs=fsc
if (fs>finf) fs=finf
if (ltstar(i)-zt(i,j)<=zmin) fs=finf
! minimum factor of safety, from Upsidedown
if (abs((inidesatt(i,j)-desatt(i,j))/inidesatt(i,j))>0.05) then
zfmin(i)=zt(i,j)
pmin(i)=pt2(i,j)
fdepth(i)=ltstar(i)-zt(i,j)
flm(j)=fs
jmin=j
end if
! end do j=1,nzst+1
1100 continue
fmn=minval(flm(:))
if (fdepth(i)==0.) then
zfmin(i)=ltstar(i)
pmin(i)=pt2(i,1)
fmn=finf
end if
fsmin(i)=fmn
! if slope failure occurs, there comes material entrainment
if (fsmin(i)<1.) then
gindx(i)=1
tempfsh(i)=fdepth(i)
tempfsrho(i)=(rhos-rhow)*cvstar+rhow
end if
! end do i=1,imx1
end do
!$omp end parallel do
end subroutine doublelayer