module cdraw_data
    implicit none
    integer::ind,ibc,iph,ivc,isc
    real(4)::baselength,grdspc
    real(4)::xmin,xmax,zmin,zmax
    real(4)::bav,bfv
    real(4)::conint,annint
    character::strx*100,strz*100
end module cdraw_data


program f90_gmt_seep
    !Drawing by GMT for seepage analysis
    use cdraw_data
    implicit none
    integer::i,j,k,l,ne
    integer::n1,n2,n3,k1,k2,k3
    character::strcom*256           !Comment
    integer::NODT                   !Number of nodes
    integer::NELT                   !Number of elements
    integer::MATEL                  !Number of material sets
    integer::KOH                    !Number of nodes with given toal head
    integer::KOQ                    !Number of nodes with given discharge
    integer::KOU                    !Number of nodes with the condition of seepage point
    integer::nt                     !Number of total degree of freedom (NODT x nhen)
    integer::nod                    !Number of nodes per element
    integer::nhen=1                 !Number of degree of freedom per node
    integer,allocatable::kakom(:,:) !Element connectivity
    integer,allocatable::matno(:)   !Material set number
    real(4),allocatable::x(:)       !x-coordinate of node
    real(4),allocatable::z(:)       !z-coordinate of node
    real(4),allocatable::hvec(:)    !Total head vector
    real(4),allocatable::qvec(:)    !Discharge vector
    real(4),allocatable::pvec(:)    !Pressure head vector
    integer,allocatable::nokh(:)
    integer,allocatable::nokq(:)
    integer,allocatable::noku(:)
    real(4),allocatable::xg(:)      !x-coordinate of centroid of element
    real(4),allocatable::zg(:)      !z-coordinate of centroid of element
    real(4),allocatable::vx(:)      !Velocity vector in x-direction
    real(4),allocatable::vz(:)      !velocity vector in z-direction
    real(4),allocatable::Akr(:)     !Ratio of permeability (unsat/sat)
    integer::iw
    real(4)::rw

    character::linebuf*1000
    character::filename*128
    character::fnameR*50

    call getarg(1,fnameR) !Input file name

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,*) strcom
    read(11,'()')
    read(11,*) nod,NODT,NELT,MATEL,KOH,KOQ,KOU
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    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(xg(1:NELT))
    allocate(zg(1:NELT))
    allocate(vx(1:NELT))
    allocate(vz(1:NELT))
    allocate(Akr(1:NELT))
    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,*) iw,(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,*) iw,x(i),z(i),hvec(i),qvec(i),pvec(i)
    end do
    read(11,'()')
    read(11,'()')
    !element,xg,zg,vx,vz,vm,Ak0,Akr,matno
    do ne=1,NELT
        read(11,*) iw,xg(ne),zg(ne),vx(ne),vz(ne),rw,rw,Akr(ne),iw
    end do
    close(11)

    open(11,file='_inp_base.txt',status='old')
    read(11,*) ind,ibc,iph,ivc,isc
    read(11,*) baselength,grdspc
    read(11,*) bav,bfv
    read(11,*) conint,annint
    read(11,*) strx
    read(11,*) strz
    close(11)

    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

    if(0<ind)call DAT_IND(kakom,x,z,NELT,NODT,nod)
    call DAT_MASK(kakom,x,z,nod,NODT,NELT)
    call DAT_MESH(kakom,matno,x,z,nokh,nokq,noku,nod,NODT,NELT,MATEL,KOH,KOQ,KOU)
    call DAT_CONT(kakom,x,z,hvec,pvec,nod,NODT,NELT)
    call DAT_VECT(xg,zg,vx,vz,NELT)
    call MAKE_BAT(MATEL,KOH,KOQ,KOU)
    call system('bash _gmt_seep.txt')

stop


contains

    real(4) function AREA(x1,y1,x2,y2,x3,y3)
        real(4),intent(in)::x1,y1,x2,y2,x3,y3
        AREA=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    end function AREA


    subroutine DAT_IND(kakom,x,z,NELT,NODT,nod)
        integer,intent(in)::NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(4),intent(in)::x(1:NODT),z(1:NODT)

        integer::ne,i,j,k,l
        real(4)::xg,zg,dts

        dts=0.001*baselength

        open (12, file='_ne_circ.txt', status='replace')
        open (13, file='_text_ne.txt', status='replace')
        select case(nod)
            case(3)
                do ne=1,NELT
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    xg=(x(i)+x(j)+x(k))/3.0
                    zg=(z(i)+z(j)+z(k))/3.0
                    write(12,*) xg,' ',zg
                    write(13,*) xg,' ',zg,' 9 0 MC ',ne
                end do
            case(4)
                do ne=1,NELT
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    l=kakom(ne,4)
                    xg=0.25*(x(i)+x(j)+x(k)+x(l))
                    zg=0.25*(z(i)+z(j)+z(k)+z(l))
                    write(12,*) xg,' ',zg
                    write(13,*) xg,' ',zg,' 9 0 MC ',ne
                end do
        end select
        close(12)
        close(13)

        open (12, file='_text_nd.txt', status='replace')
        do i=1,NODT
            write(12,*) x(i)-dts,' ',z(i)-dts,' 9 0 TR ',i
        end do
        close(12)
    end subroutine DAT_IND


    subroutine DAT_MASK(kakom,x,z,nod,NODT,NELT)
        integer::nod,NODT,NELT
        integer::kakom(1:NELT,1:nod)
        real(4)::x(1:NODT),z(1:NODT)

        integer::i,j,k,l,ne,ii,jj,nnn1,nnn2
        real(4)::a,b,area1,area2,area3,area4

        nnn1=int((xmax-xmin)/grdspc)
        nnn2=int((zmax-zmin)/grdspc)
        open (12, file='_mask_data.txt', status='replace')
        select case(nod)
            case(3)
                do ii=0,nnn1
                do jj=0,nnn2
                    a=xmin+grdspc*real(ii)
                    b=zmin+grdspc*real(jj)
                    do ne=1,NELT
                        i=kakom(ne,1)
                        j=kakom(ne,2)
                        k=kakom(ne,3)
                        area1=AREA(a,b,x(i),z(i),x(j),z(j))
                        area2=AREA(a,b,x(j),z(j),x(k),z(k))
                        area3=AREA(a,b,x(k),z(k),x(i),z(i))
                        if(0.0<=area1.and.0.0<=area2.and.0.0<=area3)then
                            write(12,*) a,' ',b
                            exit
                        end if
                    end do
                end do
                end do
            case(4)
                do ii=0,nnn1
                do jj=0,nnn2
                    a=xmin+grdspc*real(ii)
                    b=zmin+grdspc*real(jj)
                    do ne=1,NELT
                        i=kakom(ne,1)
                        j=kakom(ne,2)
                        k=kakom(ne,3)
                        l=kakom(ne,4)
                        area1=AREA(a,b,x(i),z(i),x(j),z(j))
                        area2=AREA(a,b,x(j),z(j),x(k),z(k))
                        area3=AREA(a,b,x(k),z(k),x(l),z(l))
                        area4=AREA(a,b,x(l),z(l),x(i),z(i))
                        if(0.0<=area1.and.0.0<=area2.and.0.0<=area3.and.0.0<=area4)then
                            write(12,*) a,' ',b
                            exit
                        end if
                    end do
                end do
                end do
        end select
        close (12)
    end subroutine DAT_MASK


    subroutine DAT_MESH(kakom,matno,x,z,nokh,nokq,noku,nod,NODT,NELT,MATEL,KOH,KOQ,KOU)
        integer::nod,NODT,NELT,MATEL,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

        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,'(a)') '>'
                            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,'(a)') '>'
                            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
    end subroutine DAT_MESH


    subroutine DAT_CONT(kakom,x,z,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)

        integer::i

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


    subroutine DAT_VECT(xg,zg,vx,vz,NELT)
        integer::NELT
        real(4)::xg(1:NELT),zg(1:NELT)
        real(4)::vx(1:NELT),vz(1:NELT)

        integer::ne
        real(4)::vmax,vscale

        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')
        do ne=1,NELT,3
            write(12,*) xg(ne)-0.5*vscale*vx(ne),' ',zg(ne)-0.5*vscale*vz(ne),' ',&
                        xg(ne)+0.5*vscale*vx(ne),' ',zg(ne)+0.5*vscale*vz(ne)
        end do
        close(12)
    end subroutine DAT_VECT


    subroutine MAKE_BAT(MATEL,KOH,KOQ,KOU)
        integer::MATEL,KOH,KOQ,KOU

        integer::ma
        character::srange*100,sscale*100,sba*100,sgrdspc*100,ssradi*100,filename*100
        character::ssx*100,ssz*100,scc*50,sca*50,sscld*100,spsize*50,scsize*50,scol*50,snum*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::psize,work,ds
        character(len=3)::col_r(1:MATEL),col_g(1:MATEL),col_b(1:MATEL)
        character(len=11)::col_mesh(1:MATEL)

        psize=0.015*baselength
        ds=bfv

        open (11, file='_col_mesh.txt', status='old')
        do i=1,MATEL
            read(11,*) col_r(i),col_g(i),col_b(i)
            col_mesh(i)=trim(adjustl(col_r(i)))//'/'//trim(adjustl(col_g(i)))//'/'//trim(adjustl(col_b(i)))
        end do
        close(11)

        write(strw1,'(f10.3)') xmin-ds;
        write(strw2,'(f10.3)') xmax;
        write(strw3,'(f10.3)') zmin-ds;
        write(strw4,'(f10.3)') zmax;
        srange='range='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
        if(zmax-zmin<=xmax-xmin)then
            write(strw1,'(f10.3)') baselength;
            write(strw2,'(f10.3)') baselength/(xmax-xmin+ds)*(zmax-zmin+ds)
        else
            write(strw1,'(f10.3)') baselength/(zmax-zmin+ds)*(xmax-xmin+ds);
            write(strw2,'(f10.3)') baselength;
        end if
        sscale='scale='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') bav
        write(strw2,'(f10.3)') bfv
        sba='ba=a'//trim(adjustl(strw1))//'f'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') grdspc
        sgrdspc='grdspc='//trim(adjustl(strw1))
        write(strw1,'(f11.4)') grdspc/sqrt(2.0)
        ssradi='sradi='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') conint
        scc='conint='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') annint
        sca='annint='//trim(adjustl(strw1))
        ssx='strx="'//trim(adjustl(strx))//'"'
        ssz='strz="'//trim(adjustl(strz))//'"'
        if(zmax-zmin<=xmax-xmin)then
            work=baselength/(xmax-xmin+ds)*(zmax-zmin+ds)
            write(strw1,'(f10.3)') baselength+0.5
            write(strw2,'(f10.3)') work-work*0.5/2.0-0.3
            write(strw3,'(f10.3)') work*0.5
            if(4.0<work*0.5)then
                write(strw2,'(f10.3)') work-4.0/2.0-0.3
                write(strw3,'(f10.3)') 4.0
            end if
        else
            work=baselength/(zmax-zmin+ds)*(xmax-xmin+ds)
            write(strw1,'(f10.3)') work+0.5
            write(strw2,'(f10.3)') baselength-work*0.5/2.0-0.3
            write(strw3,'(f10.3)') work*0.5
            if(4.0<work*0.5)then
                write(strw2,'(f10.3)') work-4.0/2.0-0.3
                write(strw3,'(f10.3)') 4.0
            end if
        end if
        sscld='scld='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/0.3'
        write(strw1,'(f10.3)') psize
        spsize='psize='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') psize*4.0
        scsize='csize='//trim(adjustl(strw1))

        !Writing of batch file
        open (12, file='_gmt_seep.txt', status='replace')
        write(12,'(a)') 'gmt set FONT_ANNOT_PRIMARY 10'
        write(12,'(a)') 'gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c'
        write(12,'(a)') 'gmt set FONT_LABEL 10p'
        write(12,'(a)') 'gmt set MAP_LABEL_OFFSET 0.3c'
        write(12,'(a)') 'gmt set MAP_TICK_LENGTH_PRIMARY 0.2c'
        write(12,'(a)') trim(adjustl(srange))
        write(12,'(a)') trim(adjustl(sscale))
        write(12,'(a)') trim(adjustl(sba))
        write(12,'(a)') trim(adjustl(sgrdspc))
        write(12,'(a)') trim(adjustl(ssradi))
        write(12,'(a)') trim(adjustl(scc))
        write(12,'(a)') trim(adjustl(sca))
        write(12,'(a)') trim(adjustl(ssx))
        write(12,'(a)') trim(adjustl(ssz))
        write(12,'(a)') trim(adjustl(sscld))
        write(12,'(a)') trim(adjustl(spsize))
        write(12,'(a)') trim(adjustl(scsize))
        write(12,'()')
        !Drawing of mesh
        write(12,'(a)') '# *** Drawing of mesh ***'
        write(12,'(a)') 'fig=_fig_gmt_mesh.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$strz" -BWS -P -K > $fig'
        do ma=1,MATEL
            write (strw1,'(i2.2)') ma
            snum='num='//trim(adjustl(strw1))
            scol='col='//col_mesh(ma)
            write(12,'(a)') trim(adjustl(snum))
            write(12,'(a)') trim(adjustl(scol))
            write(12,'(a)') 'gmt psxy _mesh${num}.txt -R -J -W0.5,black -G$col -K -O >> $fig'
        end do
        if(0<ind)then
            write(12,'(a)') 'gmt psxy _ne_circ.txt -R -J -SC$csize -N -K -O >> $fig'
            write(12,'(a)') 'gmt pstext _text_ne.txt -R -J -F+f+a+j -N -K -O >> $fig'
            write(12,'(a)') 'gmt pstext _text_nd.txt -R -J -F+f+a+j -N -K -O >> $fig'
        end if
        if(0<ibc)then
            if(0<KOH)write(12,'(a)') 'gmt psxy _koh.txt -R -J -SS$psize -W1,black -N -K -O >> $fig'
            if(0<KOQ)write(12,'(a)') 'gmt psxy _koq.txt -R -J -SC$psize -W1,black -N -K -O >> $fig'
            if(0<KOU)write(12,'(a)') 'gmt psxy _kou.txt -R -J -St$psize -W1,black -N -K -O >> $fig'
        end if
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        !Preparing of contour drawing
        write(12,'(a)') 'gmt grdmask _mask_data.txt -G_mask.grd -I$grdspc -R$range -S$sradi -NNan/0/0'
        write(12,'(a)') 'gmt surface _pvec.txt -G_pvec.grd -I$grdspc -R$range -T0.0'
        write(12,'(a)') 'gmt surface _hvec.txt -G_hvec.grd -I$grdspc -R$range -T0.0'
        write(12,'(a)') 'gmt grdmath _pvec.grd _mask.grd OR = _pvecm.grd'
        write(12,'(a)') 'gmt grdmath _hvec.grd _mask.grd OR = _hvecm.grd'
        !Drawing of contour for pressure head
        write(12,'(a)') '# *** Drawing of contour ***'
        write(12,'(a)') 'fig=_fig_gmt_conp.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$strz" -BWS -P -K > $fig'
        write(12,'(a)') 'gmt grdimage _pvecm.grd -J -C0_cpt.cpt -K -O >> $fig'
        write(12,'(a)') 'gmt grdcontour _pvecm.grd -C$conint -A$annint -J -K -O >> $fig'
        if(iph==1)write(12,'(a)') 'gmt grdcontour _hvecm.grd -C$conint -A$annint -J -K -O >> $fig'
        if(0<isc)write(12,'(a)') 'gmt psscale -Ba$annint -C0_cpt.cpt -D$scld -E0.3 -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        !Drawing of contour for total head
        write(12,'(a)') 'fig=_fig_gmt_conh.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$strz" -BWS -P -K > $fig'
        write(12,'(a)') 'gmt grdimage _hvecm.grd -J -C0_cpt.cpt -K -O >> $fig'
        write(12,'(a)') 'gmt grdcontour _hvecm.grd -C$conint -A$annint -J -K -O >> $fig'
        if(iph==1)write(12,'(a)') 'gmt grdcontour _pvecm.grd -C$conint -A$annint -J -K -O >> $fig'
        if(0<isc)write(12,'(a)') 'gmt psscale -Ba$annint -C0_cpt.cpt -D$scld -E0.3 -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        !Drawing of vector
        write(12,'(a)') '# *** Drawing of vector ***'
        write(12,'(a)') 'fig=_fig_gmt_vect.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$strz" -BWS -P -K > $fig'
        if(ivc==1)write(12,'(a)') 'gmt grdimage _pvecm.grd -J -C0_cpt.cpt -K -O >> $fig'
        if(ivc==2)write(12,'(a)') 'gmt grdimage _hvecm.grd -J -C0_cpt.cpt -K -O >> $fig'
        write(12,'(a)') 'gmt psxy _vect.txt -R -J -Sv0.1+s+e+a30 -G0 -W0.5 -N -K -O >> $fig'
        if(0<isc)write(12,'(a)') 'gmt psscale -Ba$annint -C0_cpt.cpt -D$scld -E0.3 -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        close(12)
    end subroutine MAKE_BAT


    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_gmt_seep
