SUB_SHUILIANGJISUAN/Fortran/SUB_SHUILIANGJISUAN.f90

!  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