diff --git a/CMakeLists.txt b/CMakeLists.txt
index be62dcd57b..9517fe2ec3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -422,6 +422,7 @@ set(LIBOPENMC_FORTRAN_SRC
   src/tallies/tally_filter_energyfunc.F90
   src/tallies/tally_filter_material.F90
   src/tallies/tally_filter_mesh.F90
+  src/tallies/tally_filter_meshsurface.F90
   src/tallies/tally_filter_mu.F90
   src/tallies/tally_filter_polar.F90
   src/tallies/tally_filter_surface.F90
diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index f1633f3f9d..070b5bd201 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -109,6 +109,7 @@ Constructing Tallies
    openmc.CellbornFilter
    openmc.SurfaceFilter
    openmc.MeshFilter
+   openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter
    openmc.MuFilter
diff --git a/examples/xml/pincell/tallies.xml b/examples/xml/pincell/tallies.xml
index 0ae36eceda..7e1e0dafe4 100644
--- a/examples/xml/pincell/tallies.xml
+++ b/examples/xml/pincell/tallies.xml
@@ -8,7 +8,7 @@
   </mesh>
 
   <filter id="1" type="mesh">
-    <bins>1</bins>
+    <bins>2</bins>
   </filter>
 
   <filter id="2" type="energy">
diff --git a/include/openmc.h b/include/openmc.h
index bb04907a52..10c0b2de8e 100644
--- a/include/openmc.h
+++ b/include/openmc.h
@@ -52,6 +52,7 @@ extern "C" {
   int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n);
   int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins);
   int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
+  int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
   int openmc_next_batch();
   int openmc_nuclide_name(int index, char** name);
   void openmc_plot_geometry();
diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py
index 5a5df4814d..691ecc09cc 100644
--- a/openmc/capi/filter.py
+++ b/openmc/capi/filter.py
@@ -55,6 +55,9 @@ _dll.openmc_material_filter_set_bins.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_mesh_filter_set_mesh.restype = c_int
 _dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
+_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
 
 
 class Filter(_FortranObjectWithID):
@@ -188,6 +191,10 @@ class MeshFilter(Filter):
     filter_type = 'mesh'
 
 
+class MeshSurfaceFilter(Filter):
+    filter_type = 'meshsurface'
+
+
 class MuFilter(Filter):
     filter_type = 'mu'
 
@@ -216,6 +223,7 @@ _FILTER_TYPE_MAP = {
     'energyfunction': EnergyFunctionFilter,
     'material': MaterialFilter,
     'mesh': MeshFilter,
+    'meshsurface': MeshSurfaceFilter,
     'mu': MuFilter,
     'polar': PolarFilter,
     'surface': SurfaceFilter,
diff --git a/openmc/filter.py b/openmc/filter.py
index 26629dafe0..c1b9f8dc92 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -15,6 +15,7 @@ import openmc.checkvalue as cv
 from .cell import Cell
 from .material import Material
 from .mixin import IDManagerMixin
+from .surface import Surface
 from .universe import Universe
 
 
@@ -22,12 +23,11 @@ _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
                  'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal',
                  'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom']
 
-_CURRENT_NAMES = {1:  'x-min out', 2:  'x-min in',
-                  3:  'x-max out', 4:  'x-max in',
-                  5:  'y-min out', 6:  'y-min in',
-                  7:  'y-max out', 8:  'y-max in',
-                  9:  'z-min out', 10: 'z-min in',
-                  11: 'z-max out', 12: 'z-max in'}
+_CURRENT_NAMES = (
+    'x-min out', 'x-min in', 'x-max out', 'x-max in',
+    'y-min out', 'y-min in', 'y-max out', 'y-max in',
+    'z-min out', 'z-min in', 'z-max out', 'z-max in'
+)
 
 
 class FilterMeta(ABCMeta):
@@ -497,9 +497,9 @@ class CellFilter(WithIDFilter):
 
     Parameters
     ----------
-    bins : openmc.Cell, Integral, or iterable thereof
-        The Cells to tally. Either openmc.Cell objects or their
-        Integral ID numbers can be used.
+    bins : openmc.Cell, int, or iterable thereof
+        The cells to tally. Either openmc.Cell objects or their ID numbers can
+        be used.
     filter_id : int
         Unique identifier for the filter
 
@@ -564,86 +564,29 @@ class CellbornFilter(WithIDFilter):
     expected_type = Cell
 
 
-class SurfaceFilter(Filter):
-    """Bins particle currents on Mesh surfaces.
+class SurfaceFilter(WithIDFilter):
+    """Filters particles by surface crossing
 
     Parameters
     ----------
-    bins : Iterable of Integral
-        Indices corresponding to which face of a mesh cell the current is
-        crossing.
+    bins : openmc.Surface, int, or iterable of Integral
+        The surfaces to tally over. Either openmc.Surface objects or their ID
+        numbers can be used.
     filter_id : int
         Unique identifier for the filter
 
     Attributes
     ----------
     bins : Iterable of Integral
-        Indices corresponding to which face of a mesh cell the current is
-        crossing.
+        The surfaces to tally over. Either openmc.Surface objects or their ID
+        numbers can be used.
     id : int
         Unique identifier for the filter
     num_bins : Integral
         The number of filter bins
 
     """
-    @Filter.bins.setter
-    def bins(self, bins):
-        # Format the bins as a 1D numpy array.
-        bins = np.atleast_1d(bins)
-
-        # Check the bin values.
-        cv.check_iterable_type('filter bins', bins, Integral)
-        for edge in bins:
-            cv.check_greater_than('filter bin', edge, 0, equality=True)
-
-        self._bins = bins
-
-    @property
-    def num_bins(self):
-        # Need to handle number of bins carefully -- for surface current
-        # tallies, the number of bins depends on the mesh, which we don't have a
-        # reference to in this filter
-        return self._num_bins
-
-    def get_pandas_dataframe(self, data_size, stride, **kwargs):
-        """Builds a Pandas DataFrame for the Filter's bins.
-
-        This method constructs a Pandas DataFrame object for the filter with
-        columns annotated by filter bin information. This is a helper method for
-        :meth:`Tally.get_pandas_dataframe`.
-
-        Parameters
-        ----------
-        data_size : int
-            The total number of bins in the tally corresponding to this filter
-        stride : int
-            Stride in memory for the filter
-
-        Returns
-        -------
-        pandas.DataFrame
-            A Pandas DataFrame with a column of strings describing which surface
-            the current is crossing and which direction it points.  The number
-            of rows in the DataFrame is the same as the total number of bins in
-            the corresponding tally, with the filter bin appropriately tiled to
-            map to the corresponding tally bins.
-
-        See also
-        --------
-        Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
-
-        """
-        # Initialize Pandas DataFrame
-        df = pd.DataFrame()
-
-        filter_bins = np.repeat(self.bins, stride)
-        tile_factor = data_size // len(filter_bins)
-        filter_bins = np.tile(filter_bins, tile_factor)
-        filter_bins = [_CURRENT_NAMES[x] for x in filter_bins]
-        df = pd.concat([df, pd.DataFrame(
-            {self.short_name.lower(): filter_bins})])
-
-        return df
+    expected_type = Surface
 
 
 class MeshFilter(Filter):
@@ -689,7 +632,6 @@ class MeshFilter(Filter):
         filter_id = int(group.name.split('/')[-1].lstrip('filter '))
 
         out = cls(mesh_obj, filter_id=filter_id)
-        out._num_bins = group['n_bins'].value
 
         return out
 
@@ -705,10 +647,7 @@ class MeshFilter(Filter):
 
     @property
     def num_bins(self):
-        try:
-            return self._num_bins
-        except AttributeError:
-            return reduce(operator.mul, self.mesh.dimension)
+        return reduce(operator.mul, self.mesh.dimension)
 
     def check_bins(self, bins):
         if not len(bins) == 1:
@@ -802,7 +741,7 @@ class MeshFilter(Filter):
         filter_dict = {}
 
         # Append Mesh ID as outermost index of multi-index
-        mesh_key = 'mesh {0}'.format(self.mesh.id)
+        mesh_key = 'mesh {}'.format(self.mesh.id)
 
         # Find mesh dimensions - use 3D indices for simplicity
         n_dim = len(self.mesh.dimension)
@@ -845,6 +784,137 @@ class MeshFilter(Filter):
         return df
 
 
+class MeshSurfaceFilter(MeshFilter):
+    """Filter events by surface crossings on a regular, rectangular mesh.
+
+    Parameters
+    ----------
+    mesh : openmc.Mesh
+        The Mesh object that events will be tallied onto
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : Integral
+        The Mesh ID
+    mesh : openmc.Mesh
+        The Mesh object that events will be tallied onto
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+
+    @property
+    def num_bins(self):
+        n_dim = len(self.mesh.dimension)
+        return 4*n_dim*reduce(operator.mul, self.mesh.dimension)
+
+    def get_bin_index(self, filter_bin):
+        # Split bin into mesh/surface parts
+        *mesh_tuple, surf = filter_bin
+
+        # Get index for mesh alone
+        mesh_index = super().get_bin_index(mesh_tuple)
+
+        # Combine surface and mesh index
+        n_dim = len(self.mesh.dimension)
+        for surf_index, name in enumerate(_CURRENT_NAMES):
+            if surf == name:
+                return 4*n_dim*mesh_index + surf_index
+        else:
+            raise ValueError("'{}' is not a valid mesh surface.".format(surf))
+
+    def get_bin(self, bin_index):
+        n_dim = len(self.mesh.dimension)
+        mesh_index, surf_index = divmod(bin_index, 4*n_dim)
+        mesh_bin = super().get_bin(mesh_index)
+        return mesh_bin + (_CURRENT_NAMES[surf_index],)
+
+    def get_pandas_dataframe(self, data_size, stride, **kwargs):
+        """Builds a Pandas DataFrame for the Filter's bins.
+
+        This method constructs a Pandas DataFrame object for the filter with
+        columns annotated by filter bin information. This is a helper method for
+        :meth:`Tally.get_pandas_dataframe`.
+
+        Parameters
+        ----------
+        data_size : int
+            The total number of bins in the tally corresponding to this filter
+        stride : int
+            Stride in memory for the filter
+
+        Returns
+        -------
+        pandas.DataFrame
+            A Pandas DataFrame with three columns describing the x,y,z mesh
+            cell indices corresponding to each filter bin.  The number of rows
+            in the DataFrame is the same as the total number of bins in the
+            corresponding tally, with the filter bin appropriately tiled to map
+            to the corresponding tally bins.
+
+        See also
+        --------
+        Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
+
+        """
+        # Initialize Pandas DataFrame
+        df = pd.DataFrame()
+
+        # Initialize dictionary to build Pandas Multi-index column
+        filter_dict = {}
+
+        # Append Mesh ID as outermost index of multi-index
+        mesh_key = 'mesh {}'.format(self.mesh.id)
+
+        # Find mesh dimensions - use 3D indices for simplicity
+        if len(self.mesh.dimension) == 3:
+            nx, ny, nz = self.mesh.dimension
+        elif len(self.mesh.dimension) == 2:
+            nx, ny = self.mesh.dimension
+            nz = 1
+        else:
+            nx = self.mesh.dimension
+            ny = nz = 1
+
+        # Generate multi-index sub-column for x-axis
+        filter_bins = np.arange(1, nx + 1)
+        repeat_factor = 12 * stride
+        filter_bins = np.repeat(filter_bins, repeat_factor)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        filter_dict[(mesh_key, 'x')] = filter_bins
+
+        # Generate multi-index sub-column for y-axis
+        filter_bins = np.arange(1, ny + 1)
+        repeat_factor = 12 * nx * stride
+        filter_bins = np.repeat(filter_bins, repeat_factor)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        filter_dict[(mesh_key, 'y')] = filter_bins
+
+        # Generate multi-index sub-column for z-axis
+        filter_bins = np.arange(1, nz + 1)
+        repeat_factor = 12 * nx * ny * stride
+        filter_bins = np.repeat(filter_bins, repeat_factor)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        filter_dict[(mesh_key, 'z')] = filter_bins
+
+        # Generate multi-index sub-column for surface
+        repeat_factor = stride
+        filter_bins = np.repeat(_CURRENT_NAMES, repeat_factor)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        filter_dict[(mesh_key, 'surf')] = filter_bins
+
+        # Initialize a Pandas DataFrame from the mesh dictionary
+        return pd.concat([df, pd.DataFrame(filter_dict)])
+
+
 class RealFilter(Filter):
     """Tally modifier that describes phase-space and other characteristics
 
diff --git a/openmc/tallies.py b/openmc/tallies.py
index d8bc2136a1..6c055d1913 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -2738,7 +2738,7 @@ class Tally(IDManagerMixin):
                 new_filter = filter_type(bins)
 
                 # Set number of bins manually for mesh/distribcell filters
-                if filter_type in (openmc.DistribcellFilter, openmc.MeshFilter):
+                if filter_type is openmc.DistribcellFilter:
                     new_filter._num_bins = find_filter._num_bins
 
                 # Replace existing filter with new one
diff --git a/src/api.F90 b/src/api.F90
index d79a32d945..df183c292f 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -75,6 +75,7 @@ module openmc_api
   public :: openmc_material_filter_get_bins
   public :: openmc_material_filter_set_bins
   public :: openmc_mesh_filter_set_mesh
+  public :: openmc_meshsurface_filter_set_mesh
   public :: openmc_next_batch
   public :: openmc_nuclide_name
   public :: openmc_plot_geometry
@@ -275,8 +276,9 @@ contains
     ! Clear active tally lists
     call active_analog_tallies % clear()
     call active_tracklength_tallies % clear()
-    call active_current_tallies % clear()
+    call active_meshsurf_tallies % clear()
     call active_collision_tallies % clear()
+    call active_surface_tallies % clear()
     call active_tallies % clear()
 
     ! Reset timers
diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index 521be5e3f2..3152f1eb2b 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -73,12 +73,10 @@ contains
     integer :: ital          ! tally object index
     integer :: ijk(3)        ! indices for mesh cell
     integer :: score_index   ! index to pull from tally object
-    integer :: i_filt        ! index in filters array
     integer :: i_filter_mesh ! index for mesh filter
     integer :: i_filter_ein  ! index for incoming energy filter
     integer :: i_filter_eout ! index for outgoing energy filter
-    integer :: i_filter_surf ! index for surface filter
-    integer :: stride_surf   ! stride for surface filter
+    integer :: i_mesh        ! flattend index for mesh
     logical :: energy_filters! energy filters present
     real(8) :: flux          ! temp variable for flux
     type(RegularMesh), pointer :: m ! pointer for mesh object
@@ -95,10 +93,10 @@ contains
 
     ! Associate tallies and mesh
     associate (t => cmfd_tallies(1) % obj)
-      i_filt = t % filter(t % find_filter(FILTER_MESH))
+      i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
     end associate
 
-    select type(filt => filters(i_filt) % obj)
+    select type(filt => filters(i_filter_mesh) % obj)
     type is (MeshFilter)
       m => meshes(filt % mesh)
     end select
@@ -115,16 +113,16 @@ contains
 
       ! Associate tallies and mesh
       associate (t => cmfd_tallies(ital) % obj)
-      i_filt = t % filter(t % find_filter(FILTER_MESH))
-      select type(filt => filters(i_filt) % obj)
-      type is (MeshFilter)
-        m => meshes(filt % mesh)
-      end select
+
+      if (ital < 3) then
+        i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
+      else
+        i_filter_mesh = t % filter(t % find_filter(FILTER_MESHSURFACE))
+      end if
 
       ! Check for energy filters
       energy_filters = (t % find_filter(FILTER_ENERGYIN) > 0)
 
-      i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
       if (energy_filters) then
         i_filter_ein  = t % filter(t % find_filter(FILTER_ENERGYIN))
         i_filter_eout = t % filter(t % find_filter(FILTER_ENERGYOUT))
@@ -247,65 +245,62 @@ contains
 
               else if (ital == 3) then
 
-                i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
-                stride_surf = t % stride(t % find_filter(FILTER_SURFACE))
-
                 ! Initialize and filter for energy
                 do l = 1, size(t % filter)
                   call filter_matches(t % filter(l)) % bins % clear()
                   call filter_matches(t % filter(l)) % bins % push_back(1)
                 end do
+
+                ! Set the bin for this mesh cell
+                i_mesh = m % get_bin_from_indices([ i, j, k ])
+                filter_matches(i_filter_mesh) % bins % data(1) = 12*(i_mesh - 1) + 1
+
+                ! Set the energy bin if needed
                 if (energy_filters) then
                   filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
                 end if
 
-                ! Get the bin for this mesh cell
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices([ i, j, k ])
-
-                score_index = 1
+                score_index = 0
                 do l = 1, size(t % filter)
-                  if (t % filter(l) == i_filter_surf) cycle
                   score_index = score_index + (filter_matches(t % filter(l)) &
                        % bins % data(1) - 1) * t % stride(l)
                 end do
 
                 ! Left surface
                 cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_LEFT - 1) * stride_surf)
+                     score_index + OUT_LEFT)
                 cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_LEFT - 1) * stride_surf)
+                     score_index + IN_LEFT)
 
                 ! Right surface
                 cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_RIGHT - 1) * stride_surf)
+                     score_index + IN_RIGHT)
                 cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_RIGHT - 1) * stride_surf)
+                     score_index + OUT_RIGHT)
 
                 ! Back surface
                 cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_BACK - 1) * stride_surf)
+                     score_index + OUT_BACK)
                 cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_BACK - 1) * stride_surf)
+                     score_index + IN_BACK)
 
                 ! Front surface
                 cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_FRONT - 1) * stride_surf)
+                     score_index + IN_FRONT)
                 cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_FRONT - 1) * stride_surf)
+                     score_index + OUT_FRONT)
 
-                ! Bottom surface
+                ! Left surface
                 cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_BOTTOM - 1) * stride_surf)
+                     score_index + OUT_BOTTOM)
                 cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_BOTTOM - 1) * stride_surf)
+                     score_index + IN_BOTTOM)
 
-                ! Top surface
+                ! Right surface
                 cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (IN_TOP - 1) * stride_surf)
+                     score_index + IN_TOP)
                 cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, &
-                     score_index + (OUT_TOP - 1) * stride_surf)
-
+                     score_index + OUT_TOP)
               end if TALLY
 
             end do OUTGROUP
diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90
index dbfabb2540..86f4e2400e 100644
--- a/src/cmfd_input.F90
+++ b/src/cmfd_input.F90
@@ -373,7 +373,7 @@ contains
 
     ! Determine number of filters
     energy_filters = check_for_node(node_mesh, "energy")
-    n = merge(5, 3, energy_filters)
+    n = merge(4, 2, energy_filters)
 
     ! Extend filters array so we can add CMFD filters
     err = openmc_extend_filters(n, i_filt_start, i_filt_end)
@@ -409,44 +409,15 @@ contains
     ! Duplicate the mesh filter for the mesh current tally since other
     ! tallies use this filter and we need to change the dimension
     i_filt = i_filt + 1
-    err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR)
+    err = openmc_filter_set_type(i_filt, C_CHAR_'meshsurface' // C_NULL_CHAR)
     call openmc_get_filter_next_id(filt_id)
     err = openmc_filter_set_id(i_filt, filt_id)
-    err = openmc_mesh_filter_set_mesh(i_filt, i_start)
+    err = openmc_meshsurface_filter_set_mesh(i_filt, i_start)
 
-    ! We need to increase the dimension by one since we also need
-    ! currents coming into and out of the boundary mesh cells.
-    filters(i_filt) % obj % n_bins = product(m % dimension + 1)
-
-    ! Set up surface filter
-    i_filt = i_filt + 1
-    allocate(SurfaceFilter :: filters(i_filt) % obj)
-    select type(filt => filters(i_filt) % obj)
-    type is(SurfaceFilter)
-      filt % id = i_filt
-      filt % n_bins = 4 * m % n_dimension
-      allocate(filt % surfaces(4 * m % n_dimension))
-      if (m % n_dimension == 2) then
-        filt % surfaces = (/ OUT_LEFT, IN_LEFT, IN_RIGHT, OUT_RIGHT, &
-             OUT_BACK, IN_BACK, IN_FRONT, OUT_FRONT /)
-      elseif (m % n_dimension == 3) then
-        filt % surfaces = (/ OUT_LEFT, IN_LEFT, IN_RIGHT, OUT_RIGHT, &
-             OUT_BACK, IN_BACK, IN_FRONT, OUT_FRONT, &
-             OUT_BOTTOM, IN_BOTTOM, IN_TOP, OUT_TOP /)
-      end if
-      filt % current = .true.
-      ! Add filter to dictionary
-      call filter_dict % set(filt % id, i_filt)
-    end select
 
     ! Initialize filters
     do i = i_filt_start, i_filt_end
-      select type (filt => filters(i) % obj)
-      type is (SurfaceFilter)
-        ! Don't do anything
-      class default
-        call filt % initialize()
-      end select
+      call filters(i) % obj % initialize()
     end do
 
     ! Allocate tallies
@@ -564,13 +535,9 @@ contains
         ! Set tally estimator to analog
         t % estimator = ESTIMATOR_ANALOG
 
-        ! Set the surface filter index in the tally find_filter array
-        n_filter = n_filter + 1
-
         ! Allocate and set filters
         allocate(filter_indices(n_filter))
-        filter_indices(1) = i_filt_end - 1
-        filter_indices(n_filter) = i_filt_end
+        filter_indices(1) = i_filt_end
         if (energy_filters) then
           filter_indices(2) = i_filt_start + 1
         end if
@@ -588,7 +555,7 @@ contains
 
         ! Set macro bins
         t % score_bins(1) = SCORE_CURRENT
-        t % type = TALLY_MESH_CURRENT
+        t % type = TALLY_MESH_SURFACE
       end if
 
       ! Make CMFD tallies active from the start
diff --git a/src/constants.F90 b/src/constants.F90
index 3334a4cc76..0e138f0ee5 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -292,7 +292,7 @@ module constants
   ! Tally type
   integer, parameter :: &
        TALLY_VOLUME          = 1, &
-       TALLY_MESH_CURRENT    = 2, &
+       TALLY_MESH_SURFACE    = 2, &
        TALLY_SURFACE         = 3
 
   ! Tally estimator types
@@ -358,7 +358,7 @@ module constants
   integer, parameter :: NO_BIN_FOUND = -1
 
   ! Tally filter and map types
-  integer, parameter :: N_FILTER_TYPES = 15
+  integer, parameter :: N_FILTER_TYPES = 16
   integer, parameter :: &
        FILTER_UNIVERSE       = 1,  &
        FILTER_MATERIAL       = 2,  &
@@ -374,7 +374,8 @@ module constants
        FILTER_AZIMUTHAL      = 12, &
        FILTER_DELAYEDGROUP   = 13, &
        FILTER_ENERGYFUNCTION = 14, &
-       FILTER_CELLFROM       = 15
+       FILTER_CELLFROM       = 15, &
+       FILTER_MESHSURFACE    = 16
 
   ! Mesh types
   integer, parameter :: &
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 918d471d3c..a29bb08e83 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -10,7 +10,7 @@ module input_xml
   use distribution_multivariate
   use distribution_univariate
   use endf,             only: reaction_name
-  use error,            only: fatal_error, warning, write_message
+  use error,            only: fatal_error, warning, write_message, openmc_err_msg
   use geometry,         only: calc_offsets, maximum_levels, count_instance, &
                               neighbor_lists
   use geometry_header
@@ -2197,7 +2197,6 @@ contains
     integer :: l             ! another loop index
     integer :: filter_id     ! user-specified identifier for filter
     integer :: i_filt        ! index in filters array
-    integer :: i_filter_mesh ! index of mesh filter
     integer :: i_elem        ! index of entry in dictionary
     integer :: n             ! size of arrays in mesh specification
     integer :: n_words       ! number of words read
@@ -2209,7 +2208,6 @@ contains
     integer :: trig_ind      ! index of triggers array for each tally
     integer :: user_trig_ind ! index of user-specified triggers for each tally
     integer :: i_start, i_end
-    integer :: i_filt_start, i_filt_end
     integer(C_INT) :: err
     real(8) :: threshold     ! trigger convergence threshold
     integer :: n_order       ! moment order requested
@@ -2358,13 +2356,13 @@ contains
            call get_node_value(node_filt, "type", temp_str)
       temp_str = to_lower(temp_str)
 
-      ! Determine number of bins
+      ! Make sure bins have been set
       select case(temp_str)
       case ("energy", "energyout", "mu", "polar", "azimuthal")
         if (.not. check_for_node(node_filt, "bins")) then
           call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
         end if
-      case ("mesh", "universe", "material", "cell", "distribcell", &
+      case ("mesh", "meshsurface", "universe", "material", "cell", "distribcell", &
             "cellborn", "cellfrom", "surface", "delayedgroup")
         if (.not. check_for_node(node_filt, "bins")) then
           call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
@@ -2373,6 +2371,7 @@ contains
 
       ! Allocate according to the filter type
       err = openmc_filter_set_type(i_start + i - 1, to_c_string(temp_str))
+      if (err /= 0) call fatal_error(to_f_string(openmc_err_msg))
 
       ! Read filter data from XML
       call f % obj % from_xml(node_filt)
@@ -2835,18 +2834,18 @@ contains
                  &t % find_filter(FILTER_CELL) > 0 .or. &
                  &t % find_filter(FILTER_CELLFROM) > 0) then
 
-              ! Check to make sure that mesh currents are not desired as well
-              if (t % find_filter(FILTER_MESH) > 0) then
-                call fatal_error("Cannot tally other mesh currents &
-                     &in the same tally as surface currents")
+              ! Check to make sure that mesh surface currents are not desired as well
+              if (t % find_filter(FILTER_MESHSURFACE) > 0) then
+                call fatal_error("Cannot tally mesh surface currents &
+                     &in the same tally as normal surface currents")
               end if
 
               t % type = TALLY_SURFACE
               t % score_bins(j) = SCORE_CURRENT
 
-            else if (t % find_filter(FILTER_MESH) > 0) then
+            else if (t % find_filter(FILTER_MESHSURFACE) > 0) then
               t % score_bins(j) = SCORE_CURRENT
-              t % type = TALLY_MESH_CURRENT
+              t % type = TALLY_MESH_SURFACE
 
               ! Check to make sure that current is the only desired response
               ! for this tally
@@ -2854,75 +2853,6 @@ contains
                 call fatal_error("Cannot tally other scores in the &
                      &same tally as surface currents")
               end if
-
-              ! Get index of mesh filter
-              i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
-
-              ! Check to make sure mesh filter was specified
-              if (i_filter_mesh == 0) then
-                call fatal_error("Cannot tally surface current without a mesh &
-                     &filter.")
-              end if
-
-              ! Get pointer to mesh
-              select type(filt => filters(i_filter_mesh) % obj)
-              type is (MeshFilter)
-                m => meshes(filt % mesh)
-              end select
-
-              ! Extend the filters array so we can add a surface filter and
-              ! mesh filter
-              err = openmc_extend_filters(2, i_filt_start, i_filt_end)
-
-              ! Duplicate the mesh filter since other tallies might use this
-              ! filter and we need to change the dimension
-              filters(i_filt_start) = filters(i_filter_mesh)
-
-              ! We need to increase the dimension by one since we also need
-              ! currents coming into and out of the boundary mesh cells.
-              filters(i_filt_start) % obj % n_bins = product(m % dimension + 1)
-
-              ! Set ID
-              call openmc_get_filter_next_id(filter_id)
-              err = openmc_filter_set_id(i_filt_start, filter_id)
-
-
-              ! Add surface filter
-              allocate(SurfaceFilter :: filters(i_filt_end) % obj)
-              select type (filt => filters(i_filt_end) % obj)
-              type is (SurfaceFilter)
-                filt % n_bins = 4 * m % n_dimension
-                allocate(filt % surfaces(4 * m % n_dimension))
-                if (m % n_dimension == 1) then
-                  filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT /)
-                elseif (m % n_dimension == 2) then
-                  filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, &
-                       OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT /)
-                elseif (m % n_dimension == 3) then
-                  filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, &
-                       OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT, OUT_BOTTOM, &
-                       IN_BOTTOM, OUT_TOP, IN_TOP /)
-                end if
-                filt % current = .true.
-
-                ! Set ID
-                call openmc_get_filter_next_id(filter_id)
-                err = openmc_filter_set_id(i_filt_end, filter_id)
-              end select
-
-              ! Copy filter indices to resized array
-              n_filter = size(t % filter)
-              allocate(temp_filter(n_filter + 1))
-              temp_filter(1:size(t % filter)) = t % filter
-              n_filter = n_filter + 1
-
-              ! Set mesh and surface filters
-              temp_filter(t % find_filter(FILTER_MESH)) = i_filt_start
-              temp_filter(n_filter) = i_filt_end
-
-              ! Set filters
-              err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter)
-              deallocate(temp_filter)
             end if
 
           case ('events')
diff --git a/src/output.F90 b/src/output.F90
index 964289126c..b6809ca82f 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -771,12 +771,6 @@ contains
         end associate
       end if
 
-      ! Handle surface current tallies separately
-      if (t % type == TALLY_MESH_CURRENT) then
-        call write_surface_current(t, unit_tally)
-        cycle
-      end if
-
       ! WARNING: Admittedly, the logic for moving for printing results is
       ! extremely confusing and took quite a bit of time to get correct. The
       ! logic is structured this way since it is not practical to have a do
@@ -938,188 +932,6 @@ contains
 
   end subroutine write_tallies
 
-!===============================================================================
-! WRITE_SURFACE_CURRENT writes out surface current tallies over a mesh to the
-! tallies.out file.
-!===============================================================================
-
-  subroutine write_surface_current(t, unit_tally)
-    type(TallyObject), intent(in) :: t
-    integer, intent(in) :: unit_tally
-
-    integer :: i                    ! mesh index
-    integer :: j                    ! loop index over tally filters
-    integer :: ijk(3)               ! indices of mesh cells
-    integer :: n_dim                ! number of mesh dimensions
-    integer :: n_cells              ! number of mesh cells
-    integer :: l                    ! index for energy
-    integer :: i_filter_mesh        ! index for mesh filter
-    integer :: i_filter_ein         ! index for incoming energy filter
-    integer :: i_filter_surf        ! index for surface filter
-    integer :: stride_surf          ! stride for surface filter
-    integer :: n                    ! number of incoming energy bins
-    integer :: filter_index         ! index in results array for filters
-    integer :: nr                   ! number of realizations
-    real(8) :: x(2)                 ! mean and standard deviation
-    logical :: print_ebin           ! should incoming energy bin be displayed?
-    logical :: energy_filters       ! energy filters present
-    character(MAX_LINE_LEN) :: string
-    type(RegularMesh), pointer :: m
-    type(TallyFilterMatch), allocatable :: matches(:)
-
-    allocate(matches(n_filters))
-
-    nr = t % n_realizations
-
-    ! Get pointer to mesh
-    i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
-    select type(filt => filters(i_filter_mesh) % obj)
-    type is (MeshFilter)
-      m => meshes(filt % mesh)
-    end select
-
-    ! Get surface filter index and stride
-    i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
-    stride_surf = t % stride(t % find_filter(FILTER_SURFACE))
-
-    ! initialize bins array
-    do j = 1, size(t % filter)
-      call matches(t % filter(j)) % bins % clear()
-      call matches(t % filter(j)) % bins % push_back(1)
-    end do
-
-    ! determine how many energy in bins there are
-    energy_filters = (t % find_filter(FILTER_ENERGYIN) > 0)
-    if (energy_filters) then
-      print_ebin = .true.
-      i_filter_ein = t % filter(t % find_filter(FILTER_ENERGYIN))
-      n = filters(i_filter_ein) % obj % n_bins
-    else
-      print_ebin = .false.
-      n = 1
-    end if
-
-    ! Get the dimensions and number of cells in the mesh
-    n_dim = m % n_dimension
-    n_cells = product(m % dimension)
-
-    ! Loop over all the mesh cells
-    do i = 1, n_cells
-
-      ! Get the indices for this cell
-      call m % get_indices_from_bin(i, ijk)
-      matches(i_filter_mesh) % bins % data(1) = i
-
-      ! Write the header for this cell
-      if (n_dim == 1) then
-        string = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
-      else if (n_dim == 2) then
-        string = "Mesh Index (" // trim(to_str(ijk(1))) // ", " &
-             // trim(to_str(ijk(2))) // ")"
-      else if (n_dim == 3) then
-        string = "Mesh Index (" // trim(to_str(ijk(1))) // ", " &
-             // trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
-      end if
-
-      write(UNIT=unit_tally, FMT='(1X,A)') trim(string)
-
-      do l = 1, n
-        if (print_ebin) then
-          ! Set incoming energy bin
-          matches(i_filter_ein) % bins % data(1) = l
-
-          ! Write incoming energy bin
-          write(UNIT=unit_tally, FMT='(3X,A)') &
-               trim(filters(i_filter_ein) % obj % text_label( &
-               matches(i_filter_ein) % bins % data(1)))
-        end if
-
-        filter_index = 1
-        do j = 1, size(t % filter)
-          if (t % filter(j) == i_filter_surf) cycle
-          filter_index = filter_index + (matches(t % filter(j)) &
-               % bins % data(1) - 1) * t % stride(j)
-        end do
-
-        associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
-
-        ! Left Surface
-        x(:) = mean_stdev(r(:, 1, filter_index + (OUT_LEFT - 1) * &
-             stride_surf), nr)
-        write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-             "Outgoing Current on Left", to_str(x(1)), trim(to_str(x(2)))
-
-        x(:) = mean_stdev(r(:, 1, filter_index + (IN_LEFT - 1) * &
-             stride_surf), nr)
-        write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-             "Incoming Current on Left", to_str(x(1)), trim(to_str(x(2)))
-
-        ! Right Surface
-        x(:) = mean_stdev(r(:, 1, filter_index + (OUT_RIGHT - 1) * &
-             stride_surf), nr)
-        write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-             "Outgoing Current on Right", to_str(x(1)), trim(to_str(x(2)))
-
-        x(:) = mean_stdev(r(:, 1, filter_index + (IN_RIGHT - 1) * &
-             stride_surf), nr)
-        write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-             "Incoming Current on Right", to_str(x(1)), trim(to_str(x(2)))
-
-        if (n_dim >= 2) then
-
-          ! Back Surface
-          x(:) = mean_stdev(r(:, 1, filter_index + (OUT_BACK - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Outgoing Current on Back", to_str(x(1)), trim(to_str(x(2)))
-
-          x(:) = mean_stdev(r(:, 1, filter_index + (IN_BACK - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Incoming Current on Back", to_str(x(1)), trim(to_str(x(2)))
-
-          ! Front Surface
-          x(:) = mean_stdev(r(:, 1, filter_index + (OUT_FRONT - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Outgoing Current on Front", to_str(x(1)), trim(to_str(x(2)))
-
-          x(:) = mean_stdev(r(:, 1, filter_index + (IN_FRONT - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Incoming Current on Front", to_str(x(1)), trim(to_str(x(2)))
-        end if
-
-        if (n_dim == 3) then
-
-          ! Bottom Surface
-          x(:) = mean_stdev(r(:, 1, filter_index + (OUT_BOTTOM - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Outgoing Current on Bottom", to_str(x(1)), trim(to_str(x(2)))
-
-          x(:) = mean_stdev(r(:, 1, filter_index + (IN_BOTTOM - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Incoming Current on Bottom", to_str(x(1)), trim(to_str(x(2)))
-
-          ! Top Surface
-          x(:) = mean_stdev(r(:, 1, filter_index + (OUT_TOP - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Outgoing Current on Top", to_str(x(1)), trim(to_str(x(2)))
-
-          x(:) = mean_stdev(r(:, 1, filter_index + (IN_TOP - 1) * &
-               stride_surf), nr)
-          write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
-               "Incoming Current on Top", to_str(x(1)), trim(to_str(x(2)))
-        end if
-        end associate
-      end do
-    end do
-
-  end subroutine write_surface_current
-
 !===============================================================================
 ! MEAN_STDEV computes the sample mean and standard deviation of the mean of a
 ! single tally score
diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc
index 5ab8971e6c..204284c480 100644
--- a/src/relaxng/tallies.rnc
+++ b/src/relaxng/tallies.rnc
@@ -35,14 +35,14 @@ element tallies {
   element filter {
     (element id { xsd:int } | attribute id { xsd:int }) &
     (
-      ( (element type { ( "cell" | "cellfrom" | "cellborn" | "material" | 
-          "universe" | "surface" | "distribcell" | "mesh" | "energy" | 
-          "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | 
-          "energyfunction") } |
-         attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" | 
-          "universe" | "surface" | "distribcell" | "mesh" | "energy" | 
-          "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | 
-          "energyfunction") }) &
+      ( (element type { ( "cell" | "cellfrom" | "cellborn" | "material" |
+          "universe" | "surface" | "distribcell" | "mesh" | "energy" |
+          "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" |
+          "energyfunction" | "meshsurface") } |
+         attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" |
+          "universe" | "surface" | "distribcell" | "mesh" | "energy" |
+          "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" |
+          "energyfunction" | "meshsurface") }) &
         (element bins { list { xsd:double+ } } |
           attribute bins { list { xsd:double+ } })
       ) |
diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng
index bad2fdce91..15b6f5b249 100644
--- a/src/relaxng/tallies.rng
+++ b/src/relaxng/tallies.rng
@@ -182,6 +182,7 @@
                     <value>azimuthal</value>
                     <value>delayedgroup</value>
                     <value>energyfunction</value>
+                    <value>meshsurface</value>
                   </choice>
                 </element>
                 <attribute name="type">
@@ -201,6 +202,7 @@
                     <value>azimuthal</value>
                     <value>delayedgroup</value>
                     <value>energyfunction</value>
+                    <value>meshsurface</value>
                   </choice>
                 </attribute>
               </choice>
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index a9a5fd8bfc..f4d12be7e1 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -3064,9 +3064,9 @@ contains
 ! tally total or partial currents between two cells
 !===============================================================================
 
-    subroutine  score_surface_tally(p)
-
-    type(Particle), intent(in)    :: p
+  subroutine score_surface_tally(p, tally_vec)
+    type(Particle), intent(in)  :: p
+    type(VectorInt), intent(in) :: tally_vec
 
     integer :: i
     integer :: i_tally
@@ -3086,9 +3086,9 @@ contains
     ! No collision, so no weight change when survival biasing
     flux = p % wgt
 
-    TALLY_LOOP: do i = 1, active_surface_tallies % size()
+    TALLY_LOOP: do i = 1, tally_vec % size()
       ! Get index of tally and pointer to tally
-      i_tally = active_surface_tallies % data(i)
+      i_tally = tally_vec % data(i)
       associate (t => tallies(i_tally) % obj)
 
       ! Find all valid bins in each filter if they have not already been found
@@ -3188,290 +3188,6 @@ contains
 
   end subroutine score_surface_tally
 
-!===============================================================================
-! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually
-! determining which mesh surfaces were crossed
-!===============================================================================
-
-  subroutine score_surface_current(p)
-
-    type(Particle), intent(in) :: p
-
-    integer :: i
-    integer :: i_tally
-    integer :: j, k                 ! loop indices
-    integer :: n_dim                ! num dimensions of the mesh
-    integer :: d1                   ! dimension index
-    integer :: d2                   ! dimension index
-    integer :: d3                   ! dimension index
-    integer :: ijk0(3)              ! indices of starting coordinates
-    integer :: ijk1(3)              ! indices of ending coordinates
-    integer :: n_cross              ! number of surface crossings
-    integer :: filter_index         ! index of scoring bin
-    integer :: i_filter_mesh        ! index of mesh filter in filters array
-    integer :: i_filter_surf        ! index of surface filter in filters
-    integer :: i_filter_energy      ! index of energy filter in filters
-    real(8) :: uvw(3)               ! cosine of angle of particle
-    real(8) :: xyz0(3)              ! starting/intermediate coordinates
-    real(8) :: xyz1(3)              ! ending coordinates of particle
-    real(8) :: xyz_cross(3)         ! coordinates of bounding surfaces
-    real(8) :: d(3)                 ! distance to each bounding surface
-    real(8) :: distance             ! actual distance traveled
-    integer :: matching_bin         ! next valid filter bin
-    logical :: start_in_mesh        ! particle's starting xyz in mesh?
-    logical :: end_in_mesh          ! particle's ending xyz in mesh?
-    logical :: cross_surface        ! whether the particle crosses a surface
-    logical :: energy_filter        ! energy filter present
-    type(RegularMesh), pointer :: m
-
-    TALLY_LOOP: do i = 1, active_current_tallies % size()
-      ! Copy starting and ending location of particle
-      xyz0 = p % last_xyz_current
-      xyz1 = p % coord(1) % xyz
-
-      ! Get pointer to tally
-      i_tally = active_current_tallies % data(i)
-      associate (t => tallies(i_tally) % obj)
-
-      ! Check for energy filter
-      energy_filter = (t % find_filter(FILTER_ENERGYIN) > 0)
-
-      ! Get index for mesh, surface, and energy filters
-      i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
-      i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
-      if (energy_filter) then
-        i_filter_energy = t % filter(t % find_filter(FILTER_ENERGYIN))
-      end if
-
-      ! Reset the matching bins arrays
-      call filter_matches(i_filter_mesh) % bins % resize(1)
-      call filter_matches(i_filter_surf) % bins % resize(1)
-      if (energy_filter) then
-        call filter_matches(i_filter_energy) % bins % resize(1)
-      end if
-
-      ! Get pointer to mesh
-      select type(filt => filters(i_filter_mesh) % obj)
-      type is (MeshFilter)
-        m => meshes(filt % mesh)
-      end select
-
-      n_dim = m % n_dimension
-
-      ! Determine indices for starting and ending location
-      call m % get_indices(xyz0, ijk0, start_in_mesh)
-      call m % get_indices(xyz1, ijk1, end_in_mesh)
-
-      ! Check to see if start or end is in mesh -- if not, check if track still
-      ! intersects with mesh
-      if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
-        if (.not. m % intersects(xyz0, xyz1)) cycle
-      end if
-
-      ! Calculate number of surface crossings
-      n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
-      if (n_cross == 0) then
-        cycle
-      end if
-
-      ! Copy particle's direction
-      uvw = p % coord(1) % uvw
-
-      ! Determine incoming energy bin.  We need to tell the energy filter this
-      ! is a tracklength tally so it uses the pre-collision energy.
-      if (energy_filter) then
-        call filter_matches(i_filter_energy) % bins % clear()
-        call filter_matches(i_filter_energy) % weights % clear()
-        call filters(i_filter_energy) % obj % get_all_bins(p, &
-             ESTIMATOR_TRACKLENGTH, filter_matches(i_filter_energy))
-        if (filter_matches(i_filter_energy) % bins % size() == 0) cycle
-        matching_bin = filter_matches(i_filter_energy) % bins % data(1)
-        filter_matches(i_filter_energy) % bins % data(1) = matching_bin
-      end if
-
-      ! Bounding coordinates
-      do d1 = 1, n_dim
-        if (uvw(d1) > 0) then
-          xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
-        else
-          xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
-        end if
-      end do
-
-      do j = 1, n_cross
-        ! Reset scoring bin index
-        filter_matches(i_filter_surf) % bins % data(1) = 0
-
-        ! Set the distances to infinity
-        d = INFINITY
-
-        ! Calculate distance to each bounding surface. We need to treat
-        ! special case where the cosine of the angle is zero since this would
-        ! result in a divide-by-zero.
-        do d1 = 1, n_dim
-          if (uvw(d1) == 0) then
-            d(d1) = INFINITY
-          else
-            d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
-          end if
-        end do
-
-        ! Determine the closest bounding surface of the mesh cell by
-        ! calculating the minimum distance. Then use the minimum distance and
-        ! direction of the particle to determine which surface was crossed.
-        distance = minval(d)
-
-        ! Loop over the dimensions
-        do d1 = 1, n_dim
-
-          ! Get the other dimensions.
-          if (d1 == 1) then
-            d2 = mod(d1, 3) + 1
-            d3 = mod(d1 + 1, 3) + 1
-          else
-            d2 = mod(d1 + 1, 3) + 1
-            d3 = mod(d1, 3) + 1
-          end if
-
-          ! Check whether distance is the shortest distance
-          if (distance == d(d1)) then
-
-            ! Check whether particle is moving in positive d1 direction
-            if (uvw(d1) > 0) then
-
-              ! Outward current on d1 max surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 1
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk0)
-                filter_index = 1
-                do k = 1, size(t % filter)
-                  filter_index = filter_index + (filter_matches(t % &
-                       filter(k)) % bins % data(1) - 1) * t % stride(k)
-                end do
-!$omp atomic
-                t % results(RESULT_VALUE, 1, filter_index) = &
-                     t % results(RESULT_VALUE, 1, filter_index) + p % wgt
-              end if
-
-              ! Inward current on d1 min surface
-              cross_surface = .false.
-              select case(n_dim)
-
-              case (1)
-                if (ijk0(d1) >= 0 .and. ijk0(d1) <  m % dimension(d1)) then
-                  cross_surface = .true.
-                end if
-
-              case (2)
-                if (ijk0(d1) >= 0 .and. ijk0(d1) <  m % dimension(d1) .and. &
-                     ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
-                  cross_surface = .true.
-                end if
-
-              case (3)
-                if (ijk0(d1) >= 0 .and. ijk0(d1) <  m % dimension(d1) .and. &
-                     ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
-                     ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
-                  cross_surface = .true.
-                end if
-              end select
-
-              ! If the particle crossed the surface, tally the current
-              if (cross_surface) then
-                ijk0(d1) = ijk0(d1) + 1
-                filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 2
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk0)
-                filter_index = 1
-                do k = 1, size(t % filter)
-                  filter_index = filter_index + (filter_matches(t % &
-                       filter(k)) % bins % data(1) - 1) * t % stride(k)
-                end do
-!$omp atomic
-                t % results(RESULT_VALUE, 1, filter_index) = &
-                     t % results(RESULT_VALUE, 1, filter_index) + p % wgt
-                ijk0(d1) = ijk0(d1) - 1
-              end if
-
-              ijk0(d1) = ijk0(d1) + 1
-              xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
-
-              ! The particle is moving in the negative d1 direction
-            else
-
-              ! Outward current on d1 min surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 3
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk0)
-                filter_index = 1
-                do k = 1, size(t % filter)
-                  filter_index = filter_index + (filter_matches(t % &
-                       filter(k)) % bins % data(1) - 1) * t % stride(k)
-                end do
-!$omp atomic
-                t % results(RESULT_VALUE, 1, filter_index) = &
-                     t % results(RESULT_VALUE, 1, filter_index) + p % wgt
-              end if
-
-              ! Inward current on d1 max surface
-              cross_surface = .false.
-              select case(n_dim)
-
-              case (1)
-                if (ijk0(d1) >  1 .and. ijk0(d1) <= m % dimension(d1) + 1) then
-                  cross_surface = .true.
-                end if
-
-              case (2)
-                if (ijk0(d1) >  1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
-                     ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
-                  cross_surface = .true.
-                end if
-
-              case (3)
-                if (ijk0(d1) >  1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
-                     ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
-                     ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
-                  cross_surface = .true.
-                end if
-              end select
-
-              ! If the particle crossed the surface, tally the current
-              if (cross_surface) then
-                ijk0(d1) = ijk0(d1) - 1
-                filter_matches(i_filter_surf) % bins % data(1) = d1 * 4
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk0)
-                filter_index = 1
-                do k = 1, size(t % filter)
-                  filter_index = filter_index + (filter_matches(t % &
-                       filter(k)) % bins % data(1) - 1) * t % stride(k)
-                end do
-!$omp atomic
-                t % results(RESULT_VALUE, 1, filter_index) = &
-                     t % results(RESULT_VALUE, 1, filter_index) + p % wgt
-                ijk0(d1) = ijk0(d1) + 1
-              end if
-
-              ijk0(d1) = ijk0(d1) - 1
-              xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
-            end if
-          end if
-        end do
-
-        ! Calculate new coordinates
-        xyz0 = xyz0 + distance * uvw
-      end do
-
-      end associate
-    end do TALLY_LOOP
-
-  end subroutine score_surface_current
-
 !===============================================================================
 ! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative
 !===============================================================================
@@ -4375,7 +4091,7 @@ contains
     call active_collision_tallies % clear()
     call active_tracklength_tallies % clear()
     call active_surface_tallies % clear()
-    call active_current_tallies % clear()
+    call active_meshsurf_tallies % clear()
 
     do i = 1, n_tallies
       associate (t => tallies(i) % obj)
@@ -4392,8 +4108,8 @@ contains
             elseif (t % estimator == ESTIMATOR_COLLISION) then
               call active_collision_tallies % push_back(i)
             end if
-          elseif (t % type == TALLY_MESH_CURRENT) then
-            call active_current_tallies % push_back(i)
+          elseif (t % type == TALLY_MESH_SURFACE) then
+            call active_meshsurf_tallies % push_back(i)
           elseif (t % type == TALLY_SURFACE) then
             call active_surface_tallies % push_back(i)
           end if
diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90
index 587f8dd902..3d6eaaef21 100644
--- a/src/tallies/tally_filter.F90
+++ b/src/tallies/tally_filter.F90
@@ -19,6 +19,7 @@ module tally_filter
   use tally_filter_energyfunc
   use tally_filter_material
   use tally_filter_mesh
+  use tally_filter_meshsurface
   use tally_filter_mu
   use tally_filter_polar
   use tally_filter_surface
@@ -67,6 +68,8 @@ contains
           type_ = 'material'
         type is (MeshFilter)
           type_ = 'mesh'
+        type is (MeshSurfaceFilter)
+          type_ = 'meshsurface'
         type is (MuFilter)
           type_ = 'mu'
         type is (PolarFilter)
@@ -136,6 +139,8 @@ contains
           allocate(MaterialFilter :: filters(index) % obj)
         case ('mesh')
           allocate(MeshFilter :: filters(index) % obj)
+        case ('meshsurface')
+          allocate(MeshSurfaceFilter :: filters(index) % obj)
         case ('mu')
           allocate(MuFilter :: filters(index) % obj)
         case ('polar')
diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90
index 2092e7e717..3c0e6250c2 100644
--- a/src/tallies/tally_filter_mesh.F90
+++ b/src/tallies/tally_filter_mesh.F90
@@ -84,7 +84,6 @@ contains
     integer :: ijk1(3)              ! indices of ending coordinates
     integer :: search_iter          ! loop count for intersection search
     integer :: bin
-    real(8) :: weight               ! weight to be pushed back
     real(8) :: uvw(3)               ! cosine of angle of particle
     real(8) :: xyz0(3)              ! starting/intermediate coordinates
     real(8) :: xyz1(3)              ! ending coordinates of particle
@@ -96,8 +95,6 @@ contains
     logical :: end_in_mesh          ! ending coordinates inside mesh?
     type(RegularMesh), pointer :: m
 
-    weight = ERROR_REAL
-
     ! Get a pointer to the mesh.
     m => meshes(this % mesh)
     n = m % n_dimension
diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90
new file mode 100644
index 0000000000..e1e8c48033
--- /dev/null
+++ b/src/tallies/tally_filter_meshsurface.F90
@@ -0,0 +1,334 @@
+module tally_filter_meshsurface
+
+  use, intrinsic :: ISO_C_BINDING
+
+  use hdf5
+
+  use constants
+  use dict_header,         only: EMPTY
+  use error
+  use mesh_header,         only: RegularMesh, meshes, n_meshes, mesh_dict
+  use hdf5_interface
+  use particle_header,     only: Particle
+  use string,              only: to_str
+  use tally_filter_header
+  use xml_interface
+
+  implicit none
+  private
+  public :: openmc_meshsurface_filter_set_mesh
+
+!===============================================================================
+! MESHFILTER indexes the location of particle events to a regular mesh.  For
+! tracklength tallies, it will produce multiple valid bins and the bin weight
+! will correspond to the fraction of the track length that lies in that bin.
+!===============================================================================
+
+  type, public, extends(TallyFilter) :: MeshSurfaceFilter
+    integer :: mesh
+  contains
+    procedure :: from_xml
+    procedure :: get_all_bins
+    procedure :: to_statepoint
+    procedure :: text_label
+  end type MeshSurfaceFilter
+
+contains
+
+  subroutine from_xml(this, node)
+    class(MeshSurfaceFilter), intent(inout) :: this
+    type(XMLNode), intent(in) :: node
+
+    integer :: i_mesh
+    integer :: id
+    integer :: n
+    integer :: n_dim
+    integer :: val
+
+    n = node_word_count(node, "bins")
+
+    if (n /= 1) call fatal_error("Only one mesh can be &
+         &specified per meshsurface filter.")
+
+    ! Determine id of mesh
+    call get_node_value(node, "bins", id)
+
+    ! Get pointer to mesh
+    val = mesh_dict % get(id)
+    if (val /= EMPTY) then
+      i_mesh = val
+    else
+      call fatal_error("Could not find mesh " // trim(to_str(id)) &
+           // " specified on filter.")
+    end if
+
+    ! Determine number of bins
+    n_dim = meshes(i_mesh) % n_dimension
+    this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension)
+
+    ! Store the index of the mesh
+    this % mesh = i_mesh
+  end subroutine from_xml
+
+  subroutine get_all_bins(this, p, estimator, match)
+    class(MeshSurfaceFilter), intent(in)  :: this
+    type(Particle),           intent(in)  :: p
+    integer,                  intent(in)  :: estimator
+    type(TallyFilterMatch), intent(inout) :: match
+
+    integer :: j                    ! loop indices
+    integer :: n_dim                ! num dimensions of the mesh
+    integer :: d1                   ! dimension index
+    integer :: ijk0(3)              ! indices of starting coordinates
+    integer :: ijk1(3)              ! indices of ending coordinates
+    integer :: n_cross              ! number of surface crossings
+    integer :: i_mesh               ! flattened mesh bin index
+    integer :: i_surf               ! surface index (1--12)
+    integer :: i_bin                ! actual index for filter
+    real(8) :: uvw(3)               ! cosine of angle of particle
+    real(8) :: xyz0(3)              ! starting/intermediate coordinates
+    real(8) :: xyz1(3)              ! ending coordinates of particle
+    real(8) :: xyz_cross(3)         ! coordinates of bounding surfaces
+    real(8) :: d(3)                 ! distance to each bounding surface
+    real(8) :: distance             ! actual distance traveled
+    logical :: start_in_mesh        ! particle's starting xyz in mesh?
+    logical :: end_in_mesh          ! particle's ending xyz in mesh?
+
+    ! Copy starting and ending location of particle
+    xyz0 = p % last_xyz_current
+    xyz1 = p % coord(1) % xyz
+
+    associate (m => meshes(this % mesh))
+      n_dim = m % n_dimension
+
+      ! Determine indices for starting and ending location
+      call m % get_indices(xyz0, ijk0, start_in_mesh)
+      call m % get_indices(xyz1, ijk1, end_in_mesh)
+
+      ! Check to see if start or end is in mesh -- if not, check if track still
+      ! intersects with mesh
+      if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
+        if (.not. m % intersects(xyz0, xyz1)) return
+      end if
+
+      ! Calculate number of surface crossings
+      n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
+      if (n_cross == 0) return
+
+      ! Copy particle's direction
+      uvw = p % coord(1) % uvw
+
+      ! Bounding coordinates
+      do d1 = 1, n_dim
+        if (uvw(d1) > 0) then
+          xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
+        else
+          xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
+        end if
+      end do
+
+      do j = 1, n_cross
+        ! Set the distances to infinity
+        d = INFINITY
+
+        ! Calculate distance to each bounding surface. We need to treat
+        ! special case where the cosine of the angle is zero since this would
+        ! result in a divide-by-zero.
+        do d1 = 1, n_dim
+          if (uvw(d1) == 0) then
+            d(d1) = INFINITY
+          else
+            d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
+          end if
+        end do
+
+        ! Determine the closest bounding surface of the mesh cell by
+        ! calculating the minimum distance. Then use the minimum distance and
+        ! direction of the particle to determine which surface was crossed.
+        distance = minval(d)
+
+        ! Loop over the dimensions
+        do d1 = 1, n_dim
+
+          ! Check whether distance is the shortest distance
+          if (distance == d(d1)) then
+
+            ! Check whether particle is moving in positive d1 direction
+            if (uvw(d1) > 0) then
+
+              ! Outward current on d1 max surface
+              if (all(ijk0(:n_dim) >= 1) .and. &
+                   all(ijk0(:n_dim) <= m % dimension)) then
+                i_surf = d1 * 4 - 1
+                i_mesh = m % get_bin_from_indices(ijk0)
+                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
+
+                call match % bins % push_back(i_bin)
+                call match % weights % push_back(ONE)
+              end if
+
+              ! Advance position
+              ijk0(d1) = ijk0(d1) + 1
+              xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
+
+              ! If the particle crossed the surface, tally the inward current on
+              ! d1 min surface
+              if (all(ijk0(:n_dim) >= 1) .and. &
+                   all(ijk0(:n_dim) <= m % dimension)) then
+                i_surf = d1 * 4 - 2
+                i_mesh = m % get_bin_from_indices(ijk0)
+                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
+
+                call match % bins % push_back(i_bin)
+                call match % weights % push_back(ONE)
+              end if
+
+            else
+              ! The particle is moving in the negative d1 direction
+
+              ! Outward current on d1 min surface
+              if (all(ijk0(:n_dim) >= 1) .and. &
+                   all(ijk0(:n_dim) <= m % dimension)) then
+                i_surf = d1 * 4 - 3
+                i_mesh = m % get_bin_from_indices(ijk0)
+                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
+
+                call match % bins % push_back(i_bin)
+                call match % weights % push_back(ONE)
+              end if
+
+              ! Advance position
+              ijk0(d1) = ijk0(d1) - 1
+              xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
+
+              ! If the particle crossed the surface, tally the inward current on
+              ! d1 max surface
+              if (all(ijk0(:n_dim) >= 1) .and. &
+                   all(ijk0(:n_dim) <= m % dimension)) then
+                i_surf = d1 * 4
+                i_mesh = m % get_bin_from_indices(ijk0)
+                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
+
+                call match % bins % push_back(i_bin)
+                call match % weights % push_back(ONE)
+              end if
+            end if
+          end if
+        end do
+
+        ! Calculate new coordinates
+        xyz0 = xyz0 + distance * uvw
+      end do
+    end associate
+
+  end subroutine get_all_bins
+
+  subroutine to_statepoint(this, filter_group)
+    class(MeshSurfaceFilter), intent(in) :: this
+    integer(HID_T),           intent(in) :: filter_group
+
+    call write_dataset(filter_group, "type", "meshsurface")
+    call write_dataset(filter_group, "n_bins", this % n_bins)
+    call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
+  end subroutine to_statepoint
+
+  function text_label(this, bin) result(label)
+    class(MeshSurfaceFilter), intent(in) :: this
+    integer,                  intent(in) :: bin
+    character(MAX_LINE_LEN)              :: label
+
+    integer :: i_mesh
+    integer :: i_surf
+    integer :: n_dim
+    integer, allocatable :: ijk(:)
+
+    associate (m => meshes(this % mesh))
+      n_dim = m % n_dimension
+      allocate(ijk(n_dim))
+
+      ! Get flattend mesh index and surface index
+      i_mesh = (bin - 1) / (4*n_dim) + 1
+      i_surf = mod(bin - 1, 4*n_dim) + 1
+
+      ! Get mesh index part of label
+      call m % get_indices_from_bin(i_mesh, ijk)
+      if (m % n_dimension == 1) then
+        label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
+      elseif (m % n_dimension == 2) then
+        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+             trim(to_str(ijk(2))) // ")"
+      elseif (m % n_dimension == 3) then
+        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+             trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
+      end if
+
+      ! Get surface part of label
+      select case (i_surf)
+      case (OUT_LEFT)
+        label = trim(label) // " Outgoing, x-min"
+      case (IN_LEFT)
+        label = trim(label) // " Incoming, x-min"
+      case (OUT_RIGHT)
+        label = trim(label) // " Outgoing, x-max"
+      case (IN_RIGHT)
+        label = trim(label) // " Incoming, x-max"
+      case (OUT_BACK)
+        label = trim(label) // " Outgoing, y-min"
+      case (IN_BACK)
+        label = trim(label) // " Incoming, y-min"
+      case (OUT_FRONT)
+        label = trim(label) // " Outgoing, y-max"
+      case (IN_FRONT)
+        label = trim(label) // " Incoming, y-max"
+      case (OUT_BOTTOM)
+        label = trim(label) // " Outgoing, z-min"
+      case (IN_BOTTOM)
+        label = trim(label) // " Incoming, z-min"
+      case (OUT_TOP)
+        label = trim(label) // " Outgoing, z-max"
+      case (IN_TOP)
+        label = trim(label) // " Incoming, z-max"
+      end select
+    end associate
+  end function text_label
+
+!===============================================================================
+!                               C API FUNCTIONS
+!===============================================================================
+
+  function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) bind(C)
+    ! Set the mesh for a mesh surface filter
+    integer(C_INT32_T), value, intent(in) :: index
+    integer(C_INT32_T), value, intent(in) :: index_mesh
+    integer(C_INT) :: err
+
+    integer :: n_dim
+
+    err = 0
+    if (index >= 1 .and. index <= n_filters) then
+      if (allocated(filters(index) % obj)) then
+        select type (f => filters(index) % obj)
+        type is (MeshSurfaceFilter)
+          if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
+            f % mesh = index_mesh
+            n_dim = meshes(index_mesh) % n_dimension
+            f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension)
+          else
+            err = E_OUT_OF_BOUNDS
+            call set_errmsg("Index in 'meshes' array is out of bounds.")
+          end if
+        class default
+          err = E_INVALID_TYPE
+          call set_errmsg("Tried to set mesh on a non-mesh filter.")
+        end select
+      else
+        err = E_ALLOCATE
+        call set_errmsg("Filter type has not been set yet.")
+      end if
+    else
+      err = E_OUT_OF_BOUNDS
+      call set_errmsg("Index in filters array out of bounds.")
+    end if
+  end function openmc_meshsurface_filter_set_mesh
+
+end module tally_filter_meshsurface
diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90
index 9604d2382a..0495ce14dc 100644
--- a/src/tallies/tally_header.F90
+++ b/src/tallies/tally_header.F90
@@ -144,7 +144,7 @@ module tally_header
   ! Active tally lists
   type(VectorInt), public :: active_analog_tallies
   type(VectorInt), public :: active_tracklength_tallies
-  type(VectorInt), public :: active_current_tallies
+  type(VectorInt), public :: active_meshsurf_tallies
   type(VectorInt), public :: active_collision_tallies
   type(VectorInt), public :: active_tallies
   type(VectorInt), public :: active_surface_tallies
@@ -336,6 +336,8 @@ contains
         j = FILTER_SURFACE
       type is (MeshFilter)
         j = FILTER_MESH
+      type is (MeshSurfaceFilter)
+        j = FILTER_MESHSURFACE
       type is (EnergyFilter)
         j = FILTER_ENERGYIN
       type is (EnergyoutFilter)
@@ -416,7 +418,7 @@ contains
     ! Deallocate tally node lists
     call active_analog_tallies % clear()
     call active_tracklength_tallies % clear()
-    call active_current_tallies % clear()
+    call active_meshsurf_tallies % clear()
     call active_collision_tallies % clear()
     call active_surface_tallies % clear()
     call active_tallies % clear()
diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90
index ee2259ed48..7dcb3d71eb 100644
--- a/src/tallies/trigger.F90
+++ b/src/tallies/trigger.F90
@@ -168,7 +168,7 @@ contains
           trigger % variance = ZERO
 
           ! Mesh current tally triggers require special treatment
-          if (t % type == TALLY_MESH_CURRENT) then
+          if (t % type == TALLY_MESH_SURFACE) then
             call compute_tally_current(t, trigger)
 
           else
diff --git a/src/tracking.F90 b/src/tracking.F90
index f5839eb8e9..3fe6a065b7 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -19,9 +19,9 @@ module tracking
   use surface_header
   use tally_header
   use tally,              only: score_analog_tally, score_tracklength_tally, &
-                                score_collision_tally, score_surface_current, &
-                                score_track_derivative, score_surface_tally, &
-                                score_collision_derivative, zero_flux_derivs
+                                score_collision_tally, score_surface_tally, &
+                                score_track_derivative, zero_flux_derivs, &
+                                score_collision_derivative
   use track_output,       only: initialize_particle_track, write_particle_track, &
                                 add_particle_track, finalize_particle_track
 
@@ -185,7 +185,8 @@ contains
           p % event = EVENT_SURFACE
         end if
         ! Score cell to cell partial currents
-        if(active_surface_tallies % size() > 0) call score_surface_tally(p)
+        if(active_surface_tallies % size() > 0) &
+             call score_surface_tally(p, active_surface_tallies)
       else
         ! ====================================================================
         ! PARTICLE HAS COLLISION
@@ -200,7 +201,8 @@ contains
         ! since the direction of the particle will change and we need to use the
         ! pre-collision direction to figure out what mesh surfaces were crossed
 
-        if (active_current_tallies % size() > 0) call score_surface_current(p)
+        if (active_meshsurf_tallies % size() > 0) &
+             call score_surface_tally(p, active_meshsurf_tallies)
 
         ! Clear surface component
         p % surface = NONE
@@ -317,12 +319,12 @@ contains
       ! forward slightly so that if the mesh boundary is on the surface, it is
       ! still processed
 
-      if (active_current_tallies % size() > 0) then
+      if (active_meshsurf_tallies % size() > 0) then
         ! TODO: Find a better solution to score surface currents than
         ! physically moving the particle forward slightly
 
         p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
-        call score_surface_current(p)
+        call score_surface_tally(p, active_meshsurf_tallies)
       end if
 
       ! Score to global leakage tally
@@ -350,10 +352,10 @@ contains
       ! particle to change -- artificially move the particle slightly back in
       ! case the surface crossing is coincident with a mesh boundary
 
-      if (active_current_tallies % size() > 0) then
+      if (active_meshsurf_tallies % size() > 0) then
         xyz = p % coord(1) % xyz
         p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
-        call score_surface_current(p)
+        call score_surface_tally(p, active_meshsurf_tallies)
         p % coord(1) % xyz = xyz
       end if
 
@@ -407,10 +409,10 @@ contains
       ! Score surface currents since reflection causes the direction of the
       ! particle to change -- artificially move the particle slightly back in
       ! case the surface crossing is coincident with a mesh boundary
-      if (active_current_tallies % size() > 0) then
+      if (active_meshsurf_tallies % size() > 0) then
         xyz = p % coord(1) % xyz
         p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
-        call score_surface_current(p)
+        call score_surface_tally(p, active_meshsurf_tallies)
         p % coord(1) % xyz = xyz
       end if
 
diff --git a/tests/regression_tests/cmfd_feed/results_true.dat b/tests/regression_tests/cmfd_feed/results_true.dat
index d5eed3d3c5..5e6750fe6b 100644
--- a/tests/regression_tests/cmfd_feed/results_true.dat
+++ b/tests/regression_tests/cmfd_feed/results_true.dat
@@ -364,822 +364,6 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
 cmfd indices
 1.000000E+01
 1.000000E+00
diff --git a/tests/regression_tests/cmfd_nofeed/results_true.dat b/tests/regression_tests/cmfd_nofeed/results_true.dat
index a3fbfb9549..f00966cf8f 100644
--- a/tests/regression_tests/cmfd_nofeed/results_true.dat
+++ b/tests/regression_tests/cmfd_nofeed/results_true.dat
@@ -364,822 +364,6 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
-0.000000E+00
 cmfd indices
 1.000000E+01
 1.000000E+00
diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat
index ca063b7ce9..511b298482 100644
--- a/tests/regression_tests/filter_mesh/inputs_true.dat
+++ b/tests/regression_tests/filter_mesh/inputs_true.dat
@@ -327,18 +327,27 @@
   <filter id="1" type="mesh">
     <bins>1</bins>
   </filter>
+  <filter id="4" type="meshsurface">
+    <bins>1</bins>
+  </filter>
   <filter id="2" type="mesh">
     <bins>2</bins>
   </filter>
+  <filter id="5" type="meshsurface">
+    <bins>2</bins>
+  </filter>
   <filter id="3" type="mesh">
     <bins>3</bins>
   </filter>
+  <filter id="6" type="meshsurface">
+    <bins>3</bins>
+  </filter>
   <tally id="1" name="tally 1">
     <filters>1</filters>
     <scores>total</scores>
   </tally>
   <tally id="2" name="tally 2">
-    <filters>1</filters>
+    <filters>4</filters>
     <scores>current</scores>
   </tally>
   <tally id="3" name="tally 3">
@@ -346,7 +355,7 @@
     <scores>total</scores>
   </tally>
   <tally id="4" name="tally 4">
-    <filters>2</filters>
+    <filters>5</filters>
     <scores>current</scores>
   </tally>
   <tally id="5" name="tally 5">
@@ -354,7 +363,7 @@
     <scores>total</scores>
   </tally>
   <tally id="6" name="tally 6">
-    <filters>3</filters>
+    <filters>6</filters>
     <scores>current</scores>
   </tally>
 </tallies>
diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat
index ac3439da50..f331238b70 100644
--- a/tests/regression_tests/filter_mesh/results_true.dat
+++ b/tests/regression_tests/filter_mesh/results_true.dat
@@ -1 +1 @@
-804d161cb8eae506d3247a533d122f44a01d3cedd566b3c65c71a0b51326dd9b97f8bbf45af7304b500476f3f854d91b10ccad94122d15f23641b05b835fada6
\ No newline at end of file
+46950c046648faaa5ff3cb7b4fdd03667ae3c6da96a7ed8121761291de452cf88481f53e967ed52407d77d90109ae24839967a22ae216531a0e1491f5051ea72
\ No newline at end of file
diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index e0c6dd0f2f..8a7f7a6fca 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -30,6 +30,9 @@ class FilterMeshTestHarness(HashedPyAPITestHarness):
         mesh_1d_filter = openmc.MeshFilter(mesh_1d)
         mesh_2d_filter = openmc.MeshFilter(mesh_2d)
         mesh_3d_filter = openmc.MeshFilter(mesh_3d)
+        meshsurf_1d_filter = openmc.MeshSurfaceFilter(mesh_1d)
+        meshsurf_2d_filter = openmc.MeshSurfaceFilter(mesh_2d)
+        meshsurf_3d_filter = openmc.MeshSurfaceFilter(mesh_3d)
 
         # Initialized the tallies
         tally = openmc.Tally(name='tally 1')
@@ -38,7 +41,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness):
         self._model.tallies.append(tally)
 
         tally = openmc.Tally(name='tally 2')
-        tally.filters = [mesh_1d_filter]
+        tally.filters = [meshsurf_1d_filter]
         tally.scores = ['current']
         self._model.tallies.append(tally)
 
@@ -48,7 +51,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness):
         self._model.tallies.append(tally)
 
         tally = openmc.Tally(name='tally 4')
-        tally.filters = [mesh_2d_filter]
+        tally.filters = [meshsurf_2d_filter]
         tally.scores = ['current']
         self._model.tallies.append(tally)
 
@@ -58,7 +61,7 @@ class FilterMeshTestHarness(HashedPyAPITestHarness):
         self._model.tallies.append(tally)
 
         tally = openmc.Tally(name='tally 6')
-        tally.filters = [mesh_3d_filter]
+        tally.filters = [meshsurf_3d_filter]
         tally.scores = ['current']
         self._model.tallies.append(tally)
 
diff --git a/tests/regression_tests/score_current/results_true.dat b/tests/regression_tests/score_current/results_true.dat
index 75d61b7188..dedc473de9 100644
--- a/tests/regression_tests/score_current/results_true.dat
+++ b/tests/regression_tests/score_current/results_true.dat
@@ -1 +1 @@
-4675d5101f4f829369c39cb33d654430836b934ab07c165777ba6e214bdf3a698b8082f4f9bb9e78f1f0e495b30ea02cf9b3d14622c59915d818d678a1e5b7b1
\ No newline at end of file
+138b312cdaa822c9b62f757b3259522004b679c4ed289a92798b29a6442c26d12c53256635be273f13e3703816ff50a1b9f52d79770eade01e482384ac0d389f
\ No newline at end of file
diff --git a/tests/regression_tests/score_current/tallies.xml b/tests/regression_tests/score_current/tallies.xml
index 3b8f60adba..a381e4c2e2 100644
--- a/tests/regression_tests/score_current/tallies.xml
+++ b/tests/regression_tests/score_current/tallies.xml
@@ -9,7 +9,7 @@
   </mesh>
 
   <filter id="1">
-    <type>mesh</type>
+    <type>meshsurface</type>
     <bins>1</bins>
   </filter>