diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp
index ff494666..9a57c9b0 100644
--- a/examples/crack_branching.cpp
+++ b/examples/crack_branching.cpp
@@ -38,6 +38,7 @@ int main( int argc, char* argv[] )
         double t_final = 43e-6;
         double dt = 5e-8;
         double output_frequency = 5;
+        bool output_reference = true;
 
         // Material constants
         double E = 72e+9;                      // [Pa]
@@ -68,7 +69,7 @@ int main( int argc, char* argv[] )
             CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
         model_type force_model( delta, K, G0 );
         CabanaPD::Inputs inputs( num_cell, low_corner, high_corner, t_final, dt,
-                                 output_frequency );
+                                 output_frequency, output_reference );
         inputs.read_args( argc, argv );
 
         // Create particles from mesh.
diff --git a/examples/elastic_wave.cpp b/examples/elastic_wave.cpp
index 2c194c18..2f4aa850 100644
--- a/examples/elastic_wave.cpp
+++ b/examples/elastic_wave.cpp
@@ -35,6 +35,7 @@ int main( int argc, char* argv[] )
         double t_final = 0.6;
         double dt = 0.01;
         int output_frequency = 5;
+        bool output_reference = true;
         double K = 1.0;
         double G = 0.5;
         double delta = 0.075;
@@ -51,7 +52,7 @@ int main( int argc, char* argv[] )
         model_type force_model( delta, K, G );
 
         CabanaPD::Inputs inputs( num_cell, low_corner, high_corner, t_final, dt,
-                                 output_frequency );
+                                 output_frequency, output_reference );
         inputs.read_args( argc, argv );
 
         // Create particles from mesh.
diff --git a/examples/kalthoff_winkler.cpp b/examples/kalthoff_winkler.cpp
index d315f03b..ce109657 100644
--- a/examples/kalthoff_winkler.cpp
+++ b/examples/kalthoff_winkler.cpp
@@ -46,6 +46,7 @@ int main( int argc, char* argv[] )
         double t_final = 70e-6;
         double dt = 0.133e-6;
         int output_frequency = 10;
+        bool output_reference = true;
 
         // Material constants
         double E = 191e+9;                           // [Pa]
@@ -82,7 +83,7 @@ int main( int argc, char* argv[] )
         //     CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
         // model_type force_model( delta, K, G, G0 );
         CabanaPD::Inputs inputs( num_cell, low_corner, high_corner, t_final, dt,
-                                 output_frequency );
+                                 output_frequency, output_reference );
         inputs.read_args( argc, argv );
 
         // Create particles from mesh.
diff --git a/src/CabanaPD_Input.cpp b/src/CabanaPD_Input.cpp
index 9a01e2c6..ade16b15 100644
--- a/src/CabanaPD_Input.cpp
+++ b/src/CabanaPD_Input.cpp
@@ -70,13 +70,14 @@ namespace CabanaPD
 // FIXME: hardcoded.
 Inputs::Inputs( const std::array<int, 3> nc, std::array<double, 3> lc,
                 std::array<double, 3> hc, const double t_f, const double dt,
-                const int of )
+                const int of, const bool ref )
     : num_cells( nc )
     , low_corner( lc )
     , high_corner( hc )
     , final_time( t_f )
     , timestep( dt )
     , output_frequency( of )
+    , output_reference( ref )
 {
     num_steps = final_time / timestep;
 }
diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp
index aaaa022b..6c7b47b1 100644
--- a/src/CabanaPD_Input.hpp
+++ b/src/CabanaPD_Input.hpp
@@ -32,12 +32,13 @@ class Inputs
     double final_time;
     double timestep;
     int output_frequency;
+    bool output_reference;
 
     bool half_neigh = false;
 
     Inputs( const std::array<int, 3> nc, std::array<double, 3> lc,
             std::array<double, 3> hc, const double t_f, const double dt,
-            const int output_freq );
+            const int output_freq, const bool output_ref );
     ~Inputs();
     void read_args( int argc, char* argv[] );
 };
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index dd38941b..b08d4e49 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -124,6 +124,7 @@ class SolverElastic
 
         num_steps = inputs->num_steps;
         output_frequency = inputs->output_frequency;
+        output_reference = inputs->output_reference;
 
         // Create integrator.
         integrator = std::make_shared<integrator_type>( inputs->timestep );
@@ -236,7 +237,7 @@ class SolverElastic
 
                 step_output( step, W );
                 particles->output( step / output_frequency,
-                                   step * inputs->timestep );
+                                   step * inputs->timestep, output_reference );
             }
             other_time += other_timer.seconds();
         }
@@ -303,6 +304,7 @@ class SolverElastic
 
     int num_steps;
     int output_frequency;
+    bool output_reference;
 
   protected:
     std::shared_ptr<particle_type> particles;
@@ -445,7 +447,7 @@ class SolverFracture
 
                 this->step_output( step, W );
                 particles->output( step / output_frequency,
-                                   step * inputs->timestep );
+                                   step * inputs->timestep, output_reference );
             }
             other_time += other_timer.seconds();
         }
@@ -456,6 +458,7 @@ class SolverFracture
 
     using base_type::num_steps;
     using base_type::output_frequency;
+    using base_type::output_reference;
 
   protected:
     using base_type::comm;