merged develop into cmfd_rbgs branch

.gitignore

@@ -1,4 +1,5 @@
 # Compiled objects and modules
+*.a
 *.o
 *.mod
 *.log
@@ -18,6 +19,7 @@ src/openmc
 
 # Documentation builds
 docs/build
+docs/source/_images/*.pdf
 
 # xml-fortran reader
 src/xml-fortran/xmlreader
@@ -26,4 +28,4 @@ src/xml-fortran/xmlreader
 src/templates/*.f90
 
 # Test results error file
-results_error.dat
+results_error.dat

LICENSE

@@ -1,4 +1,4 @@
-Copyright (c) 2011-2013 Massachusetts Institute of Technology
+Copyright (c) 2011-2014 Massachusetts Institute of Technology
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

docs/Makefile

@@ -6,13 +6,19 @@ SPHINXOPTS    =
 SPHINXBUILD   = sphinx-build
 PAPER         =
 BUILDDIR      = build
+IMAGEDIR      = source/_images
 
 # Internal variables.
 PAPEROPT_a4     = -D latex_paper_size=a4
 PAPEROPT_letter = -D latex_paper_size=letter
 ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
 
-.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest
+# SVG to PDF conversion
+SVG2PDF = inkscape
+PDFS    = $(patsubst %.svg,%.pdf,$(wildcard $(IMAGEDIR)/*.svg))
+
+
+.PHONY: help images clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest
 
 help:
 	@echo "Please use \`make <target>' where <target> is one of"
@@ -33,8 +39,16 @@ help:
 	@echo "  linkcheck  to check all external links for integrity"
 	@echo "  doctest    to run all doctests embedded in the documentation (if enabled)"
 
+# Pattern rule for converting SVG to PDF
+%.pdf: %.svg
+	$(SVG2PDF) -f $< -A $@
+
+# Rule to build PDFs
+images: $(PDFS)
+
 clean:
 	-rm -rf $(BUILDDIR)/*
+	-rm $(PDFS)
 
 html:
 	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
@@ -91,14 +105,14 @@ epub:
 	@echo
 	@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
 
-latex:
+latex: images
 	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
 	@echo
 	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
 	@echo "Run \`make' in that directory to run these through (pdf)latex" \
 	      "(use \`make latexpdf' here to do that automatically)."
 
-latexpdf:
+latexpdf: images
 	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
 	@echo "Running LaTeX files through pdflatex..."
 	make -C $(BUILDDIR)/latex all-pdf

docs/source/conf.py

@@ -39,7 +39,7 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'OpenMC'
-copyright = u'2011-2013, Massachusetts Institute of Technology'
+copyright = u'2011-2014, Massachusetts Institute of Technology'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -48,7 +48,7 @@ copyright = u'2011-2013, Massachusetts Institute of Technology'
 # The short X.Y version.
 version = "0.5"
 # The full version, including alpha/beta/rc tags.
-release = "0.5.2"
+release = "0.5.3"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -121,7 +121,7 @@ html_title = "OpenMC Documentation"
 
 # The name of an image file (relative to this directory) to place at the top
 # of the sidebar.
-html_logo = '../img/openmc.png'
+html_logo = '_images/openmc.png'
 
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
@@ -188,6 +188,8 @@ latex_documents = [
    u'Massachusetts Institute of Technology', 'manual'),
 ]
 
+latex_elements = {'preamble': '\\usepackage{enumitem}\\setlistdepth{9}'}
+
 # The name of an image file (relative to this directory) to place at the top of
 # the title page.
 #latex_logo = None

docs/source/devguide/index.rst

@@ -15,6 +15,6 @@ as debugging.
     structures
     styleguide
     workflow
-    xml-fortran
+    xml-parsing
     statepoint
     voxel

docs/source/devguide/statepoint.rst

@@ -4,6 +4,544 @@
 State Point Binary File Specifications
 ======================================
 
+-----------
+Revision 10 
+-----------
+
+**integer(4) FILETYPE_STATEPOINT**
+
+    Flags whether this file is a statepoint file or a particle restart file.
+
+**integer(4) REVISION_STATEPOINT**
+
+    Revision of the binary state point file. Any time a change is made in the
+    format of the state-point file, this integer is incremented.
+
+**integer(4) VERSION_MAJOR**
+
+    Major version number for OpenMC
+
+**integer(4) VERSION_MINOR**
+
+    Minor version number for OpenMC
+
+**integer(4) VERSION_RELEASE**
+
+    Release version number for OpenMC
+
+**character(19) time_stamp**
+
+    Date and time the state point was written.
+
+**character(255) path**
+
+    Absolute path to directory containing input files.
+
+**integer(8) seed**
+
+    Pseudo-random number generator seed.
+
+**integer(4) run_mode**
+
+    run mode used. The modes are described in constants.F90.
+
+**integer(8) n_particles**
+
+    Number of particles used per generation.
+
+**integer(4) n_batches**
+
+    Total number of batches (active + inactive).
+
+**integer(4) current_batch**
+
+    The number of batches already simulated.
+
+if (run_mode == MODE_EIGENVALUE)
+
+    **integer(4) n_inactive**
+
+        Number of inactive batches
+
+    **integer(4) gen_per_batch**
+
+        Number of generations per batch for criticality calculations
+
+    *do i = 1, current_batch \* gen_per_batch*
+
+        **real(8) k_generation(i)**
+
+             k-effective for the i-th total generation
+
+    *do i = 1, current_batch \* gen_per_batch*
+
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th total generation
+
+    **real(8) k_col_abs**
+
+        Sum of product of collision/absorption estimates of k-effective
+
+    **real(8) k_col_tra**
+
+        Sum of product of collision/track-length estimates of k-effective
+
+    **real(8) k_abs_tra**
+
+        Sum of product of absorption/track-length estimates of k-effective
+
+    **real(8) k_combined(2)**
+
+        Mean and standard deviation of a combined estimate of k-effective
+
+    **integer(4) cmfd_on**
+
+        Flag that cmfd is on
+
+    if (cmfd_on)
+
+        **integer(4) cmfd % indices**
+
+            Indices for cmfd mesh (i,j,k,g)
+
+        **real(8) cmfd % k_cmfd(1:current_batch)**
+
+            CMFD eigenvalues
+
+        **real(8) cmfd % src(1:I,1:J,1:K,1:G)**
+
+            CMFD fission source
+
+        **real(8) cmfd % entropy(1:current_batch)**
+
+            CMFD estimate of Shannon entropy
+
+        **real(8) cmfd % balance(1:current_batch)**
+
+            RMS of the residual neutron balance equation on CMFD mesh
+
+        **real(8) cmfd % dom(1:current_batch)**
+
+            CMFD estimate of dominance ratio
+
+        **real(8) cmfd % scr_cmp(1:current_batch)**
+
+            RMS comparison of difference between OpenMC and CMFD fission source
+
+**integer(4) n_meshes**
+
+    Number of meshes in tallies.xml file
+
+*do i = 1, n_meshes*
+
+    **integer(4) meshes(i) % id**
+
+        Unique ID of mesh.
+
+    **integer(4) meshes(i) % type**
+
+        Type of mesh.
+
+    **integer(4) meshes(i) % n_dimension**
+
+        Number of dimensions for mesh (2 or 3).
+
+    **integer(4) meshes(i) % dimension(:)**
+
+        Number of mesh cells in each dimension.
+
+    **real(8) meshes(i) % lower_left(:)**
+
+        Coordinates of lower-left corner of mesh.
+
+    **real(8) meshes(i) % upper_right(:)**
+
+        Coordinates of upper-right corner of mesh.
+
+    **real(8) meshes(i) % width(:)**
+
+        Width of each mesh cell in each dimension.
+
+**integer(4) n_tallies**
+
+*do i = 1, n_tallies*
+
+    **integer(4) tallies(i) % id**
+
+        Unique ID of tally.
+
+    **integer(4) tallies(i) % n_realizations**
+
+        Number of realizations for the i-th tally.
+
+    **integer(4) size(tallies(i) % scores, 1)**
+
+        Total number of score bins for the i-th tally
+
+    **integer(4) size(tallies(i) % scores, 2)**
+
+        Total number of filter bins for the i-th tally
+
+    **integer(4) tallies(i) % n_filters**
+
+    *do j = 1, tallies(i) % n_filters*
+
+        **integer(4) tallies(i) % filter(j) % type**
+
+            Type of tally filter.
+
+        **integer(4) tallies(i) % filter(j) % n_bins**
+
+            Number of bins for filter.
+
+        **integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
+
+            Value for each filter bin of this type.
+
+    **integer(4) tallies(i) % n_nuclide_bins**
+
+        Number of nuclide bins. If none are specified, this is just one.
+
+    *do j = 1, tallies(i) % n_nuclide_bins*
+
+        **integer(4) tallies(i) % nuclide_bins(j)**
+
+            Values of specified nuclide bins
+
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins.
+
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % score_bins(j)**
+
+            Values of specified scoring bins (e.g. SCORE_FLUX).
+            
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % scatt_order(j)**
+
+            Scattering Order specified scoring bins.
+            
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins without accounting for those added by
+        the scatter-pn command.
+
+**integer(4) n_realizations**
+
+    Number of realizations for global tallies.
+
+**integer(4) N_GLOBAL_TALLIES**
+
+    Number of global tally scores
+
+*do i = 1, N_GLOBAL_TALLIES*
+
+    **real(8) global_tallies(i) % sum**
+
+        Accumulated sum for the i-th global tally
+
+    **real(8) global_tallies(i) % sum_sq**
+
+        Accumulated sum of squares for the i-th global tally
+
+**integer(4) tallies_on**
+
+    Flag indicated if tallies are present in the file.
+
+if (tallies_on > 0)
+
+    *do i = 1, n_tallies*
+
+        *do k = 1, size(tallies(i) % scores, 2)*
+
+            *do j = 1, size(tallies(i) % scores, 1)*
+
+                **real(8) tallies(i) % scores(j,k) % sum**
+            
+                    Accumulated sum for the j-th score and k-th filter of the
+                    i-th tally
+
+                **real(8) tallies(i) % scores(j,k) % sum_sq**
+
+                    Accumulated sum of squares for the j-th score and k-th
+                    filter of the i-th tally
+
+if (run_mode == MODE_EIGENVALUE)
+
+    *do i = 1, n_particles*
+
+        **real(8) source_bank(i) % wgt**
+
+            Weight of the i-th source particle
+
+        **real(8) source_bank(i) % xyz(1:3)**
+
+            Coordinates of the i-th source particle.
+
+        **real(8) source_bank(i) % uvw(1:3)**
+
+            Direction of the i-th source particle
+
+        **real(8) source_bank(i) % E**
+
+            Energy of the i-th source particle.
+
+----------
+Revision 9
+----------
+
+**integer(4) FILETYPE_STATEPOINT**
+
+    Flags whether this file is a statepoint file or a particle restart file.
+
+**integer(4) REVISION_STATEPOINT**
+
+    Revision of the binary state point file. Any time a change is made in the
+    format of the state-point file, this integer is incremented.
+
+**integer(4) VERSION_MAJOR**
+
+    Major version number for OpenMC
+
+**integer(4) VERSION_MINOR**
+
+    Minor version number for OpenMC
+
+**integer(4) VERSION_RELEASE**
+
+    Release version number for OpenMC
+
+**character(19) time_stamp**
+
+    Date and time the state point was written.
+
+**character(255) path**
+
+    Absolute path to directory containing input files.
+
+**integer(8) seed**
+
+    Pseudo-random number generator seed.
+
+**integer(4) run_mode**
+
+    run mode used. The modes are described in constants.F90.
+
+**integer(8) n_particles**
+
+    Number of particles used per generation.
+
+**integer(4) n_batches**
+
+    Total number of batches (active + inactive).
+
+**integer(4) current_batch**
+
+    The number of batches already simulated.
+
+if (run_mode == MODE_EIGENVALUE)
+
+    **integer(4) n_inactive**
+
+        Number of inactive batches
+
+    **integer(4) gen_per_batch**
+
+        Number of generations per batch for criticality calculations
+
+    *do i = 1, current_batch \* gen_per_batch*
+
+        **real(8) k_generation(i)**
+
+             k-effective for the i-th total generation
+
+    *do i = 1, current_batch \* gen_per_batch*
+
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th total generation
+
+    **real(8) k_col_abs**
+
+        Sum of product of collision/absorption estimates of k-effective
+
+    **real(8) k_col_tra**
+
+        Sum of product of collision/track-length estimates of k-effective
+
+    **real(8) k_abs_tra**
+
+        Sum of product of absorption/track-length estimates of k-effective
+
+    **real(8) k_combined(2)**
+
+        Mean and standard deviation of a combined estimate of k-effective
+
+**integer(4) n_meshes**
+
+    Number of meshes in tallies.xml file
+
+*do i = 1, n_meshes*
+
+    **integer(4) meshes(i) % id**
+
+        Unique ID of mesh.
+
+    **integer(4) meshes(i) % type**
+
+        Type of mesh.
+
+    **integer(4) meshes(i) % n_dimension**
+
+        Number of dimensions for mesh (2 or 3).
+
+    **integer(4) meshes(i) % dimension(:)**
+
+        Number of mesh cells in each dimension.
+
+    **real(8) meshes(i) % lower_left(:)**
+
+        Coordinates of lower-left corner of mesh.
+
+    **real(8) meshes(i) % upper_right(:)**
+
+        Coordinates of upper-right corner of mesh.
+
+    **real(8) meshes(i) % width(:)**
+
+        Width of each mesh cell in each dimension.
+
+**integer(4) n_tallies**
+
+*do i = 1, n_tallies*
+
+    **integer(4) tallies(i) % id**
+
+        Unique ID of tally.
+
+    **integer(4) tallies(i) % n_realizations**
+
+        Number of realizations for the i-th tally.
+
+    **integer(4) size(tallies(i) % scores, 1)**
+
+        Total number of score bins for the i-th tally
+
+    **integer(4) size(tallies(i) % scores, 2)**
+
+        Total number of filter bins for the i-th tally
+
+    **integer(4) tallies(i) % n_filters**
+
+    *do j = 1, tallies(i) % n_filters*
+
+        **integer(4) tallies(i) % filter(j) % type**
+
+            Type of tally filter.
+
+        **integer(4) tallies(i) % filter(j) % n_bins**
+
+            Number of bins for filter.
+
+        **integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
+
+            Value for each filter bin of this type.
+
+    **integer(4) tallies(i) % n_nuclide_bins**
+
+        Number of nuclide bins. If none are specified, this is just one.
+
+    *do j = 1, tallies(i) % n_nuclide_bins*
+
+        **integer(4) tallies(i) % nuclide_bins(j)**
+
+            Values of specified nuclide bins
+
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins.
+
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % score_bins(j)**
+
+            Values of specified scoring bins (e.g. SCORE_FLUX).
+            
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % scatt_order(j)**
+
+            Scattering Order specified scoring bins.
+            
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins without accounting for those added by
+        the scatter-pn command.
+
+**integer(4) n_realizations**
+
+    Number of realizations for global tallies.
+
+**integer(4) N_GLOBAL_TALLIES**
+
+    Number of global tally scores
+
+*do i = 1, N_GLOBAL_TALLIES*
+
+    **real(8) global_tallies(i) % sum**
+
+        Accumulated sum for the i-th global tally
+
+    **real(8) global_tallies(i) % sum_sq**
+
+        Accumulated sum of squares for the i-th global tally
+
+**integer(4) tallies_on**
+
+    Flag indicated if tallies are present in the file.
+
+if (tallies_on > 0)
+
+    *do i = 1, n_tallies*
+
+        *do k = 1, size(tallies(i) % scores, 2)*
+
+            *do j = 1, size(tallies(i) % scores, 1)*
+
+                **real(8) tallies(i) % scores(j,k) % sum**
+            
+                    Accumulated sum for the j-th score and k-th filter of the
+                    i-th tally
+
+                **real(8) tallies(i) % scores(j,k) % sum_sq**
+
+                    Accumulated sum of squares for the j-th score and k-th
+                    filter of the i-th tally
+
+if (run_mode == MODE_EIGENVALUE)
+
+    *do i = 1, n_particles*
+
+        **real(8) source_bank(i) % wgt**
+
+            Weight of the i-th source particle
+
+        **real(8) source_bank(i) % xyz(1:3)**
+
+            Coordinates of the i-th source particle.
+
+        **real(8) source_bank(i) % uvw(1:3)**
+
+            Direction of the i-th source particle
+
+        **real(8) source_bank(i) % E**
+
+            Energy of the i-th source particle.
+
 ----------
 Revision 8
 ----------

docs/source/devguide/styleguide.rst

@@ -35,17 +35,17 @@ Don't use ``print *`` or ``write(*,*)``. If writing to a file, use a specific
 unit. Writing to standard output or standard error should be handled by the
 ``write_message`` subroutine or functionality in the error module.
 
--------------------------
-Subroutines and Functions
--------------------------
+----------
+Procedures
+----------
 
-Above each subroutine/function, include a comment block giving a brief
-description of what the subroutine or function does.
+Above each procedure, include a comment block giving a brief description of what
+the procedure does.
 
-Arguments to subroutines/functions should be explicitly specified as intent(in),
-intent(out), or intent(inout).
+Nonpointer dummy arguments to procedures should be explicitly specified as
+intent(in), intent(out), or intent(inout).
 
-Include a comment describing what each argument to a subroutine/function is.
+Include a comment describing what each argument to a procedure is.
 
 ---------
 Variables
@@ -133,9 +133,11 @@ value of f90-continuation-indent in Emacs.
 Whitespace in Expressions
 -------------------------
 
+Use a single space between arguments to procedures. 
+
 Avoid extraneous whitespace in the following situations:
 
-- In subroutine/function calls::
+- In procedure calls::
 
     Yes: call somesub(x, y(2), z)
     No:  call somesub( x, y( 2 ), z )

docs/source/devguide/workflow.rst

@@ -5,41 +5,113 @@ Development Workflow
 ====================
 
 Anyone wishing to make contributions to OpenMC should be fully acquianted and
-comfortable working with git_ and GitHub_. The primary means of modifying and
-making contributions to OpenMC is through GitHub `pull requests`_. This is
-what's known as a fork and pull development model. The steps for this are as
-follows:
+comfortable working with git_ and GitHub_. We assume here that you have git
+installed on your system, have a GitHub account, and have setup SSH keys to be
+able to create/push to repositories on GitHub. 
+
+Overview
+--------
+
+Development of OpenMC relies heavily on branching; specifically, we use a
+branching model sometimes referred to as `git flow`_. If you plan to contribute
+to OpenMC development, we highly recommend that you read the linked blog post to
+get a sense of how the branching model works. There are two main branches that
+always exist: *master* and *develop*. The *master* branch is a stable branch
+that contains the latest release of the code. The *develop* branch is where any
+ongoing development takes place prior to a release and is not guaranteed to be
+stable. When the development team decides that a release should occur, the
+*develop* branch is merged into *master*.
+
+Trivial changes to the code may be committed directly to the *develop* branch by
+a trusted developer. However, most new features should be developed on a branch
+that branches off of *develop*. When the feature is completed, a `pull request`_
+is initiated on GitHub that is then reviewed by a trusted developer. If the pull
+request is satisfactory, it is then merged into *develop*. Note that a trusted
+developer may not review their own pull request (i.e., an independent code
+review is required).
+
+Code Review Criteria
+--------------------
+
+In order to be considered suitable for inclusion in the *develop* branch, the
+following criteria must be satisfied for all proposed changes:
+
+- Changes have a clear purpose and are useful.
+- Compiles under all conditions (MPI, OpenMP, HDF5, etc.).  This is checked as
+  part of the test suite (see `test_compile.py`_).
+- Passes the regression suite.
+- If appropriate, test cases are added to regression suite.
+- No memory leaks (checked with valgrind_).
+- Conforms to the OpenMC `style guide`_.
+- No degradation of performance or greatly increased memory usage. This is not a
+  hard rule -- in certain circumstances, a performance loss might be acceptable
+  if there are compelling reasons.
+- New features/input are documented.
+- No unnecessary external software dependencies are introduced.
+
+Contributing
+------------
+
+Now that you understand the basic development workflow, let's discuss how an
+individual to contribute to development. Note that this would apply to both new
+features and bug fixes. The general steps for contributing are as follows:
 
 1. Fork the main openmc repository from `mit-crpg/openmc`_. This will create a
-repository with the same name under your personal account. As such, you can
-commit to it as you please without disrupting other developers.
+   repository with the same name under your personal account. As such, you can
+   commit to it as you please without disrupting other developers.
 
-2. Create a branch that you want merged back to `mit-crpg/openmc`_ and make
-commits that you intend to go back. If you have made other changes that should
-not be merged back, those changes should be on another branch.
+   .. image:: ../_images/fork.png
 
-3. Issue a pull request from GitHub and select the branch you want merged.
+2. Clone your fork of OpenMC and create a branch that branches off of *develop*:
 
-4. The OpenMC integration manager will review your pull request and make sure it
-conforms to the :ref:`devguide_styleguide`, compiles correctly, runs in parallel
-correctly, does not break other features in the code, etc. Any issues with the
-pull request can be discussed directly on the pull request page itself.
+   .. code-block:: sh
 
-5. After the pull request has been thoroughly vetted, it is merged back into
-`mit-crpg/openmc`_.
+       git clone git@github.com:yourusername/openmc.git
+       cd openmc
+       git checkout -b newbranch develop
+
+3. Make your changes on the new branch that you intend to have included in
+   *develop*. If you have made other changes that should not be merged back, 
+   ensure that those changes are made on a different branch.
+
+4. Issue a pull request from GitHub and select the *develop* branch of
+   mit-crpg/openmc as the target.
+
+   .. image:: ../_images/pullrequest.png
+
+   At a minimum, you should describe what the changes you've made are and why
+   you are making them. If the changes are related to an oustanding issue, make
+   sure it is cross-referenced. A wise developer would also check whether their
+   changes do indeed pass the regression test suite.
+
+5. A trusted developer will review your pull request based on the criteria
+   above. Any issues with the pull request can be discussed directly on the pull
+   request page itself.
+
+6. After the pull request has been thoroughly vetted, it is merged back into the
+   *develop* branch of mit-crpg/openmc.
+
+Private Development
+-------------------
 
 While the process above depends on the fork of the OpenMC repository being
 publicly available on GitHub, you may also wish to do development on a private
 repository for research or commercial purposes. The proper way to do this is to
 create a complete copy of the OpenMC repository (not a fork from GitHub). The
 private repository can then either be stored just locally or in conjunction with
-a private repository on Github (this requires a `paid plan`_). If you want to
-merge some changes you've made in your private repository back to
-`mit-crpg/openmc`_ repository, simply follow the steps above with an extra step
-of pulling a branch from your private repository into your public fork.
+a private repository on Github (this requires a `paid plan`_). Alternatively,
+`Bitbucket`_ offers private repositories for free. If you want to merge some
+changes you've made in your private repository back to mit-crpg/openmc
+repository, simply follow the steps above with an extra step of pulling a branch
+from your private repository into a public fork.
 
 .. _git: http://git-scm.com/
 .. _GitHub: https://github.com/
-.. _pull requests: https://help.github.com/articles/using-pull-requests
+.. _git flow: http://nvie.com/git-model
+.. _test_compile.py: https://github.com/mit-crpg/openmc/blob/develop/tests/test_compile/test_compile.py
+.. _valgrind: http://valgrind.org/
+.. _style guide: http://mit-crpg.github.io/openmc/devguide/styleguide.html
+.. _pull request: https://help.github.com/articles/using-pull-requests
 .. _mit-crpg/openmc: https://github.com/mit-crpg/openmc
 .. _paid plan: https://github.com/plans
+.. _Bitbucket: https://bitbucket.org

docs/source/devguide/xml-fortran.rst

@@ -1,40 +0,0 @@
-.. _devguide_xml-fortran:
-
-=========================
-xml-fortran Input Parsing
-=========================
-
-OpenMC relies on the xml-fortran package for reading and intrepreting the XML
-input files for geometry, materials, settings, tallies, etc. The use of an XML
-format makes writing input files considerably more flexible than would otherwise
-be possible.
-
-With the xml-fortran package, extending the user input files to include new tags
-is fairly straightforward. A "template" file exists for each diferent type of
-input file that tells xml-fortran what to expect in a file. These template files
-can be found in the src/templates directory. The steps for modifying/adding
-input are as follows:
-
-1. Add a ``<variable>``` tag to the desired template file,
-e.g. src/templates/geometry_t.xml. See the `xml-fortran documentation`_ for a
-description of the acceptable fields.
-
-2. In the input_xml module, any input given in your new tag will be read
-automatically through a call to, e.g. read_xml_file_geometry_t. Whatever
-variable name you specified should have the data available.
-
-3. Add code in the appropriate subroutine to check the variable for any possible
-errors.
-
-4. Add a variable in OpenMC to copy the temporary variable into if there are no
-errors.
-
-A set of `RELAX NG`_ schemata exists that enables real-time validation of input
-files when using the GNU Emacs text editor. You should also modify the RELAX NG
-schema for the template you changed (e.g. src/templates/geometry.rnc) so that
-those who use Emacs can confirm whether their input is valid before they
-run. You will need to be familiar with RELAX NG `compact syntax`_.
-
-.. _xml-fortran documentation: http://xml-fortran.sourceforge.net/documentation.html
-.. _RELAX NG: http://relaxng.org/
-.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html

docs/source/devguide/xml-parsing.rst

@@ -0,0 +1,38 @@
+.. _devguide_xml-parsing:
+
+=================
+XML Input Parsing
+=================
+
+OpenMC relies on the FoX_ Fortran XML library for reading and intrepreting the
+XML input files for geometry, materials, settings, tallies, etc. The use of an
+XML format makes writing input files considerably more flexible than would
+otherwise be possible.
+
+With the FoX library, extending the user input files to include new tags is
+fairly straightforward. The steps for modifying/adding input are as follows:
+
+1. Add appropriate calls to procedures from the `xml_interface module`_, such as
+``check_for_node``, ``get_node_value``, and ``get_node_array``. All input
+reading is performed in the `input_xml module`_.
+
+2. Make sure that your input can be categorized as one of the datatypes from
+`XML Schema Part 2`_ and that parsing of the data appropriately reflects
+this. For example, for a boolean_ value, true can be represented either by "true"
+or by "1".
+
+3. Add code to check the variable for any possible errors.
+
+A set of `RELAX NG`_ schemata exists that enables real-time validation of input
+files when using the GNU Emacs text editor. You should also modify the RELAX NG
+schema for the file you changed (e.g. src/relaxng/geometry.rnc) so that
+those who use Emacs can confirm whether their input is valid before they
+run. You will need to be familiar with RELAX NG `compact syntax`_.
+
+.. _FoX: https://github.com/andreww/fox
+.. _xml_interface module: https://github.com/mit-crpg/openmc/blob/develop/src/xml_interface.F90
+.. _input_xml module: https://github.com/mit-crpg/openmc/blob/develop/src/input_xml.F90
+.. _XML Schema Part 2: http://www.w3.org/TR/xmlschema-2/
+.. _boolean: http://www.w3.org/TR/xmlschema-2/#boolean
+.. _RELAX NG: http://relaxng.org/
+.. _compact syntax: http://relaxng.org/compact-tutorial-20030326.html

docs/source/license.rst

@@ -4,7 +4,7 @@
 License Agreement
 =================
 
-Copyright © 2011-2013 Massachusetts Institute of Technology
+Copyright © 2011-2014 Massachusetts Institute of Technology
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

docs/source/methods/cross_sections.rst

@@ -80,7 +80,7 @@ dashed box would need to be stored on a per-nuclide basis, and the union grid
 would need to be stored once. This method is also referred to as *double
 indexing* and is available as an option in Serpent (see paper by Leppanen_).
 
-.. figure:: ../../img/uniongrid.svg
+.. figure:: ../_images/uniongrid.*
     :width: 600px
     :align: center
     :figclass: align-center

docs/source/methods/eigenvalue.rst

@@ -108,6 +108,40 @@ at plots of :math:`k_{eff}` and the Shannon entropy. A number of methods have
 been proposed (see e.g. [Romano]_, [Ueki]_), but each of these is not without
 problems.
 
+---------------------------
+Uniform Fission Site Method
+---------------------------
+
+Generally speaking, the variance of a Monte Carlo tally will be inversely
+proportional to the number of events that score to the tally. In a reactor
+problem, this implies that regions with low relative power density will have
+higher variance that regions with high relative power density. One method to
+circumvent the uneven distribution of relative errors is the uniform fission
+site (UFS) method introduced by [Sutton]_. In this method, the portion of the
+problem containing fissionable material is subdivided into a number of cells
+(typically using a structured mesh). Rather than producing
+
+.. math::
+
+    m = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t}
+
+fission sites at each collision where :math:`w` is the weight of the neutron,
+:math:`k` is the previous-generation estimate of the neutron multiplication
+factor, :math:`\nu\Sigma_f` is the neutron production cross section, and
+:math:`\Sigma_t` is the total cross section, in the UFS method we produce
+
+.. math::
+
+    m_{UFS} = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} \frac{v_i}{s_i}
+
+fission sites at each collision where :math:`v_i` is the fraction of the total
+volume occupied by cell :math:`i` and :math:`s_i` is the fraction of the fission
+source contained in cell :math:`i`. To ensure that no bias is introduced, the
+weight of each fission site stored in the fission bank is :math:`s_i/v_i` rather
+than unity. By ensuring that the expected number of fission sites in each mesh
+cell is constant, the collision density across all cells, and hence the variance
+of tallies, is more uniform than it would be otherwise.
+
 .. _Shannon entropy: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-06-3737_entropy.pdf
 
 .. [Lieberoth] J. Lieberoth, "A Monte Carlo Technique to Solve the Static
@@ -119,5 +153,9 @@ problems.
    *Proc. International Conference on Mathematics, Computational Methods, and
    Reactor Physics*, Saratoga Springs, New York (2009).
 
+.. [Sutton] Daniel J. Kelly, Thomas M. Sutton, and Stephen C. Wilson, "MC21
+   Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark
+   Problem," *Proc. PHYSOR 2012*, Knoxville, Tennessee, Apr. 15--20 (2012).
+
 .. [Ueki] Taro Ueki, "On-the-Fly Judgments of Monte Carlo Fission Source
    Convergence," *Trans. Am. Nucl. Soc.*, **98**, 512 (2008).

docs/source/methods/geometry.rst

@@ -45,7 +45,7 @@ surface by a combination of the unique ID of the surface and a positive/negative
 sign. The following illustration shows an example of an ellipse with unique ID 1
 dividing space into two half-spaces.
 
-.. figure:: ../../img/halfspace.svg
+.. figure:: ../_images/halfspace.*
    :align: center
    :figclass: align-center
 
@@ -60,7 +60,7 @@ half-space references whose intersection defines the region. The region is then
 assigned a material defined elsewhere. The following illustration shows an
 example of a cell defined as the intersection of an ellipse and two planes.
    
-.. figure:: ../../img/union.svg
+.. figure:: ../_images/union.*
    :align: center
    :figclass: align-center
 

docs/source/methods/parallelization.rst

@@ -65,7 +65,7 @@ in the case of an eigenvalue calculation). This idea is illustrated in
 
 .. _figure-master-slave:
 
-.. figure:: ../../img/master-slave.png
+.. figure:: ../_images/master-slave.png
     :align: center
     :figclass: align-center
 
@@ -122,7 +122,7 @@ needed. This concept is illustrated in :ref:`Figure 2
 
 .. _figure-nearest-neighbor:
 
-.. figure:: ../../img/nearest-neighbor.png
+.. figure:: ../_images/nearest-neighbor.png
     :align: center
     :figclass: align-center
 
@@ -203,7 +203,7 @@ communicated between adjacent nodes.
 
 .. _figure-neighbor-example:
 
-.. figure:: ../../img/nearest-neighbor-example.png
+.. figure:: ../_images/nearest-neighbor-example.png
     :align: center
     :figclass: align-center
 

docs/source/methods/physics.rst

@@ -790,6 +790,7 @@ outgoing angle is
 
     \mu = \frac{1}{A} \ln \left ( \xi_4 e^A + (1 - \xi_4) e^{-A} \right ).
 
+.. _ace-law-61:
 
 ACE Law 61 - Correlated Energy and Angle Distribution
 +++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -952,7 +953,7 @@ as
 
     v_n \bar{\sigma} (v_n, T) = \int d\mathbf{v}_T v_r \sigma(v_r)
     M (\mathbf{v}_T)
-    
+
 where :math:`v_n` is the magnitude of the velocity of the neutron,
 :math:`\bar{\sigma}` is an effective cross section, :math:`T` is the temperature
 of the target material, :math:`\mathbf{v}_T` is the velocity of the target
@@ -1321,7 +1322,7 @@ given analytically by
 
     \mu = 1 - \frac{E_i}{E}
 
-where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. 
+where :math:`E_i` is the energy of the Bragg edge that scattered the neutron.
 
 Outgoing Angle for Incoherent Elastic Scattering
 ------------------------------------------------
@@ -1348,18 +1349,24 @@ where the interpolation factor is defined as
 Outgoing Energy and Angle for Inelastic Scattering
 --------------------------------------------------
 
-On each |sab| table, there is a correlated angle-energy secondary distribution
-for neutron thermal inelastic scattering. While the documentation for the ACE
-format implies that there are a series of equiprobable outgoing energies, the
-outgoing energies may have non-uniform probability distribution. In particular,
-if the thermal data were processed with :math:`iwt = 0` in NJOY, then the first
-and last outgoing energies have a relative probability of 1, the second and
-second to last energies have a relative probability of 4, and all other energies
-have a relative probability of 10. The procedure to determine the outgoing
-energy and angle is as such. First, the interpolation factor is determined from
-equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is sampled
-either from a uniform distribution or from the aforementioned skewed
-distribution. The outgoing energy is then interpolated between values
+Each |sab| table provides a correlated angle-energy secondary distribution for
+neutron thermal inelastic scattering.  There are three representations used
+in the ACE thermal scattering data: equiprobable discrete outgoing
+energies, non-uniform yet still discrete outgoing energies, and continuous
+outgoing energies with corresponding probability and cumulative distribution
+functions provided in tabular format.  These three representations all
+represent the angular distribution in a common format, using a series of
+discrete equiprobable outgoing cosines.
+
+Equi-Probable Outgoing Energies
++++++++++++++++++++++++++++++++
+
+If the thermal data was processed with :math:`iwt = 1` in NJOY, then the
+outgoing energy spectra is represented in the ACE data as a set of discrete and
+equiprobable outgoing energies.  The procedure to determine the outgoing energy
+and angle is as such. First, the interpolation factor is determined from
+equation :eq:`sab-interpolation-factor`.  Then, an outgoing energy bin is
+sampled from a uniform distribution and then interpolated between values
 corresponding to neighboring incoming energies:
 
 .. math::
@@ -1380,6 +1387,37 @@ uniformly and then the final cosine is interpolated on the incoming energy grid:
 where :math:`\mu_{i,j,k}` is the k-th outgoing cosine corresponding to the j-th
 outgoing energy and the i-th incoming energy.
 
+Skewed Equi-Probable Outgoing Energies
+++++++++++++++++++++++++++++++++++++++
+
+If the thermal data was processed with :math:`iwt=0` in NJOY, then the
+outgoing energy spectra is represented in the ACE data according to the
+following: the first and last outgoing energies have a relative probability of
+1, the second and second-to-last energies have a relative probability of 4, and
+all other energies have a relative probability of 10.  The procedure to
+determine the outgoing energy and angle is similar to the method discussed
+above, except that the sampled probability distribution is now skewed
+accordingly.
+
+Continuous Outgoing Energies
+++++++++++++++++++++++++++++
+
+If the thermal data was processed with :math:`iwt=2` in NJOY, then the
+outgoing energy spectra is represented by a continuous outgoing energy spectra
+in tabular form with linear-linear interpolation.  The sampling of the outgoing
+energy portion of this format is very similar to :ref:`ACE Law 61<ace-law-61>`,
+but the sampling of the correlated angle is performed as it was in the other
+two representations discussed in this sub-section.  In the Law 61 algorithm,
+we found an interpolation factor :math:`f`, statistically sampled an incoming
+energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on
+the tabulated cumulative distribution function. Once the outgoing energy has
+been determined with equation :eq:`ace-law-4-energy`, we then need to decide
+which angular distribution data to use.  Like the linear-linear interpolation
+case in Law 61, the angular distribution closest to the sampled value of the
+cumulative distribution function for the outgoing energy is utilized.  The
+actual algorithm utilized to sample the outgoing angle is shown in equation
+:eq:`inelastic-angle`.
+
 .. _probability_tables:
 
 ----------------------------------------------

docs/source/publications.rst

@@ -4,6 +4,53 @@
 Publications
 ============
 
+- Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte
+  Carlo simulations using the multipole representation," *Ann. Nucl. Energy*,
+  **64**, 78--85 (2014). `<http://dx.doi.org/10.1016/j.anucene.2013.09.043>`_
+
+- Andrew Siegel, Kord Smith, Kyle Felker, Paul Romano, Benoit Forget, and Peter
+  Beckman, "Improved cache performance in Monte Carlo transport calculations
+  using energy banding," *Comput. Phys. Commun.*
+  (2013). `<http://dx.doi.org/10.1016/j.cpc.2013.10.008>`_
+
+- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Validation of OpenMC
+  Reactor Physics Simulations with the B&W 1810 Series Benchmarks,"
+  *Trans. Am. Nucl. Soc.*, **109**, 1301--1304 (2013).
+
+- Bryan R. Herman, Benoit Forget, and Kord Smith, "Utilizing CMFD in OpenMC to
+  Estimate Dominance Ratio and Adjoint," *Trans. Am. Nucl. Soc.*, **109**,
+  1389-1392 (2013).
+
+- Timothy P. Burke, Brian C. Kiedrowski, and William R. Martin, "Flux and
+  Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*,
+  **109**, 683-686 (2013).
+
+- Paul K. Romano, Benoit Forget, Kord Smith, and Andrew Siegel, "On the use of
+  tally servers in Monte Carlo simulations of light-water reactors,"
+  *Proc. Joint International Conference on Supercomputing in Nuclear
+  Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013).
+
+- Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit
+  Forget, and Kord Smith, "OpenMC: A State-of-the-Art Monte Carlo Code for
+  Research and Development," *Proc. Joint International Conference on
+  Supercomputing in Nuclear Applications and Monte Carlo*, Paris, France,
+  Oct. 27--31 (2013).
+
+- Kyle G. Felker, Andrew R. Siegel, Kord S. Smith, Paul K. Romano, and Benoit
+  Forget, "The energy band memory server algorithm for parallel Monte Carlo
+  calculations," *Proc. Joint International Conference on Supercomputing in
+  Nuclear Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013).
+
+- John R. Tramm and Andrew R. Siegel, "Memory Bottlenecks and Memory Contention
+  in Multi-Core Monte Carlo Transport Codes," *Proc. Joint International
+  Conference on Supercomputing in Nuclear Applications and Monte Carlo*, Paris,
+  France, Oct. 27--31 (2013).
+
+- Andrew R. Siegel, Kord Smith, Paul K. Romano, Benoit Forget, and Kyle Felker,
+  "Multi-core performance studies of a Monte Carlo neutron transport code,"
+  *Int. J. High Perform. Comput. Appl.*, **28** (1), 87--96
+  (2014). `<http://dx.doi.org/10.1177/1094342013492179>`_
+
 - Paul K. Romano, Andrew R. Siegel, Benoit Forget, and Kord Smith, "Data
   decomposition of Monte Carlo particle transport simulations via tally
   servers," *J. Comput. Phys.*, **252**, 20--36

docs/source/releasenotes/index.rst

@@ -10,6 +10,7 @@ bugs fixed, and known issues for each successive release.
 .. toctree::
     :maxdepth: 1
 
+    notes_0.5.3
     notes_0.5.2
     notes_0.5.1
     notes_0.5.0

docs/source/releasenotes/notes_0.5.3.rst

@@ -4,10 +4,6 @@
 Release Notes for OpenMC 0.5.3
 ==============================
 
-.. note::
-    These release notes are for an upcoming release of OpenMC and are still
-    subject to change.
-
 -------------------
 System Requirements
 -------------------
@@ -21,12 +17,20 @@ the problem at hand (mostly on the number of nuclides in the problem).
 New Features
 ------------
 
-- Special run mode --tallies removed.
+- Output interface enhanced to allow multiple files handles to be opened
+- Particle restart file linked to output interface
 - Particle restarts and state point restarts are both identified with the -r
   command line flag.
-- New regression test suite.
-- All memory leaks fixed.
-- Shared-memory parallelism with OpenMP.
+- Particle instance no longer global, passed to all physics routines
+- Physics routines refactored to rely less on global memory, more arguments
+  passed in
+- CMFD routines refactored and now can compute dominance ratio on the fly
+- PETSc 3.4.2 or higher must be used and compiled with fortran datatype support
+- Memory leaks fixed except for ones from xml-fortran package
+- Test suite enhanced to test output with different compiler options
+- Description of OpenMC development workflow added
+- OpenMP shared-memory parallelism added
+- Special run mode --tallies removed.
 
 ---------
 Bug Fixes
@@ -35,7 +39,7 @@ Bug Fixes
 - 2b1e8a_: Normalize direction vector after reflecting particle.
 - 5853d2_: Set blank default for cross section listing alias.
 - e178c7_: Fix infinite loop with words greater than 80 characters in write_message.
-- c18a6e_: Chcek for valid secondary mode on S(a,b) tables.
+- c18a6e_: Check for valid secondary mode on S(a,b) tables.
 - 82c456_: Fix bug where last process could have zero particles.
 
 .. _2b1e8a: https://github.com/mit-crpg/openmc/commit/2b1e8a

docs/source/usersguide/input.rst

@@ -198,19 +198,26 @@ tally data, this option can significantly improve the parallel efficiency.
 --------------------
 
 The ``<output>`` element determines what output files should be written to disk
-during the run. This element has no attributes or sub-elements and should be set
-to a list of strings separated by spaces. Valid options are "summary",
-"cross-sections", and "tallies". For example, if you want the summary and cross
-sections summary file to be written, this element should be given as:
+during the run. The sub-elements are described below, where "true" will write
+out the file and "false" will not.
 
-  .. code-block:: xml
+  :cross_sections:
+    Writes out an ASCII summary file of the cross sections that were read in.
 
-      <output>summary cross_sections</output>
+    *Default*: false
 
-  .. note:: The tally results will be written to a binary/HDF5 state point file by
-            default.
+  :summary: 
+    Writes out an ASCII summary file describing all of the user input files that
+    were read in.
 
-  *Default*: "tallies"
+    *Default*: false
+
+  :tallies:
+    Write out an ASCII file of tally results.
+
+    *Default*: true
+
+  .. note:: The tally results will always be written to a binary/HDF5 state point file.
 
 ``<output_path>`` Element
 -------------------------
@@ -398,6 +405,15 @@ integers: the batch number, generation number, and particle number.
 
   *Default*: None
 
+.. _track:
+
+``<track>`` Element
+-------------------
+
+The ``<track>`` element specifies particles for which OpenMC will output binary files describing particle position at every step of its transport. This element should be followed by triplets of integers.  Each triplet describes one particle.  The integers in each triplet specify the batch number, generation number, and particle number, respectively.
+
+  *Default*: None
+
 ``<uniform_fs>`` Element
 ------------------------
 
@@ -636,9 +652,10 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
 ---------------------
 
 The ``<lattice>`` can be used to represent repeating structures (e.g. fuel pins
-in an assembly) or other geometry which naturally fits into a two-dimensional
-structured mesh. Each cell within the lattice is filled with a specified
-universe. A ``<lattice>`` accepts the following attributes or sub-elements:
+in an assembly) or other geometry which naturally fits into a two- or
+three-dimensional structured mesh. Each cell within the lattice is filled with a
+specified universe. A ``<lattice>`` accepts the following attributes or
+sub-elements:
 
   :id:
     A unique integer that can be used to identify the surface.
@@ -650,18 +667,19 @@ universe. A ``<lattice>`` accepts the following attributes or sub-elements:
     *Default*: rectangular
 
   :dimension:
-    Two integers representing the number of lattice cells in the x- and y-
-    directions, respectively.
+    Two or three integers representing the number of lattice cells in the x- and
+    y- (and z-) directions, respectively.
 
     *Default*: None
 
   :lower_left:
-    The coordinates of the lower-left corner of the lattice.
+    The coordinates of the lower-left corner of the lattice. If the lattice is
+    two-dimensional, only the x- and y-coordinates are specified.
 
     *Default*: None
 
   :width:
-    The width of the lattice cell in the x- and y- directions.
+    The width of the lattice cell in the x- and y- (and z-) directions.
 
     *Default*: None
 
@@ -1196,13 +1214,9 @@ with "false".
 ``<mesh>`` Element
 ------------------
 
-If a structured mesh is desired as a filter for a tally, it must be specified in
-a separate element with the tag name ``<mesh>``. This element has the following
+The CMFD mesh is a structured Cartesian mesh. This element has the following
 attributes/sub-elements:
 
-  :type:
-    The type of structured mesh. Only "rectangular" is currently supported.
-
   :lower_left:
     The lower-left corner of the structured mesh. If only two coordinate are
     given, it is assumed that the mesh is an x-y mesh.

docs/source/usersguide/install.rst

@@ -71,13 +71,15 @@ Prerequisites
       To compile with support for HDF5_ output (highly recommended), you will
       need to have HDF5 installed on your computer. The installed version will
       need to have been compiled with the same compiler you intend to compile
-      OpenMC with. HDF5_ must be built with parallel I/O features. An example
-      of configuring HDF5_ is listed below::
+      OpenMC with. HDF5_ must be built with parallel I/O features if you intend
+      to use HDF5_ with MPI. An example of configuring HDF5_ is listed below::
 
            FC=/opt/mpich/3.0.4-gnu/bin/mpif90 CC=/opt/mpich/3.0.4-gnu/bin/mpicc \
            ./configure --prefix=/opt/hdf5/1.8.11-gnu --enable-fortran \
                        --enable-fortran2003 --enable-parallel
 
+      You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
+
     * PETSc_ for CMFD acceleration
 
       To enable CMFD acceleration, you will need to have PETSc_ (3.4.2 or higher)
@@ -91,7 +93,7 @@ Prerequisites
                        --with-fortran-datatypes
 
       The BLAS/LAPACK library is not required to be downloaded and can be linked
-      explicitly (e.g., Intel MLK library).
+      explicitly (e.g., Intel MKL library).
 
     * git_ version control software for obtaining source code
 
@@ -356,6 +358,7 @@ OpenMC accepts the following command line flags:
 -r, --restart file     Restart a previous run from a state point or a particle
                        restart file
 -s, --threads N        Run with *N* OpenMP threads
+-t, --track            Write tracks for all particles
 -v, --version          Show version information
 
 -----------------------------------------------------

docs/source/usersguide/processing.rst

@@ -38,7 +38,7 @@ running OpenMC with the -plot or -p command-line option (See
 Plotting in 2D
 --------------
 
-.. image:: ../../img/atr.png
+.. image:: ../_images/atr.png
    :height: 200px
 
 After running OpenMC to obtain PPM files, images should be saved to another
@@ -58,7 +58,7 @@ Ubuntu: ``sudo apt-get install imagemagick``).  Images are then converted like:
 Plotting in 3D
 --------------
 
-.. image:: ../../img/3dgeomplot.png
+.. image:: ../_images/3dgeomplot.png
    :height: 200px
 
 The binary VOXEL files output by OpenMC can not be viewed directly by any
@@ -162,7 +162,7 @@ tasks will be described here in the following sections.
 Plotting in 2D
 --------------
 
-.. image:: ../../img/plotmeshtally.png
+.. image:: ../_images/plotmeshtally.png
    :height: 200px
 
 For simple viewing of 2D slices of a mesh plot, the utility plot_mesh_tally.py
@@ -170,7 +170,7 @@ is provided.  This utility provides an interactive GUI to explore and plot
 mesh tallies for any scores and filter bins.  It requires statepoint.py, as well
 as `PyQt <http://www.riverbankcomputing.com/software/pyqt>`_.
 
-.. image:: ../../img/fluxplot.png
+.. image:: ../_images/fluxplot.png
    :height: 200px
 
 Alternatively, the user can write their own Python script to manipulate the data
@@ -249,7 +249,7 @@ two heatmaps in the previous figure.
 Plotting in 3D
 --------------
 
-.. image:: ../../img/3dcore.png
+.. image:: ../_images/3dcore.png
    :height: 200px
 
 As with 3D plots of the geometry, meshtally data needs to be put into a standard
@@ -353,9 +353,89 @@ file. Note that the data contained in the output from
 ``StatePoint.extract_result`` is already in a Numpy array that can be reshaped
 and dumped to MATLAB in one step.
 
+----------------------------
+Particle Track Visualization
+----------------------------
 
+.. image:: ../_images/Tracks.png
+   :height: 200px
 
+OpenMC can dump particle tracks—the position of particles as they are
+transported through the geometry.  There are two ways to make OpenMC output
+tracks: all particle tracks through a commandline argument or specific particle
+tracks through settings.xml.
 
+Running OpenMC with the argument "-t", "-track", or "--track" will cause a track
+file to be created for every particle transported in the code.
 
+The settings.xml file can dictate that specific particle tracks are output.
+These particles are specified withen a ''track'' element.  The ''track'' element
+should contain triplets of integers specifying the batch, generation, and
+particle numbers, respectively.  For example, to output the tracks for particles
+3 and 4 of batch 1 and generation 2 the settings.xml file should contain:
 
+.. code-block:: xml
 
+      <track>
+        1 2 3
+        1 2 4
+      </track>
+
+After running OpenMC, the directory should contain a file of the form
+"track_(batch #)_(generation #)_(particle #).(binary or h5)" for each particle
+tracked.  These track files can be converted into VTK poly data files with the
+"track.py" utility.  The usage of track.py is of the form "track.py [-o OUT] IN"
+where OUT is the optional output filename and IN is one or more filenames
+describing track files.  The default output name is "track.pvtp".  A common
+usage of track.py is "track.py track*.binary" which will use the data from all
+binary track files in the directory to write a "track.pvtp" VTK output file.
+The .pvtp file can then be read and plotted by 3d visualization programs such as
+Paraview.
+
+----------------------
+Source Site Processing
+----------------------
+
+For eigenvalue problems, OpenMC will store information on the fission source
+sites in the statepoint file by default. For each source site, the weight,
+position, sampled direction, and sampled energy are stored. To extract this data
+from a statepoint file, the statepoint.py Python module can be used. Below is an
+example of an interactive ipython session using the statepoint.py Python module:
+
+.. code-block:: python
+
+    In [1]: import statepoint
+    
+    In [2]: sp = statepoint.StatePoint('statepoint.100.h5')
+    
+    In [3]: sp.read_source()
+    
+    In [4]: len(sp.source)
+    Out[4]: 1000
+    
+    In [5]: sp.source[0:10]
+    Out[5]: 
+    [<SourceSite: xyz=[  2.21980946  -8.92686048  87.93720485] at E=0.932923263566>,
+     <SourceSite: xyz=[  2.21980946  -8.92686048  87.93720485] at E=0.349240220512>,
+     <SourceSite: xyz=[-31.21542213 -30.26762771  72.10845757] at E=3.75843584486>,
+     <SourceSite: xyz=[-31.21542213 -30.26762771  72.10845757] at E=0.80550137267>,
+     <SourceSite: xyz=[   0.18805099  -69.13376508  103.67726838] at E=1.67922461097>,
+     <SourceSite: xyz=[   0.18805099  -69.13376508  103.67726838] at E=1.16304110199>,
+     <SourceSite: xyz=[ -50.42189115   -9.96571672  123.34077905] at E=0.710937974074>,
+     <SourceSite: xyz=[ -32.80427668  -15.49316628  125.26301151] at E=1.61907104162>,
+     <SourceSite: xyz=[  53.20376026  -15.38643708  120.58071044] at E=3.33962024907>,
+     <SourceSite: xyz=[  53.20376026  -15.38643708  120.58071044] at E=1.90185680329>]
+    
+    In [6]: site = sp.source[0]
+    
+    In [7]: site.weight
+    Out[7]: 1.0
+    
+    In [8]: site.xyz
+    Out[8]: array([  2.21980946,  -8.92686048,  87.93720485])
+    
+    In [9]: site.uvw
+    Out[9]: array([ 0.06740523,  0.50612814,  0.85982024])
+    
+    In [10]: site.E
+    Out[10]: 0.93292326356564159

docs/source/usersguide/troubleshoot.rst

@@ -19,25 +19,6 @@ you are using a compiler that does not support type-bound procedures from
 Fortran 2003. This affects any version of gfortran prior to 4.6. Downloading and
 installing the latest gfortran_ compiler should resolve this problem.
 
-Fatal Error: Wrong module version '4' (expected '9') for file 'xml_data_cmfd_t.mod' opened at (1)
-*************************************************************************************************
-
-The `.mod` modules files that are created by gfortran are versioned and
-sometimes are usually not backwards compatible. If gfortran is upgraded and the
-modules files for xml-fortran source files are not deleted, this error may
-occur. To fix this, clear out all module and object files with :program:`make
-distclean` and then recompiling.
-
-Fatal Error: File 'xml_data_cmfd_t.mod' opened at (1) is not a GFORTRAN module file
-***********************************************************************************
-
-When OpenMC compiles, the first thing it needs to do is compile source in the
-xml-fortran subdirectory. If you compiled everything with a compiler other than
-gfortran, performed a :program:`make clean`, and then tried to :program:`make`
-with gfortran, the xml-fortran modules would have been compiled with a different
-compiler. To fix this, try clearing out all module and object files with
-:program:`make distclean` and then recompiling.
-
 gfortran: unrecognized option '-cpp'
 ************************************
 

man/man1/openmc.1

@@ -27,6 +27,9 @@ Restart a previous run from a state point or a particle restart file named
 .BI \-s " N" "\fR,\fP \-\-threads" " N"
 Use \fIN\fP OpenMP threads.
 .TP
+.B "\-t\fR, \fP\-\-track"
+Write tracks for all particles.
+.TP
 .B "\-v\fR, \fP\-\-version"
 Show version information.
 .TP
@@ -43,7 +46,7 @@ to locate ACE format cross section libraries if the user has not specified the
 <cross_sections> tag in
 .I settings.xml\fP.
 .SH LICENSE
-Copyright \(co 2011-2013 Massachusetts Institute of Technology.
+Copyright \(co 2011-2014 Massachusetts Institute of Technology.
 .PP
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

schemas.xml

@@ -1,10 +1,10 @@
 <?xml version="1.0"?>
 <locatingRules xmlns="http://thaiopensource.com/ns/locating-rules/1.0">
-  <documentElement localName="geometry" uri="src/templates/geometry.rnc"/>
-  <documentElement localName="materials" uri="src/templates/materials.rnc"/>
-  <documentElement localName="settings" uri="src/templates/settings.rnc"/>
-  <documentElement localName="tallies" uri="src/templates/tallies.rnc"/>
-  <documentElement localName="plots" uri="src/templates/plots.rnc"/>
-  <documentElement localName="cmfd" uri="src/templates/cmfd.rnc"/>
-  <documentElement localName="cross_sections" uri="src/templates/cross_sections.rnc"/>
+  <documentElement localName="geometry" uri="src/relaxng/geometry.rnc"/>
+  <documentElement localName="materials" uri="src/relaxng/materials.rnc"/>
+  <documentElement localName="settings" uri="src/relaxng/settings.rnc"/>
+  <documentElement localName="tallies" uri="src/relaxng/tallies.rnc"/>
+  <documentElement localName="plots" uri="src/relaxng/plots.rnc"/>
+  <documentElement localName="cmfd" uri="src/relaxng/cmfd.rnc"/>
+  <documentElement localName="cross_sections" uri="src/relaxng/cross_sections.rnc"/>
 </locatingRules>

src/DEPENDENCIES

@@ -37,6 +37,7 @@ cmfd_execute.o: tally.o
 cmfd_header.o: constants.o
 
 cmfd_input.o: cmfd_header.o
+cmfd_input.o: constants.o
 cmfd_input.o: error.o
 cmfd_input.o: global.o
 cmfd_input.o: mesh_header.o
@@ -45,7 +46,7 @@ cmfd_input.o: string.o
 cmfd_input.o: tally.o
 cmfd_input.o: tally_header.o
 cmfd_input.o: tally_initialize.o
-cmfd_input.o: templates/cmfd_t.o
+cmfd_input.o: xml_interface.o
 
 cmfd_jfnk_solver.o: cmfd_loss_operator.o
 cmfd_jfnk_solver.o: cmfd_power_solver.o
@@ -218,12 +219,7 @@ input_xml.o: random_lcg.o
 input_xml.o: string.o
 input_xml.o: tally_header.o
 input_xml.o: tally_initialize.o
-input_xml.o: templates/cross_sections_t.o
-input_xml.o: templates/geometry_t.o
-input_xml.o: templates/materials_t.o
-input_xml.o: templates/plots_t.o
-input_xml.o: templates/settings_t.o
-input_xml.o: templates/tallies_t.o
+input_xml.o: xml_interface.o
 
 interpolation.o: constants.o
 interpolation.o: endf_header.o
@@ -341,6 +337,7 @@ solver_interface.o: vector_header.o
 source.o: bank_header.o
 source.o: constants.o
 source.o: error.o
+source.o: geometry.o
 source.o: geometry_header.o
 source.o: global.o
 source.o: math.o
@@ -382,6 +379,11 @@ tally_initialize.o: tally_header.o
 
 timer_header.o: constants.o
 
+track_output.o: global.o
+track_output.o: output_interface.o
+track_output.o: particle_header.o
+track_output.o: string.o
+
 tracking.o: cross_section.o
 tracking.o: error.o
 tracking.o: geometry.o
@@ -393,6 +395,10 @@ tracking.o: physics.o
 tracking.o: random_lcg.o
 tracking.o: string.o
 tracking.o: tally.o
+tracking.o: track_output.o
 
 vector_header.o: constants.o
 
+xml_interface.o: constants.o
+xml_interface.o: error.o
+xml_interface.o: global.o

src/Makefile

@@ -3,10 +3,7 @@ prefix = /usr/local
 
 source = $(wildcard *.F90)
 objects = $(source:.F90=.o)
-templates = $(wildcard templates/*.o)
-xml_fort = xml-fortran/xmlparse.o \
-           xml-fortran/read_xml_primitives.o \
-           xml-fortran/write_xml_primitives.o
+xml_lib = -Lxml/lib -lxml
 
 #===============================================================================
 # User Options
@@ -184,17 +181,17 @@ endif
 
 ifeq ($(OPENMP),yes)
   ifeq ($(COMPILER),intel)
-    F90FLAGS += -openmp -DOPENMP
+    F90FLAGS += -openmp
     LDFLAGS += -openmp
   endif
 
   ifeq ($(COMPILER),gnu)
-    F90FLAGS += -fopenmp -DOPENMP
+    F90FLAGS += -fopenmp
     LDFLAGS += -fopenmp
   endif
 
   ifeq ($(COMPILER),ibm)
-    F90FLAGS += -qsmp=omp -WF,-DOPENMP
+    F90FLAGS += -qsmp=omp
     LDFLAGS += -qsmp=omp
   endif
 endif
@@ -221,7 +218,7 @@ endif
 
 ifeq ($(MACHINE),bluegene)
   F90 = /bgsys/drivers/ppcfloor/comm/xl/bin/mpixlf2003
-  F90FLAGS = -WF,-DNO_F2008,-DMPI -O3
+  F90FLAGS = -WF,-DNO_F2008,-DMPI,-DRESTRICTED_ASSOCIATED_BUG -O3
   LDFLAGS = -lmpich.cnkf90
 endif
 
@@ -236,19 +233,18 @@ endif
 
 ifeq ($(MACHINE),bluegeneq)
   F90 = mpixlf2003
-  F90FLAGS = -WF,-DNO_F2008,-DMPI -O5
+  F90FLAGS = -WF,-DNO_F2008,-DMPI,-DRESTRICTED_ASSOCIATED_BUG -O5
 endif
 
 #===============================================================================
 # Targets
 #===============================================================================
 
-all: xml-fortran $(program)
-xml-fortran:
-	cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)"
-	cd templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)"
+all: xml $(program)
+xml:
+	cd xml; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)"
 $(program): $(objects)
-	$(F90) $(objects) $(templates) $(xml_fort) $(LDFLAGS) -o $@
+	$(F90) $(objects) $(xml_lib) $(LDFLAGS) -o $@
 install:
 	@install -D $(program) $(DESTDIR)$(prefix)/bin/$(program)
 	@install -D utils/statepoint_cmp.py $(DESTDIR)$(prefix)/bin/statepoint_cmp
@@ -264,8 +260,7 @@ uninstall:
 	@rm $(DESTDIR)$(prefix)/share/man/man1/openmc.1
 	@rm $(DESTDIR)$(prefix)/share/doc/$(program)/copyright
 distclean: clean
-	cd xml-fortran; make clean
-	cd templates; make clean
+	cd xml; make clean
 clean:
 	@rm -f *.o *.mod $(program)
 neat:
@@ -275,10 +270,11 @@ neat:
 # Rules
 #===============================================================================
 
-.PHONY: all xml-fortran install uninstall clean neat distclean 
+.SUFFIXES: .F90 .o
+.PHONY: all xml install uninstall clean neat distclean 
 
 %.o: %.F90
-	$(F90) $(F90FLAGS) -DGIT_SHA1="\"$(GIT_SHA1)\"" -Ixml-fortran -Itemplates -c $<
+	$(F90) $(F90FLAGS) -DGIT_SHA1="\"$(GIT_SHA1)\"" -Ixml/include -c $<
 
 #===============================================================================
 # Dependencies

src/ace.F90

@@ -115,9 +115,9 @@ contains
         ! search through the list of nuclides for one which has a matching zaid
         sab => sab_tables(mat % i_sab_tables(k))
 
-        ! Loop through nuclides and find  match
+        ! Loop through nuclides and find match
         FIND_NUCLIDE: do j = 1, mat % n_nuclides
-          if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then
+          if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then
             mat % i_sab_nuclides(k) = j
             exit FIND_NUCLIDE
           end if
@@ -161,16 +161,16 @@ contains
           mat % i_sab_tables(m)   = temp_table
         end do SORT_SAB
       end if
-      
+
       ! Deallocate temporary arrays for names of nuclides and S(a,b) tables
       if (allocated(mat % names)) deallocate(mat % names)
       if (allocated(mat % sab_names)) deallocate(mat % sab_names)
 
     end do MATERIAL_LOOP2
-    
+
     ! Avoid some valgrind leak errors
     call already_read % clear()
-    
+
   end subroutine read_xs
 
 !===============================================================================
@@ -244,8 +244,16 @@ contains
       ! Read first line of header
       read(UNIT=in, FMT='(A10,2G12.0,1X,A10)') name, awr, kT, date_
 
+      ! Check that correct xs was found -- if cross_sections.xml is broken, the
+      ! location of the table may be wrong
+      if(adjustl(name) /= adjustl(listing % name)) then
+        message = "XS listing entry " // trim(listing % name) // " did not &
+             &match ACE data, " // trim(name) // " found instead."
+        call fatal_error()
+      end if
+
       ! Read more header and NXS and JXS
-      read(UNIT=in, FMT=100) comment, mat, & 
+      read(UNIT=in, FMT=100) comment, mat, &
            (zaids(i), awrs(i), i=1,16), NXS, JXS
 100   format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/&
            ,8I9/8I9/8I9/8I9/8I9/8I9)
@@ -269,7 +277,7 @@ contains
            ACCESS='direct', RECL=record_length)
 
       ! Read all header information
-      read(UNIT=in, REC=location) name, awr, kT, date_, & 
+      read(UNIT=in, REC=location) name, awr, kT, date_, &
            comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS
 
       ! determine table length
@@ -325,7 +333,15 @@ contains
       sab % name = name
       sab % awr  = awr
       sab % kT   = kT
-      sab % zaid = zaids(1)
+      ! Find sab % n_zaid
+      do i = 1, 16
+        if (zaids(i) == 0) then
+          sab % n_zaid = i - 1
+          exit
+        end if
+      end do
+      allocate(sab % zaid(sab % n_zaid))
+      sab % zaid = zaids(1: sab % n_zaid)
 
       call read_thermal_data(sab)
     end select
@@ -398,7 +414,7 @@ contains
     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 :: 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
@@ -801,7 +817,7 @@ contains
 
     LED  = JXS(10)
 
-    ! Loop over all reactions 
+    ! Loop over all reactions
     do i = 1, NXS(5)
       rxn => nuc % reactions(i+1) ! skip over elastic scattering
       rxn % has_energy_dist = .true.
@@ -832,7 +848,7 @@ contains
 
     integer :: LDIS   ! location of all energy distributions
     integer :: LNW    ! location of next energy distribution if multiple
-    integer :: LAW    ! secondary energy distribution law   
+    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
@@ -1179,15 +1195,10 @@ contains
     integer :: NE_out ! number of outgoing energies
     integer :: NMU    ! number of outgoing angles
     integer :: JXS4   ! location of elastic energy table
+    integer(8), allocatable :: LOCC(:) ! Location of inelastic data
 
-    ! read secondary energy mode for inelastic scattering and check
+    ! read secondary energy mode for inelastic scattering
     table % secondary_mode = NXS(7)
-    if (table % secondary_mode /= SAB_SECONDARY_EQUAL .and. &
-         table % secondary_mode /= SAB_SECONDARY_SKEWED) then
-      message = "Unsupported secondary mode on S(a,b) table " // &
-           trim(adjustl(table % name)) // ": " // to_str(table % secondary_mode)
-      call fatal_error()
-    end if
 
     ! read number of inelastic energies and allocate arrays
     NE_in = int(XSS(JXS(1)))
@@ -1205,29 +1216,67 @@ contains
 
     ! allocate space for outgoing energy/angle for inelastic
     ! scattering
-    NE_out = NXS(4)
-    NMU = NXS(3) + 1
-    table % n_inelastic_e_out = NE_out
-    table % n_inelastic_mu = NMU
-    allocate(table % inelastic_e_out(NE_out, NE_in))
-    allocate(table % inelastic_mu(NMU, NE_out, NE_in))
+    if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. &
+         table % secondary_mode == SAB_SECONDARY_SKEWED) then
+      NMU = NXS(3) + 1
+      table % n_inelastic_mu = NMU
+      NE_out = NXS(4)
+      table % n_inelastic_e_out = NE_out
+      allocate(table % inelastic_e_out(NE_out, NE_in))
+      allocate(table % inelastic_mu(NMU, NE_out, NE_in))
+    else if (table % secondary_mode == SAB_SECONDARY_CONT) then
+      NMU = NXS(3) - 1
+      table % n_inelastic_mu = NMU
+      allocate(table % inelastic_data(NE_in))
+      allocate(LOCC(NE_in))
+      ! NE_out will be determined later
+    end if
 
     ! read outgoing energy/angle distribution for inelastic scattering
-    lc = JXS(3) - 1
-    do i = 1, NE_in
-      do j = 1, NE_out
-        ! read outgoing energy
-        table % inelastic_e_out(j,i) = XSS(lc + 1)
+    if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. &
+         table % secondary_mode == SAB_SECONDARY_SKEWED) then
+      lc = JXS(3) - 1
+      do i = 1, NE_in
+        do j = 1, NE_out
+          ! read outgoing energy
+          table % inelastic_e_out(j,i) = XSS(lc + 1)
 
-        ! read outgoing angles for this outgoing energy
-        do k = 1, NMU
-          table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
+          ! read outgoing angles for this outgoing energy
+          do k = 1, NMU
+            table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
+          end do
+
+          ! advance pointer
+          lc = lc + 1 + NMU
         end do
-
-        ! advance pointer
-        lc = lc + 1 + NMU
       end do
-    end do
+    else if (table % secondary_mode == SAB_SECONDARY_CONT) then
+      ! Get the location pointers to each Ein's DistEnergySAB data
+      LOCC = get_int(NE_in)
+      ! Get the number of outgoing energies and allocate space accordingly
+      do i = 1, NE_in
+        NE_out = int(XSS(XSS_index + i - 1))
+        table % inelastic_data(i) % n_e_out = NE_out
+        allocate(table % inelastic_data(i) % e_out (NE_out))
+        allocate(table % inelastic_data(i) % e_out_pdf (NE_out))
+        allocate(table % inelastic_data(i) % e_out_cdf (NE_out))
+        allocate(table % inelastic_data(i) % mu (NMU, NE_out))
+      end do
+
+      ! Now we can fill the inelastic_data(i) attributes
+      do i = 1, NE_in
+        XSS_index = LOCC(i)
+        NE_out = table % inelastic_data(i) % n_e_out
+        do j = 1, NE_out
+          table % inelastic_data(i) % e_out(j) = XSS(XSS_index + 1)
+          table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2)
+          table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3)
+          table % inelastic_data(i) % mu(:, j) = &
+            XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
+          XSS_index = XSS_index + 4 + NMU - 1
+        end do
+      end do
+    end if
 
     ! read number of elastic energies and allocate arrays
     JXS4 = JXS(4)

src/ace_header.F90

@@ -16,7 +16,7 @@ module ace_header
     integer, allocatable :: type(:)     ! type of distribution
     integer, allocatable :: location(:) ! location of each table
     real(8), allocatable :: data(:)     ! angular distribution data
-    
+
     ! Type-Bound procedures
     contains
       procedure :: clear => distangle_clear ! Deallocates DistAngle
@@ -32,10 +32,10 @@ module ace_header
     type(Tab1) :: p_valid             ! probability of law validity
     real(8), allocatable :: data(:)   ! energy distribution data
 
-    ! For reactions that may have multiple energy distributions such as (n.2n),
+    ! For reactions that may have multiple energy distributions such as (n,2n),
     ! this pointer allows multiple laws to be stored
     type(DistEnergy), pointer :: next => null()
-    
+
     ! Type-Bound procedures
     contains
       procedure :: clear => distenergy_clear ! Deallocates DistEnergy
@@ -57,7 +57,7 @@ module ace_header
     logical :: has_energy_dist         ! Energy distribution present?
     type(DistAngle)           :: adist ! Secondary angular distribution
     type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution
-    
+
     ! Type-Bound procedures
     contains
       procedure :: clear => reaction_clear ! Deallocates Reaction
@@ -76,7 +76,7 @@ module ace_header
     logical :: multiply_smooth ! multiply by smooth cross section?
     real(8), allocatable :: energy(:)   ! incident energies
     real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
-    
+
     ! Type-Bound procedures
     contains
       procedure :: clear => urrdata_clear ! Deallocates UrrData
@@ -137,22 +137,37 @@ module ace_header
     ! Reactions
     integer :: n_reaction ! # of reactions
     type(Reaction), pointer :: reactions(:) => null()
-    
+
     ! Type-Bound procedures
     contains
       procedure :: clear => nuclide_clear ! Deallocates Nuclide
   end type Nuclide
 
+!===============================================================================
+! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
+! thermal scattering collisions which utilize a continuous secondary energy
+! representation.
+!===============================================================================
+
+  type DistEnergySab
+    integer              :: n_e_out
+    real(8), allocatable :: e_out(:)
+    real(8), allocatable :: e_out_pdf(:)
+    real(8), allocatable :: e_out_cdf(:)
+    real(8), allocatable :: mu(:,:)
+  end type DistEnergySab
+
 !===============================================================================
 ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
 ! of light isotopes such as water, graphite, Be, etc
 !===============================================================================
-     
+
   type SAlphaBeta
-    character(10) :: name ! name of table, e.g. lwtr.10t
-    integer       :: zaid ! Z and A identifier, e.g. 6012 for Carbon-12
-    real(8)       :: awr  ! weight of nucleus in neutron masses
-    real(8)       :: kT   ! temperature in MeV (k*T)
+    character(10) :: name     ! name of table, e.g. lwtr.10t
+    real(8)       :: awr      ! weight of nucleus in neutron masses
+    real(8)       :: kT       ! temperature in MeV (k*T)
+    integer       :: n_zaid   ! Number of valid zaids
+    integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012
 
     ! threshold for S(a,b) treatment (usually ~4 eV)
     real(8) :: threshold_inelastic
@@ -162,11 +177,17 @@ module ace_header
     integer :: n_inelastic_e_in  ! # of incoming E for inelastic
     integer :: n_inelastic_e_out ! # of outgoing E for inelastic
     integer :: n_inelastic_mu    ! # of outgoing angles for inelastic
-    integer :: secondary_mode    ! secondary mode (equal/skewed)
+    integer :: secondary_mode    ! secondary mode (equal/skewed/continuous)
     real(8), allocatable :: inelastic_e_in(:)
-    real(8), allocatable :: inelastic_sigma(:) 
+    real(8), allocatable :: inelastic_sigma(:)
+    ! The following are used only if secondary_mode is 0 or 1
     real(8), allocatable :: inelastic_e_out(:,:)
     real(8), allocatable :: inelastic_mu(:,:,:)
+    ! The following is used only if secondary_mode is 3
+    ! The different implementation is necessary because the continuous
+    ! representation has a variable number of outgoing energy points for each
+    ! incoming energy
+    type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
 
     ! Elastic scattering data
     integer :: elastic_mode   ! elastic mode (discrete/exact)
@@ -214,8 +235,9 @@ module ace_header
     real(8) :: kappa_fission   ! microscopic energy-released from fission
 
     ! Information for S(a,b) use
-    integer :: index_sab   ! index in sab_tables (zero means no table)
-    real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
+    integer :: index_sab          ! index in sab_tables (zero means no table)
+    integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide
+    real(8) :: elastic_sab        ! microscopic elastic scattering on S(a,b) table
 
     ! Information for URR probability table use
     logical :: use_ptable  ! in URR range with probability tables?
@@ -239,125 +261,125 @@ module ace_header
 
 !===============================================================================
 ! DISTANGLE_CLEAR resets and deallocates data in Reaction.
-!===============================================================================    
-  
+!===============================================================================
+
     subroutine distangle_clear(this)
-      
+
       class(DistAngle), intent(inout) :: this ! The DistAngle object to clear
-      
+
       if (allocated(this % energy)) &
            deallocate(this % energy, this % type, this % location, this % data)
-      
-    end subroutine distangle_clear    
+
+    end subroutine distangle_clear
 
 !===============================================================================
 ! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
-!===============================================================================    
-  
+!===============================================================================
+
     recursive subroutine distenergy_clear(this)
-      
+
       class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
-      
+
       ! Clear p_valid
       call this % p_valid % clear()
-      
+
       if (allocated(this % data)) &
            deallocate(this % data)
-        
+
       if (associated(this % next)) then
         ! recursively clear this item
         call this % next % clear()
         deallocate(this % next)
       end if
-      
+
     end subroutine distenergy_clear
-    
+
 !===============================================================================
 ! REACTION_CLEAR resets and deallocates data in Reaction.
-!===============================================================================    
-  
+!===============================================================================
+
     subroutine reaction_clear(this)
-      
+
       class(Reaction), intent(inout) :: this ! The Reaction object to clear
-      
+
       if (allocated(this % sigma)) &
            deallocate(this % sigma)
-      
+
       if (associated(this % edist)) then
         call this % edist % clear()
         deallocate(this % edist)
       end if
-        
+
       call this % adist % clear()
-      
-    end subroutine reaction_clear    
-    
+
+    end subroutine reaction_clear
+
 !===============================================================================
 ! URRDATA_CLEAR resets and deallocates data in Reaction.
-!===============================================================================    
-  
+!===============================================================================
+
     subroutine urrdata_clear(this)
-      
+
       class(UrrData), intent(inout) :: this ! The UrrData object to clear
-      
+
       if (allocated(this % energy)) &
            deallocate(this % energy, this % prob)
-      
-    end subroutine urrdata_clear      
+
+    end subroutine urrdata_clear
 
 !===============================================================================
 ! NUCLIDE_CLEAR resets and deallocates data in Nuclide.
-!===============================================================================    
-  
+!===============================================================================
+
     subroutine nuclide_clear(this)
-      
+
       class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
-      
+
       integer :: i ! Loop counter
-      
+
       if (allocated(this % grid_index)) &
            deallocate(this % grid_index)
-      
+
       if (allocated(this % energy)) &
            deallocate(this % total, this % elastic, this % fission,  &
           this % nu_fission, this % absorption)
       if (allocated(this % heating)) &
            deallocate(this % heating)
-      
+
       if (allocated(this % index_fission)) &
            deallocate(this % index_fission)
-        
+
       if (allocated(this % nu_t_data)) &
            deallocate(this % nu_t_data)
-        
+
       if (allocated(this % nu_p_data)) &
            deallocate(this % nu_p_data)
-        
+
       if (allocated(this % nu_d_data)) &
            deallocate(this % nu_d_data)
-        
+
       if (allocated(this % nu_d_precursor_data)) &
            deallocate(this % nu_d_precursor_data)
-        
+
       if (associated(this % nu_d_edist)) then
         do i = 1, size(this % nu_d_edist)
           call this % nu_d_edist(i) % clear()
         end do
         deallocate(this % nu_d_edist)
       end if
-      
+
       if (associated(this % urr_data)) then
         call this % urr_data % clear()
         deallocate(this % urr_data)
       end if
-      
+
       if (associated(this % reactions)) then
         do i = 1, size(this % reactions)
           call this % reactions(i) % clear()
         end do
         deallocate(this % reactions)
       end if
-      
-    end subroutine nuclide_clear    
+
+    end subroutine nuclide_clear
 
 end module ace_header

src/cmfd_data.F90

@@ -767,15 +767,15 @@ contains
     use global,     only: cmfd, cmfd_hold_weights
 
     real(8) :: get_reflector_albedo ! reflector albedo
-    integer :: i                    ! iteration counter for x
-    integer :: j                    ! iteration counter for y
-    integer :: k                    ! iteration counter for z
-    integer :: g                    ! iteration counter for groups
-    integer :: l                    ! iteration counter for leakages
+    integer, intent(in) :: i ! iteration counter for x
+    integer, intent(in) :: j ! iteration counter for y
+    integer, intent(in) :: k ! iteration counter for z
+    integer, intent(in) :: g ! iteration counter for groups
+    integer, intent(in) :: l ! iteration counter for leakages
 
-    integer :: shift_idx            ! parameter to shift index by +1 or -1
-    real(8) :: current(12)          ! partial currents for all faces of mesh cell            
-    real(8) :: albedo               ! the albedo
+    integer :: shift_idx   ! parameter to shift index by +1 or -1
+    real(8) :: current(12) ! partial currents for all faces of mesh cell            
+    real(8) :: albedo      ! the albedo
 
     ! Get partial currents from object
     current = cmfd%current(:,g,i,j,k)
@@ -798,7 +798,7 @@ contains
   end function get_reflector_albedo
 
 !===============================================================================
-! FIX_NEUTRON_BALANCE
+! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance
 !===============================================================================
 
   subroutine fix_neutron_balance()
@@ -934,7 +934,7 @@ contains
   subroutine compute_effective_downscatter()
 
     use constants, only: ZERO, CMFD_NOACCEL
-    use global,    only: cmfd, cmfd_downscatter
+    use global,    only: cmfd
 
     integer :: nx                ! number of mesh cells in x direction
     integer :: ny                ! number of mesh cells in y direction

src/cmfd_execute.F90

@@ -122,9 +122,9 @@ contains
 
   subroutine process_cmfd_options()
 
+#ifdef PETSC
     use global,       only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err
 
-#ifdef PETSC
     ! Check for snes monitor
     if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", &
          "stdout", mpi_err)
@@ -143,9 +143,10 @@ contains
   subroutine calc_fission_source()
 
     use constants,  only: CMFD_NOACCEL, ZERO, TWO
-    use global,     only: cmfd, cmfd_coremap, master, mpi_err, entropy_on, &
-                          current_batch
+    use global,     only: cmfd, cmfd_coremap, master, entropy_on, current_batch
+
 #ifdef MPI
+    use global,     only: mpi_err
     use mpi
 #endif
 
@@ -259,34 +260,36 @@ contains
   end subroutine calc_fission_source
 
 !===============================================================================
-! CMFD_REWEIGHT
+! CMFD_REWEIGHT performs weighting of particles in the source bank
 !===============================================================================
 
   subroutine cmfd_reweight(new_weights)
 
     use constants,   only: ZERO, ONE
     use error,       only: warning, fatal_error
-    use global,      only: n_particles, meshes, source_bank, work,             &
-                           n_user_meshes, message, cmfd, master, mpi_err
+    use global,      only: meshes, source_bank, work, n_user_meshes, message, &
+                           cmfd, master
     use mesh_header, only: StructuredMesh
     use mesh,        only: count_bank_sites, get_mesh_indices
     use search,      only: binary_search
 
 #ifdef MPI
+    use global,      only: mpi_err
     use mpi
 #endif
 
-    integer :: nx                      ! maximum number of cells in x direction
-    integer :: ny                      ! maximum number of cells in y direction
-    integer :: nz                      ! maximum number of cells in z direction
-    integer :: ng                      ! maximum number of energy groups
-    integer :: i                       ! iteration counter
-    integer :: ijk(3)                  ! spatial bin location
-    integer :: e_bin                   ! energy bin of source particle
-    integer :: n_groups                ! number of energy groups
-    logical :: outside                 ! any source sites outside mesh
-    logical :: in_mesh                 ! source site is inside mesh
-    logical :: new_weights             ! calcualte new weights
+    logical, intent(in) :: new_weights ! calcualte new weights
+
+    integer :: nx       ! maximum number of cells in x direction
+    integer :: ny       ! maximum number of cells in y direction
+    integer :: nz       ! maximum number of cells in z direction
+    integer :: ng       ! maximum number of energy groups
+    integer :: i        ! iteration counter
+    integer :: ijk(3)   ! spatial bin location
+    integer :: e_bin    ! energy bin of source particle
+    integer :: n_groups ! number of energy groups
+    logical :: outside  ! any source sites outside mesh
+    logical :: in_mesh  ! source site is inside mesh
 
     type(StructuredMesh), pointer :: m ! point to mesh
     real(8), allocatable :: egrid(:)   ! energy grid
@@ -394,13 +397,13 @@ contains
     use global, only: cmfd, cmfd_coremap
 
     integer :: matidx ! the index location in matrix
-    integer :: i      ! current x index
-    integer :: j      ! current y index
-    integer :: k      ! current z index
-    integer :: g      ! current group index
-    integer :: nx     ! maximum number of cells in x direction
-    integer :: ny     ! maximum number of cells in y direction
-    integer :: ng     ! maximum number of energy groups
+    integer, intent(in) :: i  ! current x index
+    integer, intent(in) :: j  ! current y index
+    integer, intent(in) :: k  ! current z index
+    integer, intent(in) :: g  ! current group index
+    integer, intent(in) :: nx ! maximum number of cells in x direction
+    integer, intent(in) :: ny ! maximum number of cells in y direction
+    integer, intent(in) :: ng ! maximum number of energy groups
 
     ! Check if coremap is used
     if (cmfd_coremap) then
@@ -418,7 +421,7 @@ contains
   end function get_matrix_idx
 
 !===============================================================================
-! CMFD_TALLY_RESET
+! CMFD_TALLY_RESET resets all cmfd tallies
 !===============================================================================
 
   subroutine cmfd_tally_reset()

src/cmfd_header.F90

@@ -93,8 +93,8 @@ contains
 
   subroutine allocate_cmfd(this, n_batches)
 
-    integer :: n_batches
-    type(cmfd_type) :: this
+    integer, intent(in)            :: n_batches ! number of batches in calc
+    type(cmfd_type), intent(inout) :: this      ! cmfd instance
 
     integer :: nx  ! number of mesh cells in x direction
     integer :: ny  ! number of mesh cells in y direction
@@ -169,7 +169,7 @@ contains
 
   subroutine deallocate_cmfd(this)
 
-    type(cmfd_type) :: this
+    type(cmfd_type), intent(inout) :: this ! cmfd instance
 
     if (allocated(this % egrid))         deallocate(this % egrid)
     if (allocated(this % totalxs))       deallocate(this % totalxs)

src/cmfd_input.F90

@@ -20,7 +20,9 @@ contains
 
     use cmfd_header,  only: allocate_cmfd
 
+#ifdef PETSC
     integer :: new_comm ! new mpi communicator
+#endif
     integer :: color    ! color group of processor
 
     ! Read in cmfd input file
@@ -34,17 +36,17 @@ contains
     end if
 
     ! Split up procs
-# ifdef PETSC 
-    call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,0,new_comm,mpi_err)
-# endif
+#ifdef PETSC
+    call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, new_comm, mpi_err)
+#endif
 
     ! assign to PETSc
-# ifdef PETSC
+#ifdef PETSC
     PETSC_COMM_WORLD = new_comm
 
     ! Initialize PETSc on all procs
-    call PetscInitialize(PETSC_NULL_CHARACTER,mpi_err)
-# endif
+    call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err)
+#endif
 
     ! Initialize timers
     call time_cmfd % reset()
@@ -61,16 +63,22 @@ contains
 !===============================================================================
 
   subroutine read_cmfd_xml()
-    
+
     use error,   only: fatal_error, warning
+    use global
     use output,  only: write_message
     use string,  only: lower_case
-    use xml_data_cmfd_t
+    use xml_interface
     use, intrinsic :: ISO_FORTRAN_ENV
 
-    integer :: ng                       ! number of energy groups
-    logical :: file_exists              ! does cmfd.xml exist?
-    character(MAX_LINE_LEN) :: filename ! name of input file
+    integer :: ng
+    integer, allocatable :: iarray(:)
+    logical :: file_exists ! does cmfd.xml exist?
+    logical :: found
+    character(MAX_LINE_LEN) :: filename
+    character(MAX_LINE_LEN) :: temp_str
+    type(Node), pointer :: doc => null()
+    type(Node), pointer :: node_mesh => null()
 
     ! Read cmfd input file
     filename = trim(path_input) // "cmfd.xml"
@@ -91,16 +99,25 @@ contains
     end if
 
     ! Parse cmfd.xml file
-    call read_xml_file_cmfd_t(filename)
+    call open_xmldoc(doc, filename)
 
-    ! Set spatial dimensions in cmfd object (structed Cartesian mesh)
-    cmfd % indices(1:3) = mesh_ % dimension(1:3)
+    ! Get pointer to mesh XML node
+    call get_node_ptr(doc, "mesh", node_mesh, found = found)
 
-    ! Get number of energy groups or set to 1 group default
-    if (associated(mesh_ % energy)) then
-      ng = size(mesh_ % energy)
-      if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(ng))
-      cmfd % egrid = mesh_ % energy 
+    ! Check if mesh is there
+    if (.not.found) then
+      message = "No CMFD mesh specified in CMFD XML file."
+      call fatal_error()
+    end if
+
+    ! Set spatial dimensions in cmfd object
+    call get_node_array(node_mesh, "dimension", cmfd % indices(1:3))
+
+    ! Get number of energy groups
+    if (check_for_node(node_mesh, "energy")) then
+      ng = get_arraysize_double(node_mesh, "energy")
+      if(.not.allocated(cmfd%egrid)) allocate(cmfd%egrid(ng))
+      call get_node_array(node_mesh, "energy", cmfd%egrid)
       cmfd % indices(4) = ng - 1 ! sets energy group dimension
     else
       if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
@@ -108,90 +125,115 @@ contains
       cmfd % indices(4) = 1 ! one energy group
     end if
     
-    ! Set global albedo, these can be overwritten by coremap
-    if (associated(mesh_ % albedo)) then
-      cmfd % albedo = mesh_ % albedo
+    ! Set global albedo
+    if (check_for_node(node_mesh, "albedo")) then
+      call get_node_array(node_mesh, "albedo", cmfd % albedo)
     else
       cmfd % albedo = (/1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
     end if
 
-    ! Get core map overlay for a subset mesh for CMFD 
-    if (associated(mesh_ % map)) then
-
-      ! Allocate a core map with appropriate dimensions
+    ! Get acceleration map
+    if (check_for_node(node_mesh, "map")) then
       allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
            cmfd % indices(3)))
-
-      ! Check and make sure it is of correct size
-      if (size(mesh_ % map) /= product(cmfd % indices(1:3))) then
-        message = 'CMFD coremap not to correct dimensions'
+      if (get_arraysize_integer(node_mesh, "map") /= &
+          product(cmfd % indices(1:3))) then
+        message = 'FATAL==>CMFD coremap not to correct dimensions'
         call fatal_error() 
       end if
-
-      ! Reshape core map vector into (x,y,z) array
-      cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3)))
-
-      ! Indicate to cmfd that a core map overlay is active
+      allocate(iarray(get_arraysize_integer(node_mesh, "map")))
+      call get_node_array(node_mesh, "map", iarray)
+      cmfd % coremap = reshape(iarray,(cmfd % indices(1:3)))
       cmfd_coremap = .true.
-
-      ! Write to stdout that a core map overlay is active
-      message = "Core Map Overlay Activated"
-      call write_message(10)
-
+      deallocate(iarray)
     end if
 
-    ! Get normalization constant for source (default is 1.0 from XML) 
-    cmfd % norm = norm_
+    ! Check for normalization constant
+    if (check_for_node(doc, "norm")) then
+      call get_node_value(doc, "norm", cmfd % norm)
+    end if
 
     ! Set feedback logical
-    call lower_case(feedback_)
-    if (feedback_ == 'true' .or. feedback_ == '1') cmfd_feedback = .true.
+    if (check_for_node(doc, "feedback")) then
+      call get_node_value(doc, "feedback", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_feedback = .true.
+    end if
 
     ! Set downscatter logical
-    call lower_case(downscatter_)
-    if (downscatter_ == 'true' .or. downscatter_ == '1') &
-         cmfd_downscatter = .true.
+    if (check_for_node(doc, "downscatter")) then
+      call get_node_value(doc, "downscatter", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_downscatter = .true.
+    end if
 
-    ! Set the solver type (default power from XML)
-    cmfd_solver_type = solver_(1:10)
+    ! Set the solver type
+    if (check_for_node(doc, "solver")) &
+      call get_node_value(doc, "solver", cmfd_solver_type)
 
-    ! Set convergence monitoring 
-    call lower_case(snes_monitor_)
-    call lower_case(ksp_monitor_)
-    call lower_case(power_monitor_)
-    if (snes_monitor_ == 'true' .or. snes_monitor_ == '1') &
-         cmfd_snes_monitor = .true.
-    if (ksp_monitor_ == 'true' .or. ksp_monitor_ == '1') &
-         cmfd_ksp_monitor = .true.
-    if (power_monitor_ == 'true' .or. power_monitor_ == '1') &
-         cmfd_power_monitor = .true.
+    ! Set monitoring
+    if (check_for_node(doc, "snes_monitor")) then
+      call get_node_value(doc, "snes_monitor", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_snes_monitor = .true.
+    end if
+    if (check_for_node(doc, "ksp_monitor")) then
+      call get_node_value(doc, "ksp_monitor", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_ksp_monitor = .true.
+    end if
+    if (check_for_node(doc, "power_monitor")) then
+      call get_node_value(doc, "power_monitor", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_power_monitor = .true.
+    end if
 
     ! Output logicals
-    call lower_case(write_matrices_)
-    if (write_matrices_ == 'true' .or. write_matrices_ == '1') &
-         cmfd_write_matrices = .true.
+    if (check_for_node(doc, "write_matrices")) then
+      call get_node_value(doc, "write_matices", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_write_matrices = .true.
+    end if
 
     ! Run an adjoint calc
-    call lower_case(run_adjoint_)
-    if (run_adjoint_ == 'true' .or. run_adjoint_ == '1') &
-         cmfd_run_adjoint = .true.
+    if (check_for_node(doc, "run_adjoint")) then
+      call get_node_value(doc, "run_adjoint", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+        cmfd_run_adjoint = .true.
+    end if
 
-    ! Batch to begin cmfd (default is 1 from XML)
-    cmfd_begin = begin_
+    ! Batch to begin cmfd
+    if (check_for_node(doc, "begin")) &
+      call get_node_value(doc, "begin", cmfd_begin) 
 
-    ! Tally during inactive batches (by default we will always tally from 1)
-    call lower_case(inactive_)
-    if (inactive_ == 'false' .or. inactive_ == '0') cmfd_tally_on = .false.
+    ! Tally during inactive batches
+    if (check_for_node(doc, "inactive")) then
+      call get_node_value(doc, "inactive", temp_str)
+      call lower_case(temp_str)
+      if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') &
+        cmfd_tally_on = .false.
+    end if
 
     ! Inactive batch flush window
-    cmfd_inact_flush(1) = inactive_flush_ ! the interval of batches
-    cmfd_inact_flush(2) = num_flushes_    ! number of times to do this
+    if (check_for_node(doc, "inactive_flush")) &
+      call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1))
+    if (check_for_node(doc, "num_flushes")) &
+      call get_node_value(doc, "num_flushes", cmfd_inact_flush(2))
 
     ! Last flush before active batches
-    cmfd_act_flush = active_flush_
+    if (check_for_node(doc, "active_flush")) &
+      call get_node_value(doc, "active_flush", cmfd_act_flush)
 
     ! Get display
-    cmfd_display = display_
+    if (check_for_node(doc, "display")) &
+      call get_node_value(doc, "display", cmfd_display)
     if (trim(cmfd_display) == 'dominance' .and. &
         trim(cmfd_solver_type) /= 'power') then
       message = 'Dominance Ratio only aviable with power iteration solver'
@@ -200,15 +242,24 @@ contains
     end if
 
     ! Read in spectral radius estimate and tolerances
-    cmfd_spectral = spectral_
-    if (shift_ /= ZERO) cmfd_shift = shift_
-    cmfd_ktol = ktol_
-    cmfd_stol = stol_
-    cmfd_atoli = atoli_
-    cmfd_rtoli = rtoli_
-
+    if (check_for_node(doc, "spectral")) &
+      call get_node_value(doc, "spectral", cmfd_spectral)
+    if (check_for_node(doc, "shift")) &
+      call get_node_value(doc, "shift", cmfd_shift)
+    if (check_for_node(doc, "ktol")) &
+      call get_node_value(doc, "ktol", cmfd_ktol)
+    if (check_for_node(doc, "stol")) &
+      call get_node_value(doc, "stol", cmfd_stol)
+    if (check_for_node(doc, "atoli")) &
+      call get_node_value(doc, "atoli", cmfd_atoli)
+    if (check_for_node(doc, "rtoli")) &
+      call get_node_value(doc, "rtoli", cmfd_rtoli)
+   
     ! Create tally objects
-    call create_cmfd_tally()
+    call create_cmfd_tally(doc)
+
+    ! Close CMFD XML file
+    call close_xmldoc(doc)
 
   end subroutine read_cmfd_xml
 
@@ -221,29 +272,31 @@ contains
 !   3: Surface current
 !===============================================================================
 
-  subroutine create_cmfd_tally()
+  subroutine create_cmfd_tally(doc)
 
+    use constants,        only: MAX_LINE_LEN
     use error,            only: fatal_error, warning
     use mesh_header,      only: StructuredMesh
     use string
     use tally,            only: setup_active_cmfdtallies
     use tally_header,     only: TallyObject, TallyFilter
     use tally_initialize, only: add_tallies
-    use xml_data_cmfd_t
+    use xml_interface
 
+    type(Node), pointer :: doc ! pointer to XML doc info
+
+    character(MAX_LINE_LEN) :: temp_str ! temp string
     integer :: i           ! loop counter
     integer :: n           ! size of arrays in mesh specification
     integer :: ng          ! number of energy groups (default 1)
     integer :: n_filters   ! number of filters
     integer :: i_filter_mesh ! index for mesh filter
-    character(MAX_LINE_LEN) :: filename
+    integer :: iarray3(3) ! temp integer array
+    real(8) :: rarray3(3) ! temp double array
     type(TallyObject),    pointer :: t => null()
     type(StructuredMesh), pointer :: m => null()
     type(TallyFilter) :: filters(N_FILTER_TYPES) ! temporary filters
-
-    ! Parse cmfd.xml file
-     filename = trim(path_input) // "cmfd.xml"
-     call read_xml_file_cmfd_t(filename)
+    type(Node), pointer :: node_mesh => null()
 
     ! Set global variables if they are 0 (this can happen if there is no tally
     ! file)
@@ -259,8 +312,11 @@ contains
     ! Set mesh type to rectangular
     m % type = LATTICE_RECT
 
+    ! Get pointer to mesh XML node
+    call get_node_ptr(doc, "mesh", node_mesh)
+
     ! Determine number of dimensions for mesh
-    n = size(mesh_ % dimension)
+    n = get_arraysize_integer(node_mesh, "dimension")
     if (n /= 2 .and. n /= 3) then
        message = "Mesh must be two or three dimensions."
        call fatal_error()
@@ -274,75 +330,80 @@ contains
     allocate(m % upper_right(n))
 
     ! Check that dimensions are all greater than zero
-    if (any(mesh_ % dimension <= 0)) then
-       message = "All entries on the <dimension> element for a tally mesh &
-            &must be positive."
-       call fatal_error()
+    call get_node_array(node_mesh, "dimension", iarray3(1:n))
+    if (any(iarray3(1:n) <= 0)) then
+      message = "All entries on the <dimension> element for a tally mesh &
+           &must be positive."
+      call fatal_error()
     end if
 
     ! Read dimensions in each direction
-    m % dimension = mesh_ % dimension
+    m % dimension = iarray3(1:n)
 
     ! Read mesh lower-left corner location
-    if (m % n_dimension /= size(mesh_ % lower_left)) then
-       message = "Number of entries on <lower_left> must be the same as &
-            &the number of entries on <dimension>."
-       call fatal_error()
+    if (m % n_dimension /= get_arraysize_double(node_mesh, "lower_left")) then
+      message = "Number of entries on <lower_left> must be the same as &
+           &the number of entries on <dimension>."
+      call fatal_error()
     end if
-    m % lower_left = mesh_ % lower_left
+    call get_node_array(node_mesh, "lower_left", m % lower_left)
 
-    ! Make sure either upper-right or width was specified
-    if (associated(mesh_ % upper_right) .and. &
-         associated(mesh_ % width)) then
-       message = "Cannot specify both <upper_right> and <width> on a &
-             &tally mesh."
-       call fatal_error()
+    ! Make sure both upper-right or width were specified
+    if (check_for_node(node_mesh, "upper_right") .and. &
+        check_for_node(node_mesh, "width")) then
+      message = "Cannot specify both <upper_right> and <width> on a &
+           &tally mesh."
+      call fatal_error()
     end if
 
     ! Make sure either upper-right or width was specified
-    if (.not. associated(mesh_ % upper_right) .and. &
-         .not. associated(mesh_ % width)) then
-       message = "Must specify either <upper_right> and <width> on a &
-            &tally mesh."
-       call fatal_error()
+    if (.not.check_for_node(node_mesh, "upper_right") .and. &
+        .not.check_for_node(node_mesh, "width")) then
+      message = "Must specify either <upper_right> and <width> on a &
+           &tally mesh."
+      call fatal_error()
     end if
 
-    if (associated(mesh_ % width)) then
-       ! Check to ensure width has same dimensions
-       if (size(mesh_ % width) /= size(mesh_ % lower_left)) then
-          message = "Number of entries on <width> must be the same as the &
-               &number of entries on <lower_left>."
-          call fatal_error()
-       end if
+    if (check_for_node(node_mesh, "width")) then
+      ! Check to ensure width has same dimensions
+      if (get_arraysize_double(node_mesh, "width") /= &
+          get_arraysize_double(node_mesh, "lower_left")) then
+        message = "Number of entries on <width> must be the same as the &
+             &number of entries on <lower_left>."
+        call fatal_error()
+      end if
 
-       ! Check for negative widths
-       if (any(mesh_ % width < ZERO)) then
-          message = "Cannot have a negative <width> on a tally mesh."
-          call fatal_error()
-       end if
+      ! Check for negative widths
+      call get_node_array(node_mesh, "width", rarray3(1:n))
+      if (any(rarray3(1:n) < ZERO)) then
+        message = "Cannot have a negative <width> on a tally mesh."
+        call fatal_error()
+      end if
 
-       ! Set width and upper right coordinate
-       m % width = mesh_ % width
-       m % upper_right = m % lower_left + m % dimension * m % width
+      ! Set width and upper right coordinate
+      m % width = rarray3(1:n)
+      m % upper_right = m % lower_left + m % dimension * m % width
 
-    elseif (associated(mesh_ % upper_right)) then
-       ! Check to ensure width has same dimensions
-       if (size(mesh_ % upper_right) /= size(mesh_ % lower_left)) then
-          message = "Number of entries on <upper_right> must be the same as &
-               &the number of entries on <lower_left>."
-          call fatal_error()
-       end if
+    elseif (check_for_node(node_mesh, "upper_right")) then
+      ! Check to ensure width has same dimensions
+      if (get_arraysize_double(node_mesh, "upper_right") /= &
+          get_arraysize_double(node_mesh, "lower_left")) then
+        message = "Number of entries on <upper_right> must be the same as &
+             &the number of entries on <lower_left>."
+        call fatal_error()
+      end if
 
-       ! Check that upper-right is above lower-left
-       if (any(mesh_ % upper_right < mesh_ % lower_left)) then
-          message = "The <upper_right> coordinates must be greater than the &
-               &<lower_left> coordinates on a tally mesh."
-          call fatal_error()
-       end if
+      ! Check that upper-right is above lower-left
+      call get_node_array(node_mesh, "upper_right", rarray3(1:n))
+      if (any(rarray3(1:n) < m % lower_left)) then
+        message = "The <upper_right> coordinates must be greater than the &
+             &<lower_left> coordinates on a tally mesh."
+        call fatal_error()
+      end if
 
-       ! Set width and upper right coordinate
-       m % upper_right = mesh_ % upper_right
-       m % width = (m % upper_right - m % lower_left) / m % dimension
+      ! Set upper right coordinate and width
+      m % upper_right = rarray3(1:n)
+      m % width = (m % upper_right - m % lower_left) / real(m % dimension, 8)
     end if
 
     ! Set volume fraction
@@ -361,8 +422,12 @@ contains
       t => cmfd_tallies(i)
 
       ! Set reset property
-      call lower_case(reset_)
-      if (reset_ == 'true' .or. reset_ == '1') t % reset = .true.
+      if (check_for_node(doc, "reset")) then
+        call get_node_value(doc, "reset", temp_str)
+        call lower_case(temp_str)
+        if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
+          t % reset = .true.
+      end if
 
       ! Set up mesh filter
       n_filters = 1
@@ -373,13 +438,14 @@ contains
       t % find_filter(FILTER_MESH) = n_filters
 
       ! Read and set incoming energy mesh filter
-      if (associated(mesh_ % energy)) then
+      if (check_for_node(node_mesh, "energy")) then
         n_filters = n_filters + 1
         filters(n_filters) % type = FILTER_ENERGYIN
-        ng = size(mesh_ % energy)
+        ng = get_arraysize_double(node_mesh, "energy") 
         filters(n_filters) % n_bins = ng - 1
         allocate(filters(n_filters) % real_bins(ng))
-        filters(n_filters) % real_bins = mesh_ % energy
+        call get_node_array(node_mesh, "energy", &
+             filters(n_filters) % real_bins)
         t % find_filter(FILTER_ENERGYIN) = n_filters
       end if
 
@@ -433,14 +499,15 @@ contains
         ! Set tally type to volume
         t % type = TALLY_VOLUME
 
-        ! Read and set outgoing energy mesh filter
-        if (associated(mesh_ % energy)) then
+        ! read and set outgoing energy mesh filter
+        if (check_for_node(node_mesh, "energy")) then
           n_filters = n_filters + 1
           filters(n_filters) % type = FILTER_ENERGYOUT
-          ng = size(mesh_ % energy)
+          ng = get_arraysize_double(node_mesh, "energy")
           filters(n_filters) % n_bins = ng - 1
           allocate(filters(n_filters) % real_bins(ng))
-          filters(n_filters) % real_bins = mesh_ % energy
+          call get_node_array(node_mesh, "energy", &
+               filters(n_filters) % real_bins)
           t % find_filter(FILTER_ENERGYOUT) = n_filters
         end if
 
@@ -449,8 +516,8 @@ contains
         allocate(t % filters(n_filters))
         t % filters = filters(1:n_filters)
 
-        ! Deallocate filters bins array
-        if (associated(mesh_ % energy)) &
+        ! deallocate filters bins array
+        if (check_for_node(node_mesh, "energy")) &
              deallocate(filters(n_filters) % real_bins)
 
         ! Allocate macro reactions
@@ -518,12 +585,15 @@ contains
 
       ! Deallocate filter bins
       deallocate(filters(1) % int_bins)
-      if (associated(mesh_ % energy)) deallocate(filters(2) % real_bins)
+      if (check_for_node(node_mesh, "energy")) &
+        deallocate(filters(2) % real_bins)
 
     end do
 
     ! Put cmfd tallies into active tally array and turn tallies on
+!$omp parallel
     call setup_active_cmfdtallies()
+!$omp end parallel
     tallies_on = .true.
 
   end subroutine create_cmfd_tally

src/cmfd_jfnk_solver.F90

@@ -30,7 +30,7 @@ contains
 
     use global,  only: time_cmfdbuild, time_cmfdsolve
 
-    logical, optional :: adjoint ! adjoint calculation
+    logical, intent(in), optional :: adjoint ! adjoint calculation
 
     ! Check for adjoint
     adjoint_calc = .false.
@@ -84,8 +84,8 @@ contains
     use constants, only: ZERO, ONE
     use global,    only: cmfd, cmfd_adjoint_type, current_batch
 
-    logical :: physical_adjoint
-    integer :: n
+    logical :: physical_adjoint ! physical adjoing calculation logical
+    integer :: n ! size of matrices
 
     ! Set up all matrices
     call init_loss_matrix(loss)
@@ -153,8 +153,8 @@ contains
 
   subroutine init_jacobian_matrix()
 
-    integer :: nnz
-    integer :: n
+    integer :: nnz ! number of nonzeros in matrix
+    integer :: n   ! dimension of matrix
 
     ! Get length of matrix and number of nonzeros total in loss matrix
     nnz = loss % nnz
@@ -179,7 +179,7 @@ contains
 
     use constants,  only: ONE
 
-    type(Vector) :: x
+    type(Vector), intent(in) :: x ! solution vector
 
     integer      :: i      ! loop counter for jacobian rows
     integer      :: jjac   ! loop counter for jacobian cols
@@ -267,8 +267,8 @@ contains
 
     use global,       only: cmfd_write_matrices
 
-    type(Vector) :: x
-    type(Vector) :: res
+    type(Vector), intent(in)    :: x   ! solution vector
+    type(Vector), intent(inout) :: res ! residual vector
 
     character(len=25) :: filename
     integer           :: n

src/cmfd_loss_operator.F90

@@ -16,7 +16,7 @@ contains
 
   subroutine init_loss_matrix(loss_matrix)
 
-    type(Matrix) :: loss_matrix ! cmfd loss matrix
+    type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
 
     integer :: nx   ! maximum number of x cells
     integer :: ny   ! maximum number of y cells
@@ -26,8 +26,10 @@ contains
     integer :: nnz  ! number of nonzeros in matrix
     integer :: n_i  ! number of interior cells
     integer :: n_c  ! number of corner cells
-    integer :: n_s  ! number side cells
+    integer :: n_s  ! number of side cells
+    integer :: n_e  ! number of edge cells
     integer :: nz_c ! number of non-zero corner cells
+    integer :: nz_e ! number of non-zero edge cells
     integer :: nz_s ! number of non-zero side cells
     integer :: nz_i ! number of non-zero interior cells
 
@@ -48,13 +50,16 @@ contains
     if (cmfd_coremap) then
       nnz = preallocate_loss_matrix(nx, ny, nz, ng, n)
     else ! structured Cartesian grid
-      n_c  = 4                   ! define # of corners
-      n_s  = 2*(nx + ny) - 8     ! define # of sides
-      n_i  = nx*ny - (n_c + n_S) ! define # of interiors  
-      nz_c = ng*n_c*(3 + ng - 1) ! define # nonzero corners
-      nz_s = ng*n_s*(4 + ng - 1) ! define # nonzero sides
-      nz_i = ng*n_i*(5 + ng - 1) ! define # nonzero interiors
-      nnz  = nz_c + nz_s + nz_i
+      n_c  = 8                   ! define # of corners
+      n_e  = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges
+      n_s  = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) &
+           + 2*(ny - 2)*(nz - 2) ! define # of sides
+      n_i  = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors  
+      nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners
+      nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges
+      nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides
+      nz_i = ng*n_i*(7 + ng - 1) ! define # nonzero interiors
+      nnz  = nz_c + nz_e + nz_s + nz_i
     end if
 
     ! Configure loss matrix
@@ -68,12 +73,12 @@ contains
 
   function preallocate_loss_matrix(nx, ny, nz, ng, n) result(nnz)
 
-    integer, intent(in)    :: nx   ! maximum number of x cells
-    integer, intent(in)    :: ny   ! maximum number of y cells
-    integer, intent(in)    :: nz   ! maximum number of z cells
-    integer, intent(in)    :: ng   ! maximum number of groups
-    integer, intent(in)    :: n    ! total length of matrix
-    integer                :: nnz  ! number of nonzeros
+    integer, intent(in) :: nx  ! maximum number of x cells
+    integer, intent(in) :: ny  ! maximum number of y cells
+    integer, intent(in) :: nz  ! maximum number of z cells
+    integer, intent(in) :: ng  ! maximum number of groups
+    integer, intent(in) :: n   ! total length of matrix
+    integer             :: nnz ! number of nonzeros
 
     integer :: i             ! iteration counter for x
     integer :: j             ! iteration counter for y
@@ -181,39 +186,39 @@ contains
 
   subroutine build_loss_matrix(loss_matrix, adjoint)
 
-    type(Matrix) :: loss_matrix     ! cmfd loss matrix
-    logical, optional :: adjoint    ! set up the adjoint
+    type(Matrix), intent(inout)   :: loss_matrix ! cmfd loss matrix
+    logical, intent(in), optional :: adjoint     ! set up the adjoint
 
-    integer :: nxyz(3,2)            ! single vector containing bound. locations
-    integer :: i                    ! iteration counter for x
-    integer :: j                    ! iteration counter for y
-    integer :: k                    ! iteration counter for z
-    integer :: g                    ! iteration counter for groups
-    integer :: l                    ! iteration counter for leakages
-    integer :: h                    ! energy group when doing scattering
-    integer :: nx                   ! maximum number of x cells
-    integer :: ny                   ! maximum number of y cells
-    integer :: nz                   ! maximum number of z cells
-    integer :: ng                   ! maximum number of groups
-    integer :: neig_mat_idx         ! matrix index of neighbor cell
-    integer :: scatt_mat_idx        ! matrix index for h-->g scattering terms
-    integer :: bound(6)             ! vector for comparing when looking for bound
-    integer :: xyz_idx              ! index for determining if x,y or z leakage
-    integer :: dir_idx              ! index for determining - or + face of cell
-    integer :: neig_idx(3)          ! spatial indices of neighbour
-    integer :: shift_idx            ! parameter to shift index by +1 or -1
-    integer :: irow                 ! iteration counter over row
-    logical :: adjoint_calc         ! is this a physical adjoint calculation?
-    real(8) :: totxs                ! total macro cross section
-    real(8) :: scattxsgg            ! scattering macro cross section g-->g
-    real(8) :: scattxshg            ! scattering macro cross section h-->g
-    real(8) :: dtilde(6)            ! finite difference coupling parameter
-    real(8) :: dhat(6)              ! nonlinear coupling parameter
-    real(8) :: hxyz(3)              ! cell lengths in each direction
-    real(8) :: jn                   ! direction dependent leakage coeff to neig
-    real(8) :: jo(6)                ! leakage coeff in front of cell flux
-    real(8) :: jnet                 ! net leakage from jo
-    real(8) :: val                  ! temporary variable before saving to 
+    integer :: nxyz(3,2)     ! single vector containing bound. locations
+    integer :: i             ! iteration counter for x
+    integer :: j             ! iteration counter for y
+    integer :: k             ! iteration counter for z
+    integer :: g             ! iteration counter for groups
+    integer :: l             ! iteration counter for leakages
+    integer :: h             ! energy group when doing scattering
+    integer :: nx            ! maximum number of x cells
+    integer :: ny            ! maximum number of y cells
+    integer :: nz            ! maximum number of z cells
+    integer :: ng            ! maximum number of groups
+    integer :: neig_mat_idx  ! matrix index of neighbor cell
+    integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
+    integer :: bound(6)      ! vector for comparing when looking for bound
+    integer :: xyz_idx       ! index for determining if x,y or z leakage
+    integer :: dir_idx       ! index for determining - or + face of cell
+    integer :: neig_idx(3)   ! spatial indices of neighbour
+    integer :: shift_idx     ! parameter to shift index by +1 or -1
+    integer :: irow          ! iteration counter over row
+    logical :: adjoint_calc  ! is this a physical adjoint calculation?
+    real(8) :: totxs         ! total macro cross section
+    real(8) :: scattxsgg     ! scattering macro cross section g-->g
+    real(8) :: scattxshg     ! scattering macro cross section h-->g
+    real(8) :: dtilde(6)     ! finite difference coupling parameter
+    real(8) :: dhat(6)       ! nonlinear coupling parameter
+    real(8) :: hxyz(3)       ! cell lengths in each direction
+    real(8) :: jn            ! direction dependent leakage coeff to neig
+    real(8) :: jo(6)         ! leakage coeff in front of cell flux
+    real(8) :: jnet          ! net leakage from jo
+    real(8) :: val           ! temporary variable before saving to 
 
     ! Check for adjoint
     adjoint_calc = .false.
@@ -366,14 +371,14 @@ contains
 
   subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
 
-    integer :: matidx ! the index location in matrix
-    integer :: i      ! current x index
-    integer :: j      ! current y index
-    integer :: k      ! current z index
-    integer :: g      ! current group index
-    integer :: nx     ! maximum number of x cells
-    integer :: ny     ! maximum number of y cells
-    integer :: ng     ! maximum number of groups
+    integer, intent(out) :: matidx ! the index location in matrix
+    integer, intent(in)  :: i      ! current x index
+    integer, intent(in)  :: j      ! current y index
+    integer, intent(in)  :: k      ! current z index
+    integer, intent(in)  :: g      ! current group index
+    integer, intent(in)  :: nx     ! maximum number of x cells
+    integer, intent(in)  :: ny     ! maximum number of y cells
+    integer, intent(in)  :: ng     ! maximum number of groups
 
     ! Check if coremap is used
     if (cmfd_coremap) then
@@ -396,15 +401,15 @@ contains
 
   subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
 
-    integer :: i     ! iteration counter for x
-    integer :: j     ! iteration counter for y
-    integer :: k     ! iteration counter for z
-    integer :: g     ! iteration counter for groups
-    integer :: irow  ! iteration counter over row (0 reference)
-    integer :: nx    ! maximum number of x cells
-    integer :: ny    ! maximum number of y cells
-    integer :: nz    ! maximum number of z cells
-    integer :: ng    ! maximum number of groups
+    integer, intent(out) :: i     ! iteration counter for x
+    integer, intent(out) :: j     ! iteration counter for y
+    integer, intent(out) :: k     ! iteration counter for z
+    integer, intent(out) :: g     ! iteration counter for groups
+    integer, intent(in)  :: irow  ! iteration counter over row (0 reference)
+    integer, intent(in)  :: nx    ! maximum number of x cells
+    integer, intent(in)  :: ny    ! maximum number of y cells
+    integer, intent(in)  :: nz    ! maximum number of z cells
+    integer, intent(in)  :: ng    ! maximum number of groups
 
     ! Check for core map
     if (cmfd_coremap) then

src/cmfd_power_solver.F90

@@ -41,11 +41,11 @@ contains
 
     use global,  only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
 
-    real(8), optional :: k_tol    ! tolerance on keff
-    real(8), optional :: s_tol    ! tolerance on source
-    logical, optional :: adjoint  ! adjoint calc
+    real(8), intent(in), optional :: k_tol    ! tolerance on keff
+    real(8), intent(in), optional :: s_tol    ! tolerance on source
+    logical, intent(in), optional :: adjoint  ! adjoint calc
 
-    logical :: physical_adjoint = .false.
+    logical :: physical_adjoint = .false. ! physical adjoint default false
 
     ! Set tolerances if present
     if (present(k_tol)) ktol = k_tol
@@ -102,9 +102,11 @@ contains
   subroutine init_data(adjoint)
 
     use constants, only: ONE, ZERO
+#ifdef PETSC
     use global,    only: cmfd_write_matrices
+#endif
 
-    logical :: adjoint
+    logical, intent(in) :: adjoint ! adjoint calcualtion
 
     integer :: n      ! problem size
     real(8) :: guess  ! initial guess
@@ -238,7 +240,7 @@ contains
     use global,     only: cmfd_power_monitor, master
     use, intrinsic :: ISO_FORTRAN_ENV
 
-    integer     :: iter           ! iteration number
+    integer, intent(in) :: iter ! iteration number
 
     ! Reset convergence flag
     iconv = .false.

src/cmfd_prod_operator.F90

@@ -16,7 +16,7 @@ contains
 
   subroutine init_prod_matrix(prod_matrix)
 
-    type(Matrix) :: prod_matrix
+    type(Matrix), intent(inout) :: prod_matrix ! production matrix
 
     integer :: nx   ! maximum number of x cells
     integer :: ny   ! maximum number of y cells
@@ -50,8 +50,8 @@ contains
 
   subroutine build_prod_matrix(prod_matrix, adjoint)
 
-    type(Matrix) :: prod_matrix
-    logical, optional :: adjoint
+    type(Matrix), intent(inout)   :: prod_matrix ! production matrix
+    logical, intent(in), optional :: adjoint     ! adjoint calculation logical
 
     integer :: i                  ! iteration counter for x
     integer :: j                  ! iteration counter for y
@@ -135,14 +135,14 @@ contains
 
   subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
 
-    integer :: matidx ! the index location in matrix
-    integer :: i      ! current x index
-    integer :: j      ! current y index
-    integer :: k      ! current z index
-    integer :: g      ! current group index
-    integer :: nx     ! maximum number of x cells
-    integer :: ny     ! maximum number of y cells
-    integer :: ng     ! maximum number of groups
+    integer, intent(out) :: matidx ! the index location in matrix
+    integer, intent(in)  :: i      ! current x index
+    integer, intent(in)  :: j      ! current y index
+    integer, intent(in)  :: k      ! current z index
+    integer, intent(in)  :: g      ! current group index
+    integer, intent(in)  :: nx     ! maximum number of x cells
+    integer, intent(in)  :: ny     ! maximum number of y cells
+    integer, intent(in)  :: ng     ! maximum number of groups
 
     ! check if coremap is used
     if (cmfd_coremap) then
@@ -165,15 +165,15 @@ contains
 
   subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
 
-    integer :: i     ! iteration counter for x
-    integer :: j     ! iteration counter for y
-    integer :: k     ! iteration counter for z
-    integer :: g     ! iteration counter for groups
-    integer :: irow  ! iteration counter over row (0 reference)
-    integer :: nx    ! maximum number of x cells
-    integer :: ny    ! maximum number of y cells
-    integer :: nz    ! maximum number of z cells
-    integer :: ng    ! maximum number of groups
+    integer, intent(out) :: i    ! iteration counter for x
+    integer, intent(out) :: j    ! iteration counter for y
+    integer, intent(out) :: k    ! iteration counter for z
+    integer, intent(out) :: g    ! iteration counter for groups
+    integer, intent(in)  :: irow ! iteration counter over row (0 reference)
+    integer, intent(in)  :: nx   ! maximum number of x cells
+    integer, intent(in)  :: ny   ! maximum number of y cells
+    integer, intent(in)  :: nz   ! maximum number of z cells
+    integer, intent(in)  :: ng   ! maximum number of groups
 
     ! check for core map
     if (cmfd_coremap) then

src/constants.F90

@@ -8,7 +8,7 @@ module constants
   ! OpenMC major, minor, and release numbers
   integer, parameter :: VERSION_MAJOR   = 0
   integer, parameter :: VERSION_MINOR   = 5
-  integer, parameter :: VERSION_RELEASE = 2
+  integer, parameter :: VERSION_RELEASE = 3
 
   ! Revision numbers for binary files
   integer, parameter :: REVISION_STATEPOINT       = 10
@@ -20,7 +20,7 @@ module constants
        FILETYPE_PARTICLE_RESTART = -2
 
   ! ============================================================================
-  ! ADJUSTABLE PARAMETERS 
+  ! ADJUSTABLE PARAMETERS
 
   ! NOTE: This is the only section of the constants module that should ever be
   ! adjusted. Modifying constants in other sections may cause the code to fail.
@@ -50,6 +50,10 @@ module constants
   integer, parameter :: MAX_WORD_LEN = 150
   integer, parameter :: MAX_FILE_LEN = 255
 
+  ! Maximum number of external source spatial resamples to encounter before an
+  ! error is thrown.
+  integer, parameter :: MAX_EXTSRC_RESAMPLES = 10000
+
   ! ============================================================================
   ! PHYSICAL CONSTANTS
 
@@ -110,9 +114,9 @@ module constants
 
   ! 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_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
@@ -143,7 +147,7 @@ module constants
        ELECTRON = 3
 
   ! Angular distribution type
-  integer, parameter :: & 
+  integer, parameter :: &
        ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
        ANGLE_32_EQUI   = 2, & ! 32 equiprobable bins
        ANGLE_TABULAR   = 3    ! Tabular angular distribution
@@ -151,7 +155,8 @@ module constants
   ! Secondary energy mode for S(a,b) inelastic scattering
   integer, parameter :: &
        SAB_SECONDARY_EQUAL  = 0, & ! Equally-likely outgoing energy bins
-       SAB_SECONDARY_SKEWED = 1    ! Skewed outgoing energy bins
+       SAB_SECONDARY_SKEWED = 1, & ! Skewed outgoing energy bins
+       SAB_SECONDARY_CONT   = 2    ! Continuous, linear-linear interpolation
 
   ! Elastic mode for S(a,b) elastic scattering
   integer, parameter :: &
@@ -275,7 +280,7 @@ module constants
        SCORE_KAPPA_FISSION = -12, & ! fission energy production rate
        SCORE_CURRENT       = -13, & ! partial current
        SCORE_EVENTS        = -14    ! number of events
-       
+
   ! Maximum scattering order supported
   integer, parameter :: SCATT_ORDER_MAX = 10
   character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10"
@@ -322,7 +327,7 @@ module constants
 
   ! Source angular distribution types
   integer, parameter :: &
-       SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular 
+       SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular
        SRC_ANGLE_MONO      = 2, & ! Monodirectional source
        SRC_ANGLE_TABULAR   = 3    ! Tabular distribution
 
@@ -332,7 +337,7 @@ module constants
        SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum
        SRC_ENERGY_WATT    = 3, & ! Watt fission spectrum
        SRC_ENERGY_TABULAR = 4    ! Tabular distribution
-       
+
   ! ============================================================================
   ! MISCELLANEOUS CONSTANTS
 

src/cross_section.F90

@@ -93,6 +93,8 @@ contains
       ! Calculate microscopic cross section for this nuclide
       if (p % E /= micro_xs(i_nuclide) % last_E) then
         call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
+      else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
+        call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
       end if
 
       ! ========================================================================
@@ -235,13 +237,8 @@ contains
       end if
     end if
 
-    ! Set last evaluated energy -- if we're in S(a,b) region, force
-    ! re-calculation of cross-section
-    if (i_sab == 0) then
-      micro_xs(i_nuclide) % last_E = E
-    else
-      micro_xs(i_nuclide) % last_E = ZERO
-    end if
+    micro_xs(i_nuclide) % last_E = E
+    micro_xs(i_nuclide) % last_index_sab = i_sab
 
   end subroutine calculate_nuclide_xs
 

src/doppler.F90

@@ -87,7 +87,7 @@ contains
        ! =======================================================================
        ! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
 
-       if (k == 1 .and. a >= 4.0) then
+       if (k == 1 .and. a >= -4.0) then
           ! Since x = 0, this implies that a = -y
           F_b = F_a
           a = -y
@@ -101,7 +101,7 @@ contains
        end if
 
        ! =======================================================================
-       ! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
+       ! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
 
        k = i
        b = ZERO

src/eigenvalue.F90

@@ -166,7 +166,7 @@ contains
 
   subroutine finalize_generation()
 
-#ifdef OPENMP
+#ifdef _OPENMP
     ! Join the fission bank from each thread into one global fission bank
     call join_bank_from_threads()
 #endif
@@ -821,7 +821,7 @@ contains
 
   end subroutine replay_batch_history
 
-#ifdef OPENMP
+#ifdef _OPENMP
 !===============================================================================
 ! JOIN_BANK_FROM_THREADS
 !===============================================================================

src/global.F90

@@ -158,7 +158,7 @@ module global
   ! Source and fission bank
   type(Bank), allocatable, target :: source_bank(:)
   type(Bank), allocatable, target :: fission_bank(:)
-#ifdef OPENMP
+#ifdef _OPENMP
   type(Bank), allocatable, target :: master_fission_bank(:)
 #endif
   integer(8) :: n_bank       ! # of sites in fission bank
@@ -205,7 +205,7 @@ module global
   integer :: MPI_BANK              ! MPI datatype for fission bank
   integer :: MPI_TALLYRESULT       ! MPI datatype for TallyResult
 
-#ifdef OPENMP
+#ifdef _OPENMP
   integer :: n_threads = NONE      ! number of OpenMP threads
   integer :: thread_id             ! ID of a given thread
 #endif
@@ -283,6 +283,10 @@ module global
   integer    :: trace_gen
   integer(8) :: trace_particle
 
+  ! Particle tracks
+  logical :: write_all_tracks = .false.
+  integer, allocatable :: track_identifiers(:,:)
+
   ! Particle restart run
   logical :: particle_restart_run = .false.
 
@@ -442,7 +446,7 @@ contains
 !$omp parallel
     if (allocated(fission_bank)) deallocate(fission_bank)
 !$omp end parallel
-#ifdef OPENMP
+#ifdef _OPENMP
     if (allocated(master_fission_bank)) deallocate(master_fission_bank)
 #endif
     if (allocated(source_bank)) deallocate(source_bank)
@@ -459,6 +463,9 @@ contains
     call active_tracklength_tallies % clear()
     call active_current_tallies % clear()
     call active_tallies % clear()
+
+    ! Deallocate track_identifiers
+    if (allocated(track_identifiers)) deallocate(track_identifiers)
     
     ! Deallocate dictionaries
     call cell_dict % clear()

src/hdf5_summary.F90

@@ -100,6 +100,7 @@ contains
 
     integer          :: i, j, k, m
     integer          :: n_x, n_y, n_z
+    integer          :: length(3)
     integer, allocatable :: lattice_universes(:,:,:)
     type(Cell),     pointer :: c => null()
     type(Surface),  pointer :: s => null()
@@ -312,8 +313,8 @@ contains
           end do
         end do
       end do
-      call su % write_data(lattice_universes, "universes", &
-           length=(/n_x, n_y, n_z/), &
+      length = [n_x, n_y, n_z]
+      call su % write_data(lattice_universes, "universes", length=length, &
            group="geometry/lattices/lattice " // trim(to_str(lat % id)))
       deallocate(lattice_universes)
 

src/initialize.F90

@@ -26,7 +26,7 @@ module initialize
   use mpi
 #endif
 
-#ifdef OPENMP
+#ifdef _OPENMP
   use omp_lib
 #endif
 
@@ -372,7 +372,7 @@ contains
           ! Read number of threads
           i = i + 1
 
-#ifdef OPENMP          
+#ifdef _OPENMP          
           ! Read and set number of OpenMP threads
           n_threads = str_to_int(argv(i))
           if (n_threads < 1) then
@@ -394,6 +394,9 @@ contains
         case ('-eps_tol', '-ksp_gmres_restart')
           ! Handle options that would be based to PETSC
           i = i + 1
+        case ('-t', '-track', '--track')
+          write_all_tracks = .true.
+          i = i + 1
         case default
           message = "Unknown command line option: " // argv(i)
           call fatal_error()
@@ -827,14 +830,15 @@ contains
       call fatal_error()
     end if
 
-#ifdef OPENMP
+#ifdef _OPENMP
     ! If OpenMP is being used, each thread needs its own private fission
     ! bank. Since the private fission banks need to be combined at the end of a
     ! generation, there is also a 'master_fission_bank' that is used to collect
     ! the sites from each thread.
 
+    n_threads = omp_get_max_threads()
+
 !$omp parallel
-    n_threads = omp_get_num_threads()
     thread_id = omp_get_thread_num()
 
     if (thread_id == 0) then

src/matrix_header.F90

@@ -46,9 +46,9 @@ contains
 
   subroutine matrix_create(self, n, nnz)
 
-    integer       :: n
-    integer       :: nnz
-    class(Matrix) :: self
+    integer, intent(in) :: n             ! dimension of matrix
+    integer, intent(in) :: nnz           ! number of nonzeros
+    class(Matrix), intent(inout) :: self ! matrix instance
 
     ! Preallocate vectors
     if (.not.allocated(self % row)) allocate(self % row(n+1))
@@ -74,7 +74,7 @@ contains
 
   subroutine matrix_destroy(self)
 
-    class(Matrix) :: self
+    class(Matrix), intent(inout) :: self ! matrix instance
 
 #ifdef PETSC
     if (self % petsc_active) call MatDestroy(self % petsc_mat, petsc_err)
@@ -92,9 +92,9 @@ contains
 
   subroutine matrix_add_value(self, col, val)
 
-    integer       :: col
-    real(8)       :: val
-    class(Matrix) :: self
+    integer, intent(in) :: col           ! col location in matrix
+    real(8), intent(in) :: val           ! value to store in matrix
+    class(Matrix), intent(inout) :: self ! matrix instance
 
     ! Record the data
     self % col(self % nz_count) = col
@@ -115,7 +115,7 @@ contains
 
   subroutine matrix_new_row(self)
 
-    class(Matrix) :: self
+    class(Matrix), intent(inout) :: self ! matrix instance
 
     ! Record the current number of nonzeros
     self % row(self % n_count) = self % nz_count
@@ -135,7 +135,7 @@ contains
 
   subroutine matrix_assemble(self)
 
-    class(Matrix) :: self
+    class(Matrix), intent(inout) :: self ! matrix instance
 
     integer :: i
     integer :: first
@@ -161,12 +161,12 @@ contains
 
   recursive subroutine sort_csr(col, val, first, last)
 
-    integer :: col(:)
-    integer :: first
-    integer :: last
-    real(8) :: val(:)
+    integer :: col(:) ! column vector to sort
+    integer :: first  ! first value in sort
+    integer :: last   ! last value in sort
+    real(8) :: val(:) ! value vector to be sorted like columns
 
-    integer :: mid
+    integer :: mid ! midpoint value
 
     if (first < last) then
       call split(col, val, first, last, mid)
@@ -182,18 +182,18 @@ contains
 
   subroutine split(col, val, low, high, mid)
 
-    integer :: col(:)
-    integer :: low
-    integer :: high
-    integer :: mid
-    real(8) :: val(:)
+    integer :: col(:) ! column vector to sort
+    integer :: low    ! low index to sort
+    integer :: high   ! high index to sort
+    integer :: mid    ! middle of sort
+    real(8) :: val(:) ! value vector to be sorted like columns
 
-    integer :: left
-    integer :: right
-    integer :: iswap
-    integer :: pivot
-    real(8) :: rswap
-    real(8) :: val0
+    integer :: left  ! contains left value in sort
+    integer :: right ! contains right value in sort
+    integer :: iswap ! temporary interger swap
+    integer :: pivot ! pivotting variable for columns
+    real(8) :: rswap ! temporary real swap
+    real(8) :: val0  ! pivot for value vector
 
     left = low
     right = high
@@ -235,14 +235,14 @@ contains
   end subroutine split
 
 !===============================================================================
-! MATRIX_GET_ROW
+! MATRIX_GET_ROW is a method to get row and checks for PETSc C indexing use
 !===============================================================================
 
   function matrix_get_row(self, i) result(row)
 
-    class(Matrix) :: self
-    integer       :: i
-    integer       :: row
+    class(Matrix), intent(in) :: self ! matrix instance
+    integer, intent(in)       :: i    ! row to get
+    integer                   :: row  ! index in col where row begins
 
     row = self % row(i)
 
@@ -251,14 +251,14 @@ contains
   end function matrix_get_row
 
 !===============================================================================
-! MATRIX_GET_COL
+! MATRIX_GET_COL is a method to get column and checks for PETSc C index use
 !===============================================================================
 
   function matrix_get_col(self, i) result(col)
 
-    class(Matrix) :: self
-    integer       :: i
-    integer       :: col
+    class(Matrix), intent(in) :: self ! matrix instance
+    integer, intent(in)       :: i    ! index from row vector
+    integer                   :: col  ! the actual column
 
     col = self % col(i)
 
@@ -272,7 +272,7 @@ contains
 
   subroutine matrix_setup_petsc(self)
 
-    class(Matrix) :: self
+    class(Matrix), intent(inout) :: self ! matrix instance
 
     ! change indices to c notation
     self % row = self % row - 1
@@ -295,11 +295,11 @@ contains
 
   subroutine matrix_write_petsc_binary(self, filename)
 
-    character(*) :: filename
-    class(Matrix) :: self
+    character(*), intent(in)  :: filename ! file name to write to
+    class(Matrix), intent(in) :: self     ! matrix instance
 
 #ifdef PETSC
-    type(PetscViewer) :: viewer
+    type(PetscViewer) :: viewer ! a petsc viewer instance
 
     call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
          FILE_MODE_WRITE, viewer, petsc_err)
@@ -315,7 +315,7 @@ contains
 
   subroutine matrix_transpose(self)
 
-    class(Matrix) :: self
+    class(Matrix), intent(inout) :: self ! matrix instance
 
 #ifdef PETSC
     call MatTranspose(self % petsc_mat, MAT_REUSE_MATRIX, self % petsc_mat, &
@@ -333,12 +333,12 @@ contains
     use constants,      only: ZERO
     use vector_header,  only: Vector
 
-    class(Matrix) :: self
-    type(Vector)  :: vec_in
-    type(Vector)  :: vec_out
+    class(Matrix), intent(in)   :: self    ! matrix instance
+    type(Vector), intent(in)    :: vec_in  ! vector to multiply matrix against
+    type(Vector), intent(inout) :: vec_out ! resulting vector
 
-    integer :: i
-    integer :: j
+    integer :: i ! row iteration counter
+    integer :: j ! column iteration counter
 
     ! Begin loop around rows
     ROWS: do i = 1, self % n

src/output.F90

@@ -31,6 +31,10 @@ contains
 
   subroutine title()
 
+#ifdef _OPENMP
+    use omp_lib
+#endif
+
     write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
          '       .d88888b.                             888b     d888  .d8888b.', &
          '      d88P" "Y88b                            8888b   d8888 d88P  Y88b', &
@@ -46,25 +50,31 @@ contains
 
     ! Write version information
     write(UNIT=OUTPUT_UNIT, FMT=*) &
-         '     Copyright:     2011-2013 Massachusetts Institute of Technology'
+         '     Copyright:      2011-2014 Massachusetts Institute of Technology'
     write(UNIT=OUTPUT_UNIT, FMT=*) &
-         '     License:       http://mit-crpg.github.io/openmc/license.html'
-    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
+         '     License:        http://mit-crpg.github.io/openmc/license.html'
+    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') &
          VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
 #ifdef GIT_SHA1
-    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1
+    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1
 #endif
 
     ! Write the date and time
-    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') &
+    write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') &
          time_stamp()
 
 #ifdef MPI
     ! Write number of processors
-    write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') &
+    write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') &
          trim(to_str(n_procs))
 #endif
 
+#ifdef _OPENMP
+    ! Write number of OpenMP threads
+    write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') &
+         trim(to_str(omp_get_max_threads()))
+#endif
+
   end subroutine title
 
 !===============================================================================
@@ -126,7 +136,7 @@ contains
     ! print header based on level
     select case (header_level)
     case (1)
-      write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & 
+      write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
            repeat('=', n) // '>     ' // trim(line) // '     <' // &
            repeat('=', m), repeat('=', 75)
     case (2)
@@ -172,6 +182,7 @@ contains
       write(OUTPUT_UNIT,*) '  -r, --restart          Restart a previous run from a state point'
       write(OUTPUT_UNIT,*) '                         or a particle restart file'
       write(OUTPUT_UNIT,*) '  -s, --threads          Number of OpenMP threads'
+      write(OUTPUT_UNIT,*) '  -t, --track            Write tracks for all particles'
       write(OUTPUT_UNIT,*) '  -v, --version          Show version information'
       write(OUTPUT_UNIT,*) '  -?, --help             Show this message'
     end if
@@ -179,7 +190,7 @@ contains
   end subroutine print_usage
 
 !===============================================================================
-! WRITE_MESSAGE displays an informational message to the log file and the 
+! WRITE_MESSAGE displays an informational message to the log file and the
 ! standard output stream.
 !===============================================================================
 
@@ -214,7 +225,7 @@ contains
           ! Determine last space in current line
           last_space = index(message(i_start+1:i_start+line_wrap), &
                ' ', BACK=.true.)
-          if (last_space == 0) then 
+          if (last_space == 0) then
             i_end = min(length + 1, i_start+line_wrap) - 1
             write(ou, fmt='(1X,A)') message(i_start+1:i_end)
           else
@@ -1027,8 +1038,10 @@ contains
     type(SAlphaBeta), pointer :: sab
     integer,         optional :: unit
 
-    integer :: size_sab ! memory used by S(a,b) table
-    integer :: unit_    ! unit to write to
+    integer :: size_sab   ! memory used by S(a,b) table
+    integer :: unit_      ! unit to write to
+    integer :: i          ! Loop counter for parsing through sab % zaid
+    integer :: char_count ! Counter for the number of characters on a line
 
     ! set default unit for writing information
     if (present(unit)) then
@@ -1039,7 +1052,30 @@ contains
 
     ! Basic S(a,b) table information
     write(unit_,*) 'S(a,b) Table ' // trim(sab % name)
-    write(unit_,*) '  zaid = ' // trim(to_str(sab % zaid))
+    write(unit_,'(A)',advance="no") '   zaids = '
+    ! Initialize the counter based on the above string
+    char_count = 11
+    do i = 1, sab % n_zaid
+      ! Deal with a line thats too long
+      if (char_count >= 73) then  ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma)
+        ! End the line
+        write(unit_,*) ""
+        ! Add 11 leading blanks
+        write(unit_,'(A)', advance="no") "           "
+        ! reset the counter to 11
+        char_count = 11
+      end if
+      if (i < sab % n_zaid) then
+        ! Include a comma
+        write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", "
+        char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2
+      else
+        ! Don't include a comma, since we are all done
+        write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i)))
+      end if
+
+    end do
+    write(unit_,*) "" ! Move to next line
     write(unit_,*) '  awr = ' // trim(to_str(sab % awr))
     write(unit_,*) '  kT = ' // trim(to_str(sab % kT))
 
@@ -1241,7 +1277,7 @@ contains
       end select
     end if
     write(UNIT=ou, FMT=*)
-    
+
     write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "========="
     write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
     if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
@@ -1271,7 +1307,7 @@ contains
     if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
          entropy(overall_gen)
 
-    if (overall_gen - n_inactive*gen_per_batch > 1) then 
+    if (overall_gen - n_inactive*gen_per_batch > 1) then
       write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
            keff, keff_std
     end if
@@ -1299,7 +1335,7 @@ contains
          entropy(current_batch*gen_per_batch)
 
     ! write out accumulated k-effective if after first active batch
-    if (overall_gen - n_inactive*gen_per_batch > 1) then 
+    if (overall_gen - n_inactive*gen_per_batch > 1) then
       write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
            keff, keff_std
     else
@@ -1309,7 +1345,7 @@ contains
     ! write out cmfd keff if it is active and other display info
     if (cmfd_on) then
       write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
-         cmfd % k_cmfd(current_batch) 
+         cmfd % k_cmfd(current_batch)
       select case(trim(cmfd_display))
         case('entropy')
           write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
@@ -1348,7 +1384,7 @@ contains
 
       ! Plot id
       write(ou,100) "Plot ID:", trim(to_str(pl % id))
-      
+
       ! Plot type
       if (pl % type == PLOT_TYPE_SLICE) then
         write(ou,100) "Plot Type:", "Slice"
@@ -1386,7 +1422,7 @@ contains
              trim(to_str(pl % pixels(2)))
       else if (pl % type == PLOT_TYPE_VOXEL) then
         write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // &
-             trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3))) 
+             trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
       end if
 
       write(ou,*)
@@ -1538,7 +1574,7 @@ contains
 
     write(ou,100) 'Cell ID','No. Overlap Checks'
 
-    do i = 1, n_cells    
+    do i = 1, n_cells
       write(ou,101) cells(i) % id, overlap_check_cnt(i)
       if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1
     end do
@@ -1572,7 +1608,7 @@ contains
     integer :: k            ! loop index for scoring bins
     integer :: n            ! loop index for nuclides
     integer :: l            ! loop index for user scores
-    integer :: type         ! type of tally filter        
+    integer :: type         ! type of tally filter
     integer :: indent       ! number of spaces to preceed output
     integer :: filter_index ! index in results array for filters
     integer :: score_index  ! scoring bin index
@@ -1748,13 +1784,13 @@ contains
                 score_name = 'P' // trim(to_str(t % scatt_order(k))) // &
                   ' Scattering Moment'
               end if
-              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & 
+              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
                 repeat(" ", indent), score_name, &
                 to_str(t % results(score_index,filter_index) % sum), &
                 trim(to_str(t % results(score_index,filter_index) % sum_sq))
             else if (t % score_bins(k) == SCORE_SCATTER_PN) then
               score_name = "Scattering Rate"
-              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & 
+              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
                 repeat(" ", indent), score_name, &
                 to_str(t % results(score_index,filter_index) % sum), &
                 trim(to_str(t % results(score_index,filter_index) % sum_sq))
@@ -1762,7 +1798,7 @@ contains
                 score_index = score_index + 1
                 score_name = 'P' // trim(to_str(n_order)) // &
                   ' Scattering Moment'
-                write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & 
+                write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
                   repeat(" ", indent), score_name, &
                   to_str(t % results(score_index,filter_index) % sum), &
                   trim(to_str(t % results(score_index,filter_index) % sum_sq))
@@ -1774,7 +1810,7 @@ contains
               else
                 score_name = score_names(abs(t % score_bins(k)))
               end if
-              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & 
+              write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
                 repeat(" ", indent), score_name, &
                 to_str(t % results(score_index,filter_index) % sum), &
                 trim(to_str(t % results(score_index,filter_index) % sum_sq))
@@ -1812,8 +1848,8 @@ contains
     integer :: i_filter_ein         ! index for incoming energy filter
     integer :: i_filter_surf        ! index for surface filter
     integer :: n                    ! number of incoming energy bins
-    integer :: len1                 ! length of string 
-    integer :: len2                 ! length of string 
+    integer :: len1                 ! length of string
+    integer :: len2                 ! length of string
     integer :: filter_index         ! index in results array for filters
     logical :: print_ebin           ! should incoming energy bin be displayed?
     character(MAX_LINE_LEN) :: string
@@ -1863,14 +1899,14 @@ contains
                  mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_RIGHT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Left", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_RIGHT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Left", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
@@ -1880,14 +1916,14 @@ contains
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_RIGHT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Right", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_RIGHT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Right", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
@@ -1897,14 +1933,14 @@ contains
                  mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_FRONT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Back", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_FRONT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Back", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
@@ -1914,14 +1950,14 @@ contains
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_FRONT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Front", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_FRONT
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Front", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
@@ -1931,14 +1967,14 @@ contains
                  mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_TOP
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Bottom", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_TOP
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Bottom", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
@@ -1948,14 +1984,14 @@ contains
                  mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
             matching_bins(i_filter_surf) = IN_TOP
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Incoming Current from Top", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))
 
             matching_bins(i_filter_surf) = OUT_TOP
             filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
-            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & 
+            write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
                  "Outgoing Current to Top", &
                  to_str(t % results(1,filter_index) % sum), &
                  trim(to_str(t % results(1,filter_index) % sum_sq))

src/output_interface.F90

@@ -1644,6 +1644,12 @@ contains
     call h5dclose_f(dset, hdf5_err)
     if (present(group)) call hdf5_close_group(self % hdf5_grp)
 
+# elif MPI
+
+    ! Write out tally buffer
+    call MPI_FILE_READ(self % unit_fh, buffer, n1*n2, MPI_TALLYRESULT, &
+         MPI_STATUS_IGNORE, mpiio_err)
+
 #else
 
     ! Read tally result

src/particle_header.F90

@@ -76,6 +76,9 @@ module particle_header
     ! Statistical data
     integer    :: n_collision   ! # of collisions
 
+    ! Track output
+    logical    :: write_track = .false.
+
   contains
     procedure :: initialize => initialize_particle
     procedure :: clear => clear_particle

src/physics.F90

@@ -47,7 +47,7 @@ contains
     if (verbosity >= 10 .or. trace) then
       message = "    " // trim(reaction_name(p % event_MT)) // " with " // &
            trim(adjustl(nuclides(p % event_nuclide) % name)) // &
-           ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV." 
+           ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV."
       call write_message()
     end if
 
@@ -158,7 +158,7 @@ contains
         call fatal_error()
       end if
 
-      ! Find atom density 
+      ! Find atom density
       i_nuclide    = mat % nuclide(i)
       atom_density = mat % atom_density(i)
 
@@ -209,7 +209,7 @@ contains
       i_reaction = nuc % index_fission(1)
       return
     end if
-    
+
     ! Get grid index and interpolatoin factor and sample fission cdf
     i_grid = micro_xs(i_nuclide) % index_grid
     f      = micro_xs(i_nuclide) % interp_factor
@@ -226,7 +226,7 @@ contains
       if (i_grid < rxn % threshold) cycle
 
       ! add to cumulative probability
-      prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & 
+      prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
            + f*(rxn%sigma(i_grid - rxn%threshold + 2)))
 
       ! Create fission bank sites if fission occus
@@ -382,7 +382,7 @@ contains
         if (i_grid < rxn % threshold) cycle
 
         ! add to cumulative probability
-        prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & 
+        prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
              + f*(rxn%sigma(i_grid - rxn%threshold + 2)))
       end do
 
@@ -490,13 +490,24 @@ contains
     integer :: k            ! outgoing cosine bin
     integer :: n_energy_out ! number of outgoing energy bins
     real(8) :: f            ! interpolation factor
-    real(8) :: r            ! used for skewed sampling
+    real(8) :: r            ! used for skewed sampling & continuous
     real(8) :: E_ij         ! outgoing energy j for E_in(i)
     real(8) :: E_i1j        ! outgoing energy j for E_in(i+1)
     real(8) :: mu_ijk       ! outgoing cosine k for E_in(i) and E_out(j)
     real(8) :: mu_i1jk      ! outgoing cosine k for E_in(i+1) and E_out(j)
     real(8) :: prob         ! probability for sampling Bragg edge
     type(SAlphaBeta), pointer, save :: sab => null()
+    ! Following are needed only for SAB_SECONDARY_CONT scattering
+    integer :: l              ! sampled incoming E bin (is i or i + 1)
+    real(8) :: E_i_1, E_i_J   ! endpoints on outgoing grid i
+    real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1
+    real(8) :: E_1, E_J       ! endpoints interpolated between i and i+1
+    real(8) :: E_l_j, E_l_j1  ! adjacent E on outgoing grid l
+    real(8) :: p_l_j, p_l_j1  ! adjacent p on outgoing grid l
+    real(8) :: c_j, c_j1      ! cumulative probability
+    real(8) :: frac           ! interpolation factor on outgoing energy
+    real(8) :: r1             ! RNG for outgoing energy
+
 !$omp threadprivate(sab)
 
     ! Get pointer to S(a,b) table
@@ -513,7 +524,7 @@ contains
         f = ZERO
       else
         i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
-        f = (E - sab%elastic_e_in(i)) / & 
+        f = (E - sab%elastic_e_in(i)) / &
              (sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
       end if
 
@@ -554,8 +565,7 @@ contains
       ! Outgoing energy is same as incoming energy -- no need to do anything
 
     else
-      ! Determine number of outgoing energy and angle bins
-      n_energy_out = sab % n_inelastic_e_out
+      ! Perform inelastic calculations
 
       ! Get index and interpolation factor for inelastic grid
       if (E < sab % inelastic_e_in(1)) then
@@ -563,61 +573,150 @@ contains
         f = ZERO
       else
         i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
-        f = (E - sab%inelastic_e_in(i)) / & 
+        f = (E - sab%inelastic_e_in(i)) / &
              (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
       end if
 
       ! Now that we have an incoming energy bin, we need to determine the
       ! outgoing energy bin. This will depend on the "secondary energy
       ! mode". If the mode is 0, then the outgoing energy bin is chosen from a
-      ! set of equally-likely bins. However, if the mode is 1, then the first
+      ! set of equally-likely bins. If the mode is 1, then the first
       ! two and last two bins are skewed to have lower probabilities than the
       ! other bins (0.1 for the first and last bins and 0.4 for the second and
-      ! second to last bins, relative to a normal bin probability of 1)
+      ! second to last bins, relative to a normal bin probability of 1).
+      ! Finally, if the mode is 2, then a continuous distribution (with
+      ! accompanying PDF and CDF is utilized)
 
-      if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
-        ! All bins equally likely
-        j = 1 + int(prn() * n_energy_out)
-      elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
-        r = prn() * (n_energy_out - 3)
-        if (r > ONE) then
-          ! equally likely N-4 middle bins
-          j = int(r) + 2
-        elseif (r > 0.6) then
-          ! second to last bin has relative probability of 0.4
-          j = n_energy_out - 1
-        elseif (r > 0.5) then
-          ! last bin has relative probability of 0.1
-          j = n_energy_out
-        elseif (r > 0.1) then
-          ! second bin has relative probability of 0.4
-          j = 2
-        else
-          ! first bin has relative probability of 0.1
-          j = 1
+      if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
+          (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
+        if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
+          ! All bins equally likely
+
+          j = 1 + int(prn() * sab % n_inelastic_e_out)
+        elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
+          ! Distribution skewed away from edge points
+
+          ! Determine number of outgoing energy and angle bins
+          n_energy_out = sab % n_inelastic_e_out
+
+          r = prn() * (n_energy_out - 3)
+          if (r > ONE) then
+            ! equally likely N-4 middle bins
+            j = int(r) + 2
+          elseif (r > 0.6) then
+            ! second to last bin has relative probability of 0.4
+            j = n_energy_out - 1
+          elseif (r > 0.5) then
+            ! last bin has relative probability of 0.1
+            j = n_energy_out
+          elseif (r > 0.1) then
+            ! second bin has relative probability of 0.4
+            j = 2
+          else
+            ! first bin has relative probability of 0.1
+            j = 1
+          end if
         end if
+
+        ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
+        E_ij  = sab % inelastic_e_out(j,i)
+        E_i1j = sab % inelastic_e_out(j,i+1)
+
+        ! Outgoing energy
+        E = (1 - f)*E_ij + f*E_i1j
+
+        ! Sample outgoing cosine bin
+        k = 1 + int(prn() * sab % n_inelastic_mu)
+
+        ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
+        mu_ijk  = sab % inelastic_mu(k,j,i)
+        mu_i1jk = sab % inelastic_mu(k,j,i+1)
+
+        ! Cosine of angle between incoming and outgoing neutron
+        mu = (1 - f)*mu_ijk + f*mu_i1jk
+
+      else if (sab % secondary_mode == SAB_SECONDARY_CONT) then
+        ! Continuous secondary energy - this is to be similar to
+        ! Law 61 interpolation on outgoing energy
+
+        ! Sample between ith and (i+1)th bin
+        r = prn()
+        if (f > r) then
+          l = i + 1
+        else
+          l = i
+        end if
+
+        ! Determine endpoints on grid i
+        n_energy_out = sab % inelastic_data(i) % n_e_out
+        E_i_1 = sab % inelastic_data(i) % e_out(1)
+        E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
+
+        ! Determine endpoints on grid i + 1
+        n_energy_out = sab % inelastic_data(i + 1) % n_e_out
+        E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
+        E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
+
+        E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
+        E_J = E_i_J + f * (E_i1_J - E_i_J)
+
+        ! Determine outgoing energy bin
+        ! (First reset n_energy_out to the right value)
+        n_energy_out = sab % inelastic_data(l) % n_e_out
+        r1 = prn()
+        c_j = sab % inelastic_data(l) % e_out_cdf(1)
+        do j = 1, n_energy_out - 1
+          c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
+          if (r1 < c_j1) exit
+          c_j = c_j1
+        end do
+
+        ! check to make sure k is <= n_energy_out - 1
+        j = min(j, n_energy_out - 1)
+
+        ! Get the data to interpolate between
+        E_l_j = sab % inelastic_data(l) % e_out(j)
+        p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
+
+        ! Next part assumes linear-linear interpolation in standard
+        E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
+        p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
+
+        ! Find secondary energy (variable E)
+        frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
+        if (frac == ZERO) then
+          E = E_l_j + (r1 - c_j) / p_l_j
+        else
+          E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
+                       TWO * frac * (r1 - c_j))) - p_l_j) / frac
+        end if
+
+        ! Now interpolate between incident energy bins i and i + 1
+        if (l == i) then
+          E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
+        else
+          E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
+        end if
+
+        ! Find angular distribution for closest outgoing energy bin
+        if (r1 - c_j < c_j1 - r1) then
+          j = j
+        else
+          j = j + 1
+        end if
+
+        ! Sample outgoing cosine bin
+        k = 1 + int(prn() * sab % n_inelastic_mu)
+
+        ! Will use mu from the randomly chosen incoming and closest outgoing
+        ! energy bins
+        mu = sab % inelastic_data(l) % mu(k, j)
+
       else
         message = "Invalid secondary energy mode on S(a,b) table " // &
              trim(sab % name)
-      end if
-
-      ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
-      E_ij  = sab % inelastic_e_out(j,i)
-      E_i1j = sab % inelastic_e_out(j,i+1)
-
-      ! Outgoing energy
-      E = (1 - f)*E_ij + f*E_i1j
-
-      ! Sample outgoing cosine bin
-      k = 1 + int(prn() * sab % n_inelastic_mu)
-
-      ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
-      mu_ijk  = sab % inelastic_mu(k,j,i)
-      mu_i1jk = sab % inelastic_mu(k,j,i+1)
-
-      ! Cosine of angle between incoming and outgoing neutron
-      mu = (1 - f)*mu_ijk + f*mu_i1jk
-    end if
+      end if  ! (inelastic secondary energy treatment)
+    end if  ! (elastic or inelastic)
 
     ! change direction of particle
     uvw = rotate_angle(uvw, mu)
@@ -974,7 +1073,7 @@ contains
       E_cm = E
 
       ! determine outgoing energy in lab
-      E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & 
+      E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
            / ((A+ONE)*(A+ONE))
 
       ! determine outgoing angle in lab
@@ -1037,7 +1136,7 @@ contains
       r = ONE
     else
       i = binary_search(rxn % adist % energy, n, E)
-      r = (E - rxn % adist % energy(i)) / & 
+      r = (E - rxn % adist % energy(i)) / &
            (rxn % adist % energy(i+1) - rxn % adist % energy(i))
     end if
 
@@ -1166,7 +1265,7 @@ contains
     end if
 
   end function rotate_angle
-    
+
 !===============================================================================
 ! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary
 ! neutron from a collision or for a prompt/delayed fission neutron
@@ -1322,7 +1421,7 @@ contains
       ! =======================================================================
       ! CONTINUOUS TABULAR DISTRIBUTION
 
-      ! read number of interpolation regions and incoming energies 
+      ! read number of interpolation regions and incoming energies
       NR  = int(edist % data(1))
       NE  = int(edist % data(2 + 2*NR))
       if (NR == 1) then
@@ -1348,7 +1447,7 @@ contains
         r = ONE
       else
         i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
-        r = (E_in - edist%data(lc+i)) / & 
+        r = (E_in - edist%data(lc+i)) / &
              (edist%data(lc+i+1) - edist%data(lc+i))
       end if
 
@@ -1452,7 +1551,7 @@ contains
       ! =======================================================================
       ! MAXWELL FISSION SPECTRUM
 
-      ! read number of interpolation regions and incoming energies 
+      ! read number of interpolation regions and incoming energies
       NR = int(edist % data(1))
       NE = int(edist % data(2 + 2*NR))
 
@@ -1484,7 +1583,7 @@ contains
       ! =======================================================================
       ! EVAPORATION SPECTRUM
 
-      ! read number of interpolation regions and incoming energies 
+      ! read number of interpolation regions and incoming energies
       NR = int(edist % data(1))
       NE = int(edist % data(2 + 2*NR))
 
@@ -1565,7 +1664,7 @@ contains
         call fatal_error()
       end if
 
-      ! read number of interpolation regions and incoming energies 
+      ! read number of interpolation regions and incoming energies
       NR = int(edist % data(1))
       NE = int(edist % data(2 + 2*NR))
       if (NR > 0) then
@@ -1587,7 +1686,7 @@ contains
         r = ONE
       else
         i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
-        r = (E_in - edist%data(lc+i)) / & 
+        r = (E_in - edist%data(lc+i)) / &
              (edist%data(lc+i+1) - edist%data(lc+i))
       end if
 
@@ -1714,11 +1813,11 @@ contains
 
       if (.not. present(mu_out)) then
         ! call write_particle_restart()
-        message = "Law 44 called without giving mu_out as argument."
+        message = "Law 61 called without giving mu_out as argument."
         call fatal_error()
       end if
 
-      ! read number of interpolation regions and incoming energies 
+      ! read number of interpolation regions and incoming energies
       NR = int(edist % data(1))
       NE = int(edist % data(2 + 2*NR))
       if (NR > 0) then
@@ -1740,7 +1839,7 @@ contains
         r = ONE
       else
         i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
-        r = (E_in - edist%data(lc+i)) / & 
+        r = (E_in - edist%data(lc+i)) / &
              (edist%data(lc+i+1) - edist%data(lc+i))
       end if
 

src/random_lcg.F90

@@ -5,8 +5,13 @@ module random_lcg
   private
   save
 
+  ! Random number streams
+  integer, parameter :: N_STREAMS = 2
+  integer, parameter :: STREAM_TRACKING = 1
+  integer, parameter :: STREAM_TALLIES  = 2
+
   integer(8) :: prn_seed0  ! original seed
-  integer(8) :: prn_seed   ! current seed
+  integer(8) :: prn_seed(N_STREAMS) ! current seed
   integer(8) :: prn_mult   ! multiplication factor, g
   integer(8) :: prn_add    ! additive factor, c
   integer    :: prn_bits   ! number of bits, M
@@ -14,6 +19,7 @@ module random_lcg
   integer(8) :: prn_mask   ! 2^M - 1
   integer(8) :: prn_stride ! stride between particles
   real(8)    :: prn_norm   ! 2^(-M)
+  integer    :: stream     ! current RNG stream
 
 !$omp threadprivate(prn_seed)
 
@@ -21,6 +27,8 @@ module random_lcg
   public :: initialize_prng
   public :: set_particle_seed
   public :: prn_skip
+  public :: prn_set_stream
+  public :: STREAM_TRACKING, STREAM_TALLIES
 
 contains
 
@@ -35,12 +43,12 @@ contains
     ! This algorithm uses bit-masking to find the next integer(8) value to be
     ! used to calculate the random number
 
-    prn_seed = iand(prn_mult*prn_seed + prn_add, prn_mask)
+    prn_seed(stream) = iand(prn_mult*prn_seed(stream) + prn_add, prn_mask)
 
     ! Once the integer is calculated, we just need to divide by 2**m,
     ! represented here as multiplying by a pre-calculated factor
 
-    pseudo_rn = prn_seed * prn_norm
+    pseudo_rn = prn_seed(stream) * prn_norm
 
   end function prn
 
@@ -53,8 +61,13 @@ contains
 
     use global, only: seed
 
+    integer :: i
+
+    stream     = STREAM_TRACKING
     prn_seed0  = seed
-    prn_seed   = prn_seed
+    do i = 1, N_STREAMS
+      prn_seed(i) = prn_seed0 + i - 1
+    end do
     prn_mult   = 2806196910506780709_8
     prn_add    = 1_8
     prn_bits   = 63
@@ -74,7 +87,11 @@ contains
 
     integer(8), intent(in) :: id
 
-    prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0)
+    integer :: i
+
+    do i = 1, N_STREAMS
+      prn_seed(i) = prn_skip_ahead(id*prn_stride, prn_seed0 + i - 1)
+    end do
 
   end subroutine set_particle_seed
     
@@ -86,7 +103,7 @@ contains
 
     integer(8), intent(in) :: n ! number of seeds to skip
 
-    prn_seed = prn_skip_ahead(n, prn_seed)
+    prn_seed(stream) = prn_skip_ahead(n, prn_seed(stream))
 
   end subroutine prn_skip
 
@@ -151,4 +168,18 @@ contains
 
   end function prn_skip_ahead
 
+!===============================================================================
+! PRN_SET_STREAM changes the random number stream. If random numbers are needed
+! in routines not used directly for tracking (e.g. physics), this allows the
+! numbers to be generated without affecting reproducibility of the physics.
+!===============================================================================
+
+  subroutine prn_set_stream(i)
+
+    integer, intent(in) :: i
+
+    stream = i
+
+  end subroutine prn_set_stream
+
 end module random_lcg

src/relaxng/cross_sections.rnc

@@ -10,7 +10,9 @@ element cross_sections {
     (element temperature { xsd:double } | attribute temperature { xsd:double }) &
     (element path { xsd:string { maxLength = "255" } } |
       attribute path { xsd:string { maxLength = "255" } }) &
-    (element location { xsd:int } | attribute location { xsd:int })?
+    (element location { xsd:int } | attribute location { xsd:int })? &
+    (element filetype { ( "ascii" | "binary" ) } |
+      attribute filetype { ( "ascii" | "binary" ) })?
   }* &
 
   element directory { xsd:string { maxLength = "255" } }? &

src/relaxng/settings.rnc

@@ -40,8 +40,12 @@ element settings {
 
   element no_reduce { xsd:boolean }? &
 
-  element output { list { 
-    ( "summary" | "cross_sections" | "tallies" )+ } }? &
+  element output {
+    (element summary { xsd:boolean } | attribute summary { xsd:boolean })? &
+    (element cross_sections { xsd:boolean } |
+      attribute cross_sections { xsd:boolean })? &
+    (element tallies { xsd:boolean } | attribute tallies { xsd:boolean })?
+  }? &
 
   element output_path { xsd:string { maxLength = "255" } }? &
 
@@ -101,6 +105,8 @@ element settings {
 
   element trace { list { xsd:positiveInteger+ } }? &
 
+  element track { list { xsd:positiveInteger+ } }? &
+
   element verbosity { xsd:positiveInteger }? &
 
   element uniform_fs{

src/solver_interface.F90

@@ -16,8 +16,8 @@ module solver_interface
   ! GMRES solver type 
   type, public :: GMRESSolver 
 #ifdef PETSC
-    type(ksp) :: ksp_
-    type(pc)  :: pc_
+    type(ksp) :: ksp_ ! Krylov linear solver instance
+    type(pc)  :: pc_  ! Preconditioner instance
 #endif
    contains
 #ifdef PETSC
@@ -37,11 +37,11 @@ module solver_interface
   ! JFNK solver type 
   type, public :: JFNKSolver 
 #ifdef PETSC
-    type(ksp)  :: ksp_
-    type(pc)   :: pc_
-    type(snes) :: snes_
-    type(mat)  :: jac_mf
-    integer    :: ls
+    type(ksp)  :: ksp_   ! Krylov linear solver instance
+    type(pc)   :: pc_    ! Preconditioner instance
+    type(snes) :: snes_  ! Nonlinear solver instance
+    type(mat)  :: jac_mf ! Matrix free jacobian instance
+    integer    :: ls     ! Line search instance
 #endif
    contains
 #ifdef PETSC
@@ -56,13 +56,13 @@ module solver_interface
   abstract interface
     subroutine res_interface(x, res)
       import :: Vector
-      type(Vector) :: x
-      type(Vector) :: res
+      type(Vector), intent(in)    :: x   ! solution vector
+      type(Vector), intent(inout) :: res ! residual vector
     end subroutine res_interface
 
     subroutine jac_interface(x)
       import :: Vector
-      type(Vector) :: x
+      type(Vector), intent(in) :: x ! solution vector
     end subroutine jac_interface
   end interface
 
@@ -79,7 +79,7 @@ contains
 
   subroutine petsc_gmres_create(self)
 
-    class(GMRESSolver) :: self
+    class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
 
     integer :: ilu_levels = 5
     real(8) :: rtol = 1.0e-10_8
@@ -102,9 +102,9 @@ contains
 
   subroutine petsc_gmres_set_oper(self, prec_mat, mat_in)
 
-    class(GMRESSolver) :: self
-    type(Matrix)       :: prec_mat
-    type(Matrix)       :: mat_in
+    class(GMRESSolver), intent(inout) :: self     ! GMRES solver instanace
+    type(Matrix), intent(inout)       :: prec_mat ! preconditioner matrix
+    type(Matrix), intent(inout)       :: mat_in   ! coefficient matrix
 
     call KSPSetOperators(self % ksp_, mat_in % petsc_mat, prec_mat % petsc_mat, &
          SAME_NONZERO_PATTERN, petsc_err)
@@ -118,7 +118,7 @@ contains
 
   subroutine petsc_gmres_destroy(self)
 
-   class(GMRESSolver) :: self
+   class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
 
     call KSPDestroy(self % ksp_, petsc_err)
 
@@ -130,9 +130,9 @@ contains
 
   subroutine petsc_gmres_solve(self, b, x)
 
-    class(GMRESSolver) :: self
-    type(Vector)       :: b
-    type(Vector)       :: x
+    class(GMRESSolver), intent(inout) :: self ! GMRES solver instance
+    type(Vector), intent(inout)       :: b    ! right hand side vector
+    type(Vector), intent(inout)       :: x    ! solution vector
 
     call KSPSolve(self % ksp_, b % petsc_vec, x % petsc_vec, petsc_err)
 
@@ -144,7 +144,7 @@ contains
 
   subroutine petsc_jfnk_create(self)
 
-    class(JFNKSolver) :: self
+    class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
 
     ! Turn on mf_operator option for matrix free jacobian
     call PetscOptionsSetValue("-snes_mf_operator", "TRUE", petsc_err)
@@ -176,7 +176,7 @@ contains
 
   subroutine petsc_jfnk_destroy(self)
 
-    class(JFNKSolver) :: self
+    class(JFNKSolver), intent(inout) :: self ! JFNK solver instance
 
     call MatDestroy(self % jac_mf, petsc_err)
     call SNESDestroy(self % snes_, petsc_err)
@@ -189,10 +189,10 @@ contains
 
   subroutine petsc_jfnk_set_functions(self, ctx, res, jac_prec)
 
-    class(JFNKSolver) :: self
-    type(Jfnk_ctx)    :: ctx
-    type(Vector)      :: res
-    type(Matrix)      :: jac_prec
+    class(JFNKSolver), intent(inout) :: self     ! JFNK solver instance
+    type(Jfnk_ctx),    intent(inout) :: ctx      ! JFNK context instance
+    type(Vector),      intent(inout) :: res      ! residual vector
+    type(Matrix),      intent(inout) :: jac_prec ! preconditioner matrix
 
     ! Set residual procedure
     call SNESSetFunction(self % snes_, res % petsc_vec, &
@@ -217,8 +217,8 @@ contains
 
   subroutine petsc_jfnk_solve(self, xvec)
 
-    class(JFNKSolver) :: self
-    type(Vector)      :: xvec
+    class(JFNKSolver), intent(inout) :: self ! JFNK instance
+    type(Vector), intent(inout)      :: xvec ! solution vector
 
     call SNESSolve(self % snes_, PETSC_NULL_DOUBLE, xvec % petsc_vec, petsc_err)
 
@@ -230,14 +230,14 @@ contains
 
   subroutine petsc_jfnk_compute_residual(snes_, x, res, ctx, ierr)
 
-    type(snes)     :: snes_
-    type(vec)      :: x
-    type(vec)      :: res
-    integer        :: ierr
-    type(Jfnk_ctx) :: ctx
+    type(snes), intent(inout)     :: snes_ ! PETSc SNES object
+    type(vec), intent(inout)      :: x     ! PETSc solution vector
+    type(vec), intent(inout)      :: res   ! PETSc residual vector
+    integer, intent(inout)        :: ierr  ! error code
+    type(Jfnk_ctx), intent(inout) :: ctx ! JFNK context instance
 
-    type(Vector) :: xvec
-    type(Vector) :: resvec
+    type(Vector) :: xvec   ! solution vector
+    type(Vector) :: resvec ! residual vector
 
     ! We need to use the x and res that come from PETSc because the pointer
     ! location changes and therefore we can not use module variables
@@ -263,15 +263,15 @@ contains
   subroutine petsc_jfnk_compute_jacobian(snes_, x, jac_mf, jac_prec, flag, &
              ctx, ierr)
 
-    type(snes)     :: snes_
-    type(vec)      :: x
-    type(mat)      :: jac_mf
-    type(mat)      :: jac_prec
-    integer        :: flag 
-    type(Jfnk_ctx) :: ctx
-    integer        :: ierr
+    type(snes), intent(inout)     :: snes_     ! PETSc snes instance
+    type(vec), intent(inout)      :: x         ! PETSc solution vector
+    type(mat), intent(inout)      :: jac_mf    ! PETSc matrix free jacobian
+    type(mat), intent(inout)      :: jac_prec  ! PETSc matrix jacobian precond.
+    integer, intent(inout)        :: flag      ! unused madatory flag
+    type(Jfnk_ctx), intent(inout) :: ctx       ! JFNK context instance
+    integer, intent(inout)        :: ierr      ! error code
 
-    type(Vector) :: xvec
+    type(Vector) :: xvec ! solution vector
 
     ! Again, we use the vector that comes from Petsc to build the Jacobian
     ! matrix

src/source.F90

@@ -3,6 +3,7 @@ module source
   use bank_header,     only: Bank
   use constants
   use error,           only: fatal_error
+  use geometry,        only: find_cell
   use geometry_header, only: BASE_UNIVERSE
   use global
   use math,            only: maxwell_spectrum, watt_spectrum
@@ -73,6 +74,9 @@ contains
     real(8) :: p_max(3)   ! maximum coordinates of source
     real(8) :: a          ! Arbitrary parameter 'a'
     real(8) :: b          ! Arbitrary parameter 'b'
+    logical :: found      ! Does the source particle exist within geometry?
+    type(Particle) :: p   ! Temporary particle for using find_cell
+    integer, save :: num_resamples = 0 ! Number of resamples encountered
 
     ! Set weight to one by default
     site % wgt = ONE
@@ -80,11 +84,32 @@ contains
     ! Sample position
     select case (external_source % type_space)
     case (SRC_SPACE_BOX)
-      ! Coordinates sampled uniformly over a box
-      p_min = external_source % params_space(1:3)
-      p_max = external_source % params_space(4:6)
-      r = (/ (prn(), i = 1,3) /)
-      site % xyz = p_min + r*(p_max - p_min)
+      ! Set particle defaults
+      call p % initialize()
+      ! Repeat sampling source location until a good site has been found
+      found = .false.
+      do while (.not.found)
+        ! Coordinates sampled uniformly over a box
+        p_min = external_source % params_space(1:3)
+        p_max = external_source % params_space(4:6)
+        r = (/ (prn(), i = 1,3) /)
+        site % xyz = p_min + r*(p_max - p_min)
+
+        ! Fill p with needed data
+        p % coord0 % xyz = site % xyz
+        p % coord0 % uvw = [ ONE, ZERO, ZERO ]
+
+        ! Now search to see if location exists in geometry
+        call find_cell(p, found)
+        if (.not. found) then
+          num_resamples = num_resamples + 1
+          if (num_resamples == MAX_EXTSRC_RESAMPLES) then
+            message = "Maximum number of external source spatial resamples &
+                      &reached!"
+            call fatal_error()
+          end if
+        end if
+      end do
 
     case (SRC_SPACE_POINT)
       ! Point source
@@ -146,7 +171,7 @@ contains
   end subroutine sample_external_source
 
 !===============================================================================
-! GET_SOURCE_PARTICLE returns the next source particle 
+! GET_SOURCE_PARTICLE returns the next source particle
 !===============================================================================
 
   subroutine get_source_particle(p, index_source)
@@ -155,6 +180,7 @@ contains
     integer(8),     intent(in)    :: index_source
 
     integer(8) :: particle_seed  ! unique index for particle
+    integer :: i
     type(Bank), pointer, save :: src => null()
 !$omp threadprivate(src)
 
@@ -177,6 +203,21 @@ contains
     if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
          p % id == trace_particle) trace = .true.
 
+    ! Set particle track.
+    p % write_track = .false.
+    if (write_all_tracks) then
+      p % write_track = .true.
+    else if (allocated(track_identifiers)) then
+      do i=1, size(track_identifiers(1,:))
+        if (current_batch == track_identifiers(1,i) .and. &
+             &current_gen == track_identifiers(2,i) .and. &
+             &p % id == track_identifiers(3,i)) then
+          p % write_track = .true.
+          exit
+        end if
+      end do
+    end if
+
   end subroutine get_source_particle
 
 !===============================================================================

src/state_point.F90

@@ -464,6 +464,7 @@ contains
     character(19)           :: current_time
     integer                 :: i
     integer                 :: j
+    integer                 :: length(4)
     integer                 :: int_array(3)
     integer, allocatable    :: temp_array(:)
     real(8)                 :: real_array(3) 
@@ -474,7 +475,7 @@ contains
     call write_message(1)
 
     ! Open file for reading
-    call sp % file_open(path_state_point, 'r')
+    call sp % file_open(path_state_point, 'r', serial = .false.)
 
     ! Read filetype
     call sp % read_data(int_array(1), "filetype")
@@ -542,10 +543,9 @@ contains
         call sp % read_data(cmfd % indices, "indicies", length=4, group="cmfd")
         call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, &
              group="cmfd")
+        length = cmfd % indices([4,1,2,3])
         call sp % read_data(cmfd % cmfd_src, "cmfd_src", &
-             length=(/cmfd % indices(4), cmfd % indices(1), &
-             cmfd % indices(2), cmfd % indices(3)/), &
-             group="cmfd")
+             length=length, group="cmfd")
         call sp % read_data(cmfd % entropy, "cmfd_entropy", &
                        length=restart_batch, group="cmfd")
         call sp % read_data(cmfd % balance, "cmfd_balance", &
@@ -705,13 +705,6 @@ contains
 
         end do TALLY_RESULTS
       end if
-
-#ifdef MPI
-      ! If using MPI, file needs to be closed and reopened in parallel
-      ! If serial, we cannot close the file or we will lose our file position
-      call sp % file_close()
-# endif
-
     end if
 
     ! Read source if in eigenvalue mode 
@@ -739,14 +732,6 @@ contains
         ! Open source file 
         call sp % file_open(filename, 'r', serial = .false.)
 
-      else
-
-#ifdef MPI
-      ! Reopen statepoint file in parallel, but only if MPI
-      ! We will compute the position where the source begins
-      call sp % file_open(path_state_point, 'r', serial = .false.)
-#endif
-
       end if
 
       ! Write out source

src/templates/cmfd_t.xml

@@ -1,40 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <options rootname="cmfd" />
-
-  <typedef name="mesh_xml">
-    <component name="dimension"   type="integer-array" />
-    <component name="lower_left"  type="double-array"  />
-    <component name="upper_right" type="double-array"  />
-    <component name="width"       type="double-array"  />
-    <component name="albedo"      type="double-array"  />
-    <component name="map"         type="integer-array" />
-    <component name="energy"      type="double-array"  />
-  </typedef>
-
-  <variable name="mesh_"           tag="mesh"           type="mesh_xml"                                    />
-  <variable name="norm_"           tag="norm"           type="double"        default="1.0"                 />
-  <variable name="feedback_"       tag="feedback"       type="word"          length="5"                    />
-  <variable name="reset_"          tag="reset"          type="word"          length="5"                    />
-  <variable name="downscatter_"    tag="downscatter"    type="word"          length="5"                    />
-  <variable name="solver_"         tag="solver"         type="word"          default="'power'" length="25" />
-  <variable name="snes_monitor_"   tag="snes_monitor"   type="word"          length="5"                    />
-  <variable name="ksp_monitor_"    tag="ksp_monitor"    type="word"          length="5"                    />
-  <variable name="power_monitor_"  tag="power_monitor"  type="word"          length="5"                    />
-  <variable name="write_matrices_" tag="write_matrices" type="word"          length="5"                    />
-  <variable name="run_adjoint_"    tag="run_adjoint"    type="word"          length="5"                    />
-  <variable name="begin_"          tag="begin"          type="integer"       default="1"                   />
-  <variable name="inactive_"       tag="inactive"       type="word"          length="5"                    />
-  <variable name="active_flush_"   tag="active_flush"   type="integer"       default="0"                   />
-  <variable name="inactive_flush_" tag="inactive_flush" type="integer"       default="9999"                />
-  <variable name="num_flushes_"    tag="num_flushes"    type="integer"       default="1"                   />
-  <variable name="display_"        tag="display"        type="word"          default="''" length="25"      />
-  <variable name="spectral_"       tag="spectral"       type="double"        default="0.0"                 />
-  <variable name="shift_"          tag="shift"          type="double"        default="0.0"                 />
-  <variable name="ktol_"           tag="ktol"           type="double"        default="1.e-8"               />
-  <variable name="stol_"           tag="stol"           type="double"        default="1.e-8"               />
-  <variable name="atoli_"          tag="atoli"          type="double"        default="1.e-10"              />
-  <variable name="rtoli_"          tag="rtoli"          type="double"        default="1.e-5"               />
-
-</template>

src/templates/cross_sections_t.xml

@@ -1,26 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <!-- This is the template for reading cross section listings in OpenMC -->
-
-  <options rootname="cross_sections" />
-
-  <typedef name="ace_table_xml">
-    <component name="name" type="word" length="15" />
-    <component name="alias" type="word" length="15" default="''" />
-    <component name="type" type="word" length="10" default="'neutron'" />
-    <component name="zaid" type="integer" default="0" />
-    <component name="metastable" type="integer" default="0" />
-    <component name="awr" type="double" default="0.0" />
-    <component name="temperature" type="double" default="0.0" />
-    <component name="path" type="word" length="255" />
-    <component name="location" type="integer" default="0" />
-  </typedef>
-
-  <variable name="ace_tables_" tag="ace_table" type="ace_table_xml" dimension="1" />
-  <variable name="directory_" tag="directory" type="word" length="255" />
-  <variable name="filetype_" tag="filetype" type="word" length="255" />
-  <variable name="record_length_" tag="record_length" type="integer" />
-  <variable name="entries_" tag="entries" type="integer" />
-
-</template>

src/templates/geometry_t.xml

@@ -1,39 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <!-- This is the template for reading geometry in OpenMC -->
-
-  <options rootname="geometry" />
-
-  <typedef name="cell_xml">
-    <component name="id" type="integer" />
-    <component name="universe" type="integer" default="0" />
-    <component name="material" type="word" length="12" default="''" />
-    <component name="fill" type="integer" default="0" />
-    <component name="surfaces" type="integer-array" />
-    <component name="rotation" type="double-array" />
-    <component name="translation" type="double-array" />
-  </typedef>
-
-  <typedef name="surface_xml">
-    <component name="id" type="integer" />
-    <component name="type" type="word" length="15" />
-    <component name="coeffs" type="double-array" />
-    <component name="boundary" type="word" length="12" default="'transmit'" />
-  </typedef>
-
-  <typedef name="lattice_xml">
-    <component name="id" type="integer" />
-    <component name="type" type="word" length="12" default="'rectangular'" />
-    <component name="dimension" type="integer-array" />
-    <component name="lower_left" type="double-array" />
-    <component name="width" type="double-array" />
-    <component name="universes" type="integer-array" />
-    <component name="outside" type="integer" default="0"/>
-  </typedef>
-
-  <variable name="cell_" tag="cell" type="cell_xml" dimension="1" />
-  <variable name="surface_" tag="surface" type="surface_xml" dimension="1" />
-  <variable name="lattice_" tag="lattice" type="lattice_xml" dimension="1" />
-
-</template>

src/templates/materials_t.xml

@@ -1,44 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <!-- This is the template for reading materials in OpenMC -->
-
-  <options rootname="materials" />
-
-  <!-- Type for specifying density of a material -->
-
-  <typedef name="density_xml">
-    <component name="value" type="double" default="0.0" />
-    <component name="units" type="word" length="10" default="'atom/b-cm'" />
-  </typedef>
-
-  <!-- Type for specifying a nuclide -->
-
-  <typedef name="nuclide_xml">
-    <component name="name" type="word" length="7" default="''" />
-    <component name="xs" type="word" length="3" default="''" />
-    <component name="ao" type="double" default="0.0" />
-    <component name="wo" type="double" default="0.0" />
-  </typedef>
-
-  <!-- Type for specifying S(a,b) data -->
-
-  <typedef name="sab_xml">
-    <component name="name" type="word" length="10" />
-    <component name="xs" type="word" length="3" />
-  </typedef>
-
-  <!-- Type for specifying a material -->
-
-  <typedef name="material_xml">
-    <component name="id" type="integer" />
-    <component name="density" type="density_xml" />
-    <component name="nuclides" tag="nuclide" type="nuclide_xml" dimension="1" />
-    <component name="elements" tag="element" type="nuclide_xml" dimension="1" />
-    <component name="sab" type="sab_xml" dimension="1" />
-  </typedef>
-
-  <variable name="material_" tag="material" type="material_xml" dimension="1" />
-  <variable name="default_xs_" tag="default_xs" type="word" length="3" />
-
-</template>

src/templates/plots_t.xml

@@ -1,32 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <options rootname="plots" />
-
-  <typedef name="col_spec_xml">
-    <component name="id" type="integer" />
-    <component name="rgb" type="integer-array" />
-  </typedef>
-
-  <typedef name="mask_xml">
-    <component name="components" type="integer-array" />
-    <component name="background" type="integer-array" />
-  </typedef>
-
-  <typedef name="plot_xml">
-    <component name="id" type="integer" />
-    <component name="filename" type="word" length="50" default="'plot'" />
-    <component name="type" type="word" length="10" default="'slice'"/>
-    <component name="color" type="word" length="10" default="'cell'"/>
-    <component name="origin" type="double-array" />
-    <component name="width" type="double-array" />
-    <component name="basis" type="word" length="3" default="'xy'" />
-    <component name="pixels" type="integer-array" />
-    <component name="background" type="integer-array"/>
-    <component name="col_spec_" tag="col_spec" type="col_spec_xml" dimension="1" />
-    <component name="mask_" tag="mask" type="mask_xml" dimension="1" />
-  </typedef>
-
-  <variable name="plot_" tag="plot" type="plot_xml" dimension="1" />
-
-</template>

src/templates/settings_t.xml

@@ -1,69 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <options rootname="settings" />
-
-  <typedef name="run_parameters_xml">
-    <component name="batches" type="integer" default="0" />
-    <component name="inactive" type="integer" />
-    <component name="particles" type="word" length="25" default="''" />
-    <component name="generations_per_batch" type="integer" default="1" />
-  </typedef>
-
-  <typedef name="distribution_xml">
-    <component name="type" type="word" length="16" />
-    <component name="length" type="integer" />
-    <component name="interpolation" type="word" length="10" />
-    <component name="parameters" type="double-array" />
-  </typedef>
-
-  <typedef name="source_xml">
-    <component name="file" type="word" length="255" />
-    <component name="space" type="distribution_xml" />
-    <component name="angle" type="distribution_xml"  />
-    <component name="energy" type="distribution_xml" />
-  </typedef>
-
-  <typedef name="cutoff_xml">
-    <component name="weight" type="double" default="0.25" />
-    <component name="weight_avg" type="double" default="1.0" />
-  </typedef>
-
-  <typedef name="mesh_xml">
-    <component name="dimension" type="integer-array" />
-    <component name="lower_left" type="double-array" />
-    <component name="upper_right" type="double-array" />
-  </typedef>
-
-  <typedef name="statepoint_xml">
-    <component name="batches" type="integer-array" />
-    <component name="interval" type="integer" default="0" />
-    <component name="source_separate" type="word" length="5" default="''" />
-    <component name="source_write" type="word" length="5" default="''" />
-  </typedef>
-
-  <variable name="confidence_intervals_" tag="confidence_intervals" type="word" length="5" />
-  <variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
-  <variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
-  <variable name="eigenvalue_" tag="eigenvalue" type="run_parameters_xml" />
-  <variable name="energy_grid_" tag="energy_grid" type="word" length="7" />
-  <variable name="entropy_" tag="entropy" type="mesh_xml" dimension="1" />
-  <variable name="fixed_source_" tag="fixed_source" type="run_parameters_xml" />
-  <variable name="no_reduce_" tag="no_reduce" type="word" length="5" />
-  <variable name="output_" tag="output" type="word-array" length="20" />
-  <variable name="output_path_" tag="output_path" type="word" length="255" />
-  <variable name="ptables_" tag="ptables" type="word" length="5" />
-  <variable name="run_cmfd_" tag="run_cmfd" type="word" length="5" />
-  <variable name="seed_" tag="seed" type="integer" />
-  <variable name="source_" tag="source" type="source_xml" />
-  <variable name="state_point_" tag="state_point" type="statepoint_xml" dimension="1" />
-  <variable name="survival_" tag="survival_biasing" type="word" length="5" />
-  <variable name="threads_" tag="threads" type="integer" />
-  <variable name="trace_" tag="trace" type="integer-array" />
-  <variable name="uniform_fs_" tag="uniform_fs" type="mesh_xml" dimension="1" />
-  <variable name="verbosity_" tag="verbosity" type="integer" />
-
-  <!-- Check for old criticality tag -->
-  <variable name="criticality_" tag="criticality" type="run_parameters_xml" />
-
-</template>

src/templates/tallies_t.xml

@@ -1,37 +0,0 @@
-<?xml version="1.0"?>
-<template>
-
-  <options rootname="tallies" />
-
-  <typedef name="mesh_xml">
-    <component name="id" type="integer" />
-    <component name="type" type="word" length="12" />
-    <component name="dimension" type="integer-array" />
-    <component name="lower_left" type="double-array" />
-    <component name="upper_right" type="double-array" />
-    <component name="width" type="double-array" />
-  </typedef>
-
-  <typedef name="filter_xml">
-    <component name="type" type="word" length="20" default = "''" />
-    <component name="bins" type="word-array" length="20" />
-    <component name="groups" type="word" length="20" default = "''" />
-  </typedef>
-
-  <typedef name="tally_xml">
-    <component name="id" type="integer" />
-    <component name="label" type="word" length="52" default="''" />
-    <component name="estimator" type="word" length="12" default="''" />
-    <component name="filter" type="filter_xml" dimension="1" />
-    <component name="scores" type="word-array" length="20" />
-    <component name="nuclides" type="word-array" length="12" />
-
-    <!-- Check for old tally filter format -->
-    <component name="filters" type="filter_xml" dimension="1" />
-  </typedef>
-
-  <variable name="mesh_" tag="mesh" type="mesh_xml" dimension="1" />
-  <variable name="tally_" tag="tally" type="tally_xml" dimension="1" />
-  <variable name="separate_" tag="assume_separate" type="word" length="3" />
-
-</template>

src/track_output.F90

@@ -0,0 +1,79 @@
+!===============================================================================
+! TRACK_OUTPUT handles output of particle tracks (the paths taken by particles
+! as they are transported through the geometry).
+!===============================================================================
+
+module track_output
+
+  use global
+  use output_interface,  only: BinaryOutput
+  use particle_header,   only: Particle
+  use string,            only: to_str
+
+  implicit none
+
+  integer, private                  :: n_tracks     ! total number of tracks
+  real(8), private, allocatable     :: coords(:,:)  ! track coordinates
+!$omp threadprivate(n_tracks, coords)
+
+contains
+
+!===============================================================================
+! INITIALIZE_PARTICLE_TRACK
+!===============================================================================
+
+  subroutine initialize_particle_track()
+    n_tracks = 0
+  end subroutine initialize_particle_track
+
+!===============================================================================
+! WRITE_PARTICLE_TRACK copies particle position to an array.
+!===============================================================================
+
+  subroutine write_particle_track(p)
+    type(Particle), intent(in)  :: p
+
+    real(8), allocatable :: new_coords(:, :)
+
+    ! Add another column to coords
+    n_tracks = n_tracks + 1
+    if (allocated(coords)) then
+      allocate(new_coords(3, n_tracks))
+      new_coords(:, 1:n_tracks-1) = coords
+      call move_alloc(FROM=new_coords, TO=coords)
+    else
+      allocate(coords(3,1))
+    end if
+
+    ! Write current coordinates into the newest column.
+    coords(:, n_tracks) = p % coord0 % xyz
+  end subroutine write_particle_track
+
+!===============================================================================
+! FINALIZE_PARTICLE_TRACK writes the particle track array to disk.
+!===============================================================================
+
+  subroutine finalize_particle_track(p)
+    type(Particle), intent(in)  :: p
+
+    integer                  :: length(2)
+    character(MAX_FILE_LEN)  :: fname
+    type(BinaryOutput)       :: binout
+
+#ifdef HDF5
+    fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
+         // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
+         // '.h5'
+#else
+    fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
+         // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) &
+         // '.binary'
+#endif
+    call binout % file_create(fname)
+    length = [3, n_tracks]
+    call binout % write_data(coords, 'coordinates', length=length)
+    call binout % file_close()
+    deallocate(coords)
+  end subroutine finalize_particle_track
+
+end module track_output

src/tracking.F90

@@ -13,6 +13,8 @@ module tracking
   use string,          only: to_str
   use tally,           only: score_analog_tally, score_tracklength_tally, &
                              score_surface_current
+  use track_output,    only: initialize_particle_track, write_particle_track, &
+                             finalize_particle_track
 
 contains
 
@@ -67,8 +69,16 @@ contains
     ! Force calculation of cross-sections by setting last energy to zero 
     micro_xs % last_E = ZERO
 
+    ! Prepare to write out particle track.
+    if (p % write_track) then
+      call initialize_particle_track()
+    endif
+
     do while (p % alive)
 
+      ! Write particle track.
+      if (p % write_track) call write_particle_track(p)
+
       if (check_overlaps) call check_cell_overlap(p)
 
       ! Calculate microscopic and macroscopic cross sections -- note: if the
@@ -196,6 +206,12 @@ contains
 
     end do
 
+    ! Finish particle track output.
+    if (p % write_track) then
+      call write_particle_track(p)
+      call finalize_particle_track(p)
+    endif
+
   end subroutine transport
 
 end module tracking

src/utils/build_dependencies.py

@@ -12,10 +12,8 @@ for src in glob.iglob('*.F90'):
                    open(src, 'r').read())
     for name in d:
         if name in ['mpi', 'hdf5', 'h5lt', 'petscsys', 'petscmat', 'petscksp',
-                    'petscsnes', 'petscvec', 'omp_lib']:
+                    'petscsnes', 'petscvec', 'omp_lib', 'fox_dom']:
             continue
-        if name.startswith('xml_data_'):
-            name = name.replace('xml_data_', 'templates/')
         deps.add(name)
     if deps:
         dependencies[module] = sorted(list(deps))

src/utils/convert_xsdir.py

@@ -35,8 +35,12 @@ class Xsdir(object):
         words = line.split()
         if words:
             if words[0].lower().startswith('datapath'):
-                index = line.index('=')
-                self.datapath = line[index+1:].strip()
+                if '=' in words[0]:
+                    index = line.index('=')
+                    self.datapath = line[index+1:].strip()
+                else:
+                    if len(line.strip()) > 8:
+                        self.datapath = line[8:].strip()
             else:
                 self.f.seek(0)
 
@@ -71,7 +75,7 @@ class Xsdir(object):
             # Handle continuation lines
             while words[-1] == '+':
                 extraWords = self.f.readline().split()
-                words = words + extraWords
+                words = words[:-1] + extraWords
             assert len(words) >= 7
 
             # Create XsdirTable object and add to line

src/utils/plot_mesh_tally.py

@@ -1,262 +1,279 @@
 #!/usr/bin/env python2
 
-'''Python script to plot tally data generated by OpenMC.'''
+"""Python script to plot tally data generated by OpenMC."""
 
+import os
 import sys
-from statepoint import *
 
-# Color intensity dependent on individual score?
-
-from PyQt4.QtCore import *
-from PyQt4.QtGui import *
-import matplotlib.pyplot as plt
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg
 from matplotlib.figure import Figure
-from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
-from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar
+import matplotlib.pyplot as plt
 import numpy as np
 
-class AppForm(QMainWindow):
-    def __init__(self, parent=None):
-        QMainWindow.__init__(self, parent)
+from statepoint import *
+
+if sys.version_info[0] < 3:
+    import Tkinter as tk
+else:
+    import tkinter as tk
+import tkFileDialog
+import tkFont
+import tkMessageBox
+import ttk
+
+
+class MeshPlotter(tk.Frame):
+    def __init__(self, parent, filename):
+        tk.Frame.__init__(self, parent)
+
+        self.labels = {'cell': 'Cell:', 'cellborn': 'Cell born:',
+                       'surface': 'Surface:', 'material': 'Material:',
+                       'universe': 'Universe:', 'energyin': 'Energy in:',
+                       'energyout': 'Energy out:'}
+
+        self.filterBoxes = {}
 
         # Read data from source or leakage fraction file
-        self.get_file_data()
-        self.main_frame = QWidget() 
-        self.setCentralWidget(self.main_frame)
-       
-        # Create the Figure, Canvas, and Axes
+        self.get_file_data(filename)
+
+        # Set up top-level window
+        top = self.winfo_toplevel()
+        top.title('Mesh Tally Plotter: ' + filename)
+        top.rowconfigure(0, weight=1)
+        top.columnconfigure(0, weight=1)
+        self.grid(sticky=tk.W+tk.N)
+
+        # Create widgets and draw to screen
+        self.create_widgets()
+        self.update()
+
+    def create_widgets(self):
+        figureFrame = tk.Frame(self)
+        figureFrame.grid(row=0, column=0)
+
+        # Create the Figure and Canvas
         self.dpi = 100
-        self.fig = Figure((5.0, 15.0), dpi=self.dpi)
-        self.canvas = FigureCanvas(self.fig)
-        self.canvas.setParent(self.main_frame)
-        self.axes = self.fig.add_subplot(111)
-        
+        self.fig = Figure((5.0, 5.0), dpi=self.dpi)
+        self.canvas = FigureCanvasTkAgg(self.fig, master=figureFrame)
+        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
+
         # Create the navigation toolbar, tied to the canvas
-        self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame)
+        self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame)
+        self.mpl_toolbar.update()
+        self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1)
 
-        # Grid layout at bottom
-        self.grid = QGridLayout()
+        # Create frame for comboboxes
+        self.selectFrame = tk.Frame(self)
+        self.selectFrame.grid(row=1, column=0, sticky=tk.W+tk.E)
 
-        # Overall layout
-        self.vbox = QVBoxLayout()
-        self.vbox.addWidget(self.canvas)
-        self.vbox.addWidget(self.mpl_toolbar)
-        self.vbox.addLayout(self.grid)
-        self.main_frame.setLayout(self.vbox)
+        # Tally selection
+        labelTally = tk.Label(self.selectFrame, text='Tally:')
+        labelTally.grid(row=0, column=0, sticky=tk.W)
+        self.tallyBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.tallyBox['values'] = [self.datafile.tallies[i].id
+                                   for i in self.meshTallies]
+        self.tallyBox.current(0)
+        self.tallyBox.grid(row=0, column=1, sticky=tk.W+tk.E)
+        self.tallyBox.bind('<<ComboboxSelected>>', self.update)
 
-        # Tally selections
-        label_tally = QLabel("Tally:")
-        self.tally = QComboBox()
-        self.tally.addItems([(str(i + 1)) for i in range(self.n_tallies)])
-        self.connect(self.tally, SIGNAL('activated(int)'),
-                     self._update)
-        self.connect(self.tally, SIGNAL('activated(int)'),
-                     self.populate_boxes)
-        self.connect(self.tally, SIGNAL('activated(int)'),
-                     self.on_draw)
+        # Planar basis selection
+        labelBasis = tk.Label(self.selectFrame, text='Basis:')
+        labelBasis.grid(row=1, column=0, sticky=tk.W)
+        self.basisBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.basisBox['values'] = ('xy', 'yz', 'xz')
+        self.basisBox.current(0)
+        self.basisBox.grid(row=1, column=1, sticky=tk.W+tk.E)
+        self.basisBox.bind('<<ComboboxSelected>>', self.update)
 
-        # Planar basis
-        label_basis = QLabel("Basis:")
-        self.basis = QComboBox()
-        self.basis.addItems(['xy', 'yz', 'xz'])
+        # Axial level
+        labelAxial = tk.Label(self.selectFrame, text='Axial level:')
+        labelAxial.grid(row=2, column=0, sticky=tk.W)
+        self.axialBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.axialBox.grid(row=2, column=1, sticky=tk.W+tk.E)
+        self.axialBox.bind('<<ComboboxSelected>>', self.redraw)
 
-        # Update window when 'Basis' selection is changed
-        self.connect(self.basis, SIGNAL('activated(int)'),
-                     self._update)
-        self.connect(self.basis, SIGNAL('activated(int)'),
-                     self.populate_boxes)
-        self.connect(self.basis, SIGNAL('activated(int)'),
-                     self.on_draw)
+        # Option for mean/uncertainty
+        labelMean = tk.Label(self.selectFrame, text='Mean/Uncertainty:')
+        labelMean.grid(row=3, column=0, sticky=tk.W)
+        self.meanBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.meanBox['values'] = ('Mean', 'Absolute uncertainty',
+                                  'Relative uncertainty')
+        self.meanBox.current(0)
+        self.meanBox.grid(row=3, column=1, sticky=tk.W+tk.E)
+        self.meanBox.bind('<<ComboboxSelected>>', self.update)
 
-        # Axial level within selected basis
-        label_axial_level = QLabel("Axial Level:")
-        self.axial_level = QComboBox()
-        self.connect(self.axial_level, SIGNAL('activated(int)'),
-                     self.on_draw)
-                     
-        # Add Option to plot mean or uncertainty
-        label_mean = QLabel("Mean or Uncertainty:")
-        self.mean = QComboBox()
-        self.mean.addItems(['Mean','Absolute Uncertainty', 
-                            'Relative Uncertainty'])
-        
-        # Update window when mean selection is changed
-        self.connect(self.mean, SIGNAL('activated(int)'),
-                     self.on_draw)
-        
+        # Scores
+        labelScore = tk.Label(self.selectFrame, text='Score:')
+        labelScore.grid(row=4, column=0, sticky=tk.W)
+        self.scoreBox = ttk.Combobox(self.selectFrame, state='readonly')
+        self.scoreBox.grid(row=4, column=1, sticky=tk.W+tk.E)
+        self.scoreBox.bind('<<ComboboxSelected>>', self.redraw)
 
-        self.label_filters = QLabel("Filter options:")
+        # Filter label
+        font = tkFont.Font(weight='bold')
+        labelFilters = tk.Label(self.selectFrame, text='Filters:', font=font)
+        labelFilters.grid(row=5, column=0, sticky=tk.W)
 
-        # Labels for all possible filters
-        self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ',
-                       'surface': 'Surface: ', 'material': 'Material',
-                       'universe': 'Universe: ', 'energyin': 'Energy in: ',
-                       'energyout': 'Energy out: '}
+    def update(self, event=None):
+        if not event:
+            widget = None
+        else:
+            widget = event.widget
 
-        # Empty reusable labels
-        self.qlabels = {}
-        for j in range(8):
-            self.nextLabel = QLabel
-            self.qlabels[j] = self.nextLabel
+        tally_id = self.meshTallies[self.tallyBox.current()]
+        selectedTally = self.datafile.tallies[tally_id]
 
-        # Reusable comboboxes labelled with filter names
-        self.boxes = {}
-        for key in self.labels.keys():
-            self.nextBox = QComboBox()
-            self.connect(self.nextBox, SIGNAL('activated(int)'),
-                         self.on_draw)
-            self.boxes[key] = self.nextBox
+        # Get mesh for selected tally
+        self.mesh = self.datafile.meshes[
+            selectedTally.filters['mesh'].bins[0] - 1]
 
-        # Combobox to select among scores
-        self.score_label = QLabel("Score:")
-        self.scoreBox = QComboBox()
-        for item in self.tally_scores[0]:
-            self.scoreBox.addItems(str(item))
-        self.connect(self.scoreBox, SIGNAL('activated(int)'),
-                     self.on_draw)
+        # Get mesh dimensions
+        self.nx, self.ny, self.nz = self.mesh.dimension
 
-        # Fill layout
-        self.grid.addWidget(label_tally, 0, 0)
-        self.grid.addWidget(self.tally, 0, 1)
-        self.grid.addWidget(label_basis, 1, 0)
-        self.grid.addWidget(self.basis, 1, 1)
-        self.grid.addWidget(label_axial_level, 2, 0)
-        self.grid.addWidget(self.axial_level, 2, 1)
-        self.grid.addWidget(label_mean, 3, 0)
-        self.grid.addWidget(self.mean, 3, 1)
-        self.grid.addWidget(self.label_filters, 4, 0)
+        # Repopulate comboboxes baesd on current basis selection
+        text = self.basisBox['values'][self.basisBox.current()]
+        if text == 'xy':
+            self.axialBox['values'] = [str(i+1) for i in range(self.nz)]
+        elif text == 'yz':
+            self.axialBox['values'] = [str(i+1) for i in range(self.nx)]
+        else:
+            self.axialBox['values'] = [str(i+1) for i in range(self.ny)]
+        self.axialBox.current(0)
 
-        self._update()
-        self.populate_boxes() 
-        self.on_draw()
+        # If update() was called by a change in the basis combobox, we don't
+        # need to repopulate the filters
+        if widget == self.basisBox:
+            self.redraw()
+            return
 
-    def get_file_data(self):
-        # Get data file name from "open file" browser
-        filename = QFileDialog.getOpenFileName(self, 'Select statepoint file', '.')
+        # Update scores
+        self.scoreBox['values'] = selectedTally.scores
+        self.scoreBox.current(0)
 
-        # Create StatePoint object and read in data
-        self.datafile = StatePoint(str(filename))
-        self.datafile.read_results()
-        self.datafile.generate_stdev()
+        # Remove any filter labels/comboboxes that exist
+        for row in range(6, self.selectFrame.grid_size()[1]):
+            for w in self.selectFrame.grid_slaves(row=row):
+                w.grid_forget()
+                w.destroy()
 
-        self.setWindowTitle('Core Map Tool : ' + str(self.datafile.path))
+        # create a label/combobox for each filter in selected tally
+        count = 0
+        for filterType in selectedTally.filters:
+            if filterType == 'mesh':
+                continue
+            count += 1
 
-        # Set maximum colorbar value by maximum tally data value
-        self.maxvalue = self.datafile.tallies[0].results.max()
+            # Create label and combobox for this filter
+            label = tk.Label(self.selectFrame, text=self.labels[filterType])
+            label.grid(row=count+6, column=0, sticky=tk.W)
+            combobox = ttk.Combobox(self.selectFrame, state='readonly')
+            self.filterBoxes[filterType] = combobox
 
-        self.labelList = []
+            # Set combobox items
+            f = selectedTally.filters[filterType]
+            if filterType in ['energyin', 'energyout']:
+                combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2])
+                                      for i in range(f.length)]
+            else:
+                combobox['values'] = [str(i) for i in f.bins]
 
-        # Read mesh dimensions
-#        for mesh in self.datafile.meshes:
-#            self.nx, self.ny, self.nz = mesh.dimension
+            combobox.current(0)
+            combobox.grid(row=count+6, column=1, sticky=tk.W+tk.E)
+            combobox.bind('<<ComboboxSelected>>', self.redraw)
 
-        # Read filter types from statepoint file
-        self.n_tallies = len(self.datafile.tallies)
-        self.tally_list = []
-        for tally in self.datafile.tallies:
-            self.filter_types = []
-            for f in tally.filters:
-                self.filter_types.append(f)
-            self.tally_list.append(self.filter_types)
+        # If There are no filters, leave a 'None available' message
+        if count == 0:
+            count += 1
+            label = tk.Label(self.selectFrame, text="None Available")
+            label.grid(row=count+6, column=0, sticky=tk.W)
 
-        # Read score types from statepoint file
-        self.tally_scores = []
-        for tally in self.datafile.tallies:
-            self.score_types = []
-            for s in tally.scores:
-                self.score_types.append(s)
-            self.tally_scores.append(self.score_types)
-#        print 'self.tally_scores = ', self.tally_scores
+        self.redraw()
 
-    def on_draw(self):
-        """ Redraws the figure
-        """
+    def redraw(self, event=None):
+        basis = self.basisBox.current() + 1
+        axial_level = self.axialBox.current() + 1
+        is_mean = self.meanBox.current()
 
-#        print 'Calling on_draw...'
-        # Get selected basis, axial_level and stage
-        basis = self.basis.currentIndex() + 1
-        axial_level = self.axial_level.currentIndex() + 1
-        is_mean = self.mean.currentIndex()
+        # Get selected tally
+        tally_id = self.meshTallies[self.tallyBox.current()]
+        selectedTally = self.datafile.tallies[tally_id]
 
         # Create spec_list
         spec_list = []
-        for tally in self.datafile.tallies[self.tally.currentIndex()].filters.values():
-            if tally.type == 'mesh':
+        for f in selectedTally.filters.values():
+            if f.type == 'mesh':
                 continue
-            index = self.boxes[tally.type].currentIndex()
-            spec_list.append((tally.type, index))
-        
+            index = self.filterBoxes[f.type].current()
+            spec_list.append((f.type, index))
+
         # Take is_mean and convert it to an index of the score
         score_loc = is_mean
         if score_loc > 1:
             score_loc = 1
-        
-        if self.basis.currentText() == 'xy':
+
+        text = self.basisBox['values'][self.basisBox.current()]
+        if text == 'xy':
             matrix = np.zeros((self.nx, self.ny))
             for i in range(self.nx):
                 for j in range(self.ny):
-                    matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), 
-                        spec_list + [('mesh', (i, j, axial_level))], 
-                        self.scoreBox.currentIndex())[score_loc]
-                    # Calculate relative uncertainty from absolute, if
-                    # requested
+                    matrix[i, j] = self.datafile.get_value(tally_id,
+                        spec_list + [('mesh', (i + 1, j + 1, axial_level))],
+                        self.scoreBox.current())[score_loc]
+                    # Calculate relative uncertainty from absolute, if requested
                     if is_mean == 2:
                         # Take care to handle zero means when normalizing
-                        mean_val = self.datafile.get_value(self.tally.currentIndex(), 
-                            spec_list + [('mesh', (i, j, axial_level))], 
-                            self.scoreBox.currentIndex())[0]
+                        mean_val = self.datafile.get_value(tally_id,
+                            spec_list + [('mesh', (i + 1, j + 1, axial_level))],
+                            self.scoreBox.current())[0]
                         if mean_val > 0.0:
-                            matrix[i,j] = matrix[i,j] / mean_val
+                            matrix[i, j] = matrix[i, j] / mean_val
                         else:
-                            matrix[i,j] = 0.0
-                        
-        elif self.basis.currentText() == 'yz':
+                            matrix[i, j] = 0.0
+
+        elif text == 'yz':
             matrix = np.zeros((self.ny, self.nz))
             for i in range(self.ny):
                 for j in range(self.nz):
-                    matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), 
-                        spec_list + [('mesh', (axial_level, i, j))], 
-                        self.scoreBox.currentIndex())[score_loc]
-                    # Calculate relative uncertainty from absolute, if
-                    # requested
+                    matrix[i, j] = self.datafile.get_value(tally_id,
+                        spec_list + [('mesh', (axial_level, i + 1, j + 1))],
+                        self.scoreBox.current())[score_loc]
+                    # Calculate relative uncertainty from absolute, if requested
                     if is_mean == 2:
                         # Take care to handle zero means when normalizing
-                        mean_val = self.datafile.get_value(self.tally.currentIndex(), 
-                            spec_list + [('mesh', (axial_level, i, j))], 
-                            self.scoreBox.currentIndex())[0]
+                        mean_val = self.datafile.get_value(tally_id,
+                            spec_list + [('mesh', (axial_level, i + 1, j + 1))],
+                            self.scoreBox.current())[0]
                         if mean_val > 0.0:
-                            matrix[i,j] = matrix[i,j] / mean_val
+                            matrix[i, j] = matrix[i, j] / mean_val
                         else:
-                            matrix[i,j] = 0.0
-                   
+                            matrix[i, j] = 0.0
+
         else:
             matrix = np.zeros((self.nx, self.nz))
             for i in range(self.nx):
                 for j in range(self.nz):
-                    matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), 
-                        spec_list + [('mesh', (i, axial_level, j))], 
-                        self.scoreBox.currentIndex())[score_loc]
-                    # Calculate relative uncertainty from absolute, if
-                    # requested
+                    matrix[i, j] = self.datafile.get_value(tally_id,
+                        spec_list + [('mesh', (i + 1, axial_level, j + 1))],
+                        self.scoreBox.current())[score_loc]
+                    # Calculate relative uncertainty from absolute, if requested
                     if is_mean == 2:
                         # Take care to handle zero means when normalizing
-                        mean_val = self.datafile.get_value(self.tally.currentIndex(), 
-                            spec_list + [('mesh', (i, axial_level, j))], 
-                            self.scoreBox.currentIndex())[0]
+                        mean_val = self.datafile.get_value(tally_id,
+                            spec_list + [('mesh', (i + 1, axial_level, j + 1))],
+                            self.scoreBox.current())[0]
                         if mean_val > 0.0:
-                            matrix[i,j] = matrix[i,j] / mean_val
+                            matrix[i, j] = matrix[i, j] / mean_val
                         else:
-                            matrix[i,j] = 0.0
-
-#        print spec_list
+                            matrix[i, j] = 0.0
 
         # Clear the figure
         self.fig.clear()
 
         # Make figure, set up color bar
         self.axes = self.fig.add_subplot(111)
-        cax = self.axes.imshow(matrix, vmin=0.0, vmax=matrix.max(), 
-          interpolation="nearest")
+        cax = self.axes.imshow(matrix.transpose(), vmin=0.0, vmax=matrix.max(),
+                               interpolation='none', origin='lower')
         self.fig.colorbar(cax)
 
         self.axes.set_xticks([])
@@ -266,95 +283,43 @@ class AppForm(QMainWindow):
         # Draw canvas
         self.canvas.draw()
 
-    def _update(self):
-        '''Updates widget to display new relevant comboboxes and figure data
-        '''
-#        print 'Calling _update...'
+    def get_file_data(self, filename):
+        # Create StatePoint object and read in data
+        self.datafile = StatePoint(filename)
+        self.datafile.read_results()
+        self.datafile.generate_stdev()
 
-        self.mesh = self.datafile.meshes[
-          self.datafile.tallies[
-          self.tally.currentIndex()].filters['mesh'].bins[0] - 1]
+        # Find which tallies are mesh tallies
+        self.meshTallies = []
+        for itally, tally in enumerate(self.datafile.tallies):
+            if 'mesh' in tally.filters:
+                self.meshTallies.append(itally)
 
-        self.nx, self.ny, self.nz = self.mesh.dimension
-
-        # Clear axial level combobox
-        self.axial_level.clear()
-
-        # Repopulate axial level combobox based on current basis selection
-        if (self.basis.currentText() == 'xy'):
-            self.axial_level.addItems([str(i+1) for i in range(self.nz)])
-        elif (self.basis.currentText() == 'yz'):
-            self.axial_level.addItems([str(i+1) for i in range(self.nx)])
-        else:
-            self.axial_level.addItems([str(i+1) for i in range(self.ny)])
-
-        # Determine maximum value from current tally data set
-        self.maxvalue = self.datafile.tallies[
-          self.tally.currentIndex()].results.max()
-#        print self.maxvalue
-
-        # Clear and hide old filter labels
-        for item in self.labelList:
-            item.clear()
-
-        # Clear and hide old filter boxes
-        for j in self.labels:
-            self.boxes[j].clear()
-            self.boxes[j].setParent(None)
-
-        self.update()
-
-    def populate_boxes(self):
-#        print 'Calling populate_boxes...'
-
-        n = 5
-        labels = {'cell': 'Cell : ',
-                       'cellborn': 'Cell born: ',
-                       'surface': 'Surface: ',
-                       'material': 'Material: ',
-                       'universe': 'Universe: '}
-
-        # For each filter in newly-selected tally, name a label and fill the
-        # relevant combobox with options
-        for element in self.tally_list[self.tally.currentIndex()]:
-            nextFilter = self.datafile.tallies[
-              self.tally.currentIndex()].filters[element]
-            if element == 'mesh':
-                continue
-
-            label = QLabel(self.labels[element])
-            self.labelList.append(label)
-            combobox = self.boxes[element]
-            self.grid.addWidget(label, n, 0)
-            self.grid.addWidget(combobox, n, 1)
-            n += 1
-
-#            print element
-            if element in ['cell', 'cellborn', 'surface', 'material', 'universe']:
-                combobox.addItems([str(i) for i in nextFilter.bins])
-#                for i in nextFilter.bins:
-#                    print i
-
-            elif element == 'energyin' or element == 'energyout':
-                for i in range(nextFilter.length):
-                    text = (str(nextFilter.bins[i]) + ' to ' + 
-                      str(nextFilter.bins[i+1]))
-                    combobox.addItem(text)
-
-        self.scoreBox.clear()
-        for item in self.tally_scores[self.tally.currentIndex()]:
-            self.scoreBox.addItem(str(item))
-        self.grid.addWidget(self.score_label, n, 0)
-        self.grid.addWidget(self.scoreBox, n, 1)
+        if not self.meshTallies:
+            tkMessageBox.showerror("Invalid StatePoint File",
+                                   "File does not contain mesh tallies!")
+            sys.exit(1)
 
 
+if __name__ == '__main__':
+    # Hide root window
+    root = tk.Tk()
+    root.withdraw()
 
-def main():
-    app = QApplication(sys.argv)
-    form = AppForm()
-    form.show()
-    app.exec_()
+    # If no filename given as command-line argument, open file dialog
+    if len(sys.argv) < 2:
+        filename = tkFileDialog.askopenfilename(title='Select statepoint file',
+                                                initialdir='.')
+    else:
+        filename = sys.argv[1]
 
+    if filename:
+        # Check to make sure file exists
+        if not os.path.isfile(filename):
+            tkMessageBox.showerror("File not found",
+                                   "Could not find regular file: " + filename)
+            sys.exit(1)
 
-if __name__ == "__main__":
-    main()
+        app = MeshPlotter(root, filename)
+        root.deiconify()
+        root.mainloop()

src/utils/statepoint.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python2
 
 import struct
-from math import sqrt
 from collections import OrderedDict
 
 import numpy as np
@@ -10,10 +9,10 @@ import scipy.stats
 filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn',
                 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'}
 
-score_types = {-1: 'flux', 
+score_types = {-1: 'flux',
                 -2: 'total',
                 -3: 'scatter',
-                -4: 'nu-scatter', 
+                -4: 'nu-scatter',
                 -5: 'scatter-n',
                 -6: 'scatter-pn',
                 -7: 'transport',
@@ -276,7 +275,7 @@ class StatePoint(object):
                     f.bins = self._get_int(path=base+'bins')
                 else:
                     f.bins = self._get_int(f.length, path=base+'bins')
-            
+
             base = 'tallies/tally' + str(i+1) + '/'
 
             # Read nuclide bins
@@ -329,7 +328,8 @@ class StatePoint(object):
                 if self._hdf5:
                     path = 'tallies/tally{0}/results'.format(i+1)
                     data = self._f[path].value
-                    t.results = np.dstack((data['sum'], data['sum_sq']))
+                    t.results = np.column_stack((data['sum'], data['sum_sq']))
+                    t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
                 else:
                     t.results = np.array(self._get_double(2*n))
                     t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
@@ -378,7 +378,7 @@ class StatePoint(object):
         Calculates the sample mean and standard deviation of the mean for each
         tally bin.
         """
-        
+
         # Determine number of realizations
         n = self.n_realizations
 
@@ -386,14 +386,14 @@ class StatePoint(object):
         for i in range(len(self.global_tallies)):
             # Get sum and sum of squares
             s, s2 = self.global_tallies[i]
-                    
+
             # Calculate sample mean and replace value
             s /= n
             self.global_tallies[i,0] = s
 
             # Calculate standard deviation
             if s != 0.0:
-                self.global_tallies[i,1] = t_value*sqrt((s2/n - s*s)/(n-1))
+                self.global_tallies[i,1] = t_value*np.sqrt((s2/n - s*s)/(n-1))
 
         # Regular tallies
         for t in self.tallies:
@@ -401,14 +401,14 @@ class StatePoint(object):
                 for j in range(t.results.shape[1]):
                     # Get sum and sum of squares
                     s, s2 = t.results[i,j]
-                    
+
                     # Calculate sample mean and replace value
                     s /= n
                     t.results[i,j,0] = s
 
                     # Calculate standard deviation
                     if s != 0.0:
-                        t.results[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1))
+                        t.results[i,j,1] = t_value*np.sqrt((s2/n - s*s)/(n-1))
 
     def get_value(self, tally_index, spec_list, score_index):
         """Returns a tally score given a list of filters to satisfy.
@@ -457,7 +457,7 @@ class StatePoint(object):
                 filter_index += value*t.filters[f_type].stride
             else:
                 filter_index += f_index*t.filters[f_type].stride
-        
+
         # Return the desired result from Tally.results. This could be the sum and
         # sum of squares, or it could be mean and stdev if self.generate_stdev()
         # has been called already.
@@ -530,7 +530,7 @@ class StatePoint(object):
         for i in range(n_filters):
 
             # compute indices for filter combination
-            filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % 
+            filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) %
                    np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1
 
             # append in dictionary bin with filter
@@ -543,14 +543,14 @@ class StatePoint(object):
                 dims.reverse()
                 dims = np.asarray(dims)
                 if score_str == 'current':
-                    dims += 1 
-                meshmax[1:4] = dims 
+                    dims += 1
+                meshmax[1:4] = dims
                 mesh_bins = np.zeros((n_bins,3))
-                mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % 
+                mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) %
                             np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1
-                mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % 
+                mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) %
                             np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1
-                mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % 
+                mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) %
                             np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1
                 data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],
                             mesh_bins[:,2])})
@@ -575,7 +575,7 @@ class StatePoint(object):
     def _get_data(self, n, typeCode, size):
         return list(struct.unpack('={0}{1}'.format(n,typeCode),
                                   self._f.read(n*size)))
-    
+
     def _get_int(self, n=1, path=None):
         if self._hdf5:
             return [int(v) for v in self._f[path].value]

src/utils/track.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python2
+"""Convert binary particle track to VTK poly data.
+
+Usage information can be obtained by running 'track.py --help':
+
+    usage: track.py [-h] [-o OUT] IN [IN ...]
+
+    Convert particle track file to a .pvtp file.
+
+    positional arguments:
+      IN                 Input particle track data filename(s).
+
+    optional arguments:
+      -h, --help         show this help message and exit
+      -o OUT, --out OUT  Output VTK poly data filename.
+
+"""
+
+import os
+import argparse
+import struct
+import vtk
+
+
+def _parse_args():
+    # Create argument parser.
+    parser = argparse.ArgumentParser(
+        description='Convert particle track file to a .pvtp file.')
+    parser.add_argument('input', metavar='IN', type=str, nargs='+',
+                        help='Input particle track data filename(s).')
+    parser.add_argument('-o', '--out', metavar='OUT', type=str, dest='out',
+                        help='Output VTK poly data filename.')
+
+    # Parse and return commandline arguments.
+    return parser.parse_args()
+
+
+def main():
+    # Parse commandline arguments.
+    args = _parse_args()
+
+    # Check input file extensions.
+    for fname in args.input:
+        if not (fname.endswith('.h5') or fname.endswith('.binary')):
+            raise ValueError("Input file names must either end with '.h5' or"
+                             "'.binary'.")
+    
+    # Make sure that the output filename ends with '.pvtp'.
+    if not args.out:
+        args.out = 'tracks.pvtp'
+    elif os.path.splitext(args.out)[1] != '.pvtp':
+        args.out = ''.join([args.out, '.pvtp'])
+
+    # Import HDF library if HDF files are present
+    for fname in args.input:
+        if fname.endswith('.h5'):
+            import h5py
+            break
+
+    # Initialize data arrays and offset.
+    points = vtk.vtkPoints()
+    cells = vtk.vtkCellArray()
+    point_offset = 0
+    for fname in args.input:
+        # Write coordinate values to points array.
+        if fname.endswith('.binary'):
+            track = open(fname, 'rb').read()
+            coords = [struct.unpack("ddd", track[24*i : 24*(i+1)])
+                      for i in range(len(track)/24)]
+            n_points = len(coords)
+            for triplet in coords:
+                points.InsertNextPoint(triplet)
+        else:
+            coords = h5py.File(fname).get('coordinates')
+            n_points = coords.shape[0]
+            for i in range(n_points):
+                points.InsertNextPoint(coords[i,:])
+                
+        # Create VTK line and assign points to line.
+        line = vtk.vtkPolyLine()
+        line.GetPointIds().SetNumberOfIds(n_points)
+        for i in range(n_points):
+            line.GetPointIds().SetId(i, point_offset+i)
+        
+        cells.InsertNextCell(line)
+        point_offset += n_points
+    data = vtk.vtkPolyData()
+    data.SetPoints(points)
+    data.SetLines(cells)
+    
+    writer = vtk.vtkXMLPPolyDataWriter()
+    writer.SetInput(data)
+    writer.SetFileName(args.out)
+    writer.Write()
+
+
+
+if __name__ == '__main__':
+    main()

src/vector_header.F90

@@ -14,9 +14,9 @@ module vector_header
     real(8), allocatable :: data(:) ! where vector data is stored
     real(8), pointer :: val(:) ! pointer to vector data
 #  ifdef PETSC
-    type(vec) :: petsc_vec
+    type(vec) :: petsc_vec ! PETSc vector
 #  endif
-    logical :: petsc_active
+    logical :: petsc_active ! Logical if PETSc is being used
    contains
      procedure :: create       => vector_create
      procedure :: destroy      => vector_destroy
@@ -26,7 +26,7 @@ module vector_header
      procedure :: copy         => vector_copy
   end type Vector
 
-  integer :: petsc_err
+  integer :: petsc_err ! petsc error code
 
 contains
 
@@ -36,8 +36,8 @@ contains
 
   subroutine vector_create(self, n)
 
-    integer       :: n
-    class(Vector), target :: self
+    integer, intent(in)                    :: n    ! size of vector
+    class(Vector), intent(inout), target   :: self ! vector instance
 
     ! Preallocate vector
     if (.not.allocated(self % data)) allocate(self % data(n))
@@ -60,7 +60,7 @@ contains
 
   subroutine vector_destroy(self)
 
-    class(Vector) :: self
+    class(Vector), intent(inout) :: self ! vector instance
 
 #ifdef PETSC
     if (self % petsc_active) call VecDestroy(self % petsc_vec, petsc_err)
@@ -77,9 +77,9 @@ contains
 
   subroutine vector_add_value(self, idx, val)
 
-    integer       :: idx
-    real(8)       :: val
-    class(Vector) :: self
+    integer, intent(in)          :: idx  ! index location in vector
+    real(8), intent(in)          :: val  ! value to add
+    class(Vector), intent(inout) :: self ! vector instance
 
     self % val(idx) = val
 
@@ -91,7 +91,7 @@ contains
 
   subroutine vector_setup_petsc(self)
 
-    class(Vector) :: self
+    class(Vector), intent(inout) :: self ! vector instance
 
     ! Link to PETSc 
 #ifdef PETSC
@@ -110,11 +110,11 @@ contains
 
   subroutine vector_write_petsc_binary(self, filename)
 
-    character(*) :: filename
-    class(Vector) :: self
+    character(*), intent(in)  :: filename ! name of file to write to
+    class(Vector), intent(in) :: self     ! vector instance
 
 #ifdef PETSC
-    type(PetscViewer) :: viewer
+    type(PetscViewer) :: viewer ! PETSc viewer instance
 
     call PetscViewerBinaryOpen(PETSC_COMM_WORLD, trim(filename), &
          FILE_MODE_WRITE, viewer, petsc_err)

src/xml-fortran/read_from_buffer.inc

@@ -1,79 +0,0 @@
-! Part of XML-Fortran library:
-!
-! $Id: read_from_buffer.inc,v 1.2 2006/03/26 19:05:48 arjenmarkus Exp $
-!
-   character(len=*), intent(in)                  :: buffer
-   integer, intent(inout)                        :: ierror
-
-   integer                                       :: n
-   integer                                       :: i
-   integer                                       :: step
-   integer                                       :: ierr
-   !
-   ! First allocate an array that is surely large enough
-   ! Note:
-   ! This is not completely failsafe: with list-directed
-   ! input you can also use repeat counts (10000*1.0 for
-   ! instance).
-   !
-   allocate( work(len(buffer)/2+1) )
-
-   !
-   ! NOTE:
-   ! This is not portable!!
-   !
-   ! read( buffer, *, iostat = ierror ) (work(n), n=1,size(work))
-   !
-   ! So, use a different strategy: a binary search
-   ! First: establish that we have at least one item to read
-   ! Second: do the binary search
-   !
-!   read( buffer, *, iostat = ierr ) work(1)
-!   if ( ierr /= 0 ) then
-!      n = 0
-!   else
-      n = 1
-      do while ( n <= size(work) )
-         n = 2 * n
-      enddo
-      n    = n / 2
-      step = n / 2
-!      step = n / 2
-
-      do while ( step > 0 )
-         read( buffer, *, iostat = ierr ) (work(i), i = 1,n)
-         if ( ierr /= 0 ) then
-            ierror = ierr       ! Store the error code for later use
-            n = n - step
-         else
-            n = n + step
-         endif
-         step = step / 2
-      enddo
-!   endif
-
-   !
-   ! Then allocate an array of the actual size needed
-   ! and copy the data
-   !
-   !
-   if ( associated( var ) ) then
-      deallocate( var )
-   endif
-   !
-   ! One complication: we may have one too many
-   ! (consequence of the binary search)
-   !
-   read( buffer, *, iostat = ierr ) (work(i), i = 1,n)
-   if ( ierr < 0 ) then
-      n = n - 1
-   endif
-
-   allocate( var(n) )
-   var(1:n) = work(1:n)
-   deallocate( work )
-
-   if ( ierror .lt. 0 ) then
-      ierror = 0
-   endif
-

src/xml-fortran/read_xml_array.inc

@@ -1,54 +0,0 @@
-! Part of XML-Fortran library:
-!
-! $Id: read_xml_array.inc,v 1.3 2007/02/26 20:33:38 arjenmarkus Exp $
-!
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   logical, intent(inout)                        :: has_var
-
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: idx
-   integer                                       :: ierr
-
-   !
-   ! The big trick:
-   ! A string long enough to hold all data strings
-   !
-   character(len=nodata*(len(data(1))+1))        :: bufferd
-   integer                                       :: start
-
-   !
-   ! The value can be stored in an attribute values="..." or in
-   ! the data
-   !
-   has_var = .false.
-   idx = xml_find_attrib( attribs, noattribs, 'values', buffer )
-   if ( idx .gt. 0 ) then
-      call read_from_buffer( buffer, var, ierr )
-      if ( buffer .ne. ' ' ) then
-         has_var = .true.
-      endif
-   else
-      bufferd = ' '
-      start   = 1
-      do idx = 1,nodata
-         if ( data(idx) .ne. ' ' ) then
-            bufferd(start:) = data(idx)
-            start           = start + len(data(idx)) + 1
-         endif
-      enddo
-      call read_from_buffer( bufferd, var, ierr )
-      if ( bufferd .ne. ' ' ) then
-         has_var = .true.
-      endif
-   endif
-
-   if ( ierr .ne. 0 ) then
-      write(*,*) 'Error reading variable - tag = ', trim(tag)
-      has_var = .false.
-   endif

src/xml-fortran/read_xml_primitives.f90

@@ -1,527 +0,0 @@
-! read_xml_prims.f90 - Read routines for primitive data
-!
-! $Id: read_xml_prims.f90,v 1.7 2007/12/07 10:38:41 arjenmarkus Exp $
-!
-! Arjen Markus
-!
-! General information:
-! This module is part of the XML-Fortran library. Its
-! purpose is to help read individual items from an XML
-! file into the variables that have been connected to
-! the various tags. It is used by the code generated
-! by the make_xml_reader program.
-!
-! Because the routines differ mostly by the type of the
-! output variable, the body is included, to prevent
-! too much repeated blocks of code with all the maintenance
-! issues that causes.
-!
-module read_xml_primitives
-   use xmlparse
-   implicit none
-
-   private :: read_from_buffer
-   private :: read_from_buffer_integers
-   private :: read_from_buffer_reals
-   private :: read_from_buffer_doubles
-   private :: read_from_buffer_logicals
-   private :: read_from_buffer_words
-
-   interface read_from_buffer
-      module procedure read_from_buffer_integers
-      module procedure read_from_buffer_reals
-      module procedure read_from_buffer_doubles
-      module procedure read_from_buffer_logicals
-      module procedure read_from_buffer_words
-   end interface
-
-contains
-
-! skip_until_endtag --
-!    Routine to read the XML file until the end tag is encountered
-!
-! Arguments:
-!    info        The XML file data structure
-!    tag         The tag in question
-!    attribs     Array of attributes and their values
-!    data        Array of strings, representing the data
-!    error       Has an error occurred?
-!
-subroutine skip_until_endtag( info, tag, attribs, data, error )
-   type(XML_PARSE), intent(inout)                  :: info
-   character(len=*), intent(in)                    :: tag
-   character(len=*), dimension(:,:), intent(inout) :: attribs
-   character(len=*), dimension(:), intent(inout)   :: data
-   logical, intent(out)                            :: error
-
-   integer                                         :: noattribs
-   integer                                         :: nodata
-   integer                                         :: ierr
-   logical                                         :: endtag
-   character(len=len(tag))                         :: newtag
-
-   error = .true.
-   do
-      call xml_get( info, newtag, endtag, attribs, noattribs, &
-                    data, nodata )
-      if ( xml_error(info) ) then
-         error = .true.
-         exit
-      endif
-      if ( endtag .and. newtag == tag ) then
-         exit
-      endif
-   enddo
-end subroutine skip_until_endtag
-
-! read_xml_integer --
-!    Routine to read a single integer from the parsed data
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question (error message only)
-!    endtag      End tag found? (Dummy argument, actually)
-!    attribs     Array of attributes and their values
-!    noattribs   Number of attributes found
-!    data        Array of strings, representing the data
-!    nodata      Number of data strings
-!    var         Variable to be filled
-!    has_var     Has the variable been set?
-!
-subroutine read_xml_integer( info, tag, endtag, attribs, noattribs, data, nodata, &
-                             var, has_var )
-   integer, intent(inout)                       :: var
-
-   include 'read_xml_scalar.inc'
-
-end subroutine read_xml_integer
-
-! read_xml_line --
-!    Routine to read a single line of text from the parsed data
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question (error message only)
-!    endtag      End tag found? (Dummy argument, actually)
-!    attribs     Array of attributes and their values
-!    noattribs   Number of attributes found
-!    data        Array of strings, representing the data
-!    nodata      Number of data strings
-!    var         Variable to be filled
-!    has_var     Has the variable been set?
-!
-subroutine read_xml_line( info, tag, endtag, attribs, noattribs, data, nodata, &
-                          var, has_var )
-   type(XML_PARSE), intent(inout)               :: info
-   character(len=*), intent(in)                 :: tag
-   logical, intent(inout)                       :: endtag
-   character(len=*), dimension(:,:), intent(in) :: attribs
-   integer, intent(in)                          :: noattribs
-   character(len=*), dimension(:), intent(in)   :: data
-   integer, intent(in)                          :: nodata
-   character(len=*), intent(inout)              :: var
-   logical, intent(inout)                       :: has_var
-
-   character(len=len(attribs(1,1)))             :: buffer
-   integer                                      :: idx
-   integer                                      :: ierr
-
-   !
-   ! The value can be stored in an attribute value="..." or in
-   ! the data
-   !
-   has_var = .false.
-   idx = xml_find_attrib( attribs, noattribs, 'value', buffer )
-   if ( idx > 0 ) then
-      var     = buffer
-      has_var = .true.
-   else
-      do idx = 1,nodata
-         if ( data(idx) /= ' ' ) then
-            var = data(idx)
-            has_var = .true.
-            exit
-         endif
-      enddo
-   endif
-end subroutine read_xml_line
-
-! read_xml_real, ... --
-!    See read_xml_integer for an explanation
-!
-subroutine read_xml_real( info, tag, endtag, attribs, noattribs, data, nodata, &
-                          var, has_var )
-   real, intent(inout)                          :: var
-
-   include 'read_xml_scalar.inc'
-
-end subroutine read_xml_real
-
-subroutine read_xml_double( info, tag, endtag, attribs, noattribs, data, nodata, &
-                            var, has_var )
-   real(kind=kind(1.0d00)), intent(inout)       :: var
-
-   include 'read_xml_scalar.inc'
-
-end subroutine read_xml_double
-
-subroutine read_xml_logical( info, tag, endtag, attribs, noattribs, data, nodata, &
-                             var, has_var )
-   logical, intent(inout)       :: var
-
-   include 'read_xml_scalar.inc'
-
-end subroutine read_xml_logical
-
-subroutine read_xml_word( info, tag, endtag, attribs, noattribs, data, nodata, &
-                          var, has_var )
-   character(len=*), intent(inout)       :: var
-
-   include 'read_xml_word.inc'
-
-end subroutine read_xml_word
-
-! read_xml_integer_array --
-!    Routine to read a one-dimensional integer array from the parsed
-!    ata
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question (error message only)
-!    endtag      End tag found? (Dummy argument, actually)
-!    attribs     Array of attributes and their values
-!    noattribs   Number of attributes found
-!    data        Array of strings, representing the data
-!    nodata      Number of data strings
-!    var         Variable to be filled
-!    has_var     Has the variable been set?
-!
-subroutine read_xml_integer_array( info, tag, endtag, attribs, noattribs, data, &
-                                   nodata, var, has_var )
-   integer, dimension(:), pointer                :: var
-
-   include 'read_xml_array.inc'
-
-end subroutine read_xml_integer_array
-
-! read_xml_line_array --
-!    Routine to read an array of lines of text from the parsed data
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question (error message only)
-!    attribs     Array of attributes and their values
-!    noattribs   Number of attributes found
-!    data        Array of strings, representing the data
-!    nodata      Number of data strings
-!    var         Variable to be filled
-!    has_var     Has the variable been set?
-!
-subroutine read_xml_line_array( info, tag, endtag, attribs, noattribs, data, &
-                                nodata, var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   character(len=*), dimension(:), pointer       :: var
-   logical, intent(inout)                        :: has_var
-
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: idx
-   integer                                       :: idxv
-   integer                                       :: ierr
-   logical                                       :: started
-
-   !
-   ! The value can be stored in an attribute values="..." or in
-   ! the data
-   !
-   has_var = .false.
-   idx = xml_find_attrib( attribs, noattribs, 'values', buffer )
-   if ( idx > 0 ) then
-      allocate( var(1:1) )
-      var(1) = buffer
-      if ( buffer /= ' ' ) then
-         has_var = .true.
-      endif
-   else
-      idxv    = 0
-      started = .false.
-      do idx = 1,nodata
-         if ( data(idx) /= ' ' .or. started ) then
-            if ( .not. started ) then
-               allocate( var(1:nodata-idx+1) )
-               started = .true.
-            endif
-            idxv = idxv + 1
-            var(idxv) = data(idx)
-         endif
-      enddo
-      if ( started ) then
-         has_var = .true.
-      endif
-   endif
-end subroutine read_xml_line_array
-
-! read_xml_real_array, ... --
-!    See read_xml_integer_array for an explanation
-!
-subroutine read_xml_real_array( info, tag, endtag, attribs, noattribs, data, &
-                                nodata, var, has_var )
-   real, dimension(:), pointer :: var
-
-   include 'read_xml_array.inc'
-
-end subroutine read_xml_real_array
-
-subroutine read_xml_double_array( info, tag, endtag, attribs, noattribs, data, &
-                                  nodata, var, has_var )
-   real(kind=kind(1.0d00)), dimension(:), pointer :: var
-
-   include 'read_xml_array.inc'
-
-end subroutine read_xml_double_array
-
-subroutine read_xml_logical_array( info, tag, endtag, attribs, noattribs, data, &
-                                   nodata, var, has_var )
-   logical, dimension(:), pointer :: var
-
-   include 'read_xml_array.inc'
-
-end subroutine read_xml_logical_array
-
-subroutine read_xml_word_array( info, tag, endtag, attribs, noattribs, data, &
-                                nodata, var, has_var )
-   character(len=*), dimension(:), pointer :: var
-
-   include 'read_xml_array.inc'
-
-end subroutine read_xml_word_array
-
-! read_from_buffer_integers --
-!    Routine to read all integers from a long string
-!
-! Arguments:
-!    buffer      String containing the data
-!    var         Variable to be filled
-!    ierror      Error flag
-!
-subroutine read_from_buffer_integers( buffer, var, ierror )
-   integer, dimension(:), pointer                :: var
-   integer, dimension(:), pointer                :: work
-
-   include 'read_from_buffer.inc'
-
-end subroutine read_from_buffer_integers
-
-! read_xml_from_buffer_reals, ... -
-!    See read_xml_from_buffer_integers for an explanation
-!
-subroutine read_from_buffer_reals( buffer, var, ierror )
-   real, dimension(:), pointer                :: var
-   real, dimension(:), pointer                :: work
-
-   include 'read_from_buffer.inc'
-
-end subroutine read_from_buffer_reals
-
-subroutine read_from_buffer_doubles( buffer, var, ierror )
-   real(kind=kind(1.0d00)), dimension(:), pointer :: var
-   real(kind=kind(1.0d00)), dimension(:), pointer :: work
-
-   include 'read_from_buffer.inc'
-
-end subroutine read_from_buffer_doubles
-
-subroutine read_from_buffer_logicals( buffer, var, ierror )
-   logical, dimension(:), pointer :: var
-   logical, dimension(:), pointer :: work
-
-   include 'read_from_buffer.inc'
-
-end subroutine read_from_buffer_logicals
-
-subroutine read_from_buffer_words( buffer, var, ierror )
-   character(len=*), dimension(:), pointer :: var
-   character(len=len(var)), dimension(:), pointer :: work
-
-   include 'read_from_buffer.inc'
-
-end subroutine read_from_buffer_words
-
-! read_xml_word_1dim, ... -
-!    Read an array of "words" (or ...) but from different elements
-!
-subroutine read_xml_integer_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                                  var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   integer, dimension(:), pointer                :: var
-   logical, intent(inout)                        :: has_var
-
-   integer,dimension(:), pointer                 :: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_integer( info, tag, endtag, attribs, noattribs, data, nodata, &
-                          var(newsize), has_var )
-
-end subroutine read_xml_integer_1dim
-
-subroutine read_xml_real_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                               var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   real, dimension(:), pointer                   :: var
-   logical, intent(inout)                        :: has_var
-
-   real, dimension(:), pointer                   :: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_real( info, tag, endtag, attribs, noattribs, data, nodata, &
-                       var(newsize), has_var )
-
-end subroutine read_xml_real_1dim
-
-subroutine read_xml_double_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                                 var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   real(kind=kind(1.0d00)), dimension(:), pointer:: var
-   logical, intent(inout)                        :: has_var
-
-   real(kind=kind(1.0d00)), dimension(:), pointer:: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_double( info, tag, endtag, attribs, noattribs, data, nodata, &
-                       var(newsize), has_var )
-
-end subroutine read_xml_double_1dim
-
-subroutine read_xml_logical_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                                  var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   logical, dimension(:), pointer                :: var
-   logical, intent(inout)                        :: has_var
-
-   logical, dimension(:), pointer                :: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_logical( info, tag, endtag, attribs, noattribs, data, nodata, &
-                          var(newsize), has_var )
-
-end subroutine read_xml_logical_1dim
-
-subroutine read_xml_word_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                             var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   character(len=*), dimension(:), pointer       :: var
-   logical, intent(inout)                        :: has_var
-
-   character(len=len(var)),dimension(:), pointer :: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_word( info, tag, endtag, attribs, noattribs, data, nodata, &
-                       var(newsize), has_var )
-
-end subroutine read_xml_word_1dim
-
-subroutine read_xml_line_1dim( info, tag, endtag, attribs, noattribs, data, nodata, &
-                             var, has_var )
-   type(XML_PARSE), intent(inout)                :: info
-   character(len=*), intent(in)                  :: tag
-   logical, intent(inout)                        :: endtag
-   character(len=*), dimension(:,:), intent(in)  :: attribs
-   integer, intent(in)                           :: noattribs
-   character(len=*), dimension(:), intent(in)    :: data
-   integer, intent(in)                           :: nodata
-   character(len=*), dimension(:), pointer       :: var
-   logical, intent(inout)                        :: has_var
-
-   character(len=len(var)),dimension(:), pointer :: newvar
-   character(len=len(attribs(1,1)))              :: buffer
-   integer                                       :: newsize
-   integer                                       :: ierr
-
-   newsize = size(var) + 1
-   allocate( newvar(1:newsize) )
-   newvar(1:newsize-1) = var
-   deallocate( var )
-   var => newvar
-
-   call read_xml_line( info, tag, endtag, attribs, noattribs, data, nodata, &
-                       var(newsize), has_var )
-
-end subroutine read_xml_line_1dim
-
-
-end module read_xml_primitives

src/xml-fortran/read_xml_scalar.inc

@@ -1,40 +0,0 @@
-! Part of XML-Fortran library:
-!
-! $Id: read_xml_scalar.inc,v 1.3 2007/02/26 20:33:38 arjenmarkus Exp $
-!
-   type(XML_PARSE), intent(inout)               :: info
-   character(len=*), intent(in)                 :: tag
-   logical, intent(inout)                       :: endtag
-   character(len=*), dimension(:,:), intent(in) :: attribs
-   integer, intent(in)                          :: noattribs
-   character(len=*), dimension(:), intent(in)   :: data
-   integer, intent(in)                          :: nodata
-   logical, intent(inout)                       :: has_var
-
-   character(len=len(attribs(1,1)))             :: buffer
-   integer                                      :: idx
-   integer                                      :: ierr
-
-   !
-   ! The value can be stored in an attribute value="..." or in
-   ! the data
-   !
-   has_var = .false.
-   idx = xml_find_attrib( attribs, noattribs, 'value', buffer )
-   if ( idx .gt. 0 ) then
-      read( buffer, *, iostat=ierr ) var
-      has_var = .true.
-   else
-      do idx = 1,nodata
-         if ( data(idx) .ne. ' ' ) then
-            read( data(idx), *, iostat=ierr ) var
-            has_var = .true.
-            exit
-         endif
-      enddo
-   endif
-
-   if ( ierr .ne. 0 ) then
-      write(*,*) 'Error reading variable - tag = ', trim(tag)
-      has_var = .false.
-   endif

src/xml-fortran/read_xml_word.inc

@@ -1,38 +0,0 @@
-! Part of XML-Fortran library:
-!
-   type(XML_PARSE), intent(inout)               :: info
-   character(len=*), intent(in)                 :: tag
-   logical, intent(inout)                       :: endtag
-   character(len=*), dimension(:,:), intent(in) :: attribs
-   integer, intent(in)                          :: noattribs
-   character(len=*), dimension(:), intent(in)   :: data
-   integer, intent(in)                          :: nodata
-   logical, intent(inout)                       :: has_var
-
-   character(len=len(attribs(1,1)))             :: buffer
-   integer                                      :: idx
-   integer                                      :: ierr
-
-   !
-   ! The value can be stored in an attribute value="..." or in
-   ! the data
-   !
-   has_var = .false.
-   idx = xml_find_attrib( attribs, noattribs, 'value', buffer )
-   if ( idx .gt. 0 ) then
-      read( buffer, *, iostat=ierr ) var
-      has_var = .true.
-   else
-      do idx = 1,nodata
-         if ( data(idx) .ne. ' ' ) then
-            read( data(idx), '(A)', iostat=ierr ) var
-            has_var = .true.
-            exit
-         endif
-      enddo
-   endif
-
-   if ( ierr .ne. 0 ) then
-      write(*,*) 'Error reading variable - tag = ', trim(tag)
-      has_var = .false.
-   endif

src/xml-fortran/write_xml_primitives.f90

@@ -1,489 +0,0 @@
-! write_xml_prims.f90 - Write routines for primitive data
-!
-! $Id: write_xml_prims.f90,v 1.2 2007/12/27 05:13:59 arjenmarkus Exp $
-!
-! Arjen Markus
-!
-! General information:
-! This module is part of the XML-Fortran library. Its
-! purpose is to write individual items to an XML
-! file using the right tag. It is used by the code generated
-! by the make_xml_reader program.
-!
-module write_xml_primitives
-    use xmlparse
-    implicit none
-
-!   interface write_to_xml
-!       module procedure write_to_xml_integers
-!       module procedure write_to_xml_reals
-!       module procedure write_to_xml_doubles
-!       module procedure write_to_xml_logicals
-!       module procedure write_to_xml_words
-!   end interface
-    interface write_to_xml_word
-       module procedure write_to_xml_string
-    end interface
-    interface write_to_xml_line
-       module procedure write_to_xml_string
-    end interface
-
-contains
-
-! write_to_xml_integer --
-!    Routine to write a single integer to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_integer( info, tag, indent, value )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    integer, intent(in)                          :: value
-
-    character(len=100)                           :: indentation
-
-    indentation = ' '
-    write( info%lun, '(4a,i0,3a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>', value, '</', trim(tag), '>'
-
-end subroutine write_to_xml_integer
-
-! write_to_xml_integer_1dim --
-!    Routine to write an array of integers to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    values      Values to be written
-!
-subroutine write_to_xml_integer_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    integer, dimension(:), intent(in)            :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_integer( info, tag, indent, values(i) )
-    enddo
-
-end subroutine write_to_xml_integer_1dim
-
-! write_to_xml_real --
-!    Routine to write a single real value (single precision) to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_real( info, tag, indent, value )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real, intent(in)                             :: value
-
-    character(len=100)                           :: indentation
-    character(len=12)                            :: buffer
-
-    indentation = ' '
-    write( buffer, '(1pg12.4)' ) value
-    write( info%lun, '(8a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>', trim(adjustl(buffer)), '</', trim(tag), '>'
-
-end subroutine write_to_xml_real
-
-! write_to_xml_real_1dim --
-!    Routine to write an array of reals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    values      Values to be written
-!
-subroutine write_to_xml_real_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real, dimension(:), intent(in)               :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_real( info, tag, indent, values(i) )
-    enddo
-
-end subroutine write_to_xml_real_1dim
-
-! write_to_xml_double --
-!    Routine to write one real value (double precision) to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_double( info, tag, indent, value )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real(kind=kind(1.0d0)), intent(in)           :: value
-
-    character(len=100)                           :: indentation
-    character(len=16)                            :: buffer
-
-    indentation = ' '
-    write( buffer, '(1pg16.7)' ) value
-    write( info%lun, '(8a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>', trim(adjustl(buffer)), '</', trim(tag), '>'
-
-end subroutine write_to_xml_double
-
-! write_to_xml_double_1dim --
-!    Routine to write an array of double precision reals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    values      Values to be written
-!
-subroutine write_to_xml_double_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real(kind=kind(1.0d00)), dimension(:), intent(in) :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_double( info, tag, indent, values(i) )
-    enddo
-
-end subroutine write_to_xml_double_1dim
-
-! write_to_xml_string --
-!    Routine to write one string to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_string( info, tag, indent, value )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    character(len=*), intent(in)                 :: value
-
-    character(len=100)                           :: indentation
-
-    !
-    ! NOTE: No guards against <, >, & and " yet!
-    ! NOTE: difference needed between words and lines?
-    !
-    indentation = ' '
-    write( info%lun, '(8a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>', trim(value), '</', trim(tag), '>'
-
-end subroutine write_to_xml_string
-
-! write_to_xml_word_1dim --
-!    Routine to write an array of single words to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_word_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    character(len=*), dimension(:), intent(in)   :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_string( info, tag, indent, values(i) )
-    enddo
-end subroutine write_to_xml_word_1dim
-
-! write_to_xml_string_1dim --
-!    Routine to write an array of strings to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    values      Values to be written
-!
-subroutine write_to_xml_string_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    character(len=*), dimension(:), intent(in)   :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_string( info, tag, indent, values(i) )
-    enddo
-
-end subroutine write_to_xml_string_1dim
-
-! write_to_xml_logical --
-!    Routine to write one logical to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    value       Value to be written
-!
-subroutine write_to_xml_logical( info, tag, indent, value )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    logical, intent(in)                          :: value
-
-    character(len=100)                           :: indentation
-
-    indentation = ' '
-    if ( value ) then
-        write( info%lun, '(8a)' ) indentation(1:min(indent,100)), &
-            '<', trim(tag), '>true</', trim(tag), '>'
-    else
-        write( info%lun, '(8a)' ) indentation(1:min(indent,100)), &
-            '<', trim(tag), '>false</', trim(tag), '>'
-    endif
-
-end subroutine write_to_xml_logical
-
-! write_to_xml_logical_1dim --
-!    Routine to write an array of logicals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    values      Values to be written
-!
-subroutine write_to_xml_logical_1dim( info, tag, indent, values )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    logical, dimension(:), intent(in)            :: values
-
-    integer                                      :: i
-
-    do i = 1,size(values)
-        call write_to_xml_logical( info, tag, indent, values(i) )
-    enddo
-
-end subroutine write_to_xml_logical_1dim
-
-! write_to_xml_integer_array --
-!    Routine to write an array of integers to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_integer_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    integer, dimension(:), intent(in)            :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i, i2, j
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array),10
-        i2 = min( i + 9, size(array) )
-        write( info%lun, '(a,10i12)' ) indentation(1:min(indent+4,100)), &
-            ( array(j) ,j = i,i2 )
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_integer_array
-
-! write_to_xml_real_array --
-!    Routine to write an array of single precision reals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_real_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real, dimension(:), intent(in)               :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i, i2, j
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array),10
-        i2 = min( i + 9, size(array) )
-        write( info%lun, '(a,10g12.4)' ) indentation(1:min(indent+4,100)), &
-            ( array(j) ,j = i,i2 )
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_real_array
-
-! write_to_xml_double_array --
-!    Routine to write an array of double precision reals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_double_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    real(kind=kind(1.0d0)), dimension(:), intent(in) :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i, i2, j
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array),5
-        i2 = min( i + 4, size(array) )
-        write( info%lun, '(a,5g20.7)' ) indentation(1:min(indent+4,100)), &
-            ( array(j) ,j = i,i2 )
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_double_array
-
-! write_to_xml_logical_array --
-!    Routine to write an array of logicals to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_logical_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    logical, dimension(:), intent(in)            :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i, i2, j
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array),10
-        i2 = min( i + 9, size(array) )
-        write( info%lun, '(a,10a)' ) indentation(1:min(indent+4,100)), &
-            ( merge('true  ', 'false ', array(j)) ,j = i,i2 )
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_logical_array
-
-! write_to_xml_word_array --
-!    Routine to write an array of words to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_word_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    character(len=*), dimension(:), intent(in)   :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i, i2, j
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array),10
-        i2 = min( i + 9, size(array) )
-        write( info%lun, '(a,20a)' ) indentation(1:min(indent+4,100)), &
-            ( trim(array(j)) , ' ' ,j = i,i2 )
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_word_array
-
-! write_to_xml_line_array --
-!    Routine to write an array of lines to the XML file
-!
-! Arguments:
-!    info        XML parser structure
-!    tag         The tag in question
-!    indent      Number of spaces for indentation
-!    array       Values to be written
-!
-subroutine write_to_xml_line_array( info, tag, indent, array )
-    type(XML_PARSE), intent(in)                  :: info
-    character(len=*), intent(in)                 :: tag
-    integer, intent(in)                          :: indent
-    logical, dimension(:), intent(in)            :: array
-
-    character(len=100)                           :: indentation
-    integer                                      :: i
-
-    indentation = ' '
-
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '<', trim(tag), '>'
-    do i = 1,size(array)
-        write( info%lun, '(a)' ) indentation(1:min(indent+4,100)), &
-            array(i)
-    enddo
-    write( info%lun, '(4a)' ) indentation(1:min(indent,100)), &
-        '</', trim(tag), '>'
-
-end subroutine write_to_xml_line_array
-
-end module write_xml_primitives

src/xml-fortran/xml-fortran.LICENSE

@@ -1,37 +0,0 @@
-Note: The official xml-fortran repository on sourceforge
-<http://xml-fortran.sourceforge.net> only lists the license as "BSD License". It
-is assumed that this refers to the 3-clause BSD license ("modified" or "new")
-that is shown below.
-
-The license for xml-fortran covers all source code in the xml-fortran/
-directory.
-
-================================================================================
-
-Copyright (c) 2008 Arjen Markus
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-    * Redistributions of source code must retain the above copyright notice,
-      this list of conditions and the following disclaimer.
-
-    * Redistributions in binary form must reproduce the above copyright notice,
-      this list of conditions and the following disclaimer in the documentation
-      and/or other materials provided with the distribution.
-
-    * Neither the name of Arjen Markus nor the names of its contributors may be
-      used to endorse or promote products derived from this software without
-      specific prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
-ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-DISCLAIMED. IN NO EVENT SHALL ARJEN MARKUS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
-OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
-ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.