program f90_num_heat
    !****************************************************************
    !Optimization of node numbers for 2D heat conduction analysis
    !****************************************************************
    implicit none
    integer::i,j,k,l,ne

    character(len=50)::ftinp         !File name for temperature time history
    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 temerature
    integer::KOC                     !Number of sides with the condition of heat transfer boundary
    integer::nod=4                   !Number of nodes for 1 element
    integer::nhen=1                  !Degree of freedom of 1 node
    integer::nt                      !Total degrees of freedom (NODT x nhen)
    integer::mm                      !Degrees of freedom of reduced FE equation
    integer::ib                      !Band width of reduced FE equation
    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 node
    real(8),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
    real(8),allocatable::alphac(:)   !Heat transfer rate of specified side
    real(8),allocatable::tempe0(:)   !Initial value of nodal temperature
    real(8)::delta                   !Time increment
    real(8)::ttime                   !Time
    integer::itmax                   !Number of time steps
    integer::nout                    !Number of nodes for output of temperature time history
    integer,allocatable::noout(:)    !Node number for output of temperature time history
    integer::ntout                   !Frequency of output of all nodal temperature
    integer,allocatable::nntout(:)   !Step number for output of all nodal temperature
    real(8),allocatable::tempend(:,:)!Nodal temperature for output

    integer,allocatable::numr(:)    !Original node numbers
    integer,allocatable::jnd(:)     !Revised node numbers
    character::linebuf*1000
    character::str*20

    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

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,'(a)') ftinp
    read(11,*) NODT,NELT,MATEL,KOT,KOC,delta
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:MATEL,1:5))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(nokt(1:KOT))
    allocate(nekc(1:KOC,1:2))
    allocate(alphac(1:KOC))
    allocate(tempe0(1:nt))
    allocate(numr(1:NODT))
    allocate(jnd(1:NODT))
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3),wm(i,4),wm(i,5)
    end do
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),kakom(ne,3),kakom(ne,4),matno(ne)
    end do
    do i=1,NODT
        read(11,*) x(i),y(i),tempe0(i)
    end do
    if(0<KOT)then
        do i=1,KOT
            read(11,*) nokt(i)
        end do
    end if
    if(0<KOC)then
        do i=1,KOC
            read(11,*) nekc(i,1),nekc(i,2),alphac(i)
        end do
    end if
    read(11,*) nout
    allocate(noout(1:nout))
    read(11,*) (noout(i),i=1,nout)
    read(11,*) ntout
    if(0<ntout)then
        allocate(nntout(1:ntout))
        allocate(tempend(1:NODT,1:ntout))
        read(11,*) (nntout(i),i=1,ntout)
    end if
    close(11)

    !*****************************************************
    !Renumbering
    ! To save revised node numbers in array 'jnd(i)' where 'i' is original node number.
    ! To save original node numbers in array 'numr(i)' where 'i' is revised node number.
    !*****************************************************
    call optnum(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(12,'(a)') trim(adjustl(ftinp))
    write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOT,',',KOC,',',delta
    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)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    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
    do i=1,NODT
        write(linebuf,*) x(numr(i)),',',y(numr(i)),',',tempe0(numr(i))
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    if(0<KOT)then
        do i=1,KOT
            write(linebuf,*) jnd(nokt(i))
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(0<KOC)then
        do i=1,KOC
            write(linebuf,*) nekc(i,1),',',jnd(nekc(i,2)),',',alphac(i)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    write(linebuf,*) nout
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    if(0<nout)then
        write(linebuf,*) jnd(noout(1))
        do i=2,nout
            write(str,*) jnd(noout(i))
            linebuf=trim(linebuf)//','//str
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end if
    write(linebuf,*) ntout
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    if(0<ntout)then
        write(linebuf,*) nntout(1)
        do i=2,ntout
            write(str,*) nntout(i)
            linebuf=trim(linebuf)//','//str
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end if
    close(12)

stop

contains

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

        integer::i,j,k,ii,jj,iii,ik
        integer::idiff,ndiff,jnti,jsub,jjt,mem1,minmax,max,k4,k5
        integer::jmem(1:NODT),memjt(1:NODT*8),newjt(1:NODT),joint(1:NODT)

        idiff=NODT
        do j=1,NODT
            jmem(j)=0
        end do
        do j=1,NELT
            do i=1,nod
        print "('j=',i0,' i=',i0)",j,i
                jnti=kakom(j,i)
                if(jnti==0)exit
                jsub=(jnti-1)*8
                do ii=1,nod
                    if(ii==i)cycle
                    jjt=kakom(j,ii)
                    if(jjt==0)exit
                    mem1=jmem(jnti)
                    if(mem1/=0)then
                        do iii=1,mem1
                            if(memjt(jsub+iii)==jjt)exit
                        end do
                    end if
                    if(memjt(jsub+iii)/=jjt)then
                        jmem(jnti)=jmem(jnti)+1
                        memjt(jsub+jmem(jnti))=jjt
                        if(abs(jnti-jjt)>idiff)idiff=abs(jnti-jjt)
                    end if
                end do
            end do
        end do

        minmax=idiff
        do ik=1,NODT
            do j=1,NODT
                joint(j)=0
                newjt(j)=0
            end do
            max=0
            i=1
            newjt(1)=ik
            joint(ik)=1
            k=1
            do
                k4=jmem(newjt(i))
                if(k4/=0)then
                    jsub=(newjt(i)-1)*8
                    do jj=1,k4
                        k5=memjt(jsub+jj)
                        if(joint(k5)>0)cycle
                        k=k+1
                        newjt(k)=k5
                        joint(k5)=k
                        ndiff=abs(i-k)
                        if(ndiff>=minmax)exit
                        if(ndiff>max)max=ndiff
                    end do
                end if
                if(ndiff>=minmax)exit
                if(k==NODT)exit
                i=i+1
            end do
            if(ndiff>=minmax)cycle
            if(k==NODT)minmax=max
            do j=1,NODT
                jnd(j)=joint(j)
            end do
        end do

    end subroutine optnum


    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_heat
