program gameoflife
implicit none
include 'mpif.h'

! mpi declarations
integer :: ierror,comm,csize,myrank,status(mpi_status_size),request 

! parameters
integer, parameter :: GSIZEX=48
integer, parameter :: GSIZEY=48
integer, parameter :: MAXGREY=1

! dlength is an array with the number of nodes in each directions
! mpigrid is the cartesian communicator
! nrow and ncol are the local dimensions of the array locarray
integer :: dlength(2),nrow,ncol
integer, allocatable :: locarray(:,:),garray(:,:)
integer mpi_grid,reorder,coords(2),snd_coords(2),ii,jj,ix,iy
integer me,proc_up,proc_down,proc_left,proc_right
logical periods(0:1)
integer, allocatable :: edge_up(:),edge_down(:),edge_left(:),edge_right(:)
integer, allocatable :: scratch(:,:) 
! derived datatypes
integer mpi_block,coltype,rowtype
integer nn,dum,istep
real*8 :: x

dlength=0
coords=0
snd_coords=0
proc_up=0
proc_down=0
proc_right=0
proc_left=0

! initialize MPI environment !
call mpi_init(ierror)
call mpi_comm_size(MPI_COMM_WORLD,csize,ierror)
call mpi_comm_rank(MPI_COMM_WORLD,myrank,ierror)

! find the optimal x and y (dlength(1,2))
call mpi_dims_create(csize,2,dlength,ierror)

! create the cartesian topology
reorder=1
periods(0)=.true.
periods(1)=.true.
call mpi_cart_create(MPI_COMM_WORLD,2,dlength,periods,reorder,mpi_grid,ierror)
call mpi_cart_coords(mpi_grid,myrank,2,coords,ierror)
write(*,*)"DEBUG: cart_coords1",myrank,coords
nrow=GSIZEX/dlength(1)
ncol=GSIZEY/dlength(2)
allocate(edge_up(ncol))
allocate(edge_down(ncol))
allocate(edge_left(nrow))
allocate(edge_right(nrow))

! initialize the global and local arrays
if(myrank.eq.0)then
        allocate(garray(GSIZEX,GSIZEY))
        garray=0 
endif 

allocate(locarray(nrow,ncol))

locarray=0
if(mod(coords(1)+coords(2)+1,2)==1)then
        locarray=0
else
        locarray=MAXGREY
endif

! uncomment if you want random numbers
!do jj=1,ncol
!        do ii=1,nrow
!                call random_number(x)
!                if(x.lt.0.5d0)locarray(ii,jj)=1
!        enddo
!enddo

call mpi_barrier(MPI_COMM_WORLD,ierror)
! create a datatype
call mpi_type_vector(ncol,nrow,GSIZEY,MPI_INTEGER,mpi_block,ierror)
call mpi_type_commit(mpi_block,ierror)

call mpi_isend(locarray(1,1),nrow*ncol,MPI_INTEGER,0,10,mpi_grid,request,ierror)

if(myrank.eq.0)then
        do ii=0,csize-1
                call mpi_cart_coords(mpi_grid,ii,2,snd_coords,ierror)
                ix=nrow*snd_coords(1)+1 
                iy=ncol*snd_coords(2)+1
                call mpi_recv(garray(ix,iy),1,mpi_block,ii,10,mpi_grid,status,ierror)
        enddo
endif

if(myrank.eq.0)then
        do iy=1,GSIZEY
        write(*,fmt='(48i1)')(garray(ix,iy),ix=1,GSIZEX)
        enddo
endif
call mpi_barrier(mpi_grid,ierror)

call mpi_type_vector(ncol,1,nrow,MPI_INTEGER,rowtype,ierror)
call mpi_type_commit(rowtype,ierror)
call mpi_cart_shift(mpi_grid,1,1,proc_left,proc_right,ierror)
call mpi_cart_shift(mpi_grid,0,1,proc_up,proc_down,ierror)


call mpi_sendrecv(locarray(1,1),nrow,MPI_INTEGER,proc_left,0, &
                  edge_right,nrow,MPI_INTEGER,proc_right,0,mpi_grid,status,ierror)
call mpi_sendrecv(locarray(1,ncol),nrow,MPI_INTEGER,proc_right,0, &
                  edge_left,nrow,MPI_INTEGER,proc_left,0,mpi_grid,status,ierror)
call mpi_sendrecv(locarray(1,1),1,rowtype,proc_up,0, &
                  edge_down,ncol,MPI_INTEGER,proc_down,0,mpi_grid,status,ierror)
call mpi_sendrecv(locarray(nrow,1),1,rowtype,proc_down,0, & 
                  edge_up,ncol,MPI_INTEGER,proc_up,0,mpi_grid,status,ierror)



!! count neighbours for internal regions

allocate(scratch(0:nrow+1,0:ncol+1))
scratch=0

do istep=1,100
        scratch(1:nrow,0)=edge_left(:)
        scratch(1:nrow,ncol+1)=edge_right(:)
        scratch(0,1:ncol)=edge_up(:)
        scratch(nrow+1,1:ncol)=edge_down(:)
        scratch(1:nrow,1:ncol)=locarray(:,:)

        do ii=1,nrow
                do jj=1,ncol
                        nn=0
                        if(scratch(ii-1,jj).eq.1)nn=nn+1
                        if(scratch(ii+1,jj).eq.1)nn=nn+1
                        if(scratch(ii,jj-1).eq.1)nn=nn+1
                        if(scratch(ii,jj+1).eq.1)nn=nn+1
                        if(nn.lt.2)locarray(ii,jj)=0
                        if(nn.eq.3)locarray(ii,jj)=1
                        if(nn.gt.3)locarray(ii,jj)=0
                 enddo
        enddo

        call mpi_barrier(mpi_grid,ierror)
        call mpi_cart_shift(mpi_grid,1,1,proc_left,proc_right,ierror)
        call mpi_cart_shift(mpi_grid,0,1,proc_up,proc_down,ierror)


        call mpi_sendrecv(locarray(1,1),nrow,MPI_INTEGER,proc_left,0, &
                edge_right,nrow,MPI_INTEGER,proc_right,0,mpi_grid,status,ierror)
        call mpi_sendrecv(locarray(1,ncol),nrow,MPI_INTEGER,proc_right,0, &
                edge_left,nrow,MPI_INTEGER,proc_left,0,mpi_grid,status,ierror)
        call mpi_sendrecv(locarray(1,1),1,rowtype,proc_up,0, &
                edge_down,ncol,MPI_INTEGER,proc_down,0,mpi_grid,status,ierror)
        call mpi_sendrecv(locarray(nrow,1),1,rowtype,proc_down,0, & 
                edge_up,ncol,MPI_INTEGER,proc_up,0,mpi_grid,status,ierror)

        if(myrank.eq.0)then
                do ii=1,csize-1
                call mpi_cart_coords(mpi_grid,ii,2,snd_coords,ierror)
                ix=nrow*snd_coords(1)+1
                iy=ncol*snd_coords(2)+1
                call mpi_recv(garray(ix,iy),1,mpi_block,ii,10,mpi_grid,status,ierror)
                enddo
        else
                call mpi_send(locarray(1,1),nrow*ncol,MPI_INTEGER,0,10,mpi_grid,ierror)
        endif

        call mpi_barrier(mpi_grid,ierror)

        ! I/O
        if(myrank.eq.0)then
                do iy=1,GSIZEY
                write(50+istep,fmt='(48i1)')(garray(ix,iy),ix=1,GSIZEX)
                enddo
        endif
enddo

call mpi_type_free(mpi_block,ierror)
call mpi_type_free(rowtype,ierror)

if(myrank.eq.0)deallocate(garray)
deallocate(locarray)
call mpi_finalize(ierror)
end program