PROGRAM GameOfLife
! 2D cellular automaton ("Game of Life")
! (c) Pekka Manninen, Jussi Enkovaara, CSC - IT Center for Science Ltd
! based on the code presented in Haataja et al., Rinnakkaislaskentaa MPI:lla 
! (CSC 1997)
  IMPLICIT NONE
  INTEGER, PARAMETER :: &
      max_iter = 1000, & ! maximum number of iterations
      max_width = 1000, & ! max board width
      max_height = 1000, & ! max board height
      freq = 10           ! after how many iterations we print the board
  INTEGER, DIMENSION(:,:), POINTER :: &
      board_1, board_2     ! boards
  INTEGER :: &
      height, width, &     ! board dimensions
      iter, i, &
      t0, t1, cnt, cnt_max

  ! read in the initial configuration
  CALL read_data(iter, height, width)
  ! initialize boards
  CALL init_board(board_1, board_2, height, width)
  ! develop the automaton
  CALL SYSTEM_CLOCK(t0,cnt,cnt_max)
  DO i = 1, iter
    CALL iterate(board_1, board_2)
    IF (MOD(i,freq)==0) CALL draw(board_1, i)
  END DO
  CALL SYSTEM_CLOCK(t1,cnt,cnt_max)
  WRITE(*,'(A,I4,A,F6.2,A)') 'Completed ', iter, ' iterations in ', &
       REAL(t1-t0)/REAL(cnt), ' seconds.'

CONTAINS

  SUBROUTINE read_data(iter, height, width)
    IMPLICIT NONE
    INTEGER, INTENT(OUT) :: iter, height, width
    INTEGER :: task_cols, col_mod, i, rc

    WRITE (*,*) 'Insert the number of iterations'
    READ(*,*) iter
    IF (iter < 1 .OR. iter > max_iter) iter = max_iter
    WRITE (*,*) 'Insert the size of the board (height width)'
    READ(*,*) height, width
    IF (height < 1 .OR. height > max_height) height = max_height
    IF (width < 1 .OR. width > max_width) width = max_width
  END SUBROUTINE read_data
  
  SUBROUTINE init_board(board_1, board_2, height, width)
    IMPLICIT NONE
    INTEGER, DIMENSION(:,:), POINTER :: board_1, board_2
    INTEGER, INTENT(IN) :: height, width

    ALLOCATE(board_1(height,0:width+1))
    ALLOCATE(board_2(height,0:width+1))
    ! start from a plus-sign shaped pattern
    board_1(:,:) = 0
    board_2(:,:) = 0
    board_1(height/2,:) = 1
    board_1(:,width/2) = 1
    CALL draw(board_1, 0)
  END SUBROUTINE init_board
  
  SUBROUTINE iterate(board_1, board_2)
    IMPLICIT NONE
    INTEGER, PARAMETER :: ALIVE = 1, DEAD = 0
    INTEGER, DIMENSION(:,:), POINTER :: board_1, board_2
    INTEGER, DIMENSION(:,:), POINTER :: temp_board
    INTEGER :: hl,hu,wl, wu, rc,ip,im
    INTEGER :: i,j,sum

    hl = LBOUND(board_1,1)
    hu = UBOUND(board_1,1)
    wl = LBOUND(board_1,2)
    wu = UBOUND(board_1,2)
    
    !Copy boarder areas into extra space reserved for "ghost" layer
    !This will make the computation peridic (loops below from wl+1...wu-1)
    board_1(:,wl)=board_1(:,wu-1)
    board_1(:,wu)=board_1(:,wl+1)
    
    DO i=hl,hu
      !take periodic boundary conditions into account in another way in first dimension
      ip=i+1
      im=i-1
      IF(i==hl)  im=hu
      IF(i==hu)  ip=hl
      
      DO j=wl+1,wu-1
      ! determine the number of alive neighbor cells
        sum=board_1(im,j-1)+board_1(i,j-1)+board_1(ip,j-1)+ &
            board_1(im,j)+board_1(ip,j)+ &
            board_1(im,j+1)+board_1(i,j+1)+board_1(ip,j+1);
        ! set the status of a cell in the next iteration
        IF(sum< 2 .OR. sum> 3) THEN
          board_2(i,j)=DEAD
        ELSE IF (sum .EQ. 3) THEN
          board_2(i,j)=ALIVE
        ELSE
          board_2(i,j)=board_1(i,j)
        END IF
      END DO
    END DO
  

    temp_board => board_1
    board_1 => board_2
    board_2 => temp_board
  END SUBROUTINE iterate

  SUBROUTINE draw(board_1, iter)
    IMPLICIT NONE
    INTEGER :: print_unit = 10
    INTEGER, DIMENSION(:,:), POINTER :: board_1
    INTEGER, INTENT(IN) :: iter
    CHARACTER(LEN=15) :: filename, cnt
    INTEGER :: i, w, h, ios
    h = SIZE(board_1,1)
    w = SIZE(board_1,2)
    WRITE(cnt,'(I8)') iter
    WRITE(filename,'(A)')  'life_' // TRIM(ADJUSTL(cnt)) // '.pbm'
    OPEN(print_unit, file=filename, status='replace', &
        action='write', iostat=ios)
    WRITE(print_unit,'(A)') 'P1'
    WRITE(print_unit,'(A,I4)') '# Iteration number ', iter
    WRITE(print_unit,'(2I4)') w-2, h
    DO i = 1, h
      WRITE(print_unit,'(1000I2)') board_1(i,1:w-2)
    END DO
    CLOSE(print_unit, status='keep', iostat=ios)
  END SUBROUTINE draw
  
END PROGRAM GameOfLife