diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp
index 6e281e25e..bc314626d 100644
--- a/Source/Diffusion/ComputeTurbulentViscosity.cpp
+++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp
@@ -151,12 +151,17 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
           Real eps       = std::numeric_limits<Real>::epsilon();
           Real dtheta_dz;
           if (use_most && k==klo) {
+#ifdef ERF_EXPLICIT_MOST_STRESS
+              dtheta_dz = ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
+                          - cell_data(i,j,k  ,RhoTheta_comp)/cell_data(i,j,k  ,Rho_comp) )*dzInv;
+#else
               dtheta_dz = 0.5 * (-3 * cell_data(i,j,k  ,RhoTheta_comp)
                                     / cell_data(i,j,k  ,Rho_comp)
                                 + 4 * cell_data(i,j,k+1,RhoTheta_comp)
                                     / cell_data(i,j,k+1,Rho_comp)
                                 -     cell_data(i,j,k+2,RhoTheta_comp)
                                     / cell_data(i,j,k+2,Rho_comp) ) * dzInv;
+#endif
           } else {
               dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
                                 - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;