diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h
index 89717ffb2..e1738bc30 100644
--- a/include/cantera/oneD/Domain1D.h
+++ b/include/cantera/oneD/Domain1D.h
@@ -1,6 +1,5 @@
-/**
- *  @file Domain1D.h
- */
+ //! @file Domain1D.h
+
 /*
  *  Copyright 2002 California Institute of Technology
  */
@@ -82,10 +81,8 @@ public:
         return *m_container;
     }
 
-    /*!
-     * Specify the container object for this domain, and the
-     * position of this domain in the list.
-     */
+    //! Specify the container object for this domain, and the position of this
+    //! domain in the list.
     void setContainer(OneDim* c, size_t index) {
         m_container = c;
         m_index = index;
@@ -116,9 +113,9 @@ public:
     }
 
     /*!
-     * 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.
+     * 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.
      */
     virtual void init() {  }
 
@@ -126,9 +123,9 @@ public:
     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.
+     * 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.
      */
     virtual void resize(size_t nv, size_t np) {
         // if the number of components is being changed, then a
@@ -162,7 +159,7 @@ public:
         return m_nv;
     }
 
-    //! Check that the specified component index is in range
+    //! Check that the specified component index is in range.
     //! Throws an exception if n is greater than nComponents()-1
     void checkComponentIndex(size_t n) const {
         if (n >= m_nv) {
@@ -170,7 +167,7 @@ public:
         }
     }
 
-    //! Check that an array size is at least nComponents()
+    //! Check that an array size is at least nComponents().
     //! Throws an exception if nn is less than nComponents(). Used before calls
     //! which take an array pointer.
     void checkComponentArraySize(size_t nn) const {
@@ -184,7 +181,7 @@ public:
         return m_points;
     }
 
-    //! Check that the specified point index is in range
+    //! Check that the specified point index is in range.
     //! Throws an exception if n is greater than nPoints()-1
     void checkPointIndex(size_t n) const {
         if (n >= m_points) {
@@ -192,7 +189,7 @@ public:
         }
     }
 
-    //! Check that an array size is at least nPoints()
+    //! Check that an array size is at least nPoints().
     //! Throws an exception if nn is less than nPoints(). Used before calls
     //! which take an array pointer.
     void checkPointArraySize(size_t nn) const {
@@ -239,21 +236,21 @@ public:
 
     //! Set tolerances for time-stepping mode
     /*!
-     *  @param rtol Relative tolerance
-     *  @param atol Absolute tolerance
-     *  @param n    component index these tolerances apply to. If set to -1
-     *      (the default), these tolerances will be applied to all solution
+     * @param rtol Relative tolerance
+     * @param atol Absolute tolerance
+     * @param n    component index these tolerances apply to. If set to -1 (the
+     *      default), these tolerances will be applied to all solution
      *      components.
      */
     void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos);
 
     //! Set tolerances for steady-state mode
     /*!
-     *  @param rtol Relative tolerance
-     *  @param atol Absolute tolerance
-     *  @param n    component index these tolerances apply to. If set to -1
-     *      (the default), these tolerances will be applied to all solution
-     *      components.
+     * @param rtol Relative tolerance
+     * @param atol Absolute tolerance
+     * @param n    component index these tolerances apply to. If set to -1 (the
+     *     default), these tolerances will be applied to all solution
+     *     components.
      */
     void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos);
 
@@ -277,18 +274,18 @@ public:
         return m_min[n];
     }
 
+    //! Prepare to do time stepping with time step dt
     /*!
-     * Prepare to do time stepping with time step dt. Copy the internally-
-     * stored solution at the last time step to array x0.
+     * 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
     /*!
-     * Prepare to solve the steady-state problem. Set the internally-stored
-     * reciprocal of the time step to 0,0
+     * Set the internally-stored reciprocal of the time step to 0.0
      */
     void setSteadyMode() {
         m_rdt = 0.0;
@@ -363,15 +360,14 @@ public:
 
     //! Save the current solution for this domain into an XML_Node
     /*!
-     *  Base class version of the general domain1D save function. Derived
-     *  classes should call the base class method in addition to saving their
-     *  own data.
+     * Base class version of the general domain1D save function. Derived classes
+     * should call the base class method in addition to saving their own data.
      *
-     *  @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.
-     *  @return     XML_Node created to represent this domain
+     * @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.
+     * @return     XML_Node created to represent this domain
      */
     virtual XML_Node& save(XML_Node& o, const doublereal* const sol);
 
@@ -425,16 +421,16 @@ public:
     }
 
     /**
-     * The index of the first (i.e., left-most) grid point
-     * belonging to this domain.
+     * The index of the first (i.e., left-most) grid point belonging to this
+     * domain.
      */
     size_t firstPoint() const {
         return m_jstart;
     }
 
     /**
-     * The index of the last (i.e., right-most) grid point
-     * belonging to this domain.
+     * The index of the last (i.e., right-most) grid point belonging to this
+     * domain.
      */
     size_t lastPoint() const {
         return m_jstart + m_points - 1;
@@ -540,11 +536,10 @@ public:
     virtual void setupGrid(size_t n, const doublereal* z);
 
     /**
-     * Writes some or all initial solution values into the global
-     * solution array, beginning at the location pointed to by
-     * x. This method is called by the Sim1D constructor, and
-     * allows default values or ones that have been set locally
-     * prior to installing this domain into the container to be
+     * Writes some or all initial solution values into the global solution
+     * array, beginning at the location pointed to by x. This method is called
+     * by the Sim1D constructor, and allows default values or ones that have
+     * been set locally prior to installing this domain into the container to be
      * written to the global solution vector.
      */
     virtual void _getInitialSoln(doublereal* x);
@@ -553,14 +548,12 @@ public:
     virtual doublereal initialValue(size_t n, size_t j);
 
     /**
-     * In some cases, a domain may need to set parameters that
-     * depend on the initial solution estimate. In such cases, the
-     * parameters may be set in method _finalize. This method is
-     * called just before the Newton solver is called, and the x
-     * array is guaranteed to be the local solution vector for
-     * this domain that will be used as the initial guess. If no
-     * such parameters need to be set, then method _finalize does
-     * not need to be overloaded.
+     * In some cases, a domain may need to set parameters that depend on the
+     * initial solution estimate. In such cases, the parameters may be set in
+     * method _finalize. This method is called just before the Newton solver is
+     * called, and the x array is guaranteed to be the local solution vector for
+     * this domain that will be used as the initial guess. If no such parameters
+     * need to be set, then method _finalize does not need to be overloaded.
      */
     virtual void _finalize(const doublereal* x) {}
 
diff --git a/include/cantera/oneD/Inlet1D.h b/include/cantera/oneD/Inlet1D.h
index 3497498dd..75a4591c2 100644
--- a/include/cantera/oneD/Inlet1D.h
+++ b/include/cantera/oneD/Inlet1D.h
@@ -23,15 +23,14 @@ const int LeftInlet = 1;
 const int RightInlet = -1;
 
 /**
- * The base class for boundaries between one-dimensional spatial
- * domains. The boundary may have its own internal variables, such
- * as surface species coverages.
+ * The base class for boundaries between one-dimensional spatial domains. The
+ * boundary may have its own internal variables, such as surface species
+ * coverages.
  *
- * The boundary types are an inlet, an outlet, a symmetry plane,
- * and a surface.
+ * The boundary types are an inlet, an outlet, a symmetry plane, and a surface.
  *
- * The public methods are all virtual, and the base class
- * implementations throw exceptions.
+ * The public methods are all virtual, and the base class implementations throw
+ * exceptions.
  * @ingroup onedim
  */
 class Bdry1D : public Domain1D
@@ -198,8 +197,8 @@ public:
 };
 
 /**
- * A symmetry plane. The axial velocity u = 0, and all other
- * components have zero axial gradients.
+ * A symmetry plane. The axial velocity u = 0, and all other components have
+ * zero axial gradients.
  * @ingroup onedim
  */
 class Symm1D : public Bdry1D
@@ -225,7 +224,8 @@ public:
 
 
 /**
- *  An outlet.
+ * An outlet.
+ * @ingroup onedim
  */
 class Outlet1D : public Bdry1D
 {
@@ -291,10 +291,10 @@ protected:
 };
 
 /**
- * A non-reacting surface. The axial velocity is zero
- * (impermeable), as is the transverse velocity (no slip). The
- * temperature is specified, and a zero flux condition is imposed
- * for the species.
+ * A non-reacting surface. The axial velocity is zero (impermeable), as is the
+ * transverse velocity (no slip). The temperature is specified, and a zero flux
+ * condition is imposed for the species.
+ * @ingroup onedim
  */
 class Surf1D : public Bdry1D
 {
diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h
index 79c2af23e..d82c9c38d 100644
--- a/include/cantera/oneD/MultiJac.h
+++ b/include/cantera/oneD/MultiJac.h
@@ -1,6 +1,4 @@
-/**
- *  @file MultiJac.h
- */
+//! @file MultiJac.h
 
 /*
  *  Copyright 2002 California Institute of Technology
@@ -16,10 +14,10 @@ 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
- * 1D domains, with different variables in each domain.
+ * 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 1D domains, with different variables
+ * in each domain.
  * @ingroup onedim
  */
 class MultiJac : public BandMatrix
@@ -28,10 +26,10 @@ public:
     MultiJac(OneDim& r);
 
     /**
-     * Evaluate the Jacobian at x0. The unperturbed residual
-     * function is resid0, which must be supplied on input. The
-     * third parameter 'rdt' is the reciprocal of the time
-     * step. If zero, the steady-state Jacobian is evaluated.
+     * Evaluate the Jacobian at x0. The unperturbed residual function is resid0,
+     * which must be supplied on input. The third parameter 'rdt' is the
+     * reciprocal of the time step. If zero, the steady-state Jacobian is
+     * evaluated.
      */
     void eval(doublereal* x0, doublereal* resid0, double rdt);
 
@@ -45,8 +43,7 @@ public:
         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;
     }
@@ -70,10 +67,10 @@ public:
     void incrementDiagonal(int j, doublereal d);
 
 protected:
-    //!  Residual evaluator for this Jacobian
+    //! 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;
 
diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h
index 40d82071d..eb3d6b26b 100644
--- a/include/cantera/oneD/MultiNewton.h
+++ b/include/cantera/oneD/MultiNewton.h
@@ -1,6 +1,4 @@
-/**
- *  @file MultiNewton.h
- */
+//! @file MultiNewton.h
 
 /*
  *  Copyright 2002 California Institute of Technology
@@ -43,12 +41,11 @@ public:
                          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.
+     * 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,
@@ -59,9 +56,9 @@ public:
                      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.
+     * 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);
diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h
index 8bf3d96d2..e4f5f6229 100644
--- a/include/cantera/oneD/OneDim.h
+++ b/include/cantera/oneD/OneDim.h
@@ -57,7 +57,7 @@ public:
 
     size_t domainIndex(const std::string& name);
 
-    //! Check that the specified domain index is in range
+    //! Check that the specified domain index is in range.
     //! Throws an exception if n is greater than nDomains()-1
     void checkDomainIndex(size_t n) const {
         if (n >= m_nd) {
@@ -65,7 +65,7 @@ public:
         }
     }
 
-    //! Check that an array size is at least nDomains()
+    //! Check that an array size is at least nDomains().
     //! Throws an exception if nn is less than nDomains(). Used before calls
     //! which take an array pointer.
     void checkDomainArraySize(size_t nn) const {
@@ -99,10 +99,8 @@ public:
         return m_nvars[jg];
     }
 
-    /**
-     * Location in the solution vector of the first component of
-     *  global point jg.
-     */
+    //! Location in the solution vector of the first component of global point
+    //! jg.
     size_t loc(size_t jg) {
         return m_loc[jg];
     }
@@ -198,11 +196,12 @@ public:
     double timeStep(int nsteps, double dt, double* x,
                     double* r, int loglevel);
 
-    //! Write statistics about the number of iterations and Jacobians at each grid level
+    //! Write statistics about the number of iterations and Jacobians at each
+    //! grid level
     /*!
-     *  @param printTime  Boolean that indicates whether time should be printed out
-     *                    The default is true. It's turned off for test problems where
-     *                    we don't want to print any times
+     *  @param printTime  Boolean that indicates whether time should be printed
+     *                    out The default is true. It's turned off for test
+     *                    problems where we don't want to print any times
      */
     void writeStats(int printTime = 1);
 
@@ -233,11 +232,11 @@ public:
      * 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
+     * - 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();
 
@@ -254,21 +253,22 @@ public:
 protected:
     void evalSSJacobian(doublereal* x, doublereal* xnew);
 
-    doublereal m_tmin; // minimum timestep size
-    doublereal m_tmax; // maximum timestep size
-    doublereal m_tfactor; // factor time step is multiplied by
-    // if time stepping fails ( < 1 )
+    doublereal m_tmin; //!< minimum timestep size
+    doublereal m_tmax; //!< maximum timestep size
 
-    std::unique_ptr<MultiJac> m_jac; // Jacobian evaluator
-    std::unique_ptr<MultiNewton> m_newt; // Newton iterator
-    doublereal m_rdt; // reciprocal of time step
-    bool m_jac_ok; // if true, Jacobian is current
+    //! factor time step is multiplied by  if time stepping fails ( < 1 )
+    doublereal m_tfactor;
+
+    std::unique_ptr<MultiJac> m_jac; //!< Jacobian evaluator
+    std::unique_ptr<MultiNewton> m_newt; //!< Newton iterator
+    doublereal m_rdt; //!< reciprocal of time step
+    bool m_jac_ok; //!< if true, Jacobian is current
 
     //! number of domains
     size_t m_nd;
 
-    size_t m_bw; // Jacobian bandwidth
-    size_t m_size; // solution vector size
+    size_t m_bw; //!< Jacobian bandwidth
+    size_t m_size; //!< solution vector size
 
     std::vector<Domain1D*> m_dom, m_connect, m_bulk;
 
diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h
index 23ccb5fb4..bb73a4209 100644
--- a/include/cantera/oneD/Sim1D.h
+++ b/include/cantera/oneD/Sim1D.h
@@ -20,18 +20,17 @@ 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.
-     * The domain pointers must be entered in left-to-right order --- i.e.,
-     * the pointer to the leftmost domain is domain[0], the pointer to the
-     * domain to its right is domain[1], etc.
+     *     The domain pointers must be entered in left-to-right order --- i.e.,
+     *     the pointer to the leftmost domain is domain[0], the pointer to the
+     *     domain to its right is domain[1], etc.
      */
     Sim1D(std::vector<Domain1D*>& domains);
 
@@ -128,10 +127,10 @@ public:
 
     //! Set the minimum grid spacing in the specified domain(s).
     /*!
-     *  @param dom Domain index. If dom == -1, the specified spacing
-                   is applied to all domains.
-        @param gridmin The minimum allowable grid spacing [m]
-    */
+     * @param dom Domain index. If dom == -1, the specified spacing is applied
+     *            to all domains.
+     * @param gridmin The minimum allowable grid spacing [m]
+     */
     void setGridMin(int dom, double gridmin);
 
     //! Initialize the solution with a previously-saved solution.
@@ -169,7 +168,8 @@ private:
     /// Calls method _finalize in each domain.
     void finalize();
 
-    /*! Wrapper around the Newton solver.
+    //! Wrapper around the Newton solver
+    /*!
      * @return 0 if successful, -1 on failure
      */
     int newtonSolve(int loglevel);
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 65f6cdaa9..4f257a9ef 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -1,6 +1,5 @@
-/**
- * @file StFlow.h
- */
+//! @file StFlow.h
+
 // Copyright 2001  California Institute of Technology
 
 #ifndef CT_STFLOW_H
@@ -13,6 +12,7 @@
 
 namespace Cantera
 {
+
 //------------------------------------------
 //   constants
 //------------------------------------------
@@ -33,7 +33,7 @@ class Transport;
 
 /**
  *  This class represents 1D flow domains that satisfy the one-dimensional
- *  similarity solution for chemically-reacting, axisymmetric, flows.
+ *  similarity solution for chemically-reacting, axisymmetric flows.
  *  @ingroup onedim
  */
 class StFlow : public Domain1D
@@ -85,9 +85,8 @@ public:
         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;
     }
@@ -196,10 +195,10 @@ public:
 
     //! Set the emissivities for the boundary values
     /*!
-    *   Reads the emissivities for the left and right boundary values in the
-    *   radiative term and writes them into the variables, which are used for
-    *   the calculation.
-    */
+     * Reads the emissivities for the left and right boundary values in the
+     * radiative term and writes them into the variables, which are used for the
+     * calculation.
+     */
     void setBoundaryEmissivities(doublereal e_left, doublereal e_right) {
         if (e_left < 0 || e_left > 1) {
             throw CanteraError("setBoundaryEmissivities",
@@ -474,11 +473,10 @@ protected:
     std::vector<bool> m_do_species;
     int m_transport_option;
 
-    // flag for the radiative heat loss
+    //! flag for the radiative heat loss
     bool m_do_radiation;
 
-    // radiative heat loss vector
-    // vector which contains the values of the radiative heat loss
+    //! radiative heat loss vector
     vector_fp m_qdotRadiation;
 
     // fixed T and Y values
diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp
index 92c1f5fc6..d1a342a61 100644
--- a/src/oneD/MultiJac.cpp
+++ b/src/oneD/MultiJac.cpp
@@ -1,6 +1,4 @@
-/**
- *  @file MultiJac.cpp Implementation file for class MultiJac
- */
+//! @file MultiJac.cpp Implementation file for class MultiJac
 
 /*
  *  Copyright 2002 California Institute of Technology
diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp
index 0bc603044..3d4197baa 100644
--- a/src/oneD/MultiNewton.cpp
+++ b/src/oneD/MultiNewton.cpp
@@ -1,6 +1,4 @@
-/**
- *  @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
- */
+//! @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
 
 /*
  *  Copyright 2001 California Institute of Technology
@@ -83,8 +81,8 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
 }
 
 /**
- * This function computes the square of a weighted norm of a step
- * vector for one domain.
+ * This function computes the square of a weighted norm of a step vector for one
+ * domain.
  *
  * @param x     Solution vector for this domain.
  * @param step  Newton step vector for this domain.
@@ -99,10 +97,10 @@ doublereal bound_step(const doublereal* x, const doublereal* step,
  * \f[
  *     w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
  * \f]
- * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
- * component n, and multiplies the average magnitude of
- * solution component n in the domain. The second term,
- * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component n.
+ * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for component n,
+ * and multiplies the average magnitude of 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)
@@ -130,16 +128,12 @@ doublereal norm_square(const doublereal* x,
 
 } // end unnamed-namespace
 
-//-----------------------------------------------------------
-//                  constants
-//-----------------------------------------------------------
 
+// constants
 const doublereal DampFactor = sqrt(2.0);
 const size_t NDAMP = 7;
 
-//-----------------------------------------------------------
-//                 MultiNewton methods
-//-----------------------------------------------------------
+// ---------------- MultiNewton methods ----------------
 
 MultiNewton::MultiNewton(int sz)
     : m_maxAge(5)
@@ -237,18 +231,15 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
     // compute the multiplier to keep all components in bounds
     doublereal fbound = boundStep(x0, step0, r, loglevel-1);
 
-    // if fbound is very small, then x0 is already close to the
-    // boundary and step0 points out of the allowed domain. In
-    // this case, the Newton algorithm fails, so return an error
-    // condition.
+    // if fbound is very small, then x0 is already close to the boundary and
+    // step0 points out of the allowed domain. In this case, the Newton
+    // algorithm fails, so return an error condition.
     if (fbound < 1.e-10) {
         debuglog("\nAt limits.\n", loglevel);
         return -3;
     }
 
-    //--------------------------------------------
-    //           Attempt damped step
-    //--------------------------------------------
+    // ---------- Attempt damped step ----------
 
     // damping coefficient starts at 1.0
     doublereal damp = 1.0;
@@ -262,8 +253,7 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
             x1[j] = ff*step0[j] + x0[j];
         }
 
-        // compute the next undamped step that would result if x1
-        // is accepted
+        // compute the next undamped step that would result if x1 is accepted
         step(x1, step1, r, jac, loglevel-1);
 
         // compute the weighted norm of step1
@@ -278,20 +268,19 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
                      jac.nEvals(), jac.age(), m_maxAge);
         }
 
-        // if the norm of s1 is less than the norm of s0, then
-        // accept this damping coefficient. Also accept it if this
-        // step would result in a converged solution. Otherwise,
-        // decrease the damping coefficient and try again.
+        // if the norm of s1 is less than the norm of s0, then accept this
+        // damping coefficient. Also accept it if this step would result in a
+        // converged solution. Otherwise, decrease the damping coefficient and
+        // try again.
         if (s1 < 1.0 || s1 < s0) {
             break;
         }
         damp /= DampFactor;
     }
 
-    // If a damping coefficient was found, return 1 if the
-    // solution after stepping by the damped step would represent
-    // a converged solution, and return 0 otherwise. If no damping
-    // coefficient could be found, return -2.
+    // If a damping coefficient was found, return 1 if the solution after
+    // stepping by the damped step would represent a converged solution, and
+    // return 0 otherwise. If no damping coefficient could be found, return -2.
     if (m < NDAMP) {
         if (s1 > 1.0) {
             return 0;
@@ -352,8 +341,8 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
         }
         frst = false;
 
-        // Successful step, but not converged yet. Take the damped
-        // step, and try again.
+        // Successful step, but not converged yet. Take the damped step, and try
+        // again.
         if (m == 0) {
             copy(x1, x1 + m_n, m_x.begin());
         } else if (m == 1) {
@@ -361,9 +350,9 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
             jac.setAge(0); // for efficient sensitivity analysis
             break;
         } else if (m < 0) {
-            // If dampStep fails, first try a new Jacobian if an old
-            // one was being used. If it was a new Jacobian, then
-            // return -1 to signify failure.
+            // If dampStep fails, first try a new Jacobian if an old one was
+            // being used. If it was a new Jacobian, then return -1 to signify
+            // failure.
             if (jac.age() > 1) {
                 forceNewJac = true;
                 if (nJacReeval > 3) {
diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp
index f0ad83d5a..8a011cfec 100644
--- a/src/oneD/OneDim.cpp
+++ b/src/oneD/OneDim.cpp
@@ -70,8 +70,7 @@ void OneDim::addDomain(Domain1D* d)
         m_bulk.push_back(d);
     }
 
-    // add it also to the global domain list, and set its
-    // container and position
+    // add it also to the global domain list, and set its container and position
     m_dom.push_back(d);
     d->setContainer(this, m_nd);
     m_nd++;
@@ -265,14 +264,13 @@ void OneDim::initTimeInteg(doublereal dt, doublereal* x)
     doublereal rdt_old = m_rdt;
     m_rdt = 1.0/dt;
 
-    // if the stepsize has changed, then update the transient
-    // part of the Jacobian
+    // if the stepsize has changed, then update the transient part of the
+    // Jacobian
     if (fabs(rdt_old - m_rdt) > Tiny) {
         m_jac->updateTransient(m_rdt, m_mask.data());
     }
 
-    // iterate over all domains, preparing each one to begin
-    // time stepping
+    // iterate over all domains, preparing each one to begin time stepping
     Domain1D* d = left();
     while (d) {
         d->initTimeInteg(dt, x);
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index e4e3c2d4b..6d1fc252c 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -349,9 +349,9 @@ int Sim1D::refine(int loglevel)
                     xnew.push_back(value(n, i, m));
                 }
 
-                // now check whether a new point is needed in the
-                // interval to the right of point m, and if so, add
-                // entries to znew and xnew for this new point
+                // now check whether a new point is needed in the interval to
+                // the right of point m, and if so, add entries to znew and xnew
+                // for this new point
                 if (r.newPointNeeded(m) && m + 1 < npnow) {
                     // add new point at midpoint
                     zmid = 0.5*(d.grid(m) + d.grid(m+1));
@@ -372,10 +372,9 @@ int Sim1D::refine(int loglevel)
         dsize.push_back(znew.size() - nstart);
     }
 
-    // At this point, the new grid znew and the new solution
-    // vector xnew have been constructed, but the domains
-    // themselves have not yet been modified.  Now update each
-    // domain with the new grid.
+    // At this point, the new grid znew and the new solution vector xnew have
+    // been constructed, but the domains themselves have not yet been modified.
+    // Now update each domain with the new grid.
 
     size_t gridstart = 0, gridsize;
     for (size_t n = 0; n < m_nd; n++) {
@@ -413,7 +412,8 @@ int Sim1D::setFixedTemperature(doublereal t)
         Domain1D& d = domain(n);
         size_t comp = d.nComponents();
 
-        // loop over points in the current grid to determine where new point is needed.
+        // loop over points in the current grid to determine where new point is
+        // needed.
         FreeFlame* d_free = dynamic_cast<FreeFlame*>(&domain(n));
         size_t npnow = d.nPoints();
         size_t nstart = znew.size();
@@ -466,10 +466,9 @@ int Sim1D::setFixedTemperature(doublereal t)
         dsize.push_back(znew.size() - nstart);
     }
 
-    // At this point, the new grid znew and the new solution
-    // vector xnew have been constructed, but the domains
-    // themselves have not yet been modified.  Now update each
-    // domain with the new grid.
+    // At this point, the new grid znew and the new solution vector xnew have
+    // been constructed, but the domains themselves have not yet been modified.
+    // Now update each domain with the new grid.
     size_t gridstart = 0, gridsize;
     for (n = 0; n < m_nd; n++) {
         Domain1D& d = domain(n);
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 43978bebe..2ec35d0df 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -1,6 +1,5 @@
-/**
- * @file StFlow.cpp
- */
+//! @file StFlow.cpp
+
 // Copyright 2002  California Institute of Technology
 
 #include "cantera/oneD/StFlow.h"
@@ -218,9 +217,8 @@ void StFlow::_finalize(const doublereal* x)
 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
+    // if evaluating a Jacobian, and the global point is outside the domain of
+    // influence for this domain, then skip evaluating the residual
     if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
         return;
     }
@@ -251,9 +249,7 @@ void StFlow::eval(size_t jg, doublereal* xg,
 
     size_t j, k;
 
-    //-----------------------------------------------------
-    //              update properties
-    //-----------------------------------------------------
+    // ------------ update properties ------------
 
     updateThermo(x, j0, j1);
     // update transport properties only if a Jacobian is not being evaluated
@@ -344,25 +340,23 @@ void StFlow::eval(size_t jg, doublereal* xg,
         if (j == 0) {
             // these may be modified by a boundary object
 
-            // Continuity. This propagates information right-to-left,
-            // since rho_u at point 0 is dependent on rho_u at point 1,
-            // but not on mdot from the inlet.
+            // Continuity. This propagates information right-to-left, since
+            // rho_u at point 0 is dependent on rho_u at point 1, but not on
+            // mdot from the inlet.
             rsd[index(c_offset_U,0)] =
                 -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
                 -(density(1)*V(x,1) + density(0)*V(x,0));
 
-            // the inlet (or other) object connected to this one
-            // will modify these equations by subtracting its values
-            // for V, T, and mdot. As a result, these residual equations
-            // will force the solution variables to the values for
-            // the boundary object
+            // the inlet (or other) object connected to this one will modify
+            // these equations by subtracting its values for V, T, and mdot. As
+            // a result, these residual equations will force the solution
+            // variables to the values for the boundary object
             rsd[index(c_offset_V,0)] = V(x,0);
             rsd[index(c_offset_T,0)] = T(x,0);
             rsd[index(c_offset_L,0)] = -rho_u(x,0);
 
-            // The default boundary condition for species is zero
-            // flux. However, the boundary object may modify
-            // this.
+            // The default boundary condition for species is zero flux. However,
+            // the boundary object may modify this.
             sum = 0.0;
             for (k = 0; k < m_nsp; k++) {
                 sum += Y(x,k,0);
@@ -694,10 +688,9 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
                 soln[index(2,j)] = x[j];
             }
 
-            // For fixed-temperature simulations, use the
-            // imported temperature profile by default.  If
-            // this is not desired, call setFixedTempProfile
-            // *after* restoring the solution.
+            // For fixed-temperature simulations, use the imported temperature
+            // profile by default.  If this is not desired, call
+            // setFixedTempProfile *after* restoring the solution.
             vector_fp zz(np);
             for (size_t jj = 0; jj < np; jj++) {
                 zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp
index dd4c25606..e130ad367 100644
--- a/src/oneD/boundaries1D.cpp
+++ b/src/oneD/boundaries1D.cpp
@@ -1,6 +1,5 @@
-/**
- * @file boundaries1D.cpp
- */
+//! @file boundaries1D.cpp
+
 // Copyright 2002-3  California Institute of Technology
 
 #include "cantera/oneD/Inlet1D.h"
@@ -55,8 +54,7 @@ void Bdry1D::_init(size_t n)
         }
     }
 
-    // if this is not the last domain, see what is connected on
-    // the right
+    // if this is not the last domain, see what is connected on the right
     if (m_index + 1 < container().nDomains()) {
         Domain1D& r = container().domain(m_index+1);
         if (r.domainType() == cFlowType) {
@@ -73,9 +71,7 @@ void Bdry1D::_init(size_t n)
     }
 }
 
-//----------------------------------------------------------
-//   Inlet1D methods
-//----------------------------------------------------------
+// ---------------- Inlet1D methods ----------------
 
 void Inlet1D::setMoleFractions(const std::string& xin)
 {
@@ -120,10 +116,9 @@ void Inlet1D::init()
     setSteadyTolerances(1e-4, 1e-5);
     setTransientTolerances(1e-4, 1e-5);
 
-    // if a flow domain is present on the left, then this must be
-    // a right inlet. Note that an inlet object can only be a
-    // terminal object - it cannot have flows on both the left and
-    // right
+    // if a flow domain is present on the left, then this must be a right inlet.
+    // Note that an inlet object can only be a terminal object - it cannot have
+    // flows on both the left and right
     if (m_flow_left) {
         m_ilr = RightInlet;
         m_flow = m_flow_left;
@@ -173,21 +168,20 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
         xb = x + 2;
         rb = r + 2;
 
-        // The first flow residual is for u. This, however, is not
-        // modified by the inlet, since this is set within the flow
-        // domain from the continuity equation.
+        // The first flow residual is for u. This, however, is not modified by
+        // the inlet, since this is set within the flow domain from the
+        // continuity equation.
 
         // spreading rate. The flow domain sets this to V(0),
         // so for finite spreading rate subtract m_V0.
         rb[1] -= m_V0;
 
-        // The third flow residual is for T, where it is set to
-        // T(0).  Subtract the local temperature to hold the flow
-        // T to the inlet T.
+        // The third flow residual is for T, where it is set to T(0).  Subtract
+        // the local temperature to hold the flow T to the inlet T.
         rb[2] -= x[1];
 
-        // The flow domain sets this to -rho*u. Add mdot to
-        // specify the mass flow rate.
+        // The flow domain sets this to -rho*u. Add mdot to specify the mass
+        // flow rate.
         rb[3] += x[0];
 
         // add the convective term to the species residual equations
@@ -195,9 +189,8 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
             rb[4+k] += x[0]*m_yin[k];
         }
 
-        // if the flow is a freely-propagating flame, mdot is not
-        // specified.  Set mdot equal to rho*u, and also set
-        // lambda to zero.
+        // if the flow is a freely-propagating flame, mdot is not specified.
+        // Set mdot equal to rho*u, and also set lambda to zero.
         if (!m_flow->fixed_mdot()) {
             m_mdot = m_flow->density(0)*xb[0];
             r[0] = m_mdot - x[0];
@@ -252,9 +245,7 @@ void Inlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(2,1);
 }
 
-//--------------------------------------------------
-//      Empty1D
-//--------------------------------------------------
+// ------------- Empty1D -------------
 
 string Empty1D::componentName(size_t n) const
 {
@@ -305,9 +296,7 @@ void Empty1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-//--------------------------------------------------
-//      Symm1D
-//--------------------------------------------------
+// -------------- Symm1D --------------
 
 string Symm1D::componentName(size_t n) const
 {
@@ -384,9 +373,7 @@ void Symm1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-//--------------------------------------------------
-//      Outlet1D
-//--------------------------------------------------
+// -------- Outlet1D --------
 
 string Outlet1D::componentName(size_t n) const
 {
@@ -477,9 +464,7 @@ void Outlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-//--------------------------------------------------
-//      OutletRes1D
-//--------------------------------------------------
+// -------- OutletRes1D --------
 
 void OutletRes1D::setMoleFractions(const std::string& xres)
 {
@@ -626,9 +611,7 @@ void OutletRes1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-//-----------------------------------------------------------
-//  Surf1D
-//-----------------------------------------------------------
+// -------- Surf1D --------
 
 string Surf1D::componentName(size_t n) const
 {
@@ -701,9 +684,7 @@ void Surf1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
-//-----------------------------------------------------------
-//  ReactingSurf1D
-//-----------------------------------------------------------
+// -------- ReactingSurf1D --------
 
 string ReactingSurf1D::componentName(size_t n) const
 {