新人求助用龙库塔迭代求解水面线，程序已编完无法运行 - Fortran Coder

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

program main

    dimension x(0:20),y(0:20),z(0:20),p(0:21),r(0:20),c(0:20)

    double precision x,y,z,p,h,x0,y0,z0,p0,n,k1,k2,k3,k4,k,q,g,a,b,m,s,r,c

    k=0.009

    g=9.8

    s=0.022

    m=1.5

    a=20

    b=10

    q=45

    r=((b+m*y)*y)/(b+2*sqrt(1+m**2)*y)

    c=(1/s)*r**(1/6)

    f(x,y)=(k-q**2/(a*c*sqrt(r))**2)/(1-q**2*b/(g*a**3)

    h=30

    n=20

    x(0)=0.0

    y(0)=1.88

    r(0)=1.44

    c(0)=363.6

    do i=0,n-1

         r(i+1)=((b+m*y(i+1))*y(i+1))/(b+2*sqrt(1+m**2)*y(i+1))

         c(i+1)=(1/s)*r(i+1)**(1/6)

         x(i+1)=(i+1)*h

k1=f(x(i),y(i))

k2=f(x(i)+0.5*h,y(i)+0.5*h*k1)

k3=f(x(i)+0.5*h,y(i)+0.5*h*k2)

k4=f(x(i)+h,y(i)+h*k3)

y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6

    end do

    p(0)=0

    do i=0,n

        z(i)=p(i)+exp(-p(i))-1

        p(i+1)=p(i)+h

    end do

    open (10,file="1,txt")

    write(10,*) (x(i),y(i),z(i),i=0,20)

    end