program f90_num_seep
    !****************************************************************
    !Optimization of node number for 2D seepage flow analysis
    !****************************************************************
    implicit none
    integer::i,j,k,l,ne
    character::strcom*256           !Comment
    integer::NODT                   !Number of nodes
    integer::NELT                   !Number of elements
    integer::MATEL                  !Number of material sets
    integer::KOH                    !Number of nodes with given total head
    integer::KOQ                    !Number of nodes with given discharge
    integer::KOU                    !Number of nodes with the condition of seepage point
    integer::idan                   !Section
    integer::nod                    !Number of node of 1 element
    integer::nhen=1                 !Degree of freedom of 1 node
    integer::nt                     !Total degrees of freedom
    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 in x-direction
    real(8),allocatable::z(:)       !y-coordinate in y-direction
    integer,allocatable::nokh(:)    !Node number with given total head
    integer,allocatable::nokq(:)    !Node number with given discharge
    integer,allocatable::noku(:)    !Node number with the condition of seepage point
    real(8),allocatable::Hinp(:)    !Given total head of node
    real(8),allocatable::Qinp(:)    !Given discharge node
    real(8),allocatable::hvec(:)    !Initial value of total head
    integer,allocatable::numr(:)    !Original node number
    integer,allocatable::jnd(:)     !Revised node number
    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,KOH,KOQ,KOU,idan
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:3))
    allocate(x(1:NODT))
    allocate(z(1:NODT))
    if(0<KOH)allocate(nokh(1:KOH))
    if(0<KOQ)allocate(nokq(1:KOQ))
    if(0<KOU)allocate(noku(1:KOU))
    allocate(Hinp(1:KOH))
    allocate(Qinp(1:KOQ))
    allocate(hvec(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)
    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),z(i),hvec(i)
    end do
    if(0<KOH)then
        do i=1,KOH
            read(11,*) nokh(i),Hinp(i)
        end do
    end if
    if(0<KOQ)then
        do i=1,KOQ
            read(11,*) nokq(i),Qinp(i)
        end do
    end if
    if(0<KOU)then
        do i=1,KOU
            read(11,*) noku(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.
    !*****************************************************
    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(linebuf,*) nod,',',NODT,',',NELT,',',MATEL,',',KOH,',',KOQ,',',KOU,',',idan
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    do i=1,MATEL
        write(linebuf,*) wm(i,1),',',wm(i,2),',',wm(i,3)
        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)),',',z(numr(i)),',',hvec(numr(i))
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    if(0<KOH)then
        do i=1,KOH
            write(linebuf,*) jnd(nokh(i)),',',Hinp(i)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(0<KOQ)then
        do i=1,KOQ
            write(linebuf,*) jnd(nokq(i)),',',Qinp(i)
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    end if
    if(0<KOU)then
        do i=1,KOU
            write(linebuf,*) jnd(noku(i))
            call del_spaces(linebuf)
            write (12,'(a)') trim(linebuf)
        end do
    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_seep
