diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp
index 51d91dc88..dec2ac7f7 100644
--- a/Source/BoundaryConditions/ABLMost.cpp
+++ b/Source/BoundaryConditions/ABLMost.cpp
@@ -130,9 +130,9 @@ ABLMost::impose_most_bcs (const int& lev,
 #endif
                           MultiFab* z_phys)
 {
-    const int zlo = 0;
+    const int klo = 0;
     if (flux_type == FluxCalcType::MOENG) {
-        moeng_flux flux_comp(zlo);
+        moeng_flux flux_comp(klo);
         compute_most_bcs(lev, mfs,
 #ifdef ERF_EXPLICIT_MOST_STRESS
                          xzmom_flux, xzmom_flux,
@@ -143,7 +143,7 @@ ABLMost::impose_most_bcs (const int& lev,
 #endif
                          z_phys, m_geom[lev].CellSize(2), flux_comp);
     } else if (flux_type == FluxCalcType::DONELAN) {
-        donelan_flux flux_comp(zlo);
+        donelan_flux flux_comp(klo);
         compute_most_bcs(lev, mfs,
 #ifdef ERF_EXPLICIT_MOST_STRESS
                          xzmom_flux, xzmom_flux,
@@ -154,7 +154,7 @@ ABLMost::impose_most_bcs (const int& lev,
 #endif
                          z_phys, m_geom[lev].CellSize(2), flux_comp);
     } else {
-        custom_flux flux_comp(zlo);
+        custom_flux flux_comp(klo);
         compute_most_bcs(lev, mfs,
 #ifdef ERF_EXPLICIT_MOST_STRESS
                          xzmom_flux, xzmom_flux,
@@ -191,13 +191,13 @@ ABLMost::compute_most_bcs (const int& lev,
                            const Real& dz_no_terrain,
                            const FluxCalc& flux_comp)
 {
-    const int zlo   = 0;
+    const int klo   = 0;
     const int icomp = 0;
     for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
     {
         // TODO: No LSM lateral ghost cells, should this change?
         // Valid CC box
-        Box vbx = mfi.validbox(); vbx.makeSlab(2,zlo-1);
+        Box vbx = mfi.validbox(); vbx.makeSlab(2,klo-1);
 
         Box vbxx = surroundingNodes(vbx,0);
         Box vbxy = surroundingNodes(vbx,1);
@@ -248,36 +248,37 @@ ABLMost::compute_most_bcs (const int& lev,
 
             if (var_idx == Vars::cons) {
                 Box b2d = bx;
-                b2d.setBig(2,zlo-1);
+                b2d.setBig(2,klo-1);
                 int n = RhoTheta_comp;
                 ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
+                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
+
 #ifdef ERF_EXPLICIT_MOST_STRESS
-                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,zlo+1) - zphys_arr(i,j,zlo) ) : dz_no_terrain;
-#endif
+                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain;
 
                     // This is the _kinematic_ heat flux [K m/s]
+                    Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz, dz1,
+                                                          cons_arr, velx_arr, vely_arr,
+                                                          umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr,
+                                                          dest_arr);
+#else
                     Real Tflux = flux_comp.compute_t_flux(i, j, k, n, icomp, dz,
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                                                          dz1,
-#endif
                                                           cons_arr, velx_arr, vely_arr,
-#ifndef ERF_EXPLICIT_MOST_STRESS
                                                           eta_arr,
-#endif
                                                           umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr,
                                                           dest_arr);
+#endif
 
 
                     // TODO: make sure not to double-count surface heat flux if using a LSM
-                    int is_land = (lmask_arr) ? lmask_arr(i,j,zlo) : 1;
+                    int is_land = (lmask_arr) ? lmask_arr(i,j,klo) : 1;
                     if (is_land && lsm_flux_arr && vbx.contains(i,j,k)) {
-                        lsm_flux_arr(i,j,zlo) = Tflux;
+                        lsm_flux_arr(i,j,klo) = Tflux;
                     }
 #ifdef ERF_EXPLICIT_MOST_STRESS
-                    else if ((k == zlo-1) && vbx.contains(i,j,k)) {
-                        hfx_arr(i,j,zlo) = Tflux;
+                    else if ((k == klo-1) && vbx.contains(i,j,k)) {
+                        hfx_arr(i,j,klo) = Tflux;
                     }
 #endif
                 });
@@ -287,83 +288,74 @@ ABLMost::compute_most_bcs (const int& lev,
                     n = RhoQ1_comp;
                     ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                     {
-                        Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
+                        Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
 #ifdef ERF_EXPLICIT_MOST_STRESS
-                        Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,zlo+1) - zphys_arr(i,j,zlo) ) : dz_no_terrain;
-#endif
-
+                        Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain;
+                        Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz, dz1,
+                                                              cons_arr, velx_arr, vely_arr,
+                                                              umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr,
+                                                              dest_arr);
+#else
                         Real Qflux = flux_comp.compute_q_flux(i, j, k, n, icomp, dz,
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                                                              dz1,
-#endif
                                                               cons_arr, velx_arr, vely_arr,
-#ifndef ERF_EXPLICIT_MOST_STRESS
                                                               eta_arr,
-#endif
                                                               umm_arr, tm_arr, u_star_arr, q_star_arr, t_surf_arr,
                                                               dest_arr);
+#endif
                     });
                 }
 
             } else if (var_idx == Vars::xvel) {
 
                 Box xb2d = surroundingNodes(bx,0);
-                xb2d.setBig(2,zlo-1);
+                xb2d.setBig(2,klo-1);
 
                 ParallelFor(xb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
+                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
 #ifdef ERF_EXPLICIT_MOST_STRESS
-                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,zlo+1) - zphys_arr(i,j,zlo) ) : dz_no_terrain;
-#endif
-
+                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain;
+                    Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz, dz1,
+                                                            cons_arr, velx_arr, vely_arr,
+                                                            umm_arr, um_arr, u_star_arr,
+                                                            dest_arr);
+                    if ((k == klo-1) && vbxx.contains(i,j,k)) {
+                        t13_arr(i,j,klo) = -stressx;
+                        if (t31_arr) t31_arr(i,j,klo) = -stressx;
+                    }
+#else
                     Real stressx = flux_comp.compute_u_flux(i, j, k, icomp, dz,
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                                                            dz1,
-#endif
                                                             cons_arr, velx_arr, vely_arr,
-#ifndef ERF_EXPLICIT_MOST_STRESS
                                                             eta_arr,
-#endif
                                                             umm_arr, um_arr, u_star_arr,
                                                             dest_arr);
-
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                    if ((k == zlo-1) && vbxx.contains(i,j,k)) {
-                        t13_arr(i,j,zlo) = -stressx;
-                        if (t31_arr) t31_arr(i,j,zlo) = -stressx;
-                    }
 #endif
                 });
 
             } else if (var_idx == Vars::yvel) {
 
                 Box yb2d = surroundingNodes(bx,1);
-                yb2d.setBig(2,zlo-1);
+                yb2d.setBig(2,klo-1);
 
                 ParallelFor(yb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
                 {
-                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
+                    Real dz = (zphys_arr) ? ( zphys_arr(i,j,klo) - zphys_arr(i,j,klo-1) ) : dz_no_terrain;
 #ifdef ERF_EXPLICIT_MOST_STRESS
-                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,zlo+1) - zphys_arr(i,j,zlo) ) : dz_no_terrain;
-#endif
-
+                    Real dz1 = (zphys_arr) ? ( zphys_arr(i,j,klo+1) - zphys_arr(i,j,klo) ) : dz_no_terrain;
+                    Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz, dz1,
+                                                            cons_arr, velx_arr, vely_arr,
+                                                            umm_arr, vm_arr, u_star_arr,
+                                                            dest_arr);
+                    if ((k == klo-1) && vbxy.contains(i,j,k)) {
+                        t23_arr(i,j,klo) = -stressy;
+                        if (t32_arr) t32_arr(i,j,klo) = -stressy;
+                    }
+#else
                     Real stressy = flux_comp.compute_v_flux(i, j, k, icomp, dz,
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                                                            dz1,
-#endif
                                                             cons_arr, velx_arr, vely_arr,
-#ifndef ERF_EXPLICIT_MOST_STRESS
                                                             eta_arr,
-#endif
                                                             umm_arr, vm_arr, u_star_arr,
                                                             dest_arr);
-
-#ifdef ERF_EXPLICIT_MOST_STRESS
-                    if ((k == zlo-1) && vbxy.contains(i,j,k)) {
-                        t23_arr(i,j,zlo) = -stressy;
-                        if (t32_arr) t32_arr(i,j,zlo) = -stressy;
-                    }
 #endif
                 });
             }
diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H
index 00e6e7015..1c58f53f6 100644
--- a/Source/BoundaryConditions/MOSTStress.H
+++ b/Source/BoundaryConditions/MOSTStress.H
@@ -618,14 +618,6 @@ struct moeng_flux
     {
         amrex::Real rho, eta, qv;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         int ic, jc;
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
@@ -642,6 +634,11 @@ struct moeng_flux
 #ifdef ERF_EXPLICIT_MOST_STRESS
         amrex::Abort("Explicit MOST stress not implemented for Moeng moisture flux");
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz);
@@ -688,14 +685,6 @@ struct moeng_flux
         iy = iy > ubound(vely_arr).x   ? ubound(vely_arr).x   : iy;
         jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
         ic = ic > ubound(cons_arr).x ? ubound(cons_arr).x : ic;
@@ -724,6 +713,11 @@ struct moeng_flux
         amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz;
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         // Note: Kh = eta/rho
@@ -766,14 +760,6 @@ struct moeng_flux
         jy = j  < lbound(vely_arr).y   ? lbound(vely_arr).y   : j;
         jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
-        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
         jc = j  < lbound(cons_arr).y   ? lbound(cons_arr).y   : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -804,6 +790,11 @@ struct moeng_flux
         amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(velx - ugrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
+        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5 *(  eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
                       + eta_arr(ie  ,je,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);
@@ -844,14 +835,6 @@ struct moeng_flux
         int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
         int jxhi = j >  ubound(velx_arr).y ? ubound(velx_arr).y : j;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
-        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x   ? lbound(cons_arr).x   : i;
         jc = j  < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -882,6 +865,11 @@ struct moeng_flux
         amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(vely - vgrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
+        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5*(  eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
                      + eta_arr(ie,je  ,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);
@@ -934,14 +922,6 @@ struct donelan_flux
     {
         amrex::Real rho, eta, qv;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         int ic, jc;
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
@@ -958,6 +938,11 @@ struct donelan_flux
 #ifdef ERF_EXPLICIT_MOST_STRESS
         amrex::Abort("Explicit MOST stress not implemented for Donelan moisture flux");
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz);
@@ -993,14 +978,6 @@ struct donelan_flux
     {
         amrex::Real rho, theta, eta;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         int ic, jc;
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
@@ -1022,6 +999,11 @@ struct donelan_flux
         amrex::Real thetagrad = (cons_arr(ic,jc,zlo+1,RhoTheta_comp) - cons_arr(ic,jc,zlo,RhoTheta_comp)) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - thetagrad * deltaz;
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   =  eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         // Note: Kh = eta/rho
@@ -1064,14 +1046,6 @@ struct donelan_flux
         jy = j  < lbound(vely_arr).y   ? lbound(vely_arr).y   : j;
         jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
-        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
         jc = j  < lbound(cons_arr).y   ? lbound(cons_arr).y   : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -1102,6 +1076,11 @@ struct donelan_flux
         amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(velx - ugrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
+        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5 *(  eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
                       + eta_arr(ie  ,je,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);
@@ -1142,14 +1121,6 @@ struct donelan_flux
         int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
         int jxhi = j >  ubound(velx_arr).y ? ubound(velx_arr).y : j;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
-        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x   ? lbound(cons_arr).x   : i;
         jc = j  < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -1180,6 +1151,11 @@ struct donelan_flux
         amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(vely - vgrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
+        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5*(  eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
                      + eta_arr(ie,je  ,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);
@@ -1231,14 +1207,6 @@ struct custom_flux
     {
         amrex::Real rho, eta, qv;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         int ic, jc;
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
@@ -1255,6 +1223,11 @@ struct custom_flux
 #ifdef ERF_EXPLICIT_MOST_STRESS
         amrex::Abort("Explicit MOST stress not implemented for custom moisture flux");
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Q_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         dest_arr(i,j,k,icomp+n) = rho*(qv - moflux*rho/eta*deltaz);
@@ -1290,14 +1263,6 @@ struct custom_flux
     {
         amrex::Real rho, eta, theta;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
-        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
-        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
-        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
-#endif
-
         int ic, jc;
         ic = i  < lbound(cons_arr).x ? lbound(cons_arr).x : i;
         jc = j  < lbound(cons_arr).y ? lbound(cons_arr).y : j;
@@ -1318,6 +1283,11 @@ struct custom_flux
                                 - cons_arr(ic,jc,zlo  ,RhoTheta_comp)/cons_arr(ic,jc,zlo  ,Rho_comp)) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp+n) = cons_arr(ic,jc,zlo,RhoTheta_comp) - rho*thetagrad * deltaz;
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x ? lbound(eta_arr).x : i;
+        je = j  < lbound(eta_arr).y ? lbound(eta_arr).y : j;
+        ie = ie > ubound(eta_arr).x ? ubound(eta_arr).x : ie;
+        je = je > ubound(eta_arr).y ? ubound(eta_arr).y : je;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
         eta   = amrex::max(eta,eta_eps);
         // Note: Kh = eta/rho
@@ -1360,14 +1330,6 @@ struct custom_flux
         jy = j  < lbound(vely_arr).y   ? lbound(vely_arr).y   : j;
         jy = jy > ubound(vely_arr).y-1 ? ubound(vely_arr).y-1 : jy;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
-        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x+1 ? lbound(cons_arr).x+1 : i;
         jc = j  < lbound(cons_arr).y   ? lbound(cons_arr).y   : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -1389,6 +1351,11 @@ struct custom_flux
         amrex::Real ugrad = (velx_arr(i,j,zlo+1) - velx) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(velx - ugrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x+1 ? lbound(eta_arr).x+1 : i;
+        je = j  < lbound(eta_arr).y   ? lbound(eta_arr).y   : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5 *(  eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
                       + eta_arr(ie  ,je,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);
@@ -1429,14 +1396,6 @@ struct custom_flux
         int jxlo = j <= lbound(velx_arr).y ? lbound(velx_arr).y : j-1;
         int jxhi = j >  ubound(velx_arr).y ? ubound(velx_arr).y : j;
 
-#ifndef ERF_EXPLICIT_MOST_STRESS
-        int ie, je;
-        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
-        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
-        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
-        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
-#endif
-
         ic = i  < lbound(cons_arr).x   ? lbound(cons_arr).x   : i;
         jc = j  < lbound(cons_arr).y+1 ? lbound(cons_arr).y+1 : j;
         ic = ic > ubound(cons_arr).x   ? ubound(cons_arr).x   : ic;
@@ -1458,6 +1417,11 @@ struct custom_flux
         amrex::Real vgrad = (vely_arr(i,j,zlo+1) - vely) / (0.5*(dz+dz1));
         dest_arr(i,j,k,icomp) = rho*(vely - vgrad * deltaz);
 #else
+        int ie, je;
+        ie = i  < lbound(eta_arr).x   ? lbound(eta_arr).x   : i;
+        je = j  < lbound(eta_arr).y+1 ? lbound(eta_arr).y+1 : j;
+        ie = ie > ubound(eta_arr).x   ? ubound(eta_arr).x   : ie;
+        je = je > ubound(eta_arr).y   ? ubound(eta_arr).y   : je;
         eta   = 0.5*(  eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
                      + eta_arr(ie,je  ,zlo,EddyDiff::Mom_v) );
         eta   = amrex::max(eta,eta_eps);