      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)
      parameter (maxlp=10000)
      integer r1(2),r(2),dr(2,maxz),iflag(2)
      integer perim(2,maxperim)
      common list(2,0:maxsize)
      integer sitepl(2*maxlp-1),perimpl(2*maxlp-1)
     c,sitech(maxsize),perimch(maxperim)
c
c-----------------------------------------------------------------------
c              Function "lp" gives the plane number
c-----------------------------------------------------------------------
c
      lp(i,j)=mod(100*i+99*j,maxlp)+maxlp
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
      ncord=4
      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            Initialize planes and chains for hash coding
c-----------------------------------------------------------------------
c
      do 130 i=1,2*maxlp-1
      sitepl(i)=0
      perimpl(i)=0
  130 continue
      sitepl(lp(0,0))=1
      do 140 i=1,maxsize
      sitech(i)=0
  140 continue
      do 150 i=1,maxperim
      perimch(i)=0
  150 continue
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,ncord
      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
      lp1=lp(r1(1),r1(2))
      n1=sitepl(lp1)
      if(n1.eq.0) go to 425
      n2=n1
  427 if(r1(1).eq.list(1,n2).and.r1(2).eq.list(2,n2)) go to 400
      n2=sitech(n2)
      if(n2.eq.0) go to 425
      go to 427
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
      n3=perimpl(lp1)
      if(n3.eq.0) go to 430
      n4=n3
  437 if(r1(1).eq.perim(1,n4).and.r1(2).eq.perim(2,n4)) go to 400
      n4=perimch(n4)
      if(n4.eq.0) go to 430
      go to 437
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)
      perimpl(lp1)=m3
      perimch(m3)=n3
      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)
      sitepl(lp1)=m2
      sitech(m2)=n1
      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
