! 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 ) ! 计算断面输出水位 ! 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