program f90_num_plnt
    !****************************************************************
    !Optimization of node numbers for 2D stress analysis
    !****************************************************************
    implicit none
    integer::i,j,k,l,ne
    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::NF                         !Number of loaded nodes
    integer::nt                         !Total degrees of freedom
    integer::nod                        !Number of node of 1 element
    integer::nhen=2                     !Degrees of freedom of 1 node
    integer,allocatable::kakom(:,:)     !Element connectivity
    integer,allocatable::matno(:)       !Material set number
    real(8),allocatable::wm(:,:)        !work for material properties
    real(8),allocatable::x(:)           !x-coordinate of element
    real(8),allocatable::y(:)           !y-coordinate of element
    real(8),allocatable::deltaT(:)      !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
    real(8),allocatable::f(:)           !External load vector
    real(8),allocatable::rdis(:)        !Displacements of restricted node
    integer,allocatable::nfd(:)         !Node number of loaded node
    integer,allocatable::numr(:)        !Original node number
    integer,allocatable::jnd(:)         !Revised node number

    integer::IPR                        !Flag for outout of stress
    character :: linebuf*1000

    character::fnameR*50,fnameW*50      !Input file name, output file name

    call getarg(1,fnameR) !Input file name
    call getarg(2,fnameW) !Output file name

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:8))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(deltaT(1:NODT))
    allocate(nokx(1:KOX))
    allocate(noky(1:KOY))
    allocate(f(1:nt))
    allocate(rdis(1:nt))
    allocate(nfd(1:NF))
    allocate(numr(1:NODT))
    allocate(jnd(1:NODT))
    !材料特性（t,Em,po,gamma,gkh,gkv,alpha,ts）
    do i=1,MATEL
        read(11,*) (wm(i,j),j=1,8)
    end do
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,nod),matno(ne)
    end do
    do i=1,NODT
        read(11,*) x(i),y(i),deltaT(i)
    end do
    do i=1,nt
        rdis(i)=0.0D0
    end do
    if(0<KOX)then
        do i=1,KOX
            read(11,*) nokx(i),rdis(2*nokx(i)-1)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            read(11,*) noky(i),rdis(2*noky(i))
        end do
    end if
    do i=1,nt
        f(i)=0.0D0
    end do
    if(0<NF)then
        do i=1,NF
            read(11,*) nfd(i),f(2*nfd(i)-1),f(2*nfd(i))
        end do
    end if
    close(11)

    !*****************************************************
    !Renumbering
    ! To save revised node number in array 'jnd(i)' where 'i' is original node number.
    ! To save original node number in array 'numr(i)' where 'i' is revised node number.
    !*****************************************************
    do i=1,NODT
        jnd(i)=i
    end do

    call opt(kakom,jnd,NODT,NELT,nod)

    do i=1,NODT
        do j=1,NODT
            if(i==jnd(j))exit
        end do
        numr(i)=j
    end do

    !*****************************************************
    !To print out
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(linebuf,*) nod,',',NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',NF,',',NSTRES,',',IPR
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    do i=1,MATEL
        write(linebuf,*) wm(i,1),',',wm(i,2),',',wm(i,3),',',wm(i,4),',',wm(i,5),',',wm(i,6),',',wm(i,7),',',wm(i,8)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    select case(nod)
        case(3)
            do ne=1,NELT
                i=kakom(ne,1)
                j=kakom(ne,2)
                k=kakom(ne,3)
                kakom(ne,1)=jnd(i)
                kakom(ne,2)=jnd(j)
                kakom(ne,3)=jnd(k)
                write(linebuf,*) kakom(ne,1),',',kakom(ne,2),',',kakom(ne,3),',',matno(ne)
                call del_spaces(linebuf)
                write (12,'(a)') trim(linebuf)
            end do
        case(4)
            do ne=1,NELT
                i=kakom(ne,1)
                j=kakom(ne,2)
                k=kakom(ne,3)
                l=kakom(ne,4)
                kakom(ne,1)=jnd(i)
                kakom(ne,2)=jnd(j)
                kakom(ne,3)=jnd(k)
                kakom(ne,4)=jnd(l)
                write(linebuf,*) kakom(ne,1),',',kakom(ne,2),',',kakom(ne,3),',',kakom(ne,4),',',matno(ne)
                call del_spaces(linebuf)
                write (12,'(a)') trim(linebuf)
            end do
    end select
    do i=1,NODT
        write(linebuf,*) x(numr(i)),',',y(numr(i)),',',deltaT(numr(i))
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    if(0<KOX)then
        do i=1,KOX
            write(linebuf,*) jnd(nokx(i)),',',rdis(2*nokx(i)-1)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            write(linebuf,*) jnd(noky(i)),',',rdis(2*noky(i))
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(0<NF)then
        do i=1,NF
            write(linebuf,*) jnd(nfd(i)),',',f(2*nfd(i)-1),',',f(2*nfd(i))
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    close(12)

stop

contains

    subroutine opt(kakom,jnd,NODT,NELT,nod)
        integer::NODT,NELT,nod
        integer::kakom(1:NELT,1:nod),jnd(1:NODT)

        integer::i,j,k,l,ne

        open(12, file='fnameT1.txt', status='replace')
        write(linebuf,*) NODT,',',NELT,',',nod
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
        select case(nod)
            case(3)
                do ne=1,NELT
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    write(linebuf,*) kakom(ne,1),',',kakom(ne,2),',',kakom(ne,3)
                    call del_spaces(linebuf)
                    write (12,'(a)') trim(linebuf)
                end do
            case(4)
                do ne=1,NELT
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    l=kakom(ne,3)
                    write(linebuf,*) kakom(ne,1),',',kakom(ne,2),',',kakom(ne,3),',',kakom(ne,4)
                    call del_spaces(linebuf)
                    write (12,'(a)') trim(linebuf)
                end do
        end select
        close(12)

        call system('./f90_num')

        open(11,file='fnameT2.txt',status='old')
        do i=1,NODT
            read(11,*) k,jnd(k)
        end do
        close(11)

        call system('del fnameT1.txt')
        call system('del fnameT2.txt')

    end subroutine opt


    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_num_plnt
