program cluster
c
c ***********************************************************************
c     Program to generate percolation clusters by breadth first search
c ***********************************************************************
c
      parameter (maxsize=50000,maxz=4,maxchem=1501,maxperim=100000)
      parameter (maxlp=10000)
      common list(2,0:maxsize)
      real etime,time(2)
      integer iflag(2)
      character*12 out
   10 write(6,120)
  120 format(' SQ Cluster:  p, n, ntri, isd ?')
      read(5,*) p,n,ntri,isd
      write(6,121)
  121 format(' file name to record cluster sites ?')
      read(5,122) out
  122 format(a12)
      do 300 itri=1,ntri
      isd0=isd
      t0=etime(time)
      call bmake2(p,n,m2,iflag,isuc,isd)
      t1=etime(time)
      if(isuc.eq.1) then
         write(6,200) m2,t1-t0,iflag,isd0,isd
  200    format(1x,'size=',i6,', CPU=',1p,e10.3,', flags=',2i1,2x,2i12)
         open(1,file=out)
         rewind(1)
         do 250 i=1,m2
  250    write(1,251) list(1,i),list(2,i)
  251    format(1x,2(i6,1x))
         close(1)
         stop
      else
         write(6,520) itri,m2,iflag
  520    format(' trial#',i3,', size=',i5,', flag=',2i1)
      endif
  300 continue
      stop
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c     "bmake2" to make a percolation cluster by breadth first
c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      subroutine bmake2(p,n,m2,iflag,isuc,isd)
      parameter (maxsize=50000,maxz=4,maxperim=100000)
      integer r1(2),r(2),dr(2,maxz),iflag(2)
      integer perim(2,maxperim)
      common list(2,0:maxsize)
c
c-----------------------------------------------------------------------
c    The following block specifies lattice geometry for SQ.
c-----------------------------------------------------------------------
c
      data dr/0,1,0,-1,1,0,-1,0/
c
c-----------------------------------------------------------------------
c                   Initialize exit conditions
c-----------------------------------------------------------------------
c
      lim=maxsize
      if(n.lt.0) lim=min0(-n,maxsize)
      nok=min0(iabs(n),lim)
c
c-----------------------------------------------------------------------
c             Put origin on stack
c-----------------------------------------------------------------------
c
      list(1,1)=0
      list(2,1)=0
      list(1,0)=-1000000
      list(2,0)=-1000000
c
c-----------------------------------------------------------------------
c     m1 is the lower pointer where link check is being done and
c     m2 is the upper pointer where the cluster has grown up to.
c     m3 is the upper index for perim and
c-----------------------------------------------------------------------
c
      m1=1
      m2=1
      m3=0
c
c-----------------------------------------------------------------------
c          No neighbor yet for the origin.  Set flags to OK for now.
c-----------------------------------------------------------------------
c
      iflag(1)=0
      iflag(2)=0
c
c======================================================================
c               START GENERATION OF CLUSTER
c======================================================================
c
  300 continue
      r(1)=list(1,m1)
      r(2)=list(2,m1)
      mk2=m2
      mk3=m3
c
c
c
      do 400 icord=1,maxz
      r1(1)=r(1)+dr(1,icord)
      r1(2)=r(2)+dr(2,icord)
c
c-----------------------------------------------------------------------
c     Check if this neighbor is already on the stack
c-----------------------------------------------------------------------
c
      do 427 i=1,m2
  427 if(r1(1).eq.list(1,i).and.r1(2).eq.list(2,i)) go to 400
c
c----------------------------------------------------------------------
c     If any of the limits has been reached, then no more addition of
c     sites is done.
c----------------------------------------------------------------------
c
  425 if((iflag(1)+iflag(2)).gt.0) go to 500
c
c------------------------------------------------------------------------
c     Check if this site has already been decided to be a perimeter site.
c------------------------------------------------------------------------
c
      do 437 i=1,m3
  437 if(r1(1).eq.perim(1,i).and.r1(2).eq.perim(2,i)) go to 400
c
c-----------------------------------------------------------------------
c     Generate new site with probability p.
c-----------------------------------------------------------------------
c
  430 if(ranf(isd).le.p) go to 440
      m3=m3+1
      perim(1,m3)=r1(1)
      perim(2,m3)=r1(2)
      if(m3.eq.maxperim) iflag(2)=1
      go to 400
c
c------------------------------------------------------------------------
c     Add the new site on 'list'
c------------------------------------------------------------------------
c
  440 m2=m2+1
      list(1,m2)=r1(1)
      list(2,m2)=r1(2)
      if(m2.eq.lim) iflag(1)=1
c
c-----------------------------------------------------------------------
c     Advance the lower pointer m1 to check the next site on 'list'.
c-----------------------------------------------------------------------
c
  400 continue
      if(m1.eq.m2) go to 500
      m1=m1+1
      go to 300
c
c=======================================================================
c                    GENERATION ATTEMPT IS FINISHED
c=======================================================================
c
  500 if(m2.ge.nok) then
         isuc=1
      else
         isuc=0
      endif
c
c
c
      return
      end
c
c
c
      real function ranf(iseed)
      iseed=iseed*1566083941+1
      ranf=float(iseed)*2.328306e-10+0.5
      return
      end