diff --git a/src/constants.F90 b/src/constants.F90
index 9bd5df3dcb..20edd77e5a 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -243,7 +243,7 @@ module constants
        EVENT_FISSION =  3 
 
   ! Tally score type
-  integer, parameter :: N_SCORE_TYPES = 16
+  integer, parameter :: N_SCORE_TYPES = 17
   integer, parameter :: &
        SCORE_FLUX       = -1,  & ! flux
        SCORE_TOTAL      = -2,  & ! total reaction rate
@@ -252,15 +252,16 @@ module constants
        SCORE_SCATTER_1  = -5,  & ! first scattering moment
        SCORE_SCATTER_2  = -6,  & ! second scattering moment
        SCORE_SCATTER_3  = -7,  & ! third scattering moment
-       SCORE_DIFFUSION  = -8,  & ! diffusion coefficient
-       SCORE_N_1N       = -9,  & ! (n,1n) rate
-       SCORE_N_2N       = -10, & ! (n,2n) rate
-       SCORE_N_3N       = -11, & ! (n,3n) rate
-       SCORE_N_4N       = -12, & ! (n,4n) rate
-       SCORE_ABSORPTION = -13, & ! absorption rate
-       SCORE_FISSION    = -14, & ! fission rate
-       SCORE_NU_FISSION = -15, & ! neutron production rate
-       SCORE_CURRENT    = -16    ! partial current
+       SCORE_TRANSPORT  = -8,  & ! transport reaction rate
+       SCORE_DIFFUSION  = -9,  & ! diffusion coefficient
+       SCORE_N_1N       = -10, & ! (n,1n) rate
+       SCORE_N_2N       = -11, & ! (n,2n) rate
+       SCORE_N_3N       = -12, & ! (n,3n) rate
+       SCORE_N_4N       = -13, & ! (n,4n) rate
+       SCORE_ABSORPTION = -14, & ! absorption rate
+       SCORE_FISSION    = -15, & ! fission rate
+       SCORE_NU_FISSION = -16, & ! neutron production rate
+       SCORE_CURRENT    = -17    ! partial current
 
   ! Tally map bin finding
   integer, parameter :: NO_BIN_FOUND = -1
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 48eed7784b..2c01771985 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -1093,6 +1093,11 @@ contains
              case ('scatter-3')
                 t % score_bins(j) % scalar = SCORE_SCATTER_3
 
+                ! Set tally estimator to analog
+                t % estimator = ESTIMATOR_ANALOG
+             case('transport')
+                t % score_bins(j) % scalar = SCORE_TRANSPORT
+
                 ! Set tally estimator to analog
                 t % estimator = ESTIMATOR_ANALOG
              case ('diffusion')
diff --git a/src/tally.F90 b/src/tally.F90
index bef12c571c..b90750c4d1 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -248,6 +248,8 @@ contains
     real(8) :: last_wgt             ! pre-collision particle weight
     real(8) :: wgt                  ! post-collision particle weight
     real(8) :: mu                   ! cosine of angle of collision
+    real(8) :: macro_total          ! material macro total xs
+    real(8) :: macro_scatt          ! material macro scatt xs
     logical :: found_bin            ! scoring bin found?
     type(TallyObject), pointer :: t => null()
 
@@ -346,6 +348,20 @@ contains
 
              score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu)
 
+          case (SCORE_TRANSPORT)
+            ! Skip any event where the particle didn't scatter
+            if (p % event /= EVENT_SCATTER) cycle
+
+            ! get material macros
+            macro_total = material_xs % total
+            macro_scatt = material_xs % total - material_xs % absorption
+
+            ! Score total rate - p1 scatter rate
+            ! Note estimator needs to be adjusted since tallying is only
+            ! occuring when a scatter has happend. Effectively this means
+            ! multiplying the estimator by total/scatter macro
+            score = (macro_total - mu*macro_scatt)*(ONE/macro_scatt)
+
           case (SCORE_DIFFUSION)
              ! Skip any event where the particle didn't scatter
              if (p % event /= EVENT_SCATTER) cycle
@@ -1528,6 +1544,7 @@ contains
     score_name(abs(SCORE_SCATTER_1))  = "First Scattering Moment"
     score_name(abs(SCORE_SCATTER_2))  = "Second Scattering Moment"
     score_name(abs(SCORE_SCATTER_3))  = "Third Scattering Moment"
+    score_name(abs(SCORE_TRANSPORT))  = "Transport Rate"
     score_name(abs(SCORE_DIFFUSION))  = "Diffusion Coefficient"
     score_name(abs(SCORE_N_1N))       = "(n,1n) Rate"
     score_name(abs(SCORE_N_2N))       = "(n,2n) Rate"