module gvariables
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    real(8),parameter::g=9.80665D0     !Gravity acceleration
    integer,parameter::N=10            !Number of study cases for diameter of shaft

    character(len=256)::Title   !Title                       (input data)
    real(8)::Hg              ! Total head                    (input data)
    real(8)::Q0              ! Max.discharge                 (input data)
    real(8)::L               ! Length of pressure tunnel     (input data)
    real(8)::c               ! Head loss coefficient         (input data)
    real(8)::Cd              ! Discharge coefficient of port (input data)
    real(8)::zm              ! Max.water level               (input data) 
    real(8)::d0              ! Diameter of pressure tunnel   (input data)
    real(8)::xc              ! Chosen diameter of port
    real(8)::yc              ! Chosen diameter of shaft
    real(8)::zc              ! Max.water level of shaft with chosen diameter
    real(8)::R_DFC2          ! Min.diameter of shaft defined by dynamic stability
    real(8)::DF(1:N)         ! Diameter of shaft (array for case study)
    real(8)::zmax(1:N)       ! Max.water level (array for case study)

    real(4)::draw_xmin,draw_xmax,draw_dx
    real(4)::draw_ymin,draw_ymax,draw_dy

end module gvariables


program f90_surgeSCVF
    !*********************************************
    ! Stability and Max.Water level of surge tank
    !*********************************************
    use gvariables
    implicit none
    character(len=100)::fnameR     !Input data file
    character(len=100)::fnameF_SC  !Drawing (eps) file name for SC
    character(len=100)::fnameF_VF  !drawing (eps) file name for VF

    call getarg(1,fnameR)
    call getarg(2,fnameF_SC)
    call getarg(3,fnameF_VF)

    open(11,file=fnameR,status='old')
    read(11,'(a)') Title
    read(11,'()')
    read(11,*) Hg,Q0,L,d0,c,Cd,zm,xc,yc
    close(11)

    call CAL_SC()
    call PLOT_SC(fnameF_SC)
    call CAL_VF()
    call PLOT_VF(fnameF_VF)

    write(6,*) 'End of calculation.'

stop


contains

    subroutine CAL_SC()
        integer::iflag
        real(8)::fwa,v0
        real(8)::h0                 ! Head loss of pressure tunnel before interception (h0=c*v0*v0)
        real(8)::k0                 ! Resistance of port at interception
        real(8)::DFC1,DFC2,DFC3     ! Diameter of shaft (critical value)
        real(8)::Fsa                ! Section area of shaft
        real(8)::Fpa                ! Section area of port
        real(8)::d2,d2min,d2max,d2d ! Study range of diameter of port

        ! Stability condition
        fwa=0.25D0*pi*d0*d0
        v0=Q0/fwa
        h0=c*v0*v0

        Fsa=L*fwa/c/2.0D0/g/(Hg-zm)
        DFC2=sqrt(4.0D0*Fsa/pi)
        R_DFC2=DFC2

        d2min=0.1D0
        d2max=dble(int(DFC2+1.0))
        d2d=0.1D0

        open(12, file='temp_SC.csv', status='replace')
        write(12,'(a)') '#d2,k0-h0,DFC1,DFC2,DFC3'
        d2=d2min-d2d
        iflag=-1
        do while(d2<d2max)
            d2=d2+d2d
            Fpa=0.25D0*pi*d2*d2
            k0=0.5D0/g*Q0*Q0/Cd/Cd/Fpa/Fpa
            if(k0-h0<0.0D0)exit
 
            Fsa=h0*L*fwa/c/(h0+k0)/g/Hg
            DFC1=sqrt(4.0D0*Fsa/pi)
 
            Fsa=h0*L*fwa*fwa*fwa/2.0D0/g/c/c/k0/Q0/Q0
            DFC3=sqrt(4.0D0*Fsa/pi)
 
            write(12,*) d2,',',k0-h0,',',DFC1,',',DFC2,',',DFC3
        end do
        close(12)

        draw_xmin=0.0
        draw_xmax=real(int(d2-d2d))
        draw_dx=0.5
        draw_ymin=0.0
        draw_ymax=real(int((DFC2+1.0D0)/2.0D0))*6.0
        draw_dy=2.0

    end subroutine CAL_SC


    subroutine CAL_VF()
        integer::i,iflag,io
        real(8)::fwa,v0,md
        real(8)::Fsa                ! 水槽断面積 
        real(8)::Fpa                ! 制水口内径,断面積 
        real(8)::d1,d1min,d1d       ! 水槽内径計算範囲指定 
        real(8)::d2,d2min,d2max,d2d ! 制水口内径計算範囲指定 
        real(8)::h0                 ! 遮断前圧力水路損失水頭 h0=c*v0*v0 
        real(8)::k0                 ! 遮断時制水口抵抗 
        character(len=1000)::linebuf
        character(len=50)::str
        real(4)::work

        ! Max.water level by Vogt-Forchheimer formula
        fwa=0.25D0*pi*d0*d0
        v0=Q0/fwa
        h0=c*v0*v0

        d1min=dble(int(R_DFC2+1.0D0))
        d1d=1.0D0
        do i=1,N
            DF(i)=d1min+d1d*dble(i-1)
        end do

        open(12, file='temp_VF.csv', status='replace')
        linebuf='#d2,k0-h0'
        do i=1,N
            write (str,'(f5.1)') DF(i)
            linebuf=trim(linebuf)//','//str
        end do
        write(12,'(a)') trim(linebuf)
        d2min=0.1D0
        d2d=0.1D0
        d2=d2min-d2d
        do while(d2<draw_xmax)
            d2=d2+d2d
            if(draw_xmax<d2)exit
            Fpa=0.25D0*pi*d2*d2
            k0=0.5D0/g*Q0*Q0/Cd/Cd/Fpa/Fpa
            if(k0-h0<0.0D0)exit
            iflag=1 ! flag for log function
            do i=1,N
                d1=DF(i)
                Fsa=0.25D0*pi*d1*d1
                md=2.0D0*g*Fsa*(h0+k0)/L/fwa/v0/v0
                zmax(i)=CAL_ZM(k0,h0,md,iflag)
            end do
            if(iflag==1)then
                write (str,'(f8.3)') d2
                linebuf=trim(str)
                write (str,'(f10.3)') k0-h0
                linebuf=trim(linebuf)//','//trim(str)
                do i=1,N
                    write (str,'(f8.3)') zmax(i)
                    linebuf=trim(linebuf)//','//trim(str)
                end do
                write(12,'(a)') trim(linebuf)
            end if
        end do
        close(12)

        iflag=1
        d1=yc  !Selected diameter of shaft
        d2=xc  !Selected diameter of port
        Fsa=0.25D0*pi*d1*d1
        Fpa=0.25D0*pi*d2*d2
        k0=0.5D0/g*Q0*Q0/Cd/Cd/Fpa/Fpa
        md=2.0D0*g*Fsa*(h0+k0)/L/fwa/v0/v0
        zc=CAL_ZM(k0,h0,md,iflag) !Max.water level

        open(11,file='temp_VF.csv',status='old')
        read(11,'()')
        do 
            read(11,*,iostat=io) work,work,(zmax(i),i=1,N)
            if(io<0)exit
        end do
        close(11)

        draw_ymin=0.0
        draw_ymax=real(int(0.1*zmax(1))+1.0)*10.0
        draw_dy=10.0
        draw_ymax=100.0
    end subroutine CAL_VF


    real(8) function CAL_ZM(k0,h0,md,iflag)
        real(8),intent(in)::k0,h0,md
        integer,intent(inout)::iflag

        real(8)::ff,df,z1,z2,z
        real(8)::eps=1.0D-10

        ff=eps*10.0D0
        if(md*k0<1.0D0)then
            z1=-1.0D0/md+0.0001D0
            do while(eps<abs(ff))
                z=z1
                if(1.0D0+md*z<0.0D0)then
                    iflag=0
                    exit
                end if
                ff=((1.0D0+md*z)-log(1.0D0+md*z))-((1.0D0+md*h0)-log(1.0D0-md*k0))
                df=md*(1.0D0-1.0D0/(1.0D0+md*z))
                z2=z1-ff/df
                z1=z2
            end do
        end if
        if(1.0D0<md*k0)then
            z1=1.0D0/md+0.0001D0
            do while(eps<abs(ff))
                z=z1
                if(md*z-1.0D0<0.0D0)then
                    iflag=0
                    exit
                end if
                ff=((md*z-1.0D0)+log(md*z-1.0D0))-(log(md*k0-1.0D0)-(md*h0+1.0D0))
                df=md*(1.0D0+1.0D0/(md*z-1.0D0))
                z2=z1-ff/df
                z1=z2
            end do
        end if
        if(iflag==0)z1=0.0D0

        CAL_ZM=abs(z1)
    end  function CAL_ZM


    subroutine PLOT_SC(fnameF_SC)
        !******************************
        !Make script file for gnuplot (SC)
        !******************************
        character(len=100),intent(in)::fnameF_SC  !Drawing (eps) file name for SC
        character(len=1000)::linebuf
        character(len=100)::dummys
        character(len=50)::dummy1,dummy2,sloc
        real(4)::xx,yy
        integer::k

        open(12, file='gpl_temp_SC.txt', status='replace')
        write(12,'(a)') 'set terminal postscript eps enhanced font "Helvetica" 16'
        linebuf='set output "'//trim(adjustl(fnameF_SC))//'"'
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set datafile separator ","'
        write(12,'(a)') 'set size ratio 1.0'
        if(0<index(Title,'Load'))write(12,'(a)') 'set key left top'
        if(0<index(Title,'Input'))write(12,'(a)') 'set key right top'
        write(12,'(a)') 'set key width -2'
        write(12,'(a)') 'set key spacing 1.2'
        write(12,'(a)') 'set key box'
        write(12,'(a)') 'set xlabel "Diameter of the port (m)"'
        write(12,'(a)') 'set ylabel "Diameter of the shaft (m)"'

        !Range of axes
        write (dummy1,'(f8.3)') draw_xmin
        write (dummy2,'(f8.3)') draw_xmax
        linebuf='set xrange ['//trim(adjustl(dummy1))//':'//trim(adjustl(dummy2))//']'
        write(12,'(a)') trim(adjustl(linebuf))
        write (dummy1,'(f8.3)') draw_ymin
        write (dummy2,'(f8.3)') draw_ymax
        linebuf='set yrange ['//trim(adjustl(dummy1))//':'//trim(adjustl(dummy2))//']'
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set tmargin 2'
        !x-axis
        write (dummy1,'(f8.3)') draw_dx
        linebuf='set xtics '//trim(adjustl(dummy1))
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set format x "%.1f"'
        !y-axis
        write (dummy1,'(f8.3)') draw_dy
        linebuf='set ytics '//trim(adjustl(dummy1))
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set format y "%.0f"'

        write(12,'(a)') 'set tmargin 2'
        !Title
        k=1
        dummys=Title
        xx=0.0
        yy=1.03
        sloc='left'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set rmargin 8'

        !Show selected point
        write(12,'(a)') 'set style arrow 2 size graph 0.05,20 filled nohead linewidth 1 linetype 1'
        call MAKE_ARROW(2,draw_xmin,real(yc),real(xc),real(yc),linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
        call MAKE_ARROW(2,real(xc),draw_ymin,real(xc),real(yc),linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
        k=k+1
        write (dummy1,'(f8.1)') xc
        dummys=trim(adjustl(dummy1))//'m'
        xx=(xc-draw_xmin)/(draw_xmax-draw_xmin)+0.02
        yy=0.02
        sloc='left rotate by 90'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(linebuf)
        k=k+1
        write (dummy1,'(f8.1)') yc
        dummys=trim(adjustl(dummy1))//'m'
        xx=0.02
        yy=(yc-draw_ymin)/(draw_ymax-draw_ymin)+0.02
        sloc='left'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(linebuf)

       !Required min.diameter of shaft
        k=k+1
        write (dummy1,'(f8.3)') R_DFC2
        write (dummy2,'(f8.1)') zm
        dummys=trim(adjustl(dummy1))//'m (Target Z{/=10 m}='//trim(adjustl(dummy2))//'m)'
        xx=0.02
        yy=(R_DFC2-draw_ymin)/(draw_ymax-draw_ymin)+0.02
        sloc='left'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(linebuf)

        !Set axes
        write(12,'(a)') 'set grid xtics ytics mxtics mytics'
        !Plot command
        write(12,'(a)') 'plot \'
        write(12,'(a)') '"temp_SC.csv" using 1:3 with lines linewidth 3 linetype 4 title "D{/=10 req} > D{/=10 1}", \'
        write(12,'(a)') '"temp_SC.csv" using 1:4 with lines linewidth 3 linetype 1 title "D{/=10 req} > D{/=10 2}", \'
        write(12,'(a)') '"temp_SC.csv" using 1:5 with lines linewidth 3 linetype 2 title "D{/=10 req} < D{/=10 3}", \'
        write(12,'(a)') '"-" with points pointsize 2 pointtype 7 title "Selected point"'
        write (dummy1,'(f8.3)') xc
        write (dummy2,'(f8.3)') yc
        linebuf=dummy1//','//dummy2
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'e'
        close(12)

        !**************************
        !Execute gnuplot
        !**************************
        call system('gnuplot "gpl_temp_SC.txt"')

    end subroutine PLOT_SC


    subroutine PLOT_VF(fnameF_VF)
        !******************************
        !Make script file for gnuplot (VF)
        !******************************
        character(len=100),intent(in)::fnameF_VF  !drawing (eps) file name for VF
        character(len=1000)::linebuf
        character(len=100)::dummys
        character(len=50)::dummy1,dummy2,sloc
        real(4)::xx,yy
        integer::i,k

        open(12, file='gpl_temp_VF.txt', status='replace')
        write(12,'(a)') 'set terminal postscript eps enhanced font "Helvetica" 16'
        linebuf='set output "'//trim(adjustl(fnameF_VF))//'"'
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set datafile separator ","'
        write(12,'(a)') 'set size ratio 1.0'
        if(0<index(Title,'Load'))write(12,'(a)') 'set key left top'
        if(0<index(Title,'Input'))write(12,'(a)') 'set key right top'
        write(12,'(a)') 'set key width -4'
        write(12,'(a)') 'set key spacing 1.2'
        write(12,'(a)') 'set key box'
        write(12,'(a)') 'set xlabel "Diameter of the port (m)"'
        write(12,'(a)') 'set ylabel "Maximum WL  Z{/=10 max} (m)  or  Difference of pressure  |Z^{/Symbol \054}{/=10 rm}| (m)"'

        !Range of axes
        write (dummy1,'(f8.3)') draw_xmin
        write (dummy2,'(f8.3)') draw_xmax
        linebuf='set xrange ['//trim(adjustl(dummy1))//':'//trim(adjustl(dummy2))//']'
        write(12,'(a)') trim(adjustl(linebuf))
        write (dummy1,'(f8.3)') draw_ymin
        write (dummy2,'(f8.3)') draw_ymax
        linebuf='set yrange ['//trim(adjustl(dummy1))//':'//trim(adjustl(dummy2))//']'
        write(12,'(a)') trim(adjustl(linebuf))
        !x-axis
        write (dummy1,'(f8.3)') draw_dx
        linebuf='set xtics '//trim(adjustl(dummy1))
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set format x "%.1f"'
        !y-axis
        write (dummy1,'(f8.3)') draw_dy
        linebuf='set ytics '//trim(adjustl(dummy1))
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set format y "%.0f"'

        write(12,'(a)') 'set tmargin 2'
        !Title
        k=1
        dummys=Title
        xx=0.0
        yy=1.03
        sloc='left'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'set rmargin 8'

        !Required Max.water level
        write(12,'(a)') 'set style arrow 2 size graph 0.05,20 filled nohead linewidth 1 linetype 1'
        call MAKE_ARROW(2,draw_xmin,real(zm),draw_xmax,real(zm),linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
        k=k+1
        write (dummy1,'(f8.3)') zm
        dummys='Target Z{/=10 m}='//trim(adjustl(dummy1))//'m'
        xx=0.02
        yy=(zm-draw_ymin)/(draw_ymax-draw_ymin)+0.02
        sloc='left'
        call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
        write(12,'(a)') trim(linebuf)

        !Shaft diameter
        do i=1,N
            k=k+1
            write (dummy1,'(i2)') int(DF(i))
            dummys='D='//trim(adjustl(dummy1))//'m'
            xx=1.01
            yy=(real(zmax(i))-draw_ymin)/(draw_ymax-draw_ymin)
            sloc='left'
            call MAKE_LABEL(k,dummys,xx,yy,sloc,linebuf)
            write(12,'(a)') trim(linebuf)
        end do

        !Set axses
        write(12,'(a)') 'set grid xtics ytics mxtics mytics'
        !Plot command
        write(12,'(a)') 'plot \'
        do i=1,N
            write (dummy1,'(i2)') 2+i
            linebuf='"temp_VF.csv" using 1:'//trim(adjustl(dummy1))//' with line linewidth 3 linetype 1 notitle, \'
            write(12,'(a)') trim(adjustl(linebuf))
        end do
        write(12,'(a)') '"temp_VF.csv" using 1:2 with linespoints linewidth 3 pointtype 6 &
                         title "|Z^{/Symbol \054}{/=10 rm}| = k{/=10 0} - h{/=10 0}", \'
        write(12,'(a)') '"-" with points pointsize 2 pointtype 7 title "Selected point"'
        write (dummy1,'(f8.3)') xc
        write (dummy2,'(f8.3)') zc
        linebuf=dummy1//','//dummy2
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'e'
        close(12)

        !**************************
        !Execute gnuplot
        !**************************
        call system('gnuplot "gpl_temp_VF.txt"')

    end subroutine PLOT_VF


    subroutine MAKE_LABEL(i,dummys,xx,yy,sloc,linebuf)
        integer,intent(in)::i
        character(len=100),intent(in)::dummys
        real(4),intent(in)::xx,yy
        character(len=50),intent(in)::sloc
        character(len=1000),intent(out)::linebuf

        character(len=20)::dummyi,dummyx,dummyy

        write (dummyi,'(i2)') i
        write (dummyx,'(f8.3)') xx
        write (dummyy,'(f8.3)') yy
        linebuf='set label '//trim(adjustl(dummyi))//' "'//trim(adjustl(dummys))//'" at graph '//&
                 trim(adjustl(dummyx))//',graph '//trim(adjustl(dummyy))//' '//trim(adjustl(sloc))
    end subroutine MAKE_LABEL


    subroutine MAKE_ARROW(iarrow,x1,y1,x2,y2,linebuf)
        integer,intent(in)::iarrow
        real(4),intent(in)::x1,y1,x2,y2
        character(len=1000),intent(out)::linebuf

        character(len=20)::dummyi,dummyx1,dummyy1,dummyx2,dummyy2

        write (dummyi,'(i2)') iarrow
        write (dummyx1,'(f8.3)') x1
        write (dummyy1,'(f8.3)') y1
        write (dummyx2,'(f8.3)') x2
        write (dummyy2,'(f8.3)') y2
        linebuf='set arrow as '//trim(adjustl(dummyi))//' from '//trim(adjustl(dummyx1))//','//trim(adjustl(dummyy1))//&
                ' to '//trim(adjustl(dummyx2))//','//trim(adjustl(dummyy2))
    end subroutine MAKE_ARROW


end program f90_surgeSCVF