program f90_fig_mode
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    integer::i,j,ne,kk
    character::strcom*50                !
    integer::NODT                       !
    integer::NELT                       !
    integer::MATEL                      !
    integer::KOX                        !
    integer::KOY                        !
    integer::KOZ                        !
    integer::NF                         !
    integer,allocatable::kakom(:,:)     !
    real(8),allocatable::x(:)           !
    real(8),allocatable::y(:)           !
    integer,allocatable::nokx(:)        !
    integer,allocatable::noky(:)        !
    integer,allocatable::nokz(:)        !
    real(8),allocatable::ff(:)          !
    real(8),allocatable::dist(:)        !
    integer::nt                         !
    integer::nod=2                      !
    integer::nhen=3                     !
    character::fnameR*100,fnameW*100,fnameS*100
    integer::nlmax,itmax,nload
    integer::itera
    real(8)::lam
    real(8)::deltaS
    real(8)::Tpress

    integer::nstep
    integer,allocatable::istep(:)
    real(8),allocatable::Tpr(:)
    real(8),allocatable::fx(:,:)
    real(8),allocatable::fy(:,:)
    real(8),allocatable::disx(:,:)
    real(8),allocatable::disy(:,:)

    character :: linebuf*1000,dummy*100,filename*100
    integer::iw,k1,k2,k3,n1,n2,n3
    real(8)::rw
    real(8)::coef,ddx,ddy,xmin,xmax,ymin,ymax,ddh,fmax,xx,yy

    !Input from command line
    call getarg(1,fnameR)
    call getarg(2,fnameS)
    call getarg(3,dummy)
    read(dummy,*) nstep
    allocate(istep(1:nstep))
    do j=1,nstep
        call getarg(3+j,dummy)
        read(dummy,*) istep(j)
    end do

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,'()')
    read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ,NF,nlmax,itmax,rw
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    if(1<=KOX)allocate(nokx(1:KOX))
    if(1<=KOY)allocate(noky(1:KOY))
    if(1<=KOZ)allocate(nokz(1:KOZ))
    allocate(ff(1:nt))
    allocate(dist(1:nt))
    allocate(fx(1:NODT,1:nstep))
    allocate(fy(1:NODT,1:nstep))
    allocate(disx(1:NODT,1:nstep))
    allocate(disy(1:NODT,1:nstep))
    allocate(Tpr(1:nstep))

    read(11,'()')
    read(11,'()')
    n1=0
    n2=0
    n3=0
    do i=1,NODT
        read(11,*) iw,x(i),y(i),rw,rw,rw,k1,k2,k3
        if(0<KOX)then
            if(k1==1)then
                n1=n1+1
                nokx(n1)=i
            end if
        end if
        if(0<KOY)then
            if(k2==1)then
                n2=n2+1
                noky(n2)=i
            end if
        end if
        if(0<KOZ)then
            if(k3==1)then
                n3=n3+1
                nokz(n3)=i
            end if
        end if
    end do
    read(11,'()')
    read(11,'()')
    do ne=1,NELT
        read(11,*) iw,(kakom(ne,j),j=1,nod),rw,rw,rw,iw
    end do
    do nload=1,nlmax
        read(11,'(a)') linebuf
        if(linebuf(1:1)/='*')exit
        read(11,'()')
        read(11,*) iw,itera,Tpress,lam,deltaS,rw
        if(itmax<itera)exit
        read(11,'()')
        read(11,'()')
        do i=1,NODT
            read(11,*) iw,rw,rw,dist(3*i-2),dist(3*i-1),dist(3*i),ff(3*i-2),ff(3*i-1),ff(3*i),rw,rw,rw
        end do
        read(11,'()')
        read(11,'()')
        do ne=1,NELT
            read(11,'()')
        end do
        do j=1,nstep
            if(nload==istep(j))then
                Tpr(j)=Tpress
                do i=1,NODT
                    disx(i,j)=dist(3*i-2)
                    disy(i,j)=dist(3*i-1)
                    fx(i,j)=ff(3*i-2)
                    fy(i,j)=ff(3*i-1)
                end do
            end if
        end do
    end do
    close(11)


    do kk=1,nstep
        xmin=x(1)
        xmax=x(1)
        ymin=y(1)
        ymax=y(1)
        do i=1,NODT
            if(x(i)<xmin)xmin=x(i)
            if(y(i)<ymin)ymin=y(i)
            if(xmax<x(i))xmax=x(i)
            if(ymax<y(i))ymax=y(i)
        end do
        ddx=xmax-xmin
        ddy=ymax-ymin
        if(ddx<=ddy)coef=0.20*ddy
        if(ddy<=ddx)coef=0.20*ddx

        write (filename, '("_dis", i3.3, ".txt")') kk
        open(12, file=filename, status='replace')
        do ne=1,NELT
            i=kakom(ne,1)
            j=kakom(ne,2)
            write(12,'(a)') '>'
            write(12,*) x(i)+disx(i,kk),y(i)+disy(i,kk)
            write(12,*) x(j)+disx(j,kk),y(j)+disy(j,kk)
        end do
        close(12)

        if(ddx<=ddy)ddh=0.05*ddy
        if(ddy<=ddx)ddh=0.05*ddx
        write (filename, '("_bou", i3.3, ".txt")') kk
        open(12, file=filename, status='replace')
        if(0<KOX)then
            do i=1,KOX
                xx=x(nokx(i))+disx(nokx(i),kk)
                yy=y(nokx(i))+disy(nokx(i),kk)
                write(12,'(a)') '>'
                write(12,*) xx,yy
                write(12,*) xx-ddh*cos(30.0D0/180.0D0*pi),yy+0.5*ddh
                write(12,*) xx-ddh*cos(30.0D0/180.0D0*pi),yy-0.5*ddh
                write(12,*) xx,yy
                write(12,'(a)') '>'
                write(12,*) xx-1.5D0*ddh*cos(30.0D0/180.0D0*pi),yy+0.75*ddh
                write(12,*) xx-1.5D0*ddh*cos(30.0D0/180.0D0*pi),yy-0.75*ddh
            end do
        end if
        if(0<KOY)then
            do i=1,KOY
                xx=x(noky(i))+disx(noky(i),kk)
                yy=y(noky(i))+disy(noky(i),kk)
                write(12,'(a)') '>'
                write(12,*) xx,yy
                write(12,*) xx+0.5*ddh,yy-ddh*cos(30.0D0/180.0D0*pi)
                write(12,*) xx-0.5*ddh,yy-ddh*cos(30.0D0/180.0D0*pi)
                write(12,*) xx,yy
                write(12,'(a)') '>'
                write(12,*) xx-0.75*ddh,yy-1.5D0*ddh*cos(30.0D0/180.0D0*pi)
                write(12,*) xx+0.75*ddh,yy-1.5D0*ddh*cos(30.0D0/180.0D0*pi)
            end do
        end if
        close(12)

        fmax=0.0
        do i=1,NODT
            if(fmax<abs(fx(i,kk)))fmax=abs(fx(i,kk))
            if(fmax<abs(fy(i,kk)))fmax=abs(fy(i,kk))
        end do
        if(ddx<=ddy)coef=0.1*ddy
        if(ddy<=ddx)coef=0.1*ddx
        write (filename, '("_load", i3.3, ".txt")') kk
        open(12, file=filename, status='replace')
        do i=1,NODT
            write(12,*) x(i)+disx(i,kk)-fx(i,kk)/fmax*coef,y(i)+disy(i,kk)-fy(i,kk)/fmax*coef,&
                        x(i)+disx(i,kk),y(i)+disy(i,kk)
        end do
        close(12)
    end do

    open(12, file='_elm.txt', status='replace')
    do ne=1,NELT
        i=kakom(ne,1)
        j=kakom(ne,2)
        write(12,'(a)') '>'
        write(12,*) x(i),y(i)
        write(12,*) x(j),y(j)
    end do
    close(12)

    call GMTBAT(nstep,Tpr,fnameS,fnameW)
    dummy='bash '//fnameW
    call system(trim(adjustl(dummy)))

stop

contains

    subroutine GMTBAT(nstep,Tpr,fnameS,fnameW)
        integer,intent(in)::nstep
        real(8),intent(in)::Tpr(1:nstep)
        character,intent(in)::fnameS*100
        character,intent(out)::fnameW*100
        integer::kk
        character::linebuf*1000,filename*100,dummy*50

        linebuf='_gmt_'//fnameS//'_mode.txt'
        call del_spaces(linebuf)
        write(fnameW,'(a)') trim(linebuf)
        open(12, file=fnameW, status='replace')
        write(12,'(a)') 'gmt set MAP_FRAME_PEN 0'
        write(12,'(a)') 'gmt set MAP_VECTOR_SHAPE 0'
        write(12,'()')

        do kk=1,nstep
            write(dummy, '(i3.3)') kk
            filename='fig_'//trim(fnameS)//dummy//'.eps'
            call del_spaces(filename)
            linebuf='fig='//trim(filename)
            write(12,'(a)') trim(linebuf)
            write(12,'(a)') 'gmt psbasemap -R-2500/2500/-2500/2500 -JX10/10 -Bxg2500 -Byg2500 -P -K > $fig'
            write(12,'(a)') 'gmt psxy _elm.txt -R -J -W2,5_2:0  -O -K >> $fig'
            write(filename, '("_dis", i3.3, ".txt")') kk
            linebuf='gmt psxy '//trim(filename)//' -R -J -W2  -O -K >> $fig'
            write(12,'(a)') trim(linebuf)
            write(filename, '("_bou", i3.3, ".txt")') kk
            linebuf='gmt psxy '//trim(filename)//' -R -J -W1 -O -K >> $fig'
            write(12,'(a)') trim(linebuf)
            write(filename, '("_load", i3.3, ".txt")') kk
            linebuf='gmt psxy '//trim(filename)//' -R -J -Sv0.2+s+e+a60  -W3,gray -O -K >> $fig'
            write(12,'(a)') trim(linebuf)
            write(dummy,'(E10.3)') Tpr(kk)
            if(kk== 1)linebuf="P="//trim(dummy)
            if(kk== 2)linebuf="P="//trim(dummy)
            if(kk== 3)linebuf="P="//trim(dummy)
            if(kk== 4)linebuf="P="//trim(dummy)
            if(kk== 5)linebuf="P="//trim(dummy)
            if(kk== 6)linebuf="P="//trim(dummy)
            if(kk== 7)linebuf="P="//trim(dummy)
            if(kk== 8)linebuf="P="//trim(dummy)
            if(kk== 9)linebuf="P="//trim(dummy)
            write(12,'(a)') 'echo 0.02 0.98 10 0 TL '//trim(linebuf)//'| gmt pstext -R0/1/0/1 -JX -F+f+a+j -O >> $fig'
            write(12,'()')
        end do
        close(12)
    end subroutine GMTBAT


    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_fig_mode
