diff --git a/Cantera/src/CVode.cpp b/Cantera/src/CVode.cpp
index 14ec47083..f7d3c5d44 100755
--- a/Cantera/src/CVode.cpp
+++ b/Cantera/src/CVode.cpp
@@ -40,7 +40,7 @@ extern "C" {
         double* ydata = N_VDATA(y);
         double* ydotdata = N_VDATA(ydot);
         Cantera::FuncEval* f = (Cantera::FuncEval*)f_data;
-        f->eval(t, ydata, ydotdata);
+        f->eval(t, ydata, ydotdata, NULL);
     }
 
 
@@ -71,7 +71,7 @@ extern "C" {
             dy = 1.0/ewtdata[j];
             ydata[j] = ysave + dy;
             dy = ydata[j] - ysave;
-            func->eval(t, ydata, ydot);
+            func->eval(t, ydata, ydot, NULL);
             for (i=0; i < N; i++) {
                 col_j[i] = (ydot[i] - fydata[i])/dy;
             }
@@ -212,6 +212,8 @@ namespace Cantera {
                 m_iter, m_itol, &m_reltol,
                 &m_abstols, m_data, NULL, TRUE, m_iopt, 
                 m_ropt.begin(), NULL);
+            cout << "m_reltol = " << m_reltol << endl;
+            cout << "m_abstols = " << m_abstols << endl;
         }
 
         if (!m_cvode_mem) throw CVodeErr("CVodeMalloc failed.");
diff --git a/Cantera/src/CVodesIntegrator.cpp b/Cantera/src/CVodesIntegrator.cpp
index 56f66e482..d5350d2a1 100644
--- a/Cantera/src/CVodesIntegrator.cpp
+++ b/Cantera/src/CVodesIntegrator.cpp
@@ -1,11 +1,13 @@
 /**
- *  @file CVodeInt.cpp
+ *  @file CVodesIntegrator.cpp
  *
  */
 
 // Copyright 2001  California Institute of Technology
 
 #include "CVodesIntegrator.h"
+#include "stringUtils.h"
+
 #include <iostream>
 using namespace std;
 
@@ -24,6 +26,21 @@ inline static N_Vector nv(void* x) {
     return reinterpret_cast<N_Vector>(x);
 }
 
+namespace Cantera {
+    
+    class FuncData {
+    public:
+        FuncData(FuncEval* f, int npar = 0) {
+            m_pars.resize(npar, 1.0);
+            m_func = f;
+        }
+        virtual ~FuncData() {}
+        vector_fp m_pars;
+        FuncEval* m_func;
+    };
+}
+
+        
 extern "C" {
 
     /**
@@ -39,8 +56,18 @@ extern "C" {
         void *f_data) {
         double* ydata = NV_DATA_S(y); //N_VDATA(y);
         double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot);
-        Cantera::FuncEval* f = (Cantera::FuncEval*)f_data;
-        f->eval(t, ydata, ydotdata, NULL);
+        Cantera::FuncData* d = (Cantera::FuncData*)f_data;
+        Cantera::FuncEval* f = d->m_func;
+        //try {
+            if (d->m_pars.size() == 0)
+                f->eval(t, ydata, ydotdata, NULL);
+            else
+                f->eval(t, ydata, ydotdata, d->m_pars.begin());
+            //}
+            //catch (...) {
+            //Cantera::showErrors();
+            //Cantera::error("Teminating execution");
+            //}
     }
 
 }
@@ -64,9 +91,11 @@ namespace Cantera {
                            m_maxord(0),
                            m_reltol(1.e-9), 
                            m_abstols(1.e-15), 
+                                           m_reltolsens(1.0e-5),
+                                           m_abstolsens(1.0e-4),
                            m_nabs(0), 
                            m_hmax(0.0),
-                           m_maxsteps(20000)
+                                           m_maxsteps(20000), m_np(0)
     {
         //m_ropt.resize(OPT_SIZE,0.0);
         //m_iopt = new long[OPT_SIZE];
@@ -77,9 +106,15 @@ namespace Cantera {
     /// Destructor.
     CVodesIntegrator::~CVodesIntegrator()
     {   
-      if (m_cvode_mem) CVodeFree(m_cvode_mem);
+        if (m_cvode_mem) {
+            if (m_np > 0) 
+                CVodeSensFree(m_cvode_mem);
+            CVodeFree(m_cvode_mem);
+        }
       if (m_y) N_VDestroy_Serial(nv(m_y)); //N_VFree(nv(m_y));
       if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); //N_VFree(nv(m_abstol));
+      delete m_fdata;
+
       //delete[] m_iopt;
     }
     
@@ -109,6 +144,11 @@ namespace Cantera {
         m_abstols = abstol;
     }
 
+    void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
+        m_reltolsens = reltol; 
+        m_abstolsens = abstol;
+    }
+
     void CVodesIntegrator::setProblemType(int probtype) {
       m_type = probtype;
     }
@@ -151,6 +191,29 @@ namespace Cantera {
             throw CVodesErr("unknown iterator");
     }
 
+    void CVodesIntegrator::sensInit(double t0, FuncEval& func) {
+        m_np = func.nparams();
+        long int nv = func.neq();
+
+        doublereal* data;
+        int n, j;
+        m_yS = N_VNewVectorArray_Serial(m_np, nv);
+        for (n = 0; n < m_np; n++) {
+            data = NV_DATA_S(m_yS[n]);
+            for (j = 0; j < nv; j++) {
+                data[j] =0.0;
+            }
+        }
+        int flag;
+        flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
+        if (flag != CV_SUCCESS) 
+            throw CVodesErr("Error in CVodeSensMalloc");
+        vector_fp atol(m_np, m_abstolsens);
+        double rtol = m_reltolsens;
+        cout << "atol = " << atol[0] << " " << atol[m_np-1] << endl;
+        flag = CVodeSetSensTolerances(m_cvode_mem, CV_SS, rtol, atol.begin());
+    }
+
     void CVodesIntegrator::initialize(double t0, FuncEval& func) 
     {
         m_neq = func.neq();
@@ -166,11 +229,13 @@ namespace Cantera {
         // check abs tolerance array size
         if (m_itol == CV_SV && m_nabs < m_neq) 
             throw CVodesErr("not enough absolute tolerance values specified.");
-        func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
-
-        //m_iopt[MXSTEP] = m_maxsteps;
-        //m_iopt[MAXORD] = m_maxord;
-        //m_ropt[HMAX]   = m_hmax;
+        //try {
+            func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
+            //}
+            //catch (CanteraError) {
+            //showErrors();
+            //error("Teminating execution");
+            // }
 
         if (m_cvode_mem) CVodeFree(m_cvode_mem);
 
@@ -196,12 +261,17 @@ namespace Cantera {
             //      &m_abstols, m_data, NULL, TRUE, m_iopt, 
             //      m_ropt.begin(), NULL);
         }
-        if (flag == CV_MEM_FAIL) {
-            throw CVodesErr("Memory allocation failed.");
-        }
-        else if (flag == CV_ILL_INPUT) {
-            throw CVodesErr("Illegal value for CVodeMalloc input argument.");
+        if (flag != CV_SUCCESS) {
+            if (flag == CV_MEM_FAIL) {
+                throw CVodesErr("Memory allocation failed.");
+            }
+            else if (flag == CV_ILL_INPUT) {
+                throw CVodesErr("Illegal value for CVodeMalloc input argument.");
+            }
+            else 
+                throw CVodesErr("CVodeMalloc failed.");
         }
+        cout << "returned from CVodeMalloc. m_cvode_mem = " << m_cvode_mem << endl;
 
 
         if (m_type == DENSE + NOJAC) {
@@ -219,11 +289,20 @@ namespace Cantera {
         }
 
         // pass a pointer to func in m_data 
-        m_data = (void*)&func;
-        flag = CVodeSetFdata(m_cvode_mem, m_data);
+        m_fdata = new FuncData(&func, func.nparams());
+
+        //m_data = (void*)&func;
+
+        flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
         if (flag != CV_SUCCESS) 
             throw CVodesErr("CVodeSetFdata failed.");
 
+        if (func.nparams() > 0) {
+            sensInit(t0, func);
+            flag = CVodeSetSensParams(m_cvode_mem, m_fdata->m_pars.begin(), 
+                NULL, NULL);
+        }
+
         // set options
         if (m_maxord > 0)
             flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
@@ -237,16 +316,15 @@ namespace Cantera {
     void CVodesIntegrator::reinitialize(double t0, FuncEval& func) 
     {
         m_t0  = t0;
-        func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
+        //try {
+            func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
+            //}
+            //catch (CanteraError) {
+            //showErrors();
+            //error("Teminating execution");
+            //}
 
-        // set options
-        //        m_iopt[MXSTEP] = m_maxsteps;
-        //m_iopt[MAXORD] = m_maxord;
-        //m_ropt[HMAX]   = m_hmax;
-
-        //if (m_cvode_mem) CVodeFree(m_cvode_mem);
-
-        int result;
+        int result, flag;
         if (m_itol == CV_SV) {
             result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), 
                 m_itol, m_reltol,
@@ -257,8 +335,8 @@ namespace Cantera {
                 m_itol, m_reltol,
                 &m_abstols);
         }
-
-        if (result != 0) throw CVodesErr("CVReInit failed.");
+        cout << "returned from CVodeReInit. m_cvode_mem = " << m_cvode_mem << endl;
+        if (result != CV_SUCCESS) throw CVodesErr("CVReInit failed. result = "+int2str(result));
 
         if (m_type == DENSE + NOJAC) {
             long int N = m_neq;
@@ -274,11 +352,6 @@ namespace Cantera {
             throw CVodesErr("unsupported option");
         }
 
-        // pass a pointer to func in m_data 
-        m_data = (void*)&func;
-        long int flag = CVodeSetFdata(m_cvode_mem, m_data);
-        if (flag != CV_SUCCESS) 
-            throw CVodesErr("CVodeSetFdata failed.");
 
         // set options
         if (m_maxord > 0)
@@ -296,7 +369,10 @@ namespace Cantera {
 	flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
 	if (flag != CV_SUCCESS) 
 	  throw CVodesErr(" CVodes error encountered.");
-      }
+        if (m_np > 0) {
+            CVodeGetSens(m_cvode_mem, tout, m_yS);
+        }
+    }
 
     double CVodesIntegrator::step(double tout)
     {
@@ -310,10 +386,18 @@ namespace Cantera {
 
     int CVodesIntegrator::nEvals() const {
         long int ne;
-        return CVodeGetNumRhsEvals(m_cvode_mem, &ne);
+        CVodeGetNumRhsEvals(m_cvode_mem, &ne);
         return ne;
         //return m_iopt[NFE]; 
     }
+
+    double CVodesIntegrator::sensitivity(int k, int p) {
+        if (k < 0 || k >= m_neq) 
+            throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")");
+        if (p < 0 || p >= m_np)
+            throw CVodesErr("sensitivity: p out of range ("+int2str(p)+")");
+        return NV_Ith_S(m_yS[p],k); 
+    }
 }
 
 
diff --git a/Cantera/src/CVodesIntegrator.h b/Cantera/src/CVodesIntegrator.h
index cfb45d145..1753d2b44 100644
--- a/Cantera/src/CVodesIntegrator.h
+++ b/Cantera/src/CVodesIntegrator.h
@@ -23,8 +23,13 @@
 #include "ctexceptions.h"
 #include "ct_defs.h"
 
+#include <nvector.h>
+#include <nvector_serial.h>
+
 namespace Cantera {
 
+    class FuncData;
+
     /**
      * Exception thrown when a CVODES error is encountered.
      */
@@ -49,6 +54,7 @@ namespace Cantera {
         virtual ~CVodesIntegrator();
         virtual void setTolerances(double reltol, int n, double* abstol);
         virtual void setTolerances(double reltol, double abstol);
+        virtual void setSensitivityTolerances(double reltol, double abstol);
         virtual void setProblemType(int probtype);
         virtual void initialize(double t0, FuncEval& func);
         virtual void reinitialize(double t0, FuncEval& func);
@@ -65,8 +71,13 @@ namespace Cantera {
         virtual void setMinStepSize(double hmin);
         virtual void setMaxSteps(int nmax);
 
+        virtual int nSensParams() { return m_np; }
+        virtual double sensitivity(int k, int p);
+
     private:
 
+        void sensInit(double t0, FuncEval& func);
+
 	int m_neq;
         void* m_cvode_mem;
         double m_t0;
@@ -78,10 +89,14 @@ namespace Cantera {
         int m_maxord;
         double m_reltol;
         double m_abstols;
+        double m_reltolsens, m_abstolsens;
         int m_nabs;
         double m_hmax, m_hmin;
         int m_maxsteps;
-        void* m_data;
+        FuncData* m_fdata;
+        N_Vector*  m_yS;
+        int m_np;
+
     };
 
 }    // namespace
diff --git a/Cantera/src/ChemEquil.cpp b/Cantera/src/ChemEquil.cpp
index d2a7fb266..ec8fdeadf 100755
--- a/Cantera/src/ChemEquil.cpp
+++ b/Cantera/src/ChemEquil.cpp
@@ -352,6 +352,7 @@ namespace Cantera {
             m_p2 = new DensityCalculator<thermo_t>;
             break; 
         default:
+            endLogGroup();
             throw CanteraError("equilibrate","illegal property pair.");
         }
 
@@ -364,6 +365,7 @@ namespace Cantera {
         if (tempFixed) {
             double tfixed = s.temperature();
             if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) {
+                endLogGroup();
                 throw CanteraError("ChemEquil","Specified temperature ("
                     +fp2str(m_thermo->temperature())+" K) outside "
                     "valid range of "+fp2str(m_thermo->minTemp())+" K to "
@@ -573,6 +575,8 @@ namespace Cantera {
             fail++;
             if (fail > 3) {
                 addLogEntry("dampStep","Failed 3 times. Giving up.");
+                endLogGroup(); // iteration
+                endLogGroup(); // equilibrate
                 s.restoreState(state);
                 throw CanteraError("equilibrate",
                     "Cannot find an acceptable Newton damping coefficient.");
diff --git a/Cantera/src/EdgeKinetics.cpp b/Cantera/src/EdgeKinetics.cpp
index 146c48299..a2ca8a540 100644
--- a/Cantera/src/EdgeKinetics.cpp
+++ b/Cantera/src/EdgeKinetics.cpp
@@ -263,7 +263,10 @@ namespace Cantera {
             //cout << "i, beta = " << i << " " << m_beta[i] << endl;
             if (eamod != 0.0 && m_E[i] != 0.0) {
                 ea = GasConstant * m_E[i];
-                if (eamod + ea < 0.0) eamod = -ea;
+                if (eamod + ea < 0.0) {
+                    writelog("Warning: act energy mod too large");
+                    eamod = -ea;
+                }
                 kf[irxn] *= exp(-eamod*rrt);
             }
         }
diff --git a/Cantera/src/FalloffMgr.h b/Cantera/src/FalloffMgr.h
index ab6de82e4..fd62c9459 100755
--- a/Cantera/src/FalloffMgr.h
+++ b/Cantera/src/FalloffMgr.h
@@ -86,8 +86,10 @@ namespace Cantera {
          */
         void updateTemp(doublereal t, workPtr work) {
             int i;
-            for (i = 0; i < m_n; i++) m_falloff[i]->updateTemp(t, 
-                work + m_offset[i]);
+            for (i = 0; i < m_n; i++) {
+                m_falloff[i]->updateTemp(t, 
+                    work + m_offset[i]);
+            }
         }
 
         /**
diff --git a/Cantera/src/Func1.h b/Cantera/src/Func1.h
index bf0b7b125..5ec3ab4b8 100644
--- a/Cantera/src/Func1.h
+++ b/Cantera/src/Func1.h
@@ -32,20 +32,30 @@ namespace Cantera {
 
     /**
      * Base class for 'functor' classes that evaluate a function of
-     * one variable.
+     * one variable. 
      */
     class Func1 {
     public:
         Func1() {}
         virtual ~Func1() {}
+        /// Calls method eval to evaluate the function
         doublereal operator()(doublereal t) { return eval(t); }
+        /// Evaluate the function.
         virtual doublereal eval(doublereal t) { return 0.0; }
-    protected:
-        int m_n;
     private:
     };
 
 
+    /**
+     * A Gaussian. 
+     * \f[
+     * f(t) = A e^{-[(t - t_0)/\tau]^2}
+     * \f]
+     * where \f[ \tau = \frac{fwhm}{2\sqrt{\ln 2}} \f]
+     * @param A peak value
+     * @param t0 offset
+     * @param fwhm full width at half max  
+     */
     class Gaussian : public Func1 {
     public:
         Gaussian(double A, double t0, double fwhm) {
@@ -86,23 +96,6 @@ namespace Cantera {
             return r;
         }
 
-//         virtual string show(doublereal t) {
-//             int n;
-//             string s = "";
-//             doublereal r = m_c[m_n-1];
-//             for (n = m_n-1; n >= 0; n--) {
-//                 s += fp2str(m_c[n]);
-//                 if (n > 0) s += "*x";
-//                 if (n > 1) s += "^"+int2str(n);
-//                 if (n > 0) {
-//                     if (m_c[n] < 0.0) s += " - ";
-//                     else 
-//                 r *= t;
-//                 r += m_c[m_n - n - 1];
-//             }
-//             return r;
-//         }            
-
     protected:
         int m_n;
         vector_fp m_c;
@@ -110,7 +103,12 @@ namespace Cantera {
 
 
     /**
-     * Fourier cosine/sin series.
+     * Fourier cosine/sine series.
+     * 
+     * \f[
+     * f(t) = \frac{A_0}{2} + 
+     * \sum_{n=1}^N A_n \cos (n \omega t) + B_n \sin (n \omega t)
+     * \f]
      */
     class Fourier1 : public Func1 {
     public:
@@ -146,6 +144,9 @@ namespace Cantera {
 
     /**
      * Sum of Arrhenius terms.
+     * \f[
+     * f(T) = \sum_{n=1}^N A_n T^b_n \exp(-E_n/T)
+     * \f]
      */
     class Arrhenius1 : public Func1 {
     public:
diff --git a/Cantera/src/FuncEval.h b/Cantera/src/FuncEval.h
index f64c48ad9..0904af0a2 100755
--- a/Cantera/src/FuncEval.h
+++ b/Cantera/src/FuncEval.h
@@ -37,8 +37,9 @@ namespace Cantera {
          * @param t time. (input)
          * @param y solution vector. (input)
          * @param ydot rate of change of solution vector. (output)
+         * @param p parameter vector
          */
-	virtual void eval(double t, double* y, double* ydot)=0;
+	virtual void eval(double t, double* y, double* ydot, double* p)=0;
 
         /**
          * Fill the solution vector with the initial conditions 
@@ -49,6 +50,9 @@ namespace Cantera {
         /** Number of equations. */
         virtual int neq()=0;
 
+        /// Number of parameters.
+        virtual int nparams() { return 0; }
+
     protected:
 
     private:
diff --git a/Cantera/src/GRI_30_Kinetics.cpp b/Cantera/src/GRI_30_Kinetics.cpp
index dcc9bd893..2fc8ef8e5 100755
--- a/Cantera/src/GRI_30_Kinetics.cpp
+++ b/Cantera/src/GRI_30_Kinetics.cpp
@@ -25,10 +25,6 @@
 using namespace std;
 
 
-#ifdef HAVE_INTEL_MKL
-#include "mkl_vml.h"
-#endif
-
 namespace Cantera {
 
     /**
@@ -44,11 +40,9 @@ namespace Cantera {
             doublereal logT = log(T);
             m_kdata->m_logc_ref = m_kdata->m_logp_ref - logT;
             update_rates(T, logT, m_kdata->m_rfn.begin());
-
-            m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin()); 
+            m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin());
             m_falloff_high_rates.update(T, logT, m_kdata->m_rfn_high.begin());
             m_falloffn.updateTemp(T, m_kdata->falloff_work.begin());
-
             m_kdata->m_temp = T;
             gri30_updateKc();
             m_kdata->m_ROP_ok = false;
diff --git a/Cantera/src/ImplicitChem.cpp b/Cantera/src/ImplicitChem.cpp
index 12d6d3609..5a05e2f76 100755
--- a/Cantera/src/ImplicitChem.cpp
+++ b/Cantera/src/ImplicitChem.cpp
@@ -16,7 +16,7 @@
 #endif
 
 #include "ImplicitChem.h"
-#include "CVode.h"
+#include "Integrator.h"
 
 namespace Cantera {
 
@@ -24,7 +24,7 @@ namespace Cantera {
         : FuncEval(), m_kin(&kin), m_thermo(&therm), m_integ(0),
           m_atol(1.e-15), m_rtol(1.e-7), m_maxstep(0.0), m_energy(false)
     {
-        m_integ = new CVodeInt;
+        m_integ = newIntegrator("CVODE"); //CVodeInt;
         //m_mix = &kin.phase();
         m_wt = m_thermo->molecularWeights();
 
@@ -74,7 +74,8 @@ namespace Cantera {
     /**
      * Called by the integrator to evaluate ydot given y at time 'time'.
      */
-    void ImplicitChem::eval(doublereal time, doublereal* y, doublereal* ydot) 
+    void ImplicitChem::eval(doublereal time, doublereal* y, 
+        doublereal* ydot, doublereal* p) 
     {
         updateState(y);   // synchronize the mixture state with y
         m_thermo->setPressure(m_press);
diff --git a/Cantera/src/ImplicitChem.h b/Cantera/src/ImplicitChem.h
index e20d4f8b9..96180bdee 100755
--- a/Cantera/src/ImplicitChem.h
+++ b/Cantera/src/ImplicitChem.h
@@ -17,7 +17,7 @@
 #endif
 
 #include "FuncEval.h"
-#include "CVode.h"
+#include "Integrator.h"
 #include "Kinetics.h"
 #include "ThermoPhase.h"
 
@@ -83,7 +83,8 @@ namespace Cantera {
 
         // overloaded methods of class FuncEval
         virtual int neq() { return m_nsp; }
-	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);
 
diff --git a/Cantera/src/ImplicitSurfChem.cpp b/Cantera/src/ImplicitSurfChem.cpp
index 4e348c190..b6a1b5367 100755
--- a/Cantera/src/ImplicitSurfChem.cpp
+++ b/Cantera/src/ImplicitSurfChem.cpp
@@ -17,7 +17,8 @@
 #endif
 
 #include "ImplicitSurfChem.h"
-#include "CVode.h"
+
+#include "Integrator.h"
 
 namespace Cantera {
 
@@ -41,7 +42,7 @@ namespace Cantera {
             nt = k[n]->nTotalSpecies();
             if (nt > ntmax) ntmax = nt;
         }
-        m_integ = new CVodeInt;
+        m_integ = newIntegrator("CVODE");// CVodeInt;
 
         // use backward differencing, with a full Jacobian computed
         // numerically, and use a Newton linear iterator
@@ -88,12 +89,13 @@ namespace Cantera {
      * Called by the integrator to evaluate ydot given y at time 'time'.
      */
     void ImplicitSurfChem::eval(doublereal time, doublereal* y, 
-        doublereal* ydot) 
+        doublereal* ydot, doublereal* p) 
     {
+        int n;
         updateState(y);   // synchronize the surface state(s) with y
         doublereal rs0, sum;
         int loc, k, kstart;
-        for (int n = 0; n < m_nsurf; n++) {
+        for (n = 0; n < m_nsurf; n++) {
             rs0 = 1.0/m_surf[n]->siteDensity();
             m_kin[n]->getNetProductionRates(m_work.begin());
             kstart = m_kin[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
diff --git a/Cantera/src/ImplicitSurfChem.h b/Cantera/src/ImplicitSurfChem.h
index b5dc4b694..799f2cdf8 100755
--- a/Cantera/src/ImplicitSurfChem.h
+++ b/Cantera/src/ImplicitSurfChem.h
@@ -20,7 +20,7 @@
 #endif
 
 #include "FuncEval.h"
-#include "CVode.h"
+#include "Integrator.h"
 #include "InterfaceKinetics.h"
 #include "SurfPhase.h"
 
@@ -82,7 +82,8 @@ namespace Cantera {
 
         // 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);
 
diff --git a/Cantera/src/Integrator.h b/Cantera/src/Integrator.h
index 6a0413520..1fb48d270 100755
--- a/Cantera/src/Integrator.h
+++ b/Cantera/src/Integrator.h
@@ -85,6 +85,9 @@ namespace Cantera {
         virtual void setTolerances(doublereal reltol, doublereal abstol)
             { warn("setTolerances"); }
 
+        virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol)
+            { warn("setSensitivityTolerances"); }
+
         /**
          * Set problem type.
          */
@@ -157,6 +160,13 @@ namespace Cantera {
         virtual void setMaxSteps(int nmax) 
             { warn("setMaxStep"); }
 
+        virtual int nSensParams() 
+            { warn("nSensParams()"); return 0; }
+
+        virtual double sensitivity(int k, int p) { 
+            warn("sensitivity"); return 0.0;
+        }
+        
     private:
 
         doublereal m_dummy;
@@ -167,6 +177,9 @@ namespace Cantera {
 
     };
 
+    // defined in ODE_integrators.cpp
+    Integrator* newIntegrator(string itype);
+
 }    // namespace
 
 #endif 
diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in
index db4c7c9ac..4a3edaa16 100755
--- a/Cantera/src/Makefile.in
+++ b/Cantera/src/Makefile.in
@@ -18,6 +18,8 @@ CANTERA_LIB = @buildlib@/libcantera.a
 CXX_FLAGS   = @CXXFLAGS@ $(CXX_OPT) $(LOCAL_DEFNS)
 EXT         = ../../ext
 do_ranlib = @DO_RANLIB@
+USE_SUNDIALS = @use_sundials@
+SUNDIALS_INC = @sundials_include@
 
 #----------------------
 #  kernel components
@@ -87,8 +89,8 @@ RPATH      = $(RPATH_OBJ)
 
 
 # solvers
-SOLVERS_OBJ = CVode.o  BandMatrix.o
-SOLVERS_H   = CVode.h  BandMatrix.h Integrator.h
+SOLVERS_OBJ = ODE_integrators.o  BandMatrix.o
+SOLVERS_H   = BandMatrix.h Integrator.h
 SOLVERS    = $(SOLVERS_OBJ)
 
 # 1D flow models
@@ -183,6 +185,9 @@ ALL_H = $(BASE_H) $(THERMO_H) $(KINETICS_H) $(HETEROKIN_H) \
 .cpp.o:
 	@CXX@ -c $< $(CXX_INCLUDES) $(CXX_FLAGS) 
 
+ODE_integrators.o:
+	@CXX@ -c ODE_integrators.cpp $(CXX_INCLUDES) $(SUNDIALS_INC) $(CXX_FLAGS) 
+
 lib:  $(OBJ_LIB)
 	$(RM)     $(CANTERA_LIB)
 	@ARCHIVE@ $(CANTERA_LIB) *.o > /dev/null
diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp
index 45ba0f96e..50d0f1607 100644
--- a/Cantera/src/MultiPhase.cpp
+++ b/Cantera/src/MultiPhase.cpp
@@ -534,7 +534,6 @@ namespace Cantera {
                 "No convergence for T");
         }
         else if (XY == SP) {
-            writelog("SP\n");
             s0 = entropy();
             start = true;
             Tlow = 1.0; // m_Tmin;      // lower bound on T
@@ -610,7 +609,6 @@ namespace Cantera {
                 "No convergence for T");
         }
 
-
 //         else if (XY == SP) {
 //             if (loglevel > 0) {
 //                 addLogEntry("problem type","fixed S,P");
diff --git a/Cantera/src/NasaPoly1.h b/Cantera/src/NasaPoly1.h
index 0fc4983f3..dcfce287a 100755
--- a/Cantera/src/NasaPoly1.h
+++ b/Cantera/src/NasaPoly1.h
@@ -13,6 +13,7 @@
 #ifndef CT_NASAPOLY1_H
 #define CT_NASAPOLY1_H
 
+#include "global.h"
 #include "SpeciesThermoInterpType.h"
 
 namespace Cantera {
@@ -136,6 +137,8 @@ namespace Cantera {
             cp_R[m_index] = cp;
             h_RT[m_index] = h;
             s_R[m_index] = s;
+            //writelog("NASA1: for species "+int2str(m_index)+", h_RT = "+
+            //    fp2str(h)+"\n");
         }
 
 	/**
diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h
index 31aba5a18..ac26c6171 100755
--- a/Cantera/src/NasaThermo.h
+++ b/Cantera/src/NasaThermo.h
@@ -17,6 +17,7 @@
 #include "NasaPoly1.h"
 #include "speciesThermoTypes.h"
 #include "polyfit.h"
+#include "global.h"
 
 namespace Cantera {
 
@@ -70,6 +71,8 @@ namespace Cantera {
 			     doublereal minTemp, doublereal maxTemp,
 			     doublereal refPressure) { 
 
+            //writelog("installing NASA thermo for species "+name+"\n");
+            //writelog(" index = "+int2str(index)+"\n");
             int imid = int(c[0]);       // midpoint temp converted to integer
             int igrp = m_index[imid];   // has this value been seen before?
             if (igrp == 0) {            // if not, prepare new group
diff --git a/Cantera/src/ODE_integrators.cpp b/Cantera/src/ODE_integrators.cpp
index b3d34d5b3..fbf53d89f 100644
--- a/Cantera/src/ODE_integrators.cpp
+++ b/Cantera/src/ODE_integrators.cpp
@@ -7,19 +7,19 @@
 #include "CVode.cpp"
 #endif
 
-// namespace Cantera {
+namespace Cantera {
 
-//     Integrator* newIntegrator(string itype) {
-//         if (itype == "CVODE") {
-// #ifdef HAS_SUNDIALS
-//             return new CVodesIntegrator();
-// #else
-//             return new CVodeInt();
-// #endif
-//         }
-//         else {
-//             throw CanteraError("newIntegrator",
-//                 "unknown ODE integrator: "+itype);
-//         }
-//     }
-// }
+    Integrator* newIntegrator(string itype) {
+        if (itype == "CVODE") {
+#ifdef HAS_SUNDIALS
+            return new CVodesIntegrator();
+#else
+            return new CVodeInt();
+#endif
+        }
+        else {
+            throw CanteraError("newIntegrator",
+                "unknown ODE integrator: "+itype);
+        }
+    }
+}
diff --git a/Cantera/src/ReactionPath.cpp b/Cantera/src/ReactionPath.cpp
index 3241d1896..de6fc7808 100755
--- a/Cantera/src/ReactionPath.cpp
+++ b/Cantera/src/ReactionPath.cpp
@@ -450,7 +450,7 @@ namespace Cantera {
             }
         }
         s << " label = " << "\"" << "Scale = " 
-          << flmax << "\";" << endl; //\\l\\l created with Cantera (www.cantera.org)\\l\";"
+          << flmax << "\\l " << title << "\";" << endl; //created with Cantera (www.cantera.org)\\l\";"
         s  << " fontname = \""+m_font+"\";" << endl << "}" << endl;
     }
 
@@ -658,21 +658,72 @@ namespace Cantera {
         g.fmt(out, m_elementSymbols);
     }
 
+    void ReactionPathBuilder::findElements(Kinetics& kin) {
+
+        string ename;
+        m_enamemap.clear();
+        m_nel = 0;
+        int i, np = kin.nPhases();
+        ThermoPhase* p;
+        map<string, int> enamemap;
+        for (i = 0; i < np; i++) {
+            p = &kin.thermo(i);
+            // iterate over the elements in this phase
+            int m, nel = p->nElements();
+            for (m = 0; m < nel; m++) {
+                ename = p->elementName(m);
+
+                // if no entry is found for this element name, then
+                // it is a new element. In this case, add the name
+                // to the list of names, increment the element count, 
+                // and add an entry to the name->(index+1) map.
+                if (m_enamemap.find(ename) == m_enamemap.end()) {
+                    m_enamemap[ename] = m_nel + 1;
+                    m_elementSymbols.push_back(ename);
+                    m_nel++;
+                }
+            }
+        }
+        m_atoms.resize(kin.nTotalSpecies(), m_nel, 0.0);
+        string sym;
+        int k, ip, nsp, mlocal, kp, m;
+        // iterate over the elements
+        for (m = 0; m < m_nel; m++) {
+            sym = m_elementSymbols[m];
+            k = 0;
+            // iterate over the phases
+            for (ip = 0; ip < np; ip++) {
+                phase_t* p = &kin.thermo(ip);
+                nsp = p->nSpecies();
+                mlocal = p->elementIndex(sym);    
+                for (kp = 0; kp < nsp; kp++) {
+                    if (mlocal >= 0) {
+                        m_atoms(k, m) = p->nAtoms(kp, mlocal);
+                    }
+                    k++;
+                }
+            }
+        }
+    }
+
+
+
     int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) {
         //m_warn.clear();
         m_transfer.clear();
 
-        const Kinetics::thermo_t& ph = kin.thermo();
-
-        m_nel = ph.nElements();
-        m_ns = ph.nSpecies();
-        m_nr = kin.nReactions();
+        //const Kinetics::thermo_t& ph = kin.thermo();
 
         m_elementSymbols.clear();
+        findElements(kin);
+        //m_nel = ph.nElements();
+        m_ns = kin.nTotalSpecies(); //ph.nSpecies();
+        m_nr = kin.nReactions();
+
         int m, i;
-        for (m = 0; m < m_nel; m++) {
-            m_elementSymbols.push_back(ph.elementName(m));
-        }
+        //for (m = 0; m < m_nel; m++) {
+        //    m_elementSymbols.push_back(ph.elementName(m));
+        //}
 
         // all reactants / products, even ones appearing on both sides
         // of the reaction
@@ -736,7 +787,7 @@ namespace Cantera {
             for (n = 0; n < nrnet; n++) {
                 k = m_reac[i][n];
                 for (int m = 0; m < m_nel; m++) {
-                    m_elatoms(m,i) += ph.nAtoms(k,m);
+                    m_elatoms(m,i) += m_atoms(k,m); //ph.nAtoms(k,m);
                 }
             }            
         }
@@ -746,7 +797,7 @@ namespace Cantera {
         m_sgroup.resize(m_ns);
         int j;
         for (j = 0; j < m_ns; j++) {
-            for (int m = 0; m < m_nel; m++) comp[m] = int(ph.nAtoms(j,m));
+            for (int m = 0; m < m_nel; m++) comp[m] = int(m_atoms(j,m)); //ph.nAtoms(j,m));
             m_sgroup[j] = Group(comp);
         }
 
@@ -768,10 +819,11 @@ namespace Cantera {
                 nar = 0;
                 nap = 0;
                 for (j = 0;  j < nr; j++) {
-                    if (ph.nAtoms(m_reac[i][j],m) > 0) nar++;
+                    //                    if (ph.nAtoms(m_reac[i][j],m) > 0) nar++;
+                    if (m_atoms(m_reac[i][j],m) > 0) nar++;
                 }
                 for (j = 0;  j < np; j++) {
-                    if (ph.nAtoms(m_prod[i][j],m) > 0) nap++;
+                    if (m_atoms(m_prod[i][j],m) > 0) nap++;
                 }
                 if (nar > 1 && nap > 1) {
                     m_determinate[i] = false; break;
@@ -812,19 +864,19 @@ namespace Cantera {
         doublereal threshold = 0.0;
         bool fwd_incl, rev_incl, force_incl;
 
-        const Kinetics::thermo_t& ph = s.thermo();
-        int m = ph.elementIndex(element);
+        //        const Kinetics::thermo_t& ph = s.thermo();
+        int m = m_enamemap[element]-1; //ph.elementIndex(element);
 
         r.element = element;
         if (m < 0) return -1;
         
         //int k;
-        int kk = ph.nSpecies();
+        int kk = s.nTotalSpecies();
 
         s.getFwdRatesOfProgress(m_ropf.begin());
         s.getRevRatesOfProgress(m_ropr.begin());
 
-        ph.getMoleFractions(m_x.begin());
+        //ph.getMoleFractions(m_x.begin());
 
         //doublereal sum = 0.0;
         //for (k = 0; k < kk; k++) {
@@ -869,7 +921,7 @@ namespace Cantera {
                         revlabel = "";
                         for (l = 0; l < np; l++) {
                             if (l != kp) 
-                                revlabel += " + "+ ph.speciesName(m_prod[i][l]);
+                                revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][l]);
                         }
                         if (s.reactionType(i) == THREE_BODY_RXN)
                             revlabel += " + M "; 
@@ -882,8 +934,8 @@ namespace Cantera {
                         // element m, and both are allowed to appear in
                         // the diagram
 
-                        if ((kkr != kkp) && (ph.nAtoms(kkr,m) > 0 
-                            &&  ph.nAtoms(kkp,m) > 0)
+                        if ((kkr != kkp) && (m_atoms(kkr,m) > 0 
+                            &&  m_atoms(kkp,m) > 0)
                             && status[kkr] >= 0 && status[kkp] >= 0) 
                         {
 
@@ -894,8 +946,8 @@ namespace Cantera {
                             // reactant species was the source of a
                             // given m-atom in the product
  
-                            if ( (ph.nAtoms(kkp,m) < m_elatoms(m, i)) && 
-                                (ph.nAtoms(kkr,m) < m_elatoms(m, i)) ) 
+                            if ( (m_atoms(kkp,m) < m_elatoms(m, i)) && 
+                                (m_atoms(kkr,m) < m_elatoms(m, i)) ) 
                             {
                                 map<int, map<int, Group> >& g = m_transfer[i];
                                 if (g.empty()) {
@@ -926,7 +978,7 @@ namespace Cantera {
                             // the same expression.
 
                             else {
-                                f = ph.nAtoms(kkp,m) * ph.nAtoms(kkr,m) / m_elatoms(m, i);
+                                f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i);
                             }
 
                             fwd = ropf*f;
@@ -940,10 +992,10 @@ namespace Cantera {
                             if (fwd_incl || rev_incl) 
                             {
                                 if (!r.hasNode(kkr)) {
-                                    r.addNode(kkr, ph.speciesName(kkr), m_x[kkr]);
+                                    r.addNode(kkr, s.kineticsSpeciesName(kkr), m_x[kkr]);
                                 }
                                 if (!r.hasNode(kkp)) {
-                                    r.addNode(kkp, ph.speciesName(kkp), m_x[kkp]);
+                                    r.addNode(kkp, s.kineticsSpeciesName(kkp), m_x[kkp]);
                                 }
                             }
                             if (fwd_incl) {
diff --git a/Cantera/src/ReactionPath.h b/Cantera/src/ReactionPath.h
index dbaa397d4..09b6597a7 100755
--- a/Cantera/src/ReactionPath.h
+++ b/Cantera/src/ReactionPath.h
@@ -252,6 +252,7 @@ namespace Cantera {
         void writeGroup(ostream& out, const Group& g);
 
     protected:
+        void findElements(Kinetics& kin);
 
         int m_nr;
         int m_ns;
@@ -268,6 +269,8 @@ namespace Cantera {
         //        map<int, int> m_warn;
         map<int, map<int, map<int, Group> > >  m_transfer;
         vector<bool> m_determinate;
+        Array2D m_atoms;
+        map<string,int> m_enamemap;
     };
 
 }
diff --git a/Cantera/src/SpeciesThermoFactory.cpp b/Cantera/src/SpeciesThermoFactory.cpp
index eb5c59982..4d58ac61b 100755
--- a/Cantera/src/SpeciesThermoFactory.cpp
+++ b/Cantera/src/SpeciesThermoFactory.cpp
@@ -72,6 +72,7 @@ namespace Cantera {
 	  iother = 1;
 	}
 	if (iother) {
+            writelog("returning new GeneralSpeciesThermo");
 	  return new GeneralSpeciesThermo();
 	}
         return newSpeciesThermo(NASA*inasa
@@ -90,6 +91,7 @@ namespace Cantera {
 	  }
         }
 	if (iother) {
+            writelog("returning new GeneralSpeciesThermo");
 	  return new GeneralSpeciesThermo();
 	}
         return newSpeciesThermo(NASA*inasa
@@ -110,6 +112,7 @@ namespace Cantera {
 	  }
         }
 	if (iother) {
+            writelog("returning new GeneralSpeciesThermo");
 	  return new GeneralSpeciesThermo();
 	}
         return newSpeciesThermo(NASA*inasa
diff --git a/Cantera/src/ThermoPhase.cpp b/Cantera/src/ThermoPhase.cpp
index 2492df091..9288e5302 100644
--- a/Cantera/src/ThermoPhase.cpp
+++ b/Cantera/src/ThermoPhase.cpp
@@ -123,13 +123,14 @@ namespace Cantera {
             dt = (u - intEnergy_mass())/cv_mass();
             if (dt > 100.0) dt = 100.0;
             else if (dt < -100.0) dt = -100.0;
-            setTemperature(temperature() + dt);
+            setTemperature(temperature() + 0.5*dt);
             if (fabs(dt) < tol) {
                 return;
             }
         }
         throw CanteraError("setState_UV",
             "no convergence. dt = " + fp2str(dt)+"\n"
+            +"tol = "+fp2str(tol)+"\n"
             +"u = "+fp2str(u)+" v = "+fp2str(v)+"\n");
     }
 
diff --git a/Cantera/src/ctvector.h b/Cantera/src/ctvector.h
index 6df9de3e0..9e48565d8 100755
--- a/Cantera/src/ctvector.h
+++ b/Cantera/src/ctvector.h
@@ -22,6 +22,10 @@
 #include <iostream>
 #include "config.h"
 
+#ifdef DEBUG_MODE
+#include "ctexceptions.h"
+#endif
+
 namespace ct {
 
     /**
@@ -44,7 +48,13 @@ namespace ct {
         ctvector_fp operator=(const ctvector_fp& x);
         virtual ~ctvector_fp();
 
-        value_type operator[](size_t n) const { return _data[n]; }
+        value_type operator[](size_t n) const {
+#ifdef DEBUG_MODE
+            if (n < 0 || n >= _size)
+                throw CanteraError("ctvector_fp","index out of range");
+#endif
+            return _data[n]; 
+        }
         value_type& operator[](size_t n) { return _data[n]; }
 
         void resize(size_t n);
diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp
index 946c9174d..798cfe524 100755
--- a/Cantera/src/misc.cpp
+++ b/Cantera/src/misc.cpp
@@ -772,10 +772,11 @@ namespace Cantera {
         if (title != "" && title != __app->loggroups.back()) {
             writelog("Logfile error."
                 "\n   beginLogGroup: "+ __app->loggroups.back()+
-                "\n   endLogGroup;   "+title+"\n");
+                "\n   endLogGroup:   "+title+"\n");
+            cout << "calling write_logfile..." << endl;
             write_logfile("logerror"); 
-            __app->loggroups.clear();
-            __app->loglevels.clear();
+            //__app->loggroups.clear();
+            //__app->loglevels.clear();
         }
         else if (__app->loggroups.size() == 1) {
             write_logfile(__app->loggroups.back()+"_log"); 
@@ -796,7 +797,9 @@ namespace Cantera {
     /// file will be overwritten.  will be appended to the name.
     /// @ingroup HTML_logs
     void write_logfile(string file) {
-        if (!__app->xmllog) return;
+        if (!__app->xmllog) {
+            return;
+        }
         string::size_type idot = file.rfind('.');
         string ext = "";
         string nm = file;
@@ -827,6 +830,7 @@ namespace Cantera {
         // Now we have a file name that does not correspond to any 
         // existing file. Open it as an output stream, and dump the 
         // XML (HTML) tree to it.
+
         if (__app->xmllog) {
             ofstream f(fname.c_str());
             // go to the top of the tree, and write it all.
diff --git a/Cantera/src/transport/SolidTransport.cpp b/Cantera/src/transport/SolidTransport.cpp
index 95cceef5f..da9ba8b12 100644
--- a/Cantera/src/transport/SolidTransport.cpp
+++ b/Cantera/src/transport/SolidTransport.cpp
@@ -62,7 +62,7 @@ namespace Cantera {
 
     /**
      * Compute the mobilities of the species from the diffusion coefficients, 
-     * usind the Einstein relation.
+     * using the Einstein relation.
      */
     void SolidTransport::getMobilities(doublereal* mobil) {
         int k;
diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h
index 489ada4c6..39fc37305 100755
--- a/Cantera/src/transport/TransportBase.h
+++ b/Cantera/src/transport/TransportBase.h
@@ -34,6 +34,7 @@ namespace Cantera {
     const int cUserTransport       = 500;
     const int cFtnTransport        = 600;
 
+    // forward reference
     class XML_Writer;
 
 
@@ -143,16 +144,26 @@ namespace Cantera {
 
         /** 
          * Get the molar fluxes [kmol/m^2/s], given the thermodynamic
-         * state at two nearby points.  @param state1 Array of
-         * temperature, density, and mass fractions for state 1.
+         * state at two nearby points. 
+         * @param state1 Array of temperature, density, and mass
+         * fractions for state 1.
          * @param state2 Array of temperature, density, and mass
-         * fractions for state 2.  @param delta Distance from state 1
-         * to state 2 (m).
+         * fractions for state 2.  
+         * @param delta Distance from state 1 to state 2 (m).
          */ 
         virtual void getMolarFluxes(const doublereal* state1,
             const doublereal* state2, doublereal delta, 
             doublereal* fluxes) { err("getMolarFluxes"); }
 
+        /** 
+         * Get the mass fluxes [kg/m^2/s], given the thermodynamic
+         * state at two nearby points. 
+         * @param state1 Array of temperature, density, and mass
+         * fractions for state 1.
+         * @param state2 Array of temperature, density, and mass
+         * fractions for state 2.  
+         * @param delta Distance from state 1 to state 2 (m).
+         */ 
         virtual void getMassFluxes(const doublereal* state1,
             const doublereal* state2, doublereal delta, 
             doublereal* fluxes) { err("getMassFluxes"); }