河流模拟课程设计报告书-武汉大学水利水电学院港航

1、河流模拟课程设计报告水库一维泥沙淤积计算班级：2010级港航一班 姓名： 学号： 2010301580 指导老师： 吴卫民 二一四年四月目录一、设计任务及目的1二、计算条件与资料介绍12.1参数选取12.2地形资料22.3水沙资料22.4特征值统计2三、基本原理33.1基本方程及特征值、参数选取33.2方程的离散4四、计算步骤6五、关于库尾冲刷问题7六、计算成果96.1淤积量随时间变化96.2纵剖面形态分析106.3坝前断面形态图116.4泥沙质量守恒性分析126.5库容损失合理性分析13附：151 程序内变量说明152 Fortran计算源程序16一、设计任务及目的该课程设计给定水沙资料和断

2、面形态资料，解决水库一维泥沙淤积计算问题。以水流、泥沙方程为基础，构建恒定流条件下的河道一维水沙数学模型，并编制出完整的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计算模式为不饱和输沙，水流泥沙方程采用非耦合解。通过课程设计，初步掌握一维数学模型建立数学模型的的基本过程和计算方法，具备一定的解决实际问题的能力。二、计算条件与资料介绍资料见磁盘文件，包括：1河道地形，糙率等(.txt文件)；2日平均流量，日平均含沙量(.xlsx文件)2.1参数选取其中，推移质来沙量按悬移质的5计算。均匀沙粒

3、径为d=0.05mm,推移质d=2mm。其余参数在基本原理中说明。2.2地形资料地形资料包括水库前31个断面的断面起点距和高程，最远距坝里程为94645米。2.3水沙资料计算系列采用1991年2000年10年水沙数据，重复5次，计算50年的冲淤变化。此外，将给出的流量资料均乘以0.5后作为日平均流量。2.4特征值统计表2. 1水沙特征值统计年平均流量(m3/s)年平均含沙量(kg/m3)年输沙量(万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基本方程及特征值、参数选取水流连续方程：水流运动方程 或 泥沙连续方程 河床变形方程恢复饱和系数 推移质平衡输沙方程 水流挟沙力公式(张瑞瑾公式)： ,推移质输沙率公式采用Mayer-Peter公式(适用于宽浅河槽)： n为河床平整情况下的沙粒曼宁糙率系数。 为粒配曲线中90%较之为小的粒径MA

5、YER-PETER公式中的能坡J按均匀流曼宁公式近似计算（每个断面不同）。曼宁公式： 宽浅河槽， 3.2方程的离散方程（3.1）在恒定流情况下有离散为:Q=const方程（3.3）变形为或上式离散为方程（3.4）去掉时间项得到该方程的解析解为：由方程（3.4、3.5）可得:对2 号断面以下，上式可以离散为：移项变形为：对于进口断面，推移质不考虑，悬移质采用单点离散,方程（5）可离散为：四、计算步骤（1）输入河床地形糙率等数据（2）读入一个时段的水沙数据（特别注意，不要一次性将数据全部读入）（3）计算水面线，同时得到各断面的水力要素。计算前要注意在坝前输入水位，各断面均应对流量赋值（4）计算水流

6、挟沙力（5）计算推移质输沙率（6）计算各断面含沙量（7）计算各断面冲淤厚度（8）修改水各断面水下河床高程（9）重新进入（2）进行下一循环（10）计算50年河床变形，计算时段为一天，单位为秒（s）（11）每10年输出一次计算结果(如淤积总量等)五、关于库尾冲刷问题按均匀沙模型计算水库的淤积，在库尾计算结果可能会发生冲刷现象。这与实际情况不符，这是因为均匀沙模型与天然情况有较大差异所至。天然河流实际河床泥沙级配为非均匀沙，由于水库一般建在山区河流，属下切型河道，河床多为较大的漂石、砾石、甚至是基岩，如果是非均匀沙模型，来沙粒径与河床粒径可以分为多组，上游来沙粒径与河床粒径可以根据实际情况分别给出，

7、可以把河床粒径定义的很大，使得河床泥沙不会发生冲刷，而来沙粒径一般较细，在库尾不会淤积，直接泄向下游，因为库尾与天然情况较为接近，这种情况下库尾一般是不冲不淤。在均匀沙模型中，无论河床还是水流中的泥沙只有一个粒径，其代表性不能模拟实际情况，一般来说，所选粒径属于悬移质，这种粒径的泥沙在天然河流中可能属冲泻质，不可能在河床中存在，用这种粒径代替库尾的漂石、砾石等粒径，而库尾的流速很大，自然会发生冲刷。这种结果当然是错误的，如果发生冲刷，应该在建库前就发生了，而不是建库后。解决的办法：人为限制库尾的冲刷，即允许河道冲刷掉已经淤积的泥沙，而不允许冲刷原始河床，程序处理方法：1 用变量记忆每个断面的初

8、始深泓点高程；2 每个时段计算每个断面的冲淤面积时，如果冲淤面积为负（冲刷），此时判断该断面当前深泓高程与初始河床的关系，如果小于等于原始高程，则不再允许本河段冲刷，即令冲淤面积为0，反之则可以冲刷（此时冲走的是前期淤积物）。3上面情况发生时要特别注意泥沙的连续性（质量守恒），冲淤面积为0意味着下断面通过的泥沙数量与上断面相同，这种情况下，下断面含沙量和推移质均与上断面相同。4要检查泥沙连续性1）入库沙量减去出库沙量等于水库中泥沙淤积量，按时段检查并注意重量与体积的换算。2）如果不相等要逐一检查每个区间的连续性，上断面沙量减去下断面沙量等于区间淤积量。3）为了进行检查要用变量记忆下每个断面的累

9、计淤积面积。5 注意程序设计时除要单独计算水面线要外，外面要进行区间循环，在每个小区间内完成含沙量、推移质河床变形计算、修改河床高程等后再进行下一个循环。如果先计算所有断面含沙量后，再调整出口含沙量时，下一个区间进口含沙量就不对了。六、计算成果6.1淤积量随时间变化表6. 1水库淤积量(体积法)随时间变化统计淤积年数淤积量(亿t)淤积量(亿m3)102.37 1.79 204.46 3.36 306.19 4.67 407.43 5.61 508.18 6.18 图6. 1水库淤积量随时间变化图分析：1淤积量(质量)量级为亿t；2随着淤积时间的增加，淤积量呈增加态势；3淤积量曲线斜率逐渐减小，

10、表明淤积速度随时间的增加逐渐减缓。326.2纵剖面形态分析图6. 2 深泓线高程变化图分析：1随着时间的增加，深泓高程逐年抬升；2随着距坝里程的减小，深泓高程增加值越大，表明淤积程度越往坝前越明显。6.3坝前断面形态图图6. 3 近坝断面形态图分析：1随着时间的增加，河床断面各点高程逐年抬升，表明经历了持续的淤积过程；2淤积大概在280m高程以下，280m以上部分没有淤积。3 40年后，淤积变化不明显，河床高程变化已经很小。6.4泥沙质量守恒性分析检验方法：入库泥沙出库泥沙水库内泥沙淤积量(体积法)对于某一时段内：入库泥沙量出库泥沙量(质量)：水库内泥沙淤积量(质量)：表6. 2泥沙质量守恒性

11、分析淤积年数水库内泥沙淤积量(体积法)(亿t)入库减出库泥沙量(亿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利用体积法计算的水库泥沙淤积量与利用入库泥沙减出库泥沙的方法计算的结果误差在0.5以内，基本符合泥沙守恒原则；2 随着时间的增加，淤积量加大，两种算法结果之间的误差呈递增态势。6.5库容损失合理性分析库容变化曲线：图6. 4 水位库容关系曲线各年库容变化：表6. 3水位库容变化关系单

12、位：(万m3)水位原始库容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 随着时间

15、的增加，同一水位下库容呈逐渐减小态势，表明水库在逐年淤积； 2 时间越久，库容减小的程度越小，水库的淤积速率在逐渐下降。3 50年后，水库库容减少量巨大，应引起充分重视。附：1 程序内变量说明2 Fortran计算源程序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)!读入断面数据和坝前水位,数据从第十行开始open(10,file='section.txt',status='old')rea

17、d(10,*)read(10,*)read(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)!原始地形深泓点输入do i=1,nnzuidi(i)=minval

18、(x(1:npoint(i),2,i)enddo!原始数据输出和统计部分:深泓和27断面形态!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,*)'原始坝前断面起点距高程'do j=1,npoint(27)write

19、(60,*)x(j,1,27),x(j,2,27)enddoclose(60)!3计算初始静库容,j的1-26对应水位250-275open(90,file='0kurong.txt',status='new')write(90,*)'水位库容'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 wr

20、ite(90,*)shuiwei,kurong(j)enddoclose(90)!提示：时间节点对应数据组数10年3653,20年7306,30年10959,40年14612，50年18265nQS=18265rukushaliang=0.0chukushaliang=0.0jingliuliang=0.0zongliuliang=0.0zonghansha=0.0rujianchu=0.0yujiliang=0.0!为了看每年淤积量与来沙量之比设置yujiliang2(rujianchu与yujiliang不做每年清零处理)yujiliang2=0.0kurong(1:26)=0.0dt=2

21、4*60.0*60.0!读入流量与水沙资料open(20,file='rewriteQS.txt',status='old')read(20,*)do L=1,nQSread(20,*)year,month,day,Q(1),s(1)Q(1)=Q(1)*0.5 do i=2,nn Q(i)=Q(1)enddowrite(*,*)'该组流量和来沙量为:',Q(1),s(1)!推求水面线call level(x,rough,nn,zlevel,dx,Q,npoint,B,A,mm,fai)WRITE(*,*)'累积里程disd（自坝址断面起算

22、） 水位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!计算泥沙沉速和水流挟沙力dxuan=0.05*10*(-3.0)call FALLV(dxuan,w)write(*,*) 'w(悬)=',wk=0.124m=1.05g=9.81do i=1,nnss(i)=k*(v(i)*3.0/(g*h(i)*w)*menddo!计算能坡do i=1,nnJJ(i)=v(i)*2.0*rough(i)*2.0/h(i)*(4.0/3.0)e

23、nddo!计算推移质输沙率dtui=2*10*(-3.0)rough2=dtui*(1.0/6.0)/26.0p=1ps=2.65r=p*grs=ps*g!单位说明：公式中r单位为kN/m3,p单位未为t/m3gb(1)=Q(1)*s(1)*0.05/B(1)do i=2,nnpanduan=(rough2/rough(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)=0endifend

24、dowrite(*,*)'断面推移质输沙率gb'do i=1,nnwrite(*,*)i,gb(i)enddo!计算各断面含沙量、冲淤厚度，并修改断面高程，1断面单独进行!每个断面要进行如下判断：1加上dy是否冲到原始河床以下。2断面高程修改部分应在水位以下。3若不符合条件1，则令dy=0且调整含沙量S(除进口断面外)，若调整后S小于零则S=0afa=0.25p0=1.325*1000dt=86400.0!1断面write(*,*)'断面冲淤变化量dy'dy(1)=afa*w*(s(1)-ss(1)*dt/p0jiancha=minval(x(1:npoint(

25、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),'断面冲至原始河床以下,已进行修改dy=0'endif!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

26、(-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)/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.

27、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),'断面冲至原始河床以下,已进行修改dy=0'S(i+1)=(p0*dx(i+1)*B(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)=0e

28、ndifenddowrite(*,*)'断面含沙量s挟沙力ss'do i=1,nnwrite(*,*)i,s(i),ss(1)enddo!验算淤积守恒性及部分特征值统计结果write(*,*)'该来水来沙条件下河床地形已调整，以下验证淤积守恒'rukushaliang=rukushaliang+Q(1)*s(1)*dtchukushaliang=chukushaliang+Q(31)*S(31)*dtjingliuliang=jingliuliang+Q(1)*dtzongliuliang=zongliuliang+Q(1)zonghansha=zonghans

29、ha+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(*,*)'入库-出库泥沙淤积量'write(*,*)rujianchu,yujiliangwrite(*,*)'以上数据L=&

30、#39;,L!进行年平均流量、含沙量、输沙量、径流量统计!一年一次,1991if(L.eq.365)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1991tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukush

31、aliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1992-366if(L.eq.731)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliuliangopen(100,file='1992tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhan

32、sha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1993if(L.eq.1096)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1993tongji.txt',status='new')write(100,*)'年

33、平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1994if(L.eq.1461)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1994to

34、ngji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1995if(L.eq.1826)thenpjliuliang=zongliuliang/365pjha

35、nsha=zonghansha/zongliuliangopen(100,file='1995tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次

36、,1996-366if(L.eq.2192)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliuliangopen(100,file='1996tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingl

37、iuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1997if(L.eq.2557)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1997tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,y

38、ujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1998if(L.eq.2922)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1998tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(

39、100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,1999if(L.eq.3287)thenpjliuliang=zongliuliang/365pjhansha=zonghansha/zongliuliangopen(100,file='1999tongji.txt',status='new&#

40、39;)write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0endif!一年一次,2000-366if(L.eq.3653)thenpjliuliang=zongliuliang/366pjhansha=zonghansha/zongliulian

41、gopen(100,file='2000tongji.txt',status='new')write(100,*)'年平均流量 含沙量 输沙量 径流量 淤积量'write(100,*)pjliuliang,pjhansha,rukushaliang,jingliuliang,yujiliang2close(100)rukushaliang=0chukushaliang=0jingliuliang=0zongliuliang=0zonghansha=0yujiliang2=0write(*,*)'年统计完毕，进行第11年的计算'en

42、dif!进行冲淤后数据输出和统计部分：深泓和27断面形态，静库容计算!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'

43、;)write(60,*)'27断面高程'do j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file='10yujiliang.txt',status='new')write(80,*)'淤积量入库-出库'write(80,*)yujiliang,rujianchuclose(80)!4open(90,file='10kurong.txt',status='new')write(90,*)'库容'do j=1,2

44、6shuiwei=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!20年7306if(L.eq.7306)thendo i=1,nnzuidi2(i)=minval(x(1:npoint(i),2,i)enddo!1open(70,file='20ShenHong.t

45、xt',status='new')write(70,*)'冲淤后深泓'do i=1,nnwrite(70,*)zuidi2(i)enddoclose(70)!2open(60,file='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,*)'淤积量入库-出库&

46、#39;write(80,*)yujiliang,rujianchuclose(80)!4open(90,file='20kurong.txt',status='new')write(90,*)'库容'do j=1,26shuiwei=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,*)k

47、urong(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)enddoclose(70)!2open(60,file='30mouduanmian.txt',status='new')do

48、j=1,npoint(27)write(60,*)x(j,2,27)enddoclose(60)!3open(80,file='30yujiliang.txt',status='new')write(80,*)'淤积量入库-出库'write(80,*)yujiliang,rujianchuclose(80)!4open(90,file='30kurong.txt',status='new')write(90,*)'库容'do j=1,26shuiwei=j+249kurong(1:26)=0do i

49、=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')w

50、rite(70,*)'冲淤后深泓'do i=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,*)'淤积量入库-出库'write(80,*)yujiliang,rujia

51、nchuclose(80)!4open(90,file='40kurong.txt',status='new')write(90,*)'库容'do j=1,26shuiwei=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

52、年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