zzx/MAIN_UFRN
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
MAIN_UFRN/Fortran/SUB_QZ.f90

350 lines
7.8 KiB
Fortran
Raw Permalink Normal View History

1
2025-05-14 17:42:10 +08:00
! ==================================================================================
! <20><><EFBFBD>ĵ<EFBFBD><C4B5><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD><EFBFBD><EFBFBD>δ֪ʱ<D6AA><CAB1>ˮλZ<CEBB><5A><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Q<EFBFBD><51><EFBFBD>ӳ<EFBFBD><D3B3><EFBFBD>
! ==================================================================================
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
! ----------------------------------------------------------------------------------
! <09><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƹ<EFBFBD>ʽϵ<CABD><CFB5>L<EFBFBD><4C>M<EFBFBD><4D>P<EFBFBD><50>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
! ------------------------------------------------------
! <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӷμ<D3B6><CEBC><EFBFBD>
!
! <09><><EFBFBD><EFBFBD><EFBFBD>Ӷ<EFBFBD>M<EFBFBD>㴦<EFBFBD><E3B4A6>ˮ<EFBFBD><CBAE><EFBFBD><EFBFBD>B<EFBFBD><42><EFBFBD><EFBFBD><EFBFBD><EFBFBD>A<EFBFBD><41>ˮ<EFBFBD><CBAE><EFBFBD>뾶R<EBBEB6><52><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Q<EFBFBD><51>ˮλ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))
! ------------------------------------------------------
! <09><><EFBFBD><EFBFBD><EFBFBD>Ӷη<D3B6><CEB7><EFBFBD><EFBFBD>е<EFBFBD>ϵ<EFBFBD><CFB5>
! ------------------------------------------------------
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
! --------------------------------------------------------
! <09><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㴦<EFBFBD><E3B4A6><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͨ<EFBFBD><CDA8><EFBFBD><EFBFBD>Lc=1<><31><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Lc=2<><32><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ˮ<EFBFBD><CBAE>Lc=3<><33>
! --------------------------------------------------------
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
! ------------------------------------------------------
! <09><><EFBFBD><EFBFBD><EFBFBD>α߽<CEB1><DFBD><EFBFBD><EFBFBD><EFBFBD>Ϊˮλ<CBAE><CEBB><EFBFBD><EFBFBD><EFBFBD>ߵĵ<DFB5><C4B5><EFBFBD>ϵ<EFBFBD><CFB5>
! ------------------------------------------------------
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
! ------------------------------------------------------
! <09><><EFBFBD><EFBFBD><EFBFBD>α߽<CEB1><DFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ߵĵ<DFB5><C4B5><EFBFBD>ϵ<EFBFBD><CFB5>
! ------------------------------------------------------
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
! ----------------------------------------------------------------------------------
! <09>ش<EFBFBD><D8B4><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Q<EFBFBD><51>ˮλ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
Reference in New Issue Copy Permalink