PROGRAM life
      IMPLICIT NONE
!
!     MAXLOOP iterations of an NxN life board
!
      INTEGER, PARAMETER :: N=200, MAXLOOP=30
      INTEGER loop, i, j, cr, nx, r, c, nb
      INTEGER, DIMENSION(N,N,2) :: step
      CHARACTER*(16) name
!      
!     OMP variables
      INTEGER :: OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS
      REAL(8) :: OMP_GET_WTIME, wall1, wall2, wall
!   
!     Initialise step
!
      step = 0
      ! Ten active cells in a line
      step(N/2,((N/2)-5):((N/2)+4),1) = 1
!
! Insert parallel region, pay attention to  shared or private variables
!
! Print starting configuration
!
! Only master thread executes the three following statements
	  loop = 0
	  WRITE(name,"(A4,I4.4)") "Life",loop
      CALL INTDATA2MPGM(N,N,step(:,:,1),name,4)

!
!     Perform MAXLOOP updates
! 
	  wall1 = 0.0; wall2 = 0.0; wall = 0.0
      DO loop = 1, MAXLOOP
!$         wall1 = OMP_GET_WTIME()
         cr = 2-mod(loop,2)
         nx = 2-mod(loop+1,2)
         ! Do it in parallel
         DO r = 2, N-1
            DO c = 2, N-1
               step(r,c,nx) = step(r,c,cr)   ! Inizialize cell
               !  Per each cell
               nb = 0
               DO i = r-1,r+1
                  DO j = c-1,c+1
                      ! Count active neighbours
                  	  IF ( .NOT. (i == r .AND. j == c) ) THEN
                  	  	  nb = nb + step(i,j,cr)
                  	  END IF
                  END DO
               END DO
               IF (nb.lt.2 .or. nb.gt.3) THEN
                  step(r,c,nx) = 0
               ELSE IF (nb.eq.3) THEN
                  step(r,c,nx) = 1
               END IF
            END DO
         END DO
         
!        Write out computed state of board
!
!$	     wall2 = OMP_GET_WTIME()
	     wall = wall + wall2 - wall1
! Only master thread executes the two following statements
	       WRITE(name,"(A4,I4.4)") "Life",loop
    	   CALL INTDATA2MPGM(N,N,step(:,:,nx),name,4)

      END DO
! Only master thread executes the following statement
	  WRITE(*,*) "Computing elapsed time: ",wall

! Close parallel region
	
      END
!
subroutine IntData2Mpgm(s1,s2,idata,name,Mag)
!  Simple subroutine to dump integer data in a PGM format
   implicit none
   INTEGER, INTENT(IN) :: s1, s2, Mag
   INTEGER, DIMENSION(s1,s2), intent(in) :: idata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni, r, m
   integer :: vp, vs
   real :: rmin, rmax
   character(80) :: fname
!
!  Write on unit 700 with  PGM format
   ouni = 700
   fname = TRIM(name)//".pgm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      print*,"!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P2"
!  Dimensions
   write(ouni,*) s1*Mag, s2*Mag
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(idata); rmax = MAXVAL(idata)
   print*,"rmin, rmax = ", rmin, rmax
   vs = 0
   do i = 1, s1
      DO r = 1, Mag
      do j = 1, s2
         vp = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         DO m = 1, Mag
            write(ouni,"(i4)",advance="no") vp
            vs = vs + 1
         END DO
         if (vs >= 10 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
      END DO
   end do
   close(ouni)
   return
end subroutine IntData2Mpgm