program f90_GMTseep
    !Saturation line
    implicit none
    integer::i,j,k,l,n,ne,kk,ii,jj
    integer::n1,n2,n3,k1,k2,k3
    character::strcom*256           !書出用コメント
    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(:)    !圧力指定点流量
    integer,allocatable::nokh(:)
    integer,allocatable::nokq(:)
    integer,allocatable::noku(:)
    real(4),allocatable::vx(:)      !x方向流速ベクトル
    real(4),allocatable::vz(:)      !y方向流速ベクトル
    real(4),allocatable::Akr(:)     !透水係数比
    real(4)::QTi,QTo
    real(4)::rw
    real(4)::x1,y1,z1,x2,y2,z2,x3,y3,z3
    real(4)::xg,zg,xx,zz,yy
    real(4),allocatable::hpv(:)
    real(4),allocatable::satx(:,:)
    real(4),allocatable::satz(:,:)
    integer,allocatable::nd(:)
    real(4),allocatable::yp(:)
    integer::id,iw,ik
    integer::kkmax,kkm1,kkm2,nsat,kdraw

    character :: linebuf*1000
    character :: filename*128

    real(4)::bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4


    !***************************
    !data input
    !***************************
    open(11,file='fnameW.csv',status='old')
    read(11,*) strcom,bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4
    read(11,'()')
    read(11,*) nod,NODT,NELT,MATEL,KOH,KOQ,KOU
    !------------------------------------------------------------------------------------------
    !配列寸法宣言
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    if(nod==3)nsat=NELT
    if(nod==4)nsat=NELT*4
    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(hpv(1:nt))
    allocate(vx(1:NELT))
    allocate(vz(1:NELT))
    allocate(Akr(1:NELT))
    allocate(nd(1:kkmax))
    allocate(satx(1:kkmax,1:nsat))
    allocate(satz(1:kkmax,1:nsat))
    if(0<KOH)allocate(nokh(1:KOH))
    if(0<KOQ)allocate(nokq(1:KOQ))
    if(0<KOU)allocate(noku(1:KOU))

    read(11,'()')
    read(11,'()')
    !node,x,z,KOH,KOQ,KOU,Hinp,Qinp
    n1=0
    n2=0
    n3=0
    do i=1,NODT
        read(11,*) j,rw,rw,k1,k2,k3,rw,rw
        if(0<KOH)then
            if(k1==1)then
                n1=n1+1
                nokh(n1)=i
            end if
        end if
        if(0<KOQ)then
            if(k2==1)then
                n2=n2+1
                nokq(n2)=i
            end if
        end if
        if(0<KOU)then
            if(k3==1)then
                n3=n3+1
                noku(n3)=i
            end if
        end if
    end do
    read(11,'()')
    read(11,'()')
    !element,node-1,node-2,node-3,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
    read(11,'()')
    read(11,'()')
    !element,vx,vz,vm,Ak0,Akr,matno
    do ne=1,NELT
        read(11,*) j,vx(ne),vz(ne),rw,rw,Akr(ne),iw
    end do
    close(11)


    call DRAW_MESH(kakom,matno,x,z,nokh,nokq,noku,nod,NODT,NELT,KOH,KOQ,KOU)
    call DRAW_CONT(kakom,x,z,bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4,hvec,pvec,nod,NODT,NELT)
    call DRAW_VECT(kakom,x,z,bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4,pvec,vx,vz,nod,NODT,NELT)


stop


contains

    subroutine DRAW_MESH(kakom,matno,x,z,nokh,nokq,noku,nod,NODT,NELT,KOH,KOQ,KOU)
        integer::nod,NODT,NELT,KOH,KOQ,KOU
        integer::kakom(1:NELT,1:nod),matno(1:NELT)
        integer::nokh(1:KOH),nokq(1:KOQ),noku(1:KOU)
        real(4)::x(1:NODT),z(1:NODT)

        integer::i,j,k,l,ma,ne
        real(4)::xmin,xmax,zmin,zmax
        character::srange*100,sscale*100,sbmap*200,spsxy*200,filename*50,scol*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::baseheight=4.0

        xmin=(x(1))
        xmax=(x(1))
        zmin=(z(1))
        zmax=(z(1))
        do i=1,NODT
            if((x(i))<xmin)xmin=(x(i))
            if(xmax<(x(i)))xmax=(x(i))
            if((z(i))<zmin)zmin=(z(i))
            if(zmax<(z(i)))zmax=(z(i))
        end do

        do ma=1,MATEL
            write (filename, '("_mesh", i2.2, ".txt")') ma
            open (12, file=filename, status='replace')
            select case(nod)
                case(3)
                    do ne=1,NELT
                        if(matno(ne)==ma)then
                            i=kakom(ne,1)
                            j=kakom(ne,2)
                            k=kakom(ne,3)
                            write(12,*) '>'
                            write(12,*) (x(i)),' ',(z(i))
                            write(12,*) (x(j)),' ',(z(j))
                            write(12,*) (x(k)),' ',(z(k))
                        end if
                    end do
                case(4)
                    do ne=1,NELT
                        if(matno(ne)==ma)then
                            i=kakom(ne,1)
                            j=kakom(ne,2)
                            k=kakom(ne,3)
                            l=kakom(ne,4)
                            write(12,*) '>'
                            write(12,*) (x(i)),' ',(z(i))
                            write(12,*) (x(j)),' ',(z(j))
                            write(12,*) (x(k)),' ',(z(k))
                            write(12,*) (x(l)),' ',(z(l))
                        end if
                    end do
            end select
            close (12)
        end do

        if(0<KOH)then
            open (12, file='_koh.txt', status='replace')
            do i=1,KOH
                write(12,*) x(nokh(i)),' ',z(nokh(i))
            end do
            close(12)
        end if
        if(0<KOQ)then
            open (12, file='_koq.txt', status='replace')
            do i=1,KOQ
                write(12,*) x(nokq(i)),' ',z(nokq(i))
            end do
            close(12)
        end if
        if(0<KOU)then
            open (12, file='_kou.txt', status='replace')
            do i=1,KOU
                write(12,*) x(noku(i)),' ',z(noku(i))
            end do
            close(12)
        end if

       write(strw1,'(f10.3)') xmin;
       write(strw2,'(f10.3)') xmax;
       write(strw3,'(f10.3)') zmin;
       write(strw4,'(f10.3)') zmax;
       srange=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
       write(strw1,'(f10.3)') baseheight/(zmax-zmin)*(xmax-xmin);
       write(strw2,'(f10.3)') baseheight;
       sscale=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
       sbmap=trim(adjustl('psbasemap -R'//trim(adjustl(srange))//' -JX'//trim(adjustl(sscale))//&
              ' -Ba5f1:"Distance (m)":/a5f1:"Altitude (m)":WS -P -K > %fig%'))

        open (12, file='bat_gmt_mesh.bat', status='replace')
        write(12,'(a)') 'gmtset ANOT_FONT_SIZE 10'
        write(12,'(a)') 'gmtset ANOT_OFFSET 0.2c'
        write(12,'(a)') 'gmtset LABEL_FONT_SIZE 10'
        write(12,'(a)') 'gmtset LABEL_OFFSET 0.1c'
        write(12,'(a)') 'set fig=fig_gmt_mesh.eps'
        write(12,'(a)') trim(adjustl(sbmap))
        do ma=1,MATEL
            write (filename, '("_mesh", i2.2, ".txt")') ma
            if(ma==1)scol='255/255/204'
            if(ma==2)scol='0/255/255'
            if(ma==3)scol='255/0/255'
            if(ma==4)scol='0/255/255'
            if(ma==5)scol='255/255/204'
            spsxy='psxy '//trim(adjustl(filename))//' -R -J -W0.5/black -G'//trim(adjustl(scol))//' -K -O >> %fig%'
            write(12,'(a)') trim(adjustl(spsxy))
        end do
        if(0<KOH)then
            write(12,'(a)') 'psxy _koh.txt -R -J -SS0.15 -W2/black -N -K -O >> %fig%'
        end if
        if(0<KOQ)then
            write(12,'(a)') 'psxy _koq.txt -R -J -SC0.15 -W2/black -N -K -O >> %fig%'
        end if
        if(0<KOU)then
            write(12,'(a)') 'psxy _kou.txt -R -J -ST0.1 -W2/black -N -K -O >> %fig%'
        end if
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'(a)') 'set fig='
        write(12,'(a)') 'del .gmt*'
        close(12)

        call system('bat_gmt_mesh')
    end subroutine DRAW_MESH


    subroutine DRAW_CONT(kakom,x,z,bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4,hvec,pvec,nod,NODT,NELT)
        integer::nod,NODT,NELT
        integer::kakom(1:NELT,1:nod)
        real(4)::x(1:NODT),z(1:NODT),hvec(1:NODT),pvec(1:NODT)
        real(4)::bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4

        integer::i,j,k,l,ne
        real(4)::xmin,xmax,zmin,zmax
        character::srange*100,sscale*100,sbmap*200,spsxy*200,filename*50,scol*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::baseheight=4.0

        open (12, file='_cp.cpt', status='replace')
        write(12,'(a)') '-20 255 255 190 0 255 255 125'
        write(12,'(a)') '0 0 255 255 20 0 125 255'
        write(12,'(a)') 'N 255 255 255 N 255 255 255'
        close(12)

        open (12, file='_mask.txt', status='replace')
        write(12,*) bdx1,' ',bdz1
        write(12,*) bdx2,' ',bdz2
        write(12,*) bdx3,' ',bdz3
        write(12,*) bdx4,' ',bdz4
        close(12)

        open (12, file='_temp1.txt', status='replace')
        do i=1,NODT
            write(12,*) x(i),' ',z(i),' ',pvec(i)
        end do
        close(12)
        open (12, file='_temp2.txt', status='replace')
        do i=1,NODT
            write(12,*) x(i),' ',z(i),' ',hvec(i)
        end do
        close(12)

        xmin=(x(1))
        xmax=(x(1))
        zmin=(z(1))
        zmax=(z(1))
        do i=1,NODT
            if((x(i))<xmin)xmin=(x(i))
            if(xmax<(x(i)))xmax=(x(i))
            if((z(i))<zmin)zmin=(z(i))
            if(zmax<(z(i)))zmax=(z(i))
        end do

       write(strw1,'(f10.3)') xmin;
       write(strw2,'(f10.3)') xmax;
       write(strw3,'(f10.3)') zmin;
       write(strw4,'(f10.3)') zmax;
       srange=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
       write(strw1,'(f10.3)') baseheight/(zmax-zmin)*(xmax-xmin);
       write(strw2,'(f10.3)') baseheight;
       sscale=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
       sbmap=trim(adjustl('psbasemap -R'//trim(adjustl(srange))//' -JX'//trim(adjustl(sscale))//&
              ' -Ba5f1:"Distance (m)":/a5f1:"Altitude (m)":WS -P -K > %fig%'))

        open (12, file='bat_gmt_cont.bat', status='replace')
        write(12,'(a)') 'gmtset ANOT_FONT_SIZE 10'
        write(12,'(a)') 'gmtset ANOT_OFFSET 0.2c'
        write(12,'(a)') 'gmtset LABEL_FONT_SIZE 10'
        write(12,'(a)') 'gmtset LABEL_OFFSET 0.1c'
        write(12,'(a)') 'set fig=fig_gmt_cont.eps'
        write(12,'(a)') trim(adjustl(sbmap))
        write(12,'(a)') 'grdmask _mask.txt -G_mask.grd -I0.1 -R -NNaN/0/0'
        write(12,'(a)') 'surface _temp1.txt -G_temp.grd -I0.1 -R -T0.0'
        write(12,'(a)') 'grdmath _temp.grd _mask.grd OR = _masked_data.grd'
        write(12,'(a)') 'grdimage _masked_data.grd -J -C_cp.cpt -K -O >> %fig%'
        write(12,'(a)') 'grdcontour _masked_data.grd -C1 -A5 -J -K -O >> %fig%'
        write(12,'(a)') 'surface _temp2.txt -G_temp.grd -I0.1 -R -T0.0'
        write(12,'(a)') 'grdmath _temp.grd _mask.grd OR = _masked_data.grd'
        write(12,'(a)') 'grdcontour _masked_data.grd -C1 -A5 -J -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'(a)') 'set fig='
        write(12,'(a)') 'del .gmt*'
        close(12)

        call system('bat_gmt_cont')
    end subroutine DRAW_CONT


    subroutine DRAW_VECT(kakom,x,z,bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4,pvec,vx,vz,nod,NODT,NELT)
        integer::nod,NODT,NELT
        integer::kakom(1:NELT,1:nod)
        real(4)::x(1:NODT),z(1:NODT),pvec(1:NODT)
        real(4)::vx(1:NELT),vz(1:NELT)
        real(4)::bdx1,bdz1,bdx2,bdz2,bdx3,bdz3,bdx4,bdz4

        integer::i,j,k,l,ne
        real(4)::xmin,xmax,zmin,zmax,vmax,vscale,xx,zz
        character::srange*100,sscale*100,sbmap*200,spsxy*200,filename*50,scol*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::baseheight=4.0

        open (12, file='_cp.cpt', status='replace')
        write(12,'(a)') '-20 255 255 190 0 255 255 125'
        write(12,'(a)') '0 0 255 255 20 0 125 255'
        write(12,'(a)') 'N 255 255 255 N 255 255 255'
        close(12)

        open (12, file='_mask.txt', status='replace')
        write(12,*) bdx1,' ',bdz1
        write(12,*) bdx2,' ',bdz2
        write(12,*) bdx3,' ',bdz3
        write(12,*) bdx4,' ',bdz4
        close(12)

        open (12, file='_temp1.txt', status='replace')
        do i=1,NODT
            write(12,*) x(i),' ',z(i),' ',pvec(i)
        end do
        close(12)

        xmin=(x(1))
        xmax=(x(1))
        zmin=(z(1))
        zmax=(z(1))
        do i=1,NODT
            if((x(i))<xmin)xmin=(x(i))
            if(xmax<(x(i)))xmax=(x(i))
            if((z(i))<zmin)zmin=(z(i))
            if(zmax<(z(i)))zmax=(z(i))
        end do
        vmax=0.0
        do ne=1,NELT
            if(vmax<sqrt(vx(ne)**2.0+vz(ne)**2.0))vmax=sqrt(vx(ne)**2.0+vz(ne)**2.0)
        end do
        if(xmax-xmin<zmax-zmin)then
            vscale=(zmax-zmin)/20.0/vmax
        else
            vscale=(xmax-xmin)/20.0/vmax
        end if
            print "('vmax=',e15.7)",vmax

        open (12, file='_vect.txt', status='replace')
        select case(nod)
            case(3)
                do ne=1,NELT,3
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    xx=(x(i)+x(j)+x(k))/3.0
                    zz=(z(i)+z(j)+z(k))/3.0
                    write(12,*) xx-0.5*vscale*vx(ne),' ',zz-0.5*vscale*vz(ne),' ',&
                                xx+0.5*vscale*vx(ne),' ',zz+0.5*vscale*vz(ne)
                end do
            case(4)
                do ne=1,NELT,3
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    l=kakom(ne,4)
                    xx=0.25*(x(i)+x(j)+x(k)+x(l))
                    zz=0.25*(z(i)+z(j)+z(k)+z(l))
                    write(12,*) xx-0.5*vscale*vx(ne),' ',zz-0.5*vscale*vz(ne),' ',&
                                xx+0.5*vscale*vx(ne),' ',zz+0.5*vscale*vz(ne)
                end do
        end select
        close(12)

       write(strw1,'(f10.3)') xmin;
       write(strw2,'(f10.3)') xmax;
       write(strw3,'(f10.3)') zmin;
       write(strw4,'(f10.3)') zmax;
       srange=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
       write(strw1,'(f10.3)') baseheight/(zmax-zmin)*(xmax-xmin);
       write(strw2,'(f10.3)') baseheight;
       sscale=trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
       sbmap=trim(adjustl('psbasemap -R'//trim(adjustl(srange))//' -JX'//trim(adjustl(sscale))//&
              ' -Ba5f1:"Distance (m)":/a5f1:"Altitude (m)":WS -P -K > %fig%'))

        open (12, file='bat_gmt_vect.bat', status='replace')
        write(12,'(a)') 'gmtset ANOT_FONT_SIZE 10'
        write(12,'(a)') 'gmtset ANOT_OFFSET 0.2c'
        write(12,'(a)') 'gmtset LABEL_FONT_SIZE 10'
        write(12,'(a)') 'gmtset LABEL_OFFSET 0.1c'
        write(12,'(a)') 'gmtset VECTOR_SHAPE=1'
        write(12,'(a)') 'set fig=fig_gmt_vect.eps'
        write(12,'(a)') trim(adjustl(sbmap))
        write(12,'(a)') 'grdmask _mask.txt -G_mask.grd -I0.1 -R -NNaN/0/0'
        write(12,'(a)') 'surface _temp1.txt -G_temp.grd -I0.1 -R -T0.0'
        write(12,'(a)') 'grdmath _temp.grd _mask.grd OR = _masked_data.grd'
        write(12,'(a)') 'grdimage _masked_data.grd -J -C_cp.cpt -K -O >> %fig%'
        write(12,'(a)') 'psxy _vect.txt -R -J -Svs0.005/0.1/0.03 -G0 -N -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'(a)') 'set fig='
        write(12,'(a)') 'del .gmt*'
        close(12)

        call system('bat_gmt_vect')
    end subroutine DRAW_VECT


    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


end program f90_GMTseep