impure subroutine parallel_obj_Topology(this,is_periodic,Ng,Nb)
      !> Builds a Cartesian topolgy with MPI.
      implicit none
      class(parallel_obj), intent(inout) :: this                               !! Parallel object
      logical,             intent(in)    :: is_periodic(3)                     !! Periodicity
      integer,             intent(in)    :: Ng(3)                              !! # Grid points
      integer,optional,    intent(in)    :: Nb(3)                              !! Explicit block decomposition
      ! Work variables
      integer :: dims(3)
      integer :: coords(3)
      logical :: reorder
      integer :: ndims
      integer :: ierr

      ! Initialize the dims var
      dims  = 0
      ndims = 3
      if (present(Nb)) then
        this%np = Nb
        dims    = Nb
      else
        ! Enforce 1 directional node for 1D & 2D simulations
        if (Ng(1).eq.1) dims(1)=1
        if (Ng(2).eq.1) dims(2)=1
        if (Ng(3).eq.1) dims(3)=1
        call MPI_DIMS_CREATE(this%nproc,ndims,dims,ierr)
        this%np  = dims
      end if
      this%npx = this%np(1)
      this%npy = this%np(2)
      this%npz = this%np(3)

      ! Test if nproc is correct
      if (this%nproc .ne. this%npx*this%npy*this%npz) &
        call this%Stop('Error setting up the parallel topology: mismatch in parallel decomposition')

      ! Set MPI topology
      reorder = .true.
      call MPI_CART_CREATE(this%comm%w,ndims,dims,is_periodic,reorder,this%comm%g,ierr)

      ! Return the cartesian coordinate of this rank
      call MPI_COMM_RANK(this%comm%g,this%rank%mine,ierr)
      call MPI_CART_COORDS(this%comm%g,this%rank%mine,ndims,coords,ierr)
      this%rank%mine   = this%rank%mine + 1
      this%rank%dir(1) = coords(1) + 1
      this%rank%dir(2) = coords(2) + 1
      this%rank%dir(3) = coords(3) + 1

                                                                               !! Define a root processor at coordinates (1,1,1)
      coords = [1,1,1]
      call MPI_CART_RANK(this%comm%g,coords-1,this%RootRank,ierr)
      this%RootRank  = this%RootRank + 1

      ! Look up the ranks for the neighbors
      call MPI_CART_SHIFT(this%comm%g,0,1,this%rank%L(1),this%rank%R(1),ierr)
      this%rank%L(1)=this%rank%L(1)+1
      this%rank%R(1)=this%rank%R(1)+1

      call MPI_CART_SHIFT(this%comm%g,1, 1,this%rank%L(2),this%rank%R(2),ierr)
      this%rank%L(2)=this%rank%L(2)+1
      this%rank%R(2)=this%rank%R(2)+1

      call MPI_CART_SHIFT(this%comm%g,2, 1,this%rank%L(3),this%rank%R(3),ierr)
      this%rank%L(3)=this%rank%L(3)+1
      this%rank%R(3)=this%rank%R(3)+1

      return
    end subroutine parallel_obj_Topology