! 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)) :: wall, cpu, x
    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.

    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

! 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.

    CALL SYSTEM_CLOCK(count,rate)
    wall = count/(rate+0.0D0)
    CALL CPU_TIME(cpu)
! The PRIVATE(n) is needed.
!$OMP PARALLEL DEFAULT(NONE), SHARED(TOTAL,SIZE,STEPS), PRIVATE(n)
    DO n = 1,steps
        CALL Iterate(total)
!$OMP SECTIONS
!$OMP SECTION
        CALL Zero(total(:0,:size))
!$OMP SECTION
        CALL Zero(total(1:,:0))
!$OMP SECTION
        CALL Zero(total(size+1:,1:))
!$OMP SECTION
        CALL Zero(total(:size,size+1:))
!$OMP END SECTIONS
    END DO
!$OMP END PARALLEL
    CALL SYSTEM_CLOCK(count,rate)
    wall = count/(rate+0.0D0)-wall
    CALL CPU_TIME(x)
    cpu = x-cpu

! Print the result and time taken.

    DO n = 1,size
        PRINT format, chars(total(1:size,n))
    END DO
    PRINT *
    PRINT '("Time = ",F0.2," seconds, CPU = ",F0.2," seconds")', &
        wall,cpu
    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
! a more advanced OpenMP course!

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

! n is PRIVATE by default.
        n = UBOUND(grid,1)
!$OMP SINGLE
        DO j = 2,n-1
            DO i = 2,n-1
                m = Get(grid(i,j))
                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);           m = 17*m
CASE (3);           m = m+16
                END SELECT
                CALL Set(grid(i,j),m)
            END DO
        END DO
!$OMP END SINGLE
!$OMP DO PRIVATE(i,j)
        DO j = 1,n
            DO i = 1,n
                CALL Set(grid(i,j),ISHFT(Get(grid(i,j))+0,-4))
            END DO
        END DO
!$OMP END DO
    END SUBROUTINE Iterate

! These procedures are a gratuitously inefficient way of performing
! these tasks, but are used to ensure that there is a chance of wrong
! answers if the program has race conditions.  The obvious and more
! efficient code will work even in that case, primarily because of
! atomicity.  Despite my best attempts, some compilers will optimise
! these to uselessness anyway.  VOLATILE is not used because it's
! not designed or suitable for this purpose.

    SUBROUTINE Zero (grid)
        INTEGER(KIND=byte) :: grid(:,:)
        INTEGER :: i, j

        DO j = 1,UBOUND(grid,2)
            DO i = 1,UBOUND(grid,1)
                grid(i,j) = IEOR(grid(i,j)+0,INT(Z'3F'))
            END DO
        END DO
        DO j = 1,UBOUND(grid,2)
            DO i = 1,UBOUND(grid,1)
                CALL Set(grid(i,j),0)
            END DO
        END DO
    END SUBROUTINE Zero

    FUNCTION Get (value)
        INTEGER(KIND=byte) :: Get, value

        Get = value
        value = IEOR(value+0,INT(Z'3F'))
    END FUNCTION Get

    SUBROUTINE Set (location, value)
        INTEGER(KIND=byte) :: location
        INTEGER, INTENT(IN) :: value

        location = ISHFT(INT(location),-6)+value
    END SUBROUTINE Set

END PROGRAM Life
