PROGRAM cca
      IMPLICIT NONE
      INTEGER first_particle(1000), next_particle(1000), L, N, nx, ny,
     +       site, IR_, ncl, x, y
      COMMON /GLBLI/L, N, nx(4), ny(4), site(50,50), x(1000), y(1000)
      CALL GWopen(IR_, 0)
      CALL initial(first_particle,next_particle,ncl)
      DO WHILE(ncl .NE. 1)
         CALL move(next_particle,first_particle,ncl)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(first_particle,next_particle,ncl)
      IMPLICIT NONE
      CHARACTER buff*80
      REAL xwin, ywin
      REAL dmy, rnd
      INTEGER first_particle(1000), next_particle(1000), L,
     +       N, nx, ny, site, DATA_(8)
      INTEGER IR_, i, nn, xi, yi, ncl, x, y
      LOGICAL Ltmp1_
      COMMON /GLBLI/L, N, nx(4), ny(4), site(50,50), x(1000), y(1000)
      INTEGER DP_
      DATA DP_/1/,DATA_/1,0,-1,0,0,1,0,-1/
      dmy = rnd(-1)
      DO nn = 1, 4
         nx(nn) = DATA_(DP_)
         DP_ = DP_ + 1
         ny(nn) = DATA_(DP_)
         DP_ = DP_ + 1
      END DO
      WRITE(*,'(A,$)') 'L = '
      READ(*,*) L
      WRITE(*,'(A,$)') 'N = '
      READ(*,*) N
      CALL compute_aspect_ratio(REAL(L+1),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWrect(IR_, 1.0, 1.0, REAL(L+1), REAL(L+1))
      WRITE(buff, *) 'L = ', L, ', N = ', N
      CALL GWputtxt(IR_, 1.01, (L + 1)*1.01, buff)
      CALL GWanchor(IR_, 1)
      ncl = 0                               ! number of clusters
      DO i = 1, N
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            x(i) = int(L*rnd(0)) + 1
            y(i) = int(L*rnd(0)) + 1
            Ltmp1_ = site(x(i),y(i)) .NE. 0
         END DO
         ncl = ncl + 1
         site(x(i),y(i)) = ncl
         CALL GWsetogn(IR_, 0, (ncl+1));
         CALL GWsrect(IR_, REAL(x(i)),   REAL(y(i)), 
     +                     REAL(x(i)+1), REAL(y(i)+1), 16)
         first_particle(ncl) = i
         next_particle(i) = i
         xi = x(i)
         yi = y(i)
         CALL neighbor(xi,yi,next_particle,first_particle,ncl)
      END DO
      END
C
      SUBROUTINE neighbor(xi,yi,next_particle,first_particle,ncl)
      IMPLICIT NONE
      INTEGER next_particle(1000), first_particle(1000), L, N, nx, ny,
     +       site, pbc
      INTEGER nn, part, perim, px, py, xi, yi, ncl, x, y
      COMMON /GLBLI/L, N, nx(4), ny(4), site(50,50), x(1000), y(1000)
      DO nn = 1, 4
         px = pbc(xi + nx(nn),L)
         py = pbc(yi + ny(nn),L)
         perim = site(px,py)
         part = site(xi,yi)
         IF(perim .NE. 0 .AND. perim .NE. part) THEN
            CALL mergecl(perim,part,first_particle,next_particle,ncl)
         END IF
      END DO
      END
C
      SUBROUTINE mergecl(c1,c2,first_particle,next_particle,ncl)
      IMPLICIT NONE
      INTEGER p1next, plast, first_particle(1000), next_particle(1000),
     +       L, N, nx, ny, site, IR_
      INTEGER p, p1, p2, p2next, c1, c2, ncl, x, y
      LOGICAL Ltmp1_
      COMMON /GLBLI/L, N, nx(4), ny(4), site(50,50), x(1000), y(1000)
      p1 = first_particle(c1)
      p2 = first_particle(c2)
      p1next = next_particle(p1)
      p2next = next_particle(p2)
      next_particle(p1) = p2next
      next_particle(p2) = p1next
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         site(x(p2next),y(p2next)) = c1
         p2next = next_particle(p2next)
         Ltmp1_ = p2next .NE. p1next
      END DO
      CALL GWsetogn(IR_, (c2+1), (c1+1))
      plast = first_particle(ncl)
      IF(c2 .NE. ncl) THEN
         p = plast
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            site(x(p),y(p)) = c2
            p = next_particle(p)
            Ltmp1_ = .NOT.(p .EQ. plast)
         END DO
         first_particle(c2) = plast
         CALL GWsetogn(IR_, (ncl+1), (c2+1))
      END IF
      ncl = ncl - 1
      END
C
      SUBROUTINE move(next_particle,first_particle,ncl)
      IMPLICIT NONE
      INTEGER dx, dy, p1, next_particle(1000), first_particle(1000), L,
     +       N, nx, ny, site, pbc
      REAL rnd
      INTEGER IR_, c, direction, i, xi, yi, ncl, x, y
      LOGICAL Ltmp1_
      COMMON /GLBLI/L, N, nx(4), ny(4), site(50,50), x(1000), y(1000)
      c = int(ncl*rnd(0)) + 1
      direction = int(4*rnd(0)) + 1
      p1 = first_particle(c)
      i = p1
      dx = nx(direction)
      dy = ny(direction)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         site(x(i),y(i)) = 0
         x(i) = pbc(x(i) + dx,L)
         y(i) = pbc(y(i) + dy,L)
         i = next_particle(i)
         Ltmp1_ = i .NE. p1
      END DO
      CALL GWsetogn(IR_, 0, -(c+1))
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL GWsrect(IR_, REAL(x(i)),   REAL(y(i)), 
     +                     REAL(x(i)+1), REAL(y(i)+1), 16)
         site(x(i),y(i)) = c
         i = next_particle(i)
         Ltmp1_ = i .NE. p1
      END DO
      CALL GWflush(IR_, -(c+1))
      CALL GWsetogn(IR_, 0, 1)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         xi = x(i)
         yi = y(i)
         CALL neighbor(xi,yi,next_particle,first_particle,ncl)
         i = next_particle(i)
         Ltmp1_ = i .NE. p1
      END DO
      END
C
      FUNCTION pbc(s,L)
      IMPLICIT NONE
      INTEGER pbc, s, L
      IF(s .GT. L) THEN
         pbc = 1
      ELSE IF(s .LT. 1) THEN
         pbc = L
      ELSE
         pbc = s
      END IF
      END