diff --git a/src_common/ComputeAverages.cpp b/src_common/ComputeAverages.cpp
index 34791dcc..2b474cc5 100644
--- a/src_common/ComputeAverages.cpp
+++ b/src_common/ComputeAverages.cpp
@@ -95,8 +95,12 @@ void WriteHorizontalAverage(const MultiFab& mf_in, const int& dir, const int& in
 
 void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
                                 const int& dir, const int& incomp,
-                                const int& ncomp)
+                                const int& ncomp, int outcomp)
 {
+    if (outcomp == -1) {
+        outcomp = incomp; // default outcomp is incomp unless specified
+    }
+
     // number of points in the averaging direction
     int npts = n_cells[dir];
 
@@ -166,7 +170,7 @@ void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
         const Array4<Real> mf = mf_out.array(mfi);
 
         for (auto n=0; n<ncomp; ++n) {
-            comp = incomp+n;
+            comp = outcomp+n;
             for (auto k = lo.z; k <= hi.z; ++k) {
             for (auto j = lo.y; j <= hi.y; ++j) {
             for (auto i = lo.x; i <= hi.x; ++i) {
diff --git a/src_common/common_functions.H b/src_common/common_functions.H
index 353d38ef..7cfeba66 100644
--- a/src_common/common_functions.H
+++ b/src_common/common_functions.H
@@ -148,7 +148,7 @@ void WriteHorizontalAverage(const MultiFab& mf_in, const int& dir, const int& in
 
 void WriteHorizontalAverageToMF(const MultiFab& mf_in, MultiFab& mf_out,
                                 const int& dir, const int& incomp,
-                                const int& ncomp);
+                                const int& ncomp, int outcomp=-1);
     
 void ComputeVerticalAverage(const MultiFab & mf, MultiFab & mf_flat, const Geometry & geom,
                             const int& dir, const int& incomp, const int& ncomp,
diff --git a/src_common/common_functions.cpp b/src_common/common_functions.cpp
index e9ca0ad0..bbc4d429 100644
--- a/src_common/common_functions.cpp
+++ b/src_common/common_functions.cpp
@@ -240,6 +240,7 @@ int                        common::plot_means;
 int                        common::plot_vars;
 int                        common::plot_covars;
 int                        common::plot_cross;
+int                        common::plot_deltaY_dir;
 int                        common::particle_motion;
 
 AMREX_GPU_MANAGED amrex::Real common::turb_a;
@@ -608,6 +609,7 @@ void InitializeCommonNamespace() {
     plot_vars = 0;
     plot_covars = 0;
     plot_cross = 0;
+    plot_deltaY_dir = -1;
     particle_motion = 0;
 
     // turblent forcing parameters
@@ -1135,6 +1137,7 @@ void InitializeCommonNamespace() {
     pp.query("plot_vars",plot_vars);
     pp.query("plot_covars",plot_covars);
     pp.query("plot_cross",plot_cross);
+    pp.query("plot_deltaY_dir",plot_deltaY_dir);
     pp.query("particle_motion",particle_motion);
     pp.query("turb_a",turb_a);
     pp.query("turb_b",turb_b);
diff --git a/src_common/common_namespace.H b/src_common/common_namespace.H
index 7e91aa22..b5493d3a 100644
--- a/src_common/common_namespace.H
+++ b/src_common/common_namespace.H
@@ -314,6 +314,7 @@ namespace common {
     extern int                        plot_vars;
     extern int                        plot_covars;
     extern int                        plot_cross;
+    extern int                        plot_deltaY_dir;
     extern int                        particle_motion;
 
     // parameters for turbulent forcing example
diff --git a/src_compressible_stag/writePlotFileStag.cpp b/src_compressible_stag/writePlotFileStag.cpp
index c606b196..820c4cd8 100644
--- a/src_compressible_stag/writePlotFileStag.cpp
+++ b/src_compressible_stag/writePlotFileStag.cpp
@@ -87,7 +87,11 @@ void WritePlotFileStag(int step,
 
 	if (nspec_surfcov>0) nplot += nspec_surfcov*6;
     }
-   
+
+    if (plot_deltaY_dir != -1) {
+        nplot += nspecies;
+    }
+
     amrex::BoxArray ba = cuMeans.boxArray();
     amrex::DistributionMapping dmap = cuMeans.DistributionMap();
 
@@ -257,6 +261,16 @@ void WritePlotFileStag(int step,
 	}
     }
 
+    if (plot_deltaY_dir != -1) {
+        MultiFab Ybar(ba, dmap, nspecies, 0);
+        // Yk is component 6: in prim
+        WriteHorizontalAverageToMF(prim,Ybar,plot_deltaY_dir,6,nspecies,0);
+        Ybar.mult(-1.);
+        amrex::MultiFab::Add(Ybar,prim,6,0,nspecies,0);
+        amrex::MultiFab::Copy(plotfile,Ybar,0,cnt,nspecies,0);
+        cnt+= nspecies;
+    }
+
     // Set variable names
     cnt = 0;
 
@@ -445,6 +459,14 @@ void WritePlotFileStag(int step,
 
     }
 
+    if (plot_deltaY_dir != 1) {
+        x = "deltaYk_";
+        for (i=0; i<nspecies; i++) {
+            varNames[cnt] = x;
+            varNames[cnt++] += 48+i;
+        }
+    }
+
     AMREX_ASSERT(cnt==nplot);
 
     // write a plotfile