From 32fed991cf3fdb14a6b2a06a0cfdcb9599ca8add Mon Sep 17 00:00:00 2001
From: Dave Goodwin <dgg@caltech.edu>
Date: Tue, 13 May 2003 19:43:30 +0000
Subject: [PATCH] -

---
 Cantera/src/DASPK.h               |   5 +-
 Cantera/src/InterfaceKinetics.cpp |   9 +-
 Cantera/src/ReactionPath.cpp      |  16 +
 Cantera/src/ReactionPath.h        |  13 +-
 Cantera/src/funcs.cpp             |  23 ++
 Cantera/src/importCTML.cpp        |   2 +-
 Cantera/src/misc.cpp              |  15 +-
 Cantera/src/oneD/Inlet1D.h        | 552 +++++++++---------------------
 Cantera/src/oneD/Makefile.in      |   2 +-
 Cantera/src/oneD/MultiNewton.cpp  |   4 +-
 Cantera/src/oneD/OneDim.cpp       | 141 +++++---
 Cantera/src/oneD/OneDim.h         |  89 +++--
 Cantera/src/oneD/Resid1D.h        | 215 ++++++++++--
 Cantera/src/oneD/Sim1D.cpp        | 306 +++++++++++++++++
 Cantera/src/oneD/Sim1D.h          |  89 +++++
 Cantera/src/oneD/StFlow.cpp       | 173 +++++-----
 Cantera/src/oneD/StFlow.h         | 142 +++++---
 Cantera/src/oneD/newton_utils.cpp |  17 +-
 Cantera/src/oneD/refine.cpp       | 215 ++++++++++++
 Cantera/src/oneD/refine.h         |  47 +++
 20 files changed, 1443 insertions(+), 632 deletions(-)
 create mode 100644 Cantera/src/oneD/Sim1D.cpp
 create mode 100644 Cantera/src/oneD/Sim1D.h
 create mode 100644 Cantera/src/oneD/refine.cpp
 create mode 100644 Cantera/src/oneD/refine.h

diff --git a/Cantera/src/DASPK.h b/Cantera/src/DASPK.h
index b3d31cc73..d119420b1 100755
--- a/Cantera/src/DASPK.h
+++ b/Cantera/src/DASPK.h
@@ -3,8 +3,9 @@
  *  @file DASPK.h
  *
  *  Header file for class DASPK
- *
- *  $Author$
+ */
+
+/*  $Author$
  *  $Date$
  *  $Revision$
  *
diff --git a/Cantera/src/InterfaceKinetics.cpp b/Cantera/src/InterfaceKinetics.cpp
index f45dea5c8..95c0f16ab 100644
--- a/Cantera/src/InterfaceKinetics.cpp
+++ b/Cantera/src/InterfaceKinetics.cpp
@@ -39,6 +39,7 @@ namespace Cantera {
 
     doublereal SurfPhase::
     enthalpy_mole() const {
+        if (m_n0 <= 0.0) return 0.0; 
         _updateThermo();
         return mean_X(m_h0.begin());
     }
@@ -58,10 +59,14 @@ namespace Cantera {
     }
 
     void SurfPhase::
-    getActivityConcentrations(doublereal* c) const { getConcentrations(c); }
+    getActivityConcentrations(doublereal* c) const { 
+        getConcentrations(c); 
+    }
 
     doublereal SurfPhase::
-    standardConcentration(int k) const { return m_n0/size(k); }
+    standardConcentration(int k) const {
+        return m_n0/size(k); 
+    }
 
     doublereal SurfPhase::
     logStandardConc(int k) const {
diff --git a/Cantera/src/ReactionPath.cpp b/Cantera/src/ReactionPath.cpp
index 8c5aa3b53..87368a448 100755
--- a/Cantera/src/ReactionPath.cpp
+++ b/Cantera/src/ReactionPath.cpp
@@ -28,6 +28,22 @@
 
 namespace Cantera {
      
+    /// add a path to or from this node
+    void SpeciesNode::addPath(Path* path) { 
+        m_paths.push_back(path);
+        if (path->begin() == this) m_out += path->flow();
+        else if (path->end() == this) m_in += path->flow();
+        else throw CanteraError("addPath","path added to wrong node");
+    }
+        
+    void SpeciesNode::printPaths() {
+        for (int i = 0; i < m_paths.size(); i++) {
+            cout << m_paths[i]->begin()->name << " -->  "
+                 << m_paths[i]->end()->name << ":   "
+                 << m_paths[i]->flow() << endl;
+        }
+    }
+
 
     /** 
      * Construct a path connecting two species nodes.
diff --git a/Cantera/src/ReactionPath.h b/Cantera/src/ReactionPath.h
index 2602144e3..5e51dfca2 100755
--- a/Cantera/src/ReactionPath.h
+++ b/Cantera/src/ReactionPath.h
@@ -46,7 +46,7 @@ namespace Cantera {
 
         /// Default constructor
         SpeciesNode() : number(-1), name(""), value(0.0), 
-                        visible(false) {}
+                        visible(false), m_in(0.0), m_out(0.0) {}
 
         /// Destructor
         virtual ~SpeciesNode() {}
@@ -75,10 +75,17 @@ namespace Cantera {
         int nPaths() const { return m_paths.size(); }
 
         /// add a path to or from this node
-        void addPath(Path* path) { m_paths.push_back(path); }
+        void addPath(Path* path);
 
+        double outflow() {return m_out;}
+        double inflow() {return m_in;}
+        double netOutflow() {return m_out - m_in;}
+        
+        void printPaths();
+
+        
     protected:
-
+        double m_in, m_out;
         path_list m_paths;
     };
 
diff --git a/Cantera/src/funcs.cpp b/Cantera/src/funcs.cpp
index 9456df6f7..4e6712da4 100755
--- a/Cantera/src/funcs.cpp
+++ b/Cantera/src/funcs.cpp
@@ -23,6 +23,29 @@ extern "C" {
 
 namespace Cantera {
 
+    /**
+     * Linearly interpolate a function defined on a discrete grid.
+     * vector xpts contains a monotonic sequence of grid points, and 
+     * vector fpts contains function values defined at these points.
+     * The value returned is the linear interpolate at point x.
+     * If x is outside the range of xpts, the value of fpts at the 
+     * nearest end is returned.
+     */
+
+    doublereal linearInterp(doublereal x, const vector_fp& xpts, 
+        const vector_fp& fpts) {
+        if (x <= xpts[0]) return fpts[0];
+        if (x >= xpts.back()) return fpts.back();
+        const doublereal* loc = lower_bound(xpts.begin(), xpts.end(), x);
+        int iloc = int(loc - xpts.begin()) - 1;
+        doublereal ff = fpts[iloc] + 
+            (x - xpts[iloc])*(fpts[iloc + 1] 
+                - fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
+        return ff;
+    }
+
+
+
     doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, 
         int maxdeg, int& ndeg, doublereal eps, doublereal* r) {
         integer nn = n;
diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp
index 9b87a3e6c..4b346d217 100755
--- a/Cantera/src/importCTML.cpp
+++ b/Cantera/src/importCTML.cpp
@@ -771,7 +771,7 @@ namespace Cantera {
         int nt = th.size();
 
         // for each referenced phase, attempt to find its id among those
-        // phases that have already been built. 
+        // phases specified. 
         bool phase_ok;
 
         string phase_id;
diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp
index 1d76e5ca4..f9e2398c3 100755
--- a/Cantera/src/misc.cpp
+++ b/Cantera/src/misc.cpp
@@ -30,15 +30,16 @@ namespace Cantera {
      */
     class Application {
     public:
-        Application() : linelen(0) {}
+        Application() : linelen(0), stop_on_error(false), write_log_to_cout(true) {}
         virtual ~Application(){}
         vector<string> inputDirs;
         vector<string> errorMessage;
         vector<string> warning;
         vector<string> errorRoutine;
         string msglog;
-        bool stop_on_error;
         size_t linelen;
+        bool stop_on_error;
+        bool write_log_to_cout;
         map<string, string>     options;
     };
 
@@ -87,9 +88,9 @@ namespace Cantera {
         int i = __app->errorMessage.size();
         if (i == 0) return;
         f << endl << endl;
-        f << "**********************" << endl;
-        f << "    Cantera Error!    " << endl;
-        f << "**********************" << endl << endl;
+        f << "************************************************" << endl;
+        f << "                Cantera Error!                  " << endl;
+        f << "************************************************" << endl << endl;
         int j;
         for (j = 0; j < i; j++) {
             f << endl;
@@ -259,6 +260,10 @@ namespace Cantera {
             __app->msglog += "\n";
             __app->linelen = 0;
         }
+        if (__app->write_log_to_cout) {
+            cout << __app->msglog;
+            clearlog();
+        }
     }
     void writelog(const char* msg) {writelog(string(msg));}
     void getlog(string& s) {
diff --git a/Cantera/src/oneD/Inlet1D.h b/Cantera/src/oneD/Inlet1D.h
index 685d1a4ed..6b53e2ee5 100644
--- a/Cantera/src/oneD/Inlet1D.h
+++ b/Cantera/src/oneD/Inlet1D.h
@@ -1,13 +1,28 @@
+/**
+ * @file Inlet1D.h
+ *
+ * Boundary objects for one-dimensional simulations.
+ *
+ */
+
+/*
+ * $Author$
+ * $Revision$
+ * $Date$
+ *
+ * Copyright 2002-3  California Institute of Technology
+ */
+
 
 #ifndef CT_BDRY1D_H
 #define CT_BDRY1D_H
 
 #include "Resid1D.h"
-//#include "surfacePhase.h"
-//#include "surfKinetics.h"
+#include "../SurfPhase.h"
+#include "../InterfaceKinetics.h"
 #include "StFlow.h"
 #include "OneDim.h"
-#include "ctml.h"
+#include "../ctml.h"
 
 namespace Cantera {
 
@@ -27,56 +42,75 @@ namespace Cantera {
      */
     class Bdry1D : public Resid1D {
     public:
-        Bdry1D() : Resid1D(1, 1, 0.0) {}
+
+        Bdry1D();
+ 
         virtual ~Bdry1D() {}
 
         /// Initialize.
-        virtual void init(){err("init");}
+        virtual void init() { _init(1); }
 
         /// Set the temperature.
-        virtual void setTemperature(doublereal t){err("setTemperature");}
+        virtual void setTemperature(doublereal t){m_temp = t;}
 
         /// Temperature [K].
-        virtual doublereal temperature() {err("temperature"); return 0.0;}
+        virtual doublereal temperature() {return m_temp;}
 
         /// Set the mole fractions by specifying a string.
         virtual void setMoleFractions(string xin){err("setMoleFractions");}
 
         /// Set the mole fractions by specifying an array.
         virtual void setMoleFractions(doublereal* xin){err("setMoleFractions");}
+
         /// Mass fraction of species k.
         virtual doublereal massFraction(int k) {err("massFraction"); return 0.0;}
+
         /// Set the total mass flow rate.
-        virtual void setMdot(doublereal mdot){err("setMdot");}
+        virtual void setMdot(doublereal mdot){m_mdot = mdot;}
 
         /// The total mass flow rate [kg/m2/s].
-        virtual doublereal mdot() {err("mdot"); return 0.0;}
+        virtual doublereal mdot() {return m_mdot;}
 
     protected:
+
+        void _init(int n);
+
+        StFlow *m_flow_left, *m_flow_right;
+        int m_ilr, m_left_nv, m_right_nv;
+        int m_left_loc, m_right_loc;
+        int m_left_points;
+        int m_nv, m_left_nsp, m_right_nsp;
+        int m_sp_left, m_sp_right;
+        int m_start_left, m_start_right;
+        ThermoPhase *m_phase_left, *m_phase_right;
+        doublereal m_temp, m_mdot;
+
     private:
         void err(string method) {
-            throw CanteraError("Bdry1D::"+method, "attempt to call base class method "+method);
+            throw CanteraError("Bdry1D::"+method, 
+                "attempt to call base class method "+method);
         }
     };
 
 
+    /**
+     * An inlet.
+     */
     class Inlet1D : public Bdry1D {
 
     public:
 
-        Inlet1D(int ilr = 1) {
-            m_type = cInletType; 
-            m_flow = 0;
-            m_ilr = ilr;
+        /**
+         * Constructor. Create a new Inlet1D instance. If invoked
+         * without parameters, a left inlet (facing right) is
+         * constructed).
+         */
+        Inlet1D() : m_V0(0.0), m_nsp(0), m_flow(0) {
+            m_type = cInletType;
+            m_xstr = "";
         }
         virtual ~Inlet1D(){}
 
-        /// Set the inlet temperature
-        virtual void setTemperature(doublereal t) {
-            m_temp = t;
-            needJacUpdate();
-        }
-
         /// set spreading rate
         virtual void setSpreadRate(doublereal V0) {
             m_V0 = V0;
@@ -88,432 +122,170 @@ namespace Cantera {
             return m_V0;
         }
 
-        /// Temperature [K].
-        doublereal temperature() {return m_temp;}
-
-        virtual void setMoleFractions(string xin) {
-            m_xstr = xin;
+        virtual void showSolution(ostream& s, const doublereal* x) {
+            s << "-------------------  Inlet " << domainIndex() << " ------------------- " << endl;
+            s << "  mdot:        " << m_mdot << " kg/m^2/s" << "   " << x[0] << endl;
+            s << "  temperature: " << m_temp << " K" << "    " << x[1] << endl;
             if (m_flow) {
-                m_flow->phase().setMoleFractionsByName(xin);
-                m_flow->phase().getMassFractions(m_yin.begin());
-                needJacUpdate();
-            }
-        }
-
-        virtual void setMoleFractions(doublereal* xin) {
-            if (m_flow) {
-                m_flow->phase().setMoleFractions(xin);
-                m_flow->phase().getMassFractions(m_yin.begin());
-                needJacUpdate();
-            }
-        }
-
-        virtual doublereal massFraction(int k) {return m_yin[k];}
-
-        virtual void setMdot(doublereal mdot) { m_mdot = mdot; }
-
-        virtual string componentName(int n) const { 
-            switch (n) {
-            case 0: return "mdot"; break;
-            case 1: return "temperature"; break;
-            default: return "unknown";
-            }
-        }
-
-        virtual void init() {
-            if (m_index < 0) {
-                throw CanteraError("Inlet1D", 
-                    "install in container before calling init.");
-            }
-            resize(2,1);
-
-            // set bounds
-            const doublereal lower[2] = {-1.0e5, 200.0};
-            const doublereal upper[2] = {1.0e5, 1.e5};
-            setBounds(2, lower, 2, upper);
-
-            // set tolerances
-            vector_fp rtol(2, 1e-4);
-            vector_fp atol(2, 1.e-5);
-            setTolerances(2, rtol.begin(), 2, atol.begin());
-
-            // if a flow domain is present on the left, then this must
-            // be a right inlet
-            if (m_index > 0) {
-                Resid1D& r = container().domain(m_index-1);
-                if (r.domainType() == cFlowType) {
-                    m_ilr = RightInlet;
-                    m_flow = (StFlow*)&r;
-                }
-                else 
-                    throw CanteraError("Inlet1D::init",
-                        "Inlet domains can only be connected to a flow domain.");
-            }
-            else {
-                if (container().nDomains() > 1) {
-                    Resid1D& r = container().domain(1);
-                    if (r.domainType() == cFlowType) {
-                        m_ilr = LeftInlet;
-                        m_flow = (StFlow*)&r;
+                s << "  mass fractions: " << endl;
+                for (int k = 0; k < m_flow->phase().nSpecies(); k++) {
+                    if (m_yin[k] != 0.0) {
+                        s << "      " << m_flow->phase().speciesName(k) 
+                          << "  " << m_yin[k] << endl;
                     }
-                    else 
-                        throw CanteraError("Inlet1D::init",
-                            "An inlet domain can only be connected to a flow domain.");
                 }
-                else 
-                    throw CanteraError("Inlet1D::init",
-                        "An inlet domain must be connected to a flow domain.");
             }
-            // components = u, V, T, lambda, + mass fractions
-            m_nsp = m_flow->nComponents() - 4;
-            m_yin.resize(m_nsp, 0.0);
-            if (m_xstr != "") 
-                setMoleFractions(m_xstr);
-            else
-                m_yin[0] = 1.0;
+            s << endl;
         }
 
+        virtual void _getInitialSoln(doublereal* x) {
+            x[0] = m_mdot;
+            x[1] = m_temp;
+        }
 
+        virtual void _finalize(const doublereal* x) {
+            ; //m_mdot = x[0];
+            //m_temp = x[1];
+        }
+
+        virtual void setMoleFractions(string xin);
+        virtual void setMoleFractions(doublereal* xin);
+        virtual doublereal massFraction(int k) {return m_yin[k];}
+        virtual string componentName(int n) const;
+        virtual void init();
         virtual void eval(int jg, doublereal* xg, doublereal* rg, 
-            integer* diagg, doublereal rdt) {
-            int k;
-            if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
-
-            // start of local part of global arrays
-            doublereal* x = xg + loc();
-            doublereal* r = rg + loc();
-            integer* diag = diagg + loc();
-            doublereal *xb, *rb;
-
-            // residual equations for the two local variables
-            r[0] = m_mdot - x[0];
-            r[1] = m_temp - x[1];
-
-            // both are algebraic constraints
-            diag[0] = 0;
-            diag[1] = 0;
-
-            // if it is a left inlet, then the flow solution vector
-            // starts 2 to the right in the global solution vector
-            if (m_ilr == LeftInlet) {
-                xb = x + 2;
-                rb = r + 2;
-
-                // If the energy equation is being solved, then
-                // the flow domain set this residual to T(0). 
-                // Subtract the inlet temperature.
-                if (m_flow->doEnergy(0)) {
-                    rb[2] -= x[1];   // T
-                }
-
-                // spreading rate. Flow domain sets this to V(0),
-                // so for finite spreading rate subtract m_V0.
-                rb[1] -= m_V0;
-
-                rb[3] += x[0];       // lambda
-                for (k = 0; k < m_nsp; k++) {
-                    rb[4+k] += x[0]*m_yin[k];
-                }
-            }
-
-            // right inlet.
-            else {
-                int boffset = m_flow->nComponents();
-                xb = x - boffset;
-                rb = r - boffset;
-                rb[1] -= m_V0;
-                rb[2] -= x[1]; // T
-                xb[0] += x[0]; // u
-                for (k = 0; k < m_nsp; k++) 
-                    rb[4+k] += x[0]*m_yin[k];
-            }                
-        }
-
-        virtual void save(XML_Node& o, doublereal* soln) {
-            doublereal* s = soln + loc();
-            XML_Node& inlt = o.addChild("inlet");
-            for (int k = 0; k < 2; k++) {
-                ctml::addFloat(inlt, componentName(k), s[k], "", "",0.0, 1.0);
-            }
-        }
+            integer* diagg, doublereal rdt);
+        virtual void save(XML_Node& o, doublereal* soln);
                 
+
     protected:
 
         int m_ilr;
-        doublereal m_mdot, m_temp, m_V0;
-        StFlow *m_flow;
+        doublereal m_V0;
         int m_nsp;
         vector_fp m_yin;
         string m_xstr;
+        StFlow *m_flow;
     };
 
 
+    /**
+     * A symmetry plane. The axial velocity u = 0, and all other
+     * components have zero axial gradients.
+     */
     class Symm1D : public Bdry1D {
 
     public:
 
-        Symm1D(int ilr = 1) {
+        Symm1D() {
             m_type = cSymmType; 
-            m_flow = 0;
-            m_ilr = ilr;
         }
         virtual ~Symm1D(){}
 
-        virtual string componentName(int n) const { 
-            switch (n) {
-            case 0: return "dummy"; break;
-            default: return "<unknown>";
-            }
-        }
-
-        virtual void init() {
-            if (m_index < 0) {
-                throw CanteraError("Symm1D", 
-                    "install in container before calling init.");
-            }
-            resize(1,1);
-
-            // set bounds
-            const doublereal lower = -1.0e5;
-            const doublereal upper =  1.0e5;
-            setBounds(1, &lower, 1, &upper);
-
-            // set tolerances
-            doublereal rtol = 1e-4;
-            doublereal atol = 1.e-5;
-            setTolerances(1, &rtol, 1, &atol);
-
-            if (m_index > 0) {
-                Resid1D& r = container().domain(m_index-1);
-                if (r.domainType() == cFlowType) {
-                    m_ilr = -1;
-                    m_flow = (StFlow*)&r;
-                }
-                else 
-                    throw CanteraError("Symm1D::init",
-                        "Symmetry planes can only be connected to flow domains.");
-            }
-            else {
-                if (container().nDomains() > 1) {
-                    Resid1D& r = container().domain(1);
-                    if (r.domainType() == cFlowType) {
-                        m_ilr = 1;
-                        m_flow = (StFlow*)&r;
-                    }
-                    else 
-                        throw CanteraError("Symm1D::init",
-                            "Symmetry planes can only be connected to flow domains.");
-                }
-                else 
-                    throw CanteraError("Symm1D::init",
-                        "A symmetry plane must be connected to a flow domain.");
-            }
-            m_nsp = m_flow->nComponents() - 4;
-        }
+        virtual string componentName(int n) const;
 
+        virtual void init();
 
         virtual void eval(int jg, doublereal* xg, doublereal* rg, 
-            integer* diagg, doublereal rdt) {
-            int k;
-            if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
+            integer* diagg, doublereal rdt);
 
-            // start of local part of global arrays
-            doublereal* x = xg + loc();
-            doublereal* r = rg + loc();
-            integer* diag = diagg + loc();
-            doublereal *xb, *rb;
-            // integer *db = diag + loc();
-
-            r[0] = x[0];
-            diag[0] = 0;
-            int nc = m_flow->nComponents();
-
-            if (m_ilr == 1) {
-                xb = x + 1;
-                rb = r + 1;
-                rb[0] = xb[0];
-                rb[1] = xb[1] - xb[1 + nc];
-                if (m_flow->doEnergy(0)) {
-                    rb[2] = xb[2] - xb[2 + nc];
-                }
-                rb[3] = xb[3] - xb[3 + nc];
-                for (k = 0; k < m_nsp; k++) {
-                    rb[4+k] = xb[4+k] - xb[4+k+nc];
-                }
-            }
-            else {
-                xb = x - nc;
-                rb = r - nc;
-                rb[0] = xb[0];
-                rb[1] = xb[1] - xb[1 - nc];
-                //               if (m_flow->doEnergy(0)) {
-                rb[2] = xb[2] - xb[2 - nc];
-                    // }
-                rb[3] = xb[3] - xb[3 - nc];
-                for (k = 0; k < m_nsp; k++) {
-                    rb[4+k] = xb[4+k] - xb[4+k-nc];
-                }
-            }                
-        }
-
-        virtual void save(XML_Node& o, doublereal* soln) {
-	    // doublereal* s = soln + loc();
-            // XML_Node& symm = o.addChild("symmetry_plane");
-            (void) o.addChild("symmetry_plane");
-        }
+        virtual void save(XML_Node& o, doublereal* soln);
 
     protected:
 
-        int m_ilr;
-        StFlow *m_flow;
-        int m_nsp;
     };
 
 
-
-    /////////////////////////////////////////////////////////////
-    //
-    // surf1D
-    //
-    ////////////////////////////////////////////////////////////
-
+    /**
+     * 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.
+     */
     class Surf1D : public Bdry1D {
 
     public:
 
-        Surf1D(int ilr = 1) {
+        Surf1D() {
             m_type = cSurfType; 
-            m_flow = 0;
-            m_ilr = ilr;
         }
         virtual ~Surf1D(){}
 
-        /// Set the surface temperature
-        virtual void setTemperature(doublereal t) {
-            m_temp = t;
-            needJacUpdate();
-        }
-
-        /// Temperature [K].
-        doublereal temperature() {return m_temp;}
-
-        virtual string componentName(int n) const { 
-            switch (n) {
-            case 0: return "dummy"; break;
-            case 1: return "temperature"; break;
-            default: return "<unknown>";
-            }
-        }
-
-        virtual void init() {
-            if (m_index < 0) {
-                throw CanteraError("Surf1D", 
-                    "install in container before calling init.");
-            }
-            resize(2,1);
-
-            // set bounds
-            const doublereal lower[2] = {-1.0e5, 200.0};
-            const doublereal upper[2] = {1.0e5, 1.e5};
-            setBounds(2, lower, 2, upper);
-
-            // set tolerances
-            vector_fp rtol(2, 1e-4);
-            vector_fp atol(2, 1.e-5);
-            setTolerances(2, rtol.begin(), 2, atol.begin());
-
-            if (m_index > 0) {
-                Resid1D& r = container().domain(m_index-1);
-                if (r.domainType() == cFlowType) {
-                    m_ilr = -1;
-                    m_flow = (StFlow*)&r;
-                }
-                else 
-                    throw CanteraError("Surf1D::init",
-                        "Surface domains can only be connected to a flow domain.");
-            }
-            else {
-                if (container().nDomains() > 1) {
-                    Resid1D& r = container().domain(1);
-                    if (r.domainType() == cFlowType) {
-                        m_ilr = 1;
-                        m_flow = (StFlow*)&r;
-                    }
-                    else 
-                        throw CanteraError("Surf1D::init",
-                            "A surface domain can only be connected to a flow domain.");
-                }
-                else 
-                    throw CanteraError("Surf1D::init",
-                        "A surface domain must be connected to a flow domain.");
-            }
-            m_nsp = m_flow->nComponents() - 4;
-        }
+        virtual string componentName(int n) const;
 
+        virtual void init();
 
         virtual void eval(int jg, doublereal* xg, doublereal* rg, 
-            integer* diagg, doublereal rdt) {
-            int k;
-            if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
+            integer* diagg, doublereal rdt);
 
-            // start of local part of global arrays
-            doublereal* x = xg + loc();
-            doublereal* r = rg + loc();
-            integer* diag = diagg + loc();
-            doublereal *xb, *rb;
-            //integer *db = diag + loc();
+        virtual void save(XML_Node& o, doublereal* soln);
 
-            r[0] = x[0];
-            r[1] = m_temp - x[1];
-            diag[0] = 0;
-            diag[1] = 0;
-            int nc = m_flow->nComponents();
-
-            if (m_ilr == 1) {
-                xb = x + 2;
-                rb = r + 2;
-                rb[0] = xb[0];
-                rb[1] = xb[1];   // T
-                //                if (m_flow->doEnergy(0)) {
-                rb[2] = xb[2] - x[1];   // T
-                    //}
-                rb[3] = xb[3] - xb[3 + nc];   // lambda
-                for (k = 0; k < m_nsp; k++) {
-                    rb[4+k] += xb[4+k] - xb[4+k+nc];
-                }
-            }
-            else {
-                xb = x - nc;
-                rb = r - nc;
-                rb[0] = xb[0];
-                rb[1] = xb[1];   // T
-                //                if (m_flow->doEnergy(0)) {
-                rb[2] = xb[2] - x[1];   // T
-                    //}
-                rb[3] = xb[3] - xb[3 - nc];   // lambda
-                for (k = 0; k < m_nsp; k++) {
-                    rb[4+k] += xb[4+k] - xb[4+k-nc];
-                }
-            }
+        virtual void _getInitialSoln(doublereal* x) {
+            x[0] = m_temp;
         }
 
-        virtual void save(XML_Node& o, doublereal* soln) {
-            doublereal* s = soln + loc();
-            XML_Node& srf = o.addChild("surface");
-            for (int k = 1; k < 2; k++) {
-                ctml::addFloat(srf, componentName(k), s[k], "", "",0.0, 1.0);
-            }
+        virtual void _finalize(const doublereal* x) {
+            ; //m_temp = x[0];
+        }
+
+        virtual void showSolution(ostream& s, const doublereal* x) {
+            s << "-------------------  Surface " << domainIndex() << " ------------------- " << endl;
+            s << "  temperature: " << m_temp << " K" << "    " << x[0] << endl;
         }
-                
     protected:
 
-        int m_ilr;
-        doublereal m_temp;
-        StFlow *m_flow;
-        int m_nsp;
     };
 
-
 }
 
+// //     };
+
+
+
+// //     /////////////////////////////////////////////////////////////
+// //     //
+// //     // surf1D
+// //     //
+// //     ////////////////////////////////////////////////////////////
+// // #ifdef WITH_CHEMSURF
+
+// //     class ChemSurf1D : public Bdry1D {
+
+// //     public:
+
+// //         ChemSurf1D(InterfaceKinetics* skin = 0) {
+// //             m_type = cSurfType;
+// //             m_kin = 0;
+// //             m_sphase = 0;
+// //             m_nsp = 0;
+// //             if (skin) setKinetics(skin);
+// //         }
+// //         virtual ~ChemSurf1D(){}
+
+// //         void setKinetics(InterfaceKinetics* kin);
+
+// //         virtual string componentName(int n) const;
+// //         virtual void init();
+
+// //         virtual void eval(int jg, doublereal* xg, doublereal* rg, 
+// //             integer* diagg, doublereal rdt);
+
+// //         virtual void save(XML_Node& o, doublereal* soln);
+                
+// //     protected:
+
+// //         int m_ilr, m_nsp;
+//         InterfaceKinetics* m_kin;
+//         SurfPhase* m_sphase;
+//         vector_fp m_work;
+//         const doublereal *m_molwt_right, *m_molwt_left; 
+//         int m_start_surf;
+//         vector<ThermoPhase*> m_bulk;
+//         vector<int> m_nbulk;
+//         int m_nsurf;
+//         vector_fp m_mult;
+//         vector<bool> m_do_surf_species;
+//         vector_fp m_fixed_cov;
+//     };
+
 #endif
diff --git a/Cantera/src/oneD/Makefile.in b/Cantera/src/oneD/Makefile.in
index 6e7b1ab48..a6a4159f5 100644
--- a/Cantera/src/oneD/Makefile.in
+++ b/Cantera/src/oneD/Makefile.in
@@ -16,7 +16,7 @@ OBJDIR = .
 CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
 
 # stirred reactors
-OBJS    = MultiJac.o MultiNewton.o newton_utils.o OneDim.o StFlow.o
+OBJS    = MultiJac.o MultiNewton.o newton_utils.o OneDim.o StFlow.o boundaries1D.o refine.o Sim1D.o 
 
 CXX_INCLUDES = -I..
 ONED_LIB = @buildlib@/liboneD.a
diff --git a/Cantera/src/oneD/MultiNewton.cpp b/Cantera/src/oneD/MultiNewton.cpp
index 8cbf7106b..7e89da014 100644
--- a/Cantera/src/oneD/MultiNewton.cpp
+++ b/Cantera/src/oneD/MultiNewton.cpp
@@ -176,7 +176,7 @@ namespace Cantera {
             writelog("\n\nDamped Newton iteration:\n");
             writelog(dashedline);
 
-            sprintf(m_buf,"\n%s  %8s   %8s   %8s  %8s  %8s  %5s\n",
+            sprintf(m_buf,"\n%s  %9s   %9s   %9s  %9s  %9s  %5s\n",
                 "m","F_damp","F_bound","log10(ss)",
                 "log10(s0)","log10(s1)","N_jac");
             writelog(m_buf);
@@ -226,7 +226,7 @@ namespace Cantera {
             // write log information
             if (loglevel > 0) {
                 doublereal ss = r.ssnorm(x1,step1);
-                sprintf(m_buf,"\n%d  %8.4f   %8.4f   %8.4f   %8.4f   %8.4f  %5d ",
+                sprintf(m_buf,"\n%d  %9.5f   %9.5f   %9.5f   %9.5f   %9.5f  %5d ",
                     m,damp,fbound,log10(ss+SmallNumber),
                     log10(s0+SmallNumber),
                     log10(s1+SmallNumber), 
diff --git a/Cantera/src/oneD/OneDim.cpp b/Cantera/src/oneD/OneDim.cpp
index 918c9d80a..1a324a129 100644
--- a/Cantera/src/oneD/OneDim.cpp
+++ b/Cantera/src/oneD/OneDim.cpp
@@ -17,15 +17,16 @@ namespace Cantera {
      * Default constructor. Create an empty object.
      */
     OneDim::OneDim() 
-      : m_jac(0), m_newt(0), 
-	m_rdt(0.0), m_jac_ok(false),
-	m_nd(0), m_bw(0), m_size(0),
-	m_nflow(0), m_nbdry(0), m_init(false),
-	m_ss_jac_age(10), m_ts_jac_age(20),
-	m_nevals(0), m_evaltime(0.0)
+        : m_tmin(1.0e-16), m_tmax(0.1), m_tfactor(0.5),
+          m_jac(0), m_newt(0), 
+          m_rdt(0.0), m_jac_ok(false),
+          m_nd(0), m_bw(0), m_size(0),
+          m_init(false),
+          m_ss_jac_age(10), m_ts_jac_age(20),
+          m_nevals(0), m_evaltime(0.0)
     {
         m_newt = new MultiNewton(1);
-        m_solve_time = 0.0;
+        //m_solve_time = 0.0;
     }
 
 
@@ -33,11 +34,12 @@ namespace Cantera {
      * Construct a OneDim container for the domains pointed at by the
      * input vector of pointers.
     */ 
-    OneDim::OneDim(vector<Resid1D*> domains) : 
+    OneDim::OneDim(vector<Resid1D*> domains) :
+        m_tmin(1.0e-16), m_tmax(0.1), m_tfactor(0.5),
         m_jac(0), m_newt(0), 
 	m_rdt(0.0), m_jac_ok(false),
 	m_nd(0), m_bw(0), m_size(0),
-	m_nflow(0), m_nbdry(0), m_init(false),
+	m_init(false),
 	m_ss_jac_age(10), m_ts_jac_age(20),
 	m_nevals(0), m_evaltime(0.0)
     {
@@ -46,17 +48,25 @@ namespace Cantera {
         m_newt = new MultiNewton(1);
         int nd = domains.size();
         int i;
-        for (i = 0; i < nd; i++) addDomain(domains[i]);
+        for (i = 0; i < nd; i++) {
+            addDomain(domains[i]);
+        }
         init();
         resize();
     }
 
+
+    /**
+     * Domains are added left-to-right. 
+     */
     void OneDim::addDomain(Resid1D* d) {
 
-        // if not the first domain, link it to the last (rightmost)
-        // current domain
+        // if 'd' is not the first domain, link it to the last domain
+        // added (the rightmost one)
         int n = m_dom.size();
-        if (m_dom.size() > 0) m_dom.back()->append(d);
+        if (n > 0) m_dom.back()->append(d);
+
+        // every other domain is a connector
         if (2*(n/2) == n)
             m_connect.push_back(d);
         else
@@ -93,6 +103,18 @@ namespace Cantera {
         }
     }
 
+    /**
+     * Save statistics on function and Jacobiab evaulation, 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) {
             int nev = m_jac->nEvals();
@@ -108,16 +130,18 @@ namespace Cantera {
         }
     }
 
+
+    /**
+     * Call after one or more grids has been refined.
+     */
     void OneDim::resize() {
         int i;
         m_bw = 0;
-        m_nvars.clear();
-        m_loc.clear();
-
+        vector_int nvars, loc;
         int lc = 0;
 
+        // save the statistics for the last grid
         saveStats();
-
         m_pts = 0;
         for (i = 0; i < m_nd; i++) {
             Resid1D* d = m_dom[i];
@@ -125,8 +149,8 @@ namespace Cantera {
             int np = d->nPoints();
             int nv = d->nComponents();
             for (int n = 0; n < np; n++) {
-                m_nvars.push_back(nv);
-                m_loc.push_back(lc);
+                nvars.push_back(nv);
+                loc.push_back(lc);
                 lc += nv;
                 m_pts++;
             }
@@ -147,9 +171,13 @@ namespace Cantera {
 
             m_size = d->loc() + d->size();
         }
+        m_nvars = nvars;
+        m_loc = loc;
+
         m_newt->resize(size());
         m_mask.resize(size());
-        
+
+        // delete the current Jacobian evaluator and create a new one        
         delete m_jac;
         m_jac = new MultiJac(*this);
         m_jac_ok = false;
@@ -160,20 +188,13 @@ namespace Cantera {
 
 
     int OneDim::solve(doublereal* x, doublereal* xnew, int loglevel) {
-        clock_t t0 = clock();
-        static int iok = 0, inotok = 0;
         if (!m_jac_ok) {
-            eval(-1, x, xnew, 0.0, 0); // m_rdt);
-            m_jac->eval(x, xnew, 0.0); // m_rdt);
+            eval(-1, x, xnew, 0.0, 0);
+            m_jac->eval(x, xnew, 0.0);
             m_jac->updateTransient(m_rdt, m_mask.begin());
             m_jac_ok = true;
-            inotok++;
         }
-        else
-            iok++;
         int m = m_newt->solve(x, xnew, *this, *m_jac, loglevel);
-        clock_t t1 = clock();
-        m_solve_time += (t1 - t0)/(1.0*CLOCKS_PER_SEC);
         return m;
     }
 
@@ -202,16 +223,21 @@ namespace Cantera {
      */
     void OneDim::eval(int j, double* x, double* r, doublereal rdt, int count) {
         clock_t t0 = clock();
-        fill(r, r+m_size, 0.0);
+        fill(r, r + m_size, 0.0);
         fill(m_mask.begin(), m_mask.end(), 0);
         if (rdt < 0.0) rdt = m_rdt;
 
         vector<Resid1D*>::iterator d; 
+
+        // iterate over the bulk domains first
         for (d = m_bulk.begin(); d != m_bulk.end(); ++d)
             (*d)->eval(j, x, r, m_mask.begin(), rdt);
+
+        // then over the connector domains
         for (d = m_connect.begin(); d != m_connect.end(); ++d)
             (*d)->eval(j, x, r, m_mask.begin(), rdt);
 
+        // increment counter and time
         if (count) {
             clock_t t1 = clock();
             m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC;
@@ -220,7 +246,7 @@ namespace Cantera {
     }
 
 
-    /** 
+    /**
      * The 'infinity' (maximum magnitude) norm of the steady-state
      * residual. Used only for diagnostic output.
      */
@@ -240,9 +266,15 @@ namespace Cantera {
     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 (fabs(rdt_old - m_rdt) > Tiny) {
             m_jac->updateTransient(m_rdt, m_mask.begin());
         }
+
+        // iterate over all domains, preparing each one to begin
+        // time stepping
         Resid1D* d = left();
         while (d) {
             d->initTimeInteg(dt, x);
@@ -250,7 +282,8 @@ namespace Cantera {
         }
     }
     
-    /** 
+
+    /**
      * 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.
@@ -276,6 +309,7 @@ namespace Cantera {
         m_init = true;
     }
 
+
     /**
      * Signal that the current Jacobian is no longer valid.
      */
@@ -284,6 +318,7 @@ namespace Cantera {
             m_container->jacobian().setAge(10000);
     }
 
+
     /**
      * Take time steps using Backward Euler.
      *
@@ -294,13 +329,15 @@ namespace Cantera {
      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);
                 
          if (loglevel > 0) {
-             writelog("Begin time integration.\n\n");
-             writelog(" step    size (s)    log10(ss) ");
-             writelog("===============================");
+             //writelog("Begin time stepping.\n\n");
+             writelog("\n\n step    size (s)    log10(ss) \n");
+             writelog("===============================\n");
          }
+
          int n = 0, m;
          doublereal ss;
          char str[80];
@@ -310,32 +347,54 @@ namespace Cantera {
                  sprintf(str, " %4d  %10.4g  %10.4g" , n,dt,log10(ss));
                  writelog(str);
              }
+
+             // set up for time stepping with stepsize dt
              initTimeInteg(dt,x);
+
+             // solve the transient problem
              m = solve(x, r, loglevel-1);
+
+             // successful time step. Copy the new solution in r to 
+             // the current solution in x.
              if (m >= 0) {
                  n += 1;
                  if (loglevel > 0) writelog("\n");
                  copy(r, r + m_size, x);
+                 if (m == 100) {
+                     dt *= 1.5;
+                     cout << "m = 100, dt = " << dt << endl;
+                 }
+                 if (dt > m_tmax) dt = m_tmax;
              }
+
+             // No solution could be found with this time step. 
+             // Decrease the stepsize and try again.
              else {
-                 if (loglevel > 0) writelog("...failure.");
-                 dt *= 0.5;
-                 if (dt < 1.e-14)
+                 if (loglevel > 0) writelog("...failure.\n");
+                 dt *= m_tfactor;
+                 cout << "halved dt = " << dt << endl;
+                 if (dt < m_tmin)
                      throw CanteraError("OneDim::timeStep",
                          "Time integration failed.");
              }
          }
+
+         // Prepare to solve the steady problem.
          setSteadyMode();
-         newton().setOptions(m_ss_jac_age);        
+         newton().setOptions(m_ss_jac_age);     
+
+         // return the value of the last stepsize, which may be smaller
+         // than the initial stepsize   
          return dt;
      }
 
+
     void OneDim::save(string fname, string id, string desc, doublereal* sol) {
 
         struct tm *newtime;
         time_t aclock;
-        ::time( &aclock );                  /* Get time in seconds */
-        newtime = localtime( &aclock );   /* Convert time to struct tm form */
+        ::time( &aclock );              /* Get time in seconds */
+        newtime = localtime( &aclock ); /* Convert time to struct tm form */
 
         XML_Node root("doc");
         ifstream fin(fname.c_str());
diff --git a/Cantera/src/oneD/OneDim.h b/Cantera/src/oneD/OneDim.h
index 9685ed573..7b2ff0849 100644
--- a/Cantera/src/oneD/OneDim.h
+++ b/Cantera/src/oneD/OneDim.h
@@ -17,6 +17,7 @@ namespace Cantera {
      * represented by an instance of Resid1D.
      */
     class OneDim {
+
     public: 
 
         // Default constructor.
@@ -31,14 +32,19 @@ namespace Cantera {
         /// Add a domain.
         void addDomain(Resid1D* d);
 
-        /// Return a reference to the Jacobian.
+        /// Return a reference to the Jacobian evaluator.
         MultiJac& jacobian();
 
         /// Return a reference to the Newton iterator.
         MultiNewton& newton();
 
-        /// Solve F(x) = 0, where F(x) is the multi-domain residual.
-        int solve(doublereal* x, doublereal* xnew, int loglevel);
+        /**
+         * Solve F(x) = 0, where F(x) is the multi-domain residual function.
+         * @param x0         Starting estimate of solution.
+         * @param x1         Final solution satisfying F(x1) = 0.
+         * @param loglevel   Controls amount of diagnostic output.
+         */
+        int solve(doublereal* x0, doublereal* x1, int loglevel);
 
         /// Number of domains.
         int nDomains() const { return m_nd; }
@@ -61,8 +67,10 @@ namespace Cantera {
         /// Number of solution components at global point jg.
         int nVars(int jg) { 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.
+         */
         int loc(int jg) { return m_loc[jg]; }
 
         /// Jacobian bandwidth.
@@ -74,13 +82,16 @@ namespace Cantera {
         /// Total number of points.
         int points() { return m_pts; }
 
-        /// Staedy-state max norm of the residual.
+        /**
+         * Steady-state max norm of the residual evaluated using solution x.
+         * On return, array r contains the steady-state residual values.
+         */
         doublereal ssnorm(doublereal* x, doublereal* r);
 
         /// Reciprocal of the time step. 
         doublereal rdt() const { return m_rdt; }
         
-        /// Prepare for time stepping.
+        /// Prepare for time stepping beginning with solution x.
         void initTimeInteg(doublereal dt, doublereal* x);
 
         /// True if transient mode.
@@ -89,10 +100,25 @@ namespace Cantera {
         /// True if steady mode.
         bool steady() const { return (m_rdt == 0.0); }
 
-        /// Set steady mode
+
+        /**
+         * Set steady mode.  After invoking this method, subsequent
+         * calls to solve() will solve the steady-state problem.
+         */
         void setSteadyMode();
 
-        /// Evaluate the multi-domain residual function 
+
+        /**
+         * Evaluate the multi-domain residual function
+         * 
+         * @param j       if j > 0, only evaluate residual for points j-1, j, 
+         *                and j + 1; otherwise, evaluate at all grid points.
+         * @param x       solution vector
+         * @param r       on return, contains the residual vector
+         * @param rdt     Reciprocal of the time step. if omitted, then
+         *                  the default value is used.
+         * @param count   Set to zero to omit this call from the statistics
+         */ 
         void eval(int j, double* x, double* r, doublereal rdt=-1.0, 
             int count = 1);
 
@@ -100,7 +126,8 @@ namespace Cantera {
         Resid1D* pointDomain(int i);
 
         void resize();
-        doublereal solveTime() { return m_solve_time; }
+
+        //doublereal solveTime() { return m_solve_time; }
 
         //void setTransientMask();
         vector_int& transientMask() { return m_mask; }
@@ -113,21 +140,35 @@ namespace Cantera {
 
         void save(string fname, string id, string desc, doublereal* sol);
 
+        // options
+        void setMinTimeStep(doublereal tmin) { m_tmin = tmin; }
+        void setMaxTimeStep(doublereal tmax) { m_tmax = tmax; }
+        void setTimeStepFactor(doublereal tfactor) { m_tfactor = tfactor; }
+        void setJacAge(int ss_age, int ts_age=-1) {
+            m_ss_jac_age = ss_age;
+            if (ts_age > 0) 
+                m_ts_jac_age = ts_age;
+            else
+                m_ts_jac_age = m_ss_jac_age;
+        }
+
     protected:
 
-        MultiJac* m_jac;
-        MultiNewton* m_newt;
-        doublereal m_rdt;
-        bool m_jac_ok;
-        int m_nd;
-        int m_bw;
-        int m_size;
-        vector_int m_states;
-        vector_int m_start;
-        vector_int m_comp, m_points;
+        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 )
+
+        MultiJac* m_jac;          // Jacobian evaluator
+        MultiNewton* m_newt;      // Newton iterator
+        doublereal m_rdt;         // reciprocal of time step
+        bool m_jac_ok;            // if true, Jacobian is current
+        int m_nd;                 // number of domains
+        int m_bw;                 // Jacobian bandwidth
+        int m_size;               // solution vector size
+        
         vector<Resid1D*> m_dom, m_connect, m_bulk;
-        vector_int m_flow, m_bdry;
-        int m_nflow, m_nbdry;
+
         bool m_init;
         vector_int m_nvars;
         vector_int m_loc;
@@ -135,9 +176,10 @@ namespace Cantera {
         int m_pts;
         doublereal m_solve_time;
 
+        // options
         int m_ss_jac_age, m_ts_jac_age;
 
-        // stats
+        // statistics
         int m_nevals;
         doublereal m_evaltime;
         vector_int m_gridpts;
@@ -147,6 +189,7 @@ namespace Cantera {
         vector_fp m_funcElapsed;
 
     private:
+
     };
 
 }
diff --git a/Cantera/src/oneD/Resid1D.h b/Cantera/src/oneD/Resid1D.h
index c9d9b0b24..468f65c01 100644
--- a/Cantera/src/oneD/Resid1D.h
+++ b/Cantera/src/oneD/Resid1D.h
@@ -13,9 +13,10 @@
 #define CT_RESID1D_H
 
 
-//#include "stringUtils.h"
 #include "../ctexceptions.h"
 #include "../xml.h"
+#include "refine.h"
+
 
 namespace Cantera {
 
@@ -32,16 +33,15 @@ namespace Cantera {
 
 
     /**
-     * Base class for single-domain, one-dimensional residual function
-     * evaluators.
+     * Base class for one-dimensional domains.
      */
     class Resid1D {
     public:
 
         /**
          * Constructor.
-         * @param nv Number of variables at each grid point.
-         * @param points Number of grid points.
+         * @param nv      Number of variables at each grid point.
+         * @param points  Number of grid points.
          */
         Resid1D(int nv=1, int points=1, 
             doublereal time = 0.0) : 
@@ -52,17 +52,26 @@ namespace Cantera {
             m_iloc(0),
             m_jstart(0),
             m_left(0), 
-            m_right(0) {
+            m_right(0),
+            m_refiner(0) {
             resize(nv, points);
         }
 
-        /// Destructor.
-        virtual ~Resid1D(){}
+        /// Destructor. Does nothing
+        virtual ~Resid1D(){ delete m_refiner; }
 
         /// Domain type flag.
         const int domainType() { return m_type; }
 
-        const OneDim& container() const{ return *m_container; }
+        /**
+         * The left-to-right location of this domain.
+         */
+        const int domainIndex() { return m_index; }
+
+        /**
+         * The container holding this domain.
+         */
+        const OneDim& container() const { return *m_container; }
 
         /**
          * Specify the container object for this domain, and the
@@ -73,25 +82,40 @@ namespace Cantera {
             m_index = index;
         }
 
-        /** Initialize. Base class method does nothing, but may be
-         * overloaded.
+        /** 
+         * Initialize. Base class method does nothing, but may be
+         * overloaded. 
          */
         virtual void init(){}
 
+        virtual void setInitialState(doublereal* xlocal = 0){}
+        virtual void setState(int 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.
          */
         virtual void resize(int nv, int np) {
+            if (nv != m_nv || !m_refiner) {
+                m_nv = nv;
+                delete m_refiner;
+                m_refiner = new Refiner(*this);
+                cout << "created refiner with nv = " << m_nv << endl;
+            }
             m_nv = nv;
             m_max.resize(m_nv, 0.0);
             m_min.resize(m_nv, 0.0);
             m_rtol.resize(m_nv, 0.0);
             m_atol.resize(m_nv, 0.0);
             m_points = np;
+            m_z.resize(np, 0.0);
             m_slast.resize(m_nv * m_points, 0.0);
             locate();
         }
 
+        Refiner& refiner() { return *m_refiner; }
+
         /// Number of components at each grid point.
         int nComponents() const { return m_nv; }
 
@@ -107,40 +131,69 @@ namespace Cantera {
          */
         void setBounds(int nl, const doublereal* lower, 
             int nu, const doublereal* upper) {
-            if (nl != m_nv || nu != m_nv)
+            if (nl < m_nv || nu < m_nv)
                 throw CanteraError("Resid1D::setBounds",
-                    "wrong array size for solution bounds");
+                    "wrong array size for solution bounds. "
+                    "Size should be at least "+int2str(m_nv));
             copy(upper, upper + m_nv, m_max.begin());
             copy(lower, lower + m_nv, m_min.begin());
         }
 
         void setTolerances(int nr, const doublereal* rtol, 
             int na, const doublereal* atol) {
-            if (nr != m_nv || na != m_nv)
+            if (nr < m_nv || na < m_nv)
                 throw CanteraError("Resid1D::setTolerances",
-                    "wrong array size for solution error tolerances. Size should be "+int2str(m_nv));
+                    "wrong array size for solution error tolerances. "
+                    "Size should be at least "+int2str(m_nv));
             copy(rtol, rtol + m_nv, m_rtol.begin());
             copy(atol, atol + m_nv, m_atol.begin());
         }
 
+        /// Relative tolerance of the nth component.
         doublereal rtol(int n) { return m_rtol[n]; }
+
+        /// Absolute tolerance of the nth component.
         doublereal atol(int n) { return m_atol[n]; }
 
+        /// Upper bound on the nth component.
         doublereal upperBound(int n) const { return m_max[n]; }
+
+        /// Lower bound on the nth component
         doublereal lowerBound(int 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.
+         */
         void initTimeInteg(doublereal dt, const doublereal* x0) {
             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
+         */
         void setSteadyMode() { m_rdt = 0.0; }
 
+        /// True if in steady-state mode
         bool steady() { return (m_rdt == 0.0); }
+
+        /// 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.
+         */
         void needJacUpdate();
 
+        /**
+         * Evaluate the steady-state residual at all points, even if in 
+         * transient mode. Used to print diagnostic output.
+         */
         void evalss(doublereal* x, doublereal* r, integer* mask) {
             eval(-1,x,r,mask,0.0);
         }
@@ -155,59 +208,163 @@ namespace Cantera {
                 "residual function not defined.");
         }
 
+        /**
+         * Does nothing.
+         */
         virtual void update(doublereal* x) {}
 
-        doublereal time() { return m_time;}
+        doublereal time() const { return m_time;}
         void incrementTime(doublereal dt) { m_time += dt; }
         size_t index(int n, int j) const { return m_nv*j + n; }
+        doublereal value(doublereal* x, int n, int j) const {
+            return x[index(n,j)];
+        }
 
         virtual void setJac(MultiJac* jac){}
         virtual void save(XML_Node& o, doublereal* sol) {
             throw CanteraError("Resid1D::save","base class method called");
         }
 
-        int size() { return m_nv*m_points; }
+        int size() const { return m_nv*m_points; }
 
+        /**
+         * Find the index of the first grid point in this domain, and
+         * the start of its variables in the global solution vector.
+         */
         void locate() {
+
             if (m_left) {
+                // there is a domain on the left, so the first grid point
+                // in this domain is one more than the last one on the left
                 m_jstart = m_left->lastPoint() + 1;
+
+                // the starting location in the solution vector
                 m_iloc = m_left->loc() + m_left->size();
             }
             else {
+                // this is the left-most domain
                 m_jstart = 0;
                 m_iloc = 0;
             }
+            // if there is a domain to the right of this one, then 
+            // repeat this for it
             if (m_right) m_right->locate();
         }
 
-        virtual int loc(int j = 0) { return m_iloc; }
+        /**
+         * Location of the start of the local solution vector in the global
+         * solution vector,
+         */
+        virtual int loc(int j = 0) const { return m_iloc; }
 
-        int firstPoint() { return m_jstart; }
-        int lastPoint() { return m_jstart + m_points - 1; }
+        /**
+         * The index of the first (i.e., left-most) grid point
+         * belonging to this domain.
+         */
+        int firstPoint() const { return m_jstart; }
 
+        /**
+         * The index of the last (i.e., right-most) grid point
+         * belonging to this domain.
+         */
+        int lastPoint() const { return m_jstart + m_points - 1; }
+
+        /**
+         * 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(Resid1D* left) { 
+            m_left = left;
+            locate();
+        }
+
+        /**
+         * Set the right neighbor to domain 'right.' 
+         */
+        void linkRight(Resid1D* right) { m_right = right; }
+
+        /**
+         * Append domain 'right' to this one, and update all links.
+         */
         void append(Resid1D* right) {
             linkRight(right);
             right->linkLeft(this);
         }
 
-        void linkLeft(Resid1D* left) { 
-            m_left = left;
-            locate();
-        }
-        void linkRight(Resid1D* right) { m_right = right; }
+        /**
+         * Return a pointer to the left neighbor.
+         */
+        Resid1D* left() const { return m_left; }
 
-        Resid1D* left() { return m_left; }
-        Resid1D* right() { return m_right; }
+        /**
+         * Return a pointer to the right neighbor.
+         */
+        Resid1D* right() const { return m_right; }
 
-        double prevSoln(int n, int j) const{
+        /**
+         * Value of component n at point j in the previous solution.
+         */
+        double prevSoln(int n, int j) const {
             return m_slast[m_nv*j + n];
         }
         
+        /**
+         * Specify an identifying tag for this domain.
+         */
         void setID(const string& s) {m_id = s;}
+
+        /**
+         * Specify descriptive text for this domain.
+         */
         void setDesc(const string& s) {m_desc = s;}
 
         virtual void getTransientMask(integer* mask){}
 
+        virtual void showSolution(ostream& s, const doublereal* x) {}
+
+        doublereal z(int jlocal) const {
+            return m_z[jlocal];
+        }
+        doublereal zmin() const { return m_z[0]; }
+        doublereal zmax() const { return m_z[m_points - 1]; }
+
+
+        void setProfile(string name, doublereal* values, doublereal* soln) {
+            int n, j;
+            for (n = 0; n < m_nv; n++) {
+                if (name == componentName(n)) {
+                    for (j = 0; j < m_points; j++) {
+                        soln[index(n, j) + m_iloc] = values[j];
+                    }
+                    return;
+                }
+            }
+            throw CanteraError("Resid1D::setProfile",
+                "unknown component: "+name);
+        }
+
+        vector_fp& grid() { return m_z; }
+        const vector_fp& grid() const { return m_z; }
+        doublereal grid(int point) { return m_z[point]; }
+
+        virtual void setupGrid(int n, const doublereal* z) {}
+
+        /**
+         * Writes some or all initial solution values into array x,
+         * which is the solution vector for this domain. This allows
+         * initial values that have been set prior to installing this
+         * domain into the container to be written to the global
+         * solution vector.
+         */
+        virtual void _getInitialSoln(doublereal* x) {cout << "base class method _getInitialSoln called!" << endl;}
+
+        /**
+         * Perform any necessary domain-specific initialization using 
+         * local solution vector x.
+         */ 
+        virtual void _finalize(const doublereal* x) {cout << "base class method _finalize called!" << endl;}
+
     protected:
 
         doublereal m_rdt;
@@ -219,6 +376,7 @@ namespace Cantera {
         vector_fp m_min;
         vector_fp m_rtol;
         vector_fp m_atol;
+        vector_fp m_z;
         OneDim* m_container;
         int m_index;
         int m_type;
@@ -226,6 +384,7 @@ namespace Cantera {
         int m_jstart;
         Resid1D *m_left, *m_right;
         string m_id, m_desc;
+        Refiner* m_refiner;
 
     private:
 
diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp
new file mode 100644
index 000000000..e7ffc7fbb
--- /dev/null
+++ b/Cantera/src/oneD/Sim1D.cpp
@@ -0,0 +1,306 @@
+/**
+ * @file Sim1D.cpp
+ */
+
+#include "Sim1D.h"
+
+namespace Cantera {
+
+
+    Sim1D::Sim1D(vector<Resid1D*>& domains) : OneDim(domains) {
+
+        // 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);
+        for (int n = 0; n < m_nd; n++) {
+            domain(n)._getInitialSoln(m_x.begin() + start(n));
+        }
+
+        // set some defaults
+        m_tstep = 1.0e-5;
+        //m_maxtimestep = 10.0;
+        m_steps.push_back(1);
+        m_steps.push_back(2);
+        m_steps.push_back(5);
+        m_steps.push_back(10);
+
+    }
+
+
+    /**
+     * 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(int dom, int comp, int localPoint,  doublereal value) {
+        int iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
+        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(int dom, int comp, int localPoint) const {
+        int iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
+        return m_x[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(int dom, int comp, 
+        const vector_fp& pos, const vector_fp& values) {
+
+        Resid1D& d = domain(dom);
+        int np = d.nPoints();
+        int n;
+        doublereal z0 = d.zmin();
+        doublereal z1 = d.zmax();
+        doublereal zpt, frac, v;
+        for (n = 0; n < np; n++) {
+            zpt = d.z(n);
+            frac = (zpt - z0)/(z1 - z0);
+            v = linearInterp(frac, pos, values);
+            setValue(dom, comp, n, v);
+        }
+    }
+
+
+
+    void Sim1D::setFlatProfile(int dom, int comp, doublereal v) {
+        int np = domain(dom).nPoints();
+        int n;
+        for (n = 0; n < np; n++) { setValue(dom, comp, n, v); }
+    }
+
+
+    void Sim1D::showSolution(ostream& s) {
+        for (int n = 0; n < m_nd; n++) {
+            domain(n).showSolution(s, m_x.begin() + start(n));
+        }
+    }
+
+
+    void Sim1D::finalize() {
+        for (int n = 0; n < m_nd; n++) {
+            domain(n)._finalize(m_x.begin() + start(n));
+        }
+    }
+
+
+    void Sim1D::setTimeStep(doublereal stepsize, int n, integer* tsteps) {
+        m_tstep = stepsize;
+        m_steps.resize(n);
+        for (int i = 0; i < n; i++) m_steps[i] = tsteps[i];
+    }
+
+
+    void Sim1D::newtonSolve(int loglevel) {
+        int m = OneDim::solve(m_x.begin(), m_xnew.begin(), loglevel);
+        if (m >= 0) 
+            copy(m_xnew.begin(), m_xnew.end(), m_x.begin());
+        else if (m > -10)
+            throw CanteraError("Sim1D::newtonSolve","no solution found");
+        else {
+            cout << "ERROR: solve returned m = " << m << endl;
+            exit(-1);
+        }
+    }
+
+
+    void Sim1D::solve(int loglevel, bool refine_grid) {
+
+        int new_points = 1;
+        int istep, nsteps;
+        doublereal dt = m_tstep;
+        int soln_number = -1;
+
+        finalize();
+
+        while (new_points > 0) {
+
+            istep = 0;
+            nsteps = m_steps[istep];
+
+            bool ok = false;
+            while (!ok) {
+                    
+                try {
+                    if (loglevel > 0) {
+                        writelog("Attempt Newton solution of steady-state problem...");
+                    }
+                    newtonSolve(loglevel-1);
+                    
+                    if (loglevel > 0) {
+                        writelog("success.\n\n");
+                        //writelog("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n"); 
+                        writelog("Problem solved on [");
+                        for (int mm = 1; mm < nDomains(); mm+=2) {
+                            writelog(int2str(domain(mm).nPoints()));
+                            if (mm < nDomains() - 2) writelog(", ");
+                        }
+                        writelog("]");
+                        writelog(" point grid(s).\n\n");
+                        //writelog("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"); 
+                    }
+                    ok = true;
+                    soln_number++;
+                    
+                }
+
+                catch (CanteraError) {
+
+                    char buf[100];
+                    if (loglevel > 0) writelog("failure. \n\n");
+                    if (loglevel == 1) writelog("Take "+int2str(nsteps)+" timesteps   ");
+                    dt = timeStep(nsteps, dt, m_x.begin(), m_xnew.begin(), loglevel-1);
+                    if (loglevel == 1) {
+                        sprintf(buf, " %10.4g %10.4g \n", dt, 
+                            log10(ssnorm(m_x.begin(), m_xnew.begin())));
+                        writelog(buf);
+                    }
+                    istep++;
+                    if (istep >= int(m_steps.size())) {
+                        nsteps = m_steps.back();
+                        dt *= 2.0;
+                        cout << " doubled dt = " << dt << endl;
+                    }
+                    else {
+                        nsteps = m_steps[istep];
+                    }
+                    if (dt > m_tmax) dt = m_tmax;
+                }
+            }
+            if (loglevel > 2) showSolution(cout);
+            
+            if (refine_grid) {
+                new_points = refine(loglevel);
+            }
+            else {
+                new_points = 0;
+            }
+        }
+    }
+
+    /**
+     * Refine the grid in all domains.
+     */
+    int Sim1D::refine(int loglevel) {
+        int np = 0;
+        vector_fp znew, xnew;
+        doublereal xmid, zmid;
+        int strt, n, m, i;
+
+        for (n = 0; n < m_nd; n++) {
+            strt = znew.size();
+            Resid1D& d = domain(n);
+            Refiner& r = d.refiner();
+
+            // determine where new points are needed
+            r.analyze(d.grid().size(), d.grid().begin(), m_x.begin() + start(n));
+            if (loglevel > 0) { r.show(); }
+
+            np += r.nNewPoints();
+            int comp = d.nComponents();
+
+            // loop over points in the current grid
+            int npnow = d.nPoints();
+            for (m = 0; m < npnow; m++) {
+
+                // add the current grid point to the new grid
+                znew.push_back(d.grid(m));
+
+                // do the same for the solution at this point
+                for (i = 0; i < comp; i++) {
+                    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
+
+                if (r.newPointNeeded(m)) {
+
+                    // add new point at midpoint
+                    zmid = 0.5*(d.grid(m) + d.grid(m+1));
+                    znew.push_back(zmid);
+
+                    // for each component, linearly interpolate the solution to 
+                    // this point
+                    for (i = 0; i < comp; i++) {
+                        xmid = 0.5*(value(n, i, m) + value(n, i, m+1));
+                        xnew.push_back(xmid);
+                    }
+                }
+            }
+        }
+
+        // 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.
+
+        int gridstart = 0, gridsize;
+        for (n = 0; n < m_nd; n++) {
+            Resid1D& d = domain(n);
+            Refiner& r = d.refiner();
+            gridsize = d.nPoints() + r.nNewPoints();
+            d.setupGrid(gridsize, znew.begin() + gridstart);
+            gridstart += gridsize;
+        }
+
+        // Replace the current solution vector with the new one
+        m_x.resize(xnew.size());
+        copy(xnew.begin(), xnew.end(), m_x.begin());
+
+        // resize the work array
+        m_xnew.resize(xnew.size());
+
+        //        copy(xnew.begin(), xnew.end(), m_xnew.begin());
+
+        resize();
+        finalize();
+        return np;
+    }
+
+
+    /**
+     * 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) {
+        if (dom >= 0) {
+            Refiner& r = domain(dom).refiner();
+            r.setCriteria(ratio, slope, curve);
+        }
+        else {
+            for (int n = 0; n < m_nd; n++) {
+                Refiner& r = domain(n).refiner();
+                r.setCriteria(ratio, slope, curve);
+            }                    
+        }
+    }
+            
+}
diff --git a/Cantera/src/oneD/Sim1D.h b/Cantera/src/oneD/Sim1D.h
new file mode 100644
index 000000000..7760b6fe1
--- /dev/null
+++ b/Cantera/src/oneD/Sim1D.h
@@ -0,0 +1,89 @@
+/**
+ * @file Sim1D.h
+ */
+
+#ifndef CT_SIM1D_H
+#define CT_SIM1D_H
+
+#include "OneDim.h"
+#include "../funcs.h"
+
+namespace Cantera {
+
+    /**
+     * One-dimensional simulations.
+     */
+    class Sim1D : public OneDim {
+
+    public: 
+
+        /**
+         * Default constructor. This constructor can be used to create
+         * a dummy object if necessary, but is not usually called in
+         * user programs. Use the next constructor instead.
+         */
+        Sim1D();
+
+
+        /**
+         * Standard constructor. 
+         * @param domains A vector of pointers to the domains to be linked together.
+         * The domains must appear in left-to-right order.
+         */ 
+        Sim1D(vector<Resid1D*>& domains);
+
+        /// Destructor. Does nothing.
+        virtual ~Sim1D(){}
+
+
+        /// Set one entry in the solution vector.
+        void setValue(int dom, int comp, int localPoint,  doublereal value);
+
+        /// Get one entry in the solution vector.
+        doublereal value(int dom, int comp, int localPoint) const;
+
+        /// Specify a profile for one component of one domain.
+        void setProfile(int dom, int comp, const vector_fp& pos, 
+            const vector_fp& values);
+
+        /// Set component 'comp' of domain 'dom' to value 'v' at all points.
+        void setFlatProfile(int dom, int comp, doublereal v);
+
+        /// Print to stream s the current solution for all domains.     
+        void showSolution(ostream& s);
+
+        /// Calls method _finalize in each domain.
+        void finalize();
+
+        void setTimeStep(doublereal stepsize, int n, integer* tsteps);
+
+        //void setMaxTimeStep(doublereal tmax) { m_maxtimestep = tmax; }
+
+        void solve(int loglevel = 0, bool refine_grid = true);
+
+        int refine(int loglevel=0);
+
+        void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
+            doublereal slope = 0.8, doublereal curve = 0.8);
+
+    protected:
+
+        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
+
+    private:
+
+        void newtonSolve(int loglevel);
+
+
+    };
+
+}
+
+#endif
+
+
diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp
index d75f46945..75d05d435 100644
--- a/Cantera/src/oneD/StFlow.cpp
+++ b/Cantera/src/oneD/StFlow.cpp
@@ -101,28 +101,6 @@ namespace Cantera {
 
     //--------------------- linear interp ------------------------------
 
-    /**
-     * Linearly interpolate a function defined on a discrete grid.
-     * vector xpts contains a monotonic sequence of grid points, and 
-     * vector fpts contains function values defined at these points.
-     * The value returned is the linear interpolate at point x.
-     * If x is outside the range of xpts, the value of fpts at the 
-     * nearest end is returned.
-     */
-
-    doublereal linearInterp(doublereal x, vector_fp& xpts, vector_fp& fpts) {
-        if (x <= xpts[0]) return fpts[0];
-        if (x >= xpts.back()) return fpts.back();
-        doublereal* loc = lower_bound(xpts.begin(), xpts.end(), x);
-        int iloc = int(loc - xpts.begin()) - 1;
-        doublereal ff = fpts[iloc] + 
-            (x - xpts[iloc])*(fpts[iloc + 1] 
-                - fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
-        return ff;
-    }
-
-
-
     StFlow::StFlow(igthermo_t* ph, int nsp, int points) : 
         Resid1D(nsp+4, points),
         m_inlet_u(0.0),
@@ -146,10 +124,16 @@ namespace Cantera {
 
         m_points = points;
         m_thermo = ph;
-        m_nv = m_nsp + 4;
 
         if (ph == 0) return; // used to create a dummy object
 
+        int nsp2 = m_thermo->nSpecies();
+        if (nsp2 != m_nsp) {
+            m_nsp = nsp2;
+            Resid1D::resize(m_nsp+4, points);
+        }
+
+
         // make a local copy of the species molecular weight vector
         m_wt = m_thermo->molecularWeights();
 
@@ -170,7 +154,9 @@ namespace Cantera {
         m_surfdot.resize(m_nsp, 0.0);
         m_ybar.resize(m_nsp);
 
-        // default solution bounds
+
+        //-------------- default solution bounds --------------------
+
         vector_fp vmin(m_nv), vmax(m_nv);
         
         // no bounds on u
@@ -178,7 +164,7 @@ namespace Cantera {
         vmax[0] = 1.e20;
 
         // no negative V
-        vmin[1] = -0.01;
+        vmin[1] = -0.1;
         vmax[1] = 1.e20;
 
         // temperature bounds
@@ -187,17 +173,30 @@ namespace Cantera {
 
         // lamda should be negative
         vmin[3] = -1.e20;
-        vmax[3] = 0.001;
+        vmax[3] = 1.0;
 
         // mass fraction bounds
         int k;
         for (k = 0; k < m_nsp; k++) {
-            vmin[4+k] = -1.e-5;
+            vmin[4+k] = -1.0e-5;
             vmax[4+k] = 1.1;
         }
         setBounds(vmin.size(), vmin.begin(), vmax.size(), vmax.begin());
+
+
+        //-------------------- default error tolerances ----------------
+        vector_fp rtol(m_nv, 1.0e-8);
+        vector_fp atol(m_nv, 1.0e-15);
+        setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin());
+
+        //-------------------- grid refinement -------------------------
+        m_refiner->setActive(0, false);
+        m_refiner->setActive(1, false);
+        m_refiner->setActive(2, false);
+        m_refiner->setActive(3, false);
     }
 
+
     /**
      * Change the grid size. Called after grid refinement.
      */
@@ -258,6 +257,7 @@ namespace Cantera {
             throw CanteraError("setTransport","unknown transport model.");
     }
 
+
     /**
      * Set the gas object state to be consistent with the solution at
      * point j.
@@ -269,6 +269,11 @@ namespace Cantera {
         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,int j) {
         m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
         const doublereal* yyj = x + m_nv*j + c_offset_Y;
@@ -280,44 +285,26 @@ namespace Cantera {
     }
 
 
-//     /**
-//      * Integrate the species mass fractions at each point separately,
-//      * without the transport terms. This method is provided to
-//      * condition a poor estimate of the solution to produce a better
-//      * starting estimate for Newton iteration. It is not used by any
-//      * other method, but is available for use in user codes, if
-//      * desired.
-//      */
-//     void StFlow::integrateChem(doublereal* x,doublereal dt) {
-//         int j;
-//         if (!ready()) return;
-//         if (m_integrator == 0) {
-//              m_integrator = new ImplicitChem(*m_kin, *m_thermo);
-//             m_integrator->initialize(0.0);
-//         }
-//         for (j = 0; j < m_points; j++) {
-//             setGas(x,j);
-//             m_integrator->integrate(0.0, dt);
-//             m_thermo->getMassFractions(m_nsp, &x[index(c_offset_Y,j)]);
-//             T(x,j) = m_thermo->temperature();
-//         }
-//     }
-
-
 
     /**
-     *  Evaluate the residual function for 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.
+     *  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 AxiStagnFlow::eval(int 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 (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return;
 
+        // if evaluating a Jacobian, compute the steady-state residual
         if (jg >= 0) rdt = 0.0;
 
         // start of local part of global arrays
@@ -328,11 +315,11 @@ namespace Cantera {
         int jmin, jmax, jpt;
         jpt = jg - firstPoint();
 
-        if (jg < 0) {
+        if (jg < 0) {      // evaluate all points
             jmin = 0;
             jmax = m_points - 1;
         }
-        else {
+        else {            // evaluate points for Jacobian
             jmin = max(jpt-1, 0);
             jmax = min(jpt+1,m_points-1);
         }
@@ -349,14 +336,16 @@ namespace Cantera {
         //              update properties
         //-----------------------------------------------------
 
-        // thermodynamic properties
+        // thermodynamic properties only if a Jacobian is
+        // not being evaluated
         if (jpt < 0) updateThermo(x, j0, j1);
 
         // update transport properties only if a Jacobian is
         // not being evaluated
         if (jpt < 0) updateTransport(x, j0, j1);
 
-        // update the species diffusive mass fluxes
+        // update the species diffusive mass fluxes whether or not a
+        // Jacobian is being evaluated
         updateDiffFluxes(x, j0, j1);
 
 
@@ -380,17 +369,28 @@ namespace Cantera {
 
 #define NEW_INLET
 #ifdef NEW_INLET
-                 // continuity
+
+                // 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
                  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);
-                 //cout << "rsd: " << rsd[0] << "  " << rsd[1] << "  " << rsd[2] << endl;
+                 //cout << "density = " << density(0) << "  " << u(x,0)
+                 //      << "   " << rho_u(x,0) << endl;
 
-                 // zero flux
+                 // The default boundary condition for species is zero
+                 // flux. However, the boundary object may modify
+                 // this.
                  for (k = 0; k < m_nsp; k++) {
                      rsd[index(c_offset_Y + k, 0)] =  
                          -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
@@ -421,21 +421,18 @@ namespace Cantera {
 
 
             //----------------------------------------------
-            //         right boundary
             //
-            //  The right boundary residuals are for a nonreacting,
-            //  impermeable wall. Since domains are evaluated left to
-            //  right, the surface object may add terms to these
-            //  residual equations.
+            //         right boundary
             //
             //----------------------------------------------
 
             else if (j == m_points - 1) {
 
-                //m_boundary[1]->eval(x + index(0, j), m_rho[j], 
-                //    m_flux.begin() + m_nsp*(j-1), 
-                //    rsd + index(0, j));
-
+                 // the boundary object connected to the right of this
+                 // one may 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(0,j)] = rho_u(x,j);
                 rsd[index(1,j)] = V(x,j);
                 rsd[index(2,j)] = T(x,j);
@@ -445,10 +442,14 @@ namespace Cantera {
                     sum += Y(x,k,j);
                     rsd[index(k+4,j)] = rho_u(x,j)*Y(x,k,j) + m_flux(k,j-1);
                 }
+
+                // TODO: why is this done here, but not for the left
+                // boundary or interior?
                 rsd[index(4,j)] = 1.0 - sum;
                 diag[index(4,j)] = 0;
             }
 
+
             //------------------------------------------
             //     interior points                 
             //------------------------------------------
@@ -855,7 +856,8 @@ namespace Cantera {
 
 
 
-    void StFlow::outputTEC(ostream &s, const doublereal* x, string title, int zone) {
+    void StFlow::outputTEC(ostream &s, const doublereal* x, 
+        string title, int zone) {
         int j,k;
         s << "TITLE     = \"" + title + "\"" << endl;
         s << "VARIABLES = \"Z (m)\"" << endl;
@@ -884,17 +886,19 @@ namespace Cantera {
 
     string StFlow::componentName(int n) const {
         switch(n) {
-        case 0: return "u [m/s]";
-        case 1: return "V [1/s]";
-        case 2: return "T [K]";
+        case 0: return "u";
+        case 1: return "V";
+        case 2: return "T";
         case 3: return "lambda";
         default:
             if (n >= (int) c_offset_Y && n < (int) (c_offset_Y + m_nsp)) {
-                if (m_do_species[n - c_offset_Y]) 
-                    return m_thermo->speciesName(n - c_offset_Y)+"  ";
-                else
-                    return m_thermo->speciesName(n - c_offset_Y)+" *";
+                    return m_thermo->speciesName(n - c_offset_Y);
             }
+                 //    if (m_do_species[n - c_offset_Y]) 
+                 //   return m_thermo->speciesName(n - c_offset_Y)+"  ";
+                 // else
+                 //   return m_thermo->speciesName(n - c_offset_Y)+" *";
+                 //}
             else 
                 return "<unknown>";
         }
@@ -1016,7 +1020,6 @@ namespace Cantera {
         for (n = 0; n < nd; n++) {
             XML_Node& fa = *d[n];
             nm = fa["title"];
-            cout << "nm = " << nm << endl;
             getFloatArray(fa,x,false);
             if (nm == "u") {
                 writelog("axial velocity   ");
@@ -1134,7 +1137,7 @@ namespace Cantera {
         flow.addAttribute("id",id);
         addString(flow,"timestamp",asctime(newtime));
         addFloat(flow, "pressure", m_press, "Pa", "pressure"); 
-        addString(flow,"solve_time",fp2str(m_container->solveTime()));
+        //        addString(flow,"solve_time",fp2str(m_container->solveTime()));
         if (desc != "") addString(flow,"description",desc);
         XML_Node& gv = flow.addChild("grid_data");
         addFloatArray(gv,"z",m_z.size(),m_z.begin(),
@@ -1220,9 +1223,9 @@ namespace Cantera {
         m_jac = jac;
     }
 
-    void StFlow::requestJacUpdate() {
-        if (m_jac) m_jac->setAge(10000);
-    }
+    //void StFlow::requestJacUpdate() {
+    //    if (m_jac) m_jac->setAge(10000);
+    //}
 
     void StFlow::setEnergyFactor(doublereal efctr) {
         doublereal de = efctr - m_efctr;
diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h
index a1204a1b9..63d5f494d 100644
--- a/Cantera/src/oneD/StFlow.h
+++ b/Cantera/src/oneD/StFlow.h
@@ -15,17 +15,15 @@
 #define CT_STFLOW_H
 
 #include "../transport/TransportBase.h"
-//#include "IdealGasMix.h"
 #include "Resid1D.h"
-//#include "../ChemEquil.h"
 #include "../Array.h"
 #include "../sort.h"
-//#include "ImplicitChem.h"
 #include "../IdealGasPhase.h"
 #include "../Kinetics.h"
-
+#include "../funcs.h"
 #include "../flowBoundaries.h"
 
+
 namespace Cantera {
 
     typedef IdealGasPhase igthermo_t;
@@ -82,19 +80,18 @@ namespace Cantera {
          */
         //@{
 
-        void setupGrid(int n, const doublereal* z);
+        virtual void setupGrid(int n, const doublereal* z);
 
         //thermo_t& phase() { return *m_phase; }
         thermo_t& phase() { return *m_thermo; }
         kinetics_t& kinetics() { return *m_kin; }
 
         /**
-         * Set the thermo manager. Note that the flow equations
-         * assume the ideal gas equation.
+         * Set the thermo manager. Note that the flow equations assume
+         * the ideal gas equation.
          */
         void setThermo(igthermo_t& th) { 
             m_thermo = &th;
-            //m_phase = &th.phase();
         }
 
         /// set the kinetics manager
@@ -104,13 +101,56 @@ namespace Cantera {
         void setTransport(Transport& trans, bool withSoret = false);
 
         /// set the pressure
-        void setPressure(doublereal p) {
-            m_press = p;
-        }
+        void setPressure(doublereal p) { m_press = p; }
 
         /// Check that all required parameters have been set.
         bool ready();
 
+        virtual void setState(int point, const doublereal* state) {
+            setTemperature(point, state[2]);
+            int k;
+            for (k = 0; k < m_nsp; k++) {
+                setMassFraction(point, k, state[4+k]);
+            }
+        }
+
+
+        virtual void _getInitialSoln(doublereal* x) {
+            int k, j;
+            for (j = 0; j < m_points; j++) {
+                x[index(2,j)] = T_fixed(j);
+                for (k = 0; k < m_nsp; k++) {
+                    x[index(4+k,j)] = Y_fixed(k,j);
+                }
+            }
+        }   
+
+        virtual void _finalize(const doublereal* x) {
+            int k, j;
+            doublereal zz, tt;
+            int nz = m_zfix.size();
+            bool e = m_do_energy[0];
+            for (j = 0; j < m_points; j++) {
+                if (e || nz == 0) 
+                    setTemperature(j, T(x, j));
+                else {
+                    zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
+                    tt = linearInterp(zz, m_zfix, m_tfix);
+                    setTemperature(j, tt);
+                }   
+                for (k = 0; k < m_nsp; k++) {
+                    setMassFraction(j, k, Y(x, k, j));
+                }
+            }
+            if (e) solveEnergyEqn();
+        }
+
+
+        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
@@ -128,10 +168,17 @@ namespace Cantera {
          */
         void setMassFraction(int j, int k, doublereal y) {
             m_fixedy(k,j) = y;
-            m_do_species[k] = false;
+            m_do_species[k] = true; // false;
         }
 
+        /**
+         * The fixed temperature value at point j.
+         */
         doublereal T_fixed(int j) const {return m_fixedtemp[j];}
+
+        /**
+         * The fixed mass fraction value of species k at point j.
+         */
         doublereal Y_fixed(int k, int j) const {return m_fixedy(k,j);}
 
         virtual string componentName(int n) const;
@@ -143,7 +190,7 @@ namespace Cantera {
         void outputTEC(ostream &s, const doublereal* x, 
             string title, int zone);
 
-        void showSolution(ostream& s, const doublereal* x);
+        virtual void showSolution(ostream& s, const doublereal* x);
 
         void save(string fname, string id, string desc, doublereal* soln);
         virtual void save(XML_Node& o, doublereal* sol);
@@ -159,9 +206,12 @@ namespace Cantera {
             if (j < 0)
                 for (int i = 0; i < m_points; i++)
                     m_do_energy[i] = true;
-            else
-            m_do_energy[j] = true;
-            requestJacUpdate();
+            else 
+                m_do_energy[j] = true;
+            m_refiner->setActive(0, true);
+            m_refiner->setActive(1, true);
+            m_refiner->setActive(2, true);
+            needJacUpdate();
         }
 
         void fixTemperature(int j=-1) {
@@ -170,7 +220,10 @@ namespace Cantera {
                     m_do_energy[i] = false;
                 }
             else m_do_energy[j] = false;
-            requestJacUpdate();
+            m_refiner->setActive(0, false);
+            m_refiner->setActive(1, false);
+            m_refiner->setActive(2, false);
+            needJacUpdate();
         }
 
         bool doSpecies(int k) { return m_do_species[k]; }
@@ -182,7 +235,7 @@ namespace Cantera {
                     m_do_species[i] = true;
             }
             else m_do_species[k] = true;
-            requestJacUpdate();
+            needJacUpdate();
         }
 
         void setEnergyFactor(doublereal efctr);
@@ -193,17 +246,14 @@ namespace Cantera {
                     m_do_species[i] = false;
             }
             else m_do_species[k] = false;
-            requestJacUpdate();
+            needJacUpdate();
             //m_jac->setAge(10000);
         }
 
         void integrateChem(doublereal* x,doublereal dt);
 
-        doublereal z(int j) const {return m_z[j];}
-        doublereal zmin() const { return m_z[0]; }
-        doublereal zmax() const { return m_z[m_points - 1]; }
-
         void resize(int points);
+
         virtual void setFixedPoint(int j0, doublereal t0){}
 
 
@@ -226,14 +276,10 @@ namespace Cantera {
 
     protected:
 
-    // used to write mole fractions to plot files.
+        // used to write mass fractions to plot files.
 
         doublereal component(const doublereal* x, int i, int j) const {
             doublereal xx = x[index(i,j)];
-            //if (i >= 4) {
-            //    return xx*m_wtm[j]/m_wt[i-4];
-            //}
-            //else return xx;
             return xx;
         }
 
@@ -247,11 +293,16 @@ namespace Cantera {
 
         doublereal wdot(int k, int j) const {return m_wdot(k,j);}
 
+        /// write the net production rates at point j into array m_wdot
         void getWdot(doublereal* x,int 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.
+         */
         void updateThermo(const doublereal* x, int j0, int j1) {
             int j;
             for (j = j0; j <= j1; j++) {
@@ -262,6 +313,7 @@ namespace Cantera {
             }
         }
 
+
         //--------------------------------
         // central-differenced derivatives
         //--------------------------------
@@ -329,17 +381,14 @@ namespace Cantera {
         // differencing, assuming u(z) is negative
 
         doublereal dVdz(const doublereal* x,int j) const {
-            //return (V(x,j+1) - V(x,j-1))/(m_dz[j] + m_dz[j-1]);
             return (V(x,j) - V(x,j-1))/m_dz[j-1];
         } 
 
         doublereal dYdz(const doublereal* x,int k, int j) const {
-            //return (Y(x,k,j+1) - Y(x,k,j))/m_dz[j]; 
             return (Y(x,k,j) - Y(x,k,j-1))/m_dz[j-1]; 
         } 
 
         doublereal dTdz(const doublereal* x,int j) const {
-            // return (T(x,j+1) - T(x,j))/m_dz[j];
             return (T(x,j) - T(x,j-1))/m_dz[j-1];
         }
         
@@ -357,6 +406,12 @@ namespace Cantera {
 
         void updateDiffFluxes(const doublereal* x, int j0, int j1);
 
+        //---------------------------------------------------------
+        //
+        //             member data
+        //
+        //---------------------------------------------------------
+
         // inlet
         doublereal m_inlet_u;
         doublereal m_inlet_V;
@@ -369,8 +424,9 @@ namespace Cantera {
 
         doublereal m_press;        // pressure
 
+        // grid parameters
         vector_fp m_dz;
-        vector_fp m_z;
+        //vector_fp m_z;
 
         // mixture thermo properties
         vector_fp m_rho;
@@ -393,14 +449,11 @@ namespace Cantera {
 
         int m_nsp;
 
-        //IdealGasMix*    m_fluid;
-        //thermo_t*          m_phase;
-        igthermo_t*         m_thermo;
-        kinetics_t*       m_kin;
+        igthermo_t*     m_thermo;
+        kinetics_t*     m_kin;
         Transport*      m_trans;
-        //ImplicitChem*   m_integrator;
 
-        MultiJac*          m_jac;
+        MultiJac*       m_jac;
 
         bool m_ok;
 
@@ -410,23 +463,30 @@ namespace Cantera {
         vector<bool> m_do_species;
         int m_transport_option;
 
-        vector_fp m_zest;
-        Array2D   m_yest;
+        // solution estimate
+        //vector_fp m_zest;
+        //Array2D   m_yest;
+
+        // fixed T and Y values
         Array2D   m_fixedy;
         vector_fp m_fixedtemp;
         vector_fp m_zfix;
         vector_fp m_tfix;
 
         vector<FlowBdry::Boundary*>  m_boundary;
+
         doublereal m_efctr;
 
     private:
 
-        void requestJacUpdate();
         vector_fp m_ybar;
     };
 
 
+
+    /**
+     * A class for axisymmetric stagnation flows.
+     */
     class AxiStagnFlow : public StFlow {
     public:
         AxiStagnFlow(igthermo_t* ph = 0, int nsp = 1, int points = 1) :
diff --git a/Cantera/src/oneD/newton_utils.cpp b/Cantera/src/oneD/newton_utils.cpp
index 775e65637..8b20750c8 100644
--- a/Cantera/src/oneD/newton_utils.cpp
+++ b/Cantera/src/oneD/newton_utils.cpp
@@ -38,10 +38,11 @@ namespace Cantera {
             for (j = 0; j < np; j++) {
                 val = x[index(m,j)];
                 if (loglevel > 0) {
-                    if (val > above + Tiny || val < below - Tiny)
-                        cout << "ERROR: solution out of bounds. " 
-                             << r.componentName(m) << "(" << j << ") = " << val
-                             << " (" << below << ", " << above << ")" << endl;
+                    if (val > above + Tiny || val < below - Tiny) {
+                        sprintf(buf, "domain %d: %20s(%d) = %10.3e (%10.3e, %10.3e)\n",
+                            r.domainIndex(), r.componentName(m).c_str(), j, val, below, above);
+                        writelog(string("ERROR: solution out of bounds.\n")+buf);
+                    }
                 }
 
                 newval = val + step[index(m,j)];
@@ -57,13 +58,13 @@ namespace Cantera {
                 if (loglevel > 1 && (newval > above || newval < below)) {
                     if (!wroteTitle) { 
                         writelog("\nNewton step takes solution out of bounds.\n\n");
-                        sprintf(buf,"  %12s  %4s  %10s  %10s  %10s  %10s\n",
-                            "component","pt","value","step","min","max");
+                        sprintf(buf,"  %12s  %12s  %4s  %10s  %10s  %10s  %10s\n",
+                            "domain","component","pt","value","step","min","max");
                         wroteTitle = true;
                         writelog(buf);
                     }
-                    sprintf(buf, "  %12s  %4i  %10.3e  %10.3e  %10.3e  %10.3e\n",
-                        r.componentName(m).c_str(), j, val, 
+                    sprintf(buf, "          %4i  %12s  %4i  %10.3e  %10.3e  %10.3e  %10.3e\n",
+                        r.domainIndex(), r.componentName(m).c_str(), j, val, 
                         step[index(m,j)], below, above);
                     writelog(buf);
                 }
diff --git a/Cantera/src/oneD/refine.cpp b/Cantera/src/oneD/refine.cpp
new file mode 100644
index 000000000..44a7018f0
--- /dev/null
+++ b/Cantera/src/oneD/refine.cpp
@@ -0,0 +1,215 @@
+
+#include <map>
+#include <algorithm>
+#include "Resid1D.h"
+
+#include "refine.h"
+
+using namespace std;
+
+namespace Cantera {
+
+
+    template<class M>
+    static bool has_key(const M& m, int j) {
+        if (m.find(j) != m.end()) return true;
+        return false;
+    }
+
+
+    /**
+     * Return the square root of machine precision.
+     */
+    static doublereal eps() {
+        doublereal e = 1.0;
+        while (1.0 + e != 1.0) e *= 0.5;
+        return sqrt(e);
+    }
+
+    
+    Refiner::Refiner(Resid1D& domain) :
+        m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_min_range(0.01),
+        m_domain(&domain)
+    {
+        m_nv = m_domain->nComponents();
+        m_active.resize(m_nv, true);
+        m_thresh = eps();
+    }
+            
+
+    int Refiner::analyze(int n, const doublereal* z, 
+        const doublereal* x) {
+
+        if (m_domain->nPoints() <= 1) return 0;
+        m_nv = m_domain->nComponents();
+
+        //m_ok = false;
+
+        // check consistency
+        if (n != m_domain->nPoints()) return -1;
+
+        m_loc.clear();
+        m_c.clear();
+
+        /**
+         * find locations where cell size ratio is too large.
+         */
+        int j;
+        vector_fp dz(n-1, 0.0);
+        dz[0] = z[1] - z[0];
+        for (j = 1; j < n-1; j++) {
+            dz[j] = z[j+1] - z[j];
+            if (dz[j] > m_ratio*dz[j-1]) {
+                m_loc[j] = 1;
+                m_c["point "+int2str(j)] = 1;
+            }
+            if (dz[j] < dz[j-1]/m_ratio) {
+                m_loc[j-1] = 1;
+                m_c["point "+int2str(j-1)] = 1;                
+            }
+        }
+
+        string name;
+        doublereal vmin, vmax, smin, smax, aa, ss;
+        doublereal dmax, r;
+        vector_fp v(n), s(n-1);
+        for (int i = 0; i < m_nv; i++) {
+            //cout << i << "   " << m_nv << "  " << m_active[i] << endl;
+            if (m_active[i]) {
+                name = m_domain->componentName(i);
+            
+                // get component i at all points
+                for (j = 0; j < n; j++) v[j] = value(x, i, j); 
+
+                // slope of component i
+                for (j = 0; j < n-1; j++)
+                    s[j] = (value(x, i, j+1) - value(x, i, j))/
+                           (z[j+1] - z[j]);
+
+                // find the range of values and slopes
+            
+                vmin = *min_element(v.begin(), v.end());
+                vmax = *max_element(v.begin(), v.end());
+                smin = *min_element(s.begin(), s.end());
+                smax = *max_element(s.begin(), s.end());
+
+                // max absolute values of v and s
+                aa = fmaxx(abs(vmax), abs(vmin));
+                ss = fmaxx(abs(smax), abs(smin));
+
+
+                // refine based on component i only if the range of v is
+                // greater than a fraction 'min_range' of max |v|. This
+                // eliminates components that consist of small fluctuations
+                // on a constant background.
+
+                if ((vmax - vmin) > m_min_range*aa) {
+
+                    // maximum allowable difference in value between
+                    // adjacent points.
+                
+                    dmax = m_slope*(vmax - vmin) + m_thresh;
+                    for (j = 0; j < n-1; j++) {
+                        r = abs(v[j+1] - v[j])/dmax;
+                        if (r > 1.0) {
+                            m_loc[j] = 1;
+                            m_c[name] = 1;
+                        }
+                    }
+                }
+
+                        
+                // refine based on the slope of component i only if the
+                // range of s is greater than a fraction 'min_range' of max
+                // |s|. This eliminates components that consist of small
+                // fluctuations on a constant slope background.
+            
+                if ((smax - smin) > m_min_range*ss) {
+
+                    // maximum allowable difference in slope between
+                    // adjacent points.
+                    dmax = m_curve*(smax - smin);
+                    for (j = 0; j < n-2; j++) {
+                        r = abs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]);
+                        if (r > 1.0) {
+                            m_c[name] = 1;
+                            m_loc[j] = 1;
+                            m_loc[j+1] = 1;
+                        }
+                        //cout << "at point " << j << " slope r = "
+                        //     << r << " for " << name << endl
+                        //     << "    threshold = " << m_thresh << endl;
+                    }
+                }
+                //cout << name << "  " << m_curve << "  " << smax << "  " << smin << "   " << ss << "   " << m_min_range << endl;
+            }
+        }
+        return m_loc.size();
+    }
+
+    double Refiner::value(const double* x, int i, int j) {
+        return x[m_domain->index(i,j)];
+    }
+
+    void Refiner::show() {
+        int nnew = m_loc.size();
+        if (nnew > 0) {
+            writelog("Refining grid.  "
+                "New points inserted after grid points ");
+            map<int, int>::const_iterator b = m_loc.begin();
+            for (; b != m_loc.end(); ++b) {
+                writelog(int2str(b->first)+" ");
+            }
+            writelog("\n");
+            writelog("to resolve ");
+            map<string, int>::const_iterator bb = m_c.begin();
+            for (; bb != m_c.end(); ++bb) {
+                writelog(string(bb->first)+" ");
+            }
+            writelog("\n");
+        }
+    }
+
+
+    int Refiner::getNewGrid(int n, const doublereal* z, 
+        int nn, doublereal* zn) {
+        int j;
+        int nnew = m_loc.size();
+        if (nnew + n > nn) {
+            throw CanteraError("Refine::getNewGrid",
+                "array size too small.");
+            return -1;
+        }
+
+        int jn = 0;
+        if (m_loc.size() == 0) {
+            copy(z, z + n,  zn);
+            return 0;
+        }
+
+        for (j = 0; j < n - 1; j++) {
+            zn[jn] = z[j];
+            jn++;
+            if (has_key(m_loc, j)) {
+                zn[jn] = 0.5*(z[j] + z[j+1]);
+                jn++;
+            }
+        }
+        zn[jn] = z[n-1];
+        return 0;
+    }
+
+//         int npts = znew.size();
+//         newsoln.resize(npts*ncomp);
+//         newsoln = Numeric.zeros((npts, ncomp),'d')
+//         for i in range(ncomp):
+//             for j in range(npts):
+//                 newsoln[j,i] = interp.interp(znew[j],grid,solution[:,i])
+
+//         return (Numeric.array(znew), Numeric.array(znew), newsoln, self.ok)
+                               
+    
+                
+            
+    
+}
diff --git a/Cantera/src/oneD/refine.h b/Cantera/src/oneD/refine.h
new file mode 100644
index 000000000..363a2b5ec
--- /dev/null
+++ b/Cantera/src/oneD/refine.h
@@ -0,0 +1,47 @@
+#ifndef CT_REFINE_H
+#define CT_REFINE_H
+
+namespace Cantera {
+
+    class Resid1D;
+
+    class Refiner {
+
+    public:
+
+        Refiner(Resid1D& domain);
+        virtual ~Refiner(){}
+
+        void setCriteria(doublereal ratio = 10.0,
+            doublereal slope = 0.8, 
+            doublereal curve = 0.8) {
+            m_ratio = ratio; m_slope = slope; m_curve = curve;
+        }
+        void setActive(int comp, bool state = true) { m_active[comp] = state; }
+
+        int analyze(int n, const doublereal* z, const doublereal* x);
+        int getNewGrid(int n, const doublereal* z, int nn, doublereal* znew);
+        //int getNewSoln(int n, const doublereal* x, doublereal* xnew);
+        int nNewPoints() { return m_loc.size(); }
+        void show();
+        bool newPointNeeded(int j) { 
+            return m_loc.find(j) != m_loc.end();
+        }
+        double value(const double* x, int i, int j);
+
+    protected:
+
+        map<int, int>       m_loc;
+        map<string, int>    m_c;
+        vector<bool>        m_active;
+        doublereal          m_ratio, m_slope, m_curve;
+        doublereal          m_min_range;
+        Resid1D*            m_domain;
+        int                 m_nv;
+        doublereal          m_thresh;
+
+    };
+
+}
+
+#endif