PROGRAM eden
      IMPLICIT NONE
      INTEGER IC1_, IR_, L, status(-1:2), mass(1000)
      DATA IC1_/1/, mass/1000*0/
      CALL GWopen(IR_, 0)
      CALL initial(L,status,IC1_)
      CALL grow(L,mass,status,IC1_)
      CALL mass_dist(L,mass)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L,status,IC1_)
      IMPLICIT NONE
      REAL GWaspect, xwin, ywin
      REAL dmy, rnd
      INTEGER IR_, x, y, L, status(-1:2), IC1_
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'value of L (odd) = '
      READ(*,*) L
      CALL GWvport(IR_, 0.01*GWaspect(1), 0.01, 0.7*GWaspect(1), 0.99)
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsetpen(IR_, 13, -1, -1, 4)
      CALL GWcmbmrk(IR_, 1, 6, 1.0, 13, -1, 4)
      CALL GWsetpen(IR_, 16, -1, -1, 4)
      CALL GWcmbmrk(IR_, 2, 6, 1.0, 16, -1, 4)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWsetbrs(IR_, 0, 0, -1)
      CALL GWrect(IR_, 0.5, 0.5, L+0.5, L+0.5)
      DO y = 1, L
         DO x = 1, L
            CALL GWsetpxl(IR_, REAL(x), REAL(y), 0) ! tested, unoccupied site
         END DO
      END DO
      END
C
      SUBROUTINE grow(L,mass,status,IC1_)
*     generate single percolation cluster          
      IMPLICIT NONE
      REAL*8 dx, dy, site(0:201,0:201)
      REAL rnd
      INTEGER IR_, nx(4), ny(4), mass(1000), i, iper, nn, nper,
     +        perx(5000), pery(5000), r, x, DP_, DATA_(8),
     +        xnew, xseed, y, ynew, yseed, L, status(-1:2), IC1_
      CHARACTER boundary*80
      LOGICAL Ltmp1_
      DATA DP_/1/,DATA_/1,0,-1,0,0,1,0,-1/
      boundary = ''
*     set up boundary
      DO i = 1, L
         site(0,i) = -1
         site(L+1,i) = -1
         site(i,L+1) = -1
         site(i,0) = -1
      END DO
*     seed at center of lattice
      xseed = L/2 + 1
      yseed = xseed
      site(xseed,yseed) = 1                 ! seed site
      CALL GWputcmb(IR_, 1, REAL(xseed), REAL(yseed), 1.0, 1.0, 0)
      DO i = 1, 4
*        nx,ny direction vectors for new perimeter sites
         nx(i) = DATA_(DP_)
         DP_ = DP_ + 1
         ny(i) = DATA_(DP_)
         DP_ = DP_ + 1
*        perx,pery, positions of perimeter sites of seed
         perx(i) = xseed + nx(i)
         pery(i) = yseed + ny(i)
*        perimeter sites labeled by 2
         site(perx(i),pery(i)) = 2          ! site placed on perimeter list
         CALL GWputcmb(IR_, 2, REAL(perx(i)), REAL(pery(i)), 
     +                 1.0, 1.0, 0)
      END DO
      nper = 4                              ! initial number of perimeter sites
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         iper = int(rnd(0)*nper) + 1        ! randomly choose perimeter site
         x = perx(iper)                     ! coordinate of perimeter site
         y = pery(iper)
*        last perimeter site in array replaces newly occupied site
         perx(iper) = perx(nper)
         pery(iper) = pery(nper)
         nper = nper - 1
         site(x,y) = 1
         CALL GWputcmb(IR_, 1, REAL(x), REAL(y), 1.0, 1.0, 0)
         dx = x - xseed
         dy = y - yseed
         r = sqrt(dx*dx + dy*dy)            ! distance from seed
         mass(r) = mass(r) + 1              ! number of sites at distance r
         DO nn = 1, 4
            xnew = x + nx(nn)
            ynew = y + ny(nn)
            IF(site(xnew,ynew) .EQ. 0) THEN
               nper = nper + 1
               perx(nper) = xnew
               pery(nper) = ynew
               site(xnew,ynew) = 2          ! site placed on perimeter list
               CALL GWputcmb(IR_, 2, REAL(perx(nn)), REAL(pery(nn)), 
     +                       1.0, 1.0, 0)
            END IF
         END DO
         IF(x .EQ. L .OR. x .EQ. 1) boundary = 'on'
         IF(y .EQ. L .OR. y .EQ. 1) boundary = 'on'
         Ltmp1_ = boundary .NE. 'on'
      END DO
      CALL GWrect(IR_, 0.5, 0.5, L+0.5, L+0.5)  ! redraw box
      END
C
      SUBROUTINE mass_dist(L,mass)
*     print mass as a function of r
      IMPLICIT NONE
      REAL GWaspect
      INTEGER IR_, masstotal, mass(1000), r, L
      CALL GWvport(IR_, 0.71*GWaspect(1), 0.01, GWaspect(1), 1.0)
      masstotal = 0
      WRITE(*,'(2A9)') 'r', 'M(r)'
      WRITE(*,*)
      DO r = 1, L
         IF(mass(r) .EQ. 0) RETURN
         masstotal = masstotal + mass(r)
         WRITE(*,'(2I9)') r, masstotal
      END DO
      END