From f42fd721444034f25af1b9380cbc1222c53a1ddd Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Mon, 1 Aug 2011 22:17:36 -0400
Subject: [PATCH] Moved types to header modules. Renamed a few important
 arrays. Added OBJECTS and DEPENDENCIES for Makefile.

---
 src/DEPENDENCIES                           |  129 ++
 src/Makefile                               |   80 +-
 src/OBJECTS                                |   28 +
 src/ace.f90                                | 1221 -------------------
 src/bank_header.f90                        |   18 +
 src/constants.f90                          |  177 +++
 src/cross_section.f90                      | 1239 +++++++++++++++++++-
 src/cross_section_header.f90               |  136 +++
 src/{data_structures.f90 => datatypes.f90} |   28 +-
 src/datatypes_header.f90                   |  113 ++
 src/endf.f90                               |    2 +-
 src/energy_grid.f90                        |   32 +-
 src/error.f90                              |    4 +-
 src/fileio.f90                             |   73 +-
 src/geometry.f90                           |   26 +-
 src/geometry_header.f90                    |   77 ++
 src/global.f90                             |  210 +---
 src/logging.f90                            |    5 +-
 src/main.f90                               |   34 +-
 src/material_header.f90                    |   21 +
 src/mpi_routines.f90                       |    9 +-
 src/output.f90                             |   23 +-
 src/particle_header.f90                    |   30 +
 src/physics.f90                            |   87 +-
 src/score.f90                              |   11 +-
 src/search.f90                             |    6 +-
 src/source.f90                             |   18 +-
 src/source_header.f90                      |   15 +
 src/string.f90                             |   16 +-
 src/tally_header.f90                       |   39 +
 src/timing.f90                             |   78 +-
 src/types.f90                              |  415 -------
 src/unittest.f90                           |   92 --
 33 files changed, 2290 insertions(+), 2202 deletions(-)
 create mode 100644 src/DEPENDENCIES
 create mode 100644 src/OBJECTS
 delete mode 100644 src/ace.f90
 create mode 100644 src/bank_header.f90
 create mode 100644 src/constants.f90
 create mode 100644 src/cross_section_header.f90
 rename src/{data_structures.f90 => datatypes.f90} (98%)
 create mode 100644 src/datatypes_header.f90
 create mode 100644 src/geometry_header.f90
 create mode 100644 src/material_header.f90
 create mode 100644 src/particle_header.f90
 create mode 100644 src/source_header.f90
 create mode 100644 src/tally_header.f90
 delete mode 100644 src/types.f90
 delete mode 100644 src/unittest.f90

diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
new file mode 100644
index 0000000000..2c1afbfb66
--- /dev/null
+++ b/src/DEPENDENCIES
@@ -0,0 +1,129 @@
+cross_section.o: constants.o
+cross_section.o: cross_section_header.o
+cross_section.o: datatypes.o
+cross_section.o: datatypes_header.o
+cross_section.o: endf.o
+cross_section.o: error.o
+cross_section.o: fileio.o
+cross_section.o: global.o 
+cross_section.o: material_header.o
+cross_section.o: output.o 
+cross_section.o: string.o 
+
+datatypes.o: datatypes_header.o
+
+endf.o: constants.o 
+endf.o: string.o
+
+error.o: global.o
+
+energy_grid.o: constants.o
+energy_grid.o: datatypes.o
+energy_grid.o: datatypes_header.o
+energy_grid.o: global.o
+energy_grid.o: output.o
+
+fileio.o: constants.o
+fileio.o: cross_section_header.o
+fileio.o: datatypes.o
+fileio.o: datatypes_header.o
+fileio.o: error.o
+fileio.o: geometry_header.o
+fileio.o: global.o
+fileio.o: output.o
+fileio.o: source_header.o
+fileio.o: string.o
+
+geometry.o: constants.o
+geometry.o: datatypes.o
+geometry.o: error.o
+geometry.o: geometry_header.o
+geometry.o: global.o
+geometry.o: output.o
+geometry.o: particle_header.o
+geometry.o: string.o
+
+global.o: bank_header.o
+global.o: constants.o
+global.o: cross_section_header.o
+global.o: datatypes_header.o
+global.o: geometry_header.o
+global.o: material_header.o
+global.o: particle_header.o
+global.o: source_header.o
+global.o: tally_header.o
+global.o: timing.o
+
+logging.o: constants.o
+logging.o: global.o
+
+main.o: constants.o
+main.o: cross_section.o
+main.o: energy_grid.o
+main.o: fileio.o
+main.o: geometry.o
+main.o: geometry_header.o
+main.o: global.o
+main.o: logging.o
+main.o: mcnp_random.o
+main.o: mpi_routines.o
+main.o: output.o
+main.o: particle_header.o
+main.o: physics.o
+main.o: score.o
+main.o: source.o
+main.o: string.o
+main.o: timing.o
+
+mpi_routines.o: constants.o
+mpi_routines.o: error.o
+mpi_routines.o: global.o
+mpi_routines.o: mcnp_random.o
+mpi_routines.o: output.o
+mpi_routines.o: source.o
+
+output.o: constants.o
+output.o: datatypes.o
+output.o: endf.o
+output.o: geometry_header.o
+output.o: global.o
+output.o: string.o
+
+physics.o: constants.o
+physics.o: cross_section_header.o
+physics.o: endf.o
+physics.o: error.o
+physics.o: geometry.o
+physics.o: geometry_header.o
+physics.o: global.o
+physics.o: mcnp_random.o
+physics.o: output.o
+physics.o: particle_header.o
+physics.o: score.o
+physics.o: search.o
+physics.o: string.o
+
+score.o: constants.o
+score.o: cross_section.o
+score.o: error.o
+score.o: global.o
+score.o: output.o
+score.o: search.o
+score.o: string.o
+
+search.o: constants.o
+search.o: error.o
+
+source.o: bank_header.o
+source.o: constants.o
+source.o: cross_section_header.o
+source.o: global.o
+source.o: mcnp_random.o
+source.o: output.o
+source.o: particle_header.o
+source.o: physics.o
+
+string.o: constants.o
+string.o: error.o
+
+timing.o: constants.o
diff --git a/src/Makefile b/src/Makefile
index ddc265f69e..bfef04b267 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,33 +1,15 @@
 program = openmc
 
-main = main.f90
-test = unittest.f90
-modules = 	ace.f90 \
-		cross_section.f90 \
-		data_structures.f90 \
-		endf.f90 \
-		energy_grid.f90 \
-		error.f90 \
-		fileio.f90 \
-		geometry.f90 \
-		global.f90 \
-		logging.f90 \
-		mcnp_random.f90 \
-		mpi_routines.f90 \
-		output.f90 \
-		physics.f90 \
-		search.f90 \
-		source.f90 \
-		string.f90 \
-		score.f90 \
-		timing.f90 \
-		types.f90
+#===============================================================================
+# Object Files
+#===============================================================================
 
-main_objects = $(modules:.f90=.o) $(main:.f90=.o)
-test_objects = $(modules:.f90=.o) $(test:.f90=.o)
+include OBJECTS
 
-#--------------------------------------------------------------------
+#===============================================================================
 # Compiler Options
+#===============================================================================
+
 F90 = ifort
 F90FLAGS = -fpp -g -traceback
 
@@ -39,52 +21,30 @@ ifeq ($(USE_MPI),yes)
   F90FLAGS += -DMPI -f90=ifort
 endif
 
-#--------------------------------------------------------------------
+#===============================================================================
 # Targets
+#===============================================================================
+
 all:	$(program)
-$(program): $(main_objects)
-	$(F90) $(main_objects) -o $@
+$(program): $(objects)
+	$(F90) $(objects) -o $@
 clean:
 	@rm -f *.o *.mod $(program)
 neat:
 	@rm -f *.o *.mod
-unittest: $(test_objects)
-	$(F90) $(test_objects) -o $@
 
-#--------------------------------------------------------------------
-# Rules & misc
+#===============================================================================
+# Rules
+#===============================================================================
+
 .SUFFIXES: .f90 .o
 .PHONY: all clean neat
 
 %.o: %.f90
 	$(F90) -c $(F90FLAGS) $<
 
-#--------------------------------------------------------------------
+#===============================================================================
 # Dependencies
-ace.o: global.o output.o string.o fileio.o string.o endf.o error.o
-cross_section.o: global.o string.o data_structures.o output.o error.o
-data_structures.o: global.o
-endf.o: global.o string.o
-energy_grid.o: global.o output.o data_structures.o
-error.o: global.o
-fileio.o: types.o global.o string.o output.o data_structures.o \
-          error.o
-geometry.o: types.o global.o output.o string.o data_structures.o \
-            error.o
-global.o: types.o
-logging.o: global.o
-main.o: global.o fileio.o output.o geometry.o mcnp_random.o \
-		source.o physics.o cross_section.o data_structures.o \
-		ace.o energy_grid.o score.o logging.o string.o
-mpi_routines.o: global.o output.o types.o mcnp_random.o source.o \
-                error.o
-output.o: global.o types.o data_structures.o endf.o string.o
-physics.o: types.o global.o mcnp_random.o geometry.o output.o \
-		search.o endf.o score.o error.o string.o
-score.o: global.o output.o ace.o error.o string.o
-search.o: global.o error.o
-source.o: global.o mcnp_random.o output.o physics.o
-string.o: global.o error.o
-timing.o: global.o types.o error.o
-types.o:
-unittest.o: global.o energy_grid.o mpi_routines.o
+#===============================================================================
+
+include DEPENDENCIES
diff --git a/src/OBJECTS b/src/OBJECTS
new file mode 100644
index 0000000000..d9af5f5577
--- /dev/null
+++ b/src/OBJECTS
@@ -0,0 +1,28 @@
+objects = \
+bank_header.o \
+cross_section_header.o \
+cross_section.o \
+datatypes.o \
+datatypes_header.o \
+endf.o \
+energy_grid.o \
+error.o \
+fileio.o \
+geometry.o \
+geometry_header.o \
+global.o \
+logging.o \
+main.o \
+material_header.o \
+mcnp_random.o \
+mpi_routines.o \
+output.o \
+particle_header.o \
+physics.o \
+score.o \
+search.o \
+source.o \
+source_header.o \
+string.o \
+tally_header.o \
+timing.o
diff --git a/src/ace.f90 b/src/ace.f90
deleted file mode 100644
index d30a63ad6f..0000000000
--- a/src/ace.f90
+++ /dev/null
@@ -1,1221 +0,0 @@
-module ace
-
-  use global
-  use error,           only: fatal_error
-  use output,          only: message
-  use string,          only: lower_case
-  use fileio,          only: read_line, read_data, skip_lines
-  use string,          only: split_string, str_to_real
-  use data_structures, only: dict_create, dict_add_key, dict_has_key, &
-       &                     dict_get_key, dict_delete
-  use endf,            only: reaction_name
-
-  integer :: NXS(16)             ! Descriptors for ACE XSS tables
-  integer :: JXS(32)             ! Pointers into ACE XSS tables
-  real(8), allocatable :: XSS(:) ! Cross section data
-  integer :: XSS_index           ! current index in XSS data
-
-  private :: NXS
-  private :: JXS
-  private :: XSS
-
-contains
-
-!===============================================================================
-! READ_XS reads all the cross sections for the problem and stores them in
-! xs_continuous and xs_thermal arrays
-!===============================================================================
-
-  subroutine read_xs()
-
-    integer        :: i                ! index in materials array
-    integer        :: j                ! index over isotopes in material
-    integer        :: index            ! index in xsdatas array
-    integer        :: n                ! length of key
-    integer        :: index_continuous ! index in xs_continuous
-    integer        :: index_thermal    ! index in xs_thermal
-    character(10)  :: key              ! name of isotope, e.g. 92235.03c
-    character(max_line_len) :: msg              ! output/error message
-    type(Material),      pointer :: mat => null()
-    type(xsData),        pointer :: iso => null()
-    type(AceContinuous), pointer :: ace_cont => null()
-    type(AceThermal),    pointer :: ace_thermal => null()
-    type(DictionaryCI),  pointer :: temp_dict => null()
-
-    call dict_create(ace_dict)
-
-    ! determine how many continuous-energy tables and how many S(a,b) thermal
-    ! scattering tables there are
-    index_continuous = 0
-    index_thermal = 0
-    do i = 1, n_materials
-       mat => materials(i)
-       do j = 1, mat%n_isotopes
-          index = mat%isotopes(j)
-          key = xsdatas(index)%id
-          n = len_trim(key)
-          call lower_case(key)
-          select case (key(n:n))
-          case ('c')
-             if (.not. dict_has_key(ace_dict, key)) then
-                index_continuous = index_continuous + 1
-                call dict_add_key(ace_dict, key, index_continuous)
-                mat%table(j) = index_continuous
-             else
-                mat%table(j) = dict_get_key(ace_dict, key)
-             end if
-          case ('t')
-             ! Need same logic as for continuous tables
-             n_thermal = n_thermal + 1
-          case default
-             msg = "Unknown cross section table type: " // key
-             call fatal_error(msg)
-          end select
-       end do
-    end do
-
-    n_continuous = index_continuous
-    n_thermal = index_thermal
-
-    ! allocate arrays for ACE table storage
-    allocate(xs_continuous(n_continuous))
-    allocate(xs_thermal(n_thermal))
-
-    ! loop over all nuclides in xsdata
-    call dict_create(temp_dict)
-
-    index_continuous = 0
-    index_thermal = 0
-    do i = 1, n_materials
-       mat => materials(i)
-       do j = 1, mat%n_isotopes
-          index = mat%isotopes(j)
-          key = xsdatas(index)%id
-          n = len_trim(key)
-          call lower_case(key)
-          select case (key(n:n))
-          case ('c')
-             if (.not. dict_has_key(temp_dict, key)) then
-                index_continuous = index_continuous + 1
-                call read_ACE_continuous(index_continuous, index)
-             end if
-          case ('t')
-             n_thermal = n_thermal + 1
-          end select
-       end do
-    end do
-
-    ! delete dictionary
-    call dict_delete(temp_dict)
-       
-  end subroutine read_xs
-
-!===============================================================================
-! READ_ACE_CONTINUOUS reads in a single ACE continuous-energy neutron cross
-! section table. This routine reads the header data and then calls appropriate
-! subroutines to parse the actual data.
-!===============================================================================
-
-  subroutine read_ACE_continuous(index_table, index)
-
-    integer, intent(in) :: index_table ! index in xs_continuous array
-    integer, intent(in) :: index       ! index in xsdatas array
-
-    integer                 :: in = 7           ! unit to read from
-    integer                 :: ioError          ! error status for file access
-    integer                 :: words_per_line   ! number of words per line (data)
-    integer                 :: lines            ! number of lines (data
-    integer                 :: n                ! number of data values
-    real(8)                 :: kT               ! ACE table temperature
-    logical                 :: file_exists      ! does ACE library exist?
-    logical                 :: found_xs         ! did we find table in library?
-    character(7)            :: readable         ! is ACE library readable?
-    character(max_line_len) :: msg              ! output/error message
-    character(max_line_len) :: line             ! single line to read
-    character(max_word_len) :: words(max_words) ! words on a line
-    character(max_word_len) :: filename         ! name of ACE library file
-    character(10)           :: tablename        ! name of cross section table
-    type(AceContinuous), pointer :: table => null()
-
-    filename = xsdatas(index)%path
-    tablename = xsdatas(index)%id
-
-    table => xs_continuous(index_table)
-
-    ! Check if input file exists and is readable
-    inquire(FILE=filename, EXIST=file_exists, READ=readable)
-    if (.not. file_exists) then
-       msg = "ACE library '" // trim(filename) // "' does not exist!"
-       call fatal_error(msg)
-    elseif (readable(1:3) == 'NO') then
-       msg = "ACE library '" // trim(filename) // "' is not readable! &
-            &Change file permissions with chmod command."
-       call fatal_error(msg)
-    end if
-
-    ! display message
-    msg = "Loading ACE cross section table: " // tablename
-    call message(msg, 6)
-
-    ! open file
-    open(file=filename, unit=in, status='old', & 
-         & action='read', iostat=ioError)
-    if (ioError /= 0) then
-       msg = "Error while opening file: " // filename
-       call fatal_error(msg)
-    end if
-
-    found_xs = .false.
-    do while (.not. found_xs)
-       call read_line(in, line, ioError)
-       if (ioError < 0) then
-          msg = "Could not find ACE table " // tablename // "."
-          call fatal_error(msg)
-       end if
-       call split_string(line, words, n)
-       if (trim(words(1)) == trim(tablename)) then
-          found_xs = .true.
-          table%name = words(1)
-          table%awr = str_to_real(words(2))
-          kT = str_to_real(words(3))
-          table%temp = kT / K_BOLTZMANN
-       end if
-       
-       ! Skip 5 lines
-       call skip_lines(in, 5, ioError)
-
-       ! Read NXS data
-       lines = 2
-       words_per_line = 8
-       call read_data(in, NXS, 16, lines, words_per_line)
-
-       ! Read JXS data
-       lines = 4
-       call read_data(in, JXS, 32, lines, words_per_line)
-
-       ! Calculate how many data points and lines in the XSS array
-       n = NXS(1)
-       lines = (n + 3)/4
-
-       if (found_xs) then
-          ! allocate storage for XSS array
-          allocate(XSS(n))
-          
-          ! Read XSS
-          words_per_line = 4
-          call read_data(in, XSS, n, lines, words_per_line)
-       else
-          call skip_lines(in, lines, ioError)
-       end if
-
-    end do
-
-    call read_esz(table)
-    call read_nu_data(table)
-    call read_reactions(table)
-    call read_angular_dist(table)
-    call read_energy_dist(table)
-    call read_unr_res(table)
-
-    ! Free memory from XSS array
-    if(allocated(XSS)) deallocate(XSS)
-    if(associated(table)) nullify(table)
-
-    close(unit=in)
-
-  end subroutine read_ACE_continuous
-
-!===============================================================================
-! READ_ESZ - reads through the ESZ block. This block contains the energy grid,
-! total xs, absorption xs, elastic scattering xs, and heating numbers.
-!===============================================================================
-
-  subroutine read_esz(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: NE ! number of energy points for total and elastic cross sections
-
-    ! determine number of energy points
-    NE = NXS(3)
-    table%n_grid = NE
-
-    ! allocate storage for arrays
-    allocate(table%energy(NE))
-    allocate(table%sigma_t(NE))
-    allocate(table%sigma_a(NE))
-    allocate(table%sigma_el(NE))
-    allocate(table%heating(NE))
-
-    ! read data from XSS -- right now the total, absorption and elastic
-    ! scattering are read in to these special arrays, but in reality, it should
-    ! be necessary to only store elastic scattering and possibly total
-    ! cross-section for total material xs generation.
-    XSS_index = 1
-    table%energy = get_real(NE)
-    table%sigma_t = get_real(NE)
-    table%sigma_a = get_real(NE)
-    table%sigma_el = get_real(NE)
-    table%heating = get_real(NE)
-    
-  end subroutine read_esz
-
-!===============================================================================
-! READ_NU_DATA reads data given on the number of neutrons emitted from fission
-! as a function of the incoming energy of a neutron. This data may be broken
-! down into prompt and delayed neutrons emitted as well.
-!===============================================================================
-
-  subroutine read_nu_data(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: JXS2   ! location for fission nu data
-    integer :: JXS24  ! location for delayed neutron data
-    integer :: KNU    ! location for nu data
-    integer :: LNU    ! type of nu data (polynomial or tabular)
-    integer :: NC     ! number of polynomial coefficients
-    integer :: NR     ! number of interpolation regions
-    integer :: NE     ! number of energies 
-    integer :: NPCR   ! number of delayed neutron precursor groups
-    integer :: LED    ! location of energy distribution locators
-    integer :: LDIS   ! location of all energy distributions
-    integer :: LOCC   ! location of energy distributions for given MT
-    integer :: LNW    ! location of next energy distribution if multiple
-    integer :: LAW    ! secondary energy distribution law   
-    integer :: IDAT   ! location of first energy distribution for given MT
-    integer :: loc    ! locator
-    integer :: length             ! length of data to allocate
-    integer :: length_interp_data ! length of interpolation data
-
-    JXS2  = JXS(2)
-    JXS24 = JXS(24)
-
-    if (JXS2 == 0) then
-       ! =======================================================================
-       ! NO PROMPT/TOTAL NU DATA
-       table % nu_t_type = NU_NONE
-       table % nu_p_type = NU_NONE
-
-    elseif (XSS(JXS2) > 0) then
-       ! =======================================================================
-       ! PROMPT OR TOTAL NU DATA
-       KNU = JXS2
-       LNU = int(XSS(KNU))
-       if (LNU == 1) then
-          ! Polynomial data
-          table % nu_t_type = NU_POLYNOMIAL
-          table % nu_p_type = NU_NONE
-
-          ! allocate determine how many coefficients for polynomial
-          NC = int(XSS(KNU+1))
-          length = NC + 1
-       elseif (LNU == 2) then
-          ! Tabular data
-          table % nu_t_type = NU_TABULAR
-          table % nu_p_type = NU_NONE
-
-          ! determine number of interpolation regions and number of energies
-          NR = int(XSS(KNU+1))
-          NE = int(XSS(KNU+2+2*NR))
-          length = 2 + 2*NR + 2*NE
-       end if
-
-       ! allocate space for nu data storage
-       allocate(table % nu_t_data(length))
-
-       ! read data -- for polynomial, this is the number of coefficients and the
-       ! coefficients themselves, and for tabular, this is interpolation data
-       ! and tabular E/nu
-       XSS_index = KNU + 1
-       table % nu_t_data = get_real(length)
-
-    elseif (XSS(JXS2) < 0) then
-       ! =======================================================================
-       ! PROMPT AND TOTAL NU DATA -- read prompt data first
-       KNU = JXS2 + 1
-       LNU = XSS(KNU)
-       if (LNU == 1) then
-          ! Polynomial data
-          table % nu_p_type = NU_POLYNOMIAL
-
-          ! allocate determine how many coefficients for polynomial
-          NC = XSS(KNU+1)
-          length = NC + 1
-       elseif (LNU == 2) then
-          ! Tabular data
-          table % nu_p_type = NU_TABULAR
-
-          ! determine number of interpolation regions and number of energies
-          NR = XSS(KNU+1)
-          NE = XSS(KNU+2+2*NR)
-          length = 2 + 2*NR + 2*NE
-       end if
-
-       ! allocate space for nu data storage
-       allocate(table % nu_p_data(length))
-
-       ! read data
-       XSS_index = KNU + 1
-       table % nu_p_data = get_real(length)
-
-       ! Now read total nu data
-       KNU = JXS2 + abs(XSS(JXS2)) + 1
-       LNU = XSS(KNU)
-       if (LNU == 1) then
-          ! Polynomial data
-          table % nu_t_type = NU_POLYNOMIAL
-
-          ! allocate determine how many coefficients for polynomial
-          NC = int(XSS(KNU+1))
-          length = NC + 1
-       elseif (LNU == 2) then
-          ! Tabular data
-          table % nu_t_type = NU_TABULAR
-
-          ! determine number of interpolation regions and number of energies
-          NR = int(XSS(KNU+1))
-          NE = int(XSS(KNU+2+2*NR))
-          length = 2 + 2*NR + 2*NE
-       end if
-
-       ! allocate space for nu data storage
-       allocate(table % nu_t_data(length))
-
-       ! read data
-       XSS_index = KNU + 1
-       table % nu_t_data = get_real(length)
-    end if
-
-    if (JXS24 > 0) then
-       ! =======================================================================
-       ! DELAYED NU DATA
-
-       table % nu_d_type = NU_TABULAR
-       KNU = JXS24
-
-       ! determine size of tabular delayed nu data
-       NR = int(XSS(KNU+1))
-       NE = int(XSS(KNU+2+2*NR))
-       length = 2 + 2*NR + 2*NE
-
-       ! allocate space for delayed nu data
-       allocate(table % nu_d_data(length))
-       
-       ! read delayed nu data
-       XSS_index = KNU + 1
-       table % nu_d_data = get_real(length)
- 
-       ! =======================================================================
-       ! DELAYED NEUTRON ENERGY DISTRIBUTION
-
-       ! Allocate space for secondary energy distribution
-       NPCR = NXS(8)
-       table % n_precursor = NPCR
-       allocate(table % nu_d_edist(NPCR))
-
-       LED  = JXS(26)
-       LDIS = JXS(27)
-
-       ! Loop over all delayed neutron precursor groups
-       do i = 1, NPCR
-          ! find location of energy distribution data
-          LOCC = XSS(LED + i - 1)
-
-          LNW  = XSS(LDIS + LOCC - 1)
-          LAW  = XSS(LDIS + LOCC)
-          IDAT = XSS(LDIS + LOCC + 1)
-          NR   = XSS(LDIS + LOCC + 2)
-          table % nu_d_edist(i) % law = LAW
-
-          ! allocate space for ENDF interpolation parameters
-          if (NR > 0) then
-             allocate(table % nu_d_edist(i) % nbt(NR))
-             allocate(table % nu_d_edist(i) % int(NR))
-          end if
-
-          ! read ENDF interpolation parameters
-          XSS_index = LDIS + LOCC + 3
-          table % nu_d_edist(i) % nbt = get_real(NR)
-          table % nu_d_edist(i) % int = get_real(NR)
-
-          ! allocate space for law validity data
-          NE = XSS(LDIS + LOCC + 3 + 2*NR)
-          allocate(table % nu_d_edist(i) % energy(NE))
-          allocate(table % nu_d_edist(i) % pvalid(NE))
-
-          length_interp_data = 5 + 2*(NR + NE)
-
-          ! read law validity data
-          XSS_index = LDIS + LOCC + 4 + 2*NR
-          table % nu_d_edist(i) % energy = get_real(NE)
-          table % nu_d_edist(i) % pvalid = get_real(NE)
-
-          ! Set index to beginning of IDAT array
-          loc = LDIS + IDAT - 2
-
-          ! determine length of energy distribution
-          length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
-
-          ! allocate secondary energy distribution array
-          allocate(table % nu_d_edist(i) % data(length))
-
-          ! read secondary energy distribution
-          XSS_index = loc + 1
-          table % nu_d_edist(i) % data = get_real(length)
-       end do
-
-       ! =======================================================================
-       ! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS
-
-       ! determine length of all precursor constants/yields/interp data
-       length = 0
-       loc = JXS(25)
-       do i = 1, NPCR
-          NR = XSS(loc + length + 1)
-          NE = XSS(loc + length + 2 + 2*NR)
-          length = length + 3 + 2*NR + 2*NE
-       end do
-
-       ! allocate space for precusor data
-       allocate(table % nu_d_precursor_data(length))
-
-       ! read delayed neutron precursor data
-       XSS_index = loc
-       table % nu_d_precursor_data = get_real(length)
-
-    else
-       table % nu_d_type = NU_NONE
-    end if
-
-  end subroutine read_nu_data
-
-!===============================================================================
-! READ_REACTIONS - Get the list of reaction MTs for this cross-section
-! table. The MT values are somewhat arbitrary. Also read in Q-values, neutron
-! multiplicities, and cross-sections.
-!===============================================================================
-
-  subroutine read_reactions(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: LMT   ! index of MT list in XSS
-    integer :: NMT   ! Number of reactions
-    integer :: JXS4  ! index of Q values in XSS
-    integer :: JXS5  ! index of neutron multiplicities in XSS
-    integer :: JXS7  ! index of reactions cross-sections in XSS
-    integer :: LXS   ! location of cross-section locators
-    integer :: LOCA  ! location of cross-section for given MT
-    integer :: NE    ! number of energies for reaction
-    type(AceReaction), pointer :: rxn => null()
-    
-    LMT  = JXS(3)
-    JXS4 = JXS(4)
-    JXS5 = JXS(5)
-    LXS  = JXS(6)
-    JXS7 = JXS(7)
-    NMT  = NXS(4)
-
-    ! allocate array of reactions. Add one since we need to include an elastic
-    ! scattering channel
-    table%n_reaction = NMT + 1
-    allocate(table%reactions(NMT+1))
-
-    ! Store elastic scattering cross-section on reaction one
-    rxn => table%reactions(1)
-    rxn%MT      = 2
-    rxn%Q_value = ZERO
-    rxn%TY      = 1
-    rxn%IE      = 1
-    allocate(rxn%sigma(table%n_grid))
-    rxn%sigma = table%sigma_el
-
-    do i = 1, NMT
-       rxn => table%reactions(i+1)
-
-       ! read MT number, Q-value, and neutrons produced
-       rxn % MT      = XSS(LMT+i-1)
-       rxn % Q_value = XSS(JXS4+i-1)
-       rxn % TY      = XSS(JXS5+i-1)
-
-       ! read cross section values
-       LOCA = XSS(LXS+i-1)
-       rxn % IE = XSS(JXS7 + LOCA - 1)
-       NE = XSS(JXS7 + LOCA)
-       allocate(rxn%sigma(NE))
-       XSS_index = JXS7 + LOCA + 1
-       rxn % sigma = get_real(NE)
-
-       ! set defaults
-       rxn % has_angle_dist  = .false.
-       rxn % has_energy_dist = .false.
-    end do
-
-  end subroutine read_reactions
-
-!===============================================================================
-! READ_ANGULAR_DIST parses the angular distribution for each reaction with
-! secondary neutrons
-!===============================================================================
-
-  subroutine read_angular_dist(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: JXS8   ! location of angular distribution locators
-    integer :: JXS9   ! location of angular distributions
-    integer :: LOCB   ! location of angular distribution for given MT
-    integer :: NE     ! number of incoming energies
-    integer :: NP     ! number of points for cosine distribution
-    integer :: LC     ! locator
-    integer :: i      ! index in reactions array
-    integer :: j      ! index over incoming energies
-    integer :: length ! length of data array to allocate
-    type(AceReaction), pointer :: rxn => null()
-
-    JXS8 = JXS(8)
-    JXS9 = JXS(9)
-
-    ! loop over all reactions with secondary neutrons -- NXS(5) does not include
-    ! elastic scattering
-    do i = 1, NXS(5) + 1
-       rxn => table%reactions(i)
-
-       ! find location of angular distribution
-       LOCB = XSS(JXS8 + i - 1)
-       if (LOCB == -1) then
-          ! Angular distribution data are specified through LAWi = 44 in the DLW
-          ! block
-          cycle
-       elseif (LOCB == 0) then
-          ! No angular distribution data are given for this reaction, isotropic
-          ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0)
-          cycle
-       end if
-       rxn % has_angle_dist = .true.
-
-       ! allocate space for incoming energies and locations
-       NE = XSS(JXS9 + LOCB - 1)
-       rxn % adist % n_energy = NE
-       allocate(rxn % adist % energy(NE))
-       allocate(rxn % adist % type(NE))
-       allocate(rxn % adist % location(NE))
-
-       ! read incoming energy grid and location of tables
-       XSS_index = JXS9 + LOCB
-       rxn % adist % energy   = get_real(NE)
-       rxn % adist % location = get_int(NE)
-
-       ! determine dize of data block
-       length = 0
-       do j = 1, NE
-          LC = rxn % adist % location(j)
-          if (LC == 0) then
-             ! isotropic
-             rxn % adist % type(j) = ANGLE_ISOTROPIC
-          elseif (LC > 0) then
-             ! 32 equiprobable bins
-             rxn % adist % type(j) = ANGLE_32_EQUI
-             length = length + 33
-          elseif (LC < 0) then
-             ! tabular distribution
-             rxn % adist % type(j) = ANGLE_TABULAR
-             NP = XSS(JXS9 + abs(LC))
-             length = length + 2 + 3*NP
-          end if
-       end do
-       
-       ! allocate angular distribution data and read
-       allocate(rxn % adist % data(length))
-
-       ! read angular distribution -- currently this does not actually parse the
-       ! angular distribution tables for each incoming energy, that must be done
-       ! on-the-fly
-       LC = rxn % adist % location(1)
-       XSS_index = JXS9 + abs(LC) - 1
-       rxn % adist % data = get_real(length)
-
-       ! change location pointers since they are currently relative to JXS(9)
-       LC = abs(rxn % adist % location(1))
-       rxn % adist % location = abs(rxn % adist % location) - LC
-
-    end do
-
-  end subroutine read_angular_dist
-
-!===============================================================================
-! READ_ENERGY_DIST parses the secondary energy distribution for each reaction
-! with seconary neutrons (except elastic scattering)
-!===============================================================================
-
-  subroutine read_energy_dist(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: LED   ! location of energy distribution locators
-    integer :: LDIS  ! location of all energy distributions
-    integer :: LOCC  ! location of energy distributions for given MT
-    integer :: LNW   ! location of next energy distribution if multiple
-    integer :: LAW   ! secondary energy distribution law   
-    integer :: NR    ! number of interpolation regions
-    integer :: NE    ! number of incoming energies
-    integer :: IDAT  ! location of first energy distribution for given MT
-    integer :: loc   ! locator
-    integer :: length             ! length of data to allocate
-    integer :: length_interp_data ! length of interpolation data
-    integer :: i, j, k, l         ! indices
-    type(AceReaction), pointer :: rxn => null()
-
-    LED  = JXS(10)
-    LDIS = JXS(11)
-
-    ! Loop over all reactions 
-    do i = 1, NXS(5)
-       rxn => table % reactions(i+1) ! skip over elastic scattering
-       rxn % has_energy_dist = .true.
-
-       ! find location of energy distribution data
-       LOCC = XSS(LED + i - 1)
-       
-       LNW  = XSS(LDIS + LOCC - 1)
-       LAW  = XSS(LDIS + LOCC)
-       IDAT = XSS(LDIS + LOCC + 1)
-       NR   = XSS(LDIS + LOCC + 2)
-       rxn % edist % law = LAW
-       rxn % edist % n_interp = NR
-
-       ! allocate space for ENDF interpolation parameters
-       if (NR > 0) then
-          allocate(rxn % edist % nbt(NR))
-          allocate(rxn % edist % int(NR))
-       end if
-       
-       ! read ENDF interpolation parameters
-       XSS_index = LDIS + LOCC + 3
-       rxn % edist % nbt = get_real(NR)
-       rxn % edist % int = get_real(NR)
-
-       ! allocate space for law validity data
-       NE = XSS(LDIS + LOCC + 3 + 2*NR)
-       allocate(rxn % edist % energy(NE))
-       allocate(rxn % edist % pvalid(NE))
-       
-       length_interp_data = 5 + 2*(NR + NE)
-       
-       ! read law validity data
-       XSS_index = LDIS + LOCC + 4 + 2*NR
-       rxn % edist % energy = get_real(NE)
-       rxn % edist % pvalid = get_real(NE)
-
-       ! Set index to beginning of IDAT array
-       loc = LDIS + IDAT - 2
-
-       ! determine length of energy distribution
-       length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
-
-       ! allocate secondary energy distribution array
-       allocate(rxn % edist % data(length))
-
-       ! read secondary energy distribution
-       XSS_index = loc + 1
-       rxn % edist % data = get_real(length)
-    end do
-
-  end subroutine read_energy_dist
-
-!===============================================================================
-! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy
-! distribution array based on the secondary energy law and location in XSS
-!===============================================================================
-
-  function length_energy_dist(loc, law, LOCC, lid) result(length)
-
-    integer, intent(in) :: loc    ! location in XSS array
-    integer, intent(in) :: law    ! energy distribution law
-    integer, intent(in) :: LOCC   ! location of energy distribution
-    integer, intent(in) :: lid    ! length of interpolation data
-    integer             :: length ! length of energy distribution (LDAT)
-
-    integer :: i,j,k
-    integer :: NR    ! number of interpolation regions
-    integer :: NE    ! number of incoming energies
-    integer :: NP    ! number of points in outgoing energy distribution
-    integer :: NMU   ! number of points in outgoing cosine distribution
-    integer :: NRa   ! number of interpolation regions for Watt 'a'
-    integer :: NEa   ! number of energies for Watt 'a'
-    integer :: NRb   ! number of interpolation regions for Watt 'b'
-    integer :: NEb   ! number of energies for Watt 'b'
-
-    ! initialize length
-    length = 0
-
-    select case (law)
-    case (1)
-       ! Tabular equiprobable energy bins
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       NP = XSS(loc + 3 + 2*NR + NE)
-       length = 3 + 2*NR + NE + 3*NP*NE
-
-    case (2)
-       ! Discrete photon energy
-       length = 2
-
-    case (3)
-       ! Level scattering
-       length = 2
-
-    case (4)
-       ! Continuous tabular distribution
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       length = length + 2 + 2*NR + 2*NE
-       do i = 1,NE
-          ! determine length
-          NP = XSS(loc + length + 2)
-          length = length + 2 + 3*NP
-
-          ! adjust location for this block
-          j = loc + 2 + 2*NR + NE + i
-          XSS(j) = XSS(j) - LOCC - lid
-       end do
-
-    case (5)
-       ! General evaporation spectrum
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       NP = XSS(loc + 3 + 2*NR + 2*NE)
-       length = 3 + 2*NR + 2*NE + NP
-
-    case (7)
-       ! Maxwell fission spectrum
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       length = 3 + 2*NR + 2*NE
-
-    case (9)
-       ! Evaporation spectrum
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       length = 3 + 2*NR + 2*NE
-
-    case (11)
-       ! Watt spectrum
-       NRa = XSS(loc + 1)
-       NEa = XSS(loc + 2 + 2*NRa)
-       NRb = XSS(loc + 3 + 2*(NRa+NEa))
-       NEb = XSS(loc + 4 + 2*(NRa+NEa+NRb))
-       length = 5 + 2*(NRa + NEa + NRb + NEb)
-
-    case (44)
-       ! Kalbach-Mann correlated scattering
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       length = length + 2 + 2*NR + 2*NE
-       do i = 1,NE
-          NP = XSS(loc + length + 2)
-          length = length + 2 + 5*NP
-
-          ! adjust location for this block
-          j = loc + 2 + 2*NR + NE + i
-          XSS(j) = XSS(j) - LOCC - lid
-       end do
-
-    case (61)
-       ! Correlated energy and angle distribution
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       length = length + 2 + 2*NR + 2*NE
-       do i = 1,NE
-          ! outgoing energy distribution
-          NP = XSS(loc + length + 2)
-
-          ! adjust locators for angular distribution
-          do j = 1, NP
-             k = loc + length + 2 + 3*NP + j
-             if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
-          end do
-
-          length = length + 2 + 4*NP
-          do j = 1, NP
-             ! outgoing angle distribution -- NMU here is actually
-             ! referred to as NP in the MCNP documentation
-             NMU = XSS(loc + length + 2)
-             length = length + 2 + 3*NMU
-          end do
-
-          ! adjust locators for energy distribution
-          j = loc + 2 + 2*NR + NE + i
-          XSS(j) = XSS(j) - LOCC - lid
-       end do
-
-    case (66)
-       ! N-body phase space distribution
-       length = 2
-
-    case (67)
-       ! Laboratory energy-angle law
-       NR = XSS(loc + 1)
-       NE = XSS(loc + 2 + 2*NR)
-       NMU = XSS(loc + 4 + 2*NR + 2*NE)
-       length = 4 + 2*(NR + NE + NMU)
-
-    end select
-
-  end function length_energy_dist
-
-!===============================================================================
-! READ_UNR_RES reads in unresolved resonance probability tables if present.
-!===============================================================================
-
-  subroutine read_unr_res(table)
-
-    type(AceContinuous), pointer :: table
-
-    integer :: JXS23 ! location of URR data
-    integer :: loc   ! locator
-    integer :: N     ! # of incident energies
-    integer :: M     ! # of probabilities
-    integer :: i     ! index over incoming energies
-    integer :: j     ! index over values
-    integer :: k     ! index over cumulative probabilities
-
-    ! determine locator for URR data
-    JXS23 = JXS(23)
-
-    ! check if URR data is present
-    if (JXS23 /= 0) then
-       table % urr_present = .true.
-       allocate(table % urr_params(6))
-       loc = JXS23
-    else
-       table % urr_present = .false.
-       return
-    end if
-
-    ! read parameters
-    table % urr_params(1) = XSS(loc)     ! # of incident energies
-    table % urr_params(2) = XSS(loc + 1) ! # of probabilities
-    table % urr_params(3) = XSS(loc + 2) ! interpolation parameter
-    table % urr_params(4) = XSS(loc + 3) ! inelastic competition flag
-    table % urr_params(5) = XSS(loc + 4) ! other absorption flag
-    table % urr_params(6) = XSS(loc + 5) ! factors flag
-
-    ! allocate incident energies and probability tables
-    N = table % urr_params(1)
-    M = table % urr_params(2)
-    allocate(table % urr_energy(N))
-    allocate(table % urr_prob(N,6,M))
-
-    ! read incident energies
-    XSS_index = loc + 6
-    table % urr_energy = get_real(N)
-
-    ! read probability tables
-    do i = 1, N
-       do j = 1, 6
-          table % urr_prob(i,j,1:M) = get_real(M)
-       end do
-    end do
-
-  end subroutine read_unr_res
-
-!===============================================================================
-! READ_ACE_THERMAL reads in a single ACE thermal scattering S(a,b) table. This
-! routine in particular reads the header data and then calls read_thermal_data
-! to parse the actual data blocks.
-!===============================================================================
-
-  subroutine read_ACE_thermal(index_table, index)
-
-    integer, intent(in) :: index_table ! index in xs_thermal array
-    integer, intent(in) :: index       ! index in xsdatas array
-
-    integer                 :: in = 7           ! unit to read from
-    integer                 :: ioError          ! error status for file access
-    integer                 :: words_per_line   ! number of words per line (data)
-    integer                 :: lines            ! number of lines (data
-    integer                 :: n                ! number of data values
-    real(8)                 :: kT               ! ACE table temperature
-    logical                 :: file_exists      ! does ACE library exist?
-    logical                 :: found_xs         ! did we find table in library?
-    character(7)            :: readable         ! is ACE library readable?
-    character(max_line_len) :: msg              ! output/error message
-    character(max_line_len) :: line             ! single line to read
-    character(max_word_len) :: words(max_words) ! words on a line
-    character(max_word_len) :: filename         ! name of ACE library file
-    character(10)           :: tablename        ! name of cross section table
-    type(AceThermal), pointer :: table => null()
-
-    filename = xsdatas(index)%path
-    tablename = xsdatas(index)%id
-
-    table => xs_thermal(index_table)
-
-    ! Check if input file exists and is readable
-    inquire(FILE=filename, EXIST=file_exists, READ=readable)
-    if (.not. file_exists) then
-       msg = "ACE library '" // trim(filename) // "' does not exist!"
-       call fatal_error(msg)
-    elseif (readable(1:3) == 'NO') then
-       msg = "ACE library '" // trim(filename) // "' is not readable! &
-            &Change file permissions with chmod command."
-       call fatal_error(msg)
-    end if
-
-    ! display message
-    msg = "Loading ACE cross section table: " // tablename
-    call message(msg, 6)
-
-    ! open file
-    open(file=filename, unit=in, status='old', & 
-         & action='read', iostat=ioError)
-    if (ioError /= 0) then
-       msg = "Error while opening file: " // filename
-       call fatal_error(msg)
-    end if
-
-    found_xs = .false.
-    do while (.not. found_xs)
-       call read_line(in, line, ioError)
-       if (ioError < 0) then
-          msg = "Could not find ACE table " // tablename // "."
-          call fatal_error(msg)
-       end if
-       call split_string(line, words, n)
-       if (trim(words(1)) == trim(tablename)) then
-          found_xs = .true.
-          table%name = words(1)
-          table%awr = str_to_real(words(2))
-          kT = str_to_real(words(3))
-          table%temp = kT / K_BOLTZMANN
-       end if
-       
-       ! Skip 5 lines
-       call skip_lines(in, 5, ioError)
-
-       ! Read NXS data
-       lines = 2
-       words_per_line = 8
-       call read_data(in, NXS, 16, lines, words_per_line)
-
-       ! Read JXS data
-       lines = 4
-       call read_data(in, JXS, 32, lines, words_per_line)
-
-       ! Calculate how many data points and lines in the XSS array
-       n = NXS(1)
-       lines = (n + 3)/4
-
-       if (found_xs) then
-          ! allocate storage for XSS array
-          allocate(XSS(n))
-          
-          ! Read XSS
-          words_per_line = 4
-          call read_data(in, XSS, n, lines, words_per_line)
-       else
-          call skip_lines(in, lines, ioError)
-       end if
-
-    end do
-
-    call read_thermal_data(table)
-
-    ! Free memory from XSS array
-    if(allocated(XSS)) deallocate(XSS)
-    if(associated(table)) nullify(table)
-
-    close(unit=in)
-
-  end subroutine read_ACE_thermal
-
-!===============================================================================
-! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding
-! secondary energy/angle distributions derived from experimental S(a,b)
-! data. Namely, in MCNP language, this routine reads the ITIE, ITCE, ITXE, and
-! ITCA blocks.
-!===============================================================================
-
-  subroutine read_thermal_data(table)
-
-    type(AceThermal), pointer :: table
-
-    integer :: i      ! index for incoming energies
-    integer :: j      ! index for outgoing energies
-    integer :: k      ! index for outoging angles
-    integer :: loc    ! location in XSS array
-    integer :: NE_in  ! number of incoming energies
-    integer :: NE_out ! number of outgoing energies
-    integer :: NMU    ! number of outgoing angles
-    integer :: JXS4   ! location of elastic energy table
-
-    ! read number of inelastic energies and allocate arrays
-    NE_in = XSS(JXS(1))
-    table % n_inelastic_e_in = NE_in
-    allocate(table % inelastic_e_in(NE_in))
-    allocate(table % inelastic_sigma(NE_in))
-
-    ! read inelastic energies and cross-sections
-    XSS_index = JXS(1) + 1
-    table % inelastic_e_in = get_real(NE_in)
-    table % inelastic_sigma = get_real(NE_in)
-
-    ! allocate space for outgoing energy/angle for inelastic
-    ! scattering
-    NE_out = NXS(4)
-    NMU = NXS(3) + 1
-    allocate(table % inelastic_e_out(NE_out, NE_in))
-    allocate(table % inelastic_mu(NMU, NE_out, NE_in))
-
-    ! read outgoing energy/angle distribution for inelastic scattering
-    loc = JXS(3) - 1
-    do i = 1, NE_in
-       do j = 1, NE_out
-          ! read outgoing energy
-          table % inelastic_e_out(j,i) = XSS(loc + 1)
-
-          ! read outgoing angles for this outgoing energy
-          do k = 1, NMU
-             table % inelastic_mu(k,j,i) = XSS(loc + 1 + k)
-          end do
-
-          ! advance pointer
-          loc = loc + 1 + NMU
-       end do
-    end do
-
-    ! read number of elastic energies and allocate arrays
-    JXS4 = JXS(4)
-    if (JXS4 /= 0) then
-       NE_in = XSS(JXS4)
-       table % n_elastic_e_in = NE_in
-       allocate(table % elastic_e_in(NE_in))
-       allocate(table % elastic_P(NE_in))
-
-       ! determine whether sigma=P or sigma = P/E
-       table % n_elastic_type = NXS(5)
-    else
-       table % n_elastic_e_in = 0
-    end if
-
-    ! allocate space for outgoing energy/angle for elastic scattering
-    NMU = NXS(6) + 1
-    table % n_elastic_mu = NMU
-    if (NMU > 0) then
-       allocate(table % elastic_mu(NMU, NE_in))
-    end if
-
-    ! read equiprobable outgoing cosines for elastic scattering each
-    ! incoming energy
-    if (JXS4 /= 0 .and. NMU /= 0) then
-       loc = JXS(6) - 1
-       do i = 1, NE_in
-          do j = 1, NMU
-             table % elastic_mu(j,i) = XSS(loc + j)
-          end do
-          loc = loc + NMU
-       end do
-    end if
-    
-  end subroutine read_thermal_data
-
-!===============================================================================
-! GET_INT returns an array of integers read from the current position in the XSS
-! array
-!===============================================================================
-
-  function get_int(n_values) result(array)
-
-    integer, intent(in) :: n_values        ! number of values to read
-    integer             :: array(n_values) ! array of values
-
-    array = int(XSS(XSS_index:XSS_index + n_values - 1))
-    XSS_index = XSS_index + n_values
-
-  end function get_int
-
-!===============================================================================
-! GET_REAL returns an array of real(8)s read from the current position in the
-! XSS array
-!===============================================================================
-
-  function get_real(n_values) result(array)
-
-    integer, intent(in) :: n_values        ! number of values to read
-    real(8)             :: array(n_values) ! array of values
-
-    array = XSS(XSS_index:XSS_index + n_values - 1)
-    XSS_index = XSS_index + n_values
-
-  end function get_real
-
-!===============================================================================
-! GET_MACRO_XS
-!===============================================================================
-
-  function get_macro_xs(p, mat, MT) result(xs)
-
-    type(Particle), pointer    :: p
-    type(Material), pointer    :: mat
-    integer,        intent(in) :: MT
-    real(8)                    :: xs
-
-    integer :: i
-    integer :: n_isotopes
-    integer :: IE
-    real(8) :: density
-    real(8) :: density_i
-    real(8) :: sigma_i
-    real(8) :: f
-    type(AceContinuous), pointer :: table => null()
-    type(AceReaction),   pointer :: rxn => null()
-
-    ! initialize xs
-    xs = ZERO
-
-    ! find material atom density
-    density = cMaterial%atom_density
-
-    ! loop over all isotopes in material
-    n_isotopes = cMaterial % n_isotopes
-    do i = 1, n_isotopes
-       table => xs_continuous(cMaterial % table(i))
-
-       ! determine nuclide atom density
-       density_i = cMaterial%atom_percent(i) * density
-
-       ! search nuclide energy grid
-       IE = table%grid_index(p % IE)
-       f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
-       
-       ! handle special case of total cross section
-       if (MT == 1) then
-          xs = xs + density * (ONE-f) * table%sigma_t(IE) + & 
-               & f * (table%sigma_t(IE+1))
-          cycle
-       end if
-
-       ! loop over reactions in isotope
-       do j = 1, table % n_reaction
-          rxn => table % reactions(i)
-
-          ! check for matching MT
-          if (MT /= rxn % MT) cycle
-
-          ! if energy is below threshold for this reaction, skip it
-          if (IE < rxn % IE) cycle
-
-          ! add to cumulative probability
-          sigma_i = (ONE-f) * rxn%sigma(IE-rxn%IE+1) + & 
-               & f * (rxn%sigma(IE-rxn%IE+2))
-       end do
-
-       ! calculate nuclide macroscopic cross-section
-       xs = xs + density_i * sigma_i
-    end do
-
-  end function get_macro_xs
-
-end module ace
diff --git a/src/bank_header.f90 b/src/bank_header.f90
new file mode 100644
index 0000000000..b896382c98
--- /dev/null
+++ b/src/bank_header.f90
@@ -0,0 +1,18 @@
+module bank_header
+
+  implicit none
+
+!===============================================================================
+! BANK is used for storing fission sites in criticality calculations. Since all
+! the state information of a neutron is not needed, this type allows sites to be
+! stored with less memory
+!===============================================================================
+
+  type Bank
+     integer(8) :: uid    ! Unique ID
+     real(8)    :: xyz(3) ! location of bank particle
+     real(8)    :: uvw(3) ! diretional cosines
+     real(8)    :: E      ! energy
+  end type Bank
+
+end module bank_header
diff --git a/src/constants.f90 b/src/constants.f90
new file mode 100644
index 0000000000..fda54f9e41
--- /dev/null
+++ b/src/constants.f90
@@ -0,0 +1,177 @@
+module constants
+
+  implicit none
+
+  ! Physical constants
+  real(8), parameter ::            &
+       PI           = 3.1415926535898_8, & ! pi
+       MASS_NEUTRON = 1.0086649156,      & ! mass of a neutron
+       MASS_PROTON  = 1.00727646677,     & ! mass of a proton
+       AMU          = 1.66053873e-27,    & ! 1 amu in kg
+       N_AVOGADRO   = 0.602214179,       & ! Avogadro's number in 10^24/mol
+       K_BOLTZMANN  = 8.617342e-5,       & ! Boltzmann constant in eV/K
+       INFINITY     = huge(0.0_8),       & ! positive infinity
+       ZERO         = 0.0_8,             &
+       ONE          = 1.0_8,             &
+       TWO          = 2.0_8
+
+  ! Boundary conditions
+  integer, parameter ::    &
+       & BC_TRANSMIT = 0,  & ! Transmission boundary condition (default)
+       & BC_VACUUM   = 1,  & ! Vacuum boundary condition
+       & BC_REFLECT  = 2,  & ! Reflecting boundary condition
+       & BC_PERIODIC = 3     ! Periodic boundary condition
+
+  ! Logical operators for cell definitions
+  integer, parameter ::                &
+       & OP_LEFT_PAREN  = huge(0),     & ! Left parentheses
+       & OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
+       & OP_UNION       = huge(0) - 2, & ! Union operator
+       & OP_DIFFERENCE  = huge(0) - 3    ! Difference operator
+
+  ! Cell types
+  integer, parameter ::    &
+       & CELL_NORMAL  = 1, & ! Cell with a specified material
+       & CELL_FILL    = 2, & ! Cell filled by a separate universe
+       & CELL_LATTICE = 3    ! Cell filled with a lattice
+
+  ! Lattice types
+  integer, parameter :: &
+       & LATTICE_RECT = 1, &
+       & LATTICE_HEX  = 2
+
+  ! Surface types
+  integer, parameter ::    &
+       & SURF_PX     =  1, & ! Plane parallel to x-plane 
+       & SURF_PY     =  2, & ! Plane parallel to y-plane 
+       & SURF_PZ     =  3, & ! Plane parallel to z-plane 
+       & SURF_PLANE  =  4, & ! Arbitrary plane
+       & SURF_CYL_X  =  5, & ! Cylinder along x-axis
+       & SURF_CYL_Y  =  6, & ! Cylinder along y-axis
+       & SURF_CYL_Z  =  7, & ! Cylinder along z-axis
+       & SURF_SPHERE =  8, & ! Sphere
+       & SURF_BOX_X  =  9, & ! Box extending infinitely in x-direction
+       & SURF_BOX_Y  = 10, & ! Box extending infinitely in y-direction
+       & SURF_BOX_Z  = 11, & ! Box extending infinitely in z-direction
+       & SURF_BOX    = 12, & ! Rectangular prism
+       & SURF_GQ     = 13    ! General quadratic surface
+
+  ! Surface senses
+  integer, parameter ::      &
+       & SENSE_POSITIVE = 1, &
+       & SENSE_NEGATIVE = -1
+
+  ! Codes for read errors -- better hope these numbers are never used in an
+  ! input file!
+  integer, parameter :: ERROR_INT  = -huge(0)
+  real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8
+
+  ! Source types
+  integer, parameter ::   &
+       SRC_BOX     = 1, & ! Source in a rectangular prism
+       SRC_CELL    = 2, & ! Source in a cell
+       SRC_SURFACE = 3    ! Source on a surface
+
+  integer, parameter ::        &
+       PROB_SOURCE      = 1, & ! External source problem
+       PROB_CRITICALITY = 2    ! Criticality problem
+
+  ! Interpolation flag
+  integer, parameter ::     &
+       HISTOGRAM     = 1, & ! y is constant in x
+       LINEAR_LINEAR = 2, & ! y is linear in x
+       LINEAR_LOG    = 3, & ! y is linear in ln(x)
+       LOG_LINEAR    = 4, & ! ln(y) is linear in x
+       LOG_LOG       = 5    ! ln(y) is linear in ln(x)
+
+  ! Particle type
+  integer, parameter :: &
+       NEUTRON  = 1, &
+       PHOTON   = 2, &
+       ELECTRON = 3
+
+  ! Angular distribution type
+  integer, parameter :: & 
+       & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
+       & ANGLE_32_EQUI   = 2, & ! 32 equiprobable bins
+       & ANGLE_TABULAR   = 3    ! Tabular angular distribution
+
+  ! Reaction types
+  integer, parameter :: &
+       TOTAL_XS  = 1, &
+       ELASTIC   = 2, &
+       N_LEVEL   = 4, &
+       MISC      = 5, &
+       N_2ND     = 11, &
+       N_2N      = 16, &
+       N_3N      = 17, &
+       FISSION   = 18, &
+       N_F       = 19, &
+       N_NF      = 20, &
+       N_2NF     = 21, &
+       N_NA      = 22, &
+       N_N3A     = 23, &
+       N_2NA     = 24, &
+       N_3NA     = 25, &
+       N_NP      = 28, &
+       N_N2A     = 29, &
+       N_2N2A    = 30, &
+       N_ND      = 32, &
+       N_NT      = 33, &
+       N_N3HE    = 34, &
+       N_ND2A    = 35, &
+       N_NT2A    = 36, &
+       N_4N      = 37, &
+       N_3NF     = 38, &
+       N_2NP     = 41, &
+       N_3NP     = 42, &
+       N_N2P     = 44, &
+       N_NPA     = 45, &
+       N_N1      = 51, &
+       N_N40     = 90, &
+       N_NC      = 91, &
+       N_GAMMA   = 102, &
+       N_P       = 103, &
+       N_D       = 104, &
+       N_T       = 105, &
+       N_3HE     = 106, &
+       N_A       = 107, &
+       N_2A      = 108, &
+       N_3A      = 109, &
+       N_2P      = 111, &
+       N_PA      = 112, &
+       N_T2A     = 113, &
+       N_D2A     = 114, &
+       N_PD      = 115, &
+       N_PT      = 116, &
+       N_DA      = 117
+
+  ! Tally distinguishers
+  integer, parameter :: &
+       TALLY_FLUX = -1, & ! tally the flux only
+       TALLY_ALL  =  0, & ! tally all reactions
+       TALLY_BINS =  1, & ! tally individual (reactions, cells)
+       TALLY_SUM  =  2    ! tally sum of specified (reactions, cells)
+
+  ! Fission neutron emission (nu) type
+  integer, parameter ::     &
+       NU_NONE       = 0, & ! No nu values (non-fissionable)
+       NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
+       NU_TABULAR    = 2    ! Nu values given by tabular distribution
+
+  ! Versioning numbers
+  integer, parameter :: VERSION_MAJOR = 0
+  integer, parameter :: VERSION_MINOR = 2
+  integer, parameter :: VERSION_RELEASE = 2
+
+  ! Maximum number of words in a single line, length of line, and length of
+  ! single word
+  integer, parameter :: MAX_WORDS    = 500
+  integer, parameter :: MAX_LINE_LEN = 250
+  integer, parameter :: MAX_WORD_LEN = 150
+
+  ! Unit numbers
+  integer, parameter :: UNIT_LOG = 9 ! unit # for writing log file
+
+end module constants
+
diff --git a/src/cross_section.f90 b/src/cross_section.f90
index 6706975194..dcbe06bd22 100644
--- a/src/cross_section.f90
+++ b/src/cross_section.f90
@@ -1,16 +1,1231 @@
 module cross_section
 
+  use constants
+  use cross_section_header, only: Nuclide, Reaction, SAB_Table, xsData
+  use datatypes,            only: dict_create, dict_add_key, dict_get_key, &
+                                  dict_has_key, dict_delete
+  use datatypes_header,     only: DictionaryCI
+  use endf,                 only: reaction_name
+  use error,                only: fatal_error
+  use fileio,               only: read_line, read_data, skip_lines
   use global
-  use string,          only: split_string, str_to_int, str_to_real
-  use data_structures, only: dict_create, dict_add_key, dict_get_key
-  use output,          only: message
-  use error,           only: fatal_error
-  use types,           only: xsData
+  use material_header,      only: Material
+  use output,               only: message
+  use string,               only: split_string, str_to_int, str_to_real, &
+                                  lower_case
 
   implicit none
 
+  integer :: NXS(16)             ! Descriptors for ACE XSS tables
+  integer :: JXS(32)             ! Pointers into ACE XSS tables
+  real(8), allocatable :: XSS(:) ! Cross section data
+  integer :: XSS_index           ! current index in XSS data
+
+  private :: NXS
+  private :: JXS
+  private :: XSS
+
 contains
 
+!===============================================================================
+! READ_XS reads all the cross sections for the problem and stores them in
+! nuclides and sab_tables arrays
+!===============================================================================
+
+  subroutine read_xs()
+
+    integer        :: i                ! index in materials array
+    integer        :: j                ! index over isotopes in material
+    integer        :: index            ! index in xsdatas array
+    integer        :: n                ! length of key
+    integer        :: index_nuclides   ! index in nuclides
+    integer        :: index_sab        ! index in sab_tables
+    character(10)  :: key              ! name of isotope, e.g. 92235.03c
+    character(MAX_LINE_LEN) :: msg     ! output/error message
+    type(Material),     pointer :: mat => null()
+    type(xsData),       pointer :: iso => null()
+    type(Nuclide),      pointer :: ace_cont => null()
+    type(SAB_Table),    pointer :: ace_thermal => null()
+    type(DictionaryCI), pointer :: temp_dict => null()
+
+    call dict_create(nuclide_dict)
+
+    ! determine how many continuous-energy tables and how many S(a,b) thermal
+    ! scattering tables there are
+    index_nuclides = 0
+    index_sab = 0
+    do i = 1, n_materials
+       mat => materials(i)
+       do j = 1, mat%n_isotopes
+          index = mat%isotopes(j)
+          key = xsdatas(index)%id
+          n = len_trim(key)
+          call lower_case(key)
+          select case (key(n:n))
+          case ('c')
+             if (.not. dict_has_key(nuclide_dict, key)) then
+                index_nuclides = index_nuclides + 1
+                call dict_add_key(nuclide_dict, key, index_nuclides)
+                mat%table(j) = index_nuclides
+             else
+                mat%table(j) = dict_get_key(nuclide_dict, key)
+             end if
+          case ('t')
+             ! Need same logic as for continuous tables
+             n_sab_tables = n_sab_tables + 1
+          case default
+             msg = "Unknown cross section table type: " // key
+             call fatal_error(msg)
+          end select
+       end do
+    end do
+
+    n_nuclides   = index_nuclides
+    n_sab_tables = index_sab
+
+    ! allocate arrays for ACE table storage
+    allocate(nuclides(n_nuclides))
+    allocate(sab_tables(n_sab_tables))
+
+    ! loop over all nuclides in xsdata
+    call dict_create(temp_dict)
+
+    index_nuclides = 0
+    index_sab = 0
+    do i = 1, n_materials
+       mat => materials(i)
+       do j = 1, mat%n_isotopes
+          index = mat%isotopes(j)
+          key = xsdatas(index)%id
+          n = len_trim(key)
+          call lower_case(key)
+          select case (key(n:n))
+          case ('c')
+             if (.not. dict_has_key(temp_dict, key)) then
+                index_nuclides = index_nuclides + 1
+                call read_ACE_continuous(index_nuclides, index)
+             end if
+          case ('t')
+             n_sab_tables = n_sab_tables + 1
+          end select
+       end do
+    end do
+
+    ! delete dictionary
+    call dict_delete(temp_dict)
+       
+  end subroutine read_xs
+
+!===============================================================================
+! READ_ACE_CONTINUOUS reads in a single ACE continuous-energy neutron cross
+! section table. This routine reads the header data and then calls appropriate
+! subroutines to parse the actual data.
+!===============================================================================
+
+  subroutine read_ACE_continuous(index_table, index)
+
+    integer, intent(in) :: index_table ! index in nuclides array
+    integer, intent(in) :: index       ! index in xsdatas array
+
+    integer                 :: in = 7           ! unit to read from
+    integer                 :: ioError          ! error status for file access
+    integer                 :: words_per_line   ! number of words per line (data)
+    integer                 :: lines            ! number of lines (data
+    integer                 :: n                ! number of data values
+    real(8)                 :: kT               ! ACE table temperature
+    logical                 :: file_exists      ! does ACE library exist?
+    logical                 :: found_xs         ! did we find table in library?
+    character(7)            :: readable         ! is ACE library readable?
+    character(MAX_LINE_LEN) :: msg              ! output/error message
+    character(MAX_LINE_LEN) :: line             ! single line to read
+    character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line
+    character(MAX_WORD_LEN) :: filename         ! name of ACE library file
+    character(10)           :: tablename        ! name of cross section table
+    type(Nuclide), pointer :: table => null()
+
+    filename = xsdatas(index)%path
+    tablename = xsdatas(index)%id
+
+    table => nuclides(index_table)
+
+    ! Check if input file exists and is readable
+    inquire(FILE=filename, EXIST=file_exists, READ=readable)
+    if (.not. file_exists) then
+       msg = "ACE library '" // trim(filename) // "' does not exist!"
+       call fatal_error(msg)
+    elseif (readable(1:3) == 'NO') then
+       msg = "ACE library '" // trim(filename) // "' is not readable! &
+            &Change file permissions with chmod command."
+       call fatal_error(msg)
+    end if
+
+    ! display message
+    msg = "Loading ACE cross section table: " // tablename
+    call message(msg, 6)
+
+    ! open file
+    open(file=filename, unit=in, status='old', & 
+         & action='read', iostat=ioError)
+    if (ioError /= 0) then
+       msg = "Error while opening file: " // filename
+       call fatal_error(msg)
+    end if
+
+    found_xs = .false.
+    do while (.not. found_xs)
+       call read_line(in, line, ioError)
+       if (ioError < 0) then
+          msg = "Could not find ACE table " // tablename // "."
+          call fatal_error(msg)
+       end if
+       call split_string(line, words, n)
+       if (trim(words(1)) == trim(tablename)) then
+          found_xs = .true.
+          table%name = words(1)
+          table%awr = str_to_real(words(2))
+          kT = str_to_real(words(3))
+          table%temp = kT / K_BOLTZMANN
+       end if
+       
+       ! Skip 5 lines
+       call skip_lines(in, 5, ioError)
+
+       ! Read NXS data
+       lines = 2
+       words_per_line = 8
+       call read_data(in, NXS, 16, lines, words_per_line)
+
+       ! Read JXS data
+       lines = 4
+       call read_data(in, JXS, 32, lines, words_per_line)
+
+       ! Calculate how many data points and lines in the XSS array
+       n = NXS(1)
+       lines = (n + 3)/4
+
+       if (found_xs) then
+          ! allocate storage for XSS array
+          allocate(XSS(n))
+          
+          ! Read XSS
+          words_per_line = 4
+          call read_data(in, XSS, n, lines, words_per_line)
+       else
+          call skip_lines(in, lines, ioError)
+       end if
+
+    end do
+
+    call read_esz(table)
+    call read_nu_data(table)
+    call read_reactions(table)
+    call read_angular_dist(table)
+    call read_energy_dist(table)
+    call read_unr_res(table)
+
+    ! Free memory from XSS array
+    if(allocated(XSS)) deallocate(XSS)
+    if(associated(table)) nullify(table)
+
+    close(unit=in)
+
+  end subroutine read_ACE_continuous
+
+!===============================================================================
+! READ_ESZ - reads through the ESZ block. This block contains the energy grid,
+! total xs, absorption xs, elastic scattering xs, and heating numbers.
+!===============================================================================
+
+  subroutine read_esz(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: NE ! number of energy points for total and elastic cross sections
+
+    ! determine number of energy points
+    NE = NXS(3)
+    table%n_grid = NE
+
+    ! allocate storage for arrays
+    allocate(table%energy(NE))
+    allocate(table%sigma_t(NE))
+    allocate(table%sigma_a(NE))
+    allocate(table%sigma_el(NE))
+    allocate(table%heating(NE))
+
+    ! read data from XSS -- right now the total, absorption and elastic
+    ! scattering are read in to these special arrays, but in reality, it should
+    ! be necessary to only store elastic scattering and possibly total
+    ! cross-section for total material xs generation.
+    XSS_index = 1
+    table%energy = get_real(NE)
+    table%sigma_t = get_real(NE)
+    table%sigma_a = get_real(NE)
+    table%sigma_el = get_real(NE)
+    table%heating = get_real(NE)
+    
+  end subroutine read_esz
+
+!===============================================================================
+! READ_NU_DATA reads data given on the number of neutrons emitted from fission
+! as a function of the incoming energy of a neutron. This data may be broken
+! down into prompt and delayed neutrons emitted as well.
+!===============================================================================
+
+  subroutine read_nu_data(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: i, j   ! loop indices
+    integer :: JXS2   ! location for fission nu data
+    integer :: JXS24  ! location for delayed neutron data
+    integer :: KNU    ! location for nu data
+    integer :: LNU    ! type of nu data (polynomial or tabular)
+    integer :: NC     ! number of polynomial coefficients
+    integer :: NR     ! number of interpolation regions
+    integer :: NE     ! number of energies 
+    integer :: NPCR   ! number of delayed neutron precursor groups
+    integer :: LED    ! location of energy distribution locators
+    integer :: LDIS   ! location of all energy distributions
+    integer :: LOCC   ! location of energy distributions for given MT
+    integer :: LNW    ! location of next energy distribution if multiple
+    integer :: LAW    ! secondary energy distribution law   
+    integer :: IDAT   ! location of first energy distribution for given MT
+    integer :: loc    ! locator
+    integer :: length             ! length of data to allocate
+    integer :: length_interp_data ! length of interpolation data
+
+    JXS2  = JXS(2)
+    JXS24 = JXS(24)
+
+    if (JXS2 == 0) then
+       ! =======================================================================
+       ! NO PROMPT/TOTAL NU DATA
+       table % nu_t_type = NU_NONE
+       table % nu_p_type = NU_NONE
+
+    elseif (XSS(JXS2) > 0) then
+       ! =======================================================================
+       ! PROMPT OR TOTAL NU DATA
+       KNU = JXS2
+       LNU = int(XSS(KNU))
+       if (LNU == 1) then
+          ! Polynomial data
+          table % nu_t_type = NU_POLYNOMIAL
+          table % nu_p_type = NU_NONE
+
+          ! allocate determine how many coefficients for polynomial
+          NC = int(XSS(KNU+1))
+          length = NC + 1
+       elseif (LNU == 2) then
+          ! Tabular data
+          table % nu_t_type = NU_TABULAR
+          table % nu_p_type = NU_NONE
+
+          ! determine number of interpolation regions and number of energies
+          NR = int(XSS(KNU+1))
+          NE = int(XSS(KNU+2+2*NR))
+          length = 2 + 2*NR + 2*NE
+       end if
+
+       ! allocate space for nu data storage
+       allocate(table % nu_t_data(length))
+
+       ! read data -- for polynomial, this is the number of coefficients and the
+       ! coefficients themselves, and for tabular, this is interpolation data
+       ! and tabular E/nu
+       XSS_index = KNU + 1
+       table % nu_t_data = get_real(length)
+
+    elseif (XSS(JXS2) < 0) then
+       ! =======================================================================
+       ! PROMPT AND TOTAL NU DATA -- read prompt data first
+       KNU = JXS2 + 1
+       LNU = XSS(KNU)
+       if (LNU == 1) then
+          ! Polynomial data
+          table % nu_p_type = NU_POLYNOMIAL
+
+          ! allocate determine how many coefficients for polynomial
+          NC = XSS(KNU+1)
+          length = NC + 1
+       elseif (LNU == 2) then
+          ! Tabular data
+          table % nu_p_type = NU_TABULAR
+
+          ! determine number of interpolation regions and number of energies
+          NR = XSS(KNU+1)
+          NE = XSS(KNU+2+2*NR)
+          length = 2 + 2*NR + 2*NE
+       end if
+
+       ! allocate space for nu data storage
+       allocate(table % nu_p_data(length))
+
+       ! read data
+       XSS_index = KNU + 1
+       table % nu_p_data = get_real(length)
+
+       ! Now read total nu data
+       KNU = JXS2 + abs(XSS(JXS2)) + 1
+       LNU = XSS(KNU)
+       if (LNU == 1) then
+          ! Polynomial data
+          table % nu_t_type = NU_POLYNOMIAL
+
+          ! allocate determine how many coefficients for polynomial
+          NC = int(XSS(KNU+1))
+          length = NC + 1
+       elseif (LNU == 2) then
+          ! Tabular data
+          table % nu_t_type = NU_TABULAR
+
+          ! determine number of interpolation regions and number of energies
+          NR = int(XSS(KNU+1))
+          NE = int(XSS(KNU+2+2*NR))
+          length = 2 + 2*NR + 2*NE
+       end if
+
+       ! allocate space for nu data storage
+       allocate(table % nu_t_data(length))
+
+       ! read data
+       XSS_index = KNU + 1
+       table % nu_t_data = get_real(length)
+    end if
+
+    if (JXS24 > 0) then
+       ! =======================================================================
+       ! DELAYED NU DATA
+
+       table % nu_d_type = NU_TABULAR
+       KNU = JXS24
+
+       ! determine size of tabular delayed nu data
+       NR = int(XSS(KNU+1))
+       NE = int(XSS(KNU+2+2*NR))
+       length = 2 + 2*NR + 2*NE
+
+       ! allocate space for delayed nu data
+       allocate(table % nu_d_data(length))
+       
+       ! read delayed nu data
+       XSS_index = KNU + 1
+       table % nu_d_data = get_real(length)
+ 
+       ! =======================================================================
+       ! DELAYED NEUTRON ENERGY DISTRIBUTION
+
+       ! Allocate space for secondary energy distribution
+       NPCR = NXS(8)
+       table % n_precursor = NPCR
+       allocate(table % nu_d_edist(NPCR))
+
+       LED  = JXS(26)
+       LDIS = JXS(27)
+
+       ! Loop over all delayed neutron precursor groups
+       do i = 1, NPCR
+          ! find location of energy distribution data
+          LOCC = XSS(LED + i - 1)
+
+          LNW  = XSS(LDIS + LOCC - 1)
+          LAW  = XSS(LDIS + LOCC)
+          IDAT = XSS(LDIS + LOCC + 1)
+          NR   = XSS(LDIS + LOCC + 2)
+          table % nu_d_edist(i) % law = LAW
+
+          ! allocate space for ENDF interpolation parameters
+          if (NR > 0) then
+             allocate(table % nu_d_edist(i) % nbt(NR))
+             allocate(table % nu_d_edist(i) % int(NR))
+          end if
+
+          ! read ENDF interpolation parameters
+          XSS_index = LDIS + LOCC + 3
+          table % nu_d_edist(i) % nbt = get_real(NR)
+          table % nu_d_edist(i) % int = get_real(NR)
+
+          ! allocate space for law validity data
+          NE = XSS(LDIS + LOCC + 3 + 2*NR)
+          allocate(table % nu_d_edist(i) % energy(NE))
+          allocate(table % nu_d_edist(i) % pvalid(NE))
+
+          length_interp_data = 5 + 2*(NR + NE)
+
+          ! read law validity data
+          XSS_index = LDIS + LOCC + 4 + 2*NR
+          table % nu_d_edist(i) % energy = get_real(NE)
+          table % nu_d_edist(i) % pvalid = get_real(NE)
+
+          ! Set index to beginning of IDAT array
+          loc = LDIS + IDAT - 2
+
+          ! determine length of energy distribution
+          length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
+
+          ! allocate secondary energy distribution array
+          allocate(table % nu_d_edist(i) % data(length))
+
+          ! read secondary energy distribution
+          XSS_index = loc + 1
+          table % nu_d_edist(i) % data = get_real(length)
+       end do
+
+       ! =======================================================================
+       ! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS
+
+       ! determine length of all precursor constants/yields/interp data
+       length = 0
+       loc = JXS(25)
+       do i = 1, NPCR
+          NR = XSS(loc + length + 1)
+          NE = XSS(loc + length + 2 + 2*NR)
+          length = length + 3 + 2*NR + 2*NE
+       end do
+
+       ! allocate space for precusor data
+       allocate(table % nu_d_precursor_data(length))
+
+       ! read delayed neutron precursor data
+       XSS_index = loc
+       table % nu_d_precursor_data = get_real(length)
+
+    else
+       table % nu_d_type = NU_NONE
+    end if
+
+  end subroutine read_nu_data
+
+!===============================================================================
+! READ_REACTIONS - Get the list of reaction MTs for this cross-section
+! table. The MT values are somewhat arbitrary. Also read in Q-values, neutron
+! multiplicities, and cross-sections.
+!===============================================================================
+
+  subroutine read_reactions(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: i     ! loop indices
+    integer :: LMT   ! index of MT list in XSS
+    integer :: NMT   ! Number of reactions
+    integer :: JXS4  ! index of Q values in XSS
+    integer :: JXS5  ! index of neutron multiplicities in XSS
+    integer :: JXS7  ! index of reactions cross-sections in XSS
+    integer :: LXS   ! location of cross-section locators
+    integer :: LOCA  ! location of cross-section for given MT
+    integer :: NE    ! number of energies for reaction
+    type(Reaction), pointer :: rxn => null()
+    
+    LMT  = JXS(3)
+    JXS4 = JXS(4)
+    JXS5 = JXS(5)
+    LXS  = JXS(6)
+    JXS7 = JXS(7)
+    NMT  = NXS(4)
+
+    ! allocate array of reactions. Add one since we need to include an elastic
+    ! scattering channel
+    table%n_reaction = NMT + 1
+    allocate(table%reactions(NMT+1))
+
+    ! Store elastic scattering cross-section on reaction one
+    rxn => table%reactions(1)
+    rxn%MT      = 2
+    rxn%Q_value = ZERO
+    rxn%TY      = 1
+    rxn%IE      = 1
+    allocate(rxn%sigma(table%n_grid))
+    rxn%sigma = table%sigma_el
+
+    do i = 1, NMT
+       rxn => table%reactions(i+1)
+
+       ! read MT number, Q-value, and neutrons produced
+       rxn % MT      = XSS(LMT+i-1)
+       rxn % Q_value = XSS(JXS4+i-1)
+       rxn % TY      = XSS(JXS5+i-1)
+
+       ! read cross section values
+       LOCA = XSS(LXS+i-1)
+       rxn % IE = XSS(JXS7 + LOCA - 1)
+       NE = XSS(JXS7 + LOCA)
+       allocate(rxn%sigma(NE))
+       XSS_index = JXS7 + LOCA + 1
+       rxn % sigma = get_real(NE)
+
+       ! set defaults
+       rxn % has_angle_dist  = .false.
+       rxn % has_energy_dist = .false.
+    end do
+
+  end subroutine read_reactions
+
+!===============================================================================
+! READ_ANGULAR_DIST parses the angular distribution for each reaction with
+! secondary neutrons
+!===============================================================================
+
+  subroutine read_angular_dist(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: JXS8   ! location of angular distribution locators
+    integer :: JXS9   ! location of angular distributions
+    integer :: LOCB   ! location of angular distribution for given MT
+    integer :: NE     ! number of incoming energies
+    integer :: NP     ! number of points for cosine distribution
+    integer :: LC     ! locator
+    integer :: i      ! index in reactions array
+    integer :: j      ! index over incoming energies
+    integer :: length ! length of data array to allocate
+    type(Reaction), pointer :: rxn => null()
+
+    JXS8 = JXS(8)
+    JXS9 = JXS(9)
+
+    ! loop over all reactions with secondary neutrons -- NXS(5) does not include
+    ! elastic scattering
+    do i = 1, NXS(5) + 1
+       rxn => table%reactions(i)
+
+       ! find location of angular distribution
+       LOCB = XSS(JXS8 + i - 1)
+       if (LOCB == -1) then
+          ! Angular distribution data are specified through LAWi = 44 in the DLW
+          ! block
+          cycle
+       elseif (LOCB == 0) then
+          ! No angular distribution data are given for this reaction, isotropic
+          ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0)
+          cycle
+       end if
+       rxn % has_angle_dist = .true.
+
+       ! allocate space for incoming energies and locations
+       NE = XSS(JXS9 + LOCB - 1)
+       rxn % adist % n_energy = NE
+       allocate(rxn % adist % energy(NE))
+       allocate(rxn % adist % type(NE))
+       allocate(rxn % adist % location(NE))
+
+       ! read incoming energy grid and location of tables
+       XSS_index = JXS9 + LOCB
+       rxn % adist % energy   = get_real(NE)
+       rxn % adist % location = get_int(NE)
+
+       ! determine dize of data block
+       length = 0
+       do j = 1, NE
+          LC = rxn % adist % location(j)
+          if (LC == 0) then
+             ! isotropic
+             rxn % adist % type(j) = ANGLE_ISOTROPIC
+          elseif (LC > 0) then
+             ! 32 equiprobable bins
+             rxn % adist % type(j) = ANGLE_32_EQUI
+             length = length + 33
+          elseif (LC < 0) then
+             ! tabular distribution
+             rxn % adist % type(j) = ANGLE_TABULAR
+             NP = XSS(JXS9 + abs(LC))
+             length = length + 2 + 3*NP
+          end if
+       end do
+       
+       ! allocate angular distribution data and read
+       allocate(rxn % adist % data(length))
+
+       ! read angular distribution -- currently this does not actually parse the
+       ! angular distribution tables for each incoming energy, that must be done
+       ! on-the-fly
+       LC = rxn % adist % location(1)
+       XSS_index = JXS9 + abs(LC) - 1
+       rxn % adist % data = get_real(length)
+
+       ! change location pointers since they are currently relative to JXS(9)
+       LC = abs(rxn % adist % location(1))
+       rxn % adist % location = abs(rxn % adist % location) - LC
+
+    end do
+
+  end subroutine read_angular_dist
+
+!===============================================================================
+! READ_ENERGY_DIST parses the secondary energy distribution for each reaction
+! with seconary neutrons (except elastic scattering)
+!===============================================================================
+
+  subroutine read_energy_dist(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: LED   ! location of energy distribution locators
+    integer :: LDIS  ! location of all energy distributions
+    integer :: LOCC  ! location of energy distributions for given MT
+    integer :: LNW   ! location of next energy distribution if multiple
+    integer :: LAW   ! secondary energy distribution law   
+    integer :: NR    ! number of interpolation regions
+    integer :: NE    ! number of incoming energies
+    integer :: IDAT  ! location of first energy distribution for given MT
+    integer :: loc   ! locator
+    integer :: length             ! length of data to allocate
+    integer :: length_interp_data ! length of interpolation data
+    integer :: i, j, k, l         ! indices
+    type(Reaction), pointer :: rxn => null()
+
+    LED  = JXS(10)
+    LDIS = JXS(11)
+
+    ! Loop over all reactions 
+    do i = 1, NXS(5)
+       rxn => table % reactions(i+1) ! skip over elastic scattering
+       rxn % has_energy_dist = .true.
+
+       ! find location of energy distribution data
+       LOCC = XSS(LED + i - 1)
+       
+       LNW  = XSS(LDIS + LOCC - 1)
+       LAW  = XSS(LDIS + LOCC)
+       IDAT = XSS(LDIS + LOCC + 1)
+       NR   = XSS(LDIS + LOCC + 2)
+       rxn % edist % law = LAW
+       rxn % edist % n_interp = NR
+
+       ! allocate space for ENDF interpolation parameters
+       if (NR > 0) then
+          allocate(rxn % edist % nbt(NR))
+          allocate(rxn % edist % int(NR))
+       end if
+       
+       ! read ENDF interpolation parameters
+       XSS_index = LDIS + LOCC + 3
+       rxn % edist % nbt = get_real(NR)
+       rxn % edist % int = get_real(NR)
+
+       ! allocate space for law validity data
+       NE = XSS(LDIS + LOCC + 3 + 2*NR)
+       allocate(rxn % edist % energy(NE))
+       allocate(rxn % edist % pvalid(NE))
+       
+       length_interp_data = 5 + 2*(NR + NE)
+       
+       ! read law validity data
+       XSS_index = LDIS + LOCC + 4 + 2*NR
+       rxn % edist % energy = get_real(NE)
+       rxn % edist % pvalid = get_real(NE)
+
+       ! Set index to beginning of IDAT array
+       loc = LDIS + IDAT - 2
+
+       ! determine length of energy distribution
+       length = length_energy_dist(loc, LAW, LOCC, length_interp_data)
+
+       ! allocate secondary energy distribution array
+       allocate(rxn % edist % data(length))
+
+       ! read secondary energy distribution
+       XSS_index = loc + 1
+       rxn % edist % data = get_real(length)
+    end do
+
+  end subroutine read_energy_dist
+
+!===============================================================================
+! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy
+! distribution array based on the secondary energy law and location in XSS
+!===============================================================================
+
+  function length_energy_dist(loc, law, LOCC, lid) result(length)
+
+    integer, intent(in) :: loc    ! location in XSS array
+    integer, intent(in) :: law    ! energy distribution law
+    integer, intent(in) :: LOCC   ! location of energy distribution
+    integer, intent(in) :: lid    ! length of interpolation data
+    integer             :: length ! length of energy distribution (LDAT)
+
+    integer :: i,j,k
+    integer :: NR    ! number of interpolation regions
+    integer :: NE    ! number of incoming energies
+    integer :: NP    ! number of points in outgoing energy distribution
+    integer :: NMU   ! number of points in outgoing cosine distribution
+    integer :: NRa   ! number of interpolation regions for Watt 'a'
+    integer :: NEa   ! number of energies for Watt 'a'
+    integer :: NRb   ! number of interpolation regions for Watt 'b'
+    integer :: NEb   ! number of energies for Watt 'b'
+
+    ! initialize length
+    length = 0
+
+    select case (law)
+    case (1)
+       ! Tabular equiprobable energy bins
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       NP = XSS(loc + 3 + 2*NR + NE)
+       length = 3 + 2*NR + NE + 3*NP*NE
+
+    case (2)
+       ! Discrete photon energy
+       length = 2
+
+    case (3)
+       ! Level scattering
+       length = 2
+
+    case (4)
+       ! Continuous tabular distribution
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       length = length + 2 + 2*NR + 2*NE
+       do i = 1,NE
+          ! determine length
+          NP = XSS(loc + length + 2)
+          length = length + 2 + 3*NP
+
+          ! adjust location for this block
+          j = loc + 2 + 2*NR + NE + i
+          XSS(j) = XSS(j) - LOCC - lid
+       end do
+
+    case (5)
+       ! General evaporation spectrum
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       NP = XSS(loc + 3 + 2*NR + 2*NE)
+       length = 3 + 2*NR + 2*NE + NP
+
+    case (7)
+       ! Maxwell fission spectrum
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       length = 3 + 2*NR + 2*NE
+
+    case (9)
+       ! Evaporation spectrum
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       length = 3 + 2*NR + 2*NE
+
+    case (11)
+       ! Watt spectrum
+       NRa = XSS(loc + 1)
+       NEa = XSS(loc + 2 + 2*NRa)
+       NRb = XSS(loc + 3 + 2*(NRa+NEa))
+       NEb = XSS(loc + 4 + 2*(NRa+NEa+NRb))
+       length = 5 + 2*(NRa + NEa + NRb + NEb)
+
+    case (44)
+       ! Kalbach-Mann correlated scattering
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       length = length + 2 + 2*NR + 2*NE
+       do i = 1,NE
+          NP = XSS(loc + length + 2)
+          length = length + 2 + 5*NP
+
+          ! adjust location for this block
+          j = loc + 2 + 2*NR + NE + i
+          XSS(j) = XSS(j) - LOCC - lid
+       end do
+
+    case (61)
+       ! Correlated energy and angle distribution
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       length = length + 2 + 2*NR + 2*NE
+       do i = 1,NE
+          ! outgoing energy distribution
+          NP = XSS(loc + length + 2)
+
+          ! adjust locators for angular distribution
+          do j = 1, NP
+             k = loc + length + 2 + 3*NP + j
+             if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
+          end do
+
+          length = length + 2 + 4*NP
+          do j = 1, NP
+             ! outgoing angle distribution -- NMU here is actually
+             ! referred to as NP in the MCNP documentation
+             NMU = XSS(loc + length + 2)
+             length = length + 2 + 3*NMU
+          end do
+
+          ! adjust locators for energy distribution
+          j = loc + 2 + 2*NR + NE + i
+          XSS(j) = XSS(j) - LOCC - lid
+       end do
+
+    case (66)
+       ! N-body phase space distribution
+       length = 2
+
+    case (67)
+       ! Laboratory energy-angle law
+       NR = XSS(loc + 1)
+       NE = XSS(loc + 2 + 2*NR)
+       NMU = XSS(loc + 4 + 2*NR + 2*NE)
+       length = 4 + 2*(NR + NE + NMU)
+
+    end select
+
+  end function length_energy_dist
+
+!===============================================================================
+! READ_UNR_RES reads in unresolved resonance probability tables if present.
+!===============================================================================
+
+  subroutine read_unr_res(table)
+
+    type(Nuclide), pointer :: table
+
+    integer :: JXS23 ! location of URR data
+    integer :: loc   ! locator
+    integer :: N     ! # of incident energies
+    integer :: M     ! # of probabilities
+    integer :: i     ! index over incoming energies
+    integer :: j     ! index over values
+    integer :: k     ! index over cumulative probabilities
+
+    ! determine locator for URR data
+    JXS23 = JXS(23)
+
+    ! check if URR data is present
+    if (JXS23 /= 0) then
+       table % urr_present = .true.
+       allocate(table % urr_params(6))
+       loc = JXS23
+    else
+       table % urr_present = .false.
+       return
+    end if
+
+    ! read parameters
+    table % urr_params(1) = XSS(loc)     ! # of incident energies
+    table % urr_params(2) = XSS(loc + 1) ! # of probabilities
+    table % urr_params(3) = XSS(loc + 2) ! interpolation parameter
+    table % urr_params(4) = XSS(loc + 3) ! inelastic competition flag
+    table % urr_params(5) = XSS(loc + 4) ! other absorption flag
+    table % urr_params(6) = XSS(loc + 5) ! factors flag
+
+    ! allocate incident energies and probability tables
+    N = table % urr_params(1)
+    M = table % urr_params(2)
+    allocate(table % urr_energy(N))
+    allocate(table % urr_prob(N,6,M))
+
+    ! read incident energies
+    XSS_index = loc + 6
+    table % urr_energy = get_real(N)
+
+    ! read probability tables
+    do i = 1, N
+       do j = 1, 6
+          table % urr_prob(i,j,1:M) = get_real(M)
+       end do
+    end do
+
+  end subroutine read_unr_res
+
+!===============================================================================
+! READ_ACE_THERMAL reads in a single ACE thermal scattering S(a,b) table. This
+! routine in particular reads the header data and then calls read_thermal_data
+! to parse the actual data blocks.
+!===============================================================================
+
+  subroutine read_ACE_thermal(index_table, index)
+
+    integer, intent(in) :: index_table ! index in sab_tables array
+    integer, intent(in) :: index       ! index in xsdatas array
+
+    integer                 :: in = 7           ! unit to read from
+    integer                 :: ioError          ! error status for file access
+    integer                 :: words_per_line   ! number of words per line (data)
+    integer                 :: lines            ! number of lines (data
+    integer                 :: n                ! number of data values
+    real(8)                 :: kT               ! ACE table temperature
+    logical                 :: file_exists      ! does ACE library exist?
+    logical                 :: found_xs         ! did we find table in library?
+    character(7)            :: readable         ! is ACE library readable?
+    character(MAX_LINE_LEN) :: msg              ! output/error message
+    character(MAX_LINE_LEN) :: line             ! single line to read
+    character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line
+    character(MAX_WORD_LEN) :: filename         ! name of ACE library file
+    character(10)           :: tablename        ! name of cross section table
+    type(SAB_Table), pointer :: table => null()
+
+    filename = xsdatas(index)%path
+    tablename = xsdatas(index)%id
+
+    table => sab_tables(index_table)
+
+    ! Check if input file exists and is readable
+    inquire(FILE=filename, EXIST=file_exists, READ=readable)
+    if (.not. file_exists) then
+       msg = "ACE library '" // trim(filename) // "' does not exist!"
+       call fatal_error(msg)
+    elseif (readable(1:3) == 'NO') then
+       msg = "ACE library '" // trim(filename) // "' is not readable! &
+            &Change file permissions with chmod command."
+       call fatal_error(msg)
+    end if
+
+    ! display message
+    msg = "Loading ACE cross section table: " // tablename
+    call message(msg, 6)
+
+    ! open file
+    open(file=filename, unit=in, status='old', & 
+         & action='read', iostat=ioError)
+    if (ioError /= 0) then
+       msg = "Error while opening file: " // filename
+       call fatal_error(msg)
+    end if
+
+    found_xs = .false.
+    do while (.not. found_xs)
+       call read_line(in, line, ioError)
+       if (ioError < 0) then
+          msg = "Could not find ACE table " // tablename // "."
+          call fatal_error(msg)
+       end if
+       call split_string(line, words, n)
+       if (trim(words(1)) == trim(tablename)) then
+          found_xs = .true.
+          table%name = words(1)
+          table%awr = str_to_real(words(2))
+          kT = str_to_real(words(3))
+          table%temp = kT / K_BOLTZMANN
+       end if
+       
+       ! Skip 5 lines
+       call skip_lines(in, 5, ioError)
+
+       ! Read NXS data
+       lines = 2
+       words_per_line = 8
+       call read_data(in, NXS, 16, lines, words_per_line)
+
+       ! Read JXS data
+       lines = 4
+       call read_data(in, JXS, 32, lines, words_per_line)
+
+       ! Calculate how many data points and lines in the XSS array
+       n = NXS(1)
+       lines = (n + 3)/4
+
+       if (found_xs) then
+          ! allocate storage for XSS array
+          allocate(XSS(n))
+          
+          ! Read XSS
+          words_per_line = 4
+          call read_data(in, XSS, n, lines, words_per_line)
+       else
+          call skip_lines(in, lines, ioError)
+       end if
+
+    end do
+
+    call read_thermal_data(table)
+
+    ! Free memory from XSS array
+    if(allocated(XSS)) deallocate(XSS)
+    if(associated(table)) nullify(table)
+
+    close(unit=in)
+
+  end subroutine read_ACE_thermal
+
+!===============================================================================
+! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding
+! secondary energy/angle distributions derived from experimental S(a,b)
+! data. Namely, in MCNP language, this routine reads the ITIE, ITCE, ITXE, and
+! ITCA blocks.
+!===============================================================================
+
+  subroutine read_thermal_data(table)
+
+    type(SAB_Table), pointer :: table
+
+    integer :: i      ! index for incoming energies
+    integer :: j      ! index for outgoing energies
+    integer :: k      ! index for outoging angles
+    integer :: loc    ! location in XSS array
+    integer :: NE_in  ! number of incoming energies
+    integer :: NE_out ! number of outgoing energies
+    integer :: NMU    ! number of outgoing angles
+    integer :: JXS4   ! location of elastic energy table
+
+    ! read number of inelastic energies and allocate arrays
+    NE_in = XSS(JXS(1))
+    table % n_inelastic_e_in = NE_in
+    allocate(table % inelastic_e_in(NE_in))
+    allocate(table % inelastic_sigma(NE_in))
+
+    ! read inelastic energies and cross-sections
+    XSS_index = JXS(1) + 1
+    table % inelastic_e_in = get_real(NE_in)
+    table % inelastic_sigma = get_real(NE_in)
+
+    ! allocate space for outgoing energy/angle for inelastic
+    ! scattering
+    NE_out = NXS(4)
+    NMU = NXS(3) + 1
+    allocate(table % inelastic_e_out(NE_out, NE_in))
+    allocate(table % inelastic_mu(NMU, NE_out, NE_in))
+
+    ! read outgoing energy/angle distribution for inelastic scattering
+    loc = JXS(3) - 1
+    do i = 1, NE_in
+       do j = 1, NE_out
+          ! read outgoing energy
+          table % inelastic_e_out(j,i) = XSS(loc + 1)
+
+          ! read outgoing angles for this outgoing energy
+          do k = 1, NMU
+             table % inelastic_mu(k,j,i) = XSS(loc + 1 + k)
+          end do
+
+          ! advance pointer
+          loc = loc + 1 + NMU
+       end do
+    end do
+
+    ! read number of elastic energies and allocate arrays
+    JXS4 = JXS(4)
+    if (JXS4 /= 0) then
+       NE_in = XSS(JXS4)
+       table % n_elastic_e_in = NE_in
+       allocate(table % elastic_e_in(NE_in))
+       allocate(table % elastic_P(NE_in))
+
+       ! determine whether sigma=P or sigma = P/E
+       table % n_elastic_type = NXS(5)
+    else
+       table % n_elastic_e_in = 0
+    end if
+
+    ! allocate space for outgoing energy/angle for elastic scattering
+    NMU = NXS(6) + 1
+    table % n_elastic_mu = NMU
+    if (NMU > 0) then
+       allocate(table % elastic_mu(NMU, NE_in))
+    end if
+
+    ! read equiprobable outgoing cosines for elastic scattering each
+    ! incoming energy
+    if (JXS4 /= 0 .and. NMU /= 0) then
+       loc = JXS(6) - 1
+       do i = 1, NE_in
+          do j = 1, NMU
+             table % elastic_mu(j,i) = XSS(loc + j)
+          end do
+          loc = loc + NMU
+       end do
+    end if
+    
+  end subroutine read_thermal_data
+
+!===============================================================================
+! GET_INT returns an array of integers read from the current position in the XSS
+! array
+!===============================================================================
+
+  function get_int(n_values) result(array)
+
+    integer, intent(in) :: n_values        ! number of values to read
+    integer             :: array(n_values) ! array of values
+
+    array = int(XSS(XSS_index:XSS_index + n_values - 1))
+    XSS_index = XSS_index + n_values
+
+  end function get_int
+
+!===============================================================================
+! GET_REAL returns an array of real(8)s read from the current position in the
+! XSS array
+!===============================================================================
+
+  function get_real(n_values) result(array)
+
+    integer, intent(in) :: n_values        ! number of values to read
+    real(8)             :: array(n_values) ! array of values
+
+    array = XSS(XSS_index:XSS_index + n_values - 1)
+    XSS_index = XSS_index + n_values
+
+  end function get_real
+
+!===============================================================================
+! GET_MACRO_XS
+!===============================================================================
+
+  function get_macro_xs(p, mat, MT) result(xs)
+
+    type(Particle), pointer    :: p
+    type(Material), pointer    :: mat
+    integer,        intent(in) :: MT
+    real(8)                    :: xs
+
+    integer :: i, j
+    integer :: n_isotopes
+    integer :: IE
+    real(8) :: density
+    real(8) :: density_i
+    real(8) :: sigma_i
+    real(8) :: f
+    type(Nuclide), pointer :: table => null()
+    type(Reaction),   pointer :: rxn => null()
+
+    ! initialize xs
+    xs = ZERO
+
+    ! find material atom density
+    density = cMaterial%atom_density
+
+    ! loop over all isotopes in material
+    n_isotopes = cMaterial % n_isotopes
+    do i = 1, n_isotopes
+       table => nuclides(cMaterial % table(i))
+
+       ! determine nuclide atom density
+       density_i = cMaterial%atom_percent(i) * density
+
+       ! search nuclide energy grid
+       IE = table%grid_index(p % IE)
+       f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
+       
+       ! handle special case of total cross section
+       if (MT == 1) then
+          xs = xs + density * (ONE-f) * table%sigma_t(IE) + & 
+               & f * (table%sigma_t(IE+1))
+          cycle
+       end if
+
+       ! loop over reactions in isotope
+       do j = 1, table % n_reaction
+          rxn => table % reactions(i)
+
+          ! check for matching MT
+          if (MT /= rxn % MT) cycle
+
+          ! if energy is below threshold for this reaction, skip it
+          if (IE < rxn % IE) cycle
+
+          ! add to cumulative probability
+          sigma_i = (ONE-f) * rxn%sigma(IE-rxn%IE+1) + & 
+               & f * (rxn%sigma(IE-rxn%IE+2))
+       end do
+
+       ! calculate nuclide macroscopic cross-section
+       xs = xs + density_i * sigma_i
+    end do
+
+  end function get_macro_xs
+
 !===============================================================================
 ! READ_XSDATA reads the data in a SERPENT xsdata file and builds a dictionary to
 ! find cross-section information later on.
@@ -21,10 +1236,10 @@ contains
     character(*), intent(in) :: path
 
     type(xsData), pointer    :: iso => null()
-    character(max_line_len)  :: line
-    character(max_line_len)  :: msg
-    character(max_word_len)  :: words(max_words)
-    character(max_word_len)  :: filename
+    character(MAX_LINE_LEN)  :: line
+    character(MAX_LINE_LEN)  :: msg
+    character(MAX_WORD_LEN)  :: words(MAX_WORDS)
+    character(MAX_WORD_LEN)  :: filename
     integer                  :: i, n
     integer                  :: in = 7
     logical                  :: file_exists
@@ -128,7 +1343,7 @@ contains
   subroutine material_total_xs()
 
     type(Material),      pointer :: mat => null()
-    type(AceContinuous), pointer :: table => null()
+    type(Nuclide), pointer :: table => null()
     integer :: i, j, k
     integer :: index
     integer :: IE
@@ -141,7 +1356,7 @@ contains
     real(8) :: E_i1
     real(8) :: sigma_i
     real(8) :: sigma_i1
-    character(max_line_len) :: msg
+    character(MAX_LINE_LEN) :: msg
 
     msg = "Creating material total cross-sections..."
     call message(msg, 4)
@@ -160,7 +1375,7 @@ contains
        ! loop over each isotope
        do j = 1, mat%n_isotopes
           index = mat % table(j)
-          table => xs_continuous(index)
+          table => nuclides(index)
 
           ! find atom density of isotope
           density_i = density * mat % atom_percent(j)
diff --git a/src/cross_section_header.f90 b/src/cross_section_header.f90
new file mode 100644
index 0000000000..fcd9c57996
--- /dev/null
+++ b/src/cross_section_header.f90
@@ -0,0 +1,136 @@
+module cross_section_header
+
+  implicit none
+
+!===============================================================================
+! DISTANGLE contains data for a tabular secondary angle distribution whether it
+! be tabular or 32 equiprobable cosine bins
+!===============================================================================
+
+  type DistAngle
+     integer              :: n_energy    ! # of incoming energies
+     real(8), allocatable :: energy(:)   ! incoming energy grid
+     integer, allocatable :: type(:)     ! type of distribution
+     integer, allocatable :: location(:) ! location of each table
+     real(8), allocatable :: data(:)     ! angular distribution data
+  end type DistAngle
+
+!===============================================================================
+! DISTENERGY contains data for a secondary energy distribution for all
+! scattering laws
+!===============================================================================
+
+  type DistEnergy
+     integer :: law                    ! secondary distribution law
+     integer :: n_interp               ! # of interpolation regions
+     integer, allocatable :: nbt(:)    ! ENDF interpolation parameters
+     integer, allocatable :: int(:)    ! ''
+     integer :: n_energy               ! # of energies for law validity
+     real(8), allocatable :: energy(:) ! energy grid for law validity
+     real(8), allocatable :: pvalid(:) ! probability of law validity
+     real(8), allocatable :: data(:)   ! energy distribution data
+  end type DistEnergy
+
+!===============================================================================
+! REACTION contains the cross-section and secondary energy and angle
+! distributions for a single reaction in a continuous-energy ACE-format table
+!===============================================================================
+
+  type Reaction
+     integer :: MT                     ! ENDF MT value
+     real(8) :: Q_value                ! Reaction Q value
+     integer :: TY                     ! Number of neutrons released
+     integer :: IE                     ! Starting energy grid index
+     real(8), allocatable :: sigma(:)  ! Cross section values
+     logical :: has_angle_dist         ! Angle distribution present?
+     logical :: has_energy_dist        ! Energy distribution present?
+     type(DistAngle)  :: adist      ! Secondary angular distribution
+     type(DistEnergy) :: edist      ! Secondary energy distribution
+  end type Reaction
+
+!===============================================================================
+! NUCLIDE contains all the data for an ACE-format continuous-energy cross
+! section. The ACE format (A Compact ENDF format) is used in MCNP and several
+! other Monte Carlo codes.
+!===============================================================================
+
+  type Nuclide
+     character(20) :: name
+     real(8) :: awr
+     real(8) :: temp
+     integer :: n_grid
+     integer, allocatable :: grid_index(:)
+     real(8), allocatable :: energy(:)
+     real(8), allocatable :: sigma_t(:)
+     real(8), allocatable :: sigma_a(:)
+     real(8), allocatable :: sigma_el(:)
+     real(8), allocatable :: heating(:)
+
+     ! Total fission neutron emission
+     integer :: nu_t_type
+     real(8), allocatable :: nu_t_data(:)
+
+     ! Prompt fission neutron emission
+     integer :: nu_p_type
+     real(8), allocatable :: nu_p_data(:)
+     
+     ! Delayed fission neutron emission
+     integer :: nu_d_type
+     integer :: n_precursor
+     real(8), allocatable :: nu_d_data(:)
+     real(8), allocatable :: nu_d_precursor_data(:)
+     type(DistEnergy), allocatable :: nu_d_edist(:)
+
+     ! Unresolved resonance data
+     logical :: urr_present
+     integer, allocatable :: urr_params(:)
+     real(8), allocatable :: urr_energy(:)
+     real(8), allocatable :: urr_prob(:,:,:)
+
+     ! Reactions
+     integer :: n_reaction
+     type(Reaction), pointer :: reactions(:) => null()
+
+  end type Nuclide
+
+!===============================================================================
+! SAB_TABLE contains S(a,b) data for thermal neutron scattering, typically off
+! of light isotopes such as water, graphite, Be, etc
+!===============================================================================
+     
+  type SAB_Table
+     character(20) :: name
+     real(8) :: awr
+     real(8) :: temp
+     integer :: n_inelastic_e_in
+     integer :: n_inelastic_e_out
+     integer :: n_inelastic_mu
+     real(8), allocatable :: inelastic_e_in(:)
+     real(8), allocatable :: inelastic_sigma(:) 
+     real(8), allocatable :: inelastic_e_out(:,:)
+     real(8), allocatable :: inelastic_mu(:,:,:)
+     integer :: n_elastic_e_in
+     integer :: n_elastic_type
+     integer :: n_elastic_mu
+     real(8), allocatable :: elastic_e_in(:)
+     real(8), allocatable :: elastic_P(:)
+     real(8), allocatable :: elastic_mu(:,:)
+  end type SAB_Table
+
+!===============================================================================
+! XSDATA contains data read in from a SERPENT xsdata file
+!===============================================================================
+
+  type xsData
+     character(10) :: alias
+     character(10) :: id
+     integer :: type
+     integer :: zaid
+     integer :: isomeric
+     real(8) :: awr
+     real(8) :: temp
+     integer :: binary
+     character(150) :: path
+  end type xsData
+
+end module cross_section_header
diff --git a/src/data_structures.f90 b/src/datatypes.f90
similarity index 98%
rename from src/data_structures.f90
rename to src/datatypes.f90
index de4f66d21e..b7739aa081 100644
--- a/src/data_structures.f90
+++ b/src/datatypes.f90
@@ -1,7 +1,7 @@
-module data_structures
+module datatypes
 
 !===============================================================================
-! DATA_STRUCTURES module
+! DATATYPES module
 !
 ! This module implements a dictionary that has (key,value) pairs. This data
 ! structure is used to provide lookup features, e.g. cells and surfaces by name.
@@ -12,13 +12,13 @@ module data_structures
 ! types (changes made to ListData, dict_create, dict_add_key, and dict_get_key).
 !===============================================================================
 
-  use types
+  use datatypes_header
 
   implicit none
 
-  integer, parameter :: hash_size  = 4993
-  integer, parameter :: multiplier = 31
-  integer, parameter :: DICT_NULL = 0
+  integer, parameter :: HASH_SIZE       = 4993
+  integer, parameter :: HASH_MULTIPLIER = 31
+  integer, parameter :: DICT_NULL       = 0
 
 !===============================================================================
 ! LIST Interfaces -- these allow one to use a single subroutine or function for
@@ -642,9 +642,9 @@ contains
     integer :: i
 
     allocate(dict)
-    allocate(dict%table(hash_size))
+    allocate(dict%table(HASH_SIZE))
 
-    do i = 1,hash_size
+    do i = 1,HASH_SIZE
        dict%table(i)%list => null()
     enddo
 
@@ -661,9 +661,9 @@ contains
     integer :: i
 
     allocate(dict)
-    allocate(dict%table(hash_size))
+    allocate(dict%table(HASH_SIZE))
 
-    do i = 1,hash_size
+    do i = 1,HASH_SIZE
        dict%table(i)%list => null()
     enddo
 
@@ -998,12 +998,12 @@ contains
     val = 0
 
     do i = 1,len(key)
-       val = multiplier * val + ichar(key(i:i))
+       val = HASH_MULTIPLIER * val + ichar(key(i:i))
     enddo
 
     ! Added the absolute val on val-1 since the sum in the do loop is
     ! susceptible to integer overflow
-    val = 1 + mod(abs(val-1), hash_size)
+    val = 1 + mod(abs(val-1), HASH_SIZE)
 
   end function dict_ci_hashkey
 
@@ -1025,7 +1025,7 @@ contains
 
     ! Added the absolute val on val-1 since the sum in the do loop is
     ! susceptible to integer overflow
-    val = 1 + mod(abs(key-1), hash_size)
+    val = 1 + mod(abs(key-1), HASH_SIZE)
 
   end function dict_ii_hashkey
 
@@ -1093,4 +1093,4 @@ contains
 
   end function dict_ii_keys
 
-end module data_structures
+end module datatypes
diff --git a/src/datatypes_header.f90 b/src/datatypes_header.f90
new file mode 100644
index 0000000000..7d332cea63
--- /dev/null
+++ b/src/datatypes_header.f90
@@ -0,0 +1,113 @@
+module datatypes_header
+
+!===============================================================================
+! DATATYPES_HEADER module contains the derived types for (key,value) pairs,
+! lists, and dictionaries. It would be nice to add sets as well. This module
+! could be made a LOT cleaner by using unlimited polymorphism, but we will wait
+! on that since compiler support is still limited (namely gfortran which doesn't
+! support it yet)
+!===============================================================================
+
+  implicit none
+
+!===============================================================================
+! KEYVALUECI stores the (key,value) pair for a dictionary where the key is a
+! string and the value is an integer. Note that we need to store the key in
+! addition to the value for collision resolution.
+!===============================================================================
+
+  ! Key length for dictionary
+  integer, parameter :: DICT_KEY_LENGTH = 20
+
+  type KeyValueCI
+     character(len=DICT_KEY_LENGTH) :: key
+     integer                        :: value
+  end type KeyValueCI
+
+!===============================================================================
+! KEYVALUEII stores the (key,value) pair for a dictionary where the key is an
+! integer and the value is an integer. Note that we need to store the key in
+! addition to the value for collision resolution.
+!===============================================================================
+
+  type KeyValueII
+     integer :: key
+     integer :: value
+  end type KeyValueII
+
+!===============================================================================
+! LISTKEYVALUECI stores a linked list of (key,value) pairs where the key is a
+! character and the value is an integer
+!===============================================================================
+
+  type ListKeyValueCI
+     type(ListKeyValueCI), pointer :: next => null()
+     type(KeyValueCI)              :: data
+  end type ListKeyValueCI
+
+!===============================================================================
+! LISTKEYVALUEII stores a linked list of (key,value) pairs where the key is a
+! character and the value is an integer
+!===============================================================================
+
+  type ListKeyValueII
+     type(ListKeyValueII), pointer :: next => null()
+     type(KeyValueII)              :: data
+  end type ListKeyValueII
+
+!===============================================================================
+! LISTREAL stores a linked list of real values. This is used for constructing a
+! unionized energy grid.
+!===============================================================================
+
+  type ListReal
+     type(ListReal), pointer :: next => null()
+     real(8)                 :: data
+  end type ListReal
+
+!===============================================================================
+! LISTINT stores a linked list of integer values.
+!===============================================================================
+
+  type ListInt
+     type(ListInt), pointer :: next => null()
+     integer                :: data
+  end type ListInt
+
+!===============================================================================
+! HASHLISTCI - Since it's not possible to directly do an array of pointers, this
+! derived type provides a pointer
+!===============================================================================
+
+  type HashListCI
+     type(ListKeyValueCI), pointer :: list => null()
+  end type HashListCI
+
+!===============================================================================
+! HASHLISTII - Since it's not possible to directly do an array of pointers, this
+! derived type provides a pointer
+!===============================================================================
+
+  type HashListII
+     type(ListKeyValueII), pointer :: list => null()
+  end type HashListII
+
+!===============================================================================
+! DICTIONARYCI provides a dictionary data structure of (key,value) pairs where
+! the keys are strings and values are integers.
+!===============================================================================
+
+  type DictionaryCI
+     type(HashListCI), pointer :: table(:) => null()
+  end type DictionaryCI
+
+!===============================================================================
+! DICTIONARYII provides a dictionary data structure of (key,value) pairs where
+! the keys and values are both integers.
+!===============================================================================
+
+  type DictionaryII
+     type(HashListII), pointer :: table(:) => null()
+  end type DictionaryII
+
+end module datatypes_header
diff --git a/src/endf.f90 b/src/endf.f90
index d87e98d670..8232e8f5a4 100644
--- a/src/endf.f90
+++ b/src/endf.f90
@@ -1,6 +1,6 @@
 module endf
 
-  use global
+  use constants
   use string, only: int_to_str
 
 contains
diff --git a/src/energy_grid.f90 b/src/energy_grid.f90
index 5b96c072ed..fa31ae929d 100644
--- a/src/energy_grid.f90
+++ b/src/energy_grid.f90
@@ -1,9 +1,10 @@
 module energy_grid
 
+  use constants,        only: MAX_LINE_LEN
+  use datatypes,        only: list_insert, list_size, list_delete
+  use datatypes_header, only: ListReal
   use global
-  use output, only: message
-  use data_structures, only: list_insert, list_size, list_delete, &
-       &                     list_size
+  use output,           only: message
 
 contains
 
@@ -19,13 +20,12 @@ contains
 
     type(ListReal), pointer :: list => null()
     type(ListReal), pointer :: current => null()
-
-    type(Material),       pointer :: mat => null()
-    type(AceContinuous),  pointer :: table => null()
-    type(AceReaction),    pointer :: rxn => null()
+    type(Material), pointer :: mat => null()
+    type(Nuclide),  pointer :: table => null()
+    type(Reaction), pointer :: rxn => null()
     integer :: i, j
     integer :: n
-    character(max_line_len) :: msg
+    character(MAX_LINE_LEN) :: msg
 
     msg = "Creating unionized energy grid..."
     call message(msg, 5)
@@ -36,7 +36,7 @@ contains
        
        ! loop over all isotopes
        do j = 1, mat%n_isotopes
-          table => xs_continuous(mat%table(j))
+          table => nuclides(mat%table(j))
 
           ! loop over energy points
           n = size(table%energy)
@@ -159,26 +159,26 @@ contains
 
     integer :: i, j
     integer :: index
-    type(AceContinuous), pointer :: acecont
+    type(Nuclide), pointer :: nuc
 
     real(8) :: union_energy
     real(8) :: energy
     
 
-    do i = 1, n_continuous
-       acecont => xs_continuous(i)
-       allocate(acecont%grid_index(n_grid))
+    do i = 1, n_nuclides
+       nuc => nuclides(i)
+       allocate(nuc%grid_index(n_grid))
 
        index = 1
-       energy = acecont%energy(index)
+       energy = nuc%energy(index)
 
        do j = 1, n_grid
           union_energy = e_grid(j)
           if (union_energy >= energy) then
              index = index + 1
-             energy = acecont%energy(index)
+             energy = nuc%energy(index)
           end if
-          acecont%grid_index(j) = index-1
+          nuc%grid_index(j) = index-1
        end do
 
     end do
diff --git a/src/error.f90 b/src/error.f90
index 410691e1d7..83d95f695a 100644
--- a/src/error.f90
+++ b/src/error.f90
@@ -1,10 +1,12 @@
 module error
 
   use ISO_FORTRAN_ENV
-  use global
+
+  use global, only: master, free_memory
 
   implicit none
 
+  ! Short names for output and error units
   integer :: ou = OUTPUT_UNIT
   integer :: eu = ERROR_UNIT
 
diff --git a/src/fileio.f90 b/src/fileio.f90
index 1f96d83ef0..ebc32ce2f6 100644
--- a/src/fileio.f90
+++ b/src/fileio.f90
@@ -1,13 +1,18 @@
 module fileio
 
+  use constants
+  use cross_section_header, only: xsData
+  use datatypes,            only: dict_create, dict_add_key, dict_get_key, &
+                                  dict_has_key, DICT_NULL, dict_keys, list_size
+  use datatypes_header,     only: DictionaryII, ListKeyValueII, ListKeyValueCI
+  use error,                only: fatal_error, warning, ERROR_INT, ERROR_REAL
+  use geometry_header,      only: Cell, Surface, BASE_UNIVERSE
   use global
-  use types,  only: Cell, Surface, ExtSource, ListKeyValueCI
-  use string, only: split_string_wl, lower_case, split_string, concatenate, &
-                    is_number, str_to_int, int_to_str, str_to_real
-  use error,  only: fatal_error, warning
-  use output, only: message
-  use data_structures, only: dict_create, dict_add_key, dict_get_key, &
-                             dict_has_key, DICT_NULL, dict_keys, list_size
+  use output,               only: message
+  use source_header,        only: ExtSource
+  use string,               only: split_string_wl, lower_case, split_string, &
+                                  concatenate, is_number, str_to_int, & 
+                                  int_to_str, str_to_real
 
   implicit none
 
@@ -37,7 +42,7 @@ contains
 
     integer                 :: argc        ! number of command line arguments
     logical                 :: file_exists ! does specified input file exist?
-    character(max_line_len) :: msg         ! error message
+    character(MAX_LINE_LEN) :: msg         ! error message
     character(7)            :: readable    ! is input file readable?
 
     argc = COMMAND_ARGUMENT_COUNT()
@@ -79,9 +84,9 @@ contains
     integer                 :: n                ! number of words on a line
     integer                 :: count            ! number of cells in a universe
     integer                 :: universe_num     ! user-specified universe #
-    character(max_line_len) :: line             ! a line of words in input file
-    character(max_line_len) :: msg              ! output/error message
-    character(max_word_len) :: words(max_words) ! words on a line
+    character(MAX_LINE_LEN) :: line             ! a line of words in input file
+    character(MAX_LINE_LEN) :: msg              ! output/error message
+    character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line
     type(ListKeyValueII), pointer :: key_list => null()
     type(Universe),       pointer :: univ => null()
 
@@ -230,9 +235,9 @@ contains
     integer :: index_material          ! index in materials array
     integer :: index_source            ! index in source array (?)
     integer :: index_tally             ! index in tally array
-    character(max_line_len) :: line             ! a line of words
-    character(max_line_len) :: msg              ! output/error message
-    character(max_word_len) :: words(max_words) ! words on a single line
+    character(MAX_LINE_LEN) :: line             ! a line of words
+    character(MAX_LINE_LEN) :: msg              ! output/error message
+    character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a single line
 
     msg = "Second pass through input file..."
     call message(msg, 5)
@@ -335,7 +340,7 @@ contains
     integer                 :: index    ! index in surfaces/materials array 
     integer                 :: surf_num ! user-specified surface number
     integer                 :: bc       ! boundary condition
-    character(max_line_len) :: msg      ! output/error message
+    character(MAX_LINE_LEN) :: msg      ! output/error message
     type(Cell),           pointer :: c => null()
     type(Surface),        pointer :: surf => null()
     type(ListKeyValueII), pointer :: key_list => null()
@@ -469,8 +474,8 @@ contains
     integer                 :: i            ! index for surface list in a cell
     integer                 :: universe_num ! user-specified universe number
     integer                 :: n_surfaces   ! number of surfaces in a cell
-    character(max_line_len) :: msg          ! output/error message
-    character(max_word_len) :: word         ! single word
+    character(MAX_LINE_LEN) :: msg          ! output/error message
+    character(MAX_WORD_LEN) :: word         ! single word
     type(Cell), pointer     :: c => null()
 
     c => cells(index)
@@ -562,8 +567,8 @@ contains
     integer                 :: ioError     ! error status for file access
     integer                 :: i           ! index for surface coefficients
     integer                 :: coeffs_reqd ! number of coefficients are required
-    character(max_line_len) :: msg         ! output/error message
-    character(max_word_len) :: word        ! single word
+    character(MAX_LINE_LEN) :: msg         ! output/error message
+    character(MAX_WORD_LEN) :: word        ! single word
     type(Surface), pointer  :: surf => null()
 
     surf => surfaces(index)
@@ -651,8 +656,8 @@ contains
 
     integer                 :: surface_uid ! User-specified uid of surface
     integer                 :: bc          ! Boundary condition
-    character(max_word_len) :: word        ! Boundary condition (in input file)
-    character(max_line_len) :: msg         ! Output/error message
+    character(MAX_WORD_LEN) :: word        ! Boundary condition (in input file)
+    character(MAX_LINE_LEN) :: msg         ! Output/error message
 
     ! Read surface identifier
     surface_uid = str_to_int(words(2))
@@ -696,8 +701,8 @@ contains
     integer                 :: n_y          ! number of lattice cells in y direction
     integer                 :: i,j          ! loop indices for Lattice % universes
     integer                 :: index_word   ! index in words array
-    character(max_line_len) :: msg          ! output/error/message
-    character(max_word_len) :: word         ! single word
+    character(MAX_LINE_LEN) :: msg          ! output/error/message
+    character(MAX_WORD_LEN) :: word         ! single word
     type(Lattice), pointer  :: lat => null()
     
     lat => lattices(index)
@@ -779,8 +784,8 @@ contains
     integer                 :: i           ! index in values list
     integer                 :: ioError     ! error status for file access
     integer                 :: values_reqd ! # of values required to specify source
-    character(max_line_len) :: msg  ! output/error message
-    character(max_word_len) :: word ! single word
+    character(MAX_LINE_LEN) :: msg  ! output/error message
+    character(MAX_WORD_LEN) :: word ! single word
 
     ! Read source type
     word = words(2)
@@ -826,8 +831,8 @@ contains
     integer :: cell_uid
     integer :: r_bins, c_bins, e_bins
     real(8) :: E
-    character(max_word_len) :: word
-    character(max_line_len) :: msg
+    character(MAX_WORD_LEN) :: word
+    character(MAX_LINE_LEN) :: msg
     type(Tally), pointer :: t => null()
 
     t => tallies(index)
@@ -977,7 +982,7 @@ contains
     integer                 :: i          ! index over isotopes
     integer                 :: ioError    ! error status for file access
     integer                 :: n_isotopes ! number of isotopes in material
-    character(max_line_len) :: msg        ! output/error message
+    character(MAX_LINE_LEN) :: msg        ! output/error message
     type(Material), pointer :: mat => null()
 
     ! Check for correct number of arguments
@@ -1032,7 +1037,7 @@ contains
     logical        :: percent_in_atom ! isotopes specified in atom percent?
     logical        :: density_in_atom ! density specified in atom/b-cm?
     character(10)  :: key             ! name of isotopes, e.g. 92235.03c
-    character(max_line_len) :: msg    ! output/error message
+    character(MAX_LINE_LEN) :: msg    ! output/error message
     type(xsData),   pointer :: iso => null()
     type(Material), pointer :: mat => null()
     
@@ -1122,7 +1127,7 @@ contains
     character(*), intent(in) :: words(n_words) ! words on criticality card
     integer,      intent(in) :: n_words        ! number of words
 
-    character(max_line_len) :: msg ! output/error message
+    character(MAX_LINE_LEN) :: msg ! output/error message
 
     ! Set problem type to criticality
     problem_type = PROB_CRITICALITY
@@ -1173,12 +1178,12 @@ contains
   subroutine get_next_line(unit, words, n, ioError)
 
     integer,      intent(in)  :: unit             ! unit to read from
-    character(*), intent(out) :: words(max_words) ! words read
+    character(*), intent(out) :: words(MAX_WORDS) ! words read
     integer,      intent(out) :: n                ! number of words
     integer,      intent(out) :: ioError          ! error status
 
-    character(max_line_len) :: line                   ! single line
-    character(max_word_len) :: local_words(max_words) ! words on one line
+    character(MAX_LINE_LEN) :: line                   ! single line
+    character(MAX_WORD_LEN) :: local_words(MAX_WORDS) ! words on one line
     integer                 :: index                  ! index of words
 
     index = 0
@@ -1222,7 +1227,7 @@ contains
     integer, intent(out) :: ioError ! error status 
 
     integer                 :: i        ! index for number of lines
-    character(max_line_len) :: tmp      ! single line
+    character(MAX_LINE_LEN) :: tmp      ! single line
 
     do i = 1, n_lines
        read(UNIT=unit, FMT='(A)', IOSTAT=ioError) tmp
diff --git a/src/geometry.f90 b/src/geometry.f90
index e0b996db2b..5d4b4321ae 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -1,11 +1,13 @@
 module geometry
 
+  use constants
+  use datatypes,       only: dict_get_key
+  use error,           only: fatal_error
+  use geometry_header, only: Cell, Surface, Universe, Lattice
   use global
-  use types,  only: Cell, Surface
-  use output, only: message
-  use error,  only: fatal_error
-  use data_structures, only: dict_get_key
-  use string, only: int_to_str
+  use output,          only: message
+  use particle_header, only: Particle
+  use string,          only: int_to_str
 
   implicit none
      
@@ -28,7 +30,7 @@ contains
     integer :: i               ! index of surfaces in cell
     integer :: surf_num        ! index in surfaces array (with sign)
     integer :: current_surface ! current surface of particle (with sign)
-    character(max_line_len) :: msg      ! output/error message
+    character(MAX_LINE_LEN) :: msg      ! output/error message
     type(Surface), pointer  :: surf => null()
 
     current_surface = p%surface
@@ -94,7 +96,7 @@ contains
     type(Particle), pointer :: p          ! pointer to particle
     logical,  intent(inout) :: found      ! particle found?
 
-    character(max_line_len) :: msg        ! error message
+    character(MAX_LINE_LEN) :: msg        ! error message
     integer                 :: i          ! index over cells
     integer                 :: x, y
     type(Cell),     pointer :: c          ! pointer to cell
@@ -178,7 +180,7 @@ contains
     integer                 :: i          ! index of neighbors
     integer                 :: index_cell ! index in cells array
     logical                 :: found      ! particle found in universe?
-    character(max_line_len) :: msg        ! output/error message?
+    character(MAX_LINE_LEN) :: msg        ! output/error message?
     type(Surface),  pointer :: surf
     type(Cell),     pointer :: c
     type(Universe), pointer :: lower_univ => null()
@@ -288,7 +290,7 @@ contains
     real(8)        :: x0       ! half the width of lattice element
     real(8)        :: y0       ! half the height of lattice element
     logical        :: found    ! particle found in cell?
-    character(max_line_len) :: msg  ! output/error message
+    character(MAX_LINE_LEN) :: msg  ! output/error message
     type(Lattice),  pointer :: lat
     type(Universe), pointer :: univ
 
@@ -407,11 +409,11 @@ contains
     real(8) :: a,b,c,k      ! quadratic equation coefficients
     real(8) :: quad         ! discriminant of quadratic equation
     logical :: on_surface   ! is particle on surface?
-    character(max_line_len) :: msg   ! output/error message
+    character(MAX_LINE_LEN) :: msg   ! output/error message
     type(Cell),    pointer  :: cell_p => null()
     type(Cell),    pointer  :: parent_p => null()
     type(Surface), pointer  :: surf_p => null()
-    type(LatticE), pointer  :: lat => null()
+    type(Lattice), pointer  :: lat => null()
 
     cell_p => cells(p%cell)
 
@@ -922,7 +924,7 @@ contains
     integer, allocatable :: count_positive(:) ! # of cells on positive side
     integer, allocatable :: count_negative(:) ! # of cells on negative side
     logical :: positive   ! positive side specified in surface list
-    character(max_line_len) :: msg ! output/error message
+    character(MAX_LINE_LEN) :: msg ! output/error message
     type(Cell),    pointer  :: c
     type(Surface), pointer  :: surf
 
diff --git a/src/geometry_header.f90 b/src/geometry_header.f90
new file mode 100644
index 0000000000..7b44c0a570
--- /dev/null
+++ b/src/geometry_header.f90
@@ -0,0 +1,77 @@
+module geometry_header
+
+  implicit none
+
+!===============================================================================
+! UNIVERSE defines a geometry that fills all phase space
+!===============================================================================
+
+  type Universe
+     integer :: uid                   ! Unique ID
+     integer :: type                  ! Type
+     integer :: level                 ! Level of universe (0=base)
+     integer :: n_cells               ! # of cells within
+     integer, allocatable :: cells(:) ! List of cells within
+     real(8) :: x0                    ! Translation in x-coordinate
+     real(8) :: y0                    ! Translation in y-coordinate
+     real(8) :: z0                    ! Translation in z-coordinate
+     integer, allocatable :: tallies(:)
+  end type Universe
+
+!===============================================================================
+! LATTICE is an ordered array of elements (either rectangular, hexagonal, or
+! triangular)
+!===============================================================================
+
+  type Lattice
+     integer :: uid      ! Universe number for lattice
+     integer :: type     ! Type of lattice (rectangular, hex, etc)
+     integer :: level    ! Level of lattice
+     integer :: n_x      ! number of lattice cells in x-direction
+     integer :: n_y      ! number of lattice cells in y-direction
+     real(8) :: x0       ! x-coordinate of lattice origin
+     real(8) :: y0       ! y-coordinate of lattice origin
+     real(8) :: width_x  ! width of lattice cell 
+     real(8) :: width_y  ! width of lattice cell
+     integer, allocatable :: element(:,:) ! specified universes
+     integer, allocatable :: tallies(:)
+  end type Lattice
+
+!===============================================================================
+! SURFACE type defines a first- or second-order surface that can be used to
+! construct closed volumes (cells)
+!===============================================================================
+
+  type Surface
+     integer :: uid                    ! Unique ID
+     integer :: type                   ! Type of surface
+     real(8), allocatable :: coeffs(:) ! Definition of surface
+     integer, allocatable :: & 
+          & neighbor_pos(:), &         ! List of cells on positive side
+          & neighbor_neg(:)            ! List of cells on negative side
+     integer :: bc                     ! Boundary condition
+  end type Surface
+
+!===============================================================================
+! CELL defines a closed volume by its bounding surfaces
+!===============================================================================
+
+  type Cell
+     integer :: uid        ! Unique ID
+     integer :: type       ! Type of cell (normal, universe, lattice)
+     integer :: universe   ! universe # this cell is in
+     integer :: fill       ! universe # filling this cell
+     integer :: parent     ! cell within which this cell resides
+     integer :: material   ! Material within cell (0 for universe)
+     integer :: n_surfaces ! Number of surfaces within
+     integer, allocatable :: & 
+          & surfaces(:)    ! List of surfaces bounding cell -- note that
+                           ! parentheses, union, etc operators will be listed
+                           ! here too
+     integer, allocatable :: tallies(:)
+  end type Cell
+
+  ! array index of universe 0
+  integer :: BASE_UNIVERSE
+
+end module geometry_header
diff --git a/src/global.f90 b/src/global.f90
index f3e30fa6aa..465d0a2007 100644
--- a/src/global.f90
+++ b/src/global.f90
@@ -1,6 +1,15 @@
 module global
 
-  use types
+  use bank_header,          only: Bank
+  use constants
+  use cross_section_header, only: Nuclide, SAB_Table, xsData
+  use datatypes_header,     only: DictionaryII, DictionaryCI
+  use geometry_header,      only: Cell, Universe, Lattice, Surface
+  use material_header,      only: Material
+  use particle_header,      only: Particle
+  use source_header,        only: ExtSource
+  use tally_header,         only: Tally
+  use timing,               only: Timer
 
 #ifdef MPI
   use mpi
@@ -8,12 +17,6 @@ module global
 
   implicit none
 
-  ! Maximum number of words in a single line, length of line, and length of
-  ! single word
-  integer, parameter :: max_words    = 500
-  integer, parameter :: max_line_len = 250
-  integer, parameter :: max_word_len = 150
-
   ! Main arrays for cells, surfaces, materials
   type(Cell),     allocatable, target :: cells(:)
   type(Universe), allocatable, target :: universes(:)
@@ -37,16 +40,15 @@ module global
   type(DictionaryII), pointer :: lattice_dict
   type(DictionaryII), pointer :: surface_dict
   type(DictionaryII), pointer :: material_dict
-  type(DictionaryII), pointer :: isotope_dict
   type(DictionaryII), pointer :: tally_dict
   type(DictionaryCI), pointer :: xsdata_dict
-  type(DictionaryCI), pointer :: ace_dict
+  type(DictionaryCI), pointer :: nuclide_dict
 
   ! Cross section arrays
-  type(AceContinuous), allocatable, target :: xs_continuous(:)
-  type(AceThermal),    allocatable, target :: xs_thermal(:)
-  integer :: n_continuous
-  integer :: n_thermal
+  type(Nuclide),   allocatable, target :: nuclides(:)
+  type(SAB_Table), allocatable, target :: sab_tables(:)
+  integer :: n_nuclides
+  integer :: n_sab_tables
 
   ! Current cell, surface, material
   type(Cell),     pointer :: cCell
@@ -87,191 +89,23 @@ module global
   logical :: mpi_enabled ! is MPI in use and initialized?
 
   ! Timing variables
-  type(TimerObj) :: time_total       ! timer for total run
-  type(TimerObj) :: time_init        ! timer for initialization
-  type(TimerObj) :: time_intercycle  ! timer for intercycle synchronization
-  type(TimerObj) :: time_compute     ! timer for computation
+  type(Timer) :: time_total       ! timer for total run
+  type(Timer) :: time_init        ! timer for initialization
+  type(Timer) :: time_intercycle  ! timer for intercycle synchronization
+  type(Timer) :: time_compute     ! timer for computation
 
   ! Paths to input file, cross section data, etc
-  character(max_word_len) :: & 
+  character(MAX_WORD_LEN) :: & 
        & path_input,         &
        & path_xsdata
-  integer, parameter :: UNIT_LOG = 9 ! unit # for writing log file
-
-  ! Physical constants
-  real(8), parameter ::            &
-       & PI           = 3.1415926535898_8, & ! pi
-       & MASS_NEUTRON = 1.0086649156,      & ! mass of a neutron
-       & MASS_PROTON  = 1.00727646677,     & ! mass of a proton
-       & AMU          = 1.66053873e-27,    & ! 1 amu in kg
-       & N_AVOGADRO   = 0.602214179,       & ! Avogadro's number in 10^24/mol
-       & K_BOLTZMANN  = 8.617342e-5,       & ! Boltzmann constant in eV/K
-       & INFINITY     = huge(0.0_8),       & ! positive infinity
-       & ZERO         = 0.0_8,             &
-       & ONE          = 1.0_8,             &
-       & TWO          = 2.0_8
-
-  ! Boundary conditions
-  integer, parameter ::    &
-       & BC_TRANSMIT = 0,  & ! Transmission boundary condition (default)
-       & BC_VACUUM   = 1,  & ! Vacuum boundary condition
-       & BC_REFLECT  = 2,  & ! Reflecting boundary condition
-       & BC_PERIODIC = 3     ! Periodic boundary condition
-
-  ! Logical operators for cell definitions
-  integer, parameter ::                &
-       & OP_LEFT_PAREN  = huge(0),     & ! Left parentheses
-       & OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
-       & OP_UNION       = huge(0) - 2, & ! Union operator
-       & OP_DIFFERENCE  = huge(0) - 3    ! Difference operator
-
-  ! Cell types
-  integer, parameter ::    &
-       & CELL_NORMAL  = 1, & ! Cell with a specified material
-       & CELL_FILL    = 2, & ! Cell filled by a separate universe
-       & CELL_LATTICE = 3    ! Cell filled with a lattice
-
-  ! Lattice types
-  integer, parameter :: &
-       & LATTICE_RECT = 1, &
-       & LATTICE_HEX  = 2
-
-  ! array index of universe 0
-  integer :: BASE_UNIVERSE
-
-  ! Surface types
-  integer, parameter ::    &
-       & SURF_PX     =  1, & ! Plane parallel to x-plane 
-       & SURF_PY     =  2, & ! Plane parallel to y-plane 
-       & SURF_PZ     =  3, & ! Plane parallel to z-plane 
-       & SURF_PLANE  =  4, & ! Arbitrary plane
-       & SURF_CYL_X  =  5, & ! Cylinder along x-axis
-       & SURF_CYL_Y  =  6, & ! Cylinder along y-axis
-       & SURF_CYL_Z  =  7, & ! Cylinder along z-axis
-       & SURF_SPHERE =  8, & ! Sphere
-       & SURF_BOX_X  =  9, & ! Box extending infinitely in x-direction
-       & SURF_BOX_Y  = 10, & ! Box extending infinitely in y-direction
-       & SURF_BOX_Z  = 11, & ! Box extending infinitely in z-direction
-       & SURF_BOX    = 12, & ! Rectangular prism
-       & SURF_GQ     = 13    ! General quadratic surface
-
-  ! Surface senses
-  integer, parameter ::      &
-       & SENSE_POSITIVE = 1, &
-       & SENSE_NEGATIVE = -1
-
-  ! Source types
-  integer, parameter ::   &
-       & SRC_BOX     = 1, & ! Source in a rectangular prism
-       & SRC_CELL    = 2, & ! Source in a cell
-       & SRC_SURFACE = 3    ! Source on a surface
 
   ! Problem type
   integer :: problem_type
-  integer, parameter ::        &
-       & PROB_SOURCE      = 1, & ! External source problem
-       & PROB_CRITICALITY = 2    ! Criticality problem
-
-  ! Angular distribution type
-  integer, parameter :: & 
-       & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
-       & ANGLE_32_EQUI   = 2, & ! 32 equiprobable bins
-       & ANGLE_TABULAR   = 3    ! Tabular angular distribution
-
-  ! Interpolation flag
-  integer, parameter ::     &
-       & HISTOGRAM     = 1, & ! y is constant in x
-       & LINEAR_LINEAR = 2, & ! y is linear in x
-       & LINEAR_LOG    = 3, & ! y is linear in ln(x)
-       & LOG_LINEAR    = 4, & ! ln(y) is linear in x
-       & LOG_LOG       = 5    ! ln(y) is linear in ln(x)
-
-  ! Particle type
-  integer, parameter :: &
-       & NEUTRON  = 1, &
-       & PHOTON   = 2, &
-       & ELECTRON = 3
-
-  ! Reaction types
-  integer, parameter :: &
-       & TOTAL_XS  = 1, &
-       & ELASTIC   = 2, &
-       & N_LEVEL   = 4, &
-       & MISC      = 5, &
-       & N_2ND     = 11, &
-       & N_2N      = 16, &
-       & N_3N      = 17, &
-       & FISSION   = 18, &
-       & N_F       = 19, &
-       & N_NF      = 20, &
-       & N_2NF     = 21, &
-       & N_NA      = 22, &
-       & N_N3A     = 23, &
-       & N_2NA     = 24, &
-       & N_3NA     = 25, &
-       & N_NP      = 28, &
-       & N_N2A     = 29, &
-       & N_2N2A    = 30, &
-       & N_ND      = 32, &
-       & N_NT      = 33, &
-       & N_N3HE    = 34, &
-       & N_ND2A    = 35, &
-       & N_NT2A    = 36, &
-       & N_4N      = 37, &
-       & N_3NF     = 38, &
-       & N_2NP     = 41, &
-       & N_3NP     = 42, &
-       & N_N2P     = 44, &
-       & N_NPA     = 45, &
-       & N_N1      = 51, &
-       & N_N40     = 90, &
-       & N_NC      = 91, &
-       & N_GAMMA   = 102, &
-       & N_P       = 103, &
-       & N_D       = 104, &
-       & N_T       = 105, &
-       & N_3HE     = 106, &
-       & N_A       = 107, &
-       & N_2A      = 108, &
-       & N_3A      = 109, &
-       & N_2P      = 111, &
-       & N_PA      = 112, &
-       & N_T2A     = 113, &
-       & N_D2A     = 114, &
-       & N_PD      = 115, &
-       & N_PT      = 116, &
-       & N_DA      = 117
-
-  ! Tally distinguishers
-  integer, parameter :: &
-       & TALLY_FLUX = -1, & ! tally the flux only
-       & TALLY_ALL  =  0, & ! tally all reactions
-       & TALLY_BINS =  1, & ! tally individual (reactions, cells)
-       & TALLY_SUM  =  2    ! tally sum of specified (reactions, cells)
-
-  ! Fission neutron emission (nu) type
-  integer, parameter ::     &
-       & NU_NONE       = 0, & ! No nu values (non-fissionable)
-       & NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
-       & NU_TABULAR    = 2    ! Nu values given by tabular distribution
-
-  ! Integer code for read error -- better hope this number is never
-  ! used in an input file!
-  integer, parameter :: &
-       & ERROR_INT = -huge(0)
-
-  real(8), parameter :: &
-       & ERROR_REAL = -huge(0.0_8) * 0.917826354_8
 
   ! The verbosity controls how much information will be printed to the
   ! screen and in logs
   integer :: verbosity
 
-  ! Versioning numbers
-  integer, parameter :: VERSION_MAJOR = 0
-  integer, parameter :: VERSION_MINOR = 2
-  integer, parameter :: VERSION_RELEASE = 2
-
 contains
 
 !===============================================================================
@@ -311,8 +145,8 @@ contains
 
     ! Deallocate cross section data
     if (allocated(xsdatas)) deallocate(xsdatas)
-    if (allocated(xs_continuous)) deallocate(xs_continuous)
-    if (allocated(xs_thermal)) deallocate(xs_thermal)
+    if (allocated(nuclides)) deallocate(nuclides)
+    if (allocated(sab_tables)) deallocate(sab_tables)
 
     ! Deallocate energy grid
     if (allocated(e_grid)) deallocate(e_grid)
diff --git a/src/logging.f90 b/src/logging.f90
index 0d951471dc..5c8279e228 100644
--- a/src/logging.f90
+++ b/src/logging.f90
@@ -1,6 +1,7 @@
 module logging
 
-  use global
+  use constants, only: MAX_WORD_LEN, UNIT_LOG
+  use global,    only: path_input
 
   implicit none
 
@@ -12,7 +13,7 @@ contains
 
   subroutine create_log()
 
-    character(max_word_len) :: path_log     ! path of log file
+    character(MAX_WORD_LEN) :: path_log     ! path of log file
     logical                 :: file_exists  ! does log file already exist?
     integer                 :: ioError      ! error status for file access
 
diff --git a/src/main.f90 b/src/main.f90
index 71773b42ed..06d89900da 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -1,22 +1,24 @@
 program main
 
+  use constants
+  use cross_section,   only: read_xs, read_xsdata, material_total_xs
+  use energy_grid,     only: unionized_grid, original_indices
+  use fileio,          only: read_input, read_command_line, read_count, &
+                             normalize_ao, build_universe
+  use geometry,        only: neighbor_lists
+  use geometry_header, only: Universe, BASE_UNIVERSE
   use global
-  use fileio,        only: read_input, read_command_line, read_count, &
-                           normalize_ao, build_universe
-  use output,        only: title, echo_input, message, print_summary, &
-                           print_particle, header, print_runtime
-  use geometry,      only: neighbor_lists
-  use mcnp_random,   only: RN_init_problem, RN_init_particle
-  use source,        only: init_source, get_source_particle
-  use physics,       only: transport
-  use cross_section, only: read_xsdata, material_total_xs
-  use ace,           only: read_xs
-  use energy_grid,   only: unionized_grid, original_indices
-  use mpi_routines,  only: setup_mpi, synchronize_bank
-  use score,         only: calculate_keff
-  use logging,       only: create_log
-  use timing,        only: timer_start, timer_stop
-  use string,        only: int_to_str
+  use logging,         only: create_log
+  use mcnp_random,     only: RN_init_problem, RN_init_particle
+  use mpi_routines,    only: setup_mpi, synchronize_bank
+  use output,          only: title, echo_input, message, print_summary, &
+                             print_particle, header, print_runtime
+  use particle_header, only: Particle
+  use physics,         only: transport
+  use score,           only: calculate_keff
+  use source,          only: init_source, get_source_particle
+  use string,          only: int_to_str
+  use timing,          only: timer_start, timer_stop
 
 #ifdef MPI
   use mpi
diff --git a/src/material_header.f90 b/src/material_header.f90
new file mode 100644
index 0000000000..1d539a3533
--- /dev/null
+++ b/src/material_header.f90
@@ -0,0 +1,21 @@
+module material_header
+
+  implicit none
+
+!===============================================================================
+! MATERIAL describes a material by its constituent isotopes
+!===============================================================================
+
+  type Material
+     integer :: uid
+     integer :: n_isotopes
+     character(10), allocatable :: names(:)  ! isotope names
+     integer, allocatable :: isotopes(:)     ! index in xsdata list
+     integer, allocatable :: table(:)        ! index in nuclides array
+     real(8)              :: atom_density    ! total atom density in atom/b-cm
+     real(8), allocatable :: atom_percent(:)
+     real(8), allocatable :: total_xs(:)     ! macroscopic cross-section
+     integer :: sab_table
+  end type Material
+
+end module material_header
diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90
index f5a7a0c826..c7a2902889 100644
--- a/src/mpi_routines.f90
+++ b/src/mpi_routines.f90
@@ -1,9 +1,10 @@
 module mpi_routines
 
-  use global
+  use constants,   only: max_line_len
   use error,       only: fatal_error
-  use output,      only: message
+  use global
   use mcnp_random, only: rang, RN_init_particle, RN_skip
+  use output,      only: message
   use source,      only: copy_from_bank, source_index
 
 #ifdef MPI
@@ -33,7 +34,7 @@ contains
     integer        :: bank_types(4)  ! Datatypes
     integer(MPI_ADDRESS_KIND) :: bank_disp(4)   ! Displacements
     integer(MPI_ADDRESS_KIND) :: base
-    character(max_line_len) :: msg            ! Error message
+    character(MAX_LINE_LEN) :: msg            ! Error message
     type(Bank)     :: b
 
     mpi_enabled = .true.
@@ -129,7 +130,7 @@ contains
          & temp_sites(:),      & ! local array of extra sites on each node
          & left_bank(:),       & ! bank sites to send/recv to or from left node
          & right_bank(:)         ! bank sites to send/recv to or fram right node
-    character(max_line_len) :: msg
+    character(MAX_LINE_LEN) :: msg
 
     msg = "Collecting number of fission sites..."
     call message(msg, 8)
diff --git a/src/output.f90 b/src/output.f90
index 4eba340092..ba3eb4f70f 100644
--- a/src/output.f90
+++ b/src/output.f90
@@ -1,14 +1,17 @@
 module output
 
   use ISO_FORTRAN_ENV
-  use global
-  use types,           only: Cell, Universe, Surface
-  use data_structures, only: dict_get_key
+
+  use constants
+  use datatypes,       only: dict_get_key
   use endf,            only: reaction_name
+  use geometry_header, only: Cell, Universe, Surface
+  use global
   use string,          only: upper_case, int_to_str, real_to_str
 
   implicit none
 
+  ! Short names for output and error units
   integer :: ou = OUTPUT_UNIT
   integer :: eu = ERROR_UNIT
 
@@ -273,7 +276,7 @@ contains
 
   subroutine print_reaction(rxn)
 
-    type(AceReaction), pointer :: rxn
+    type(Reaction), pointer :: rxn
 
     write(ou,*) 'Reaction ' // reaction_name(rxn % MT)
     write(ou,*) '    MT = ' // int_to_str(rxn % MT)
@@ -474,17 +477,17 @@ contains
 
     type(Material), pointer :: mat
 
-    integer                      :: i
-    integer                      :: n_lines
-    real(8)                      :: density
-    character(max_line_len)      :: string
-    type(AceContinuous), pointer :: table => null()
+    integer                 :: i
+    integer                 :: n_lines
+    real(8)                 :: density
+    character(max_line_len) :: string
+    type(Nuclide),  pointer :: table => null()
 
     write(ou,*) 'Material ' // int_to_str(mat % uid)
     write(ou,*) '    Atom Density = ' // trim(real_to_str(mat % atom_density)) &
          & // ' atom/b-cm'
     do i = 1, mat % n_isotopes
-       table => xs_continuous(mat % table(i))
+       table => nuclides(mat % table(i))
        density = mat % atom_density * mat % atom_percent(i)
        string = '    ' // trim(table % name) // ' = ' // &
             & trim(real_to_str(density)) // ' atom/b-cm'
diff --git a/src/particle_header.f90 b/src/particle_header.f90
new file mode 100644
index 0000000000..fc2650eb89
--- /dev/null
+++ b/src/particle_header.f90
@@ -0,0 +1,30 @@
+module particle_header
+
+  implicit none
+
+!===============================================================================
+! PARTICLE describes the state of a particle being transported through the
+! geometry
+!===============================================================================
+
+  type Particle
+    integer(8) :: uid          ! Unique ID
+    integer    :: type         ! Particle type (n, p, e, etc)
+    real(8)    :: xyz(3)       ! location
+    real(8)    :: xyz_local(3) ! local location (after transformations)
+    real(8)    :: uvw(3)       ! directional cosines
+    real(8)    :: wgt          ! particle weight
+    real(8)    :: E            ! energy
+    integer    :: IE           ! index on energy grid
+    real(8)    :: interp       ! interpolation factor for energy grid
+    integer    :: cell         ! current cell
+    integer    :: universe     ! current universe
+    integer    :: lattice      ! current lattice
+    integer    :: surface      ! current surface
+    integer    :: index_x      ! lattice index for x direction
+    integer    :: index_y      ! lattice index for y direction
+    logical    :: alive        ! is particle alive?
+    integer    :: n_coll       ! # of collisions
+  end type Particle
+
+end module particle_header
diff --git a/src/physics.f90 b/src/physics.f90
index 51dd328883..4c3df4583f 100644
--- a/src/physics.f90
+++ b/src/physics.f90
@@ -1,16 +1,19 @@
 module physics
 
+  use constants
+  use cross_section_header, only: Nuclide, Reaction, DistEnergy
+  use endf,                 only: reaction_name
+  use error,                only: fatal_error, warning
+  use geometry,             only: find_cell, dist_to_boundary, cross_surface, &
+                                  cross_lattice
+  use geometry_header,      only: Universe, BASE_UNIVERSE
   use global
-  use geometry,    only: find_cell, dist_to_boundary, cross_surface, &
-       &                 cross_lattice
-  use types,       only: Particle
-  use mcnp_random, only: rang
-  use output,      only: message, print_particle
-  use error,       only: fatal_error, warning
-  use search,      only: binary_search
-  use endf,        only: reaction_name
-  use score,       only: score_tally
-  use string,      only: int_to_str
+  use mcnp_random,          only: rang
+  use output,               only: message, print_particle
+  use particle_header,      only: Particle
+  use score,                only: score_tally
+  use search,               only: binary_search
+  use string,               only: int_to_str
 
   implicit none
 
@@ -33,7 +36,7 @@ contains
     real(8)        :: f              ! interpolation factor
     logical        :: found_cell     ! found cell which particle is in?
     logical        :: in_lattice     ! is surface crossing in lattice?
-    character(max_line_len) :: msg   ! output/error message
+    character(MAX_LINE_LEN) :: msg   ! output/error message
     type(Universe), pointer :: univ
 
     if (p % cell == 0) then
@@ -159,9 +162,9 @@ contains
     real(8) :: r1         ! random number
     real(8) :: flux       ! collision estimator of flux
     real(8), allocatable :: Sigma_rxn(:) ! macroscopic xs for each nuclide
-    character(max_line_len)       :: msg ! output/error message
-    type(AceContinuous), pointer :: table
-    type(AceReaction),   pointer :: rxn
+    character(MAX_LINE_LEN)       :: msg ! output/error message
+    type(Nuclide),  pointer :: table
+    type(Reaction), pointer :: rxn
 
     density = cMaterial%atom_density
 
@@ -170,7 +173,7 @@ contains
     n_isotopes = cMaterial%n_isotopes
     allocate(Sigma_rxn(n_isotopes))
     do i = 1, n_isotopes
-       table => xs_continuous(cMaterial%table(i))
+       table => nuclides(cMaterial%table(i))
 
        ! determine nuclide atom density
        density_i = cMaterial%atom_percent(i) * density
@@ -195,7 +198,7 @@ contains
 
     ! Get table, total xs, interpolation factor
     density_i = cMaterial%atom_percent(i)*density
-    table => xs_continuous(cMaterial%table(i))
+    table => nuclides(cMaterial%table(i))
     sigma = Sigma_rxn(i) / density_i
     IE = table%grid_index(p % IE)
     f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
@@ -273,9 +276,9 @@ contains
 
   subroutine elastic_scatter(p, table, rxn)
 
-    type(Particle),      pointer :: p
-    type(AceContinuous), pointer :: table
-    type(AceReaction),   pointer :: rxn
+    type(Particle), pointer :: p
+    type(Nuclide),  pointer :: table
+    type(Reaction), pointer :: rxn
 
     real(8) :: awr ! atomic weight ratio of target
     real(8) :: phi ! azimuthal angle
@@ -348,9 +351,9 @@ contains
 
   subroutine n_fission(p, table, rxn)
 
-    type(Particle),      pointer :: p
-    type(AceContinuous), pointer :: table
-    type(AceReaction),   pointer :: rxn
+    type(Particle), pointer :: p
+    type(Nuclide),  pointer :: table
+    type(Reaction), pointer :: rxn
 
     integer :: i           ! loop index
     integer :: j           ! index on nu energy grid / precursor group
@@ -374,7 +377,7 @@ contains
     real(8) :: xi          ! random number
     real(8) :: yield       ! delayed neutron precursor yield
     real(8) :: prob        ! cumulative probability
-    character(max_line_len) :: msg  ! error message
+    character(MAX_LINE_LEN) :: msg  ! error message
 
     ! copy energy of neutron
     E = p % E
@@ -638,9 +641,9 @@ contains
 
   subroutine inelastic_scatter(p, table, rxn)
 
-    type(Particle),      pointer :: p
-    type(AceContinuous), pointer :: table
-    type(AceReaction),   pointer :: rxn
+    type(Particle), pointer :: p
+    type(Nuclide),  pointer :: table
+    type(Reaction), pointer :: rxn
 
     integer :: law
     real(8) :: A          ! atomic weight ratio of nuclide
@@ -650,7 +653,7 @@ contains
     real(8) :: E_cm       ! outgoing energy in center-of-mass
     real(8) :: u,v,w      ! direction cosines
     real(8) :: Q          ! Q-value of reaction
-    character(max_line_len) :: msg ! error message
+    character(MAX_LINE_LEN) :: msg ! error message
     
     ! copy energy of neutron
     E_in = p % E
@@ -707,9 +710,9 @@ contains
 
   subroutine n_xn(p, table, rxn)
 
-    type(Particle),      pointer :: p
-    type(AceContinuous), pointer :: table
-    type(AceReaction),   pointer :: rxn
+    type(Particle), pointer :: p
+    type(Nuclide),  pointer :: table
+    type(Reaction), pointer :: rxn
 
     integer :: i           ! loop index
     integer :: n_secondary ! number of secondary particles
@@ -721,7 +724,7 @@ contains
     real(8) :: E_cm        ! outgoing energy in center-of-mass
     real(8) :: u,v,w       ! direction cosines
     real(8) :: Q           ! Q-value of reaction
-    character(max_line_len) :: msg  ! error message
+    character(MAX_LINE_LEN) :: msg  ! error message
     
     ! copy energy of neutron
     E_in = p % E
@@ -788,7 +791,7 @@ contains
     type(Particle), pointer :: p
 
     integer                 :: cell_num ! user-specified cell number
-    character(max_line_len) :: msg      ! output/error message
+    character(MAX_LINE_LEN) :: msg      ! output/error message
 
     p % alive = .false.
     if (verbosity >= 10) then
@@ -807,8 +810,8 @@ contains
 
   function sample_angle(rxn, E) result(mu)
 
-    type(AceReaction), pointer    :: rxn ! reaction
-    real(8),           intent(in) :: E   ! incoming energy
+    type(Reaction), pointer    :: rxn ! reaction
+    real(8),        intent(in) :: E   ! incoming energy
 
     real(8)        :: xi      ! random number on [0,1)
     integer        :: interp  ! type of interpolation
@@ -826,7 +829,7 @@ contains
     real(8)        :: c_k     ! cumulative frequency at k
     real(8)        :: c_k1    ! cumulative frequency at k+1
     real(8)        :: p0,p1   ! probability distribution
-    character(max_line_len) :: msg     ! error message
+    character(MAX_LINE_LEN) :: msg     ! error message
 
     ! check if reaction has angular distribution -- if not, sample outgoing
     ! angle isotropically
@@ -974,12 +977,12 @@ contains
 
   subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
 
-    type(AceDistEnergy), intent(inout) :: edist
-    real(8), intent(in)                :: E_in
-    real(8), intent(out)               :: E_out
-    real(8), intent(inout), optional   :: mu_out
-    real(8), intent(in),    optional   :: A
-    real(8), intent(in),    optional   :: Q
+    type(DistEnergy),  intent(inout) :: edist
+    real(8), intent(in)              :: E_in
+    real(8), intent(out)             :: E_out
+    real(8), intent(inout), optional :: mu_out
+    real(8), intent(in),    optional :: A
+    real(8), intent(in),    optional :: Q
 
     integer :: i           ! index on incoming energy grid
     integer :: k           ! sampled index on outgoing grid
@@ -1025,7 +1028,7 @@ contains
     real(8) :: E_max       ! parameter for n-body dist
     real(8) :: x, y, v     ! intermediate variables for n-body dist
     real(8) :: r1, r2, r3, r4, r5, r6
-    character(max_line_len) :: msg  ! error message
+    character(MAX_LINE_LEN) :: msg  ! error message
 
     ! TODO: If there are multiple scattering laws, sample scattering law
 
diff --git a/src/score.f90 b/src/score.f90
index 68c2e6a947..096a249bd9 100644
--- a/src/score.f90
+++ b/src/score.f90
@@ -1,11 +1,12 @@
 module score
 
+  use constants
+  use cross_section, only: get_macro_xs
+  use error,         only: fatal_error
   use global
-  use output, only: message
-  use error,  only: fatal_error
-  use ace,    only: get_macro_xs
-  use search, only: binary_search
-  use string, only: int_to_str
+  use output,        only: message
+  use search,        only: binary_search
+  use string,        only: int_to_str
 
 #ifdef MPI
   use mpi
diff --git a/src/search.f90 b/src/search.f90
index b039a0f2b9..b114b04db4 100644
--- a/src/search.f90
+++ b/src/search.f90
@@ -1,7 +1,7 @@
 module search
 
-  use global
-  use error, only: fatal_error
+  use constants, only: ONE, MAX_LINE_LEN
+  use error,     only: fatal_error
 
 contains
 
@@ -20,7 +20,7 @@ contains
     integer :: L
     integer :: R
     real(8) :: testval
-    character(max_line_len) :: msg
+    character(MAX_LINE_LEN) :: msg
 
     L = 1
     R = n
diff --git a/src/source.f90 b/src/source.f90
index de6f4bb03e..8c69b1565a 100644
--- a/src/source.f90
+++ b/src/source.f90
@@ -1,9 +1,13 @@
 module source
 
+  use bank_header,          only: Bank
+  use constants,            only: ONE, MAX_LINE_LEN
+  use cross_section_header, only: Nuclide
   use global
-  use mcnp_random, only: rang, RN_init_particle
-  use output, only: message
-  use physics, only: watt_spectrum
+  use mcnp_random,          only: rang, RN_init_particle
+  use output,               only: message
+  use particle_header,      only: Particle
+  use physics,              only: watt_spectrum
 
   implicit none
 
@@ -28,7 +32,7 @@ contains
     real(8)    :: E          ! outgoing energy
     real(8)    :: p_min(3)   ! minimum coordinates of source
     real(8)    :: p_max(3)   ! maximum coordinates of source
-    character(max_line_len) :: msg    ! error message
+    character(MAX_LINE_LEN) :: msg    ! error message
 
     msg = 'Initializing source particles...'
     call message(msg, 6)
@@ -135,9 +139,9 @@ contains
 
   subroutine add_bank_sites(p, table, n)
 
-    type(Particle),      pointer    :: p
-    type(AceContinuous), pointer    :: table
-    integer,             intent(in) :: n
+    type(Particle), pointer    :: p
+    type(Nuclide),  pointer    :: table
+    integer,        intent(in) :: n
 
     integer :: i
 
diff --git a/src/source_header.f90 b/src/source_header.f90
new file mode 100644
index 0000000000..077e6338d8
--- /dev/null
+++ b/src/source_header.f90
@@ -0,0 +1,15 @@
+module source_header
+
+  implicit none
+
+!===============================================================================
+! EXTSOURCE describes an external source of neutrons for a fixed-source problem
+! or for the starting source in a criticality problem
+!===============================================================================
+
+  type ExtSource
+     integer :: type                    ! type of source, e.g. 'box' or 'cell'
+     real(8), allocatable :: values(:)  ! values for particular source type
+  end type ExtSource
+
+end module source_header
diff --git a/src/string.f90 b/src/string.f90
index 7c176789c5..3c65cd0732 100644
--- a/src/string.f90
+++ b/src/string.f90
@@ -1,7 +1,7 @@
 module string
 
-  use global
-  use error, only: warning
+  use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL
+  use error,     only: warning
 
   implicit none
 
@@ -24,14 +24,14 @@ contains
   subroutine split_string(string, words, n)
 
     character(*), intent(in)  :: string
-    character(*), intent(out) :: words(max_words)
+    character(*), intent(out) :: words(MAX_WORDS)
     integer,      intent(out) :: n
 
     character(1)  :: char    ! current character
     integer       :: i       ! current index
     integer       :: i_start ! starting index of word
     integer       :: i_end   ! ending index of word
-    character(max_line_len) :: msg
+    character(MAX_LINE_LEN) :: msg
 
     i_start = 0
     i_end = 0
@@ -77,7 +77,7 @@ contains
   subroutine split_string_wl(string, words, n)
 
     character(*), intent(in)  :: string
-    character(*), intent(out) :: words(max_words)
+    character(*), intent(out) :: words(MAX_WORDS)
     integer,      intent(out) :: n
 
     character(1)  :: char    ! current character
@@ -136,7 +136,7 @@ contains
 
     character(*),   intent(in)  :: words(n_words)
     integer,        intent(in)  :: n_words
-    character(max_line_len)     :: string
+    character(MAX_LINE_LEN)     :: string
 
     integer :: i ! index
 
@@ -279,7 +279,7 @@ contains
     integer :: ioError
 
     character(8)            :: fmt ! format for reading string
-    character(max_line_len) :: msg ! error message
+    character(MAX_LINE_LEN) :: msg ! error message
 
     ! Determine total field width
     w = len_trim(string)
@@ -336,7 +336,7 @@ contains
 
     ! Create format specifier for writing character
     num2 = abs(num)
-    if (num2 == ZERO) then
+    if (num2 == 0.0_8) then
        write(fmt, '("(F",I2,".",I2,")")') width, 1
     elseif (num2 < 1.0e-1_8) then
        write(fmt, '("(ES",I2,".",I2,")")') width, decimal - 1
diff --git a/src/tally_header.f90 b/src/tally_header.f90
new file mode 100644
index 0000000000..df49d5e4f9
--- /dev/null
+++ b/src/tally_header.f90
@@ -0,0 +1,39 @@
+module tally_header
+
+  implicit none
+
+!===============================================================================
+! TALLYSCORE
+!===============================================================================
+
+  type TallyScore
+     integer :: n_events
+     real(8) :: val
+     real(8) :: val_sq
+  end type TallyScore
+
+!===============================================================================
+! TALLY
+!===============================================================================
+
+  type Tally
+     integer :: uid
+     integer :: type
+     real(8) :: volume
+     integer :: cell_type
+     integer :: reaction_type
+     integer :: material_type
+     integer, allocatable :: reactions(:)
+     integer, allocatable :: cells(:)
+     integer, allocatable :: materials(:)
+     integer, allocatable :: universes(:)
+     real(8), allocatable :: energies(:)
+     real(8) :: xyz_min(3)
+     real(8) :: xyz_max(3)
+     integer :: n_x
+     integer :: n_y
+     integer :: n_z
+     type(TallyScore), allocatable :: score(:,:,:)
+  end type Tally
+
+end module tally_header
diff --git a/src/timing.f90 b/src/timing.f90
index 90d7e1d79d..92bcb571c8 100644
--- a/src/timing.f90
+++ b/src/timing.f90
@@ -1,33 +1,39 @@
 module timing
 
-  use global
-  use types, only: TimerObj
-  use error, only: warning
+  use constants, only: ZERO
 
   implicit none
 
+!===============================================================================
+! TIMER represents a timer that can be started and stopped to measure how long
+! different routines run. The intrinsic routine system_clock is used to measure
+! time rather than cpu_time.
+!===============================================================================
+
+  type Timer
+     logical :: running      = .false. ! is timer running?
+     integer :: start_counts = 0       ! counts when started
+     real(8) :: elapsed      = 0.      ! total time elapsed in seconds
+   contains
+     procedure :: start     => timer_start
+     procedure :: get_value => timer_get_value
+     procedure :: stop      => timer_stop
+     procedure :: reset     => timer_reset
+  end type Timer
+
 contains
 
 !===============================================================================
 ! TIMER_START
 !===============================================================================
 
-  subroutine timer_start(timer)
+  subroutine timer_start(self)
 
-    type(TimerObj), intent(inout) :: timer
-
-    character(max_line_len) :: msg ! warning message
-
-    ! Check if timer is already running
-    if (timer % running) then
-       msg = "Tried to start a timer that was already running!"
-       call warning(msg)
-       return
-    end if
+    class(Timer), intent(inout) :: self
 
     ! Turn timer on and measure starting time
-    timer % running = .true.
-    call system_clock(timer % start_counts)
+    self % running = .true.
+    call system_clock(self % start_counts)
 
   end subroutine timer_start
 
@@ -35,21 +41,21 @@ contains
 ! TIMER_GET_VALUE
 !===============================================================================
 
-  function timer_get_value(timer) result(elapsed)
+  function timer_get_value(self) result(elapsed)
 
-    type(TimerObj), intent(in) :: timer   ! the timer
-    real(8)                    :: elapsed ! total elapsed time
+    class(Timer), intent(in) :: self   ! the timer
+    real(8)                 :: elapsed ! total elapsed time
 
     integer :: end_counts   ! current number of counts
     integer :: count_rate   ! system-dependent counting rate
     real    :: elapsed_time ! elapsed time since last start
 
-    if (timer % running) then
+    if (self % running) then
        call system_clock(end_counts, count_rate)
-       elapsed_time = real(end_counts - timer % start_counts)/real(count_rate)
-       elapsed = timer % elapsed + elapsed_time
+       elapsed_time = real(end_counts - self % start_counts)/real(count_rate)
+       elapsed = self % elapsed + elapsed_time
     else
-       elapsed = timer % elapsed
+       elapsed = self % elapsed
     end if
 
   end function timer_get_value
@@ -58,22 +64,16 @@ contains
 ! TIMER_STOP
 !===============================================================================
 
-  subroutine timer_stop(timer)
+  subroutine timer_stop(self)
 
-    type(TimerObj), intent(inout) :: timer
-
-    character(max_line_len) :: msg
+    class(Timer), intent(inout) :: self
 
     ! Check to make sure timer was running
-    if (.not. timer % running) then
-       msg = "Tried to stop a timer that was not running!"
-       call warning(msg)
-       return
-    end if
+    if (.not. self % running) return
 
     ! Stop timer and add time
-    timer % elapsed = timer_get_value(timer)
-    timer % running = .false.
+    self % elapsed = self % get_value()
+    self % running = .false.
 
   end subroutine timer_stop
 
@@ -81,13 +81,13 @@ contains
 ! TIMER_RESET
 !===============================================================================
 
-  subroutine timer_reset(timer)
+  subroutine timer_reset(self)
 
-    type(TimerObj), intent(inout) :: timer
+    class(Timer), intent(inout) :: self
 
-    timer % running      = .false.
-    timer % start_counts = 0
-    timer % elapsed      = ZERO
+    self % running      = .false.
+    self % start_counts = 0
+    self % elapsed      = ZERO
 
   end subroutine timer_reset
 
diff --git a/src/types.f90 b/src/types.f90
deleted file mode 100644
index 4b9fc4bc39..0000000000
--- a/src/types.f90
+++ /dev/null
@@ -1,415 +0,0 @@
-module types
-
-  implicit none
-
-!===============================================================================
-! UNIVERSE defines a geometry that fills all phase space
-!===============================================================================
-
-  type Universe
-     integer :: uid                   ! Unique ID
-     integer :: type                  ! Type
-     integer :: level                 ! Level of universe (0=base)
-     integer :: n_cells               ! # of cells within
-     integer, allocatable :: cells(:) ! List of cells within
-     real(8) :: x0                    ! Translation in x-coordinate
-     real(8) :: y0                    ! Translation in y-coordinate
-     real(8) :: z0                    ! Translation in z-coordinate
-     integer, allocatable :: tallies(:)
-  end type Universe
-
-!===============================================================================
-! LATTICE is an ordered array of elements (either rectangular, hexagonal, or
-! triangular)
-!===============================================================================
-
-  type Lattice
-     integer :: uid      ! Universe number for lattice
-     integer :: type     ! Type of lattice (rectangular, hex, etc)
-     integer :: level    ! Level of lattice
-     integer :: n_x      ! number of lattice cells in x-direction
-     integer :: n_y      ! number of lattice cells in y-direction
-     real(8) :: x0       ! x-coordinate of lattice origin
-     real(8) :: y0       ! y-coordinate of lattice origin
-     real(8) :: width_x  ! width of lattice cell 
-     real(8) :: width_y  ! width of lattice cell
-     integer, allocatable :: element(:,:) ! specified universes
-     integer, allocatable :: tallies(:)
-  end type Lattice
-
-!===============================================================================
-! SURFACE type defines a first- or second-order surface that can be used to
-! construct closed volumes (cells)
-!===============================================================================
-
-  type Surface
-     integer :: uid                    ! Unique ID
-     integer :: type                   ! Type of surface
-     real(8), allocatable :: coeffs(:) ! Definition of surface
-     integer, allocatable :: & 
-          & neighbor_pos(:), &         ! List of cells on positive side
-          & neighbor_neg(:)            ! List of cells on negative side
-     integer :: bc                     ! Boundary condition
-  end type Surface
-
-!===============================================================================
-! CELL defines a closed volume by its bounding surfaces
-!===============================================================================
-
-  type Cell
-     integer :: uid        ! Unique ID
-     integer :: type       ! Type of cell (normal, universe, lattice)
-     integer :: universe   ! universe # this cell is in
-     integer :: fill       ! universe # filling this cell
-     integer :: parent     ! cell within which this cell resides
-     integer :: material   ! Material within cell (0 for universe)
-     integer :: n_surfaces ! Number of surfaces within
-     integer, allocatable :: & 
-          & surfaces(:)    ! List of surfaces bounding cell -- note that
-                           ! parentheses, union, etc operators will be listed
-                           ! here too
-     integer, allocatable :: tallies(:)
-  end type Cell
-
-!===============================================================================
-! PARTICLE describes the state of a particle being transported through the
-! geometry
-!===============================================================================
-
-  type Particle
-    integer(8) :: uid          ! Unique ID
-    integer    :: type         ! Particle type (n, p, e, etc)
-    real(8)    :: xyz(3)       ! location
-    real(8)    :: xyz_local(3) ! local location (after transformations)
-    real(8)    :: uvw(3)       ! directional cosines
-    real(8)    :: wgt          ! particle weight
-    real(8)    :: E            ! energy
-    integer    :: IE           ! index on energy grid
-    real(8)    :: interp       ! interpolation factor for energy grid
-    integer    :: cell         ! current cell
-    integer    :: universe     ! current universe
-    integer    :: lattice      ! current lattice
-    integer    :: surface      ! current surface
-    integer    :: index_x      ! lattice index for x direction
-    integer    :: index_y      ! lattice index for y direction
-    logical    :: alive        ! is particle alive?
-    integer    :: n_coll       ! # of collisions
-  end type Particle
-
-!===============================================================================
-! BANK is used for storing fission sites in criticality calculations. Since all
-! the state information of a neutron is not needed, this type allows sites to be
-! stored with less memory
-!===============================================================================
-
-  type Bank
-     integer(8) :: uid    ! Unique ID
-     real(8)    :: xyz(3) ! location of bank particle
-     real(8)    :: uvw(3) ! diretional cosines
-     real(8)    :: E      ! energy
-  end type Bank
-
-!===============================================================================
-! TALLYSCORE
-!===============================================================================
-
-  type TallyScore
-     integer :: n_events
-     real(8) :: val
-     real(8) :: val_sq
-  end type TallyScore
-
-!===============================================================================
-! TALLY
-!===============================================================================
-
-  type Tally
-     integer :: uid
-     integer :: type
-     real(8) :: volume
-     integer :: cell_type
-     integer :: reaction_type
-     integer :: material_type
-     integer, allocatable :: reactions(:)
-     integer, allocatable :: cells(:)
-     integer, allocatable :: materials(:)
-     integer, allocatable :: universes(:)
-     real(8), allocatable :: energies(:)
-     real(8) :: xyz_min(3)
-     real(8) :: xyz_max(3)
-     integer :: n_x
-     integer :: n_y
-     integer :: n_z
-     type(TallyScore), allocatable :: score(:,:,:)
-  end type Tally
-
-!===============================================================================
-! MATERIAL describes a material by its constituent isotopes
-!===============================================================================
-
-  type Material
-     integer :: uid
-     integer :: n_isotopes
-     character(10), allocatable :: names(:)  ! isotope names
-     integer, allocatable :: isotopes(:)     ! index in xsdata list
-     integer, allocatable :: table(:)        ! index in xs array
-     real(8)              :: atom_density    ! total atom density in atom/b-cm
-     real(8), allocatable :: atom_percent(:)
-     real(8), allocatable :: total_xs(:)     ! macroscopic cross-section
-     integer :: sab_table
-  end type Material
-
-!===============================================================================
-! EXTSOURCE describes an external source of neutrons for a fixed-source problem
-! or for the starting source in a criticality problem
-!===============================================================================
-
-  type ExtSource
-     integer :: type                    ! type of source, e.g. 'box' or 'cell'
-     real(8), allocatable :: values(:)  ! values for particular source type
-  end type ExtSource
-
-!===============================================================================
-! ACEDISTANGLE contains data for a tabular secondary angle distribution whether
-! it be tabular or 32 equiprobable cosine bins
-!===============================================================================
-
-  type AceDistAngle
-     integer              :: n_energy    ! # of incoming energies
-     real(8), allocatable :: energy(:)   ! incoming energy grid
-     integer, allocatable :: type(:)     ! type of distribution
-     integer, allocatable :: location(:) ! location of each table
-     real(8), allocatable :: data(:)     ! angular distribution data
-  end type AceDistAngle
-
-!===============================================================================
-! ACEDISTENERGY contains data for a secondary energy distribution for all
-! scattering laws
-!===============================================================================
-
-  type AceDistEnergy
-     integer :: law                    ! secondary distribution law
-     integer :: n_interp               ! # of interpolation regions
-     integer, allocatable :: nbt(:)    ! ENDF interpolation parameters
-     integer, allocatable :: int(:)    ! ''
-     integer :: n_energy               ! # of energies for law validity
-     real(8), allocatable :: energy(:) ! energy grid for law validity
-     real(8), allocatable :: pvalid(:) ! probability of law validity
-     real(8), allocatable :: data(:)   ! energy distribution data
-  end type AceDistEnergy
-
-!===============================================================================
-! ACEREACTION contains the cross-section and secondary energy and angle
-! distributions for a single reaction in a continuous-energy ACE-format table
-!===============================================================================
-
-  type AceReaction
-     integer :: MT                     ! ENDF MT value
-     real(8) :: Q_value                ! Reaction Q value
-     integer :: TY                     ! Number of neutrons released
-     integer :: IE                     ! Starting energy grid index
-     real(8), allocatable :: sigma(:)  ! Cross section values
-     logical :: has_angle_dist         ! Angle distribution present?
-     logical :: has_energy_dist        ! Energy distribution present?
-     type(AceDistAngle)  :: adist      ! Secondary angular distribution
-     type(AceDistEnergy) :: edist      ! Secondary energy distribution
-  end type AceReaction
-
-!===============================================================================
-! ACECONTINUOUS contains all the data for an ACE-format continuous-energy cross
-! section. The ACE format (A Compact ENDF format) is used in MCNP and several
-! other Monte Carlo codes.
-!===============================================================================
-
-  type AceContinuous
-     character(20) :: name
-     real(8) :: awr
-     real(8) :: temp
-     integer :: n_grid
-     integer, allocatable :: grid_index(:)
-     real(8), allocatable :: energy(:)
-     real(8), allocatable :: sigma_t(:)
-     real(8), allocatable :: sigma_a(:)
-     real(8), allocatable :: sigma_el(:)
-     real(8), allocatable :: heating(:)
-
-     ! Total fission neutron emission
-     integer :: nu_t_type
-     real(8), allocatable :: nu_t_data(:)
-
-     ! Prompt fission neutron emission
-     integer :: nu_p_type
-     real(8), allocatable :: nu_p_data(:)
-     
-     ! Delayed fission neutron emission
-     integer :: nu_d_type
-     integer :: n_precursor
-     real(8), allocatable :: nu_d_data(:)
-     real(8), allocatable :: nu_d_precursor_data(:)
-     type(AceDistEnergy), allocatable :: nu_d_edist(:)
-
-     ! Unresolved resonance data
-     logical :: urr_present
-     integer, allocatable :: urr_params(:)
-     real(8), allocatable :: urr_energy(:)
-     real(8), allocatable :: urr_prob(:,:,:)
-
-     ! Reactions
-     integer :: n_reaction
-     type(AceReaction), pointer :: reactions(:) => null()
-
-  end type AceContinuous
-
-!===============================================================================
-! ACETHERMAL contains S(a,b) data for thermal neutron scattering, typically off
-! of light isotopes such as water, graphite, Be, etc
-!===============================================================================
-     
-  type AceThermal
-     character(20) :: name
-     real(8) :: awr
-     real(8) :: temp
-     integer :: n_inelastic_e_in
-     integer :: n_inelastic_e_out
-     integer :: n_inelastic_mu
-     real(8), allocatable :: inelastic_e_in(:)
-     real(8), allocatable :: inelastic_sigma(:) 
-     real(8), allocatable :: inelastic_e_out(:,:)
-     real(8), allocatable :: inelastic_mu(:,:,:)
-     integer :: n_elastic_e_in
-     integer :: n_elastic_type
-     integer :: n_elastic_mu
-     real(8), allocatable :: elastic_e_in(:)
-     real(8), allocatable :: elastic_P(:)
-     real(8), allocatable :: elastic_mu(:,:)
-  end type AceThermal
-
-!===============================================================================
-! XSDATA contains data read in from a SERPENT xsdata file
-!===============================================================================
-
-  type xsData
-     character(10) :: alias
-     character(10) :: id
-     integer :: type
-     integer :: zaid
-     integer :: isomeric
-     real(8) :: awr
-     real(8) :: temp
-     integer :: binary
-     character(150) :: path
-  end type xsData
-
-!===============================================================================
-! TIMEROBJ represents a timer that can be started and stopped to measure how
-! long different routines run. The intrinsic routine system_clock is used to
-! measure time rather than cpu_time.
-!===============================================================================
-
-  type TimerObj
-     logical :: running      = .false. ! is timer running?
-     integer :: start_counts = 0       ! counts when started
-     real(8) :: elapsed      = 0.      ! total time elapsed in seconds
-  end type TimerObj
-
-!===============================================================================
-! KEYVALUECI stores the (key,value) pair for a dictionary where the key is a
-! string and the value is an integer. Note that we need to store the key in
-! addition to the value for collision resolution.
-!===============================================================================
-
-  ! Key length for dictionary
-  integer, parameter :: DICT_KEY_LENGTH = 20
-
-  type KeyValueCI
-     character(len=DICT_KEY_LENGTH) :: key
-     integer                        :: value
-  end type KeyValueCI
-
-!===============================================================================
-! KEYVALUEII stores the (key,value) pair for a dictionary where the key is an
-! integer and the value is an integer. Note that we need to store the key in
-! addition to the value for collision resolution.
-!===============================================================================
-
-  type KeyValueII
-     integer :: key
-     integer :: value
-  end type KeyValueII
-
-!===============================================================================
-! LISTKEYVALUECI stores a linked list of (key,value) pairs where the key is a
-! character and the value is an integer
-!===============================================================================
-
-  type ListKeyValueCI
-     type(ListKeyValueCI), pointer :: next => null()
-     type(KeyValueCI)              :: data
-  end type ListKeyValueCI
-
-!===============================================================================
-! LISTKEYVALUEII stores a linked list of (key,value) pairs where the key is a
-! character and the value is an integer
-!===============================================================================
-
-  type ListKeyValueII
-     type(ListKeyValueII), pointer :: next => null()
-     type(KeyValueII)              :: data
-  end type ListKeyValueII
-
-!===============================================================================
-! LISTREAL stores a linked list of real values. This is used for constructing a
-! unionized energy grid.
-!===============================================================================
-
-  type ListReal
-     type(ListReal), pointer :: next => null()
-     real(8)                 :: data
-  end type ListReal
-
-!===============================================================================
-! LISTINT stores a linked list of integer values.
-!===============================================================================
-
-  type ListInt
-     type(ListInt), pointer :: next => null()
-     integer                :: data
-  end type ListInt
-
-!===============================================================================
-! HASHLISTCI - Since it's not possible to directly do an array of pointers, this
-! derived type provides a pointer
-!===============================================================================
-
-  type HashListCI
-     type(ListKeyValueCI), pointer :: list => null()
-  end type HashListCI
-
-!===============================================================================
-! HASHLISTII - Since it's not possible to directly do an array of pointers, this
-! derived type provides a pointer
-!===============================================================================
-
-  type HashListII
-     type(ListKeyValueII), pointer :: list => null()
-  end type HashListII
-
-!===============================================================================
-! DICTIONARYCI provides a dictionary data structure of (key,value) pairs where
-! the keys are strings and values are integers.
-!===============================================================================
-
-  type DictionaryCI
-     type(HashListCI), pointer :: table(:) => null()
-  end type DictionaryCI
-
-!===============================================================================
-! DICTIONARYII provides a dictionary data structure of (key,value) pairs where
-! the keys and values are both integers.
-!===============================================================================
-
-  type DictionaryII
-     type(HashListII), pointer :: table(:) => null()
-  end type DictionaryII
-
-end module types
diff --git a/src/unittest.f90 b/src/unittest.f90
deleted file mode 100644
index dcfb144288..0000000000
--- a/src/unittest.f90
+++ /dev/null
@@ -1,92 +0,0 @@
-program unittest
-
-  use global
-  use types,       only: Bank
-  use energy_grid, only: add_grid_points
-  use mpi_routines, only: setup_mpi, MPI_BANK
-  use mpi
-
-  call test_mpi_bank()
-
-contains
-
-!===============================================================================
-! TEST_ENERGY_GRID is a test of the energy grid unionization subroutine on a
-! simplified set of data
-!===============================================================================
-
-  subroutine test_energy_grid()
-
-    type(ListReal), pointer :: list => null()
-    type(ListReal), pointer :: current => null()
-    real(8), allocatable :: energy(:)
-    integer :: n, i
-    integer :: n_list
-
-    n = 10
-    allocate(energy(n))
-    do i = 1,n
-       energy(i) = 1.1248_8*i
-    end do
-
-    ! create list
-    allocate(list)
-    current => list
-    current%data = 4.23_8
-    allocate(current%next)
-    current => current%next
-    current%data = 7.23_8
-
-    call add_grid_points(list, energy, n)
-
-    current => list
-    do while(associated(current))
-       print *,current%data
-       current => current%next
-    end do
-
-    deallocate(energy)
-
-  end subroutine test_energy_grid
-
-!===============================================================================
-! TEST_MPI_BANK is a test of the derived datatype MPI_BANK that is used to send
-! fission bank sites during synchronization
-!===============================================================================
-
-  subroutine test_mpi_bank()
-
-    integer :: ierr
-    type(Bank) :: temp_bank(2)
-
-    call setup_mpi()
-
-    if (rank == 0) then
-       temp_bank(1) % uid = 123
-       temp_bank(1) % xyz = (/ 13.0_8, 9.0_8, 8.0_8 /)
-       temp_bank(1) % uvw = (/ 7.5_8, 6.0_8, 2.0_8 /)
-       temp_bank(1) % E   = 4.592834_8
-
-       temp_bank(2) % uid = 321
-       temp_bank(2) % xyz = (/ 23.0_8, 29.0_8, 28.0_8 /)
-       temp_bank(2) % uvw = (/ 1.0_8, 9.0_8, 11.0_8 /)
-       temp_bank(2) % E   = 0.327438_8
-       call MPI_SEND(temp_bank, 2, MPI_BANK, 1, 3, MPI_COMM_WORLD, ierr)
-    elseif (rank == 1) then
-       call MPI_RECV(temp_bank, 2, MPI_BANK, 0, 3, MPI_COMM_WORLD, &
-            & MPI_STATUS_IGNORE, ierr)
-       print *, temp_bank(1) % uid
-       print *, temp_bank(1) % xyz
-       print *, temp_bank(1) % uvw
-       print *, temp_bank(1) % E
-       print *, temp_bank(2) % uid
-       print *, temp_bank(2) % xyz
-       print *, temp_bank(2) % uvw
-       print *, temp_bank(2) % E
-    end if
-
-    call MPI_FINALIZE(ierr)
-    
-  end subroutine test_mpi_bank
-
-end program unittest