program f90_pq
    !Saturation line
    implicit none
    integer::i,j,k,n,ne,kk
    character::strcom*50            !Comment
    integer::NODT                   !Number of nodes
    integer::NELT                   !Number of elements
    integer::MATEL                  !Number of material sets
    integer::KOH                    !Number of given total head nodes
    integer::KOQ                    !Number of given discharge nodes
    integer::KOU                    !Number of given seepage point nodes
    integer::nt                     !Total degrees of freedom
    integer::nod                    !Number of nodes per element
    integer::nhen=1                 !Degree of freedom pwe node
    integer,allocatable::kakom(:,:) !Element-Nodes relationship
    integer,allocatable::matno(:)   !kind of material set
    real(4),allocatable::x(:)       !x-coordinate of node
    real(4),allocatable::z(:)       !y-coordinate of node
    real(4),allocatable::hvec(:)    !Nodal head vector
    real(4),allocatable::qvec(:)    !Nodal discharge vector
    real(4),allocatable::pvec(:)    !Discharge at given head nodes
    real(4)::QTi,QTo,Qmax
    real(4)::rw,el
    real(4)::xmin,xmax,dx,ymin,ymax,dy,y2min,y2max,dy2
    real(4),allocatable::xs(:)
    real(4),allocatable::zs(:)
    real(4),allocatable::qs(:)
    real(4),allocatable::ps(:)
    real(4),allocatable::xw(:)
    real(4),allocatable::zw(:)
    real(4),allocatable::qw(:)
    integer::nd
    integer::id,iw,ik
    real(4)::x1u,z1u,x2u,z2u,x1d,z1d,x2d,z2d
    real(4)::a1,b1,a2,b2,xx,zz
    real(4)::ds=0.01
    character::s*20,sline*100
    character::fnameR*50,fnameF*50,strx1*20,stry1*20,strx2*20,stry2*20

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

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,*) strcom,x1u,z1u,x2u,z2u,x1d,z1d,x2d,z2d
    read(11,'()')
    read(11,*) nod,NODT,NELT,MATEL,KOH,KOQ,KOU
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(x(1:NODT))
    allocate(z(1:NODT))
    allocate(hvec(1:nt))
    allocate(qvec(1:nt))
    allocate(pvec(1:nt))
    allocate(xs(1:NODT))
    allocate(zs(1:NODT))
    allocate(qs(1:NODT))
    allocate(ps(1:NODT))
    allocate(xw(1:NODT))
    allocate(zw(1:NODT))
    allocate(qw(1:NODT))

    read(11,'()')
    read(11,'()')
    !node,x,z,KOH,KOQ,KOU,Hinp,Qinp
    k=0
    do i=1,NODT
        read(11,*) iw,rw,rw,ik,iw,iw,rw,rw
    end do
    read(11,'()')
    read(11,'()')
    !element,node-1,node-2,node-3,node-4,Ak0,alpha,em,matno
    do ne=1,NELT
        read(11,*) i,(kakom(ne,j),j=1,nod),rw,rw,rw,matno(ne)
    end do
    read(11,'()')
    read(11,'()')
    !node,coord-x,coord-z,hvec,qvec,pvec
    do i=1,NODT
        read(11,*) j,x(i),z(i),hvec(i),qvec(i),pvec(i)
    end do
    close(11)


    !Search nodes on the dam surface
    a1=(z2u-z1u)/(x2u-x1u)
    b1=(x2u*z1u-x1u*z2u)/(x2u-x1u)

    a2=(z2d-z1d)/(x2d-x1d)
    b2=(x2d*z1d-x1d*z2d)/(x2d-x1d)

    k=0
    do i=1,NODT
        xx=x(i)
        zz=z(i)
        id=0
        if(a1*xx+b1-ds<zz)id=1
        if(a2*xx+b2-ds<zz)id=1
        if((x2u+ds<xx.and.xx<x1d-ds).and.z2u-ds<zz)id=1
        if((x1u+ds<xx.and.xx<x2d-ds).and.zz<ds)id=1
        if(id==1)then
            k=k+1
            xs(k)=x(i)
            zs(k)=z(i)
            qs(k)=qvec(i)
            ps(k)=pvec(i)
        end if
    end do
    nd=k
    !Sort od nodes
    do i=1,nd-1
        do j=i+1,nd
            if(xs(i)>xs(j))then
                rw=xs(i)
                xs(i)=xs(j)
                xs(j)=rw

                rw=zs(i)
                zs(i)=zs(j)
                zs(j)=rw

                rw=qs(i)
                qs(i)=qs(j)
                qs(j)=rw

                rw=ps(i)
                ps(i)=ps(j)
                ps(j)=rw
            end if
        end do
    end do

    QTi=0.0
    QTo=0.0
    do i=1,nd
        if(0.0<qs(i))QTi=QTi+qs(i)
        if(qs(i)<0.0)QTo=QTo+qs(i)
    end do

    !Plot data of Discharge
    k=0
    k=k+1
    xw(k)=xs(1)
    qw(k)=qs(1)
    do i=1,nd
        if(ds<zs(i))then
            k=k+1
            xw(k)=xs(i)
            qw(k)=qs(i)
        end if
    end do
    k=k+1
    xw(k)=xs(nd)
    qw(k)=qs(nd)
    el=0.5*(xw(2)-xw(1))
    zw(1)=xw(1)
    qw(1)=qw(1)/el
    do i=2,k-1
        el=0.5*(xw(i)-xw(i-1))+0.5*(xw(i+1)-xw(i))
        zw(i)=0.5*(xw(i-1)+xw(i))
        qw(i)=qw(i)/el
    end do
    el=0.5*(xw(k)-xw(k-1))
    zw(k)=xw(k)
    qw(k)=qw(k)/el
    Qmax=abs(qw(1))
    do i=1,k
        if(Qmax<abs(qw(i)))Qmax=abs(qw(i))
    end do

    !Drawing range
    xmin=0.0
    xmax=x2d
    if(xmax<=10.0)dx=1.0
    if(10.0<xmax.and.xmax<=30.0)dx=2.0
    if(30.0<xmax.and.xmax<=60.0)dx=5.0
    if(60.0<xmax)dx=10.0

    ymax=z2u+5.0
    ymin=-ymax
    dy=5.0

    if(5.5e-4<Qmax)y2max=1.0e-3
    if(1.0e-4<Qmax.and.Qmax<=5.5e-4)y2max=5.0e-4
    if(5.0e-5<Qmax.and.Qmax<=1.0e-4)y2max=1.0e-4
    if(1.0e-5<Qmax.and.Qmax<=5.0e-5)y2max=5.0e-5
    if(Qmax<=1.0e-5)y2max=1.0e-5
    y2min=-y2max
    dy2=y2max/5.0

    open(12, file='_gpl_pq.txt', status='replace')
    write(12,'(a)') 'set terminal postscript eps enhanced size 15cm,10cm font "Helvetica" 16'
    write(12,'(a)') 'set output "'//trim(adjustl(fnameF))//'"'
    write(12,'(a)') '#'
    write(strx1,*) xmin
    write(strx2,*) xmax
    write(stry1,*) ymin
    write(stry2,*) ymax
    write(12,'(a)') 'set xrange ['//trim(adjustl(strx1))//':'//trim(adjustl(strx2))//']'
    write(12,'(a)') 'set yrange ['//trim(adjustl(stry1))//':'//trim(adjustl(stry2))//']'
    write(stry1,*) y2min
    write(stry2,*) y2max
    write(12,'(a)') 'set y2range ['//trim(adjustl(stry1))//':'//trim(adjustl(stry2))//']'
    write(strx1,*) dx
    write(stry1,*) dy
    write(stry2,*) dy2
    write(12,'(a)') 'set xtics '//trim(adjustl(strx1))
    write(12,'(a)') 'set ytics '//trim(adjustl(stry1))
    write(12,'(a)') 'set y2tics '//trim(adjustl(stry2))
    write(12,'(a)') 'set format x "%.0f"'
    write(12,'(a)') 'set format y "%.0f"'
    write(12,'(a)') 'set format y2 "%9.1E"'
    write(12,'(a)') 'set grid xtics ytics mxtics mytics'
    write(12,'(a)') 'set xlabel "Horizontal distance from upstream edge (m)"'
    write(12,'(a)') 'set ylabel "Pressure head  H (m)"'
    write(12,'(a)') 'set y2label "Seepage discharge  Q (m^3/s/m^2)"'
    write(12,'(a)') 'set key box'
    write(12,'(a)') 'set key left bottom'
    write(12,'(a)') 'set xzeroaxis lt -1'
    write(12,'(a)') '#'
    write (s,'(e15.7)') QTi
    sline='set label "QTi='//trim(s)//'" at graph 0.95,graph 0.95 right'
    write(12,'(a)') trim(sline)
    write (s,'(e15.7)') QTo
    sline='set label "QTo='//trim(s)//'" at graph 0.95,graph 0.90 right'
    write(12,'(a)') trim(sline)
    sline='set label "'//trim(adjustl(strcom))//'" at graph 0.02,graph 0.95 left'
    write(12,'(a)') trim(sline)
    write(12,'(a)') '#'
    write(12,'(a)') 'plot \'
    write(12,'(a)') '"-" axes x1y1 with lp pt 4 ps 1 lt 1 lw 1 title "H of Dam surface",\'
    write(12,'(a)') '"-" axes x1y1 with lp pt 8 ps 1 lt 1 lw 1 title "H of Bottom line",\'
    write(12,'(a)') '"-" axes x1y2 with steps lt 1 lw 3 title "Q of Dam surface"'
    !Pressure head on the dam surface
    write(12,*) xs(1),' ',ps(1)
    do i=1,nd
        if(ds<zs(i))write(12,*) xs(i),' ',ps(i)
    end do
    write(12,*) xs(nd),' ',ps(nd)
    write(12,'(a)') 'e'
    !Pressure head on the dam foundation
    do i=1,nd
        if(zs(i)<ds)write(12,*) xs(i),' ',ps(i)
    end do
    write(12,'(a)') 'e'
    !Discharge from dam surface
    do i=1,k
        write(12,*) zw(i),' ',qw(i)
    end do
    write(12,'(a)') 'e'
    write(12,'(a)') 'quit'
    close(12)

    !gnuplot
    call system('gnuplot _gpl_pq.txt')

stop

end program f90_pq
