! ================================================================================== ! 用四点隐格式计算未知时层水位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