diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index aee3c5c25..dc8681aee 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -568,6 +568,10 @@ contains
     ! count how many fuel assemblies exist
     cmfd % mat_dim = sum(cmfd % coremap - 1)
 
+    ! allocate indexmap
+    if (.not. allocated(cmfd % indexmap)) allocate                             &
+   &                                      (cmfd % indexmap(cmfd % mat_dim,3))
+
     ! begin loops over spatial indices
     ZLOOP: do k = 1,nz
 
@@ -585,6 +589,9 @@ contains
 
             ! must be a fuel --> give unique id number
             cmfd % coremap(i,j,k) = kount
+            cmfd % indexmap(kount,1) = i
+            cmfd % indexmap(kount,2) = j
+            cmfd % indexmap(kount,3) = k 
             kount = kount + 1
 
           end if
diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index 65d1d93d4..bf9bed0e7 100644
--- a/src/cmfd_execute.F90
+++ b/src/cmfd_execute.F90
@@ -42,7 +42,7 @@ contains
 
       ! begin timer
 !     call timer_start(time_cmfd)
-print *,'Setting up data'
+
       ! set up cmfd
       if(.not. cmfd_only) call set_up_cmfd()
 
@@ -52,7 +52,6 @@ print *,'Setting up data'
     call cmfd_bcast()
 
 #ifdef PETSC
-print *,'Executing SLEPC'
       ! execute snes solver
       call cmfd_slepc_execute()
 #endif
@@ -80,6 +79,7 @@ print *,'Executing SLEPC'
   subroutine cmfd_bcast()
 
     use cmfd_header, only: allocate_cmfd
+    use global,      only: cmfd_coremap
 
     integer :: nx  ! number of mesh cells in x direction
     integer :: ny  ! number of mesh cells in y direction
@@ -111,6 +111,15 @@ print *,'Executing SLEPC'
     call MPI_BCAST(cmfd%hxyz,3*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
     call MPI_BCAST(cmfd%current,12*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
 
+    ! broadcast coremap info
+    if (cmfd_coremap) then
+      call MPI_BCAST(cmfd%coremap,nx*ny*nz,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
+      call MPI_BCAST(cmfd%mat_dim,1,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
+      if (.not. allocated(cmfd % indexmap)) allocate                           &
+     &           (cmfd % indexmap(cmfd % mat_dim,3))
+      call MPI_BCAST(cmfd%indexmap,cmfd%mat_dim*3,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
+    end if
+
     ! sync up procs
     call MPI_Barrier(MPI_COMM_WORLD,mpi_err)
 
diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90
index b9ac7ff84..58df1ee38 100644
--- a/src/cmfd_header.F90
+++ b/src/cmfd_header.F90
@@ -50,6 +50,7 @@ module cmfd_header
 
     ! core map for no reflector accel
     integer, allocatable :: coremap(:,:,:)
+    integer, allocatable :: indexmap(:,:)
     integer :: mat_dim
 
     ! eigenvector/eigenvalue from cmfd run
@@ -125,12 +126,12 @@ contains
     if (allocated(this % dhat))          deallocate(this % dhat)
     if (allocated(this % hxyz))          deallocate(this % hxyz)
     if (allocated(this % coremap))       deallocate(this % coremap)
+    if (allocated(this % indexmap))      deallocate(this % indexmap)
     if (allocated(this % phi))           deallocate(this % phi)
     if (allocated(this % sourcepdf))     deallocate(this % sourcepdf)
     if (allocated(this % sourcecounts))  deallocate(this % sourcecounts)
     if (allocated(this % weightfactors)) deallocate(this % weightfactors)
 
-
   end subroutine deallocate_cmfd
 
 end module cmfd_header
diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90
index 08c7f984b..d8aeb9519 100644
--- a/src/cmfd_loss_operator.F90
+++ b/src/cmfd_loss_operator.F90
@@ -281,17 +281,32 @@ contains
 
   subroutine matrix_to_indices(irow,g,i,j,k)
 
+    use global, only: cmfd,cmfd_coremap
+
     integer :: i                    ! iteration counter for x
     integer :: j                    ! iteration counter for y
     integer :: k                    ! iteration counter for z
     integer :: g                    ! iteration counter for groups
     integer :: irow                 ! iteration counter over row (0 reference)
 
-    ! compute indices
-    k = irow/(nx*ny*nz) + 1
-    j = mod(irow,nx*ny*nz)/(ny*nz) + 1
-    i = mod(irow,ny*nz)/nz + 1
-    g = mod(irow,nz) + 1 
+    ! check for core map
+    if (cmfd_coremap) then
+
+      ! get indices from indexmap
+      g = mod(irow,ng) + 1
+      i = cmfd % indexmap(irow+1,1)
+      j = cmfd % indexmap(irow+1,2)
+      k = cmfd % indexmap(irow+1,3)
+
+    else
+
+      ! compute indices
+      g = mod(irow,ng) + 1 
+      i = mod(irow,ng*nx)/ng + 1
+      j = mod(irow,ng*nx*ny)/(ng*nx)+ 1
+      k = mod(irow,ng*nx*ny*nz)/(ng*nx*ny) + 1 
+
+    end if
 
   end subroutine matrix_to_indices
 
diff --git a/src/cmfd_prod_operator.F90 b/src/cmfd_prod_operator.F90
index 33a120f6a..f5d8c45b5 100644
--- a/src/cmfd_prod_operator.F90
+++ b/src/cmfd_prod_operator.F90
@@ -175,17 +175,32 @@ contains
 
   subroutine matrix_to_indices(irow,g,i,j,k)
 
+    use global, only: cmfd,cmfd_coremap
+
     integer :: i                    ! iteration counter for x
     integer :: j                    ! iteration counter for y
     integer :: k                    ! iteration counter for z
     integer :: g                    ! iteration counter for groups
     integer :: irow                 ! iteration counter over row (0 reference)
 
-    ! compute indices
-    k = irow/(nx*ny*nz) + 1
-    j = mod(irow,nx*ny*nz)/(ny*nz) + 1
-    i = mod(irow,ny*nz)/nz + 1
-    g = mod(irow,nz) + 1
+    ! check for core map
+    if (cmfd_coremap) then
+
+      ! get indices from indexmap
+      g = mod(irow,ng) + 1
+      i = cmfd % indexmap(irow+1,1)
+      j = cmfd % indexmap(irow+1,2)
+      k = cmfd % indexmap(irow+1,3)
+
+    else
+
+      ! compute indices
+      g = mod(irow,ng) + 1 
+      i = mod(irow,ng*nx)/ng + 1
+      j = mod(irow,ng*nx*ny)/(ng*nx)+ 1
+      k = mod(irow,ng*nx*ny*nz)/(ng*nx*ny) + 1
+
+    end if
 
   end subroutine matrix_to_indices