*  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 fish4

      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 = 66

      integer 	fish(MAXG,MAXG), direction(MAXG,MAXG), move(MAXG,MAXG)
CMF$  layout    fish(:news,:news),direction(:news,:news), move(:news,:news)

      integer 	d, i, j, 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
          call CM_timer_start(1)
*
* Pick a random direction ( 1=E, 2=N, 3=W, 4=S, 0= no move)
          call cmf_random(direction,5)

*          print *,nbrcnt(2:d-1,2:d-1)
          call updateboundaries(fish, direction, d)
          call MoveFishEast(fish,direction,move,d)
          call MoveFishNorth(fish,direction,move,d)
          call MoveFishWest(fish,direction,move,d)
          call MoveFishSouth(fish,direction,move,d)
*          print *,fish(1:d,1:d)
          call CM_timer_stop(1)
*
* Display output every dis_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 MoveFishEast(fish, direction, move, d)
      integer fish(:,:), direction(:,:), move(:,:), d

      where ((fish(2:d+1,1:d+1) .eq. 1) .and. (direction(2:d+1,1:d+1) .eq. 1) .and. (fish(2:d+1,2:d+2) .eq. 0))
          move(2:d+1,1:d+1) = 1
      elsewhere
          move(2:d+1,1:d+1) = 0
      end where
      where ( ((fish(2:d+1,2:d+1) .eq. 1) .and. (move(2:d+1,2:d+1) .eq. 0)) .or. (move(2:d+1,1:d) .eq. 1))
          fish(2:d+1,2:d+1) = 1
      elsewhere
          fish(2:d+1,2:d+1) = 0
      end where
      call updateboundaries(fish, direction, d)
      end 

      subroutine MoveFishNorth(fish, direction, move, d)
      integer fish(:,:), direction(:,:), move(:,:), d

      where ((fish(2:d+2,2:d+1) .eq. 1) .and. (direction(2:d+2,2:d+1) .eq. 2) .and. (fish(1:d+1,2:d+1) .eq. 0))
          move(2:d+2,2:d+1) = 1
      elsewhere
          move(2:d+2,2:d+1) = 0
      end where
      where ( ((fish(2:d+1,2:d+1) .eq. 1) .and. (move(2:d+1,2:d+1) .eq. 0)) .or. (move(3:d+2,2:d+1) .eq. 1))
          fish(2:d+1,2:d+1) = 1
      elsewhere
          fish(2:d+1,2:d+1) = 0
      end where
      call updateboundaries(fish, direction, d)
      end 

      subroutine MoveFishWest(fish, direction, move, d)
      integer fish(:,:), direction(:,:), move(:,:), d

      where ((fish(2:d+1,2:d+2) .eq. 1) .and. (direction(2:d+1,2:d+2) .eq. 3) .and. (fish(2:d+1,1:d+1) .eq. 0))
          move(2:d+1,2:d+2) = 1
      elsewhere
          move(2:d+1,2:d+2) = 0
      end where
      where ( ((fish(2:d+1,2:d+1) .eq. 1) .and. (move(2:d+1,2:d+1) .eq. 0)) .or. (move(2:d+1,3:d+2) .eq. 1))
          fish(2:d+1,2:d+1) = 1
      elsewhere
          fish(2:d+1,2:d+1) = 0
      end where
      call updateboundaries(fish, direction, d)
      end 

      subroutine MoveFishSouth(fish, direction, move, d)
      integer fish(:,:), direction(:,:), move(:,:), d

      where ((fish(1:d+1,2:d+1) .eq. 1) .and. (direction(1:d+1,2:d+1) .eq. 4) .and. (fish(2:d+2,2:d+1) .eq. 0))
          move(1:d+1,2:d+1) = 1
      elsewhere
          move(1:d+1,2:d+1) = 0
      end where
      where ( ((fish(2:d+1,2:d+1) .eq. 1) .and. (move(2:d+1,2:d+1) .eq. 0)) .or. (move(1:d,2:d+1) .eq. 1))
          fish(2:d+1,2:d+1) = 1
      elsewhere
          fish(2:d+1,2:d+1) = 0
      end where
      call updateboundaries(fish, direction, d)
      end 

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

      fish(1,:) = fish(d+1,:)
      fish(d+2,:) = fish(2,:)
      fish(:,1) = fish(:,d+1)
      fish(:,d+2) = fish(:,2)
      direction(1,:) = direction(d+1,:)
      direction(d+2,:) = direction(2,:)
      direction(:,1) = direction(:,d+1)
      direction(:,d+2) = direction(:,2)
      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 = stride/2

      print *,'limit = ',limit
      print *,'stride = ',stride
      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