      PROGRAM ECO
      IMPLICIT NONE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C THIS FILE CONTAINS FORTRAN CODE FOR MODELLING A SIMPLE PREDATOR-PREY
C MODEL.  THE TWO ANIMAL POPULATIONS MODELLED ARE RABBITS AND FOXES THAT
C LIVE ON SOME PIECE OF LAND.  THE TWO ANIMAL POPULATIONS MODELLED ARE
C RABBITS AND FOXES THAT LIVE ON SOME PIECE OF LAND. THE ANIMAL
C POPULATIONS ARE REPRESENTED BY TWO TWO-DIMENSIONAL ARRAYS OF SIZE
C NS_Size BY WE_Size. PRECISELY, THE NUMBER OF FOXES THAT LIVE IN THE
C I,J-TH STRETCH OF LAND (WITH 1<=I<=NS_Size AND 1<=J<=WE_Size) IS STORED
C AS FOX (I, J). THE CORRESPONDING FIGURE FOR RABBITS IS STORED IN RABBIT
C (I, J). BOUNDARY CONDITIONS ARE ENFORCED SO THAT THE LAND IS ACTUALLY
C A TORUS AND CAN BE THOUGHT TO EXTEND PERIODICALLY BEYOND THE GIVEN
C BOUNDS NS_Size AND WE_Size.  THE POPULATIONS OF EACH STRETCH IS UPDATED
C AT EACH GENERATION ACCORDING TO SOME SIMPLE RULES AND THE TOTAL
C POPULATIONS ARE SUMMED AT REGULAR INTERVALS OF TIME.
C 
C YOUR FIRST MAJOR TASK WHEN PARALLELIZING THIS CODE WILL BE TO DEFINE A
C GEOMETRIC DATA DECOMPOSITION. COMMENTS HAVE BEEN ADDED TO ASSIST YOU.
C THEY ASSUME THAT THE PROCESSES HAVE BEEN ORGANISED INTO A TWO
C DIMENSIONAL CARTESIAN TOPOLOGY AS SHOWN IN THE DIAGRAM BELOW. ARROWS
C POINT TOWARD THE DIRECTION OF INCREASING INDICES.
C
C        +-----J----->       +-----Y----->
C        |                   |              
C        |                   |                  +--> EAST
C        I   A(I,J)          X  PROC(X,Y)       |        
C        |                   |                  V 
C        |                   |                SOUTH
C        V                   V               
C
C PROCEDURES AND THEIR ROLES:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C SetMesh: DEFINES THE PROCESS TOPOLOGY TO BE A 2-DIMENSIONAL
C          TORUS.
C SetLand: DEFINES A GEOMETRIC PARTITIONING OF THE PROBLEM.
C          AND INITIALISES THE LOCAL RABBIT AND FOX POPULATIONS.
C SetComm: DEFINES THE MPI DATA TYPES FOR EXCHANGING ROWS AND
C          COLUMNS AMONG NEAREST NEIGHBOUR PROCESSES.
C Evolve:  IS CALLED REPEATEDLY TO COMPUTE THE NEXT GENERATIONS
C          OF BOTH FOXES AND RABBITS. Evolve CALLS THE SUBROUTINE
C          FillBorder TO EXCHANGE HALO DATA BETWEEN PROCESSES.
C GetPopulation: COMPUTES THE TOTAL POPULATION OF THE SPECIES
C          IT IS GIVEN AS ARGUMENT.
C FillBorder: DOES THE ACTUAL HALO DATA SWAPS FOR THE SPECIES IT
C          IS GIVEN HAS ARGUMENT.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ARRAY DIMENSIONS AND RELATED PARAMETERS
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      INCLUDE 'mpi.inc'
      INCLUDE  'param.inc'
      INTEGER NITER, PERIOD
      PARAMETER (NITER=50, PERIOD=1)
C
C ARRAYS
C
C The arrays storing the number of animals in each stretch of land
C have borders that are used as buffers for data from nearest neighbour
C processors and to simplify the coding of periodic conditions. The
C data assigned to the current process is stored starting from column 1 
C row 1 of the two arrays below.
C 
      REAL Rabbit(0:LocalMaxX+1,0:LocalMaxY+1),
     &     Fox(0:LocalMaxX+1,0:LocalMaxY+1)
C
C The array named model is used to pass the model parameters
C across procedures.
C
      REAL model(2,3)
C
C The local indices of a stretch of land are the indices used by the
C process assigned to that stretch of land. The global indices of the
C same stretch of land are its local indices when only one process is
C involved in the computation, i.e. when the rabbit and fox populations
C are kept in one big array each. The array lbnds provides a mapping
C between local and global indices, specifically:
C 
C The north-eastern most stretch of land assigned to the current process
C has global indices (lbnds (NS, 1), lbnds (WE, 1) ) and
C the south-western most stretch of land assigned to the current process
C has global indices (lbnds (NS, 2), lbnds (WE, 2) ).
C
      INTEGER lbnds(NADIM,2)
C
C The ranks of the current process and of its nearest neighbours in
C the 2-dimensional mesh are kept in the next array. Constants
C HERE, NORTH, SOUTH, EAST and WEST are defined in the file param.h
C to help use this array.
C 
      INTEGER neigh(0:NBNGH)
C
C FUNCTION & INTRINSIC
C
      INTEGER NINT, MOD
      INTRINSIC NINT, MOD, MAX
C
C LOOP INDICES AND COUNTERS
C
      INTEGER K
      INTEGER nbrab, nbfox
C
C MPI PARAMETERS AND VARIABLES
C
      INTEGER comm, ierror, CHK
      INTEGER types(NADIM)
C
C CALL MPI
C
      CALL MPI_Init(ierror)
C
C INITIALISE THE PROBLEM
C
      CALL SetMesh(comm,ierror)
      CALL SetLand(Rabbit,Fox,lbnds,model,comm,ierror)
      CALL SetComm(types,neigh,lbnds,comm,ierror)
C
C ITERATE
C
      DO K = 1, NITER

         CALL Evolve(Rabbit,Fox,model,types,neigh,lbnds,
     &      comm,ierror)
         IF( MOD(K,PERIOD).EQ.0 ) THEN
            CALL GetPopulation(Rabbit,lbnds,nbrab,0,comm,ierror)
            CALL GetPopulation(Fox,lbnds,nbfox,0,comm,ierror)
            IF( neigh(HERE).EQ.0 ) WRITE (*, 100) 'Year ',
     &         K, ' : ', nbrab, ' rabbits and ', nbfox, ' foxes.'
         END IF
      END DO
 100  FORMAT(A5,I3,A3,I10,A13,I10,A7)
C
C CLOSE MPI
C
      CALL MPI_Finalize(ierror)

      STOP
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C SET UP A VIRTUAL TOPOLOGY FOR A 2-D MESH OF PROCESSES.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SetMesh(comm,ierror)
      IMPLICIT NONE
#include <mpif.h>
      INCLUDE 'param.inc'
      INTEGER comm, ierror
      INTEGER npes            ! THE TOTAL NUMBER OF PROCESSES
      INTEGER dims(NPDIM)     ! THE NUMBER OF PROCESSES IN EACH DIMENSION
      INTEGER i, CHK
      LOGICAL periods(NPDIM), reorder
C
C FIND OUT HOW MANY PROCESSES THERE ARE OUT THERE.
C
      CALL MPI_Comm_size(MPI_COMM_WORLD,npes,ierror)
C
C SET THE DIMENSIONS OF THE PROCESS MESH. IN THE FIRST INSTANCE, MIMIC
C A 1 X npes RING BY SETTING dims[NS]=1 and dims[WE]=npes.
C
      DO i = 1, NPDIM
         dims(i) = 0
      END DO
      CALL MPI_Dims_create(npes,NPDIM,dims,ierror)
C
C THE PROCESS MESH IS TO BE PERIODIC IN ALL DIMENSIONS.
C
      DO i = 1, NPDIM
         periods(i) = .TRUE.
      END DO
      reorder = .TRUE.
C
C DEFINE THE MESH COMMUNICATOR.
C
      CALL MPI_Cart_create(MPI_COMM_WORLD,NPDIM,dims,
     &   periods,reorder,comm,ierror)

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C INITIALISE THE POPULATION OF RABBITS AND FOXES.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SetLand(Rabbit,Fox,lbnds,model,comm,ierror)
      IMPLICIT NONE
#ifdef T3D
      INCLUDE 'mpi.inc'
#endif
      INCLUDE 'param.inc'
      REAL Rabbit(0:LocalMaxX+1,0:LocalMaxY+1),
     &   Fox(0:LocalMaxX+1,0:LocalMaxY+1), model(2,3)
      INTEGER comm, lbnds(NADIM,2), ierror, CHK

      INTEGER gi, gj, li, lj
      INTEGER dx, dy, fx, fy, lx, ly
      INTEGER lc, lr
      INTEGER dims(NPDIM)   ! SEE THE COMMENTS IN SetComm ABOVE
      INTEGER coords(NPDIM) ! THE CURRENT PROCESS HAS COORDINATES
                            ! (coords(NS),coords(WE)).
      LOGICAL periods(NPDIM)
C
C GET PROCESS INFO.
C
      CALL MPI_Cart_get(comm,NPDIM,dims,periods,coords,ierror)
C
C SET THE PARAMETERS OF THE PREDATOR-PREY MODEL
C
       model(RAB,SAME) = -0.2
       model(RAB,OTHER) = 0.6
       model(RAB,MIGRANT) = 0.01
       model(FXX,OTHER) = 0.6
       model(FXX,SAME) = -1.8
       model(FXX,MIGRANT) = 0.02
C
C USE A GEOMETRIC DECOMPOSITION AND COMPUTE THE GLOBAL INDICES OF THE
C FIRST AND LAST LOCAL ROW OF THE CURRENT PROCESS.
C 
      dx = (NS_Size+dims(NS)-1)/dims(NS)
      fx = coords(NS)*dx + 1
      lx = fx + dx -1
      IF( fx .GT. NS_Size ) THEN
         lx = fx - 1
      ELSE IF( lx .GT. NS_Size ) THEN
         lx = NS_Size
      END IF
C
C DO THE SAME FOR COLUMNS INSTEAD OF ROWS.
C
      dy = (WE_Size+dims(WE)-1)/dims(WE)
      fy = coords(WE)*dy + 1
      ly = fy + dy -1
      IF( fy .GT. WE_Size ) THEN
         ly = fy - 1
      ELSE IF( ly .GT. WE_Size ) THEN
         ly = WE_Size
      END IF
C
C COMPUTE THE INDICES OF THE LAST ROW AND COLUMN OF DATA IN THE 
C LOCAL ARRAYS Fox AND Rabbit.
C 
      lr = lx - fx + 1
      lc = ly - fy + 1
      IF(lc.GT.LocalMaxY .OR. lr.GT.LocalMaxX) THEN
         CALL MPI_Abort(MPI_COMM_WORLD,-1)
      END IF

C
C FILL-IN THE LOCAL ARRAYS FOR FOXES AND RABBITS.
C
      DO gj = fy, ly
         lj = gj - fy + 1
         DO gi = fx, lx
             li = gi - fx + 1
             Rabbit(li,lj) =
     &          128.0*(gi-1)*(NS_Size-gi)*(gj-1)*(WE_Size-gj) /
     &          REAL(NS_Size*NS_Size*WE_Size*WE_Size)
             Fox(li,lj) =
     &          8.0*(gi/REAL(NS_Size)-0.5)*(gi/REAL(NS_Size)-0.5)+ 
     &          8.0*(gj/REAL(WE_Size)-0.5)*(gj/REAL(WE_Size)-0.5)
         END DO
      END DO
C
C SAVE THE LOWER AND UPPER BOUNDS OF GLOBAL ARRAY INDICES FOR LATER USE.
C
      lbnds(NS,1) = fx
      lbnds(NS,2) = lx
      lbnds(WE,1) = fy
      lbnds(WE,2) = ly

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C INITIALISE MPI TYPES AND COMPUTE NEAREST NEIGHBOURS RANKS
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SetComm(types,neigh,lbnds,comm,ierror)
      IMPLICIT NONE
#include <mpif.h>
      INCLUDE 'param.inc'
      INTEGER types(NADIM), neigh(0:NBNGH), lbnds(NADIM,2)
      INTEGER comm, ierror
      INTEGER lc, lr, CHK
C
C FIND THE RANK OF THE CURRENT PROCESS
C
      CALL MPI_Comm_rank(comm,neigh(HERE),ierror)
C
C FIND THE RANKS OF THE NEAREST NEIGHBOURS. IF YOU HAVE SET dims(WE)=npes
C IN SetComm IT MAY BE SAFER TO SPECIFY neigh[SOUTH]=neigh[NORTH]=
C neigh[HERE] AT THIS POINT.
C
      CALL MPI_Cart_shift(comm,0,1,neigh(NORTH),neigh(SOUTH),
     &   ierror)
      CALL MPI_Cart_shift(comm,1,1,neigh(EAST),neigh(WEST),
     &   ierror)
C
C COMPUTE THE INDICES OF THE LAST ROW AND COLUMN OF DATA IN THE 
C LOCAL ARRAYS Fox AND Rabbit.
C 
      lr = lbnds(NS,2) - lbnds(NS,1) + 1
      lc = lbnds(WE,2) - lbnds(WE,1) + 1
C
C DEFINE A MPI DATA TYPE FOR MOVING ONE LOCAL COLUMN AMONG NEAREST
C NEIGHBOUR PROCESSES.
C
      CALL MPI_Type_vector(1,lr,LocalMaxX+2,MPI_REAL,
     &   types(COLUMN), ierror)
      CALL MPI_Type_commit(types(COLUMN),ierror)
C
C DEFINE A MPI DATA TYPE FOR MOVING ONE LOCAL ROW AMONG NEAREST
C NEIGHBOUR PROCESSES. THAT DATA TYPE NEED NOT BE USED IF 
C dims(WE)=npes.
C
      CALL MPI_Type_vector(lc,1,LocalMaxX+2,MPI_REAL,
     &   types(ROW), ierror)
      CALL MPI_Type_commit(types(ROW),ierror)

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C COMPUTE THE NEXT GENERATIONS
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE Evolve(Rabbit,Fox,model,types,neigh,lbnds,comm,
     & ierror)
      IMPLICIT NONE
      INCLUDE 'param.inc'
      REAL Rabbit(0:LocalMaxX+1,0:LocalMaxY+1),
     &   Fox(0:LocalMaxX+1,0:LocalMaxY+1), model(2,3)
      INTEGER types(NADIM),neigh(0:NBNGH),lbnds(NADIM,2), comm
      INTEGER ierror

      INTEGER li, lj
      INTEGER lc, lr, CHK

      REAL TRabbit(0:LocalMaxX+1,0:LocalMaxY+1),
     &   TFox(0:LocalMaxX+1,0:LocalMaxY+1)
      REAL AlR, BtR, MuR, AlF, BtF, MuF

      AlR = model(RAB,SAME)
      BtR = model(RAB,OTHER)
      MuR  = model(RAB,MIGRANT)
      BtF = model(FXX,SAME)
      AlF = model(FXX,OTHER)
      MuF  = model(FXX,MIGRANT)
C
C FILL-IN THE BORDERS OF THE LOCAL Rabbit and Fox ARRAYS.
C
      CALL FillBorder(Rabbit,types,neigh,lbnds,comm,ierror)
      CALL FillBorder(Fox,types,neigh,lbnds,comm,ierror)
C
C COMPUTE THE INDICES OF THE LAST ROW AND COLUMN OF DATA IN THE 
C LOCAL ARRAYS Fox AND Rabbit.
C 
      lr = lbnds(NS,2) - lbnds(NS,1) + 1
      lc = lbnds(WE,2) - lbnds(WE,1) + 1
C
C UPDATE THE LOCAL POPULATION DATA.
C
      DO lj = 1, lc
         DO li = 1, lr
            TRabbit(li,lj) = (1.0+AlR-4.0*MuR)*Rabbit(li,lj) +
     &         BtR*Fox(li,lj) +
     &         MuR*(Rabbit(li,lj-1)+Rabbit(li,lj+1)+
     &             Rabbit(li-1,lj)+Rabbit(li+1,lj))

            TFox(li,lj) = AlF*Rabbit(li,lj) +
     &         (1.0+BtF-4.0*MuF)*Fox(li,lj) +
     &         MuF*(Fox(li,lj-1)+Fox(li,lj+1)+
     &             Fox(li-1,lj)+Fox(li+1,lj))
         END DO
      END DO
C
C ENSURE THE NUMBERS IN Rabbit and Fox ARE NON-NEGATIVE.
C
      DO lj = 1, lc
         DO li = 1, lr
            Rabbit(li,lj) = MAX(0.0, TRabbit(li,lj))
            Fox(li,lj) = MAX(0.0, TFox(li,lj))
         END DO
      END DO

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Set the margin of a data array
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE FillBorder(Animal,types,neigh,lbnds,comm,ierror)
      IMPLICIT NONE
#include <mpif.h>
      INCLUDE 'param.inc'
      REAL Animal(0:LocalMaxX+1,0:LocalMaxY+1)
      INTEGER types(NADIM),neigh(0:NBNGH),lbnds(NADIM,2)
      INTEGER comm, ierror, status(MPI_STATUS_SIZE)
      INTEGER lc, lr, CHK
      INTEGER i, j
C
C COMPUTE THE INDICES OF THE LAST ROW AND COLUMN OF DATA IN THE 
C LOCAL ARRAYS Fox AND Rabbit.
C 
      lr = lbnds(NS,2) - lbnds(NS,1) + 1
      lc = lbnds(WE,2) - lbnds(WE,1) + 1
C
C MOVE BORDERS TO THE EAST AND TO THE WEST.
C
      IF( lc .EQ. WE_Size ) THEN
C THE CURRENT NODE HAS NO PROPER EAST AND WEST NEIGHBOURS.
         DO i = 1, lr
            Animal(i,0) = Animal(i,lc)
            Animal(i,lc+1) = Animal(i,1)
         END DO
      ELSE IF( lc .GT. 0 ) THEN
C MOVE LOCAL BORDER TO THE EASTERN NEIGHBOUR.
         CALL MPI_Send(Animal(1,1),1,types(COLUMN),neigh(EAST),
     &      EAST,comm,ierror)
         CALL MPI_Recv(Animal(1,lc+1),1,types(COLUMN),neigh(WEST),
     &      EAST,comm,status,ierror)
C MOVE LOCAL BORDER TO THE WESTERN NEIGHBOUR.
         CALL MPI_Send(Animal(1,lc),1,types(COLUMN),neigh(WEST),
     &      WEST,comm,ierror)
         CALL MPI_Recv(Animal(1,0),1,types(COLUMN),neigh(EAST),
     &      WEST,comm,status,ierror)
      ELSE
C SHIFT DATA RECEIVED FROM NEAREST NEIGHBOURS.
C EASTWARD FIRST
         CALL MPI_Recv(Animal(1,lc+1),1,types(COLUMN),neigh(WEST),
     &      EAST,comm,status,ierror)
         CALL MPI_Send(Animal(1,1),1,types(COLUMN),neigh(EAST),
     &      EAST,comm,ierror)
C WESTWARD NEXT
         CALL MPI_Recv(Animal(1,0),1,types(COLUMN),neigh(EAST),
     &      WEST,comm,status,ierror)
         CALL MPI_Send(Animal(1,lc),1,types(COLUMN),neigh(WEST),
     &      WEST,comm,ierror)
      END IF    
C
C MOVE BORDERS TO THE NORTH AND TO THE SOUTH.
C
      IF( lr .EQ. NS_Size ) THEN
C THE CURRENT PROCESS HAS NO PROPER NORTH AND SOUTH NEIGHBOURS.
         DO j = 1, lc
            Animal(0,j) = Animal(lr,j)
            Animal(lr+1,j) = Animal(1,j)
         END DO
      ELSE IF( lr .GT. 0 ) THEN
C MOVE LOCAL BORDER TO NEAREST NEIGHBOUR TO THE NORTH.
         CALL MPI_Send(Animal(1,1),1,types(ROW),neigh(NORTH),
     &      NORTH,comm,ierror)
         CALL MPI_Recv(Animal(lr+1,1),1,types(ROW),neigh(SOUTH),
     &      NORTH,comm,status,ierror)
C MOVE LOCAL BORDER TO NEAREST NEIGHBOUR TO THE SOUTH.
         CALL MPI_Send(Animal(lr,1),1,types(ROW),neigh(SOUTH),
     &      SOUTH,comm,ierror)
         CALL MPI_Recv(Animal(0,1),1,types(ROW),neigh(NORTH),
     &      SOUTH,comm,status,ierror)
      ELSE
C THE CURRENT PROCESS HAS NO ASSIGNED PIECE OF LAND
C SHIFT DATA NORTHWARD FIRST
         CALL MPI_Recv(Animal(lr+1,1),1,types(ROW),neigh(SOUTH),
     &      NORTH,comm,status,ierror)
         CALL MPI_Send(Animal(1,1),1,types(ROW),neigh(NORTH),
     &      NORTH,comm,ierror)
C SHIFT DATA SOUTHWARD NEXT
         CALL MPI_Recv(Animal(0,1),1,types(ROW),neigh(NORTH),
     &      SOUTH,comm,status,ierror)
         CALL MPI_Send(Animal(lr,1),1,types(ROW),neigh(SOUTH),
     &      SOUTH,comm,ierror)
      END IF

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Compute the number of individuals in a population.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE GetPopulation(Animal,lbnds,count,root,comm,ierror)
      IMPLICIT NONE
#include <mpif.h>
      INCLUDE 'param.inc'
      REAL Animal(0:LocalMaxX+1,0:LocalMaxY+1)
      INTEGER lbnds(NADIM,2), count, root, comm, ierror
      INTEGER i, j, lr, lc, CHK
      REAL p, q

      INTRINSIC NINT

C
C COMPUTE THE INDICES OF THE LAST ROW AND COLUMN OF DATA IN THE 
C LOCAL ARRAYS Fox AND Rabbit.
C 
      lr = lbnds(NS,2) - lbnds(NS,1) + 1
      lc = lbnds(WE,2) - lbnds(WE,1) + 1
C
C COMPUTE THE LOCAL POPULATION.
C
      p = 0.0
      DO j = 1, lc
         DO i = 1, lr
            p = p + Animal(i,j)
         END DO
      END DO
C
C COMPUTE THE GLOBAL POPULATION.
C
      CALL MPI_Reduce(p,q,1,MPI_REAL,MPI_SUM,root,comm,ierror)

      count = NINT(q)

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PRINTEM(A,lbnds,comm,ierror)
      INCLUDE 'param.inc'
      REAL A(0:LocalMaxX+1,0:LocalMaxY+1)
      INTEGER lbnds(NADIM,2), comm, ierror
      INTEGER gi, gj, li, lj
      INTEGER fx, fy, lx, ly

      fx = lbnds(NS,1)
      lx = lbnds(NS,2)
      fy = lbnds(WE,1)
      ly = lbnds(WE,2)

      do gi = 1,NS_Size
         do gj = 1, WE_Size
            CALL MPI_barrier(comm,ierror)
            if((gi.ge.fx .and. gi.le.lx) .and.
     &         (gj.ge.fy .and. gj.le.ly))
     &            print *, gi, gj, ':',A(gi-fx+1,gj-fy+1)
         end do
      end do
      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION CHK(ierror)
#include <mpif.h>
      INTEGER CHK, rlen,  jerror
      CHARACTER*(MPI_MAX_ERROR_STRING) msg
      IF(ierror.NE.0)THEN
         CALL MPI_Error_string(ierror,msg,rlen,jerror)
         PRINT *, msg
         CALL MPI_Abort(MPI_COMM_WORLD,jerror)
      END IF

      CHK = ierror
      RETURN
      END

