      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
C INCLUDE THE MPI HEADER FILE
      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 IN COLUMNS 1 THROUGH
C LOCALMAXY AND IN ROWS 1 THROUGH LOCALMAXX OF THE TWO ARRAYS BELOW.
C THE HALO DATA IS STORED IN ROWS 0 AND NS_Size+1 AND IN COLUMNS 0 AND
C WE_Size+1 OF THOSE ARRAYS.
C 
      REAL Rabbit(0:LocalMaxX+1,0:LocalMaxY+1),
     &     Fox(0:LocalMaxX+1,0:LocalMaxY+1)
C
C THE PARAMETERS OF THE PREDATOR-PREY MODEL ARE STORED IN THE ARRAY
C NAMED model.
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, PICTORIALLY:
C
C                      lbnds(WE, 1)   lbnds(WE, 2)
C                           |              |
C                           V              V
C                           ----------------
C            lbnds(NS, 1)-> |              |
C                           |  LOCAL DATA  |
C                           |  OF CURRENT  |
C                           |   PROCESS    |
C            lbnds(NS, 2)-> |              |
C                           ----------------
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.inc
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
      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
C INCLUDE THE MPI HEADER FILE
      INCLUDE 'param.inc'
      INTEGER comm, ierror
      INTEGER npes            ! THE TOTAL NUMBER OF PROCESSES
      INTEGER dims(NPDIM)     ! THE NUMBER OF PROCESSES IN EACH DIMENSION
C
C FIND OUT HOW MANY PROCESSES THERE ARE OUT THERE.
C
      npes = 1
C
C SET THE DIMENSIONS OF THE PROCESS MESH. IN THE FIRST INSTANCE, MIMIC
C SAY A 1 X npes RING BY SETTING dims(NS)=1 and dims(WE)=npes.
C

      dims(NS) = 1
      dims(WE) = 1

C
C THE PROCESS MESH IS TO BE PERIODIC IN ALL DIMENSIONS.
C




C
C DEFINE THE MESH COMMUNICATOR.
C



      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C INITIALISE THE POPULATION OF RABBITS AND FOXES.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SetLand(Rabbit,Fox,lbnds,model,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 comm, lbnds(NADIM,2), ierror

      INTEGER gi, gj, li, lj
      INTEGER dx, dy        ! THE MAXIMUM NUMBERS OF ROWS AND
                            ! COLUMNS ASSIGNED TO A PROCESSOR
      INTEGER fx, lx        ! GLOBAL INDICES OF THE FIRST AND
                            ! LAST ROW ASSIGNED TO THIS PROCESS
      INTEGER fy, ly        ! SIMILAR AS ABOVE
      INTEGER lc, lr        ! LOCAL INDICES OF THE LAST COLUMN
                            ! AND ROW ASSIGNED TO THIS PROCESS
C
C GET PROCESS INFO.
C

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
      fx = 1
      lx = NS_Size
C
C     CHECK THE VALIDITY OF GLOBAL INDICES
C
      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
      fy = 1
      ly = WE_Size
C
C     CHECK THE VALIDITY OF GLOBAL INDICES
C
      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. CHECK THAT LOCAL ARRAY BOUNDS FOR 
C Rabbit and Fox ARE NOT EXCEEDED.
C 
      lr = NS_Size
      lc = WE_Size
      IF(lc.GT.LocalMaxY .OR. lr.GT.LocalMaxX) THEN
         PRINT *, 'LOCAL ARRAY BOUNDS EXCEEDED IN SUBROUTINE SetLand'
      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
C INCLUDE THE MPI HEADER FILE
      INCLUDE 'param.inc'
      INTEGER types(NADIM), neigh(0:NBNGH), lbnds(NADIM,2)
      INTEGER comm, ierror
      INTEGER lc, lr
C
C FIND THE RANK OF THE CURRENT PROCESS
C
      neigh(HERE) = 0
C
C FIND THE RANKS OF THE NEAREST NEIGHBOURS.
C
      neigh(NORTH) = 0
      neigh(SOUTH) = 0
      neigh(EAST) = 0
      neigh(WEST) = 0
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

      types(COLUMN) = 0

C
C DEFINE A MPI DATA TYPE FOR MOVING ONE LOCAL ROW AMONG NEAREST
C NEIGHBOUR PROCESSES. 
C

      types(ROW) = 0


      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

      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 one of the local data array.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE FillBorder(Animal,types,neigh,lbnds,comm,ierror)
      IMPLICIT NONE
C INCLUDE THE MPI HEADER FILE
      INCLUDE 'param.inc'
      REAL Animal(0:LocalMaxX+1,0:LocalMaxY+1)
      INTEGER types(NADIM),neigh(0:NBNGH),lbnds(NADIM,2)
      INTEGER comm, ierror
      INTEGER lc, lr
      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 
C MOVE LOCAL BORDER TO THE EASTERN NEIGHBOUR.




C MOVE LOCAL BORDER TO THE WESTERN NEIGHBOUR.




      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 
C MOVE LOCAL BORDER TO NEAREST NEIGHBOUR TO THE NORTH.




C MOVE LOCAL BORDER TO NEAREST NEIGHBOUR TO THE SOUTH.




      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
C INCLUDE THE MPI HEADER FILE
      INCLUDE 'param.inc'
      REAL Animal(0:LocalMaxX+1,0:LocalMaxY+1)
      INTEGER lbnds(NADIM,2), count, root, comm, ierror
      INTEGER i, j, lr, lc
      REAL p

      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


      count = NINT(p)

      RETURN
      END

