module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi


module cdraw_data
    implicit none
    integer::ind,ibc,ifc,isc
    real(4)::baselength,grdspc
    real(4)::xmin,xmax,ymin,ymax
    real(4)::bav,bfv
    real(4)::conint,annint
    character::strx*100,stry*100
    integer::ncut,numcut(10)

    real(4)::fmax,fact0=1.0E-10
end module cdraw_data


program f90_GMTplnt
    !Drawing by GMT for 2D Stress Analysis
    use defpi
    use cdraw_data
    implicit none
    integer::i,j,k,l,kk,ne,n,ia,ja,k1,k2,n1,n2,n3
    character::strcom*50                !Comment
    integer::NSTRES                     !Stress state (0: plane strain, 1: plane stress)
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of restricted nodes in x-direction
    integer::KOY                        !Number of restricted nodes in y-direction
    integer::KXY                        !Number of restricted nodes in x and y-directions
    integer::NF                         !Number of loaded nodes
    integer,allocatable::kakom(:,:)     !Element connectivity
    integer,allocatable::matno(:)       !Material set number
    real(4),allocatable::x(:)           !x-coordinate of node
    real(4),allocatable::y(:)           !y-coordinate of node
    real(4),allocatable::grdspcT(:)     !Temperature change of node (positive: temperature rise)
    integer,allocatable::nokx(:)        !Restricted node number in x-direction
    integer,allocatable::noky(:)        !Restricted node number in y-direction
    integer,allocatable::noxy(:)        !Restricted node number in x and y-directions
    real(4),allocatable::f(:)           !External load vector
    real(4),allocatable::dis(:)         !Displacement vector
    real(4),allocatable::strsave(:,:,:) !Stresses of element
    integer,allocatable::noten(:,:)     !Stress state (notension or not)
    integer::nt                         !Total degree of freedom (NODT x nhen)
    integer::nod                        !Number of nodes per element
    integer::nhen=2                     !Number of degree of freedom per node
    integer::kkmax                      !Number of integration point
    real(4),allocatable::xg(:,:)        !x-coordinate of centroid of element
    real(4),allocatable::yg(:,:)        !y-coordinate of centroid of element
    real(4)::rw
    integer::iw

    integer::IPR
    character::linebuf*1000
    character::fnameR*50

    call getarg(1,fnameR) !Input file name

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,'()')
    read(11,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    if(nod==3)IPR=0
    nt=NODT*nhen
    if(nod==3)kkmax=1
    if(nod==4)kkmax=4
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(grdspcT(1:NODT))
    allocate(nokx(1:KOX))
    allocate(noky(1:KOY))
    allocate(noxy(1:KOX))
    allocate(f(1:nt))
    allocate(dis(1:nt))
    allocate(strsave(1:NELT,1:kkmax,1:3))
    allocate(noten(1:NELT,1:kkmax))
    allocate(xg(1:NELT,1:kkmax))
    allocate(yg(1:NELT,1:kkmax))
    read(11,'()')
    read(11,'()')
    n1=0
    n2=0
    n3=0
    do i=1,NODT
        read(11,*) iw,rw,rw,f(2*i-1),f(2*i),k1,k2,rw,rw,grdspcT(i)
        if(0<KOX)then
            if(k1==1.and.k2==0)then
                n1=n1+1
                nokx(n1)=i
            end if
        end if
        if(0<KOY)then
            if(k1==0.and.k2==1)then
                n2=n2+1
                noky(n2)=i
            end if
        end if
        if(0<KOX.and.0<KOY)then
            if(k1==1.and.k2==1)then
                n3=n3+1
                noxy(n3)=i
            end if
        end if
    end do
    KOX=n1
    KOY=n2
    KXY=n3
    read(11,'()')
    read(11,'()')
    do ne=1,NELT
        read(11,*) iw,(kakom(ne,j),j=1,nod),rw,rw,rw,rw,rw,rw,rw,rw,matno(ne)
    end do
    read(11,'()')
    read(11,'()')
    do i=1,NODT
        read(11,*) iw,x(i),y(i),dis(2*i-1),dis(2*i),rw,rw,rw,rw
    end do
    read(11,'()')
    read(11,'()')
    do ne=1,NELT
        select case(IPR)
            case(0)
                do kk=1,kkmax
                    read(11,*) iw,iw,xg(ne,kk),yg(ne,kk),(strsave(ne,kk,j),j=1,3),rw,rw,rw,noten(ne,kk),matno(ne)
                end do
            case(1)
                kk=1
                read(11,*) iw,iw,xg(ne,kk),yg(ne,kk),(strsave(ne,kk,j),j=1,3),rw,rw,rw,noten(ne,kk),matno(ne)
        end select
    end do
    close(11)

    open(11,file='_inp_base.txt',status='old')
    read(11,*) ind,ibc,ifc,isc
    read(11,*) baselength,grdspc
    read(11,*) bav,bfv
    read(11,*) conint,annint
    read(11,*) strx
    read(11,*) stry
    read(11,*) ncut,(numcut(i),i=1,ncut)
    close(11)

    xmin=(x(1))
    xmax=(x(1))
    ymin=(y(1))
    ymax=(y(1))
    do i=1,NODT
        if((x(i))<xmin)xmin=(x(i))
        if(xmax<(x(i)))xmax=(x(i))
        if((y(i))<ymin)ymin=(y(i))
        if(ymax<(y(i)))ymax=(y(i))
    end do

    if(0<ind)call DAT_IND(kakom,x,y,NELT,NODT,nod)
    call DAT_MASK(kakom,x,y,nod,NODT,NELT)
    call DAT_MESH(kakom,matno,x,y,nokx,noky,noxy,f,nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt)
    call DAT_PS(kakom,matno,x,y,xg,yg,strsave,nod,NODT,NELT,MATEL,IPR,nt,kkmax)
    call MAKE_BAT(MATEL,KOX,KOY,KXY)
    call system('bat_gmt_plnt')

stop

contains

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

        integer::ne,i,j,k,l
        real(4)::xg,yg,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
                    yg=(y(i)+y(j)+y(k))/3.0
                    write(12,*) xg,' ',yg
                    write(13,*) xg,' ',yg,' 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))
                    yg=0.25*(y(i)+y(j)+y(k)+y(l))
                    write(12,*) xg,' ',yg
                    write(13,*) xg,' ',yg,' 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,' ',y(i)-dts,' 9 0 0 TR ',i
        end do
        close(12)
    end subroutine DAT_IND


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

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

        nnn1=int((xmax-xmin)/grdspc)
        nnn2=int((ymax-ymin)/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=ymin+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),y(i),x(j),y(j))
                        area2=AREA(a,b,x(j),y(j),x(k),y(k))
                        area3=AREA(a,b,x(k),y(k),x(i),y(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=ymin+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),y(i),x(j),y(j))
                        area2=AREA(a,b,x(j),y(j),x(k),y(k))
                        area3=AREA(a,b,x(k),y(k),x(l),y(l))
                        area4=AREA(a,b,x(l),y(l),x(i),y(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,y,nokx,noky,noxy,f,nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt)
        integer::nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt
        integer::kakom(1:NELT,1:nod),matno(1:NELT)
        integer::nokx(1:KOX),noky(1:KOY),noxy(1:KOX)
        real(4)::x(1:NODT),y(1:NODT),f(1:nt)

        integer::i,j,k,l,ma,ne
        real(4)::fscale
        character::filename*50

        if(ymax-ymin<xmax-xmin)then
            fscale=0.05*(xmax-xmin)
        else
            fscale=0.05*(ymax-ymin)
        end if

        fmax=0.0
        do i=1,NODT
            if(fmax<sqrt(f(2*i-1)**2.0+f(2*i)**2.0))fmax=sqrt(f(2*i-1)**2.0+f(2*i)**2.0)
        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)),' ',(y(i))
                            write(12,*) (x(j)),' ',(y(j))
                            write(12,*) (x(k)),' ',(y(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)),' ',(y(i))
                            write(12,*) (x(j)),' ',(y(j))
                            write(12,*) (x(k)),' ',(y(k))
                            write(12,*) (x(l)),' ',(y(l))
                        end if
                    end do
            end select
            close (12)
        end do

        if(0<KOX)then
            open (12, file='_kox.txt', status='replace')
            do i=1,KOX
                write(12,*) x(nokx(i)),' ',y(nokx(i))
            end do
            close(12)
        end if
        if(0<KOY)then
            open (12, file='_koy.txt', status='replace')
            do i=1,KOY
                write(12,*) x(noky(i)),' ',y(noky(i))
            end do
            close(12)
        end if
        if(0<KXY)then
            open (12, file='_kxy.txt', status='replace')
            do i=1,KXY
                write(12,*) x(noxy(i)),' ',y(noxy(i))
            end do
            close(12)
        end if

        if(fact0<fmax)then
            open (12, file='_force.txt', status='replace')
            do i=1,NODT
                if(fact0<sqrt(f(2*i-1)**2.0+f(2*i)**2.0))then
                    write(12,*) x(i),' ',y(i),' ',x(i)+f(2*i-1)/fmax*fscale,' ',y(i)+f(2*i)/fmax*fscale
                end if
            end do
            close(12)
        end if
    end subroutine DAT_MESH


    subroutine DAT_PS(kakom,matno,x,y,xg,yg,strsave,nod,NODT,NELT,MATEL,IPR,nt,kkmax)
        integer::nod,NODT,NELT,MATEL,IPR,nt,kkmax
        integer::kakom(1:NELT,1:nod),matno(1:NELT)
        real(4)::x(1:NODT),y(1:NODT)
        real(4)::xg(1:NELT,1:kkmax),yg(1:NELT,1:kkmax)
        real(4)::strsave(1:NELT,1:kkmax,1:3)

        integer::ne,nc,kk,kkk
        real(4)::sigx,sigy,tauxy,ps1,ps2,ang
        real(4)::fact0=1.0E-10              !0判定

        if(nod==3)kkk=1
        if(nod==4.and.IPR==1)kkk=1
        if(nod==4.and.IPR==0)kkk=4

        open (12, file='_ps1.txt', status='replace')
        open (13, file='_ps2.txt', status='replace')
        if(ncut==0)then
            do ne=1,NELT
                do kk=1,kkk
                    sigx =strsave(ne,kk,1)
                    sigy =strsave(ne,kk,2)
                    tauxy=strsave(ne,kk,3)
                    call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                    write(12,*) xg(ne,kk),' ',yg(ne,kk),' ',ps1
                    write(13,*) xg(ne,kk),' ',yg(ne,kk),' ',ps2
                end do
            end do
        else
            do ne=1,NELT
                do nc=1,ncut
                    if(matno(ne)/=numcut(nc))then
                        do kk=1,kkk
                            sigx =strsave(ne,kk,1)
                            sigy =strsave(ne,kk,2)
                            tauxy=strsave(ne,kk,3)
                            call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
                            write(12,*) xg(ne,kk),' ',yg(ne,kk),' ',ps1
                            write(13,*) xg(ne,kk),' ',yg(ne,kk),' ',ps2
                        end do
                    end if
                end do
            end do
        end if
        close(12)
        close(13)

        if(0<ncut)then
            open (12, file='_cutm.txt', status='replace')
            select case(nod)
                case(3)
                    do ne=1,NELT
                        do nc=1,ncut
                            if(matno(ne)==numcut(nc))then
                                i=kakom(ne,1)
                                j=kakom(ne,2)
                                k=kakom(ne,3)
                                write(12,*) '>'
                                write(12,*) (x(i)),' ',(y(i))
                                write(12,*) (x(j)),' ',(y(j))
                                write(12,*) (x(k)),' ',(y(k))
                            end if
                        end do
                    end do
                case(4)
                    do ne=1,NELT
                        do nc=1,ncut
                            if(matno(ne)==numcut(nc))then
                                i=kakom(ne,1)
                                j=kakom(ne,2)
                                k=kakom(ne,3)
                                l=kakom(ne,4)
                                write(12,*) '>'
                                write(12,*) (x(i)),' ',(y(i))
                                write(12,*) (x(j)),' ',(y(j))
                                write(12,*) (x(k)),' ',(y(k))
                                write(12,*) (x(l)),' ',(y(l))
                            end if
                        end do
                    end do
            end select
            close (12)
        end if
    end subroutine DAT_PS


    subroutine MAKE_BAT(MATEL,KOX,KOY,KXY)
        integer::MATEL,KOX,KOY,KXY

        integer::ma
        character::srange*100,sscale*100,sba*100,sgrdspc*100,ssradi*100,filename*100
        character::ssx*100,ssy*100,scc*50,sca*50,sscld*100,scol*50,snum*50,spsize*50,scsize*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::psize,ds
        real(4)::work
        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=bav

        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)') ymin-ds;
        write(strw4,'(f10.3)') ymax;
        srange='set range='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
        if(ymax-ymin<=xmax-xmin)then
            write(strw1,'(f10.3)') baselength;
            write(strw2,'(f10.3)') baselength/(xmax-xmin+ds)*(ymax-ymin+ds)
        else
            write(strw1,'(f10.3)') baselength/(ymax-ymin+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))//'"'
        ssy='set stry="'//trim(adjustl(stry))//'"'
        if(ymax-ymin<=xmax-xmin)then
            work=baselength/(xmax-xmin)*(ymax-ymin)
            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/(ymax-ymin+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_plnt.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(ssy))
        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%:%stry%: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<KOX)write(12,'(a)') 'psxy _kox.txt -R -J -Sd%psize% -W3/black -N -K -O >> %fig%'
            if(0<KOY)write(12,'(a)') 'psxy _koy.txt -R -J -St%psize% -W3/black -N -K -O >> %fig%'
            if(0<KXY)write(12,'(a)') 'psxy _kxy.txt -R -J -SC%psize% -W3/black -N -K -O >> %fig%'
        end if
        if(0<ifc.and.fact0<fmax)write(12,'(a)') 'psxy _force.txt -R -J -Svs0.02/0.2/0.1 -G0/0/255 -N -K -O >> %fig%'
        write(12,'(a)') 'echo 0 0 | psxy -R -J -Sp -O >> %fig%'
        write(12,'()')
        !drawing of contour of ps1
        write(12,'(a)') 'rem *** Drawing of contour of ps1 ***'
        write(12,'(a)') 'set fig=fig_gmt_ps1.eps'
        write(12,'(a)') 'grdmask _mask_data.txt -G_mask.grd -I%grdspc% -R%range% -S%sradi% -NNan/0/0'
        write(12,'(a)') 'surface _ps1.txt -G_ps1.grd -I%grdspc% -R -T0.0'
        write(12,'(a)') 'grdmath _ps1.grd _mask.grd OR = _ps1m.grd'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%stry%:WS -P -K > %fig%'
        write(12,'(a)') 'grdimage _ps1m.grd -J -C0_cpt.cpt -K -O >> %fig%'
        if(0<ncut)write(12,'(a)') 'psxy _cutm.txt -R -J -G0 -K -O >> %fig%'
!        write(12,'(a)') 'grdcontour _ps1m.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 contour of ps2
        write(12,'(a)') 'rem *** Drawing of contour of ps1 ***'
        write(12,'(a)') 'set fig=fig_gmt_ps2.eps'
        write(12,'(a)') 'surface _ps2.txt -G_ps2.grd -I%grdspc% -R -T0.0'
        write(12,'(a)') 'grdmath _ps2.grd _mask.grd OR = _ps2m.grd'
        write(12,'(a)') 'psbasemap -R%range% -JX%scale% -B%ba%:%strx%:/%ba%:%stry%:WS -P -K > %fig%'
        write(12,'(a)') 'grdimage _ps2m.grd -J -C0_cpt.cpt -K -O >> %fig%'
        if(0<ncut)write(12,'(a)') 'psxy _cutm.txt -R -J -G0 -K -O >> %fig%'
!        write(12,'(a)') 'grdcontour _ps2m.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,'()')
        !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 PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        !主応力計算
        real(4),intent(in)::sigx,sigy,tauxy
        real(4),intent(out)::ps1,ps2,ang

        ps1=0.5*(sigx+sigy)+sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        ps2=0.5*(sigx+sigy)-sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        if(sigx==sigy)then
            if(tauxy>0.0)ang=45.0
            if(tauxy<0.0)ang=135.0
            if(tauxy==0.0)ang=0.0
        else
            ang=0.5*atan(2.0*tauxy/(sigx-sigy))
            ang=180.0/pi*ang
            if(sigx>sigy.and.tauxy>=0.0)ang=ang
            if(sigx>sigy.and.tauxy<0.0)ang=ang+180.0
            if(sigx<sigy)ang=ang+90.0
        end if
    end subroutine PST_PL


    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 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_GMTplnt