PROGRAM dla
      IMPLICIT NONE
      INTEGER L, site(-100:100,-100:100), IR_
      CALL GWopen(IR_, 0)
      CALL initial(site,L)
      CALL grow_cluster(site,L)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(site,L)
      IMPLICIT NONE
      REAL xwin, ywin
      REAL dmy, rnd
      INTEGER x, y, site(-100:100,-100:100), L, IR_, Itmp1_, Itmp2_
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'L = '
      READ(*,*) L
      CALL compute_aspect_ratio(REAL(L+2),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWrect(IR_, -L-0.5, -L-0.5, L+0.5, L+0.5)
      DO Itmp1_ = -100, 100
         DO Itmp2_ = -100, 100
            site(Itmp1_, Itmp2_) = 0
         END DO
      END DO
      DO y = -L, L
         DO x = -L, L
            CALL GWsetpxl(IR_, REAL(x), REAL(y), 0)
         END DO
      END DO
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWsetmrk(IR_, 6, 1.0, 13, -1, 4);
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE grow_cluster(site,L)
      IMPLICIT NONE
      REAL*8 pi, theta
      PARAMETER(pi = 3.141592653589793D0)
      REAL rnd
      INTEGER IR_, site(-100:100,-100:100), L, N, R0, x, y
      LOGICAL TBkeyinput
      R0 = 3                                ! start walker at distance R0 from origin
      site(0,0) = 1                         ! seed site
      N = 1                                 ! number of particles in cluster
      CALL GWputmrk(IR_, 0.0, 0.0)
      DO WHILE(.NOT.TBkeyinput())
*        find random initial position of new walker
         theta = 2*pi*rnd(0)
         x = int(R0*cos(theta))
         y = int(R0*sin(theta))
         CALL walk(site,L,x,y,R0,N)         ! random walk 
      END DO
      END
C
      SUBROUTINE walk(site,L,x,y,R0,N)
*     walk until on a perimeter site of cluster  
*     or walker strays too far from cluster
      IMPLICIT NONE
      REAL*8 r
      REAL rnd, random
      INTEGER step, x, y, R0, N, site(-100:100,-100:100), L
      CHARACTER onperimeter*80
      LOGICAL Ltmp1_
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         onperimeter = 'no'
         r = sqrt(DBLE(x*x + y*y))
         IF(r .LT. R0 + 1) THEN
*           test if walker on perimeter
            CALL test(site,L,x,y,r,R0,N,onperimeter)
         END IF
         step = int(r - R0) - 1             ! big step 
         IF(step .LT. 1) step = 1
         IF(onperimeter .EQ. 'no') THEN
            random = rnd(0)
            IF(random .LT. 0.25) THEN
               x = x + step
            ELSE IF(random .LT. 0.5) THEN
               x = x - step
            ELSE IF(random .LT. 0.75) THEN
               y = y + step
            ELSE
               y = y - step
            END IF
         END IF
         Ltmp1_ = .NOT.(onperimeter .EQ. 'yes' .OR. r .GT. 2*R0)
      END DO
      END
C
      SUBROUTINE test(site,L,x,y,r,R0,N,onperimeter)
      IMPLICIT NONE
      REAL*8 sum, r
      INTEGER IR_, site(-100:100,-100:100), L, x, y, R0, N
      CHARACTER*80 onperimeter, Stmp1_
      sum = site(x,y+1) + site(x,y-1) + site(x+1,y) + site(x-1,y)
      IF(sum .GT. 0) THEN                   ! walker on perimeter site 
         site(x,y) = 1
         onperimeter = 'yes'
         IF(abs(x) .LE. L .AND. abs(y) .LE. L) THEN
            CALL GWputmrk(IR_, REAL(x), REAL(y))
            N = N + 1
            CALL GWsetogn(IR_, 0, -2)
            WRITE(Stmp1_, *) 'N = ',N
            CALL GWputtxt(IR_, REAL(-L), REAL(L+1), Stmp1_);
            CALL GWflush(IR_, -2)
            CALL GWsetogn(IR_, 0, 1);
         ELSE
            CALL GWquit(IR_);               ! cluster outside box
            STOP
         END IF
         IF(r .GE. R0 - 1) R0 = int(r+2)
      END IF
      END