      PROGRAM site2
*     draw site percolation configurations           
      IMPLICIT NONE
      REAL*8 delta_p, p, rsite(64,64)
      INTEGER IR_, L, site(0:64,0:64), np(0:1000), mrk
      DATA site/4225*0/, p/0.0D0/
      CALL GWopen(IR_, 0)
*     define lattice and screen parameters
      CALL initial(rsite,L,p,delta_p,mrk)
*     occupy sites for given probability p
      CALL configuration(rsite,L,p,delta_p,site,np,mrk)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(rsite,L,p,delta_p,mrk)
      IMPLICIT NONE
      CHARACTER Stmp1_*80
      REAL*8 rsite(64,64), p, delta_p
      REAL xwin, ywin
      REAL dmy, rnd
      INTEGER IR_, mrk, x, y, L
      dmy = rnd(-1)
      CALL GWinput(IR_, 'linear dimension of square lattice = ', Stmp1_)
      READ(Stmp1_,*) L
      CALL GWinput(IR_, 'initial probability p = ', Stmp1_)
      READ(Stmp1_,*) p
*     desired increase in p upon mouse click
      CALL GWinput(IR_, 'delta p = ', Stmp1_)
      READ(Stmp1_,*) delta_p
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, -0.1*xwin, -0.1*ywin, 1.1*xwin, 1.1*ywin)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWsetbrs(IR_, -1, 0, -1)
      CALL GWcmbmrk(mrk, 0, 5, 1.0, 15, -1, 4)
      CALL GWanchor(IR_, 1);
      DO y = 1, L
         DO x = 1, L
*           assign random number to each lattice site
            rsite(x,y) = rnd(0)
         END DO
      END DO
      END
C
      SUBROUTINE configuration(rsite,L,p,delta_p,site,np,mrk)
      IMPLICIT NONE
      CHARACTER Stmp1_*80
      REAL*8 rsite(64,64), p, delta_p
      REAL xs, ys
      INTEGER IR_, L, vspan, x, y, site(0:64,0:64), np(0:1000), mrk
      CALL GWsetogn(IR_, 0, -2);
      DO WHILE(p .LE. 1)
         WRITE(Stmp1_,*) 'p =', p
         CALL GWsettxt(IR_, L/25.0, 0.0, 1, 0, -1, '')
         CALL GWputtxt(IR_, 0.0, REAL(L+1), Stmp1_)
         CALL GWsettxt(IR_, 0.6, 0.0, 0, 0, -1, '')
         DO y = 1, L
            DO x = 1, L
               IF(rsite(x,y) .LT. p) THEN
                  site(x,y) = -1
                  CALL GWputcmb(IR_, mrk, REAL(x), REAL(y), 0.0,0.0,0)
               END IF
            END DO
         END DO
         CALL GWrect(IR_, 0.5, 0.5, L+0.5, L+0.5)
                                            ! redraw lines
         CALL assign(site, np, L)
         CALL show_label(site, np, L, 'yes')
         CALL span(site, vspan, L)
         CALL GWflush(IR_, -2)
         CALL GWcappnt(IR_, xs, ys, '')
         IF(xs .LT. 0.0) RETURN
         p = p + delta_p
      END DO
      END
