武汉大学《水力学》梯形渠道恒定非均匀渐变流水面线计算程序WTSFFOR.docx

moduleNrType ! 确定常数implicit noneinteger,parameter::I2B=selected_int_kind(4)integer,parameter::I4B=selected_int_kind(9)integer,parameter::SP=kind(1.0) ! 单精度integer,parameter::DP=kind(1.0d0) ! 双精度real(DP),parameter::ERR=1.0d-6 ! 两实数差值的绝对值小于此常数可认为相等end module NrTypemodule wtsfuse NrTypeimplicit nonesave ! 模块内全局变量均为静态变量integer(I2B),parameter::NOUT=2 ! 输出文件号integer(I2B)::dr ! 水面线计算方向real(DP),parameter::G=9.8d0 ! 重力加速度real(DP)::q,i,n,m,b,csn,srm,alfa! 以上分别为流量,底坡比降,糙率,边坡系数,底宽,底坡余弦值,边坡余割值,动能修正系数real(DP)::v2,j1,j2,es1,es2,ds! 流速,水力比降,断面比能,断面间距real(DP),allocatable::h(:),v(:),s(:) ! 水深,流速,距起始(控制)断面的距离! 外部函数real(DP),external::ErFenFa ! 二分法real(DP),external::zBrent ! 布伦特法containssubroutine ReadData() ! 读取数据! 使用的全局变量有:q,i,n,m,b,NOUTwrite(*,*)键盘输入:流量Q,底坡比降i,糙率n,边坡系数m,底宽b(单位均为国际单位制)read(*,*)q,i,n,m,b ! 输入流量,底坡比降,糙率,边坡系数,底宽if(qERR.or.nERR.or.m0.0d0.or.b0.0d0)then ! m,b可以为0,i可正可负write(*,*)错误:Q,n必须均为正数!m,b必须均为非负数!stopend ifwrite(*,(5x,Q=,f8.2,(m**3/s),4x,i=,f8.5,4x,n=,f6.4,4x,m=,f4.2,4x,b=,f7.2,(m)))q,i,n,m,bwrite(NOUT,(5x,Q=,f8.2,(m**3/s),4x,i=,f8.5,4x,n=,f6.4,4x,m=,f4.2,4x,b=,f7.2,(m)))q,i,n,m,bend subroutine ReadDatasubroutine Cal() ! 计算! 全局变量均使用real(DP)::a,r,hk,h0,hd,hb ! 过水面积,水力半径,临界水深,正常水深,控制水深,区间端点水深integer(I2B)::l,ns,j ! ns为区间总数,j为计算方法alfa=1.05d0 ! 动能修正系数csn=sqrt(1.0d0-i*i) ! 底坡余弦值srm=sqrt(1.0d0+m*m) ! 边坡余割值call SubHk(hk) ! 调用子程序SubHk求临界水深write(*,(5x,临界水深Hk=,f10.6,(m)))hkwrite(NOUT,(5x,临界水深Hk=,f10.6,(m)))hk! 计算正常水深H0if(iERR)then ! i=0即平坡或逆坡h0=100.0d0 ! 正常水深H0设为一大值write(*,*)i=0,没有正常水深write(NOUT,*)i=0,没有正常水深 else call SubH0(h0) ! 调用子程序SubHk求正常水深write(*,(5x,正常水深H0=,f10.6,(m)))h0write(NOUT,(5x,正常水深H0=,f10.6,(m)))h0endifwrite(*,*)键盘输入:控制水深Hd,断面间距Ds,区间总数Ns(断面总数为Ns+1)read(*,*)hd,ds,ns ! 输入控制水深,断面间距,区间总数if(hdERR.or.dsERR.or.ns=0)thenwrite(*,*)错误:Hd,Ds必须均为非负数!Ns必须为正整数!stopend if! 断面总数为ns+1allocate(h(ns+1),v(ns+1),s(ns+1)) ! 动态数组分配内存h=0.0d0;v=0.0d0;s=0.0d0 ! 赋初值! 计算起始(控制)断面参数h(1)=hd ! 第一断面为控制断面a=(b+m*h(1))*h(1) ! 过水面积r=a/(b+2.0d0*h(1)*srm) ! 水力半径v(1)=q/a ! 流速j1=(n*v(1))**2/r**(4.0d0/3.0