module defgpi
    implicit none
    real(8),parameter::g=9.80665D0
    real(8),parameter::pi=3.14159265358979323846D0
end module defgpi
!-----------------------------------------------------------
program f90_STEEP
! Water Surface Profile
!-----------------------------------------------------------
    use defgpi
    implicit none
    real(8),parameter::n=0.013
    real(8),parameter::Q=730.0D0
    real(8),parameter::B=16.0D0
    real(8)::x
    real(8)::y
    real(8)::h
    real(8)::theta0,theta1,theta2
    real(8)::dx
    real(8)::hc
    real(8)::h0
    real(8)::ddx,k1,k2,k3,k4

    write(6,'(6(A10))') 'x','y','theta','h','h0','hc'
    x=0.0D0
    y=FLOOR_Y(x)
    theta0=0.0D0
    hc=(Q*Q/g/B/B/cos(theta0))**(1.0D0/3.0D0)
    h0=-1.0D0
    write(6,'(6(f10.3))') x,y,theta0/pi*180.0D0,hc,h0,hc

    x=0.0D0
    dx=0.001D0
    ddx=0.01D0*dx
    do
        x=x+dx
        theta0=atan((FLOOR_Y(x-ddx)-FLOOR_Y(x))/ddx)
        hc=(Q*Q/g/B/B/cos(theta0))**(1.0D0/3.0D0)
        h0=CALH0(n,Q,B,theta0,hc)
        if(0.0D0<h0)exit
    end do
    y=FLOOR_Y(x)
    write(6,'(6(f10.3))') x,y,theta0/pi*180.0D0,hc,h0,hc

    x=0.05D0
    dx=0.05D0
    ddx=0.01D0*dx
    theta0=atan((FLOOR_Y(x-ddx)-FLOOR_Y(x))/ddx)
    hc=(Q*Q/g/B/B/cos(theta0))**(1.0D0/3.0D0)
    h=hc-0.001
    y=FLOOR_Y(x)
    h0=CALH0(n,Q,B,theta0,hc)
    write(6,'(6(f10.3))') x,y,theta0/pi*180.0D0,hc,h0,hc
    do
        x=x+dx
        y=FLOOR_Y(x)
        theta0=atan((FLOOR_Y(x-ddx)-FLOOR_Y(x))/ddx)
        theta1=atan((FLOOR_Y(x+0.5D0*dx-ddx)-FLOOR_Y(x+0.5D0*dx))/ddx)
        theta2=atan((FLOOR_Y(x+1.0D0*dx-ddx)-FLOOR_Y(x+1.0D0*dx))/ddx)
        k1=FUNC(n,Q,B,x,h,theta0)
        k2=FUNC(n,Q,B,x+0.5D0*dx,h+0.5D0*dx*k1,theta1)
        k3=FUNC(n,Q,B,x+0.5D0*dx,h+0.5D0*dx*k2,theta1)
        k4=FUNC(n,Q,B,x+1.0D0*dx,h+1.0D0*dx*k3,theta2)
        h=h+dx/6.0D0*(k1+2.0D0*k2+2.0D0*k3+k4)
        hc=(Q*Q/g/B/B/cos(theta2))**(1.0D0/3.0D0)
        h0=CALH0(n,Q,B,theta2,hc)
        write(6,'(6(f10.3))') x,y,theta2/pi*180.0D0,h,h0,hc
        if(110.0D0<=x)exit
    end do

stop

contains
!-----------------------------------------------------------
    real(8) function FUNC(n,Q,B,x,h,theta)
    ! Function for Runge-Kutta method
        real(8),intent(in)::n,Q,B,x,h,theta
        real(8)::A,R
        A=B*h
        R=A/(B+2.0D0*h)
        FUNC=(sin(theta)-n*n*Q*Q/R**(4.0D0/3.0D0)/A/A)/(cos(theta)-Q*Q*B/g/A/A/A)
    end function FUNC
!-----------------------------------------------------------
    real(8) function FLOOR_Y(x)
    ! Propilr of Invert
        real(8),intent(in)::x
        real(8)::y
        if(0.0D0<=x.and.x<27.606D0)y=691.500D0-x*x/44.1696D0
        if(27.606D0<=x.and.x<83.199D0)y=708.754D0-1.25D0*x
        if(83.199D0<=x.and.x<97.254D0)y=616.000D0-sqrt(18.0D0*18.0D0-(x-97.254D0)*(x-97.254D0))
        if(97.254D0<=x)y=598.000D0
        FLOOR_Y=y
    end function FLOOR_Y
!-----------------------------------------------------------
    real(8) function CALH0(n,Q,B,theta,hc)
    ! Calculation of normal depth (Bisection method)
        real(8),intent(in)::n,Q,B,theta,hc
        real(8)::f1,f2,fi,h1,h2,hi
        h1=1.0D-6        ; f1=CALF(n,Q,B,h1,theta)
        h2=hc            ; f2=CALF(n,Q,B,h2,theta)
        hi=0.5D0*(h1+h2) ; fi=CALF(n,Q,B,hi,theta)
        do
            if(f1*fi<0.0D0)h2=hi
            if(fi*f2<0.0D0)h1=hi
            if(0.0D0<f1*fi.and.0.0D0<fi*f2)h1=-1.0D0
            if(f1*fi<0.0D0.and.fi*f2<0.0D0)h2=-1.0D0
            if(h1<0.0D0.or.h2<0.0D0)exit
            hi=0.5D0*(h1+h2)
            f1=CALF(n,Q,B,h1,theta)
            f2=CALF(n,Q,B,h2,theta)
            fi=CALF(n,Q,B,hi,theta)
            if(abs(fi)<1.0D-3)exit
        end do
        if(h1<0.0D0.or.h2<0.0D0)hi=-1.0D0
        CALH0=hi
    end function CALH0
!-----------------------------------------------------------
    real(8) function CALF(n,Q,B,h,theta)
    ! Function for Bisection method
        real(8),intent(in)::n,Q,B,h,theta
        real(8)::A,R
        A=B*h
        R=A/(B+2.0D0*h)
        CALF=sin(theta)-n*n*Q*Q/R**(4.0D0/3.0D0)/A/A
    end function CALF
!-----------------------------------------------------------
end program f90_STEEP
!-----------------------------------------------------------