diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h
index 0c3dad0d9..5276c45ee 100644
--- a/include/cantera/equil/MultiPhase.h
+++ b/include/cantera/equil/MultiPhase.h
@@ -426,7 +426,7 @@ public:
 
     //! Number of phases.
     size_t nPhases() const {
-        return m_np;
+        return m_phase.size();
     }
 
     //! Return true is species \a kGlob is a species in a multicomponent
@@ -617,9 +617,6 @@ private:
      */
     std::map<std::string, size_t> m_enamemap;
 
-    //! Number of phases in the MultiPhase object
-    size_t m_np;
-
     //! Current value of the temperature (kelvin)
     doublereal m_temp;
 
diff --git a/include/cantera/equil/MultiPhaseEquil.h b/include/cantera/equil/MultiPhaseEquil.h
index 2b502a3cc..5d1d69c93 100644
--- a/include/cantera/equil/MultiPhaseEquil.h
+++ b/include/cantera/equil/MultiPhaseEquil.h
@@ -161,7 +161,7 @@ protected:
         return (m_nsp > m_nel) ? m_nsp - m_nel : 0;
     }
 
-    size_t m_nel_mix, m_nsp_mix, m_np;
+    size_t m_nel_mix, m_nsp_mix;
     size_t m_nel, m_nsp;
     size_t m_eloc;
     int m_iter;
diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h
index 45c1801e3..a598b230e 100644
--- a/include/cantera/kinetics/GasKinetics.h
+++ b/include/cantera/kinetics/GasKinetics.h
@@ -70,8 +70,6 @@ public:
     virtual void update_rates_C();
 
 protected:
-    size_t m_nfall;
-
     //! Reaction index of each falloff reaction
     std::vector<size_t> m_fallindx;
 
diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h
index f2cc1757d..2f39f2aac 100644
--- a/include/cantera/kinetics/ImplicitSurfChem.h
+++ b/include/cantera/kinetics/ImplicitSurfChem.h
@@ -227,22 +227,12 @@ protected:
 
     std::vector<size_t> m_specStartIndex;
 
-    //! Total number of surface phases.
-    /*!
-     * This is also equal to the number of InterfaceKinetics objects
-     * as there is a 1-1 correspondence between InterfaceKinetics objects
-     * and surface phases.
-     */
-    size_t m_nsurf;
-
     //! Total number of surface species in all surface phases
     /*!
      * This is the total number of unknowns in m_mode 0 problem
      */
     size_t m_nv;
 
-    size_t m_numBulkPhases;
-    std::vector<size_t> m_nspBulkPhases;
     size_t m_numTotalBulkSpecies;
     size_t m_numTotalSpecies;
 
diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h
index bf013b916..f7509d99b 100644
--- a/include/cantera/kinetics/InterfaceKinetics.h
+++ b/include/cantera/kinetics/InterfaceKinetics.h
@@ -404,12 +404,6 @@ protected:
      */
     std::vector<size_t> m_irrev;
 
-    //! Number of irreversible reactions in the mechanism
-    size_t m_nirrev;
-
-    //! Number of reversible reactions in the mechanism
-    size_t m_nrev;
-
     //! Array of concentrations for each species in the kinetics mechanism
     /*!
      * An array of generalized concentrations \f$ C_k \f$ that are defined
diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h
index 358feb94a..88166713a 100644
--- a/include/cantera/kinetics/RxnRates.h
+++ b/include/cantera/kinetics/RxnRates.h
@@ -129,12 +129,12 @@ public:
         m_mcov = 0.0;
         size_t k;
         doublereal th;
-        for (size_t n = 0; n < m_ncov; n++) {
+        for (size_t n = 0; n < m_ac.size(); n++) {
             k = m_sp[n];
             m_acov += m_ac[n] * theta[k];
             m_ecov += m_ec[n] * theta[k];
         }
-        for (size_t n = 0; n < m_nmcov; n++) {
+        for (size_t n = 0; n < m_mc.size(); n++) {
             k = m_msp[n];
             th = std::max(theta[k], Tiny);
             m_mcov += m_mc[n]*std::log(th);
@@ -177,7 +177,6 @@ protected:
     doublereal m_acov, m_ecov, m_mcov;
     std::vector<size_t> m_sp, m_msp;
     vector_fp m_ac, m_ec, m_mc;
-    size_t m_ncov, m_nmcov;
 };
 
 
diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h
index 41c087220..49aa7f322 100644
--- a/include/cantera/numerics/Func1.h
+++ b/include/cantera/numerics/Func1.h
@@ -820,9 +820,8 @@ class Poly1 : public Func1
 public:
     Poly1(size_t n, doublereal* c) :
         Func1() {
-        m_n = n+1;
         m_cpoly.resize(n+1);
-        std::copy(c, c+m_n, m_cpoly.begin());
+        std::copy(c, c+m_cpoly.size(), m_cpoly.begin());
     }
 
     Poly1(const Poly1& b) :
@@ -836,7 +835,6 @@ public:
         }
         Func1::operator=(right);
         m_cpoly = right.m_cpoly;
-        m_n = right.m_n;
         m_parent = 0;
         return *this;
     }
@@ -847,16 +845,15 @@ public:
     }
 
     virtual doublereal eval(doublereal t) const {
-        doublereal r = m_cpoly[m_n-1];
-        for (size_t n = 1; n < m_n; n++) {
+        doublereal r = m_cpoly[m_cpoly.size()-1];
+        for (size_t n = 1; n < m_cpoly.size(); n++) {
             r *= t;
-            r += m_cpoly[m_n - n - 1];
+            r += m_cpoly[m_cpoly.size() - n - 1];
         }
         return r;
     }
 
 protected:
-    size_t m_n;
     vector_fp m_cpoly;
 };
 
@@ -875,7 +872,6 @@ public:
     Fourier1(size_t n, doublereal omega, doublereal a0,
              doublereal* a, doublereal* b) :
         Func1() {
-        m_n = n;
         m_omega = omega;
         m_a0_2 = 0.5*a0;
         m_ccos.resize(n);
@@ -898,7 +894,6 @@ public:
         m_a0_2 = right.m_a0_2;
         m_ccos = right.m_ccos;
         m_csin = right.m_csin;
-        m_n = right.m_n;
         m_parent = 0;
         return *this;
     }
@@ -911,7 +906,7 @@ public:
     virtual doublereal eval(doublereal t) const {
         size_t n, nn;
         doublereal sum = m_a0_2;
-        for (n = 0; n < m_n; n++) {
+        for (n = 0; n < m_ccos.size(); n++) {
             nn = n + 1;
             sum += m_ccos[n]*std::cos(m_omega*nn*t)
                    + m_csin[n]*std::sin(m_omega*nn*t);
@@ -920,7 +915,6 @@ public:
     }
 
 protected:
-    size_t m_n;
     doublereal m_omega, m_a0_2;
     vector_fp m_ccos, m_csin;
 };
@@ -937,7 +931,6 @@ class Arrhenius1 : public Func1
 public:
     Arrhenius1(size_t n, doublereal* c) :
         Func1() {
-        m_n = n;
         m_A.resize(n);
         m_b.resize(n);
         m_E.resize(n);
@@ -959,7 +952,6 @@ public:
             return *this;
         }
         Func1::operator=(right);
-        m_n = right.m_n;
         m_A = right.m_A;
         m_b = right.m_b;
         m_E = right.m_E;
@@ -974,14 +966,13 @@ public:
 
     virtual doublereal eval(doublereal t) const {
         doublereal sum = 0.0;
-        for (size_t n = 0; n < m_n; n++) {
+        for (size_t n = 0; n < m_A.size(); n++) {
             sum += m_A[n]*std::pow(t,m_b[n])*std::exp(-m_E[n]/t);
         }
         return sum;
     }
 
 protected:
-    size_t m_n;
     vector_fp m_A, m_b, m_E;
 };
 
diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h
index e4f5f6229..80572977f 100644
--- a/include/cantera/oneD/OneDim.h
+++ b/include/cantera/oneD/OneDim.h
@@ -47,7 +47,7 @@ public:
 
     /// Number of domains.
     size_t nDomains() const {
-        return m_nd;
+        return m_dom.size();
     }
 
     /// Return a reference to domain i.
@@ -60,8 +60,8 @@ public:
     //! Check that the specified domain index is in range.
     //! Throws an exception if n is greater than nDomains()-1
     void checkDomainIndex(size_t n) const {
-        if (n >= m_nd) {
-            throw IndexError("checkDomainIndex", "domains", n, m_nd-1);
+        if (n >= m_dom.size()) {
+            throw IndexError("checkDomainIndex", "domains", n, m_dom.size()-1);
         }
     }
 
@@ -69,8 +69,8 @@ public:
     //! Throws an exception if nn is less than nDomains(). Used before calls
     //! which take an array pointer.
     void checkDomainArraySize(size_t nn) const {
-        if (m_nd > nn) {
-            throw ArraySizeError("checkDomainArraySize", nn, m_nd);
+        if (m_dom.size() > nn) {
+            throw ArraySizeError("checkDomainArraySize", nn, m_dom.size());
         }
     }
 
@@ -264,9 +264,6 @@ protected:
     doublereal m_rdt; //!< reciprocal of time step
     bool m_jac_ok; //!< if true, Jacobian is current
 
-    //! number of domains
-    size_t m_nd;
-
     size_t m_bw; //!< Jacobian bandwidth
     size_t m_size; //!< solution vector size
 
diff --git a/include/cantera/thermo/AdsorbateThermo.h b/include/cantera/thermo/AdsorbateThermo.h
index acf18a31b..31edaf83c 100644
--- a/include/cantera/thermo/AdsorbateThermo.h
+++ b/include/cantera/thermo/AdsorbateThermo.h
@@ -30,9 +30,7 @@ class Adsorbate : public SpeciesThermoInterpType
 {
 public:
     //! Empty constructor
-    Adsorbate() :
-        m_nFreqs(0) {
-    }
+    Adsorbate() {}
 
     //! Full Constructor
     /*!
@@ -43,10 +41,9 @@ public:
     Adsorbate(double tlow, double thigh, double pref, const double* coeffs)
         : SpeciesThermoInterpType(tlow, thigh, pref)
     {
-        m_nFreqs = int(coeffs[0]);
+        m_freq.resize(int(coeffs[0]));
         m_be = coeffs[1];
-        m_freq.resize(m_nFreqs);
-        std::copy(coeffs+2, coeffs + 2 + m_nFreqs, m_freq.begin());
+        std::copy(coeffs+2, coeffs + 2 + m_freq.size(), m_freq.begin());
     }
 
     virtual SpeciesThermoInterpType*
@@ -76,15 +73,14 @@ public:
         tlow = m_lowT;
         thigh = m_highT;
         pref = m_Pref;
-        coeffs[0] = static_cast<double>(m_nFreqs);
+        coeffs[0] = static_cast<double>(m_freq.size());
         coeffs[1] = m_be;
-        for (size_t i = 2; i < m_nFreqs+2; i++) {
+        for (size_t i = 2; i < m_freq.size()+2; i++) {
             coeffs[i] = m_freq[i-2];
         }
     }
 
 protected:
-    size_t m_nFreqs;
     //! array of vib frequencies
     vector_fp m_freq;
     doublereal m_be;
@@ -92,7 +88,7 @@ protected:
     doublereal _energy_RT(double T) const {
         doublereal x, hnu_kt, hnu, sum = 0.0;
         doublereal kt = T*Boltzmann;
-        for (size_t i = 0; i < m_nFreqs; i++) {
+        for (size_t i = 0; i < m_freq.size(); i++) {
             hnu = Planck * m_freq[i];
             hnu_kt = hnu/kt;
             x = exp(-hnu_kt);
@@ -104,7 +100,7 @@ protected:
     doublereal _free_energy_RT(double T) const {
         doublereal x, hnu_kt, sum = 0.0;
         doublereal kt = T*Boltzmann;
-        for (size_t i = 0; i < m_nFreqs; i++) {
+        for (size_t i = 0; i < m_freq.size(); i++) {
             hnu_kt = Planck * m_freq[i] / kt;
             x = exp(-hnu_kt);
             sum += log(1.0 - x);
diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h
index 41e354296..5022bab65 100644
--- a/include/cantera/thermo/LatticeSolidPhase.h
+++ b/include/cantera/thermo/LatticeSolidPhase.h
@@ -453,9 +453,6 @@ protected:
     //! Current value of the molar density
     doublereal m_molar_density;
 
-    //! Number of sublattice phases
-    size_t m_nlattice;
-
     //! Vector of sublattic ThermoPhase objects
     std::vector<LatticePhase*> m_lattice;
 
diff --git a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h
index f08ca16ad..e7e900fe7 100644
--- a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h
+++ b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h
@@ -98,9 +98,6 @@ public:
     virtual void modifyParameters(doublereal* coeffs);
 
 protected:
-    //! Number of temperature regions
-    size_t m_numTempRegions;
-
     //! Lower boundaries of each temperature regions
     vector_fp m_lowerTempBounds;
 
diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp
index e7256522e..f5029750d 100644
--- a/src/equil/MultiPhase.cpp
+++ b/src/equil/MultiPhase.cpp
@@ -16,7 +16,6 @@ namespace Cantera
 {
 
 MultiPhase::MultiPhase() :
-    m_np(0),
     m_temp(298.15),
     m_press(OneBar),
     m_nel(0),
@@ -29,7 +28,6 @@ MultiPhase::MultiPhase() :
 }
 
 MultiPhase::MultiPhase(const MultiPhase& right) :
-    m_np(0),
     m_temp(298.15),
     m_press(OneBar),
     m_nel(0),
@@ -54,7 +52,6 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right)
         m_spstart = right.m_spstart;
         m_enames = right.m_enames;
         m_enamemap = right.m_enamemap;
-        m_np = right.m_np;
         m_temp = right.m_temp;
         m_press = right.m_press;
         m_nel = right.m_nel;
@@ -72,7 +69,7 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right)
 void MultiPhase::addPhases(MultiPhase& mix)
 {
     size_t n;
-    for (n = 0; n < mix.m_np; n++) {
+    for (n = 0; n < mix.nPhases(); n++) {
         addPhase(mix.m_phase[n], mix.m_moles[n]);
     }
 }
@@ -102,9 +99,7 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
     m_moles.push_back(moles);
     m_temp_OK.push_back(true);
 
-    // update the number of phases and the total number of
-    // species
-    m_np = m_phase.size();
+    // update the total number of species
     m_nsp += p->nSpecies();
 
     // determine if this phase has new elements for each new element, add an
@@ -174,7 +169,7 @@ void MultiPhase::init()
         sym = m_enames[m];
         k = 0;
         // iterate over the phases
-        for (ip = 0; ip < m_np; ip++) {
+        for (ip = 0; ip < nPhases(); ip++) {
             ThermoPhase* p = m_phase[ip];
             nsp = p->nSpecies();
             mlocal = p->elementIndex(sym);
@@ -248,7 +243,7 @@ doublereal MultiPhase::elementMoles(size_t m) const
 {
     doublereal sum = 0.0, phasesum;
     size_t i, k = 0, ik, nsp;
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         phasesum = 0.0;
         nsp = m_phase[i]->nSpecies();
         for (ik = 0; ik < nsp; ik++) {
@@ -264,7 +259,7 @@ doublereal MultiPhase::charge() const
 {
     doublereal sum = 0.0;
     size_t i;
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         sum += phaseCharge(i);
     }
     return sum;
@@ -301,7 +296,7 @@ void MultiPhase::getChemPotentials(doublereal* mu) const
 {
     size_t i, loc = 0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         m_phase[i]->getChemPotentials(mu + loc);
         loc += m_phase[i]->nSpecies();
     }
@@ -313,7 +308,7 @@ void MultiPhase::getValidChemPotentials(doublereal not_mu,
     size_t i, loc = 0;
     updatePhases();
     // iterate over the phases
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
             if (!standard) {
                 m_phase[i]->getChemPotentials(mu + loc);
@@ -341,7 +336,7 @@ doublereal MultiPhase::gibbs() const
     size_t i;
     doublereal sum = 0.0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (m_moles[i] > 0.0) {
             sum += m_phase[i]->gibbs_mole() * m_moles[i];
         }
@@ -354,7 +349,7 @@ doublereal MultiPhase::enthalpy() const
     size_t i;
     doublereal sum = 0.0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (m_moles[i] > 0.0) {
             sum += m_phase[i]->enthalpy_mole() * m_moles[i];
         }
@@ -367,7 +362,7 @@ doublereal MultiPhase::IntEnergy() const
     size_t i;
     doublereal sum = 0.0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (m_moles[i] > 0.0) {
             sum += m_phase[i]->intEnergy_mole() * m_moles[i];
         }
@@ -380,7 +375,7 @@ doublereal MultiPhase::entropy() const
     size_t i;
     doublereal sum = 0.0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (m_moles[i] > 0.0) {
             sum += m_phase[i]->entropy_mole() * m_moles[i];
         }
@@ -393,7 +388,7 @@ doublereal MultiPhase::cp() const
     size_t i;
     doublereal sum = 0.0;
     updatePhases();
-    for (i = 0; i < m_np; i++) {
+    for (i = 0; i < nPhases(); i++) {
         if (m_moles[i] > 0.0) {
             sum += m_phase[i]->cp_mole() * m_moles[i];
         }
@@ -437,7 +432,7 @@ void MultiPhase::getMoles(doublereal* molNum) const
     copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
     size_t ik;
     doublereal* dtmp = molNum;
-    for (size_t ip = 0; ip < m_np; ip++) {
+    for (size_t ip = 0; ip < nPhases(); ip++) {
         doublereal phasemoles = m_moles[ip];
         ThermoPhase* p = m_phase[ip];
         size_t nsp = p->nSpecies();
@@ -455,7 +450,7 @@ void MultiPhase::setMoles(const doublereal* n)
     size_t ip, loc = 0;
     size_t ik, k = 0, nsp;
     doublereal phasemoles;
-    for (ip = 0; ip < m_np; ip++) {
+    for (ip = 0; ip < nPhases(); ip++) {
         ThermoPhase* p = m_phase[ip];
         nsp = p->nSpecies();
         phasemoles = 0.0;
@@ -523,7 +518,7 @@ void MultiPhase::calcElemAbundances() const
     for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
         m_elemAbundances[eGlobal] = 0.0;
     }
-    for (size_t ip = 0; ip < m_np; ip++) {
+    for (size_t ip = 0; ip < nPhases(); ip++) {
         ThermoPhase* p = m_phase[ip];
         size_t nspPhase = p->nSpecies();
         doublereal phasemoles = m_moles[ip];
@@ -542,7 +537,7 @@ doublereal MultiPhase::volume() const
 {
     int i;
     doublereal sum = 0;
-    for (i = 0; i < int(m_np); i++) {
+    for (i = 0; i < int(nPhases()); i++) {
         double vol = 1.0/m_phase[i]->molarDensity();
         sum += m_moles[i] * vol;
     }
@@ -858,7 +853,7 @@ std::string MultiPhase::phaseName(const size_t iph) const
 int MultiPhase::phaseIndex(const std::string& pName) const
 {
     std::string tmp;
-    for (int iph = 0; iph < (int) m_np; iph++) {
+    for (int iph = 0; iph < (int) nPhases(); iph++) {
         const ThermoPhase* tptr = m_phase[iph];
         tmp = tptr->id();
         if (tmp == pName) {
@@ -896,7 +891,7 @@ bool MultiPhase::tempOK(const size_t p) const
 void MultiPhase::uploadMoleFractionsFromPhases()
 {
     size_t ip, loc = 0;
-    for (ip = 0; ip < m_np; ip++) {
+    for (ip = 0; ip < nPhases(); ip++) {
         ThermoPhase* p = m_phase[ip];
         p->getMoleFractions(&m_moleFractions[loc]);
         loc += p->nSpecies();
@@ -907,7 +902,7 @@ void MultiPhase::uploadMoleFractionsFromPhases()
 void MultiPhase::updatePhases() const
 {
     size_t p, nsp, loc = 0;
-    for (p = 0; p < m_np; p++) {
+    for (p = 0; p < nPhases(); p++) {
         nsp = m_phase[p]->nSpecies();
         m_phase[p]->setState_TPX(m_temp, m_press, &m_moleFractions[loc]);
         loc += nsp;
diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp
index d813fcce6..f8c52efbb 100644
--- a/src/equil/MultiPhaseEquil.cpp
+++ b/src/equil/MultiPhaseEquil.cpp
@@ -17,7 +17,6 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
     // store some mixture parameters locally
     m_nel_mix = mix->nElements();
     m_nsp_mix = mix->nSpecies();
-    m_np = mix->nPhases();
     m_press = mix->pressure();
     m_temp = mix->temperature();
 
@@ -580,7 +579,7 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
 
             // sum over solution phases
             doublereal sum = 0.0, psum;
-            for (ip = 0; ip < m_np; ip++) {
+            for (ip = 0; ip < m_mix->nPhases(); ip++) {
                 ThermoPhase& p = m_mix->phase(ip);
                 if (p.nSpecies() > 1) {
                     psum = 0.0;
@@ -691,7 +690,6 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
     size_t nSpecies;
     double vol = 0.0;
     string sName;
-    size_t nphase = m_np;
     FILE* FP = fopen(reportFile.c_str(), "w");
     if (!FP) {
         throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file");
@@ -708,7 +706,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
     vector_fp molalities;
 
     vol = 0.0;
-    for (size_t iphase = 0; iphase < nphase; iphase++) {
+    for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
         istart = m_mix->speciesIndex(0, iphase);
         ThermoPhase& tref = m_mix->phase(iphase);
         nSpecies = tref.nSpecies();
@@ -730,7 +728,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
     fprintf(FP,"Pressure     = %11.5g Pascal\n", pres);
     fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
 
-    for (size_t iphase = 0; iphase < nphase; iphase++) {
+    for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
         istart = m_mix->speciesIndex(0, iphase);
         ThermoPhase& tref = m_mix->phase(iphase);
         ThermoPhase* tp = &tref;
diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp
index bd8b9769d..5286745ae 100644
--- a/src/kinetics/GasKinetics.cpp
+++ b/src/kinetics/GasKinetics.cpp
@@ -12,7 +12,6 @@ namespace Cantera
 {
 GasKinetics::GasKinetics(thermo_t* thermo) :
     BulkKinetics(thermo),
-    m_nfall(0),
     m_logp_ref(0.0),
     m_logc_ref(0.0),
     m_logStandConc(0.0),
@@ -139,7 +138,7 @@ void GasKinetics::processFalloffReactions()
     // use m_ropr for temporary storage of reduced pressure
     vector_fp& pr = m_ropr;
 
-    for (size_t i = 0; i < m_nfall; i++) {
+    for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
         pr[i] = concm_falloff_values[i] * m_rfn_low[i] / (m_rfn_high[i] + SmallNumber);
         AssertFinite(pr[i], "GasKinetics::processFalloffReactions",
                      "pr[{}] is not finite.", i);
@@ -147,7 +146,7 @@ void GasKinetics::processFalloffReactions()
 
     m_falloffn.pr_to_falloff(pr.data(), falloff_work.data());
 
-    for (size_t i = 0; i < m_nfall; i++) {
+    for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
         if (reactionType(m_fallindx[i]) == FALLOFF_RXN) {
             pr[i] *= m_rfn_high[i];
         } else { // CHEMACT_RXN
@@ -155,7 +154,7 @@ void GasKinetics::processFalloffReactions()
         }
     }
 
-    scatter_copy(pr.begin(), pr.begin() + m_nfall,
+    scatter_copy(pr.begin(), pr.begin() + m_falloff_low_rates.nReactions(),
                  m_ropf.begin(), m_fallindx.begin());
 }
 
@@ -175,7 +174,7 @@ void GasKinetics::updateROP()
         m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
     }
 
-    if (m_nfall) {
+    if (m_falloff_high_rates.nReactions()) {
         processFalloffReactions();
     }
 
@@ -223,7 +222,7 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd)
         m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
     }
 
-    if (m_nfall) {
+    if (m_falloff_high_rates.nReactions()) {
         processFalloffReactions();
     }
 
@@ -271,14 +270,15 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r)
 {
     // install high and low rate coeff calculators and extend the high and low
     // rate coeff value vectors
-    m_falloff_high_rates.install(m_nfall, r.high_rate);
+    size_t nfall = m_falloff_high_rates.nReactions();
+    m_falloff_high_rates.install(nfall, r.high_rate);
     m_rfn_high.push_back(0.0);
-    m_falloff_low_rates.install(m_nfall, r.low_rate);
+    m_falloff_low_rates.install(nfall, r.low_rate);
     m_rfn_low.push_back(0.0);
 
     // add this reaction number to the list of falloff reactions
     m_fallindx.push_back(nReactions()-1);
-    m_rfallindx[nReactions()-1] = m_nfall;
+    m_rfallindx[nReactions()-1] = nfall;
 
     // install the enhanced third-body concentration calculator
     map<size_t, double> efficiencies;
@@ -292,14 +292,11 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r)
                 "' while adding reaction '" + r.equation() + "'");
         }
     }
-    m_falloff_concm.install(m_nfall, efficiencies,
+    m_falloff_concm.install(nfall, efficiencies,
                             r.third_body.default_efficiency);
 
     // install the falloff function calculator for this reaction
-    m_falloffn.install(m_nfall, r.reaction_type, r.falloff);
-
-    // increment the falloff reaction counter
-    ++m_nfall;
+    m_falloffn.install(nfall, r.reaction_type, r.falloff);
 }
 
 void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r)
diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp
index 267ba4b26..3559d0e42 100644
--- a/src/kinetics/ImplicitSurfChem.cpp
+++ b/src/kinetics/ImplicitSurfChem.cpp
@@ -16,9 +16,7 @@ namespace Cantera
 {
 
 ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
-    m_nsurf(0),
     m_nv(0),
-    m_numBulkPhases(0),
     m_numTotalBulkSpecies(0),
     m_numTotalSpecies(0),
     m_atol(1.e-14),
@@ -30,12 +28,11 @@ ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
     m_commonTempPressForPhases(true),
     m_ioFlag(0)
 {
-    m_nsurf = k.size();
     size_t ns, nsp;
     size_t nt, ntmax = 0;
     size_t kinSpIndex = 0;
     // Loop over the number of surface kinetics objects
-    for (size_t n = 0; n < m_nsurf; n++) {
+    for (size_t n = 0; n < k.size(); n++) {
         InterfaceKinetics* kinPtr = k[n];
         m_vecKinPtrs.push_back(kinPtr);
         ns = k[n]->surfacePhaseIndex();
@@ -60,9 +57,7 @@ ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
                 ThermoPhase* thPtr = & kinPtr->thermo(ip);
                 if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
                     m_bulkPhases.push_back(thPtr);
-                    m_numBulkPhases++;
                     nsp = thPtr->nSpecies();
-                    m_nspBulkPhases.push_back(nsp);
                     m_numTotalBulkSpecies += nsp;
                     imatch = m_bulkPhases.size() - 1;
                 }
@@ -110,7 +105,7 @@ void ImplicitSurfChem::getInitialConditions(doublereal t0, size_t lenc,
 void ImplicitSurfChem::getState(doublereal* c)
 {
     size_t loc = 0;
-    for (size_t n = 0; n < m_nsurf; n++) {
+    for (size_t n = 0; n < m_surf.size(); n++) {
         m_surf[n]->getCoverages(c + loc);
         loc += m_nsp[n];
     }
@@ -139,7 +134,7 @@ void ImplicitSurfChem::integrate0(doublereal t0, doublereal t1)
 void ImplicitSurfChem::updateState(doublereal* c)
 {
     size_t loc = 0;
-    for (size_t n = 0; n < m_nsurf; n++) {
+    for (size_t n = 0; n < m_surf.size(); n++) {
         m_surf[n]->setCoverages(c + loc);
         loc += m_nsp[n];
     }
@@ -151,7 +146,7 @@ void ImplicitSurfChem::eval(doublereal time, doublereal* y,
     updateState(y); // synchronize the surface state(s) with y
     doublereal rs0, sum;
     size_t loc = 0, kstart;
-    for (size_t n = 0; n < m_nsurf; n++) {
+    for (size_t n = 0; n < m_surf.size(); n++) {
         rs0 = 1.0/m_surf[n]->siteDensity();
         m_vecKinPtrs[n]->getNetProductionRates(m_work.data());
         kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
@@ -244,13 +239,13 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
 void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
 {
     size_t kstart;
-    for (size_t ip = 0; ip < m_nsurf; ip++) {
+    for (size_t ip = 0; ip < m_surf.size(); ip++) {
         ThermoPhase* TP_ptr = m_surf[ip];
         kstart = m_specStartIndex[ip];
         TP_ptr->getConcentrations(vecConcSpecies + kstart);
     }
     kstart = m_nv;
-    for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
+    for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
         ThermoPhase* TP_ptr = m_bulkPhases[ip];
         TP_ptr->getConcentrations(vecConcSpecies + kstart);
         kstart += TP_ptr->nSpecies();
@@ -260,13 +255,13 @@ void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
 void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
 {
     size_t kstart;
-    for (size_t ip = 0; ip < m_nsurf; ip++) {
+    for (size_t ip = 0; ip < m_surf.size(); ip++) {
         ThermoPhase* TP_ptr = m_surf[ip];
         kstart = m_specStartIndex[ip];
         TP_ptr->setConcentrations(vecConcSpecies + kstart);
     }
     kstart = m_nv;
-    for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
+    for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
         ThermoPhase* TP_ptr = m_bulkPhases[ip];
         TP_ptr->setConcentrations(vecConcSpecies + kstart);
         kstart += TP_ptr->nSpecies();
@@ -275,11 +270,11 @@ void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
 
 void ImplicitSurfChem::setCommonState_TP(doublereal TKelvin, doublereal PresPa)
 {
-    for (size_t ip = 0; ip < m_nsurf; ip++) {
+    for (size_t ip = 0; ip < m_surf.size(); ip++) {
         ThermoPhase* TP_ptr = m_surf[ip];
         TP_ptr->setState_TP(TKelvin, PresPa);
     }
-    for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
+    for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
         ThermoPhase* TP_ptr = m_bulkPhases[ip];
         TP_ptr->setState_TP(TKelvin, PresPa);
     }
diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp
index 33759b197..82839eebb 100644
--- a/src/kinetics/InterfaceKinetics.cpp
+++ b/src/kinetics/InterfaceKinetics.cpp
@@ -19,8 +19,6 @@ namespace Cantera
 
 InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
     m_redo_rates(false),
-    m_nirrev(0),
-    m_nrev(0),
     m_surf(0),
     m_integrator(0),
     m_logp0(0.0),
@@ -65,8 +63,6 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
     m_rates = right.m_rates;
     m_redo_rates = right.m_redo_rates;
     m_irrev = right.m_irrev;
-    m_nirrev = right.m_nirrev;
-    m_nrev = right.m_nrev;
     m_conc = right.m_conc;
     m_actConc = right.m_actConc;
     m_mu0 = right.m_mu0;
@@ -198,7 +194,7 @@ void InterfaceKinetics::updateKc()
 {
     fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
 
-    if (m_nrev > 0) {
+    if (m_revindex.size() > 0) {
         /*
          * Get the vector of standard state electrochemical potentials for
          * species in the Interfacial kinetics object and store it in m_mu0[]
@@ -210,7 +206,7 @@ void InterfaceKinetics::updateKc()
         // compute Delta mu^0 for all reversible reactions
         getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
 
-        for (size_t i = 0; i < m_nrev; i++) {
+        for (size_t i = 0; i < m_revindex.size(); i++) {
             size_t irxn = m_revindex[i];
             if (irxn == npos || irxn >= nReactions()) {
                 throw CanteraError("InterfaceKinetics", "illegal value: irxn = {}", irxn);
@@ -218,7 +214,7 @@ void InterfaceKinetics::updateKc()
             // WARNING this may overflow HKM
             m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
         }
-        for (size_t i = 0; i != m_nirrev; ++i) {
+        for (size_t i = 0; i != m_irrev.size(); ++i) {
             m_rkcn[ m_irrev[i] ] = 0.0;
         }
     }
@@ -250,7 +246,7 @@ void InterfaceKinetics::checkPartialEquil()
 
     vector_fp dmu(nTotalSpecies(), 0.0);
     vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
-    if (m_nrev > 0) {
+    if (m_revindex.size() > 0) {
         cout << "T = " << thermo(0).temperature() << " " << thermo(0).RT() << endl;
         size_t nsp, ik=0;
         doublereal delta;
@@ -267,7 +263,7 @@ void InterfaceKinetics::checkPartialEquil()
         // compute Delta mu^ for all reversible reactions
         getRevReactionDelta(dmu.data(), rmu.data());
         updateROP();
-        for (size_t i = 0; i < m_nrev; i++) {
+        for (size_t i = 0; i < m_revindex.size(); i++) {
             size_t irxn = m_revindex[i];
             writelog("Reaction {} {}\n",
                      reactionString(irxn), rmu[irxn]/thermo(0).RT());
@@ -695,10 +691,8 @@ bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
 
     if (r.reversible) {
         m_revindex.push_back(i);
-        m_nrev++;
     } else {
         m_irrev.push_back(i);
-        m_nirrev++;
     }
 
     m_rxnPhaseIsReactant.emplace_back(nPhases(), false);
diff --git a/src/kinetics/RxnRates.cpp b/src/kinetics/RxnRates.cpp
index ec7f947db..5126b6db5 100644
--- a/src/kinetics/RxnRates.cpp
+++ b/src/kinetics/RxnRates.cpp
@@ -32,8 +32,6 @@ SurfaceArrhenius::SurfaceArrhenius()
     , m_acov(0.0)
     , m_ecov(0.0)
     , m_mcov(0.0)
-    , m_ncov(0)
-    , m_nmcov(0)
 {
 }
 
@@ -44,22 +42,18 @@ SurfaceArrhenius::SurfaceArrhenius(double A, double b, double Ta)
     , m_acov(0.0)
     , m_ecov(0.0)
     , m_mcov(0.0)
-    , m_ncov(0)
-    , m_nmcov(0)
 {
 }
 
 void SurfaceArrhenius::addCoverageDependence(size_t k, doublereal a,
                                doublereal m, doublereal e)
 {
-    m_ncov++;
     m_sp.push_back(k);
     m_ac.push_back(a);
     m_ec.push_back(e);
     if (m != 0.0) {
         m_msp.push_back(k);
         m_mc.push_back(m);
-        m_nmcov++;
     }
 }
 
diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp
index 8a011cfec..afbe58517 100644
--- a/src/oneD/OneDim.cpp
+++ b/src/oneD/OneDim.cpp
@@ -15,7 +15,7 @@ namespace Cantera
 OneDim::OneDim()
     : m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
       m_rdt(0.0), m_jac_ok(false),
-      m_nd(0), m_bw(0), m_size(0),
+      m_bw(0), m_size(0),
       m_init(false), m_pts(0), m_solve_time(0.0),
       m_ss_jac_age(10), m_ts_jac_age(20),
       m_interrupt(0), m_nevals(0), m_evaltime(0.0)
@@ -26,7 +26,7 @@ OneDim::OneDim()
 OneDim::OneDim(vector<Domain1D*> domains) :
     m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
     m_rdt(0.0), m_jac_ok(false),
-    m_nd(0), m_bw(0), m_size(0),
+    m_bw(0), m_size(0),
     m_init(false), m_solve_time(0.0),
     m_ss_jac_age(10), m_ts_jac_age(20),
     m_interrupt(0), m_nevals(0), m_evaltime(0.0)
@@ -46,7 +46,7 @@ OneDim::~OneDim()
 
 size_t OneDim::domainIndex(const std::string& name)
 {
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < m_dom.size(); n++) {
         if (domain(n).id() == name) {
             return n;
         }
@@ -72,8 +72,7 @@ void OneDim::addDomain(Domain1D* d)
 
     // add it also to the global domain list, and set its container and position
     m_dom.push_back(d);
-    d->setContainer(this, m_nd);
-    m_nd++;
+    d->setContainer(this, m_dom.size()-1);
     resize();
 }
 
@@ -140,7 +139,7 @@ void OneDim::resize()
     // save the statistics for the last grid
     saveStats();
     m_pts = 0;
-    for (size_t i = 0; i < m_nd; i++) {
+    for (size_t i = 0; i < nDomains(); i++) {
         Domain1D* d = m_dom[i];
 
         size_t np = d->nPoints();
@@ -181,7 +180,7 @@ void OneDim::resize()
     m_jac.reset(new MultiJac(*this));
     m_jac_ok = false;
 
-    for (size_t i = 0; i < m_nd; i++) {
+    for (size_t i = 0; i < nDomains(); i++) {
         m_dom[i]->setJac(m_jac.get());
     }
 }
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index 6c6cba4a2..804834f2e 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -22,7 +22,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
     // domain-specific initialization of the solution vector.
     m_x.resize(size(), 0.0);
     m_xnew.resize(size(), 0.0);
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         domain(n)._getInitialSoln(&m_x[start(n)]);
     }
 
@@ -33,7 +33,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
 
 void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vector_fp& vals)
 {
-    for (size_t dom=0; dom<m_nd; dom++) {
+    for (size_t dom=0; dom<nDomains(); dom++) {
         Domain1D& d = domain(dom);
         size_t ncomp = d.nComponents();
         for (size_t comp=0; comp<ncomp; comp++) {
@@ -116,13 +116,13 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
     }
 
     vector<XML_Node*> xd = f->getChildren("domain");
-    if (xd.size() != m_nd) {
+    if (xd.size() != nDomains()) {
         throw CanteraError("Sim1D::restore", "Solution does not contain the "
             " correct number of domains. Found {} expected {}.\n",
-            xd.size(), m_nd);
+            xd.size(), nDomains());
     }
     size_t sz = 0;
-    for (size_t m = 0; m < m_nd; m++) {
+    for (size_t m = 0; m < nDomains(); m++) {
         if (loglevel > 0 && xd[m]->attrib("id") != domain(m).id()) {
             writelog("Warning: domain names do not match: '" +
                      (*xd[m])["id"] + + "' and '" + domain(m).id() + "'\n");
@@ -131,7 +131,7 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
     }
     m_x.resize(sz);
     m_xnew.resize(sz);
-    for (size_t m = 0; m < m_nd; m++) {
+    for (size_t m = 0; m < nDomains(); m++) {
         domain(m).restore(*xd[m], &m_x[domain(m).loc()], loglevel);
     }
     resize();
@@ -149,7 +149,7 @@ void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
 
 void Sim1D::showSolution(ostream& s)
 {
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         if (domain(n).domainType() != cEmptyType) {
             domain(n).showSolution_s(s, &m_x[start(n)]);
         }
@@ -158,7 +158,7 @@ void Sim1D::showSolution(ostream& s)
 
 void Sim1D::showSolution()
 {
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         if (domain(n).domainType() != cEmptyType) {
             writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
                      +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
@@ -169,14 +169,14 @@ void Sim1D::showSolution()
 
 void Sim1D::getInitialSoln()
 {
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         domain(n)._getInitialSoln(&m_x[start(n)]);
     }
 }
 
 void Sim1D::finalize()
 {
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         domain(n)._finalize(&m_x[start(n)]);
     }
 }
@@ -319,7 +319,7 @@ int Sim1D::refine(int loglevel)
     doublereal xmid, zmid;
     std::vector<size_t> dsize;
 
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         Domain1D& d = domain(n);
         Refiner& r = d.refiner();
 
@@ -377,7 +377,7 @@ int Sim1D::refine(int loglevel)
     // Now update each domain with the new grid.
 
     size_t gridstart = 0, gridsize;
-    for (size_t n = 0; n < m_nd; n++) {
+    for (size_t n = 0; n < nDomains(); n++) {
         Domain1D& d = domain(n);
         gridsize = dsize[n];
         d.setupGrid(gridsize, &znew[gridstart]);
@@ -406,7 +406,7 @@ int Sim1D::setFixedTemperature(doublereal t)
     size_t m1 = 0;
     std::vector<size_t> dsize;
 
-    for (n = 0; n < m_nd; n++) {
+    for (n = 0; n < nDomains(); n++) {
         bool addnewpt=false;
         Domain1D& d = domain(n);
         size_t comp = d.nComponents();
@@ -469,7 +469,7 @@ int Sim1D::setFixedTemperature(doublereal t)
     // been constructed, but the domains themselves have not yet been modified.
     // Now update each domain with the new grid.
     size_t gridstart = 0, gridsize;
-    for (n = 0; n < m_nd; n++) {
+    for (n = 0; n < nDomains(); n++) {
         Domain1D& d = domain(n);
         gridsize = dsize[n];
         d.setupGrid(gridsize, &znew[gridstart]);
@@ -493,7 +493,7 @@ void Sim1D::setRefineCriteria(int dom, doublereal ratio,
         Refiner& r = domain(dom).refiner();
         r.setCriteria(ratio, slope, curve, prune);
     } else {
-        for (size_t n = 0; n < m_nd; n++) {
+        for (size_t n = 0; n < nDomains(); n++) {
             Refiner& r = domain(n).refiner();
             r.setCriteria(ratio, slope, curve, prune);
         }
@@ -506,7 +506,7 @@ void Sim1D::setGridMin(int dom, double gridmin)
         Refiner& r = domain(dom).refiner();
         r.setGridMin(gridmin);
     } else {
-        for (size_t n = 0; n < m_nd; n++) {
+        for (size_t n = 0; n < nDomains(); n++) {
             Refiner& r = domain(n).refiner();
             r.setGridMin(gridmin);
         }
@@ -519,7 +519,7 @@ void Sim1D::setMaxGridPoints(int dom, int npoints)
         Refiner& r = domain(dom).refiner();
         r.setMaxPoints(npoints);
     } else {
-        for (size_t n = 0; n < m_nd; n++) {
+        for (size_t n = 0; n < nDomains(); n++) {
             Refiner& r = domain(n).refiner();
             r.setMaxPoints(npoints);
         }
diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp
index 1477852eb..847d755c1 100644
--- a/src/thermo/LatticeSolidPhase.cpp
+++ b/src/thermo/LatticeSolidPhase.cpp
@@ -20,15 +20,13 @@ namespace Cantera
 {
 LatticeSolidPhase::LatticeSolidPhase() :
     m_press(-1.0),
-    m_molar_density(0.0),
-    m_nlattice(0)
+    m_molar_density(0.0)
 {
 }
 
 LatticeSolidPhase::LatticeSolidPhase(const LatticeSolidPhase& right) :
     m_press(-1.0),
-    m_molar_density(0.0),
-    m_nlattice(0)
+    m_molar_density(0.0)
 {
     *this = right;
 }
@@ -40,7 +38,6 @@ LatticeSolidPhase& LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
         m_tlast = right.m_tlast;
         m_press = right.m_press;
         m_molar_density = right.m_molar_density;
-        m_nlattice = right.m_nlattice;
         deepStdVectorPointerCopy<LatticePhase>(right.m_lattice, m_lattice);
         m_x = right.m_x;
         theta_ = right.theta_;
@@ -52,7 +49,7 @@ LatticeSolidPhase& LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
 LatticeSolidPhase::~LatticeSolidPhase()
 {
     // We own the sublattices. So we have to delete the sublattices
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         delete m_lattice[n];
         m_lattice[n] = 0;
     }
@@ -66,14 +63,14 @@ ThermoPhase* LatticeSolidPhase::duplMyselfAsThermoPhase() const
 doublereal LatticeSolidPhase::minTemp(size_t k) const
 {
     if (k != npos) {
-        for (size_t n = 0; n < m_nlattice; n++) {
+        for (size_t n = 0; n < m_lattice.size(); n++) {
             if (lkstart_[n+1] < k) {
                 return m_lattice[n]->minTemp(k-lkstart_[n]);
             }
         }
     }
     doublereal mm = 1.0E300;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         double ml = m_lattice[n]->minTemp();
         mm = std::min(mm, ml);
     }
@@ -83,14 +80,14 @@ doublereal LatticeSolidPhase::minTemp(size_t k) const
 doublereal LatticeSolidPhase::maxTemp(size_t k) const
 {
     if (k != npos) {
-        for (size_t n = 0; n < m_nlattice; n++) {
+        for (size_t n = 0; n < m_lattice.size(); n++) {
             if (lkstart_[n+1] < k) {
                 return (m_lattice[n])->maxTemp(k - lkstart_[n]);
             }
         }
     }
     doublereal mm = -1.0E300;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         double ml = m_lattice[n]->maxTemp();
         mm = std::max(mm, ml);
     }
@@ -106,7 +103,7 @@ doublereal LatticeSolidPhase::enthalpy_mole() const
 {
     _updateThermo();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->enthalpy_mole();
     }
     return sum;
@@ -116,7 +113,7 @@ doublereal LatticeSolidPhase::intEnergy_mole() const
 {
     _updateThermo();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->intEnergy_mole();
     }
     return sum;
@@ -126,7 +123,7 @@ doublereal LatticeSolidPhase::entropy_mole() const
 {
     _updateThermo();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->entropy_mole();
     }
     return sum;
@@ -136,7 +133,7 @@ doublereal LatticeSolidPhase::gibbs_mole() const
 {
     _updateThermo();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->gibbs_mole();
     }
     return sum;
@@ -146,7 +143,7 @@ doublereal LatticeSolidPhase::cp_mole() const
 {
     _updateThermo();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->cp_mole();
     }
     return sum;
@@ -156,7 +153,7 @@ void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         m_lattice[n]->getMoleFractions(c+strt);
         strt += m_lattice[n]->nSpecies();
     }
@@ -182,7 +179,7 @@ doublereal LatticeSolidPhase::logStandardConc(size_t k) const
 void LatticeSolidPhase::setPressure(doublereal p)
 {
     m_press = p;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         m_lattice[n]->setPressure(m_press);
     }
     calcDensity();
@@ -191,7 +188,7 @@ void LatticeSolidPhase::setPressure(doublereal p)
 doublereal LatticeSolidPhase::calcDensity()
 {
     double sum = 0.0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         sum += theta_[n] * m_lattice[n]->density();
     }
     Phase::setDensity(sum);
@@ -201,13 +198,13 @@ doublereal LatticeSolidPhase::calcDensity()
 void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
 {
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nsp = m_lattice[n]->nSpecies();
         m_lattice[n]->setMoleFractions(x + strt);
         strt += nsp;
     }
     for (size_t k = 0; k < strt; k++) {
-        m_x[k] = x[k] / m_nlattice;
+        m_x[k] = x[k] / m_lattice.size();
     }
     Phase::setMoleFractions(m_x.data());
     calcDensity();
@@ -218,7 +215,7 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
     size_t strt = 0;
     // the ifdef block should be the way we calculate this.!!!!!
     Phase::getMoleFractions(x);
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nsp = m_lattice[n]->nSpecies();
         double sum = 0.0;
         for (size_t k = 0; k < nsp; k++) {
@@ -245,7 +242,7 @@ void LatticeSolidPhase::getChemPotentials(doublereal* mu) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nlsp = m_lattice[n]->nSpecies();
         m_lattice[n]->getChemPotentials(mu+strt);
         strt += nlsp;
@@ -256,7 +253,7 @@ void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nlsp = m_lattice[n]->nSpecies();
         m_lattice[n]->getPartialMolarEnthalpies(hbar + strt);
         strt += nlsp;
@@ -267,7 +264,7 @@ void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nlsp = m_lattice[n]->nSpecies();
         m_lattice[n]->getPartialMolarEntropies(sbar + strt);
         strt += nlsp;
@@ -278,7 +275,7 @@ void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nlsp = m_lattice[n]->nSpecies();
         m_lattice[n]->getPartialMolarCp(cpbar + strt);
         strt += nlsp;
@@ -289,7 +286,7 @@ void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nlsp = m_lattice[n]->nSpecies();
         m_lattice[n]->getPartialMolarVolumes(vbar + strt);
         strt += nlsp;
@@ -300,7 +297,7 @@ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
 {
     _updateThermo();
     size_t strt = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         m_lattice[n]->getStandardChemPotentials(mu0+strt);
         strt += m_lattice[n]->nSpecies();
     }
@@ -309,7 +306,7 @@ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
 void LatticeSolidPhase::getGibbs_RT_ref(doublereal* grt) const
 {
     _updateThermo();
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
     }
 }
@@ -330,7 +327,7 @@ void LatticeSolidPhase::installSlavePhases(XML_Node* phaseNode)
 
     XML_Node& la = phaseNode->child("thermo").child("LatticeArray");
     std::vector<XML_Node*> lattices = la.getChildren("phase");
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         LatticePhase* lp = m_lattice[n];
         vector_fp constArr(lp->nElements());
         const vector_fp& aws = lp->atomicWeights();
@@ -365,11 +362,11 @@ void LatticeSolidPhase::initThermo()
 {
     initLengths();
     size_t loc = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nsp = m_lattice[n]->nSpecies();
         lkstart_[n] = loc;
         for (size_t k = 0; k < nsp; k++) {
-            m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_nlattice;
+            m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_lattice.size();
             loc++;
         }
         lkstart_[n+1] = loc;
@@ -380,8 +377,8 @@ void LatticeSolidPhase::initThermo()
 
 void LatticeSolidPhase::initLengths()
 {
-    theta_.resize(m_nlattice,0);
-    lkstart_.resize(m_nlattice+1);
+    theta_.resize(m_lattice.size(), 0);
+    lkstart_.resize(m_lattice.size() + 1);
     m_x.resize(m_kk, 0.0);
     tmpV_.resize(m_kk, 0.0);
 }
@@ -392,7 +389,7 @@ void LatticeSolidPhase::_updateThermo() const
     if (m_tlast != tnow) {
         getMoleFractions(m_x.data());
         size_t strt = 0;
-        for (size_t n = 0; n < m_nlattice; n++) {
+        for (size_t n = 0; n < m_lattice.size(); n++) {
             m_lattice[n]->setTemperature(tnow);
             m_lattice[n]->setMoleFractions(&m_x[strt]);
             m_lattice[n]->setPressure(m_press);
@@ -406,7 +403,7 @@ void LatticeSolidPhase::setLatticeMoleFractionsByName(int nn, const std::string&
 {
     m_lattice[nn]->setMoleFractionsByName(x);
     size_t loc = 0;
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nsp = m_lattice[n]->nSpecies();
         double ndens = m_lattice[n]->molarDensity();
         for (size_t k = 0; k < nsp; k++) {
@@ -422,18 +419,17 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
     eosdata._require("model","LatticeSolid");
     XML_Node& la = eosdata.child("LatticeArray");
     std::vector<XML_Node*> lattices = la.getChildren("phase");
-    m_nlattice = lattices.size();
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < lattices.size(); n++) {
         m_lattice.push_back((LatticePhase*)newPhase(*lattices[n]));
     }
     std::vector<string> pnam;
     std::vector<string> pval;
     int np = getPairs(eosdata.child("LatticeStoichiometry"), pnam, pval);
-    theta_.resize(m_nlattice);
+    theta_.resize(m_lattice.size());
     for (int i = 0; i < np; i++) {
         double val = fpValueCheck(pval[i]);
         bool found = false;
-        for (size_t j = 0; j < m_nlattice; j++) {
+        for (size_t j = 0; j < m_lattice.size(); j++) {
             ThermoPhase& tp = *m_lattice[j];
             string idj = tp.id();
             if (idj == pnam[i]) {
@@ -450,7 +446,7 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
 
 void LatticeSolidPhase::modifyOneHf298SS(const size_t k, const doublereal Hf298New)
 {
-    for (size_t n = 0; n < m_nlattice; n++) {
+    for (size_t n = 0; n < m_lattice.size(); n++) {
         if (lkstart_[n+1] < k) {
             size_t kk = k-lkstart_[n];
             SpeciesThermo& l_spthermo = m_lattice[n]->speciesThermo();
diff --git a/src/thermo/Nasa9PolyMultiTempRegion.cpp b/src/thermo/Nasa9PolyMultiTempRegion.cpp
index d070d8d81..088b66412 100644
--- a/src/thermo/Nasa9PolyMultiTempRegion.cpp
+++ b/src/thermo/Nasa9PolyMultiTempRegion.cpp
@@ -20,25 +20,22 @@ using namespace std;
 namespace Cantera
 {
 Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion() :
-    m_numTempRegions(0),
     m_currRegion(0)
 {
 }
 
 Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector<Nasa9Poly1*>& regionPts) :
-    m_numTempRegions(0),
     m_currRegion(0)
 {
-    m_numTempRegions = regionPts.size();
     // From now on, we own these pointers
     for (Nasa9Poly1* region : regionPts) {
         m_regionPts.emplace_back(region);
     }
-    m_lowerTempBounds.resize(m_numTempRegions);
+    m_lowerTempBounds.resize(regionPts.size());
     m_lowT = m_regionPts[0]->minTemp();
-    m_highT = m_regionPts[m_numTempRegions-1]->maxTemp();
+    m_highT = m_regionPts[m_regionPts.size()-1]->maxTemp();
     m_Pref = m_regionPts[0]->refPressure();
-    for (size_t i = 0; i < m_numTempRegions; i++) {
+    for (size_t i = 0; i < m_regionPts.size(); i++) {
         m_lowerTempBounds[i] = m_regionPts[i]->minTemp();
         if (fabs(m_regionPts[i]->refPressure() - m_Pref) > 0.0001) {
             throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
@@ -59,12 +56,11 @@ Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector<Nasa9Poly1*>& regionPt
 
 Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(const Nasa9PolyMultiTempRegion& b) :
     SpeciesThermoInterpType(b),
-    m_numTempRegions(b.m_numTempRegions),
     m_lowerTempBounds(b.m_lowerTempBounds),
     m_currRegion(b.m_currRegion)
 {
-    m_regionPts.resize(m_numTempRegions);
-    for (size_t i = 0; i < m_numTempRegions; i++) {
+    m_regionPts.resize(b.m_regionPts.size());
+    for (size_t i = 0; i < m_regionPts.size(); i++) {
         m_regionPts[i].reset(new Nasa9Poly1(*b.m_regionPts[i]));
     }
 }
@@ -74,11 +70,10 @@ Nasa9PolyMultiTempRegion::operator=(const Nasa9PolyMultiTempRegion& b)
 {
     if (&b != this) {
         SpeciesThermoInterpType::operator=(b);
-        m_numTempRegions = b.m_numTempRegions;
         m_lowerTempBounds = b.m_lowerTempBounds;
         m_currRegion = b.m_currRegion;
-        m_regionPts.resize(m_numTempRegions);
-        for (size_t i = 0; i < m_numTempRegions; i++) {
+        m_regionPts.resize(b.m_regionPts.size());
+        for (size_t i = 0; i < m_regionPts.size(); i++) {
             m_regionPts[i].reset(new Nasa9Poly1(*b.m_regionPts[i]));
         }
     }
@@ -117,7 +112,7 @@ void Nasa9PolyMultiTempRegion::updateProperties(const doublereal* tt,
         doublereal* s_R) const
 {
     m_currRegion = 0;
-    for (size_t i = 1; i < m_numTempRegions; i++) {
+    for (size_t i = 1; i < m_regionPts.size(); i++) {
         if (tt[0] < m_lowerTempBounds[i]) {
             break;
         }
@@ -133,7 +128,7 @@ void Nasa9PolyMultiTempRegion::updatePropertiesTemp(const doublereal temp,
 {
     // Now find the region
     m_currRegion = 0;
-    for (size_t i = 1; i < m_numTempRegions; i++) {
+    for (size_t i = 1; i < m_regionPts.size(); i++) {
         if (temp < m_lowerTempBounds[i]) {
             break;
         }
@@ -154,12 +149,12 @@ void Nasa9PolyMultiTempRegion::reportParameters(size_t& n, int& type,
     thigh = m_highT;
     pref = m_Pref;
     double ctmp[12];
-    coeffs[0] = double(m_numTempRegions);
+    coeffs[0] = double(m_regionPts.size());
     int index = 1;
     size_t n_tmp = 0;
     int type_tmp = 0;
     double pref_tmp = 0.0;
-    for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
+    for (size_t iReg = 0; iReg < m_regionPts.size(); iReg++) {
         m_regionPts[iReg]->reportParameters(n_tmp, type_tmp,
                                             coeffs[index], coeffs[index+1],
                                             pref_tmp, ctmp);
@@ -173,7 +168,7 @@ void Nasa9PolyMultiTempRegion::reportParameters(size_t& n, int& type,
 void Nasa9PolyMultiTempRegion::modifyParameters(doublereal* coeffs)
 {
     int index = 3;
-    for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
+    for (size_t iReg = 0; iReg < m_regionPts.size(); iReg++) {
         m_regionPts[iReg]->modifyParameters(coeffs + index);
         index += 11;
     }