diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp
index a40bb292..d05f68da 100644
--- a/src/convection/incflo_compute_MAC_projected_velocities.cpp
+++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp
@@ -208,6 +208,12 @@ incflo::compute_MAC_projected_velocities (
             velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
         }
 
+        // On a direction_dependent face IncfloVelFill writes the prescribed velocity
+        // into the ghost cell only where it points into the domain; on the outflow
+        // part it copies the interior cell, which is 0 here because
+        // time_dep_inflow_vel was zeroed above.  The inflow tests below must
+        // therefore be strict, otherwise the outflow part of the face is treated as
+        // zero inflow and the extrapolated outflow velocity is lost.
         for (MFIter mfi(time_dep_inflow_vel,false); mfi.isValid(); ++mfi)
         {
             Box const& bx = mfi.validbox();
@@ -222,11 +228,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 0;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (i == dlo.x && ( bc.lo(0) == BCType::ext_dir ||
-                                   (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) >= Real(0.0)) ) ) {
+                                   (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) >  Real(0.0)) ) ) {
                     umac_arr(i,j,k) = cc_arr(i-1,j,k,0);
                 }
                 if (i == dhi.x && ( bc.hi(0) == BCType::ext_dir ||
-                                   (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) <= Real(0.0)) ) ) {
+                                   (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) <  Real(0.0)) ) ) {
                     umac_arr(i+1,j,k) = cc_arr(i+1,j,k,0);
                 }
             });
@@ -236,11 +242,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 1;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (j == dlo.y && ( bc.lo(1) == BCType::ext_dir||
-                                   (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) >= Real(0.0)) ) ) {
+                                   (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) >  Real(0.0)) ) ) {
                     vmac_arr(i,j,k) = cc_arr(i,j-1,k,1);
                 }
                 if (j == dhi.y && ( bc.hi(1) == BCType::ext_dir ||
-                                   (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) <= Real(0.0)) ) ) {
+                                   (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) <  Real(0.0)) ) ) {
                     vmac_arr(i,j+1,k) = cc_arr(i,j+1,k,1);
                 }
             });
@@ -251,11 +257,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 2;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (k == dlo.z && ( bc.lo(2) == BCType::ext_dir ||
-                                   (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) >= Real(0.0)) ) ) {
+                                   (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) >  Real(0.0)) ) ) {
                     wmac_arr(i,j,k) = cc_arr(i,j,k-1,2);
                 }
                 if (k == dhi.z && ( bc.hi(2) == BCType::ext_dir ||
-                                   (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) <= Real(0.0)) ) ) {
+                                   (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) <  Real(0.0)) ) ) {
                     wmac_arr(i,j,k+1) = cc_arr(i,j,k+1,2);
                 }
             });
diff --git a/src/incflo.cpp b/src/incflo.cpp
index 7a29256d..ed47f4e5 100644
--- a/src/incflo.cpp
+++ b/src/incflo.cpp
@@ -170,6 +170,13 @@ void incflo::Evolve()
         {
             if (m_verbose > 0) amrex::Print() << "Regridding...\n";
             regrid(0, m_cur_time);
+#ifdef INCFLO_USE_PARTICLES
+            // This must be done after regrid() returns: inside RemakeLevel /
+            // MakeNewLevelFromCoarse the ParGDB still sees the old BoxArray,
+            // DistributionMapping and finest_level, so a Redistribute there would
+            // map the particles onto the pre-regrid layout.
+            particleData.Redistribute();
+#endif
             if (m_verbose > 0 && ParallelDescriptor::IOProcessor()) {
                 printGridSummary(amrex::OutStream(), 0, finest_level);
             }
diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp
index 114e3a7e..5f1d0847 100644
--- a/src/incflo_apply_predictor.cpp
+++ b/src/incflo_apply_predictor.cpp
@@ -229,7 +229,10 @@ void incflo::ApplyPredictor (bool incremental_projection)
     // **************************************************************************************
     // Update the particle positions
     // **************************************************************************************
-    if (m_advection_type != "MOL") {
+    // Not during the initial pressure iterations: the fields are reset after each
+    // of those, but the particle positions would not be, so the particles would
+    // start the run m_initial_iterations*dt ahead of the fluid.
+    if (m_advection_type != "MOL" && !incremental_projection) {
         evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac),
                                            GetVecOfConstPtrs(w_mac)));
     }
diff --git a/src/incflo_compute_dt.cpp b/src/incflo_compute_dt.cpp
index b4b5f1ca..f96269ed 100644
--- a/src/incflo_compute_dt.cpp
+++ b/src/incflo_compute_dt.cpp
@@ -52,6 +52,51 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
 
        compute_vel_forces_on_level (lev, vel_forces, vel, rho, tra_o, tra);
 
+       // Explicit-diffusion bound: the largest diffusivity that is applied
+       // explicitly, divided by rho.  We must use the strain-rate dependent
+       // viscosity here (not just the constant m_mu) and we must cover the tracer
+       // and temperature diffusivities as well, since those are updated explicitly
+       // with the same m_diff_type switch.  Covered cells are set to zero.
+       MultiFab nu;
+       if (explicit_diffusion) {
+           nu.define(grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev));
+           compute_viscosity_at_level(lev, &nu, &m_leveldata[lev]->density,
+                                      &m_leveldata[lev]->velocity, geom[lev],
+                                      m_cur_time, 0);
+           Real mu_s_max = Real(0.0);
+           if (m_advect_tracer) {
+               for (int n = 0; n < m_ntrac; ++n) {
+                   mu_s_max = amrex::max(mu_s_max, m_mu_s[n]);
+               }
+           }
+           Real mu_T_eff = m_use_temperature ? m_mu_T / m_cp : Real(0.0);
+#ifdef AMREX_USE_EB
+           auto const& dt_flags = EBFactory(lev).getMultiEBCellFlagFab();
+#endif
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+           for (MFIter mfi(nu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+               Box const& bx = mfi.tilebox();
+               Array4<Real> const& nu_a = nu.array(mfi);
+               Array4<Real const> const& r = rho.const_array(mfi);
+#ifdef AMREX_USE_EB
+               Array4<EBCellFlag const> const& f = dt_flags.const_array(mfi);
+#endif
+               ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+               {
+#ifdef AMREX_USE_EB
+                   if (f(i,j,k).isCovered()) { nu_a(i,j,k) = Real(0.0); return; }
+#endif
+                   Real rinv = Real(1.0)/r(i,j,k);
+                   // Velocity and temperature diffuse with eta/rho and mu_T/(rho cp);
+                   // a non-conservative tracer diffuses with mu_s itself.
+                   nu_a(i,j,k) = amrex::max(amrex::max(nu_a(i,j,k), mu_T_eff)*rinv,
+                                            mu_s_max*amrex::max(Real(1.0), rinv));
+               });
+           }
+       }
+
 #ifdef AMREX_USE_EB
         if (!vel.isAllRegular()) {
             auto const& flag = EBFactory(lev).getMultiEBCellFlagFab();
@@ -72,22 +117,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
                            return mx;
                        });
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, flag, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r,
-                                                      Array4<EBCellFlag const> const& f) -> Real
-                          {
-                              Real mx = Real(-1.0);
-                              amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                              {
-                                  if (!f(i,j,k).isCovered()) {
-                                      Real rho_inv = Real(1.0)/r(i,j,k);
-                                      mx = amrex::max(rho_inv, mx);
-                                  }
-                              });
-                              return mx;
-                          });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing
@@ -131,19 +161,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
                        });
 
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r) -> Real
-                           {
-                               Real mx = Real(-1.0);
-                               amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                               {
-                                   Real rho_inv = Real(1.0)/r(i,j,k);
-                                   mx = amrex::max(rho_inv, mx);
-                               });
-                               return mx;
-                           });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing
diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp
index 11ef44c7..37adfc7d 100644
--- a/src/incflo_regrid.cpp
+++ b/src/incflo_regrid.cpp
@@ -59,10 +59,6 @@ void incflo::MakeNewLevelFromCoarse (int lev,
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,lev));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Remake an existing level using provided BoxArray and DistributionMapping and
@@ -120,10 +116,6 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,finest_level));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Rebuild macproj on demand.  Note that this must not be done inside ClearLevel:
diff --git a/src/prob/prob_bc.H b/src/prob/prob_bc.H
index 29382e68..d85af581 100644
--- a/src/prob/prob_bc.H
+++ b/src/prob/prob_bc.H
@@ -5,6 +5,95 @@
 #include <AMReX_Geometry.H>
 #include <AMReX_PhysBCFunct.H>
 
+// Normal velocity prescribed on a domain face, including the probtype-specific
+// inflow profiles.  Every fill functor uses this so that the direction_dependent
+// inflow/outflow decision is the same for velocity, density, tracer and
+// temperature.  (Probtype 16 is a moving wall, not an inflow, and so is handled
+// in IncfloVelFill directly.)
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real incflo_bc_normal_velocity (int probtype, amrex::Orientation face,
+                                       int i, int j, int k,
+                                       amrex::Box const& domain_box, amrex::Real time,
+                                       amrex::GpuArray<amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>, AMREX_SPACEDIM*2> const& bcv_vel)
+{
+    amrex::ignore_unused(i,j,k,time);
+    const int dir = face.coordDir();
+    amrex::Real norm_vel = bcv_vel[face][dir];
+
+    if (dir == 0 && face.isLow())
+    {
+        if (42 == probtype)
+        {
+            norm_vel = time;
+        }
+        else if (31 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1));
+            norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y);
+        }
+        else if (43 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1);
+        }
+#if (AMREX_SPACEDIM == 3)
+        else if (311 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
+            norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z);
+        }
+        else if (41 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
+            norm_vel = amrex::Real(0.5) * z;
+        }
+#endif
+    }
+    else if (dir == 0 && face.isHigh())
+    {
+        if (42 == probtype)
+        {
+            norm_vel = time;
+        }
+        else if (43 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1);
+        }
+    }
+    else if (dir == 1 && face.isLow())
+    {
+#if (AMREX_SPACEDIM == 3)
+        if (32 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2));
+            norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z);
+        }
+#endif
+        if (322 == probtype)
+        {
+            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+            norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x);
+        }
+    }
+#if (AMREX_SPACEDIM == 3)
+    else if (dir == 2 && face.isLow())
+    {
+        if (33 == probtype)
+        {
+            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+            norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x);
+        }
+        else if (333 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y);
+        }
+    }
+#endif
+    return norm_vel;
+}
+
 struct IncfloVelFill
 {
     int probtype;
@@ -44,35 +133,9 @@ struct IncfloVelFill
             if (i < domain_box.smallEnd(0))
             {
                 int dir = 0;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (42 == probtype)
-                {
-                    norm_vel = time;
-                }
-                else if (31 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1));
-                    norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y);
-                }
-                else if (43 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1);
-                }
-#if (AMREX_SPACEDIM == 3)
-                else if (311 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
-                    norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z);
-                }
-                else if (41 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
-                    norm_vel = amrex::Real(0.5) * z;
-                }
-#endif
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::x,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1101 == probtype)
@@ -103,18 +166,9 @@ struct IncfloVelFill
             if (i > domain_box.bigEnd(0))
             {
                 int dir = 0;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (42 == probtype)
-                {
-                    norm_vel = time;
-                }
-                else if (43 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::x,amrex::Orientation::high),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1101 == probtype)
@@ -145,21 +199,9 @@ struct IncfloVelFill
             if (j < domain_box.smallEnd(1))
             {
                 int dir = 1;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-#if (AMREX_SPACEDIM == 3)
-                if (32 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2));
-                    norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z);
-                }
-#endif
-                if (322 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::y,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 if ( (bc.lo(dir) == amrex::BCType::ext_dir) ||
                      (bc.lo(dir) == amrex::BCType::direction_dependent && norm_vel >= 0.) )
@@ -182,14 +224,9 @@ struct IncfloVelFill
             if (j > domain_box.bigEnd(1))
             {
                 int dir = 1;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (16 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::y,amrex::Orientation::high),
+                    i, j, k, domain_box, time, bcv_vel);
 
 #if (AMREX_SPACEDIM == 3)
                 if (1102 == probtype)
@@ -207,6 +244,12 @@ struct IncfloVelFill
                     {
                         if ( orig_comp+nc == dir ) {
                             vel(i,j,k,dcomp+nc) = norm_vel;
+                        } else if (16 == probtype && orig_comp+nc == 0) {
+                            // Burggraf lid: the prescribed profile is the *tangential*
+                            // velocity u = 16 x^2 (1-x)^2 on the high-y wall; the normal
+                            // velocity keeps the wall value.
+                            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+                            vel(i,j,k,dcomp+nc) = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x);
                         } else {
                             vel(i,j,k,dcomp+nc) =  bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][orig_comp+nc];
                         }
@@ -225,19 +268,9 @@ struct IncfloVelFill
             if (k < domain_box.smallEnd(2))
             {
                 int dir = 2;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (33 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x);
-                }
-                else if (333 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::z,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1100 == probtype)
@@ -270,7 +303,9 @@ struct IncfloVelFill
             if (k > domain_box.bigEnd(2))
             {
                 int dir = 2;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][dir];
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::z,amrex::Orientation::high),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1100 == probtype)
@@ -314,7 +349,7 @@ struct IncfloDenFill
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& rho,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -381,48 +416,75 @@ struct IncfloDenFill
         else if ( (i < domain_box.smallEnd(0)) &&
                   ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                     (bc.lo(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
         }
         else if ( (i > domain_box.bigEnd(0)) &&
                   ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                     (bc.hi(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)];
         }
+        // Outflow part of a direction_dependent face: extrapolate the interior cell.
+        // Nobody else fills these ghost cells -- filcc has no case for
+        // BCType::direction_dependent -- so without this they keep stale data.
+        else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(domain_box.smallEnd(0),j,k);
+        }
+        else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(domain_box.bigEnd(0),j,k);
+        }
 
         if ( (j < domain_box.smallEnd(1)) &&
               ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                 (bc.lo(1) == amrex::BCType::direction_dependent &&
-                 bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                 incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)];
         }
         else if ( (j > domain_box.bigEnd(1)) &&
                   ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                     (bc.hi(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)];
         }
+        else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,domain_box.smallEnd(1),k);
+        }
+        else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,domain_box.bigEnd(1),k);
+        }
 #if (AMREX_SPACEDIM == 3)
          if ( (k < domain_box.smallEnd(2)) &&
               ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                 (bc.lo(2) == amrex::BCType::direction_dependent &&
-                 bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                 incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)];
         }
+        else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,j,domain_box.smallEnd(2));
+        }
 
         if ( (k > domain_box.bigEnd(2)) &&
                   ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                     (bc.hi(2) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)];
         }
+        else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,j,domain_box.bigEnd(2));
+        }
 #endif
     }
 };
@@ -442,7 +504,7 @@ struct IncfloTracFill
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& tracer,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -511,47 +573,74 @@ struct IncfloTracFill
             else if ( (i < domain_box.smallEnd(0)) &&
                       ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                         (bc.lo(0) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][n];
             }
             else if ( (i > domain_box.bigEnd(0)) &&
                       ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                         (bc.hi(0) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][n];
             }
+            // Outflow part of a direction_dependent face: extrapolate the interior cell.
+            // Nobody else fills these ghost cells -- filcc has no case for
+            // BCType::direction_dependent -- so without this they keep stale data.
+            else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(domain_box.smallEnd(0),j,k,n);
+            }
+            else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(domain_box.bigEnd(0),j,k,n);
+            }
 
             if ( (j < domain_box.smallEnd(1)) &&
                   ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                     (bc.lo(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][n];
             }
             else if ( (j > domain_box.bigEnd(1)) &&
                       ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                         (bc.hi(1) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][n];
             }
+            else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,domain_box.smallEnd(1),k,n);
+            }
+            else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,domain_box.bigEnd(1),k,n);
+            }
 #if (AMREX_SPACEDIM == 3)
             if ( (k < domain_box.smallEnd(2)) &&
                  ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                    (bc.lo(2) == amrex::BCType::direction_dependent &&
-                    bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                    incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][n];
             }
             else if ( (k > domain_box.bigEnd(2)) &&
                       ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                         (bc.hi(2) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][n];
             }
+            else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,j,domain_box.smallEnd(2),n);
+            }
+            else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,j,domain_box.bigEnd(2),n);
+            }
 #endif
         }
     }
@@ -572,7 +661,7 @@ struct IncfloTempFill
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& temperature,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -639,47 +728,74 @@ struct IncfloTempFill
         else if ( (i < domain_box.smallEnd(0)) &&
                   ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                     (bc.lo(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
         }
         else if ( (i > domain_box.bigEnd(0)) &&
                   ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                     (bc.hi(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)];
         }
+        // Outflow part of a direction_dependent face: extrapolate the interior cell.
+        // Nobody else fills these ghost cells -- filcc has no case for
+        // BCType::direction_dependent -- so without this they keep stale data.
+        else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(domain_box.smallEnd(0),j,k);
+        }
+        else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(domain_box.bigEnd(0),j,k);
+        }
 
         if ( (j < domain_box.smallEnd(1)) &&
              ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                (bc.lo(1) == amrex::BCType::direction_dependent &&
-                bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)];
         }
         else if ( (j > domain_box.bigEnd(1)) &&
                   ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                     (bc.hi(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)];
         }
+        else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,domain_box.smallEnd(1),k);
+        }
+        else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,domain_box.bigEnd(1),k);
+        }
 #if (AMREX_SPACEDIM == 3)
         if ( (k < domain_box.smallEnd(2)) &&
              ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                (bc.lo(2) == amrex::BCType::direction_dependent &&
-                bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)];
         }
         else if ( (k > domain_box.bigEnd(2)) &&
                   ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                     (bc.hi(2) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)];
         }
+        else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,j,domain_box.smallEnd(2));
+        }
+        else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,j,domain_box.bigEnd(2));
+        }
 #endif
     }
 };
diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp
index 5762b87c..b2aef967 100644
--- a/src/prob/prob_init_fluid.cpp
+++ b/src/prob/prob_init_fluid.cpp
@@ -1035,7 +1035,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
             if (nt > 2 && i <= dhi.x*3/4) tracer(i,j,k,2) = 3.0;
         });
     }
-    else if (42 == m_probtype || 43 == m_probtype)
+    else if (42 == m_probtype)
     {
         ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
@@ -1046,6 +1046,25 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
             }
         });
     }
+    else if (43 == m_probtype)
+    {
+        // Start from the zero-net-flux in/out profile that IncfloVelFill imposes on
+        // the x faces (prob_bc.H), so that the outflow part of the direction_dependent
+        // boundary has a non-zero velocity when the initial projection enforces
+        // in/out solvability.
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            Real y = Real(j+0.5)*dyinv;
+            AMREX_D_TERM(vel(i,j,k,0) = Real(6.0) * y * (Real(1.0)-y) - Real(1.0);,
+                         vel(i,j,k,1) = Real(0.0);,
+                         vel(i,j,k,2) = Real(0.0););
+
+            const int nt = tracer.nComp();
+            for (int n = 0; n < nt; ++n) {
+                tracer(i,j,k,n) = 0.0;
+            }
+        });
+    }
     else if (32 == m_probtype)
     {
         Real v = m_ic_v;