MAIN_UFRN/Fortran/MAIN_UFRN.f90

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)
    
    INTEGER::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