program f90_floodr
    implicit none
    integer,parameter::nnmax=1000
    integer,parameter::ndmax=1000
    integer,parameter::imax=10
    real(8)::rh(1:nnmax),rv(1:nnmax),vcoef   ! Reservoir H-V data
    real(8)::ti(1:ndmax),q_in(1:ndmax),ELini,outlet ! Inflow data
    real(8)::EL,VOL,q_out
    real(8)::hh,dh,eps,q1,q2,Qin,Qout,R,elv,hmax,ELcr,wlmax
    integer::i,j,nn,nd,iUD,io,iii
    character(len=50)::fnameR1,fnameR2,fnameW

    call getarg(1,fnameR1)
    call getarg(2,fnameR2)
    call getarg(3,fnameW)

    ! Input H-V data
    open(11,file=fnameR1,status='old')
        read(11,*) vcoef
        do i=1,nnmax
            read(11,*,iostat=io) rh(i),rv(i)
            rv(i)=rv(i)*vcoef
            if(io<0)exit
        end do
    close(11)
    nn=i-1

    ! Input time sequence of inflow
    open(11,file=fnameR2,status='old')
        read(11,*) ELini,outlet
        do i=1,ndmax
            read(11,*,iostat=io) ti(i),q_in(i)
            if(io<0)exit
        end do
    close(11)
    nd=i-1

    ELcr=ELini
    hmax=-99.9D0
    do iii=1,imax
        ! Initial outflow
        elv=ELini
        EL=elv
        VOL=RET_V(nn,rh,rv,EL)
        elv=RET_H(nn,rh,rv,VOL)
        q_out=SPQ(ELcr,elv,hmax)+outlet
        i=1
        iUD=0
        open(12,file=fnameW,status='replace')
        write(12,'(2(a5),6(a15))') 'i','iUD','time','EL','EL-elv','VOL','q_in','q_out'
        write(12,'(2(i5),6(e15.7))') i,iUD,ti(i),EL,EL-elv,VOL,q_in(i),q_out

        ! Iterative calculation
        wlmax=-99.9D0
        iUD=1
        dh=0.0001
        eps=0.0001
        do i=1,nd-1
            Qin=0.5D0*(q_in(i)+q_in(i+1))*(ti(i+1)-ti(i))*3600.0D0
            hh=0.0D0
            j=0
            do
                j=j+1;
                hh=hh+dble(iUD)*dh
                elv=EL    ;q1=SPQ(ELcr,elv,hmax)+outlet
                elv=EL+hh ;q2=SPQ(ELcr,elv,hmax)+outlet
                Qout=0.5D0*(q1+q2)*(ti(i+1)-ti(i))*3600.0D0
                R=VOL+Qin-Qout
                elv=RET_H(nn,rh,rv,R)
                if(iUD==1.and.j==10)then
                    if(EL+hh>elv)then
                        iUD=-1
                        hh=0.0D0
                        j=0
                    end if
                end if
                if(iUD==-1.and.j==10)then
                    if(EL+hh<elv)then
                        iUD=1
                        hh=0.0D0
                        j=0
                    end if
                end if
                if(abs(EL+hh-elv)<eps)exit
            end do
            VOL=R    ! Cumulative volume
            EL=EL+hh ! Elevation
            q_out=q2 ! Outflow
            write(12,'(2(i5),6(e15.7))') i+1,iUD,ti(i+1),EL,EL-elv,VOL,q_in(i+1),q_out
            if(wlmax<=EL)wlmax=EL
        end do
        close(12)
        write(6,*) iii,hmax,wlmax
        if(abs(ELcr+hmax-wlmax)<0.001D0)exit
        hmax=wlmax-ELcr
    end do

stop

contains

real(8) function SPQ(ELcr,elv,hmax)
    ! To obtain the discharge of overflow crest from the water level
    real(8),intent(in)::ELcr,elv,hmax
    real(8)::b,x,Cd,a,h,q,c,kp,ka,Hd,bx
    integer::n

    b=20.5D0
    n=4
    x=0.1387D0
    Cd=1.971D0+0.498D0*x+6.63D0*x*x
    a=(1.6D0-Cd)/(Cd-3.2D0)
    h=elv-ELcr
    q=0.0D0
    if(0.0<h)then
        if(hmax<0.0D0)Hd=h
        if(0.0D0<hmax)Hd=hmax
        c=1.6D0*(1.0D0+2.0D0*a*h/Hd)/(1.0D0+a*h/Hd)
        kp=0.112703962703963D0*(h/Hd)**2D0-0.237885780885781D0*(h/Hd)+0.126496503496503D0
        ka=-0.0813160731673404D0*(h/Hd)**2D0+0.255143813465017D0*(h/Hd)
        bx=b*dble(n)-2.0D0*(dble(n-1)*kp+ka)*h
        q=c*bx*h**1.5D0
    end if
    SPQ=q
end function SPQ

real(8) function RET_V(nn,rh,rv,elv)
    ! To obtain the cumulative volume from the water level
    integer,intent(in)::nn
    real(8),intent(in)::rh(1:nnmax),rv(1:nnmax),elv
    integer::i
    real(8)::x1,y1,x2,y2,a,b,v

    do i=1,nn-1
        if(rh(i)<=elv.and.elv<=rh(i+1))exit
    end do
    if(rh(nn)<elv)i=nn-1
    x1=rv(i)
    y1=rh(i)
    x2=rv(i+1)
    y2=rh(i+1)
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    v=(elv-b)/a
    RET_V=v
end function RET_V

real(8) function RET_H(nn,rh,rv,v)
    ! To obtain the water level from cumulative volume
    integer,intent(in)::nn
    real(8),intent(in)::rh(1:nnmax),rv(1:nnmax),v
    integer::i
    real(8)::x1,y1,x2,y2,a,b,elv

    do i=1,nn-1
        if(rv(i)<=v.and.v<=rv(i+1))exit
    end do
    if(rv(nn)<v)i=nn-1
    x1=rv(i)
    y1=rh(i)
    x2=rv(i+1)
    y2=rh(i+1)
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*v+b
    RET_H=elv
end function RET_H

end program f90_floodr
