module gvariables
    implicit none

    real(8),parameter::pi=3.14159265358979323846D0
    real(8),parameter::g=9.80665D0     !Gravity acceleration
    integer,parameter::NUMS=10         !Limit number of section definition points of Surge Tank
    integer,parameter::NUMQ=20         !Limit number of discharge definition points

    character(len=256)::Title   !Title
    integer::ICT                !Calculation case (1: normal, 2: AFC)
    real(8)::AFCA               !Discharge of half amplitude of AFC(m3/s)
    real(8)::AFCT               !Period of AFC(s)

    !Definition of times
    real(8)::TMAX            !End time for calculation
    real(8)::dt              !Time increment for calculation
    real(8)::DTWR            !Time increment for print out

    !Definition of properties of orifice (port)
    real(8)::PAA             !Section area of port
    real(8)::PCI             !Discharge coefficient of port for in-flow
    real(8)::PCO             !Discharge coefficient of port for out-flow

    !Definition of propertoes of pressure tunnel & reservior
    real(8)::RWL             !Water level at reservoir (EL)
    real(8)::TNL             !Length of pressure tunnel
    real(8)::TNA             !Section area of presure tunnel
    real(8)::TNC             !Head loss coefficient of pressure tunnel

    !Definition of properties of shaft
    integer::NST                   !Number of sections
    real(8)::SEL(NUMS+2)           !Level definition of sections of shaft (EL)
    real(8)::SAA(NUMS+2)           !Section area of specified level in the shaft
    character(len=50)::SLB(NUMS+1) !Label for section (text: for drawing)

    !Difinition of discharge
    integer::NQT             !Number of discharge input points
    real(8)::QTQ(NUMQ+1)     !Discharge at specified time point
    real(8)::QTI(NUMQ+1)     !Time at specified time point

    !Others
    real(8)::AFCS            !Starting time of interception during AFC operation
    real(8)::QI              !Initial discharge

    character(len=50)::fnameR,fnameW !Input & output file name

end module gvariables



program f90_surge
    use gvariables
    implicit none
    integer::iret

    call getarg(1,fnameR)
    call getarg(2,fnameW)

    call DINP()       !Data input
    call PRE()        !Pre-treatment
    call CALC(iret)   !Execution of calculation
    call PLTDATA()    !DATA write for drawing

    print "('End of calculation')"
    if(iret==0)print "(i0,': Succsessfull termination')",iret
    if(iret==1)print "(i0,': Water Level is greater than Top of the shaft')",iret
    if(iret==2)print "(i0,': Water Level is less than Bottom of the shaft')",iret

    stop


contains

    subroutine DINP()
        integer::i

        open(11,file=fnameR,status='old')
        read(11,'(a)') Title
        read(11,*) ICT,AFCA,AFCT
        read(11,*) TMAX,dt,DTWR
        read(11,*) PAA,PCI,PCO
        read(11,*) RWL,TNL,TNA,TNC
        read(11,*) NST
        do i=1,NST
            read(11,*) SAA(i),SEL(i),SLB(i)
        end do
        read(11,*) NQT
        do i=1,NQT
            read(11,*) QTQ(i),QTI(i)
        end do
        close(11)
    end subroutine DINP


    subroutine PRE()
        integer::i,j
        real(8)::wa,wb
        character(len=50)::wc

        !QI : Initial discharge
        !AFCS : Starting time of interception (normal condition: AFCS=-1.0)
        !In normal condition, starting time of interception is set on QTI(1).
        !In AFC operation, starting time of interception is set on QTI(2).

        QI=QTQ(1)
        select case(ICT)
            case(1)
                AFCS=-1.0D0 !Normal
            case(2)
                AFCS=QTI(2) !AFC
        end select

        !To verify the order of section levels of shaft (from bottom to top)
        do i=1,NST
            do j=i+1,NST
                if(SEL(i)>=SEL(j))then  !To sort section numbers in small order of elevation
                    wa=SEL(j)
                    wb=SAA(j)
                    wc=SLB(j)
                    SEL(j)=SEL(i)
                    SAA(j)=SAA(i)
                    SLB(j)=SLB(i)
                    SEL(i)=wa
                    SAA(i)=wb
                    SLB(i)=wc
                end if
            end do
        end do
        SAA(NST+1)=SAA(NST)!Dummy for judgement of top level (section area)
        SEL(NST+1)=SEL(NST)!Dummy for judgement of top level (EL)
    end subroutine PRE


    subroutine CALC(iret)
        integer:: j,ic,np,npe,npn,npw,iw
        real(8) dto,qn,vi,zi,zo,tmi,tmx,zmi,zmx,dts,fa,dv,dz,zw,vw,dti,dtw
        integer,intent(out)::iret

        iret=0
        open(12, file=fnameW, status='replace')
        write (12,'(a)') '#T(s),WL(m),v(m/s),Q(m3/s),k(m)'

        np=int(DTWR/dt+0.00001D0) !Step number of print out
        dto=dt                    !Time increment
        qn=QI                     !Initial discharge
        vi=qn/TNA                 !Initial velocity of tunnel discharge
        zi=TNC*vi*vi              !Initial head loss due to initial velocity
        if(vi<0.0D0)zi=-TNC*vi*vi !Definition of sign of initial head loss (direction of flow)
        zo=RWL-zi                 !Initial water level in Surge Tank

        !To print initial values (time=0)
        dts=0.0D0
        zw=zo
        vw=vi
        zi=zi
        qn=qn
        call WOUT(dts,zw,vw,qn)

        !To search the section number with initial water level
        do j=NST,1,-1
            if(SEL(j)<zo)exit
        end do

        npe=j      !Number of section with initial water level
        npw=0      !Initialization of counter for print out
        dts=0.0D0  !Initialization of elapsed time

        do while(dts<=TMAX)
            dts=dts+dto  !The time found solution
            npw=npw+1
            fa=SAA(npe)
            call RUNGE(qn,vi,zi,fa,dto,dv,dz,dts)
            zw=RWL-zi-dz

            iw=npe+1
            if(iw>NST)iw=NST
            if(zw<SEL(npe).or.SEL(iw)<zw)then
                dti=0.1D0*dto
                do ic=1,10
                    dtw=dts-dto+dti*dble(ic)
                    zw=RWL-zi
                    if(SEL(npe+1)<zw)npe=npe+1
                    if(zw<SEL(npe))npe=npe-1
                    if(NST<npe)iret=1 !water level greater than USWL
                    if(npe<1)iret=2   !water level less than DSWL
                    fa=SAA(npe)
                    call RUNGE(qn,vi,zi,fa,dti,dv,dz,dtw)
                    vi=vi+dv
                    zi=zi+dz
                end do
            else
                vi=vi+dv
                zi=zi+dz
            end if
            !To print results
            if(npw>=np)then
                npw=0
                zw=RWL-zi
                vw=vi
                call WOUT(dts,zw,vw,qn)
            end if
        end do
        iret=iret
        close(12)
    end subroutine CALC


    subroutine RUNGE(qo,vo,zo,fa,dtt,dv,dz,dts)
        real(8),intent(inout)::qo
        real(8),intent(out)::dv,dz
        real(8),intent(in)::vo,zo,fa,dtt,dts
        real(8)::wrk,dz0,dv0,cf,dtw,v1,z1,dv1,dz1,dv2,dz2,dv3,dz3

        wrk=0.0D0
        dz0=dtt*FZ(vo,qo,fa)
        cf=PCI*PAA                !Water level rising (in-flow from port)
        if(dz0>0.0D0)cf=PCO*PAA   !Water level drop (out-flow from port)
        wrk=FK(vo,qo,cf)
        dv0=dtt*FV(zo,vo,wrk)

        v1=vo+0.5D0*dv0
        z1=zo+0.5D0*dz0
        dtw=dts-dtt*0.5D0
        qo=QNCAL(dtw)
        wrk=FK(v1,qo,cf)
        dv1=dtt*FV(z1,v1,wrk)
        dz1=dtt*FZ(v1,qo,fa)

        v1=vo+0.5D0*dv1
        z1=zo+0.5D0*dz1
        wrk=FK(v1,qo,cf)
        dv2=dtt*FV(z1,v1,wrk)
        dz2=dtt*FZ(v1,qo,fa)

        v1=vo+dv2
        z1=zo+dz2
        dtw=dts
        qo=QNCAL(dtw)
        wrk=FK(v1,qo,cf)
        dv3=dtt*FV(z1,v1,wrk)
        dz3=dtt*FZ(v1,qo,fa)

        dv=(dv0+2.0D0*dv1+2.0D0*dv2+dv3)/6.0D0
        dz=(dz0+2.0D0*dz1+2.0D0*dz2+dz3)/6.0D0
    end subroutine RUNGE


    real(8) function FZ(v,q,fa)
        ! Calculation of dz 
        real(8),intent(in)::v,q,fa
        FZ=(q-TNA*v)/fa
    end function FZ


    real(8) function FK(v,q,cf)
        ! Calculation of k 
        real(8),intent(in)::v,q,cf
        FK=abs((TNA*v-q)/cf)*(TNA*v-q)/cf/2.0D0/g
    end function FK


    real(8) function FV(z,v,wrk)
        ! Calculation of dv 
        real(8),intent(in)::z,v,wrk
        FV=(z-TNC*abs(v)*v-wrk)/TNL*g
    end function FV


    real(8) function QNCAL(tt)
        ! Calculation of discharge by interpolation method 
        real(8),intent(in)::tt
        real(8)::qn
        integer::i

        if(tt>=AFCS)then
            do i=2,NQT
                if(QTI(i-1)<tt.and.tt<=QTI(i))then
                    qn=QTQ(i-1)+(QTQ(i)-QTQ(i-1))/(QTI(i)-QTI(i-1))*(tt-QTI(i-1))
                    exit
                end if
            end do
        else
            qn=QI+AFCA*sin(2.0D0*pi/AFCT*tt)
        end if
        QNCAL=qn
    end function QNCAL


    subroutine WOUT(dts,zw,vw,qn)
        real(8),intent(in)::dts,zw,vw,qn
        real(8)::k,cf
        character(len=1000)::linebuf
        character(len=200)::fmt1

        fmt1="(e15.7,',',e15.7,',',e15.7,',',e15.7,',',e15.7)"

        cf=0.5D0*(PCI+PCO)*PAA
        k=FK(vw,qn,cf)

        write(linebuf,FMT=fmt1) dts,zw,vw,qn,k
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end subroutine WOUT


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


    subroutine PLTDATA()
        integer,parameter::nn=10000
        integer::i,io
        real(8)::ti,ti_max,ti_min
        real(8)::wl,wl_max,wl_min

        open(11, file=fnameW, status='old')
        read(11,'()')
        wl_max=-1.0D5
        wl_min=1.0D5
        do i=1,nn
            read(11,*,iostat=io) ti,wl
            if(wl_max<wl)then
                wl_max=wl;ti_max=ti
            end if
            if(wl<wl_min)then
                wl_min=wl;ti_min=ti
            end if
            if(io<0)exit
        end do
        close(11)

        open(12,file='_draw.txt',status='replace')
        write(12,'(a)') Title
        write(12,*) 0,',',int(TMAX),',',60
        write(12,*) int((SEL(1)-10.0D0)/10.0D0)*10,',',int((SEL(NST)+20.0D0)/10.0D0)*10,',',10
        write(12,*) ti_max,',',wl_max
        write(12,*) ti_min,',',wl_min
        write(12,*) RWL,',','RWL'
        do i=1,NST
            write(12,*) SEL(i),',',SLB(i)
        end do
        close(12)
    end subroutine PLTDATA

end program f90_surge
