diff --git a/Mechanisms/Aromatic_KrNara/mechanism.H b/Mechanisms/Aromatic_KrNara/mechanism.H
index 9c1dc4296..6459411a9 100644
--- a/Mechanisms/Aromatic_KrNara/mechanism.H
+++ b/Mechanisms/Aromatic_KrNara/mechanism.H
@@ -175,6 +175,7 @@
 #define NUM_GAS_ELEMENTS 5     // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 158    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 1049 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0  // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/Aromatic_KrNara/mechanism.cpp b/Mechanisms/Aromatic_KrNara/mechanism.cpp
index bd7cf76bc..30f3bb432 100644
--- a/Mechanisms/Aromatic_KrNara/mechanism.cpp
+++ b/Mechanisms/Aromatic_KrNara/mechanism.cpp
@@ -877,7 +877,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -887,12 +887,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[158]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 158; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -900,7 +900,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 1049; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -915,7 +915,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[158];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -947,7 +947,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 158; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -1627,7 +1627,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "H";
   ename[2] = "O";
@@ -1639,7 +1639,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(158);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "O2";
   kname[2] = "O";
diff --git a/Mechanisms/BurkeDryer/mechanism.H b/Mechanisms/BurkeDryer/mechanism.H
index e3963ebdf..c38425bb3 100644
--- a/Mechanisms/BurkeDryer/mechanism.H
+++ b/Mechanisms/BurkeDryer/mechanism.H
@@ -28,9 +28,10 @@
 #define CO_ID 11
 #define CO2_ID 12
 
-#define NUM_GAS_ELEMENTS 6   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 13   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 27 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 6    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 13    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 27  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/BurkeDryer/mechanism.cpp b/Mechanisms/BurkeDryer/mechanism.cpp
index 8e1cbc018..bb8f24668 100644
--- a/Mechanisms/BurkeDryer/mechanism.cpp
+++ b/Mechanisms/BurkeDryer/mechanism.cpp
@@ -48,7 +48,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -58,12 +58,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[13]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 13; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -71,7 +71,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 27; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -86,7 +86,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[13];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -119,7 +119,7 @@ CKNCF(int* ncf)
 {
   int kd = 6;
   // Zero ncf
-  for (int id = 0; id < kd * 13; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -173,7 +173,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(6);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "N";
@@ -186,7 +186,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(13);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "H2";
   kname[2] = "O";
diff --git a/Mechanisms/C1-C2-NO/mechanism.H b/Mechanisms/C1-C2-NO/mechanism.H
index fb72e3874..96ac4b1e8 100644
--- a/Mechanisms/C1-C2-NO/mechanism.H
+++ b/Mechanisms/C1-C2-NO/mechanism.H
@@ -61,6 +61,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 44    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 269 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/C1-C2-NO/mechanism.cpp b/Mechanisms/C1-C2-NO/mechanism.cpp
index 3488e9a65..4bbcef37e 100644
--- a/Mechanisms/C1-C2-NO/mechanism.cpp
+++ b/Mechanisms/C1-C2-NO/mechanism.cpp
@@ -199,7 +199,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -209,12 +209,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[44]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 44; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -222,7 +222,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 269; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -237,7 +237,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[44];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -269,7 +269,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 44; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -460,7 +460,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -472,7 +472,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(44);
+  kname.resize(NUM_SPECIES);
   kname[0] = "AR";
   kname[1] = "N2";
   kname[2] = "H";
diff --git a/Mechanisms/C1-C2-NO_qss/mechanism.H b/Mechanisms/C1-C2-NO_qss/mechanism.H
index 552156769..b1c5924f2 100644
--- a/Mechanisms/C1-C2-NO_qss/mechanism.H
+++ b/Mechanisms/C1-C2-NO_qss/mechanism.H
@@ -53,9 +53,10 @@
 #define NCO_ID 37
 #define CH3O2_ID 38
 
-#define NUM_GAS_ELEMENTS 5      // Elements in the homogeneous phase
-#define NUM_QSSA_GAS_SPECIES 39 // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 269   // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 39    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 269 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 5 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
@@ -64,7 +65,7 @@
 #define NUM_SURFACE_REACTIONS 0 // Reactions in the heterogeneous phase
 
 #define NUM_ELEMENTS (NUM_GAS_ELEMENTS + NUM_SURFACE_ELEMENTS)
-#define NUM_SPECIES (NUM_QSSA_GAS_SPECIES + NUM_SURFACE_SPECIES)
+#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)
 #define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)
 
 #define NUM_IONS 0
diff --git a/Mechanisms/C1-C2-NO_qss/mechanism.cpp b/Mechanisms/C1-C2-NO_qss/mechanism.cpp
index cd03f71c7..471f5a97a 100644
--- a/Mechanisms/C1-C2-NO_qss/mechanism.cpp
+++ b/Mechanisms/C1-C2-NO_qss/mechanism.cpp
@@ -199,7 +199,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -209,12 +209,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[39]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 39; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -222,7 +222,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 269; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -237,9 +237,9 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[39];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
-  amrex::Real g_RT_qss[5];
+  amrex::Real g_RT_qss[NUM_GAS_SPECIES];
   gibbs_qss(g_RT_qss, T);
 
   amrex::Real sc_qss[5];
@@ -276,7 +276,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 39; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -448,7 +448,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "Ar";
   ename[1] = "N";
   ename[2] = "H";
@@ -460,7 +460,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(39);
+  kname.resize(NUM_SPECIES);
   kname[0] = "AR";
   kname[1] = "N2";
   kname[2] = "H";
diff --git a/Mechanisms/CH4_lean/mechanism.H b/Mechanisms/CH4_lean/mechanism.H
index 68d4557f5..346e474fd 100644
--- a/Mechanisms/CH4_lean/mechanism.H
+++ b/Mechanisms/CH4_lean/mechanism.H
@@ -31,9 +31,10 @@
 #define CH2OH_ID 15
 #define N2_ID 16
 
-#define NUM_GAS_ELEMENTS 5   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 17   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 73 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 17    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 73  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/CH4_lean/mechanism.cpp b/Mechanisms/CH4_lean/mechanism.cpp
index 06e75113f..80ae4cd86 100644
--- a/Mechanisms/CH4_lean/mechanism.cpp
+++ b/Mechanisms/CH4_lean/mechanism.cpp
@@ -79,7 +79,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -89,12 +89,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[17]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 17; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -102,7 +102,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 73; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -117,7 +117,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[17];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -149,7 +149,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 17; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -224,7 +224,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -236,7 +236,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(17);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/CH4_lean_qss/mechanism.H b/Mechanisms/CH4_lean_qss/mechanism.H
index d30d57655..280edc7ad 100644
--- a/Mechanisms/CH4_lean_qss/mechanism.H
+++ b/Mechanisms/CH4_lean_qss/mechanism.H
@@ -26,9 +26,10 @@
 #define CH2O_ID 11
 #define N2_ID 12
 
-#define NUM_GAS_ELEMENTS 4      // Elements in the homogeneous phase
-#define NUM_QSSA_GAS_SPECIES 13 // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 73    // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 13    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 73  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 4 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
@@ -37,7 +38,7 @@
 #define NUM_SURFACE_REACTIONS 0 // Reactions in the heterogeneous phase
 
 #define NUM_ELEMENTS (NUM_GAS_ELEMENTS + NUM_SURFACE_ELEMENTS)
-#define NUM_SPECIES (NUM_QSSA_GAS_SPECIES + NUM_SURFACE_SPECIES)
+#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)
 #define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)
 
 #define NUM_IONS 0
diff --git a/Mechanisms/CH4_lean_qss/mechanism.cpp b/Mechanisms/CH4_lean_qss/mechanism.cpp
index e223a62ff..11ccb17e1 100644
--- a/Mechanisms/CH4_lean_qss/mechanism.cpp
+++ b/Mechanisms/CH4_lean_qss/mechanism.cpp
@@ -79,7 +79,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -89,12 +89,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[13]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 13; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -102,7 +102,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 73; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -117,9 +117,9 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[13];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
-  amrex::Real g_RT_qss[4];
+  amrex::Real g_RT_qss[NUM_GAS_SPECIES];
   gibbs_qss(g_RT_qss, T);
 
   amrex::Real sc_qss[4];
@@ -155,7 +155,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 13; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -212,7 +212,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "C";
@@ -223,7 +223,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(13);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/Davis/mechanism.H b/Mechanisms/Davis/mechanism.H
index 9c6525933..a3c6412ba 100644
--- a/Mechanisms/Davis/mechanism.H
+++ b/Mechanisms/Davis/mechanism.H
@@ -29,9 +29,10 @@
 #define H2O2_ID 12
 #define CO2_ID 13
 
-#define NUM_GAS_ELEMENTS 6   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 14   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 38 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 6    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 14    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 38  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/Davis/mechanism.cpp b/Mechanisms/Davis/mechanism.cpp
index de9d095cd..ac30cb6db 100644
--- a/Mechanisms/Davis/mechanism.cpp
+++ b/Mechanisms/Davis/mechanism.cpp
@@ -54,7 +54,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -64,12 +64,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[14]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 14; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -77,7 +77,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 38; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -92,7 +92,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[14];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -125,7 +125,7 @@ CKNCF(int* ncf)
 {
   int kd = 6;
   // Zero ncf
-  for (int id = 0; id < kd * 14; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -184,7 +184,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(6);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -197,7 +197,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(14);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "AR";
diff --git a/Mechanisms/FFCM1_Red/mechanism.H b/Mechanisms/FFCM1_Red/mechanism.H
index afbd4933f..02d88a066 100644
--- a/Mechanisms/FFCM1_Red/mechanism.H
+++ b/Mechanisms/FFCM1_Red/mechanism.H
@@ -38,6 +38,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 21    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 124 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/FFCM1_Red/mechanism.cpp b/Mechanisms/FFCM1_Red/mechanism.cpp
index dfbaaa59d..d5c602fb8 100644
--- a/Mechanisms/FFCM1_Red/mechanism.cpp
+++ b/Mechanisms/FFCM1_Red/mechanism.cpp
@@ -110,7 +110,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -120,12 +120,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[21]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 21; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -133,7 +133,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 124; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -148,7 +148,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[21];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -180,7 +180,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 21; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -272,7 +272,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -284,7 +284,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(21);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "H2";
   kname[2] = "H";
diff --git a/Mechanisms/H2-CO-CO2-3spec/mechanism.H b/Mechanisms/H2-CO-CO2-3spec/mechanism.H
index d9c7d11c1..8983388ad 100644
--- a/Mechanisms/H2-CO-CO2-3spec/mechanism.H
+++ b/Mechanisms/H2-CO-CO2-3spec/mechanism.H
@@ -15,9 +15,10 @@
 #define CO_ID 1
 #define CO2_ID 2
 
-#define NUM_GAS_ELEMENTS 3  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 3   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 3    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 3     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/H2-CO-CO2-3spec/mechanism.cpp b/Mechanisms/H2-CO-CO2-3spec/mechanism.cpp
index 320ae8330..2d46c0e40 100644
--- a/Mechanisms/H2-CO-CO2-3spec/mechanism.cpp
+++ b/Mechanisms/H2-CO-CO2-3spec/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[3]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 3; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -76,7 +76,7 @@ CKNCF(int* ncf)
 {
   int kd = 3;
   // Zero ncf
-  for (int id = 0; id < kd * 3; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -96,7 +96,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(3);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "C";
   ename[2] = "O";
@@ -106,7 +106,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(3);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "CO";
   kname[2] = "CO2";
diff --git a/Mechanisms/IonizedAir/mechanism.H b/Mechanisms/IonizedAir/mechanism.H
index 5c6d9364f..84a57930d 100644
--- a/Mechanisms/IonizedAir/mechanism.H
+++ b/Mechanisms/IonizedAir/mechanism.H
@@ -17,9 +17,10 @@
 #define On_ID 3
 #define O2n_ID 4
 
-#define NUM_GAS_ELEMENTS 3  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 5   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 3    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 5     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/IonizedAir/mechanism.cpp b/Mechanisms/IonizedAir/mechanism.cpp
index 73d0ad8dd..8150d47cc 100644
--- a/Mechanisms/IonizedAir/mechanism.cpp
+++ b/Mechanisms/IonizedAir/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[5]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 5; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -76,7 +76,7 @@ CKNCF(int* ncf)
 {
   int kd = 3;
   // Zero ncf
-  for (int id = 0; id < kd * 5; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -102,7 +102,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(3);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "N";
   ename[2] = "E";
@@ -112,7 +112,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(5);
+  kname.resize(NUM_SPECIES);
   kname[0] = "O2";
   kname[1] = "N2";
   kname[2] = "E";
diff --git a/Mechanisms/JL4/mechanism.H b/Mechanisms/JL4/mechanism.H
index 9595b7185..5370cedd5 100644
--- a/Mechanisms/JL4/mechanism.H
+++ b/Mechanisms/JL4/mechanism.H
@@ -20,9 +20,10 @@
 #define CO2_ID 5
 #define H2_ID 6
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 7   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 4 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 7     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 4   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/JL4/mechanism.cpp b/Mechanisms/JL4/mechanism.cpp
index 81183c58e..cb32edfde 100644
--- a/Mechanisms/JL4/mechanism.cpp
+++ b/Mechanisms/JL4/mechanism.cpp
@@ -36,7 +36,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -46,12 +46,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[7]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 7; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -59,7 +59,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 4; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -74,7 +74,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[7];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -105,7 +105,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 7; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -139,7 +139,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "O";
   ename[2] = "H";
@@ -150,7 +150,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(7);
+  kname.resize(NUM_SPECIES);
   kname[0] = "CH4";
   kname[1] = "O2";
   kname[2] = "H2O";
diff --git a/Mechanisms/Kolla/mechanism.H b/Mechanisms/Kolla/mechanism.H
index 6063885f1..459a5b21f 100644
--- a/Mechanisms/Kolla/mechanism.H
+++ b/Mechanisms/Kolla/mechanism.H
@@ -25,9 +25,10 @@
 #define HCO_ID 10
 #define N2_ID 11
 
-#define NUM_GAS_ELEMENTS 4   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 12   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 29 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 12    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 29  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/Kolla/mechanism.cpp b/Mechanisms/Kolla/mechanism.cpp
index df49868d8..87794f928 100644
--- a/Mechanisms/Kolla/mechanism.cpp
+++ b/Mechanisms/Kolla/mechanism.cpp
@@ -49,7 +49,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -59,12 +59,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[12]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 12; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -72,7 +72,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 29; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -87,7 +87,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[12];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -118,7 +118,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 12; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -171,7 +171,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -182,7 +182,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(12);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "O2";
   kname[2] = "O";
diff --git a/Mechanisms/LiDryer/mechanism.H b/Mechanisms/LiDryer/mechanism.H
index 526639c41..0a257f0d7 100644
--- a/Mechanisms/LiDryer/mechanism.H
+++ b/Mechanisms/LiDryer/mechanism.H
@@ -21,9 +21,10 @@
 #define H2O2_ID 7
 #define N2_ID 8
 
-#define NUM_GAS_ELEMENTS 3   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 9    // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 21 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 3    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 9     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 21  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/LiDryer/mechanism.cpp b/Mechanisms/LiDryer/mechanism.cpp
index ea90c0929..d0051d1a8 100644
--- a/Mechanisms/LiDryer/mechanism.cpp
+++ b/Mechanisms/LiDryer/mechanism.cpp
@@ -45,7 +45,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -55,12 +55,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[9]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 9; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -68,7 +68,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 21; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -83,7 +83,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[9];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -113,7 +113,7 @@ CKNCF(int* ncf)
 {
   int kd = 3;
   // Zero ncf
-  for (int id = 0; id < kd * 9; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -153,7 +153,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(3);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "N";
@@ -163,7 +163,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(9);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "O2";
   kname[2] = "H2O";
diff --git a/Mechanisms/LuDME/mechanism.H b/Mechanisms/LuDME/mechanism.H
index b33eaacfb..2f372088d 100644
--- a/Mechanisms/LuDME/mechanism.H
+++ b/Mechanisms/LuDME/mechanism.H
@@ -55,6 +55,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 39    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 175 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/LuDME/mechanism.cpp b/Mechanisms/LuDME/mechanism.cpp
index 503b62886..cb5eae7aa 100644
--- a/Mechanisms/LuDME/mechanism.cpp
+++ b/Mechanisms/LuDME/mechanism.cpp
@@ -141,7 +141,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -151,12 +151,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[39]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 39; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -164,7 +164,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 175; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -179,7 +179,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[39];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -210,7 +210,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 39; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -389,7 +389,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -400,7 +400,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(39);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "H2";
   kname[2] = "CH2";
diff --git a/Mechanisms/LuEthylene/mechanism.H b/Mechanisms/LuEthylene/mechanism.H
index 911aa44d4..4edb624b4 100644
--- a/Mechanisms/LuEthylene/mechanism.H
+++ b/Mechanisms/LuEthylene/mechanism.H
@@ -49,6 +49,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 32    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 206 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/LuEthylene/mechanism.cpp b/Mechanisms/LuEthylene/mechanism.cpp
index 7b4ac4b7e..4f0c1138f 100644
--- a/Mechanisms/LuEthylene/mechanism.cpp
+++ b/Mechanisms/LuEthylene/mechanism.cpp
@@ -161,7 +161,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -171,12 +171,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[32]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 32; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -184,7 +184,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 206; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -199,7 +199,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[32];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -231,7 +231,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 32; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -370,7 +370,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -382,7 +382,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(32);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/LuEthylene_qss/mechanism.H b/Mechanisms/LuEthylene_qss/mechanism.H
index b79d5638f..6b87db283 100644
--- a/Mechanisms/LuEthylene_qss/mechanism.H
+++ b/Mechanisms/LuEthylene_qss/mechanism.H
@@ -35,9 +35,10 @@
 #define C3H6_ID 20
 #define N2_ID 21
 
-#define NUM_GAS_ELEMENTS 4      // Elements in the homogeneous phase
-#define NUM_QSSA_GAS_SPECIES 22 // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 206   // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4     // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 22     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 206  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 10 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
@@ -46,7 +47,7 @@
 #define NUM_SURFACE_REACTIONS 0 // Reactions in the heterogeneous phase
 
 #define NUM_ELEMENTS (NUM_GAS_ELEMENTS + NUM_SURFACE_ELEMENTS)
-#define NUM_SPECIES (NUM_QSSA_GAS_SPECIES + NUM_SURFACE_SPECIES)
+#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)
 #define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)
 
 #define NUM_IONS 0
diff --git a/Mechanisms/LuEthylene_qss/mechanism.cpp b/Mechanisms/LuEthylene_qss/mechanism.cpp
index 3d048e92c..00e934cf2 100644
--- a/Mechanisms/LuEthylene_qss/mechanism.cpp
+++ b/Mechanisms/LuEthylene_qss/mechanism.cpp
@@ -161,7 +161,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -171,12 +171,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[22]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 22; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -184,7 +184,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 206; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -199,9 +199,9 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[22];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
-  amrex::Real g_RT_qss[10];
+  amrex::Real g_RT_qss[NUM_GAS_SPECIES];
   gibbs_qss(g_RT_qss, T);
 
   amrex::Real sc_qss[10];
@@ -237,7 +237,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 22; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -333,7 +333,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "C";
@@ -344,7 +344,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(22);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/NUIGalway/mechanism.H b/Mechanisms/NUIGalway/mechanism.H
index eecb5f66b..7c47f5d4f 100644
--- a/Mechanisms/NUIGalway/mechanism.H
+++ b/Mechanisms/NUIGalway/mechanism.H
@@ -136,6 +136,7 @@
 #define NUM_GAS_ELEMENTS 6     // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 118    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 1246 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0  // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/NUIGalway/mechanism.cpp b/Mechanisms/NUIGalway/mechanism.cpp
index 7ec74d5c1..edfdce3a3 100644
--- a/Mechanisms/NUIGalway/mechanism.cpp
+++ b/Mechanisms/NUIGalway/mechanism.cpp
@@ -908,7 +908,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -918,12 +918,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[118]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 118; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -931,7 +931,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 1246; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -946,7 +946,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[118];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -979,7 +979,7 @@ CKNCF(int* ncf)
 {
   int kd = 6;
   // Zero ncf
-  for (int id = 0; id < kd * 118; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -1534,7 +1534,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(6);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "N";
@@ -1547,7 +1547,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(118);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "H2";
   kname[2] = "C";
diff --git a/Mechanisms/SootReaction/mechanism.H b/Mechanisms/SootReaction/mechanism.H
index 91b0b1490..d0e03dde4 100644
--- a/Mechanisms/SootReaction/mechanism.H
+++ b/Mechanisms/SootReaction/mechanism.H
@@ -63,6 +63,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 47    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 224 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/SootReaction/mechanism.cpp b/Mechanisms/SootReaction/mechanism.cpp
index a9d1b643b..e5d845e47 100644
--- a/Mechanisms/SootReaction/mechanism.cpp
+++ b/Mechanisms/SootReaction/mechanism.cpp
@@ -200,7 +200,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -210,12 +210,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[47]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 47; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -223,7 +223,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 224; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -238,7 +238,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[47];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -269,7 +269,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 47; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -465,7 +465,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "C";
   ename[2] = "H";
@@ -476,7 +476,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(47);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "S-CH2";
   kname[2] = "T-CH2";
diff --git a/Mechanisms/air/mechanism.H b/Mechanisms/air/mechanism.H
index 29be62ea2..a2f4b4de0 100644
--- a/Mechanisms/air/mechanism.H
+++ b/Mechanisms/air/mechanism.H
@@ -13,9 +13,10 @@
 #define O2_ID 0
 #define N2_ID 1
 
-#define NUM_GAS_ELEMENTS 2  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 2   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 2    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 2     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/air/mechanism.cpp b/Mechanisms/air/mechanism.cpp
index 8b2aa66d9..278a0da23 100644
--- a/Mechanisms/air/mechanism.cpp
+++ b/Mechanisms/air/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[2]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 2; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -75,7 +75,7 @@ CKNCF(int* ncf)
 {
   int kd = 2;
   // Zero ncf
-  for (int id = 0; id < kd * 2; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -90,7 +90,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(2);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "N";
 }
@@ -99,7 +99,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(2);
+  kname.resize(NUM_SPECIES);
   kname[0] = "O2";
   kname[1] = "N2";
 }
diff --git a/Mechanisms/alzeta/mechanism.H b/Mechanisms/alzeta/mechanism.H
index ba3098163..2f9b972a4 100644
--- a/Mechanisms/alzeta/mechanism.H
+++ b/Mechanisms/alzeta/mechanism.H
@@ -90,6 +90,7 @@
 #define NUM_GAS_ELEMENTS 6    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 72    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 518 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/alzeta/mechanism.cpp b/Mechanisms/alzeta/mechanism.cpp
index 9ab11d840..37d2f19c0 100644
--- a/Mechanisms/alzeta/mechanism.cpp
+++ b/Mechanisms/alzeta/mechanism.cpp
@@ -357,7 +357,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -367,12 +367,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[72]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 72; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -380,7 +380,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 518; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -395,7 +395,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[72];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -428,7 +428,7 @@ CKNCF(int* ncf)
 {
   int kd = 6;
   // Zero ncf
-  for (int id = 0; id < kd * 72; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -747,7 +747,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(6);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -760,7 +760,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(72);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/chem-CH4-2step/mechanism.H b/Mechanisms/chem-CH4-2step/mechanism.H
index 4d875c157..f65c7d902 100644
--- a/Mechanisms/chem-CH4-2step/mechanism.H
+++ b/Mechanisms/chem-CH4-2step/mechanism.H
@@ -19,9 +19,10 @@
 #define CO2_ID 4
 #define N2_ID 5
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 6   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 3 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 6     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 3   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/chem-CH4-2step/mechanism.cpp b/Mechanisms/chem-CH4-2step/mechanism.cpp
index fda87060c..b507fc32a 100644
--- a/Mechanisms/chem-CH4-2step/mechanism.cpp
+++ b/Mechanisms/chem-CH4-2step/mechanism.cpp
@@ -35,7 +35,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -45,12 +45,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[6]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 6; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -58,7 +58,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 3; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -73,7 +73,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[6];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -104,7 +104,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 6; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -135,7 +135,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -146,7 +146,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(6);
+  kname.resize(NUM_SPECIES);
   kname[0] = "O2";
   kname[1] = "H2O";
   kname[2] = "CH4";
diff --git a/Mechanisms/chem-H/mechanism.H b/Mechanisms/chem-H/mechanism.H
index 889329d33..050980593 100644
--- a/Mechanisms/chem-H/mechanism.H
+++ b/Mechanisms/chem-H/mechanism.H
@@ -21,9 +21,10 @@
 #define H2O2_ID 7
 #define N2_ID 8
 
-#define NUM_GAS_ELEMENTS 3   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 9    // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 27 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 3    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 9     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 27  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/chem-H/mechanism.cpp b/Mechanisms/chem-H/mechanism.cpp
index d5202ab62..558e587fb 100644
--- a/Mechanisms/chem-H/mechanism.cpp
+++ b/Mechanisms/chem-H/mechanism.cpp
@@ -48,7 +48,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -58,12 +58,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[9]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 9; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -71,7 +71,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 27; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -86,7 +86,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[9];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -116,7 +116,7 @@ CKNCF(int* ncf)
 {
   int kd = 3;
   // Zero ncf
-  for (int id = 0; id < kd * 9; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -156,7 +156,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(3);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "N";
@@ -166,7 +166,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(9);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/decane_3sp/mechanism.H b/Mechanisms/decane_3sp/mechanism.H
index e745dedf7..7d2be8861 100644
--- a/Mechanisms/decane_3sp/mechanism.H
+++ b/Mechanisms/decane_3sp/mechanism.H
@@ -16,9 +16,10 @@
 #define O2_ID 1
 #define N2_ID 2
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 3   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 3     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/decane_3sp/mechanism.cpp b/Mechanisms/decane_3sp/mechanism.cpp
index 37bd152ff..04578e68e 100644
--- a/Mechanisms/decane_3sp/mechanism.cpp
+++ b/Mechanisms/decane_3sp/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[3]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 3; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -77,7 +77,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 3; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -96,7 +96,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -107,7 +107,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(3);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC10H22";
   kname[1] = "O2";
   kname[2] = "N2";
diff --git a/Mechanisms/dodecane_lu/mechanism.H b/Mechanisms/dodecane_lu/mechanism.H
index b1b85d5b4..15f5f38e6 100644
--- a/Mechanisms/dodecane_lu/mechanism.H
+++ b/Mechanisms/dodecane_lu/mechanism.H
@@ -69,6 +69,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 53    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 268 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/dodecane_lu/mechanism.cpp b/Mechanisms/dodecane_lu/mechanism.cpp
index 6ead27ac1..fc16e9996 100644
--- a/Mechanisms/dodecane_lu/mechanism.cpp
+++ b/Mechanisms/dodecane_lu/mechanism.cpp
@@ -201,7 +201,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -211,12 +211,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[53]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 53; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -224,7 +224,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 268; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -239,7 +239,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[53];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -270,7 +270,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 53; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -495,7 +495,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -506,7 +506,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(53);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC12H26";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/dodecane_lu_qss/mechanism.H b/Mechanisms/dodecane_lu_qss/mechanism.H
index 5360d90c1..677c92d7e 100644
--- a/Mechanisms/dodecane_lu_qss/mechanism.H
+++ b/Mechanisms/dodecane_lu_qss/mechanism.H
@@ -48,9 +48,10 @@
 #define OC12H23OOH_ID 33
 #define N2_ID 34
 
-#define NUM_GAS_ELEMENTS 4      // Elements in the homogeneous phase
-#define NUM_QSSA_GAS_SPECIES 35 // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 268   // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4     // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 35     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 268  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 18 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
@@ -59,7 +60,7 @@
 #define NUM_SURFACE_REACTIONS 0 // Reactions in the heterogeneous phase
 
 #define NUM_ELEMENTS (NUM_GAS_ELEMENTS + NUM_SURFACE_ELEMENTS)
-#define NUM_SPECIES (NUM_QSSA_GAS_SPECIES + NUM_SURFACE_SPECIES)
+#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)
 #define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)
 
 #define NUM_IONS 0
diff --git a/Mechanisms/dodecane_lu_qss/mechanism.cpp b/Mechanisms/dodecane_lu_qss/mechanism.cpp
index d25f0d0a4..f968c74bb 100644
--- a/Mechanisms/dodecane_lu_qss/mechanism.cpp
+++ b/Mechanisms/dodecane_lu_qss/mechanism.cpp
@@ -201,7 +201,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -211,12 +211,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[35]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 35; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -224,7 +224,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 268; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -239,9 +239,9 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[35];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
-  amrex::Real g_RT_qss[18];
+  amrex::Real g_RT_qss[NUM_GAS_SPECIES];
   gibbs_qss(g_RT_qss, T);
 
   amrex::Real sc_qss[18];
@@ -277,7 +277,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 35; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -426,7 +426,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -437,7 +437,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(35);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC12H26";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/dodecane_wang/mechanism.H b/Mechanisms/dodecane_wang/mechanism.H
index ad1f17ff2..3b4d16f4f 100644
--- a/Mechanisms/dodecane_wang/mechanism.H
+++ b/Mechanisms/dodecane_wang/mechanism.H
@@ -72,6 +72,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 56    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 289 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/dodecane_wang/mechanism.cpp b/Mechanisms/dodecane_wang/mechanism.cpp
index ca582fb41..3a58236a6 100644
--- a/Mechanisms/dodecane_wang/mechanism.cpp
+++ b/Mechanisms/dodecane_wang/mechanism.cpp
@@ -241,7 +241,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -251,12 +251,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[56]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 56; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -264,7 +264,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 289; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -279,7 +279,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[56];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -310,7 +310,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 56; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -549,7 +549,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -560,7 +560,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(56);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "O";
   kname[2] = "OH";
diff --git a/Mechanisms/dodmethair_4sp/mechanism.H b/Mechanisms/dodmethair_4sp/mechanism.H
index 206450251..88ce13316 100644
--- a/Mechanisms/dodmethair_4sp/mechanism.H
+++ b/Mechanisms/dodmethair_4sp/mechanism.H
@@ -17,9 +17,10 @@
 #define O2_ID 2
 #define N2_ID 3
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 4   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 4     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/dodmethair_4sp/mechanism.cpp b/Mechanisms/dodmethair_4sp/mechanism.cpp
index 87cdbc66b..e85528f62 100644
--- a/Mechanisms/dodmethair_4sp/mechanism.cpp
+++ b/Mechanisms/dodmethair_4sp/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[4]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 4; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -77,7 +77,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 4; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -100,7 +100,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -111,7 +111,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(4);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC12H26";
   kname[1] = "CH4";
   kname[2] = "O2";
diff --git a/Mechanisms/drm19/mechanism.H b/Mechanisms/drm19/mechanism.H
index 2d5b9f66c..273bad1f5 100644
--- a/Mechanisms/drm19/mechanism.H
+++ b/Mechanisms/drm19/mechanism.H
@@ -35,9 +35,10 @@
 #define N2_ID 19
 #define AR_ID 20
 
-#define NUM_GAS_ELEMENTS 5   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 21   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 84 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 21    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 84  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/drm19/mechanism.cpp b/Mechanisms/drm19/mechanism.cpp
index b99a99cc2..c5525e196 100644
--- a/Mechanisms/drm19/mechanism.cpp
+++ b/Mechanisms/drm19/mechanism.cpp
@@ -85,7 +85,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -95,12 +95,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[21]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 21; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -108,7 +108,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 84; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -123,7 +123,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[21];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -155,7 +155,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 21; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -245,7 +245,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -257,7 +257,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(21);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/ethylene_af/mechanism.H b/Mechanisms/ethylene_af/mechanism.H
index 27b57ac00..4d05b4a6c 100644
--- a/Mechanisms/ethylene_af/mechanism.H
+++ b/Mechanisms/ethylene_af/mechanism.H
@@ -45,6 +45,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 29    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 172 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/ethylene_af/mechanism.cpp b/Mechanisms/ethylene_af/mechanism.cpp
index 693e932d5..2af77b6e5 100644
--- a/Mechanisms/ethylene_af/mechanism.cpp
+++ b/Mechanisms/ethylene_af/mechanism.cpp
@@ -139,7 +139,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -149,12 +149,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[29]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 29; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -162,7 +162,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 172; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -177,7 +177,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[29];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -208,7 +208,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 29; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -334,7 +334,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "H";
   ename[2] = "O";
@@ -345,7 +345,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(29);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "H";
   kname[2] = "H2";
diff --git a/Mechanisms/grimech12/mechanism.H b/Mechanisms/grimech12/mechanism.H
index 91633180e..8249e66c7 100644
--- a/Mechanisms/grimech12/mechanism.H
+++ b/Mechanisms/grimech12/mechanism.H
@@ -49,6 +49,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 32    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 177 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/grimech12/mechanism.cpp b/Mechanisms/grimech12/mechanism.cpp
index 9057b1a98..d1577d0a8 100644
--- a/Mechanisms/grimech12/mechanism.cpp
+++ b/Mechanisms/grimech12/mechanism.cpp
@@ -144,7 +144,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -154,12 +154,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[32]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 32; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -167,7 +167,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 177; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -182,7 +182,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[32];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -214,7 +214,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 32; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -352,7 +352,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -364,7 +364,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(32);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/grimech30-noArN/mechanism.H b/Mechanisms/grimech30-noArN/mechanism.H
index 47c95d4a9..50b56451a 100644
--- a/Mechanisms/grimech30-noArN/mechanism.H
+++ b/Mechanisms/grimech30-noArN/mechanism.H
@@ -51,6 +51,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 35    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 217 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/grimech30-noArN/mechanism.cpp b/Mechanisms/grimech30-noArN/mechanism.cpp
index 753adbb61..d81ebaa88 100644
--- a/Mechanisms/grimech30-noArN/mechanism.cpp
+++ b/Mechanisms/grimech30-noArN/mechanism.cpp
@@ -168,7 +168,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -178,12 +178,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[35]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 35; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -191,7 +191,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 217; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -206,7 +206,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[35];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -237,7 +237,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 35; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -390,7 +390,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -401,7 +401,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(35);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/grimech30/mechanism.H b/Mechanisms/grimech30/mechanism.H
index 4e035c297..4035d0954 100644
--- a/Mechanisms/grimech30/mechanism.H
+++ b/Mechanisms/grimech30/mechanism.H
@@ -70,6 +70,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 53    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 325 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/grimech30/mechanism.cpp b/Mechanisms/grimech30/mechanism.cpp
index d54dc8493..5ec07826d 100644
--- a/Mechanisms/grimech30/mechanism.cpp
+++ b/Mechanisms/grimech30/mechanism.cpp
@@ -237,7 +237,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -247,12 +247,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[53]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 53; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -260,7 +260,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 325; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -275,7 +275,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[53];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -307,7 +307,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 53; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -541,7 +541,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -553,7 +553,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(53);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H2";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/heptane_3sp/mechanism.H b/Mechanisms/heptane_3sp/mechanism.H
index aec578456..e8a2a9eb2 100644
--- a/Mechanisms/heptane_3sp/mechanism.H
+++ b/Mechanisms/heptane_3sp/mechanism.H
@@ -16,9 +16,10 @@
 #define O2_ID 1
 #define N2_ID 2
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 3   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 3     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/heptane_3sp/mechanism.cpp b/Mechanisms/heptane_3sp/mechanism.cpp
index 135bf36ac..23d27a900 100644
--- a/Mechanisms/heptane_3sp/mechanism.cpp
+++ b/Mechanisms/heptane_3sp/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[3]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 3; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -77,7 +77,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 3; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -96,7 +96,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -107,7 +107,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(3);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC7H16";
   kname[1] = "O2";
   kname[2] = "N2";
diff --git a/Mechanisms/heptane_fc/mechanism.H b/Mechanisms/heptane_fc/mechanism.H
index 172076875..372cb1279 100644
--- a/Mechanisms/heptane_fc/mechanism.H
+++ b/Mechanisms/heptane_fc/mechanism.H
@@ -68,6 +68,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 52    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 218 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/heptane_fc/mechanism.cpp b/Mechanisms/heptane_fc/mechanism.cpp
index 94ce4a4da..7f8318eee 100644
--- a/Mechanisms/heptane_fc/mechanism.cpp
+++ b/Mechanisms/heptane_fc/mechanism.cpp
@@ -168,7 +168,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -178,12 +178,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[52]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 52; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -191,7 +191,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 218; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -206,7 +206,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[52];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -237,7 +237,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 52; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -464,7 +464,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "O";
   ename[2] = "H";
@@ -475,7 +475,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(52);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "O";
   kname[2] = "H2";
diff --git a/Mechanisms/heptane_lu_88sk/mechanism.H b/Mechanisms/heptane_lu_88sk/mechanism.H
index 456f521dd..edb2c01cb 100644
--- a/Mechanisms/heptane_lu_88sk/mechanism.H
+++ b/Mechanisms/heptane_lu_88sk/mechanism.H
@@ -104,6 +104,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 88    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 729 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/heptane_lu_88sk/mechanism.cpp b/Mechanisms/heptane_lu_88sk/mechanism.cpp
index 856e946e4..fa8af8d84 100644
--- a/Mechanisms/heptane_lu_88sk/mechanism.cpp
+++ b/Mechanisms/heptane_lu_88sk/mechanism.cpp
@@ -488,7 +488,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -498,12 +498,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[88]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 88; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -511,7 +511,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 729; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -526,7 +526,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[88];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -557,7 +557,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 88; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -963,7 +963,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "C";
   ename[2] = "O";
@@ -974,7 +974,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(88);
+  kname.resize(NUM_SPECIES);
   kname[0] = "h";
   kname[1] = "h2";
   kname[2] = "o";
diff --git a/Mechanisms/heptane_lu_qss/mechanism.H b/Mechanisms/heptane_lu_qss/mechanism.H
index b51508687..c3c14743b 100644
--- a/Mechanisms/heptane_lu_qss/mechanism.H
+++ b/Mechanisms/heptane_lu_qss/mechanism.H
@@ -68,9 +68,10 @@
 #define c7h14o24_ID 53
 #define n2_ID 54
 
-#define NUM_GAS_ELEMENTS 4      // Elements in the homogeneous phase
-#define NUM_QSSA_GAS_SPECIES 55 // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 724   // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4     // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 55     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 724  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 33 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
@@ -79,7 +80,7 @@
 #define NUM_SURFACE_REACTIONS 0 // Reactions in the heterogeneous phase
 
 #define NUM_ELEMENTS (NUM_GAS_ELEMENTS + NUM_SURFACE_ELEMENTS)
-#define NUM_SPECIES (NUM_QSSA_GAS_SPECIES + NUM_SURFACE_SPECIES)
+#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)
 #define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)
 
 #define NUM_IONS 0
diff --git a/Mechanisms/heptane_lu_qss/mechanism.cpp b/Mechanisms/heptane_lu_qss/mechanism.cpp
index d896fce50..056fb18da 100644
--- a/Mechanisms/heptane_lu_qss/mechanism.cpp
+++ b/Mechanisms/heptane_lu_qss/mechanism.cpp
@@ -484,7 +484,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -494,12 +494,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[55]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 55; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -507,7 +507,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 724; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -522,9 +522,9 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[55];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
-  amrex::Real g_RT_qss[33];
+  amrex::Real g_RT_qss[NUM_GAS_SPECIES];
   gibbs_qss(g_RT_qss, T);
 
   amrex::Real sc_qss[33];
@@ -560,7 +560,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 55; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -810,7 +810,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "C";
@@ -821,7 +821,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(55);
+  kname.resize(NUM_SPECIES);
   kname[0] = "h";
   kname[1] = "h2";
   kname[2] = "o";
diff --git a/Mechanisms/isooctane_lu/mechanism.H b/Mechanisms/isooctane_lu/mechanism.H
index 1f358d9c4..5d9b1e1ce 100644
--- a/Mechanisms/isooctane_lu/mechanism.H
+++ b/Mechanisms/isooctane_lu/mechanism.H
@@ -159,6 +159,7 @@
 #define NUM_GAS_ELEMENTS 4     // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 143    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 1178 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0  // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/isooctane_lu/mechanism.cpp b/Mechanisms/isooctane_lu/mechanism.cpp
index a02aa1a91..ad49c8c64 100644
--- a/Mechanisms/isooctane_lu/mechanism.cpp
+++ b/Mechanisms/isooctane_lu/mechanism.cpp
@@ -861,7 +861,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -871,12 +871,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[143]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 143; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -884,7 +884,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 1178; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -899,7 +899,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[143];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -930,7 +930,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 143; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -1589,7 +1589,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "C";
   ename[1] = "H";
   ename[2] = "O";
@@ -1600,7 +1600,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(143);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "H2";
   kname[2] = "O";
diff --git a/Mechanisms/methaneIons_diRenzo/mechanism.H b/Mechanisms/methaneIons_diRenzo/mechanism.H
index f41e7fed5..ecb0a0ecc 100644
--- a/Mechanisms/methaneIons_diRenzo/mechanism.H
+++ b/Mechanisms/methaneIons_diRenzo/mechanism.H
@@ -43,6 +43,7 @@
 #define NUM_GAS_ELEMENTS 5    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 26    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 134 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/methaneIons_diRenzo/mechanism.cpp b/Mechanisms/methaneIons_diRenzo/mechanism.cpp
index f2068f409..ff29efab1 100644
--- a/Mechanisms/methaneIons_diRenzo/mechanism.cpp
+++ b/Mechanisms/methaneIons_diRenzo/mechanism.cpp
@@ -117,7 +117,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -127,12 +127,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[26]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 26; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -140,7 +140,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 134; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -155,7 +155,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[26];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -187,7 +187,7 @@ CKNCF(int* ncf)
 {
   int kd = 5;
   // Zero ncf
-  for (int id = 0; id < kd * 26; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -301,7 +301,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(5);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "H";
   ename[2] = "O";
@@ -313,7 +313,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(26);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "H2";
   kname[2] = "H";
diff --git a/Mechanisms/ndodecane_35/mechanism.H b/Mechanisms/ndodecane_35/mechanism.H
index 6d664ce25..0d2fb1e46 100644
--- a/Mechanisms/ndodecane_35/mechanism.H
+++ b/Mechanisms/ndodecane_35/mechanism.H
@@ -48,9 +48,10 @@
 #define OC12H23OOH_ID 33
 #define N2_ID 34
 
-#define NUM_GAS_ELEMENTS 4  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 35  // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 35    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/ndodecane_35/mechanism.cpp b/Mechanisms/ndodecane_35/mechanism.cpp
index 906ea3184..49bcc9c56 100644
--- a/Mechanisms/ndodecane_35/mechanism.cpp
+++ b/Mechanisms/ndodecane_35/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[35]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 35; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -77,7 +77,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 35; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -226,7 +226,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "O";
   ename[1] = "H";
   ename[2] = "C";
@@ -237,7 +237,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(35);
+  kname.resize(NUM_SPECIES);
   kname[0] = "NC12H26";
   kname[1] = "H";
   kname[2] = "O";
diff --git a/Mechanisms/nitrogens/mechanism.H b/Mechanisms/nitrogens/mechanism.H
index 3ad0c49cb..d5c2ad52e 100644
--- a/Mechanisms/nitrogens/mechanism.H
+++ b/Mechanisms/nitrogens/mechanism.H
@@ -12,9 +12,10 @@
 #define N2_ID 0
 #define N2a_ID 1
 
-#define NUM_GAS_ELEMENTS 1  // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 2   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 0 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 1    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 2     // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 0   // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/nitrogens/mechanism.cpp b/Mechanisms/nitrogens/mechanism.cpp
index c1389ad7d..e044d77c1 100644
--- a/Mechanisms/nitrogens/mechanism.cpp
+++ b/Mechanisms/nitrogens/mechanism.cpp
@@ -19,7 +19,7 @@ CKINU(const int i, int& nspec, int* /*ki*/, int* /*nu*/)
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -29,12 +29,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[2]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 2; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -74,7 +74,7 @@ CKNCF(int* ncf)
 {
   int kd = 1;
   // Zero ncf
-  for (int id = 0; id < kd * 2; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -89,7 +89,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(1);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
 }
 
@@ -97,7 +97,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(2);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "N2a";
 }
diff --git a/Mechanisms/propane_fc/mechanism.H b/Mechanisms/propane_fc/mechanism.H
index d54a811e4..d1690b343 100644
--- a/Mechanisms/propane_fc/mechanism.H
+++ b/Mechanisms/propane_fc/mechanism.H
@@ -50,6 +50,7 @@
 #define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
 #define NUM_GAS_SPECIES 34    // Species in the homogeneous phase
 #define NUM_GAS_REACTIONS 167 // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/propane_fc/mechanism.cpp b/Mechanisms/propane_fc/mechanism.cpp
index 67e791b59..205e68739 100644
--- a/Mechanisms/propane_fc/mechanism.cpp
+++ b/Mechanisms/propane_fc/mechanism.cpp
@@ -137,7 +137,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -147,12 +147,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[34]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 34; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -160,7 +160,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 167; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -175,7 +175,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[34];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -206,7 +206,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 34; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -356,7 +356,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "N";
   ename[1] = "O";
   ename[2] = "H";
@@ -367,7 +367,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(34);
+  kname.resize(NUM_SPECIES);
   kname[0] = "N2";
   kname[1] = "O";
   kname[2] = "H";
diff --git a/Mechanisms/sCO2/mechanism.H b/Mechanisms/sCO2/mechanism.H
index 45956f6c6..2c55ca049 100644
--- a/Mechanisms/sCO2/mechanism.H
+++ b/Mechanisms/sCO2/mechanism.H
@@ -29,9 +29,10 @@
 #define CH2O_ID 14
 #define HCO_ID 15
 
-#define NUM_GAS_ELEMENTS 4   // Elements in the homogeneous phase
-#define NUM_GAS_SPECIES 16   // Species in the homogeneous phase
-#define NUM_GAS_REACTIONS 60 // Reactions in the homogeneous phase
+#define NUM_GAS_ELEMENTS 4    // Elements in the homogeneous phase
+#define NUM_GAS_SPECIES 16    // Species in the homogeneous phase
+#define NUM_GAS_REACTIONS 60  // Reactions in the homogeneous phase
+#define NUM_QSS_GAS_SPECIES 0 // QSS species in the homogeneous phase
 
 #define SITE_DENSITY 0.000000E+00 // mol/cm^2
 
diff --git a/Mechanisms/sCO2/mechanism.cpp b/Mechanisms/sCO2/mechanism.cpp
index d67606b08..c7430573c 100644
--- a/Mechanisms/sCO2/mechanism.cpp
+++ b/Mechanisms/sCO2/mechanism.cpp
@@ -71,7 +71,7 @@ CKINU(const int i, int& nspec, int ki[], int nu[])
   }
 }
 
-// Returns the progress rates of each reactions
+// Returns the progress rates of each reaction
 // Given P, T, and mole fractions
 void
 CKKFKR(
@@ -81,12 +81,12 @@ CKKFKR(
   amrex::Real q_f[],
   amrex::Real q_r[])
 {
-  amrex::Real c[16]; // temporary storage
+  amrex::Real c[NUM_SPECIES]; // temporary storage
   amrex::Real PORT =
-    1e6 * P / (8.31446261815324e+07 * T); // 1e6 * P/RT so c goes to SI units
+    1e6 * P / (8.31446261815324e+07 * T); // convert to SI (mol/cm^3 to mol/m^3)
 
   // Compute conversion, see Eq 10
-  for (int id = 0; id < 16; ++id) {
+  for (int id = 0; id < NUM_GAS_SPECIES; ++id) {
     c[id] = x[id] * PORT;
   }
 
@@ -94,7 +94,7 @@ CKKFKR(
   progressRateFR(q_f, q_r, c, T);
 
   // convert to chemkin units
-  for (int id = 0; id < 60; ++id) {
+  for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {
     q_f[id] *= 1.0e-6;
     q_r[id] *= 1.0e-6;
   }
@@ -109,7 +109,7 @@ progressRateFR(
   const amrex::Real invT = 1.0 / T;
   const amrex::Real logT = log(T);
   // compute the Gibbs free energy
-  amrex::Real g_RT[16];
+  amrex::Real g_RT[NUM_SPECIES];
   gibbs(g_RT, T);
 
   amrex::Real sc_qss[1];
@@ -140,7 +140,7 @@ CKNCF(int* ncf)
 {
   int kd = 4;
   // Zero ncf
-  for (int id = 0; id < kd * 16; ++id) {
+  for (int id = 0; id < kd * NUM_GAS_SPECIES; ++id) {
     ncf[id] = 0;
   }
 
@@ -215,7 +215,7 @@ CKNCF(int* ncf)
 void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
-  ename.resize(4);
+  ename.resize(NUM_ELEMENTS);
   ename[0] = "H";
   ename[1] = "O";
   ename[2] = "N";
@@ -226,7 +226,7 @@ CKSYME_STR(amrex::Vector<std::string>& ename)
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
-  kname.resize(16);
+  kname.resize(NUM_SPECIES);
   kname[0] = "H";
   kname[1] = "O2";
   kname[2] = "O";
diff --git a/Support/ceptr/ceptr/ck.py b/Support/ceptr/ceptr/ck.py
index 025de2282..230c941e7 100644
--- a/Support/ceptr/ceptr/ck.py
+++ b/Support/ceptr/ceptr/ck.py
@@ -2,6 +2,7 @@
 
 import ceptr.constants as cc
 import ceptr.thermo as cth
+import ceptr.utilities as cu
 import ceptr.writer as cw
 
 
@@ -15,31 +16,42 @@ def ckawt(fstream, mechanism):
     cw.writer(fstream, "}")
 
 
-def ckncf(fstream, mechanism, species_info):
-    """Write ckncf."""
+def ckncf(fstream, mechanism, species_info, write_sk=False):
+    """Write ckncf/skncf."""
     n_elements = mechanism.n_elements
-    n_species = species_info.n_species
+    sp_list = (
+        species_info.surface_species_list
+        if write_sk
+        else species_info.nonqssa_species_list
+    )
 
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("Returns the elemental composition "))
     cw.writer(fstream, cw.comment("of the speciesi (mdim is num of elements)"))
-    cw.writer(fstream, "void CKNCF" + cc.sym + "(int * ncf)")
+    cw.writer(
+        fstream, f"void {cu.get_function_prefix(write_sk)}KNCF" + cc.sym + "(int * ncf)"
+    )
     cw.writer(fstream, "{")
     cw.writer(fstream, f"int kd = {n_elements}; ")
     cw.writer(fstream, cw.comment("Zero ncf"))
-    cw.writer(fstream, f"for (int id = 0; id < kd * {n_species}; ++ id) {{")
+    cw.writer(
+        fstream,
+        f"for (int id = 0; id < kd * NUM_{cu.get_phase(write_sk).upper()}_SPECIES; ++"
+        " id) {",
+    )
     cw.writer(fstream, " ncf[id] = 0; ")
     cw.writer(fstream, "}")
 
+    array_idx_correction = species_info.n_gas_species if write_sk else 0
     cw.writer(fstream)
-    for sp in species_info.nonqssa_species_list:
+    for sp in sp_list:
         spec_idx = species_info.ordered_idx_map[sp]
-        species = species_info.nonqssa_species[spec_idx]
+        species = species_info.all_species[spec_idx]
         cw.writer(fstream, cw.comment(f"{species.name}"))
         for elem, coef in mechanism.species(sp).composition.items():
             cw.writer(
                 fstream,
-                f"ncf[ {species_info.ordered_idx_map[sp]} * kd +"
+                f"ncf[ {species_info.ordered_idx_map[sp] - array_idx_correction} * kd +"
                 f" {mechanism.element_index(elem)} ] = {int(coef)}; "
                 + cw.comment(f"{elem}"),
             )
@@ -47,9 +59,10 @@ def ckncf(fstream, mechanism, species_info):
     cw.writer(fstream, "}")
 
 
-def cksyme_str(fstream, mechanism, species_info):
+def cksyme_str(fstream, mechanism, interface):
     """Write cksyme."""
-    n_elements = mechanism.n_elements
+    element_names = cu.get_element_names(mechanism, interface)
+
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("Returns the vector of strings of element names"))
     cw.writer(
@@ -57,18 +70,18 @@ def cksyme_str(fstream, mechanism, species_info):
         "void CKSYME_STR" + cc.sym + "(amrex::Vector<std::string>& ename)",
     )
     cw.writer(fstream, "{")
-    cw.writer(fstream, f"ename.resize({n_elements});")
-    for elem in mechanism.element_names:
+    cw.writer(fstream, "ename.resize(NUM_ELEMENTS);")
+    for elem in element_names:
+        idx = cu.get_element_id(mechanism, interface, elem)
         cw.writer(
             fstream,
-            f'ename[{mechanism.element_index(elem)}] = "{elem}";',
+            f'ename[{idx}] = "{elem}";',
         )
     cw.writer(fstream, "}")
 
 
-def cksyms_str(fstream, mechanism, species_info):
+def cksyms_str(fstream, species_info, mech_is_heterogeneous):
     """Write cksyms."""
-    n_species = species_info.n_species
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("Returns the vector of strings of species names"))
     cw.writer(
@@ -76,12 +89,18 @@ def cksyms_str(fstream, mechanism, species_info):
         "void CKSYMS_STR" + cc.sym + "(amrex::Vector<std::string>& kname)",
     )
     cw.writer(fstream, "{")
-    cw.writer(fstream, f"kname.resize({n_species});")
+    cw.writer(fstream, "kname.resize(NUM_SPECIES);")
     for species in species_info.nonqssa_species_list:
         cw.writer(
             fstream,
             f'kname[{species_info.ordered_idx_map[species]}] = "{species}";',
         )
+    if mech_is_heterogeneous:
+        for species in species_info.surface_species_list:
+            cw.writer(
+                fstream,
+                f'kname[{species_info.ordered_idx_map[species]}] = "{species}";',
+            )
     cw.writer(fstream, "}")
 
 
@@ -2424,8 +2443,7 @@ def ckinu(fstream, mechanism, species_info, reaction_info, write_sk=False):
     """Write ckinu/skinu."""
     n_reactions = mechanism.n_reactions
     n_gas_reactions = reaction_info.n_reactions
-    phase = "surface" if write_sk else "gas"
-    function_prefix = "S" if write_sk else "C"
+    phase = cu.get_phase(is_heterogeneous=write_sk)
     function_args = (
         "int* /*ki*/, int* /*nu*/" if n_reactions == 0 else "int ki[], int nu[]"
     )
@@ -2486,7 +2504,7 @@ def ckinu(fstream, mechanism, species_info, reaction_info, write_sk=False):
     cw.writer(fstream, cw.comment("and stoichiometric coefficients. (Eq 50)"))
     cw.writer(
         fstream,
-        f"void {function_prefix}KINU"
+        f"void {cu.get_function_prefix(write_sk)}KINU"
         + cc.sym
         + f"(const int i, int& nspec, {function_args})",
     )
@@ -2536,14 +2554,13 @@ def ckinu(fstream, mechanism, species_info, reaction_info, write_sk=False):
     cw.writer(fstream, "}")
 
 
-def ckkfkr(fstream, mechanism, species_info):
+def ckkfkr(fstream, mech_is_heterogeneous, mech_is_reacting):
     """Write ckkfkr."""
-    n_species = species_info.n_species
-    n_reactions = mechanism.n_reactions
+    comment_str = ", surface coverages" if mech_is_heterogeneous else ""
 
     cw.writer(fstream)
-    cw.writer(fstream, cw.comment("Returns the progress rates of each reactions"))
-    cw.writer(fstream, cw.comment("Given P, T, and mole fractions"))
+    cw.writer(fstream, cw.comment("Returns the progress rates of each reaction"))
+    cw.writer(fstream, cw.comment(f"Given P, T, and mole fractions{comment_str}"))
     cw.writer(
         fstream,
         "void CKKFKR"
@@ -2555,7 +2572,7 @@ def ckkfkr(fstream, mechanism, species_info):
 
     cw.writer(
         fstream,
-        f"amrex::Real c[{n_species}]; " + cw.comment("temporary storage"),
+        "amrex::Real c[NUM_SPECIES]; " + cw.comment("temporary storage"),
     )
     cw.writer(
         fstream,
@@ -2564,13 +2581,36 @@ def ckkfkr(fstream, mechanism, species_info):
         " T); " + cw.comment("1e6 * P/RT so c goes to SI units"),
     )
 
+    if mech_is_heterogeneous:
+        cw.writer(
+            fstream,
+            "amrex::Real site_density = 1e4 * SITE_DENSITY;"
+            + cw.comment("convert to SI (mol/cm^2 to mol/m^2)"),
+        )
+
     # now compute conversion
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("Compute conversion, see Eq 10"))
-    cw.writer(fstream, f"for (int id = 0; id < {n_species}; ++id) {{")
+    cw.writer(fstream, "for (int id = 0; id < NUM_GAS_SPECIES; ++id) {")
     cw.writer(fstream, "c[id] = x[id]*PORT;")
     cw.writer(fstream, "}")
 
+    if mech_is_heterogeneous:
+        cw.writer(fstream)
+        cw.writer(
+            fstream,
+            cw.comment(
+                "Compute surface species concentrations (aka surface coverages)"
+            ),
+        )
+        cw.writer(
+            fstream,
+            cw.comment("Assuming unit site occupancy number for all surface species"),
+        )
+        cw.writer(fstream, "for (int id = NUM_GAS_SPECIES; id < NUM_SPECIES; ++id) {")
+        cw.writer(fstream, "c[id] = x[id]*site_density;")
+        cw.writer(fstream, "}")
+
     # call progressRateFR
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("convert to chemkin units"))
@@ -2578,11 +2618,21 @@ def ckkfkr(fstream, mechanism, species_info):
 
     # convert qdot to chemkin units
     cw.writer(fstream)
-    if n_reactions > 0:
+    if mech_is_reacting:
         cw.writer(fstream, cw.comment("convert to chemkin units"))
-        cw.writer(fstream, f"for (int id = 0; id < {n_reactions}; ++id) {{")
+        cw.writer(fstream, "for (int id = 0; id < NUM_GAS_REACTIONS; ++id) {")
         cw.writer(fstream, "q_f[id] *= 1.0e-6;")
         cw.writer(fstream, "q_r[id] *= 1.0e-6;")
         cw.writer(fstream, "}")
 
+        if mech_is_heterogeneous:
+            cw.writer(fstream)
+            cw.writer(fstream, cw.comment("convert surface qf/qr to chemkin units"))
+            cw.writer(
+                fstream, "for (int id = NUM_GAS_SPECIES; id < NUM_SPECIES; ++id) {"
+            )
+            cw.writer(fstream, "q_f[id] *= 1.0e-4;")
+            cw.writer(fstream, "q_r[id] *= 1.0e-4;")
+            cw.writer(fstream, "}")
+
     cw.writer(fstream, "}")
diff --git a/Support/ceptr/ceptr/converter.py b/Support/ceptr/ceptr/converter.py
index fbdc8d926..b7b099484 100644
--- a/Support/ceptr/ceptr/converter.py
+++ b/Support/ceptr/ceptr/converter.py
@@ -18,6 +18,7 @@
 import ceptr.symbolic_math as csm
 import ceptr.thermo as cth
 import ceptr.transport as ctr
+import ceptr.utilities as cu
 import ceptr.writer as cw
 
 
@@ -33,7 +34,7 @@ def __init__(
         qss_format_input=None,
         qss_symbolic_jacobian=False,
     ):
-        self.mechIsAHetMech = chemistry == "heterogeneous"
+        self.mech_is_heterogeneous = chemistry == "heterogeneous"
 
         self.mechanism = mechanism
         self.interface = interface
@@ -45,7 +46,7 @@ def __init__(
 
         self.mechpath = (
             pathlib.Path(self.interface.source)
-            if self.mechIsAHetMech
+            if self.mech_is_heterogeneous
             else pathlib.Path(self.mechanism.source)
         )
 
@@ -115,7 +116,7 @@ def set_species(self):
 
         self.species_info.n_all_species = (
             self.species_info.n_species + self.interface.n_species
-            if self.mechIsAHetMech
+            if self.mech_is_heterogeneous
             else self.species_info.n_species
         )
 
@@ -196,7 +197,7 @@ def set_species(self):
             "d",
         )
 
-        if self.mechIsAHetMech:
+        if self.mech_is_heterogeneous:
             # Initialize gas-solid interface species
             self.species_info.n_surface_species = self.interface.n_species
             for id, species in enumerate(self.interface.species()):
@@ -235,17 +236,17 @@ def writer(self):
             cri.rmap(cpp, self.reaction_info)
             cri.get_rmap(cpp, self.reaction_info)
             cck.ckinu(cpp, self.mechanism, self.species_info, self.reaction_info)
-            cck.ckkfkr(cpp, self.mechanism, self.species_info)
-            cp.progress_rate_fr(
-                cpp, self.mechanism, self.species_info, self.reaction_info
+            cck.ckkfkr(
+                cpp, self.mech_is_heterogeneous, self.reaction_info.n_reactions > 0
             )
+            cp.progress_rate_fr(cpp, self.species_info, self.reaction_info)
             self.atomic_weight(cpp)
             cck.ckawt(cpp, self.mechanism)
             cck.ckncf(cpp, self.mechanism, self.species_info)
-            cck.cksyme_str(cpp, self.mechanism, self.species_info)
-            cck.cksyms_str(cpp, self.mechanism, self.species_info)
+            cck.cksyme_str(cpp, self.mechanism, self.interface)
+            cck.cksyms_str(cpp, self.species_info, self.mech_is_heterogeneous)
             csp.sparsity(cpp, self.species_info)
-            if self.interface is not None:
+            if self.mech_is_heterogeneous:
                 cck.ckinu(
                     cpp,
                     self.interface,
@@ -253,6 +254,7 @@ def writer(self):
                     self.reaction_info,
                     write_sk=True,
                 )
+                cck.ckncf(cpp, self.interface, self.species_info, write_sk=True)
 
             # This is for the header file
             cw.writer(hdr, "#ifndef MECHANISM_H")
@@ -516,13 +518,15 @@ def formatter(self):
 
     def atomic_weight(self, fstream):
         """Write the atomic weight."""
+        element_names = cu.get_element_names(self.mechanism, self.interface)
+
         cw.writer(fstream)
         cw.writer(fstream, cw.comment("save atomic weights into array"))
         cw.writer(fstream, "void atomicWeight(amrex::Real *  awt)")
         cw.writer(fstream, "{")
-        for elem in self.mechanism.element_names:
-            idx = self.mechanism.element_index(elem)
-            aw = self.mechanism.atomic_weight(elem)
+        for elem in element_names:
+            idx = cu.get_element_id(self.mechanism, self.interface, elem)
+            aw = cu.get_atomic_weight(self.mechanism, self.interface, elem)
             cw.writer(fstream, f"awt[{idx}] = {aw:f}; " + cw.comment(f"{elem}"))
         cw.writer(fstream, "}")
 
@@ -700,7 +704,7 @@ def mechanism_header_includes(self, fstream):
         for elem in self.mechanism.element_names:
             cw.writer(fstream, f"{self.mechanism.element_index(elem)}  {elem}")
 
-        if self.interface is not None:
+        if self.mech_is_heterogeneous:
             n_het_species = self.interface.n_species
             n_het_reactions = self.interface.n_reactions
             all_species_list += self.interface.species_names
@@ -722,8 +726,6 @@ def mechanism_header_includes(self, fstream):
             if s[-1] == "n" or s[-1] == "p" or s == "E":
                 nb_ions += 1
 
-        qssa_str = "QSSA_" if self.species_info.n_qssa_species > 0 else ""
-
         cw.writer(fstream)
         cw.writer(
             fstream,
@@ -732,7 +734,7 @@ def mechanism_header_includes(self, fstream):
         )
         cw.writer(
             fstream,
-            f"#define NUM_{qssa_str}GAS_SPECIES {n_hom_species}"
+            f"#define NUM_GAS_SPECIES {n_hom_species}"
             + cw.comment("Species in the homogeneous phase"),
         )
         cw.writer(
@@ -741,7 +743,13 @@ def mechanism_header_includes(self, fstream):
             + cw.comment("Reactions in the homogeneous phase"),
         )
 
-        if not isinstance(self.interface, type(None)):
+        cw.writer(
+            fstream,
+            f"#define NUM_QSS_GAS_SPECIES {self.species_info.n_qssa_species}"
+            + cw.comment("QSS species in the homogeneous phase"),
+        )
+
+        if self.mech_is_heterogeneous:
             site_density = 0.1 * self.interface.site_density  # Kmol/m**2 to mol/cm**2
 
         cw.writer(fstream)
@@ -771,7 +779,7 @@ def mechanism_header_includes(self, fstream):
         )
         cw.writer(
             fstream,
-            f"#define NUM_SPECIES (NUM_{qssa_str}GAS_SPECIES + NUM_SURFACE_SPECIES)",
+            "#define NUM_SPECIES (NUM_GAS_SPECIES + NUM_SURFACE_SPECIES)",
         )
         cw.writer(
             fstream, "#define NUM_REACTIONS (NUM_GAS_REACTIONS + NUM_SURFACE_REACTIONS)"
diff --git a/Support/ceptr/ceptr/production.py b/Support/ceptr/ceptr/production.py
index d59f241a9..366b5e731 100644
--- a/Support/ceptr/ceptr/production.py
+++ b/Support/ceptr/ceptr/production.py
@@ -1435,16 +1435,14 @@ def production_rate_light(fstream, mechanism, species_info, reaction_info):
     cw.writer(fstream)
 
 
-def progress_rate_fr(fstream, mechanism, species_info, reaction_info):
+def progress_rate_fr(fstream, species_info, reaction_info):
     """Write progress rates."""
-    n_reactions = mechanism.n_reactions
-
-    assert len(reaction_info.index) == 7
+    assert len(reaction_info.index) == 7 or len(reaction_info.index) == 9
 
     cw.writer(fstream)
     cw.writer(fstream, cw.comment("compute the progress rate for each reaction"))
     cw.writer(fstream, cw.comment("USES progressRate : todo switch to GPU"))
-    if n_reactions > 0:
+    if reaction_info.n_reactions > 0:
         cw.writer(
             fstream,
             "void progressRateFR"
@@ -1461,17 +1459,17 @@ def progress_rate_fr(fstream, mechanism, species_info, reaction_info):
 
     cw.writer(fstream, "{")
 
-    if n_reactions > 0:
+    if reaction_info.n_reactions > 0:
         cw.writer(fstream, "const amrex::Real invT = 1.0 / T;")
         cw.writer(fstream, "const amrex::Real logT = log(T);")
 
         cw.writer(fstream, cw.comment("compute the Gibbs free energy"))
-        cw.writer(fstream, f"amrex::Real g_RT[{species_info.n_species}];")
+        cw.writer(fstream, "amrex::Real g_RT[NUM_SPECIES];")
         cw.writer(fstream, "gibbs(g_RT, T);")
         if species_info.n_qssa_species > 0:
             cw.writer(
                 fstream,
-                f"amrex::Real g_RT_qss[{species_info.n_qssa_species}];",
+                "amrex::Real g_RT_qss[NUM_GAS_SPECIES];",
             )
             cw.writer(fstream, "gibbs_qss(g_RT_qss, T);")
 
diff --git a/Support/ceptr/ceptr/utilities.py b/Support/ceptr/ceptr/utilities.py
index 555dfab75..df5e567a5 100644
--- a/Support/ceptr/ceptr/utilities.py
+++ b/Support/ceptr/ceptr/utilities.py
@@ -451,3 +451,41 @@ def enhancement_d(mechanism, species_info, reaction, syms=None):
         return " + ".join(alpha).replace("+ -", "- "), enhancement_smp
     else:
         return " + ".join(alpha).replace("+ -", "- ")
+
+
+def get_function_prefix(is_heterogeneous):
+    """Get the function prefix for homogeneous/heterogeneous mechanisms."""
+    return "S" if is_heterogeneous else "C"
+
+
+def get_phase(is_heterogeneous):
+    """Get the phase for homogeneous/heterogeneous mechanisms."""
+    return "surface" if is_heterogeneous else "gas"
+
+
+def get_element_id(mechanism, interface, element):
+    """Get the element id for elements associated with homogeneous/heterogeneous species."""
+    return (
+        mechanism.element_index(element)
+        if element in mechanism.element_names
+        else interface.element_index(element) + len(mechanism.element_names)
+    )
+
+
+def get_atomic_weight(mechanism, interface, element):
+    """Get atomic weight of element associated with homogeneous/heterogeneous species."""
+    return (
+        mechanism.atomic_weight(element)
+        if element in mechanism.element_names
+        else interface.atomic_weight(element)
+    )
+
+
+def get_element_names(mechanism, interface):
+    """Get the names of all elements associated with homogeneous/heterogeneous species."""
+    surface_elements_set = (
+        set(interface.element_names) - set(mechanism.element_names)
+        if interface is not None
+        else set()
+    )
+    return mechanism.element_names + list(surface_elements_set)