编辑代码

module global!定义模块名
implicit none
integer,parameter::NN=200
integer,parameter::MM=200
character(len=4)bianhao(MM)!储存和操作文本数据，字符长度为4,编号MM=200
real z1(NN,NN),ze(NN,NN)!实数z1?ze?
real zz(MM)
real,parameter::fai=0.5!定义实数类型的常量fai为0.5。parameter表示该常量是一个参数real rough
real dx
real q
end module

!主程序
program main
use global
implicit none

integer i,j,k,1,m,n,p,ni,nj,duanmianshu!指数代表的意思?
real dl,aa,bc,zzm,z1,z2,areal,xx1,rr1,bb1!指数代表的意思?
real juli(MM)

!读入断面数
open(800,file='canshu.txt')!参数文件
read(800,*)
read(800,*)duanmianshu

!读入地形数据
open(500,file='duanmiandixing.txt')!断面地形
do i=1,duanmianshu
read(500,*)bianhao(i)
do j=1,161
read(500,*)zl(i,j),ze(i,j)
end do
end do

!读入参数
read(800,*)
read(800,*)rough!粗糙系数
read(800,*)
read(800,*)dx!断面间距
read(800,*)
read(800,*)q!进口流量
read(800,*)
read(800,*)zz(duanmianshu)!出口水位

!计算每个断面的水位
do i=duanmianshu,2,-1!从断面21开始以编号数步长为1的递减
call level(i,zz(i))
end do

!输出结果
open(600,file='shuimianxian.txt')
write(600,*)'断面编号 距进口距离 水位’
write(*,*)'断面编号  距进口距离  水位’
do i=1,duanmianshu
juli(i)=(i-1)*dx
write(600,*)bianhao(i),juli(i),zz(i)
write(*,*)bianhao(i),juli(i),zz(i)
end do
stop
end

!水面线计算子程序
subroutine level(duanmian,shuiwei)!调用时必须用call命令
use globalimplicit none
real shuiwei,zmin,dz,aa,bb,xx,a,b,x,fmin,dz1,flevel,zmax
real fmax,r,rr
integer duanmian,duanmianl
  dz=0.10
  dz1=0.01
  zmin=shuiwei+dz1duanmianl=duanmian-1
  call area(duanmian,shuiwei,aa,rr)call area(duanmianl,zmin,a,r)
  fmin=flevel(shuiwei,zmin,aa,rr,a,r)
  if(fmin>0) then
    zmax=zmin+dz
    call area(duanmian-1,zmax,a,r)
    fmax=flevel(shuiwei,zmax,aa,rr,a,r)
    do while(fmax>0)
    zmax=zmax+dz
    call area(duanmian-1,zmax,a,r)
    fmax=flevel(shuiwei,zmax,aa,rr,a,r)
    end do
  end if
  call bisec(shuiwei,zmin,zmax,fmin,a,r,aa,rr,duanmian)
  return
end

!水面线函数计算程序
function flevel(z0,zlevel,alower,rlower,au,ru)!必须返回计算结果，伯努利能量方程
  use global
  implicit none
  real z0,zlevel,alower,blower,au,ru,hlower,h
  real aa,b1,b2,bb,c1,c2,cc,flevel,rlower
  !hlower=alower/blower
  !h=au/bu
  aa=z0-zlevel
  bl=fai*(q/alower)**2.0/rlower**(4/3.0)
  b2=(1-fai)*(q/au)**2.0/ru**(4/3.0)
  bb=dx*rough**2.0*(b1+b2)
  cl=(q/alower)**2.0
  c2=(q/au)**2.0
  cc=(cl-c2)/(2*9.8)
  flevel=aa+bb+cc
  return
    end

!二分法求根子程序
subroutine bisec(shuiwei,zmin,zmax,fmin,a,r,aa,rr,duanmian)
use global
implicit none
real err,ddz,ff,zmin,zmax,fmin,a,r,aa,rr,zlevel,f,shuiwei,flevel
integer duanmian

err=0.0001 !精度控制，可以修改。
ddz=zmax-zmin
do while(ddz>=err)
zlevel=0.5*(zmin+zmax)
  call area(duanmian-1,zlevel,a,r)
  f=flevel(shuiwei,zlevel,aa,rr,a,r)
  ff=f*fmin
  if (ff>0)then
    zmin=zlevel
  else
    zmax=zlevel
  end if
  ddz=zmax-zmin
end do
zz(duanmian-1)=zlevel
return
end

!过水断面面积计算子程序
subroutine area(ii,zz1,areal,rr1)
use global
implicit none
integer ii,i
real zz1,areal,bb1,xx1,rr1
real dl,z1,z2,bc
areal=0.0
xx1=0.0
bb1=0.0
do i=1,160
dl=zl(ii,i+1)-zl(ii,i)
z1=max(ze(ii,i),ze(ii,i+1))
z2=min(ze(ii,i),ze(ii,i+1))
if((zl-zz1)*(z2-zz1)>=0)then
  if(z2>=zz1)then
  areal=areal+0.0
  xxl=xx1+0.0
  bbl=bbl+0.0
  else if(zl<=zz1)then
  areal=areal+dl*(zz1-0.5*(ze(ii,i)+ze(ii,i+1)))
  xx1=xxl+sqrt((z1-z2)**2.0+dl**2.0)
  bbl=bbl+dl
  end if
end if
if((z1-zz1)*(z2-zz1)<0)then
  if(zzl<ze(ii,i))then
  bc=dl*(zz1-ze(ii,i+1))/(ze(ii,i)-ze(ii,i+1))
  areal=areal+bc*0.5*(zzl-ze(ii,i+1))
  xxl=xxl+sqrt(bc**2.0+(zzl-ze(ii,i+1))**2.0)
  bbl=bbl+bc
  else if(zz1>ze(ii,i))then
  bc=dl*(zz1-ze(ii,i))/(ze(ii,i+1)-ze(ii,i))
  areal=areal+bc*0.5*(zz1-ze(ii,i))
  xxl=xxl+sqrt(bc**2.0+(zz1-ze(ii,i))**2.0)
  bbl=bbl+bc
  end if
end if
end do
rrl=areal/xx1
return
end