      program clusterb
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)
      real p
      integer iflag(4)
      common /blk1/ link(maxz,maxsize)
      common /blk2/ list(2,0:maxsize),mult(maxsize)
      character*12 out
   10 write(6,120)
  120 format(' SQ Cluster:  p, n, lmax, ntri, isd ?')
      read(5,*) p,n,lmax1,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
      call pmake2(p,n,m2,lmax1,lmax,iflag,isuc,isd)
      if(isuc.eq.1) then
         write(6,200) m2,lmax,iflag,isd0,isd
  200    format(1x,'size=',i6,', gen=',i4,', flags=',4i1,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=',i2)
      endif
  300 continue
      stop
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c     subrouttine "pmake2" to make a percolation cluster on SQ
c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      subroutine pmake2(p,n,m2,lmax1,lmax,iflag,isuc,isd)
      parameter (maxsize=50000,maxz=4,maxchem=1501,maxperim=100000)
      parameter (maxlp=10000)
      real p,ranf
      integer r1(2),r(2),dr(2,maxz),iflag(4)
      integer iperim(2,maxperim),lchem(maxchem)
      common /blk1/ link(maxz,maxsize)
      common /blk2/ list(2,0:maxsize),mult(maxsize)
      integer sitepl(2*maxlp-1),perimpl(2*maxlp-1)
     c,sitech(maxsize),perimch(maxperim)
c
c-----------------------------------------------------------------------
c              Funcordtion "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 only.
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)
      llim=maxchem
      if(lmax1.lt.0) llim=min0(-lmax1,maxchem)
      lok=min0(iabs(lmax1),llim)
c
c-----------------------------------------------------------------------
c             Put origin on stack and point links to itself
c-----------------------------------------------------------------------
c
      list(1,1)=0
      list(2,1)=0
      do 120 icord=1,ncord
  120 link(icord,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     m4 is the upper index for the current generation.
c-----------------------------------------------------------------------
c
      m1=1
      m2=1
      m3=0
      m4=1
c
c-----------------------------------------------------------------------
c     lmax is the generation being checked and lchem(lmax) stores the
c     top index for the lmax-th generation.
c-----------------------------------------------------------------------
c
      lmax=1
      lchem(1)=1
c
c-----------------------------------------------------------------------
c          No neighbor yet for the origin.  Set flags to OK for now.
c-----------------------------------------------------------------------
c
      mult(1)=0
      iflag(1)=0
      iflag(2)=0
      iflag(3)=0
      iflag(4)=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)
      mcord=2*mod(icord,2)+icord-1
c
c     Check if this neighbor is among those already linked.
c
      do 422 j=1,ncord
      jj=link(j,m1)
      if(list(1,jj).eq.r1(1).and.list(2,jj).eq.r1(2)) go to 400
  422 continue
c
c-----------------------------------------------------------------------
c     Check if this neighbor is already on the list but not yet linked.
c-----------------------------------------------------------------------
c
      lp1=lp(r1(1),r1(2))
      n1=sitepl(lp1)
      if(mk2.eq.m1) go to 425
      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 426
      n2=sitech(n2)
      if(n2.eq.0) go to 425
      go to 427
  426 mult(m1)=mult(m1)+1
      link(icord,m1)=n2
      mult(n2)=mult(n2)+1
      link(mcord,n2)=m1
      go to 400
c
c----------------------------------------------------------------------
c     If any of the limits has been reached, then no more addition of
c     sites is done.  Only the linking must be finished up.
c----------------------------------------------------------------------
c
  425 if((iflag(1)+iflag(2)+iflag(3)+iflag(4)).gt.0) go to 400
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.iperim(1,n4).and.r1(2).eq.iperim(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
      iperim(1,m3)=r1(1)
      iperim(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' and link it.
c------------------------------------------------------------------------
c
  440 m2=m2+1
      list(1,m2)=r1(1)
      list(2,m2)=r1(2)
      mult(m1)=mult(m1)+1
      link(icord,m1)=m2
      mult(m2)=1
      do 450 jcord=1,ncord
  450 link(jcord,m2)=0
      link(mcord,m2)=m1
      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
      if(m1.eq.m4) then
         m4=m2
         lmax=lmax+1
         lchem(lmax)=m4
         if(lmax.eq.llim) iflag(3)=1
      endif
      m1=m1+1
      go to 300
c
c=======================================================================
c                    GENERATION ATTEMPT IS FINISHED
c=======================================================================
c
  500 if(m2.ge.nok.and.lmax.ge.lok) 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
