diff --git a/Source/MacProj.cpp b/Source/MacProj.cpp
index 81b047e4..c69bb01a 100644
--- a/Source/MacProj.cpp
+++ b/Source/MacProj.cpp
@@ -561,6 +561,14 @@ MacProj::mac_sync_compute (int                   level,
         forcing_term = std::make_unique<MultiFab>(grids, dmap, num_state_comps, NavierStokesBase::nghost_force());
         divu_fp.reset(ns_level.getDivCond(NavierStokesBase::nghost_force(),prev_time));
 
+        // Get divu to time n+1/2, exactly as velocity_advection and
+        // scalar_advection do before predicting the edge states, so that the
+        // sync rebuilds the edge states the way the advance made them.
+        {
+            std::unique_ptr<MultiFab> dsdt(ns_level.getDsdt(NavierStokesBase::nghost_force(),prev_time));
+            MultiFab::Saxpy(*divu_fp, 0.5*dt, *dsdt, 0, 0, 1, NavierStokesBase::nghost_force());
+        }
+
         MultiFab& Gp = ns_level.get_old_data(Gradp_Type);
 
         visc_terms.setVal(0.0); // Initialize to make calls below safe
@@ -770,12 +778,20 @@ MacProj::mac_sync_compute (int                    level,
     bool do_crse_add = false;
     bool do_fine_add = update_fluxreg;
 
+    //
+    // ComputeAofs always builds an Array4 from the state MultiFab, even though
+    // it does not use the values when the edge states are known. So it must be
+    // handed a defined MultiFab; an alias of Sync has the right BoxArray and
+    // DistributionMap.
+    //
+    MultiFab Smf(Sync, amrex::make_alias, Sync_indx, ncomp);
+
     //
     // Compute the mac sync correction.
     //
     ns_level.ComputeAofs(Sync, /*Ssync_comp*/ Sync_indx, state_comp,
                          ncomp,
-                         /*State*/ MultiFab(), /*S_comp*/ int(),//not used when known_edgestates
+                         /*State*/ Smf, /*S_comp*/ 0,           //not used when known_edgestates
                          /*forcing*/ nullptr, /*f_comp*/ int(), //not used when known_edgestates
                          /*constraint divU*/ nullptr,           //not used when known_edgestates
                          fluxes, /*flux_comp*/ 0,
@@ -988,6 +1004,11 @@ MacProj::test_umac_periodic (int       level,
     Vector<IntVect>         pshifts(27);
     std::vector< std::pair<int,Box> > isects;
 
+    //
+    // MultiFabCopyDescriptor and the FArrayBox comparison below both operate
+    // on the host, so give them host-accessible copies of u_mac.
+    //
+    Array<MultiFab,AMREX_SPACEDIM> h_umac;
 
     for (int dim = 0; dim < AMREX_SPACEDIM; dim++)
     {
@@ -995,7 +1016,13 @@ MacProj::test_umac_periodic (int       level,
         {
             Box eDomain = amrex::surroundingNodes(geom.Domain(),dim);
 
-            mfid[dim] = mfcd.RegisterMultiFab(&u_mac[dim]);
+            h_umac[dim].define(u_mac[dim].boxArray(), u_mac[dim].DistributionMap(),
+                               u_mac[dim].nComp(), u_mac[dim].nGrowVect(),
+                               MFInfo().SetArena(The_Pinned_Arena()));
+            amrex::dtoh_memcpy(h_umac[dim], u_mac[dim]);
+            Gpu::streamSynchronize();
+
+            mfid[dim] = mfcd.RegisterMultiFab(&h_umac[dim]);
 
             // How to combine pirm into one global pirm?
             // don't think std::vector::push_back() is thread safe
@@ -1061,11 +1088,11 @@ MacProj::test_umac_periodic (int       level,
         AMREX_ASSERT(pirm_i.m_srcBox.sameSize(pirm_i.m_dstBox));
         AMREX_ASSERT(u_mac[dim].DistributionMap()[pirm_i.m_idx] == ParallelDescriptor::MyProc());
 
-        diff.resize(pirm_i.m_srcBox, 1);
+        diff.resize(pirm_i.m_srcBox, 1, The_Pinned_Arena());
 
         mfcd.FillFab(mfid[dim], pirm_i.m_fbid, diff);
 
-        diff.minus<RunOn::Host>(u_mac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);
+        diff.minus<RunOn::Host>(h_umac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);
 
         const Real max_norm = diff.norm<RunOn::Host>(0);
 
diff --git a/Source/NS_derive.cpp b/Source/NS_derive.cpp
index ca60ed8d..e7bf02f0 100644
--- a/Source/NS_derive.cpp
+++ b/Source/NS_derive.cpp
@@ -131,25 +131,51 @@ namespace derive_functions
         // Need to check if there are covered cells in neighbours --
         // -- if so, use one-sided difference computation (but still quadratic)
         if (!flag_fab(i,j,k).isConnected( 1,0,0)) {
-          vx = - (c0 * dat_arr(i  ,j,k,1)
-              + c1 * dat_arr(i-1,j,k,1)
-              + c2 * dat_arr(i-2,j,k,1)) * idx;
+          if (flag_fab(i,j,k).isConnected(-1,0,0)) {
+            if (flag_fab(i-1,j,k).isConnected(-1,0,0)) {
+              vx = - (c0 * dat_arr(i  ,j,k,1)
+                  + c1 * dat_arr(i-1,j,k,1)
+                  + c2 * dat_arr(i-2,j,k,1)) * idx;
+            } else {
+              // Only one uncovered cell that way, drop to linear
+              vx = (dat_arr(i  ,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
+            }
+          }
+          // Covered on both sides, leave the derivative at zero
         } else if (!flag_fab(i,j,k).isConnected(-1,0,0)) {
-          vx = (c0 * dat_arr(i  ,j,k,1)
-            + c1 * dat_arr(i+1,j,k,1)
-            + c2 * dat_arr(i+2,j,k,1)) * idx;
+          if (flag_fab(i+1,j,k).isConnected( 1,0,0)) {
+            vx = (c0 * dat_arr(i  ,j,k,1)
+              + c1 * dat_arr(i+1,j,k,1)
+              + c2 * dat_arr(i+2,j,k,1)) * idx;
+          } else {
+            // Only one uncovered cell that way, drop to linear
+            vx = (dat_arr(i+1,j,k,1) - dat_arr(i  ,j,k,1)) * idx;
+          }
         } else {
           vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
         }
         // Do the same in y-direction
         if (!flag_fab(i,j,k).isConnected( 0,1,0)) {
-          uy = - (c0 * dat_arr(i,j  ,k,0)
-              + c1 * dat_arr(i,j-1,k,0)
-              + c2 * dat_arr(i,j-2,k,0)) * idy;
+          if (flag_fab(i,j,k).isConnected(0,-1,0)) {
+            if (flag_fab(i,j-1,k).isConnected(0,-1,0)) {
+              uy = - (c0 * dat_arr(i,j  ,k,0)
+                  + c1 * dat_arr(i,j-1,k,0)
+                  + c2 * dat_arr(i,j-2,k,0)) * idy;
+            } else {
+              // Only one uncovered cell that way, drop to linear
+              uy = (dat_arr(i,j  ,k,0) - dat_arr(i,j-1,k,0)) * idy;
+            }
+          }
+          // Covered on both sides, leave the derivative at zero
         } else if (!flag_fab(i,j,k).isConnected(0,-1,0)) {
-          uy = (c0 * dat_arr(i,j  ,k,0)
-            + c1 * dat_arr(i,j+1,k,0)
-            + c2 * dat_arr(i,j+2,k,0)) * idy;
+          if (flag_fab(i,j+1,k).isConnected( 0,1,0)) {
+            uy = (c0 * dat_arr(i,j  ,k,0)
+              + c1 * dat_arr(i,j+1,k,0)
+              + c2 * dat_arr(i,j+2,k,0)) * idy;
+          } else {
+            // Only one uncovered cell that way, drop to linear
+            uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j  ,k,0)) * idy;
+          }
         } else {
           uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy;
         }
@@ -167,20 +193,35 @@ namespace derive_functions
         // -- if so, use one-sided difference computation (but still quadratic)
         if (!flag_fab(i,j,k).isConnected( 1,0,0)) {
           // Covered cell to the right, go fish left
-          vx = - (c0 * dat_arr(i  ,j,k,1)
-              + c1 * dat_arr(i-1,j,k,1)
-              + c2 * dat_arr(i-2,j,k,1)) * idx;
-          wx = - (c0 * dat_arr(i  ,j,k,2)
-              + c1 * dat_arr(i-1,j,k,2)
-              + c2 * dat_arr(i-2,j,k,2)) * idx;
+          if (flag_fab(i,j,k).isConnected(-1,0,0)) {
+            if (flag_fab(i-1,j,k).isConnected(-1,0,0)) {
+              vx = - (c0 * dat_arr(i  ,j,k,1)
+                  + c1 * dat_arr(i-1,j,k,1)
+                  + c2 * dat_arr(i-2,j,k,1)) * idx;
+              wx = - (c0 * dat_arr(i  ,j,k,2)
+                  + c1 * dat_arr(i-1,j,k,2)
+                  + c2 * dat_arr(i-2,j,k,2)) * idx;
+            } else {
+              // Only one uncovered cell that way, drop to linear
+              vx = (dat_arr(i  ,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
+              wx = (dat_arr(i  ,j,k,2) - dat_arr(i-1,j,k,2)) * idx;
+            }
+          }
+          // Covered on both sides, leave the derivative at zero
         } else if (!flag_fab(i,j,k).isConnected(-1,0,0)) {
           // Covered cell to the left, go fish right
-          vx = (c0 * dat_arr(i  ,j,k,1)
-            + c1 * dat_arr(i+1,j,k,1)
-            + c2 * dat_arr(i+2,j,k,1)) * idx;
-          wx = (c0 * dat_arr(i  ,j,k,2)
-            + c1 * dat_arr(i+1,j,k,2)
-            + c2 * dat_arr(i+2,j,k,2)) * idx;
+          if (flag_fab(i+1,j,k).isConnected( 1,0,0)) {
+            vx = (c0 * dat_arr(i  ,j,k,1)
+              + c1 * dat_arr(i+1,j,k,1)
+              + c2 * dat_arr(i+2,j,k,1)) * idx;
+            wx = (c0 * dat_arr(i  ,j,k,2)
+              + c1 * dat_arr(i+1,j,k,2)
+              + c2 * dat_arr(i+2,j,k,2)) * idx;
+          } else {
+            // Only one uncovered cell that way, drop to linear
+            vx = (dat_arr(i+1,j,k,1) - dat_arr(i  ,j,k,1)) * idx;
+            wx = (dat_arr(i+1,j,k,2) - dat_arr(i  ,j,k,2)) * idx;
+          }
         } else {
           // No covered cells right or left, use standard stencil
           vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
@@ -188,38 +229,68 @@ namespace derive_functions
         }
         // Do the same in y-direction
         if (!flag_fab(i,j,k).isConnected(0, 1,0)) {
-          uy = - (c0 * dat_arr(i,j  ,k,0)
-              + c1 * dat_arr(i,j-1,k,0)
-              + c2 * dat_arr(i,j-2,k,0)) * idy;
-          wy = - (c0 * dat_arr(i,j  ,k,2)
-              + c1 * dat_arr(i,j-1,k,2)
-              + c2 * dat_arr(i,j-2,k,2)) * idy;
+          if (flag_fab(i,j,k).isConnected(0,-1,0)) {
+            if (flag_fab(i,j-1,k).isConnected(0,-1,0)) {
+              uy = - (c0 * dat_arr(i,j  ,k,0)
+                  + c1 * dat_arr(i,j-1,k,0)
+                  + c2 * dat_arr(i,j-2,k,0)) * idy;
+              wy = - (c0 * dat_arr(i,j  ,k,2)
+                  + c1 * dat_arr(i,j-1,k,2)
+                  + c2 * dat_arr(i,j-2,k,2)) * idy;
+            } else {
+              // Only one uncovered cell that way, drop to linear
+              uy = (dat_arr(i,j  ,k,0) - dat_arr(i,j-1,k,0)) * idy;
+              wy = (dat_arr(i,j  ,k,2) - dat_arr(i,j-1,k,2)) * idy;
+            }
+          }
+          // Covered on both sides, leave the derivative at zero
         } else if (!flag_fab(i,j,k).isConnected(0,-1,0)) {
-          uy = (c0 * dat_arr(i,j  ,k,0)
-            + c1 * dat_arr(i,j+1,k,0)
-            + c2 * dat_arr(i,j+2,k,0)) * idy;
-          wy = (c0 * dat_arr(i,j  ,k,2)
-            + c1 * dat_arr(i,j+1,k,2)
-            + c2 * dat_arr(i,j+2,k,2)) * idy;
+          if (flag_fab(i,j+1,k).isConnected(0, 1,0)) {
+            uy = (c0 * dat_arr(i,j  ,k,0)
+              + c1 * dat_arr(i,j+1,k,0)
+              + c2 * dat_arr(i,j+2,k,0)) * idy;
+            wy = (c0 * dat_arr(i,j  ,k,2)
+              + c1 * dat_arr(i,j+1,k,2)
+              + c2 * dat_arr(i,j+2,k,2)) * idy;
+          } else {
+            // Only one uncovered cell that way, drop to linear
+            uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j  ,k,0)) * idy;
+            wy = (dat_arr(i,j+1,k,2) - dat_arr(i,j  ,k,2)) * idy;
+          }
         } else {
           uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy;
           wy = 0.5 * (dat_arr(i,j+1,k,2) - dat_arr(i,j-1,k,2)) * idy;
         }
         // Do the same in z-direction
         if (!flag_fab(i,j,k).isConnected(0,0, 1)) {
-          uz = - (c0 * dat_arr(i,j,k  ,0)
-              + c1 * dat_arr(i,j,k-1,0)
-              + c2 * dat_arr(i,j,k-2,0)) * idz;
-          vz = - (c0 * dat_arr(i,j,k  ,1)
-              + c1 * dat_arr(i,j,k-1,1)
-              + c2 * dat_arr(i,j,k-2,1)) * idz;
+          if (flag_fab(i,j,k).isConnected(0,0,-1)) {
+            if (flag_fab(i,j,k-1).isConnected(0,0,-1)) {
+              uz = - (c0 * dat_arr(i,j,k  ,0)
+                  + c1 * dat_arr(i,j,k-1,0)
+                  + c2 * dat_arr(i,j,k-2,0)) * idz;
+              vz = - (c0 * dat_arr(i,j,k  ,1)
+                  + c1 * dat_arr(i,j,k-1,1)
+                  + c2 * dat_arr(i,j,k-2,1)) * idz;
+            } else {
+              // Only one uncovered cell that way, drop to linear
+              uz = (dat_arr(i,j,k  ,0) - dat_arr(i,j,k-1,0)) * idz;
+              vz = (dat_arr(i,j,k  ,1) - dat_arr(i,j,k-1,1)) * idz;
+            }
+          }
+          // Covered on both sides, leave the derivative at zero
         } else if (!flag_fab(i,j,k).isConnected(0,0,-1)) {
-          uz = (c0 * dat_arr(i,j,k  ,0)
-            + c1 * dat_arr(i,j,k+1,0)
-            + c2 * dat_arr(i,j,k+2,0)) * idz;
-          vz = (c0 * dat_arr(i,j,k  ,1)
-            + c1 * dat_arr(i,j,k+1,1)
-            + c2 * dat_arr(i,j,k+2,1)) * idz;
+          if (flag_fab(i,j,k+1).isConnected(0,0, 1)) {
+            uz = (c0 * dat_arr(i,j,k  ,0)
+              + c1 * dat_arr(i,j,k+1,0)
+              + c2 * dat_arr(i,j,k+2,0)) * idz;
+            vz = (c0 * dat_arr(i,j,k  ,1)
+              + c1 * dat_arr(i,j,k+1,1)
+              + c2 * dat_arr(i,j,k+2,1)) * idz;
+          } else {
+            // Only one uncovered cell that way, drop to linear
+            uz = (dat_arr(i,j,k+1,0) - dat_arr(i,j,k  ,0)) * idz;
+            vz = (dat_arr(i,j,k+1,1) - dat_arr(i,j,k  ,1)) * idz;
+          }
         } else {
           uz = 0.5 * (dat_arr(i,j,k+1,0) - dat_arr(i,j,k-1,0)) * idz;
           vz = 0.5 * (dat_arr(i,j,k+1,1) - dat_arr(i,j,k-1,1)) * idz;
diff --git a/Source/NS_getForce.cpp b/Source/NS_getForce.cpp
index e66bee13..6d8cb027 100644
--- a/Source/NS_getForce.cpp
+++ b/Source/NS_getForce.cpp
@@ -92,10 +92,13 @@ NavierStokesBase::getForce (FArrayBox&       force,
 #endif
 
        // Compute min/max
-       for (int n=0; n<ncomp; n++) {
-           amrex::Print() << "State comp " << scomp+n << " min/max "
-                          << State.min<RunOn::Gpu>(scomp+n) << " / "
-                          << State.max<RunOn::Gpu>(scomp+n) << '\n';
+       // State may not hold all of the state components (scalar callers pass
+       // a velocity-only FAB), so index it by its own layout rather than by
+       // the absolute state component scomp+n.
+       for (int n=0; n<State.nComp(); n++) {
+           amrex::Print() << "State comp " << n << " min/max "
+                          << State.min<RunOn::Gpu>(n) << " / "
+                          << State.max<RunOn::Gpu>(n) << '\n';
        }
        for (int n=auxScomp; n<Aux.nComp(); n++) {
            amrex::Print() << "aux comp " << n << " min/max "
@@ -178,10 +181,14 @@ NavierStokesBase::getForce (FArrayBox&       force,
 
    if (ParallelDescriptor::IOProcessor() && getForceVerbose) {
        // Compute min/max
+       // For scalar-only calls force holds ncomp components starting at
+       // component 0 (see scomp_scal above), so index force by its own layout
+       // while still reporting the absolute state component in the label.
        for (int n=0; n<ncomp; n++) {
+           const int fcomp = ( scomp<AMREX_SPACEDIM ) ? scomp+n : n;
            amrex::Print() << "Force comp " << scomp+n << " min/max "
-                          << force.min<RunOn::Gpu>(scomp+n) << " / "
-                          << force.max<RunOn::Gpu>(scomp+n) << '\n';
+                          << force.min<RunOn::Gpu>(fcomp) << " / "
+                          << force.max<RunOn::Gpu>(fcomp) << '\n';
        }
 
       amrex::Print() << "NavierStokesBase::getForce(): Leaving..."
diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp
index 5096470c..6d71f67b 100644
--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -2774,17 +2774,15 @@ NavierStokesBase::scalar_advection_update (Real dt,
             // Must do FillPatch here instead of MF iterator because we need the
             // boundary values in the old data (especially at inflow)
             //
-            const int index_new_s   = Density;
-            const int index_new_rho = Density;
-            const int index_old_s   = index_new_s   - Density;
-            const int index_old_rho = index_new_rho - Density;
-
             FillPatchIterator S_fpi(*this,S_old,1,prev_time,State_Type,Density,1);
             MultiFab& Smf=S_fpi.get_mf();
 
-            ConservativeScalMinMax(S_new, index_new_s, index_new_rho,
-                                   Smf,   index_old_s, index_old_rho);
-
+            //
+            // Clamp the new density to the min/max of the old-time density in
+            // the neighborhood. (The conservative s/rho form with s = rho is
+            // an identity, i.e. a no-op.)
+            //
+            ConvectiveScalMinMax(S_new, Density, Smf, 0);
         }
 
         ++sComp;
@@ -4277,8 +4275,13 @@ NavierStokesBase::ConservativeScalMinMax ( amrex::MultiFab&       Snew, const in
         amrex::ParallelFor(bx, [=]
         AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
+#ifdef AMREX_USE_EB
+            // Leave covered cells alone; their stencil holds no valid data.
+            if ( vfrac(i,j,k) == 0. ) { return; }
+#endif
+
             Real smn = std::numeric_limits<Real>::max();
-            Real smx = std::numeric_limits<Real>::min();
+            Real smx = std::numeric_limits<Real>::lowest();
 
 #if (AMREX_SPACEDIM==3)
             int ks = -1;
@@ -4335,8 +4338,13 @@ NavierStokesBase::ConvectiveScalMinMax ( amrex::MultiFab&       Snew, const int
         amrex::ParallelFor(bx, [=]
         AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
+#ifdef AMREX_USE_EB
+            // Leave covered cells alone; their stencil holds no valid data.
+            if ( vfrac(i,j,k) == 0. ) { return; }
+#endif
+
             Real smn = std::numeric_limits<Real>::max();
-            Real smx = std::numeric_limits<Real>::min();
+            Real smx = std::numeric_limits<Real>::lowest();
 
 #if (AMREX_SPACEDIM==3)
             int ks = -1;
diff --git a/Source/main.cpp b/Source/main.cpp
index 3967b3ec..bcd5ae6b 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -105,7 +105,8 @@ main (int   argc,
     //
     if (Amr::RegridOnRestart())
     {
-        if (    (amrptr->levelSteps(0) >= max_step ) ||
+        if (    ( (max_step >= 0) &&
+                  (amrptr->levelSteps(0) >= max_step) ) ||
                 ( (stop_time >= 0.0) &&
                   (amrptr->cumTime() >= stop_time)  )    )
         {
diff --git a/Source/prob/prob_init.cpp b/Source/prob/prob_init.cpp
index d7f023df..fe3ffb21 100644
--- a/Source/prob/prob_init.cpp
+++ b/Source/prob/prob_init.cpp
@@ -630,6 +630,14 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx,
 {
   const auto domlo = amrex::lbound(domain);
 
+  //
+  // amrex::Abort() cannot be called from within the device lambda below, so
+  // validate the mean flow direction here.
+  //
+  if ( IC.meanFlowDir < -3 || IC.meanFlowDir > 3 ) {
+    amrex::Abort("\n   init_ConvectedVortex: prob.meanFlowDir must be 0 (no mean flow) or +/-1, +/-2, +/-3\n   in the inputs file.");
+  }
+
   amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
   {
     AMREX_D_TERM(Real x = problo[0] + (i - domlo.x + 0.5)*dx[0];,
@@ -661,13 +669,13 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx,
                       vel(i,j,k,2) = w_vort);
          break;
       case 2 :
-         AMREX_D_TERM(vel(i,j,k,0) = v_vort;,
-                      vel(i,j,k,1) = IC.meanFlowMag + u_vort;,
+         AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
+                      vel(i,j,k,1) = IC.meanFlowMag + v_vort;,
                       vel(i,j,k,2) = w_vort);
          break;
       case -2 :
-         AMREX_D_TERM(vel(i,j,k,0) = v_vort;,
-                      vel(i,j,k,1) = -IC.meanFlowMag + u_vort;,
+         AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
+                      vel(i,j,k,1) = -IC.meanFlowMag + v_vort;,
                       vel(i,j,k,2) = w_vort);
          break;
       case 3 :
@@ -680,7 +688,8 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx,
                       vel(i,j,k,1) = -IC.meanFlowMag + v_vort;,
                       vel(i,j,k,2) = w_vort);
          break;
-      default : // no mean flow, i.e. the vortex alone
+      case 0 :  // no mean flow, i.e. the vortex alone
+      default :
          AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
                       vel(i,j,k,1) = v_vort;,
                       vel(i,j,k,2) = w_vort);
diff --git a/Tutorials/HIT/NS_getForce.cpp b/Tutorials/HIT/NS_getForce.cpp
index b08f466a..b25379a7 100644
--- a/Tutorials/HIT/NS_getForce.cpp
+++ b/Tutorials/HIT/NS_getForce.cpp
@@ -318,9 +318,10 @@ NavierStokesBase::getForce (FArrayBox&       force,
      int ff_jlo = jlo/TurbulentForcing::ff_factor;
      int ff_klo = klo/TurbulentForcing::ff_factor;
 
-     int ff_ihi = (ihi+1)/TurbulentForcing::ff_factor;
-     int ff_jhi = (jhi+1)/TurbulentForcing::ff_factor;
-     int ff_khi = (khi+1)/TurbulentForcing::ff_factor;
+     // +1 so that the (ff_i+1) node used by the interpolation below always exists
+     int ff_ihi = ihi/TurbulentForcing::ff_factor + 1;
+     int ff_jhi = jhi/TurbulentForcing::ff_factor + 1;
+     int ff_khi = khi/TurbulentForcing::ff_factor + 1;
 
      // adjust for ghost cells
      if (ilo < (ff_ilo*TurbulentForcing::ff_factor)) {
@@ -332,20 +333,13 @@ NavierStokesBase::getForce (FArrayBox&       force,
      if (klo < (ff_klo*TurbulentForcing::ff_factor)) {
          ff_klo=ff_klo-1;
      }
-     if (ihi == (ff_ihi*TurbulentForcing::ff_factor)) {
-         ff_ihi=ff_ihi+1;
-     }
-     if (jhi == (ff_jhi*TurbulentForcing::ff_factor)) {
-         ff_jhi=ff_jhi+1;
-     }
-     if (khi == (ff_khi*TurbulentForcing::ff_factor)) {
-         ff_khi=ff_khi+1;
-     }
 
      // allocate coarse force array
      Box ffbx(IntVect(ff_ilo, ff_jlo, ff_klo), IntVect(ff_ihi, ff_jhi, ff_khi));
-     // not sure if want elixir, gpu::sync, or async_arena here...
      FArrayBox ff_force(ffbx,AMREX_SPACEDIM);
+     // The Elixir defers freeing ff_force's device memory until the
+     // asynchronous kernels below that read it have completed. No-op on CPU.
+     Elixir ff_force_i = ff_force.elixir();
      const auto& ffarr = ff_force.array();
 
      // Construct node-based coarse forcing
diff --git a/Tutorials/HIT/TurbulentForcing_def.H b/Tutorials/HIT/TurbulentForcing_def.H
index d9a4325b..429015f0 100644
--- a/Tutorials/HIT/TurbulentForcing_def.H
+++ b/Tutorials/HIT/TurbulentForcing_def.H
@@ -105,13 +105,18 @@ TurbulentForcing::init_turbulent_forcing (const amrex::GpuArray<amrex::Real,AMRE
     int nymodes = nmodes*(int)(0.5+Ly/Lmin);
     int nzmodes = nmodes*(int)(0.5+Lz/Lmin);
 
+    // Mode indices run 0:n?modes, so they must fit in the 0:array_size-1
+    // extent of the forcing arrays.
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::max(nxmodes,std::max(nymodes,nzmodes)) < array_size,
+                                     "TurbulentForcing: too many modes; reduce turb.nmodes or increase array_size");
+
     if (verbose)
     {
         Print() << "Lmin = " << Lmin << '\n';
         Print() << "kappaMax = " << kappaMax << '\n';
         Print() << "nxmodes = " << nxmodes << '\n';
-        Print() << "nymodes = " << nxmodes << '\n';
-        Print() << "nzmodes = " << nxmodes << '\n';
+        Print() << "nymodes = " << nymodes << '\n';
+        Print() << "nzmodes = " << nzmodes << '\n';
     }
 
     Real freqMin = 1.0/forcing_time_scale_max;
@@ -355,7 +360,9 @@ TurbulentForcing::init_turbulent_forcing (const amrex::GpuArray<amrex::Real,AMRE
     if (Gpu::inLaunchRegion())
     {
         forcedata = static_cast<Real*>(The_Arena()->alloc(tmp_size*sizeof(Real)));
-        Gpu::htod_memcpy_async(forcedata, tmp, tmp_size*sizeof(Real));
+        // Synchronous copy: tmp is a stack array that dies when this
+        // function returns, so the copy must complete before then.
+        Gpu::htod_memcpy(forcedata, tmp, tmp_size*sizeof(Real));
     }
     else
 #endif
diff --git a/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp b/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp
index 847afc4a..8fa08944 100644
--- a/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp
+++ b/Tutorials/TaylorGreen/benchmarks/ViscBench.cpp
@@ -145,7 +145,7 @@ main (int   argc,
     Vector<MultiFab*> error(finestLevel+1);
     Vector<MultiFab*> dataE(finestLevel+1);
 
-    std::cout << "Level Delta L"<< norm << " norm of Error in Each Component" << '\n'
+    amrex::Print() << "Level Delta L"<< norm << " norm of Error in Each Component" << '\n'
          << "-----------------------------------------------" << '\n';
 
     for (int iLevel = 0; iLevel <= finestLevel; ++iLevel)
@@ -162,18 +162,20 @@ main (int   argc,
 
         MultiFab dataI(baI, dmap, nComp, 0);
 
-        for (int iGrid=0; iGrid<baI.size(); ++iGrid)
+        //
+        // Fill the data from the plotfile. Use the distributed MultiFab
+        // overload of FillVar so that each rank fills only the boxes it owns;
+        // the FArrayBox overload fills on a single rank only.
+        //
+        for (int iComp=0; iComp<nComp; ++iComp)
         {
-            const Box& dataGrid = baI[iGrid];
-
-            for (int iComp=0; iComp<nComp; ++iComp)
-            {
-                FArrayBox tmpFab(dataGrid,1);
+            amrDataI.FillVar(dataI, iLevel, derives[iComp], iComp);
+        }
 
-                amrDataI.FillVar(&tmpFab, dataGrid,
-                                 iLevel, derives[iComp], 0);
-                dataI[iGrid].copy(tmpFab,0,iComp,1);
-            }
+        for (MFIter mfi(*dataE[iLevel]); mfi.isValid(); ++mfi)
+        {
+            const int iGrid = mfi.index();
+            const Box& dataGrid = mfi.validbox();
 
             //
             // Fill exact solution
@@ -186,14 +188,14 @@ main (int   argc,
 #if (AMREX_SPACEDIM == 3)
             FORT_VISCBENCH(&time, &mu, &unifDir,
                            lo, hi, &nComp,
-                           ((*dataE[iLevel])[iGrid]).dataPtr(),
+                           ((*dataE[iLevel])[mfi]).dataPtr(),
                            AMREX_ARLIM(lo), AMREX_ARLIM(hi),
                            delI.dataPtr(),
                            xlo.dataPtr(), xhi.dataPtr());
 #else
             FORT_VISCBENCH(&time, &mu,
                            lo, hi, &nComp,
-                           ((*dataE[iLevel])[iGrid]).dataPtr(),
+                           ((*dataE[iLevel])[mfi]).dataPtr(),
                            AMREX_ARLIM(lo), AMREX_ARLIM(hi),
                            delI.dataPtr(),
                            xlo.dataPtr(), xhi.dataPtr());
@@ -213,7 +215,7 @@ main (int   argc,
 
         Real delAvg = pow(cellvol, (1.0/BL_SPACEDIM));
 
-        std::cout << "  " << iLevel << " " << delAvg << "    ";
+        amrex::Print() << "  " << iLevel << " " << delAvg << "    ";
         Real Ln = 0.0;
 
         for (int iComp=0; iComp<nComp; ++iComp)
@@ -230,9 +232,9 @@ main (int   argc,
                 // L_2 norm
                 Ln = (*error[iLevel]).norm2(iComp)*sqrt(cellvol);
             }
-            std::cout << Ln << "  ";
+            amrex::Print() << Ln << "  ";
         }
-        std::cout << '\n';
+        amrex::Print() << '\n';
     }
 
     //