求助！！！防波堤设计计算，波浪折射算角度的程序 - Fortran Coder

[Fortran] 纯文本查看 复制代码
program main

implicit none

integer::r,i,j,m,n,k,l,h

real::T0,s,g,pi,SL,dalpha1,dalpha2,depth,L1,L0,c0,c1,p,q,beta,gamma,Lsy,theta1,theta2,theta0

!!!!!导入水深矩阵

real,dimension(0:16,0:8)::d

real,dimension(0:1000)::x,y

real,dimension(1:1000)::ZZJG

open(1,file='matrix5.dat')

do j=0,8

    do i=0,16

  read(1,*),d(i,j)

  write (*,*) d(i,j)

enddo

enddo

close(1)

!!!!

!!!起点各种参数

print *,"please enter:T0,start position m=，n=,入射角theta0="

read *,T0,m,n,s

g=9.81

pi=3.1415926

SL=200.0

L0=g*(T0**2)/(2*pi)

c0=L0/T0

x(0)=m*400.0

y(0)=n*600.0  

theta0=s

!!!从起点算至工程点

do r=0,1000

    theta0=theta0

    x(r+1)=x(r)+cosd(theta0)*SL                                                        !!!!!!!!!!!!!!!theta是角度

    y(r+1)=y(r)+sind(theta0)*SL

    p=x(r+1)/400.0

    q=y(r+1)/600.0

!!!!!判断区间

k=0

l=0

    do i=0,16

      k=k

   beta=p-k

    if (beta>m) then

       k=i+1

    else 

    endif

    enddo

    do j=0,8

      l=l

   gamma=q-l

     if (gamma>n) then

       l=j+1

     else

  endif

    enddo

!!!!!

         dalpha1=abs((d(k-1,l)-d(k-1,l-1)))*(q-(l-1))+min(d(k-1,l-1),d(k-1,l))

   dalpha2=abs((d(k,l)-d(k,l-1)))*(q-(l-1))+min(d(k,l-1),d(k,l))

   depth=abs(dalpha1-dalpha2)*(p-(k-1))+min(dalpha1,dalpha2)

   if (depth<=1.2) then 

   go to 33

      else if (depth>1.2) then

!!!!!水深插值

!!!!!

                         L1=Lsy(depth)

                   c1=L1/T0

        !角度

                theta1=theta0+((y(r+1)-y(r))*2/(c1+c0))*((c1-c0)/(x(r+1)-x(r))-(1/tan(theta0))*((c1-c0)/(y(r+1)-y(r)))) !角度

                   ZZJG(r+1)=theta1 !角度

                theta0=theta1

       c0=c1 

   ENDIF

   

!!!!!角度计算

enddo

33 do h=1,1000

print *,"各次角度",ZZJG(h)

enddo

endprogram



function Lsy(dd) result(finals)

implicit none

integer::t

real::eps,dd,temp,finals,L0,g,T0,pi

real,dimension(0:1000)::Lr

eps=0.01

Lr(0)=L0

do t=0,1000

Lr(t+1)=g*(T0**2)*tanh(2*pi*dd/L0)/(2*pi)

temp=Lr(t+1)-Lr(t)

if ((abs(temp))>=eps) then

L0=Lr(t+1)

else if (abs(temp)<eps) then

finals=Lr(t+1)

endif

end do

endfunction