      SUBROUTINE assign(site,np,L)
*     assign cluster numbers to occupied sites
      IMPLICIT NONE
      INTEGER down, left, ncluster, x, y, site(0:64,0:64), np(0:1000), 
     +        L, proper, Itmp1_
      DO Itmp1_ = 0, 1000
         np(Itmp1_) = 0
      END DO
      ncluster = 0                          ! cluster number
      DO y = 1, L
         DO x = 1, L
            IF(site(x,y) .LT. 0) THEN
               down = y - 1                 ! square lattice
               left = x - 1
               IF(site(x,down) + site(left,y) .EQ. 0) THEN
                  ncluster = ncluster + 1   ! new cluster
                  site(x,y) = ncluster
                  np(ncluster) = ncluster   ! proper label
               ELSE
                  CALL neighbor(site,np,x,y)
               END IF
            END IF
         END DO
      END DO
*     assign proper labels to cluster array
      DO y = 1, L
         DO x = 1, L
            site(x,y) = proper(np,site(x,y))
         END DO
      END DO
      END
C
      SUBROUTINE neighbor(site,np,x,y)
      IMPLICIT NONE
      INTEGER down, left, site(0:64,0:64), np(0:1000), x, y
      down = y - 1
      left = x - 1
      IF(site(x,down)*site(left,y) .GT. 0) THEN
*        both neighbors occupied
         CALL label_min(site,np,x,y,left,down)
      ELSE IF(site(x,down) .GT. 0) THEN
         site(x,y) = site(x,down)
      ELSE
         site(x,y) = site(left,y)
      END IF
      END
C
      SUBROUTINE label_min(site,np,x,y,left,down)
*     both neighbors occupied, determine minimum cluster number
      IMPLICIT NONE
      INTEGER cl_down, cl_left, nmin, nmax, site(0:64,0:64), 
     +        np(0:1000), x, y, left, down, proper
      IF(site(left,y) .EQ. site(x,down)) THEN
*        both neighbors have same cluster label
         site(x,y) = site(left,y)
      ELSE
*        determine minimum cluster label
         cl_left = proper(np,site(left,y))
         cl_down = proper(np,site(x,down))
         nmax = max(cl_left,cl_down)
         nmin = min(cl_left,cl_down)
         site(x,y) = nmin
         IF(nmin .NE. nmax) THEN
            np(nmax) = nmin                 ! set improper label nmax = nmin
         END IF
      END IF
      END
