SUBROUTINE MAIN_UFRN(nriver ,& !河道数量 nsect ,& !河道计算断面数最大值 krc ,& !河道汊点数最大值 ndata ,& !计算时刻数 Ns ,& !河道断面数 Pc ,& !河道汊点数 Nrc ,& !河道汊点微段号 Lc ,& !汊点类型 Dric ,& !流动方向 XJ ,& !流量过程 Asave ,& !水库面积 ds ,& !河道几何数据 bd ,& !梯形断面底宽 zd ,& !梯形断面底高程 sm ,& !梯形断面边坡系数 rough ,& !梯形断面糙率 ZZ0 ,& !计算区域初始水位 Zctr ,& !河道末端挡潮闸运行控制水位 Qp ,& !计算区域入河流量过程 UB1 ,& !上游边界条件类型1(1——水位边界,2——流量边界,3——水位流量关系) UB2 ,& !上游边界条件类型2(1——河网外边界,2——河网内边界) DB1 ,& !下游边界条件类型1(1——水位边界,2——流量边界,3——水位流量关系) DB2 ,& !下游边界条件类型2(1——河网外边界,2——河网内边界) UBV ,& !上游外边界断面处的水位,或流量过程 NUB ,& !上游内边界断面所连接的河道编号、微段编号 DBV ,& !下游外边界断面处的水位,或流量,或闸下潮位过程 Gate ,& !Gate(1,i) = 闸净宽;Gate(2,i) = 最大过流量;Gate(3,i) = 流量公式C1;Gate(4,i) = 流量公式C2 (流量公式Q=C1*B*e*dZ^C2中的C1和C2) NDB ,& !下游内边界断面所连接的河道编号、微段编号 Aphi ,& !下游流量内边界(管涵)的过流断面面积、流速系数 period ,& !计算小时数 dt ,& !计算时间步长 sita ,& !水量隐格式权系数 sorz ,& !水位迭代松弛因子 sorq ,& !流量迭代松弛因子 epsz ,& !水位迭代控制精度 epsq ,& !流量迭代控制精度 Bsor1 ,& !水闸外边界计算中的松弛因子 Bsor2 ,& !管涵内边界计算中的松弛因子 Z ,& !输出河道水位 Q )& !输出河道流量 BIND(C, NAME="MAIN_UFRN") !DEC$ ATTRIBUTES DLLEXPORT::MAIN_UFRN IMPLICIT NONE INTEGER::nriver ! 河道数量 INTEGER::nsect ! 河道计算断面数最大值 INTEGER::krc ! 河道汊点数最大值 INTEGER::ndata ! 计算时刻数(时间步数据点数量) INTEGER::MRIVER integer River,tstep INTEGER::I INTEGER::II ! 河道组态数据 INTEGER::Ns(nriver) ! 各河道断面数(维度: NRIVER) INTEGER::Pc(nriver) ! 各河道汊点数(维度: NRIVER) INTEGER::Nrc(krc,nriver) ! 汊点类型(1-普通,2-集中入流,3-调蓄水库)(维度: KRC, NRIVER) INTEGER::Lc(krc,nriver) ! 汊点连接的微段编号(维度: KRC, NRIVER) !边界条件 INTEGER::UB1(nriver) ! 上游边界类型(1-水位,2-流量,3-水位流量关系)(维度: NRIVER) INTEGER::UB2(nriver) ! 上游边界位置(1-外边界,2-内边界)(维度: NRIVER) INTEGER::DB1(nriver) ! 下游边界类型(同上)(维度: NRIVER) INTEGER::DB2(nriver) ! 下游边界位置(同上)(维度: NRIVER) INTEGER::NUB(2,nriver) ! 上游内边界连接的河道和微段(维度: 2, NRIVER) INTEGER::NDB(2,nriver) ! 下游内边界连接的河道和微段(维度: 2, NRIVER) REAL::ZZ0 ! 计算区域初始水位(m) REAL::Zctr ! 河道末端挡潮闸控制水位(m) REAL::period ! 计算总时长(小时) REAL::dt ! 计算时间步长(秒) REAL::sita ! 隐格式权系数(0.5~1.0) REAL::sorz ! 水位迭代松弛因子(0~1) REAL::sorq ! 流量迭代松弛因子(0~1) REAL::epsz ! 水位迭代收敛精度(m) REAL::epsq ! 流量迭代收敛精度(m?/s) REAL::Bsor1 ! 水闸边界松弛因子 REAL::Bsor2 ! 管涵边界松弛因子 REAL::ql REAL::condu,condd(3) REAL::Bs,As,Rs,Cs !河道汊点数据 REAL::Dric(krc,nriver) ! 汊点流动方向(流入1,流出-1)(维度: KRC, NRIVER) REAL::Qj(ndata,krc,nriver) ! 集中入流流量过程(维度: NDATA, KRC, NRIVER) REAL::Asave(krc,nriver) ! 调蓄水库面积(m?)(维度: KRC, NRIVER) !河道几何数据 REAL::ds(nsect - 1,nriver) ! 微段长度(m)(维度: nsect - 1, NRIVER) REAL::bd(nsect,nriver) ! 断面底宽(m)(维度: NSECT, NRIVER) REAL::zd(nsect,nriver) ! 断面底高程(m)(维度: NSECT, NRIVER) REAL::sm(nsect,nriver) ! 边坡系数(-)(维度: NSECT, NRIVER) REAL::rough(nsect,nriver) ! 糙率(曼宁系数)(维度: NSECT, NRIVER) !边界过程数据 REAL::UBV(ndata,nriver) ! 上游外边界水位/流量过程(维度: NDATA, NRIVER) REAL::DBV(ndata,nriver) ! 下游外边界水位/流量过程(维度: NDATA, NRIVER) REAL::Gate(4,nriver) ! 闸门参数(净宽、最大流量、C1、C2)(维度: 4, NRIVER) REAL::Aphi(2,nriver) ! 管涵过流参数(面积、流速系数)(维度: 2, NRIVER) REAL::Qp(ndata) ! 计算区域入河总流量过程(m?/s)(维度: NDATA) REAL::Z0(nsect,nriver) REAL::Q0(nsect,nriver) REAL::Z(nsect,nriver) REAL::Q(nsect,nriver) REAL::V(nsect,nriver) REAL::Zc(nsect,nriver),Qc(nsect,nriver) REAL::XJ(ndata,nriver*krc) REAL::maxtstep REAL::maxiter REAL::ttime REAL::fc REAL::Scanal REAL::dzmax REAL::dqmax REAL::dz REAL::dq REAL::maxz_r REAL::maxz_s REAL::maxq_r REAL::maxq_s INTEGER::Is INTEGER::iter INTEGER::J,K MRIVER = nsect - 1 Scanal = 0.0 ! Scanal——河网所有河道的合计长度,用于计算旁侧入流量ql DO I = 1,ndata DO J = 1,KRC DO K = 1,NRIVER Qj(i,j,k) = xj(i,j*k) end do end do end do DO J = 1,nriver DO I = 1,Ns(j)-1 Scanal = Scanal + ds(i,j) ! 每条河道各微段的长度,普通汊点、集中入流微段及调蓄水库的长度取0,每条河道的首末微段长度必须取一个非零值。 END DO END DO ! 计算时间步长 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 DO tstep=1,maxtstep ! 步数累加器 ttime=dt*tstep/3600.0 ! 真实时刻 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 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 exit 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 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 END SUBROUTINE MAIN_UFRN