module cdraw_data
    implicit none
    real(4)::xmin,xmax,ymin,ymax,zmin,zmax
    integer::ind,ibc
    real(4)::baselength,baseheight
    real(4)::bavx,bfvx
    real(4)::bavz,bfvz
    character::strx*100,stry*100
end module cdraw_data

program f90_gmt_heat
    !Drawing by GMT
    !2-dimensional unsteady heat conduction analysis
    use cdraw_data
    implicit none
    integer::i,j,n,ne

    character(len=100)::strcom       !Comment
    integer::NODT                    !Number of nodes
    integer::NELT                    !Number of elements
    integer::MATEL                   !Number of material sets
    integer::KOT                     !Number of nodes with given temperature
    integer::KOC                     !Number of sides with the condition of heat transfer boundary
    integer::nod=4                   !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::y(:)        !y-coordinate of node
    integer,allocatable::nokt(:)     !Node number with given temperature
    integer,allocatable::nekc(:,:)   !Element number and side number with the condition of heat transfer boundary
    integer::kchen                   !work for side number with heat transfer boundary
    real(4)::delta                   !Time increment
    integer::itmax                   !Number of time history steps
    integer::nout                    !Number of of nodes for all time history steps
    integer::ntout                   !Frequency of output of all node's temperatures
    real(4),allocatable::ttime(:)    !Time
    integer,allocatable::noout(:)    !Node number for output of temperature of all time history steps
    integer,allocatable::nntout(:)   !Step number for all node's temperarture output
    real(4),allocatable::tempe_n(:,:)!Nodal temperature for output of all time history steps
    real(4),allocatable::tempe_s(:,:)!Nodal temperature for output of all node's temperatures
    integer::io
    real(4)::rw,alp(1:4)
    integer::iw
    character,allocatable::str(:)*20
    character::strdum*100,str1*20
    real(4)::fact0=1e-10

    character::linebuf*1000
    character::fnameR*50

    call getarg(1,fnameR) !Input fule name

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,'()')
    read(11,*) NODT,NELT,MATEL,KOT,KOC,delta,nout,ntout,itmax
    allocate(kakom(1:NELT,1:4))
    allocate(matno(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(nokt(1:KOT))
    allocate(nekc(1:KOC,1:2))
    allocate(noout(1:nout))
    allocate(nntout(1:ntout+1))
    allocate(ttime(1:itmax+1))
    allocate(tempe_n(1:itmax+1,1:nout))
    allocate(tempe_s(1:NODT,1:ntout+1))
    if(ntout<=nout)allocate(str(1:nout))
    if(nout<=ntout)allocate(str(1:ntout))
    read(11,'()')
    read(11,'()')
    j=0
    do i=1,NODT
        read(11,*) iw,x(i),y(i),n,rw
        if(n==1)then
            j=j+1
            nokt(j)=i
        end if
    end do
    read(11,'()')
    read(11,'()')
    i=0
    do ne=1,NELT
        read(11,*) iw,(kakom(ne,j),j=1,nod),(alp(j),j=1,4),rw,rw,rw,rw,rw,matno(ne)
        do j=1,4
            if(fact0<alp(j))then
                i=i+1
                nekc(i,1)=ne
                nekc(i,2)=j
            end if
        end do
    end do
    read(11,'()')
    read(11,*) strdum,strdum,(str(i),i=1,nout)
    str1='node-'
    do i=1,nout
        noout(i)=TAKENUM(str(i),str1)
    end do
    do i=1,itmax+1
        read(11,*) iw,ttime(i),(tempe_n(i,j),j=1,nout) ! i=1:initial temperature
    end do
    read(11,'()')
    if(0<ntout)then
        read(11,*) strdum,strdum,(str(i),i=1,ntout)
        str1='step-'
        nntout(1)=0
        do i=2,ntout+1
            nntout(i)=TAKENUM(str(i-1),str1)
        end do
        do i=1,NODT
            read(11,*) iw,(tempe_s(i,j),j=1,ntout+1) ! j=1:initial temperature
        end do
    end if
    close(11)

    open(11,file='_inp_base.txt',status='old')
    read(11,*) ind,ibc
    read(11,*) baselength,baseheight
    read(11,*) bavx,bfvx
    read(11,*) bavz,bfvz
    read(11,*) strx
    read(11,*) stry
    close(11)

    xmin=x(1)
    xmax=x(1)
    ymin=y(1)
    ymax=y(1)
    zmin=tempe_s(1,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))
        do j=1,ntout+1
            if(tempe_s(i,j)<zmin)zmin=tempe_s(i,j)
            if(zmax<tempe_s(i,j))zmax=tempe_s(i,j)
        end do
    end do

    call DAT_MESH(kakom,matno,x,y,nokt,nekc,nod,NODT,NELT,MATEL,KOT,KOC)
    if(0<ind)call DAT_IND(kakom,x,y,NELT,NODT,nod)
    if(0<ntout)call DAT_XYZ(kakom,x,y,tempe_s,NELT,NODT,nod,ntout)
    call MAKE_BAT(MATEL,KOT,KOC,delta,nntout,ntout)
    call system('bash _gmt_xyz.txt')

stop

Contains

    integer function TAKENUM(str0,str1)
        character,intent(in)::str0*20,str1*20

        integer::n1,n2,num
        character::str2*20

        n1=len_trim(str1)+1
        n2=len_trim(str0)
        str2=str0(n1:n2)
        read(str2,*) num
        TAKENUM=num
    end function TAKENUM


    subroutine DAT_MESH(kakom,matno,x,y,nokt,nekc,nod,NODT,NELT,MATEL,KOT,KOC)
        integer::nod,NODT,NELT,MATEL,KOT,KOC
        integer::kakom(1:NELT,1:nod),matno(1:NELT)
        integer::nokt(1:KOT),nekc(1:KOC,1:2)
        real(4)::x(1:NODT),y(1:NODT)

        integer::i,j,k,l,ma,ne,ia,ja
        character::filename*50

        do ma=1,MATEL
            write (filename, '("_mesh", i2.2, ".txt")') ma
            open (12, file=filename, status='replace')
            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)),' ',(y(i))
                    write(12,*) (x(j)),' ',(y(j))
                    write(12,*) (x(k)),' ',(y(k))
                    write(12,*) (x(l)),' ',(y(l))
                end if
            end do
            close (12)
        end do

        if(0<KOT)then
            open (12, file='_kot.txt', status='replace')
            do i=1,KOT
                write(12,*) x(nokt(i)),' ',y(nokt(i))
            end do
            close(12)
        end if
        if(0<KOC)then
            open (12, file='_koc.txt', status='replace')
            do k=1,KOC
                ne=nekc(k,1)
                ia=nekc(k,2)
                ja=nekc(k,2)+1
                if(ia==4)ja=1
                write(12,*) x(kakom(ne,ia)),' ',y(kakom(ne,ia))
                write(12,*) x(kakom(ne,ja)),' ',y(kakom(ne,ja))
            end do
            close(12)
        end if
    end subroutine DAT_MESH


    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')
        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,' 9p 0 MC ',ne
        end do
        close(12)
        close(13)

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


    subroutine DAT_XYZ(kakom,x,y,tempe_s,NELT,NODT,nod,ntout)
        integer,intent(in)::NELT,NODT,nod,ntout
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(4),intent(in)::x(1:NODT),y(1:NODT),tempe_s(1:NODT,1:ntout+1)

        integer::ne,i,j,k,l,n
        character::filename*50

        do n=1,ntout+1
            write (filename,'("_inp_xyz_",i2.2,".txt")') n
            open(13,file=filename,status='replace')
            do ne=1,NELT
                i=kakom(ne,1)
                j=kakom(ne,2)
                k=kakom(ne,3)
                l=kakom(ne,4)
                write(13,'(a)') '>'
                write(13,*) x(i),y(i),tempe_s(i,n)
                write(13,*) x(j),y(j),tempe_s(j,n)
                write(13,*) x(k),y(k),tempe_s(k,n)
                write(13,*) x(l),y(l),tempe_s(l,n)
                write(13,*) x(i),y(i),tempe_s(i,n)
            end do
            close(13)
        end do
    end subroutine DAT_XYZ


    subroutine MAKE_BAT(MATEL,KOT,KOC,delta,nntout,ntout)
        integer,intent(in)::MATEL,KOT,KOC,ntout
        integer,intent(in)::nntout(1:ntout+1)
        real(4),intent(in)::delta

        integer::i,n,ma
        character::srangex*100,srangez*100,sscalex*100,sscalez*100,sstep*50,filename*50
        character::sbx*50,sbz*50,ssx*100,ssy*100,ssz*100,sinpf*50
        character::scol*50,snum*50,spsize*50,scsize*50
        character::strw1*50,strw2*50,strw3*50,strw4*50,strw5*50,strw6*50
        real(4)::psize,ds,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=bfvx

        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;
        srangex='rangex='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))&
                        //'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
        write(strw1,'(f10.3)') xmin;
        write(strw2,'(f10.3)') xmax;
        write(strw3,'(f10.3)') ymin;
        write(strw4,'(f10.3)') ymax;
        write(strw5,'(f10.3)') zmin;
        write(strw6,'(f10.3)') zmax;
        srangez='rangez='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))&
                        //'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))&
                        //'/'//trim(adjustl(strw5))//'/'//trim(adjustl(strw6))
        write(strw1,'(f10.3)') baselength/(xmax-xmin)*(xmax-xmin)
        write(strw2,'(f10.3)') baselength/(xmax-xmin)*(ymax-ymin)
        write(strw3,'(f10.3)') baseheight
        sscalex='scalex='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
        sscalez='scalez='//trim(adjustl(strw3))
        write(strw1,'(f10.3)') bavx
        write(strw2,'(f10.3)') bfvx
        sbx='bax=a'//trim(adjustl(strw1))//'f'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') bavz
        write(strw2,'(f10.3)') bfvz
        sbz='baz=a'//trim(adjustl(strw1))//'f'//trim(adjustl(strw2))
        ssx='strx="'//trim(adjustl(strx))//'"'
        ssy='stry="'//trim(adjustl(stry))//'"'
        ssz='strz="Temperature (@%12%\260@%%C)"'
        write(strw1,'(f10.3)') psize
        spsize='psize='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') psize*4.0
        scsize='csize='//trim(adjustl(strw1))

        open (12, file='_gmt_xyz.txt', status='replace')
        write(12,'(a)') 'gmt set FONT_ANNOT_PRIMARY 12'
        write(12,'(a)') 'gmt set MAP_ANNOT_OFFSET_PRIMARY 0.5c'
        write(12,'(a)') 'gmt set FONT_LABEL 12p'
        write(12,'(a)') 'gmt set MAP_LABEL_OFFSET 0.5c'
        write(12,'(a)') trim(adjustl(srangex))
        write(12,'(a)') trim(adjustl(srangez))
        write(12,'(a)') trim(adjustl(sscalex))
        write(12,'(a)') trim(adjustl(sscalez))
        write(12,'(a)') trim(adjustl(sbx))
        write(12,'(a)') trim(adjustl(sbz))
        write(12,'(a)') trim(adjustl(ssx))
        write(12,'(a)') trim(adjustl(ssy))
        write(12,'(a)') trim(adjustl(ssz))
        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$rangex -JX$scalex -Bx$bax+l"$strx" -By$bax+l"$stry" -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<KOT)write(12,'(a)') 'gmt psxy _kot.txt -R -J -SC$psize -W2,black -N -K -O >> $fig'
            if(0<KOC)write(12,'(a)') 'gmt psxy _koc.txt -R -J -St$psize -W2,black -N -K -O >> $fig'
        end if
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        if(0<ntout)then
            do n=1,ntout+1
                write(12,'()')
                write (filename,'("fig=_fig_gmt_xyz_",i2.2,".eps")') n
                write (sstep,'("step=",f6.1,"hr")') real(nntout(n))*delta
                call del_spaces(sstep)
                write (sinpf,'("inpf=_inp_xyz_",i2.2,".txt")') n
                write(12,'(a)') trim(adjustl(filename))
                write(12,'(a)') trim(adjustl(sstep))
                write(12,'(a)') trim(adjustl(sinpf))
                write(12,'(a)') 'gmt psbasemap -R$rangez -JX$scalex -JZ$scalez &
                        -Bx$bax+l"$strx" -By$bax+l"$stry" -Bz$baz+l"$strz" -BESZnw -p150/30 -P -K > $fig'
                write(12,'(a)') 'gmt psxyz $inpf -JX -JZ -R -W1 -p150/30 -K -O >> $fig'
                write(12,'(a)') 'echo 0.90 0.90 14 0 TL t=$step | gmt pstext -R0/1/0/1 -JX -F+f+a+j -N -O >> $fig'
            end do
        end if

        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_heat
