!**********************************************************************
!   Virtual topologies and communicators
!     
!****************************************************************************
program communicators

  include 'mpif.h'
  integer, parameter :: DP=kind(0.0D0)
  
  integer :: myid,nproc,rank,i,j,subrank,data,mydata
  integer :: stat,ierr,proc_grid, proc_row,reorder,ndim
  integer,dimension(2)::coords,pos,dims,periods
  
  

  ! Initialize MPI 
  call MPI_INIT( ierr )
  call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, nproc, ierr )
  reorder=1;ndim=2;dims=0;periods=0




  ! Create a virtual 2D-grid topology
  call MPI_Dims_create(nproc,ndim,dims,ierr)
  call MPI_Cart_create(MPI_COMM_WORLD,ndim,dims,periods,reorder,proc_grid,ierr)
  call MPI_Comm_rank(proc_grid,myid,ierr)    ! Note: use proc_grid




  ! Create a communicator for each row
  call MPI_Cart_coords(proc_grid,myid,ndim,coords,ierr)
  call MPI_Comm_split(proc_grid,coords(1),coords(2),proc_row,ierr)
  call MPI_Comm_rank(proc_row,subrank,ierr)




  ! Broadcast within a row
  if (subrank==0) then
     mydata=coords(1)
  end if
  call MPI_Bcast(mydata,1,MPI_INTEGER,0,proc_row,ierr);



  ! Check the result of Broadcast
  do i=0,dims(1)-1
     do j=0,dims(2)-1
      pos(1)=i; pos(2)=j;
      call MPI_Cart_rank(proc_grid,pos,rank,ierr);
	if (rank==myid) then
           write(*,*) 'Rank: ',myid,'  Coords: ',i,' ',j,'   Data: ',mydata
        end if
     end do
  end do



  ! Exit and clean up MPI variables
  call MPI_Comm_free(proc_row,ierr);
  call MPI_Comm_free(proc_grid,ierr);

  call MPI_FINALIZE(ierr)
  
end program communicators