571 lines
15 KiB
Fortran
571 lines
15 KiB
Fortran
! SUB_SHUILIANGJISUAN.f90
|
||
!
|
||
! FUNCTIONS/SUBROUTINES exported from SUB_SHUILIANGJISUAN.dll:
|
||
! SUB_SHUILIANGJISUAN - subroutine
|
||
!
|
||
subroutine SUB_SHUILIANGJISUAN( NDATA ,& ! 计算时刻数 // 输入参数
|
||
LONG ,& ! 河道长度
|
||
INbd ,& ! 断面底宽
|
||
INzd ,& ! 底高程
|
||
INsm ,& ! 边坡系数
|
||
INrough ,& ! 糙率
|
||
ZZ0 ,& ! 计算区域初始水位
|
||
DB ,& ! 下游边界条件类型 1.水位边界 2.流量边界
|
||
UQ ,& ! 上游流量过程
|
||
DQH ,& ! 上游水位/流量过程
|
||
dt ,& ! 计算步长（ 秒）
|
||
OUTT ,& ! 计算断面输出时刻
|
||
OUTQ ,& ! 计算断面输出流量
|
||
OUTH )& ! 计算断面输出水位
|
||
BIND(C, NAME="SUB_SHUILIANGJISUAN")
|
||
! Expose subroutine SUB_SHUILIANGJISUAN to users of this DLL
|
||
!
|
||
!DEC$ ATTRIBUTES DLLEXPORT::SUB_SHUILIANGJISUAN
|
||
|
||
! Variables
|
||
|
||
INTEGER::FHD
|
||
|
||
INTEGER::NRIVER
|
||
INTEGER::NSECT
|
||
INTEGER::MRIVER
|
||
INTEGER::KRC
|
||
INTEGER::NDATA
|
||
|
||
REAL::ZZ0
|
||
|
||
INTEGER::DB
|
||
|
||
REAL::UQ(NDATA)
|
||
REAL::DQH(NDATA)
|
||
|
||
REAL::OUTQ(1000)
|
||
REAL::OUTH(1000)
|
||
REAL::OUTT(1000)
|
||
|
||
REAL::DDGC ! 计算断面堤顶高程
|
||
|
||
REAL::MAXH
|
||
|
||
INTEGER::MY ! 1.漫溢 2.未漫溢
|
||
|
||
integer River,tstep
|
||
! node
|
||
INTEGER,ALLOCATABLE::Ns(:),Pc(:),Nrc(:,:),Lc(:,:)
|
||
REAL,ALLOCATABLE::Dric(:,:),Qj(:,:,:),Asave(:,:)
|
||
! section
|
||
REAL,ALLOCATABLE::ds(:,:),bd(:,:),zd(:,:),sm(:,:), rough(:,:)
|
||
! boundary
|
||
INTEGER,ALLOCATABLE::UB1(:),UB2(:),DB1(:),DB2(:),NUB(:,:),NDB(:,:)
|
||
REAL,ALLOCATABLE::UBV(:,:),Aphi(:,:),DBV(:,:),Gate(:,:)
|
||
REAL::ql,Zctr
|
||
! calcu
|
||
REAL::dt,sita,Bsor1,Bsor2
|
||
! zzqq
|
||
REAL,ALLOCATABLE::Z0(:,:),Q0(:,:),Z(:,:),Q(:,:),V(:,:)
|
||
! r_bv
|
||
REAL::condu,condd(3)
|
||
! BARC
|
||
REAL::Bs,As,Rs,Cs
|
||
|
||
REAL,ALLOCATABLE::Qp(:)
|
||
REAL,ALLOCATABLE::Zc(:,:),Qc(:,:)
|
||
|
||
INTEGER::I
|
||
INTEGER::II
|
||
|
||
REAL::LONG
|
||
|
||
REAL::INbd
|
||
REAL::INzd
|
||
REAL::INsm
|
||
REAL::INrough
|
||
|
||
INTEGER::TEMP
|
||
|
||
!OPEN(1,FILE='OUTFHDL.TXT')
|
||
!WRITE(1,*)' ! 计算时刻数 ', NDATA ,&
|
||
! ' ! 河道长度 ', LONG ,&
|
||
! '! 断面底宽 ', INbd ,&
|
||
! '! 底高程 ' , INzd ,&
|
||
! '! 边坡系数 ', INsm ,&
|
||
! ' ! 糙率 ' , INrough ,&
|
||
! ' ! 计算区域初始水位 ', ZZ0 ,&
|
||
! ' ! 下游边界条件类型 1.水位边界 2.流量边界 ', DB ,&
|
||
! ' ! 上游流量过程 ' , UQ ,&
|
||
! ' ! 上游水位/流量过程 ', DQH ,&
|
||
! ' ! 计算步长（ 秒） ', dt ,&
|
||
! '! 计算断面堤顶高程 ', DDGC
|
||
FHD = 0
|
||
TEMP = 0
|
||
|
||
! 1.读入河道组态数据
|
||
!open(1,file='河道组态数据.TXT')
|
||
|
||
! 读入：1.河道数量 2.河道计算断面数最大值 3.河道汊点数最大值 4.计算时刻数
|
||
!CALL SUB_GETNXT(1)
|
||
!read(1,*)
|
||
nriver = 1
|
||
nsect = 6
|
||
krc = 0
|
||
|
||
mriver=nsect - 1 ! 这是我写的，微段数不就等于断面数-1吗
|
||
|
||
ALLOCATE(Ns(nriver),Pc(nriver),Nrc(krc,nriver),Lc(krc,nriver))
|
||
ALLOCATE(Dric(krc,nriver),Qj(ndata,krc,nriver),Asave(krc,nriver))
|
||
ALLOCATE(ds(mriver,nriver),bd(nsect,nriver),zd(nsect,nriver),&
|
||
sm(nsect,nriver), rough(nsect,nriver))
|
||
ALLOCATE(UB1(nriver),UB2(nriver),DB1(nriver),DB2(nriver),&
|
||
NUB(2,nriver),NDB(2,nriver))
|
||
ALLOCATE(UBV(ndata,nriver),Aphi(2,nriver),DBV(ndata,nriver),&
|
||
Gate(4,nriver))
|
||
ALLOCATE(Z0(nsect,nriver),Q0(nsect,nriver),Z(nsect,nriver),&
|
||
Q(nsect,nriver),V(nsect,nriver))
|
||
ALLOCATE(Qp(ndata))
|
||
ALLOCATE(Zc(nsect,nriver),Qc(nsect,nriver))
|
||
|
||
!读入： 1.河道序号 2.河道断面数 3.河道汊点数
|
||
! （所有汊点均变成三汊口处理）
|
||
! CALL SUB_GETNXT(1)
|
||
! DO I = 1,nriver
|
||
!read(1,*)II,Ns(i),Pc(i)
|
||
! END DO
|
||
Ns(1) = 6
|
||
Pc(1) = 0
|
||
|
||
!读入：! 1.河道序号 2.河道汊点微段号 3.汊点类型 4.流动方向
|
||
! 汊点类型：1——普通汊点，2——集中入流，3——调蓄水库
|
||
! 与汊点相连的河道或集中入流的水流方向，流入汊点：1，流出汊点：-1，调蓄水库任意给定
|
||
Nrc = 0
|
||
Lc = 0
|
||
Dric = 0
|
||
! DO J = 1,nriver
|
||
! DO I = 1,Pc(J)
|
||
! CALL SUB_GETNXT(1)
|
||
!read(1,*)II,Nrc(i,j),Lc(i,j),Dric(i,j)
|
||
! END DO
|
||
! END DO
|
||
|
||
! 集中入流的汊点读入流量过程
|
||
Qj = 0.0
|
||
!DO k=1,nriver
|
||
! DO j=1,Pc(k)
|
||
! CALL SUB_GETNXT(1)
|
||
! if(Lc(j,k).eq.2)then
|
||
! read(1,*)II,II,(Qj(i,j,k),i=1,ndata)
|
||
! end if
|
||
! END DO
|
||
! END DO
|
||
|
||
! 连接调蓄水库的汊点读入水库面积（m2）
|
||
Asave = 0.0
|
||
!DO k=1,nriver
|
||
! DO j=1,Pc(k)
|
||
! if(Lc(j,k).eq.3)then
|
||
! CALL SUB_GETNXT(1)
|
||
! read(1,*)Asave(j,k)
|
||
! end if
|
||
! END DO
|
||
! END DO
|
||
|
||
! 2.读入河道几何数据
|
||
! CALL SUB_GETNXT(1)
|
||
! ds = 0.0
|
||
! DO J = 1,nriver
|
||
! DO I = 1,Ns(j)-1
|
||
!read(1,*)II,II,ds(i,j) ! 每条河道各微段的长度，普通汊点、集中入流微段及调蓄水库的长度取0，每条河道的首末微段长度必须取一个非零值。
|
||
! END DO
|
||
! END DO
|
||
|
||
ds(:,1) = LONG/5
|
||
|
||
! CALL SUB_GETNXT(1)
|
||
! bd = 0.0
|
||
! zd = 0.0
|
||
! sm = 0.0
|
||
! rough = 0.0
|
||
! DO j=1,nriver
|
||
! DO I =1,Ns(j)
|
||
!read(1,*)II,II,bd(i,j),zd(i,j),sm(i,j),rough(i,j) ! 梯形断面底宽、底高程、边坡系数、糙率
|
||
! END DO
|
||
! END DO
|
||
|
||
bd = INbd
|
||
zd = INzd
|
||
sm = INsm
|
||
rough = INrough
|
||
|
||
! 3.读入河道边界数据
|
||
! CALL SUB_GETNXT(1)
|
||
! read(1,*)ZZ0,Zctr ! ZZ0——计算区域初始水位,Zctr——河道末端挡潮闸运行控制水位,
|
||
|
||
Zctr = 0.0
|
||
|
||
Scanal = 0.0 ! Scanal——河网所有河道的合计长度，用于计算旁侧入流量ql
|
||
|
||
DO J = 1,nriver
|
||
DO I = 1,Ns(j)-1
|
||
Scanal = Scanal + ds(i,j) ! 每条河道各微段的长度，普通汊点、集中入流微段及调蓄水库的长度取0，每条河道的首末微段长度必须取一个非零值。
|
||
END DO
|
||
END DO
|
||
|
||
!CALL SUB_GETNXT(1)
|
||
!DO I = 1,ndata
|
||
! read(1,*)II,Qp(i) ! Qp——计算区域入河流量过程
|
||
! END DO
|
||
|
||
Qp = 0.0
|
||
|
||
! UB1上游边界条件类型1（1——水位边界，2——流量边界,3——水位流量关系）
|
||
! UB2上游边界条件类型2（1——河网外边界，2——河网内边界）
|
||
! DB1下游边界条件类型1（1——水位边界，2——流量边界,3——水位流量关系）
|
||
! DB2下游边界条件类型2（1——河网外边界，2——河网内边界）
|
||
!CALL SUB_GETNXT(1)
|
||
! DO I = 1,nriver
|
||
! read(1,*)II,UB1(i),UB2(i),DB1(i),DB2(i)
|
||
!END DO
|
||
|
||
|
||
|
||
UB1(1) = 2
|
||
UB2(1) = 1
|
||
|
||
DB1(1) = DB
|
||
DB2(1) = 1
|
||
|
||
|
||
!CALL SUB_GETNXT(1)
|
||
UBV = 0.0
|
||
NUB = 0
|
||
!DO i=1,nriver
|
||
! if(UB2(i).eq.1)then ! 如果是上游河网外部边界
|
||
! read(1,*)II,II,II,(UBV(j,i),j=1,ndata) ! 读入上游外边界断面处的水位，或流量过程
|
||
! elseif(UB2(i).eq.2)then ! 如果是上游河网内部边界
|
||
! read(1,*)II,II,II,(NUB(j,i),j=1,2) ! 读入上游内边界断面所连接的河道编号、微段编号
|
||
! endif
|
||
! END DO
|
||
|
||
DO I = 1,ndata
|
||
UBV(I,1) = UQ(I)
|
||
END DO
|
||
|
||
!CALL SUB_GETNXT(1)
|
||
!DO i=1,nriver
|
||
! if(DB2(i).eq.1)then ! 下游外部
|
||
! if(DB1(i).eq.3)then ! 如果是闸下潮位过程
|
||
! read(1,*)II,II,II,(DBV(j,i),j=1,ndata),(Gate(j,i),j=1,4) ! 读入 Gate(1,i) = 闸净宽；Gate(2,i) = 最大过流量；Gate(3,i) = 流量公式C1；Gate(4,i) = 流量公式C2 （流量公式Q=C1*B*e*dZ^C2中的C1和C2）
|
||
! ELSE
|
||
! read(1,*)II,II,II,(DBV(j,i),j=1,ndata) ! 读入下游外边界断面处的水位，或流量，或闸下潮位过程
|
||
! end if
|
||
! elseif(DB2(i).eq.2)then ! 下游内部
|
||
! if(DB1(i).eq.2)then
|
||
! read(1,*)II,II,II,(NDB(j,i),j=1,2),(Aphi(j,i),j=1,2) ! 下游流量内边界（管涵）的过流断面面积、流速系数 下游内边界一般用水位，如果用流量，一般是管涵，需要提供过流面积和流速系数
|
||
! ELSE
|
||
! read(1,*)II,II,II,(NDB(j,i),j=1,2) ! 读入下游内边界断面所连接的河道编号、微段编号
|
||
! end if
|
||
! endif
|
||
!END DO
|
||
|
||
DO I = 1,ndata
|
||
DBV(I,1) = DQH(I)
|
||
END DO
|
||
|
||
! 4.读入河道计算条件
|
||
!CALL SUB_GETNXT(1)
|
||
!read(1,*)period,dt,sita ! period——计算小时数，dt——计算时间步长（秒）,sita——水量隐格式权系数
|
||
!CALL SUB_GETNXT(1)
|
||
!read(1,*)sorz,sorq,epsz,epsq ! sorz、sorq——水位、流量迭代松弛因子,epsz、epsq——水位、流量迭代控制精度
|
||
!CALL SUB_GETNXT(1)
|
||
!read(1,*)Bsor1,Bsor2 ! Bsor1、Bsor2——水闸外边界、管涵内边界计算中的松弛因子
|
||
!close(1)
|
||
|
||
period =REAL( NDATA - 1 )
|
||
|
||
! read(1,*)dt
|
||
|
||
sita = 0.7
|
||
|
||
sorz = 0.3
|
||
sorq = 0.3
|
||
epsz = 0.01
|
||
epsq = 0.1
|
||
|
||
Bsor1 = 0.3
|
||
Bsor2 = 0.1
|
||
|
||
|
||
! 计算时间步长
|
||
maxtstep=period*3600/dt
|
||
maxiter=1000 ! 仿真模拟步数
|
||
|
||
! 给河道水位、流量赋初值
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
Z0(Is,river)=ZZ0
|
||
Q0(Is,river)=0.0
|
||
END DO
|
||
END DO
|
||
|
||
! 流量、水位输出
|
||
!open(1,file='河道1各断面水位变化.TXT')
|
||
!open(2,file='河道2各断面水位变化.TXT')
|
||
!open(3,file='河道3各断面水位变化.TXT')
|
||
!open(4,file='河道4各断面水位变化.TXT')
|
||
!
|
||
!open(11,file='河道1各断面流量变化.TXT')
|
||
!open(12,file='河道2各断面流量变化.TXT')
|
||
!open(13,file='河道3各断面流量变化.TXT')
|
||
!open(14,file='河道4各断面流量变化.TXT')
|
||
|
||
! ----------------------------------------------------------------------------------
|
||
! 迭代求解每个时刻的水位、流量
|
||
! ----------------------------------------------------------------------------------
|
||
DO tstep=1,maxtstep ! 步数累加器
|
||
|
||
TEMP = TEMP + 1
|
||
ttime=dt*tstep/3600.0 ! 真实时刻
|
||
! write(*,15)ttime
|
||
|
||
iter=0
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
Z(Is,river)=Z0(Is,river)
|
||
Q(Is,river)=Q0(Is,river)
|
||
END DO
|
||
END DO
|
||
|
||
! 计算当前时刻入河流量大小
|
||
call int_a( ttime ,& ! 真实时刻 // 输入参数
|
||
ndata ,& ! 入河流量变化过程时段数
|
||
Qp ,& ! 入河流量变化过程线
|
||
fc ) ! 当前时刻入河流量大小 // 输出参数
|
||
|
||
ql=fc/Scanal
|
||
|
||
25 iter=iter+1
|
||
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
Zc(Is,river)=Z(Is,river)
|
||
Qc(Is,river)=Q(Is,river)
|
||
END DO
|
||
END DO
|
||
|
||
DO river=1,nriver
|
||
! 确定计算时刻河段上下游边界条件值
|
||
call sub_bound( NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river ,&
|
||
tstep ,&
|
||
Ns ,&
|
||
Pc ,&
|
||
Nrc ,&
|
||
Lc ,&
|
||
Dric ,&
|
||
Qj ,&
|
||
Asave ,&
|
||
UB1 ,&
|
||
UB2 ,&
|
||
DB1 ,&
|
||
DB2 ,&
|
||
NUB ,&
|
||
NDB ,&
|
||
UBV ,&
|
||
Aphi ,&
|
||
DBV ,&
|
||
Gate ,&
|
||
ql ,&
|
||
Zctr ,&
|
||
dt ,&
|
||
sita ,&
|
||
Bsor1 ,&
|
||
Bsor2 ,&
|
||
Z0 ,&
|
||
Q0 ,&
|
||
Z ,&
|
||
Q ,&
|
||
V ,&
|
||
condu ,&
|
||
condd )
|
||
|
||
! 用四点隐格式计算未知时层水位Z和流量Q
|
||
call 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 )
|
||
|
||
END DO
|
||
|
||
dzmax=0.0
|
||
dqmax=0.0
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
dz=abs(Z(Is,river)-Zc(Is,river))
|
||
dq=abs(Q(Is,river)-Qc(Is,river))
|
||
if(dz.gt.dzmax)then
|
||
dzmax=dz
|
||
maxz_r=river
|
||
maxz_s=Is
|
||
end if
|
||
if(dq.gt.dqmax)then
|
||
dqmax=dq
|
||
maxq_r=river
|
||
maxq_s=Is
|
||
end if
|
||
END DO
|
||
END DO
|
||
|
||
if(dzmax.gt.epsz.or.dqmax.gt.epsq)then
|
||
! write(*,50)iter,maxz_r,maxz_s,dzmax,maxq_r,maxq_s,dqmax
|
||
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
Z(Is,river)=(1.0-sorz)*Zc(Is,river)+&
|
||
sorz*Z(Is,river)
|
||
Q(Is,river)=(1.0-sorq)*Qc(Is,river)+&
|
||
sorq*Q(Is,river)
|
||
END DO
|
||
END DO
|
||
|
||
if(iter.le.maxiter)then
|
||
goto 25
|
||
else
|
||
! write(*,55)
|
||
goto 600
|
||
end if
|
||
|
||
else
|
||
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
! 计算断面几何要素
|
||
call sub_sect(NRIVER ,&
|
||
NSECT ,&
|
||
MRIVER ,&
|
||
KRC ,&
|
||
NDATA ,&
|
||
river,Is,Z(Is,river),&
|
||
ds ,&
|
||
bd ,&
|
||
zd ,&
|
||
sm ,&
|
||
rough ,&
|
||
Bs ,&
|
||
As ,&
|
||
Rs ,&
|
||
Cs )
|
||
|
||
|
||
V(Is,river)=Q(Is,river)/As
|
||
END DO
|
||
END DO
|
||
|
||
! 输出几个典型河道的水位流量变化过程
|
||
DO jj=1,nriver
|
||
!write(jj,80)ttime,(Z(ii,jj),ii=1,Ns(jj))
|
||
END DO
|
||
DO jj=1,nriver
|
||
jj1=jj+10
|
||
!write(jj1,80)ttime,(Q(ii,jj),ii=1,Ns(jj))
|
||
END DO
|
||
|
||
IF(FHD==1)THEN
|
||
OUTQ(TEMP ) = Q(5,1)*0.9
|
||
OUTH(TEMP ) = Z(5,1)*0.95
|
||
ELSE
|
||
OUTQ(TEMP ) = Q(5,1)
|
||
OUTH(TEMP ) = Z(5,1)
|
||
END IF
|
||
|
||
OUTT(TEMP ) = ttime
|
||
|
||
DO river=1,nriver
|
||
DO Is=1,Ns(river)
|
||
Z0(Is,river)=Z(Is,river)
|
||
Q0(Is,river)=Q(Is,river)
|
||
END DO
|
||
END DO
|
||
|
||
end if
|
||
|
||
END DO
|
||
|
||
600 close(100)
|
||
close(2)
|
||
close(3)
|
||
close(4)
|
||
|
||
close(11)
|
||
close(12)
|
||
close(13)
|
||
close(14)
|
||
|
||
15 format(55x,'time=',f10.2)
|
||
50 format(55x,'iter=',i5/1x,'maxz_r=',i5,5x,'maxz_s=',i5,&
|
||
5x,'dzmax=',e10.4/1x,'maxq_r=',i5,5x,'maxq_s=',i5,5x,'dqmax=',e10.4)
|
||
55 format(1x,'迭代发散')
|
||
80 format(1x,40f12.4)
|
||
|
||
|
||
!MAXH = MAXVAL(OUTH)
|
||
!IF(MAXH.GT.DDGC)THEN
|
||
! MY = 1
|
||
!ELSE
|
||
! MY = 0
|
||
!END IF
|
||
|
||
!WRITE(1,*)' ! 1.漫溢 0.未漫溢 ' , MY ,&
|
||
! '! 计算断面输出时刻 ', OUTT ,&
|
||
! '! 计算断面输出流量 ', OUTQ ,&
|
||
! '! 计算断面输出水位 ', OUTH
|
||
|
||
CLOSE(1)
|
||
|
||
END SUBROUTINE SUB_SHUILIANGJISUAN
|
||
|