diff --git a/src/geometry.F90 b/src/geometry.F90
index 2d9f3355a8..c79f0e6119 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -11,7 +11,7 @@ module geometry
   use tally,                  only: score_surface_current
 
   implicit none
- 
+
 contains
 
 !===============================================================================
@@ -265,13 +265,13 @@ contains
                 lattice_edge = .true.
               end if
             end if
-            
+
             if (lattice_edge) then
-              
+
               ! In this case the neutron is leaving the lattice, so we move it
               ! out, remove all lower coordinate levels and then search from
               ! universe 0.
-            
+
               p % coord => p % coord0
               call deallocate_coord(p % coord % next)
 
@@ -288,7 +288,7 @@ contains
               p % last_material = p % material
               p % material = c % material
 
-              ! We'll still make a new coordinate for the particle, as 
+              ! We'll still make a new coordinate for the particle, as
               ! distance_to_boundary will still need to track through lattice
               ! widths even though there's nothing in them but this material
 
@@ -407,10 +407,9 @@ contains
 
       ! Score to global leakage tally
       if (tallies_on) then
-!$omp critical
+!$omp atomic
         global_tallies(LEAKAGE) % value = &
            global_tallies(LEAKAGE) % value + p % wgt
-!$omp end critical
       end if
 
       ! Display message
@@ -645,7 +644,7 @@ contains
         return
       end if
     end if
-       
+
   end subroutine cross_surface
 
 !===============================================================================
@@ -907,7 +906,7 @@ contains
             if (quad < ZERO) then
               ! no intersection with cylinder
 
-              d = INFINITY 
+              d = INFINITY
 
             elseif (on_surface) then
               ! particle is on the cylinder, thus one distance is
@@ -956,7 +955,7 @@ contains
             if (quad < ZERO) then
               ! no intersection with cylinder
 
-              d = INFINITY 
+              d = INFINITY
 
             elseif (on_surface) then
               ! particle is on the cylinder, thus one distance is
@@ -1005,7 +1004,7 @@ contains
             if (quad < ZERO) then
               ! no intersection with cylinder
 
-              d = INFINITY 
+              d = INFINITY
 
             elseif (on_surface) then
               ! particle is on the cylinder, thus one distance is
@@ -1052,7 +1051,7 @@ contains
           if (quad < ZERO) then
             ! no intersection with sphere
 
-            d = INFINITY 
+            d = INFINITY
 
           elseif (on_surface) then
             ! particle is on the sphere, thus one distance is
@@ -1099,7 +1098,7 @@ contains
           if (quad < ZERO) then
             ! no intersection with cone
 
-            d = INFINITY 
+            d = INFINITY
 
           elseif (on_surface) then
             ! particle is on the cone, thus one distance is positive/negative
@@ -1118,7 +1117,7 @@ contains
             d = (-k - quad)/a
             b = (-k + quad)/a
 
-            ! determine the smallest positive solution 
+            ! determine the smallest positive solution
             if (d < ZERO) then
               if (b > ZERO) then
                 d = b
@@ -1148,7 +1147,7 @@ contains
           if (quad < ZERO) then
             ! no intersection with cone
 
-            d = INFINITY 
+            d = INFINITY
 
           elseif (on_surface) then
             ! particle is on the cone, thus one distance is positive/negative
@@ -1167,7 +1166,7 @@ contains
             d = (-k - quad)/a
             b = (-k + quad)/a
 
-            ! determine the smallest positive solution 
+            ! determine the smallest positive solution
             if (d < ZERO) then
               if (b > ZERO) then
                 d = b
@@ -1197,7 +1196,7 @@ contains
           if (quad < ZERO) then
             ! no intersection with cone
 
-            d = INFINITY 
+            d = INFINITY
 
           elseif (on_surface) then
             ! particle is on the cone, thus one distance is positive/negative
@@ -1216,7 +1215,7 @@ contains
             d = (-k - quad)/a
             b = (-k + quad)/a
 
-            ! determine the smallest positive solution 
+            ! determine the smallest positive solution
             if (d < ZERO) then
               if (b > ZERO) then
                 d = b
@@ -1274,7 +1273,7 @@ contains
           ! logic here checks whether the relative difference is within floating
           ! point precision.
 
-          if (d < dist) then 
+          if (d < dist) then
             if (abs(d - dist)/dist >= FP_REL_PRECISION) then
               dist = d
               if (u > 0) then
@@ -1566,9 +1565,8 @@ contains
 
     ! Increment number of lost particles
     p % alive = .false.
-!$omp critical
+!$omp atomic
     n_lost_particles = n_lost_particles + 1
-!$omp end critical
 
     ! Abort the simulation if the maximum number of lost particles has been
     ! reached
diff --git a/src/global.F90 b/src/global.F90
index f4c50b8330..874b8d0d1a 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -62,7 +62,7 @@ module global
   ! Cross section arrays
   type(Nuclide),    allocatable, target :: nuclides(:)    ! Nuclide cross-sections
   type(SAlphaBeta), allocatable, target :: sab_tables(:)  ! S(a,b) tables
-  type(XsListing),  allocatable, target :: xs_listings(:) ! cross_sections.xml listings 
+  type(XsListing),  allocatable, target :: xs_listings(:) ! cross_sections.xml listings
 
   ! Cross section caches
   type(NuclideMicroXS), allocatable :: micro_xs(:)  ! Cache for each nuclide
@@ -119,7 +119,7 @@ module global
   !   2) track-length estimate of k-eff
   !   3) leakage fraction
 
-  type(TallyResult), target :: global_tallies(N_GLOBAL_TALLIES)
+  type(TallyResult), allocatable, target :: global_tallies(:)
 
   ! Tally map structure
   type(TallyMap), allocatable :: tally_maps(:)
@@ -300,7 +300,7 @@ module global
   logical :: write_initial_source = .false.
 
   ! ============================================================================
-  ! CMFD VARIABLES 
+  ! CMFD VARIABLES
 
   ! Main object
   type(cmfd_type) :: cmfd
@@ -310,11 +310,11 @@ module global
 
   ! CMFD communicator
   integer :: cmfd_comm
- 
+
   ! Timing objects
   type(Timer) :: time_cmfd      ! timer for whole cmfd calculation
   type(Timer) :: time_cmfdbuild ! timer for matrix build
-  type(Timer) :: time_cmfdsolve ! timer for solver 
+  type(Timer) :: time_cmfdsolve ! timer for solver
 
   ! Flag for active core map
   logical :: cmfd_coremap = .false.
@@ -390,7 +390,7 @@ module global
   ! RESONANCE SCATTERING VARIABLES
 
   logical :: treat_res_scat = .false. ! is resonance scattering treated?
-  integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers 
+  integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers
   type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info
 
 !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
@@ -399,14 +399,14 @@ module global
 contains
 
 !===============================================================================
-! FREE_MEMORY deallocates and clears  all global allocatable arrays in the 
+! FREE_MEMORY deallocates and clears  all global allocatable arrays in the
 ! program
 !===============================================================================
 
   subroutine free_memory()
-    
+
     integer :: i ! Loop Index
-    
+
     ! Deallocate cells, surfaces, materials
     if (allocated(cells)) deallocate(cells)
     if (allocated(universes)) deallocate(universes)
@@ -449,6 +449,7 @@ contains
     if (allocated(entropy_p)) deallocate(entropy_p)
 
     ! Deallocate tally-related arrays
+    if (allocated(global_tallies)) deallocate(global_tallies)
     if (allocated(meshes)) deallocate(meshes)
     if (allocated(tallies)) then
     ! First call the clear routines
@@ -488,7 +489,7 @@ contains
 
     ! Deallocate track_identifiers
     if (allocated(track_identifiers)) deallocate(track_identifiers)
-    
+
     ! Deallocate dictionaries
     call cell_dict % clear()
     call universe_dict % clear()
@@ -525,7 +526,7 @@ contains
         if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width)
         deallocate(ufs_mesh)
     end if
-    
+
   end subroutine free_memory
 
 end module global
diff --git a/src/physics.F90 b/src/physics.F90
index 8fa0575cdc..7988970ef4 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -256,22 +256,19 @@ contains
       p % last_wgt = p % wgt
 
       ! Score implicit absorption estimate of keff
-!$omp critical
+!$omp atomic
       global_tallies(K_ABSORPTION) % value = &
            global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
            micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
-!$omp end critical
-
     else
       ! See if disappearance reaction happens
       if (micro_xs(i_nuclide) % absorption > &
            prn() * micro_xs(i_nuclide) % total) then
         ! Score absorption estimate of keff
-!$omp critical
+!$omp atomic
         global_tallies(K_ABSORPTION) % value = &
              global_tallies(K_ABSORPTION) % value + p % wgt * &
              micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
-!$omp end critical
 
         p % alive = .false.
         p % event = EVENT_ABSORB
@@ -797,7 +794,7 @@ contains
         sampling_scheme = 'cxs'
       end if
 
-    ! otherwise, use free gas model  
+    ! otherwise, use free gas model
     else
       if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then
         v_target = ZERO
@@ -859,7 +856,7 @@ contains
       m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
        & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
       xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
-      
+
       ! get max 0K xs value over range of practical relative energies
       xs_max = max(xs_low, &
         & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
@@ -972,7 +969,7 @@ contains
     case default
       call fatal_error("Not a recognized resonance scattering treatment!")
     end select
-    
+
   end subroutine sample_target_velocity
 
 !===============================================================================
diff --git a/src/tally.F90 b/src/tally.F90
index 87f5a2cab7..7c7747ff51 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -281,10 +281,9 @@ contains
               ! get the score and tally it
               score = last_wgt * calc_pn(n, mu)
 
-!$omp critical
+!$omp atomic
               t % results(score_index, filter_index) % value = &
                 t % results(score_index, filter_index) % value + score
-!$omp end critical
             end do
             j = j + t % moment_order(j)
             cycle SCORE_LOOP
@@ -347,10 +346,9 @@ contains
               ! get the score and tally it
               score = wgt * calc_pn(n, mu)
 
-!$omp critical
+!$omp atomic
               t % results(score_index, filter_index) % value = &
                 t % results(score_index, filter_index) % value + score
-!$omp end critical
             end do
             j = j + t % moment_order(j)
             cycle SCORE_LOOP
@@ -542,10 +540,9 @@ contains
           end select
 
           ! Add score to tally
-!$omp critical
+!$omp atomic
           t % results(score_index, filter_index) % value = &
                t % results(score_index, filter_index) % value + score
-!$omp end critical
 
         end do SCORE_LOOP
 
@@ -617,10 +614,9 @@ contains
       i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
       ! Add score to tally
-!$omp critical
+!$omp atomic
       t % results(i_score, i_filter) % value = &
            t % results(i_score, i_filter) % value + score
-!$omp end critical
     end do
 
     ! reset outgoing energy bin and score index
@@ -1015,10 +1011,9 @@ contains
             end if
 
             ! Add score to tally
-!$omp critical
+!$omp atomic
             t % results(score_index, filter_index) % value = &
                  t % results(score_index, filter_index) % value + score
-!$omp end critical
 
           end do SCORE_LOOP
 
@@ -1214,10 +1209,9 @@ contains
         end select
 
         ! Add score to tally
-!$omp critical
+!$omp atomic
         t % results(score_index, filter_index) % value = &
              t % results(score_index, filter_index) % value + score
-!$omp end critical
 
       end do SCORE_LOOP
 
@@ -1366,10 +1360,9 @@ contains
       end select
 
       ! Add score to tally
-!$omp critical
+!$omp atomic
       t % results(score_index, filter_index) % value = &
            t % results(score_index, filter_index) % value + score
-!$omp end critical
 
     end do MATERIAL_SCORE_LOOP
 
@@ -1787,10 +1780,9 @@ contains
               end if
 
               ! Add score to tally
-!$omp critical
+!$omp atomic
               t % results(score_index, filter_index) % value = &
                    t % results(score_index, filter_index) % value + score
-!$omp end critical
 
             end do SCORE_LOOP
 
@@ -2021,10 +2013,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         else
@@ -2035,10 +2026,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         end if
@@ -2053,10 +2043,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         else
@@ -2067,10 +2056,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         end if
@@ -2085,10 +2073,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         else
@@ -2099,10 +2086,9 @@ contains
               matching_bins(i_filter_mesh) = &
                    mesh_indices_to_bin(m, ijk0 + 1, .true.)
               filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-!$omp critical
+!$omp atomic
               t % results(1, filter_index) % value = &
                    t % results(1, filter_index) % value + p % wgt
-!$omp end critical
             end if
           end do
         end if
@@ -2224,10 +2210,9 @@ contains
           end if
 
           ! Add to surface current tally
-!$omp critical
+!$omp atomic
           t % results(1, filter_index) % value = &
                t % results(1, filter_index) % value + p % wgt
-!$omp end critical
         end if
 
         ! Calculate new coordinates
diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90
index 23ab5b5634..b7dc5ed98a 100644
--- a/src/tally_initialize.F90
+++ b/src/tally_initialize.F90
@@ -19,6 +19,9 @@ contains
 
   subroutine configure_tallies()
 
+    ! Allocate global tallies
+    allocate(global_tallies(N_GLOBAL_TALLIES))
+
     call setup_tally_arrays()
     call setup_tally_maps()
 
diff --git a/src/tracking.F90 b/src/tracking.F90
index d108325fa0..8eca3871c7 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -62,9 +62,8 @@ contains
     n_event = 0
 
     ! Add paricle's starting weight to count for normalizing tallies later
-!$omp critical
+!$omp atomic
     total_weight = total_weight + p % wgt
-!$omp end critical
 
     ! Force calculation of cross-sections by setting last energy to zero
     micro_xs % last_E = ZERO
@@ -112,11 +111,10 @@ contains
            call score_tracklength_tally(p, distance)
 
       ! Score track-length estimate of k-eff
-!$omp critical
+!$omp atomic
       global_tallies(K_TRACKLENGTH) % value = &
            global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
            material_xs % nu_fission
-!$omp end critical
 
       if (d_collision > d_boundary) then
         ! ====================================================================
@@ -140,11 +138,10 @@ contains
         ! PARTICLE HAS COLLISION
 
         ! Score collision estimate of keff
-!$omp critical
+!$omp atomic
         global_tallies(K_COLLISION) % value = &
              global_tallies(K_COLLISION) % value + p % wgt * &
              material_xs % nu_fission / material_xs % total
-!$omp end critical
 
         ! score surface current tallies -- this has to be done before the collision
         ! since the direction of the particle will change and we need to use the