! Implement the game of Life, as invented by J.H. Conway.

PROGRAM Life
    IMPLICIT NONE
    INTEGER, PARAMETER :: unit = 42
    INTEGER, PARAMETER ::  byte = SELECTED_INT_KIND(2)
    INTEGER :: steps, size, count, rate, n
    REAL(KIND=KIND(0.0D0)) :: stamp
    CHARACTER(LEN=100) :: filename
    CHARACTER(LEN=20) :: format
    INTEGER(KIND=byte), ALLOCATABLE :: total(:,:)
    CHARACTER, PARAMETER :: chars(0:2) = (/ '.', '*', '?' /)
    CHARACTER, ALLOCATABLE :: line(:)

! Read in the parameters.

    CALL SYSTEM_CLOCK(count,rate)
    stamp = count/(rate+0.0D0)
    READ '(A)', filename
    READ *, steps
    IF (steps < 1 .OR. steps > 1000000) THEN
        PRINT *, 'Invalid number of steps'
        STOP 1
    END IF

! Read in the size and allocate the grid, including a boundary in
! parallel mode, to simplify later code.

    OPEN (unit,file=filename,STATUS='OLD',ACTION='READ')
    READ (unit,*) size
    IF (size < 1 .OR. size > 1000) THEN
        PRINT *, 'Invalid grid size'
        STOP 1
    END IF
    ALLOCATE (total(0:size+1,0:size+1))
    total(:0,:) = 0
    total(:,:0) = 0
    total(size+1:,:) = 0
    total(:,size+1:) = 0

! Read in the data, and display the initial grid.  This is not very
! efficient in memory use, but so what?  For real programs, this
! often needs to be handled more subtly, which can be tricky.  Also,
! a loop may be more efficient than WHERE, but the latter is clearer.

    ALLOCATE (line(size))
    WRITE (format,'("(",I0,"A1)")') size
    DO n = 1,size
        READ (unit,format) line
        WHERE (line == '.')
            total(1:size,n) = 0
        ELSEWHERE (line == '*')
            total(1:size,n) = 1
        ELSEWHERE
            total(1:size,n) = 2
        END WHERE
        PRINT format, chars(total(1:size,n))
    END DO
    PRINT *
    DEALLOCATE(line)

! Check the input grid for validity.

    IF (MAXVAL(total) >= 2) THEN
        PRINT *, 'Invalid input grid'
        STOP 1
    END IF

! Now run the game.  We need to pass down most of the variables, so we
! rely on host association.  We need to keep the boundary empty for
! repeatable results.

    DO n = 1,steps
        CALL Iterate(total)
        total(:0,:) = 0
        total(:,:0) = 0
        total(size+1:,:) = 0
        total(:,size+1:) = 0
    END DO

! Print the result and time taken.

    DO n = 1,size
        PRINT format, chars(total(1:size,n))
    END DO
    PRINT *
    CALL SYSTEM_CLOCK(count,rate)
    PRINT '(A,F0.1," seconds")', 'Time taken was about ',    &
        count/(rate+0.0D0)-stamp
    DEALLOCATE(total)

CONTAINS

! This performs one step in the game on the local grid.  It assumes
! that the grid is valid, both in shape and content.  It would be
! possible to do this more subtly using WHERE, but there is little
! point.  One of the higher bits is used to store what the value will
! be next step (note that the neighbour count can reach 8).  Using a
! second grid would be needed when using a bitmap and provide the
! compiler more scope for autoparallelisation, but that's an issue for
! an OpenMP course, not a MPI one!

    SUBROUTINE Iterate (grid)
        INTEGER(KIND=byte) :: grid(:,:)
        INTEGER :: i, j, k, n

        n = UBOUND(grid,1)
        DO j = 2,n-1
            DO i = 2,n-1
                k = grid(i-1,j-1)+grid(i-1,j+1)+grid(i+1,j-1)+    &
                        grid(i+1,j+1)+grid(i-1,j)+grid(i+1,j)+    &
                        grid(i,j-1)+grid(i,j+1)
                SELECT CASE (IAND(k,15))
CASE (2);             grid(i,j) = 17*grid(i,j)
CASE (3);             grid(i,j) = grid(i,j)+16
                END SELECT
            END DO
        END DO
        grid = ISHFT(grid,-4)
    END SUBROUTINE Iterate

END PROGRAM Life
