      PROGRAM site
*     draw site percolation configurations           
      IMPLICIT NONE
      REAL*8 delta_p, p, rsite(64,64)
      INTEGER IR_, L, occup
      CALL GWopen(IR_, 0)
*     define lattice and screen parameters
      CALL initial(rsite,L,p,delta_p, occup)
*     occupy sites for given probability p
      CALL configuration(rsite,L,p,delta_p, occup)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(rsite,L,p,delta_p,occup)
      IMPLICIT NONE
      REAL*8 rsite(64,64), p, delta_p
      REAL r, xwin, ywin, pix, piy
      REAL dmy, rnd
      INTEGER IR_, x, y, L, occup
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'linear dimension of square lattice = '
      READ(*,*) L
      WRITE(*,'(A,$)') 'initial probability p = '
      READ(*,*) p
*     desired increase in p upon mouse click
      WRITE(*,'(A,$)') 'delta p = '
      READ(*,*) delta_p
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, -0.1*xwin, -0.1*ywin, xwin, ywin)
*     represent occupied site as blue circle with unit diameter
      CALL GWsetpen(IR_, 16, -1, -1, 4)
      CALL GWpxlsiz(IR_, pix, piy)
      CALL GWcmbmrk(occup, 1, 6, 1 + 2*piy, 16, -1, 4)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      r = 0.5D0
      CALL GWrect(IR_, 1-r, 1-r, L+r, L+r)
      DO y = 1, L
         DO x = 1, L
*           assign random number to each lattice site
            rsite(x,y) = rnd(0)
            CALL GWsetpxl(IR_, REAL(x), REAL(y), 0)
         END DO
      END DO
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE configuration(rsite,L,p,delta_p, occup)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 rsite(64,64), p, delta_p
      REAL xs, ys
      INTEGER IR_, L, x, y, occup
      LOGICAL Ltmp1_
      DO WHILE(p .LE. 1)
         CALL GWsetogn(IR_, 0,-10)
         WRITE(Stmp1_,*) 'p =', p
         CALL GWputtxt(IR_, 0.0, REAL(L+1), Stmp1_)
        ! represent occupied site as blue circle with unit diameter
         DO y = 1, L
            DO x = 1, L
               IF(rsite(x,y) .LT. p) THEN
                  CALL GWputcmb(IR_, occup, REAL(x), REAL(y), 
     +                          0.0, 0.0, 0)
               END IF
            END DO
         END DO
         CALL GWflush(IR_, -10)
         CALL GWsetogn(IR_, 0, 2)
         CALL GWsetbrs(IR_, 13, 1, -1)
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            CALL GWcappnt(IR_, xs, ys, '')
            IF(xs .LT. 0.0) RETURN
*           click on occupied site to see cluster
*           click mouse to right of lattice to increase p by delta_p
            IF(xs .LT. L+0.5 .AND. xs .GT.-0.5
     +       .AND. ys .LT. L+0.5 .AND. ys .GT. -0.5) THEN
               CALL GWflood(IR_, REAL(xs), REAL(ys))
            END IF
            Ltmp1_ = xs .LT. L+0.5
         END DO
         p = p + delta_p
         CALL GWerase(IR_, 2, -1)
      END DO
      END
