350 lines
7.8 KiB
Fortran
350 lines
7.8 KiB
Fortran
! ==================================================================================
|
||
! 用四点隐格式计算未知时层水位Z和流量Q的子程序
|
||
! ==================================================================================
|
||
subroutine sub_QZ( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river ,&
|
||
tstep ,&
|
||
Ns ,&
|
||
Pc ,&
|
||
Nrc ,&
|
||
Lc ,&
|
||
Dric ,&
|
||
Qj ,&
|
||
Asave ,&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
UB1 ,&
|
||
UB2 ,&
|
||
DB1 ,&
|
||
DB2 ,&
|
||
NUB ,&
|
||
NDB ,&
|
||
UBV ,&
|
||
Aphi ,&
|
||
DBV ,&
|
||
Gate ,&
|
||
ql ,&
|
||
Zctr ,&
|
||
dt ,&
|
||
sita ,&
|
||
Bsor1 ,&
|
||
Bsor2 ,&
|
||
Z0 ,&
|
||
Q0 ,&
|
||
Z ,&
|
||
Q ,&
|
||
V ,&
|
||
condu ,&
|
||
condd ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
|
||
|
||
INTEGER::NRIVER
|
||
INTEGER::NSECT
|
||
INTEGER::MRIVER
|
||
INTEGER::KRC
|
||
INTEGER::NDATA
|
||
integer River,tstep
|
||
REAL Qjj(ndata)
|
||
real L(nsect),M(nsect),P(nsect),R(nsect)
|
||
|
||
! node
|
||
INTEGER::Ns(nriver),Pc(nriver),Nrc(krc,nriver),Lc(krc,nriver)
|
||
REAL::Dric(krc,nriver),Qj(ndata,krc,nriver),Asave(krc,nriver)
|
||
! section
|
||
REAL::ds(mriver,nriver),bd(nsect,nriver),zd(nsect,nriver)
|
||
REAL::sm(nsect,nriver), rough(nsect,nriver)
|
||
! boundary
|
||
INTEGER::UB1(nriver),UB2(nriver),DB1(nriver),DB2(nriver)
|
||
INTEGER::NUB(2,nriver),NDB(2,nriver)
|
||
REAL::UBV(ndata,nriver),Aphi(2,nriver),DBV(ndata,nriver)
|
||
REAL::Gate(4,nriver),ql,Zctr
|
||
! calcu
|
||
REAL::dt,sita,Bsor1,Bsor2
|
||
! zzqq
|
||
REAL::Z0(nsect,nriver),Q0(nsect,nriver),Z(nsect,nriver)
|
||
REAL::Q(nsect,nriver),V(nsect,nriver)
|
||
! r_bv
|
||
REAL::condu,condd(3)
|
||
! BARC
|
||
REAL::Bs,As,Rs,Cs
|
||
|
||
|
||
f(t10,t1,t20,t2)=0.5*(sita*(t2+t1)+(1.0-sita)*(t20+t10))
|
||
tc=dt*float(tstep)/3600.0
|
||
|
||
if(UB1(river).eq.1)Z(1,river)=condu
|
||
if(UB1(river).eq.2)Q(1,river)=condu
|
||
|
||
! ----------------------------------------------------------------------------------
|
||
! 计算递推公式系数L、M、P、R
|
||
! ----------------------------------------------------------------------------------
|
||
P(1)=condu
|
||
R(1)=0.0
|
||
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,1,Z0(1,river),&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
BB10=Bs
|
||
AA10=As
|
||
RR10=Rs
|
||
CC10=Cs
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,1,Z(1,river),&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
BB1=Bs
|
||
AA1=As
|
||
RR1=Rs
|
||
CC1=Cs
|
||
do 100 k=1,Ns(river)-1
|
||
Ik=0
|
||
do 20 i=1,Pc(river)
|
||
if(Nrc(i,river).eq.k)then
|
||
Ic=i
|
||
Ik=1
|
||
end if
|
||
20 continue
|
||
|
||
if(Ik.ne.1)then
|
||
! ------------------------------------------------------
|
||
! 非汊点河段计算
|
||
!
|
||
! 计算河段M点处的水面宽B、面积A、水力半径R、流量Q和水位Z
|
||
! ------------------------------------------------------
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,k+1,Z0(k+1,river),&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
BB20=Bs
|
||
AA20=As
|
||
RR20=Rs
|
||
CC20=Cs
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,k+1,Z(k+1,river),&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
BB2=Bs
|
||
AA2=As
|
||
RR2=Rs
|
||
CC2=Cs
|
||
Bm=f(BB10,BB1,BB20,BB2)
|
||
Am=f(AA10,AA1,AA20,AA2)
|
||
Rm=f(RR10,RR1,RR20,RR2)
|
||
Cm=f(CC10,CC1,CC20,CC2)
|
||
Qm=f(Q0(k,river),Q(k,river),Q0(k+1,river),Q(k+1,river))
|
||
Zm=f(Z0(k,river),Z(k,river),Z0(k+1,river),Z(k+1,river))
|
||
|
||
! ------------------------------------------------------
|
||
! 计算河段方程中的系数
|
||
! ------------------------------------------------------
|
||
a1=1.0
|
||
c1=2.0*sita*dt/(ds(k,river)*Bm)
|
||
e1=Z0(k,river)+Z0(k+1,river)+(1.0-sita)/sita*c1*(Q0(k,river)&
|
||
-Q0(k+1,river))+2.0*dt*ql/Bm
|
||
a2=2.0*sita*dt/ds(k,river)*((Qm/Am)**2.0*Bm-9.8*Am)
|
||
c2=1.0-4.0*sita*dt/ds(k,river)*Qm/Am
|
||
d2=1.0+4.0*sita*dt/ds(k,river)*Qm/Am
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,k,Zm,&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
AZm1=As
|
||
call sub_sect( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,k+1,Zm,&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
AZm2=As
|
||
e2=(1.0-sita)/sita*a2*(Z0(k+1,river)-Z0(k,river))+(1.0-4.0*&
|
||
(1.0-sita)*dt/ds(k,river)*Qm/Am)*Q0(k+1,river)+(1.0+4.0*&
|
||
(1.0-sita)*dt/ds(k,river)*Qm/Am)*Q0(k,river)+2.0*dt*&
|
||
(Qm/Am)**2.0*(AZm2-AZm1)/ds(k,river)-2.0*dt*9.8*abs(Qm)*Qm&
|
||
/(Am*Cm*Cm*Rm)
|
||
else
|
||
|
||
! --------------------------------------------------------
|
||
! 汊点计算处理（普通汊点Lc=1，集中入流Lc=2，调蓄水库Lc=3）
|
||
! --------------------------------------------------------
|
||
if(Lc(Ic,river).eq.1)then
|
||
do 35 j=1,nriver
|
||
if(NUB(1,j).eq.river.and.NUB(2,j).eq.Nrc(Ic,river))then
|
||
a1=0.0
|
||
c1=1.0
|
||
e1=Dric(Ic,river)*Q(1,j)
|
||
a2=1.0
|
||
c2=0.0
|
||
d2=0.0
|
||
e2=0.0
|
||
elseif(NDB(1,j).eq.river.and.NDB(2,j).eq.Nrc(Ic,river))then
|
||
a1=0.0
|
||
c1=1.0
|
||
e1=Dric(Ic,river)*Q(Ns(j),j)
|
||
a2=1.0
|
||
c2=0.0
|
||
d2=0.0
|
||
e2=0.0
|
||
end if
|
||
35 continue
|
||
elseif(Lc(Ic,river).eq.2)then
|
||
do 55 ii=1,ndata
|
||
Qjj(ii)=Qj(ii,Ic,river)
|
||
55 continue
|
||
call int_a(tc,ndata,Qjj,fc)
|
||
a1=0.0
|
||
c1=1.0
|
||
e1=Dric(Ic,river)*fc
|
||
a2=1.0
|
||
c2=0.0
|
||
d2=0.0
|
||
e2=0.0
|
||
elseif(Lc(Ic,river).eq.3)then
|
||
a1=0.5*Asave(Ic,river)/dt
|
||
c1=1.0
|
||
e1=Asave(Ic,river)*Z0(k,river)/dt
|
||
a2=1.0
|
||
c2=0.0
|
||
d2=0.0
|
||
e2=0.0
|
||
end if
|
||
|
||
|
||
end if
|
||
|
||
! ------------------------------------------------------
|
||
! 求上游边界条件为水位过程线的递推系数
|
||
! ------------------------------------------------------
|
||
if(UB1(river).eq.1)then
|
||
Y1=a1*R(k)-c1
|
||
Y2=a2*R(k)+c2
|
||
Y3=e1-a1*P(k)
|
||
Y4=e2-a2*P(k)
|
||
Y5=a1*Y2+a2*Y1
|
||
L(k+1)=(a2*Y3+a1*Y4)/Y5
|
||
M(k+1)=-(a1*d2+a2*c1)/Y5
|
||
P(k+1)=(Y2*Y3-Y1*Y4)/Y5
|
||
R(k+1)=(d2*Y1-c1*Y2)/Y5
|
||
end if
|
||
! ------------------------------------------------------
|
||
! 求上游边界条件为流量过程线的递推系数
|
||
! ------------------------------------------------------
|
||
if(UB1(river).eq.2)then
|
||
Y1=a1-c1*R(k)
|
||
Y2=a2+c2*R(k)
|
||
Y3=e1+c1*P(k)
|
||
Y4=e2-c2*P(k)
|
||
Y5=d2*Y1-c1*Y2
|
||
L(k+1)=(d2*Y3-c1*Y4)/Y5
|
||
M(k+1)=-(a1*d2+a2*c1)/Y5
|
||
P(k+1)=(Y1*Y4-Y2*Y3)/Y5
|
||
R(k+1)=(a1*Y2+a2*Y1)/Y5
|
||
end if
|
||
|
||
BB10=BB20
|
||
AA10=AA20
|
||
RR10=RR20
|
||
CC10=CC20
|
||
BB1=BB2
|
||
AA1=AA2
|
||
RR1=RR2
|
||
CC1=CC2
|
||
100 continue
|
||
! ----------------------------------------------------------------------------------
|
||
! 回代求出流量Q和水位Z
|
||
! ----------------------------------------------------------------------------------
|
||
if(UB1(river).eq.1)then
|
||
Q(Ns(river),river)=(condd(3)-condd(1)*P(Ns(river)))&
|
||
/(condd(1)*R(Ns(river))+condd(2))
|
||
do 120 i=Ns(river),2,-1
|
||
Z(i,river)=P(i)+R(i)*Q(i,river)
|
||
Q(i-1,river)=L(i)+M(i)*Q(i,river)
|
||
120 continue
|
||
end if
|
||
|
||
if(UB1(river).eq.2)then
|
||
Z(Ns(river),river)=(condd(3)-condd(2)*P(Ns(river)))&
|
||
/(condd(2)*R(Ns(river))+condd(1))
|
||
do 130 i=Ns(river),2,-1
|
||
Q(i,river)=P(i)+R(i)*Z(i,river)
|
||
Z(i-1,river)=L(i)+M(i)*Z(i,river)
|
||
130 continue
|
||
end if
|
||
|
||
return
|
||
end |