diff --git a/src_MFsurfchem/MFsurfchem_functions.cpp b/src_MFsurfchem/MFsurfchem_functions.cpp
index 31eb492d..104d0a84 100644
--- a/src_MFsurfchem/MFsurfchem_functions.cpp
+++ b/src_MFsurfchem/MFsurfchem_functions.cpp
@@ -12,8 +12,8 @@ AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_SPECIES> MFsurfchem::des_rate;
 AMREX_GPU_MANAGED int MFsurfchem::stoch_surfcov0;
 AMREX_GPU_MANAGED int MFsurfchem::stoch_MFsurfchem;
 
-// temperature correction exponent
-#define BETA    0.5
+AMREX_GPU_MANAGED amrex::Real MFsurfchem::k_beta;
+AMREX_GPU_MANAGED amrex::Real MFsurfchem::e_beta;
 
 void InitializeMFSurfchemNamespace()
 {
@@ -63,6 +63,12 @@ void InitializeMFSurfchemNamespace()
 
     stoch_MFsurfchem = 1; // default value
     pp.query("stoch_MFsurfchem",stoch_MFsurfchem);
+
+    k_beta = 0.5; // default value
+    pp.query("k_beta",k_beta);
+
+    e_beta = 0.5; // default value
+    pp.query("e_beta",e_beta);
     return;
 }
 
@@ -146,7 +152,7 @@ void sample_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
 
                     amrex::Real theta = surfcov_arr(i,j,k,m);
 
-                    amrex::Real meanNads = ads_rate_const[m]*dens*(1-sumtheta)*Ntot*dt*pow(tempratio,BETA);
+  		    amrex::Real meanNads = ads_rate_const[m]*dens*(1-sumtheta)*Ntot*dt*pow(tempratio,k_beta);
                     amrex::Real meanNdes = des_rate[m]*theta*Ntot*dt;
 
                     amrex::Real Nads;
@@ -193,8 +199,7 @@ void update_MFsurfchem(MultiFab& cu, MultiFab& prim, MultiFab& surfcov, MultiFab
                 for (int m=0;m<n_ads_spec;m++) {
                     amrex::Real dN = dNadsdes_arr(i,j,k,m);
                     amrex::Real factor1 = molmass[m]/AVONUM/(dx[0]*dx[1]*dx[2]);
-                    amrex::Real factor2 = (BETA*k_B*T_inst+(e0[m]+hcv[m]*T_inst)*molmass[m]/AVONUM)/(dx[0]*dx[1]*dx[2]);
-
+		    amrex::Real factor2 = (e_beta*k_B*T_inst+(e0[m]+hcv[m]*T_inst)*molmass[m]/AVONUM)/(dx[0]*dx[1]*dx[2]);
                     surfcov_arr(i,j,k,m) += dN/Ntot;
                     cu_arr(i,j,k,0) -= factor1*dN;
                     cu_arr(i,j,k,5+m) -= factor1*dN;
diff --git a/src_MFsurfchem/MFsurfchem_namespace.H b/src_MFsurfchem/MFsurfchem_namespace.H
index f09c1ace..1449319b 100644
--- a/src_MFsurfchem/MFsurfchem_namespace.H
+++ b/src_MFsurfchem/MFsurfchem_namespace.H
@@ -11,4 +11,7 @@ namespace MFsurfchem {
     extern AMREX_GPU_MANAGED int stoch_surfcov0;
     extern AMREX_GPU_MANAGED int stoch_MFsurfchem;
 
+    extern AMREX_GPU_MANAGED amrex::Real k_beta;
+    extern AMREX_GPU_MANAGED amrex::Real e_beta;
+
 }
diff --git a/src_compressible_stag/Checkpoint.cpp b/src_compressible_stag/Checkpoint.cpp
index f03802ba..caf1a722 100644
--- a/src_compressible_stag/Checkpoint.cpp
+++ b/src_compressible_stag/Checkpoint.cpp
@@ -45,6 +45,7 @@ void WriteCheckPoint3D(int step,
                        const amrex::MultiFab& surfcov,
                        const amrex::MultiFab& surfcovMeans,
                        const amrex::MultiFab& surfcovVars,
+		       const amrex::MultiFab& surfcovcoVars,
                        const Vector<Real>& spatialCross, int ncross,
                        TurbForcingComp& turbforce)
 
@@ -230,6 +231,8 @@ void WriteCheckPoint3D(int step,
                      amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "surfcovMeans"));
         VisMF::Write(surfcovVars,
                      amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "surfcovVars"));
+	VisMF::Write(surfcovcoVars,
+		     amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "surfcovcoVars"));
     }
 }
 
@@ -612,6 +615,7 @@ void ReadCheckPoint3D(int& step,
                      amrex::MultiFab& surfcov,
                      amrex::MultiFab& surfcovMeans,
                      amrex::MultiFab& surfcovVars,
+		     amrex::MultiFab& surfcovcoVars,
                      Vector<Real>& spatialCross,
                      int ncross,
                      TurbForcingComp& turbforce,
@@ -733,6 +737,7 @@ void ReadCheckPoint3D(int& step,
             surfcov.define(ba,dmap,n_ads_spec,0);
             surfcovMeans.define(ba,dmap,n_ads_spec,0);
             surfcovVars.define(ba,dmap,n_ads_spec,0);
+	    surfcovcoVars.define(ba,dmap,n_ads_spec*6,0);
         }
 
     }
@@ -830,6 +835,7 @@ void ReadCheckPoint3D(int& step,
             for (int m=0;m<n_ads_spec;m++) {
                 surfcovMeans.setVal(0.0);
                 surfcovVars.setVal(0.0);
+		surfcovcoVars.setVal(0.0);
             }
         }
     }
@@ -859,6 +865,7 @@ void ReadCheckPoint3D(int& step,
         if (n_ads_spec>0) {
             Read_Copy_MF_Checkpoint(surfcovMeans,"surfcovMeans",checkpointname,ba_old,dmap_old,n_ads_spec,0);
             Read_Copy_MF_Checkpoint(surfcovVars,"surfcovVars",checkpointname,ba_old,dmap_old,n_ads_spec,0);
+	    Read_Copy_MF_Checkpoint(surfcovcoVars,"surfcovcoVars",checkpointname,ba_old,dmap_old,n_ads_spec*6,0);
         }
     }
 
diff --git a/src_compressible_stag/compressible_functions_stag.H b/src_compressible_stag/compressible_functions_stag.H
index 1ebaeca2..10143477 100644
--- a/src_compressible_stag/compressible_functions_stag.H
+++ b/src_compressible_stag/compressible_functions_stag.H
@@ -31,6 +31,7 @@ void WritePlotFileStag(int step,
                        const MultiFab& surfcov,
                        const MultiFab& surfcovMeans,
                        const MultiFab& surfcovVars,
+		       const MultiFab& surfcovcoVars,
                        const MultiFab& eta, 
                        const MultiFab& kappa,
                        const MultiFab& zeta);
@@ -194,6 +195,7 @@ void WriteCheckPoint3D(int step,
                        const amrex::MultiFab& surfcov,
                        const amrex::MultiFab& surfcovMeans,
                        const amrex::MultiFab& surfcovVars,
+		       const amrex::MultiFab& surfcovcoVars,
                        const Vector<Real>& spatialCross, int ncross,
                        TurbForcingComp& turbforce);
 
@@ -218,6 +220,7 @@ void ReadCheckPoint3D(int& step,
                       amrex::MultiFab& surfcov,
                       amrex::MultiFab& surfcovMeans,
                       amrex::MultiFab& surfcovVars,
+		      amrex::MultiFab& surfcovcoVars,
                       Vector<Real>& spatialCross, int ncross,
                       TurbForcingComp& turbforce,
                       BoxArray& ba, DistributionMapping& dmap);
@@ -322,7 +325,7 @@ void evaluateStatsStag3D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          Vector<Real>& dataSliceMeans_xcross,
                          Vector<Real>& spatialCross3D, const int ncross,
                          const amrex::Box& domain,
@@ -338,7 +341,7 @@ void evaluateStatsStag2D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          MultiFab& spatialCross2D, const int ncross,
                          const int steps,
                          const Geometry& geom);
@@ -352,7 +355,7 @@ void evaluateStatsStag1D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          MultiFab& spatialCross1D, const int ncross,
                          const int steps,
                          const Geometry& geom);
@@ -374,7 +377,7 @@ void EvaluateVarsCoVars(const MultiFab& cons, const MultiFab& consMean, MultiFab
                         const std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                         std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                         MultiFab& coVar,
-                        const MultiFab& theta, const MultiFab& thetaMean, MultiFab& thetaVar,
+                        const MultiFab& theta, const MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                         const int steps);
 
 void GetSliceAverageCross(Vector<Real>& dataAvMeans_x,
diff --git a/src_compressible_stag/main_driver.cpp b/src_compressible_stag/main_driver.cpp
index 1a1a0904..100d4979 100644
--- a/src_compressible_stag/main_driver.cpp
+++ b/src_compressible_stag/main_driver.cpp
@@ -198,6 +198,7 @@ void main_driver(const char* argv)
 
     MultiFab surfcovMeans;  // used in either MFsurfchem or surfchem_mui
     MultiFab surfcovVars;   // used in either MFsurfchem or surfchem_mui
+    MultiFab surfcovcoVars; // used in either MFsurfchem or surfchem_mui
 
     std::array< MultiFab, AMREX_SPACEDIM > velMeans;
     std::array< MultiFab, AMREX_SPACEDIM > velVars;
@@ -482,7 +483,7 @@ void main_driver(const char* argv)
         else {
             ReadCheckPoint3D(step_start, time, statsCount, geom, domain, cu, cuMeans, cuVars, prim,
                              primMeans, primVars, cumom, cumomMeans, cumomVars, 
-                             vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, spatialCross3D, ncross, turbforce, ba, dmap);
+                             vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, surfcovcoVars, spatialCross3D, ncross, turbforce, ba, dmap);
         }
 
         if (reset_stats == 1) statsCount = 1;
@@ -626,8 +627,10 @@ void main_driver(const char* argv)
         if (nspec_surfcov>0) {
             surfcovMeans.define(ba,dmap,nspec_surfcov,0);
             surfcovVars.define(ba,dmap,nspec_surfcov,0);
+	    surfcovcoVars.define(ba,dmap,nspec_surfcov*6,0);
             surfcovMeans.setVal(0.0);
             surfcovVars.setVal(0.0);
+	    surfcovcoVars.setVal(0.0);
         }
 
         for (int d=0; d<AMREX_SPACEDIM; d++) {
@@ -717,7 +720,7 @@ void main_driver(const char* argv)
         
         if (plot_int > 0) {
             WritePlotFileStag(0, 0.0, geom, cu, cuMeans, cuVars, cumom, cumomMeans, cumomVars, 
-                          prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, eta, kappa, zeta);
+                          prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, surfcovcoVars, eta, kappa, zeta);
 #if defined(TURB)
             if (turbForcing > 0) {
                 EvaluateWritePlotFileVelGrad(0, 0.0, geom, vel, vel_decomp);
@@ -1079,6 +1082,7 @@ void main_driver(const char* argv)
             if (nspec_surfcov>0) {
                 surfcovMeans.setVal(0.0);
                 surfcovVars.setVal(0.0);
+		surfcovcoVars.setVal(0.0);
             }
 
             if (do_1D) {
@@ -1101,19 +1105,19 @@ void main_driver(const char* argv)
         if (do_1D) {
             evaluateStatsStag1D(cu, cuMeans, cuVars, prim, primMeans, primVars, vel, 
                                 velMeans, velVars, cumom, cumomMeans, cumomVars, coVars,
-                                surfcov, surfcovMeans, surfcovVars,
+                                surfcov, surfcovMeans, surfcovVars, surfcovcoVars,
                                 spatialCross1D, ncross, statsCount, geom);
         }
         else if (do_2D) {
             evaluateStatsStag2D(cu, cuMeans, cuVars, prim, primMeans, primVars, vel, 
                                 velMeans, velVars, cumom, cumomMeans, cumomVars, coVars,
-                                surfcov, surfcovMeans, surfcovVars,
+                                surfcov, surfcovMeans, surfcovVars, surfcovcoVars,
                                 spatialCross2D, ncross, statsCount, geom);
         }
         else {
             evaluateStatsStag3D(cu, cuMeans, cuVars, prim, primMeans, primVars, vel, 
                                 velMeans, velVars, cumom, cumomMeans, cumomVars, coVars,
-                                surfcov, surfcovMeans, surfcovVars,
+                                surfcov, surfcovMeans, surfcovVars, surfcovcoVars,
                                 dataSliceMeans_xcross, spatialCross3D, ncross, domain,
                                 statsCount, geom);
         }
@@ -1145,7 +1149,7 @@ void main_driver(const char* argv)
             //yzAverage(cuMeans, cuVars, primMeans, primVars, spatialCross,
             //          cuMeansAv, cuVarsAv, primMeansAv, primVarsAv, spatialCrossAv);
             WritePlotFileStag(step, time, geom, cu, cuMeans, cuVars, cumom, cumomMeans, cumomVars,
-                              prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, eta, kappa, zeta);
+                              prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, surfcovcoVars, eta, kappa, zeta);
             
 #if defined(TURB)
             if (turbForcing > 0) {
@@ -1485,7 +1489,7 @@ void main_driver(const char* argv)
             else {
                 WriteCheckPoint3D(step, time, statsCount, geom, cu, cuMeans, cuVars, prim,
                                   primMeans, primVars, cumom, cumomMeans, cumomVars, 
-                                  vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, spatialCross3D, ncross, turbforce);
+                                  vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, surfcovcoVars, spatialCross3D, ncross, turbforce);
             }
         }
 
diff --git a/src_compressible_stag/statsStag.cpp b/src_compressible_stag/statsStag.cpp
index 33063f7f..f6491113 100644
--- a/src_compressible_stag/statsStag.cpp
+++ b/src_compressible_stag/statsStag.cpp
@@ -16,7 +16,7 @@ void evaluateStatsStag3D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          Vector<Real>& dataSliceMeans_xcross,
                          Vector<Real>& spatialCross3D, const int ncross,
                          const amrex::Box& domain,
@@ -35,7 +35,7 @@ void evaluateStatsStag3D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
     //// Evaluate Variances and Covariances
     if ((plot_vars) or (plot_covars)) {
         EvaluateVarsCoVars(cons,consMean,consVar,prim_in,primMean,primVar,velMean,velVar,
-                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,steps);
+                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,thetacoVar,steps);
         consVar.FillBoundary(geom.periodicity());
         primVar.FillBoundary(geom.periodicity());
     }
@@ -81,7 +81,7 @@ void evaluateStatsStag2D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          MultiFab& /*spatialCross2D*/, const int /*ncross*/,
                          const int steps,
                          const Geometry& geom)
@@ -98,7 +98,7 @@ void evaluateStatsStag2D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
     //// Evaluate Variances and Covariances
     if ((plot_vars) or (plot_covars)) {
         EvaluateVarsCoVars(cons,consMean,consVar,prim_in,primMean,primVar,velMean,velVar,
-                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,steps);
+                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,thetacoVar,steps);
         consVar.FillBoundary(geom.periodicity());
         primVar.FillBoundary(geom.periodicity());
     }
@@ -130,7 +130,7 @@ void evaluateStatsStag1D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                          std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                          MultiFab& coVar, 
-                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar,
+                         MultiFab& theta, MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                          MultiFab& spatialCross1D, const int ncross,
                          const int steps,
                          const Geometry& geom)
@@ -147,7 +147,7 @@ void evaluateStatsStag1D(MultiFab& cons, MultiFab& consMean, MultiFab& consVar,
     //// Evaluate Variances and Covariances
     if ((plot_vars) or (plot_covars)) {
         EvaluateVarsCoVars(cons,consMean,consVar,prim_in,primMean,primVar,velMean,velVar,
-                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,steps);
+                           cumom,cumomMean,cumomVar,coVar,theta,thetaMean,thetaVar,thetacoVar,steps);
         consVar.FillBoundary(geom.periodicity());
         primVar.FillBoundary(geom.periodicity());
     }
@@ -362,7 +362,7 @@ void EvaluateVarsCoVars(const MultiFab& cons, const MultiFab& consMean, MultiFab
                         const std::array<MultiFab, AMREX_SPACEDIM>& cumomMean,
                         std::array<MultiFab, AMREX_SPACEDIM>& cumomVar,
                         MultiFab& coVar,
-                        const MultiFab& theta, const MultiFab& thetaMean, MultiFab& thetaVar,
+                        const MultiFab& theta, const MultiFab& thetaMean, MultiFab& thetaVar, MultiFab& thetacoVar,
                         const int steps)
 {
     BL_PROFILE_VAR("EvaluateVarsCoVars()",EvaluateVarsCoVars);
@@ -397,6 +397,11 @@ void EvaluateVarsCoVars(const MultiFab& cons, const MultiFab& consMean, MultiFab
         const Array4<      Real> momyvars  = cumomVar[1].array(mfi);
         const Array4<      Real> momzvars  = cumomVar[2].array(mfi);
 
+        const Array4<const Real> surfcov      = (nspec_surfcov>0) ? theta.array(mfi) : cons.array(mfi);
+        const Array4<const Real> surfcovmeans = (nspec_surfcov>0) ? thetaMean.array(mfi) : consMean.array(mfi);
+        const Array4<      Real> surfcovvars  = (nspec_surfcov>0) ? thetaVar.array(mfi) : consVar.array(mfi);
+        const Array4<      Real> surfcovcoVars= (nspec_surfcov>0) ? thetacoVar.array(mfi) : consVar.array(mfi);
+
         // update momentum and velocity variances
         amrex::ParallelFor(tbx, tby, tbz,
         [=] AMREX_GPU_DEVICE (int i, int j, int k) 
@@ -460,6 +465,7 @@ void EvaluateVarsCoVars(const MultiFab& cons, const MultiFab& consMean, MultiFab
         const Array4<const Real> surfcov      = (nspec_surfcov>0) ? theta.array(mfi) : cons.array(mfi);
         const Array4<const Real> surfcovmeans = (nspec_surfcov>0) ? thetaMean.array(mfi) : consMean.array(mfi);
         const Array4<      Real> surfcovvars  = (nspec_surfcov>0) ? thetaVar.array(mfi) : consVar.array(mfi);
+        const Array4<      Real> surfcovcoVars= (nspec_surfcov>0) ? thetacoVar.array(mfi) : consVar.array(mfi);
         
         // other cell-centered variances and covariances (temperature fluctuation and variance)
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -570,6 +576,12 @@ void EvaluateVarsCoVars(const MultiFab& cons, const MultiFab& consMean, MultiFab
                 for (int m=0; m<nspec_surfcov; ++m) {
                     Real delsurfcov = surfcov(i,j,k,m) - surfcovmeans(i,j,k,m);
                     surfcovvars(i,j,k,m) = (surfcovvars(i,j,k,m)*stepsminusone + delsurfcov*delsurfcov)*stepsinv;
+		    surfcovcoVars(i,j,k,6*m) = (surfcovcoVars(i,j,k,6*m)*stepsminusone + delrhoYk[m]*deltemp)*stepsinv;        // <rhoYk T>
+		    surfcovcoVars(i,j,k,6*m+1) = (surfcovcoVars(i,j,k,6*m+1)*stepsminusone + delsurfcov*delrhoYk[m])*stepsinv; // <theta rhoYk>
+		    surfcovcoVars(i,j,k,6*m+2) = (surfcovcoVars(i,j,k,6*m+2)*stepsminusone + delsurfcov*delvelx)*stepsinv;     // <theta vx>
+		    surfcovcoVars(i,j,k,6*m+3) = (surfcovcoVars(i,j,k,6*m+3)*stepsminusone + delsurfcov*delvely)*stepsinv;     // <theta vy>
+		    surfcovcoVars(i,j,k,6*m+4) = (surfcovcoVars(i,j,k,6*m+4)*stepsminusone + delsurfcov*delvelz)*stepsinv;     // <theta vz>
+		    surfcovcoVars(i,j,k,6*m+5) = (surfcovcoVars(i,j,k,6*m+5)*stepsminusone + delsurfcov*deltemp)*stepsinv;     // <theta T>
                 }
             }
         });
diff --git a/src_compressible_stag/writePlotFileStag.cpp b/src_compressible_stag/writePlotFileStag.cpp
index 2be58ab8..c606b196 100644
--- a/src_compressible_stag/writePlotFileStag.cpp
+++ b/src_compressible_stag/writePlotFileStag.cpp
@@ -23,6 +23,7 @@ void WritePlotFileStag(int step,
                        const amrex::MultiFab& surfcov,
                        const amrex::MultiFab& surfcovMeans,
                        const amrex::MultiFab& surfcovVars,
+		       const amrex::MultiFab& surfcovcoVars,
                        const amrex::MultiFab& eta,
                        const amrex::MultiFab& kappa,
                        const amrex::MultiFab& zeta)
@@ -83,6 +84,8 @@ void WritePlotFileStag(int step,
     // co-variances -- see the list in main_driver
     if (plot_covars == 1) {
         nplot += 26;
+
+	if (nspec_surfcov>0) nplot += nspec_surfcov*6;
     }
    
     amrex::BoxArray ba = cuMeans.boxArray();
@@ -246,6 +249,12 @@ void WritePlotFileStag(int step,
         numvars = 26;
         amrex::MultiFab::Copy(plotfile,coVars,0,cnt,numvars,0);
         cnt+=numvars;
+
+	if (nspec_surfcov>0) {
+	    numvars = nspec_surfcov*6;
+	    amrex::MultiFab::Copy(plotfile,surfcovcoVars,0,cnt,numvars,0);
+	    cnt+=numvars;
+	}
     }
 
     // Set variable names
@@ -422,6 +431,18 @@ void WritePlotFileStag(int step,
         varNames[cnt++] = "YkH-vx";
         varNames[cnt++] = "rhoYkL-vx";
         varNames[cnt++] = "rhoYkH-vx";
+
+	if (nspec_surfcov>0) {
+	    for (i=0; i<nspec_surfcov; i++) {
+		varNames[cnt++] = "rhoYk-T";
+	        varNames[cnt++] = "theta-rhoYk";
+	        varNames[cnt++] = "theta-vx";
+	        varNames[cnt++] = "theta-vy";
+	        varNames[cnt++] = "theta-vz";
+	        varNames[cnt++] = "theta-T";
+	    }
+	}
+
     }
 
     AMREX_ASSERT(cnt==nplot);