!!$!!!!!!!!!!!!!!!!!!!!!!!!
!!$!!!! Copyright 2002 Lorenzo Giada and Matteo Marsili
!!$!!!! lgiada@mpikg-golm.mpg.de    marsili@sissa.it
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$  This program is free software; you can redistribute it and/or modify
!!$  it under the terms of the GNU General Public License as published by
!!$  the Free Software Foundation; either version 2 of the License, or
!!$  (at your option) any later version.
!!$
!!$  This program is distributed in the hope that it will be useful,
!!$  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!$  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!$  GNU General Public License for more details.
!!$
!!$  You should have received a copy of the GNU General Public License
!!$  along with this program; if not, write to the Free Software
!!$  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module mod_regrpsan
  
  integer :: N, N0, k0
  integer :: be_in,be_fin,t_wri,t_temp,sg,nr_f
  integer :: idum,rns,clnm
  integer, allocatable, dimension(:)   :: a,m,ags,m_gs,a0,m0,a_t
  integer, allocatable, dimension(:,:) :: indx

  real, allocatable, dimension(:)   :: x,e,x_gs,e_gs
  real, allocatable, dimension(:,:) :: C,C0
  real    :: egs, cf, t_wrir

  character(12) :: filfin,filclu,fildat,filgs,filout
  character(9)  :: filord
  character(3)  :: filxxx
  character(7)  :: filin
  character(12), allocatable, dimension(:) :: cells
  character(1)  :: u

end module mod_regrpsan
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program group_fit

!!$ calculates the gs by simulated annealing:
!!$ 1. sweep of the lattice (random or sequential)
!!$ 2. merging of 2 clusters chosen at random
!!$ 3. splitting of a cluster by a deterministic procedure (max covariance)
!!$ 4. do 1. -> 3. a few times, then increase or decrease the beta
!!$ 5. go up and down in beta a certain number of times
!!$ 6. last sweeps at very low temperature
!!$ 7. SAVES the gs config and write a readable table of clusters
!!$ it allows for RECLUSTERING depending on flag k0
!!$ REQUIRES the files:
!!$ regrpsan.par
!!$ covariance matrix
!!$ names of series (optional) (written as 2 columns: number, name)
!!$ names of initial configuration (optional)
!!$ WRITES the files:
!!$ .ent.gsK the energy at a given temperature and other run-time info
!!$ .now.gsK the GROUND STATE energy up to now. (it is never initialized)
!!$ .cnf.gsK the ground state configuration (at the end)
!!$ .clt.gsK the table with readable clusters (at the end)
!!$ .cls.gsK index of sets ordered per cluster (at the end)

  use mod_regrpsan
  implicit none

  integer :: nr,sw,ugh,out
  integer, allocatable, dimension(:) :: indi
  integer :: mi_new,mf_new,ai,af,nupd

  real    :: xi_new,xf_new,ei_new,ef_new,de,beta
  integer :: i,j,k,l,ik,nt,t,ibeta
  real, allocatable, dimension(:)    :: ee
  real    :: etm
  real    :: rand
  external rand

!!$ READ data and parameters and name of files from regrpsan.par
  call ini(ugh)
  if (ugh.eq.1) goto 998
  allocate(ee(N))
  nr=1
  sg=1   !+1 for decreasing T, -1 for increasing T
  egs=10000.0
  t=0
  call gs(t,0)

  open(20,file=filout,status='unknown')
!!$ TEMPERATURE CICLE
55 call cic_temp(t)
  cf=1./cf
  sw=be_in
  be_in=be_fin
  be_fin=sw
!!$ decides whether cicles are finished
  if (sg.eq.-1) then
     if (nr.eq.nr_f) go to 16
     nr=nr+1
  end if
!!$ if not, invert direction
  sg=-sg
  go to 55   ! another temperature cicle
16 continue

!!$ last Cycle at Low Temperature starting from Ground State
  do i=1,N
     a(i)=ags(i)
  end do
  call cond
  be_in=be_fin*.8           !80% of T_f
  be_fin=be_fin*1.4         !+40% of T_f
  sg=1
  call cic_temp(t)
  call gs(t,0)
!!$ now SORT the resulting clusters in order of internal energy
  do i=1,N
     if(m_gs(i).gt.0) then
        ee(i)=e_gs(i)/real(m_gs(i))
     else
        ee(i)=1.0
     end if
  end do
  a=ags ; e=e_gs ; x=x_gs ; m=m_gs
  allocate(indi(N))
  call indicize(N,e,indi,out)
  if (out>0) then
     print *, 'nothing done...'
     goto 998
  end if
  do i=1,N
     e_gs(i)=e(indi(i))
     do j=1,N
        if(a(i).eq.indi(j)) ags(i)=j
     end do
     x_gs(i)=x(indi(i))
     m_gs(i)=m(indi(i))
  end do
  a=ags
  call cond
  deallocate(ee)

!!$ Saves the GS in file filfin (...cnf...)
  open(12,file=filfin,status='unknown')
  do i=1,N0
     if (k0.eq.0) then
        write(12,*) ags(i)
     else    ! case for reclustering
        write(12,*) ags(a0(i))
     end if
  end do
  close(12)
  k=0
  do i=1,N
     if(m(i).gt.0) k=k+1
  end do
  write(6,*) 'there are ',k,'clusters, and the GS energy is ',sum(e)/real(N)
!!$ Makes file with readable clusters SORTED in file filclu (...clt...)
!!$ and only index of series in clusters in file fildat (...cls...)
  open(12,file=filclu)
  open(13,file=fildat)
  write(12,*) 'there are ',k,'clusters, and the GS energy is ',sum(e)/real(N)
  allocate (ee(N))
  ee=0.0
  do j=1,N
     if (m(j).gt.0)then
        if (m(j).eq.1)then
           write(12,333) j,m(j),  0., x(j)/real(m(j)),e(j)/real(m(j))
        else
           write(12,333) j,m(j),&
                &(x(j)-m(j))/(m(j)**2-x(j)),x(j)/real(m(j)),e(j)/real(m(j))
        end if
333     format (' Cluster', I5, ', size=',I5, ', g=', F10.5,&
             &', c/n=',F10.5, ', e/n=',F10.5)
!!$ calculate extraction energies
        do i=1,m(j)
           xi_new=x(j)
           do k=1,m(j)
              xi_new=xi_new-C(indx(j,i),indx(j,k))
           end do
           call energy(m(j)-1,xi_new,etm)
           ee(i)=e(j)-etm
        end do
!!$ sort for minimal of these extraction energies
        do l=1,m(j)
           ef_new=1000.0
           do i=1,m(j)
              if (ef_new.gt.ee(i)) then
                 k=i
                 ef_new=ee(i)
              end if
           end do
           if (k0.eq.0) then
              if (clnm.eq.1) then
                 write(12,*) indx(j,k),cells(indx(j,k)), ee(k)
                 write(13,*) indx(j,k),cells(indx(j,k)), ee(k)
              else
                 write(12,*) indx(j,k), ee(k)
                 write(13,*) indx(j,k), ee(k)
              end if
           else  ! case of reclustering
              do i=1,N0
                 if (a0(i).eq.indx(j,k))then
                    if (clnm.eq.1) then
                       write(12,*) i,indx(j,k),cells(i), ee(k)
                       write(13,*) i,indx(j,k),cells(i), ee(k)
                    else
                       write(12,*) i,indx(j,k), ee(k)
                       write(13,*) i,indx(j,k), ee(k)
                    end if
                 end if
              end do
           end if
           ee(k)=1000.0
        end do
        write(12,*) '  '
     end if
  end do
  close(12)
  close(13)
  close(20)

998 continue
end program group_fit
!!$cccccccCCCCCCCCcccccccCCCCCCcccccccCCCCCccccccc
subroutine sweep(beta,t,nupd)
!!$ SWEEP of the lattice
  use mod_regrpsan
  implicit none

  real      :: beta
  integer   :: nupd
  integer   :: i,k,j,l,ai,af,t,t0,ncls
  integer   :: mi_new,mf_new
  real      :: de,dxf
  real      :: xi_new,ei_new,xf_new,ef_new
  real      :: rand,et1,et2
  external rand

  do i=1,N
     if (rns.eq.1) then     !random sweep
        l=int(rand(idum)*N)+1
     else
        l=i                 !sequential sweep
     end if
     ai=a(l)
11   af=N*rand(idum)+1
     if (af.eq.ai) go to 11 ! no change of cluster: RETRY
     if (m(af).eq.0.and.m(ai).eq.1) go to 11 ! just a relabelling: RETRY
     if (m(ai).eq.1) then
        mi_new=0
        xi_new=0.0
        ei_new=0.0
     else
        xi_new=x(ai)
        do k=1,m(ai)
           j=indx(ai,k)
           xi_new=xi_new-C(l,j)
        end do
        mi_new=m(ai)-1
        call energy(mi_new,xi_new,ei_new)
     end if
     mf_new=m(af)+1
     dxf=0.0
     do k=1,m(af)
        j=indx(af,k)
        dxf=dxf+C(l,j)
     end do
     xf_new=x(af)+1.+dxf
     call energy(mf_new,xf_new,ef_new)
     de=ei_new+ef_new-e(ai)-e(af)
!!$ METROPOLIS
!!$   update index: 1. remove the index i from cluster a(i)
     if (de.lt.0) then  ! accept move
        do k=1,m(ai)-1
           if (indx(ai,k).eq.l) indx(ai,k)=indx(ai,m(ai))
        end do
        indx(af,mf_new)=l   ! 2. add new index
!!$    update cluster variables
        m(ai)=mi_new
        x(ai)=xi_new
        e(ai)=ei_new
        m(af)=mf_new
        x(af)=xf_new
        e(af)=ef_new
!!$    update Potts variable
        a(l)=af
        nupd=nupd+1
        call gs(beta,t)
     else if (rand(idum).le.exp(-beta*de)) then ! accept move
        do k=1,m(ai)-1
           if (indx(ai,k).eq.l) indx(ai,k)=indx(ai,m(ai))
        end do
        indx(af,mf_new)=l   ! 2. add new index
!!$    update cluster variables
        m(ai)=mi_new
        x(ai)=xi_new
        e(ai)=ei_new
        m(af)=mf_new
        x(af)=xf_new
        e(af)=ef_new
!!$    update Potts variable
        a(l)=af
        nupd=nupd+1
     end if
  end do                    !end lattice sweep

  return
end subroutine sweep
!!$cccccccCCCCCCCCccccccCCCCCccccccccCCCCcccCCCCCc
subroutine cic_temp(t)
  use mod_regrpsan
  implicit none

  integer :: nupd,nmupd,nsupd
  integer :: i,k,j,l,t,ai,af,t0,ncls
  real    :: beta,etot,de,dxf,mav
  real    :: xi_new,ei_new,xf_new,ef_new
  integer :: mi_new,mf_new
  real    :: rand
  external rand
!!$ccc CICLE IN TEMPERATURE
  beta=float(be_in)
  do
     t0=t
     nupd=0
     nmupd=0
     nsupd=0
!!$ccc Lattice sweep         
10   call sweep(beta,t,nupd)

!!$cccCCCcccCCC MERGING cccCCC
     call merge(beta,t,nmupd)

!!$cccCCCcccCCC SPLITTING cccCCC
     call split(beta,t,nsupd)
     t=t+1
     if (mod(t,t_wri).eq.0) then
        etot=0.0
        mav=0.0
        ncls=0
        etot=sum(e)
        do i=1,N
           if (m(i).gt.0) ncls=ncls+1
           mav=mav+m(i)
        end do
!!$ WRITES time, beta, energy, avg. # updates, # clusters,
!!$ avg. occupation of cluster, # merging, # splitting
!!$ in file filout (xxx.ent.gsK)
        write(20,990) real(t)/t_wrir,beta,etot/float(N),&
             &real(nupd)/(t_wrir*real(N)),ncls,&
             &real(mav)/real(ncls),real(nmupd)/t_wrir,real(nsupd)/t_wrir
        nupd=0
        nmupd=0
        nsupd=0
990     format(f8.2,f11.4,f14.8,f11.5,i6,f10.4,f9.4,f9.4)
     end if
     if (t-t0.lt.t_temp) go to 10
     if (sg.eq.1) then
        if (int(beta).ge.be_fin) go to 15
     else
        if (int(beta).le.be_fin) go to 15
     end if
     beta=beta/cf
  end do
15 return
end subroutine cic_temp
!!$ccccccCCCCCCcccccccccccccCCCCCCCCcccccccCCCCCC
subroutine merge(beta,t,nmupd)
!!$cccc   TRY the MERGING of two CLUSTERS
  use mod_regrpsan
  implicit none

  integer :: i,k,j,l,ai,af,nm,t,nmupd
  integer :: mi_new,mf_new

  real    :: xi_new,ei_new,de,beta,xf_new,ef_new
  real    :: rand
  external rand

!!$ Check if there are at least 2 clusters with at least 2 elements
  nm = 0
  do i=1,N
     if (m(i).ge.2) then
        nm = nm + 1
        if (nm.eq.2) go to 20
     endif
  end do
  return
!!$ picks at random 2 clusters with at least 2 elements
20 ai=N*rand(idum)+1
  if (m(ai).lt.2) go to 20
21 af=N*rand(idum)+1
  if ((m(af).lt.2).or.(af.eq.ai)) go to 21

!!$c New parameters for cluster ai
  mi_new=0
  xi_new=0.0
  ei_new=0.0
       
!!$ New parameters for cluster af
  mf_new=m(af)+m(ai)
  do k=1,m(ai)
     indx(af, m(af) + k) = indx(ai, k)
  end do
  xf_new=x(af)
  do k=m(af)+1,mf_new
     do l=1,k
        xf_new=xf_new+C(indx(af,k),indx(af,l))
     end do
  end do
  call energy(mf_new,xf_new,ef_new)

!!$  Metropolis
  de=ef_new-e(ai)-e(af)
  if (rand(idum).le.exp(-beta*de)) then ! accept move
!!$    update Potts variable
     do k=1,m(ai)
        a(indx(ai, k)) = af
     end do

!!$     update cluster variables
     m(ai)=mi_new
     x(ai)=xi_new
     e(ai)=ei_new
     m(af)=mf_new
     x(af)=xf_new
     e(af)=ef_new

     nmupd=nmupd+1
     if (de.lt.0) then
        call gs(beta,t)
     end if
  end if
  return
end subroutine merge
!!$cccccccccCCCCCCCCCCCCCccccccccccccCCCCCCCCCCC
subroutine split(beta,t,nsupd)
!!$ DETERMINISTIC SPLITTING OF A CLUSTER !!!!!!!!!
  use mod_regrpsan
  implicit none

  integer :: i,k,j,ai,t
  integer :: mi_new,nsupd
  integer, dimension(2) :: aff,ii,mff

  real    :: xi_new,ei_new,de,beta
  real    :: xff(2),eff(2),cmax
  real    :: rand
  external rand

!!$ picks at random a cluster with more than 3 elements
  do i=1,N
     if (m(i).gt.3) go to 6
  end do
  return
6 continue
  ai=N*rand(idum)+1
  if (m(ai).lt.4) go to 6
!!$ picks 2 empty clusters
  k=1
  do i=1,N
     if (m(i).eq.0) then
        aff(k)=i
        k=k+1
        if (k.gt.2) goto 7
     end if
  end do
!!$ picks the 2 most distant elements based on C_i,j
7 continue
  cmax=2.
  do i=1,m(ai)-1
     do j=i+1,m(ai)
        if(cmax.gt.C(indx(ai,i),indx(ai,j)))then
           cmax=C(indx(ai,i),indx(ai,j))
           ii(1)=indx(ai,i)
           ii(2)=indx(ai,j)
           indx(aff(1),1)=ii(1)
           indx(aff(2),1)=ii(2)
           mff(1)=1
           mff(2)=1
        end if
     end do
  end do
  xff(1)=1.
  xff(2)=1.
!!$ divides the cluster according to distance (C_i,j)
  do i=1,m(ai)
     j=indx(ai,i)
     if (j.ne.ii(1).and.j.ne.ii(2)) then
        if (C(ii(1),j).gt.C(ii(2),j)) then
           mff(1)=mff(1)+1
           indx(aff(1),mff(1))=j
           do k=1,mff(1)
              xff(1)=xff(1)+C(indx(aff(1),k),j)
           end do
        else
           mff(2)=mff(2)+1
           indx(aff(2),mff(2))=j
           do k=1,mff(2)
              xff(2)=xff(2)+C(indx(aff(2),k),j)
           end do
        end if
     end if
  end do
  mi_new=0
  xi_new=0.0
  ei_new=0.0
  do i=1,2
     call energy(mff(i),xff(i),eff(i))
  end do
!!$   METROPOLIS
  de=eff(1)+eff(2)-e(ai)
  if (rand(idum).le.exp(-beta*de)) then ! accept move
!!$ update Potts variables
     do i=1,2
        do k=1,mff(i)
           a(indx(aff(i),k))=aff(i)
        end do
     end do
!!$ update Cluster variables
     m(ai)=mi_new
     x(ai)=xi_new
     e(ai)=ei_new
     do i=1,2
        m(aff(i))=mff(i)
        x(aff(i))=xff(i)
        e(aff(i))=eff(i)
     end do
     nsupd=nsupd+1
     if (de.lt.0) then
        call gs(beta,t)
     end if
  end if
  return
end subroutine split
!!$cccccccCCCCCCCCCcccccccccCCCCCCCccccccCCCCCCCccccc
subroutine cond
  use mod_regrpsan
  implicit none
  integer :: i,j,k

!!$ Calculates parameters for a given cluster configuration,
!!$ which is in the array a(N)
  m=0
  x=0.0
  e=0.0
  do i=1,N
     m(a(i))=m(a(i))+1
  end do
  do i=1,N
     k=1
     do j=1,N
        if (a(j).eq.i) then
           indx(i,k)=j
           k=k+1
        end if
     end do
  end do
  do i=1,N
     if (m(i).eq.1) x(i)=1.
     if (m(i).gt.1) then
        do j=1,m(i)
           do k=j,m(i)
              x(i)=x(i)+C(indx(i,j),indx(i,k))
           end do
        end do
     end if
  end do
  do i=1,N
     if (m(i).gt.1) then
        call energy(m(i),x(i),e(i))
     end if
  end do
  return
end subroutine cond
!!$ccccccCCCCCCccccccccccCCCCCCCcccccccccCCCCCCCCCccccc
subroutine ini(ugh)
  use mod_regrpsan
  implicit none

  integer :: i,j,k,ij,ugh, flag
  real    :: rand
  external rand

!!$ READ parameters
  open(10,file='grpsan.par',status='old',iostat=ugh)
  if(ugh.ne.0) then
     write(6,*) 'configuration file grpsan.par is missing. will STOP'
     goto 999
  end if
  read(10,*) N0
  read(10,*) filxxx
  read(10,*) k0
  read(10,*) clnm
  read(10,*) be_in
  read(10,*) be_fin
  read(10,*) t_wri
  read(10,*) t_temp
  read(10,*) nr_f
  read(10,*) cf
  read(10,*) idum
  read(10,*) flag
  read(10,*) rns
  close(10)

  allocate(C0(N0,N0),a0(N0),m0(N0),cells(N0))
  if(k0.eq.0) then
     N=N0
  else     ! compute N for number of clusters
     write(u,'(I1)') k0-1     
     open(20,file=filxxx//'.cnf.gs'//u,status='old',iostat=ugh)
     if (ugh.ne.0) then
        write(6,*) 'cannot find initial configuration for reclustering.'
        write(6,*) 'check variables k0 in file regrpsan.par. will STOP'
        goto 999
     end if
     do i=1,N0
        read(20,*) a0(i)
     end do
     close(20)
     m0=0
     do i=1,N0
        m0(a0(i))=m0(a0(i))+1
     end do
     N=0
     do i=1,N0
        if (m0(i).ne.0) N=N+1
     end do
     write(6,*)'N=', N
  end if
  write(u,'(I1)') k0
  allocate(indx(N,N),C(N,N))
  allocate(ags(N),x(N),e(N),x_gs(N),m_gs(N),e_gs(N))

  filin(1:3)=filxxx
  filin(4:7)='.cij'
  filout(1:3)=filxxx
  filout(4:12)='.ent.gs'//u
  filfin(1:3)=filxxx
  filfin(4:12)='.cnf.gs'//u
  filclu(1:3)=filxxx
  filclu(4:12)='.clt.gs'//u
  fildat(1:3)=filxxx
  fildat(4:12)='.cls.gs'//u
  filord(1:3)=filxxx
  filord(4:9)='.names'
  filgs(1:3)=filxxx
  filgs(4:12)='.cnf.gs'//u
  t_wrir=real(t_wri)

!!$ LOAD INITIAL CORRELATION MATRIX
  open(10,file=filin,status='old', iostat=ugh)
  if(ugh.ne.0) then
     write(6,*) 'cannot find file ',filxxx,'.cij with correlation matrix'
     write(6,*) 'please check variable filxxx in input file regrpsan.par'
     write(6,*) 'now I have to STOP'
     go to 999
  end if
  do i=1,N0-1
     do j=i+1,N0
        read(10,*) C0(i,j)
        C0(i,j)=2*C0(i,j)
        C0(j,i)=C0(i,j)
     end do
     C0(i,i)=1.0
  end do
  C0(N0,N0)=1.0
  close(10)

!!$ load name of series
  if(clnm.eq.1) then
     open (12,file=filord,status='old', iostat=ugh)
     if(ugh.ne.0) then
        write(6,*) 'cannot find file with name of series'
        write(6,*) 'please check variable clnm in input file grpsan.par'
        write(6,*) 'will continue without names'
        clnm=0
        goto 888
     end if
     do i=1,N0
        read(12,*) j,cells(i)
     end do
  end if
888 continue
  if (k0.ne.0) then   ! compute new c_ij for RECLUSTERING
     C=0.
     do i=1,N
        do j=i,N
           do k=1,N0
              if (a0(k).eq.i) then
                 do ij=1,N0
                    if (a0(ij).eq.j) then
                       C(i,j)=C(i,j)+C0(k,ij)
                    end if
                 end do
              end if
           end do
        end do
     end do
     do i=1,N-1
        do j=i+1,N
           if (C(i,j).ne.0) then
              C(i,j)=C(i,j)/sqrt(C(i,i)*C(j,j))
              C(j,i)=C(i,j)
           end if
        end do
        C(i,i)=1.
     end do
     C(N,N)=1.
  else       ! k0=0 case
     C=C0
  end if
  deallocate(C0)

!!$ initial condition
  allocate(a(N),m(N))
  if (flag.eq.1)then
     open(20,file=filgs,status='old',iostat=ugh)
     if (ugh.ne.0) then
        write(6,*) 'cannot find initial configuration. check variables'
        write(6,*) 'filgs and flag in file grpsan.par. will STOP'
        goto 999
     end if
     allocate(a_t(N0))
     do i=1,N0
        read(20,*) a_t(i)
     end do
     close(20)
     if(k0.ne.0) then
        do i=1,N
           do j=1,N0
              if (a0(j).eq.i) then
                 a(i)=a_t(j)
                 go to 50
              end if
           end do
50         continue
        end do
     else
        a=a_t
     end if
     deallocate(a_t)
  else if (flag.eq.0) then
     do i=1,N
        a(i)=i
     end do
  else if (flag.eq.2) then
     do i=1,N
        a(i)=int(rand(idum)*N)+1
     end do
  else 
     write(6,*) 'I do not understand the value of flag'
     write(6,*) 'in file grpsan.par: 0=uniform initial condition'
     write(6,*) '1=read from file in filgs'
     write(6,*) '2=random initial conditions'
     write(6,*) 'will STOP'
     go to 999
  end if
  call cond
  call gs(0,0)
  ugh=0
  return

999 ugh=1
  return
end subroutine ini
!!$cccccccCCCCCCCccccccccCCCCCcccccCCCCCCCCC
subroutine gs(beta,t)
!!$ updates the ground-state
  use mod_regrpsan   
  implicit none

  integer :: k,t
  real    :: etot,beta
  
  etot=sum(e)
  if (etot.lt.egs) then
     ags=a
     x_gs=x
     m_gs=m
     e_gs=e
     egs=etot
!!$ SAVE the GROUND STATE energy
     open(12,file=filxxx//'.now.gs'//u,position='append')
     write(12,*) egs/float(N),beta,t/t_wrir
     close(12)
  end if
  return
end subroutine gs
!!$cccccccccCCCCCCCCCCCCCCcccccccccccccCCCCCCCc
subroutine energy(m,x,e)
!!$ calculate the energy of a cluster
  implicit none
  integer :: m
  real    :: x,e,chi
  if (m.le.1) then
     e=0.0
!!$cCccc controllo che c_s<n_s
  else if (x.le.m) then
     e=0.0
  else
     chi=x/float(m)
     e=log(chi)+(real(m)-1)*log((real(m)-chi)/(real(m)-1.0))
  end if
  return
end subroutine energy
!!$cccccccccCCCCCCCCCCcccccccCCCCCcccccccccccc
SUBROUTINE indicize(n,arr,indx,out)
  implicit none
  INTEGER :: n,indx(n),M,STACK,out
  REAL    :: arr(n)
  PARAMETER (M=7,STACK=510)
  INTEGER :: i,indxt,ir,itemp,j,js,k,l,istack(STACK)
  REAL    :: a
  out=0
  DO j=1,n
     indx(j)=j
  END DO
  js=0
  l=1
  ir=n
1 IF(ir-l.LT.M)THEN
     DO j=l+1,ir
        indxt=indx(j)
        a=arr(indxt)
        DO i=j-1,1,-1
           IF(arr(indx(i)).LE.a)GOTO 2
           indx(i+1)=indx(i)
        END DO
        i=0
2       indx(i+1)=indxt
     END DO
     IF(js.EQ.0)RETURN
     ir=istack(js)
     l=istack(js-1)
     js=js-2
  ELSE
     k=(l+ir)/2
     itemp=indx(k)
     indx(k)=indx(l+1)
     indx(l+1)=itemp
     IF(arr(indx(l+1)).GT.arr(indx(ir)))THEN
        itemp=indx(l+1)
        indx(l+1)=indx(ir)
        indx(ir)=itemp
     ENDIF
     IF(arr(indx(l)).GT.arr(indx(ir)))THEN
        itemp=indx(l)
        indx(l)=indx(ir)
        indx(ir)=itemp
     ENDIF
     IF(arr(indx(l+1)).GT.arr(indx(l)))THEN
        itemp=indx(l+1)
        indx(l+1)=indx(l)
        indx(l)=itemp
     ENDIF

     indxt=indx(l)
     a=arr(indxt)

     do i=l+2,n
        IF(arr(indx(i))<a)GOTO 3
        do j=ir-1,1,-1
           IF(arr(indx(j))<=a)exit
        end do
        IF(j<i)GOTO 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
3       continue
     end do

5    indx(l)=indx(j)
     indx(j)=indxt
     js=js+2
     IF(js>STACK)then 
        write(6,*) 'STACK too small in indicize'
        out=10
        return
     end IF
     IF(ir-i+1 >= j-l)THEN
        istack(js)=ir
        istack(js-1)=i
        ir=j-1
     ELSE
        istack(js)=j-1
        istack(js-1)=l
        l=i
     ENDIF
  ENDIF
  GOTO 1
END SUBROUTINE indicize
!!$cccccccccCCCCCCCCCCCcccccccccccccCCCCCCCCCCC

FUNCTION rand(dummy)
  implicit none
  INTEGER dummy,MI1,IM2,IMM1,A1,A2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
  REAL rand,AM,EPS,RNMX
  PARAMETER (MI1=2147483563,IM2=2147483399,AM=1./MI1,IMM1=MI1-1,&
       &A1=40014,A2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
       &NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
  INTEGER dummy2,j,k,iv(NTAB),iy
  SAVE iv,iy,dummy2
  DATA dummy2/123456789/, iv/NTAB*0/, iy/0/
  if (dummy <= 0) then
     dummy=max(-dummy,1)
     dummy2=dummy
     do j=NTAB+8,1,-1
        k=dummy/IQ1
        dummy=A1*(dummy-k*IQ1)-k*IR1
        if (dummy < 0) dummy=dummy+MI1
        if (j <= NTAB) iv(j)=dummy
     end do
     iy=iv(1)
  endif
  k=dummy/IQ1
  dummy=A1*(dummy-k*IQ1)-k*IR1
  if (dummy < 0) dummy=dummy+MI1
  k=dummy2/IQ2
  dummy2=A2*(dummy2-k*IQ2)-k*IR2
  if (dummy2 < 0) dummy2=dummy2+IM2
  j=1+iy/NDIV
  iy=iv(j)-dummy2
  iv(j)=dummy
  if(iy.lt.1)iy=iy+IMM1
  rand=min(AM*iy,RNMX)
  return
END FUNCTION rand
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$ END OF CODE