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_GMTseep
    !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('bat_gmt_seep')

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 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 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 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,*) '>'
                            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
    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='set 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='set scale='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') bav
        write(strw2,'(f10.3)') bfv
        sba='set ba=a'//trim(adjustl(strw1))//'f'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') grdspc
        sgrdspc='set grdspc='//trim(adjustl(strw1))
        write(strw1,'(f11.4)') grdspc/sqrt(2.0)
        ssradi='set sradi='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') conint
        scc='set conint='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') annint
        sca='set annint='//trim(adjustl(strw1))
        ssx='set strx="'//trim(adjustl(strx))//'"'
        ssz='set 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='set scld='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/0.3'
        write(strw1,'(f10.3)') psize
        spsize='set psize='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') psize*4.0
        scsize='set csize='//trim(adjustl(strw1))

        !Writing of batch file
        open (12, file='bat_gmt_seep.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)') 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)') 'rem *** Drawing of mesh ***'
        write(12,'(a)') 'set fig=fig_gmt_mesh.eps'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%strz%:WS -P -K > %fig%'
        do ma=1,MATEL
            write (strw1,'(i2.2)') ma
            snum='set num='//trim(adjustl(strw1))
            scol='set col='//col_mesh(ma)
            write(12,'(a)') trim(adjustl(snum))
            write(12,'(a)') trim(adjustl(scol))
            write(12,'(a)') 'psxy _mesh%num%.txt -R -J -W0.5/black -G%col% -K -O >> %fig%'
        end do
        if(0<ind)then
            write(12,'(a)') 'psxy _ne_circ.txt -R -J -SC%csize% -N -K -O >> %fig%'
            write(12,'(a)') 'pstext _text_ne.txt -R -J -N -K -O >> %fig%'
            write(12,'(a)') 'pstext _text_nd.txt -R -J -N -K -O >> %fig%'
        end if
        if(0<ibc)then
            if(0<KOH)write(12,'(a)') 'psxy _koh.txt -R -J -SS%psize% -W2/black -N -K -O >> %fig%'
            if(0<KOQ)write(12,'(a)') 'psxy _koq.txt -R -J -SC%psize% -W2/black -N -K -O >> %fig%'
            if(0<KOU)write(12,'(a)') 'psxy _kou.txt -R -J -St%psize% -W2/black -N -K -O >> %fig%'
        end if
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'()')
        !Preparing of contour drawing
        write(12,'(a)') 'grdmask _mask_data.txt -G_mask.grd -I%grdspc% -R%range% -S%sradi% -NNan/0/0'
        write(12,'(a)') 'surface _pvec.txt -G_pvec.grd -I%grdspc% -R%range% -T0.0'
        write(12,'(a)') 'surface _hvec.txt -G_hvec.grd -I%grdspc% -R%range% -T0.0'
        write(12,'(a)') 'grdmath _pvec.grd _mask.grd OR = _pvecm.grd'
        write(12,'(a)') 'grdmath _hvec.grd _mask.grd OR = _hvecm.grd'
        !Drawing of contour for pressure head
        write(12,'(a)') 'rem *** Drawing of contour ***'
        write(12,'(a)') 'set fig=fig_gmt_conp.eps'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%strz%:WS -P -K > %fig%'
        write(12,'(a)') 'grdimage _pvecm.grd -J -C0_cpt.cpt -K -O >> %fig%'
        write(12,'(a)') 'grdcontour _pvecm.grd -C%conint% -A%annint% -J -K -O >> %fig%'
        if(iph==1)write(12,'(a)') 'grdcontour _hvecm.grd -C%conint% -A%annint% -J -K -O >> %fig%'
        if(0<isc)write(12,'(a)') 'psscale -Ba%annint% -C0_cpt.cpt -D%scld% -E0.3 -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        !Drawing of contour for total head
        write(12,'(a)') 'set fig=fig_gmt_conh.eps'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%strz%:WS -P -K > %fig%'
        write(12,'(a)') 'grdimage _hvecm.grd -J -C0_cpt.cpt -K -O >> %fig%'
        write(12,'(a)') 'grdcontour _hvecm.grd -C%conint% -A%annint% -J -K -O >> %fig%'
        if(iph==1)write(12,'(a)') 'grdcontour _pvecm.grd -C%conint% -A%annint% -J -K -O >> %fig%'
        if(0<isc)write(12,'(a)') 'psscale -Ba%annint% -C0_cpt.cpt -D%scld% -E0.3 -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'()')
        !Drawing of vector
        write(12,'(a)') 'rem *** Drawing of vector ***'
        write(12,'(a)') 'set fig=fig_gmt_vect.eps'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%strz%:WS -P -K > %fig%'
        if(ivc==1)write(12,'(a)') 'grdimage _pvecm.grd -J -C0_cpt.cpt -K -O >> %fig%'
        if(ivc==2)write(12,'(a)') 'grdimage _hvecm.grd -J -C0_cpt.cpt -K -O >> %fig%'
        write(12,'(a)') 'psxy _vect.txt -R -J -Svs0.005/0.1/0.03 -G0 -N -K -O >> %fig%'
        if(0<isc)write(12,'(a)') 'psscale -Ba%annint% -C0_cpt.cpt -D%scld% -E0.3 -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'()')
        !reset of variables
        write(12,'(a)') 'set fig='
        write(12,'(a)') 'set range='
        write(12,'(a)') 'set scale='
        write(12,'(a)') 'set ba='
        write(12,'(a)') 'set grdspc='
        write(12,'(a)') 'set sradi='
        write(12,'(a)') 'set conint='
        write(12,'(a)') 'set annint='
        write(12,'(a)') 'set strx='
        write(12,'(a)') 'set strz='
        write(12,'(a)') 'set scld='
        write(12,'(a)') 'set psize='
        write(12,'(a)') 'set csize='
        write(12,'(a)') 'del .gmt*'
        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_GMTseep