module defpara
    implicit none
    integer::INTDIV
    integer,parameter::MAXDIV=20
    real(8),parameter::ADDMAX=0.0D0
! INTDIV : number of pertial division of time increment for nemerical integration
! MAXDIV : number of division from zero to maximum run-off
! ADDMAX : coefficient of additional number at maximum run-off
!          (arrange the weight for regression analysis)
end module
!**************************************************/
program f90_SFM
!**************************************************/
!---------------------------------------------------------------------------*/
!   Storage Function Mehod
!---------------------------------------------------------------------------*/
    use defpara
    implicit none
    integer,parameter::NDmax=1000
    integer::i,j,io,nnn,nn
!---------------------------------------------------------------------------*/
! nd        : Number of observed data
! md        : Divided number of time increment
! n1        : Starting point of direct run-off
! n2        : Ending point of direct run-off
! CA        : Catchment area
! dt        : Time increment of observed data
! ddt       : Time increment for calculation
! R_org()   : Observed rainfall
! Q_org()   : Observed discharge
! RO()      : Rainfall for calculation
! QO()      : Discharge for calculation
! RE        : Effective rainfall
! QC()      : Direct run-off
! QB()      : Base run-off
! SS()      : Storage function
! QE()      : Estimated flow
! ff        : Run-off percentage
! kk, pp    : SS=kk*QC**pp (kk,pp : coefficient)
! LT        : Number of delay time (delay time = LT*dt)
!
! subroutine QSPRT()    : print out of q-s relationship
! subroutine DEFQCQB()  : definition of QC and QB
! subroutine CALC()     : control rooutine for calculation of LT, ff, pp, kk
! subroutine LSM()      : regression analysis for calculation of pp, kk
! subroutine CALAREA()  : calculation of LT, ff, pp, kk
! subroutine RUNGE()    : Runge-Kutta numerical integration
! real(8) function FQ() : function for Runge-Kutta method
! subroutine SPL3()     : spline interpolation
! subroutine STEGJ()    : solver of simultaneous equations for spline interpolation
!---------------------------------------------------------------------------*/
    integer::nd,nt,n1,n2,md,mm,num1,num2
    real(8)::CA,dt,ddt
    real(8),allocatable::R_org(:),Q_org(:)
    real(8),allocatable::RO(:),QO(:)
    real(8),allocatable::QC(:),QB(:),SS(:),QE(:)
    real(8),allocatable::datax(:),datay(:)
    real(8)::ff,kk,pp
    real(8)::Qmax,SQ,sumQO,sumQE
    real(8)::Qres1,Qres2
    integer::LT,iQmax
    real(8)::x1,y1,x2,y2,xx,w1,w2,qq,xmax

    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,*) n1      !starting point of direct run-off
    call getarg(4,dummy1)
    read(dummy1,*) mm      !mm=0: base flow=QO(n1)
    call getarg(5,dummy1)
    read(dummy1,*) INTDIV  !INTDIV: pertial division of time increment

    ! Data input
    allocate(R_org(1:NDmax),Q_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,dt     !Time increment (hour)
    read(11,'()')
    do i=1,NDmax
            read(11,*,iostat=io) dummy1,dummy1,R_org(i),Q_org(i)
            if(io<0)exit
    end do
    close(11)
    nd=i-1
    md=INTDIV
    allocate(datax(1:nd),datay(1:nd))

    n1=(n1-1)*md+1
    nt=(nd-1)*md+1
    ddt=dt/dble(md)
    allocate(RO(1:nt),QO(1:nt))
    do i=0,nd-2
        do j=1,md
            RO(i*md+j+1)=R_org(i+2)
        end do
    end do
    do i=1,nd
        datax(i)=dt*dble(i)
        datay(i)=Q_org(i)
    end do
    call SPL3(nd,datax,datay,md)
    open(11,file='_spline.txt',status='old')
    do i=1,nt
        read(11,*) dummy1,QO(i)
        if(QO(i)<0.0D0)QO(i)=0.0D0
    end do
    close(11)
    RO(1)=R_org(1)
    QO(1)=Q_org(1)
    RO(nt)=R_org(nd)
    QO(nt)=Q_org(nd)
    deallocate(R_org,Q_org)
    deallocate(datax,datay)
    allocate(QC(1:nt),QB(1:nt),SS(1:nt),QE(1:nt))

    Qmax=QO(1)
    do i=1,nt
        if(Qmax<QO(i))then
            Qmax=QO(i)
            iQmax=i
        end if
    end do

    ! Search optimal n2
    if(mm/=0)then
        num1=0
        num2=0
        Qres1=1.0D30
        Qres2=1.0D30
        do nnn=0,nt/3
            n2=nt+1-nnn
            call DEFQCQB(nt,n1,n2,CA,QO,QC,QB)
            call CALC(fnameW,nt,md,n1,n2,ddt,RO,QC,SS,ff,kk,pp,LT,0)
            if(LT<0)exit
            call RUNGE(nt,n1,LT,ddt,ff,kk,pp,RO,QE)
            SQ=0.0D0
            sumQO=0.0D0
            sumQE=0.0D0
            do i=1,nt
                if(QE(i)<0.0D0)QE(i)=0.0D0
                QE(i)=(QE(i)+QB(i))*CA/3.6D0
                QB(i)=QB(i)*CA/3.6D0
                if(n1<=i.and.i<=n2)then
                    SQ=SQ+abs(QE(i)*QE(i)-QO(i)*QO(i))/Qmax/(n2-n1)
                    sumQO=sumQO+QO(i)
                    sumQE=sumQE+QE(i)
                end if
            end do
            if(SQ<Qres1)then
                Qres1=SQ
                num1=nnn
            end if
            if(abs(sumQO-sumQE)<Qres2)then
                Qres2=abs(sumQO-sumQE)
                num2=nnn
            end if
!            write(6,*) '*',nnn,Qres2,num2
        end do
        nnn=min(num1,num2)
    else
        nnn=0
    end if

    ! Finalize
    n2=nt+1-nnn
    call DEFQCQB(nt,n1,n2,CA,QO,QC,QB)
    call CALC(fnameW,nt,md,n1,n2,ddt,RO,QC,SS,ff,kk,pp,LT,1)
    call RUNGE(nt,n1,LT,ddt,ff,kk,pp,RO,QE)
    do i=1,nt
        if(QE(i)<0.0D0)QE(i)=0.0D0
        QE(i)=(QE(i)+QB(i))*CA/3.6D0
        QB(i)=QB(i)*CA/3.6D0
    end do

    ! Print out of result
    open(12,file=fnameW,status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(a)') 'nd,n1,n2,md,LT,time,ff,pp,kk'
    write(12,'(i5,4(",",i5),4(",",e15.7))') (nt-1)/md+1,(n1-1)/md+1,(n2-1)/md+1,md,LT,ddt*LT,ff,pp,kk
    write(12,'(a)') 'i,Time,RO,QO,QE,QB,QE-QB'
    do i=1,nt,md
        write(12,'(i5,",",f5.1,5(",",f7.1))') (i-1)/md+1,ddt*(i-1),RO(i),QO(i),QE(i),QB(i),QE(i)-QB(i)
    end do
    close(12)

    ! print out of Q-S relationship
    dummy1=fnameW(1:len_trim(fnameW)-4)
    dummy2='0final'
    call QSPRT(dummy1,dummy2,LT,0.0D0,ddt,ff,pp,kk,nt,md,QC,SS)

stop

contains
!---------------------------------------------------------------------------*/
!   Print out of Q-S relationship
!---------------------------------------------------------------------------*/
    subroutine QSPRT(dummy1,dummy2,iii,area,dt,ff,pp,kk,nt,md,QC,SS)
        character(len=50),intent(in)::dummy1,dummy2
        integer,intent(in)::iii,nt,md
        real(8),intent(in)::QC(1:nt),SS(1:nt)
        real(8),intent(in)::area,dt,ff,pp,kk

        integer::k,nn
        real(8)::xmax
        character(len=100)::filename

        ! Original data plot
        filename='_'//trim(adjustl(dummy1))//'_qcss0_'//trim(adjustl(dummy2))//'.txt'
        open(13,file=filename,status='replace')
        write(13,*) iii,area
        do k=1,nt,md
            write(13,*) QC(k),SS(k)
        end do
        close(13)
        ! Estimated data plot
        xmax=0.0D0
        do k=1,nt
            if(xmax<QC(k))xmax=QC(k)
        end do
        nn=50
        filename='_'//trim(adjustl(dummy1))//'_qcss1_'//trim(adjustl(dummy2))//'.txt'
        open(13,file=filename,status='replace')
        write(13,*) dble(iii)*dt,ff,pp,kk
        do k=0,nn
            qq=xmax/dble(nn)*dble(k)
            write(13,*) qq,kk*qq**pp
        end do
        close(13)
    end subroutine QSPRT
!---------------------------------------------------------------------------*/
!   Define QC, QB
!---------------------------------------------------------------------------*/
    subroutine DEFQCQB(nt,n1,n2,CA,QO,QC,QB)
        integer,intent(in)::nt,n1,n2
        real(8),intent(in)::CA
        real(8),intent(in)::QO(1:nt)
        real(8),intent(out)::QC(1:nt),QB(1:nt)
        real(8)::q1,q2,qq
        integer::i

        q1=3.6D0*QO(n1)/CA
        if(nt<n2)then
            do i=1,nt
                qq=3.6D0*QO(i)/CA
                QB(i)=qq
                if(q1<=qq)QB(i)=q1 
                QC(i)=qq-QB(i)
            end do
        else
            q2=3.6D0*QO(n2)/CA
            do i=1,nt
                qq=3.6D0*QO(i)/CA
                QB(i)=qq
                if(n1<=i.and.i<=n2)QB(i)=q1+dble(i-n1)*(q2-q1)/dble(n2-n1)
                QC(i)=qq-QB(i)
            end do
        end if
        do i=1,nt
            if(QC(i)<=0.0D0)QC(i)=0.0D0
            if(QB(i)<=0.0D0)QB(i)=0.0D0
        end do

    end subroutine DEFQCQB
!---------------------------------------------------------------------------*/
!   Calculation of ff, pp, kk, LT
!---------------------------------------------------------------------------*/
    subroutine CALC(fnameW,nt,md,n1,n2,dt,RO,QC,SS,ff,kk,pp,LT,iprt)
        integer,intent(in)::nt,md,n1,n2,iprt
        real(8),intent(in)::dt
        real(8),intent(in)::RO(1:nt),QC(1:nt)
        character(len=50),intent(in)::fnameW
        real(8),intent(out)::SS(1:nt)
        real(8),intent(out)::ff,kk,pp
        integer,intent(out)::LT

        integer::i,j,k,nlsm,nn,nn2
        real(8)::vecx(1:nt),vecy(1:nt),SW(1:nt)
        real(8)::area1,area2,ff1,ff2,pp1,kk1,pp2,kk2
        real(8)::rf,qq,xmax,w0,w1
        character(len=50)::dummy1,dummy2

        nn2=n2
        if(nt<n2)nn2=nt
        ff1=1.0D0
        pp1=0.0D0
        kk1=0.0D0
        area1=1.0D30
        do i=0,nn2
            do j=1,nt
                SS(j)=0.0D0
            end do
            w0=0.0D0
            w1=0.0D0
            do j=n1,nn2
                w0=w0+QC(j)
                if(1<=j-i)then
                    w1=w1+RO(j-i)
                else
                    w1=w1
                end if
            end do
            ff2=w0/w1
            do j=n1,nn2-1
                if(1<=j+1-i)then
                    rf=ff2*RO(j+1-i)
                else
                    rf=0.0D0
                end if
                qq=0.5D0*(QC(j)+QC(j+1))
                SS(j+1)=SS(j)+(rf-qq)*dt
            end do
            xmax=0.0D0
            do j=1,nt
                if(SS(j)<0.0D0)SS(j)=0.0D0
                if(xmax<QC(j))xmax=QC(j)
                vecx(j)=QC(j)
                vecy(j)=SS(j)
            end do
            call CALAREA(nt,md,vecx,vecy,nlsm,pp2,kk2,area2)
            if(pp2<=0.0D0)exit
            if(nlsm<=1)exit
            if(iprt==1.and.mod(i,md)==0)then
                dummy1=fnameW(1:len_trim(fnameW)-4)
                write(dummy2,'(i2.2)') i/md
                call QSPRT(dummy1,dummy2,i,area2,dt,ff2,pp2,kk2,nt,md,QC,SS)
            end if
            if(1<=i.and.area1<area2)exit
            if(pp2<=0.0D0.or.kk2<=0.0D0)exit
            if(1<=i.and.0.3D0<abs(ff1-ff2)/ff1)exit
            ff1=ff2
            pp1=pp2
            kk1=kk2
            area1=area2
            do k=1,nt
                SW(k)=SS(k)
            end do
        end do
        ff=ff1
        pp=pp1
        kk=kk1
        LT=i-1
        do i=1,nt
            SS(i)=SW(i)
        end do
        if(nlsm<=1)LT=-99
    end subroutine CALC
!---------------------------------------------------------------------------*/
! Regression (y=aa*x+bb)
!---------------------------------------------------------------------------*/
    subroutine LSM(nn,vx,vy,pp,kk)
        integer,intent(in)::nn
        real(8),intent(in)::vx(1:nn),vy(1:nn)
        real(8),intent(out)::pp,kk
        integer::i
        real(8)::x1,y1,x2,y2,xy,xm,ym,aa,bb
    
        x1=0.0D0
        y1=0.0D0
        x2=0.0D0
        xy=0.0D0
        do i=1,nn
            x1=x1+vx(i)
            y1=y1+vy(i)
            x2=x2+vx(i)*vx(i)
            xy=xy+vx(i)*vy(i)
        end do
        xm=x1/dble(nn)
        ym=y1/dble(nn)
        aa=(dble(nn)*xy-x1*y1)/(dble(nn)*x2-x1*x1)
        bb=(x2*y1-x1*xy)/(dble(nn)*x2-x1*x1)
        pp=aa
        kk=exp(bb)

        if(1.0D0<pp)then
            x2=0.0D0
            xy=0.0D0
            do i=1,nn
                x2=x2+exp(vx(i))*exp(vx(i))
                xy=xy+exp(vx(i))*exp(vy(i))
            end do
            pp=1.0D0
            kk=xy/x2
        end if

    end subroutine LSM
!---------------------------------------------------------------------------*/
! Calculation of nlsm,pp,kk,area
!---------------------------------------------------------------------------*/
    subroutine CALAREA(nt,md,vecx,vecy,nlsm,pp,kk,area)
        integer,intent(in)::nt,md
        real(8),intent(in)::vecx(1:nt),vecy(1:nt)
        real(8),intent(out)::pp,kk,area
        integer,intent(out)::nlsm

        integer::nnd
        integer::i,ixmax,j,k,m,nk,nadd
        real(8)::dumx(1:nt),dumy(1:nt)
        real(8)::dx,x1,x2,q1,s1,q2,s2,qqm,ssm
        real(8),allocatable::datax(:),datay(:)
        real(8),allocatable::datx(:),daty(:)
        character(len=50)::filename

        nnd=MAXDIV
        allocate(datax(1:2*nnd),datay(1:2*nnd))

        ixmax=1
        qqm=0.0D0
        do i=1,nt
            dumx(i)=vecx(i)
            dumy(i)=vecy(i)
            if(qqm<dumx(i))then
                qqm=dumx(i)
                ixmax=i
            end if
        end do
        qqm=dumx(ixmax)
        ssm=dumy(ixmax)

        dx=qqm/dble(nnd)
        k=0
        do i=1,nnd
            x1=dx*dble(i-1)
            x2=dx*dble(i)
            q1=0.0D0
            q2=0.0D0
            s1=0.0D0
            s2=1.0D10
            m=0
            do j=1,nt
                if(x1<dumx(j).and.dumx(j)<=x2)then
                    m=m+1
                    if(s1<dumy(j))then
                        q1=dumx(j)
                        s1=dumy(j)
                    end if
                    if(dumy(j)<s2)then
                        q2=dumx(j)
                        s2=dumy(j)
                    end if
                end if
            end do
            if(0<m)then
                k=k+1
                datax(k)=q1
                datay(k)=s1
                if(q1/=q2.or.s1/=s2)then
                k=k+1
                datax(k)=q2
                datay(k)=s2
                end if
            end if
        end do
        nk=k
        k=0
        do i=1,nk
            if(0.0D0<datax(i).and.0.0D0<datay(i))then
                k=k+1
                dumx(k)=log(datax(i))
                dumy(k)=log(datay(i))
            end if
        end do
        nlsm=k
        if(2<=nlsm)then
            nadd=int(ADDMAX*nlsm)
            nk=nlsm+nadd
            allocate(datx(1:nk),daty(1:nk))
            do i=1,nlsm
                datx(i)=dumx(i)
                daty(i)=dumy(i)
            end do
            do i=nlsm+1,nk
                datx(i)=log(qqm)
                daty(i)=log(ssm)
            end do
            call LSM(nk,datx,daty,pp,kk)
        else
            pp=-99.0D0
            kk=-99.0D0
        end if

        area=0.0D0
        if(0.0D0<pp)then
            do i=1,nt
                area=area+(vecy(i)-kk*vecx(i)**pp)**2.0D0
            end do
        end if

    end subroutine CALAREA
!---------------------------------------------------------------------------*/
! 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.5D0*(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
!---------------------------------------------------------------------------*/
! Spline interpolation
!---------------------------------------------------------------------------*/
    subroutine SPL3(nd,datax,datay,md)
        integer,intent(in)::nd,md
        real(8),intent(in)::datax(1:nd),datay(1:nd)

        integer::i,j,nn,mm
        real(8)::a(1:nd),b(1:nd),c(1:nd),d(1:nd)
        real(8)::h0,h1,hj
        real(8)::dx,xx,sj
        real(8),allocatable::ary(:,:),vec(:)

        nn=nd-2
        mm=2
        allocate(ary(1:nn,1:mm),vec(1:nn))
        do i=1,nn
            do j=1,mm
                ary(i,j)=0.0D0
            end do
            vec(i)=0.0D0
        end do
        do i=1,nd
            a(i)=datay(i)
            b(i)=0.0D0
            c(i)=0.0D0
            d(i)=0.0D0
        end do
        do i=2,nd-1
            h0=datax(i)-datax(i-1)
            h1=datax(i+1)-datax(i)
            ary(i-1,1)=2.0D0*(h0+h1)
            ary(i-1,2)=h1
            vec(i-1)=3.0D0*(a(i+1)-a(i))/h1-3.0D0*(a(i)-a(i-1))/h0
        end do

        call CHBAND(nn,mm,ary,vec)

        c(1)=0.0D0
        do i=2,nd-1
            c(i)=vec(i-1)
        end do
        c(nd)=0.0D0
        do i=1,nd-1
            hj=datax(i+1)-datax(i)
            d(i)=(c(i+1)-c(i))/(3.0D0*hj)
            b(i)=(a(i+1)-a(i))/hj-(2.0D0*c(i)+c(i+1))*hj/3.0D0
        end do

        open(12,file='_spline.txt',status='replace')
        do i=1,nd-1
            dx=(datax(i+1)-datax(i))/dble(md)
            do j=1,md
                xx=dble(j-1)*dx
                sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0
                write(12,*) datax(i)+xx,sj
            end do
        end do
        i=nd-1
        j=md+1
        xx=dble(j-1)*dx
        sj=a(i)+b(i)*xx+c(i)*xx**2.0D0+d(i)*xx**3.0D0
        write(12,*) datax(i)+xx,sj
        close(12)

    end subroutine SPL3


    subroutine CHBAND(n1,m1,array,vector)
        integer,intent(in)::n1,m1
        real(8),intent(inout)::array(1:n1,1:m1)
        real(8),intent(inout)::vector(1:n1)
        integer :: i,j,k,mm,mn
        real(8) :: z

        array(1,1)=sqrt(array(1,1))
        do j=2,m1
            array(1,j)=array(1,j)/array(1,1)
        end do
        do i=2,n1
            if(n1-i+1>m1)then
                mm=m1
            else
                mm=n1-i+1
            end if
            do j=1,mm
                z=array(i,j)
                if(j/=m1)then
                    if(m1-j+1<i)then
                        mn=m1-j+1
                    else
                        mn=i
                    end if
                    do k=2,mn
                        z=z-array(i-k+1,k)*array(i-k+1,j+k-1)
                    end do
                    if(j>1)then
                        array(i,j)=z/array(i,1)
                    else
                        array(i,1)=sqrt(z)
                    end if
                else
                    array(i,j)=z/array(i,1)
                end if
            end do
        end do

        vector(1)=vector(1)/array(1,1)
        do i=2,n1
            z=vector(i)
            if(i>m1)then
                mm=m1
            else
                mm=i
            end if
            do j=2,mm
                z=z-array(i-j+1,j)*vector(i-j+1)
            end do
            vector(i)=z/array(i,1)
        end do

        vector(n1)=vector(n1)/array(n1,1)
        do i=n1-1,1,-1
            z=vector(i)
            do j=2,m1
                if(i+j-1>n1)exit
                z=z-array(i,j)*vector(i+j-1)
            end do
            vector(i)=z/array(i,1)
        end do
    end subroutine CHBAND

!---------------------------------------------------------------------------*/
end program f90_SFM
!---------------------------------------------------------------------------*/