diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index ed0b7c4a3c..c9b543d6f9 100644
--- a/src/cmfd_execute.F90
+++ b/src/cmfd_execute.F90
@@ -25,7 +25,6 @@ contains
 
     use cmfd_data,              only: set_up_cmfd
     use cmfd_message_passing,   only: petsc_init_mpi, cmfd_bcast
-    use cmfd_output,            only: write_hdf5_output
     use cmfd_power_solver,      only: cmfd_power_execute
     use cmfd_snes_solver,       only: cmfd_snes_execute
     use error,                  only: warning, fatal_error 
@@ -120,9 +119,6 @@ contains
     ! calculate weight factors
     if (cmfd_feedback) call cmfd_reweight(.true.)
 
-    ! write output
-    if (cmfd_write_hdf5 .and. master) call write_hdf5_output()
-
     ! stop cmfd timer
     if (master) then
       call time_cmfd % stop()
diff --git a/src/cmfd_output.F90 b/src/cmfd_output.F90
index 7f9b30391b..cb93428b35 100644
--- a/src/cmfd_output.F90
+++ b/src/cmfd_output.F90
@@ -6,7 +6,7 @@ module cmfd_output
 
   implicit none
   private
-  public :: finalize_cmfd, write_hdf5_output 
+  public :: finalize_cmfd
 
 contains
 
@@ -51,81 +51,6 @@ contains
 
   end subroutine write_neutron_balance
 
-!===============================================================================
-! WRITE_HDF5_OUTPUT 
-!===============================================================================
-
-  subroutine write_hdf5_output()
-
-# ifdef HDF5
-
-    use global,       only: cmfd, current_batch
-    use hdf5
-    use string,       only: to_str
-
-    character(100) :: grpname 
-    character(100) :: path
-    integer(HID_T) :: group_id
-    integer(HID_T) :: dataspace_id         ! Data space identifier
-    integer(HID_T) :: dataset_id           ! Dataset identifier
-    integer(HSIZE_T), dimension(1) :: dim1 ! scalar for hdf5 dimensions
-    integer(HSIZE_T), dimension(4) :: dim4 ! vector for hdf5 dimensions
-    integer(Fortran_Integer) :: hdf5_err   ! error
-    integer :: nx         ! number of mesh cells in x direction
-    integer :: ny         ! number of mesh cells in y direction
-    integer :: nz         ! number of mesh cells in z direction
-    integer :: ng         ! number of energy groups
-    integer(hid_t) :: dataType
-    integer :: r8 = 8
-
-    ! get hdf5 data type
-    dataType = h5kind_to_type(r8,H5_REAL_KIND)
-
-    ! extract spatial and energy indices from object
-    nx = cmfd % indices(1)
-    ny = cmfd % indices(2)
-    nz = cmfd % indices(3)
-    ng = cmfd % indices(4)
-
-    ! set group name
-    grpname = 'batch'//trim(to_str(current_batch))
-    ! create a new group for this active batch
-!   call h5gcreate_f(cmfd%file_id,grpname,group_id,hdf5_err)
-
-    ! write out openmc source
-    path = trim(grpname)//'/openmc_src'
-    dim4 = (/ng,nx,ny,nz/)
-    call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err)
-!   call h5dcreate_f(cmfd%file_id,trim(path),dataType,dataspace_id,   &
-!                    dataset_id,hdf5_err)
-    call h5dwrite_f(dataset_id,dataType,cmfd%openmc_src,dim4,hdf5_err)
-    call h5dclose_f(dataset_id,hdf5_err)
-
-    ! write out cmfd source
-    path = trim(grpname)//'/cmfd_src'
-    dim4 = (/ng,nx,ny,nz/)
-    call h5screate_simple_f(4,dim4,dataspace_id,hdf5_err)
-!   call h5dcreate_f(cmfd%file_id,trim(path),dataType,dataspace_id,   &
-!                    dataset_id,hdf5_err)
-    call h5dwrite_f(dataset_id,dataType,cmfd%cmfd_src,dim4,hdf5_err)
-    call h5dclose_f(dataset_id,hdf5_err)
-
-    ! write out keff
-    path = trim(grpname)//'/cmfd_keff'
-    dim1 = (/1/)
-    call h5screate_simple_f(1,dim1,dataspace_id,hdf5_err)
-!   call h5dcreate_f(cmfd%file_id,trim(path),dataType,dataspace_id,   &
-!                    dataset_id,hdf5_err)
-    call h5dwrite_f(dataset_id,dataType,cmfd%keff,dim1,hdf5_err)
-    call h5dclose_f(dataset_id,hdf5_err)
-
-    ! close the group
-    call h5gclose_f(group_id,hdf5_err)
-
-# endif
-
-  end subroutine write_hdf5_output
-
 # endif
 
 end module cmfd_output
diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90
index 1785971e4f..43b7a3ca7f 100644
--- a/src/hdf5_interface.F90
+++ b/src/hdf5_interface.F90
@@ -830,10 +830,12 @@ contains
     integer              :: i, j             ! loop indices
     integer(HSIZE_T)     :: dims(1)          ! dimensions of 1D arrays
     integer(HSIZE_T)     :: dims2(2)         ! dimensions of 2D arrays
+    integer(HSIZE_T)     :: dims4(4)         ! dimensions of 4D arrays
     integer(HID_T)       :: hdf5_state_point ! identifier for state point file
     integer(HID_T)       :: tallies_group    ! "tallies" group
     integer(HID_T)       :: temp_group       ! group for i-th tally or mesh
     integer(HID_T)       :: filter_group     ! group for i-th filter
+    integer(HID_T)       :: cmfd_group       ! group for cmfd output
     integer(HID_T)       :: dspace           ! identifier for dataspace
     integer(HID_T)       :: dset             ! identifier for dataset
     integer, allocatable :: temp_array(:)    ! nuclide bin array
@@ -1113,6 +1115,30 @@ contains
     call h5dclose_f(dset, hdf5_err)
     call h5sclose_f(dspace, hdf5_err)
 
+    ! Write out CMFD info if active
+    if (cmfd_on) then
+
+      ! Create CMFD group
+      call h5gcreate_f(hdf5_state_point, "cmfd", cmfd_group, hdf5_err)
+
+      ! write out openmc source
+      dims4 = shape(cmfd % openmc_src) 
+      call h5ltmake_dataset_double_f(cmfd_group, "openmc_src", 4, &
+           dims4, cmfd % openmc_src, hdf5_err)
+
+      ! write out cmfd source
+      dims4 = shape(cmfd % cmfd_src) 
+      call h5ltmake_dataset_double_f(cmfd_group, "cmfd_src", 4, &
+           dims4, cmfd % cmfd_src, hdf5_err)
+
+      ! write out keff
+      call hdf5_write_double(cmfd_group, "cmfd_keff", cmfd % keff) 
+
+      ! Close CMFD group
+      call h5gclose_f(tallies_group, hdf5_err)
+
+    end if
+
     ! Close HDF5 state point file
     call h5fclose_f(hdf5_state_point, hdf5_err)