【導(dǎo)讀】的解決實(shí)際問(wèn)題的能力。的計(jì)算程序，并以某個(gè)水庫(kù)為實(shí)例，進(jìn)行水庫(kù)泥沙淤積計(jì)算。算模式為不飽和輸沙，水流泥沙方程采用非耦合解。水流挾沙力公式采用張瑞瑾公式，推移質(zhì)輸沙率公式采用Mayer-_Peter公式，MAYER-PETER公式中的能坡J按均勻流曼寧公式近似計(jì)算。由方程（4-5）可得????

　　

【正文】 n) IF(ndisp==1)WRITE(*,*)39。INTO SGB39。 DO I=1,nn 一維泥沙淤積計(jì)算 ~ 22 ~ H=A(I)/B(I) U=Q(I)/A(I) VA1=U** VA2=*H*W VA3=(VA1/VA2)** Sx(I)=VA3* !懸移質(zhì)挾沙力 !推移質(zhì)輸沙率，斷面輸沙率 Gb(I)=FGb(rough(I),H,Q(I),A(I),,0)*B(i) ENDDO END FUNCTION FGb(rn,r,H,Q,A,rs,d,rou,rous,ndisp) IF(ndisp==1)WRITE(*,*)39。INTO FGb39。 A1=(rn**)*(Q**) A2=(A**)*(H**(4/)) RJ=A1/A2 rnn=(d**(1/))/ AA=((rnn/rn)**)*rs*H*RJ BB=*(rsr)*d CC=((rousrou)/rous)**(rou**)* IF((AABB)0)THEN FGB=0 ELSE FGB=((AABB)**)/CC ENDIF RETURN END 一維泥沙淤積計(jì)算 ~ 23 ~ !***斷面沖淤厚度計(jì)算 *** SUBROUTINE DEPO(nn,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alow0,npoint,faidy,ndisp) DIMENSION x(mm,2,nn),npoint(nn),B(nn),Gb(nn),dy(nn),dx(nn),alow(nn),alow0(nn),Q(nn),Sx(nn),S(nn) IF(ndisp==1)WRITE(*,*)39。INTO DEPO39。 rou0= dy(1)=*w*(S(1)Sx(1))*dt/rou0 !進(jìn)口斷面河床變形方程單點(diǎn)離散 !**上游第二斷面開始離散求沖淤量 ** DO J=2,NN qq=Q(J1)/((B(J1)+B(J))*) FAC=(*w*dx(j1))/qq SA=(S(J1)Sx(J1))*EXP(FAC) SB=(Sx(J1)Sx(J))*(1EXP(FAC))/FAC S(J)=Sx(J)+SASB dGb=(Gb(j)Gb(j1))/dx(j1) dQS=(Q(J)*S(J)Q(J1)*S(J1))/dx(j1) dy1=(dGb+dQS)*dt/rou0 dy2=B(J1)*dy(J1)*(1faidy) dy3=B(J)*faidy dy(J)=(dy1+dy2)/dy3 !****判斷庫(kù)尾沖刷問(wèn)題 **** alow(J)=x(1,2,J) DO L=2,npoint(J) if(alow(J)x(L,2,J)) alow(J)=x(L,2,J) ENDDO IF((alow(J)+dy(J))alow0(J))THEN dy(J)=0 temp=rou0*B(j1)*dy(J1)*(1faidy)/dt S(J)=(temp+dGb)*dx(J1)/Q(J1)+S(J1) 一維泥沙淤積計(jì)算 ~ 24 ~ ENDIF ENDDO END !*******檢驗(yàn)?zāi)嗌呈欠袷睾?********* SUBROUTINE CHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp) DIMENSION x(60,2,npxt),npoint(npxt),B(npxt),Gb(npxt),dy(npxt),dx(npxt),Q(npxt),S(npxt),Zlevel(npxt) IF(ndisp==1) WRITE(*,*)39。INTO CHECK39。 DVV= DO M=2,npxt DVV=DVV+(B(M1)*dy(M1)+B(M)*dy(M))*dx(M1)/ TOL1=DVV*1300 TOL2=(Q(1)*(S(1)S(M))+Gb(1)Gb(M))*dt IF(ABS(TOL1TOL2)ABS(*TOL1))THEN WRITE(*,*)39。MASS NOT EQUAL39。,M,TOL1,TOL2 PAUSE ENDIF ENDDO TOL1=DVV*1300/1000. TOL2=(Q(1)*(S(1)S(npxt))+(Gb(1)Gb(npxt)))*dt/1000. IF(ABS(TOL1TOL2)ABS(*TOL1))THEN WRITE(*,*)39。MASS NOT EQUAL39。,M,TOL1,TOL2 PAUSE ENDIF DO N=1,npxt DO K=1,npoint(N) if(X(k,2,n)Zlevel(N))THEN x(K,2,N)=x(K,2,N)+dy(N) 一維泥沙淤積計(jì)算 ~ 25 ~ IF(X(K,2,N)Zlevel(N))THEN x(K,2,N)=Zlevel(N) ENDIF ENDIF ENDDO ENDDO END !************文件輸出函數(shù) ************* SUBROUTINE FILE(n2,n3,n4,n5,n6,nyear,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dWQST,AVQ,dWQS,dWGb,ndisp) dimension x(mm,2,nn),npoint(nn),Zlevel(nn),alow(nn),dxa(nn),dx(nn) IF(ndisp==1)WRITE(*,*)39。INTO FILE39。 IF(nyear==0)then write(12,90) 90 format(1x,39。年 39。,2x,39。懸移質(zhì)輸沙量 (萬(wàn)噸 )39。,1x,39。年均流量 (萬(wàn)噸 )39。,1x,39?？傆俜e量 (億立方米 )39。) write(13,91) 91 format(1x,39。年 39。,6x,39。距壩里程 (m)39。,3x,39。深泓 (m)39。) write(14,92) 92 format(1x,39。年 39。,6x,39。距壩里程 (m)39。,3x,39。水面線 (m)39。) write(15,93) 93 format(1x,39。年 39。,6x,39。起點(diǎn)距 (m)39。,3x,39。高程 (m)39。) write(16,94) 94 format(1x,39。年 39。,7x,39。水位 (m)39。,7x,39。庫(kù)容 (億立方米 )39。) ENDIF IF(MOD(nyear,2)==0)then npxt=nn !***壩前斷面 *** DO k=1,npoint(15) write(n5,100)nyear,x(k,1,15),x(k,2,15) 一維泥沙淤積計(jì)算 ~ 26 ~ 100 format(1x,I2,4X,4X,) ENDDO !***斷面深泓 *** DO K=1,npxt alow(k)=x(1,2,k) DO l=2,npoint(k) if(alow(k)x(l,2,k))then alow(k)=x(l,2,k) endif ENDDO WRITE(N3,101)nyear,dxa(k),alow(k) 101 format(1x,I2,6X,5X,) ENDDO !***水面線 *** DO K=1,nn write(n4,102)nyear,dxa(k),zlevel(k) 102 format(1x,I2,6X,5X,) ENDDO !每?jī)赡晁粠?kù)容關(guān)系 DO k=273,240,2 ZL=K* CALL FZV(npoint,nn,x,ZL,ZV,dx,ndisp) WRITE(n6,103)nyear,ZL,ZV 103 format(1x,I2,4X,4X,) ENDDO ENDIF !****淤積總量 **** WRITE(n2,104)nyear,dWQST,AVQ,dWQS+dWGb 104 format(1x,I2,4X,4X,4X,) END 一維泥沙淤積計(jì)算 ~ 27 ~ !*******庫(kù)容計(jì)算 ******* SUBROUTINE FZV(NPOINT,NSEC,X,ZLEVEL,ZV,dx,ndisp) dimension X(60,2,NSEC),NPOINT(NSEC),AA(NSEC),dx(NSEC) IF(ndisp==1)WRITE(*,*)39。INTO FZV39。 ZV= DO I=1,NSEC CALL AREA(NPOINT(I),X(1,1,I),X(1,2,I),B,A,Xw,ZLEVEL,ndisp) AA(I)=A ENDDO DO I=1,NSEC1 ZV=ZV+(AA(I)+AA(I+1))*dx(I)* ENDDO END