program f90_MARSTON

    implicit none
    integer::i
    real(8),parameter::pi=3.14159265358979323846D0
    real(8)::K
    real(8)::mu
    real(8)::phi
    real(8)::H
    real(8)::Dc1,Dc2
    real(8)::w
    real(8)::WV1,WV2
    real(8)::eps
    real(8)::rsd
    real(8)::P

    phi=30.0D0
    K=(1.0D0-sin(phi/180.0D0*pi))/(1.0D0+sin(phi/180.0D0*pi))
    mu=tan(phi/180.0D0*pi)
    w=1.8D0
    rsd=-0.1D0
    P=1.0D0
    eps=1.0D-6

    Dc1=1.547D0*2.0D0
    Dc2=0.736D0*2.0D0
    do i=1,50
        H=dble(i)*0.1D0
        call MAIN(K,mu,H,Dc1,w,rsd,P,eps,WV1)
        call MAIN(K,mu,H,Dc2,w,rsd,P,eps,WV2)
        write(6,*) H,w*H,WV1,WV2
    end do

stop


contains

    subroutine MAIN(K,mu,H,Dc,w,rsd,P,eps,WV)
        real(8),intent(in)::K,mu,H,Dc,w,rsd,P,eps
        real(8),intent(out)::WV
        real(8)::x1,x2,xx,f1,f2,ff
        integer::i,iflg,nmax
        real(8)::Cc,He

        !Bisection method
        nmax=100
        iflg=1
        x1=0.01D0*H
        x2=H
        do i=1,nmax
            xx=0.5*(x1+x2)
            f1=Func(K,mu,Dc,H,x1,rsd,P)
            f2=Func(K,mu,Dc,H,x2,rsd,P)
            ff=Func(K,mu,Dc,H,xx,rsd,P)
            if(f1*ff<0.0D0)x2=xx
            if(ff*f2<0.0D0)x1=xx
            if(f1*ff<0.0D0.and.f2*ff<0.0D0)exit
            if(f1*ff>0.0D0.and.f2*ff>0.0D0)exit
            if(abs(f1)<eps)iflg=0
            if(abs(f2)<eps)iflg=0
            if(abs(ff)<eps)iflg=0
            if(iflg==0)exit
!            write(6,*) f1,f2,ff
        end do
        if(iflg==0)then
            He=xx
        else
            He=H
        end if

        Cc=0.5D0*(1.0D0-exp(-2.0D0*K*mu*H/Dc))/K/mu
        if(H>He)then
            Cc=0.5D0*(1.0D0-exp(-2.0D0*K*mu*He/Dc))/K/mu+(H/Dc-He/Dc)*exp(-2.0D0*K*mu*He/Dc)
        end if

        WV=Cc*w*Dc   

    end subroutine MAIN


    real(8) function Func(K,mu,Dc,H,He,rsd,P)
        Real(8),intent(in)::K,mu,Dc,H,He,rsd,P
        Func=0.5D0/K/mu*(1.0D0-exp(-2.0D0*K*mu*He/Dc))*(0.5D0/K/mu-(H/Dc-He/Dc)-rsd*P/3.0D0)-0.5D0*(He/Dc)**2.0D0&
        -rsd*P/3.0D0*(H/Dc-He/Dc)*exp(-2.0D0*K*mu*He/Dc)-0.5/K/mu*He/Dc+H/Dc*He/Dc+(rsd*P*H/Dc)
    end function Func


end program f90_MARSTON