program f90_PQ4
    !Saturation line
    implicit none
    integer::i,j,k,n,ne,kk
    character::strcom*50            !書出用コメント
    integer::NODT                   !節点総数
    integer::NELT                   !要素総数
    integer::MATEL                  !材料種類数
    integer::KOH                    !全水頭指定節点数
    integer::KOQ                    !流量指定節点数
    integer::KOU                    !浸潤境界指定節点数
    integer::nt                     !全自由度（総節点数×１[１節点自由度]）
    integer::nod                    !1要素節点数
    integer::nhen=1                 !1節点自由度
    integer,allocatable::kakom(:,:) !要素構成節点番号
    integer,allocatable::matno(:)   !材料種別No
    real(4),allocatable::x(:)       !節点x座標
    real(4),allocatable::z(:)       !節点y座標
    real(4),allocatable::hvec(:)    !全体全水頭ベクトル
    real(4),allocatable::qvec(:)    !全体節点流量ベクトル
    real(4),allocatable::pvec(:)    !圧力指定点流量
    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


    !***************************
    !data input
    !***************************
    open(11,file='fnameW.csv',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)


    !堤体表面の節点の検索
    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
    !x座標の小さい順の並び替え
    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

    !流量プロットデータ作成
    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

    !作図範囲設定
    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.0e-4<Qmax)y2max=1.0e-3
    if(1.0e-4<Qmax.and.Qmax<=5.0e-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,*) 'set terminal postscript eps enhanced font "Helvetica" 16'
    write(12,*) 'set output "fig_pq.eps"'
    write(12,*) '#'
    write(12,*) 'set xrange [',xmin,':',xmax,']'
    write(12,*) 'set yrange [',ymin,':',ymax,']'
    write(12,*) 'set y2range [',y2min,':',y2max,']'
    write(12,*) 'set xtics ',dx
    write(12,*) 'set ytics ',dy
    write(12,*) 'set y2tics ',dy2
    write(12,*) 'set format x "%.0f"'
    write(12,*) 'set format y "%.0f"'
    write(12,*) 'set format y2 "%9.1E"'
    write(12,*) 'set grid xtics ytics mxtics mytics'
    write(12,*) 'set xlabel "Distance from upstream edge (m)"'
    write(12,*) 'set ylabel "Pressure head  H (m)"'
    write(12,*) 'set y2label "Seepage discharge  Q (m^3/s/m^2)"'
    write(12,*) 'set key box'
    write(12,*) 'set key left bottom'
    write(12,*) 'set xzeroaxis lt -1'
    write(12,*) '#'
    write (s,'(e15.7)') QTi
    sline='set label "QTi='//trim(s)//'" at graph 0.95,graph 0.95 right'
    write(12,*) trim(sline)
    write (s,'(e15.7)') QTo
    sline='set label "QTo='//trim(s)//'" at graph 0.95,graph 0.90 right'
    write(12,*) trim(sline)
    sline='set label "'//trim(adjustl(strcom))//'" at graph 0.02,graph 0.95 left'
    write(12,*) trim(sline)
    write(12,*) '#'
    write(12,*) 'plot \'
    write(12,*) '"-" axes x1y1 with lp pt 4 ps 1 lt 1 lw 1 title "H of Dam surface",\'
    write(12,*) '"-" axes x1y1 with lp pt 8 ps 1 lt 1 lw 1 title "H of Bottom line",\'
    write(12,*) '"-" axes x1y2 with steps lt 1 lw 3 title "Q of 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,*) 'e'
    !堤体基礎圧力水頭プロットデータ
    do i=1,nd
        if(zs(i)<ds)write(12,*) xs(i),' ',ps(i)
    end do
    write(12,*) 'e'
    !堤体表面流量プロットデータ
    do i=1,k
        write(12,*) zw(i),' ',qw(i)
    end do
    write(12,*) 'e'
    write(12,*) 'quit'
    close(12)

    !gnuplot 実行
    call system('gnuplot gpl_pq.txt')

stop

end program f90_PQ4