diff --git a/Cantera/src/zeroD/Makefile.in b/Cantera/src/zeroD/Makefile.in
index ba0114d0d..400b85cf4 100644
--- a/Cantera/src/zeroD/Makefile.in
+++ b/Cantera/src/zeroD/Makefile.in
@@ -18,9 +18,9 @@ do_ranlib = @DO_RANLIB@
 CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
 
 # stirred reactors
-OBJS    = Reactor.o ReactorBase.o FlowDevice.o Wall.o ReactorNet.o
+OBJS    = Reactor.o ReactorBase.o FlowDevice.o Wall.o ReactorNet.o FlowReactor.o
 ZEROD_H = Reactor.h ReactorBase.h FlowDevice.h Wall.h ReactorNet.h \
-          flowControllers.h PID_Controller.h Reservoir.h
+          flowControllers.h PID_Controller.h Reservoir.h FlowReactor.h
 
 CXX_INCLUDES = -I..
 ZEROD_LIB = @buildlib@/libzeroD.a
diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp
index bb5c9f139..bd13941af 100644
--- a/Cantera/src/zeroD/Reactor.cpp
+++ b/Cantera/src/zeroD/Reactor.cpp
@@ -13,7 +13,7 @@
 #endif
 
 #include "Reactor.h"
-#include "../CVode.h"
+//#include "../CVode.h"
 #include "FlowDevice.h"
 #include "Wall.h"
 #include "../InterfaceKinetics.h"
@@ -26,7 +26,9 @@ namespace CanteraZeroD {
     doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y);
 
     Reactor::Reactor() : ReactorBase(), 
+#ifdef INCL_REACTOR_INTEG
                          FuncEval(),
+#endif
                          m_kin(0),
                          m_temp_atol(1.e-11), 
                          m_maxstep(0.0),
@@ -34,7 +36,7 @@ namespace CanteraZeroD {
                          m_Q(0.0), 
                          m_rtol(1.e-9),
                          m_chem(true),
-                         m_energy(true)
+                         m_energy(true), m_nsens(-1)
     {
 #ifdef INCL_REACTOR_INTEG
         m_integ = new CVodeInt;
@@ -47,7 +49,6 @@ namespace CanteraZeroD {
 #endif
     }
 
-
     // overloaded method of FuncEval. Called by the integrator to
     // get the initial conditions.
     void Reactor::getInitialConditions(double t0, size_t leny, double* y) 
@@ -89,7 +90,6 @@ namespace CanteraZeroD {
         }
     }
 
-
     /*
      *  Must be called before calling method 'advance'
      */
@@ -111,8 +111,8 @@ namespace CanteraZeroD {
         m_pressure = m_thermo->pressure();
         m_intEnergy = m_thermo->intEnergy_mass();
 
-        int nt = 0, maxnt = 0;
-        for (int m = 0; m < m_nwalls; m++) {
+        int m, nt = 0, maxnt = 0;
+        for (m = 0; m < m_nwalls; m++) {
             if (m_wall[m]->kinetics(m_lr[m])) {
                 nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
                 if (nt > maxnt) maxnt = nt;
@@ -127,10 +127,23 @@ namespace CanteraZeroD {
             }   
         }
         m_work.resize(maxnt);
-
         m_init = true;
     }
     
+    int Reactor::nSensParams() {
+        if (m_nsens < 0) {
+            // determine the number of sensitivity parameters
+            int m, ns;
+            m_nsens = m_pnum.size();
+            for (m = 0; m < m_nwalls; m++) {
+                ns = m_wall[m]->nSensParams(m_lr[m]);
+                m_nsens_wall.push_back(ns);
+                m_nsens += ns;
+            }
+        }
+        return m_nsens;
+    }
+
     void Reactor::updateState(doublereal* y) {
 
         phase_t& mix = *m_mix;  // define for readability
@@ -140,6 +153,7 @@ namespace CanteraZeroD {
 
         // Set the mass fractions and  density of the mixture.
 
+
         doublereal u   = y[0];
         m_vol          = y[1];
         doublereal* mss = y + 2;
@@ -176,22 +190,46 @@ namespace CanteraZeroD {
         m_mix->saveState(m_state);
     }
 
-
+#ifdef INCL_REACTOR_INTEG
     void Reactor::eval(doublereal time, doublereal* y, doublereal* ydot) 
     {
         updateState(y);          // synchronize the reactor state with y
         evalEqs(time, y, ydot);
     }
+#endif
 
     /*
      * Called by the integrator to evaluate ydot given y at time 'time'.
      */
-    void Reactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot) 
+    void Reactor::evalEqs(doublereal time, doublereal* y, 
+        doublereal* ydot, doublereal* params) 
     {
         int i, k, nk;
         m_time = time;
         m_mix->restoreState(m_state);
 
+        Kinetics* kin;
+        int m, n, npar, ploc;
+        double mult;
+        // process sensitivity parameters
+        if (params) {
+            
+            npar = m_pnum.size();
+            for (n = 0; n < npar; n++) {
+                //m_mult_save[n] = m_kin->multiplier(m_pnum[n]);
+                mult = m_kin->multiplier(m_pnum[n]);
+                m_kin->setMultiplier(m_pnum[n], mult*params[n]);
+                //                m_kin->setMultiplier(m_pnum[n], m_mult_save[n]*params[n]);
+            }
+            ploc = npar;
+            for (m = 0; m < m_nwalls; m++) {
+                if (m_nsens_wall[m] > 0) {
+                    m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
+                    ploc += m_nsens_wall[m];
+                }
+            }
+        }
+
         //        updateState(y);          // synchronize the reactor state with y
 
         m_vdot = 0.0;
@@ -199,7 +237,7 @@ namespace CanteraZeroD {
 
         // compute wall terms
         doublereal vdot, rs0, sum, wallarea;
-        Kinetics* kin;
+        //        Kinetics* kin;
         SurfPhase* surf;
         int lr, ns, loc = m_nsp+2, surfloc;
         fill(m_sdot.begin(), m_sdot.end(), 0.0);
@@ -242,8 +280,6 @@ namespace CanteraZeroD {
          *             - \dot m_{out} Y_{k} + A \dot s_k.
          */
         const doublereal* mw = m_mix->molecularWeights().begin();
-
-        int n;
         if (m_chem) {
             m_kin->getNetProductionRates(ydot+2);   // "omega dot"
         }
@@ -268,7 +304,7 @@ namespace CanteraZeroD {
             ydot[0] = - m_thermo->pressure() * m_vdot - m_Q;
         }
         else {
-            ydot[0] = 0.0;    
+            ydot[0] = 0.0;
         }
 
         // add terms for open system
@@ -287,10 +323,7 @@ namespace CanteraZeroD {
                     ydot[2+n] -= mdot_out * mf[n];
                 }
                 if (m_energy) {
-                    //                    cout << "before = " << ydot[0] << endl;
                     ydot[0] -= mdot_out * enthalpy;
-                    //cout << mdot_out << " " << enthalpy << endl;
-                    //cout << "after = " << ydot[0] << endl;
                 }
             }
 
@@ -308,5 +341,59 @@ namespace CanteraZeroD {
                 }
             }
         }
+
+        // reset sensitivity parameters
+        if (params) {
+            npar = m_pnum.size();
+            for (n = 0; n < npar; n++) {
+                mult = m_kin->multiplier(m_pnum[n]);
+                m_kin->setMultiplier(m_pnum[n], mult/params[n]);
+                //m_kin->setMultiplier(m_pnum[n], m_mult_save[n]);
+            }
+            ploc = npar;
+            for (m = 0; m < m_nwalls; m++) {
+                if (m_nsens_wall[m] > 0) {
+                    m_wall[m]->resetSensitivityParameters(m_lr[m]);
+                    ploc += m_nsens_wall[m];
+                }
+            }
+        }
     }
+
+    void Reactor::addSensitivityReaction(int rxn) {
+        m_pnum.push_back(rxn);
+        m_pname.push_back(name()+": "+m_kin->reactionString(rxn));
+        m_mult_save.push_back(1.0);
+        if (rxn < 0 || rxn >= m_kin->nReactions()) 
+            throw CanteraError("Reactor::addSensitivityReaction",
+                "Reaction number out of range ("+int2str(rxn)+")");
+    }
+
+
+    int Reactor::componentIndex(string nm) const {
+        if (nm == "U") return 0;
+        if (nm == "V") return 1;
+        // check for a gas species name
+        int k = m_mix->speciesIndex(nm);
+        if (k >= 0) return k + 2;
+
+        // check for a wall species
+        int walloffset = 0, kp = 0;
+        thermo_t* th;
+        for (int m = 0; m < m_nwalls; m++) {
+            if (m_wall[m]->kinetics(m_lr[m])) {
+                kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
+                th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
+                k = th->speciesIndex(nm);
+                if (k >= 0) {
+                    return k + 2 + m_nsp + walloffset;
+                }
+                else {
+                    walloffset += th->nSpecies();
+                }
+            }   
+        }
+        return -1;
+    }
+
 }
diff --git a/Cantera/src/zeroD/Reactor.h b/Cantera/src/zeroD/Reactor.h
index 470e95516..542a0ceb7 100644
--- a/Cantera/src/zeroD/Reactor.h
+++ b/Cantera/src/zeroD/Reactor.h
@@ -7,7 +7,7 @@
  */
 
 // Copyright 2001  California Institute of Technology
-
+ 
 #ifndef CT_REACTOR_H
 #define CT_REACTOR_H
 
@@ -18,10 +18,10 @@
 
 #include "ReactorBase.h"
 #include "../FuncEval.h"
-#include "../CVode.h"
+#include "../Integrator.h"
 #include "../Kinetics.h"
 
-#define INCL_REACTOR_INTEG
+#undef INCL_REACTOR_INTEG
 
 namespace CanteraZeroD {
 
@@ -57,9 +57,11 @@ namespace CanteraZeroD {
      * flow rate. Class FuncEval is the class used to define a system
      * of ODE's to be integrated.
      */
-
+#ifdef INCL_REACTOR_INTEG
     class Reactor : public ReactorBase, public FuncEval {
-
+#else
+    class Reactor : public ReactorBase {
+#endif
     public:
 
         /**
@@ -67,7 +69,6 @@ namespace CanteraZeroD {
          */
         Reactor();
 
-
         /**
          * Destructor. Deletes the integrator.
          */
@@ -137,12 +138,14 @@ namespace CanteraZeroD {
             if (m_kin->nReactions() == 0) disableChemistry();
         }
 
+#ifdef INCL_REACTOR_INTEG
         /**
          * Set the maximum step size for integration.
          */
         void setMaxStep(doublereal maxstep) {
             m_maxstep = maxstep;
         }
+#endif
 
         void disableChemistry() { m_chem = false; }
         void enableChemistry() { m_chem = true; }
@@ -153,29 +156,33 @@ namespace CanteraZeroD {
             else m_energy = false;
         } 
 
+
         //-----------------------------------------------------
 
         /** @name References to internal objects */
         //@{
-
+#ifdef INCL_REACTOR_INTEG
         /// Return a reference to the integrator.
         Integrator& integrator() { return *m_integ; }
 
         //@}
-
+#endif
 
         //-----------------------------------------------------
 
         // overloaded methods of class FuncEval
         virtual int neq() { return m_nv; }
+#ifdef INCL_REACTOR_INTEG
 	virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
+#endif
         virtual void getInitialConditions(doublereal t0, size_t leny, 
             doublereal* y);
 
+
         //-----------------------------------------------------
 
         virtual void initialize(doublereal t0 = 0.0);
-	void evalEqs(doublereal t, doublereal* y, doublereal* ydot);
+	virtual void evalEqs(doublereal t, doublereal* y, doublereal* ydot, doublereal* params);
 
         /**
          * Set the mixture to a state consistent with solution
@@ -184,11 +191,24 @@ namespace CanteraZeroD {
 
         virtual void updateState(doublereal* y);
 
+        //        virtual void addSensitivityParam(int stype, int i);
+
+
+
+        virtual int nSensParams();
+        virtual void addSensitivityReaction(int rxn);
+
+        virtual string sensParamID(int p) { return m_pname[p]; }
+
+        //        virtual string component(int k) const;
+        virtual int componentIndex(string nm) const;
+
     protected:
         
         Kinetics*   m_kin;
-
+#ifdef INCL_REACTOR_INTEG
         Integrator* m_integ;         // pointer to integrator
+#endif
         doublereal m_temp_atol;      // tolerance on T
         doublereal m_maxstep;        // max step size
         doublereal m_vdot, m_Q;
@@ -200,6 +220,12 @@ namespace CanteraZeroD {
         bool m_energy;
         int m_nv;
 
+        int m_nsens;
+        vector_int m_pnum;
+        vector<string> m_pname;
+        vector_int m_nsens_wall;
+        vector_fp m_mult_save;
+
     private:
     };
 }
diff --git a/Cantera/src/zeroD/ReactorBase.h b/Cantera/src/zeroD/ReactorBase.h
index d72f0f003..2a1d654ef 100644
--- a/Cantera/src/zeroD/ReactorBase.h
+++ b/Cantera/src/zeroD/ReactorBase.h
@@ -32,6 +32,7 @@ namespace CanteraZeroD {
 
     const int ReactorType = 1;
     const int ReservoirType = 2;
+    const int FlowReactorType = 3;
 
     /**
      * Base class for stirred reactors. 
@@ -145,6 +146,8 @@ namespace CanteraZeroD {
             return 1;
         }
 
+        //    virtual void addSensitivityParam(int stype, int i) {}
+
     protected:
  
         int m_nsp;
diff --git a/Cantera/src/zeroD/ReactorNet.cpp b/Cantera/src/zeroD/ReactorNet.cpp
index 7676b38d2..b97c020ee 100644
--- a/Cantera/src/zeroD/ReactorNet.cpp
+++ b/Cantera/src/zeroD/ReactorNet.cpp
@@ -1,15 +1,16 @@
 #include "ReactorNet.h"
-#include "../CVode.h"
+#include "Integrator.h"
 
 namespace CanteraZeroD {
 
     ReactorNet::ReactorNet() : FuncEval(), m_nr(0), m_nreactors(0),
                                m_integ(0), m_time(0.0), m_init(false), 
-                               m_nv(0), m_rtol(1.0e-9), m_atols(1.0e-15),
+                               m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4), 
+                               m_atols(1.0e-15), m_atolsens(1.0e-4),
                                m_maxstep(-1.0),
-                               m_verbose(false)
+                               m_verbose(false), m_ntotpar(0)
     {
-        m_integ = new CVodeInt;
+        m_integ = newIntegrator("CVODE");// CVodeInt;
 
         // use backward differencing, with a full Jacobian computed
         // numerically, and use a Newton linear iterator
@@ -29,23 +30,34 @@ namespace CanteraZeroD {
             writelog("Initializing reactor network.\n");
         }
         for (n = 0; n < m_nr; n++) {
-            if (m_r[n]->type() == ReactorType) {
+            if (m_r[n]->type() >= ReactorType) {
                 m_r[n]->initialize(t0);
                 Reactor* r = (Reactor*)m_r[n];
                 m_reactors.push_back(r);
                 nv = r->neq();
                 m_size.push_back(nv);
+                m_nparams.push_back(r->nSensParams());
+                m_ntotpar += r->nSensParams();
                 m_nv += nv;
                 m_nreactors++;
                 if (m_verbose) {
                     sprintf(buf,"Reactor %d: %d variables.\n",n,nv);
                     writelog(buf);
+                    sprintf(buf,"            %d sensitivity params.\n",
+                        r->nSensParams());
+                    writelog(buf);
+                }
+                if (m_r[n]->type() == FlowReactorType && m_nr > 1) {
+                    throw CanteraError("ReactorNet::initialize",
+                        "FlowReactors must be used alone.");
                 }
             }
         }
+
         m_atol.resize(neq());
         fill(m_atol.begin(), m_atol.end(), m_atols);
         m_integ->setTolerances(m_rtol, neq(), m_atol.begin());
+        m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
         m_integ->setMaxStepSize(m_maxstep);
         if (m_verbose) {
             sprintf(buf, "Number of equations: %d\n", neq());
@@ -79,20 +91,39 @@ namespace CanteraZeroD {
         return m_time;
     }
 
-    void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot) {
+//     void ReactorNet::addSensitivityParam(int n, int stype, int i) {
+//         m_reactors[n]->addSensitivityParam(int stype, int i);
+//         m_sensreactor.push_back(n);
+//         m_nSenseParams++;
+//     }
+
+//     void ReactorNet::setParameters(int np, double* p) {
+//         int n, nr;
+//         for (n = 0; n < np; n++) {
+//             if (n < m_nSenseParams) {
+//                 nr = m_sensreactor[n];
+//                 m_reactors[nr]->setParameter(n, p[n]);
+//             }
+//         }
+//     }
+        
+    void ReactorNet::eval(doublereal t, doublereal* y, 
+        doublereal* ydot, doublereal* p) {
         int n;
         int start = 0;
-
+        int pstart = 0;
         // use a try... catch block, since exceptions are not passed
         // through CVODE, since it is C code
         try {
             updateState(y);
             for (n = 0; n < m_nreactors; n++) {
-                m_reactors[n]->evalEqs(t, y + start, ydot + start);
+                m_reactors[n]->evalEqs(t, y + start, 
+                    ydot + start, p + pstart);
                 start += m_size[n];
+                pstart += m_nparams[n];
             }
         }
-        catch (CanteraError) {
+        catch (...) {
             showErrors();
             error("Terminating execution.");
         }
@@ -116,5 +147,13 @@ namespace CanteraZeroD {
             start += m_size[n];
         }
     }
+
+    int ReactorNet::globalComponentIndex(string species, int reactor) {
+        int start = 0;
+        int n;
+        for (n = 0; n < reactor; n++) start += m_size[n];
+        return start + m_reactors[n]->componentIndex(species);
+    }
+
 }
 
diff --git a/Cantera/src/zeroD/ReactorNet.h b/Cantera/src/zeroD/ReactorNet.h
index ae663f486..2f652cc7b 100644
--- a/Cantera/src/zeroD/ReactorNet.h
+++ b/Cantera/src/zeroD/ReactorNet.h
@@ -20,7 +20,7 @@
 
 #include "Reactor.h"
 #include "../FuncEval.h"
-#include "../CVode.h"
+#include "../Integrator.h"
 
 namespace CanteraZeroD {
 
@@ -59,6 +59,12 @@ namespace CanteraZeroD {
             m_init = false;
         }
 
+        void setSensitivityTolerances(doublereal rtol, doublereal atol) {
+            if (rtol >= 0.0) m_rtolsens = rtol;
+            if (atol >= 0.0) m_atolsens = atol;
+            m_init = false;
+        }
+
         /// Current value of the simulation time.
         doublereal time() { return m_time; }
 
@@ -92,20 +98,33 @@ namespace CanteraZeroD {
             return *m_r[n];
         }
 
+        void setVerbose(bool v = true) { m_verbose = v; }
 
         /// Return a reference to the integrator.
         Integrator& integrator() { return *m_integ; }
 
         void updateState(doublereal* y);
 
+        double sensitivity(int k, int p) {
+            return m_integ->sensitivity(k, p)/m_integ->solution(k);
+        }
+
+        double sensitivity(string species, int p, int reactor=0) {
+            int k = globalComponentIndex(species, reactor);
+            return sensitivity(k, p);
+        }
+
         //-----------------------------------------------------
 
         // overloaded methods of class FuncEval
         virtual int neq() { return m_nv; }
-	virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
+	virtual void eval(doublereal t, doublereal* y, 
+            doublereal* ydot, doublereal* p);
         virtual void getInitialConditions(doublereal t0, size_t leny, 
             doublereal* y);
+        virtual int nparams() { return m_ntotpar; }
 
+        int globalComponentIndex(string species, int reactor=0);
 
     protected:
 
@@ -119,10 +138,12 @@ namespace CanteraZeroD {
         int m_nv;
         vector_int m_size;
         vector_fp m_atol;
-        doublereal m_rtol;
-        doublereal m_atols;
+        doublereal m_rtol, m_rtolsens;
+        doublereal m_atols, m_atolsens;
         doublereal m_maxstep;
         bool m_verbose;
+        int m_ntotpar;
+        vector_int m_nparams;
 
     private:
 
diff --git a/Cantera/src/zeroD/Wall.cpp b/Cantera/src/zeroD/Wall.cpp
index 4719b8ddc..c76885cb6 100644
--- a/Cantera/src/zeroD/Wall.cpp
+++ b/Cantera/src/zeroD/Wall.cpp
@@ -114,4 +114,58 @@ namespace CanteraZeroD {
         else
             m_surf[1]->setCoverages(m_rightcov.begin());
     }
-}        
+
+    void Wall::addSensitivityReaction(int leftright, int rxn) {
+        if (rxn < 0 || rxn >= m_chem[leftright]->nReactions()) 
+            throw CanteraError("Wall::addSensitivityReaction",
+                "Reaction number out of range ("+int2str(rxn)+")");
+        if (leftright == 0) {
+            m_pleft.push_back(rxn);
+            m_leftmult_save.push_back(1.0);
+            m_pname_left.push_back(m_chem[0]->reactionString(rxn));
+        }
+        else {
+            m_pright.push_back(rxn);
+            m_rightmult_save.push_back(1.0);
+            m_pname_right.push_back(m_chem[1]->reactionString(rxn));
+        }
+    }
+
+    void Wall::setSensitivityParameters(int lr, double* params) {
+        // process sensitivity parameters
+        int n, npar;
+        if (lr == 0) {
+            npar = m_pleft.size();
+            for (n = 0; n < npar; n++) {
+                m_leftmult_save[n] = m_chem[0]->multiplier(m_pleft[n]);
+                m_chem[0]->setMultiplier(m_pleft[n], 
+                    m_leftmult_save[n]*params[n]);
+            }
+        }
+        else {
+            npar = m_pright.size();
+            for (n = 0; n < npar; n++) {
+                m_rightmult_save[n] = m_chem[1]->multiplier(m_pright[n]);
+                m_chem[1]->setMultiplier(m_pright[n], 
+                    m_rightmult_save[n]*params[n]);
+            }            
+        }
+    }
+
+    void Wall::resetSensitivityParameters(int lr) {
+        int n, npar;
+        if (lr == 0) {
+            npar = m_pleft.size();
+            for (n = 0; n < npar; n++) {
+                m_chem[0]->setMultiplier(m_pleft[n], m_leftmult_save[n]);
+            }
+        }
+        else {
+            npar = m_pright.size();
+            for (n = 0; n < npar; n++) {
+                m_chem[1]->setMultiplier(m_pright[n], 
+                    m_rightmult_save[n]);
+            }            
+        }
+    }
+}
diff --git a/Cantera/src/zeroD/Wall.h b/Cantera/src/zeroD/Wall.h
index d56a3ab43..af25b44d7 100644
--- a/Cantera/src/zeroD/Wall.h
+++ b/Cantera/src/zeroD/Wall.h
@@ -135,6 +135,24 @@ namespace CanteraZeroD {
         void syncCoverages(int leftright);
 
 
+        int nSensParams(int lr) const { 
+            if (lr == 0) 
+                return m_pleft.size();
+            else
+                return m_pright.size();
+        }
+        void addSensitivityReaction(int leftright, int rxn);
+        string sensitivityParamID(int leftright, int p) { 
+            if (leftright == 0)
+                return m_pname_left[p]; 
+            else 
+                return m_pname_right[p];
+        }
+        void setSensitivityParameters(int lr, double* params);
+        void resetSensitivityParameters(int lr);
+        
+        //        int componentIndex(string nm) const;
+
     protected:
 
         //vector_fp m_coeffs;
@@ -150,6 +168,10 @@ namespace CanteraZeroD {
         Func1 *m_qf;
         vector_fp m_leftcov, m_rightcov;
 
+        vector_int m_pleft, m_pright;
+        vector_fp m_leftmult_save, m_rightmult_save;
+        vector<string> m_pname_left, m_pname_right;
+
     private:
 
     };