diff --git a/src/ace.F90 b/src/ace.F90
index e4bc3a2f90..82aa705baa 100644
--- a/src/ace.F90
+++ b/src/ace.F90
@@ -19,7 +19,6 @@ module ace
   integer                   :: NXS(16)   ! Descriptors for ACE XSS tables
   real(8), allocatable      :: XSS(:)    ! Cross section data
   integer                   :: XSS_index ! Current index in XSS data
-  character(2*MAX_LINE_LEN) :: message   ! Message to output unit
 
   private :: JXS
   private :: NXS
@@ -152,9 +151,9 @@ contains
 
         ! Check to make sure S(a,b) table matched a nuclide
         if (mat % i_sab_nuclides(k) == NONE) then
-          message = "S(a,b) table " // trim(mat % sab_names(k)) // " did not &
-               &match any nuclide on material " // trim(to_str(mat % id))
-          call fatal_error(message)
+          call fatal_error("S(a,b) table " // trim(mat % sab_names(k)) &
+               &// " did not match any nuclide on material " &
+               &// trim(to_str(mat % id)))
         end if
       end do ASSIGN_SAB
 
@@ -265,17 +264,14 @@ contains
     ! Check if ACE library exists and is readable
     inquire(FILE=filename, EXIST=file_exists, READ=readable)
     if (.not. file_exists) then
-      message = "ACE library '" // trim(filename) // "' does not exist!"
-      call fatal_error(message)
+      call fatal_error("ACE library '" // trim(filename) // "' does not exist!")
     elseif (readable(1:3) == 'NO') then
-      message = "ACE library '" // trim(filename) // "' is not readable! &
-           &Change file permissions with chmod command."
-      call fatal_error(message)
+      call fatal_error("ACE library '" // trim(filename) // "' is not readable!&
+           &Change file permissions with chmod command.")
     end if
 
     ! display message
-    message = "Loading ACE cross section table: " // listing % name
-    call write_message(message, 6)
+    call write_message("Loading ACE cross section table: " // listing % name, 6)
 
     if (filetype == ASCII) then
       ! =======================================================================
@@ -294,9 +290,8 @@ contains
       ! 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(message)
+        call fatal_error("XS listing entry " // trim(listing % name) // " did &
+             &not match ACE data, " // trim(name) // " found instead.")
       end if
 
       ! Read more header and NXS and JXS
@@ -378,8 +373,8 @@ contains
       ! if any fissionable material is found in a fixed source calculation,
       ! abort the run.
       if (run_mode == MODE_FIXEDSOURCE .and. nuc % fissionable) then
-        message = "Cannot have fissionable material in a fixed source run."
-        call fatal_error(message)
+        call fatal_error("Cannot have fissionable material in a fixed source &
+             &run.")
       end if
 
       ! for fissionable nuclides, precalculate microscopic nu-fission cross
@@ -1316,9 +1311,8 @@ contains
 
       ! Abort if no corresponding inelastic reaction was found
       if (nuc % urr_inelastic == NONE) then
-        message = "Could not find inelastic reaction specified on " &
-             // "unresolved resonance probability table."
-        call fatal_error(message)
+        call fatal_error("Could not find inelastic reaction specified on &
+             &unresolved resonance probability table.")
       end if
     end if
 
@@ -1344,9 +1338,8 @@ contains
 
     ! Check for negative values
     if (any(nuc % urr_data % prob < ZERO)) then
-      message = "Negative value(s) found on probability table for nuclide " &
-           // nuc % name
-      if (master) call warning(message)
+      if (master) call warning("Negative value(s) found on probability table &
+           &for nuclide " // nuc % name)
     end if
 
   end subroutine read_unr_res
diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index c49fc7808c..b392e86a28 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -11,8 +11,6 @@ module cmfd_data
   private
   public :: set_up_cmfd, neutron_balance
 
-  character(2*MAX_LINE_LEN) :: message
-
 contains
 
 !==============================================================================
@@ -162,10 +160,9 @@ contains
 
                 ! Detect zero flux, abort if located
                 if ((flux - ZERO) < TINY_BIT) then
-                  message = 'Detected zero flux without coremap overlay at: (' &
-                          // to_str(i) // ',' // to_str(j) // ',' // to_str(k) &
-                          // ') in group ' // to_str(h)
-                  call fatal_error(message)
+                  call fatal_error('Detected zero flux without coremap overlay &
+                       &at: (' // to_str(i) // ',' // to_str(j) // ',' // &
+                       &to_str(k) // ') in group ' // to_str(h))
                 end if
 
                 ! Get total rr and convert to total xs
@@ -767,8 +764,7 @@ contains
 
     ! write that dhats are zero
     if (dhat_reset) then
-      message = 'Dhats reset to zero.'
-      call write_message(message, 1)
+      call write_message('Dhats reset to zero.', 1)
     end if
 
   end subroutine compute_dhat
diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90
index 32797a616f..32866b52ec 100644
--- a/src/cmfd_input.F90
+++ b/src/cmfd_input.F90
@@ -10,8 +10,6 @@ module cmfd_input
   private
   public :: configure_cmfd
 
-  character(2*MAX_LINE_LEN) :: message                                      
-
 contains
 
 !===============================================================================
@@ -87,15 +85,14 @@ contains
     if (.not. file_exists) then
       ! CMFD is optional unless it is in on from settings
       if (cmfd_on) then
-        message = "No CMFD XML file, '" // trim(filename) // "' does not exist!"
-        call fatal_error(message)
+        call fatal_error("No CMFD XML file, '" // trim(filename) // "' does not&
+             & exist!")
       end if
       return
     else
 
       ! Tell user
-      message = "Reading CMFD XML file..."
-      call write_message(message, 5)
+      call write_message("Reading CMFD XML file...", 5)
 
     end if
 
@@ -107,8 +104,7 @@ contains
 
     ! Check if mesh is there
     if (.not.found) then
-      message = "No CMFD mesh specified in CMFD XML file."
-      call fatal_error(message)
+      call fatal_error("No CMFD mesh specified in CMFD XML file.")
     end if
 
     ! Set spatial dimensions in cmfd object
@@ -139,8 +135,7 @@ contains
            cmfd % indices(3)))
       if (get_arraysize_integer(node_mesh, "map") /= &
           product(cmfd % indices(1:3))) then
-        message = 'FATAL==>CMFD coremap not to correct dimensions'
-        call fatal_error(message)
+        call fatal_error('FATAL==>CMFD coremap not to correct dimensions')
       end if
       allocate(iarray(get_arraysize_integer(node_mesh, "map")))
       call get_node_array(node_mesh, "map", iarray)
@@ -216,8 +211,7 @@ contains
       temp_str = to_lower(temp_str)
       if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
 #ifndef PETSC
-        message = 'Must use PETSc when running adjoint option.'
-        call fatal_error(message)
+        call fatal_error('Must use PETSc when running adjoint option.')
 #endif
         cmfd_run_adjoint = .true.
     end if
@@ -246,8 +240,8 @@ contains
          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'
-      if (master) call warning(message)
+      if (master) call warning('Dominance Ratio only aviable with power &
+           &iteration solver')
       cmfd_display = ''
     end if
 
@@ -263,9 +257,8 @@ contains
     if (check_for_node(doc, "gauss_seidel_tolerance")) then
       n_params = get_arraysize_double(doc, "gauss_seidel_tolerance")
       if (n_params /= 2) then
-        message = 'Gauss Seidel tolerance is not 2 parameters &
-                   &(absolute, relative).'
-        call fatal_error(message)
+        call fatal_error('Gauss Seidel tolerance is not 2 parameters &
+                   &(absolute, relative).')
       end if
       call get_node_array(doc, "gauss_seidel_tolerance", gs_tol)
       cmfd_atoli = gs_tol(1)
@@ -335,8 +328,7 @@ contains
     ! Determine number of dimensions for mesh
     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(message)
+       call fatal_error("Mesh must be two or three dimensions.")
     end if
     m % n_dimension = n
 
@@ -349,9 +341,8 @@ contains
     ! Check that dimensions are all greater than zero
     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(message)
+      call fatal_error("All entries on the <dimension> element for a tally mesh&
+           &must be positive.")
     end if
 
     ! Read dimensions in each direction
@@ -359,42 +350,37 @@ contains
 
     ! Read mesh lower-left corner location
     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(message)
+      call fatal_error("Number of entries on <lower_left> must be the same as &
+           &the number of entries on <dimension>.")
     end if
     call get_node_array(node_mesh, "lower_left", m % lower_left)
 
     ! 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(message)
+      call fatal_error("Cannot specify both <upper_right> and <width> on a &
+           &tally mesh.")
     end if
 
     ! Make sure either upper-right or width was specified
     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(message)
+      call fatal_error("Must specify either <upper_right> and <width> on a &
+           &tally mesh.")
     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(message)
+        call fatal_error("Number of entries on <width> must be the same as the &
+             &number of entries on <lower_left>.")
       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(message)
+        call fatal_error("Cannot have a negative <width> on a tally mesh.")
       end if
 
       ! Set width and upper right coordinate
@@ -405,17 +391,15 @@ contains
       ! 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(message)
+        call fatal_error("Number of entries on <upper_right> must be the same &
+             &as the number of entries on <lower_left>.")
       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(message)
+        call fatal_error("The <upper_right> coordinates must be greater than &
+             &the <lower_left> coordinates on a tally mesh.")
       end if
 
       ! Set upper right coordinate and width
diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90
index d884514ce1..e660c92470 100644
--- a/src/cmfd_solver.F90
+++ b/src/cmfd_solver.F90
@@ -101,7 +101,6 @@ contains
   subroutine init_data(adjoint)
 
     use constants, only: ONE, ZERO
-    use error,     only: fatal_error
     use global,    only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, &
                          cmfd_write_matrices
 
@@ -235,8 +234,8 @@ contains
 
       ! Check if reached iteration 10000
       if (i == 10000) then
-        call fatal_error('Reached maximum iterations in CMFD power iteration ' &
-            & // 'solver.')
+        call fatal_error('Reached maximum iterations in CMFD power iteration &
+             &solver.')
       end if
 
       ! Compute source vector
diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90
index 8605acb0f5..c46f90ab17 100644
--- a/src/eigenvalue.F90
+++ b/src/eigenvalue.F90
@@ -29,7 +29,6 @@ module eigenvalue
   real(8)                   :: keff_generation ! Single-generation k on each
                                                ! processor
   real(8)                   :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
-  character(2*MAX_LINE_LEN) :: message         ! Message to output unit
 
 contains
 
@@ -116,8 +115,8 @@ contains
 
   subroutine initialize_batch()
 
-    message = "Simulating batch " // trim(to_str(current_batch)) // "..."
-    call write_message(message, 8)
+    call write_message("Simulating batch " // trim(to_str(current_batch)) &
+         &// "...", 8)
 
     ! Reset total starting particle weight used for normalizing tallies
     total_weight = ZERO
@@ -303,8 +302,7 @@ contains
     ! runs enough particles to avoid this in the first place.
 
     if (n_bank == 0) then
-      message = "No fission sites banked on processor " // to_str(rank)
-      call fatal_error(message)
+      call fatal_error("No fission sites banked on processor " // to_str(rank))
     end if
 
     ! Make sure all processors start at the same point for random sampling. Then
@@ -555,8 +553,7 @@ contains
 
     ! display warning message if there were sites outside entropy box
     if (sites_outside) then
-      message = "Fission source site(s) outside of entropy box."
-      if (master) call warning(message)
+      if (master) call warning("Fission source site(s) outside of entropy box.")
     end if
 
     ! sum values to obtain shannon entropy
@@ -774,8 +771,7 @@ contains
 
       ! Check for sites outside of the mesh
       if (master .and. sites_outside) then
-        message = "Source sites outside of the UFS mesh!"
-        call fatal_error(message)
+        call fatal_error("Source sites outside of the UFS mesh!")
       end if
 
 #ifdef MPI
@@ -805,8 +801,7 @@ contains
 
     ! Write message at beginning
     if (current_batch == 1) then
-      message = "Replaying history from state point..."
-      call write_message(message, 1)
+      call write_message("Replaying history from state point...", 1)
     end if
 
     do current_gen = 1, gen_per_batch
@@ -823,8 +818,7 @@ contains
 
     ! Write message at end
     if (current_batch == restart_batch) then
-      message = "Resuming simulation..."
-      call write_message(message, 1)
+      call write_message("Resuming simulation...", 1)
     end if
 
   end subroutine replay_batch_history
diff --git a/src/geometry.F90 b/src/geometry.F90
index 1476fa2e39..043ded05b5 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -11,9 +11,7 @@ module geometry
   use tally,                  only: score_surface_current
 
   implicit none
-
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-     
+ 
 contains
 
 !===============================================================================
@@ -102,11 +100,10 @@ contains
         if (simple_cell_contains(c, p)) then
           ! the particle should only be contained in one cell per level
           if (index_cell /= coord % cell) then
-            message = "Overlapping cells detected: " //               &
-                      trim(to_str(cells(index_cell) % id)) // ", " // &
-                      trim(to_str(cells(coord % cell) % id)) //       &
-                      " on universe " // trim(to_str(univ % id))
-            call fatal_error(message)
+            call fatal_error("Overlapping cells detected: " &
+                 &// trim(to_str(cells(index_cell) % id)) // ", " &
+                 &// trim(to_str(cells(coord % cell) % id)) &
+                 &// " on universe " // trim(to_str(univ % id)))
           end if
 
           overlap_check_cnt(index_cell) = overlap_check_cnt(index_cell) + 1
@@ -181,8 +178,7 @@ contains
 
         ! Show cell information on trace
         if (verbosity >= 10 .or. trace) then
-          message = "    Entering cell " // trim(to_str(c % id))
-          call write_message(message)
+          call write_message("    Entering cell " // trim(to_str(c % id)))
         end if
 
         if (c % type == CELL_NORMAL) then
@@ -386,8 +382,7 @@ contains
     i_surface = abs(p % surface)
     surf => surfaces(i_surface)
     if (verbosity >= 10 .or. trace) then
-      message = "    Crossing surface " // trim(to_str(surf % id))
-      call write_message(message)
+      call write_message("    Crossing surface " // trim(to_str(surf % id)))
     end if
 
     if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
@@ -419,8 +414,8 @@ contains
 
       ! Display message
       if (verbosity >= 10 .or. trace) then
-        message = "    Leaked out of surface " // trim(to_str(surf % id))
-        call write_message(message)
+        call write_message("    Leaked out of surface " &
+             &// trim(to_str(surf % id)))
       end if
       return
 
@@ -430,9 +425,8 @@ contains
 
       ! Do not handle reflective boundary conditions on lower universes
       if (.not. associated(p % coord, p % coord0)) then
-        message = "Cannot reflect particle " // trim(to_str(p % id)) // &
-             " off surface in a lower universe."
-        call handle_lost_particle(p)
+        call handle_lost_particle(p, "Cannot reflect particle " &
+             &// trim(to_str(p % id)) // " off surface in a lower universe.")
         return
       end if
 
@@ -564,9 +558,8 @@ contains
         w = w + 2*dot_prod*R*z
 
       case default
-        message = "Reflection not supported for surface " // &
-             trim(to_str(surf % id))
-        call fatal_error(message)
+        call fatal_error("Reflection not supported for surface " &
+             &// trim(to_str(surf % id)))
       end select
 
       ! Set new particle direction
@@ -585,8 +578,8 @@ contains
         call deallocate_coord(p % coord0 % next)
         call find_cell(p, found)
         if (.not. found) then
-          message = "Couldn't find particle after reflecting from surface."
-          call handle_lost_particle(p)
+          call handle_lost_particle(p, "Couldn't find particle after reflecting&
+               & from surface.")
           return
         end if
       end if
@@ -596,8 +589,8 @@ contains
 
       ! Diagnostic message
       if (verbosity >= 10 .or. trace) then
-        message = "    Reflected from surface " // trim(to_str(surf%id))
-        call write_message(message)
+        call write_message("    Reflected from surface " &
+             &// trim(to_str(surf%id)))
       end if
       return
     end if
@@ -645,10 +638,9 @@ contains
       ! undefined region in the geometry.
 
       if (.not. found) then
-        message = "After particle " // trim(to_str(p % id)) // " crossed surface " &
-             // trim(to_str(surfaces(i_surface) % id)) // " it could not be &
-             &located in any cell and it did not leak."
-        call handle_lost_particle(p)
+        call handle_lost_particle(p, "After particle " // trim(to_str(p % id)) &
+             &// " crossed surface " // trim(to_str(surfaces(i_surface) % id)) &
+             &// " it could not be located in any cell and it did not leak.")
         return
       end if
     end if
@@ -674,11 +666,10 @@ contains
     lat => lattices(p % coord % lattice)
 
     if (verbosity >= 10 .or. trace) then
-      message = "    Crossing lattice " // trim(to_str(lat % id)) // &
-           ". Current position (" // trim(to_str(p % coord % lattice_x)) &
-           // "," // trim(to_str(p % coord % lattice_y)) // "," // &
-           trim(to_str(p % coord % lattice_z)) // ")"
-      call write_message(message)
+      call write_message("    Crossing lattice " // trim(to_str(lat % id)) &
+           &// ". Current position (" // trim(to_str(p % coord % lattice_x)) &
+           &// "," // trim(to_str(p % coord % lattice_y)) // "," &
+           &// trim(to_str(p % coord % lattice_z)) // ")")
     end if
 
     if (lat % type == LATTICE_RECT) then
@@ -741,9 +732,8 @@ contains
       ! Search for particle
       call find_cell(p, found)
       if (.not. found) then
-        message = "Could not locate particle " // trim(to_str(p % id)) // &
-             " after crossing a lattice boundary."
-        call handle_lost_particle(p)
+        call handle_lost_particle(p, "Could not locate particle " &
+             &// trim(to_str(p % id)) // " after crossing a lattice boundary.")
         return
       end if
     else
@@ -764,9 +754,9 @@ contains
         ! Search for particle
         call find_cell(p, found)
         if (.not. found) then
-          message = "Could not locate particle " // trim(to_str(p % id)) // &
-               " after crossing a lattice boundary."
-          call handle_lost_particle(p)
+          call handle_lost_particle(p, "Could not locate particle " &
+               &// trim(to_str(p % id)) &
+               &// " after crossing a lattice boundary.")
           return
         end if
       end if
@@ -1496,8 +1486,8 @@ contains
     type(Cell),    pointer  :: c
     type(Surface), pointer  :: surf
 
-    message = "Building neighboring cells lists for each surface..."
-    call write_message(message, 4)
+    call write_message("Building neighboring cells lists for each surface...", &
+         &4)
 
     allocate(count_positive(n_surfaces))
     allocate(count_negative(n_surfaces))
@@ -1564,9 +1554,10 @@ contains
 ! HANDLE_LOST_PARTICLE
 !===============================================================================
 
-  subroutine handle_lost_particle(p)
+  subroutine handle_lost_particle(p, message)
 
     type(Particle), intent(inout) :: p
+    character(*)                  :: message
 
     ! Print warning and write lost particle file
     call warning(message)
@@ -1581,8 +1572,7 @@ contains
     ! Abort the simulation if the maximum number of lost particles has been
     ! reached
     if (n_lost_particles == MAX_LOST_PARTICLES) then
-      message = "Maximum number of lost particles has been reached."
-      call fatal_error(message)
+      call fatal_error("Maximum number of lost particles has been reached.")
     end if
 
   end subroutine handle_lost_particle
diff --git a/src/initialize.F90 b/src/initialize.F90
index dd3b284c20..e0a5bf8b0c 100644
--- a/src/initialize.F90
+++ b/src/initialize.F90
@@ -38,8 +38,6 @@ module initialize
 
   implicit none
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 contains
 
 !===============================================================================
@@ -158,10 +156,8 @@ contains
 
     ! Warn if overlap checking is on
     if (master .and. check_overlaps) then
-      message = ""
-      call write_message(message)
-      message = "Cell overlap checking is ON"
-      if (master) call warning(message)
+      call write_message("")
+      call warning("Cell overlap checking is ON")
     end if
 
     ! Stop initialization timer
@@ -344,9 +340,8 @@ contains
 
           ! Check that number specified was valid
           if (n_particles == ERROR_INT) then
-            message = "Must specify integer after " // trim(argv(i-1)) // &
-                 " command-line flag."
-            call fatal_error(message)
+            call fatal_error("Must specify integer after " // trim(argv(i-1)) &
+                 &// " command-line flag.")
           end if
         case ('-r', '-restart', '--restart')
           ! Read path for state point/particle restart
@@ -366,8 +361,7 @@ contains
             path_particle_restart = argv(i)
             particle_restart_run = .true.
           case default
-            message = "Unrecognized file after restart flag."
-            call fatal_error(message)
+            call fatal_error("Unrecognized file after restart flag.")
           end select
 
           ! If its a restart run check for additional source file
@@ -385,8 +379,8 @@ contains
               call sp % read_data(filetype, 'filetype')
               call sp % file_close()
               if (filetype /= FILETYPE_SOURCE) then
-                message = "Second file after restart flag must be a source file"
-                call fatal_error(message)
+                call fatal_error("Second file after restart flag must be a &
+                     &source file")
               end if
 
               ! It is a source file
@@ -420,13 +414,13 @@ contains
           ! Read and set number of OpenMP threads
           n_threads = int(str_to_int(argv(i)), 4)
           if (n_threads < 1) then
-            message = "Invalid number of threads specified on command line."
-            call fatal_error(message)
+            call fatal_error("Invalid number of threads specified on command &
+                 &line.")
           end if
           call omp_set_num_threads(n_threads)
 #else
-          message = "Ignoring number of threads specified on command line."
-          if (master) call warning(message)
+          if (master) call warning("Ignoring number of threads specified on &
+               &command line.")
 #endif
 
         case ('-?', '-h', '-help', '--help')
@@ -442,8 +436,7 @@ contains
           write_all_tracks = .true.
           i = i + 1
         case default
-          message = "Unknown command line option: " // argv(i)
-          call fatal_error(message)
+          call fatal_error("Unknown command line option: " // argv(i))
         end select
 
         last_flag = i
@@ -579,9 +572,8 @@ contains
             i_array = surface_dict % get_key(abs(id))
             c % surfaces(j) = sign(i_array, id)
           else
-            message = "Could not find surface " // trim(to_str(abs(id))) // &
-                 " specified on cell " // trim(to_str(c % id))
-            call fatal_error(message)
+            call fatal_error("Could not find surface " // trim(to_str(abs(id)))&
+                 &// " specified on cell " // trim(to_str(c % id)))
           end if
         end if
       end do
@@ -593,9 +585,8 @@ contains
       if (universe_dict % has_key(id)) then
         c % universe = universe_dict % get_key(id)
       else
-        message = "Could not find universe " // trim(to_str(id)) // &
-             " specified on cell " // trim(to_str(c % id))
-        call fatal_error(message)
+        call fatal_error("Could not find universe " // trim(to_str(id)) &
+             &// " specified on cell " // trim(to_str(c % id)))
       end if
 
       ! =======================================================================
@@ -609,9 +600,8 @@ contains
           c % type = CELL_NORMAL
           c % material = material_dict % get_key(id)
         else
-          message = "Could not find material " // trim(to_str(id)) // &
-               " specified on cell " // trim(to_str(c % id))
-          call fatal_error(message)
+          call fatal_error("Could not find material " // trim(to_str(id)) &
+               &// " specified on cell " // trim(to_str(c % id)))
         end if
       else
         id = c % fill
@@ -628,15 +618,14 @@ contains
           else if (material_dict % has_key(mid)) then
             c % material = material_dict % get_key(mid)
           else
-            message = "Could not find material " // trim(to_str(mid)) // &
-               " specified on lattice " // trim(to_str(lid))
-            call fatal_error(message)
+            call fatal_error("Could not find material " // trim(to_str(mid)) &
+                 &// " specified on lattice " // trim(to_str(lid)))
           end if
 
         else
-          message = "Specified fill " // trim(to_str(id)) // " on cell " // &
-               trim(to_str(c % id)) // " is neither a universe nor a lattice."
-          call fatal_error(message)
+          call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
+               &// trim(to_str(c % id)) // " is neither a universe nor a &
+               &lattice.")
         end if
       end if
     end do
@@ -661,9 +650,8 @@ contains
             if (universe_dict % has_key(id)) then
               lat % universes(j,k,m) = universe_dict % get_key(id)
             else
-              message = "Invalid universe number " // trim(to_str(id)) &
-                   // " specified on lattice " // trim(to_str(lat % id))
-              call fatal_error(message)
+              call fatal_error("Invalid universe number " // trim(to_str(id)) &
+                   &// " specified on lattice " // trim(to_str(lat % id)))
             end if
           end do
         end do
@@ -687,9 +675,8 @@ contains
             if (cell_dict % has_key(id)) then
               t % filters(j) % int_bins(k) = cell_dict % get_key(id)
             else
-              message = "Could not find cell " // trim(to_str(id)) // &
-                   " specified on tally " // trim(to_str(t % id))
-              call fatal_error(message)
+              call fatal_error("Could not find cell " // trim(to_str(id)) &
+                   &// " specified on tally " // trim(to_str(t % id)))
             end if
           end do
 
@@ -703,9 +690,8 @@ contains
             if (surface_dict % has_key(id)) then
               t % filters(j) % int_bins(k) = surface_dict % get_key(id)
             else
-              message = "Could not find surface " // trim(to_str(id)) // &
-                   " specified on tally " // trim(to_str(t % id))
-              call fatal_error(message)
+              call fatal_error("Could not find surface " // trim(to_str(id)) &
+                   &// " specified on tally " // trim(to_str(t % id)))
             end if
           end do
 
@@ -716,9 +702,8 @@ contains
             if (universe_dict % has_key(id)) then
               t % filters(j) % int_bins(k) = universe_dict % get_key(id)
             else
-              message = "Could not find universe " // trim(to_str(id)) // &
-                   " specified on tally " // trim(to_str(t % id))
-              call fatal_error(message)
+              call fatal_error("Could not find universe " // trim(to_str(id)) &
+                   &// " specified on tally " // trim(to_str(t % id)))
             end if
           end do
 
@@ -729,9 +714,8 @@ contains
             if (material_dict % has_key(id)) then
               t % filters(j) % int_bins(k) = material_dict % get_key(id)
             else
-              message = "Could not find material " // trim(to_str(id)) // &
-                   " specified on tally " // trim(to_str(t % id))
-              call fatal_error(message)
+              call fatal_error("Could not find material " // trim(to_str(id)) &
+                   &// " specified on tally " // trim(to_str(t % id)))
             end if
           end do
 
@@ -866,8 +850,7 @@ contains
 
     ! Check for allocation errors
     if (alloc_err /= 0) then
-      message = "Failed to allocate source bank."
-      call fatal_error(message)
+      call fatal_error("Failed to allocate source bank.")
     end if
 
 #ifdef _OPENMP
@@ -894,8 +877,7 @@ contains
 
     ! Check for allocation errors
     if (alloc_err /= 0) then
-      message = "Failed to allocate fission bank."
-      call fatal_error(message)
+      call fatal_error("Failed to allocate fission bank.")
     end if
 
   end subroutine allocate_banks
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 7025c4b89b..869f7752c7 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -22,7 +22,6 @@ module input_xml
 
   type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
                                          ! universe contains
-  character(2*MAX_LINE_LEN) :: message   ! Message to output unit
 
 contains
 
@@ -77,19 +76,17 @@ contains
     type(NodeList), pointer :: node_scat_list => null()
 
     ! Display output message
-    message = "Reading settings XML file..."
-    call write_message(message, 5)
+    call write_message("Reading settings XML file...", 5)
 
     ! Check if settings.xml exists
     filename = trim(path_input) // "settings.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
-      message = "Settings XML file '" // trim(filename) // "' does not exist! &
-           &In order to run OpenMC, you first need a set of input files; at a &
-           &minimum, this includes settings.xml, geometry.xml, and &
+      call fatal_error("Settings XML file '" // trim(filename) // "' does not &
+           &exist! In order to run OpenMC, you first need a set of input files;&
+           & at a minimum, this includes settings.xml, geometry.xml, and &
            &materials.xml. Please consult the user's guide at &
-           &http://mit-crpg.github.io/openmc for further information."
-      call fatal_error(message)
+           &http://mit-crpg.github.io/openmc for further information.")
     end if
 
     ! Parse settings.xml file
@@ -105,13 +102,12 @@ contains
         ! environment variable
         call get_environment_variable("CROSS_SECTIONS", env_variable)
         if (len_trim(env_variable) == 0) then
-          message = "No cross_sections.xml file was specified in settings.xml &
-               &or in the CROSS_SECTIONS environment variable. OpenMC needs a &
-               &cross_sections.xml file to identify where to find ACE cross &
-               &section libraries. Please consult the user's guide at &
-               &http://mit-crpg.github.io/openmc for information on how to set &
-               &up ACE cross section libraries."
-          call fatal_error(message)
+          call fatal_error("No cross_sections.xml file was specified in &
+               &settings.xml or in the CROSS_SECTIONS environment variable. &
+               &OpenMC needs a cross_sections.xml file to identify where to &
+               &find ACE cross section libraries. Please consult the user's &
+               &guide at http://mit-crpg.github.io/openmc for information on &
+               &how to set up ACE cross section libraries.")
         else
           path_cross_sections = trim(env_variable)
         end if
@@ -131,8 +127,7 @@ contains
     ! Make sure that either eigenvalue or fixed source was specified
     if (.not.check_for_node(doc, "eigenvalue") .and. &
          .not.check_for_node(doc, "fixed_source")) then
-      message = "<eigenvalue> or <fixed_source> not specified."
-      call fatal_error(message)
+      call fatal_error("<eigenvalue> or <fixed_source> not specified.")
     end if
 
     ! Eigenvalue information
@@ -145,8 +140,7 @@ contains
 
       ! Check number of particles
       if (.not.check_for_node(node_mode, "particles")) then
-        message = "Need to specify number of particles per generation."
-        call fatal_error(message)
+        call fatal_error("Need to specify number of particles per generation.")
       end if
 
       ! Get number of particles
@@ -180,8 +174,7 @@ contains
 
       ! Check number of particles
       if (.not.check_for_node(node_mode, "particles")) then
-        message = "Need to specify number of particles per batch."
-        call fatal_error(message)
+        call fatal_error("Need to specify number of particles per batch.")
       end if
 
       ! Get number of particles
@@ -200,14 +193,11 @@ contains
 
     ! Check number of active batches, inactive batches, and particles
     if (n_active <= 0) then
-      message = "Number of active batches must be greater than zero."
-      call fatal_error(message)
+      call fatal_error("Number of active batches must be greater than zero.")
     elseif (n_inactive < 0) then
-      message = "Number of inactive batches must be non-negative."
-      call fatal_error(message)
+      call fatal_error("Number of inactive batches must be non-negative.")
     elseif (n_particles <= 0) then
-      message = "Number of particles must be greater than zero."
-      call fatal_error(message)
+      call fatal_error("Number of particles must be greater than zero.")
     end if
 
     ! Copy random number seed if specified
@@ -225,11 +215,9 @@ contains
     case ('union')
       grid_method = GRID_UNION
     case ('lethargy')
-      message = "Lethargy mapped energy grid not yet supported."
-      call fatal_error(message)
+      call fatal_error("Lethargy mapped energy grid not yet supported.")
     case default
-      message = "Unknown energy grid method: " // trim(temp_str)
-      call fatal_error(message)
+      call fatal_error("Unknown energy grid method: " // trim(temp_str))
     end select
 
     ! Verbosity
@@ -244,14 +232,12 @@ contains
       if (n_threads == NONE) then
         call get_node_value(doc, "threads", n_threads)
         if (n_threads < 1) then
-          message = "Invalid number of threads: " // to_str(n_threads)
-          call fatal_error(message)
+          call fatal_error("Invalid number of threads: " // to_str(n_threads))
         end if
         call omp_set_num_threads(n_threads)
       end if
 #else
-      message = "Ignoring number of threads."
-      if (master) call warning(message)
+      if (master) call warning("Ignoring number of threads.")
 #endif
     end if
 
@@ -262,8 +248,7 @@ contains
     if (check_for_node(doc, "source")) then
       call get_node_ptr(doc, "source", node_source)
     else
-      message = "No source specified in settings XML file."
-      call fatal_error(message)
+      call fatal_error("No source specified in settings XML file.")
     end if
 
     ! Check if we want to write out source
@@ -282,9 +267,8 @@ contains
       ! Check if source file exists
       inquire(FILE=path_source, EXIST=file_exists)
       if (.not. file_exists) then
-        message = "Binary source file '" // trim(path_source) // &
-             "' does not exist!"
-        call fatal_error(message)
+        call fatal_error("Binary source file '" // trim(path_source) &
+             &// "' does not exist!")
       end if
 
     else
@@ -310,9 +294,8 @@ contains
           external_source % type_space = SRC_SPACE_POINT
           coeffs_reqd = 3
         case default
-          message = "Invalid spatial distribution for external source: " &
-              // trim(type)
-          call fatal_error(message)
+          call fatal_error("Invalid spatial distribution for external source: "&
+               &// trim(type))
         end select
 
         ! Determine number of parameters specified
@@ -324,21 +307,19 @@ contains
 
         ! Read parameters for spatial distribution
         if (n < coeffs_reqd) then
-          message = "Not enough parameters specified for spatial &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Not enough parameters specified for spatial &
+               &distribution of external source.")
         elseif (n > coeffs_reqd) then
-          message = "Too many parameters specified for spatial &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Too many parameters specified for spatial &
+               &distribution of external source.")
         elseif (n > 0) then
           allocate(external_source % params_space(n))
           call get_node_array(node_dist, "parameters", &
                external_source % params_space)
         end if
       else
-        message = "No spatial distribution specified for external source."
-        call fatal_error(message)
+        call fatal_error("No spatial distribution specified for external &
+             &source.")
       end if
 
       ! Determine external source angular distribution
@@ -361,9 +342,8 @@ contains
         case ('tabular')
           external_source % type_angle = SRC_ANGLE_TABULAR
         case default
-          message = "Invalid angular distribution for external source: " &
-               // trim(type)
-          call fatal_error(message)
+          call fatal_error("Invalid angular distribution for external source: "&
+               &// trim(type))
         end select
 
         ! Determine number of parameters specified
@@ -375,13 +355,11 @@ contains
 
         ! Read parameters for angle distribution
         if (n < coeffs_reqd) then
-          message = "Not enough parameters specified for angle &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Not enough parameters specified for angle &
+               &distribution of external source.")
         elseif (n > coeffs_reqd) then
-          message = "Too many parameters specified for angle &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Too many parameters specified for angle &
+               &distribution of external source.")
         elseif (n > 0) then
           allocate(external_source % params_angle(n))
           call get_node_array(node_dist, "parameters", &
@@ -415,9 +393,8 @@ contains
         case ('tabular')
           external_source % type_energy = SRC_ENERGY_TABULAR
         case default
-          message = "Invalid energy distribution for external source: " &
-               // trim(type)
-          call fatal_error(message)
+          call fatal_error("Invalid energy distribution for external source: " &
+               &// trim(type))
         end select
 
         ! Determine number of parameters specified
@@ -429,13 +406,11 @@ contains
 
         ! Read parameters for energy distribution
         if (n < coeffs_reqd) then
-          message = "Not enough parameters specified for energy &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Not enough parameters specified for energy &
+               &distribution of external source.")
         elseif (n > coeffs_reqd) then
-          message = "Too many parameters specified for energy &
-               &distribution of external source."
-          call fatal_error(message)
+          call fatal_error("Too many parameters specified for energy &
+               &distribution of external source.")
         elseif (n > 0) then
           allocate(external_source % params_energy(n))
           call get_node_array(node_dist, "parameters", &
@@ -485,9 +460,9 @@ contains
       ! Make sure that there are three values per particle
       n_tracks = get_arraysize_integer(doc, "track")
       if (mod(n_tracks, 3) /= 0) then
-        message = "Number of integers specified in 'track' is not divisible &
-             &by 3.  Please provide 3 integers per particle to be tracked."
-        call fatal_error(message)
+        call fatal_error("Number of integers specified in 'track' is not &
+             &divisible by 3.  Please provide 3 integers per particle to be &
+             &tracked.")
       end if
 
       ! Allocate space and get list of tracks
@@ -507,11 +482,11 @@ contains
 
       ! Check to make sure enough values were supplied
       if (get_arraysize_double(node_entropy, "lower_left") /= 3) then
-        message = "Need to specify (x,y,z) coordinates of lower-left corner &
-             &of Shannon entropy mesh."
+        call fatal_error("Need to specify (x,y,z) coordinates of lower-left &
+             &corner of Shannon entropy mesh.")
       elseif (get_arraysize_double(node_entropy, "upper_right") /= 3) then
-        message = "Need to specify (x,y,z) coordinates of upper-right corner &
-             &of Shannon entropy mesh."
+        call fatal_error("Need to specify (x,y,z) coordinates of upper-right &
+             &corner of Shannon entropy mesh.")
       end if
 
       ! Allocate mesh object and coordinates on mesh
@@ -527,10 +502,10 @@ contains
            entropy_mesh % upper_right)
 
       ! Check on values provided
-      if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then
-        message = "Upper-right coordinate must be greater than lower-left &
-             &coordinate for Shannon entropy mesh."
-        call fatal_error(message)
+      if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) &
+           &then
+        call fatal_error("Upper-right coordinate must be greater than &
+             &lower-left coordinate for Shannon entropy mesh.")
       end if
 
       ! Check if dimensions were specified -- if not, they will be calculated
@@ -539,9 +514,8 @@ contains
 
         ! If so, make sure proper number of values were given
         if (get_arraysize_integer(node_entropy, "dimension") /= 3) then
-          message = "Dimension of entropy mesh must be given as three &
-               &integers."
-          call fatal_error(message)
+          call fatal_error("Dimension of entropy mesh must be given as three &
+               &integers.")
         end if
 
         ! Allocate dimensions
@@ -569,14 +543,14 @@ contains
 
       ! Check to make sure enough values were supplied
       if (get_arraysize_double(node_ufs, "lower_left") /= 3) then
-        message = "Need to specify (x,y,z) coordinates of lower-left corner &
-             &of UFS mesh."
+        call fatal_error("Need to specify (x,y,z) coordinates of lower-left &
+             &corner of UFS mesh.")
       elseif (get_arraysize_double(node_ufs, "upper_right") /= 3) then
-        message = "Need to specify (x,y,z) coordinates of upper-right corner &
-             &of UFS mesh."
+        call fatal_error("Need to specify (x,y,z) coordinates of upper-right &
+             &corner of UFS mesh.")
       elseif (get_arraysize_integer(node_ufs, "dimension") /= 3) then
-        message = "Dimension of UFS mesh must be given as three integers."
-        call fatal_error(message)
+        call fatal_error("Dimension of UFS mesh must be given as three &
+             &integers.")
       end if
 
       ! Allocate mesh object and coordinates on mesh
@@ -598,9 +572,8 @@ contains
 
       ! Check on values provided
       if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then
-        message = "Upper-right coordinate must be greater than lower-left &
-             &coordinate for UFS mesh."
-        call fatal_error(message)
+        call fatal_error("Upper-right coordinate must be greater than &
+             &lower-left coordinate for UFS mesh.")
       end if
 
       ! Calculate width
@@ -733,9 +706,8 @@ contains
       do i = 1, n_source_points
         if (.not. statepoint_batch % contains(sourcepoint_batch % &
             get_item(i))) then
-          message = 'Sourcepoint batches are not a subset&
-                    & of statepoint batches.'
-          call fatal_error(message)
+          call fatal_error('Sourcepoint batches are not a subset&
+               & of statepoint batches.')
         end if
       end do
     end if
@@ -815,9 +787,8 @@ contains
 
           ! check to make sure a nuclide is specified
           if (.not. check_for_node(node_scatterer, "nuclide")) then
-            message = "No nuclide specified for scatterer " // trim(to_str(i)) &
-              // " in settings.xml file!"
-            call fatal_error(message)
+            call fatal_error("No nuclide specified for scatterer " &
+                 &// trim(to_str(i)) // " in settings.xml file!")
           end if
           call get_node_value(node_scatterer, "nuclide", &
             nuclides_0K(i) % nuclide)
@@ -829,18 +800,17 @@ contains
 
           ! check to make sure xs name for which method is applied is given
           if (.not. check_for_node(node_scatterer, "xs_label")) then
-            message = "Must specify the temperature dependent name of " // '' &
-              //"scatterer " // trim(to_str(i)) // " given in cross_sections.xml"
-            call fatal_error(message)
+            call fatal_error("Must specify the temperature dependent name of &
+                 &scatterer " // trim(to_str(i)) &
+                 &// " given in cross_sections.xml")
           end if
           call get_node_value(node_scatterer, "xs_label", &
             nuclides_0K(i) % name)
 
           ! check to make sure 0K xs name for which method is applied is given
           if (.not. check_for_node(node_scatterer, "xs_label_0K")) then
-            message = "Must specify the 0K name of " // '' &
-              //"scatterer "// trim(to_str(i)) // " given in cross_sections.xml"
-            call fatal_error(message)
+            call fatal_error("Must specify the 0K name of scatterer " &
+                 &// trim(to_str(i)) // " given in cross_sections.xml")
           end if
           call get_node_value(node_scatterer, "xs_label_0K", &
             nuclides_0K(i) % name_0K)
@@ -852,8 +822,8 @@ contains
 
           ! check that E_min is non-negative
           if (nuclides_0K(i) % E_min < ZERO) then
-            message = "Lower resonance scattering energy bound is negative"
-            call fatal_error(message)
+            call fatal_error("Lower resonance scattering energy bound is &
+                 &negative")
           end if
 
           if (check_for_node(node_scatterer, "E_max")) then
@@ -863,8 +833,8 @@ contains
 
           ! check that E_max is not less than E_min
           if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then
-            message = "Lower resonance scattering energy bound exceeds upper"
-            call fatal_error(message)
+            call fatal_error("Lower resonance scattering energy bound exceeds &
+                 &upper")
           end if
 
           nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide)
@@ -873,9 +843,8 @@ contains
           nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K)
         end do
       else
-        message = "No resonant scatterers are specified within the " // "" &
-          // "resonance_scattering element in settings.xml"
-        call fatal_error(message)
+        call fatal_error("No resonant scatterers are specified within the &
+             &resonance_scattering element in settings.xml")
       end if
     end if
 
@@ -900,8 +869,8 @@ contains
       case ('jendl-4.0')
         default_expand = JENDL_40
       case default
-        message = "Unknown natural element expansion option: " // trim(temp_str)
-        call fatal_error(message)
+        call fatal_error("Unknown natural element expansion option: " &
+             &// trim(temp_str))
       end select
     end if
 
@@ -943,8 +912,7 @@ contains
     type(NodeList), pointer :: node_lat_list => null()
 
     ! Display output message
-    message = "Reading geometry XML file..."
-    call write_message(message, 5)
+    call write_message("Reading geometry XML file...", 5)
 
     ! ==========================================================================
     ! READ CELLS FROM GEOMETRY.XML
@@ -953,8 +921,8 @@ contains
     filename = trim(path_input) // "geometry.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
-      message = "Geometry XML file '" // trim(filename) // "' does not exist!"
-      call fatal_error(message)
+      call fatal_error("Geometry XML file '" // trim(filename) // "' does not &
+           &exist!")
     end if
 
     ! Parse geometry.xml file
@@ -968,8 +936,7 @@ contains
 
     ! Check for no cells
     if (n_cells == 0) then
-      message = "No cells found in geometry.xml!"
-      call fatal_error(message)
+      call fatal_error("No cells found in geometry.xml!")
     end if
 
     ! Allocate cells array
@@ -991,8 +958,7 @@ contains
       if (check_for_node(node_cell, "id")) then
         call get_node_value(node_cell, "id", c % id)
       else
-        message = "Must specify id of cell in geometry XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify id of cell in geometry XML file.")
       end if
       if (check_for_node(node_cell, "universe")) then
         call get_node_value(node_cell, "universe", c % universe)
@@ -1007,8 +973,8 @@ contains
 
       ! Check to make sure 'id' hasn't been used
       if (cell_dict % has_key(c % id)) then
-        message = "Two or more cells use the same unique ID: " // to_str(c % id)
-        call fatal_error(message)
+        call fatal_error("Two or more cells use the same unique ID: " &
+             &// to_str(c % id))
       end if
 
       ! Read material
@@ -1028,30 +994,27 @@ contains
 
         ! Check for error
         if (c % material == ERROR_INT) then
-          message = "Invalid material specified on cell " // to_str(c % id)
-          call fatal_error(message)
+          call fatal_error("Invalid material specified on cell " &
+               &// to_str(c % id))
         end if
       end select
 
       ! Check to make sure that either material or fill was specified
       if (c % material == NONE .and. c % fill == NONE) then
-        message = "Neither material nor fill was specified for cell " // &
-             trim(to_str(c % id))
-        call fatal_error(message)
+        call fatal_error("Neither material nor fill was specified for cell " &
+             &// trim(to_str(c % id)))
       end if
 
       ! Check to make sure that both material and fill haven't been
       ! specified simultaneously
       if (c % material /= NONE .and. c % fill /= NONE) then
-        message = "Cannot specify material and fill simultaneously"
-        call fatal_error(message)
+        call fatal_error("Cannot specify material and fill simultaneously")
       end if
 
       ! Check to make sure that surfaces were specified
       if (.not. check_for_node(node_cell, "surfaces")) then
-        message = "No surfaces specified for cell " // &
-             trim(to_str(c % id))
-        call fatal_error(message)
+        call fatal_error("No surfaces specified for cell " &
+             &// trim(to_str(c % id)))
       end if
 
       ! Allocate array for surfaces and copy
@@ -1065,17 +1028,15 @@ contains
         ! Rotations can only be applied to cells that are being filled with
         ! another universe
         if (c % fill == NONE) then
-          message = "Cannot apply a rotation to cell " // trim(to_str(&
-               c % id)) // " because it is not filled with another universe"
-          call fatal_error(message)
+          call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
+               &c % id)) // " because it is not filled with another universe")
         end if
 
         ! Read number of rotation parameters
         n = get_arraysize_double(node_cell, "rotation")
         if (n /= 3) then
-          message = "Incorrect number of rotation parameters on cell " // &
-               to_str(c % id)
-          call fatal_error(message)
+          call fatal_error("Incorrect number of rotation parameters on cell " &
+               &// to_str(c % id))
         end if
 
         ! Copy rotation angles in x,y,z directions
@@ -1101,17 +1062,16 @@ contains
         ! Translations can only be applied to cells that are being filled with
         ! another universe
         if (c % fill == NONE) then
-          message = "Cannot apply a translation to cell " // trim(to_str(&
-               c % id)) // " because it is not filled with another universe"
-          call fatal_error(message)
+          call fatal_error("Cannot apply a translation to cell " &
+               &// trim(to_str(c % id)) // " because it is not filled with &
+               &another universe")
         end if
 
         ! Read number of translation parameters
         n = get_arraysize_double(node_cell, "translation")
         if (n /= 3) then
-          message = "Incorrect number of translation parameters on cell " &
-               // to_str(c % id)
-          call fatal_error(message)
+          call fatal_error("Incorrect number of translation parameters on &
+               &cell " // to_str(c % id))
         end if
 
         ! Copy translation vector
@@ -1152,8 +1112,7 @@ contains
 
     ! Check for no surfaces
     if (n_surfaces == 0) then
-      message = "No surfaces found in geometry.xml!"
-      call fatal_error(message)
+      call fatal_error("No surfaces found in geometry.xml!")
     end if
 
     ! Allocate cells array
@@ -1169,15 +1128,13 @@ contains
       if (check_for_node(node_surf, "id")) then
         call get_node_value(node_surf, "id", s % id)
       else
-        message = "Must specify id of surface in geometry XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify id of surface in geometry XML file.")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (surface_dict % has_key(s % id)) then
-        message = "Two or more surfaces use the same unique ID: " // &
-             to_str(s % id)
-        call fatal_error(message)
+        call fatal_error("Two or more surfaces use the same unique ID: " &
+             &// to_str(s % id))
       end if
 
       ! Copy and interpret surface type
@@ -1219,8 +1176,7 @@ contains
         s % type = SURF_CONE_Z
         coeffs_reqd  = 4
       case default
-        message = "Invalid surface type: " // trim(word)
-        call fatal_error(message)
+        call fatal_error("Invalid surface type: " // trim(word))
       end select
 
       ! Check to make sure that the proper number of coefficients
@@ -1229,13 +1185,11 @@ contains
 
       n = get_arraysize_double(node_surf, "coeffs")
       if (n < coeffs_reqd) then
-        message = "Not enough coefficients specified for surface: " // &
-             trim(to_str(s % id))
-        call fatal_error(message)
+        call fatal_error("Not enough coefficients specified for surface: " &
+             &// trim(to_str(s % id)))
       elseif (n > coeffs_reqd) then
-        message = "Too many coefficients specified for surface: " // &
-             trim(to_str(s % id))
-        call fatal_error(message)
+        call fatal_error("Too many coefficients specified for surface: " &
+             &// trim(to_str(s % id)))
       else
         allocate(s % coeffs(n))
         call get_node_array(node_surf, "coeffs", s % coeffs)
@@ -1255,9 +1209,8 @@ contains
         s % bc = BC_REFLECT
         boundary_exists = .true.
       case default
-        message = "Unknown boundary condition '" // trim(word) // &
-             "' specified on surface " // trim(to_str(s % id))
-        call fatal_error(message)
+        call fatal_error("Unknown boundary condition '" // trim(word) // &
+             &"' specified on surface " // trim(to_str(s % id)))
       end select
 
       ! Add surface to dictionary
@@ -1268,8 +1221,7 @@ contains
     ! Check to make sure a boundary condition was applied to at least one
     ! surface
     if (.not. boundary_exists) then
-      message = "No boundary conditions were applied to any surfaces!"
-      call fatal_error(message)
+      call fatal_error("No boundary conditions were applied to any surfaces!")
     end if
 
     ! ==========================================================================
@@ -1292,15 +1244,13 @@ contains
       if (check_for_node(node_lat, "id")) then
         call get_node_value(node_lat, "id", lat % id)
       else
-        message = "Must specify id of lattice in geometry XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify id of lattice in geometry XML file.")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (lattice_dict % has_key(lat % id)) then
-        message = "Two or more lattices use the same unique ID: " // &
-             to_str(lat % id)
-        call fatal_error(message)
+        call fatal_error("Two or more lattices use the same unique ID: " &
+             &// to_str(lat % id))
       end if
 
       ! Read lattice type
@@ -1313,15 +1263,13 @@ contains
       case ('hex', 'hexagon', 'hexagonal')
         lat % type = LATTICE_HEX
       case default
-        message = "Invalid lattice type: " // trim(word)
-        call fatal_error(message)
+        call fatal_error("Invalid lattice type: " // trim(word))
       end select
 
       ! Read number of lattice cells in each dimension
       n = get_arraysize_integer(node_lat, "dimension")
       if (n /= 2 .and. n /= 3) then
-        message = "Lattice must be two or three dimensions."
-        call fatal_error(message)
+        call fatal_error("Lattice must be two or three dimensions.")
       end if
 
       lat % n_dimension = n
@@ -1331,9 +1279,8 @@ contains
       ! Read lattice lower-left location
       if (size(lat % dimension) /= &
           get_arraysize_double(node_lat, "lower_left")) then
-        message = "Number of entries on <lower_left> must be the same as &
-             &the number of entries on <dimension>."
-        call fatal_error(message)
+        call fatal_error("Number of entries on <lower_left> must be the same &
+             &as the number of entries on <dimension>.")
       end if
 
       allocate(lat % lower_left(n))
@@ -1342,9 +1289,8 @@ contains
       ! Read lattice widths
       if (size(lat % dimension) /= &
           get_arraysize_double(node_lat, "width")) then
-        message = "Number of entries on <width> must be the same as &
-             &the number of entries on <lower_left>."
-        call fatal_error(message)
+        call fatal_error("Number of entries on <width> must be the same as &
+             &the number of entries on <lower_left>.")
       end if
 
       allocate(lat % width(n))
@@ -1363,9 +1309,8 @@ contains
       ! Check that number of universes matches size
       n = get_arraysize_integer(node_lat, "universes")
       if (n /= n_x*n_y*n_z) then
-        message = "Number of universes on <universes> does not match size of &
-             &lattice " // trim(to_str(lat % id)) // "."
-        call fatal_error(message)
+        call fatal_error("Number of universes on <universes> does not match &
+             &size of lattice " // trim(to_str(lat % id)) // ".")
       end if
 
       allocate(temp_int_array(n))
@@ -1441,15 +1386,14 @@ contains
     type(NodeList), pointer :: node_sab_list => null()
 
     ! Display output message
-    message = "Reading materials XML file..."
-    call write_message(message, 5)
+    call write_message("Reading materials XML file...", 5)
 
     ! Check is materials.xml exists
     filename = trim(path_input) // "materials.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
-      message = "Material XML file '" // trim(filename) // "' does not exist!"
-      call fatal_error(message)
+      call fatal_error("Material XML file '" // trim(filename) // "' does not &
+           &exist!")
     end if
 
     ! Initialize default cross section variable
@@ -1483,15 +1427,13 @@ contains
       if (check_for_node(node_mat, "id")) then
         call get_node_value(node_mat, "id", mat % id)
       else
-        message = "Must specify id of material in materials XML file"
-        call fatal_error(message)
+        call fatal_error("Must specify id of material in materials XML file")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (material_dict % has_key(mat % id)) then
-        message = "Two or more materials use the same unique ID: " // &
-             to_str(mat % id)
-        call fatal_error(message)
+        call fatal_error("Two or more materials use the same unique ID: " &
+             &// to_str(mat % id))
       end if
 
       if (run_mode == MODE_PLOTTING) then
@@ -1507,9 +1449,8 @@ contains
       if (check_for_node(node_mat, "density")) then
         call get_node_ptr(node_mat, "density", node_dens)
       else
-        message = "Must specify density element in material " // &
-                  trim(to_str(mat % id))
-        call fatal_error(message)
+        call fatal_error("Must specify density element in material " &
+             &// trim(to_str(mat % id)))
       end if
 
       ! Initialize value to zero
@@ -1532,9 +1473,8 @@ contains
         ! Check for erroneous density
         sum_density = .false.
         if (val <= ZERO) then
-          message = "Need to specify a positive density on material " // &
-               trim(to_str(mat % id)) // "."
-          call fatal_error(message)
+          call fatal_error("Need to specify a positive density on material " &
+               &// trim(to_str(mat % id)) // ".")
         end if
 
         ! Adjust material density based on specified units
@@ -1548,9 +1488,8 @@ contains
         case ('atom/cm3', 'atom/cc')
           mat % density = 1.0e-24 * val
         case default
-          message = "Unkwown units '" // trim(units) &
-               // "' specified on material " // trim(to_str(mat % id))
-          call fatal_error(message)
+          call fatal_error("Unkwown units '" // trim(units) &
+               &// "' specified on material " // trim(to_str(mat % id)))
         end select
       end if
 
@@ -1560,9 +1499,8 @@ contains
       ! Check to ensure material has at least one nuclide
       if (.not. check_for_node(node_mat, "nuclide") .and. &
            .not. check_for_node(node_mat, "element")) then
-        message = "No nuclides or natural elements specified on material " // &
-             trim(to_str(mat % id))
-        call fatal_error(message)
+        call fatal_error("No nuclides or natural elements specified on &
+             &material " // trim(to_str(mat % id)))
       end if
 
       ! Get pointer list of XML <nuclide>
@@ -1575,17 +1513,15 @@ contains
 
         ! Check for empty name on nuclide
         if (.not.check_for_node(node_nuc, "name")) then
-          message = "No name specified on nuclide in material " // &
-               trim(to_str(mat % id))
-          call fatal_error(message)
+          call fatal_error("No name specified on nuclide in material " &
+               &// trim(to_str(mat % id)))
         end if
 
         ! Check for cross section
         if (.not.check_for_node(node_nuc, "xs")) then
           if (default_xs == '') then
-            message = "No cross section specified for nuclide in material " &
-                 // trim(to_str(mat % id))
-            call fatal_error(message)
+            call fatal_error("No cross section specified for nuclide in &
+                 &material " // trim(to_str(mat % id)))
           else
             name = trim(default_xs)
           end if
@@ -1604,14 +1540,12 @@ contains
         ! weight percents were specified
         if (.not.check_for_node(node_nuc, "ao") .and. &
             .not.check_for_node(node_nuc, "wo")) then
-          message = "No atom or weight percent specified for nuclide " // &
-               trim(name)
-          call fatal_error(message)
+          call fatal_error("No atom or weight percent specified for nuclide " &
+               &// trim(name))
         elseif (check_for_node(node_nuc, "ao") .and. &
                 check_for_node(node_nuc, "wo")) then
-          message = "Cannot specify both atom and weight percents for a &
-               &nuclide: " // trim(name)
-          call fatal_error(message)
+          call fatal_error("Cannot specify both atom and weight percents for a &
+               &nuclide: " // trim(name))
         end if
 
         ! Copy atom/weight percents
@@ -1635,9 +1569,8 @@ contains
 
         ! Check for empty name on natural element
         if (.not.check_for_node(node_ele, "name")) then
-          message = "No name specified on nuclide in material " // &
-               trim(to_str(mat % id))
-          call fatal_error(message)
+          call fatal_error("No name specified on nuclide in material " &
+               &// trim(to_str(mat % id)))
         end if
         call get_node_value(node_ele, "name", name)
 
@@ -1646,9 +1579,8 @@ contains
           call get_node_value(node_ele, "xs", temp_str)
         else
           if (default_xs == '') then
-            message = "No cross section specified for nuclide in material " &
-                 // trim(to_str(mat % id))
-            call fatal_error(message)
+            call fatal_error("No cross section specified for nuclide in &
+                 &material " // trim(to_str(mat % id)))
           else
             temp_str = trim(default_xs)
           end if
@@ -1658,14 +1590,12 @@ contains
         ! weight percents were specified
         if (.not.check_for_node(node_ele, "ao") .and. &
             .not.check_for_node(node_ele, "wo")) then
-          message = "No atom or weight percent specified for element " // &
-               trim(name)
-          call fatal_error(message)
+          call fatal_error("No atom or weight percent specified for element " &
+               &// trim(name))
         elseif (check_for_node(node_ele, "ao") .and. &
                 check_for_node(node_ele, "wo")) then
-          message = "Cannot specify both atom and weight percents for a &
-               &element: " // trim(name)
-          call fatal_error(message)
+          call fatal_error("Cannot specify both atom and weight percents for &
+               &element: " // trim(name))
         end if
 
         ! Expand element into naturally-occurring isotopes
@@ -1674,9 +1604,8 @@ contains
           call expand_natural_element(name, temp_str, temp_dble, &
                list_names, list_density)
         else
-          message = "The ability to expand a natural element based on weight &
-               &percentage is not yet supported."
-          call fatal_error(message)
+          call fatal_error("The ability to expand a natural element based on &
+               &weight percentage is not yet supported.")
         end if
       end do NATURAL_ELEMENTS
 
@@ -1694,17 +1623,15 @@ contains
         ! Check that this nuclide is listed in the cross_sections.xml file
         name = trim(list_names % get_item(j))
         if (.not. xs_listing_dict % has_key(to_lower(name))) then
-          message = "Could not find nuclide " // trim(name) // &
-               " in cross_sections.xml file!"
-          call fatal_error(message)
+          call fatal_error("Could not find nuclide " // trim(name) &
+               &// " in cross_sections.xml file!")
         end if
 
         ! Check to make sure cross-section is continuous energy neutron table
         n = len_trim(name)
         if (name(n:n) /= 'c') then
-          message = "Cross-section table " // trim(name) // &
-               " is not a continuous-energy neutron table."
-          call fatal_error(message)
+          call fatal_error("Cross-section table " // trim(name) &
+               &// " is not a continuous-energy neutron table.")
         end if
 
         ! Find xs_listing and set the name/alias according to the listing
@@ -1733,9 +1660,8 @@ contains
       ! given
       if (.not. (all(mat % atom_density >= ZERO) .or. &
            all(mat % atom_density <= ZERO))) then
-        message = "Cannot mix atom and weight percents in material " // &
-             to_str(mat % id)
-        call fatal_error(message)
+        call fatal_error("Cannot mix atom and weight percents in material " &
+             &// to_str(mat % id))
       end if
 
       ! Determine density if it is a sum value
@@ -1771,8 +1697,8 @@ contains
           ! Determine name of S(a,b) table
           if (.not.check_for_node(node_sab, "name") .or. &
               .not.check_for_node(node_sab, "xs")) then
-            message = "Need to specify <name> and <xs> for S(a,b) table."
-            call fatal_error(message)
+            call fatal_error("Need to specify <name> and <xs> for S(a,b) &
+                 &table.")
           end if
           call get_node_value(node_sab, "name", name)
           call get_node_value(node_sab, "xs", temp_str)
@@ -1781,9 +1707,8 @@ contains
 
           ! Check that this nuclide is listed in the cross_sections.xml file
           if (.not. xs_listing_dict % has_key(to_lower(name))) then
-            message = "Could not find S(a,b) table " // trim(name) // &
-                 " in cross_sections.xml file!"
-            call fatal_error(message)
+            call fatal_error("Could not find S(a,b) table " // trim(name) &
+                 &// " in cross_sections.xml file!")
           end if
 
           ! Find index in xs_listing and set the name and alias according to the
@@ -1868,8 +1793,7 @@ contains
     end if
 
     ! Display output message
-    message = "Reading tallies XML file..."
-    call write_message(message, 5)
+    call write_message("Reading tallies XML file...", 5)
 
     ! Parse tallies.xml file
     call open_xmldoc(doc, filename)
@@ -1897,8 +1821,7 @@ contains
     ! Check for user tallies
     n_user_tallies = get_list_size(node_tal_list)
     if (n_user_tallies == 0) then
-      message = "No tallies present in tallies.xml file!"
-      if (master) call warning(message)
+      if (master) call warning("No tallies present in tallies.xml file!")
     end if
 
     ! Allocate tally array
@@ -1927,15 +1850,13 @@ contains
       if (check_for_node(node_mesh, "id")) then
         call get_node_value(node_mesh, "id", m % id)
       else
-        message = "Must specify id for mesh in tally XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify id for mesh in tally XML file.")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (mesh_dict % has_key(m % id)) then
-        message = "Two or more meshes use the same unique ID: " // &
-             to_str(m % id)
-        call fatal_error(message)
+        call fatal_error("Two or more meshes use the same unique ID: " &
+             &// to_str(m % id))
       end if
 
       ! Read mesh type
@@ -1948,15 +1869,13 @@ contains
       case ('hex', 'hexagon', 'hexagonal')
         m % type = LATTICE_HEX
       case default
-        message = "Invalid mesh type: " // trim(temp_str)
-        call fatal_error(message)
+        call fatal_error("Invalid mesh type: " // trim(temp_str))
       end select
 
       ! Determine number of dimensions for mesh
       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(message)
+        call fatal_error("Mesh must be two or three dimensions.")
       end if
       m % n_dimension = n
 
@@ -1969,9 +1888,8 @@ contains
       ! Check that dimensions are all greater than zero
       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(message)
+        call fatal_error("All entries on the <dimension> element for a tally &
+             &mesh must be positive.")
       end if
 
       ! Read dimensions in each direction
@@ -1979,42 +1897,37 @@ contains
 
       ! Read mesh lower-left corner location
       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(message)
+        call fatal_error("Number of entries on <lower_left> must be the same &
+             &as the number of entries on <dimension>.")
       end if
       call get_node_array(node_mesh, "lower_left", m % lower_left)
 
       ! 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(message)
+        call fatal_error("Cannot specify both <upper_right> and <width> on a &
+             &tally mesh.")
       end if
 
       ! Make sure either upper-right or width was specified
       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(message)
+        call fatal_error("Must specify either <upper_right> and <width> on a &
+             &tally mesh.")
       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(message)
+          call fatal_error("Number of entries on <width> must be the same as &
+               &the number of entries on <lower_left>.")
         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(message)
+          call fatal_error("Cannot have a negative <width> on a tally mesh.")
         end if
 
         ! Set width and upper right coordinate
@@ -2025,17 +1938,15 @@ contains
         ! 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(message)
+          call fatal_error("Number of entries on <upper_right> must be the &
+               &same as the number of entries on <lower_left>.")
         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(message)
+          call fatal_error("The <upper_right> coordinates must be greater than &
+               &the <lower_left> coordinates on a tally mesh.")
         end if
 
         ! Set width and upper right coordinate
@@ -2078,15 +1989,13 @@ contains
       if (check_for_node(node_tal, "id")) then
         call get_node_value(node_tal, "id", t % id)
       else
-        message = "Must specify id for tally in tally XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify id for tally in tally XML file.")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (tally_dict % has_key(t % id)) then
-        message = "Two or more tallies use the same unique ID: " // &
-             to_str(t % id)
-        call fatal_error(message)
+        call fatal_error("Two or more tallies use the same unique ID: " &
+             &// to_str(t % id))
       end if
 
       ! Copy tally label
@@ -2102,9 +2011,9 @@ contains
       ! the old format and if it is present, raises an error
 
 !     if (get_number_nodes(node_tal, "filters") > 0) then
-!       message = "Tally filters should be specified with multiple <filter> &
-!            &elements. Did you forget to change your <filters> element?"
-!       call fatal_error(message)
+!       call fatal_error("Tally filters should be specified with multiple &
+!            &<filter> elements. Did you forget to change your <filters> &
+!            &element?")
 !     end if
 
       ! Get pointer list to XML <filter> and get number of filters
@@ -2136,8 +2045,8 @@ contains
               n_words = get_arraysize_integer(node_filt, "bins")
             end if
           else
-            message = "Bins not set in filter on tally " // trim(to_str(t % id))
-            call fatal_error(message)
+            call fatal_error("Bins not set in filter on tally " &
+                 &// trim(to_str(t % id)))
           end if
 
           ! Determine type of filter
@@ -2187,8 +2096,7 @@ contains
             call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
 
           case ('surface')
-            message = "Surface filter is not yet supported!"
-            call fatal_error(message)
+            call fatal_error("Surface filter is not yet supported!")
 
             ! Set type of filter
             t % filters(j) % type = FILTER_SURFACE
@@ -2206,8 +2114,7 @@ contains
 
             ! Check to make sure multiple meshes weren't given
             if (n_words /= 1) then
-              message = "Can only have one mesh filter specified."
-              call fatal_error(message)
+              call fatal_error("Can only have one mesh filter specified.")
             end if
 
             ! Determine id of mesh
@@ -2218,9 +2125,8 @@ contains
               i_mesh = mesh_dict % get_key(id)
               m => meshes(i_mesh)
             else
-              message = "Could not find mesh " // trim(to_str(id)) // &
-                   " specified on tally " // trim(to_str(t % id))
-              call fatal_error(message)
+              call fatal_error("Could not find mesh " // trim(to_str(id)) &
+                   &// " specified on tally " // trim(to_str(t % id)))
             end if
 
             ! Determine number of bins -- this is assuming that the tally is
@@ -2259,10 +2165,9 @@ contains
 
           case default
             ! Specified tally filter is invalid, raise error
-            message = "Unknown filter type '" // &
-                 trim(temp_str) // "' on tally " // &
-                 trim(to_str(t % id)) // "."
-            call fatal_error(message)
+            call fatal_error("Unknown filter type '" &
+                 &// trim(temp_str) // "' on tally " &
+                 &// trim(to_str(t % id)) // ".")
 
           end select
 
@@ -2276,9 +2181,8 @@ contains
         ! Check that both cell and surface weren't specified
         if (t % find_filter(FILTER_CELL) > 0 .and. &
              t % find_filter(FILTER_SURFACE) > 0) then
-          message = "Cannot specify both cell and surface filters for tally " &
-               // trim(to_str(t % id))
-          call fatal_error(message)
+          call fatal_error("Cannot specify both cell and surface filters for &
+               &tally " // trim(to_str(t % id)))
         end if
 
       else
@@ -2340,10 +2244,9 @@ contains
 
                 ! Check if no nuclide was found
                 if (.not. associated(pair_list)) then
-                  message = "Could not find the nuclide " // trim(&
-                       sarray(j)) // " specified in tally " &
-                       // trim(to_str(t % id)) // " in any material."
-                  call fatal_error(message)
+                  call fatal_error("Could not find the nuclide " &
+                       &// trim(sarray(j)) // " specified in tally " &
+                       &// trim(to_str(t % id)) // " in any material.")
                 end if
                 deallocate(pair_list)
               else
@@ -2354,9 +2257,9 @@ contains
 
             ! Check to make sure nuclide specified is in problem
             if (.not. nuclide_dict % has_key(to_lower(word))) then
-              message = "The nuclide " // trim(word) // " from tally " // &
-                   trim(to_str(t % id)) // " is not present in any material."
-              call fatal_error(message)
+              call fatal_error("The nuclide " // trim(word) // " from tally " &
+                   &// trim(to_str(t % id)) &
+                   &// " is not present in any material.")
             end if
 
             ! Set bin to index in nuclides array
@@ -2405,13 +2308,13 @@ contains
                 ! User requested too many orders; throw a warning and set to the
                 ! maximum order.
                 ! The above scheme will essentially take the absolute value
-                message = "Invalid scattering order of " // trim(to_str(n_order)) // &
-                  " requested. Setting to the maximum permissible value, " // &
-                  trim(to_str(MAX_ANG_ORDER))
-                if (master) call warning(message)
+                if (master) call warning("Invalid scattering order of " &
+                     &// trim(to_str(n_order)) // " requested. Setting to the &
+                     &maximum permissible value, " &
+                     &// trim(to_str(MAX_ANG_ORDER)))
                 n_order = MAX_ANG_ORDER
-                sarray(j) = trim(MOMENT_STRS(imomstr)) // &
-                  trim(to_str(MAX_ANG_ORDER))
+                sarray(j) = trim(MOMENT_STRS(imomstr)) &
+                     &// trim(to_str(MAX_ANG_ORDER))
               end if
               ! Find total number of bins for this case
               if (imomstr >= YN_LOC) then
@@ -2436,7 +2339,8 @@ contains
           j = j + 1
           ! Get the input string in scores(l) but if score is one of the moment
           ! scores then strip off the n and store it as an integer to be used
-          ! later. Then perform the select case on this modified (number removed) string
+          ! later. Then perform the select case on this modified (number 
+          ! removed) string
           score_name = sarray(l)
           do imomstr = 1, size(MOMENT_STRS)
             if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
@@ -2447,10 +2351,6 @@ contains
                 ! User requested too many orders; throw a warning and set to the
                 ! maximum order.
                 ! The above scheme will essentially take the absolute value
-                message = "Invalid scattering order of " // trim(to_str(n_order)) // &
-                  " requested. Setting to the maximum permissible value, " // &
-                  trim(to_str(MAX_ANG_ORDER))
-                if (master) call warning(message)
                 n_order = MAX_ANG_ORDER
               end if
               score_name = trim(MOMENT_STRS(imomstr)) // "n"
@@ -2474,10 +2374,10 @@ contains
                   ! User requested too many orders; throw a warning and set to the
                   ! maximum order.
                   ! The above scheme will essentially take the absolute value
-                  message = "Invalid scattering order of " // trim(to_str(n_order)) // &
-                    " requested. Setting to the maximum permissible value, " // &
-                    trim(to_str(MAX_ANG_ORDER))
-                  if (master) call warning(message)
+                  if (master) call warning("Invalid scattering order of " &
+                       &// trim(to_str(n_order)) // " requested. Setting to &
+                       &the maximum permissible value, " &
+                       &// trim(to_str(MAX_ANG_ORDER)))
                   n_order = MAX_ANG_ORDER
                 end if
                 score_name = trim(MOMENT_N_STRS(imomstr)) // "n"
@@ -2491,26 +2391,24 @@ contains
             ! Prohibit user from tallying flux for an individual nuclide
             if (.not. (t % n_nuclide_bins == 1 .and. &
                  t % nuclide_bins(1) == -1)) then
-              message = "Cannot tally flux for an individual nuclide."
-              call fatal_error(message)
+              call fatal_error("Cannot tally flux for an individual nuclide.")
             end if
 
             t % score_bins(j) = SCORE_FLUX
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally flux with an outgoing energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally flux with an outgoing energy &
+                   &filter.")
             end if
           case ('flux-yn')
             ! Prohibit user from tallying flux for an individual nuclide
             if (.not. (t % n_nuclide_bins == 1 .and. &
                  t % nuclide_bins(1) == -1)) then
-              message = "Cannot tally flux for an individual nuclide."
-              call fatal_error(message)
+              call fatal_error("Cannot tally flux for an individual nuclide.")
             end if
 
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally flux with an outgoing energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally flux with an outgoing energy &
+                   &filter.")
             end if
 
             t % score_bins(j : j + n_bins - 1) = SCORE_FLUX_YN
@@ -2520,16 +2418,14 @@ contains
           case ('total')
             t % score_bins(j) = SCORE_TOTAL
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally total reaction rate with an &
-                   &outgoing energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally total reaction rate with an &
+                   &outgoing energy filter.")
             end if
 
           case ('total-yn')
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally total reaction rate with an &
-                   &outgoing energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally total reaction rate with an &
+                   &outgoing energy filter.")
             end if
 
             t % score_bins(j : j + n_bins - 1) = SCORE_TOTAL_YN
@@ -2598,9 +2494,8 @@ contains
             ! Set tally estimator to analog
             t % estimator = ESTIMATOR_ANALOG
           case ('diffusion')
-            message = "Diffusion score no longer supported for tallies, &
-                      &please remove"
-            call fatal_error(message)
+            call fatal_error("Diffusion score no longer supported for tallies, &
+                 &please remove")
           case ('n1n')
             t % score_bins(j) = SCORE_N_1N
 
@@ -2618,16 +2513,14 @@ contains
           case ('absorption')
             t % score_bins(j) = SCORE_ABSORPTION
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally absorption rate with an outgoing &
-                   &energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally absorption rate with an outgoing &
+                   &energy filter.")
             end if
           case ('fission')
             t % score_bins(j) = SCORE_FISSION
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              message = "Cannot tally fission rate with an outgoing &
-                   &energy filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally fission rate with an outgoing &
+                   &energy filter.")
             end if
           case ('nu-fission')
             t % score_bins(j) = SCORE_NU_FISSION
@@ -2644,10 +2537,9 @@ contains
             ! Check to make sure that current is the only desired response
             ! for this tally
             if (n_words > 1) then
-              message = "Cannot tally other scoring functions in the same &
-                   &tally as surface currents. Separate other scoring &
-                   &functions into a distinct tally."
-              call fatal_error(message)
+              call fatal_error("Cannot tally other scoring functions in the &
+                   &same tally as surface currents. Separate other scoring &
+                   &functions into a distinct tally.")
             end if
 
             ! Since the number of bins for the mesh filter was already set
@@ -2659,8 +2551,8 @@ contains
 
             ! Check to make sure mesh filter was specified
             if (k == 0) then
-              message = "Cannot tally surface current without a mesh filter."
-              call fatal_error(message)
+              call fatal_error("Cannot tally surface current without a mesh &
+                   &filter.")
             end if
 
             ! Get pointer to mesh
@@ -2706,16 +2598,14 @@ contains
               if (MT > 1) then
                 t % score_bins(j) = MT
               else
-                message = "Invalid MT on <scores>: " // &
-                     trim(sarray(l))
-                call fatal_error(message)
+                call fatal_error("Invalid MT on <scores>: " &
+                     &// trim(sarray(l)))
               end if
 
             else
               ! Specified score was not an integer
-              message = "Unknown scoring function: " // &
-                   trim(sarray(l))
-              call fatal_error(message)
+              call fatal_error("Unknown scoring function: " &
+                   &// trim(sarray(l)))
             end if
 
           end select
@@ -2726,9 +2616,8 @@ contains
         ! Deallocate temporary string array of scores
         deallocate(sarray)
       else
-        message = "No <scores> specified on tally " // trim(to_str(t % id)) &
-             // "."
-        call fatal_error(message)
+        call fatal_error("No <scores> specified on tally " &
+             &// trim(to_str(t % id)) // ".")
       end if
 
       ! =======================================================================
@@ -2746,18 +2635,16 @@ contains
           ! If the estimator was set to an analog estimator, this means the
           ! tally needs post-collision information
           if (t % estimator == ESTIMATOR_ANALOG) then
-            message = "Cannot use track-length estimator for tally " &
-                 // to_str(t % id)
-            call fatal_error(message)
+            call fatal_error("Cannot use track-length estimator for tally " &
+                 &// to_str(t % id))
           end if
 
           ! Set estimator to track-length estimator
           t % estimator = ESTIMATOR_TRACKLENGTH
 
         case default
-          message = "Invalid estimator '" // trim(temp_str) &
-               // "' on tally " // to_str(t % id)
-          call fatal_error(message)
+          call fatal_error("Invalid estimator '" // trim(temp_str) &
+               &// "' on tally " // to_str(t % id))
         end select
       end if
 
@@ -2801,13 +2688,12 @@ contains
     filename = trim(path_input) // "plots.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
-      message = "Plots XML file '" // trim(filename) // "' does not exist!"
-      call fatal_error(message)
+      call fatal_error("Plots XML file '" // trim(filename) &
+           &// "' does not exist!")
     end if
 
     ! Display output message
-    message = "Reading plot XML file..."
-    call write_message(message, 5)
+    call write_message("Reading plot XML file...", 5)
 
     ! Parse plots.xml file
     call open_xmldoc(doc, filename)
@@ -2829,15 +2715,13 @@ contains
       if (check_for_node(node_plot, "id")) then
         call get_node_value(node_plot, "id", pl % id)
       else
-        message = "Must specify plot id in plots XML file."
-        call fatal_error(message)
+        call fatal_error("Must specify plot id in plots XML file.")
       end if
 
       ! Check to make sure 'id' hasn't been used
       if (plot_dict % has_key(pl % id)) then
-        message = "Two or more plots use the same unique ID: " // &
-             to_str(pl % id)
-        call fatal_error(message)
+        call fatal_error("Two or more plots use the same unique ID: " &
+             &// to_str(pl % id))
       end if
 
       ! Copy plot type
@@ -2851,9 +2735,8 @@ contains
       case ("voxel")
         pl % type = PLOT_TYPE_VOXEL
       case default
-        message = "Unsupported plot type '" // trim(temp_str) &
-             // "' in plot " // trim(to_str(pl % id))
-        call fatal_error(message)
+        call fatal_error("Unsupported plot type '" // trim(temp_str) &
+             &// "' in plot " // trim(to_str(pl % id)))
       end select
 
       ! Set output file path
@@ -2874,33 +2757,29 @@ contains
         if (get_arraysize_integer(node_plot, "pixels") == 2) then
           call get_node_array(node_plot, "pixels", pl % pixels(1:2))
         else
-          message = "<pixels> must be length 2 in slice plot " // &
-                    trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("<pixels> must be length 2 in slice plot " &
+               &// trim(to_str(pl % id)))
         end if
       else if (pl % type == PLOT_TYPE_VOXEL) then
         if (get_arraysize_integer(node_plot, "pixels") == 3) then
           call get_node_array(node_plot, "pixels", pl % pixels(1:3))
         else
-          message = "<pixels> must be length 3 in voxel plot " // &
-                    trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("<pixels> must be length 3 in voxel plot " &
+               &// trim(to_str(pl % id)))
         end if
       end if
 
       ! Copy plot background color
       if (check_for_node(node_plot, "background")) then
         if (pl % type == PLOT_TYPE_VOXEL) then
-          message = "Background color ignored in voxel plot " // &
-                     trim(to_str(pl % id))
-          if (master) call warning(message)
+          if (master) call warning("Background color ignored in voxel plot " &
+               &// trim(to_str(pl % id)))
         end if
         if (get_arraysize_integer(node_plot, "background") == 3) then
           call get_node_array(node_plot, "background", pl % not_found % rgb)
         else
-          message = "Bad background RGB " &
-               // "in plot " // trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("Bad background RGB in plot " &
+               &// trim(to_str(pl % id)))
         end if
       else
         pl % not_found % rgb = (/ 255, 255, 255 /)
@@ -2920,9 +2799,8 @@ contains
         case ("yz")
           pl % basis = PLOT_BASIS_YZ
         case default
-          message = "Unsupported plot basis '" // trim(temp_str) &
-               // "' in plot " // trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("Unsupported plot basis '" // trim(temp_str) &
+               &// "' in plot " // trim(to_str(pl % id)))
         end select
       end if
 
@@ -2930,9 +2808,8 @@ contains
       if (get_arraysize_double(node_plot, "origin") == 3) then
         call get_node_array(node_plot, "origin", pl % origin)
       else
-        message = "Origin must be length 3 " &
-             // "in plot " // trim(to_str(pl % id))
-        call fatal_error(message)
+        call fatal_error("Origin must be length 3 in plot " &
+             &// trim(to_str(pl % id)))
       end if
 
       ! Copy plotting width
@@ -2940,17 +2817,15 @@ contains
         if (get_arraysize_double(node_plot, "width") == 2) then
           call get_node_array(node_plot, "width", pl % width(1:2))
         else
-          message = "<width> must be length 2 in slice plot " // &
-                    trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("<width> must be length 2 in slice plot " &
+               &// trim(to_str(pl % id)))
         end if
       else if (pl % type == PLOT_TYPE_VOXEL) then
         if (get_arraysize_double(node_plot, "width") == 3) then
           call get_node_array(node_plot, "width", pl % width(1:3))
         else
-          message = "<width> must be length 3 in voxel plot " // &
-                    trim(to_str(pl % id))
-          call fatal_error(message)
+          call fatal_error("<width> must be length 3 in voxel plot " &
+               &// trim(to_str(pl % id)))
         end if
       end if
 
@@ -2981,9 +2856,8 @@ contains
         end do
 
       case default
-        message = "Unsupported plot color type '" // trim(temp_str) &
-             // "' in plot " // trim(to_str(pl % id))
-        call fatal_error(message)
+        call fatal_error("Unsupported plot color type '" // trim(temp_str) &
+             &// "' in plot " // trim(to_str(pl % id)))
       end select
 
       ! Get the number of <col_spec> nodes and get a list of them
@@ -2994,9 +2868,8 @@ contains
       if (n_cols /= 0) then
 
         if (pl % type == PLOT_TYPE_VOXEL) then
-          message = "Color specifications ignored in voxel plot " // &
-                     trim(to_str(pl % id))
-          if (master) call warning(message)
+          if (master) call warning("Color specifications ignored in voxel &
+               &plot " // trim(to_str(pl % id)))
         end if
 
         do j = 1, n_cols
@@ -3006,18 +2879,16 @@ contains
 
           ! Check and make sure 3 values are specified for RGB
           if (get_arraysize_double(node_col, "rgb") /= 3) then
-            message = "Bad RGB " &
-                 // "in plot " // trim(to_str(pl % id))
-            call fatal_error(message)
+            call fatal_error("Bad RGB in plot " &
+                 &// trim(to_str(pl % id)))
           end if
 
           ! Ensure that there is an id for this color specification
           if (check_for_node(node_col, "id")) then
             call get_node_value(node_col, "id", col_id)
           else
-            message = "Must specify id for color specification in plot " // &
-                      trim(to_str(pl % id))
-            call fatal_error(message)
+            call fatal_error("Must specify id for color specification in &
+                 &plot " // trim(to_str(pl % id)))
           end if
 
           ! Add RGB
@@ -3027,9 +2898,8 @@ contains
               col_id = cell_dict % get_key(col_id)
               call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
             else
-              message = "Could not find cell " // trim(to_str(col_id)) // &
-                   " specified in plot " // trim(to_str(pl % id))
-              call fatal_error(message)
+              call fatal_error("Could not find cell " // trim(to_str(col_id)) &
+                   &// " specified in plot " // trim(to_str(pl % id)))
             end if
 
           else if (pl % color_by == PLOT_COLOR_MATS) then
@@ -3038,9 +2908,9 @@ contains
               col_id = material_dict % get_key(col_id)
               call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
             else
-              message = "Could not find material " // trim(to_str(col_id)) // &
-                   " specified in plot " // trim(to_str(pl % id))
-              call fatal_error(message)
+              call fatal_error("Could not find material " &
+                   &// trim(to_str(col_id)) // " specified in plot " &
+                   &// trim(to_str(pl % id)))
             end if
 
           end if
@@ -3053,9 +2923,8 @@ contains
       if (n_meshlines /= 0) then
 
         if (pl % type == PLOT_TYPE_VOXEL) then
-          message = "Meshlines ignored in voxel plot " // & 
-                     trim(to_str(pl % id))
-          call warning(message)
+          call warning("Meshlines ignored in voxel plot " &
+               &// trim(to_str(pl % id)))
         end if
         
         select case(n_meshlines)
@@ -3070,9 +2939,8 @@ contains
             if (check_for_node(node_meshlines, "meshtype")) then
               call get_node_value(node_meshlines, "meshtype", meshtype)
             else
-              message = "Must specify a meshtype for meshlines " // &
-                        "specification in plot " // trim(to_str(pl % id))
-              call fatal_error(message)
+              call fatal_error("Must specify a meshtype for meshlines &
+                   &specification in plot " // trim(to_str(pl % id)))
             end if
             
             ! Ensure that there is a linewidth for this meshlines specification
@@ -3080,9 +2948,8 @@ contains
               call get_node_value(node_meshlines, "linewidth", &
                   pl % meshlines_width)
             else
-              message = "Must specify a linewidth for meshlines " // &
-                        "specification in plot " // trim(to_str(pl % id))
-              call fatal_error(message)
+              call fatal_error("Must specify a linewidth for meshlines &
+                   &specification in plot " // trim(to_str(pl % id)))
             end if
 
             ! Check for color
@@ -3090,9 +2957,8 @@ contains
               
               ! Check and make sure 3 values are specified for RGB
               if (get_arraysize_double(node_meshlines, "color") /= 3) then
-                message = "Bad RGB for meshlines color " // &
-                          "in plot " // trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("Bad RGB for meshlines color in plot " &
+                     &// trim(to_str(pl % id)))
               end if
               
               call get_node_array(node_meshlines, "color", &
@@ -3108,19 +2974,17 @@ contains
             case ('ufs')
 
               if (.not. associated(ufs_mesh)) then
-                message = "No UFS mesh for meshlines on plot " // &
-                          trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("No UFS mesh for meshlines on plot " &
+                     &// trim(to_str(pl % id)))
               end if
-              
+ 
               pl % meshlines_mesh => ufs_mesh
 
             case ('cmfd')
 
               if (.not. cmfd_run) then
-                message = "Need CMFD run to plot CMFD mesh for meshlines " // &
-                          "on plot " // trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("Need CMFD run to plot CMFD mesh for &
+                     &meshlines on plot " // trim(to_str(pl % id)))
               end if
 
               i_mesh = cmfd_tallies(1) % &
@@ -3129,19 +2993,17 @@ contains
               pl % meshlines_mesh => meshes(i_mesh)
 
             case ('entropy')
-            
+ 
               if (.not. associated(entropy_mesh)) then
-                message = "No entropy mesh for meshlines on plot " // &
-                          trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("No entropy mesh for meshlines on plot " &
+                     &// trim(to_str(pl % id)))
               end if
-              
+ 
               if (.not. allocated(entropy_mesh % dimension)) then
-                message = "No dimension specified on entropy mesh for " // &
-                          "meshlines on plot " // trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("No dimension specified on entropy mesh &
+                     &for meshlines on plot " // trim(to_str(pl % id)))
               end if
-              
+ 
               pl % meshlines_mesh => entropy_mesh
 
             case ('tally')
@@ -3150,37 +3012,31 @@ contains
               if (check_for_node(node_meshlines, "id")) then
                 call get_node_value(node_meshlines, "id", meshid)
               else
-                message = "Must specify a mesh id for meshlines tally mesh" // &
-                          "specification in plot " // trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("Must specify a mesh id for meshlines tally &
+                     &mesh specification in plot " // trim(to_str(pl % id)))
               end if
 
               ! Check if the specified tally mesh exists
               if (mesh_dict % has_key(meshid)) then
                 pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid))
                 if (meshes(meshid) % type /= LATTICE_RECT) then
-                  message = "Non-rectangular mesh specified in meshlines " // &
-                            "for plot " // trim(to_str(pl % id))
-                  call fatal_error(message)
+                  call fatal_error("Non-rectangular mesh specified in &
+                       &meshlines for plot " // trim(to_str(pl % id)))
                 end if
               else
-                message = "Could not find mesh " // &
-                          trim(to_str(meshid)) // &
-                          " specified in meshlines for plot " // &
-                          trim(to_str(pl % id))
-                call fatal_error(message)
+                call fatal_error("Could not find mesh " &
+                     &// trim(to_str(meshid)) // " specified in meshlines for &
+                     &plot " // trim(to_str(pl % id)))
               end if
 
             case default
-              message = "Invalid type for meshlines on plot " // &
-                        trim(to_str(pl % id)) // ": " // trim(meshtype)
-              call fatal_error(message)
+              call fatal_error("Invalid type for meshlines on plot " &
+                    &// trim(to_str(pl % id)) // ": " // trim(meshtype))
             end select
 
           case default
-            message = "Mutliple meshlines" // &
-                 " specified in plot " // trim(to_str(pl % id))
-            call fatal_error(message)
+            call fatal_error("Mutliple meshlines specified in plot " &
+                 &// trim(to_str(pl % id)))
         end select
         
       end if
@@ -3191,16 +3047,14 @@ contains
       if (n_masks /= 0) then
 
         if (pl % type == PLOT_TYPE_VOXEL) then
-          message = "Mask ignored in voxel plot " // &
-                     trim(to_str(pl % id))
-          if (master) call warning(message)
+          if (master) call warning("Mask ignored in voxel plot " &
+               &// trim(to_str(pl % id)))
         end if
 
         select case(n_masks)
           case default
-            message = "Mutliple masks" // &
-                 " specified in plot " // trim(to_str(pl % id))
-            call fatal_error(message)
+            call fatal_error("Mutliple masks specified in plot " &
+                 &// trim(to_str(pl % id)))
           case (1)
 
             ! Get pointer to mask
@@ -3210,9 +3064,8 @@ contains
             n_comp = 0
             n_comp = get_arraysize_integer(node_mask, "components")
             if (n_comp == 0) then
-              message = "Missing <components> in mask of plot " // &
-                        trim(to_str(pl % id))
-              call fatal_error(message)
+              call fatal_error("Missing <components> in mask of plot " &
+                   &// trim(to_str(pl % id)))
             end if
             allocate(iarray(n_comp))
             call get_node_array(node_mask, "components", iarray)
@@ -3227,9 +3080,9 @@ contains
                 if (cell_dict % has_key(col_id)) then
                   iarray(j) = cell_dict % get_key(col_id)
                 else
-                  message = "Could not find cell " // trim(to_str(col_id)) // &
-                       " specified in the mask in plot " // trim(to_str(pl % id))
-                  call fatal_error(message)
+                  call fatal_error("Could not find cell " &
+                       &// trim(to_str(col_id)) // " specified in the mask in &
+                       &plot " // trim(to_str(pl % id)))
                 end if
 
               else if (pl % color_by == PLOT_COLOR_MATS) then
@@ -3237,9 +3090,9 @@ contains
                 if (material_dict % has_key(col_id)) then
                   iarray(j) = material_dict % get_key(col_id)
                 else
-                  message = "Could not find material " // trim(to_str(col_id)) // &
-                       " specified in the mask in plot " // trim(to_str(pl % id))
-                  call fatal_error(message)
+                  call fatal_error("Could not find material " &
+                       &// trim(to_str(col_id)) // " specified in the mask in &
+                       &plot " // trim(to_str(pl % id)))
                 end if
 
               end if
@@ -3251,9 +3104,8 @@ contains
                 if (check_for_node(node_mask, "background")) then
                   call get_node_array(node_mask, "background", pl % colors(j) % rgb)
                 else
-                  message = "Missing <background> in mask of plot " // &
-                            trim(to_str(pl % id))
-                  call fatal_error(message)
+                  call fatal_error("Missing <background> in mask of plot " &
+                       &// trim(to_str(pl % id)))
                 end if
               end if
             end do
@@ -3297,13 +3149,11 @@ contains
     inquire(FILE=path_cross_sections, EXIST=file_exists)
     if (.not. file_exists) then
        ! Could not find cross_sections.xml file
-       message = "Cross sections XML file '" // trim(path_cross_sections) // &
-            "' does not exist!"
-       call fatal_error(message)
+       call fatal_error("Cross sections XML file '" &
+            &// trim(path_cross_sections) // "' does not exist!")
     end if
 
-    message = "Reading cross sections XML file..."
-    call write_message(message, 5)
+    call write_message("Reading cross sections XML file...", 5)
 
     ! Parse cross_sections.xml file
     call open_xmldoc(doc, path_cross_sections)
@@ -3329,8 +3179,8 @@ contains
     elseif (len_trim(temp_str) == 0) then
        filetype = ASCII
     else
-       message = "Unknown filetype in cross_sections.xml: " // trim(temp_str)
-       call fatal_error(message)
+       call fatal_error("Unknown filetype in cross_sections.xml: " &
+            &// trim(temp_str))
     end if
 
     ! copy default record length and entries for binary files
@@ -3345,8 +3195,8 @@ contains
 
     ! Allocate xs_listings array
     if (n_listings == 0) then
-       message = "No ACE table listings present in cross_sections.xml file!"
-       call fatal_error(message)
+       call fatal_error("No ACE table listings present in cross_sections.xml &
+            &file!")
     else
        allocate(xs_listings(n_listings))
     end if
@@ -3404,8 +3254,7 @@ contains
        if (check_for_node(node_ace, "path")) then
          call get_node_value(node_ace, "path", temp_str)
        else
-         message = "Path missing for isotope " // listing % name
-         call fatal_error(message)
+         call fatal_error("Path missing for isotope " // listing % name)
        end if
 
        if (starts_with(temp_str, '/')) then
@@ -3428,9 +3277,9 @@ contains
     ! Check that 0K nuclides are listed in the cross_sections.xml file
     do i = 1, n_res_scatterers_total
       if (.not. xs_listing_dict % has_key(trim(nuclides_0K(i) % name_0K))) then
-        message = "Could not find nuclide " // trim(nuclides_0K(i) % name_0K) // &
-          " in cross_sections.xml file!"
-        call fatal_error(message)
+        call fatal_error("Could not find nuclide " &
+             &// trim(nuclides_0K(i) % name_0K) &
+             &// " in cross_sections.xml file!")
       end if
     end do
 
@@ -4269,8 +4118,7 @@ contains
       call list_density % append(density * 0.992742_8)
 
     case default
-      message = "Cannot expand element: " // name
-      call fatal_error(message)
+      call fatal_error("Cannot expand element: " // name)
 
     end select
 
diff --git a/src/output.F90 b/src/output.F90
index 5dcd668b44..4d0ef7edd4 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -1555,7 +1555,6 @@ contains
 
   subroutine print_results()
 
-    character(2*MAX_LINE_LEN) :: message
     real(8) :: alpha   ! significance level for CI
     real(8) :: t_value ! t-value for confidence intervals
 
@@ -1589,8 +1588,8 @@ contains
       write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, &
            global_tallies(LEAKAGE) % sum_sq
     else
-      message = "Could not compute uncertainties -- only one active batch simulated!"
-      if (master) call warning(message)
+      if (master) call warning("Could not compute uncertainties -- only one &
+           &active batch simulated!")
 
       write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum
       write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH)  % sum
diff --git a/src/particle_restart.F90 b/src/particle_restart.F90
index dd68d5a896..6ff5d63c2d 100644
--- a/src/particle_restart.F90
+++ b/src/particle_restart.F90
@@ -16,7 +16,6 @@ module particle_restart
   private
   public ::  run_particle_restart
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
   type(BinaryOutput)        :: pr      ! Binary file
 
 contains
@@ -73,9 +72,8 @@ contains
     type(Particle), intent(inout) :: p
 
     ! Write meessage
-    message = "Loading particle restart file " // trim(path_particle_restart) &
-              // "..."
-    call write_message(message, 1)
+    call write_message("Loading particle restart file " &
+         &// trim(path_particle_restart) // "...", 1)
 
     ! Open file
     call pr % file_open(path_particle_restart, 'r')
diff --git a/src/physics.F90 b/src/physics.F90
index 353b7fd344..8fa0575cdc 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -20,8 +20,6 @@ module physics
 
   implicit none
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 ! TODO: Figure out how to write particle restart files in sample_angle,
 ! sample_energy, etc.
 
@@ -49,17 +47,15 @@ contains
 
     ! Display information about collision
     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."
-      call write_message(message)
+      call write_message("    " // trim(reaction_name(p % event_MT)) &
+           &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) &
+           &// ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV.")
     end if
 
     ! check for very low energy
     if (p % E < 1.0e-100_8) then
       p % alive = .false.
-      message = "Killing neutron with extremely low energy"
-      if (master) call warning(message)
+      if (master) call warning("Killing neutron with extremely low energy")
     end if
 
   end subroutine collision
@@ -161,8 +157,7 @@ contains
       ! Check to make sure that a nuclide was sampled
       if (i > mat % n_nuclides) then
         call write_particle_restart(p)
-        message = "Did not sample any nuclide during collision."
-        call fatal_error(message)
+        call fatal_error("Did not sample any nuclide during collision.")
       end if
 
       ! Find atom density
@@ -370,9 +365,8 @@ contains
         ! Check to make sure inelastic scattering reaction sampled
         if (i > nuc % n_reaction) then
           call write_particle_restart(p)
-          message = "Did not sample any reaction for nuclide " // &
-               trim(nuc % name)
-          call fatal_error(message)
+          call fatal_error("Did not sample any reaction for nuclide " &
+               &// trim(nuc % name))
         end if
 
         rxn => nuc % reactions(i)
@@ -723,8 +717,8 @@ contains
         mu = sab % inelastic_data(l) % mu(k, j)
 
       else
-        message = "Invalid secondary energy mode on S(a,b) table " // &
-             trim(sab % name)
+        call fatal_error("Invalid secondary energy mode on S(a,b) table " &
+             &// trim(sab % name))
       end if  ! (inelastic secondary energy treatment)
     end if  ! (elastic or inelastic)
 
@@ -976,8 +970,7 @@ contains
       end do
 
     case default
-      message = "Not a recognized resonance scattering treatment!"
-      call fatal_error(message)
+      call fatal_error("Not a recognized resonance scattering treatment!")
     end select
     
   end subroutine sample_target_velocity
@@ -1095,8 +1088,7 @@ contains
       call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh)
       if (.not. in_mesh) then
         call write_particle_restart(p)
-        message = "Source site outside UFS mesh!"
-        call fatal_error(message)
+        call fatal_error("Source site outside UFS mesh!")
       end if
 
       if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
@@ -1121,10 +1113,9 @@ contains
 
     ! Check for fission bank size getting hit
     if (n_bank + nu > size(fission_bank)) then
-      message = "Maximum number of sites in fission bank reached. This can &
-           &result in irreproducible results using different numbers of &
-           &processes/threads."
-      if (master) call warning(message)
+      if (master) call warning("Maximum number of sites in fission bank &
+           &reached. This can result in irreproducible results using different &
+           &numbers of processes/threads.")
     end if
 
     ! Bank source neutrons
@@ -1249,9 +1240,8 @@ contains
         n_sample = n_sample + 1
         if (n_sample == MAX_SAMPLE) then
           ! call write_particle_restart(p)
-          message = "Resampled energy distribution maximum number of " // &
-               "times for nuclide " // nuc % name
-          call fatal_error(message)
+          call fatal_error("Resampled energy distribution maximum number of " &
+               &// "times for nuclide " // nuc % name)
         end if
       end do
 
@@ -1276,9 +1266,8 @@ contains
         n_sample = n_sample + 1
         if (n_sample == MAX_SAMPLE) then
           ! call write_particle_restart(p)
-          message = "Resampled energy distribution maximum number of " // &
-               "times for nuclide " // nuc % name
-          call fatal_error(message)
+          call fatal_error("Resampled energy distribution maximum number of " &
+               &// "times for nuclide " // nuc % name)
         end if
       end do
 
@@ -1459,8 +1448,7 @@ contains
         end if
       else
         ! call write_particle_restart(p)
-        message = "Unknown interpolation type: " // trim(to_str(interp))
-        call fatal_error(message)
+        call fatal_error("Unknown interpolation type: " // trim(to_str(interp)))
       end if
 
       ! Because of floating-point roundoff, it may be possible for mu to be
@@ -1471,8 +1459,8 @@ contains
 
     else
       ! call write_particle_restart(p)
-      message = "Unknown angular distribution type: " // trim(to_str(type))
-      call fatal_error(message)
+      call fatal_error("Unknown angular distribution type: " &
+           &// trim(to_str(type)))
     end if
 
   end function sample_angle
@@ -1619,9 +1607,8 @@ contains
       NET = int(edist % data(3 + 2*NR + NE))
       if (NR > 0) then
         ! call write_particle_restart(p)
-        message = "Multiple interpolation regions not supported while &
-             &attempting to sample equiprobable energy bins."
-        call fatal_error(message)
+        call fatal_error("Multiple interpolation regions not supported while &
+             &attempting to sample equiprobable energy bins.")
       end if
 
       ! determine index on incoming energy grid and interpolation factor
@@ -1684,14 +1671,12 @@ contains
       NR  = int(edist % data(1))
       NE  = int(edist % data(2 + 2*NR))
       if (NR == 1) then
-        message = "Assuming linear-linear interpolation when sampling &
-             &continuous tabular distribution"
-        if (master) call warning(message)
+        if (master) call warning("Assuming linear-linear interpolation when &
+             &sampling continuous tabular distribution")
       else if (NR > 1) then
         ! call write_particle_restart(p)
-        message = "Multiple interpolation regions not supported while &
-             &attempting to sample continuous tabular distribution."
-        call fatal_error(message)
+        call fatal_error("Multiple interpolation regions not supported while &
+             &attempting to sample continuous tabular distribution.")
       end if
 
       ! find energy bin and calculate interpolation factor -- if the energy is
@@ -1749,9 +1734,8 @@ contains
       if (ND > 0) then
         ! discrete lines present
         ! call write_particle_restart(p)
-        message = "Discrete lines in continuous tabular distributed not &
-             &yet supported"
-        call fatal_error(message)
+        call fatal_error("Discrete lines in continuous tabular distributed not &
+             &yet supported")
       end if
 
       ! determine outgoing energy bin
@@ -1791,8 +1775,7 @@ contains
         end if
       else
         ! call write_particle_restart(p)
-        message = "Unknown interpolation type: " // trim(to_str(INTT))
-        call fatal_error(message)
+        call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
       end if
 
       ! Now interpolate between incident energy bins i and i + 1
@@ -1833,8 +1816,7 @@ contains
         n_sample = n_sample + 1
         if (n_sample == MAX_SAMPLE) then
           ! call write_particle_restart(p)
-          message = "Too many rejections on Maxwell fission spectrum."
-          call fatal_error(message)
+          call fatal_error("Too many rejections on Maxwell fission spectrum.")
         end if
       end do
 
@@ -1866,8 +1848,7 @@ contains
         n_sample = n_sample + 1
         if (n_sample == MAX_SAMPLE) then
           ! call write_particle_restart(p)
-          message = "Too many rejections on evaporation spectrum."
-          call fatal_error(message)
+          call fatal_error("Too many rejections on evaporation spectrum.")
         end if
       end do
 
@@ -1908,8 +1889,7 @@ contains
         n_sample = n_sample + 1
         if (n_sample == MAX_SAMPLE) then
           ! call write_particle_restart(p)
-          message = "Too many rejections on Watt spectrum."
-          call fatal_error(message)
+          call fatal_error("Too many rejections on Watt spectrum.")
         end if
       end do
 
@@ -1919,8 +1899,7 @@ contains
 
       if (.not. present(mu_out)) then
         ! call write_particle_restart(p)
-        message = "Law 44 called without giving mu_out as argument."
-        call fatal_error(message)
+        call fatal_error("Law 44 called without giving mu_out as argument.")
       end if
 
       ! read number of interpolation regions and incoming energies
@@ -1928,9 +1907,8 @@ contains
       NE = int(edist % data(2 + 2*NR))
       if (NR > 0) then
         ! call write_particle_restart(p)
-        message = "Multiple interpolation regions not supported while &
-             &attempting to sample Kalbach-Mann distribution."
-        call fatal_error(message)
+        call fatal_error("Multiple interpolation regions not supported while &
+             &attempting to sample Kalbach-Mann distribution.")
       end if
 
       ! find energy bin and calculate interpolation factor -- if the energy is
@@ -1989,9 +1967,8 @@ contains
       if (ND > 0) then
         ! discrete lines present
         ! call write_particle_restart(p)
-        message = "Discrete lines in continuous tabular distributed not &
-             &yet supported"
-        call fatal_error(message)
+        call fatal_error("Discrete lines in continuous tabular distributed not &
+             &yet supported")
       end if
 
       ! determine outgoing energy bin
@@ -2045,8 +2022,7 @@ contains
         KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
       else
         ! call write_particle_restart()
-        message = "Unknown interpolation type: " // trim(to_str(INTT))
-        call fatal_error(message)
+        call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
       end if
 
       ! Now interpolate between incident energy bins i and i + 1
@@ -2072,8 +2048,7 @@ contains
 
       if (.not. present(mu_out)) then
         ! call write_particle_restart()
-        message = "Law 61 called without giving mu_out as argument."
-        call fatal_error(message)
+        call fatal_error("Law 61 called without giving mu_out as argument.")
       end if
 
       ! read number of interpolation regions and incoming energies
@@ -2081,9 +2056,8 @@ contains
       NE = int(edist % data(2 + 2*NR))
       if (NR > 0) then
         ! call write_particle_restart()
-        message = "Multiple interpolation regions not supported while &
-             &attempting to sample correlated energy-angle distribution."
-        call fatal_error(message)
+        call fatal_error("Multiple interpolation regions not supported while &
+             &attempting to sample correlated energy-angle distribution.")
       end if
 
       ! find energy bin and calculate interpolation factor -- if the energy is
@@ -2142,9 +2116,8 @@ contains
       if (ND > 0) then
         ! discrete lines present
         ! call write_particle_restart()
-        message = "Discrete lines in continuous tabular distributed not &
-             &yet supported"
-        call fatal_error(message)
+        call fatal_error("Discrete lines in continuous tabular distributed not &
+             &yet supported")
       end if
 
       ! determine outgoing energy bin
@@ -2185,8 +2158,7 @@ contains
         end if
       else
         ! call write_particle_restart()
-        message = "Unknown interpolation type: " // trim(to_str(INTT))
-        call fatal_error(message)
+        call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
       end if
 
       ! Now interpolate between incident energy bins i and i + 1
@@ -2249,8 +2221,7 @@ contains
         end if
       else
         ! call write_particle_restart()
-        message = "Unknown interpolation type: " // trim(to_str(JJ))
-        call fatal_error(message)
+        call fatal_error("Unknown interpolation type: " // trim(to_str(JJ)))
       end if
 
     case (66)
diff --git a/src/search.F90 b/src/search.F90
index acd1c9ccbc..c8fc71d365 100644
--- a/src/search.F90
+++ b/src/search.F90
@@ -59,8 +59,8 @@ contains
       ! check for large number of iterations
       n_iteration = n_iteration + 1
       if (n_iteration == MAX_ITERATION) then
-        call fatal_error("Reached maximum number of iterations on binary " &
-            &// "search.")
+        call fatal_error("Reached maximum number of iterations on binary &
+            &search.")
       end if
     end do
 
@@ -111,8 +111,8 @@ contains
       ! check for large number of iterations
       n_iteration = n_iteration + 1
       if (n_iteration == MAX_ITERATION) then
-        call fatal_error("Reached maximum number of iterations on binary " &
-            &// "search.")
+        call fatal_error("Reached maximum number of iterations on binary &
+            &search.")
       end if
     end do
 
@@ -163,8 +163,8 @@ contains
       ! check for large number of iterations
       n_iteration = n_iteration + 1
       if (n_iteration == MAX_ITERATION) then
-        call fatal_error("Reached maximum number of iterations on binary " &
-            &// "search.")
+        call fatal_error("Reached maximum number of iterations on binary &
+            &search.")
       end if
     end do
 
diff --git a/src/source.F90 b/src/source.F90
index 8b8542fd6a..07c1bfb411 100644
--- a/src/source.F90
+++ b/src/source.F90
@@ -19,8 +19,6 @@ module source
 
   implicit none
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 contains
 
 !===============================================================================
@@ -36,15 +34,14 @@ contains
     type(Bank), pointer :: src => null() ! source bank site
     type(BinaryOutput) :: sp ! statepoint/source binary file
 
-    message = "Initializing source particles..."
-    call write_message(message, 6)
+    call write_message("Initializing source particles...", 6)
 
     if (path_source /= '') then
       ! Read the source from a binary file instead of sampling from some
       ! assumed source distribution
 
-      message = 'Reading source file from ' // trim(path_source) // '...'
-      call write_message(message, 6)
+      call write_message('Reading source file from ' // trim(path_source) &
+           &// '...', 6)
 
       ! Open the binary file
       call sp % file_open(path_source, 'r', serial = .false.)
@@ -54,8 +51,8 @@ contains
 
       ! Check to make sure this is a source file
       if (itmp /= FILETYPE_SOURCE) then
-        message = "Specified starting source file not a source file type."
-        call fatal_error(message)
+        call fatal_error("Specified starting source file not a source file &
+             &type.")
       end if
 
       ! Read in the source bank
@@ -81,8 +78,7 @@ contains
 
     ! Write out initial source
     if (write_initial_source) then
-      message = 'Writing out initial source guess...'
-      call write_message(message, 1)
+      call write_message('Writing out initial source guess...', 1)
 #ifdef HDF5
       filename = trim(path_output) // 'initial_source.h5'
 #else
@@ -144,9 +140,8 @@ contains
         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(message)
+            call fatal_error("Maximum number of external source spatial &
+                 &resamples reached!")
           end if
         end if
       end do
@@ -174,9 +169,8 @@ contains
         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(message)
+            call fatal_error("Maximum number of external source spatial &
+                 &resamples reached!")
           end if
           cycle
         end if
@@ -209,8 +203,7 @@ contains
       site % uvw = external_source % params_angle
 
     case default
-      message = "No angle distribution specified for external source!"
-      call fatal_error(message)
+      call fatal_error("No angle distribution specified for external source!")
     end select
 
     ! Sample energy distribution
@@ -241,8 +234,7 @@ contains
       end do
 
     case default
-      message = "No energy distribution specified for external source!"
-      call fatal_error(message)
+      call fatal_error("No energy distribution specified for external source!")
     end select
 
     ! Set the random number generator back to the tracking stream.
diff --git a/src/state_point.F90 b/src/state_point.F90
index c3332ed740..5d8b370cd2 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -26,7 +26,6 @@ module state_point
 
   implicit none
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
   type(BinaryOutput)        :: sp      ! Statepoint/source output file
 
 contains
@@ -55,8 +54,7 @@ contains
 #endif
 
     ! Write message
-    message = "Creating state point " // trim(filename) // "..."
-    call write_message(message, 1)
+    call write_message("Creating state point " // trim(filename) // "...", 1)
 
     if (master) then
       ! Create statepoint file
@@ -316,8 +314,8 @@ contains
 #endif
 
         ! Write message for new file creation
-        message = "Creating source file " // trim(filename) // "..."
-        call write_message(message, 1)
+        call write_message("Creating source file " // trim(filename) // "...", &
+             &1)
 
         ! Create separate source file
         call sp % file_create(filename, serial = .false.)
@@ -361,8 +359,7 @@ contains
 #endif
 
       ! Write message for new file creation
-      message = "Creating source file " // trim(filename) // "..."
-      call write_message(message, 1)
+      call write_message("Creating source file " // trim(filename) // "...", 1)
 
       ! Always create this file because it will be overwritten
       call sp % file_create(filename, serial = .false.)
@@ -532,8 +529,8 @@ contains
     type(TallyObject), pointer :: t => null()
 
     ! Write message
-    message = "Loading state point " // trim(path_state_point) // "..."
-    call write_message(message, 1)
+    call write_message("Loading state point " // trim(path_state_point) &
+         &// "...", 1)
 
     ! Open file for reading
     call sp % file_open(path_state_point, 'r', serial = .false.)
@@ -545,9 +542,8 @@ contains
     ! current version
     call sp % read_data(int_array(1), "revision")
     if (int_array(1) /= REVISION_STATEPOINT) then
-      message = "State point version does not match current version " &
-                // "in OpenMC."
-      call fatal_error(message)
+      call fatal_error("State point version does not match current version &
+           &in OpenMC.")
     end if
 
     ! Read OpenMC version
@@ -556,9 +552,8 @@ contains
     call sp % read_data(int_array(3), "version_release")
     if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR &
         .or. int_array(3) /= VERSION_RELEASE) then
-      message = "State point file was created with a different version " &
-                // "of OpenMC."
-      if (master) call warning(message)
+      if (master) call warning("State point file was created with a different &
+           &version of OpenMC.")
     end if
 
     ! Read date and time
@@ -668,8 +663,8 @@ contains
       ! Check size of tally results array
       if (int_array(1) /= t % total_score_bins .and. &
           int_array(2) /= t % total_filter_bins) then
-        message = "Input file tally structure is different from restart."
-        call fatal_error(message)
+        call fatal_error("Input file tally structure is different from &
+             &restart.")
       end if
 
       ! Read number of filters
@@ -742,8 +737,8 @@ contains
 
     ! Check to make sure source bank is present
     if (path_source_point == path_state_point .and. .not. source_present) then
-      message = "Source bank must be contained in statepoint restart file"
-      call fatal_error(message)
+      call fatal_error("Source bank must be contained in statepoint restart &
+           &file")
     end if
 
     ! Read tallies to master
@@ -755,8 +750,8 @@ contains
       ! Read number of global tallies
       call sp % read_data(int_array(1), "n_global_tallies", collect=.false.)
       if (int_array(1) /= N_GLOBAL_TALLIES) then
-        message = "Number of global tallies does not match in state point."
-        call fatal_error(message)
+        call fatal_error("Number of global tallies does not match in state &
+             &point.")
       end if
 
       ! Read global tally data
@@ -792,8 +787,8 @@ contains
         call sp % file_close()
 
         ! Write message
-        message = "Loading source file " // trim(path_source_point) // "..."
-        call write_message(message, 1)
+        call write_message("Loading source file " // trim(path_source_point) &
+             &// "...", 1)
 
         ! Open source file
         call sp % file_open(path_source_point, 'r', serial = .false.)
diff --git a/src/string.F90 b/src/string.F90
index 7b3689f0eb..cff69bb16e 100644
--- a/src/string.F90
+++ b/src/string.F90
@@ -10,8 +10,6 @@ module string
      module procedure int4_to_str, int8_to_str, real_to_str
   end interface
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 contains
 
 !===============================================================================
@@ -51,9 +49,8 @@ contains
         if (i_end > 0) then
           n = n + 1
           if (i_end - i_start + 1 > len(words(n))) then
-            message = "The word '" // string(i_start:i_end) // &
-                 "' is longer than the space allocated for it."
-            if (master) call warning(message)
+            if (master) call warning("The word '" // string(i_start:i_end) &
+                &// "' is longer than the space allocated for it.")
           end if
           words(n) = string(i_start:i_end)
           ! reset indices
@@ -212,8 +209,8 @@ function zero_padded(num, n_digits) result(str)
   ! Make sure n_digits is reasonable. 10 digits is the maximum needed for the
   ! largest integer(4).
   if (n_digits > 10) then
-    message = 'zero_padded called with an unreasonably large n_digits (>10)'
-    call fatal_error(message)
+    call fatal_error('zero_padded called with an unreasonably large &
+         &n_digits (>10)')
   end if
 
   ! Write a format string of the form '(In.m)' where n is the max width and
diff --git a/src/tally.F90 b/src/tally.F90
index db3894379b..87f5a2cab7 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -21,7 +21,6 @@ module tally
 
   implicit none
 
-  character(2*MAX_LINE_LEN) :: message        ! Message to output unit
   integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
 !$omp threadprivate(position)
 
@@ -871,8 +870,8 @@ contains
                   end do REACTION_LOOP
 
                 else
-                  message = "Invalid score type on tally " // to_str(t % id) // "."
-                  call fatal_error(message)
+                  call fatal_error("Invalid score type on tally " &
+                       &// to_str(t % id) // ".")
                 end if
               end select
 
@@ -1009,8 +1008,8 @@ contains
                   end do
 
                 else
-                  message = "Invalid score type on tally " // to_str(t % id) // "."
-                  call fatal_error(message)
+                  call fatal_error("Invalid score type on tally " &
+                       &// to_str(t % id) // ".")
                 end if
               end select
             end if
@@ -1209,8 +1208,8 @@ contains
             end do REACTION_LOOP
 
           else
-            message = "Invalid score type on tally " // to_str(t % id) // "."
-            call fatal_error(message)
+            call fatal_error("Invalid score type on tally " &
+                 &// to_str(t % id) // ".")
           end if
         end select
 
@@ -1361,8 +1360,8 @@ contains
           end do
 
         else
-          message = "Invalid score type on tally " // to_str(t % id) // "."
-          call fatal_error(message)
+          call fatal_error("Invalid score type on tally " &
+               &// to_str(t % id) // ".")
         end if
       end select
 
@@ -1707,9 +1706,8 @@ contains
                 case (SCORE_EVENTS)
                   score = ONE
                 case default
-                  message = "Invalid score type on tally " // &
-                       to_str(t % id) // "."
-                  call fatal_error(message)
+                  call fatal_error("Invalid score type on tally " &
+                       &// to_str(t % id) // ".")
                 end select
 
               else
@@ -1783,9 +1781,8 @@ contains
                 case (SCORE_EVENTS)
                   score = ONE
                 case default
-                  message = "Invalid score type on tally " // &
-                       to_str(t % id) // "."
-                  call fatal_error(message)
+                  call fatal_error("Invalid score type on tally " &
+                       &// to_str(t % id) // ".")
                 end select
               end if
 
@@ -2223,8 +2220,7 @@ contains
           ! Check for errors
           if (filter_index <= 0 .or. filter_index > &
                t % total_filter_bins) then
-            message = "Score index outside range."
-            call fatal_error(message)
+            call fatal_error("Score index outside range.")
           end if
 
           ! Add to surface current tally
@@ -2563,18 +2559,16 @@ contains
 
     ! check to see if any of the active tally lists has been allocated
     if (active_tallies % size() > 0) then
-      message = "Active tallies should not exist before CMFD tallies!"
-      call fatal_error(message)
+      call fatal_error("Active tallies should not exist before CMFD tallies!")
     else if (active_analog_tallies % size() > 0) then
-      message = 'Active analog tallies should not exist before CMFD tallies!'
-      call fatal_error(message)
+      call fatal_error('Active analog tallies should not exist before CMFD &
+           &tallies!')
     else if (active_tracklength_tallies % size() > 0) then
-      message = "Active tracklength tallies should not exist before CMFD &
-           &tallies!"
-      call fatal_error(message)
+      call fatal_error("Active tracklength tallies should not exist before &
+           &CMFD tallies!")
     else if (active_current_tallies % size() > 0) then
-      message = "Active current tallies should not exist before CMFD tallies!"
-      call fatal_error(message)
+      call fatal_error("Active current tallies should not exist before CMFD &
+           &tallies!")
     end if
 
     do i = 1, n_cmfd_tallies
diff --git a/src/tracking.F90 b/src/tracking.F90
index eb94016310..9545d41ee6 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -16,8 +16,6 @@ module tracking
   use track_output,    only: initialize_particle_track, write_particle_track, &
                              finalize_particle_track
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 contains
 
 !===============================================================================
@@ -41,8 +39,7 @@ contains
 
     ! Display message if high verbosity or trace is on
     if (verbosity >= 9 .or. trace) then
-      message = "Simulating Particle " // trim(to_str(p % id))
-      call write_message(message)
+      call write_message("Simulating Particle " // trim(to_str(p % id)))
     end if
 
     ! If the cell hasn't been determined based on the particle's location,
@@ -52,8 +49,7 @@ contains
 
       ! Particle couldn't be located
       if (.not. found_cell) then
-        message = "Could not locate particle " // trim(to_str(p % id))
-        call fatal_error(message)
+        call fatal_error("Could not locate particle " // trim(to_str(p % id)))
       end if
 
       ! set birth cell attribute
@@ -200,9 +196,8 @@ contains
       ! If particle has too many events, display warning and kill it
       n_event = n_event + 1
       if (n_event == MAX_EVENTS) then
-        message = "Particle " // trim(to_str(p%id)) // " underwent maximum &
-             &number of events."
-        if (master) call warning(message)
+        if (master) call warning("Particle " // trim(to_str(p%id)) &
+             &// " underwent maximum number of events.")
         p % alive = .false.
       end if
 
diff --git a/src/xml_interface.F90 b/src/xml_interface.F90
index 907de77d96..1a24bc14b3 100644
--- a/src/xml_interface.F90
+++ b/src/xml_interface.F90
@@ -37,8 +37,6 @@ module xml_interface
     module procedure get_node_array_string
   end interface get_node_array
 
-  character(2*MAX_LINE_LEN) :: message ! Message to output unit
-
 contains
 
 !===============================================================================
@@ -202,9 +200,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
 
     ! Extract value
@@ -235,9 +232,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
       
     ! Extract value
@@ -268,9 +264,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
       
     ! Extract value
@@ -301,9 +296,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
       
     ! Extract value
@@ -334,9 +328,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
       
     ! Extract value
@@ -367,9 +360,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
       
     ! Extract value
@@ -400,9 +392,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then 
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " // &
+                getNodeName(ptr) // ".")
     end if
  
     ! Extract value
@@ -433,9 +424,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
 
     ! Get the size
@@ -462,9 +452,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " &
+           &// getNodeName(ptr) // ".")
     end if
 
     ! Get the size
@@ -491,9 +480,8 @@ contains
 
     ! Leave if it was not found
     if (.not. found) then
-      message = "Node " // node_name // " not part of Node " // &
-                getNodeName(ptr) // "."
-      call fatal_error(message)
+      call fatal_error("Node " // node_name // " not part of Node " // &
+                getNodeName(ptr) // ".")
     end if
 
     ! Get the size