diff --git a/Cantera/src/oneD/Domain1D.cpp b/Cantera/src/oneD/Domain1D.cpp
new file mode 100644
index 000000000..1e628d6a8
--- /dev/null
+++ b/Cantera/src/oneD/Domain1D.cpp
@@ -0,0 +1,190 @@
+/**
+ * @file Domain1D.cpp
+ *
+ */
+
+#include "Domain1D.h"
+
+namespace Cantera {
+
+    void Domain1D::
+    setTolerances(int nr, const doublereal* rtol, 
+        int na, const doublereal* atol, int ts) {
+        if (nr < m_nv || na < m_nv)
+            throw CanteraError("Domain1D::setTolerances",
+                "wrong array size for solution error tolerances. "
+                "Size should be at least "+int2str(m_nv));
+        if (ts >= 0) {
+            copy(rtol, rtol + m_nv, m_rtol_ss.begin());
+            copy(atol, atol + m_nv, m_atol_ss.begin());
+        }
+        if (ts <= 0) {
+            copy(rtol, rtol + m_nv, m_rtol_ts.begin());
+            copy(atol, atol + m_nv, m_atol_ts.begin());
+        }
+    }
+
+    void Domain1D::
+    setTolerances(int n, doublereal rtol, doublereal atol, int ts) {
+        if (ts >= 0) {
+            m_rtol_ss[n] = rtol;
+            m_atol_ss[n] = atol;
+        }
+        if (ts <= 0) {
+            m_rtol_ts[n] = rtol;
+            m_atol_ts[n] = atol;
+        }
+    }
+
+    void Domain1D::
+    setTolerances(doublereal rtol, doublereal atol,int ts) {
+        for (int n = 0; n < m_nv; n++){
+            if(ts >= 0) {
+                m_rtol_ss[n] = rtol;
+                m_atol_ss[n] = atol;
+            }
+            if (ts <= 0) {
+                m_rtol_ts[n] = rtol;
+                m_atol_ts[n] = atol;
+            }
+        }
+    }
+
+    void Domain1D::
+    setTolerancesTS(doublereal rtol, doublereal atol) {
+        for (int n = 0; n < m_nv; n++){
+            m_rtol_ts[n] = rtol;
+            m_atol_ts[n] = atol;
+        }
+    }
+
+    void Domain1D::
+    setTolerancesSS(doublereal rtol, doublereal atol) {
+        for (int n = 0; n < m_nv; n++){
+            m_rtol_ss[n] = rtol;
+            m_atol_ss[n] = atol;
+        }
+    }
+
+    void Domain1D::
+    eval(int jg, doublereal* xg, doublereal* rg, 
+        integer* mask, doublereal rdt) {
+
+        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
+        doublereal* x = xg + loc();
+        doublereal* rsd = rg + loc();
+        integer* diag = mask + loc();
+        
+        int jmin, jmax, jpt,  j, i;
+        jpt = jg - firstPoint();
+        
+        if (jg < 0) {      // evaluate all points
+            jmin = 0;
+            jmax = m_points - 1;
+        }
+        else {            // evaluate points for Jacobian
+            jmin = max(jpt-1, 0);
+            jmax = min(jpt+1,m_points-1);
+        }
+
+        for (j = jmin; j <= jmax; j++) {
+            if (j == 0 || j == m_points - 1) {
+                for (i = 0; i < m_nv; i++) {
+                    rsd[index(i,j)] = residual(x,i,j);
+                    diag[index(i,j)] = 0;
+                }
+            }
+            else {
+                for (i = 0; i < m_nv; i++) {
+                    rsd[index(i,j)] = residual(x,i,j) 
+                                      - timeDerivativeFlag(i)*rdt*(value(x,i,j) - prevSoln(i,j));
+                    diag[index(i,j)] = timeDerivativeFlag(i);
+                }
+            }
+        }
+    }
+
+
+    // called to set up initial grid, and after grid refinement
+    void Domain1D::setupGrid(int n, const doublereal* z) {
+        resize(m_nv, n);
+        int j;
+        for (j = 0; j < m_points; j++) m_z[j] = z[j];
+    }
+
+
+    void drawline() {
+        writelog("\n-------------------------------------"
+            "------------------------------------------");
+    }
+
+    /**
+     * Print the solution.
+     */    
+    void Domain1D::showSolution(const doublereal* x) {
+        int nn = m_nv/5;
+        int i, j, n;
+        //char* buf = new char[100];
+        char buf[100];
+        doublereal v;
+        for (i = 0; i < nn; i++) {
+            drawline();
+            sprintf(buf, "\n        z   ");
+            writelog(buf);
+            for (n = 0; n < 5; n++) { 
+                sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
+                writelog(buf);
+            }
+            drawline();
+            for (j = 0; j < m_points; j++) {
+                sprintf(buf, "\n %10.4g ",m_z[j]);
+                writelog(buf);
+                for (n = 0; n < 5; n++) { 
+                    v = value(x, i*5+n, j);
+                    sprintf(buf, " %10.4g ",v);
+                    writelog(buf);
+                }
+            }
+            writelog("\n");
+        }
+        int nrem = m_nv - 5*nn;
+        drawline();
+        sprintf(buf, "\n        z   ");
+        writelog(buf);
+        for (n = 0; n < nrem; n++) {
+            sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
+            writelog(buf);
+        }
+        drawline();
+        for (j = 0; j < m_points; j++) {
+            sprintf(buf, "\n %10.4g ",m_z[j]);
+            writelog(buf);
+            for (n = 0; n < nrem; n++) { 
+                v = value(x, nn*5+n, j);
+                sprintf(buf, " %10.4g ", v);
+                writelog(buf);
+            }
+        }
+        writelog("\n");
+    }
+
+    // initial solution
+    void Domain1D::_getInitialSoln(doublereal* x) {
+        for (int j = 0; j < m_points; j++) {
+            for (int n = 0; n < m_nv; n++) {
+                x[index(n,j)] = initialValue(n,j);
+            }
+        }
+    }
+
+    doublereal Domain1D::initialValue(int n, int j) { 
+        throw CanteraError("Domain1D::initialValue",
+            "base class method called!");
+    }
+
+} // namespace
diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h
index c7388f68f..44f7e225c 100644
--- a/Cantera/src/oneD/Domain1D.h
+++ b/Cantera/src/oneD/Domain1D.h
@@ -102,12 +102,16 @@ namespace Cantera {
          * actions required to resize the domain.
          */
         virtual void resize(int nv, int np) {
+            // if the number of components is being changed, then a
+            // new grid refiner is required.
             if (nv != m_nv || !m_refiner) {
                 m_nv = nv;
                 delete m_refiner;
                 m_refiner = new Refiner(*this);
             }
             m_nv = nv;
+            m_td.resize(m_nv, 1);
+            m_name.resize(m_nv,"");
             m_max.resize(m_nv, 0.0);
             m_min.resize(m_nv, 0.0);
             m_rtol_ss.resize(m_nv, 1.0e-8);
@@ -120,6 +124,7 @@ namespace Cantera {
             locate();
         }
 
+        /// Return a reference to the grid refiner.
         Refiner& refiner() { return *m_refiner; }
 
         /// Number of components at each grid point.
@@ -130,9 +135,19 @@ namespace Cantera {
 
         /// Name of the nth component. May be overloaded.
         virtual string componentName(int n) const { 
-            return "component " + int2str(n); 
+            if (m_name[n] != "") return m_name[n];
+            else return "component " + int2str(n); 
         }
 
+        void setComponentName(int n, string name) {
+            m_name[n] = name;
+        }
+
+        void setComponentType(int n, int ctype) {
+            if (ctype == 0) setAlgebraic(n);
+        }
+
+        /// index of component with name \a name.
         int componentIndex(string name) {
             int nc = nComponents();
             for (int n = 0; n < nc; n++) {
@@ -160,60 +175,23 @@ namespace Cantera {
             m_max[n] = upper;
         }
 
+        /// set the error tolerances for all solution components. 
         void setTolerances(int nr, const doublereal* rtol, 
-            int na, const doublereal* atol, int ts = 0) {
-            if (nr < m_nv || na < m_nv)
-                throw CanteraError("Domain1D::setTolerances",
-                    "wrong array size for solution error tolerances. "
-                    "Size should be at least "+int2str(m_nv));
-            if (ts >= 0) {
-                copy(rtol, rtol + m_nv, m_rtol_ss.begin());
-                copy(atol, atol + m_nv, m_atol_ss.begin());
-            }
-            if (ts <= 0) {
-                copy(rtol, rtol + m_nv, m_rtol_ts.begin());
-                copy(atol, atol + m_nv, m_atol_ts.begin());
-            }
-        }
+            int na, const doublereal* atol, int ts = 0);
 
-        void setTolerances(int n, doublereal rtol, doublereal atol, int ts = 0) {
-            if (ts >= 0) {
-                m_rtol_ss[n] = rtol;
-                m_atol_ss[n] = atol;
-            }
-            if (ts <= 0) {
-                m_rtol_ts[n] = rtol;
-                m_atol_ts[n] = atol;
-            }
-        }
+        /// set the error tolerances for solution component \a n.
+        void setTolerances(int n, doublereal rtol, doublereal atol, int ts = 0);
 
         //added by Karl Meredith
-        void setTolerances(doublereal rtol, doublereal atol,int ts=0) {
-            for (int n=0;n<m_nv;n++){
-                if(ts>=0) {
-                    m_rtol_ss[n] = rtol;
-                    m_atol_ss[n] = atol;
-                }
-                if (ts <= 0) {
-                    m_rtol_ts[n] = rtol;
-                    m_atol_ts[n] = atol;
-                }
-            }
-        }
+        /// set scalar error tolerances. All solution components will
+        /// have the same relative and absolute error tolerances.
+        void setTolerances(doublereal rtol, doublereal atol,int ts=0);
+
         //added by Karl Meredith
-        void setTolerancesTS(doublereal rtol, doublereal atol) {
-            for (int n=0;n<m_nv;n++){
-                m_rtol_ts[n] = rtol;
-                m_atol_ts[n] = atol;
-            }
-        }
+        void setTolerancesTS(doublereal rtol, doublereal atol);
+
         //added by Karl Meredith
-        void setTolerancesSS(doublereal rtol, doublereal atol) {
-            for (int n=0;n<m_nv;n++){
-                m_rtol_ss[n] = rtol;
-                m_atol_ss[n] = atol;
-            }
-        }
+        void setTolerancesSS(doublereal rtol, doublereal atol);
 				
         /// Relative tolerance of the nth component.
         doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); }
@@ -271,11 +249,15 @@ namespace Cantera {
          * evaluate the residual function at all points.
          */         
         virtual void eval(int j, doublereal* x, doublereal* r, 
-            integer* mask, doublereal rdt=0.0) {
-            throw CanteraError("Domain1D::eval",
-                "residual function not defined.");
+            integer* mask, doublereal rdt=0.0);
+
+        virtual doublereal residual(doublereal* x, int n, int j) {
+            throw CanteraError("Domain1D::residual","residual function must be overloaded in derived class");
         }
 
+        int timeDerivativeFlag(int n) { return m_td[n];}
+        void setAlgebraic(int n) { m_td[n] = 0; }
+
         /**
          * Does nothing.
          */
@@ -284,7 +266,7 @@ namespace Cantera {
         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 {
+        doublereal value(const doublereal* x, int n, int j) const {
             return x[index(n,j)];
         }
 
@@ -396,7 +378,7 @@ namespace Cantera {
         virtual void getTransientMask(integer* mask){}
 
         virtual void showSolution(ostream& s, const doublereal* x) {}
-        virtual void showSolution(const doublereal* x) {}
+        virtual void showSolution(const doublereal* x);
 
         virtual void restore(const XML_Node& dom, doublereal* soln) {}
 
@@ -425,7 +407,8 @@ namespace Cantera {
         const vector_fp& grid() const { return m_z; }
         doublereal grid(int point) { return m_z[point]; }
 
-        virtual void setupGrid(int n, const doublereal* z) {}
+        virtual void setupGrid(int n, const doublereal* z);
+
         void setGrid(int n, const doublereal* z);
 
         /**
@@ -436,10 +419,12 @@ namespace Cantera {
          * prior to installing this domain into the container to be
          * written to the global solution vector.
          */
-        virtual void _getInitialSoln(doublereal* x) {
-            throw CanteraError("Domain1D::_getInitialSoln",
-                "base class method _getInitialSoln called!");
-        }
+        virtual void _getInitialSoln(doublereal* x);
+
+        /**
+         * Initial value of solution component \a n at grid point \a j.
+         */
+        virtual doublereal initialValue(int n, int j);
 
         /**
          * In some cases, a domain may need to set parameters that
@@ -479,8 +464,8 @@ namespace Cantera {
         Domain1D *m_left, *m_right;
         string m_id, m_desc;
         Refiner* m_refiner;
-
-
+        vector_int m_td;
+        vector<string> m_name;
 
     private:
 
diff --git a/Cantera/src/oneD/Inlet1D.h b/Cantera/src/oneD/Inlet1D.h
index fbe80ac28..4071b34d5 100644
--- a/Cantera/src/oneD/Inlet1D.h
+++ b/Cantera/src/oneD/Inlet1D.h
@@ -202,6 +202,7 @@ namespace Cantera {
         virtual ~Empty1D(){}
 
         virtual string componentName(int n) const;
+        virtual void showSolution(const doublereal* x) {}
 
         virtual void init();
 
diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp
index b91662e28..d7dbdcee9 100644
--- a/Cantera/src/oneD/Sim1D.cpp
+++ b/Cantera/src/oneD/Sim1D.cpp
@@ -186,15 +186,18 @@ namespace Cantera {
 
     void Sim1D::showSolution(ostream& s) {
         for (int n = 0; n < m_nd; n++) {
-            domain(n).showSolution(s, m_x.begin() + start(n));
+            if (domain(n).domainType() != cEmptyType)
+                domain(n).showSolution(s, m_x.begin() + start(n));
         }
     }
 
     void Sim1D::showSolution() {
         for (int n = 0; n < m_nd; n++) {
-            writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
-                       +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
-            domain(n).showSolution(m_x.begin() + start(n));
+            if (domain(n).domainType() != cEmptyType) {
+                writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
+                    +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
+                domain(n).showSolution(m_x.begin() + start(n));
+            }
         }
     }
 
diff --git a/Cantera/src/oneD/oneD_files.cpp b/Cantera/src/oneD/oneD_files.cpp
index 1013494a0..6a4ce3623 100644
--- a/Cantera/src/oneD/oneD_files.cpp
+++ b/Cantera/src/oneD/oneD_files.cpp
@@ -6,3 +6,4 @@
 #include "boundaries1D.cpp"
 #include "refine.cpp" 
 #include "Sim1D.cpp" 
+#include "Domain1D.cpp"