module defpara
    implicit none
    integer::INTDIV
! INTDIV : number of pertial division of time increment for nemerical integration
end module
!**************************************************/
program f90_SFMQ
!**************************************************/
!---------------------------------------------------------------------------*/
!   Storage Function Mehod
!---------------------------------------------------------------------------*/
! nd        : Number of observed data
! md        : Divided number of time increment
! n1        : Starting point of direct run-off
! CA        : Catchment area
! dt        : Time increment of observed data
! ddt       : Time increment for calculation
! R_org()   : Observed rainfall
! RO()      : Rainfall for calculation
! RE        : Effective rainfall
! Qbase     : Base run-off (constant value)
! QB()      : Base run-off (constant value)
! QE()      : Estimated flow
! ff        : Run-off percentage
! kk, pp    : SS=kk*QC**pp (kk,pp : coefficient)
! LT        : Number of delay time (delay time TL = LT*dt)
!
! subroutine RUNGE()    : Numerical integration by Runge-Kutta method
! real(8) function FQ() : Function for Runge-Kutta method
!---------------------------------------------------------------------------*/
    use defpara
    implicit none
    integer,parameter::NDmax=1000
    integer::i,j,io,nnn,nn,iii
    integer::nd,nt,md,nh
    real(8)::CA,dt,ddt,Qbase,TL
    real(8),allocatable::R_org(:)
    real(8),allocatable::RO(:)
    real(8),allocatable::QB(:),QE(:)
    real(8),allocatable::QQ(:,:)
    real(8)::ff,kk,pp,ffe
    integer::n1,LT

    character(len=50)::fnameR,fnameW
    character(len=100)::strcom
    character(len=100)::dat
    character(len=50)::dummy1,dummy2

    ! Reading from command line
    call getarg(1,fnameR) !Input data file
    call getarg(2,fnameW) !Output data file
    call getarg(3,dummy1)
    read(dummy1,*) INTDIV  !INTDIV: pertial division of time increment

    ! Data input
    allocate(R_org(1:NDmax))
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom   !River, location, date
    read(11,*) dummy1,CA     !Catchment Area (km2)
    read(11,*) dummy1,Qbase  !Constant base flow (m3/s)
    read(11,*) dummy1,dt     !Time increment (hour)
    read(11,*) dummy1,n1,TL,ff,pp,kk !Parameter
    read(11,'()')
    do i=1,NDmax
            read(11,*,iostat=io) dummy1,dummy1,R_org(i)
            if(io<0)exit
    end do
    close(11)
    nd=i-1
    md=INTDIV

    nt=(nd-1)*md+1
    n1=(n1-1)*md+1
    ddt=dt/dble(md)
    allocate(RO(1:nt))
    do i=0,nd-2
        do j=1,md
            RO(i*md+j+1)=R_org(i+2)
        end do
    end do
    RO(1)=R_org(1)
    RO(nt)=R_org(nd)
    allocate(QB(1:nt),QE(1:nt))
    do i=1,nt
        QB(i)=Qbase*3.6D0/CA
    end do

    nh=2
    allocate(QQ(1:nh,1:nd))

    LT=int(TL/ddt)
    do iii=1,nh
        if(iii==1)call RUNGE(nt,n1,LT,ddt,ff,kk,pp,RO,QE)
        if(iii==2)call RUNGE(nt,1,0,ddt,ff,kk,pp,RO,QE)
        do i=1,nt
            if(QE(i)<0.0D0)QE(i)=0.0D0
        end do
        ffe=ROF(nt,n1,nt,LT,RO,QE)
        do i=1,nt
            QE(i)=(QE(i)+QB(i))*CA/3.6D0
        end do
        do i=1,nd
            QQ(iii,i)=QE(md*(i-1)+1)
        end do
    end do
    open(12,file=fnameW,status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(a)') 'nd,n1,md,time,pp,kk,ff,ffe'
    write(12,'(i5,2(",",i5),5(",",e15.7))') (nt-1)/md+1,(n1-1)/md+1,md,ddt*LT,pp,kk,ff,ffe
    write(12,'(a)') 'i,Time,RO,QLT,Q0'
    do i=1,nd
        write(12,'(i5,",",f5.1,",",f7.1,2(",",f7.1))') i,dt*(i-1),R_org(i),(QQ(iii,i),iii=1,nh)
    end do
    close(12)




stop

contains
!---------------------------------------------------------------------------*/
! Ratio of Flow
!---------------------------------------------------------------------------*/
    real(8) function ROF(nt,n1,n2,LT,RO,QE)
        integer,intent(in)::nt,n1,n2,LT
        real(8),intent(in)::RO(1:nt),QE(1:nt)

        integer::j
        real(8)::w0,w1

        w0=0.0D0
        w1=0.0D0
        do j=n1,n2
            w0=w0+QE(j)
            if(1<=j-LT)then
                w1=w1+RO(j-LT)
            else
                w1=w1
            end if
        end do
write(6,*) w0,w1
        ROF=w0/w1
    end function ROF
!---------------------------------------------------------------------------*/
! Numerical integration by Runge-Kutta method
!---------------------------------------------------------------------------*/
    subroutine RUNGE(nt,n1,LT,dt,ff,kk,pp,RO,QE)
        integer,intent(in)::nt,n1,LT
        real(8),intent(in)::dt,ff,kk,pp,RO(1:nt)
        real(8),intent(out)::QE(1:nd)
        integer::i
        real(8)::fk1,fk2,fk3,fk4,rf,qq,q0
        real(8)::rf1,rf2

        do i=1,nt
            QE(i)=0.0D0
        end do
        do i=1,nt-1
            if(i<LT)then
                rf1=0.0D0
                rf2=0.0D0
            end if
            if(LT==i)then
                rf1=0.0D0
                rf2=ff*RO(1)
            end if
            if(LT+1<=i)then
                rf1=ff*RO(i-LT)
                rf2=ff*RO(i+1-LT)
            end if
            rf=rf1
            q0=QE(i)
            if(i==n1)q0=0.01D0
            qq=q0
                fk1=FQ(kk,pp,rf,qq)
            rf=0.5*(rf1+rf2)
            qq=q0+0.5D0*dt*fk1
                fk2=FQ(kk,pp,rf,qq)
            rf=0.5D0*(rf1+rf2)
            qq=q0+0.5D0*dt*fk2
                fk3=FQ(kk,pp,rf,qq)
            rf=rf2
            qq=q0+dt*fk3
                fk4=FQ(kk,pp,rf,qq)
            QE(i+1)=q0+1.0D0/6.0D0*dt*(fk1+2.0D0*fk2+2.0D0*fk3+fk4)
        end do

    end subroutine RUNGE


    real(8) function FQ(kk,pp,rf,qq)
        real(8),intent(in)::kk,pp,rf,qq
        if(1.0D-10<qq)then
            FQ=(rf-qq)*qq**(1.0D0-pp)/kk/pp
        else
            FQ=0.0D0
        end if
    end function FQ
!---------------------------------------------------------------------------*/
end program f90_SFMQ
!---------------------------------------------------------------------------*/