! modified from:
!http://www.physics.nau.edu/~hart/classes/550_spring_2004/hw_schedule/solutions/perc.f90

module modperc

implicit none
public :: perccheck

contains

function perccheck(grid, L) result(per)
integer :: per
integer, intent(in) :: grid(:,:), L
integer :: i
per = 0

do i = 1, maxval(grid)
   
!   if( any(grid(1,:)==i) .and. &
!       any(grid(L,:)==i) .and. &
!       any(grid(:,1)==i) .and. &
!       any(grid(:,L)==i) )then ! ith cluster touches all sides--it percolates
   
   if( any(grid(1,:)==i) .and. &
       any(grid(L,:)==i) ) then  ! ith cluster touches top-bot --it percolates

!    if( (any(grid(1,:)==i)   .and. &
!          any(grid(L,:)==i)) .or.  &
!        (any(grid(:,1)==i)   .and. &
!          any(grid(:,L)==i)) )then ! ith cluster touches top-bot or sx-dc it percolates

   per = i
      exit
   endif
enddo
end function perccheck

end module modperc

program percolation
use modperc
implicit  none

integer, parameter :: L=40  ! size of grid
integer, dimension(1:L+1,1:L+1) :: site = 0 ! grid of sites

real :: p, pmin, pmax, deltap ! occupation ratio
real :: rnd           ! random numbers
real :: invL          ! inverse box edge length

integer :: Pcount, loop, Nloop  ! probability counter and general loop counters
integer :: kx, ky               ! x, y positions (integer)
integer :: iclust = 1           ! current cluster number
integer :: imax, imin           ! cluster indexes for linking of crystal columns

call random_seed() ! randomly set the random seed
print*, "Enter the probability range pmin, pmax, with deltap=0.01"
read*, pmin, pmax

deltap = 0.01 !Probability step
Nloop  = 500  !Number of simulation with fixed probability
p = pmin
invL = 1/real(L)

open(unit=11,file="percavg.dat",status="replace",action="write")

do !Loop on probability
   pcount=0 

   do loop = 1, Nloop !Loop with fixed probability
      site=0
      iclust=1

      do kx=1,L  !
         do ky=1,L
            call random_number(rnd)
            if(rnd .lt. p) then
               site(kx,ky) = iclust
               iclust=iclust+1
               endif
         enddo
      enddo

      do kx=1,L !grouping percolated columns
         do ky=1,L
            if( site(kx,ky).ne.0 ) then
               if(site(kx,ky+1).ne.0) site(kx,ky+1) = site(kx,ky)
            end if
         enddo
      enddo

      do kx=1,L !"linking" percolated columns that touch themselves
         do ky=1,L
            if( (site(kx,ky).ne.0).and.(site(kx+1,ky).ne.0).and.(site(kx,ky).ne.site(kx+1,ky)) ) then
               imax=max(site(kx,ky),site(kx+1,ky))
               imin=min(site(kx,ky),site(kx+1,ky))
               where(site==imax) site=imin               
            end if
         enddo
      enddo
   
      if(perccheck(site,L).ne.0) pcount = pcount+1 ! if percolated the counter increase
   enddo
   
   write(11,'(3F12.8)') invL, p, float(pcount)/Nloop
   
   if(p.gt.pmax) exit
   p = p+deltap
   
enddo


close(11)
end program percolation