diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index 0d873d518c..6a56d25eab 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -27,6 +27,7 @@ cmfd_execute.o: cmfd_snes_solver.o
 cmfd_input.o: datatypes.o
 cmfd_input.o: error.o
 cmfd_input.o: global.o
+cmfd_input.o: hdf5_interface.o
 cmfd_input.o: mesh_header.o
 cmfd_input.o: string.o
 cmfd_input.o: tally_header.o
diff --git a/src/Makefile b/src/Makefile
index 499004aee8..ae778a336e 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -28,7 +28,7 @@ USE_HDF5 = yes
 
 ifeq ($(COMPILER),gnu)
   F90 = gfortran
-  F90FLAGS := -cpp -fbacktrace
+  F90FLAGS := -cpp -fbacktrace -DN2008
   LDFLAGS =
 
   # Debugging
diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index 01e7b75238..6555e24d79 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -65,7 +65,7 @@ contains
     integer :: h                 ! iteration counter for outgoing groups
     integer :: ijk(3)            ! indices for mesh cell
     integer :: score_index       ! index to pull from tally object
-    integer :: bins(TALLY_TYPES) ! bins for filters
+    integer :: bins(N_FILTER_TYPES) ! bins for filters
 
     real(8) :: flux   ! temp variable for flux
 
@@ -103,10 +103,10 @@ contains
             ijk = (/ i, j, k /)
 
             ! get bin number for mesh indices
-            bins(T_MESH) = mesh_indices_to_bin(m,ijk)
+            bins(FILTER_MESH) = mesh_indices_to_bin(m,ijk)
 
             ! apply energy in filter
-            bins(T_ENERGYIN) = ng - h + 1
+            bins(FILTER_ENERGYIN) = ng - h + 1
 
             ! calculate score index from bins
             score_index = sum((bins - 1) * t%stride) + 1
@@ -148,7 +148,7 @@ contains
               t => tallies(2)
 
               ! set energy out bin
-              bins(T_ENERGYOUT) = ng - g + 1
+              bins(FILTER_ENERGYOUT) = ng - g + 1
 
               ! calculate score index from bins
               score_index = sum((bins - 1) * t%stride) + 1
@@ -166,59 +166,59 @@ contains
 
             ! initialize and filter for energy
             bins = 1
-            bins(TS_ENERGYIN) = ng - h + 1
+            bins(SURF_FILTER_ENERGYIN) = ng - h + 1
 
             ! left surface
             bins(1:3) = (/ i-1, j, k /) + 1
-            bins(TS_SURFACE) = IN_RIGHT
+            bins(SURF_FILTER_SURFACE) = IN_RIGHT
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
             cmfd % current(1,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_RIGHT
+            bins(SURF_FILTER_SURFACE) = OUT_RIGHT
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(2,h,i,j,k) = t % scores(score_index,1) % val
 
             ! right surface
             bins(1:3) = (/ i, j, k /) + 1
-            bins(TS_SURFACE) = IN_RIGHT
+            bins(SURF_FILTER_SURFACE) = IN_RIGHT
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(3,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_RIGHT
+            bins(SURF_FILTER_SURFACE) = OUT_RIGHT
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing 
             cmfd % current(4,h,i,j,k) = t % scores(score_index,1) % val
 
             ! back surface
             bins(1:3) = (/ i, j-1, k /) + 1
-            bins(TS_SURFACE) = IN_FRONT
+            bins(SURF_FILTER_SURFACE) = IN_FRONT
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
             cmfd % current(5,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_FRONT
+            bins(SURF_FILTER_SURFACE) = OUT_FRONT
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(6,h,i,j,k) = t % scores(score_index,1) % val
 
             ! front surface
             bins(1:3) = (/ i, j, k /) + 1
-            bins(TS_SURFACE) = IN_FRONT
+            bins(SURF_FILTER_SURFACE) = IN_FRONT
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(7,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_FRONT
+            bins(SURF_FILTER_SURFACE) = OUT_FRONT
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing 
             cmfd % current(8,h,i,j,k) = t % scores(score_index,1) % val
 
             ! bottom surface
             bins(1:3) = (/ i, j, k-1 /) + 1
-            bins(TS_SURFACE) = IN_TOP
+            bins(SURF_FILTER_SURFACE) = IN_TOP
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
             cmfd % current(9,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_TOP
+            bins(SURF_FILTER_SURFACE) = OUT_TOP
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(10,h,i,j,k) = t % scores(score_index,1) % val
 
             ! top surface
             bins(1:3) = (/ i, j, k /) + 1
-            bins(TS_SURFACE) = IN_TOP
+            bins(SURF_FILTER_SURFACE) = IN_TOP
             score_index = sum((bins - 1) * t % stride) + 1 ! incoming 
             cmfd % current(11,h,i,j,k) = t % scores(score_index,1) % val
-            bins(TS_SURFACE) = OUT_TOP
+            bins(SURF_FILTER_SURFACE) = OUT_TOP
             score_index = sum((bins - 1) * t % stride) + 1 ! outgoing 
             cmfd % current(12,h,i,j,k) = t % scores(score_index,1) % val
 
diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90
index 27bbf9c018..f048134a1a 100644
--- a/src/cmfd_input.F90
+++ b/src/cmfd_input.F90
@@ -87,6 +87,8 @@ contains
     integer :: n           ! size of arrays in mesh specification
     integer :: ng=1        ! number of energy groups (default 1)
     integer :: n_words     ! number of words read
+    integer :: n_filters   ! number of filters
+    integer :: filters(N_FILTER_TYPES) ! temp list of filters
     character(MAX_LINE_LEN) :: filename
     character(MAX_WORD_LEN) :: words(MAX_WORDS)
     type(TallyObject),    pointer :: t => null()
@@ -152,14 +154,20 @@ contains
 
     ! begin loop around tallies
     do i = 1,n_tallies
+
+      ! set n filters to 0
+      n_filters = 0
+      filters = 0
+
+      ! point t to tally variable
       t => tallies(i)
 
       ! allocate arrays for number of bins and stride in scores array
-      allocate(t % n_bins(TALLY_TYPES))
-      allocate(t % stride(TALLY_TYPES))
+      allocate(t % n_filter_bins(N_FILTER_TYPES))
+      allocate(t % stride(N_FILTER_TYPES))
 
       ! initialize number of bins and stride
-      t % n_bins = 0
+      t % n_filter_bins = 0
       t % stride = 0
 
       ! record tally id which is equivalent to loop number
@@ -168,7 +176,10 @@ contains
       ! set mesh filter mesh id = 1
       t % mesh = 1
       m => meshes(1)
-      t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
+      t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) +        &
+     &                                 product(m % dimension)
+      n_filters = n_filters + 1
+      filters(n_filters) = FILTER_MESH
 
       ! read and set incoming energy mesh filter
       if (len_trim(mesh_ % energy) > 0) then
@@ -178,22 +189,41 @@ contains
         do j = 1,n_words
           t % energy_in(j) = str_to_real(words(j))
         end do
-        t % n_bins(T_ENERGYIN) = n_words - 1
+        t % n_filter_bins(FILTER_ENERGYIN) = n_words - 1
+        n_filters = n_filters + 1
+        filters(n_filters) = FILTER_ENERGYIN
       end if
 
       if (i == 1) then
 
+        ! set tally estimator to analog
+        t % estimator = ESTIMATOR_ANALOG
+
+        ! set tally type to volume
+        t % type = TALLY_VOLUME
+
+        ! allocate and set filters
+        t % n_filters = n_filters
+        allocate(t % filters(n_filters))
+        t % filters = filters(1:n_filters)
+
         ! allocate macro reactions
-        allocate(t % macro_bins(3))
-        t % n_macro_bins = 3
+        allocate(t % score_bins(3))
+        t % n_score_bins = 3
 
         ! set macro_bins
-        t % macro_bins(1) % scalar = MACRO_FLUX
-        t % macro_bins(2) % scalar = MACRO_TOTAL
-        t % macro_bins(3) % scalar = MACRO_SCATTER_1
+        t % score_bins(1) % scalar = SCORE_FLUX
+        t % score_bins(2) % scalar = SCORE_TOTAL
+        t % score_bins(3) % scalar = SCORE_SCATTER_1
 
       else if (i == 2) then
 
+        ! set tally estimator to analog
+        t % estimator = ESTIMATOR_ANALOG
+
+        ! set tally type to volume
+        t % type = TALLY_VOLUME
+
         ! read and set outgoing energy mesh filter
         if (len_trim(mesh_ % energy) > 0) then
           call split_string(mesh_ % energy, words, n_words)
@@ -201,31 +231,47 @@ contains
           do j = 1, n_words
             t % energy_out(j) = str_to_real(words(j))
           end do
-          t % n_bins(T_ENERGYOUT) = n_words - 1
+          t % n_filter_bins(FILTER_ENERGYOUT) = n_words - 1
+          n_filters = n_filters + 1
+          filters(n_filters) = FILTER_ENERGYOUT
         end if
 
+        ! allocate and set filters
+        t % n_filters = n_filters
+        allocate(t % filters(n_filters))
+        t % filters = filters(1:n_filters)
+
         ! allocate macro reactions
-        allocate(t % macro_bins(2))
-        t % n_macro_bins = 2
+        allocate(t % score_bins(2))
+        t % n_score_bins = 2
 
         ! set macro_bins
-        t % macro_bins(1) % scalar = MACRO_NU_SCATTER
-        t % macro_bins(2) % scalar = MACRO_NU_FISSION
+        t % score_bins(1) % scalar = SCORE_NU_SCATTER
+        t % score_bins(2) % scalar = SCORE_NU_FISSION
 
       else if (i == 3) then
 
+        ! set tally estimator to analog
+        t % estimator = ESTIMATOR_ANALOG
+
+        ! allocate and set filters
+        t % n_filters = n_filters
+        allocate(t % filters(n_filters))
+        t % filters = filters(1:n_filters)
+
         ! allocate macro reactions
-        allocate(t % macro_bins(1))
-        t % n_macro_bins = 1
+        allocate(t % score_bins(1))
+        t % n_score_bins = 1
 
         ! set macro bins
-        t % macro_bins(1) % scalar = MACRO_CURRENT
-        t % surface_current = .true.
+        t % score_bins(1) % scalar = SCORE_CURRENT
+        t % type = TALLY_SURFACE_CURRENT
 
         ! since the number of bins for the mesh filter was already set
         ! assuming it was a flux tally, we need to adjust the number of
         ! bins
-        t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
+        t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH)        &
+       &                                 - product(m % dimension)
 
         ! get pointer to mesh
         id = t % mesh
@@ -234,12 +280,12 @@ contains
 
         ! we need to increase the dimension by one since we also need
         ! currents coming into and out of the boundary mesh cells.
-        if (size(m % dimension) == 2) then
-          t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
-       &                       product(m % dimension + 1) * 4
+        if (size(m % dimension) == 2) then  ! these lines need changing
+          t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) +    &
+         &                                 product(m % dimension + 1) * 4
         elseif (size(m % dimension) == 3) then
-          t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
-                               product(m % dimension + 1) * 6
+          t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) +    &
+         &                                 product(m % dimension + 1) * 6
         end if
 
       end if