PROGRAM invasion
*     generate invasion percolation cluster
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER IR_, Lx, Ly, perx(5000), pery(5000), nper
      CALL GWopen(IR_, 0)
      CALL initial(Lx,Ly)
      CALL assign(site,perx,pery,nper,Lx,Ly)
      CALL invade(site,perx,pery,nper,Lx,Ly)
      CALL average(site,Lx,Ly)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(Lx,Ly)
      IMPLICIT NONE
      REAL xwin, ywin, GWaspect
      REAL dmy, rnd
      INTEGER IR_, x, y, Lx, Ly
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'length in y direction = '
      READ(*,*) Ly
      Lx = 2*Ly
      CALL GWvport(IR_, 0.01*GWaspect(1), 0.2, 0.65*GWaspect(1), 0.99)
      CALL compute_aspect_ratio(REAL(Lx),xwin,ywin)
      CALL GWindow(IR_, -0.5, -0.5, xwin-0.5, ywin-0.5)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWrect(IR_, 1.0, 1.0, Lx+1.0, Ly+1.0)
      DO y = 1, Ly
         DO x = 1, Lx
            CALL GWsetpxl(IR_, REAL(x)+0.5, REAL(y)+0.5, 0)
         END DO
      END DO
      END
C
      SUBROUTINE assign(site,perx,pery,nper,Lx,Ly)
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      REAL rnd
      INTEGER IR_, x, y, nper, perx(5000), pery(5000), Lx, Ly
      DO y = 1, Ly
         site(1,y) = 1                      ! occupy first column
         CALL GWsrect(IR_, 1.0, REAL(y), 
     +                     2.0, REAL(y)+1.0, 16)
      END DO
*     assign random numbers to remaining sites
      DO y = 1, Ly
         DO x = 2, Lx
            site(x,y) = rnd(0)
         END DO
      END DO
*     sites in second column are initial perimeter sites
*     site(x,y) greater than 2 if perimeter site
      x = 2
      nper = 0
      DO y = 1, Ly
         site(x,y) = 2 + site(x,y)
         nper = nper + 1                    ! number of perimeter sites
*        sort perimeter sites
*        order perimeter list
         CALL insert(site,perx,pery,nper,x,y)
      END DO
      END
C
      SUBROUTINE insert(site,perx,pery,nper,x,y)
*     call linear or binary sort
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER ilist, ninsert, nper, x, y, perx(5000), pery(5000)
      CALL binary_search(site,perx,pery,nper,x,y,ninsert)
*     move sites with smaller random numbers to next higher array index
      DO ilist = nper, ninsert + 1, - 1
         perx(ilist) = perx(ilist-1)
         pery(ilist) = pery(ilist-1)
      END DO
      perx(ninsert) = x                     ! new site inserted in list
      pery(ninsert) = y
      END
C
      SUBROUTINE binary_search(site,perx,pery,nper,x,y,ninsert)
*     divide list in half and determine half in which random # belongs
*     continue this division until position of number is determined
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER nmid, xmid, ymid, nper, x, y, ninsert, nfirst, nlast, 
     +        perx(5000), pery(5000)
      LOGICAL Ltmp1_
      nfirst = 1                            ! beginning of list
      nlast = nper - 1                      ! end of list
      IF(nlast .LT. 1) nlast = 1
      nmid = (nfirst + nlast)/2             ! middle of list
*     determine which half of list new number is located
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         IF(nlast - nfirst .LE. 1) THEN
            ninsert = nlast
            RETURN
         END IF
         xmid = perx(nmid)
         ymid = pery(nmid)
         IF(site(x,y) .GT. site(xmid,ymid)) THEN
            nlast = nmid                    ! search upper half
         ELSE
            nfirst = nmid                   ! search lower half
         END IF
         nmid = (nfirst + nlast)/2
      END DO
      END
C
      SUBROUTINE linear_search(site,perx,pery,nper,x,y,ninsert)
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER iper, xperim, yperim, nper, x, y, ninsert,
     +        perx(5000), pery(5000)
      IF(nper .EQ. 1) THEN
         ninsert = 1
      ELSE
         DO iper = 1, nper - 1
            xperim = perx(iper)
            yperim = pery(iper)
            IF(site(x,y) .GT. site(xperim,yperim)) THEN
               ninsert = iper               ! insert new site
               RETURN
            END IF
         END DO
      END IF
      ninsert = nper
      END
C
      SUBROUTINE invade(site,perx,pery,nper,Lx,Ly)
*     nx and ny are components of vectors pointing to nearest neighbors
      IMPLICIT NONE
      REAL*8 site(0:200,0:100)
      INTEGER IR_, nx(4), ny(4), perx(5000), pery(5000),
     +       Lx, Ly, i, x, xnew, y, ynew, nper, DP_, DATA_(8)
      LOGICAL Ltmp1_
      DATA DP_/1/, DATA_/1,0,-1,0,0,1,0,-1/
      DO i = 1, 4
         nx(i) = DATA_(DP_)
         DP_ = DP_ + 1
         ny(i) = DATA_(DP_)
         DP_ = DP_ + 1
      END DO
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         x = perx(nper)
         y = pery(nper)
         nper = nper - 1
*        mark site occupied and no longer perimeter site
         site(x,y) = site(x,y) - 1
         CALL GWsrect(IR_, REAL(x), REAL(y), 
     +                     REAL(x)+1.0, REAL(y)+1.0, 16)
         DO i = 1, 4
            xnew = x + nx(i)
            ynew = y + ny(i)
            IF(ynew .GT. Ly) THEN
               ynew = 1
            ELSE IF(ynew .LT. 1) THEN
               ynew = Ly
            END IF
            IF(site(xnew,ynew) .LT. 1) THEN
               site(xnew,ynew) = site(xnew,ynew) + 2
               nper = nper + 1
               CALL insert(site,perx,pery,nper,xnew,ynew)
            END IF
         END DO
         Ltmp1_ = .NOT.(x .GE. Lx)
      END DO
      END
C
      SUBROUTINE average(site,Lx,Ly)
*     compute probability density P(r)
      IMPLICIT NONE
      REAL*8 site(0:200,0:100), dr, r, P(0:20)
      INTEGER Lmax, Lmin, n, nbin, nr(0:20), occupied,
     +       ibin, x, y, Lx, Ly
      DATA occupied/0/, P/21*0.0D0/, nr/21*0/
      Lmin = Lx/3
      Lmax = 2*Lmin
      n = (Lmax - Lmin + 1)*Ly              ! # sites in middle half of lattice
      dr = 0.05D0
      nbin = 1/dr
      DO x = Lmin, Lmax
         DO y = 1, Ly
            ibin = nbin*(mod(site(x,y),1))
            nr(ibin) = nr(ibin) + 1
            IF(site(x,y) .GE. 1 .AND. site(x,y) .LT. 2) THEN
               occupied = occupied + 1      ! total # of occupied sites
               P(ibin) = P(ibin) + 1
            END IF
         END DO
      END DO
      WRITE(*,*) '# occupied sites =', occupied
      WRITE(*,'(2A13)') ' r', '  P(r)'
      WRITE(*,*)
      DO ibin = 0, nbin
         r = dr*ibin
         IF(nr(ibin) .GT. 0) WRITE(*,'(2F13.6)') r, P(ibin)/nr(ibin)
      END DO
      WRITE(*,*)
      END