program f90_floodr
    implicit none
    integer,parameter::nnmax=1000
    integer,parameter::ndmax=1000
    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
    integer::i,j,nn,nd,iUD,io
    character(len=50)::fnameR1,fnameR2

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

    ! 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

    ! Initial outflow
    elv=ELini
    EL=elv
    VOL=RET_V(nn,rh,rv,EL)
    elv=RET_H(nn,rh,rv,VOL)
    q_out=SPQ(elv)+outlet
    i=1
    iUD=0
    write(6,'(2(a5),6(a15))') 'i','iUD','time','EL','EL-elv','VOL','q_in','q_out'
    write(6,'(2(i5),6(e15.7))') i,iUD,ti(i),EL,EL-elv,VOL,q_in(i),q_out

    ! Iterative calculation
    iUD=1
    dh=0.0001
    eps=0.0005
    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(elv)+outlet
            elv=EL+hh ;q2=SPQ(elv)+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(6,'(2(i5),6(e15.7))') i+1,iUD,ti(i+1),EL,EL-elv,VOL,q_in(i+1),q_out
    end do

stop

contains

real(8) function SPQ(elv)
    ! To obtain the discharge of overflow crest from the water level
    real(8),intent(in)::elv
    real(8)::q,c,b,h,ELcr

    ELcr=120.0D0
    b=180.0D0
    c=2.0D0
    h=elv-ELcr
    q=0.0D0
    if(0.0D0<h)q=c*b*h**(3.0D0/2.0D0)
    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