diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h
index 0e426ffb1..9adb825a3 100644
--- a/include/cantera/oneD/Domain1D.h
+++ b/include/cantera/oneD/Domain1D.h
@@ -8,13 +8,11 @@
 #ifndef CT_DOMAIN1D_H
 #define CT_DOMAIN1D_H
 
-
 #include "cantera/base/xml.h"
 #include "cantera/base/stringUtils.h"
 #include "cantera/base/ctexceptions.h"
 #include "refine.h"
 
-
 namespace Cantera
 {
 
@@ -32,18 +30,17 @@ const int cPorousType       = 109;
 class MultiJac;
 class OneDim;
 
-
 /**
  * Base class for one-dimensional domains.
  */
 class Domain1D
 {
 public:
-
     /**
      * Constructor.
      * @param nv      Number of variables at each grid point.
      * @param points  Number of grid points.
+     * @param time    (unused)
      */
     Domain1D(size_t nv=1, size_t points=1,
              doublereal time = 0.0) :
@@ -62,38 +59,31 @@ public:
         resize(nv, points);
     }
 
-    /// Destructor. Does nothing
     virtual ~Domain1D() {
         delete m_refiner;
     }
 
-    /// Domain type flag.
+    //! Domain type flag.
     int domainType() {
         return m_type;
     }
 
-    /**
-     * The left-to-right location of this domain.
-     */
+    //! The left-to-right location of this domain.
     size_t domainIndex() {
         return m_index;
     }
 
-    /**
-     * True if the domain is a connector domain.
-     */
+    //! True if the domain is a connector domain.
     bool isConnector() {
         return (m_type >= cConnectorType);
     }
 
-    /**
-     * The container holding this domain.
-     */
+    //! The container holding this domain.
     const OneDim& container() const {
         return *m_container;
     }
 
-    /**
+    /*!
      * Specify the container object for this domain, and the
      * position of this domain in the list.
      */
@@ -102,35 +92,31 @@ public:
         m_index = index;
     }
 
-    /*
-     * Set the Jacobian bandwidth. See the discussion of method bandwidth.
-     */
+    //! Set the Jacobian bandwidth. See the discussion of method bandwidth().
     void setBandwidth(int bw = -1) {
         m_bw = bw;
     }
 
+    //! Set the Jacobian bandwidth for this domain.
     /**
-     * Set the Jacobian bandwidth for this domain.  When class
-     * OneDim computes the bandwidth of the overall multi-domain
-     * problem (in OneDim::resize()), it calls this method for the
-     * bandwidth of each domain. If setBandwidth has not been
-     * called, then a negative bandwidth is returned, in which
-     * case OneDim assumes that this domain is dense -- that is,
-     * at each point, all components depend on the value of all
-     * other components at that point. In this case, the bandwidth
-     * is bw = 2*nComponents() - 1. However, if this domain
-     * contains some components that are uncoupled from other
-     * components at the same point, then this default bandwidth
-     * may greatly overestimate the true bandwidth, with a
-     * substantial penalty in performance. For such domains, use
-     * method setBandwidth to specify the bandwidth before passing
-     * this domain to the Sim1D or OneDim constructor.
+     * When class OneDim computes the bandwidth of the overall multi-domain
+     * problem (in OneDim::resize()), it calls this method for the bandwidth
+     * of each domain. If setBandwidth has not been called, then a negative
+     * bandwidth is returned, in which case OneDim assumes that this domain is
+     * dense -- that is, at each point, all components depend on the value of
+     * all other components at that point. In this case, the bandwidth is bw =
+     * 2*nComponents() - 1. However, if this domain contains some components
+     * that are uncoupled from other components at the same point, then this
+     * default bandwidth may greatly overestimate the true bandwidth, with a
+     * substantial penalty in performance. For such domains, use method
+     * setBandwidth to specify the bandwidth before passing this domain to the
+     * Sim1D or OneDim constructor.
      */
     size_t bandwidth() {
         return m_bw;
     }
 
-    /**
+    /*!
      * Initialize. This method is called by OneDim::init() for
      * each domain once at the beginning of a simulation. Base
      * class method does nothing, but may be overloaded.
@@ -140,7 +126,7 @@ public:
     virtual void setInitialState(doublereal* xlocal = 0) {}
     virtual void setState(size_t point, const doublereal* state, doublereal* x) {}
 
-    /**
+    /*!
      * Resize the domain to have nv components and np grid points.
      * This method is virtual so that subclasses can perform other
      * actions required to resize the domain.
@@ -168,12 +154,12 @@ public:
         locate();
     }
 
-    /// Return a reference to the grid refiner.
+    //! Return a reference to the grid refiner.
     Refiner& refiner() {
         return *m_refiner;
     }
 
-    /// Number of components at each grid point.
+    //! Number of components at each grid point.
     size_t nComponents() const {
         return m_nv;
     }
@@ -195,7 +181,7 @@ public:
         }
     }
 
-    /// Number of grid points in this domain.
+    //! Number of grid points in this domain.
     size_t nPoints() const {
         return m_points;
     }
@@ -217,7 +203,7 @@ public:
         }
     }
 
-    /// Name of the nth component. May be overloaded.
+    //! Name of the nth component. May be overloaded.
     virtual std::string componentName(size_t n) const {
         if (m_name[n] != "") {
             return m_name[n];
@@ -236,7 +222,7 @@ public:
         }
     }
 
-    /// index of component with name \a name.
+    //! index of component with name \a name.
     size_t componentIndex(const std::string& name) const {
         size_t nc = nComponents();
         for (size_t n = 0; n < nc; n++) {
@@ -248,9 +234,7 @@ public:
                            "no component named "+name);
     }
 
-    /**
-     * Set the lower and upper bounds for each solution component.
-     */
+    //! Set the lower and upper bounds for each solution component.
     void setBounds(size_t nl, const doublereal* lower,
                    size_t nu, const doublereal* upper) {
         if (nl < m_nv || nu < m_nv)
@@ -266,71 +250,69 @@ public:
         m_max[n] = upper;
     }
 
-    /// set the error tolerances for all solution components.
+    //! set the error tolerances for all solution components.
     void setTolerances(size_t nr, const doublereal* rtol,
                        size_t na, const doublereal* atol, int ts = 0);
 
-    /// set the error tolerances for solution component \a n.
+    //! set the error tolerances for solution component \a n.
     void setTolerances(size_t n, doublereal rtol, doublereal atol, int ts = 0);
 
-    /// set scalar error tolerances. All solution components will
-    /// have the same relative and absolute error tolerances.
+    //! set scalar error tolerances. All solution components will
+    //! have the same relative and absolute error tolerances.
     void setTolerances(doublereal rtol, doublereal atol,int ts=0);
 
     void setTolerancesTS(doublereal rtol, doublereal atol);
 
     void setTolerancesSS(doublereal rtol, doublereal atol);
 
-    /// Relative tolerance of the nth component.
+    //! Relative tolerance of the nth component.
     doublereal rtol(size_t n) {
         return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]);
     }
 
-    /// Absolute tolerance of the nth component.
+    //! Absolute tolerance of the nth component.
     doublereal atol(size_t n) {
         return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]);
     }
 
-    /// Upper bound on the nth component.
+    //! Upper bound on the nth component.
     doublereal upperBound(size_t n) const {
         return m_max[n];
     }
 
-    /// Lower bound on the nth component
+    //! Lower bound on the nth component
     doublereal lowerBound(size_t n) const {
         return m_min[n];
     }
 
-
-    /**
-     * Prepare to do time stepping with time step dt. Copy the
-     * internally-stored solution at the last time step to array
-     * x0.
+    /*!
+     * Prepare to do time stepping with time step dt. Copy the internally-
+     * stored solution at the last time step to array x0.
      */
     void initTimeInteg(doublereal dt, const doublereal* x0) {
         std::copy(x0 + loc(), x0 + loc() + size(), m_slast.begin());
         m_rdt = 1.0/dt;
     }
 
-    /**
-     * Prepare to solve the steady-state problem.
-     * Set the internally-stored reciprocal of the time step to 0,0
+    /*!
+     * Prepare to solve the steady-state problem. Set the internally-stored
+     * reciprocal of the time step to 0,0
      */
     void setSteadyMode() {
         m_rdt = 0.0;
     }
 
-    /// True if in steady-state mode
+    //! True if in steady-state mode
     bool steady() {
         return (m_rdt == 0.0);
     }
 
-    /// True if not in steady-state mode
+    //! True if not in steady-state mode
     bool transient() {
         return (m_rdt != 0.0);
     }
 
-    /**
+    /*!
      * Set this if something has changed in the governing
      * equations (e.g. the value of a constant has been changed,
      * so that the last-computed Jacobian is no longer valid.
@@ -338,7 +320,7 @@ public:
      */
     void needJacUpdate();
 
-    /**
+    /*!
      * Evaluate the steady-state residual at all points, even if in
      * transient mode. Used only to print diagnostic output.
      */
@@ -349,9 +331,13 @@ public:
     //! Evaluate the residual function at point j. If j == npos,
     //! evaluate the residual function at all points.
     /*!
-     *   @param j   Grid point j
-     *   @param x   Soln vector. This is the input.
-     *   @param r   residual this is the output.
+     *  @param j  Grid point at which to update the residual
+     *  @param[in] x  State vector
+     *  @param[out] r  residual vector
+     *  @param[out] mask  Boolean mask indicating whether each solution
+     *      component has a time derivative (1) or not (0).
+     *  @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-
+     *  state.)
      */
     virtual void eval(size_t j, doublereal* x, doublereal* r,
                       integer* mask, doublereal rdt=0.0);
@@ -367,9 +353,6 @@ public:
         m_td[n] = 0;
     }
 
-    /**
-     * Does nothing.
-     */
     virtual void update(doublereal* x) {}
 
     doublereal time() const {
@@ -407,7 +390,7 @@ public:
      * classes should call the base class method in addition to restoring
      * their own data.
      *
-     * @param o XML_Node for this domain
+     * @param dom XML_Node for this domain
      * @param soln Current value of the solution vector, local to this object.
      * @param loglevel 0 to suppress all output; 1 to show warnings; 2 for
      *      verbose output
@@ -468,54 +451,41 @@ public:
     }
 
     /**
-     * Set the left neighbor to domain 'left.' Method 'locate' is
-     * called to update the global positions of this domain and
-     * all those to its right.
+     * Set the left neighbor to domain 'left.' Method 'locate' is called to
+     * update the global positions of this domain and all those to its right.
      */
     void linkLeft(Domain1D* left) {
         m_left = left;
         locate();
     }
 
-    /**
-     * Set the right neighbor to domain 'right.'
-     */
+    //! Set the right neighbor to domain 'right.'
     void linkRight(Domain1D* right) {
         m_right = right;
     }
 
-    /**
-     * Append domain 'right' to this one, and update all links.
-     */
+    //! Append domain 'right' to this one, and update all links.
     void append(Domain1D* right) {
         linkRight(right);
         right->linkLeft(this);
     }
 
-    /**
-     * Return a pointer to the left neighbor.
-     */
+    //! Return a pointer to the left neighbor.
     Domain1D* left() const {
         return m_left;
     }
 
-    /**
-     * Return a pointer to the right neighbor.
-     */
+    //! Return a pointer to the right neighbor.
     Domain1D* right() const {
         return m_right;
     }
 
-    /**
-     * Value of component n at point j in the previous solution.
-     */
+    //! Value of component n at point j in the previous solution.
     double prevSoln(size_t n, size_t j) const {
         return m_slast[m_nv*j + n];
     }
 
-    /**
-     * Specify an identifying tag for this domain.
-     */
+    //! Specify an identifying tag for this domain.
     void setID(const std::string& s) {
         m_id = s;
     }
@@ -528,9 +498,7 @@ public:
         }
     }
 
-    /**
-     * Specify descriptive text for this domain.
-     */
+    //! Specify descriptive text for this domain.
     void setDesc(const std::string& s) {
         m_desc = s;
     }
@@ -541,6 +509,8 @@ public:
     virtual void getTransientMask(integer* mask) {}
 
     virtual void showSolution_s(std::ostream& s, const doublereal* x) {}
+
+    //! Print the solution.
     virtual void showSolution(const doublereal* x);
 
     doublereal z(size_t jlocal) const {
@@ -553,7 +523,6 @@ public:
         return m_z[m_points - 1];
     }
 
-
     void setProfile(const std::string& name, doublereal* values, doublereal* soln) {
         for (size_t n = 0; n < m_nv; n++) {
             if (name == componentName(n)) {
@@ -577,6 +546,7 @@ public:
         return m_z[point];
     }
 
+    //! called to set up initial grid, and after grid refinement
     virtual void setupGrid(size_t n, const doublereal* z);
 
     void setGrid(size_t n, const doublereal* z);
@@ -591,9 +561,7 @@ public:
      */
     virtual void _getInitialSoln(doublereal* x);
 
-    /**
-     * Initial value of solution component \a n at grid point \a j.
-     */
+    //! Initial value of solution component \a n at grid point \a j.
     virtual doublereal initialValue(size_t n, size_t j);
 
     /**
@@ -614,7 +582,6 @@ public:
     bool m_adiabatic;
 
 protected:
-
     doublereal m_rdt;
     size_t m_nv;
     size_t m_points;
@@ -629,11 +596,10 @@ protected:
     size_t m_index;
     int m_type;
 
-    //! Starting location within the solution vector for unknowns
-    //! that correspond to this domain
+    //! Starting location within the solution vector for unknowns that
+    //! correspond to this domain
     /*!
-     * Remember there may be multiple domains associated with
-     * this problem
+     * Remember there may be multiple domains associated with this problem
      */
     size_t m_iloc;
 
@@ -648,9 +614,6 @@ protected:
     vector_int m_td;
     std::vector<std::string> m_name;
     int m_bw;
-
-private:
-
 };
 }
 
diff --git a/include/cantera/oneD/Inlet1D.h b/include/cantera/oneD/Inlet1D.h
index 372312f37..8e0b8cc82 100644
--- a/include/cantera/oneD/Inlet1D.h
+++ b/include/cantera/oneD/Inlet1D.h
@@ -38,10 +38,8 @@ const int RightInlet = -1;
 class Bdry1D : public Domain1D
 {
 public:
-
     Bdry1D();
 
-    /// Initialize.
     virtual void init() {
         _init(1);
     }
@@ -93,7 +91,6 @@ public:
     virtual void setupGrid(size_t n, const doublereal* z) {}
 
 protected:
-
     void _init(size_t n);
 
     StFlow* m_flow_left, *m_flow_right;
@@ -119,14 +116,7 @@ private:
  */
 class Inlet1D : public Bdry1D
 {
-
 public:
-
-    /**
-     * Constructor. Create a new Inlet1D instance. If invoked
-     * without parameters, a left inlet (facing right) is
-     * constructed).
-     */
     Inlet1D() : Bdry1D(), m_V0(0.0), m_nsp(0), m_flow(0) {
         m_type = cInletType;
         m_xstr = "";
@@ -143,7 +133,6 @@ public:
         return m_V0;
     }
 
-
     virtual void showSolution(const doublereal* x) {
         char buf[80];
         sprintf(buf, "    Mass Flux:   %10.4g kg/m^2/s \n", m_mdot);
@@ -187,7 +176,6 @@ public:
     virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
 
 protected:
-
     int m_ilr;
     doublereal m_V0;
     size_t m_nsp;
@@ -196,15 +184,12 @@ protected:
     StFlow* m_flow;
 };
 
-
 /**
  * A terminator that does nothing.
  */
 class Empty1D : public Domain1D
 {
-
 public:
-
     Empty1D() : Domain1D() {
         m_type = cEmptyType;
     }
@@ -223,9 +208,6 @@ public:
     virtual void _getInitialSoln(doublereal* x) {
         x[0] = 0.0;
     }
-
-protected:
-
 };
 
 /**
@@ -234,7 +216,6 @@ protected:
  */
 class Symm1D : public Bdry1D
 {
-
 public:
 
     Symm1D() : Bdry1D() {
@@ -256,19 +237,15 @@ public:
     virtual void _getInitialSoln(doublereal* x) {
         x[0] = m_temp;
     }
-
-protected:
-
 };
 
 
 /**
+ *  An outlet.
  */
 class Outlet1D : public Bdry1D
 {
-
 public:
-
     Outlet1D() : Bdry1D() {
         m_type = cOutletType;
     }
@@ -288,23 +265,15 @@ public:
     virtual void _getInitialSoln(doublereal* x) {
         x[0] = m_temp;
     }
-protected:
-
 };
 
 
-
 /**
  * An outlet with specified composition.
  */
 class OutletRes1D : public Bdry1D
 {
-
 public:
-
-    /**
-     * Constructor.
-     */
     OutletRes1D() : Bdry1D(), m_nsp(0), m_flow(0) {
         m_type = cOutletResType;
         m_xstr = "";
@@ -337,14 +306,12 @@ public:
     virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
 
 protected:
-
     size_t m_nsp;
     vector_fp m_yres;
     std::string m_xstr;
     StFlow* m_flow;
 };
 
-
 /**
  * A non-reacting surface. The axial velocity is zero
  * (impermeable), as is the transverse velocity (no slip). The
@@ -353,9 +320,7 @@ protected:
  */
 class Surf1D : public Bdry1D
 {
-
 public:
-
     Surf1D() : Bdry1D() {
         m_type = cSurfType;
     }
@@ -389,21 +354,14 @@ public:
         writelog(buf);
         writelog("\n");
     }
-
-protected:
-
 };
 
-
 /**
  * A reacting surface.
- *
  */
 class ReactingSurf1D : public Bdry1D
 {
-
 public:
-
     ReactingSurf1D() : Bdry1D(),
         m_kin(0), m_surfindex(0), m_nsp(0) {
         m_type = cSurfType;
@@ -455,7 +413,6 @@ public:
     }
 
 protected:
-
     InterfaceKinetics* m_kin;
     SurfPhase* m_sphase;
     size_t m_surfindex, m_nsp;
diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h
index 94bc8ea70..5f68cb27f 100644
--- a/include/cantera/oneD/MultiJac.h
+++ b/include/cantera/oneD/MultiJac.h
@@ -19,17 +19,12 @@ namespace Cantera
 /**
  * Class MultiJac evaluates the Jacobian of a system of equations
  * defined by a residual function supplied by an instance of class
- * 'OneDim.' The residual function may consist of several linked
+ * OneDim. The residual function may consist of several linked
  * 1D domains, with different variables in each domain.
  */
 class MultiJac : public BandMatrix
 {
-
 public:
-
-    /**
-     * Constructor.
-     */
     MultiJac(OneDim& r);
 
     /**
@@ -40,36 +35,30 @@ public:
      */
     void eval(doublereal* x0, doublereal* resid0, double rdt);
 
-    /**
-     * Elapsed CPU time spent computing the Jacobian.
-     */
+    //! Elapsed CPU time spent computing the Jacobian.
     doublereal elapsedTime() const {
         return m_elapsed;
     }
 
-    /// Number of Jacobian evaluations.
+    //! Number of Jacobian evaluations.
     int nEvals() const {
         return m_nevals;
     }
 
-    /**
-     * Number of times 'incrementAge' has been called since the
-     * last evaluation
-     */
+    //! Number of times 'incrementAge' has been called since the last
+    //! evaluation
     int age() const {
         return m_age;
     }
 
-    /**
-     * Increment the Jacobian age.
-     */
+    //! Increment the Jacobian age.
     void incrementAge() {
         m_age++;
     }
 
     void updateTransient(doublereal rdt, integer* mask);
 
-    /// Set the age.
+    //! Set the Jacobian age.
     void setAge(int age) {
         m_age = age;
     }
@@ -81,11 +70,10 @@ public:
     void incrementDiagonal(int j, doublereal d);
 
 protected:
-
     //!  Residual evaluator for this jacobian
     /*!
-     *  This is a pointer to the residual evaluator. This
-     *  object isn't owned by this jacobian object.
+     *  This is a pointer to the residual evaluator. This object isn't owned
+     *  by this jacobian object.
      */
     OneDim* m_resid;
 
@@ -102,5 +90,3 @@ protected:
 }
 
 #endif
-
-
diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h
index 9d06911cf..8ebfd2496 100644
--- a/include/cantera/oneD/MultiNewton.h
+++ b/include/cantera/oneD/MultiNewton.h
@@ -20,9 +20,7 @@ namespace Cantera
  */
 class MultiNewton
 {
-
 public:
-
     MultiNewton(int sz);
     virtual ~MultiNewton();
 
@@ -30,21 +28,40 @@ public:
         return m_n;
     }
 
-    /// Compute undamped step
+    //! Compute the undamped Newton step.  The residual function is evaluated
+    //! at `x`, but the Jacobian is not recomputed.
     void step(doublereal* x, doublereal* step,
               OneDim& r, MultiJac& jac, int loglevel);
 
-    /// Compute factor to keep all components in bounds.
+    /**
+     * Return the factor by which the undamped Newton step 'step0'
+     * must be multiplied in order to keep all solution components in
+     * all domains between their specified lower and upper bounds.
+     */
     doublereal boundStep(const doublereal* x0, const doublereal* step0,
                          const OneDim& r, int loglevel);
 
+    /**
+     * On entry, step0 must contain an undamped Newton step for the
+     * solution x0. This method attempts to find a damping coefficient
+     * such that the next undamped step would have a norm smaller than
+     * that of step0. If successful, the new solution after taking the
+     * damped step is returned in x1, and the undamped step at x1 is
+     * returned in step1.
+     */
     int dampStep(const doublereal* x0, const doublereal* step0,
                  doublereal* x1, doublereal* step1, doublereal& s1,
                  OneDim& r, MultiJac& jac, int loglevel, bool writetitle);
 
+    //! Compute the weighted 2-norm of `step`.
     doublereal norm2(const doublereal* x, const doublereal* step,
                      OneDim& r) const;
 
+    /**
+     * Find the solution to F(X) = 0 by damped Newton iteration.  On
+     * entry, x0 contains an initial estimate of the solution.  On
+     * successful return, x1 contains the converged solution.
+     */
     int solve(doublereal* x0, doublereal* x1, OneDim& r, MultiJac& jac,
               int loglevel);
 
@@ -56,22 +73,22 @@ public:
     /// Change the problem size.
     void resize(size_t points);
 
-
 protected:
-
+    //! Get a pointer to an array of length m_n for temporary work space.
     doublereal* getWorkArray();
+
+    //! Release a work array by pushing its pointer onto the stack of
+    //! available arrays.
     void releaseWorkArray(doublereal* work);
+
     std::vector<doublereal*> m_workarrays;
     int m_maxAge;
     size_t m_nv, m_np, m_n;
     doublereal m_elapsed;
 
 private:
-
     char m_buf[100];
 };
 }
 
 #endif
-
-
diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h
index 220b13cba..5e191fcbb 100644
--- a/include/cantera/oneD/OneDim.h
+++ b/include/cantera/oneD/OneDim.h
@@ -20,19 +20,14 @@ class Func1;
  */
 class OneDim
 {
-
 public:
-
-    // Default constructor.
     OneDim();
 
-    // Constructor.
+    //! Construct a OneDim container for the domains in the list *domains*.
     OneDim(std::vector<Domain1D*> domains);
-
-    /// Destructor.
     virtual ~OneDim();
 
-    /// Add a domain.
+    /// Add a domain. Domains are added left-to-right.
     void addDomain(Domain1D* d);
 
     //! Return a reference to the Jacobian evaluator.
@@ -116,7 +111,11 @@ public:
         return m_bw;
     }
 
-    /// Initialize.
+    /*!
+     * Initialize all domains. On the first call, this methods calls the init
+     * method of each domain, proceeding from left to right. Subsequent calls
+     * do nothing.
+     */
     void init();
 
     /// Total number of points.
@@ -125,8 +124,9 @@ public:
     }
 
     /**
-     * Steady-state max norm of the residual evaluated using solution x.
-     * On return, array r contains the steady-state residual values.
+     * Steady-state max norm (infinity norm) of the residual evaluated using
+     * solution x. On return, array r contains the steady-state residual
+     * values. Used only for diagnostic output.
      */
     doublereal ssnorm(doublereal* x, doublereal* r);
 
@@ -135,7 +135,7 @@ public:
         return m_rdt;
     }
 
-    /// Prepare for time stepping beginning with solution x.
+    //! Prepare for time stepping beginning with solution *x* and timestep *dt*.
     void initTimeInteg(doublereal dt, doublereal* x);
 
     /// True if transient mode.
@@ -148,14 +148,14 @@ public:
         return (m_rdt == 0.0);
     }
 
-
-    /**
-     * Set steady mode.  After invoking this method, subsequent
-     * calls to solve() will solve the steady-state problem.
+    /*!
+     * Prepare to solve the steady-state problem. After invoking this method,
+     * subsequent calls to solve() will solve the steady-state problem. Sets
+     * the reciprocal of the time step to zero, and, if it was previously non-
+     * zero, signals that a new Jacobian will be needed.
      */
     void setSteadyMode();
 
-
     /**
      * Evaluate the multi-domain residual function
      *
@@ -170,18 +170,31 @@ public:
     void eval(size_t j, double* x, double* r, doublereal rdt=-1.0,
               int count = 1);
 
-    /// Pointer to the domain global point i belongs to.
+    //! Return a pointer to the domain global point *i* belongs to.
+    /*!
+     * The domains are scanned right-to-left, and the first one with starting
+     * location less or equal to i is returned.
+     */
     Domain1D* pointDomain(size_t i);
 
+    //! Call after one or more grids has been refined.
     void resize();
 
-    //doublereal solveTime() { return m_solve_time; }
-
     //void setTransientMask();
     vector_int& transientMask() {
         return m_mask;
     }
 
+    /*!
+     * Take time steps using Backward Euler.
+     *
+     * @param nsteps number of steps
+     * @param dt initial step size
+     * @param x current solution vector
+     * @param r solution vector after time stepping
+     * @param loglevel controls amount of printed diagnostics
+     * @returns size of last timestep taken
+     */
     double timeStep(int nsteps, double dt, double* x,
                     double* r, int loglevel);
 
@@ -214,6 +227,18 @@ public:
             m_ts_jac_age = m_ss_jac_age;
         }
     }
+
+    /**
+     * Save statistics on function and Jacobian evaluation, and reset the
+     * counters. Statistics are saved only if the number of Jacobian
+     * evaluations is greater than zero. The statistics saved are:
+     *
+     *    - number of grid points
+     *    - number of Jacobian evaluations
+     *    - CPU time spent evaluating Jacobians
+     *    - number of non-Jacobian function evaluations
+     *    - CPU time spent evaluating functions
+     */
     void saveStats();
 
     //! Set a function that will be called every time #eval is called.
@@ -224,7 +249,6 @@ public:
     }
 
 protected:
-
     void evalSSJacobian(doublereal* x, doublereal* xnew);
 
     doublereal m_tmin;        // minimum timestep size
@@ -259,7 +283,6 @@ protected:
     Func1* m_interrupt;
 
 private:
-
     // statistics
     int m_nevals;
     doublereal m_evaltime;
@@ -268,12 +291,8 @@ private:
     vector_fp m_jacElapsed;
     vector_int m_funcEvals;
     vector_fp m_funcElapsed;
-
-
 };
 
 }
 
 #endif
-
-
diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h
index 7d2e8ea8e..c162e9937 100644
--- a/include/cantera/oneD/Sim1D.h
+++ b/include/cantera/oneD/Sim1D.h
@@ -12,25 +12,21 @@ namespace Cantera
 {
 
 /**
- * One-dimensional simulations. Class Sim1D extends class OneDim
- * by storing the solution vector, and by adding a hybrid
- * Newton/time-stepping solver.
+ * One-dimensional simulations. Class Sim1D extends class OneDim by storing
+ * the solution vector, and by adding a hybrid Newton/time-stepping solver.
  */
 class Sim1D : public OneDim
 {
-
 public:
 
-
     //! Default constructor.
     /*!
-     *  This constructor is provided to make
-     *  the class default-constructible, but is not meant to be
-     *  used in most applications.  Use the next constructor
+     *  This constructor is provided to make the class default-constructible,
+     *  but is not meant to be used in most applications.  Use the next
+     *  constructor
      */
     Sim1D();
 
-
     /**
      * Standard constructor.
      * @param domains A vector of pointers to the domains to be linked together.
@@ -52,15 +48,41 @@ public:
     void setInitialGuess(const std::string& component, vector_fp& locs,
                          vector_fp& vals);
 
-    /// Set one entry in the solution vector.
+    /**
+     * Set a single value in the solution vector.
+     * @param dom domain number, beginning with 0 for the leftmost domain.
+     * @param comp component number
+     * @param localPoint grid point within the domain, beginning with 0 for
+     *     the leftmost grid point in the domain.
+     * @param value the value.
+     */
     void setValue(size_t dom, size_t comp, size_t localPoint,  doublereal value);
 
-    /// Get one entry in the solution vector.
+    /**
+     * Get one entry in the solution vector.
+     * @param dom domain number, beginning with 0 for the leftmost domain.
+     * @param comp component number
+     * @param localPoint grid point within the domain, beginning with 0 for
+     *     the leftmost grid point in the domain.
+     */
     doublereal value(size_t dom, size_t comp, size_t localPoint) const;
 
     doublereal workValue(size_t dom, size_t comp, size_t localPoint) const;
 
-    /// Specify a profile for one component of one domain.
+    /**
+     * Specify a profile for one component of one domain.
+     * @param dom domain number, beginning with 0 for the leftmost domain.
+     * @param comp component number
+     * @param pos A vector of relative positions, beginning with 0.0 at the
+     *     left of the domain, and ending with 1.0 at the right of the domain.
+     * @param values A vector of values corresponding to the relative position
+     *     locations.
+     *
+     * Note that the vector pos and values can have lengths different than the
+     * number of grid points, but their lengths must be equal. The values at
+     * the grid points will be linearly interpolated based on the (pos,
+     * values) specification.
+     */
     void setProfile(size_t dom, size_t comp, const vector_fp& pos,
                     const vector_fp& values);
 
@@ -96,10 +118,16 @@ public:
     /// Refine the grid in all domains.
     int refine(int loglevel=0);
 
+    //! Add node for fixed temperature point of freely propagating flame
     int setFixedTemperature(doublereal t);
+
     void setAdiabaticFlame(void);
 
-    /// Set the criteria for grid refinement.
+    /**
+     * Set grid refinement criteria. If dom >= 0, then the settings
+     * apply only to the specified domain.  If dom < 0, the settings
+     * are applied to each domain.  @see Refiner::setCriteria.
+     */
     void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
                            doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
     void setMaxGridPoints(int dom = -1, int npoints = 300);
@@ -112,7 +140,9 @@ public:
     */
     void setGridMin(int dom, double gridmin);
 
+    //! Initialize the solution with a previously-saved solution.
     void restore(const std::string& fname, const std::string& id, int loglevel=2);
+
     void getInitialSoln();
 
     void setSolution(const doublereal* soln) {
@@ -128,17 +158,20 @@ public:
     void evalSSJacobian();
 
 protected:
+    //! the solution vector
+    vector_fp m_x;
 
-    vector_fp m_x;          // the solution vector
-    vector_fp m_xnew;       // a work array used to hold the residual
-    //      or the new solution
-    doublereal m_tstep;     // timestep
-    vector_int m_steps;     // array of number of steps to take before
-    //      re-attempting the steady-state solution
+    //! a work array used to hold the residual or the new solution
+    vector_fp m_xnew;
 
+    //! timestep
+    doublereal m_tstep;
+
+    //! array of number of steps to take before re-attempting the steady-state
+    //! solution
+    vector_int m_steps;
 
 private:
-
     /// Calls method _finalize in each domain.
     void finalize();
 
@@ -146,11 +179,7 @@ private:
      * @return 0 if successful, -1 on failure
      */
     int newtonSolve(int loglevel);
-
-
 };
 
 }
 #endif
-
-
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 5d6058f98..f902e504b 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -1,6 +1,5 @@
 /**
  * @file StFlow.h
- *
  */
 // Copyright 2001  California Institute of Technology
 
@@ -34,36 +33,26 @@ const int c_Mixav_Transport = 0;
 const int c_Multi_Transport = 1;
 const int c_Soret = 2;
 
-//-----------------------------------------------------------
-//  Class StFlow
-//-----------------------------------------------------------
-
-
 /**
- *  This class represents 1D flow domains that satisfy the
- *  one-dimensional similarity solution for chemically-reacting,
- *  axisymmetric, flows.
+ *  This class represents 1D flow domains that satisfy the one-dimensional
+ *  similarity solution for chemically-reacting, axisymmetric, flows.
  */
 class StFlow : public Domain1D
 {
-
 public:
-
     //--------------------------------
     // construction and destruction
     //--------------------------------
 
-    /// Constructor. Create a new flow domain.
-    /// @param ph Object representing the gas phase. This object
-    /// will be used to evaluate all thermodynamic, kinetic, and transport
-    /// properties.
-    /// @param nsp Number of species.
+    //! Create a new flow domain.
+    //! @param ph Object representing the gas phase. This object will be used
+    //!     to evaluate all thermodynamic, kinetic, and transport properties.
+    //! @param nsp Number of species.
+    //! @param points Initial number of grid points
     StFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1);
 
-    /**
-     * @name Problem Specification
-     */
-    //@{
+    //! @name Problem Specification
+    //! @{
 
     virtual void setupGrid(size_t n, const doublereal* z);
 
@@ -85,31 +74,31 @@ public:
         m_thermo = &th;
     }
 
-    /// Set the kinetics manager. The kinetics manager must
+    //! Set the kinetics manager. The kinetics manager must
     void setKinetics(Kinetics& kin) {
         m_kin = &kin;
     }
 
-    /// set the transport manager
+    //! set the transport manager
     void setTransport(Transport& trans, bool withSoret = false);
     void enableSoret(bool withSoret);
     bool withSoret() const {
         return m_do_soret;
     }
 
-    /// Set the pressure. Since the flow equations are for the limit of
-    /// small Mach number, the pressure is very nearly constant
-    /// throughout the flow.
+    //! Set the pressure. Since the flow equations are for the limit of
+    //! small Mach number, the pressure is very nearly constant
+    //! throughout the flow.
     void setPressure(doublereal p) {
         m_press = p;
     }
 
-    /// The current pressure [Pa].
+    //! The current pressure [Pa].
     doublereal pressure() const {
         return m_press;
     }
 
-    /// @todo remove? may be unused
+    //! @todo remove? may be unused
     virtual void setState(size_t point, const doublereal* state,
                           doublereal* x) {
         setTemperature(point, state[2]);
@@ -118,7 +107,7 @@ public:
         }
     }
 
-    /// Write the initial solution estimate into array x.
+    //! Write the initial solution estimate into array x.
     virtual void _getInitialSoln(doublereal* x) {
         for (size_t j = 0; j < m_points; j++) {
             T(x,j) = m_thermo->temperature();
@@ -128,64 +117,58 @@ public:
 
     virtual void _finalize(const doublereal* x);
 
-
-    /// Sometimes it is desired to carry out the simulation
-    /// using a specified temperature profile, rather than
-    /// computing it by solving the energy equation. This
-    /// method specifies this profile.
+    //! Sometimes it is desired to carry out the simulation using a specified
+    //! temperature profile, rather than computing it by solving the energy
+    //! equation. This method specifies this profile.
     void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) {
         m_zfix = zfixed;
         m_tfix = tfixed;
     }
 
-    /**
-     * Set the temperature fixed point at grid point j, and
-     * disable the energy equation so that the solution will be
-     * held to this value.
+    /*!
+     * Set the temperature fixed point at grid point j, and disable the energy
+     * equation so that the solution will be held to this value.
      */
     void setTemperature(size_t j, doublereal t) {
         m_fixedtemp[j] = t;
         m_do_energy[j] = false;
     }
 
-    /**
-     * Set the mass fraction fixed point for species k at grid
-     * point j, and disable the species equation so that the
-     * solution will be held to this value.
-     * note: in practice, the species are hardly ever held fixed.
+    /*!
+     * Set the mass fraction fixed point for species k at grid point j, and
+     * disable the species equation so that the solution will be held to this
+     * value. Note: in practice, the species are hardly ever held fixed.
      */
     void setMassFraction(size_t j, size_t k, doublereal y) {
         m_fixedy(k,j) = y;
         m_do_species[k] = true; // false;
     }
 
-
-    /// The fixed temperature value at point j.
+    //! The fixed temperature value at point j.
     doublereal T_fixed(size_t j) const {
         return m_fixedtemp[j];
     }
 
-
-    /// The fixed mass fraction value of species k at point j.
+    //! The fixed mass fraction value of species k at point j.
     doublereal Y_fixed(size_t k, size_t j) const {
         return m_fixedy(k,j);
     }
 
+    // @}
 
     virtual std::string componentName(size_t n) const;
 
     size_t componentIndex(const std::string& name) const;
 
-
+    //! Print the solution.
     virtual void showSolution(const doublereal* x);
 
     //! Save the current solution for this domain into an XML_Node
     /*!
-     *
      *  @param o    XML_Node to save the solution to.
-     *  @param sol  Current value of the solution vector.
-     *              The object will pick out which part of the solution
-     *              vector pertains to this object.
+     *  @param sol  Current value of the solution vector. The object will pick
+     *              out which part of the solution vector pertains to this
+     *              object.
      */
     virtual XML_Node& save(XML_Node& o, const doublereal* const sol);
 
@@ -256,13 +239,18 @@ public:
 
     void integrateChem(doublereal* x,doublereal dt);
 
+    //! Change the grid size. Called after grid refinement.
     void resize(size_t components, size_t points);
 
     virtual void setFixedPoint(int j0, doublereal t0) {}
 
-
     void setJac(MultiJac* jac);
+
+    //! Set the gas object state to be consistent with the solution at point j.
     void setGas(const doublereal* x, size_t j);
+
+    //! Set the gas state to be consistent with the solution at the midpoint
+    //! between j and j + 1.
     void setGasAtMidpoint(const doublereal* x, size_t j);
 
     doublereal density(size_t j) const {
@@ -276,6 +264,13 @@ public:
         m_dovisc = dovisc;
     }
 
+    /*!
+     *  Evaluate the residual function for axisymmetric stagnation flow. If
+     *  jpt is less than zero, the residual function is evaluated at all grid
+     *  points. If jpt >= 0, then the residual function is only evaluated at
+     *  grid points jpt-1, jpt, and jpt+1. This option is used to efficiently
+     *  evaluate the Jacobian numerically.
+     */
     virtual void eval(size_t j, doublereal* x, doublereal* r,
                       integer* mask, doublereal rdt);
 
@@ -289,7 +284,6 @@ public:
                                 integer* diag, doublereal rdt) = 0;
 
 protected:
-
     doublereal component(const doublereal* x, size_t i, size_t j) const {
         return x[index(i,j)];
     }
@@ -306,15 +300,15 @@ protected:
         return m_wdot(k,j);
     }
 
-    /// write the net production rates at point j into array m_wdot
+    //! Write the net production rates at point `j` into array `m_wdot`
     void getWdot(doublereal* x, size_t j) {
         setGas(x,j);
         m_kin->getNetProductionRates(&m_wdot(0,j));
     }
 
     /**
-     * update the thermodynamic properties from point
-     * j0 to point j1 (inclusive), based on solution x.
+     * Update the thermodynamic properties from point j0 to point j1
+     * (inclusive), based on solution x.
      */
     void updateThermo(const doublereal* x, size_t j0, size_t j1) {
         for (size_t j = j0; j <= j1; j++) {
@@ -325,7 +319,6 @@ protected:
         }
     }
 
-
     //--------------------------------
     // central-differenced derivatives
     //--------------------------------
@@ -338,10 +331,8 @@ protected:
     }
 
 
-    //--------------------------------
-    //      solution components
-    //--------------------------------
-
+    //! @name Solution components
+    //! @{
 
     doublereal T(const doublereal* x, size_t j) const {
         return x[index(c_offset_T, j)];
@@ -391,11 +382,11 @@ protected:
     doublereal flux(size_t k, size_t j) const {
         return m_flux(k, j);
     }
+    //! @}
 
-
-    // convective spatial derivatives. These use upwind
-    // differencing, assuming u(z) is negative
-
+    //! @name convective spatial derivatives.
+    //! These use upwind differencing, assuming u(z) is negative
+    //! @{
     doublereal dVdz(const doublereal* x, size_t j) const {
         size_t jloc = (u(x,j) > 0.0 ? j : j + 1);
         return (V(x,jloc) - V(x,jloc-1))/m_dz[jloc-1];
@@ -410,6 +401,7 @@ protected:
         size_t jloc = (u(x,j) > 0.0 ? j : j + 1);
         return (T(x,jloc) - T(x,jloc-1))/m_dz[jloc-1];
     }
+    //! @}
 
     doublereal shear(const doublereal* x, size_t j) const {
         doublereal c1 = m_visc[j-1]*(V(x,j) - V(x,j-1));
@@ -427,13 +419,11 @@ protected:
         return m*m_nsp*m_nsp + m_nsp*j + k;
     }
 
+    //! Update the diffusive mass fluxes.
     void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1);
 
-
     //---------------------------------------------------------
-    //
     //             member data
-    //
     //---------------------------------------------------------
 
     // inlet
@@ -500,14 +490,15 @@ protected:
     vector_fp m_tfix;
 
     bool m_dovisc;
+
+    //! Update the transport properties at grid points in the range from `j0`
+    //! to `j1`, based on solution `x`.
     void updateTransport(doublereal* x, size_t j0, size_t j1);
 
 private:
     vector_fp m_ybar;
-
 };
 
-
 /**
  * A class for axisymmetric stagnation flows.
  */
@@ -554,26 +545,14 @@ public:
     }
 };
 
-
-/*
-class OneDFlow : public StFlow {
-public:
-    OneDFlow(igthermo_t* ph = 0, int nsp = 1, int points = 1) :
-        StFlow(ph, nsp, points) {
-    }
-    virtual ~OneDFlow() {}
-    virtual void eval(int j, doublereal* x, doublereal* r,
-        integer* mask, doublereal rdt);
-    virtual std::string flowType() { return "OneDFlow"; }
-    doublereal mdot(doublereal* x, int j) {
-        return x[index(c_offset_L,j)];
-    }
-
-private:
-    void updateTransport(doublereal* x,int j0, int j1);
-};
-*/
-
+/**
+ * Import a previous solution to use as an initial estimate. The
+ * previous solution may have been computed using a different
+ * reaction mechanism. Species in the old and new mechanisms are
+ * matched by name, and any species in the new mechanism that were
+ * not in the old one are set to zero. The new solution is created
+ * with the same number of grid points as in the old solution.
+ */
 void importSolution(size_t points, doublereal* oldSoln, IdealGasPhase& oldmech,
                     size_t size_new, doublereal* newSoln, IdealGasPhase& newmech);
 
diff --git a/include/cantera/oneD/refine.h b/include/cantera/oneD/refine.h
index 0c22c1f69..3310dece2 100644
--- a/include/cantera/oneD/refine.h
+++ b/include/cantera/oneD/refine.h
@@ -8,12 +8,25 @@ namespace Cantera
 
 class Domain1D;
 
+//! Refine Domain1D grids so that profiles satisfy adaptation tolerances
 class Refiner
 {
 public:
     Refiner(Domain1D& domain);
     virtual ~Refiner() {}
 
+    //! Set grid refinement criteria
+    /*!
+     *  @param ratio Maximum ratio between grid spacing at adjacent intervals.
+     *      E.g. `(x[j+1] - x[j]) / (x[j] - x[j-1]) < ratio`
+     *  @param slope Maximum fractional change in the value of each solution
+     *      component between adjacent grid points
+     *  @param curve Maximum fractional change in the derivative of each
+     *      solution component between adjacent grid points.
+     *  @param prune Threshold for removing unnecessary grid points. `prune`
+     *      should be smaller than both `slope` and `curve`. Set `prune <= 0`
+     *      to disable pruning.
+     */
     void setCriteria(doublereal ratio = 10.0,
                      doublereal slope = 0.8,
                      doublereal curve = 0.8,
@@ -26,6 +39,8 @@ public:
     void setActive(int comp, bool state = true) {
         m_active[comp] = state;
     }
+
+    //! Set the maximum number of points allowed in the domain
     void setMaxPoints(int npmax) {
         m_npmax = npmax;
     }
@@ -54,6 +69,7 @@ public:
         return (m_keep[j] != -1);
     }
     double value(const double* x, size_t i, size_t j);
+
     double maxRatio() {
         return m_ratio;
     }
diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp
index 585488a0d..5114dbdbb 100644
--- a/src/oneD/Domain1D.cpp
+++ b/src/oneD/Domain1D.cpp
@@ -1,6 +1,5 @@
 /**
  * @file Domain1D.cpp
- *
  */
 
 #include "cantera/oneD/Domain1D.h"
@@ -165,7 +164,6 @@ void Domain1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     }
 }
 
-// called to set up initial grid, and after grid refinement
 void Domain1D::setupGrid(size_t n, const doublereal* z)
 {
     if (n > 1) {
@@ -176,17 +174,12 @@ void Domain1D::setupGrid(size_t n, const doublereal* z)
     }
 }
 
-
 void drawline()
 {
     writelog("\n-------------------------------------"
              "------------------------------------------");
 }
 
-
-/**
- * Print the solution.
- */
 void Domain1D::showSolution(const doublereal* x)
 {
     size_t nn = m_nv/5;
@@ -235,8 +228,6 @@ void Domain1D::showSolution(const doublereal* x)
     writelog("\n");
 }
 
-
-// initial solution
 void Domain1D::_getInitialSoln(doublereal* x)
 {
     for (size_t j = 0; j < m_points; j++) {
@@ -253,5 +244,4 @@ doublereal Domain1D::initialValue(size_t n, size_t j)
     return 0.0;
 }
 
-
 } // namespace
diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp
index 7d542f499..20cf58403 100644
--- a/src/oneD/MultiJac.cpp
+++ b/src/oneD/MultiJac.cpp
@@ -1,7 +1,5 @@
 /**
- *  @file MultiJac.cpp
- *
- *  Implementation file for class MultiJac
+ *  @file MultiJac.cpp Implementation file for class MultiJac
  */
 
 /*
@@ -47,10 +45,6 @@ void MultiJac::incrementDiagonal(int j, doublereal d)
     value(j,j) = m_ssdiag[j];
 }
 
-/**
- * Evaluate the Jacobian at x0. The array of residual values at x0
- * is supplied as an input.
- */
 void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
 {
     m_nevals++;
@@ -99,5 +93,3 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
 }
 
 } // namespace
-
-// $Log: MultiJac.cpp,v
diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp
index eb22e448b..e987bb93c 100644
--- a/src/oneD/MultiNewton.cpp
+++ b/src/oneD/MultiNewton.cpp
@@ -1,7 +1,5 @@
 /**
- *  @file MultiNewton.cpp
- *
- *  Damped Newton solver for 1D multi-domain problems
+ *  @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
  */
 
 /*
@@ -9,7 +7,6 @@
  */
 
 #include <vector>
-using namespace std;
 
 #include "cantera/oneD/MultiNewton.h"
 
@@ -40,7 +37,6 @@ public:
     }
 };
 
-
 /**
  * Return a damping coefficient that keeps the solution after taking one
  * Newton step between specified lower and upper bounds. This function only
@@ -49,7 +45,6 @@ public:
 doublereal bound_step(const doublereal* x, const doublereal* step,
                       Domain1D& r, int loglevel)
 {
-
     char buf[100];
     size_t np = r.nPoints();
     size_t nv = r.nComponents();
@@ -102,7 +97,6 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
     return fbound;
 }
 
-
 /**
  * This function computes the square of a weighted norm of a step
  * vector for one domain.
@@ -125,7 +119,6 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
  * solution component n in the domain. The second term,
  * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
  * n.
- *
  */
 doublereal norm_square(const doublereal* x,
                        const doublereal* step, Domain1D& r)
@@ -165,13 +158,10 @@ const string dashedline =
 const doublereal DampFactor = sqrt(2.0);
 const size_t NDAMP = 7;
 
-
-
 //-----------------------------------------------------------
 //                 MultiNewton methods
 //-----------------------------------------------------------
 
-
 MultiNewton::MultiNewton(int sz)
     : m_maxAge(5)
 {
@@ -186,9 +176,6 @@ MultiNewton::~MultiNewton()
     }
 }
 
-/**
- * Prepare for a new solution vector length.
- */
 void MultiNewton::resize(size_t sz)
 {
     m_n = sz;
@@ -198,10 +185,6 @@ void MultiNewton::resize(size_t sz)
     m_workarrays.clear();
 }
 
-
-/**
- * Compute the weighted 2-norm of 'step'.
- */
 doublereal MultiNewton::norm2(const doublereal* x,
                               const doublereal* step, OneDim& r) const
 {
@@ -216,11 +199,6 @@ doublereal MultiNewton::norm2(const doublereal* x,
     return sqrt(sum);
 }
 
-
-/**
- * Compute the undamped Newton step.  The residual function is
- * evaluated at x, but the Jacobian is not recomputed.
- */
 void MultiNewton::step(doublereal* x, doublereal* step,
                        OneDim& r, MultiJac& jac, int loglevel)
 {
@@ -281,12 +259,6 @@ void MultiNewton::step(doublereal* x, doublereal* step,
 #endif
 }
 
-
-/**
- * Return the factor by which the undamped Newton step 'step0'
- * must be multiplied in order to keep all solution components in
- * all domains between their specified lower and upper bounds.
- */
 doublereal MultiNewton::boundStep(const doublereal* x0,
                                   const doublereal* step0, const OneDim& r, int loglevel)
 {
@@ -299,20 +271,10 @@ doublereal MultiNewton::boundStep(const doublereal* x0,
     return fbound;
 }
 
-
-/**
- * On entry, step0 must contain an undamped Newton step for the
- * solution x0. This method attempts to find a damping coefficient
- * such that the next undamped step would have a norm smaller than
- * that of step0. If successful, the new solution after taking the
- * damped step is returned in x1, and the undamped step at x1 is
- * returned in step1.
- */
 int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
                           doublereal* x1, doublereal* step1, doublereal& s1,
                           OneDim& r, MultiJac& jac, int loglevel, bool writetitle)
 {
-
     // write header
     if (loglevel > 0 && writetitle) {
         writelog("\n\nDamped Newton iteration:\n");
@@ -404,12 +366,6 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
     }
 }
 
-
-/**
- * Find the solution to F(X) = 0 by damped Newton iteration.  On
- * entry, x0 contains an initial estimate of the solution.  On
- * successful return, x1 contains the converged solution.
- */
 int MultiNewton::solve(doublereal* x0, doublereal* x1,
                        OneDim& r, MultiJac& jac, int loglevel)
 {
@@ -510,11 +466,6 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
     return m;
 }
 
-
-/**
- * Get a pointer to an array of length m_n for temporary work
- * space.
- */
 doublereal* MultiNewton::getWorkArray()
 {
     doublereal* w = 0;
@@ -528,10 +479,6 @@ doublereal* MultiNewton::getWorkArray()
     return w;
 }
 
-/**
- * Release a work array by pushing its pointer onto the stack of
- * available arrays.
- */
 void MultiNewton::releaseWorkArray(doublereal* work)
 {
     m_workarrays.push_back(work);
diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp
index b707b1805..af7276c04 100644
--- a/src/oneD/OneDim.cpp
+++ b/src/oneD/OneDim.cpp
@@ -13,9 +13,6 @@ using namespace std;
 namespace Cantera
 {
 
-/**
- * Default constructor. Create an empty object.
- */
 OneDim::OneDim()
     : m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
       m_jac(0), m_newt(0),
@@ -30,11 +27,6 @@ OneDim::OneDim()
     //m_solve_time = 0.0;
 }
 
-
-/**
- * Construct a OneDim container for the domains pointed at by the
- * input vector of pointers.
-*/
 OneDim::OneDim(vector<Domain1D*> domains) :
     m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
     m_jac(0), m_newt(0),
@@ -57,7 +49,6 @@ OneDim::OneDim(vector<Domain1D*> domains) :
     resize();
 }
 
-
 size_t OneDim::domainIndex(const std::string& name)
 {
     for (size_t n = 0; n < m_nd; n++) {
@@ -69,13 +60,8 @@ size_t OneDim::domainIndex(const std::string& name)
     return npos;
 }
 
-
-/**
- * Domains are added left-to-right.
- */
 void OneDim::addDomain(Domain1D* d)
 {
-
     // if 'd' is not the first domain, link it to the last domain
     // added (the rightmost one)
     int n = static_cast<int>(m_dom.size());
@@ -98,7 +84,6 @@ void OneDim::addDomain(Domain1D* d)
     resize();
 }
 
-
 OneDim::~OneDim()
 {
     delete m_jac;
@@ -114,7 +99,6 @@ MultiNewton& OneDim::newton()
     return *m_newt;
 }
 
-//==============================================================================================================
 void OneDim::writeStats(int printTime)
 {
     saveStats();
@@ -134,20 +118,7 @@ void OneDim::writeStats(int printTime)
         writelog(buf);
     }
 }
-//==============================================================================================================
 
-/**
- * Save statistics on function and Jacobian evaluation, and reset
- * the counters. Statistics are saved only if the number of
- * Jacobian evaluations is greater than zero. The statistics saved
- * are
- *
- *    - number of grid points
- *    - number of Jacobian evaluations
- *    - CPU time spent evaluating Jacobians
- *    - number of non-Jacobian function evaluations
- *    - CPU time spent evaluating functions
- */
 void OneDim::saveStats()
 {
     if (m_jac) {
@@ -164,10 +135,6 @@ void OneDim::saveStats()
     }
 }
 
-
-/**
- * Call after one or more grids has been refined.
- */
 void OneDim::resize()
 {
     m_bw = 0;
@@ -232,7 +199,6 @@ void OneDim::resize()
     }
 }
 
-
 int OneDim::solve(doublereal* x, doublereal* xnew, int loglevel)
 {
     if (!m_jac_ok) {
@@ -254,15 +220,6 @@ void OneDim::evalSSJacobian(doublereal* x, doublereal* xnew)
     m_rdt = rdt_save;
 }
 
-/**
- * Return a pointer to the domain that contains component i of the
- * global solution vector. The domains are scanned right-to-left,
- * and the first one with starting location less or equal to i is
- * returned.
- *
- * 8/26/02 changed '<' to '<='  DGG
- *
- */
 Domain1D* OneDim::pointDomain(size_t i)
 {
     Domain1D* d = right();
@@ -275,11 +232,6 @@ Domain1D* OneDim::pointDomain(size_t i)
     return 0;
 }
 
-
-/**
- * Evaluate the multi-domain residual function, and return the
- * result in array r.
- */
 void OneDim::eval(size_t j, double* x, double* r, doublereal rdt, int count)
 {
     clock_t t0 = clock();
@@ -312,11 +264,6 @@ void OneDim::eval(size_t j, double* x, double* r, doublereal rdt, int count)
     }
 }
 
-
-/**
- * The 'infinity' (maximum magnitude) norm of the steady-state
- * residual. Used only for diagnostic output.
- */
 doublereal OneDim::ssnorm(doublereal* x, doublereal* r)
 {
     eval(npos, x, r, 0.0, 0);
@@ -327,10 +274,6 @@ doublereal OneDim::ssnorm(doublereal* x, doublereal* r)
     return ss;
 }
 
-
-/**
- * Prepare for time stepping with timestep dt.
- */
 void OneDim::initTimeInteg(doublereal dt, doublereal* x)
 {
     doublereal rdt_old = m_rdt;
@@ -351,12 +294,6 @@ void OneDim::initTimeInteg(doublereal dt, doublereal* x)
     }
 }
 
-
-/**
- * Prepare to solve the steady-state problem.  Set the reciprocal
- * of the time step to zero, and, if it was previously non-zero,
- * signal that a new Jacobian will be needed.
- */
 void OneDim::setSteadyMode()
 {
     m_rdt = 0.0;
@@ -370,11 +307,6 @@ void OneDim::setSteadyMode()
     }
 }
 
-/**
- * Initialize all domains. On the first call, this methods calls
- * the init method of each domain, proceeding from left to right.
- * Subsequent calls do nothing.
- */
 void OneDim::init()
 {
     if (!m_init) {
@@ -387,10 +319,6 @@ void OneDim::init()
     m_init = true;
 }
 
-
-/**
- * Signal that the current Jacobian is no longer valid.
- */
 void Domain1D::needJacUpdate()
 {
     if (m_container) {
@@ -399,18 +327,9 @@ void Domain1D::needJacUpdate()
     }
 }
 
-
-/**
- * Take time steps using Backward Euler.
- *
- *  nsteps   -- number of steps
- *  dt       -- initial step size
- *  loglevel -- controls amount of printed diagnostics
- */
 doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
                             doublereal* r, int loglevel)
 {
-
     // set the Jacobian age parameter to the transient value
     newton().setOptions(m_ts_jac_age);
 
@@ -468,12 +387,10 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
     return dt;
 }
 
-
 void OneDim::save(const std::string& fname, std::string id,
                   const std::string& desc, doublereal* sol,
                   int loglevel)
 {
-
     struct tm* newtime;
     time_t aclock;
     ::time(&aclock);                /* Get time in seconds */
@@ -519,7 +436,6 @@ void OneDim::save(const std::string& fname, std::string id,
     writelog("Solution saved to file "+fname+" as solution "+id+".\n", loglevel);
 }
 
-
 void Domain1D::setGrid(size_t n, const doublereal* z)
 {
     m_z.resize(n);
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index a3a9e5dbe..25cc95c99 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -20,20 +20,18 @@ static void sim1D_drawline()
     s += '\n';
     writelog(s.c_str());
 }
-//====================================================================================================================
+
 Sim1D::Sim1D() :
     OneDim()
 {
     //writelog("Sim1D default constructor\n");
 }
-//====================================================================================================================
+
 Sim1D::Sim1D(vector<Domain1D*>& domains) :
     OneDim(domains)
 {
-
-    // resize the internal solution vector and the wprk array,
-    // and perform domain-specific initialization of the
-    // solution vector.
+    // resize the internal solution vector and the wprk array, and perform
+    // domain-specific initialization of the solution vector.
 
     m_x.resize(size(), 0.0);
     m_xnew.resize(size(), 0.0);
@@ -49,13 +47,10 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
     m_steps.push_back(2);
     m_steps.push_back(5);
     m_steps.push_back(10);
-
 }
-//====================================================================================================================
 
 void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vector_fp& vals)
 {
-
     for (size_t dom=0; dom<m_nd; dom++) {
         Domain1D& d = domain(dom);
         size_t ncomp = d.nComponents();
@@ -67,15 +62,6 @@ void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vecto
     }
 }
 
-
-/**
- * Set a single value in the solution vector.
- * @param dom domain number, beginning with 0 for the leftmost domain.
- * @param comp component number
- * @param localPoint grid point within the domain, beginning with 0 for
- * the leftmost grid point in the domain.
- * @param value the value.
- */
 void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint,  doublereal value)
 {
     size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
@@ -85,13 +71,6 @@ void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint,  doublereal val
     m_x[iloc] = value;
 }
 
-
-/**
- * @param dom domain number, beginning with 0 for the leftmost domain.
- * @param comp component number
- * @param localPoint grid point within the domain, beginning with 0 for
- * the leftmost grid point in the domain.
- */
 doublereal Sim1D::value(size_t dom, size_t comp, size_t localPoint) const
 {
     size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
@@ -110,21 +89,6 @@ doublereal Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const
     return m_xnew[iloc];
 }
 
-
-/**
- * @param dom domain number, beginning with 0 for the leftmost domain.
- * @param comp component number
- * @param pos A vector of relative positions, beginning with 0.0 at the
- * left of the domain, and ending with 1.0 at the right of the domain.
- * @param values A vector of values corresponding to the relative position
- * locations.
- *
- * Note that the vector pos and values can have lengths
- * different than the number of grid points, but their lengths
- * must be equal. The values at the grid points will be
- * linearly interpolated based on the (pos, values)
- * specification.
- */
 void Sim1D::setProfile(size_t dom, size_t comp,
                        const vector_fp& pos, const vector_fp& values)
 {
@@ -141,7 +105,6 @@ void Sim1D::setProfile(size_t dom, size_t comp,
     }
 }
 
-
 void Sim1D::save(const std::string& fname, const std::string& id,
                  const std::string& desc, int loglevel)
 {
@@ -156,9 +119,6 @@ void Sim1D::saveResidual(const std::string& fname, const std::string& id,
     OneDim::save(fname, id, desc, &res[0], loglevel);
 }
 
-/**
- * Initialize the solution with a previously-saved solution.
- */
 void Sim1D::restore(const std::string& fname, const std::string& id,
                     int loglevel)
 {
@@ -203,7 +163,6 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
     finalize();
 }
 
-
 void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
 {
     size_t np = domain(dom).nPoints();
@@ -213,7 +172,6 @@ void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
     }
 }
 
-
 void Sim1D::showSolution(ostream& s)
 {
     for (size_t n = 0; n < m_nd; n++) {
@@ -248,7 +206,6 @@ void Sim1D::finalize()
     }
 }
 
-
 void Sim1D::setTimeStep(doublereal stepsize, size_t n, integer* tsteps)
 {
     m_tstep = stepsize;
@@ -258,7 +215,6 @@ void Sim1D::setTimeStep(doublereal stepsize, size_t n, integer* tsteps)
     }
 }
 
-
 int Sim1D::newtonSolve(int loglevel)
 {
     int m = OneDim::solve(DATA_PTR(m_x), DATA_PTR(m_xnew), loglevel);
@@ -273,7 +229,6 @@ int Sim1D::newtonSolve(int loglevel)
     }
 }
 
-
 void Sim1D::solve(int loglevel, bool refine_grid)
 {
     int new_points = 1;
@@ -388,10 +343,6 @@ void Sim1D::solve(int loglevel, bool refine_grid)
     }
 }
 
-
-/**
- * Refine the grid in all domains.
- */
 int Sim1D::refine(int loglevel)
 {
     int ianalyze, np = 0;
@@ -485,10 +436,6 @@ int Sim1D::refine(int loglevel)
     return np;
 }
 
-
-/**
- * Add node for fixed temperature point of freely propagating flame
- */
 int Sim1D::setFixedTemperature(doublereal t)
 {
     int np = 0;
@@ -598,11 +545,6 @@ void Sim1D::setAdiabaticFlame(void)
     }
 }
 
-/**
- * Set grid refinement criteria. If dom >= 0, then the settings
- * apply only to the specified domain.  If dom < 0, the settings
- * are applied to each domain.  @see Refiner::setCriteria.
- */
 void Sim1D::setRefineCriteria(int dom, doublereal ratio,
                               doublereal slope, doublereal curve, doublereal prune)
 {
@@ -630,7 +572,6 @@ void Sim1D::setGridMin(int dom, double gridmin)
     }
 }
 
-
 void Sim1D::setMaxGridPoints(int dom, int npoints)
 {
     if (dom >= 0) {
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index c0787518b..5714483bf 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -18,22 +18,10 @@ using namespace std;
 namespace Cantera
 {
 
-
-//-------------------  importSolution ------------------------
-
-/**
- * Import a previous solution to use as an initial estimate. The
- * previous solution may have been computed using a different
- * reaction mechanism. Species in the old and new mechanisms are
- * matched by name, and any species in the new mechanism that were
- * not in the old one are set to zero. The new solution is created
- * with the same number of grid points as in the old solution.
- */
 void importSolution(size_t points,
                     doublereal* oldSoln, IdealGasPhase& oldmech,
                     size_t size_new, doublereal* newSoln, IdealGasPhase& newmech)
 {
-
     // Number of components in old and new solutions
     size_t nv_old = oldmech.nSpecies() + 4;
     size_t nv_new = newmech.nSpecies() + 4;
@@ -81,7 +69,6 @@ void importSolution(size_t points,
     }
 }
 
-
 static void st_drawline()
 {
     writelog("\n-------------------------------------"
@@ -193,10 +180,6 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     setID("stagnation flow");
 }
 
-
-/**
- * Change the grid size. Called after grid refinement.
- */
 void StFlow::resize(size_t ncomponents, size_t points)
 {
     Domain1D::resize(ncomponents, points);
@@ -225,7 +208,6 @@ void StFlow::resize(size_t ncomponents, size_t points)
     m_z.resize(m_points);
 }
 
-
 void StFlow::setupGrid(size_t n, const doublereal* z)
 {
     resize(m_nv, n);
@@ -238,10 +220,6 @@ void StFlow::setupGrid(size_t n, const doublereal* z)
     }
 }
 
-
-/**
- * Install a transport manager.
- */
 void StFlow::setTransport(Transport& trans, bool withSoret)
 {
     m_trans = &trans;
@@ -275,11 +253,6 @@ void StFlow::enableSoret(bool withSoret)
     }
 }
 
-
-/**
- * Set the gas object state to be consistent with the solution at
- * point j.
- */
 void StFlow::setGas(const doublereal* x, size_t j)
 {
     m_thermo->setTemperature(T(x,j));
@@ -288,11 +261,6 @@ void StFlow::setGas(const doublereal* x, size_t j)
     m_thermo->setPressure(m_press);
 }
 
-
-/**
- * Set the gas state to be consistent with the solution at the
- * midpoint between j and j + 1.
- */
 void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
 {
     m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
@@ -305,7 +273,6 @@ void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
     m_thermo->setPressure(m_press);
 }
 
-
 void StFlow::_finalize(const doublereal* x)
 {
     size_t k, j;
@@ -353,23 +320,9 @@ void StFlow::_finalize(const doublereal* x)
     }
 }
 
-
-//------------------------------------------------------
-
-/**
- *  Evaluate the residual function for axisymmetric stagnation
- *  flow. If jpt is less than zero, the residual function is
- *  evaluated at all grid points. If jpt >= 0, then the residual
- *  function is only evaluated at grid points jpt-1, jpt, and
- *  jpt+1. This option is used to efficiently evaluate the
- *  Jacobian numerically.
- *
- */
-
 void StFlow::eval(size_t jg, doublereal* xg,
                   doublereal* rg, integer* diagg, doublereal rdt)
 {
-
     // if evaluating a Jacobian, and the global point is outside
     // the domain of influence for this domain, then skip
     // evaluating the residual
@@ -555,10 +508,6 @@ void StFlow::eval(size_t jg, doublereal* xg,
     }
 }
 
-/**
- * Update the transport properties at grid points in the range
- * from j0 to j1, based on solution x.
- */
 void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
 {
     if (m_transport_option == c_Mixav_Transport) {
@@ -602,13 +551,6 @@ void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
     }
 }
 
-
-
-
-
-/**
- * Print the solution.
- */
 void StFlow::showSolution(const doublereal* x)
 {
     size_t nn = m_nv/5;
@@ -660,10 +602,6 @@ void StFlow::showSolution(const doublereal* x)
     writelog("\n");
 }
 
-
-/**
- * Update the diffusive mass fluxes.
- */
 void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
 {
     size_t j, k, m;
@@ -706,7 +644,6 @@ void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
     }
 }
 
-
 string StFlow::componentName(size_t n) const
 {
     switch (n) {
@@ -727,11 +664,8 @@ string StFlow::componentName(size_t n) const
     }
 }
 
-
 size_t StFlow::componentIndex(const std::string& name) const
 {
-
-
     if (name=="u") {
         return 0;
     } else if (name=="V") {
@@ -751,7 +685,6 @@ size_t StFlow::componentIndex(const std::string& name) const
     return npos;
 }
 
-
 void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
 {
     Domain1D::restore(dom, soln, loglevel);
@@ -922,7 +855,6 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     }
 }
 
-
 XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
 {
     size_t k;
@@ -1061,7 +993,6 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
     diag[index(4,j)] = 0;
 }
 
-
 void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
                                integer* diag, doublereal rdt)
 {
diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp
index 3f850ad79..82fe8d345 100644
--- a/src/oneD/boundaries1D.cpp
+++ b/src/oneD/boundaries1D.cpp
@@ -23,7 +23,6 @@ Bdry1D::Bdry1D() : Domain1D(1, 1, 0.0),
     m_type = cConnectorType;
 }
 
-
 void Bdry1D::
 _init(size_t n)
 {
@@ -73,16 +72,10 @@ _init(size_t n)
     }
 }
 
-
-
-
 //----------------------------------------------------------
-//
 //   Inlet1D methods
-//
 //----------------------------------------------------------
 
-
 void Inlet1D::
 setMoleFractions(const std::string& xin)
 {
@@ -121,7 +114,6 @@ componentName(size_t n) const
 void Inlet1D::
 init()
 {
-
     _init(2);
 
     // set bounds (mdot, T)
@@ -158,7 +150,6 @@ init()
     }
 }
 
-
 void Inlet1D::
 eval(size_t jg, doublereal* xg, doublereal* rg,
      integer* diagg, doublereal rdt)
@@ -234,8 +225,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
             rb[4+k] += x[0]*(m_yin[k]);
         }
     }
-
-
 }
 
 XML_Node& Inlet1D::
@@ -275,8 +264,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(2,1);
 }
 
-
-
 //--------------------------------------------------
 //      Empty1D
 //--------------------------------------------------
@@ -339,8 +326,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-
-
 //--------------------------------------------------
 //      Symm1D
 //--------------------------------------------------
@@ -413,7 +398,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
     }
 }
 
-
 XML_Node& Symm1D::
 save(XML_Node& o, const doublereal* const soln)
 {
@@ -429,7 +413,6 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-
 //--------------------------------------------------
 //      Outlet1D
 //--------------------------------------------------
@@ -466,7 +449,6 @@ init()
     }
 }
 
-
 void Outlet1D::
 eval(size_t jg, doublereal* xg, doublereal* rg,
      integer* diagg, doublereal rdt)
@@ -522,7 +504,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
     }
 }
 
-
 XML_Node& Outlet1D::
 save(XML_Node& o, const doublereal* const soln)
 {
@@ -538,14 +519,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-
-
-
 //--------------------------------------------------
 //      OutletRes1D
 //--------------------------------------------------
 
-
 void OutletRes1D::
 setMoleFractions(const std::string& xres)
 {
@@ -609,7 +586,6 @@ init()
     }
 }
 
-
 void OutletRes1D::
 eval(size_t jg, doublereal* xg, doublereal* rg,
      integer* diagg, doublereal rdt)
@@ -670,7 +646,6 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
     }
 }
 
-
 XML_Node& OutletRes1D::
 save(XML_Node& o, const doublereal* const soln)
 {
@@ -704,15 +679,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-
 //-----------------------------------------------------------
-//
 //  Surf1D
-//
 //-----------------------------------------------------------
 
-
-
 string Surf1D::componentName(size_t n) const
 {
     switch (n) {
@@ -739,7 +709,6 @@ init()
     setTolerances(1, &rtol, 1, &atol);
 }
 
-
 void Surf1D::
 eval(size_t jg, doublereal* xg, doublereal* rg,
      integer* diagg, doublereal rdt)
@@ -793,17 +762,10 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-
-
-
 //-----------------------------------------------------------
-//
 //  ReactingSurf1D
-//
 //-----------------------------------------------------------
 
-
-
 string ReactingSurf1D::componentName(size_t n) const
 {
     if (n == 0) {
@@ -842,7 +804,6 @@ init()
     setTolerances(m_nv, DATA_PTR(rtol), m_nv, DATA_PTR(atol));
 }
 
-
 void ReactingSurf1D::
 eval(size_t jg, doublereal* xg, doublereal* rg,
      integer* diagg, doublereal rdt)
@@ -967,4 +928,3 @@ restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(m_nsp+1,1);
 }
 }
-
diff --git a/src/oneD/refine.cpp b/src/oneD/refine.cpp
index 48044c1ac..20eb00278 100644
--- a/src/oneD/refine.cpp
+++ b/src/oneD/refine.cpp
@@ -8,7 +8,6 @@ using namespace std;
 
 namespace Cantera
 {
-
 static void r_drawline()
 {
     string s(78,'#');
@@ -26,11 +25,9 @@ Refiner::Refiner(Domain1D& domain) :
     m_thresh = std::sqrt(std::numeric_limits<double>::epsilon());
 }
 
-
 int Refiner::analyze(size_t n, const doublereal* z,
                      const doublereal* x)
 {
-
     if (n >= m_npmax) {
         writelog("max number of grid points reached ("+int2str(m_npmax)+".\n");
         return -2;
@@ -214,7 +211,6 @@ void Refiner::show()
     }
 }
 
-
 int Refiner::getNewGrid(int n, const doublereal* z,
                         int nn, doublereal* zn)
 {