program cluster
c
c ***********************************************************************
c     Program to generate percolation clusters only
c ***********************************************************************
c
      parameter (maxedge=2000,maxnn=20000,maxsize=1000000)
      common /blk0/ nn(maxnn),list(2,maxsize)
      dimension xiav(2)
      character*12 out,out0
   10 write(6,120)
  120 format(' SQ Cluster:  p, nedge, ntri, isd ?')
      read(5,*) p,nedge,ntri,isd
      if(nedge.gt.maxedge) go to 10
      write(6,121)
  121 format(' file name to record cluster sites ?')
      read(5,122) out
  122 format(a12)
      write(6,123)
  123 format(' file name to record cluster size distribution ?')
      read(5,122) out0
      do 300 itri=1,ntri
      isd0=isd
      call make2(p,nmax,nedge,isuc,isd,nsize)
      if(isuc.eq.1) then
         open(2,file=out0)
         rewind(2)
         call dist(nmax,xiav,px,nc,ng)
         conc=float(ng)/float(nedge**2)
         write(2,5) p,isd0,isd,nedge,ng,px,xiav(2),xiav(1),conc,nc
         write(6,5) p,isd0,isd,nedge,ng,px,xiav(2),xiav(1),conc,nc
    5    format(f9.4,6x,2i11,i5,i8,1p,4e10.3,i6)
         close(2)
         open(1,file=out)
         rewind(1)
         do 410 i=1,nsize
         write(1,411) (list(j,i),j=1,2)
  410    continue
  411    format(1x,2(i6,1x))
         close(1)
         stop
      else
         write(6,520) itri,nmax,isuc
  520    format(' trial#',i3,', maxindex=',i5,', flag=',i2)
      endif
  300 continue
      stop
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c            subroutine "make" to make a percolation cluster
c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      subroutine make2(p,nmax,nedge,isuc,isd,nsize)
c
c  *********************************************************************
c  *  This subroutine generates a percolation cluster on the square    *
c  *  lattice.                                                         *
c  *********************************************************************
c
      parameter (maxl=1000000,maxz=4)
      parameter (maxedge=2000)
      parameter (maxs=15000,maxnn=20000,maxsize=1000000)
      real ranf,p
      common /blk0/ nn(maxnn),list(2,maxsize)
      dimension latt(maxedge,maxedge),label(0:maxl)
      dimension ndr(2,maxz)
      data ndr/0,1,1,0,0,-1,-1,0/
c
c
c
      do 6 i=0,maxl
    6 label(i)=0
c
c==============================================================================
c   Start generating clusters in nedge x nedge grid.
c==============================================================================
c
c   Generate the first column from top to bottom.
c------------------------------------------------------------------------------
c   In this program, column refers to (i,.) for fixed i even though this is
c   actually the opposite of the order of matrix indices.  This is simply
c   because it is easier to interpret indices as (x,y).
c------------------------------------------------------------------------------
c
      ic=0
      if(ranf(isd).le.p) then
         ic=ic+1
         if(ic.gt.maxl) then
            isuc=-2
            return
         endif
         latt(1,1)=ic
         label(ic)=1
      else
         latt(1,1)=0
      endif
      do 10 j=2,nedge
      if(ranf(isd).le.p) then
         if(latt(1,j-1).gt.0) then
            latt(1,j)=latt(1,j-1)
            label(latt(1,j))=label(latt(1,j))+1
         else
            ic=ic+1
            if(ic.gt.maxl) then
               isuc=-2
               return
            endif
            latt(1,j)=ic
            label(ic)=1
         endif
      else
         latt(1,j)=0
      endif
   10 continue
c
c   Generate column 2 through column "nedge"; each one from top to bottom.
c
      do 20 i=2,nedge
c
c   The top site of each column.
c
      if(ranf(isd).le.p) then
         if(latt(i-1,1).gt.0) then
            latt(i,1)=latt(i-1,1)
            ip=latt(i,1)
            if(label(ip).lt.0) ip=-label(ip)
            label(ip)=label(ip)+1
         else
            ic=ic+1
            if(ic.gt.maxl) then
               isuc=-2
               return
            endif
            latt(i,1)=ic
            label(ic)=1
         endif
      else
         latt(i,1)=0
      endif
c
c   The rest of the column.
c
      do 30 j=2,nedge
      if(ranf(isd).le.p) then
         imx=max0(latt(i-1,j),latt(i,j-1))
         if(float(latt(i,j-1))*float(latt(i-1,j)).gt.0
     c   .and. latt(i,j-1).ne.latt(i-1,j)) then
            latt(i,j)=min0(latt(i-1,j),latt(i,j-1))
            ip=imx
            if(label(ip).lt.0) ip=-label(ip)
            iq=latt(i,j)
            if(label(iq).lt.0) iq=-label(iq)
            if(ip.eq.iq) then
               label(ip)=label(ip)+1
            else
               jmin=min0(ip,iq)
               jmax=max0(ip,iq)
               label(jmin)=label(jmin)+label(jmax)+1
               label(jmax)=-jmin
               do 13 l=jmax+1,ic
   13          if(label(l).eq.-jmax) label(l)=-jmin
            endif
         elseif(latt(i,j-1)+latt(i-1,j).gt.0) then
            latt(i,j)=imx
            ip=imx
            if(label(ip).lt.0) ip=-label(ip)
            label(ip)=label(ip)+1
         else
            ic=ic+1
            if(ic.gt.maxl) then
               isuc=-2
               return
            endif
            latt(i,j)=ic
            label(ic)=1
         endif
      else
         latt(i,j)=0
      endif
   30 continue
c
c
c
   20 continue
c
c==============================================================================
c   Finished generation.  Find the cluster number of the largest cluster.
c==============================================================================
c
      jmax=0
      numb=0
      id=0
      do 24 i=1,ic
      if(label(i).le.0) go to 24
      numb=numb+label(i)
      if(label(i).le.jmax) go to 24
      jmax=label(i)
      id=i
   24 continue
c     
c==============================================================================
c   Set up list array for the maximal cluster.
c==============================================================================
c
c   First assign site number to each site of the cluster
c
      nsize=0
      do 70 i=1,nedge
      do 70 j=1,nedge
      if(label(latt(i,j)).eq.-id.or.latt(i,j).eq.id) then
         nsize=nsize+1
         latt(i,j)=nsize
      else
         latt(i,j)=0
      endif
   70 continue
c
c
c
      do 60 i=1,nedge
      do 60 j=1,nedge
      if(latt(i,j).ne.0) then
         n=latt(i,j)
         list(1,n)=i
         list(2,n)=j	
      endif
   60 continue
c
c==============================================================================
c   Finished setting up the output array.
c==============================================================================
c==============================================================================
c   Set up "nn".
c==============================================================================
c
      ns1=maxs
      do 40 j=1,maxnn
   40 nn(j)=0
      do 50 i=1,ic
         ns=label(i)
         if(ns.gt.maxs) then
            ns1=ns1+1
            nn(ns1)=ns
         elseif(ns.gt.0) then
            nn(ns)=nn(ns)+1
         endif
   50 continue
      nmax=ns1
c
c==============================================================================
c   Finished setting up the output arrays.
c==============================================================================
c
      isuc=1
      return
      end
c
c
c
      subroutine dist(nmax,xiav,px,nc,ng)
      parameter (maxs=15000,maxnn=20000)
      dimension nna(8),nnb(8),nr(9),xiav(2)
      common /blk0/ nn(maxnn)
      do 10 i=1,8
      nna(i)=-1
   10 nnb(i)=-1
    8 ng=0
      nc=0
      js=0
      xx=0
      do 816  i=1,maxs
      if(nn(i).eq.0)go to 816
      js=js+1
      nna(js)=i
      nc=nc+nn(i)
      ng=ng+i*nn(i)
      xx=xx+float(i)**2*nn(i)
      nnb(js)=nn(i)
      if (js.ne.8)go to816
  815 write(2,50)(nna(j),nnb(j),j=1,8)
      write(6,50)(nna(j),nnb(j),j=1,8)
   50 format(8(i8,i8))
      js=0
  816 continue
      if(nmax.le.maxs)go to 62
      do 61 i=(maxs+1),nmax
      js=js+1
   31 nnb(js)=1
      nna(js)=nn(i)
      xx=xx+float(nn(i))**2
      nc=nc+1
      ng=ng+nn(i)
      if(js.ne.8)go to 61
      write(2,50)(nna(j),nnb(j),j=1,8)
      write(6,50)(nna(j),nnb(j),j=1,8)
      js=0
   61 continue
   62 if(js.eq.0)go to 63
      js1=js+1
      do 814 kk=js1,8
  814 nna(kk)=-1
      write(2,50)(nna(j),nnb(j),j=1,8)
      write(6,50)(nna(j),nnb(j),j=1,8)
   63 if(js.eq.0)js=8
      px=nna(js)
      xiav(1)=0
      xiav(2)=0
      if(ng.eq.0) go to 77
      xiav(1)=xx/float(ng)
      xiav(2)=(xx-px*px)/float(ng)
      px=px/float(ng)
   77 nna(1)=-1
      write(2,50)(nna(j),nnb(j),j=1,8)
      write(6,50)(nna(j),nnb(j),j=1,8)
      return
      end
c
c
c
      real function ranf(iseed)
      iseed=iseed*1566083941+1
      ranf=float(iseed)*2.328306e-10+0.5
      return
      end