SUB_SHUILIANGJISUAN/Fortran/SUB_QZ.f90

!	==================================================================================
!	 用四点隐格式计算未知时层水位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),&
		sm(nsect,nriver), rough(nsect,nriver)
    ! boundary
	INTEGER::UB1(nriver),UB2(nriver),DB1(nriver),DB2(nriver),&
		NUB(2,nriver),NDB(2,nriver)
    REAL::UBV(ndata,nriver),Aphi(2,nriver),DBV(ndata,nriver),&
		Gate(4,nriver),ql,Zctr
    ! calcu
	REAL::dt,sita,Bsor1,Bsor2
    ! zzqq
    REAL::Z0(nsect,nriver),Q0(nsect,nriver),Z(nsect,nriver),&
		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