program f90_num_mesh
    implicit none
    integer,allocatable::kakom(:,:)
    integer,allocatable::jnd(:),numr(:)
    real(8),allocatable::x(:),y(:)
    integer::NODT,NELT,nod
    integer::ne,i,j,k1,k2,k3,k4
    character :: linebuf*1000
    character(len=50)::fnameR,fnameW,dummy

    call getarg(1,dummy); read(dummy,*) nod
    call getarg(2,fnameR) !Input file name
    call getarg(3,fnameW) !Output file name

    open(11,file=fnameR,status='old')
    read(11,*) NODT,NELT
    allocate(kakom(1:NELT,1:nod))
    allocate(x(1:NODT),y(1:NODT))
    allocate(jnd(1:NODT),numr(1:NODT))
    do ne=1,NELT
        read(11,*) (kakom(ne,j),j=1,nod)
    end do
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    close(11)

    print "('NODT=',i0,' NELT=',i0,' nod=',i0)",NODT,NELT,nod

    !Minimization of band width
    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

    do i=1,NODT
        print "('original=',i0,' new=',i0)",i,jnd(i)
    end do

    open(12,file=fnameW,status='replace')
    write(12,'(2i6)') NODT,NELT
    select case(nod)
        case(3)
            do ne=1,NELT
                k1=kakom(ne,1)
                k2=kakom(ne,2)
                k3=kakom(ne,3)
                kakom(ne,1)=jnd(k1)
                kakom(ne,2)=jnd(k2)
                kakom(ne,3)=jnd(k3)
                write(12,'(3i6)') (kakom(ne,j),j=1,nod)
            end do
        case(4)
            do ne=1,NELT
                k1=kakom(ne,1)
                k2=kakom(ne,2)
                k3=kakom(ne,3)
                k4=kakom(ne,4)
                kakom(ne,1)=jnd(k1)
                kakom(ne,2)=jnd(k2)
                kakom(ne,3)=jnd(k3)
                kakom(ne,4)=jnd(k4)
                write(12,'(4i6)') (kakom(ne,j),j=1,nod)
            end do
    end select
    do i=1,NODT
        write(12,'(2e15.7)') x(numr(i)),y(numr(i))
    end do
    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

end program f90_num_mesh
