河流模擬課設

1、program mainparameter (nn=31,mm=80,nd=3653,ny=10,fai=1.0,ndisp=0,npxt=31,fai1=0.5)parameter (rs=2650*9.8,r=9800,w=0.0012,rou0=1325.0,dt=86400.0,d=0.002,rou=1000,rous=2650)dimension x(mm,2,nn),rough(nn),dx(nn),zlevel(nn),Q(nn),Npoint(nn),b(nn),Q0(nn),QQ(nd)dimension a(nn),xw(nn),dxa(nn),s(nn),alow0(n

2、n),alow(nn),h(nn)dimension u(nn),sx(nn),gb(nn),dy(nn),nday(ny)data Nday /365,731,1096,1461,1826,2192,2557,2922,3287,3653/open(10,file=地形.txt,status=old)open(12,file=深泓.txt,status=unknown)write(12,91) 91 format(3x,斷面號,6x距壩里程(千米),3x,深泓(米)open(13,file=水面線.txt,status=unknown)write(13,92) 92 format(3x,斷面

3、號,6x距壩里程(千米),3x,水面線高程(米)open(14,file=初始水位庫容.txt,status=unknown)write(14,93) 93 format(3x,年,6x水位(米),3x,庫容(億立方米)open(15,file=年淤積總量.txt,status=unknown)write(15,94) 94 format(3x年,3x淤積總量(萬立方米),3x年徑流量（萬方）,3x累計輸沙量（萬噸）,3x年均流量（萬方）)open(16,file=深泓（十年）.txt,status=unknown)write(16,95) 95 format(3x年,3x斷面號,3x深泓(米

4、)open(17,file=壩前斷面（十年）.txt,status=unknown)write(17,96) 96 format(3x年,3x起點距(米),3x高程(米)open(18,file=水位庫容（十年）.txt,status=unknown)write(18,97) 97 format(3x,年,6x水位(米),3x,庫容(億立方米)!=!=讀入地形數(shù)據(jù)=!=do i=1,npxtread(10,*)read(10,*)n1,npoint(i),dx(i),d1,dxa(i),rough(i)read(10,*)(x(j,1,i),j=1,npoint(i)read(10,*)(x(

5、j,2,i),j=1,npoint(i) end do!=!=讀入糙率=!=do i=1,npxtrough(i)=0.06enddo!=!=斷面間距計算=!=do i=1,30dx(i)=dxa(i)-dxa(i+1)end do!=!=判斷讀入數(shù)據(jù)的準確性=!=do i=1,npxtdo j=2,npoint(i)if(x(j,1,i)-x(j-1,1,i)0)thenwrite(*,*)mistake!endifenddoenddo!=!=計算初始深泓并輸出=!=do i=1,npxtalow0(i)=x(1,2,i)do j=2,npoint(i)if(x(j,2,i)-alow0(i

6、)0)thenalow0(i)=x(j,2,i)endifenddoenddodo i=1,npxtalow(i)=alow0(i)enddodo i=1,npxtwrite(12,*) i,dxa(i)/1000.0,alow(i)enddo!=!=計算初始時刻水位庫容并輸出=!=計算初始時刻水面線并輸出=!=計算初始時刻壩前斷面并輸出=!=do Z0=275,250,-1 zL=Z0*1.0 call FZV(NPOINT,npxt,X,ZL,ZVV,dx,ndisp,MM,NN) write(14,*) 0,ZL,ZVV enddodWQSt=0.0 dWQS=0.0 dWGb=0.0

7、dWQ=0.0wsout=0.0zvv=0wsin=0.0volu=0.0qq=0 zlevel(npxt)=275do i=1,npxtQ0(i)=377*0.5enddocall level(x,rough,npxt,zlevel,dx,Q0,npoint,b,a,xw,nn,mm,fai1,ndisp)do i=1,npxtwrite(13,*)0,i,dxa(i)/1000.0,zlevel(i) enddodo j=1,npoint(npxt)write(17,*)0,x(j,1,npxt),x(j,2,npxt)enddoQQ(1)=0!=!=讀入水沙數(shù)據(jù)=!=do k=1,5op

8、en(11,file=水沙.txt,status=old)do k1=1,ndread(11,*) Q(1),s(1) Q(1)=Q(1)*0.5do j=2,npxtQ(j)=Q(1)enddo!write(*,*)Q!pause!=!=利用張瑞瑾公式計算每個斷面的挾沙力Sx=!=利用mayer-peter公式計算推移質(zhì)輸沙率Gb=!=do i=1,npxth(i)=a(i)/b(i)u(i)=Q(i)/a(i)AA1=u(i)*3.0AA2=9.8*h(i)*wAA3=(AA1/AA2)*1.05Sx(i)=AA3*0.124enddogb(1)=s(1)*Q(1)*0.05do i=2,

9、npxtGb(i)=fGb(rough(i),r,h(i),Q(i),A(i),rs,d,rou,rous,ndisp)*b(i)!FUNCTION fGb(rough,r,h,Q,A,rs,d,rou,rous,ndisp)!write(*,*)gb(i)!pauseenddo!=!=計算各斷面含沙量和沖淤厚度=!=dy(1)=0.25*w*(s(1)-sx(1)*dt/rou0alow(1)=alow(1)+dy(1)if(alow(1)alow0(1)thenalow(1)=alow(1)-dy(1)dy(1)=0endif!write(*,*)dy(1)!pausedo i=2,npx

10、txishu=0.25*w*dx(i-1)*b(i)/q(1)s(i)=sx(i)+(s(i-1)-sx(i-1)*exp(-xishu)+(sx(i-1)-sx(i)/xishu*(1-exp(-xishu)if(s(i)0) thens(i)=0endifdGb=(Gb(i)-Gb(i-1)/dx(i-1)dQS=(Q(i)*S(i)-Q(i-1)*S(i-1)/dx(i-1)dy1=-(dGb+dQS)*dt/rou0dy2=-b(i-1)*dy(i-1)*(1-fai)dy3=(b(i)*fai)dy(i)=(dy1+dy2)/dy3!=!=判斷庫尾沖刷問題=!=alow(i)=alo

11、w(i)+dy(i)if(alow(i)alow0(i)thenalow(i)=alow(i)-dy(i)dy(i)=0gb(i)=gb(i-1)temp=rou0*b(i-1)*dy(i-1)*(fai-1)/dt*dx(i-1)s(i)=(temp-gb(i)+gb(i-1)/Q(i-1)+s(i-1)endif!if(k1.eq.28 .and.i.eq.9)then!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1),i!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1)!write(*,*)Q(i

12、-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)tt1,tt2,tt1,tt2!pause! endifif(s(i)0.05*abs(tt1)then!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1),i!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1)!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)Q(i-

13、1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)tt1,tt2,tt1,tt2!pause!endif!endif end dovolu=volu+dvvtotal1=dVV*rou0 !檢驗哪一個時段質(zhì)量不守恒 total2=(Q(1)*(s(1)-s(npxt)+(Gb(1)-Gb(npxt)*dt!write(*,*)dy,total1,total2!pauseif(abs(total1-total2)0.05*abs(total1)thenwrite(*,*)mass not equal,i,k1,total1,total2write(*,*)s(1

14、),s(npxt),Gb(1),Gb(npxt)pauseelseend if!=!=計算懸移質(zhì)推移質(zhì)及年均流量=!=WQS1=Q(1)*S(1)WQSi=Q(npxt)*S(npxt)dWQSt=dWQST+WQS1*dt/10000000 ! 為累計懸移質(zhì)輸沙量,單位為萬噸dWQS=dWQS+(WQS1-WQSi)*dt/(1.325*1000)/1E04! 為累計懸移質(zhì)沖淤量,單位為萬立方米dWGb=dwGb+(Gb(1)-Gb(npxt)*dt/(1.325*1000)/1E04! 為累計推移質(zhì)淤積總量,單位為萬立方米AVQ=AVQ+Q(1)*dt/10000 ! 為累計年徑流量 單位

15、萬方dWQ=dWQS+dWGbwsin=wsin+Gb(1)*dt/10000000!wsout=wsout+WQSi*dt/10000000!=!=水下斷面修改斷面點高程=!=do i=1,npxtdo j=1,Npoint(i)if(x(j,2,i)=zlevel(i) thenx(j,2,i)=zlevel(i)endifendifenddoenddoQQ(k1)=QQ(k1)+Q(1)QQ(0)=0!=!=輸出年淤積總量=!=do lk=1,nyif(k1=Nday(lk) thenif(lk1)thenQa=(QQ(k1)-QQ(k1-Nday(lk-1)/(Nday(lk)-Nda

16、y(lk-1)else Qa=QQ(k1)/Nday(lk)endifwrite(15,*) lk+10*(k-1),dWQ,avq,dWQSt+wsin,Qawrite(*,*) lk+10*(k-1),QaAVQ=0.0endifenddo!=!=十年一輸出的水面線變化=!=十年一輸出的深泓變化=!=十年一輸出壩前斷面變化=!=十年一輸出水位庫容關系=!=call level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,fai1,ndisp)if(k1=nday(10)thendo i=1,npxtwrite(13,*) 10*k,i,dxa(

17、i)/1000.0,zlevel(i)write(16,*)10*k,i,alow(i)enddodo j=1,npoint(npxt)write(17,*)10*k,x(j,1,npxt),x(j,2,npxt)enddodo Z0=275,250,-1zL=Z0*1.0call FZV(NPOINT,npxt,X,zl,ZVV,dx,ndisp,MM,NN) write(18,*)10*k,zl,zvvenddoendifenddoclose(11)enddoend!=!=計算斷面面積程序=!= SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)

18、! THIS SUBROUTINE IS USED TO CACULATE AREA FOR IRREGULAR CROSS-SECTION DIMENSION X(NPOINT),Z(NPOINT) IF(NDISP.LE.-1)WRITE(*,*)INTO AREA IF(NPOINT.LE.2)THEN WRITE(*,*)IN AREA THE INPUT DATA ARE FALSE N=,Npoint STOP ENDIF IF(Z(1).LT.ZLEVEL.OR.Z(Npoint).LT.ZLEVEL)THEN WRITE(*,*)IN AREA THE WATER LEVEL

19、IS TOO HIGH OR THE HEIGHT& OF THE SECTION AT EDGE IS TOO LOW STOP ENDIF A=0. B=0. XW=0. DO 10 I=1,NPOINT-1 ZMIN=AMIN1(Z(I),Z(I+1) IF(ZMIN.LT.ZLEVEL)THEN ZMAX=AMAX1(Z(I),Z(I+1) IF(ZMAX.LT.ZLEVEL)THEN DB=X(I+1)-X(I) DH=ZLEVEL-0.5*(Z(I)+Z(I+1) DX=(Z(I)-Z(I+1)*2.+DB*DB)*0.5 ELSE DB=(ZLEVEL-ZMIN)/(ZMAX-Z

20、MIN)*(X(I+1)-X(I) DH=0.5*(ZLEVEL-ZMIN) DX=(2*DH)*2.0+DB*DB)*0.5 ENDIF IF(DB.LT.0.0)THEN WRITE(*,*)THE DISTANCE FOR NODE ,I,AND,I+1,ARE FALSE WRITE(*,*)IN AREA,I,I+1,X(I),X(I+1) STOP ENDIF DA=DB*DH B=B+DB A=A+DA XW=XW+DX ENDIF10 CONTINUE return END!=!=恒定非均勻緩流水面線計算程序=!= subroutine level(x,rough,npxt,z

21、level,dx,q,npoint,b,a,xw,nn,mm,fai,NDISP) dimension x(mm,2,nn),rough(npxt),dx(npxt),zlevel(npxt),& q(npxt),npoint(npxt),b(npxt),a(npxt),xw(npxt) IF(NDISP.LE.-1)WRITE(*,*)INTO LEVEL nc=100 dz=0.1 dz1=0.03 dz2=0.1 zlevel(npxt)=267.0!write(*,*)npoint(npxt),x(1,1,npxt),x(1,2,npxt),b(npxt),&!a(npxt),xw(n

22、pxt),zlevel(npxt),ndisp !pause call area(npoint(npxt),x(1,1,npxt),x(1,2,npxt),b(npxt),& a(npxt),xw(npxt),zlevel(npxt),ndisp)do 10 ip=npxt-1,1,-1 NR=0 ! write(*,*)ip=,ip zmin=zlevel(ip+1)+dz ! WRITE(*,*)ZMIN1,ZMIN30 call area(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),& a(ip),xw(npxt),zmin,ndisp) if(a(ip).

23、LE.0)then zmin=zmin+0.173 goto 30 endif fmin=flevel(zlevel(ip+1),zmin,dx(ip),q(ip+1),q(ip),&rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai1,ndisp) !if(ip.eq.3.and.kl.eq.89)then!write(*,*)npoint(4),b(4),a(4),xw(4),zlevel(4),ndisp!write(*,*)(x(ikj,2,3),ikj=1,npoint(3)!write(*,*)fmin,zlevel(ip+1),zmin,dx(ip

24、),q(ip+1),q(ip)!write(*,*)rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai,ndisp!pause!endif if(fmin.gt.0)then fr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip) if(fr.lt.1)then zmax=zmin+dz20 call area(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),& a(ip),xw(ip),zmax,ndisp) fmax=flevel(zlevel(ip+1),zmax,dx(ip),q(ip+1),q(ip),& r

25、ough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai1,ndisp) if(fmax.gt.0.0)then zmax=zmax+dz goto 20 endif else zmin=zmin+dz2 NR=NR+1 IF(NR.GT.NC)THEN zmin=zmin-dz1 write(*,*)the loop is death,pause WRITE(*,*),nr,ip,zmin,fmin,nr,ip,zmin,fmin read(*,*) endif goto 30 endif else zmin=zmin-dz1 goto 30 endif call b

26、isec(zmin,zmax,fmin,zlevel(ip),npoint(ip),X(1,1,ip),& x(1,2,ip),b(ip),a(ip),xw(ip),q(ip),ndisp,zlevel(ip+1),b(ip+1),& a(ip+1),xw(ip+1),q(ip+1),dx(ip),rough(ip),fai) 10 continue return end !=二分法求根程序= subroutine bisec(zmin,zmax,fmin,zlevel,npoint,x,z,b,a,xw,q,& ndisp,zlelo,blower,alower,xwlower,qlower

27、,dx,rough,fai) dimension x(npoint),z(npoint) IF(NDISP.EQ.-1)WRITE(*,*)INTO BISEC ERR=0.00110 ddz=zmax-zmin zlevel=0.5*(zmin+zmax) call area(npoint,x,z,b,a,xw,zlevel,ndisp) f=flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,fai,ndisp) if(ddz.gt.err)then ff=f*fmin if(ff.gt.0)then zmin=zlevel el

28、se zmax=zlevel endif ddz=zmax-zmin goto 10 endif return end !=水面線函數(shù)計算程序= FUNCTION flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,fai,ndisp)!IF(NDISP.LE.-1) !WRITE(*,*)INTO FLEVEL HLOWER=ALOWER/BLOWER H=A/B AA=ZLELO-ZLEVEL B1=FAI*(QLOWER/BLOWER)*2.0/HLOWER*(10.0/3.0) B2=(1-FAI)*(Q/B)*2.0/H*(10