diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index 9943379784..a552b9d000 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -57,7 +57,8 @@ contains
                             FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT,      &
                             OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, ONE
     use error,        only: fatal_error
-    use global,       only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes
+    use global,       only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,&
+                            matching_bins
     use mesh,         only: mesh_indices_to_bin
     use mesh_header,  only: StructuredMesh
     use tally_header, only: TallyObject
@@ -137,21 +138,21 @@ contains
               TALLY: if (ital == 1) then
 
                 ! reset all bins to 1
-                t % matching_bins = 1
+                matching_bins(1:t%n_filters) = 1
 
                 ! set ijk as mesh indices
                 ijk = (/ i, j, k /)
 
                 ! get bin number for mesh indices
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
 
                 ! apply energy in filter
                 if (i_filter_ein > 0) then
-                  t % matching_bins(i_filter_ein) = ng - h + 1
+                  matching_bins(i_filter_ein) = ng - h + 1
                 end if
 
                 ! calculate score index from bins
-                score_index = sum((t % matching_bins - 1) * t%stride) + 1
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
 
                 ! get flux
                 flux = t % results(1,score_index) % sum
@@ -180,24 +181,24 @@ contains
                 INGROUP: do g = 1, ng
 
                   ! reset all bins to 1
-                  t % matching_bins = 1
+                  matching_bins(1:t%n_filters) = 1
 
                   ! set ijk as mesh indices
                   ijk = (/ i, j, k /)
 
                   ! get bin number for mesh indices
-                  t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
+                  matching_bins(i_filter_mesh) = mesh_indices_to_bin(m,ijk)
 
                   if (i_filter_ein > 0) then
                     ! apply energy in filter
-                    t % matching_bins(i_filter_ein) = ng - h + 1
+                    matching_bins(i_filter_ein) = ng - h + 1
 
                     ! set energy out bin
-                    t % matching_bins(i_filter_eout) = ng - g + 1
+                    matching_bins(i_filter_eout) = ng - g + 1
                   end if
 
                   ! calculate score index from bins
-                  score_index = sum((t % matching_bins - 1) * t%stride) + 1
+                  score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
 
                   ! get scattering
                   cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
@@ -216,69 +217,69 @@ contains
               else if (ital == 3) then
 
                 ! initialize and filter for energy
-                t % matching_bins = 1
+                matching_bins(1:t%n_filters) = 1
                 if (i_filter_ein > 0) then
-                  t % matching_bins(i_filter_ein) = ng - h + 1
+                  matching_bins(i_filter_ein) = ng - h + 1
                 end if
 
                 ! left surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i-1, j, k /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_RIGHT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = IN_RIGHT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_RIGHT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = OUT_RIGHT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
 
                 ! right surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i, j, k /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_RIGHT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = IN_RIGHT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_RIGHT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = OUT_RIGHT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
 
                 ! back surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i, j-1, k /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_FRONT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = IN_FRONT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_FRONT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = OUT_FRONT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
 
                 ! front surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i, j, k /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_FRONT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = IN_FRONT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_FRONT
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = OUT_FRONT
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
 
                 ! bottom surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i, j, k-1 /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_TOP
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = IN_TOP
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_TOP
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = OUT_TOP
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
 
                 ! top surface
-                t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
+                matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
                      (/ i, j, k /) + 1, .true.)
-                t % matching_bins(i_filter_surf) = IN_TOP
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
+                matching_bins(i_filter_surf) = IN_TOP
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
                 cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
-                t % matching_bins(i_filter_surf) = OUT_TOP
-                score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
+                matching_bins(i_filter_surf) = OUT_TOP
+                score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
                 cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
 
               end if TALLY
diff --git a/src/global.F90 b/src/global.F90
index a84e9b00b8..dd436fefc2 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -97,6 +97,7 @@ module global
 
   type(StructuredMesh), allocatable, target :: meshes(:)
   type(TallyObject),    allocatable, target :: tallies(:)
+  integer, allocatable :: matching_bins(:)
 
   ! Pointers for different tallies
   type(TallyObject), pointer :: user_tallies(:) => null()
@@ -384,7 +385,7 @@ module global
   logical :: output_tallies = .true.
 
 !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, &
-!$omp&              trace, thread_id, current_work)
+!$omp&              trace, thread_id, current_work, matching_bins)
 
 contains
 
diff --git a/src/output.F90 b/src/output.F90
index 4e1d903b8e..472cb5574b 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -1644,7 +1644,7 @@ contains
       ! to be used for a given tally.
 
       ! Initialize bins, filter level, and indentation
-      t % matching_bins = 0
+      matching_bins(1:t%n_filters) = 0
       j = 1
       indent = 0
 
@@ -1654,16 +1654,16 @@ contains
           if (t % n_filters == 0) exit find_bin
 
           ! Increment bin combination
-          t % matching_bins(j) = t % matching_bins(j) + 1
+          matching_bins(j) = matching_bins(j) + 1
 
           ! =================================================================
           ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
 
-          if (t % matching_bins(j) > t % filters(j) % n_bins) then
+          if (matching_bins(j) > t % filters(j) % n_bins) then
             ! If this is the first filter, then exit
             if (j == 1) exit print_bin
 
-            t % matching_bins(j) = 0
+            matching_bins(j) = 0
             j = j - 1
             indent = indent - 2
 
@@ -1696,7 +1696,7 @@ contains
         ! bins below the lowest filter level will be zeros
 
         if (t % n_filters > 0) then
-          filter_index = sum((max(t % matching_bins,1) - 1) * t % stride) + 1
+          filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1
         else
           filter_index = 1
         end if
@@ -1805,7 +1805,7 @@ contains
     m => meshes(t % filters(i_filter_mesh) % int_bins(1))
 
     ! initialize bins array
-    t % matching_bins = 1
+    matching_bins(1:t%n_filters) = 1
 
     ! determine how many energy in bins there are
     i_filter_ein = t % find_filter(FILTER_ENERGYIN)
@@ -1831,7 +1831,7 @@ contains
           do l = 1, n
             if (print_ebin) then
               ! Set incoming energy bin
-              t % matching_bins(i_filter_ein) = l
+              matching_bins(i_filter_ein) = l
 
               ! Write incoming energy bin
               write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
@@ -1839,102 +1839,102 @@ contains
             end if
 
             ! Left Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_RIGHT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_RIGHT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Left", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_RIGHT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_RIGHT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Left", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             ! Right Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_RIGHT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_RIGHT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Right", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_RIGHT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_RIGHT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Right", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             ! Back Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_FRONT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_FRONT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Back", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_FRONT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_FRONT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Back", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             ! Front Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_FRONT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_FRONT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Front", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_FRONT
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_FRONT
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Front", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             ! Bottom Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_TOP
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_TOP
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Bottom", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_TOP
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_TOP
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Bottom", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             ! Top Surface
-            t % matching_bins(i_filter_mesh) = &
+            matching_bins(i_filter_mesh) = &
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
-            t % matching_bins(i_filter_surf) = IN_TOP
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = IN_TOP
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Incoming Current from Top", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
-            t % matching_bins(i_filter_surf) = OUT_TOP
-            filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+            matching_bins(i_filter_surf) = OUT_TOP
+            filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
             write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
                  "Outgoing Current to Top", &
                  to_str(t % results(1,filter_index) % sum), &
@@ -1965,7 +1965,7 @@ contains
     real(8)              :: E1     ! upper bound for energy bin
     type(StructuredMesh), pointer :: m => null()
 
-    bin = t % matching_bins(i_filter)
+    bin = matching_bins(i_filter)
 
     select case(t % filters(i_filter) % type)
     case (FILTER_UNIVERSE)
diff --git a/src/tally.F90 b/src/tally.F90
index 2fafacf320..e1fcc2ff06 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -85,7 +85,7 @@ contains
       ! be accumulating the tally values
 
       ! Determine scoring index for this filter combination
-      filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+      filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
       ! Check for nuclide bins
       k = 0
@@ -388,7 +388,7 @@ contains
 
     ! save original outgoing energy bin and score index
     i = t % find_filter(FILTER_ENERGYOUT)
-    bin_energyout = t % matching_bins(i)
+    bin_energyout = matching_bins(i)
 
     ! Get number of energies on filter
     n = size(t % filters(i) % real_bins)
@@ -411,10 +411,10 @@ contains
            E_out > t % filters(i) % real_bins(n)) cycle
 
       ! change outgoing energy bin
-      t % matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
+      matching_bins(i) = binary_search(t % filters(i) % real_bins, n, E_out)
 
       ! determine scoring index
-      i_filter = sum((t % matching_bins - 1) * t % stride) + 1
+      i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
       ! Add score to tally
 !$omp critical
@@ -424,7 +424,7 @@ contains
     end do
 
     ! reset outgoing energy bin and score index
-    t % matching_bins(i) = bin_energyout
+    matching_bins(i) = bin_energyout
 
   end subroutine score_fission_eout
 
@@ -495,7 +495,7 @@ contains
       ! be accumulating the tally values
 
       ! Determine scoring index for this filter combination
-      filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+      filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
       if (t % all_nuclides) then
         call score_all_nuclides(p, i_tally, flux, filter_index)
@@ -1005,7 +1005,7 @@ contains
 !$omp threadprivate(t, m, mat, coord)
 
     t => tallies(i_tally)
-    t % matching_bins = 1
+    matching_bins(1:t%n_filters) = 1
 
     ! ==========================================================================
     ! CHECK IF THIS TRACK INTERSECTS THE MESH
@@ -1045,11 +1045,11 @@ contains
       case (FILTER_UNIVERSE)
         ! determine next universe bin
         ! TODO: Account for multiple universes when performing this filter
-        t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
+        matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
              p % coord % universe, i_tally)
 
       case (FILTER_MATERIAL)
-        t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
+        matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
              p % material, i_tally)
 
       case (FILTER_CELL)
@@ -1057,21 +1057,21 @@ contains
         coord => p % coord0
         do while(associated(coord))
           position(FILTER_CELL) = 0
-          t % matching_bins(i) = get_next_bin(FILTER_CELL, &
+          matching_bins(i) = get_next_bin(FILTER_CELL, &
                coord % cell, i_tally)
-          if (t % matching_bins(i) /= NO_BIN_FOUND) exit
+          if (matching_bins(i) /= NO_BIN_FOUND) exit
           coord => coord % next
         end do
         nullify(coord)
 
       case (FILTER_CELLBORN)
         ! determine next cellborn bin
-        t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
+        matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
              p % cell_born, i_tally)
 
       case (FILTER_SURFACE)
         ! determine next surface bin
-        t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
+        matching_bins(i) = get_next_bin(FILTER_SURFACE, &
              p % surface, i_tally)
 
       case (FILTER_ENERGYIN)
@@ -1081,17 +1081,17 @@ contains
         ! check if energy of the particle is within energy bins
         if (p % E < t % filters(i) % real_bins(1) .or. &
              p % E > t % filters(i) % real_bins(k + 1)) then
-          t % matching_bins(i) = NO_BIN_FOUND
+          matching_bins(i) = NO_BIN_FOUND
         else
           ! search to find incoming energy bin
-          t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
+          matching_bins(i) = binary_search(t % filters(i) % real_bins, &
                k + 1, p % E)
         end if
 
       end select
 
       ! Check if no matching bin was found
-      if (t % matching_bins(i) == NO_BIN_FOUND) return
+      if (matching_bins(i) == NO_BIN_FOUND) return
 
     end do FILTER_LOOP
 
@@ -1161,10 +1161,10 @@ contains
         flux = p % wgt * distance
 
         ! Determine mesh bin
-        t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
+        matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, ijk_cross)
 
         ! Determining scoring index
-        filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+        filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
         if (t % all_nuclides) then
           ! Score reaction rates for each nuclide in material
@@ -1301,7 +1301,7 @@ contains
 
     found_bin = .true.
     t => tallies(i_tally)
-    t % matching_bins = 1
+    matching_bins(1:t%n_filters) = 1
 
     FILTER_LOOP: do i = 1, t % n_filters
 
@@ -1311,16 +1311,16 @@ contains
         m => meshes(t % filters(i) % int_bins(1))
 
         ! Determine if we're in the mesh first
-        call get_mesh_bin(m, p % coord0 % xyz, t % matching_bins(i))
+        call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i))
 
       case (FILTER_UNIVERSE)
         ! determine next universe bin
         ! TODO: Account for multiple universes when performing this filter
-        t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
+        matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
              p % coord % universe, i_tally)
 
       case (FILTER_MATERIAL)
-        t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
+        matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
              p % material, i_tally)
 
       case (FILTER_CELL)
@@ -1328,21 +1328,21 @@ contains
         coord => p % coord0
         do while(associated(coord))
           position(FILTER_CELL) = 0
-          t % matching_bins(i) = get_next_bin(FILTER_CELL, &
+          matching_bins(i) = get_next_bin(FILTER_CELL, &
                coord % cell, i_tally)
-          if (t % matching_bins(i) /= NO_BIN_FOUND) exit
+          if (matching_bins(i) /= NO_BIN_FOUND) exit
           coord => coord % next
         end do
         nullify(coord)
 
       case (FILTER_CELLBORN)
         ! determine next cellborn bin
-        t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
+        matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
              p % cell_born, i_tally)
 
       case (FILTER_SURFACE)
         ! determine next surface bin
-        t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
+        matching_bins(i) = get_next_bin(FILTER_SURFACE, &
              p % surface, i_tally)
 
       case (FILTER_ENERGYIN)
@@ -1359,10 +1359,10 @@ contains
         ! check if energy of the particle is within energy bins
         if (E < t % filters(i) % real_bins(1) .or. &
              E > t % filters(i) % real_bins(n + 1)) then
-          t % matching_bins(i) = NO_BIN_FOUND
+          matching_bins(i) = NO_BIN_FOUND
         else
           ! search to find incoming energy bin
-          t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
+          matching_bins(i) = binary_search(t % filters(i) % real_bins, &
                n + 1, E)
         end if
 
@@ -1373,17 +1373,17 @@ contains
         ! check if energy of the particle is within energy bins
         if (p % E < t % filters(i) % real_bins(1) .or. &
              p % E > t % filters(i) % real_bins(n + 1)) then
-          t % matching_bins(i) = NO_BIN_FOUND
+          matching_bins(i) = NO_BIN_FOUND
         else
           ! search to find incoming energy bin
-          t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
+          matching_bins(i) = binary_search(t % filters(i) % real_bins, &
                n + 1, p % E)
         end if
 
       end select
 
       ! If the current filter didn't match, exit this subroutine
-      if (t % matching_bins(i) == NO_BIN_FOUND) then
+      if (matching_bins(i) == NO_BIN_FOUND) then
         found_bin = .false.
         return
       end if
@@ -1475,7 +1475,7 @@ contains
         end if
 
         ! search to find incoming energy bin
-        t % matching_bins(j) = binary_search(t % filters(j) % real_bins, &
+        matching_bins(j) = binary_search(t % filters(j) % real_bins, &
              n + 1, p % E)
       end if
 
@@ -1492,10 +1492,10 @@ contains
           do j = ijk0(3), ijk1(3) - 1
             ijk0(3) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_TOP
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_TOP
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1506,10 +1506,10 @@ contains
           do j = ijk0(3) - 1, ijk1(3), -1
             ijk0(3) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_TOP
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_TOP
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1524,10 +1524,10 @@ contains
           do j = ijk0(2), ijk1(2) - 1
             ijk0(2) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_FRONT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_FRONT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1538,10 +1538,10 @@ contains
           do j = ijk0(2) - 1, ijk1(2), -1
             ijk0(2) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_FRONT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_FRONT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1556,10 +1556,10 @@ contains
           do j = ijk0(1), ijk1(1) - 1
             ijk0(1) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_RIGHT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_RIGHT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1570,10 +1570,10 @@ contains
           do j = ijk0(1) - 1, ijk1(1), -1
             ijk0(1) = j
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_RIGHT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_RIGHT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
-              filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+              filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 !$omp critical
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
@@ -1598,7 +1598,7 @@ contains
 
       do k = 1, n_cross
         ! Reset scoring bin index
-        t % matching_bins(i_filter_surf) = 0
+        matching_bins(i_filter_surf) = 0
 
         ! Calculate distance to each bounding surface. We need to treat
         ! special case where the cosine of the angle is zero since this would
@@ -1625,8 +1625,8 @@ contains
             ! Crossing into right mesh cell -- this is treated as outgoing
             ! current from (i,j,k)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_RIGHT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_RIGHT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
             ijk0(1) = ijk0(1) + 1
@@ -1637,8 +1637,8 @@ contains
             ijk0(1) = ijk0(1) - 1
             xyz_cross(1) = xyz_cross(1) - m % width(1)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_RIGHT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_RIGHT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
           end if
@@ -1647,8 +1647,8 @@ contains
             ! Crossing into front mesh cell -- this is treated as outgoing
             ! current in (i,j,k)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_FRONT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_FRONT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
             ijk0(2) = ijk0(2) + 1
@@ -1659,8 +1659,8 @@ contains
             ijk0(2) = ijk0(2) - 1
             xyz_cross(2) = xyz_cross(2) - m % width(2)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_FRONT
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_FRONT
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
           end if
@@ -1669,8 +1669,8 @@ contains
             ! Crossing into top mesh cell -- this is treated as outgoing
             ! current in (i,j,k)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = OUT_TOP
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = OUT_TOP
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
             ijk0(3) = ijk0(3) + 1
@@ -1681,16 +1681,16 @@ contains
             ijk0(3) = ijk0(3) - 1
             xyz_cross(3) = xyz_cross(3) - m % width(3)
             if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
-              t % matching_bins(i_filter_surf) = IN_TOP
-              t % matching_bins(i_filter_mesh) = &
+              matching_bins(i_filter_surf) = IN_TOP
+              matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
             end if
           end if
         end if
 
         ! Determine scoring index
-        if (t % matching_bins(i_filter_surf) > 0) then
-          filter_index = sum((t % matching_bins - 1) * t % stride) + 1
+        if (matching_bins(i_filter_surf) > 0) then
+          filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
           ! Check for errors
           if (filter_index <= 0 .or. filter_index > &
diff --git a/src/tally_header.F90 b/src/tally_header.F90
index 41c5d79715..17cb1cb0de 100644
--- a/src/tally_header.F90
+++ b/src/tally_header.F90
@@ -86,7 +86,6 @@ module tally_header
     ! mapped onto one dimension in the results array, the stride attribute gives
     ! the stride for a given filter type within the results array
 
-    integer, allocatable :: matching_bins(:)
     integer, allocatable :: stride(:)
 
     ! This array provides a way to lookup what index in the filters array a
@@ -169,8 +168,6 @@ module tally_header
         deallocate(this % filters)
       end if
       
-      if (allocated(this % matching_bins)) &
-           deallocate(this % matching_bins)
       if (allocated(this % stride)) &
            deallocate(this % stride)
       
diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90
index 04928499e7..23ab5b5634 100644
--- a/src/tally_initialize.F90
+++ b/src/tally_initialize.F90
@@ -31,9 +31,10 @@ contains
 
   subroutine setup_tally_arrays()
 
-    integer :: i ! loop index for tallies
-    integer :: j ! loop index for filters
-    integer :: n ! temporary stride
+    integer :: i                 ! loop index for tallies
+    integer :: j                 ! loop index for filters
+    integer :: n                 ! temporary stride
+    integer :: max_n_filters = 0 ! maximum number of filters
     type(TallyObject), pointer :: t => null()
 
     TALLY_LOOP: do i = 1, n_tallies
@@ -42,7 +43,7 @@ contains
 
       ! Allocate stride and matching_bins arrays
       allocate(t % stride(t % n_filters))
-      allocate(t % matching_bins(t % n_filters))
+      max_n_filters = max(max_n_filters, t % n_filters)
 
       ! The filters are traversed in opposite order so that the last filter has
       ! the shortest stride in memory and the first filter has the largest
@@ -63,6 +64,11 @@ contains
 
     end do TALLY_LOOP
 
+    ! Allocate array for matching filter bins
+!$omp parallel
+    allocate(matching_bins(max_n_filters))
+!$omp end parallel
+
   end subroutine setup_tally_arrays
 
 !===============================================================================