diff --git a/include/cantera/zeroD/IdealGasReactor.h b/include/cantera/zeroD/IdealGasReactor.h
new file mode 100644
index 000000000..5452f517b
--- /dev/null
+++ b/include/cantera/zeroD/IdealGasReactor.h
@@ -0,0 +1,50 @@
+/**
+ *  @file IdealGasReactor.h
+ */
+
+// Copyright 2001  California Institute of Technology
+
+#ifndef CT_IDEALGASREACTOR_H
+#define CT_IDEALGASREACTOR_H
+
+#include "Reactor.h"
+#include "cantera/kinetics/Kinetics.h"
+
+namespace Cantera
+{
+
+/**
+ * Class IdealGasReactor is a class for stirred reactors that is specifically
+ * optimized for ideal gases. In this formulation, temperature replaces the
+ * total internal energy as a state variable.
+ */
+class IdealGasReactor : public Reactor
+{
+public:
+    IdealGasReactor() {}
+
+    virtual int type() const {
+        return IdealGasReactorType;
+    }
+
+    virtual void setThermoMgr(ThermoPhase& thermo);
+
+    virtual void getInitialConditions(doublereal t0, size_t leny,
+                                      doublereal* y);
+
+    virtual void initialize(doublereal t0 = 0.0);
+
+    virtual void evalEqs(doublereal t, doublereal* y,
+                         doublereal* ydot, doublereal* params);
+
+    virtual void updateState(doublereal* y);
+
+    virtual size_t componentIndex(const std::string& nm) const;
+
+protected:
+    vector_fp m_uk; //!< Species molar internal energies
+};
+
+}
+
+#endif
diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h
index 7112717b5..e6aa7efae 100644
--- a/include/cantera/zeroD/ReactorBase.h
+++ b/include/cantera/zeroD/ReactorBase.h
@@ -64,7 +64,7 @@ public:
      * a pointer to this substance is stored, and as the integration
      * proceeds, the state of the substance is modified.
      */
-    void setThermoMgr(thermo_t& thermo);
+    virtual void setThermoMgr(thermo_t& thermo);
 
     //! Connect an inlet FlowDevice to this reactor
     void addInlet(FlowDevice& inlet);
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index aa8b5ed64..2740ed352 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -228,6 +228,11 @@ cdef class ConstPressureReactor(Reactor):
     reactor_type = "ConstPressureReactor"
 
 
+cdef class IdealGasReactor(Reactor):
+    """ A constant volume, zero-dimensional reactor for ideal gas mixtures. """
+    reactor_type = "IdealGasReactor"
+
+
 cdef class FlowReactor(Reactor):
     """
     A steady-state plug flow reactor with constant cross sectional area.
diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py
index 1e40d661c..d577a3442 100644
--- a/interfaces/cython/cantera/test/test_reactor.py
+++ b/interfaces/cython/cantera/test/test_reactor.py
@@ -7,6 +7,8 @@ from . import utilities
 
 
 class TestReactor(utilities.CanteraTest):
+    reactorClass = ct.Reactor
+
     def make_reactors(self, independent=True, n_reactors=2,
                       T1=300, P1=101325, X1='O2:1.0',
                       T2=300, P2=101325, X2='O2:1.0'):
@@ -15,7 +17,7 @@ class TestReactor(utilities.CanteraTest):
 
         self.gas1 = ct.Solution('h2o2.xml')
         self.gas1.TPX = T1, P1, X1
-        self.r1 = ct.Reactor(self.gas1)
+        self.r1 = self.reactorClass(self.gas1)
         self.net.add_reactor(self.r1)
 
         if independent:
@@ -25,7 +27,7 @@ class TestReactor(utilities.CanteraTest):
 
         if n_reactors >= 2:
             self.gas2.TPX = T2, P2, X2
-            self.r2 = ct.Reactor(self.gas2)
+            self.r2 = self.reactorClass(self.gas2)
             self.net.add_reactor(self.r2)
 
     def add_wall(self, **kwargs):
@@ -33,7 +35,7 @@ class TestReactor(utilities.CanteraTest):
         return self.w
 
     def test_insert(self):
-        R = ct.Reactor()
+        R = self.reactorClass()
         f1 = lambda r: r.T
         f2 = lambda r: r.kinetics.net_production_rates
         self.assertRaises(Exception, f1, R)
@@ -470,6 +472,10 @@ class TestReactor(utilities.CanteraTest):
         self.assertNear(p1a, p1b)
         self.assertNear(p2a, p2b)
 
+class TestIdealGasReactor(TestReactor):
+    reactorClass = ct.IdealGasReactor
+
+
 class TestWellStirredReactorIgnition(utilities.CanteraTest):
     """ Ignition (or not) of a well-stirred reactor """
     def setup(self, T0, P0, mdot_fuel, mdot_ox):
@@ -486,7 +492,7 @@ class TestWellStirredReactorIgnition(utilities.CanteraTest):
 
         # reactor, initially filled with N2
         self.gas.TPX = T0, P0, "N2:1.0"
-        self.combustor = ct.Reactor(self.gas)
+        self.combustor = ct.IdealGasReactor(self.gas)
         self.combustor.volume = 1.0
 
         # outlet
@@ -587,7 +593,7 @@ class TestConstPressureReactor(utilities.CanteraTest):
         self.gas1.TPX = T0, P0, X0
         self.gas2.TPX = T0, P0, X0
 
-        self.r1 = ct.Reactor(self.gas1)
+        self.r1 = ct.IdealGasReactor(self.gas1)
         self.r2 = ct.ConstPressureReactor(self.gas2)
 
         self.r1.volume = 0.2
@@ -710,11 +716,11 @@ class TestWallKinetics(utilities.CanteraTest):
         self.solid = ct.Solution('diamond.xml', 'diamond')
         self.interface = ct.Interface('diamond.xml', 'diamond_100',
                                       (self.gas, self.solid))
-        self.r1 = ct.Reactor(self.gas)
+        self.r1 = ct.IdealGasReactor(self.gas)
         self.r1.volume = 0.01
         self.net.add_reactor(self.r1)
 
-        self.r2 = ct.Reactor(self.gas)
+        self.r2 = ct.IdealGasReactor(self.gas)
         self.r2.volume = 0.01
         self.net.add_reactor(self.r2)
 
@@ -809,7 +815,7 @@ class TestReactorSensitivities(utilities.CanteraTest):
         net = ct.ReactorNet()
         gas = ct.Solution('gri30.xml')
         gas.TPX = 1300, 20*101325, 'CO:1.0, H2:0.1, CH4:0.1, H2O:0.5'
-        r1 = ct.Reactor(gas)
+        r1 = ct.IdealGasReactor(gas)
         net.add_reactor(r1)
 
         self.assertEqual(net.n_sensitivity_params, 0)
@@ -831,14 +837,14 @@ class TestReactorSensitivities(utilities.CanteraTest):
         solid = ct.Solution('diamond.xml', 'diamond')
         interface = ct.Interface('diamond.xml', 'diamond_100',
                                  (gas1, solid))
-        r1 = ct.Reactor(gas1)
+        r1 = ct.IdealGasReactor(gas1)
         net.add_reactor(r1)
         net.atol_sensitivity = 1e-10
         net.rtol_sensitivity = 1e-8
 
         gas2 = ct.Solution('h2o2.xml')
         gas2.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
-        r2 = ct.Reactor(gas2)
+        r2 = ct.IdealGasReactor(gas2)
         net.add_reactor(r2)
 
         w = ct.Wall(r1, r2)
@@ -881,7 +887,7 @@ class TestReactorSensitivities(utilities.CanteraTest):
             net = ct.ReactorNet()
             gas.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
 
-            r = ct.Reactor(gas)
+            r = ct.IdealGasReactor(gas)
             net.add_reactor(r)
             return r, net
 
@@ -920,11 +926,11 @@ class TestReactorSensitivities(utilities.CanteraTest):
             net = ct.ReactorNet()
             gas1 = ct.Solution('h2o2.xml')
             gas1.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
-            rA = ct.Reactor(gas1)
+            rA = ct.IdealGasReactor(gas1)
 
             gas2 = ct.Solution('h2o2.xml')
             gas2.TPX = 920, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:0.5'
-            rB = ct.Reactor(gas2)
+            rB = ct.IdealGasReactor(gas2)
             if reverse:
                 net.add_reactor(rB)
                 net.add_reactor(rA)
@@ -986,8 +992,8 @@ class TestReactorSensitivities(utilities.CanteraTest):
             gas1.TPX = 1200, 1e3, 'H:0.002, H2:1, CH4:0.01, CH3:0.0002'
             gas2.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
             net = ct.ReactorNet()
-            rA = ct.Reactor(gas1)
-            rB = ct.Reactor(gas2)
+            rA = ct.IdealGasReactor(gas1)
+            rB = ct.IdealGasReactor(gas2)
 
             if order % 2 == 0:
                 wA = ct.Wall(rA, rB)
diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp
new file mode 100644
index 000000000..9d6e399f6
--- /dev/null
+++ b/src/zeroD/IdealGasReactor.cpp
@@ -0,0 +1,290 @@
+/**
+ *  @file IdealGasReactor.cpp A zero-dimensional reactor
+ */
+
+#include "cantera/zeroD/IdealGasReactor.h"
+#include "cantera/zeroD/FlowDevice.h"
+#include "cantera/zeroD/Wall.h"
+#include "cantera/kinetics/InterfaceKinetics.h"
+#include "cantera/thermo/SurfPhase.h"
+#include "cantera/zeroD/ReactorNet.h"
+
+#include <cfloat>
+
+using namespace std;
+
+namespace Cantera
+{
+
+void IdealGasReactor::setThermoMgr(ThermoPhase& thermo)
+{
+    //! @TODO: Add a method to ThermoPhase that indicates whether a given
+    //! subclass is compatible with this reactor model
+    if (thermo.eosType() != cIdealGas) {
+        throw CanteraError("IdealGasReactor::setThermoMgr",
+                           "Incompatible phase type provided");
+    }
+    Reactor::setThermoMgr(thermo);
+}
+
+void IdealGasReactor::getInitialConditions(double t0, size_t leny, double* y)
+{
+    m_init = true;
+    if (m_thermo == 0) {
+        cout << "Error: reactor is empty." << endl;
+        return;
+    }
+    m_thermo->restoreState(m_state);
+
+    // set the first component to the total mass
+    m_mass = m_thermo->density() * m_vol;
+    y[0] = m_mass;
+
+    // set the second component to the total volume
+    y[1] = m_vol;
+
+    // Set the third component to the temperature
+    y[2] = m_thermo->temperature();
+
+    // set components y+3 ... y+K+2 to the mass fractions of each species
+    m_thermo->getMassFractions(y+3);
+
+    // set the remaining components to the surface species
+    // coverages on the walls
+    size_t loc = m_nsp + 3;
+    SurfPhase* surf;
+    for (size_t m = 0; m < m_nwalls; m++) {
+        surf = m_wall[m]->surface(m_lr[m]);
+        if (surf) {
+            m_wall[m]->getCoverages(m_lr[m], y + loc);
+            loc += surf->nSpecies();
+        }
+    }
+}
+
+void IdealGasReactor::initialize(doublereal t0)
+{
+    m_thermo->restoreState(m_state);
+    m_sdot.resize(m_nsp, 0.0);
+    m_wdot.resize(m_nsp, 0.0);
+    m_uk.resize(m_nsp, 0.0);
+    m_nv = m_nsp + 3;
+    for (size_t w = 0; w < m_nwalls; w++)
+        if (m_wall[w]->surface(m_lr[w])) {
+            m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
+        }
+
+    m_enthalpy = m_thermo->enthalpy_mass();
+    m_pressure = m_thermo->pressure();
+    m_intEnergy = m_thermo->intEnergy_mass();
+
+    size_t nt = 0, maxnt = 0;
+    for (size_t m = 0; m < m_nwalls; m++) {
+        m_wall[m]->initialize();
+        if (m_wall[m]->kinetics(m_lr[m])) {
+            nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
+            if (nt > maxnt) {
+                maxnt = nt;
+            }
+            if (m_wall[m]->kinetics(m_lr[m])) {
+                if (&m_kin->thermo(0) !=
+                        &m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
+                    throw CanteraError("IdealGasReactor::initialize",
+                                       "First phase of all kinetics managers must be"
+                                       " the gas.");
+                }
+            }
+        }
+    }
+    m_work.resize(maxnt);
+    std::sort(m_pnum.begin(), m_pnum.end());
+    m_init = true;
+}
+
+void IdealGasReactor::updateState(doublereal* y)
+{
+    for (size_t i = 0; i < m_nv; i++) {
+        AssertFinite(y[i], "IdealGasReactor::updateState",
+                     "y[" + int2str(i) + "] is not finite");
+    }
+
+    // The components of y are [0] the total mass, [1] the total volume,
+    // [2] the temperature, [3...K+3] are the mass fractions of each species,
+    // and [K+3...] are the coverages of surface species on each wall.
+    m_mass = y[0];
+    m_vol = y[1];
+
+    m_thermo->setMassFractions_NoNorm(y+3);
+    m_thermo->setState_TR(y[2], m_mass / m_vol);
+
+    size_t loc = m_nsp + 3;
+    SurfPhase* surf;
+    for (size_t m = 0; m < m_nwalls; m++) {
+        surf = m_wall[m]->surface(m_lr[m]);
+        if (surf) {
+            m_wall[m]->setCoverages(m_lr[m], y+loc);
+            loc += surf->nSpecies();
+        }
+    }
+
+    // save parameters needed by other connected reactors
+    m_enthalpy = m_thermo->enthalpy_mass();
+    m_pressure = m_thermo->pressure();
+    m_intEnergy = m_thermo->intEnergy_mass();
+    m_thermo->saveState(m_state);
+}
+
+void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
+                      doublereal* ydot, doublereal* params)
+{
+    m_thermo->restoreState(m_state);
+
+    // process sensitivity parameters
+    if (params) {
+        size_t npar = m_pnum.size();
+        for (size_t n = 0; n < npar; n++) {
+            double mult = m_kin->multiplier(m_pnum[n]);
+            m_kin->setMultiplier(m_pnum[n], mult*params[n]);
+        }
+        size_t ploc = npar;
+        for (size_t 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];
+            }
+        }
+    }
+
+    m_vdot = 0.0;
+    m_Q    = 0.0;
+    double mcvdTdt = 0.0; // m * c_v * dT/dt
+    double dmdt = 0.0; // dm/dt (gas phase)
+    double* dYdt = ydot + 3;
+
+    m_thermo->getPartialMolarIntEnergies(&m_uk[0]);
+
+    // compute wall terms
+    size_t loc = m_nsp+3;
+    fill(m_sdot.begin(), m_sdot.end(), 0.0);
+    for (size_t i = 0; i < m_nwalls; i++) {
+        int lr = 1 - 2*m_lr[i];
+        double vdot = lr*m_wall[i]->vdot(time);
+        m_vdot += vdot;
+        m_Q += lr*m_wall[i]->Q(time);
+        Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
+        SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
+        if (surf && kin) {
+            double rs0 = 1.0/surf->siteDensity();
+            size_t nk = surf->nSpecies();
+            double sum = 0.0;
+            surf->setTemperature(m_state[0]);
+            m_wall[i]->syncCoverages(m_lr[i]);
+            kin->getNetProductionRates(DATA_PTR(m_work));
+            size_t ns = kin->surfacePhaseIndex();
+            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
+            for (size_t k = 1; k < nk; k++) {
+                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
+                sum -= ydot[loc + k];
+            }
+            ydot[loc] = sum;
+            loc += nk;
+
+            double wallarea = m_wall[i]->area();
+            for (size_t k = 0; k < m_nsp; k++) {
+                m_sdot[k] += m_work[k]*wallarea;
+            }
+        }
+    }
+
+    const vector_fp& mw = m_thermo->molecularWeights();
+    const doublereal* Y = m_thermo->massFractions();
+
+    if (m_chem) {
+        m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
+    }
+
+    double mdot_surf = 0.0; // net mass flux from surfaces
+    for (size_t k = 0; k < m_nsp; k++) {
+        // production in gas phase and from surfaces
+        dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
+        mdot_surf += m_sdot[k] * mw[k];
+    }
+    dmdt += mdot_surf;
+
+    // compression work and external heat transfer
+    mcvdTdt += - m_pressure * m_vdot - m_Q;
+
+    for (size_t n = 0; n < m_nsp; n++) {
+        // heat release from gas phase and surface reations
+        mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol;
+        mcvdTdt -= m_sdot[n] * m_uk[n];
+        // dilution by net surface mass flux
+        dYdt[n] -= Y[n] * mdot_surf / m_mass;
+    }
+
+    // add terms for open system
+    if (m_open) {
+        // outlets
+        for (size_t i = 0; i < m_nOutlets; i++) {
+            double mdot_out = m_outlet[i]->massFlowRate(time);
+            dmdt -= mdot_out; // mass flow out of system
+            mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work
+        }
+
+        // inlets
+        for (size_t i = 0; i < m_nInlets; i++) {
+            double mdot_in = m_inlet[i]->massFlowRate(time);
+            dmdt += mdot_in; // mass flow into system
+            mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
+            for (size_t n = 0; n < m_nsp; n++) {
+                double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+                // flow of species into system and dilution by other species
+                dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
+
+                // In combintion with h_in*mdot_in, flow work plus thermal
+                // energy carried with the species
+                mcvdTdt -= m_uk[n] / mw[n] * mdot_spec;
+            }
+        }
+    }
+
+    ydot[0] = dmdt;
+    ydot[1] = m_vdot;
+    if (m_energy) {
+        ydot[2] = mcvdTdt / (m_mass * m_thermo->cv_mass());
+    } else {
+        ydot[2] = 0;
+    }
+
+    for (size_t i = 0; i < m_nv; i++) {
+        AssertFinite(ydot[i], "IdealGasReactor::evalEqs",
+                     "ydot[" + int2str(i) + "] is not finite");
+    }
+
+    // reset sensitivity parameters
+    if (params) {
+        size_t npar = m_pnum.size();
+        for (size_t n = 0; n < npar; n++) {
+            double mult = m_kin->multiplier(m_pnum[n]);
+            m_kin->setMultiplier(m_pnum[n], mult/params[n]);
+        }
+        size_t ploc = npar;
+        for (size_t m = 0; m < m_nwalls; m++) {
+            if (m_nsens_wall[m] > 0) {
+                m_wall[m]->resetSensitivityParameters(m_lr[m]);
+                ploc += m_nsens_wall[m];
+            }
+        }
+    }
+}
+
+size_t IdealGasReactor::componentIndex(const string& nm) const
+{
+    if (nm == "T") {
+        return 2;
+    } else {
+        return Reactor::componentIndex(nm);
+    }
+}
+
+}