Initial structure and allocation of surface current tallies.

src/constants.f90

@@ -153,7 +153,7 @@ module constants
        N_DA        = 117
 
   ! Tally macro reactions
-  integer, parameter :: N_MACRO_TYPES = 14
+  integer, parameter :: N_MACRO_TYPES = 15
   integer, parameter :: &
        MACRO_FLUX       = -1,  & ! flux
        MACRO_TOTAL      = -2,  & ! total reaction rate
@@ -168,7 +168,8 @@ module constants
        MACRO_N_4N       = -11, & ! (n,4n) rate
        MACRO_ABSORPTION = -12, & ! absorption rate
        MACRO_FISSION    = -13, & ! fission rate
-       MACRO_NU_FISSION = -14    ! neutron production rate
+       MACRO_NU_FISSION = -14, & ! neutron production rate
+       MACRO_CURRENT    = -15    ! partial current
 
   ! Tally map bin finding
   integer, parameter :: NO_BIN_FOUND = -1

src/input_xml.f90

@@ -790,6 +790,46 @@ contains
                 end if
              case ('nu-fission')
                 t % macro_bins(j) % scalar = MACRO_NU_FISSION
+             case ('current')
+                t % macro_bins(j) % scalar = MACRO_CURRENT
+                t % surface_current = .true.
+
+                ! Check to make sure that current is the only desired response
+                ! for this tally
+                if (n_words > 1) then
+                   msg = "Cannot tally other macro reactions in the same " &
+                        // "tally as surface currents. Separate other macro " &
+                        // "reactions into a distinct tally."
+                   call fatal_error(msg)
+                end if
+
+                ! Check to make sure that only the mesh filter was specified
+                if (t % mesh == 0 .or. t % n_bins(T_MESH) /= & 
+                     product(t % n_bins, t % n_bins > 0)) then
+                   msg = "Surface currents must be used with a mesh filter only."
+                   call fatal_error(msg)
+                end if
+
+                ! Since the number of bins for the mesh filter was already set
+                ! assuming it was a flux tally, we need to adjust the number of
+                ! bins
+                t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
+
+                ! Get pointer to mesh
+                uid = t % mesh
+                index = dict_get_key(mesh_dict, uid)
+                m => meshes(index)
+
+                ! We need to increase the dimension by one since we also need
+                ! currents coming into and out of the boundary mesh cells.
+                if (size(m % dimension) == 2) then
+                   t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
+                        product(m % dimension + 1) * 4
+                elseif (size(m % dimension) == 3) then
+                   t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
+                        product(m % dimension + 1) * 6
+                end if
+
              case default
                 msg = "Unknown macro reaction: " // trim(words(j))
                 call fatal_error(msg)

src/tally.f90

@@ -274,6 +274,7 @@ contains
     integer :: score_index
     integer :: score_index0
     integer :: macro_bin
+    integer :: mesh_bin
     real(8) :: score
     real(8) :: last_wgt
     real(8) :: wgt
@@ -343,8 +344,17 @@ contains
        ! determine mesh bin
        if (t % n_bins(T_MESH) > 0) then
           m => meshes(t % mesh)
-          call get_mesh_bin(m, p % xyz, bins(T_MESH), in_mesh)
+
+          ! Determine if we're in the mesh first
+          call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh)
           if (.not. in_mesh) cycle
+
+          if (t % surface_current) then
+             msg = "Surface current mesh tally not yet implemented."
+             call fatal_error(msg)
+          else
+             bins(T_MESH) = mesh_bin
+          end if
        else
           bins(T_MESH) = 1
        end if

src/tally_header.f90

@@ -61,6 +61,7 @@ module tally_header
      integer :: uid
      integer :: type
      real(8) :: volume
+     logical :: surface_current = .false.
 
      ! Tally bin specifications