*  Move particles (fish) in periodic box according to external force field (current)
*  Inputs :
*    d = grid dimension
*    fish = Indicates fish presence at the grid points
*    tsteps = number of time steps
*  Outputs:
*    plot of fish positions after each time step
*
*
      program fish3

      include '/usr/cm5/cmf/include/cm/timer-fort.h'
      include '/usr/cm5/cmf/include/cm/CMF_defs.h'

*     Define the size of the display window
      integer,  parameter :: m = 256
      integer, parameter :: disp_interval = 1

*     Max. Grid Size
      integer, parameter :: MAXG = 64 	

      integer 	fish(MAXG,MAXG), nbrcnt(MAXG,MAXG)
CMF$  layout    fish(:news,:news),nbrcnt(:news,:news)

      integer 	d, i, tsteps, stride, limit, off
      real      etime0, etime1

      call CMF_DESCRIBE_ARRAY(fish)
      print*, 'Input...'
      print*, 'Give grid size (max 64) : '
      read (*,*) d
      call cmf_randomize(100)
      call cmf_random(fish,2)
      tsteps = 200
*
* Initialize the timer
      call CM_timer_clear(0)
      call CM_timer_clear(1)
      call CM_timer_start(0)

      call disp_init(m, d, stride, limit, off)
*
* Loop over time steps
      do i=1, tsteps
          print *,'i = ',i
*
* Get neighbor count
          call CM_timer_start(1)
	  nbrcnt(2:d-1,2:d-1) = fish(1:d-2,2:d-1) + fish(1:d-2,1:d-2) +
     !                          fish(1:d-2,3:d) + fish(2:d-1,1:d-2) +
     !                          fish(2:d-1,3:d) + fish(3:d,2:d-1) +
     !                          fish(3:d,1:d-2) + fish(3:d,3:d)

*          print *,nbrcnt(2:d-1,2:d-1)
          where ((nbrcnt(2:d-1,2:d-1) .eq. 3) .or. ( (fish(2:d-1,2:d-1) 
     !		  .eq. 1) .and. (nbrcnt(2:d-1,2:d-1) .eq. 2)))
                  fish(2:d-1,2:d-1) = 1
          elsewhere
                  fish(2:d-1,2:d-1) = 0
          end where
          call updateboundaries(fish, d)
*          print *,fish(1:d,1:d)
          call CM_timer_stop(1)
*
* Display output every disp_interval steps
          if (mod(i-1,disp_interval) .eq. 0) then
             call show_fish(fish, d, stride, limit, off, m)
          endif
      enddo
*
* Read the timer
      call CM_timer_stop(0)
      etime0 = CM_timer_read_elapsed(0)
      etime1 = CM_timer_read_elapsed(1)
      print *, ' Computation time =', etime1, ' sec'
      print *, ' Total Elapsed time =', etime0, ' sec'

      stop
      end

      subroutine disp_init(m, d, stride, limit, off)
      integer m, d, stride, limit, off
*
* Set parameters for display
      stride = m/d
      limit = d*stride
      off = 0
      off = stride/2

*      print *,'m = ',m
*      print *,'d = ',d
*      print *,'limit = ',limit
*      print *,'stride = ',stride
*      print *,'off = ',off
      end

      subroutine show_fish(fish, d, stride, limit, off, m)
      integer fish(:,:), d, stride, limit, off
      real show(m,m)

      show = 0
*
*     Map fish grid onto show
      show(off+1:off+limit:stride,off+1:off+limit:stride) = fish(1:d,1:d)
*
* The following enlarges each pixel for better visibility
      show(:,1:m-1) = show(:,1:m-1) + show(:,2:m)
      show(1:m-1,:) = show(1:m-1,:) + show(2:m,:)
      call display (show,m,m,MINVAL(show),MAXVAL(show)+0.0001)
      end

      subroutine updateboundaries(fish, d)
      integer fish(:,:), d

      fish(1,:) = fish(d-1,:)
      fish(d,:) = fish(2,:)
      fish(:,1) = fish(:,d-1)
      fish(:,d) = fish(:,2)
      end