      PROGRAM rg
      IMPLICIT NONE
      REAL*8 r(32,32)
      INTEGER L, b, IC1_, IC2_, IC3_, IC4_, IR_
      DATA IC1_/1/
      DATA IC2_/2/
      DATA IC3_/3/
      DATA IC4_/4/
      CALL GWopen(IR_, 0)
      CALL initial(r,L,b,IC1_,IC2_,IC3_,IC4_)
      CALL configuration(r,L,b,IC1_,IC2_,IC3_,IC4_)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(r,L,b,IC1_,IC2_,IC3_,IC4_)
      IMPLICIT NONE
      REAL*8 r(32,32)
      REAL GWaspect
      REAL dmy, rnd
      INTEGER IR_, x, y, L, b, IC1_, IC2_, IC3_, IC4_
      dmy = rnd(-1)
      L = 32
*     program must be modified for b <> 2
      b = 2
*     assign random number to each site
      DO y = 1, L
         DO x = 1, L
            r(x,y) = rnd(0)
         END DO
      END DO
      CALL GWvport(IR_, 0.0, 0.5, 0.5*GWaspect(1), 1.0)
      CALL lattice(L)                       ! draw original lattice in window #1
      CALL GWsavevp(IR_, IC1_)
      CALL GWvport(IR_, 0.5*GWaspect(1), 0.5, GWaspect(1), 1.0)
      CALL lattice(L/b)
      CALL GWsavevp(IR_, IC2_)
      CALL GWvport(IR_, 0.0, 0.0, 0.5*GWaspect(1), 0.5)
      CALL lattice(L/b**2)
      CALL GWsavevp(IR_, IC3_)
      CALL GWvport(IR_, REAL(0.5*GWaspect(1)), 0.0, GWaspect(1), 0.5)
      CALL lattice(L/b**3)
      CALL GWsavevp(IR_, IC4_)
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE lattice(Ln)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL xwin, ywin
      INTEGER IR_, x, y, Ln
      CALL compute_aspect_ratio(REAL(Ln),xwin,ywin)
      CALL GWindow(IR_, -0.07*xwin, -0.07*ywin, 1.07*xwin, 1.07*ywin)
      CALL GWrect(IR_, 0.0, 0.0, REAL(Ln), REAL(Ln))
      WRITE(Stmp1_,*) 'L =', Ln
      CALL GWputtxt(IR_, 0.0, ywin, Stmp1_)
      DO y = 1, Ln
         DO x = 1, Ln
            CALL GWsrect(IR_, REAL(x - 0.45), REAL(y - 0.45), 
     +                        REAL(x - 0.55), REAL(y - 0.55), 16)
         END DO
      END DO
      END
C
      SUBROUTINE configuration(r,L,b,IC1_,IC2_,IC3_,IC4_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 p, s(32,32), s1(32,16), s2(32,8), s3(32,4), r(32,32)
      INTEGER IR_, L, b, x, y, IC1_, IC2_, IC3_, IC4_
      DATA p/0.0D0/
      DO WHILE(p .LE.1)
         CALL GWselvp(IR_, IC1_)
         CALL GWinput(IR_, 'p = ', Stmp1_)
         IF(IR_ .EQ. 0) RETURN
         READ(Stmp1_, *) p
         CALL GWsetogn(IR_, 0, -IC1_)
         WRITE(Stmp1_, *) 'p =', p
         CALL GWputtxt(IR_, L*1.1, L*1.1, Stmp1_)
         DO y = 1, L
            DO x = 1, L
               IF(r(x,y) .LT. p) THEN
                  CALL GWsrect(IR_, REAL(x - 1), REAL(y - 1), 
     +                              REAL(x), REAL(y), 13)
                  s(x,y) = 1
               ELSE
                  s(x,y) = 0
               END IF
            END DO
         END DO
         CALL GWflush(IR_, -IC1_)           ! erase previous value of p
         CALL block(IC2_,L/b,s,s1)
         CALL block(IC3_,L/b**2,s1,s2)
         CALL block(IC4_,L/b**3,s2,s3)
         CALL GWsetogn(IR_, 0, 1)
      END DO
      END
C
      SUBROUTINE block(IC9_,Ln,w,wr)
      IMPLICIT NONE
      REAL*8 w(32,*), wr(32,*)
      INTEGER IR_, Ln, x, xc, y, yc, IC9_
      CALL GWselvp(IR_, IC9_)
      CALL GWsetogn(IR_, 0, -IC9_)
      DO y = 1, Ln
         yc = 2*y - 1
         DO x = 1, Ln
            xc = 2*x - 1
*           cell spans vertically -> renormalized site occupied
            IF(w(xc,yc)*w(xc,yc+1) .EQ. 1 .OR. 
     +         w(xc+1,yc)*w(xc+1,yc+1) .EQ. 1) THEN
               wr(x,y) = 1
               CALL GWsrect(IR_, REAL(x - 1), REAL(y - 1), 
     +                           REAL(x), REAL(y), 13)
            ELSE
               wr(x,y) = 0
            END IF
         END DO
      END DO
      CALL GWflush(IR_, -IC9_)
      END
