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

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