      subroutine clustercounting(sizemax,clusternofinal)
      implicit none
      common/var/s(300)
      common/dim/lx
      common/v1/buffer(300)
      common/v2/buffersize,clusterno
      common/v3/lat(300),pt(300,2)
      integer lat,pt,buffer,lx,ly,s
      integer clusterno,buffersize,size(300),sizemax
      integer i,j,n,clusternofinal
ccccccccccccccccccccccccc
cc  open result file
      open(unit=30,file='clusters',
     &      form='formatted',status='unknown')
cccccccccc initialization:     
      do 1 i=1,lx
        lat(i)=0
 1    continue
cc  read data: 
      n=0
      do 2 i=1,lx
        lat(i)=(s(i)+1)/2
	if(lat(i).eq.1) then
	  n=n+1
	  lat(i)=n
	  pt(n,1)=i
	  pt(n,2)=0
	endif
 2    continue
      if(n.eq.0) then
       write(*,*) 'Empty pattern'
       stop
      endif
cc Start:
      buffer(1)=1
      pt(1,2)=1
      clusterno=1 
      buffersize=1
      call clustexp
      do 3 i=2,n
        if(pt(i,2).eq.0) then
	 clusterno=clusterno+1
	 pt(i,2)=clusterno
         buffer(1)=i
	 buffersize=1
	 call clustexp
	endif
 3    continue
cccccccccc
 1000 continue
      close(30) 
ccccccccccccccccccccccccccc
ccc finding largest cluster:
      open(unit=30,file='clusters',
     &      form='formatted',status='old')
      do 1001 n=1,clusterno
       read(30,*) i,size(n)
 1001 continue
      close(30)
      sizemax=0
      do 1002 n=1,clusterno
       if(size(n).gt.sizemax) sizemax=size(n)
 1002 continue
      clusternofinal=clusterno
      end 
      
      subroutine clustexp
      implicit none
      common/dim/lx
      common/v1/buffer(300)
      common/v2/buffersize,clusterno
      common/v3/lat(300),pt(300,2)
      integer buffer,buffersize,clusterno,lx,ly
      integer lat,pt
      integer clustersize,n
      clustersize=0
      do 11 n=1,lx
       clustersize=clustersize+1
       call vecinos(buffer(1))
       call stack
       if(buffersize.eq.0) goto 12
 11   continue
 12   continue
      write(30,*) clusterno, clustersize
      end
      
      subroutine stack
      implicit none
      common/v1/buffer(300)
      common/v2/buffersize,clusterno
      integer buffer,buffersize,clusterno,k
      do 13 k=1, buffersize
       buffer(k)=buffer(k+1)
 13   continue
      buffersize=buffersize-1
      end
      
      subroutine vecinos(i)
      implicit none
      common/v1/buffer(300)
      common/v2/buffersize,clusterno
      common/v3/lat(300),pt(300,2)
      common/bond/statusbond(300,300),nbneigh(300),
     & ineigh(300,300) 
      integer buffer
      integer buffersize,clusterno,lat,pt,i,ii,k
      integer statusbond,nbneigh,ineigh
      do 14 k=1,nbneigh(pt(i,1))
       ii=ineigh(pt(i,1),k)
       if(statusbond(pt(i,1),k).eq.1.
     &      AND.pt(lat(ii),2).eq.0) then
        buffersize=buffersize+1
	buffer(buffersize)=lat(ii)
	pt(lat(ii),2)=clusterno
       endif
 15    continue
 14   continue
      end
