河流模擬課程設(shè)計(jì)報(bào)告書(shū)-武漢大學(xué)水利水電學(xué)院港航

1、河流模擬課程設(shè)計(jì)報(bào)告水庫(kù)一維泥沙淤積計(jì)算班級(jí)：2010級(jí)港航一班 姓名： 學(xué)號(hào)： 2010301580 指導(dǎo)老師： 吳衛(wèi)民 二一四年四月目錄一、設(shè)計(jì)任務(wù)及目的1二、計(jì)算條件與資料介紹12.1參數(shù)選取12.2地形資料22.3水沙資料22.4特征值統(tǒng)計(jì)2三、基本原理33.1基本方程及特征值、參數(shù)選取33.2方程的離散4四、計(jì)算步驟6五、關(guān)于庫(kù)尾沖刷問(wèn)題7六、計(jì)算成果96.1淤積量隨時(shí)間變化96.2縱剖面形態(tài)分析106.3壩前斷面形態(tài)圖116.4泥沙質(zhì)量守恒性分析126.5庫(kù)容損失合理性分析13附：151 程序內(nèi)變量說(shuō)明152 Fortran計(jì)算源程序16一、設(shè)計(jì)任務(wù)及目的該課程設(shè)計(jì)給定水沙資料和斷

2、面形態(tài)資料，解決水庫(kù)一維泥沙淤積計(jì)算問(wèn)題。以水流、泥沙方程為基礎(chǔ)，構(gòu)建恒定流條件下的河道一維水沙數(shù)學(xué)模型，并編制出完整的計(jì)算程序，并以某個(gè)水庫(kù)為實(shí)例，進(jìn)行水庫(kù)泥沙淤積計(jì)算。水流條件：恒定非均勻流。泥沙條件：包括懸移質(zhì)，推移質(zhì)的均勻沙模型，推移質(zhì)計(jì)算模式為飽和輸沙，懸移質(zhì)計(jì)算模式為不飽和輸沙，水流泥沙方程采用非耦合解。通過(guò)課程設(shè)計(jì)，初步掌握一維數(shù)學(xué)模型建立數(shù)學(xué)模型的的基本過(guò)程和計(jì)算方法，具備一定的解決實(shí)際問(wèn)題的能力。二、計(jì)算條件與資料介紹資料見(jiàn)磁盤(pán)文件，包括：1河道地形，糙率等(.txt文件)；2日平均流量，日平均含沙量(.xlsx文件)2.1參數(shù)選取其中，推移質(zhì)來(lái)沙量按懸移質(zhì)的5計(jì)算。均勻沙粒

3、徑為d=0.05mm,推移質(zhì)d=2mm。其余參數(shù)在基本原理中說(shuō)明。2.2地形資料地形資料包括水庫(kù)前31個(gè)斷面的斷面起點(diǎn)距和高程，最遠(yuǎn)距壩里程為94645米。2.3水沙資料計(jì)算系列采用1991年2000年10年水沙數(shù)據(jù)，重復(fù)5次，計(jì)算50年的沖淤變化。此外，將給出的流量資料均乘以0.5后作為日平均流量。2.4特征值統(tǒng)計(jì)表2. 1水沙特征值統(tǒng)計(jì)年平均流量(m3/s)年平均含沙量(kg/m3)年輸沙量(萬(wàn)t)年徑流量(億m3)1991年768 1.50 3624 242 1992年566 1.36 2437 179 1993年723 1.46 3339 228 1994年853 1.12 3017

4、269 1995年898 0.96 2724 283 1996年782 0.88 2186 247 1997年1071 1.21 4095 338 1998年707 0.69 1540 223 1999年804 0.68 1729 254 2000年846 0.54 1437 267 三、基本原理3.1基本方程及特征值、參數(shù)選取水流連續(xù)方程：水流運(yùn)動(dòng)方程 或 泥沙連續(xù)方程 河床變形方程恢復(fù)飽和系數(shù) 推移質(zhì)平衡輸沙方程 水流挾沙力公式(張瑞瑾公式)： ,推移質(zhì)輸沙率公式采用Mayer-Peter公式(適用于寬淺河槽)： n為河床平整情況下的沙粒曼寧糙率系數(shù)。 為粒配曲線中90%較之為小的粒徑MA

5、YER-PETER公式中的能坡J按均勻流曼寧公式近似計(jì)算（每個(gè)斷面不同）。曼寧公式： 寬淺河槽， 3.2方程的離散方程（3.1）在恒定流情況下有離散為:Q=const方程（3.3）變形為或上式離散為方程（3.4）去掉時(shí)間項(xiàng)得到該方程的解析解為：由方程（3.4、3.5）可得:對(duì)2 號(hào)斷面以下，上式可以離散為：移項(xiàng)變形為：對(duì)于進(jìn)口斷面，推移質(zhì)不考慮，懸移質(zhì)采用單點(diǎn)離散,方程（5）可離散為：四、計(jì)算步驟（1）輸入河床地形糙率等數(shù)據(jù)（2）讀入一個(gè)時(shí)段的水沙數(shù)據(jù)（特別注意，不要一次性將數(shù)據(jù)全部讀入）（3）計(jì)算水面線，同時(shí)得到各斷面的水力要素。計(jì)算前要注意在壩前輸入水位，各斷面均應(yīng)對(duì)流量賦值（4）計(jì)算水流

6、挾沙力（5）計(jì)算推移質(zhì)輸沙率（6）計(jì)算各斷面含沙量（7）計(jì)算各斷面沖淤厚度（8）修改水各斷面水下河床高程（9）重新進(jìn)入（2）進(jìn)行下一循環(huán)（10）計(jì)算50年河床變形，計(jì)算時(shí)段為一天，單位為秒（s）（11）每10年輸出一次計(jì)算結(jié)果(如淤積總量等)五、關(guān)于庫(kù)尾沖刷問(wèn)題按均勻沙模型計(jì)算水庫(kù)的淤積，在庫(kù)尾計(jì)算結(jié)果可能會(huì)發(fā)生沖刷現(xiàn)象。這與實(shí)際情況不符，這是因?yàn)榫鶆蛏衬Ｐ团c天然情況有較大差異所至。天然河流實(shí)際河床泥沙級(jí)配為非均勻沙，由于水庫(kù)一般建在山區(qū)河流，屬下切型河道，河床多為較大的漂石、礫石、甚至是基巖，如果是非均勻沙模型，來(lái)沙粒徑與河床粒徑可以分為多組，上游來(lái)沙粒徑與河床粒徑可以根據(jù)實(shí)際情況分別給出，

7、可以把河床粒徑定義的很大，使得河床泥沙不會(huì)發(fā)生沖刷，而來(lái)沙粒徑一般較細(xì)，在庫(kù)尾不會(huì)淤積，直接泄向下游，因?yàn)閹?kù)尾與天然情況較為接近，這種情況下庫(kù)尾一般是不沖不淤。在均勻沙模型中，無(wú)論河床還是水流中的泥沙只有一個(gè)粒徑，其代表性不能模擬實(shí)際情況，一般來(lái)說(shuō)，所選粒徑屬于懸移質(zhì)，這種粒徑的泥沙在天然河流中可能屬?zèng)_瀉質(zhì)，不可能在河床中存在，用這種粒徑代替庫(kù)尾的漂石、礫石等粒徑，而庫(kù)尾的流速很大，自然會(huì)發(fā)生沖刷。這種結(jié)果當(dāng)然是錯(cuò)誤的，如果發(fā)生沖刷，應(yīng)該在建庫(kù)前就發(fā)生了，而不是建庫(kù)后。解決的辦法：人為限制庫(kù)尾的沖刷，即允許河道沖刷掉已經(jīng)淤積的泥沙，而不允許沖刷原始河床，程序處理方法：1 用變量記憶每個(gè)斷面的初

8、始深泓點(diǎn)高程；2 每個(gè)時(shí)段計(jì)算每個(gè)斷面的沖淤面積時(shí)，如果沖淤面積為負(fù)（沖刷），此時(shí)判斷該斷面當(dāng)前深泓高程與初始河床的關(guān)系，如果小于等于原始高程，則不再允許本河段沖刷，即令沖淤面積為0，反之則可以沖刷（此時(shí)沖走的是前期淤積物）。3上面情況發(fā)生時(shí)要特別注意泥沙的連續(xù)性（質(zhì)量守恒），沖淤面積為0意味著下斷面通過(guò)的泥沙數(shù)量與上斷面相同，這種情況下，下斷面含沙量和推移質(zhì)均與上斷面相同。4要檢查泥沙連續(xù)性1）入庫(kù)沙量減去出庫(kù)沙量等于水庫(kù)中泥沙淤積量，按時(shí)段檢查并注意重量與體積的換算。2）如果不相等要逐一檢查每個(gè)區(qū)間的連續(xù)性，上斷面沙量減去下斷面沙量等于區(qū)間淤積量。3）為了進(jìn)行檢查要用變量記憶下每個(gè)斷面的累

9、計(jì)淤積面積。5 注意程序設(shè)計(jì)時(shí)除要單獨(dú)計(jì)算水面線要外，外面要進(jìn)行區(qū)間循環(huán)，在每個(gè)小區(qū)間內(nèi)完成含沙量、推移質(zhì)河床變形計(jì)算、修改河床高程等后再進(jìn)行下一個(gè)循環(huán)。如果先計(jì)算所有斷面含沙量后，再調(diào)整出口含沙量時(shí)，下一個(gè)區(qū)間進(jìn)口含沙量就不對(duì)了。六、計(jì)算成果6.1淤積量隨時(shí)間變化表6. 1水庫(kù)淤積量(體積法)隨時(shí)間變化統(tǒng)計(jì)淤積年數(shù)淤積量(億t)淤積量(億m3)102.37 1.79 204.46 3.36 306.19 4.67 407.43 5.61 508.18 6.18 圖6. 1水庫(kù)淤積量隨時(shí)間變化圖分析：1淤積量(質(zhì)量)量級(jí)為億t；2隨著淤積時(shí)間的增加，淤積量呈增加態(tài)勢(shì)；3淤積量曲線斜率逐漸減小，

10、表明淤積速度隨時(shí)間的增加逐漸減緩。326.2縱剖面形態(tài)分析圖6. 2 深泓線高程變化圖分析：1隨著時(shí)間的增加，深泓高程逐年抬升；2隨著距壩里程的減小，深泓高程增加值越大，表明淤積程度越往壩前越明顯。6.3壩前斷面形態(tài)圖圖6. 3 近壩斷面形態(tài)圖分析：1隨著時(shí)間的增加，河床斷面各點(diǎn)高程逐年抬升，表明經(jīng)歷了持續(xù)的淤積過(guò)程；2淤積大概在280m高程以下，280m以上部分沒(méi)有淤積。3 40年后，淤積變化不明顯，河床高程變化已經(jīng)很小。6.4泥沙質(zhì)量守恒性分析檢驗(yàn)方法：入庫(kù)泥沙出庫(kù)泥沙水庫(kù)內(nèi)泥沙淤積量(體積法)對(duì)于某一時(shí)段內(nèi)：入庫(kù)泥沙量出庫(kù)泥沙量(質(zhì)量)：水庫(kù)內(nèi)泥沙淤積量(質(zhì)量)：表6. 2泥沙質(zhì)量守恒性

11、分析淤積年數(shù)水庫(kù)內(nèi)泥沙淤積量(體積法)(億t)入庫(kù)減出庫(kù)泥沙量(億t)誤差()102.3656 2.3657 0.0697 204.4568 4.4573 0.1209 306.1868 6.1879 0.1710 407.4330 7.4346 0.2069 508.1826 8.1846 0.2407 分析：1利用體積法計(jì)算的水庫(kù)泥沙淤積量與利用入庫(kù)泥沙減出庫(kù)泥沙的方法計(jì)算的結(jié)果誤差在0.5以內(nèi)，基本符合泥沙守恒原則；2 隨著時(shí)間的增加，淤積量加大，兩種算法結(jié)果之間的誤差呈遞增態(tài)勢(shì)。6.5庫(kù)容損失合理性分析庫(kù)容變化曲線：圖6. 4 水位庫(kù)容關(guān)系曲線各年庫(kù)容變化：表6. 3水位庫(kù)容變化關(guān)系單

12、位：(萬(wàn)m3)水位原始庫(kù)容10年20年30年40年50年25116357 11154 5704 1379 126 48 25217241 11762 6108 1593 168 50 25318173 12394 6525 1827 225 54 25419147 13049 6955 2083 295 61 25520152 13732 7398 2349 376 77 25621185 14439 7856 2627 469 104 25722249 15168 8330 2917 576 140 25823354 15923 8821 3220 698 184 25924523 1670

13、6 9333 3538 838 240 26025773 17519 9867 3872 992 312 26127071 18364 10425 4219 1164 401 26228427 19237 11012 4580 1357 505 26329828 20139 11623 4961 1569 621 26431291 21078 12262 5367 1803 751 26532838 22061 12934 5799 2067 902 26634452 23093 13640 6262 2356 1073 26736162 24165 14382 6759 2670 1266

14、26837948 25301 15160 7295 3009 1483 26939819 26517 15982 7866 3370 1733 27041778 27813 16858 8475 3766 2017 27143807 29184 17808 9128 4202 2324 27245882 30674 18840 9834 4681 2659 27348001 32299 19974 10610 5203 3039 27450165 34062 21198 11484 5785 3473 27552374 35934 22543 12474 6443 3959 分析：1 隨著時(shí)間

15、的增加，同一水位下庫(kù)容呈逐漸減小態(tài)勢(shì)，表明水庫(kù)在逐年淤積； 2 時(shí)間越久，庫(kù)容減小的程度越小，水庫(kù)的淤積速率在逐漸下降。3 50年后，水庫(kù)庫(kù)容減少量巨大，應(yīng)引起充分重視。附：1 程序內(nèi)變量說(shuō)明2 Fortran計(jì)算源程序program mainparameter(mm=80,nn=31,fai=0.5,fai2=0.8)real zlevel,Q,s,k,m,g,dxuan,w,JJ,h,dtui,rough2,p,ps,r,rs,afa,p0,dt,fai2,mianji,rujianchu,panduan,dy,jiancha,x,zuidi,zuidi2,kurong,jingliuli

16、anginteger year,month,day,number,nQSdimension x(mm,2,nn),number(nn),npoint(nn),dx(nn),disu(nn),disd(nn),rough(nn),Q(nn),s(nn),zlevel(nn),B(nn),A(nn),v(nn),ss(nn),h(nn),JJ(nn),gb(nn),dy(nn),zuidi(nn),zuidi2(nn),kurong(26)!讀入斷面數(shù)據(jù)和壩前水位,數(shù)據(jù)從第十行開(kāi)始o(jì)pen(10,file=section.txt,status=old)read(10,*)read(10,*)rea

17、d(10,*)read(10,*)read(10,*)read(10,*)read(10,*)read(10,*)read(10,*)do i=1,nnread(10,*)read(10,*)number(i),npoint(i),dx(i),disu(i),disd(i),rough(i)read(10,*)(x(j,1,i),j=1,npoint(i)read(10,*)(x(j,2,i),j=1,npoint(i)enddoread(10,*)zlevel(nn)close(10)!原始地形深泓點(diǎn)輸入do i=1,nnzuidi(i)=minval(x(1:npoint(i),2,i)e

18、nddo!原始數(shù)據(jù)輸出和統(tǒng)計(jì)部分:深泓和27斷面形態(tài)!1open(70,file=0ShenHong.txt,status=new)write(70,*)原始地形深泓do i=1,nnwrite(70,*)disd(i),zuidi(i)enddoclose(70)!2open(60,file=0mouduanmian.txt,status=new)write(60,*)原始?jí)吻皵嗝嫫瘘c(diǎn)距高程do j=1,npoint(27)write(60,*)x(j,1,27),x(j,2,27)enddoclose(60)!3計(jì)算初始靜庫(kù)容,j的1-26對(duì)應(yīng)水位250-275open(90,file=0

19、kurong.txt,status=new)write(90,*)水位庫(kù)容do j=1,26shuiwei=j+249.0do i=1,nn-1call area(npoint(i),x(1,1,i),x(1,2,i),shuiwei,B(i),A(i)kurong(j)=kurong(j)+(A(i)*(1-fai2)+A(i+1)*fai2)*dx(i+1)enddo write(90,*)shuiwei,kurong(j)enddoclose(90)!提示：時(shí)間節(jié)點(diǎn)對(duì)應(yīng)數(shù)據(jù)組數(shù)10年3653,20年7306,30年10959,40年14612，50年18265nQS=18265rukus

20、haliang=0.0chukushaliang=0.0jingliuliang=0.0zongliuliang=0.0zonghansha=0.0rujianchu=0.0yujiliang=0.0!為了看每年淤積量與來(lái)沙量之比設(shè)置yujiliang2(rujianchu與yujiliang不做每年清零處理)yujiliang2=0.0kurong(1:26)=0.0dt=24*60.0*60.0!讀入流量與水沙資料open(20,file=rewriteQS.txt,status=old)read(20,*)do L=1,nQSread(20,*)year,month,day,Q(1),s

21、(1)Q(1)=Q(1)*0.5 do i=2,nn Q(i)=Q(1)enddowrite(*,*)該組流量和來(lái)沙量為:,Q(1),s(1)!推求水面線call level(x,rough,nn,zlevel,dx,Q,npoint,B,A,mm,fai)WRITE(*,*)累積里程disd（自壩址斷面起算） 水位ZLEVEL do i=1,nn WRITE(*,*)disd(i),ZLEVEL(i)enddodo i=1,nnv(i)=Q(i)/A(i)h(i)=A(i)/B(i)enddo!計(jì)算泥沙沉速和水流挾沙力dxuan=0.05*10*(-3.0)call FALLV(dxuan,

22、w)write(*,*) w(懸)=,wk=0.124m=1.05g=9.81do i=1,nnss(i)=k*(v(i)*3.0/(g*h(i)*w)*menddo!計(jì)算能坡do i=1,nnJJ(i)=v(i)*2.0*rough(i)*2.0/h(i)*(4.0/3.0)enddo!計(jì)算推移質(zhì)輸沙率dtui=2*10*(-3.0)rough2=dtui*(1.0/6.0)/26.0p=1ps=2.65r=p*grs=ps*g!單位說(shuō)明：公式中r單位為kN/m3,p單位未為t/m3gb(1)=Q(1)*s(1)*0.05/B(1)do i=2,nnpanduan=(rough2/rough

23、(i)*(3.0/2.0)*r*h(i)*JJ(i)-0.047*(rs-r)*dtuiif(panduan.gt.0)then gb(i)=panduan*(3.0/2.0)/(0.125*p*0.5*(ps-p)/ps*g)*10*3.0elsegb(i)=0endifenddowrite(*,*)斷面推移質(zhì)輸沙率gbdo i=1,nnwrite(*,*)i,gb(i)enddo!計(jì)算各斷面含沙量、沖淤厚度，并修改斷面高程，1斷面單獨(dú)進(jìn)行!每個(gè)斷面要進(jìn)行如下判斷：1加上dy是否沖到原始河床以下。2斷面高程修改部分應(yīng)在水位以下。3若不符合條件1，則令dy=0且調(diào)整含沙量S(除進(jìn)口斷面外)，若

24、調(diào)整后S小于零則S=0afa=0.25p0=1.325*1000dt=86400.0!1斷面write(*,*)斷面沖淤變化量dydy(1)=afa*w*(s(1)-ss(1)*dt/p0jiancha=minval(x(1:npoint(1),2,1)panduan=jiancha+dy(1)if(panduan.ge.zuidi(1)thenwrite(*,*)1,dy(1)do j=1,npoint(1)if(x(j,2,1).le.ZLEVEL(1)then x(j,2,1)=x(j,2,1)+dy(1)endifenddoelsedy(1)=0write(*,*)1,dy(1),斷面

25、沖至原始河床以下,已進(jìn)行修改dy=0endif!2斷面到31斷面do i=1,nn-1qaver=Q(i)*2.0/(B(i)+B(i+1)s(i+1)=ss(i+1)+(s(i)-ss(i)*exp(-afa*w*dx(i+1)/qaver)+(ss(i)-ss(i+1)*qaver/afa/w/dx(i+1)*(1-exp(-afa*w*dx(i+1)/qaver)ds1=(B(i+1)*gb(i+1)-B(i)*gb(i)/dx(i+1)ds2=(Q(i+1)*s(i+1)-Q(i)*s(i)/dx(i+1)dv1=B(i)*(1-fai2)*dy(i)dy(i+1)=(ds1+ds2)

26、/p0*dt+dv1)/(-B(i+1)*fai2)jiancha=minval(x(1:npoint(i+1),2,i+1)panduan=jiancha+dy(i+1)if(panduan.ge.zuidi(i+1)thenwrite(*,*)i+1,dy(i+1)do j=1,npoint(i+1)if(x(j,2,i+1).le.ZLEVEL(i+1)then x(j,2,i+1)=x(j,2,i+1)+dy(i+1)endifenddoelsedy(i+1)=0write(*,*)i+1,dy(i+1),斷面沖至原始河床以下,已進(jìn)行修改dy=0S(i+1)=(p0*dx(i+1)*B

27、(i)*dy(i)*(1-fai2)/dt+B(i+1)*gb(i+1)-B(i)*gb(i)-Q(i)*s(i)/(-Q(i+1)if(S(i+1).lt.0)S(i+1)=0endifenddowrite(*,*)斷面含沙量s挾沙力ssdo i=1,nnwrite(*,*)i,s(i),ss(1)enddo!驗(yàn)算淤積守恒性及部分特征值統(tǒng)計(jì)結(jié)果write(*,*)該來(lái)水來(lái)沙條件下河床地形已調(diào)整，以下驗(yàn)證淤積守恒rukushaliang=rukushaliang+Q(1)*s(1)*dtchukushaliang=chukushaliang+Q(31)*S(31)*dtjingliuliang

28、=jingliuliang+Q(1)*dtzongliuliang=zongliuliang+Q(1)zonghansha=zonghansha+Q(1)*S(1)rujianchu=rujianchu+(B(1)*gb(1)-B(31)*gb(31)+Q(1)*S(1)-Q(31)*S(31)*dtdo i=1,nn-1mianji=B(i)*dy(i)*(1-fai2)+B(i+1)*dy(i+1)*fai2yujiliang=yujiliang+p0*dx(i+1)*mianjiyujiliang2=yujiliang2+p0*dx(i+1)*mianjienddowrite(*,*)入

29、庫(kù)-出庫(kù)泥沙淤積量write(*,*)rujianchu,yujiliangwrite(*,*)以上數(shù)據(jù)L=,L!進(jìn)行年平均流量、含沙量、輸沙量、徑流量統(tǒng)計(jì)!一年一次,1991if(L.eq.365)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file=1991#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujili

30、ang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1992-366if(L.eq.731)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliuliangopen(100,file=1992#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukush

31、aliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1993if(L.eq.1096)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file=1993#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjl

32、iuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1994if(L.eq.1461)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file=1994#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量

33、 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1995if(L.eq.1826)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file=1995#txt,status=new)wr

34、ite(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1996-366if(L.eq.2192)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliuliangopen(100,file=1

35、996#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1997if(L.eq.2557)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zong

36、liuliangopen(100,file=1997#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1998if(L.eq.2922)thenpjliuliang=zongliuliang/365p

37、jhansha=zonghansha/zongliuliangopen(100,file=1998#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1999if(L.eq.3287)thenpjliu

38、liang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file=1999#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,200

39、0-366if(L.eq.3653)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliuliangopen(100,file=2000#txt,status=new)write(100,*)年平均流量 含沙量 輸沙量 徑流量 淤積量write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=

40、0yujiliang2=0write(*,*)年統(tǒng)計(jì)完畢，進(jìn)行第11年的計(jì)算endif!進(jìn)行沖淤后數(shù)據(jù)輸出和統(tǒng)計(jì)部分：深泓和27斷面形態(tài)，靜庫(kù)容計(jì)算!10年3653if(L.eq.3653)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file=10ShenHong.txt,status=new)write(70,*)沖淤后深泓do i=1,nnwrite(70,*)zuidi2(i)enddoclose(70)!2open(60,file=10mouduanmian.txt,status=new)write(60

41、,*)27斷面高程do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file=10yujiliang.txt,status=new)write(80,*)淤積量入庫(kù)-出庫(kù)write(80,*)yujiliang,rujianchuclose(80)!4open(90,file=10kurong.txt,status=new)write(90,*)庫(kù)容do j=1,26shuiwei=j+249kurong(1:26)=0do i=1,nn-1call area(npoint(i),x(1,1,i),x(1,2,i),shu

42、iwei,B(i),A(i)kurong(j)=kurong(j)+(A(i)*(1-fai2)+A(i+1)*fai2)*dx(i+1)enddo write(90,*)kurong(j)enddoclose(90)endif!20年7306if(L.eq.7306)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file=20ShenHong.txt,status=new)write(70,*)沖淤后深泓do i=1,nnwrite(70,*)zuidi2(i)enddoclose(70)!2open(60,fi

43、le=20mouduanmian.txt,status=new)do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file=20yujiliang.txt,status=new)write(80,*)淤積量入庫(kù)-出庫(kù)write(80,*)yujiliang,rujianchuclose(80)!4open(90,file=20kurong.txt,status=new)write(90,*)庫(kù)容do j=1,26shuiwei=j+249kurong(1:26)=0do i=1,nn-1call area(npoint(i

44、),x(1,1,i),x(1,2,i),shuiwei,B(i),A(i)kurong(j)=kurong(j)+(A(i)*(1-fai2)+A(i+1)*fai2)*dx(i+1)enddo write(90,*)kurong(j)enddoclose(90)endif!30年10959if(L.eq.10959)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file=30ShenHong.txt,status=new)write(70,*)沖淤后深泓do i=1,nnwrite(70,*)zuidi2(i)e

45、nddoclose(70)!2open(60,file=30mouduanmian.txt,status=new)do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file=30yujiliang.txt,status=new)write(80,*)淤積量入庫(kù)-出庫(kù)write(80,*)yujiliang,rujianchuclose(80)!4open(90,file=30kurong.txt,status=new)write(90,*)庫(kù)容do j=1,26shuiwei=j+249kurong(1:26)=0do i

46、=1,nn-1call area(npoint(i),x(1,1,i),x(1,2,i),shuiwei,B(i),A(i)kurong(j)=kurong(j)+(A(i)*(1-fai2)+A(i+1)*fai2)*dx(i+1)enddo write(90,*)kurong(j)enddoclose(90)endif!40年14612if(L.eq.14612)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file=40ShenHong.txt,status=new)write(70,*)沖淤后深泓do i=

47、1,nnwrite(70,*)zuidi2(i)enddoclose(70)!2open(60,file=40mouduanmian.txt,status=new)do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file=40yujiliang.txt,status=new)write(80,*)淤積量入庫(kù)-出庫(kù)write(80,*)yujiliang,rujianchuclose(80)!4open(90,file=40kurong.txt,status=new)write(90,*)庫(kù)容do j=1,26shuiwe

48、i=j+249kurong(1:26)=0do i=1,nn-1call area(npoint(i),x(1,1,i),x(1,2,i),shuiwei,B(i),A(i)kurong(j)=kurong(j)+(A(i)*(1-fai2)+A(i+1)*fai2)*dx(i+1)enddo write(90,*)kurong(j)enddoclose(90)endif!50年18265if(L.eq.18265)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file=50ShenHong.txt,status=new)write(70,*)沖淤后深泓do i=1,nnwrite(70,*)zuidi2(i)enddoclose(70)!2open(60,file=50mouduanmian.txt,status=new)do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file=50yujiliang.t