SUB assign(site(,),np(),L)
    ! assign cluster numbers to occupied sites
    DECLARE DEF proper
    MAT np = 0
    LET ncluster = 0              ! cluster number
    FOR y = 1 to L
        FOR x = 1 to L
            IF site(x,y) < 0 then      ! site occupied
               LET down = y - 1      ! square lattice
               LET left = x - 1
               IF site(x,down) + site(left,y) = 0 then
                  LET ncluster = ncluster + 1    ! new cluster
                  LET site(x,y) = ncluster
                  LET np(ncluster) = ncluster    ! proper label
               ELSE
                  CALL neighbor(site(,),np(),x,y)
               END IF
            END IF
        NEXT x
    NEXT y
    ! assign proper labels to cluster array
    FOR y = 1 to L
        FOR x = 1 to L
            LET site(x,y) = proper(np(),site(x,y))
        NEXT x
    NEXT y
END SUB

SUB neighbor(site(,),np(),x,y)    ! determine occupancy of neighbors 
    LET down = y - 1
    LET left = x - 1
    IF site(x,down)*site(left,y) > 0 then
       ! both neighbors occupied
       CALL label_min(site(,),np(),x,y,left,down)
    ELSE IF site(x,down) > 0 then    ! down neighbor occupied
       LET site(x,y) = site(x,down)
    ELSE                          ! left neighbor occupied
       LET site(x,y) = site(left,y)
    END IF
END SUB

SUB label_min(site(,),np(),x,y,left,down)
    ! both neighbors occupied, determine minimum cluster number
    DECLARE DEF proper
    IF site(left,y) = site(x,down) then
       ! both neighbors have same cluster label
       LET site(x,y) = site(left,y)
    ELSE
       ! determine minimum cluster label
       LET cl_left = proper(np(),site(left,y))
       LET cl_down = proper(np(),site(x,down))
       LET nmax = max(cl_left,cl_down)
       LET nmin = min(cl_left,cl_down)
       LET site(x,y) = nmin
       IF nmin <> nmax then
          LET np(nmax) = nmin     ! set improper label nmax = nmin
       END IF
    END IF
END SUB

FUNCTION proper(np(),label)
    ! recursive function
    IF np(label) = label then
       LET proper = label
    ELSE
       LET proper = proper(np(),np(label))
    END IF
END DEF