diff --git a/src/ace_header.F90 b/src/ace_header.F90
index 7753f4734a..209b752bf8 100644
--- a/src/ace_header.F90
+++ b/src/ace_header.F90
@@ -1,6 +1,6 @@
 module ace_header
 
-  use constants,   only: MAX_FILE_LEN
+  use constants,   only: MAX_FILE_LEN, ZERO
   use endf_header, only: Tab1
   use list_header, only: ListInt
 
@@ -167,12 +167,12 @@ module ace_header
 
   type Nuclide0K
 
-    character(10) :: nuclide            ! name of nuclide, e.g. U-238
-    character(16) :: scheme = 'ares'    ! target velocity sampling scheme
-    character(10) :: name               ! name of nuclide, e.g. 92235.03c
-    character(10) :: name_0K            ! name of 0K nuclide, e.g. 92235.00c
-    real(8)       :: E_min = 0.01e-6    ! lower cutoff energy for res scattering
-    real(8)       :: E_max = 1000.0e-6  ! upper cutoff energy for res scattering
+    character(10) :: nuclide             ! name of nuclide, e.g. U-238
+    character(16) :: scheme = 'ares'     ! target velocity sampling scheme
+    character(10) :: name                ! name of nuclide, e.g. 92235.03c
+    character(10) :: name_0K             ! name of 0K nuclide, e.g. 92235.00c
+    real(8)       :: E_min = 0.01e-6_8   ! lower cutoff energy for res scattering
+    real(8)       :: E_max = 1000.0e-6_8 ! upper cutoff energy for res scattering
 
   end type Nuclide0K
 
@@ -204,7 +204,7 @@ module ace_header
 
     ! threshold for S(a,b) treatment (usually ~4 eV)
     real(8) :: threshold_inelastic
-    real(8) :: threshold_elastic = 0.0
+    real(8) :: threshold_elastic = ZERO
 
     ! Inelastic scattering data
     integer :: n_inelastic_e_in  ! # of incoming E for inelastic
@@ -258,7 +258,7 @@ module ace_header
   type NuclideMicroXS
     integer :: index_grid      ! index on nuclide energy grid
     integer :: index_temp      ! temperature index for nuclide
-    real(8) :: last_E = 0.0    ! last evaluated energy
+    real(8) :: last_E = ZERO   ! last evaluated energy
     real(8) :: interp_factor   ! interpolation factor on nuc. energy grid
     real(8) :: total           ! microscropic total xs
     real(8) :: elastic         ! microscopic elastic scattering xs
diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90
index b362dcf99a..169c763957 100644
--- a/src/cmfd_input.F90
+++ b/src/cmfd_input.F90
@@ -44,6 +44,7 @@ contains
 
   subroutine read_cmfd_xml()
 
+    use constants, only: ZERO, ONE
     use error,   only: fatal_error, warning
     use global
     use output,  only: write_message
@@ -103,7 +104,7 @@ contains
       cmfd % indices(4) = ng - 1 ! sets energy group dimension
     else
       if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
-      cmfd % egrid = (/0.0_8,20.0_8/)
+      cmfd % egrid = [ ZERO, 20.0_8 ]
       cmfd % indices(4) = 1 ! one energy group
     end if
 
@@ -111,7 +112,7 @@ contains
     if (check_for_node(node_mesh, "albedo")) then
       call get_node_array(node_mesh, "albedo", cmfd % albedo)
     else
-      cmfd % albedo = (/1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
+      cmfd % albedo = [ ONE, ONE, ONE, ONE, ONE, ONE ]
     end if
 
     ! Get acceleration map
diff --git a/src/constants.F90 b/src/constants.F90
index 0c9d15f6a4..b1af171ef1 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -27,7 +27,7 @@ module constants
   ! adjusted. Modifying constants in other sections may cause the code to fail.
 
   ! Monoatomic ideal-gas scattering treatment threshold
-  real(8), parameter :: FREE_GAS_THRESHOLD = 400.0
+  real(8), parameter :: FREE_GAS_THRESHOLD = 400.0_8
 
   ! Significance level for confidence intervals
   real(8), parameter :: CONFIDENCE_LEVEL = 0.95_8
@@ -63,15 +63,18 @@ module constants
 
   real(8), parameter ::            &
        PI           = 3.1415926535898_8, & ! pi
-       MASS_NEUTRON = 1.008664916,       & ! mass of a neutron in amu
-       MASS_PROTON  = 1.007276466812,    & ! mass of a proton in amu
-       AMU          = 1.660538921e-27,   & ! 1 amu in kg
-       N_AVOGADRO   = 0.602214129,       & ! Avogadro's number in 10^24/mol
-       K_BOLTZMANN  = 8.6173324e-11,     & ! Boltzmann constant in MeV/K
+       MASS_NEUTRON = 1.008664916_8,     & ! mass of a neutron in amu
+       MASS_PROTON  = 1.007276466812_8,  & ! mass of a proton in amu
+       AMU          = 1.660538921e-27_8, & ! 1 amu in kg
+       N_AVOGADRO   = 0.602214129_8,     & ! Avogadro's number in 10^24/mol
+       K_BOLTZMANN  = 8.6173324e-11_8,   & ! Boltzmann constant in MeV/K
        INFINITY     = huge(0.0_8),       & ! positive infinity
        ZERO         = 0.0_8,             &
+       HALF         = 0.5_8,             &
        ONE          = 1.0_8,             &
-       TWO          = 2.0_8
+       TWO          = 2.0_8,             &
+       THREE        = 3.0_8,             &
+       FOUR         = 4.0_8
 
   ! ============================================================================
   ! GEOMETRY-RELATED CONSTANTS
@@ -317,13 +320,13 @@ module constants
        OUT_FRONT = 4,   &
        IN_TOP    = 5,   &
        OUT_TOP   = 6
-  
+
   ! Tally trigger types and threshold
   integer, parameter :: &
        VARIANCE           = 1, &
        RELATIVE_ERROR     = 2, &
-       STANDARD_DEVIATION = 3 
-       
+       STANDARD_DEVIATION = 3
+
   ! Global tallY parameters
   integer, parameter :: N_GLOBAL_TALLIES = 4
   integer, parameter :: &
diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90
index 3b7f149a87..a51c6fc698 100644
--- a/src/eigenvalue.F90
+++ b/src/eigenvalue.F90
@@ -96,7 +96,7 @@ contains
       end do GENERATION_LOOP
 
       call finalize_batch()
-      
+
       if (satisfy_triggers) exit BATCH_LOOP
 
     end do BATCH_LOOP
@@ -220,12 +220,12 @@ contains
 
     ! Calculate combined estimate of k-effective
     if (master) call calculate_combined_keff()
-    
+
     ! Check_triggers
     if (master) call check_triggers()
 #ifdef MPI
     call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
-         MPI_COMM_WORLD, mpi_err)              
+         MPI_COMM_WORLD, mpi_err)
 #endif
     if (satisfy_triggers .or. &
          (trigger_on .and. current_batch == n_max_batches)) then
@@ -544,7 +544,7 @@ contains
         ! If the user did not specify how many mesh cells are to be used in
         ! each direction, we automatically determine an appropriate number of
         ! cells
-        n = ceiling((n_particles/20)**(1.0/3.0))
+        n = ceiling((n_particles/20)**(ONE/THREE))
 
         ! copy dimensions
         m % n_dimension = 3
diff --git a/src/geometry.F90 b/src/geometry.F90
index ba429566ea..d88265c537 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -1141,8 +1141,8 @@ contains
           z = coord % xyz(3)
 
           ! determine oncoming edge
-          x0 = sign(lat % pitch(1) * 0.5_8, u)
-          y0 = sign(lat % pitch(2) * 0.5_8, v)
+          x0 = sign(lat % pitch(1) * HALF, u)
+          y0 = sign(lat % pitch(2) * HALF, v)
 
           ! left and right sides
           if (abs(x - x0) < FP_PRECISION) then
@@ -1179,7 +1179,7 @@ contains
           end if
 
           if (lat % is_3d) then
-            z0 = sign(lat % pitch(3) * 0.5_8, w)
+            z0 = sign(lat % pitch(3) * HALF, w)
 
             ! top and bottom sides
             if (abs(z - z0) < FP_PRECISION) then
@@ -1212,8 +1212,8 @@ contains
           end do
 
           ! Compute velocities along the hexagonal axes.
-          beta_dir = u*sqrt(3.0_8)/2.0_8 + v/2.0_8
-          gama_dir = u*sqrt(3.0_8)/2.0_8 - v/2.0_8
+          beta_dir = u*sqrt(THREE)/TWO + v/TWO
+          gama_dir = u*sqrt(THREE)/TWO - v/TWO
 
           ! Note that hexagonal lattice distance calculations are performed
           ! using the particle's coordinates relative to the neighbor lattice
@@ -1223,13 +1223,13 @@ contains
           ! of hex lattices.
 
           ! Upper right and lower left sides.
-          edge = -sign(lat % pitch(1)/2.0_8, beta_dir)  ! Oncoming edge
-          if (beta_dir > 0.0) then
+          edge = -sign(lat % pitch(1)/TWO, beta_dir)  ! Oncoming edge
+          if (beta_dir > ZERO) then
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, 0, 0])
           else
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 0, 0])
           end if
-          beta = xyz_t(1)*sqrt(3.0_8)/2.0_8 + xyz_t(2)/2.0_8
+          beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO
           if (abs(beta - edge) < FP_PRECISION) then
             d = INFINITY
           else if (beta_dir == ZERO) then
@@ -1246,13 +1246,13 @@ contains
           end if
 
           ! Lower right and upper left sides.
-          edge = -sign(lat % pitch(1)/2.0_8, gama_dir)  ! Oncoming edge
-          if (gama_dir > 0.0) then
+          edge = -sign(lat % pitch(1)/TWO, gama_dir)  ! Oncoming edge
+          if (gama_dir > ZERO) then
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, -1, 0])
           else
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 1, 0])
           end if
-          gama = xyz_t(1)*sqrt(3.0_8)/2.0_8 - xyz_t(2)/2.0_8
+          gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO
           if (abs(gama - edge) < FP_PRECISION) then
             d = INFINITY
           else if (gama_dir == ZERO) then
@@ -1271,8 +1271,8 @@ contains
           end if
 
           ! Upper and lower sides.
-          edge = -sign(lat % pitch(1)/2.0_8, v)  ! Oncoming edge
-          if (v > 0.0) then
+          edge = -sign(lat % pitch(1)/TWO, v)  ! Oncoming edge
+          if (v > ZERO) then
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, 1, 0])
           else
             xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, -1, 0])
@@ -1296,7 +1296,7 @@ contains
 
           ! Top and bottom sides.
           if (lat % is_3d) then
-            z0 = sign(lat % pitch(2) * 0.5_8, w)
+            z0 = sign(lat % pitch(2) * HALF, w)
 
             if (abs(z - z0) < FP_PRECISION) then
               d = INFINITY
@@ -1317,10 +1317,10 @@ contains
           end if
         end select LAT_TYPE
 
-        if (d_lat < 0.0) then
+        if (d_lat < ZERO) then
           call handle_lost_particle(p, "Particle " // trim(to_str(p % id)) &
-               &//" had a negative distance to a lattice boundary. d = " &
-               &//trim(to_str(d_lat)))
+               //" had a negative distance to a lattice boundary. d = " &
+               //trim(to_str(d_lat)))
         end if
       end if LAT_COORD
 
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 9e1f1d06ce..9a5778c173 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -1,5 +1,7 @@
 module geometry_header
 
+  use constants, only: HALF, TWO, THREE
+
   implicit none
 
 !===============================================================================
@@ -29,7 +31,7 @@ module geometry_header
     integer              :: outer            ! universe to tile outside the lat
     logical              :: is_3d            ! Lattice has cells on z axis
     integer, allocatable :: offset(:,:,:,:)  ! Distribcell offsets
-    
+
     contains
 
     procedure(are_valid_indices_), deferred :: are_valid_indices
@@ -121,7 +123,7 @@ module geometry_header
      character(len=52) :: name = ""   ! User-defined name
      integer :: type                   ! Type of surface
      real(8), allocatable :: coeffs(:) ! Definition of surface
-     integer, allocatable :: & 
+     integer, allocatable :: &
           neighbor_pos(:), &           ! List of cells on positive side
           neighbor_neg(:)              ! List of cells on negative side
      integer :: bc                     ! Boundary condition
@@ -141,7 +143,7 @@ module geometry_header
      integer :: material                    ! Material within cell (0 for universe)
      integer :: n_surfaces                  ! Number of surfaces within
      integer, allocatable :: offset (:)     ! Distribcell offset for tally counter
-     integer, allocatable :: & 
+     integer, allocatable :: &
           & surfaces(:)                     ! List of surfaces bounding cell -- note that
                                             ! parentheses, union, etc operators will be listed
                                             ! here too
@@ -189,7 +191,7 @@ contains
     real(8),            intent(in) :: global_xyz(3)
     integer                        :: i_xyz(3)
 
-    real(8) :: xyz(3)  ! global_xyz alias 
+    real(8) :: xyz(3)  ! global_xyz alias
 
     xyz = global_xyz
 
@@ -209,7 +211,7 @@ contains
     real(8),           intent(in) :: global_xyz(3)
     integer                       :: i_xyz(3)
 
-    real(8) :: xyz(3)    ! global_xyz alias 
+    real(8) :: xyz(3)    ! global_xyz alias
     real(8) :: alpha     ! Skewed coord axis
     real(8) :: xyz_t(3)  ! Local xyz
     real(8) :: dists(4)  ! Squared distances from cell centers
@@ -220,16 +222,16 @@ contains
 
     ! Index z direction.
     if (this % is_3d) then
-      i_xyz(3) = ceiling((xyz(3) - this % center(3))/this % pitch(2) + 0.5_8)&
-           &+ this % n_axial/2
+      i_xyz(3) = ceiling((xyz(3) - this % center(3))/this % pitch(2) + HALF)&
+           + this % n_axial/2
     else
       i_xyz(3) = 1
     end if
 
     ! Convert coordinates into skewed bases.  The (x, alpha) basis is used to
     ! find the index of the global coordinates to within 4 cells.
-    alpha = xyz(2) - xyz(1) / sqrt(3.0_8)
-    i_xyz(1) = floor(xyz(1) / (sqrt(3.0_8) / 2.0_8 * this % pitch(1)))
+    alpha = xyz(2) - xyz(1) / sqrt(THREE)
+    i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1)))
     i_xyz(2) = floor(alpha / this % pitch(1))
 
     ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
@@ -279,12 +281,12 @@ contains
     xyz = global_xyz
 
     local_xyz(1) = xyz(1) - (this % lower_left(1) + &
-         &(i_xyz(1) - 0.5_8)*this % pitch(1))
+         (i_xyz(1) - HALF)*this % pitch(1))
     local_xyz(2) = xyz(2) - (this % lower_left(2) + &
-         &(i_xyz(2) - 0.5_8)*this % pitch(2))
+         (i_xyz(2) - HALF)*this % pitch(2))
     if (this % is_3d) then
       local_xyz(3) = xyz(3) - (this % lower_left(3) + &
-           &(i_xyz(3) - 0.5_8)*this % pitch(3))
+           (i_xyz(3) - HALF)*this % pitch(3))
     else
       local_xyz(3) = xyz(3)
     end if
@@ -304,14 +306,14 @@ contains
 
     ! x_l = x_g - (center + pitch_x*cos(30)*index_x)
     local_xyz(1) = xyz(1) - (this % center(1) + &
-         &sqrt(3.0_8) / 2.0_8 * (i_xyz(1) - this % n_rings) * this % pitch(1))
+         sqrt(THREE) / TWO * (i_xyz(1) - this % n_rings) * this % pitch(1))
     ! y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
     local_xyz(2) = xyz(2) - (this % center(2) + &
-         &(i_xyz(2) - this % n_rings) * this % pitch(1) + &
-         &(i_xyz(1) - this % n_rings) * this % pitch(1) / 2.0_8)
+         (i_xyz(2) - this % n_rings) * this % pitch(1) + &
+         (i_xyz(1) - this % n_rings) * this % pitch(1) / TWO)
     if (this % is_3d) then
       local_xyz(3) = xyz(3) - this % center(3) &
-           &+ (this % n_axial/2 - i_xyz(3) + 1) * this % pitch(2)
+           + (this % n_axial/2 - i_xyz(3) + 1) * this % pitch(2)
     else
       local_xyz(3) = xyz(3)
     end if
diff --git a/src/global.F90 b/src/global.F90
index 2ec111dc08..982743c544 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -243,8 +243,8 @@ module global
   ! VARIANCE REDUCTION VARIABLES
 
   logical :: survival_biasing = .false.
-  real(8) :: weight_cutoff = 0.25
-  real(8) :: weight_survive = 1.0
+  real(8) :: weight_cutoff = 0.25_8
+  real(8) :: weight_survive = ONE
 
   ! ============================================================================
   ! HDF5 VARIABLES
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 53e4567bf7..9305313511 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -1796,11 +1796,11 @@ contains
         case ('g/cc', 'g/cm3')
           mat % density = -val
         case ('kg/m3')
-          mat % density = -0.001 * val
+          mat % density = -0.001_8 * val
         case ('atom/b-cm')
           mat % density = val
         case ('atom/cm3', 'atom/cc')
-          mat % density = 1.0e-24 * val
+          mat % density = 1.0e-24_8 * val
         case default
           call fatal_error("Unkwown units '" // trim(units) &
                &// "' specified on material " // trim(to_str(mat % id)))
diff --git a/src/math.F90 b/src/math.F90
index b3046ccc4a..58406bda18 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -1,6 +1,6 @@
 module math
 
-  use constants,  only: PI, ONE, TWO, ZERO
+  use constants
   use random_lcg, only: prn
 
   implicit none
@@ -41,27 +41,27 @@ contains
 
       q = sqrt(-TWO*log(p))
       z = (((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
-           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
+           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE)
 
-    elseif (p <= 1. - p_low) then
+    elseif (p <= ONE - p_low) then
       ! Rational approximation for central region
 
-      q = p - 0.5
+      q = p - HALF
       r = q*q
       z = (((((a(1)*r + a(2))*r + a(3))*r + a(4))*r + a(5))*r + a(6))*q / &
-           (((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + 1.)
+           (((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + ONE)
 
     else
       ! Rational approximation for upper region
 
-      q = sqrt(-2*log(1. - p))
+      q = sqrt(-TWO*log(ONE - p))
       z = -(((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
-           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + 1.)
+           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE)
     endif
 
     ! Refinement based on Newton's method
 #ifndef NO_F2008
-    z = z - (0.5 * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(0.5*z*z)
+    z = z - (HALF * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(HALF*z*z)
 #endif
 
   end function normal_percentile
@@ -86,13 +86,13 @@ contains
       ! For one degree of freedom, the t-distribution becomes a Cauchy
       ! distribution whose cdf we can invert directly
 
-      t = tan(PI*(p - 0.5))
+      t = tan(PI*(p - HALF))
 
     elseif (df == 2) then
       ! For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
       ! 2)). This can be directly inverted to yield the solution below
 
-      t = TWO*sqrt(TWO)*(p - 0.5)/sqrt(ONE - 4.*(p - 0.5)**2)
+      t = TWO*sqrt(TWO)*(p - HALF)/sqrt(ONE - FOUR*(p - HALF)**2)
 
     else
 
@@ -102,12 +102,12 @@ contains
       ! 16 (4), pp. 1123-1132 (1987).
 
       n = real(df,8)
-      k = 1./(n - 2.)
+      k = ONE/(n - TWO)
       z = normal_percentile(p)
       z2 = z * z
-      t = sqrt(n*k) * (z + (z2 - 3.)*z*k/4. + ((5.*z2 - 56.)*z2 + &
-           75.)*z*k*k/96. + (((z2 - 27.)*3.*z2 + 417.)*z2 - 315.) &
-           *z*k*k*k/384.)
+      t = sqrt(n*k) * (z + (z2 - THREE)*z*k/FOUR + ((5._8*z2 - 56._8)*z2 + &
+           75._8)*z*k*k/96._8 + (((z2 - 27._8)*THREE*z2 + 417._8)*z2 - 315._8) &
+           *z*k*k*k/384._8)
 
     end if
 
@@ -134,7 +134,7 @@ contains
     case(1)
       pnx = x
     case(2)
-      pnx = 1.5_8 * x * x - 0.5_8
+      pnx = 1.5_8 * x * x - HALF
     case(3)
       pnx = 2.5_8 * x * x * x - 1.5_8 * x
     case(4)
@@ -199,37 +199,37 @@ contains
       rn(3) = ONE*sqrt(w2m1) * cos(phi)
     case (2)
       ! l = 2, m = -2
-      rn(1) = 0.288675134594813_8 * (-3.0_8 * w**2 + 3.0_8) * sin(TWO*phi)
+      rn(1) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * sin(TWO*phi)
       ! l = 2, m = -1
       rn(2) = 1.73205080756888_8 * w*sqrt(w2m1) * sin(phi)
       ! l = 2, m = 0
-      rn(3) = 1.5_8 * w**2 - 0.5_8
+      rn(3) = 1.5_8 * w**2 - HALF
       ! l = 2, m = 1
       rn(4) = 1.73205080756888_8 * w*sqrt(w2m1) * cos(phi)
       ! l = 2, m = 2
-      rn(5) = 0.288675134594813_8 * (-3.0_8 * w**2 + 3.0_8) * cos(TWO*phi)
+      rn(5) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * cos(TWO*phi)
     case (3)
       ! l = 3, m = -3
-      rn(1) = 0.790569415042095_8 * (w2m1)**(3.0_8/TWO) * sin(3.0_8 * phi)
+      rn(1) = 0.790569415042095_8 * (w2m1)**(THREE/TWO) * sin(THREE * phi)
       ! l = 3, m = -2
       rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi)
       ! l = 3, m = -1
-      rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - 3.0_8/TWO) * &
+      rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
         sin(phi)
       ! l = 3, m = 0
       rn(4) = 2.5_8 * w**3 - 1.5_8 * w
       ! l = 3, m = 1
-      rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - 3.0_8/TWO) * &
+      rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
         cos(phi)
       ! l = 3, m = 2
       rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi)
       ! l = 3, m = 3
-      rn(7) = 0.790569415042095_8 * (w2m1)**(3.0_8/TWO) * cos(3.0_8* phi)
+      rn(7) = 0.790569415042095_8 * (w2m1)**(THREE/TWO) * cos(THREE* phi)
     case (4)
       ! l = 4, m = -4
       rn(1) = 0.739509972887452_8 * (w2m1)**2 * sin(4.0_8*phi)
       ! l = 4, m = -3
-      rn(2) = 2.09165006633519_8 * w*(w2m1)**(3.0_8/TWO) * sin(3.0_8* phi)
+      rn(2) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi)
       ! l = 4, m = -2
       rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
         sin(TWO*phi)
@@ -245,7 +245,7 @@ contains
       rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
         cos(TWO*phi)
       ! l = 4, m = 3
-      rn(8) = 2.09165006633519_8 * w*(w2m1)**(3.0_8/TWO) * cos(3.0_8* phi)
+      rn(8) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi)
       ! l = 4, m = 4
       rn(9) = 0.739509972887452_8 * (w2m1)**2 * cos(4.0_8*phi)
     case (5)
@@ -254,8 +254,8 @@ contains
       ! l = 5, m = -4
       rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi)
       ! l = 5, m = -3
-      rn(3) = 0.00996023841111995_8 * (w2m1)**(3.0_8/TWO)* &
-        ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(3.0_8*phi)
+      rn(3) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
+        ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)
       ! l = 5, m = -2
       rn(4) = 0.0487950036474267_8 * (w2m1)*((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * &
         sin(TWO*phi)
@@ -271,8 +271,8 @@ contains
       rn(8) = 0.0487950036474267_8 * (w2m1)* &
         ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
       ! l = 5, m = 3
-      rn(9) = 0.00996023841111995_8 * (w2m1)**(3.0_8/TWO)* &
-        ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(3.0_8*phi)
+      rn(9) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
+        ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)
       ! l = 5, m = 4
       rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi)
       ! l = 5, m = 5
@@ -286,8 +286,8 @@ contains
       rn(3) = 0.00104990131391452_8 * (w2m1)**2 * &
         ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
       ! l = 6, m = -3
-      rn(4) = 0.00575054632785295_8 * (w2m1)**(3.0_8/TWO) * &
-        ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(3.0_8*phi)
+      rn(4) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
+        ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)
       ! l = 6, m = -2
       rn(5) = 0.0345032779671177_8 * (w2m1) * &
         ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * sin(TWO*phi)
@@ -303,8 +303,8 @@ contains
       rn(9) = 0.0345032779671177_8 * (w2m1) * &
         ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * cos(TWO*phi)
       ! l = 6, m = 3
-      rn(10) = 0.00575054632785295_8 * (w2m1)**(3.0_8/TWO) * &
-        ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(3.0_8*phi)
+      rn(10) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
+        ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)
       ! l = 6, m = 4
       rn(11) = 0.00104990131391452_8 * (w2m1)**2 * &
         ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
@@ -324,9 +324,9 @@ contains
       rn(4) = 0.000548293079133141_8 * (w2m1)**2* &
         ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
       ! l = 7, m = -3
-      rn(5) = 0.00363696483726654_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(5) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
         ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
-        sin(3.0_8*phi)
+        sin(THREE*phi)
       ! l = 7, m = -2
       rn(6) = 0.025717224993682_8 * (w2m1)* &
         ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
@@ -346,9 +346,9 @@ contains
         ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
         cos(TWO*phi)
       ! l = 7, m = 3
-      rn(11) = 0.00363696483726654_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(11) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
       ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
-      cos(3.0_8*phi)
+      cos(THREE*phi)
       ! l = 7, m = 4
       rn(12) = 0.000548293079133141_8 * (w2m1)**2 * &
         ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
@@ -374,8 +374,8 @@ contains
       rn(5) = 0.000316557156832328_8 * (w2m1)**2* &
         ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * sin(4.0_8*phi)
       ! l = 8, m = -3
-      rn(6) = 0.00245204119306875_8 * (w2m1)**(3.0_8/TWO)* &
-        ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(3.0_8*phi)
+      rn(6) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
+        ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(THREE*phi)
       ! l = 8, m = -2
       rn(7) = 0.0199204768222399_8 * (w2m1)* &
         ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
@@ -396,9 +396,9 @@ contains
         45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
         315.0_8/16.0_8) * cos(TWO*phi)
       ! l = 8, m = 3
-      rn(12) = 0.00245204119306875_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(12) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
         ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
-        (10395.0_8/8.0_8)*w) * cos(3.0_8*phi)
+        (10395.0_8/8.0_8)*w) * cos(THREE*phi)
       ! l = 8, m = 4
       rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - &
         135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
@@ -431,9 +431,9 @@ contains
       rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
         675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
       ! l = 9, m = -3
-      rn(7) = 0.00173385495536766_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(7) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)* &
         ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
-        (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(3.0_8*phi)
+        (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)
       ! l = 9, m = -2
       rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- &
         135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/16.0_8 * w)*  &
@@ -452,9 +452,9 @@ contains
         135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/ 16.0_8 * w) * &
         cos(TWO*phi)
       ! l = 9, m = 3
-      rn(13) = 0.00173385495536766_8 * (w2m1)**(3.0_8/TWO)*((765765.0_8/16.0_8)*w**6 - &
+      rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)*w**6 - &
         675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8)* &
-        cos(3.0_8*phi)
+        cos(THREE*phi)
       ! l = 9, m = 4
       rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
         675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
@@ -494,9 +494,9 @@ contains
         11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
         45045.0_8/16.0_8) * sin(4.0_8*phi)
       ! l = 10, m = -3
-      rn(8) = 0.00127230170115096_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(8) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
         ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
-        (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(3.0_8*phi)
+        (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)
       ! l = 10, m = -2
       rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
         765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
@@ -517,9 +517,9 @@ contains
         765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
         45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
       ! l = 10, m = 3
-      rn(14) = 0.00127230170115096_8 * (w2m1)**(3.0_8/TWO)* &
+      rn(14) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
         ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
-        (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(3.0_8*phi)
+        (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)
       ! l = 10, m = 4
       rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
         11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
diff --git a/src/physics.F90 b/src/physics.F90
index 4e94ea2261..04be42fd99 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -597,13 +597,13 @@ contains
           if (r > ONE) then
             ! equally likely N-4 middle bins
             j = int(r) + 2
-          elseif (r > 0.6) then
+          elseif (r > 0.6_8) then
             ! second to last bin has relative probability of 0.4
             j = n_energy_out - 1
-          elseif (r > 0.5) then
+          elseif (r > HALF) then
             ! last bin has relative probability of 0.1
             j = n_energy_out
-          elseif (r > 0.1) then
+          elseif (r > 0.1_8) then
             ! second bin has relative probability of 0.4
             j = 2
           else
@@ -818,8 +818,8 @@ contains
 
     case ('dbrc')
       E_red = sqrt((awr * E) / kT)
-      E_low = (((E_red - 4.0_8)**2) * kT) / awr
-      E_up  = (((E_red + 4.0_8)**2) * kT) / awr
+      E_low = (((E_red - FOUR)**2) * kT) / awr
+      E_up  = (((E_red + FOUR)**2) * kT) / awr
 
       ! find lower and upper energy bound indices
       ! lower index
@@ -872,8 +872,8 @@ contains
 
     case ('ares')
       E_red = sqrt((awr * E) / kT)
-      E_low = (((E_red - 4.0_8)**2) * kT) / awr
-      E_up  = (((E_red + 4.0_8)**2) * kT) / awr
+      E_low = (((E_red - FOUR)**2) * kT) / awr
+      E_up  = (((E_red + FOUR)**2) * kT) / awr
 
       ! find lower and upper energy bound indices
       ! lower index
diff --git a/src/plot.F90 b/src/plot.F90
index d9f6b55e2b..01f4a60213 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -46,7 +46,7 @@ contains
   end subroutine run_plot
 
 !===============================================================================
-! POSITION_RGB computes the red/green/blue values for a given plot with the 
+! POSITION_RGB computes the red/green/blue values for a given plot with the
 ! current particle's position
 !===============================================================================
 
@@ -56,12 +56,12 @@ contains
     type(ObjectPlot), pointer, intent(in) :: pl
     integer, intent(out)                  :: rgb(3)
     integer, intent(out)                  :: id
-    
+
     logical :: found_cell
     integer :: level
     type(Cell),       pointer :: c => null()
     type(LocalCoord), pointer :: coord => null()
-    
+
     call deallocate_coord(p % coord0 % next)
     p % coord => p % coord0
 
@@ -77,7 +77,7 @@ contains
       coord => coord % next
       level = level + 1
     end do
-    
+
     if (.not. found_cell) then
       ! If no cell, revert to default color
       rgb = pl % not_found % rgb
@@ -107,7 +107,7 @@ contains
         id = -1
       end if
     end if
-    
+
   end subroutine position_rgb
 
 !===============================================================================
@@ -141,31 +141,31 @@ contains
     if (pl % basis == PLOT_BASIS_XY) then
       in_i  = 1
       out_i = 2
-      xyz(1) = pl % origin(1) - pl % width(1) / 2.0
-      xyz(2) = pl % origin(2) + pl % width(2) / 2.0
+      xyz(1) = pl % origin(1) - pl % width(1) / TWO
+      xyz(2) = pl % origin(2) + pl % width(2) / TWO
       xyz(3) = pl % origin(3)
     else if (pl % basis == PLOT_BASIS_XZ) then
       in_i  = 1
       out_i = 3
-      xyz(1) = pl % origin(1) - pl % width(1) / 2.0
+      xyz(1) = pl % origin(1) - pl % width(1) / TWO
       xyz(2) = pl % origin(2)
-      xyz(3) = pl % origin(3) + pl % width(2) / 2.0
+      xyz(3) = pl % origin(3) + pl % width(2) / TWO
     else if (pl % basis == PLOT_BASIS_YZ) then
       in_i  = 2
       out_i = 3
       xyz(1) = pl % origin(1)
-      xyz(2) = pl % origin(2) - pl % width(1) / 2.0
-      xyz(3) = pl % origin(3) + pl % width(2) / 2.0
+      xyz(2) = pl % origin(2) - pl % width(1) / TWO
+      xyz(3) = pl % origin(3) + pl % width(2) / TWO
     end if
 
     ! allocate and initialize particle
     call p % initialize()
     p % coord % xyz = xyz
-    p % coord % uvw = (/ 0.5, 0.5, 0.5 /)
+    p % coord % uvw = [ HALF, HALF, HALF ]
     p % coord % universe = BASE_UNIVERSE
 
     do y = 1, img % height
-      call progress % set_value(dble(y)/dble(img % height)*100.)
+      call progress % set_value(dble(y)/dble(img % height)*100)
       do x = 1, img % width
 
         ! get pixel color
@@ -201,10 +201,10 @@ contains
 ! DRAW_MESH_LINES draws mesh line boundaries on an image
 !===============================================================================
   subroutine draw_mesh_lines(pl, img)
-  
+
     type(ObjectPlot), pointer, intent(in)    :: pl
     type(Image),               intent(inout) :: img
-    
+
     logical :: in_mesh
     integer :: out_, in_  ! pixel location
     integer :: r, g, b    ! RGB color for meshlines pixels
@@ -221,13 +221,13 @@ contains
     real(8) :: xyz_ll(3)  ! lower left xyz
     real(8) :: xyz_ur(3)  ! upper right xyz
     type(StructuredMesh), pointer :: m => null()
-  
+
     m => pl % meshlines_mesh
-    
+
     r = pl % meshlines_color % rgb(1)
     g = pl % meshlines_color % rgb(2)
     b = pl % meshlines_color % rgb(3)
-    
+
     select case (pl % basis)
       case(PLOT_BASIS_XY)
         outer = 1
@@ -243,10 +243,10 @@ contains
     xyz_ll_plot = pl % origin
     xyz_ur_plot = pl % origin
 
-    xyz_ll_plot(outer) = pl % origin(1) - pl % width(1) / 2.0
-    xyz_ll_plot(inner) = pl % origin(2) - pl % width(2) / 2.0
-    xyz_ur_plot(outer) = pl % origin(1) + pl % width(1) / 2.0
-    xyz_ur_plot(inner) = pl % origin(2) + pl % width(2) / 2.0
+    xyz_ll_plot(outer) = pl % origin(1) - pl % width(1) / TWO
+    xyz_ll_plot(inner) = pl % origin(2) - pl % width(2) / TWO
+    xyz_ur_plot(outer) = pl % origin(1) + pl % width(1) / TWO
+    xyz_ur_plot(inner) = pl % origin(2) + pl % width(2) / TWO
 
     width = xyz_ur_plot - xyz_ll_plot
 
@@ -259,24 +259,24 @@ contains
         ! check if we're in the mesh for this ijk
         if (i > 0 .and. i <= m % dimension(outer) .and. &
             j > 0 .and. j <= m % dimension(inner)) then
-          
+
           ! get xyz's of lower left and upper right of this mesh cell
           xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
           xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
           xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
           xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
-          
+
           ! map the xyz ranges to pixel ranges
-          
+
           frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
           outrange(1) = int(frac * real(img % width, 8))
           frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
           outrange(2) = int(frac * real(img % width, 8))
 
           frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(1) = int((1. - frac) * real(img % height, 8))
+          inrange(1) = int((ONE - frac) * real(img % height, 8))
           frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(2) = int((1. - frac) * real(img % height, 8))
+          inrange(2) = int((ONE - frac) * real(img % height, 8))
 
           ! draw lines
           do out_ = outrange(1), outrange(2)
@@ -295,11 +295,11 @@ contains
               call set_pixel(img, outrange(2) - plus, in_, r, g, b)
             end do
           end do
-          
+
         end if
       end do
     end do
-    
+
   end subroutine draw_mesh_lines
 
 !===============================================================================
@@ -350,7 +350,7 @@ contains
   subroutine create_3d_dump(pl)
 
     type(ObjectPlot), pointer :: pl
-    
+
     integer :: x, y, z      ! voxel location indices
     integer :: rgb(3)       ! colors (red, green, blue) from 0-255
     integer :: id           ! id of cell or material
@@ -361,14 +361,14 @@ contains
 
     ! compute voxel widths in each direction
     vox = pl % width/dble(pl % pixels)
-    
+
     ! initial particle position
-    ll = pl % origin - pl % width / 2.0
+    ll = pl % origin - pl % width / TWO
 
     ! allocate and initialize particle
     call p % initialize()
     p % coord0 % xyz = ll
-    p % coord0 % uvw = (/ 0.5, 0.5, 0.5 /)
+    p % coord0 % uvw = [ HALF, HALF, HALF ]
     p % coord0 % universe = BASE_UNIVERSE
 
     ! Open binary plot file for writing
@@ -378,11 +378,11 @@ contains
     ! write plot header info
     write(UNIT_PLOT) pl % pixels, vox, ll
 
-    ! move to center of voxels    
-    ll = ll + vox / 2.0
+    ! move to center of voxels
+    ll = ll + vox / TWO
 
     do x = 1, pl % pixels(1)
-      call progress % set_value(dble(x)/dble(pl % pixels(1))*100.)
+      call progress % set_value(dble(x)/dble(pl % pixels(1))*100)
       do y = 1, pl % pixels(2)
         do z = 1, pl % pixels(3)
 
@@ -394,20 +394,20 @@ contains
 
           ! advance particle in z direction
           p % coord0 % xyz(3) = p % coord0 % xyz(3) + vox(3)
-          
+
         end do
-        
+
         ! advance particle in y direction
         p % coord0 % xyz(2) = p % coord0 % xyz(2) + vox(2)
         p % coord0 % xyz(3) = ll(3)
-        
+
       end do
-      
+
       ! advance particle in y direction
       p % coord0 % xyz(1) = p % coord0 % xyz(1) + vox(1)
       p % coord0 % xyz(2) = ll(2)
       p % coord0 % xyz(3) = ll(3)
-      
+
     end do
 
     close(UNIT_PLOT)
diff --git a/src/trigger.F90 b/src/trigger.F90
index d5efe96d54..d485e057d6 100644
--- a/src/trigger.F90
+++ b/src/trigger.F90
@@ -4,6 +4,7 @@ module trigger
   use mpi
 #endif
 
+  use constants
   use global
   use string,           only: to_str
   use output,           only: warning, write_message
@@ -30,7 +31,7 @@ contains
     character(len=52)  :: name            ! "eigenvalue" or tally score
 
     integer    :: n_pred_batches  ! predicted # batches to satisfy all triggers
-    
+
     ! Checks if current_batch is one for which the triggers must be checked
     if (current_batch < n_batches .or. (.not. trigger_on)) return
     if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. &
@@ -39,24 +40,24 @@ contains
     ! Check the trigger and output the result
     call check_tally_triggers(max_ratio, tally_id, name)
 
-    ! When trigger threshold is reached, write information 
+    ! When trigger threshold is reached, write information
     if (satisfy_triggers) then
       call write_message("Triggers satisfied for batch " // &
           trim(to_str(current_batch)))
-    
+
     ! When trigger is not reached write convergence info for user
     elseif (name == "eigenvalue") then
       call write_message("Triggers unsatisfied, max unc./thresh. is " // &
            trim(to_str(max_ratio)) //  " for " // trim(name))
     else
-      call write_message("Triggers unsatisfied, max unc./thresh. is " // & 
+      call write_message("Triggers unsatisfied, max unc./thresh. is " // &
            trim(to_str(max_ratio)) // " for " // trim(name) // &
            " in tally " // trim(to_str(tally_id)))
-    end if 
+    end if
 
     ! If batch_interval is not set, estimate batches till triggers are satisfied
     if (pred_batches .and. .not. satisfy_triggers) then
-         
+
       ! Estimate the number of remaining batches to convergence
       ! The prediction uses the fact that tally variances are proportional
       ! to 1/N where N is the number of the batches/particles
@@ -65,9 +66,9 @@ contains
       n_pred_batches = n_batch_interval + n_batches
 
       ! Write the predicted number of batches for the user
-      if (n_pred_batches > n_max_batches) then 
+      if (n_pred_batches > n_max_batches) then
         call warning("The estimated number of batches is " // &
-             trim(to_str(n_pred_batches)) // & 
+             trim(to_str(n_pred_batches)) // &
              " --  greater than max batches. ")
       else
         call write_message("The estimated number of batches is " // &
@@ -98,8 +99,8 @@ contains
     integer :: n_order        ! loop index for moment orders
     integer :: nm_order       ! loop index for Ynm moment orders
     real(8) :: uncertainty    ! trigger uncertainty
-    real(8) :: std_dev = 0.0  ! trigger standard deviation
-    real(8) :: rel_err = 0.0  ! trigger relative error
+    real(8) :: std_dev = ZERO ! trigger standard deviation
+    real(8) :: rel_err = ZERO ! trigger relative error
     real(8) :: ratio          ! ratio of the uncertainty/trigger threshold
     type(TallyObject), pointer     :: t               ! tally pointer
     type(TriggerObject), pointer   :: trigger         ! tally trigger
@@ -115,8 +116,8 @@ contains
       ! Check eigenvalue trigger
       if (run_mode == MODE_EIGENVALUE) then
         if (keff_trigger % trigger_type /= 0) then
-          select case (keff_trigger % trigger_type)        
-          case(VARIANCE) 
+          select case (keff_trigger % trigger_type)
+          case(VARIANCE)
             uncertainty = k_combined(2) ** 2
           case(STANDARD_DEVIATION)
             uncertainty = k_combined(2)
@@ -124,7 +125,7 @@ contains
             uncertainty = k_combined(2) / k_combined(1)
           end select
 
-          ! If uncertainty is above threshold, store uncertainty ratio         
+          ! If uncertainty is above threshold, store uncertainty ratio
           if (uncertainty > keff_trigger % threshold) then
             satisfy_triggers = .false.
             if (keff_trigger % trigger_type == VARIANCE) then
@@ -133,11 +134,11 @@ contains
               ratio = uncertainty / keff_trigger % threshold
             end if
             if (max_ratio < ratio) then
-              max_ratio = ratio             
+              max_ratio = ratio
               name = "eigenvalue"
-            end if 
-          end if 
-        end if   
+            end if
+          end if
+        end if
       end if
 
       ! Compute uncertainties for all tallies, scores with triggers
@@ -153,10 +154,10 @@ contains
           trigger => t % triggers(s)
 
           ! Initialize trigger uncertainties to zero
-          trigger % std_dev = 0.
-          trigger % rel_err = 0.
-          trigger % variance = 0.
-        
+          trigger % std_dev = ZERO
+          trigger % rel_err = ZERO
+          trigger % variance = ZERO
+
           ! Surface current tally triggers require special treatment
           if (t % type == TALLY_SURFACE_CURRENT) then
             call compute_tally_current(t, trigger)
@@ -178,9 +179,9 @@ contains
                   j = j - 1
                 else
                   if (j == t % n_filters) exit find_bin
-                end if 
+                end if
               end do find_bin
-        
+
               if (t % n_filters > 0) then
                 filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * &
                      t % stride) + 1
@@ -193,13 +194,13 @@ contains
 
               ! Initialize score bin index
               NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins
-                
+
                 select case(t % score_bins(trigger % score_index))
 
                 case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
 
                   score_index = score_index - 1
-             
+
                   do n_order = 0, t % moment_order(trigger % score_index)
                     score_index = score_index + 1
 
@@ -207,17 +208,17 @@ contains
                          score_index, filter_index, t)
 
                     if (trigger % variance < variance) then
-                      trigger % variance = std_dev ** 2 
+                      trigger % variance = std_dev ** 2
                     end if
-                    if (trigger % std_dev < std_dev) then 
+                    if (trigger % std_dev < std_dev) then
                       trigger % std_dev = std_dev
                     end if
-                    if (trigger % rel_err < rel_err) then 
+                    if (trigger % rel_err < rel_err) then
                       trigger % rel_err = rel_err
                     end if
 
                   end do
-              
+
                 case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
                      SCORE_TOTAL_YN)
 
@@ -242,32 +243,32 @@ contains
 
                     end do
                   end do
-           
+
                 case default
                   call get_trigger_uncertainty(std_dev, rel_err, &
                        score_index, filter_index, t)
 
                   if (trigger % variance < variance) then
-                    trigger % variance = std_dev ** 2 
+                    trigger % variance = std_dev ** 2
                   end if
-                  if (trigger % std_dev < std_dev) then 
+                  if (trigger % std_dev < std_dev) then
                     trigger % std_dev = std_dev
                   end if
-                  if (trigger % rel_err < rel_err) then 
+                  if (trigger % rel_err < rel_err) then
                     trigger % rel_err = rel_err
                   end if
 
                 end select
- 
-                select case (t % triggers(s) % type)        
-                case(VARIANCE) 
+
+                select case (t % triggers(s) % type)
+                case(VARIANCE)
                   uncertainty = trigger % variance
                 case(STANDARD_DEVIATION)
                   uncertainty = trigger % std_dev
                 case default
                   uncertainty = trigger % rel_err
                 end select
-              
+
                 if (uncertainty > t % triggers(s) % threshold) then
                   satisfy_triggers = .false.
 
@@ -275,8 +276,8 @@ contains
                     ratio = sqrt(uncertainty / t % triggers(s) % threshold)
                   else
                     ratio = uncertainty / t % triggers(s) % threshold
-                  end if 
-                  
+                  end if
+
                   if (max_ratio < ratio) then
                     max_ratio = ratio
                     name  = t % triggers(s) % score_name
@@ -297,7 +298,7 @@ contains
 ! COMPUTE_TALLY_CURRENT computes the current for a surface current tally with
 ! precision trigger(s).
 !===============================================================================
- 
+
  subroutine compute_tally_current(t, trigger)
 
     integer :: i                    ! mesh index for x
@@ -310,8 +311,8 @@ contains
     integer :: n                    ! number of incoming energy bins
     integer :: filter_index         ! index in results array for filters
     logical :: print_ebin           ! should incoming energy bin be displayed?
-    real(8) :: rel_err  = 0.0         ! temporary relative error of result
-    real(8) :: std_dev  = 0.0         ! temporary standard deviration of result
+    real(8) :: rel_err  = ZERO      ! temporary relative error of result
+    real(8) :: std_dev  = ZERO      ! temporary standard deviration of result
     type(TallyObject), pointer    :: t        ! surface current tally
     type(TriggerObject)           :: trigger  ! surface current tally trigger
     type(StructuredMesh), pointer :: m        ! surface current mesh
@@ -350,10 +351,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = std_dev**2
@@ -362,10 +363,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -377,22 +378,22 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
-            
+
             matching_bins(i_filter_surf) = OUT_RIGHT
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -404,10 +405,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -417,10 +418,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -432,10 +433,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -444,10 +445,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -459,10 +460,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -471,10 +472,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -486,10 +487,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -498,10 +499,10 @@ contains
             filter_index = &
                  sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
             call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-            if (trigger % std_dev < std_dev) then 
+            if (trigger % std_dev < std_dev) then
               trigger % std_dev = std_dev
             end if
-            if (trigger % rel_err < rel_err) then 
+            if (trigger % rel_err < rel_err) then
               trigger % rel_err = rel_err
             end if
             trigger % variance = trigger % std_dev**2
@@ -539,11 +540,11 @@ contains
     std_dev = sqrt((tally_result % sum_sq / n - mean * mean) / (n - 1))
 
     ! Compute the relative error if the mean is non-zero
-    if (mean == 0.) then
-      rel_err = 0.
+    if (mean == ZERO) then
+      rel_err = ZERO
     else
       rel_err = std_dev / mean
-    end if  
+    end if
 
   end subroutine get_trigger_uncertainty
 
diff --git a/tests/test_cmfd_feed/results_true.dat b/tests/test_cmfd_feed/results_true.dat
index 95dec91dd4..a1dd0d815e 100644
--- a/tests/test_cmfd_feed/results_true.dat
+++ b/tests/test_cmfd_feed/results_true.dat
@@ -14,7 +14,7 @@ tally 1:
 3.894180E+01
 1.517824E+02
 3.528006E+01
-1.246309E+02
+1.246308E+02
 2.863448E+01
 8.222321E+01
 2.125384E+01
diff --git a/tests/test_cmfd_nofeed/results_true.dat b/tests/test_cmfd_nofeed/results_true.dat
index 761a4749da..da9b802e7e 100644
--- a/tests/test_cmfd_nofeed/results_true.dat
+++ b/tests/test_cmfd_nofeed/results_true.dat
@@ -1,5 +1,5 @@
 k-combined:
-1.170519E+00 8.422959E-03
+1.170519E+00 8.422960E-03
 tally 1:
 1.078122E+01
 1.170828E+01
@@ -47,7 +47,7 @@ tally 2:
 6.213554E+00
 1.946061E+00
 7.376629E+01
-2.729167E+02
+2.729166E+02
 5.253400E+01
 1.385018E+02
 6.438590E+00
diff --git a/tests/test_density_atomcm3/results_true.dat b/tests/test_density_atomcm3/results_true.dat
index 718d2bc9f7..3feb35ba1e 100644
--- a/tests/test_density_atomcm3/results_true.dat
+++ b/tests/test_density_atomcm3/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-1.092203E+00 1.990175E-02
+1.092203E+00 1.990176E-02
diff --git a/tests/test_density_kgm3/results_true.dat b/tests/test_density_kgm3/results_true.dat
index 4fbd4e4eb9..c5942e1a23 100644
--- a/tests/test_density_kgm3/results_true.dat
+++ b/tests/test_density_kgm3/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-8.085745E-01 9.674582E-03
+8.085745E-01 9.674599E-03
diff --git a/tests/test_filter_distribcell/case-3/results_true.dat b/tests/test_filter_distribcell/case-3/results_true.dat
index 1378f9cf75..f5f85d29f5 100644
--- a/tests/test_filter_distribcell/case-3/results_true.dat
+++ b/tests/test_filter_distribcell/case-3/results_true.dat
@@ -1 +1 @@
-f5873c6ca03a2c75f5094b7297059285788bb93c450466606914805cdb24fd2664c3872ed42a8c6665a272070ab6b08b0a53b02612bcd993e10329aec2f1d313
\ No newline at end of file
+6008cf2ba8eecaaa5a600fa337cf54cef018e98bdba8e3bd26c6f44587376a838d5bc5e86301b2e308f9eb248e3efafd45a5336f4023d962d7921d158a621e0c
\ No newline at end of file
diff --git a/tests/test_filter_mesh_2d/results_true.dat b/tests/test_filter_mesh_2d/results_true.dat
index 126a6151fb..01aea5afe2 100644
--- a/tests/test_filter_mesh_2d/results_true.dat
+++ b/tests/test_filter_mesh_2d/results_true.dat
@@ -264,7 +264,7 @@ tally 1:
 5.320317E-01
 2.267827E-01
 9.971891E-01
-5.070903E-01
+5.070904E-01
 7.088786E-02
 5.025089E-03
 3.601108E-02
@@ -298,7 +298,7 @@ tally 1:
 1.145124E+00
 4.571704E-01
 1.135664E+00
-4.319757E-01
+4.319756E-01
 8.875810E-02
 7.878000E-03
 7.494206E-02
diff --git a/tests/test_filter_mesh_3d/results_true.dat b/tests/test_filter_mesh_3d/results_true.dat
index e43a5c4d95..a67c77c034 100644
--- a/tests/test_filter_mesh_3d/results_true.dat
+++ b/tests/test_filter_mesh_3d/results_true.dat
@@ -1687,7 +1687,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.206688E-02
+4.206687E-02
 1.769622E-03
 2.760907E-02
 7.622607E-04
@@ -1810,7 +1810,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 8.032141E-03
-6.451528E-05
+6.451529E-05
 1.897802E-01
 3.601652E-02
 0.000000E+00
@@ -2232,7 +2232,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 1.103939E-01
-1.218682E-02
+1.218681E-02
 1.828427E-01
 1.294521E-02
 0.000000E+00
@@ -2577,8 +2577,8 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.905315E-02
-6.249401E-03
+7.905316E-02
+6.249402E-03
 6.679360E-01
 1.650894E-01
 7.539474E-02
@@ -2616,7 +2616,7 @@ tally 1:
 1.107212E+00
 3.413570E-01
 7.306165E-02
-5.338005E-03
+5.338004E-03
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -2840,9 +2840,9 @@ tally 1:
 0.000000E+00
 0.000000E+00
 1.935296E-02
-3.745372E-04
-4.914513E-02
-2.415244E-03
+3.745370E-04
+4.914514E-02
+2.415245E-03
 3.076980E-01
 3.868177E-02
 2.022942E-01
@@ -2875,10 +2875,10 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.698117E-02
-5.106927E-03
+7.698116E-02
+5.106926E-03
 6.397704E-01
-2.487692E-01
+2.487693E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -3124,7 +3124,7 @@ tally 1:
 8.227568E-02
 6.769287E-03
 2.170015E-02
-4.708967E-04
+4.708966E-04
 1.342835E-01
 1.803205E-02
 0.000000E+00
@@ -3139,8 +3139,8 @@ tally 1:
 0.000000E+00
 1.408666E-01
 1.984341E-02
-2.912092E-02
-8.480282E-04
+2.912093E-02
+8.480284E-04
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -3193,8 +3193,8 @@ tally 1:
 0.000000E+00
 2.698309E-01
 5.659690E-02
-7.539544E-03
-5.684472E-05
+7.539545E-03
+5.684475E-05
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -3357,7 +3357,7 @@ tally 1:
 2.499121E-01
 5.858178E-01
 1.482815E-01
-4.753060E-02
+4.753061E-02
 2.259158E-03
 0.000000E+00
 0.000000E+00
@@ -3594,7 +3594,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 1.275913E-02
-8.900080E-05
+8.900081E-05
 7.974141E-02
 3.206666E-03
 0.000000E+00
@@ -3701,10 +3701,10 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.224856E-01
+2.224855E-01
 3.374440E-02
-9.939442E-02
-9.879251E-03
+9.939443E-02
+9.879253E-03
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -3769,10 +3769,10 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.394141E-02
-5.731913E-04
+2.394140E-02
+5.731908E-04
 2.225881E-01
-4.954544E-02
+4.954545E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -3797,8 +3797,8 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.232219E-03
-2.737612E-05
+5.232218E-03
+2.737611E-05
 2.699331E-01
 7.286389E-02
 0.000000E+00
@@ -3887,8 +3887,8 @@ tally 1:
 3.750508E-02
 4.494863E-01
 2.020380E-01
-1.524812E-02
-2.325052E-04
+1.524811E-02
+2.325049E-04
 0.000000E+00
 0.000000E+00
 5.591708E-01
@@ -4281,7 +4281,7 @@ tally 1:
 0.000000E+00
 3.933218E-02
 1.547021E-03
-4.824740E-02
+4.824739E-02
 2.327811E-03
 0.000000E+00
 0.000000E+00
@@ -4465,8 +4465,8 @@ tally 1:
 6.344565E-02
 4.457955E-01
 9.943168E-02
-4.346781E-03
-1.889451E-05
+4.346780E-03
+1.889450E-05
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -4585,7 +4585,7 @@ tally 1:
 1.026482E-01
 5.456703E-01
 1.139540E-01
-5.551968E-01
+5.551969E-01
 1.452833E-01
 2.630569E-01
 2.654789E-02
@@ -5166,7 +5166,7 @@ tally 1:
 4.295071E-01
 9.355140E-02
 1.311425E-03
-1.719837E-06
+1.719836E-06
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -5612,7 +5612,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 2.913120E-03
-8.486267E-06
+8.486266E-06
 9.098084E-02
 8.277514E-03
 0.000000E+00
@@ -6108,9 +6108,9 @@ tally 1:
 1.475983E-02
 2.178526E-04
 1.747711E+00
-8.535010E-01
+8.535011E-01
 7.461108E-02
-5.566813E-03
+5.566814E-03
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -6453,7 +6453,7 @@ tally 1:
 0.000000E+00
 1.194798E-01
 1.427541E-02
-2.732160E-03
+2.732161E-03
 7.464701E-06
 0.000000E+00
 0.000000E+00
@@ -6604,11 +6604,11 @@ tally 1:
 0.000000E+00
 0.000000E+00
 2.362808E-01
-5.091358E-02
-7.208809E-03
-5.196693E-05
+5.091357E-02
+7.208810E-03
+5.196695E-05
 5.757630E-01
-2.687940E-01
+2.687941E-01
 9.778317E-02
 4.820693E-03
 0.000000E+00
@@ -7349,8 +7349,8 @@ tally 1:
 1.915958E-04
 2.060666E-01
 4.246344E-02
-9.954400E-04
-9.909008E-07
+9.954401E-04
+9.909010E-07
 7.940105E-02
 6.304526E-03
 0.000000E+00
diff --git a/tests/test_lattice/results_true.dat b/tests/test_lattice/results_true.dat
index 4e02678c63..cf51dd5d73 100644
--- a/tests/test_lattice/results_true.dat
+++ b/tests/test_lattice/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-1.042387E+00 1.575316E-01
+1.042388E+00 1.575316E-01
diff --git a/tests/test_lattice_hex/results_true.dat b/tests/test_lattice_hex/results_true.dat
index e4e53357d2..b88285ff2e 100644
--- a/tests/test_lattice_hex/results_true.dat
+++ b/tests/test_lattice_hex/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-2.831017E-01 2.269866E-02
+2.831014E-01 2.269849E-02
diff --git a/tests/test_lattice_mixed/results_true.dat b/tests/test_lattice_mixed/results_true.dat
index f3c94aec22..d23feb91f0 100644
--- a/tests/test_lattice_mixed/results_true.dat
+++ b/tests/test_lattice_mixed/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-9.670751E-01 1.903257E-02
+9.605085E-01 8.180822E-03
diff --git a/tests/test_ptables_off/results_true.dat b/tests/test_ptables_off/results_true.dat
index 4566281596..17d6e91453 100644
--- a/tests/test_ptables_off/results_true.dat
+++ b/tests/test_ptables_off/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-2.984064E-01 3.464414E-03
+2.984064E-01 3.464413E-03
diff --git a/tests/test_reflective_cylinder/results_true.dat b/tests/test_reflective_cylinder/results_true.dat
index 8b58143354..f98c517889 100644
--- a/tests/test_reflective_cylinder/results_true.dat
+++ b/tests/test_reflective_cylinder/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-2.276055E+00 3.186525E-03
+2.276055E+00 3.186526E-03
diff --git a/tests/test_reflective_sphere/results_true.dat b/tests/test_reflective_sphere/results_true.dat
index 46a6c278c9..47aa5c763e 100644
--- a/tests/test_reflective_sphere/results_true.dat
+++ b/tests/test_reflective_sphere/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-2.276020E+00 2.831786E-03
+2.280812E+00 5.087968E-03
diff --git a/tests/test_score_fission/results_true.dat b/tests/test_score_fission/results_true.dat
index a79228e577..1a3d333838 100644
--- a/tests/test_score_fission/results_true.dat
+++ b/tests/test_score_fission/results_true.dat
@@ -8,7 +8,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 9.094323E-01
-1.773089E-01
+1.773090E-01
 tally 2:
 1.081920E+00
 2.404883E-01
diff --git a/tests/test_score_flux_yn/results_true.dat b/tests/test_score_flux_yn/results_true.dat
index 9a2d5f3230..8c04c7c898 100644
--- a/tests/test_score_flux_yn/results_true.dat
+++ b/tests/test_score_flux_yn/results_true.dat
@@ -46,7 +46,7 @@ tally 2:
 5.069422E-02
 -2.045060E-01
 2.197889E-01
-1.307603E-01
+1.307602E-01
 5.938649E-02
 -3.317004E-01
 1.359233E-01
@@ -54,7 +54,7 @@ tally 2:
 7.528253E-02
 7.685982E-02
 4.559715E-02
--4.668783E-02
+-4.668784E-02
 4.065614E-02
 1.447872E-01
 5.538956E-02
@@ -67,7 +67,7 @@ tally 2:
 -3.028859E-01
 3.382353E-02
 -2.320613E-01
-2.845905E-02
+2.845904E-02
 8.038034E-02
 3.466746E-02
 1.766484E-01
@@ -90,17 +90,17 @@ tally 2:
 2.704386E+01
 -1.955516E-01
 9.789273E-03
-4.070150E-02
+4.070149E-02
 3.530892E-02
 -1.446334E-01
 3.374966E-02
--4.072015E-02
+-4.072014E-02
 1.407320E-02
 1.056190E-01
 2.144988E-02
 -1.491728E-01
 1.699068E-02
-3.130539E-02
+3.130540E-02
 4.945783E-03
 -1.706770E-01
 6.943395E-03
@@ -121,17 +121,17 @@ tally 2:
 9.919028E-02
 1.278789E-02
 -1.274511E-01
-9.357942E-03
+9.357943E-03
 -9.037073E-02
 4.686665E-03
 -7.895199E-03
 8.526809E-03
--2.071828E-02
-5.844892E-03
+-2.071827E-02
+5.844891E-03
 -1.422877E-02
 2.008170E-03
 -7.787500E-02
-9.642363E-03
+9.642362E-03
 -6.105504E-02
 7.600775E-03
 -1.249117E-01
@@ -141,20 +141,20 @@ tally 2:
 -9.171222E-02
 4.220465E-03
 8.798390E-02
-5.037482E-03
+5.037483E-03
 4.366496E-02
 8.411743E-03
 1.605678E-02
-6.994512E-03
+6.994511E-03
 -2.179875E-02
 1.352098E-03
 1.116755E-01
-7.354423E-03
+7.354424E-03
 9.995739E-02
-9.710025E-03
+9.710026E-03
 3.010324E-02
 1.179938E-02
-2.175816E-02
+2.175815E-02
 1.281441E-03
 2.374886E-02
 6.943432E-03
@@ -188,7 +188,7 @@ tally 2:
 4.811523E-01
 -4.932440E-01
 1.789106E-01
--3.705747E-02
+-3.705746E-02
 5.517840E-01
 7.049764E-01
 4.243943E-01
@@ -196,7 +196,7 @@ tally 2:
 2.014189E-01
 -5.420278E-01
 1.201572E-01
-7.606098E-02
+7.606099E-02
 1.451071E-01
 2.244026E-01
 7.628680E-02
@@ -206,7 +206,7 @@ tally 2:
 1.958724E-01
 -3.799191E-01
 2.258678E-01
--5.909625E-01
+-5.909624E-01
 2.940322E-01
 -6.336238E-01
 1.272552E-01
@@ -232,7 +232,7 @@ tally 2:
 2.482914E-01
 2.924945E+01
 1.782032E+02
-5.986907E-01
+5.986908E-01
 4.163713E-01
 3.376452E-01
 1.744277E-01
@@ -258,7 +258,7 @@ tally 2:
 5.960686E-02
 -9.968550E-02
 4.039270E-02
--1.235292E-02
+-1.235291E-02
 2.805714E-01
 -7.494860E-02
 2.838999E-02
@@ -281,7 +281,7 @@ tally 2:
 7.284620E-01
 1.871262E-01
 1.360134E-01
-9.437336E-03
+9.437335E-03
 -3.088199E-01
 2.292832E-02
 -2.543226E-01
@@ -306,39 +306,39 @@ tally 2:
 2.067319E+01
 2.381967E-01
 4.786935E-02
--1.924307E-02
+-1.924308E-02
 4.029344E-02
 3.179175E-02
 2.098033E-02
--1.652307E-02
+-1.652308E-02
 9.391525E-03
--5.733037E-02
-1.788394E-02
+-5.733039E-02
+1.788395E-02
 4.942155E-02
 1.206623E-03
-1.874859E-02
+1.874861E-02
 1.917059E-02
 -1.179838E-01
 2.185564E-02
--9.270518E-03
+-9.270508E-03
 1.659002E-02
-7.496510E-02
+7.496511E-02
 6.849931E-03
 -1.004441E-01
 4.570271E-03
--9.899603E-02
+-9.899601E-02
 1.027435E-02
--9.963859E-02
+-9.963860E-02
 4.686519E-03
 -1.740638E-02
 3.936514E-02
-2.520348E-02
+2.520349E-02
 5.900520E-03
 1.227791E-03
 6.173783E-03
 -3.068441E-02
-6.867403E-03
--7.765653E-02
+6.867402E-03
+-7.765654E-02
 1.399383E-02
 7.333551E-02
 4.139444E-03
@@ -348,29 +348,29 @@ tally 2:
 4.918888E-03
 -1.106116E-01
 1.143785E-02
-7.254492E-03
-1.337966E-03
+7.254477E-03
+1.337965E-03
 2.365948E-01
 2.473485E-02
 5.176646E-02
 1.514750E-03
--8.723609E-02
-2.395099E-03
+-8.723610E-02
+2.395100E-03
 -1.010207E-01
 1.081070E-02
 -5.113298E-02
 6.053279E-03
 -1.601115E-01
 8.304230E-03
--5.105346E-03
+-5.105345E-03
 5.079654E-03
 -8.894549E-02
 2.505026E-03
 -1.023874E-01
 4.311776E-03
-7.972835E-02
+7.972837E-02
 1.178305E-02
--2.560943E-02
+-2.560942E-02
 4.101793E-03
 2.312203E-02
 3.036209E-03
@@ -382,7 +382,7 @@ tally 2:
 1.400732E+00
 4.934620E-01
 7.673322E-01
--2.964444E-02
+-2.964442E-02
 2.262500E-01
 -2.225956E-01
 4.872938E-01
@@ -418,15 +418,15 @@ tally 2:
 2.484316E-01
 -2.258513E-01
 8.839687E-02
--9.307566E-01
+-9.307565E-01
 4.802462E-01
 1.108577E-01
-4.451362E-02
+4.451363E-02
 1.118965E+00
 5.173979E-01
 4.314600E-01
 5.176751E-02
--2.864060E-01
+-2.864059E-01
 9.159895E-02
 -4.471918E-01
 1.743046E-01
diff --git a/tests/test_score_total_yn/results_true.dat b/tests/test_score_total_yn/results_true.dat
index 36d7040bd9..b2dcd39e12 100644
--- a/tests/test_score_total_yn/results_true.dat
+++ b/tests/test_score_total_yn/results_true.dat
@@ -116,7 +116,7 @@ tally 2:
 5.215322E-04
 1.571095E-02
 9.166057E-04
-7.226989E-03
+7.226990E-03
 3.683499E-04
 -2.456165E-02
 1.571441E-04
@@ -154,7 +154,7 @@ tally 2:
 6.914145E-05
 -5.873760E-03
 2.229738E-04
-6.666510E-03
+6.666511E-03
 1.394147E-04
 -3.245528E-03
 1.550876E-04
@@ -192,7 +192,7 @@ tally 2:
 6.688120E-03
 -8.231311E-02
 2.743044E-02
-1.172791E-02
+1.172790E-02
 1.089457E-02
 -1.545389E-01
 2.655167E-02
@@ -201,7 +201,7 @@ tally 2:
 2.773254E-02
 1.070478E-02
 -2.522703E-02
-1.136506E-02
+1.136507E-02
 -1.287564E-02
 4.969443E-03
 3.471629E-02
@@ -277,34 +277,34 @@ tally 2:
 5.659008E-03
 4.488733E-04
 -1.625829E-02
-2.325229E-04
+2.325228E-04
 2.554220E-02
 7.666018E-04
 -2.251672E-02
 7.703442E-04
--8.302473E-03
+-8.302474E-03
 2.748216E-04
 -7.692439E-02
 2.167893E-03
-7.267773E-03
-7.781908E-04
+7.267774E-03
+7.781907E-04
 6.821338E-03
-9.518921E-04
+9.518920E-04
 -2.203951E-02
 2.214610E-03
 2.052514E-02
 3.756195E-04
--4.158574E-02
-9.330261E-04
+-4.158575E-02
+9.330262E-04
 -3.000932E-02
 4.584914E-04
 -2.247072E-02
 2.519845E-04
 1.589401E-02
 2.302726E-04
--7.793790E-04
+-7.793787E-04
 1.588955E-04
--4.634908E-03
+-4.634907E-03
 5.395803E-04
 -1.193817E-02
 1.674419E-04
@@ -368,7 +368,7 @@ tally 2:
 4.167797E-01
 -4.347270E-01
 3.702158E-01
--7.344612E-02
+-7.344613E-02
 5.940699E-02
 1.081217E+00
 4.070812E-01
@@ -377,7 +377,7 @@ tally 2:
 2.656115E-01
 8.561866E-02
 -1.286933E-01
-8.347413E-02
+8.347414E-02
 1.258195E-01
 3.329577E-02
 -3.751107E-01
@@ -406,7 +406,7 @@ tally 2:
 5.017106E-02
 2.266133E-01
 1.131899E-01
-7.161674E-04
+7.161662E-04
 3.345992E-02
 -1.987041E-01
 1.381432E-01
diff --git a/tests/test_source_angle_mono/results_true.dat b/tests/test_source_angle_mono/results_true.dat
index e8ca8b9ea3..1aca9a5057 100644
--- a/tests/test_source_angle_mono/results_true.dat
+++ b/tests/test_source_angle_mono/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-3.034005E-01 5.499078E-03
+3.034005E-01 5.499077E-03
diff --git a/tests/test_sourcepoint_restart/results_true.dat b/tests/test_sourcepoint_restart/results_true.dat
index 0b4587d206..917fc0cdfb 100644
--- a/tests/test_sourcepoint_restart/results_true.dat
+++ b/tests/test_sourcepoint_restart/results_true.dat
@@ -610,7 +610,7 @@ tally 1:
 1.264629E-03
 5.949953E-07
 4.409634E-03
-5.032697E-06
+5.032698E-06
 0.000000E+00
 0.000000E+00
 0.000000E+00
diff --git a/tests/test_statepoint_batch/results_true.dat b/tests/test_statepoint_batch/results_true.dat
index 45b045ecf5..786fd55461 100644
--- a/tests/test_statepoint_batch/results_true.dat
+++ b/tests/test_statepoint_batch/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-3.109090E-01 5.807934E-03
+3.109090E-01 5.807935E-03
diff --git a/tests/test_statepoint_restart/results_true.dat b/tests/test_statepoint_restart/results_true.dat
index ea320003a9..378d07639f 100644
--- a/tests/test_statepoint_restart/results_true.dat
+++ b/tests/test_statepoint_restart/results_true.dat
@@ -729,7 +729,7 @@ tally 1:
 1.477216E-05
 2.698876E-03
 3.931578E-06
-7.470783E-03
+7.470784E-03
 2.887854E-05
 0.000000E+00
 0.000000E+00
@@ -1210,7 +1210,7 @@ tally 1:
 3.278961E-03
 7.369638E-06
 7.462474E-03
-2.983181E-05
+2.983182E-05
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -1690,7 +1690,7 @@ tally 1:
 9.790439E-04
 2.045088E-06
 5.377302E-03
-1.508922E-05
+1.508923E-05
 0.000000E+00
 0.000000E+00
 0.000000E+00
diff --git a/tests/test_trigger_batch_interval/results_true.dat b/tests/test_trigger_batch_interval/results_true.dat
index 46e3d115ae..65fa51bb4f 100644
--- a/tests/test_trigger_batch_interval/results_true.dat
+++ b/tests/test_trigger_batch_interval/results_true.dat
@@ -1,5 +1,5 @@
 k-combined:
-9.851940E-01 4.283951E-03
+9.851940E-01 4.283950E-03
 tally 1:
 1.972289E+01
 2.779778E+01
diff --git a/tests/test_trigger_no_batch_interval/results_true.dat b/tests/test_trigger_no_batch_interval/results_true.dat
index 46e3d115ae..65fa51bb4f 100644
--- a/tests/test_trigger_no_batch_interval/results_true.dat
+++ b/tests/test_trigger_no_batch_interval/results_true.dat
@@ -1,5 +1,5 @@
 k-combined:
-9.851940E-01 4.283951E-03
+9.851940E-01 4.283950E-03
 tally 1:
 1.972289E+01
 2.779778E+01
diff --git a/tests/test_trigger_tallies/results_true.dat b/tests/test_trigger_tallies/results_true.dat
index 418856e065..86b1c14d5c 100644
--- a/tests/test_trigger_tallies/results_true.dat
+++ b/tests/test_trigger_tallies/results_true.dat
@@ -1,5 +1,5 @@
 k-combined:
-9.824135E-01 6.844128E-03
+9.824135E-01 6.844127E-03
 tally 1:
 1.408810E+01
 1.985836E+01