program main
      implicit none

      integer MAX_X, MAX_Y, STEPS_MAX, UNCHANGED_MAX
      parameter (MAX_X = 500, MAX_Y = 500)
      parameter (STEPS_MAX = 200, UNCHANGED_MAX = 10)

      integer grid(0:MAX_X+1,0:MAX_Y+1)
      integer nx, ny
      integer max_steps, max_unchanged
      integer vegies
      integer nsteps, nsims
      integer ndied, nunsettled, nstable
      real*8  tot_steps_stable, tot_veg_stable
      real*8  prob
      integer seed, seed0
      integer i
      integer game_o_life

 10   continue
      write (6,*) ' Enter X, Y extent of wilderness: '
      read (5,*) nx, ny
      if (nx .gt. MAX_X .or. nx .gt. MAX_Y) goto 10
      
      write (6,*) ' Enter population probability: '
      read (5,*) prob
      
      write (6,*) 'Enter number of simulations: '
      read (5,*) nsims
      
      write (6, *) 'Enter random number seed: '
      read (5, *) seed0
      
      ndied = 0
      nunsettled = 0
      nstable = 0
      tot_steps_stable = 0
      tot_veg_stable = 0

      do i = 1,nsims

c Initialize the grid values using the given probability.

        seed = seed0*i
        call initialize_grid(grid, nx, ny, seed, prob)
        
c Now run the game of life simulation, returning number of steps.

        max_steps = STEPS_MAX
        max_unchanged = UNCHANGED_MAX
        nsteps = game_o_life(grid, nx, ny, max_steps, max_unchanged,
     $       vegies)

        write (6,*) ' time steps, total veggies =', nsteps, vegies

        if (vegies .eq. 0) then
          ndied = ndied + 1
        else if (nsteps .ge. max_steps) then
          nunsettled = nunsettled + 1
        else
          nstable = nstable + 1
          tot_steps_stable = tot_steps_stable + nsteps
          tot_veg_stable = tot_veg_stable + vegies
        endif
      enddo

      if (nstable .gt. 0) then
        tot_steps_stable = tot_steps_stable / nstable
        tot_veg_stable = tot_veg_stable / nstable
      endif

      write (6,*)' Percentage which died out:',100.0*ndied/nsims
      write (6,*)' Percentage unsettled:     ',100.0*nunsettled/nsims
      write (6,*)' Percentage stablized:     ',100.0*nstable/nsims
      write (6,*)'   Of which:'
      write (6,*)'   Average steps:          ',tot_steps_stable
      write (6,*)'   Average vegetation:     ',tot_veg_stable

      end


      subroutine initialize_grid(grid, nx, ny, seed, prob)
      implicit none
      integer MAX_X, MAX_Y
      parameter (MAX_X = 500, MAX_Y = 500)
      integer grid(0:MAX_X+1,0:MAX_Y+1)
      integer nx, ny
      integer seed
      real*8  prob
      integer i,j
      integer index, new_seed
      real*8 rand1

      do i = 1,nx
	do j = 1,ny
	  index = ny*i + j
          new_seed = seed + index
	  if (rand1(new_seed) .gt. prob) then
            grid(i,j) = 0
          else
            grid(i,j) = 1
          endif
        enddo
      enddo

      return
      end


      integer function game_o_life(grid, nx, ny, max_steps,
     $     max_unchanged, vegies)
      implicit none
      integer MAX_X, MAX_Y, STEPS_MAX, UNCHANGED_MAX
      parameter (MAX_X = 500, MAX_Y = 500)
      parameter (STEPS_MAX = 200, UNCHANGED_MAX = 10)
      integer grid(0:MAX_X+1,0:MAX_Y+1)
      integer nx, ny
      integer max_steps, max_unchanged
      integer step
      integer converged
      integer n_unchanged
      integer vegies, old_vegies, old2_vegies, old3_vegies
      integer neighbors
      integer temp_grid(0:MAX_X+1,0:MAX_Y+1)
      integer i,j

      step = 1
      vegies = 1
      old_vegies = -1
      old2_vegies = -1
      old3_vegies = -1
      n_unchanged = 0
      converged = 0

 10   if (converged .ne. 1 .and. vegies .gt. 0 .and.
     $     step .lt. max_steps) then

c Count the total amount of vegetation. 

        vegies = 0
        do i = 1,nx
          do j = 1,ny
            vegies = vegies + grid(i,j)
          enddo
        enddo
        
        if (vegies .eq. old_vegies .or.
     $       vegies .eq. old2_vegies .or.
     $       vegies .eq. old3_vegies) then
          n_unchanged = n_unchanged + 1
          if (n_unchanged .ge. max_unchanged) converged = 1
        else
          n_unchanged = 0
        endif

        old3_vegies = old2_vegies
        old2_vegies = old_vegies
        old_vegies = vegies
        
c       write(*,*) ' step, vegies = ', step, vegies

        if (converged .ne. 1) then

c  Copy the sides of the grid to make torus simple.

          do i = 1,nx
            grid(i,0) = grid(i,ny)
            grid(i,ny+1) = grid(i,1)
          enddo

          do j = 0,ny+1
            grid(0,j) = grid(nx,j)
            grid(nx+1,j) = grid(1,j)
          enddo

c Now run one time step, putting result in temp_grid.
          
          do i = 1,nx
            do j = 1,ny
              neighbors = grid(i-1,j-1) + grid(i-1,j) +
     $             grid(i-1,j+1) + grid(i,j-1) + grid(i,j+1) +
     $             grid(i+1,j-1) + grid(i+1,j) + grid(i+1,j+1)
              temp_grid(i,j) = grid(i,j)
              if (neighbors .ge. 25 .or. neighbors .le. 3) then
                temp_grid(i,j) = temp_grid(i,j) - 1
                if (temp_grid(i,j) .lt. 0) temp_grid(i,j) = 0
              else if (neighbors .le. 15) then
                temp_grid(i,j) = temp_grid(i,j) + 1
                if (temp_grid(i,j) .gt. 10) temp_grid(i,j) = 10
              endif
            enddo
          enddo
	  
c Now copy temp_grid back to grid.

          do i = 1,nx
            do j = 1,ny
              grid(i,j) = temp_grid(i,j)
            enddo
          enddo

          step = step + 1

        endif
        goto 10
      endif

      game_o_life = step
      
      return
      end


C Park/Miller RNG
 
      double precision function rand1(iseed)
      double precision aa,mm,sseed
      parameter (aa=16807.0D0,mm=2147483647.0D0)
      integer i
      
      do i = 1,5
        sseed = iseed
        sseed = mod(aa*sseed,mm)
        iseed = sseed
      enddo
      rand1 = sseed/mm

      return
      end