diff --git a/src/incflo_advance.cpp b/src/incflo_advance.cpp
index 431ba8ca..ecfd858c 100644
--- a/src/incflo_advance.cpp
+++ b/src/incflo_advance.cpp
@@ -108,15 +108,15 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol)
 
     amrex::Print() << "Pred:Change in domain volume  " << (vol - prev_vol) << std::endl;
     // For plane geoms:
-    // amrex::Print() << "Pred:Change in domain mass   " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
-    // amrex::Print() << "Pred:Change speed   " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
-    // amrex::Print() << "Pred:Change error   " << (sum - prev_mass - vol + prev_vol)/ (vol-prev_vol) << std::endl;
+    amrex::Print() << "Pred:Change in domain mass   " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
+    amrex::Print() << "Pred:Change speed   " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
+    amrex::Print() << "Pred:Change error   " << (sum - prev_mass - vol + prev_vol)/ (vol-prev_vol) << std::endl;
     
     // amrex::Print() << "Pred:Change over time in last time step divided by time and area " << (sum - prev_mass) / m_dt / 1.6 <<
     //                   " using " << sum << " " << prev_mass << " " << m_dt << std::endl;
 
     // For cyl :
-    amrex::Print() << "Pred:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl;
+//    amrex::Print() << "Pred:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl;
     // amrex::Print() << "Pred:Change speed   " << (sum - prev_mass - (vol-prev_vol)*(2.45+1.625)/2.) / ((2.45-1.)* dx[1]*64 *m_dt) << std::endl;
 
     // For sphere :
@@ -157,30 +157,31 @@ void incflo::Advance(Real orig_mass, Real& prev_mass, Real& prev_vol)
     amrex::Print() << "Corr:Change in domain volume " << (vol - prev_vol) <<std::endl;
     static Real verr = 0.;
     verr += (vol-prev_vol);
-    amrex::Print() << "Corr:Sum volume error " << verr << std::endl;
+    amrex::Print() << "Corr:Sum (volume error) " << verr << std::endl;
     //amrex::Print() << "Corr:Change in domain mass   " << (sum - prev_mass) << std::endl;
-    // amrex::Print() << "Corr:Change in domain mass   " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
-    // amrex::Print() << "Corr:Change speed   " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
-
-
-    amrex::Print() << "Corr:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl;
-    // This only works if we start at time 0
-    // Real m_calc = (orig_mass + (2.45-1.)*dx[1]*ncelly*m_t_new[0]*1.);
-    // amrex::Print() << "Corr:Mass conservation error 2D " << (sum - m_calc) << std::endl;
-    Real m_calc = (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.);
-    static Real err = 0.;
-    err += (sum - m_calc);
-    amrex::Print() << "Corr:Sum conservation error " << err << std::endl;
-    //if ( std::abs( (vol - prev_vol) -  (sum - prev_mass) ) > 1e-14 ) Abort("Vol diff and dens diff are off");
-    // amrex::Print() << "Corr:Sum of mass at time = " << m_cur_time+m_dt << " " << sum << " " << std::endl;
-//  amrex::Print() << "Change over time divided by time and area " << (sum - orig_mass) / (m_cur_time+m_dt) / 1.6 <<
-//                    " using " << sum << " " << orig_mass << " " << m_cur_time+m_dt << std::endl;
-    // amrex::Print() << "Corr:Change over time in last time step divided by time and area " << (sum - prev_mass) / m_dt / 1.6 <<
-    //                   " using " << sum << " " << prev_mass << " " << m_dt << std::endl;
+    amrex::Print() << "Corr:Change in domain mass   " << (sum - (prev_mass - 1.*dx[1]*16*m_dt*2.)) << std::endl;
+    amrex::Print() << "Corr:Change speed   " << (sum - prev_mass) / (1.* dx[1]*16 *m_dt) << std::endl;
+
+
+//     amrex::Print() << "Corr:Change in domain mass 2D " << (sum - (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.)) << std::endl;
+//     // This only works if we start at time 0
+//     // Real m_calc = (orig_mass + (2.45-1.)*dx[1]*ncelly*m_t_new[0]*1.);
+//     // amrex::Print() << "Corr:Mass conservation error 2D " << (sum - m_calc) << std::endl;
+//     Real m_calc = (prev_mass + (2.45-1.)*dx[1]*ncelly*m_dt*1.);
+//     static Real err = 0.;
+//     err += (sum - m_calc);
+//     amrex::Print() << "Corr:Sum conservation error " << err << std::endl;
+//     //if ( std::abs( (vol - prev_vol) -  (sum - prev_mass) ) > 1e-14 ) Abort("Vol diff and dens diff are off");
+//     // amrex::Print() << "Corr:Sum of mass at time = " << m_cur_time+m_dt << " " << sum << " " << std::endl;
+// //  amrex::Print() << "Change over time divided by time and area " << (sum - orig_mass) / (m_cur_time+m_dt) / 1.6 <<
+// //                    " using " << sum << " " << orig_mass << " " << m_cur_time+m_dt << std::endl;
+//     // amrex::Print() << "Corr:Change over time in last time step divided by time and area " << (sum - prev_mass) / m_dt / 1.6 <<
+//     //                   " using " << sum << " " << prev_mass << " " << m_dt << std::endl;
+
+#endif
 
     prev_mass = sum;
     prev_vol = vol;
-#endif
 
     // Stop timing current time step
     Real end_step = static_cast<Real>(ParallelDescriptor::second()) - strt_step;
diff --git a/src/redistribution/hydro_redistribution.cpp b/src/redistribution/hydro_redistribution.cpp
index 4d28ea17..c9d725e5 100644
--- a/src/redistribution/hydro_redistribution.cpp
+++ b/src/redistribution/hydro_redistribution.cpp
@@ -323,39 +323,40 @@ void Redistribution::Apply ( Box const& bx, int ncomp,
 
             Real eps = 1.e-14;
 
-            // Initialize Ubar ...
-            // FIXME? perhaps this fits better in state_utils
-            amrex::ParallelFor(Box(scratch),
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-            {
-                // Only include cells that have EB
-                if (vfrac_old(i,j,k) < 1.0 || vfrac_new(i,j,k) < 1.0)
-                {
-                    Vbar(i,j,k) = vfrac_old(i,j,k);
-                    //Vbar(i,j,k) = alpha(i,j,k,0)*vfrac_old(i,j,k);
-                    // Real Vbar = Real(1.0);
-                    if (vel_eb_new){
-                        Vbar(i,j,k) += vfrac_new(i,j,k);
-                    }
-
-                    for ( int n = 0; n < ncomp; n++){
-                        ubar(i,j,k,n) = vfrac_old(i,j,k)*U_in(i,j,k,n);
-                        //ubar(i,j,k,n) = alpha(i,j,k,0)*vfrac_old(i,j,k)*U_in(i,j,k,n);
-                        // ubar(i,j,k,n) = U_in(i,j,k,n);
-
-                        if (vel_eb_new){
-                            ubar(i,j,k,n) += vfrac_new(i,j,k)*out(i,j,k,n);
-                        }
-                    }
-                }
-                else
-                {
-                    Vbar(i,j,k) = 0.;
-                    for ( int n = 0; n < ncomp; n++){
-                        ubar(i,j,k,n) = 0.;
-                    }
-                }
-            });
+// Now try Ubar = Qhat for U_in (time n)
+            // // Initialize Ubar ...
+            // // FIXME? perhaps this fits better in state_utils
+            // amrex::ParallelFor(Box(scratch),
+            // [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            // {
+            //     // Only include cells that have EB
+            //     if (vfrac_old(i,j,k) < 1.0 || vfrac_new(i,j,k) < 1.0)
+            //     {
+            //         Vbar(i,j,k) = vfrac_old(i,j,k);
+            //         //Vbar(i,j,k) = alpha(i,j,k,0)*vfrac_old(i,j,k);
+            //         // Real Vbar = Real(1.0);
+            //         if (vel_eb_new){
+            //             Vbar(i,j,k) += vfrac_new(i,j,k);
+            //         }
+
+            //         for ( int n = 0; n < ncomp; n++){
+            //             ubar(i,j,k,n) = vfrac_old(i,j,k)*U_in(i,j,k,n);
+            //             //ubar(i,j,k,n) = alpha(i,j,k,0)*vfrac_old(i,j,k)*U_in(i,j,k,n);
+            //             // ubar(i,j,k,n) = U_in(i,j,k,n);
+
+            //             if (vel_eb_new){
+            //                 ubar(i,j,k,n) += vfrac_new(i,j,k)*out(i,j,k,n);
+            //             }
+            //         }
+            //     }
+            //     else
+            //     {
+            //         Vbar(i,j,k) = 0.;
+            //         for ( int n = 0; n < ncomp; n++){
+            //             ubar(i,j,k,n) = 0.;
+            //         }
+            //     }
+            // });
 
             amrex::ParallelFor(Box(scratch), ncomp,
             [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
@@ -387,6 +388,26 @@ void Redistribution::Apply ( Box const& bx, int ncomp,
 
                         delta_divU = (delta_vol - Ueb_dot_an); //* U_in(i,j,k,n);
                     }
+                        if ( i==15 && j==12){
+                        Print()<<"NC DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==16 && j==12){
+                        Print()<<"NB DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==15 && j==13){
+                        Print()<<"13 C DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==16 && j==13){
+                        Print()<<"13 R DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<delta_vol<<std::endl;
+                    }
 
                     // FIXME - need to think through NU case here. How to handle NU advective flux
                     // and kappa, as both were just added to NB before...
@@ -413,21 +434,37 @@ void Redistribution::Apply ( Box const& bx, int ncomp,
                             // + out(i,j,k,n) * (delta_vol - Ueb_dot_an_new));
                         }
 
-                    // if ( i==9 && j==8){
-                    //     Print()<<"NU DELTA DIVU "<<delta_divU
-                    //            <<" "<<Ueb_dot_an
-                    //            <<" "<<Ueb_dot_an_new
-                    //            <<" "<<delta_vol<<std::endl;
-                    // }
-                    // if ( i==10 && j==8){
-                    //     Print()<<"NUN DELTA DIVU "<<delta_divU
-                    //            <<" "<<Ueb_dot_an
-                    //            <<" "<<Ueb_dot_an_new
-                    //            <<" "<<delta_vol<<std::endl;
-                    //     Print()<<"U "<<U_in(i,j,k,n)
-                    //            <<" "<<out(i,j,k,n)
-                    //            <<std::endl;
-                    // }
+                        if ( i==15 && j==12){
+                        Print()<<"NC DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<Ueb_dot_an_new
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==16 && j==12){
+                        Print()<<"NB DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<Ueb_dot_an_new
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==15 && j==13){
+                        Print()<<"13 C DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<Ueb_dot_an_new
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    if ( i==16 && j==13){
+                        Print()<<"13 R DELTA DIVU "<<delta_divU
+                               <<" "<<Ueb_dot_an
+                               <<" "<<Ueb_dot_an_new
+                               <<" "<<delta_vol<<std::endl;
+                    }
+                    }
+
+                    if ( i==15 && j==12){
+                        Print()<<"covered U "<<Box(scratch)
+                               <<" U_old, hat "<<U_in(i,j,k,n)
+                               <<" "<<U_in(i,j,k,n)+dt*dUdt_in(i,j,k,n)
+                               <<std::endl;
                     }
 
                     // This will undo volume scaling that happens later in forming q-hat
@@ -445,101 +482,103 @@ void Redistribution::Apply ( Box const& bx, int ncomp,
                         scratch(i,j,k,n) = U_in(i,j,k,n) + dt * dUdt_in(i,j,k,n);
                     }
                     kappa(i,j,k) = dt * delta_divU;
-
-                    // Now add to Ubar if cell has EB
-                    if (vfrac_old(i,j,k) < 1.0 || vfrac_new(i,j,k) < 1.0)
-                    {
-                        for (int i_nbor = 1; i_nbor <= itr(i,j,k,0); i_nbor++)
-                        {
-                            // FIXME - check to make sure the box contains rst?
-                            int r = i+map[0][itr(i,j,k,i_nbor)];
-                            int s = j+map[1][itr(i,j,k,i_nbor)];
-                            int t = k+((AMREX_SPACEDIM < 3) ? 0 : map[2][itr(i,j,k,i_nbor)]);
-
-                            // Add to my Ubar
-                            // Now add me to my nbs Ubar
-                            // FIXME - does this double weight if nbhds reciprocal?
-                            if ( n == 0 ){
-                                Vbar(i,j,k) += vfrac_old(r,s,t);
-                                Vbar(r,s,t) += vfrac_old(i,j,k);
-                                // Vbar(i,j,k) += alpha(i,j,k,1)*vfrac_old(r,s,t)/nrs(r,s,t);
-                                if (vel_eb_new){
-                                    Vbar(i,j,k) += vfrac_new(r,s,t);
-                                    Vbar(r,s,t) += vfrac_new(i,j,k);
+                    ubar(i,j,k,n) = U_in(i,j,k,n);
+
+//                     // Now add to Ubar if cell has EB
+//                     if (vfrac_old(i,j,k) < 1.0 || vfrac_new(i,j,k) < 1.0)
+//                     {
+//                         for (int i_nbor = 1; i_nbor <= itr(i,j,k,0); i_nbor++)
+//                         {
+//                             // FIXME - check to make sure the box contains rst?
+//                             int r = i+map[0][itr(i,j,k,i_nbor)];
+//                             int s = j+map[1][itr(i,j,k,i_nbor)];
+//                             int t = k+((AMREX_SPACEDIM < 3) ? 0 : map[2][itr(i,j,k,i_nbor)]);
+
+//                             // Add to my Ubar
+//                             // Now add me to my nbs Ubar
+//                             // FIXME - does this double weight if nbhds reciprocal?
+//                             if ( n == 0 ){
+//                                 Vbar(i,j,k) += vfrac_old(r,s,t);
+//                                 Vbar(r,s,t) += vfrac_old(i,j,k);
+//                                 // Vbar(i,j,k) += alpha(i,j,k,1)*vfrac_old(r,s,t)/nrs(r,s,t);
+//                                 if (vel_eb_new){
+//                                     Vbar(i,j,k) += vfrac_new(r,s,t);
+//                                     Vbar(r,s,t) += vfrac_new(i,j,k);
                                     
-                                    if ( ((i==15 || i==16) && (j==10)) ||
-                                         ((r==15 || r==16) && (s==10)) ){
-                                        AllPrint()<<i<<"i Vbar "<<Vbar(i,j,k)<<std::endl;
-                                        AllPrint()<<r<<"r Vbar "<<Vbar(r,s,t)<<std::endl;
-                                    }
+//                                     // if ( ((i==15 || i==16) && (j==10)) ||
+//                                     //      ((r==15 || r==16) && (s==10)) ){
+//                                     //     AllPrint()<<i<<"i Vbar "<<Vbar(i,j,k)<<std::endl;
+//                                     //     AllPrint()<<r<<"r Vbar "<<Vbar(r,s,t)<<std::endl;
+//                                     // }
                                     
-                                }
-                            }
-
-// already doing n with ParallelFor
-                            //for ( int n = 0; n < ncomp; n++){
-                            ubar(i,j,k,n) += vfrac_old(r,s,t)*U_in(r,s,t,n);
-                            ubar(r,s,t,n) += vfrac_old(i,j,k)*U_in(i,j,k,n);
-                            // ubar(i,j,k,n) += alpha(i,j,k,1)*vfrac_old(r,s,t)*U_in(r,s,t,n)/nrs(r,s,t);
-                            // Vbar += Real(1.);
-                            // ubar(i,j,k,n) += U_in(r,s,t,n);
-                            //}
-
-                            if (vel_eb_new) {
-                                ubar(i,j,k,n) += vfrac_new(r,s,t)*out(r,s,t,n);
-                                ubar(r,s,t,n) += vfrac_new(i,j,k)*out(i,j,k,n);
-
-                                if ( ((i==15 || i==16) && (j==10)) ||
-                                     ((r==15 || r==16) && (s==10)) ){
-                                    AllPrint()<<i<<"i ubar "<<ubar(i,j,k)<<std::endl;
-                                    AllPrint()<<r<<"r ubar "<<ubar(r,s,t)<<std::endl;
-                                }
-
-                            }
-                        }
-                    }
+//                                 }
+//                             }
+
+// // already doing n with ParallelFor
+//                             //for ( int n = 0; n < ncomp; n++){
+//                             ubar(i,j,k,n) += vfrac_old(r,s,t)*U_in(r,s,t,n);
+//                             ubar(r,s,t,n) += vfrac_old(i,j,k)*U_in(i,j,k,n);
+//                             // ubar(i,j,k,n) += alpha(i,j,k,1)*vfrac_old(r,s,t)*U_in(r,s,t,n)/nrs(r,s,t);
+//                             // Vbar += Real(1.);
+//                             // ubar(i,j,k,n) += U_in(r,s,t,n);
+//                             //}
+
+//                             if (vel_eb_new) {
+//                                 ubar(i,j,k,n) += vfrac_new(r,s,t)*out(r,s,t,n);
+//                                 ubar(r,s,t,n) += vfrac_new(i,j,k)*out(i,j,k,n);
+
+//                                 if ( ((i==15 || i==16) && (j==10)) ||
+//                                      ((r==15 || r==16) && (s==10)) ){
+//                                     AllPrint()<<i<<"i ubar "<<ubar(i,j,k)<<std::endl;
+//                                     AllPrint()<<r<<"r ubar "<<ubar(r,s,t)<<std::endl;
+//                                 }
+
+//                             }
+                    //      }
+                    // }
                 }
                 else
                 {
                     scratch(i,j,k,n) = U_in(i,j,k,n) + dt * dUdt_in(i,j,k,n);
                     kappa(i,j,k) = 0.;
-                    // Already initialized now: ubar(i,j,k,n) = U_in(i,j,k,n);
+                    // NOt Already initialized now:
+                    ubar(i,j,k,n) = U_in(i,j,k,n);
                 }
             });
 
-            // Finish Ubar
-            // FIXME? perhaps this fits better in state_utils
-            amrex::ParallelFor(Box(scratch),
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-            {
-                for (int n = 0; n < ncomp; n++){
-                    if ( Vbar(i,j,k) > 0. )
-                    {
-                    if ((i==15 || i==16) && (j==10)){
-                        AllPrint()<<i<<" pieces "<<ubar(i,j,k,n)<<" "<<Vbar(i,j,k)<<std::endl;
-                        AllPrint()<<i<<"F Ubar "<<ubar(i,j,k,n)/Vbar(i,j,k)<<std::endl;
-                    }
-
-                        ubar(i,j,k,n) /= Vbar(i,j,k);
-                    }
-
-                    // if (i==15 && (j==10||j==9 || j==12)){
-                    //     Print()<<j<<" Ubar "<<ubar(i,j,k,n)<<std::endl;
-                    //     Print()<<j<<" kappa "<<kappa(i,j,k)<<std::endl;
-                    // }
-                    // if ((i==8 || i==9) && j==8)
-                        // {
-                        //     Print()<<Dim3{r,s,t}<<"alpha, beta, N : "<<alpha(r,s,t,0)<<" "<<alpha(r,s,t,1)
-                        //            <<" "<<nrs(r,s,t)<<std::endl;
-                        // }
-
-                        //FIXME - just using my U^n is just as good as any of the averages I tried
-                        // for nbhd kappa. Averaged Ubar better if using U^n+1 based kappa.
-                        // better than using alpha, beta (for nbhd kappa, not looked at for U^n+1)...
-                //ubar(i,j,k,n) = Real(2.); //U_in(i,j,k,n);
-                    //}
-                }
-            });
+            // // Finish Ubar
+            // // FIXME? perhaps this fits better in state_utils
+            // amrex::ParallelFor(Box(scratch),
+            // [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            // {
+            //     for (int n = 0; n < ncomp; n++){
+            //         if ( Vbar(i,j,k) > 0. )
+            //         {
+            //         // if ((i==15 || i==16) && (j==10)){
+            //         //     AllPrint()<<i<<" pieces "<<ubar(i,j,k,n)<<" "<<Vbar(i,j,k)<<std::endl;
+            //         //     AllPrint()<<i<<"F Ubar "<<ubar(i,j,k,n)/Vbar(i,j,k)<<std::endl;
+            //         // }
+
+            //             ubar(i,j,k,n) /= Vbar(i,j,k);
+            //         }
+
+            //         // if (i==15 && (j==10||j==9 || j==12)){
+            //         //     Print()<<j<<" Ubar "<<ubar(i,j,k,n)<<std::endl;
+            //         //     Print()<<j<<" kappa "<<kappa(i,j,k)<<std::endl;
+            //         // }
+            //         // if ((i==8 || i==9) && j==8)
+            //             // {
+            //             //     Print()<<Dim3{r,s,t}<<"alpha, beta, N : "<<alpha(r,s,t,0)<<" "<<alpha(r,s,t,1)
+            //             //            <<" "<<nrs(r,s,t)<<std::endl;
+            //             // }
+
+            //             //FIXME - just using my U^n is just as good as any of the averages I tried
+            //             // for nbhd kappa. Averaged Ubar better if using U^n+1 based kappa.
+            //             // better than using alpha, beta (for nbhd kappa, not looked at for U^n+1)...
+            //         //ubar(i,j,k,n) = Real(2.); //U_in(i,j,k,n);
+            //         //}
+            //     }
+            // });
         }
 
         //FIXME - For now, use the data in out, need to zero here
diff --git a/src/redistribution/hydro_state_redistribute.cpp b/src/redistribution/hydro_state_redistribute.cpp
index 81415d71..e1833951 100644
--- a/src/redistribution/hydro_state_redistribute.cpp
+++ b/src/redistribution/hydro_state_redistribute.cpp
@@ -103,15 +103,98 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
     if (is_periodic_z) domain_per_grown.grow(2,2);
 #endif
 
+    Box dpg3 = domain;
+    if (is_periodic_x) dpg3.grow(0,3);
+    if (is_periodic_y) dpg3.grow(1,3);
+#if (AMREX_SPACEDIM == 3)
+    if (is_periodic_z) dpg3.grow(2,3);
+#endif
+
     // Solution at the centroid of my nbhd
     FArrayBox    soln_hat_fab (bxg3,ncomp,The_Async_Arena());
     Array4<Real> soln_hat = soln_hat_fab.array();
 
     FArrayBox    kappa_hat_fab (bxg3,ncomp,The_Async_Arena());
     Array4<Real> kappa_hat = kappa_hat_fab.array();
+    FArrayBox    ubar_hat_fab (bxg3,ncomp,The_Async_Arena());
+    Array4<Real> ubar_hat = ubar_hat_fab.array();
 
     FArrayBox    kappa_out_fab (bxg3,ncomp,The_Async_Arena());
     Array4<Real> kappa_out = kappa_out_fab.array();
+    FArrayBox    ubar_out_fab (bxg3,ncomp,The_Async_Arena());
+    Array4<Real> ubar_out = ubar_out_fab.array();
+    FArrayBox    ubar_vhat_fab (bxg3,1,The_Async_Arena());
+    Array4<Real> ubar_vhat = ubar_vhat_fab.array();
+
+    // FIXME - In preparation for 3D case where we need alpha, beta to depend on the state...
+    // I *think* this could work for more than 1 NU cell in the NBHD
+    FArrayBox alpha_meb_fab(bxg3,ncomp,The_Async_Arena());
+    Array4<Real> alpha_meb = alpha_meb_fab.array();
+    FArrayBox beta_meb_fab(bxg3,ncomp,The_Async_Arena());
+    Array4<Real> beta_meb = beta_meb_fab.array();
+    FArrayBox Vhat_meb_fab(bxg3,ncomp,The_Async_Arena());
+    Array4<Real> Vhat_meb = Vhat_meb_fab.array();
+    Real eps = 1.e-14; // This really needs to be different for single v double...
+
+    amrex::ParallelFor(bxg3,
+    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+        // initialize the new alpha, beta
+        for ( int n = 0; n < ncomp; n++){
+            alpha_meb(i,j,k,n) = alpha(i,j,k,0);
+            beta_meb(i,j,k,n) = alpha(i,j,k,1);
+            Vhat_meb(i,j,k,n) = nbhd_vol(i,j,k);
+        }
+
+        ubar_vhat(i,j,k) = alpha(i,j,k,0)*vfrac_old(i,j,k);
+
+        // This loops over the neighbors of (i,j,k), and doesn't include (i,j,k) itself
+        for (int i_nbor = 1; i_nbor <= itracker(i,j,k,0); i_nbor++)
+        {
+            int r = i+imap[itracker(i,j,k,i_nbor)];
+            int s = j+jmap[itracker(i,j,k,i_nbor)];
+            int t = k+kmap[itracker(i,j,k,i_nbor)];
+            amrex::Gpu::Atomic::Add(&ubar_vhat(i,j,k),alpha(i,j,k,1) * vfrac_old(r,s,t) / nrs(r,s,t));
+        }
+    });
+
+// THis works but doesn't really change anything in the way of results...
+    // // Fix beta for newly covered cells so that the update really has 0 in the covered cell
+    // amrex::ParallelFor(bxg3 & dpg3, ncomp,
+    // [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+    // {
+    //     if (vfrac_old(i,j,k) > 0.0 && vfrac_new(i,j,k) == 0.0)
+    //     {
+    //         if ( std::abs(U_in(i,j,k,n)) > eps ){
+    //             Real nbs = 0.;
+    //             // This loops over the neighbors of (i,j,k), and doesn't include (i,j,k) itself
+    //             for (int i_nbor = 1; i_nbor <= itracker(i,j,k,0); i_nbor++)
+    //             {
+    //                 int r = i+imap[itracker(i,j,k,i_nbor)];
+    //                 int s = j+jmap[itracker(i,j,k,i_nbor)];
+    //                 int t = k+kmap[itracker(i,j,k,i_nbor)];
+    //                 nbs += (vfrac_old(r,s,t)*U_in(r,s,t,n) + kappa(r,s,t)*ubar(r,s,t,n)) / nrs(r,s,t);
+    //             }
+    //             if ( std::abs(nbs) > eps ){
+    //                 beta_meb(i,j,k,n) = -(vfrac_old(i,j,k)*U_in(i,j,k,n) + kappa(i,j,k)*ubar(i,j,k,n))/nbs;
+
+    //                 // Also need to fix alpha, Vhat of my neighbors
+    //                 for (int i_nbor = 1; i_nbor <= itracker(i,j,k,0); i_nbor++)
+    //                 {
+    //                     int r = i+imap[itracker(i,j,k,i_nbor)];
+    //                     int s = j+jmap[itracker(i,j,k,i_nbor)];
+    //                     int t = k+kmap[itracker(i,j,k,i_nbor)];
+    //                     amrex::Gpu::Atomic::Add(&alpha_meb(r,s,t,n), (alpha(i,j,k,1) - beta_meb(i,j,k,n))/nrs(r,s,t));
+    //                     amrex::Gpu::Atomic::Add(&Vhat_meb(r,s,t,n), (alpha_meb(r,s,t,n) - alpha(r,s,t,0))*vfrac_new(r,s,t));
+    //                 }
+    //             }
+    //             else {
+    //                 std::printf("Newly covered (%i,%i,%i) %i %f\n", i, j, k, itracker(i,j,k,0), U_in(i,j,k,n));
+    //                 Abort();
+    //             }
+    //         }
+    //     }
+    // });
 
     // Define Qhat (from Berger and Guliani) - the weighted solution average
     // Here we initialize soln_hat to equal U_in on all cells in bxg3 so that
@@ -123,6 +206,8 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
     {
         for (int n = 0; n < ncomp; n++){
             soln_hat(i,j,k,n) = U_in(i,j,k,n);
+            ubar_hat(i,j,k,n) = ubar(i,j,k,n);
+            ubar_out(i,j,k,n) = 0.;
             kappa_hat(i,j,k,n) = 0.;
             kappa_out(i,j,k,n) = 0.;
         }
@@ -133,8 +218,10 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
         {
             // Start with U_in(i,j,k) itself
             for (int n = 0; n < ncomp; n++){
-                soln_hat(i,j,k,n) = alpha(i,j,k,0) * U_in(i,j,k,n) * vfrac_old(i,j,k);
-                kappa_hat(i,j,k,n) = alpha(i,j,k,0) * kappa(i,j,k);
+                soln_hat(i,j,k,n) = alpha_meb(i,j,k,n) * U_in(i,j,k,n) * vfrac_old(i,j,k);
+                ubar_hat(i,j,k,n) = alpha_meb(i,j,k,n) * ubar(i,j,k,n) * vfrac_old(i,j,k);
+                kappa_hat(i,j,k,n) = alpha_meb(i,j,k,n) * kappa(i,j,k) * ubar(i,j,k,n);
+                //kappa_hat(i,j,k,n) = kappa(i,j,k);
             }
 
             // This loops over the neighbors of (i,j,k), and doesn't include (i,j,k) itself
@@ -147,12 +234,14 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
                 if (domain_per_grown.contains(IntVect(AMREX_D_DECL(r,s,t))))
                 {
                     for (int n = 0; n < ncomp; n++){
-                        soln_hat(i,j,k,n) += alpha(i,j,k,1) * U_in(r,s,t,n)*vfrac_old(r,s,t) / nrs(r,s,t);
-                        kappa_hat(i,j,k,n) += alpha(i,j,k,1) * kappa(r,s,t) / nrs(r,s,t);
+                        soln_hat(i,j,k,n) += beta_meb(i,j,k,n) * U_in(r,s,t,n)*vfrac_old(r,s,t) / nrs(r,s,t);
+                        ubar_hat(i,j,k,n) += beta_meb(i,j,k,n) * ubar(r,s,t,n)*vfrac_old(r,s,t) / nrs(r,s,t);
+                        kappa_hat(i,j,k,n) += beta_meb(i,j,k,n) * kappa(r,s,t) * ubar(r,s,t,n)/ nrs(r,s,t);
+                        //kappa_hat(i,j,k,n) += kappa(r,s,t);
                     }
-
                 }
             }
+            // FIXME - do we need this check? Change for Vhat ...
             if (nbhd_vol(i,j,k) < 1e-14 &&
                 !(std::abs(alpha(i,j,k,0)) < 1e-14 && std::abs(alpha(i,j,k,1)) < 1e-14) ) //&&
                 // vfrac_new(i,j,k) > 0. ) // Don't form NC nbhd like NU just yet
@@ -163,38 +252,20 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
             for (int n = 0; n < ncomp; n++)  {
                 if (nbhd_vol(i,j,k) < 1e-14 ){
                     // If we're here alpha = beta = 0, and could also use Qhat = vfold*U_in/vfnew
-                    // see comment below...
+                    // see comment below... FIXME??
                     soln_hat(i,j,k,n) = Real(0.0);
+                    ubar_hat(i,j,k,n) = Real(0.0);
                     // FIXME - think about if check needed here...
                     kappa_hat(i,j,k,n) = Real(0.0);
                 } else {
-                    soln_hat(i,j,k,n) /= nbhd_vol(i,j,k);
-                    kappa_hat(i,j,k,n) /= nbhd_vol(i,j,k);
+                    soln_hat(i,j,k,n) /= Vhat_meb(i,j,k,n);
+                    ubar_hat(i,j,k,n) /= ubar_vhat(i,j,k);
+                    kappa_hat(i,j,k,n) /= Vhat_meb(i,j,k,n);
                 }
             }
         }
     });
 
-    // FIXME - In preparation for 3D case where we need alpha, beta to depend on the state...
-    // I *think* this could work for more than 1 NU cell in the NBHD
-    FArrayBox alpha_meb_fab(bxg3,ncomp,The_Async_Arena());
-    Array4<Real> alpha_meb = alpha_meb_fab.array();
-    FArrayBox beta_meb_fab(bxg3,ncomp,The_Async_Arena());
-    Array4<Real> beta_meb = beta_meb_fab.array();
-    // FArrayBox Vhat_meb_fab(bxg3,ncomp,The_Async_Arena());
-    // Array4<Real> Vhat_meb = Vhat_meb_fab.array();
-
-    amrex::ParallelFor(bxg3,
-    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-    {
-        // initialize the new alpha, beta
-        for ( int n = 0; n < ncomp; n++){
-            alpha_meb(i,j,k,n) = alpha(i,j,k,0);
-            beta_meb(i,j,k,n) = alpha(i,j,k,1);
-            // Vhat_meb(i,j,k,n) = nbhd_vol(i,j,k);
-        }
-    });
-
     // FIXME - hack for NU cell, where it's merging partner winds up with Q-hat=0
     // For now, we use the NU Q-hat for making slopes...
     // An equally valid solution may be to set Q-hat=(V^n/V^n+1)U-hat above
@@ -202,6 +273,7 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
     FArrayBox    ssoln_hat_fab (bxg4,ncomp,The_Async_Arena());
     Array4<Real> slope_soln_hat = ssoln_hat_fab.array();
 
+    //FIXME -- need to think about ubar here...
     // FIXME - this breaks with OOB error with bxg3...
     // above loop only computes full soln_hat on bxg2 ...
     amrex::ParallelFor(bxg2 & domain_per_grown,
@@ -227,7 +299,7 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
                     // garantee the MP slope_soln_hat will be correct because maybe this comes
                     // first and then the loop hits MP and slope_soln_hat gets reset, or maybe
                     // the order is opposite...
-                    soln_hat(r,s,t,n) += alpha(r,s,t,1)*U_in(i,j,k,n)/nrs(i,j,k)/nbhd_vol(r,s,t);
+                    soln_hat(r,s,t,n) += beta_meb(r,s,t,n)*U_in(i,j,k,n)/nrs(i,j,k)/Vhat_meb(r,s,t,n);
 
                     slope_soln_hat(i,j,k,n) = soln_hat(r,s,t,n);
                 }
@@ -258,7 +330,9 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
                         // amrex::Gpu::Atomic::Add(&U_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*
                         //                         (soln_hat(i,j,k,n) + kappa_hat(i,j,k,n)*ubar(i,j,k,n)));
                         amrex::Gpu::Atomic::Add(&U_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*soln_hat(i,j,k,n));
-                        amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*kappa_hat(i,j,k,n));
+                        amrex::Gpu::Atomic::Add(&ubar_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*ubar_hat(i,j,k,n));
+                        amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*kappa_hat(i,j,k,n));//*ubar(i,j,k,n));
+                        //amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),kappa_hat(i,j,k,n)*ubar_hat(i,j,k,n));
                     }
                 }
             }
@@ -382,7 +456,9 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
                         // amrex::Gpu::Atomic::Add(&U_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*
                         //                         (update + kappa_hat(i,j,k,n)*ubar(i,j,k,n)));
                         amrex::Gpu::Atomic::Add(&U_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*update);
-                        amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*kappa_hat(i,j,k,n));
+                        amrex::Gpu::Atomic::Add(&ubar_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*ubar_hat(i,j,k,n));
+                        amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),alpha_meb(i,j,k,n)*nrs(i,j,k)*kappa_hat(i,j,k,n)); //*ubar(i,j,k,n));
+                        //amrex::Gpu::Atomic::Add(&kappa_out(i,j,k,n),kappa_hat(i,j,k,n)*ubar_hat(i,j,k,n));
                     } // if bx contains
 
                     // This loops over the neighbors of (i,j,k), and doesn't include (i,j,k) itself
@@ -401,7 +477,9 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
                             // amrex::Gpu::Atomic::Add(&U_out(r,s,t,n),beta_meb(i,j,k,n)*
                             //                         (update + kappa_hat(i,j,k,n)*ubar(r,s,t,n)));
                             amrex::Gpu::Atomic::Add(&U_out(r,s,t,n),beta_meb(i,j,k,n)*update);
-                            amrex::Gpu::Atomic::Add(&kappa_out(r,s,t,n),beta_meb(i,j,k,n)*kappa_hat(i,j,k,n));
+                            amrex::Gpu::Atomic::Add(&ubar_out(r,s,t,n),beta_meb(i,j,k,n)*ubar_hat(i,j,k,n));
+                            amrex::Gpu::Atomic::Add(&kappa_out(r,s,t,n),beta_meb(i,j,k,n)*kappa_hat(i,j,k,n)); //*ubar(r,s,t,n));
+                            //amrex::Gpu::Atomic::Add(&kappa_out(r,s,t,n),kappa(i,j,k)*ubar_hat(r,s,t,n));
                         } // if bx contains
                     } // i_nbor
                 } // n
@@ -416,11 +494,81 @@ Redistribution::StateRedistribute ( Box const& bx, int ncomp,
         {
             // This seems to help with a compiler issue ...
             Real denom = 1. / (nrs(i,j,k) + 1.e-40);
-            U_out(i,j,k,n) += kappa_out(i,j,k,n)*ubar(i,j,k,n);
+            //U_out(i,j,k,n) += kappa_out(i,j,k,n);//*ubar_out(i,j,k,n);
             U_out(i,j,k,n) *= denom;
+            //U_out(i,j,k,n) += kappa_out(i,j,k,n);
+            
+            if ( i==15 && j==12){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<kappa_hat(i,j,k,n)
+                       <<" "<<beta_meb(i,j,k,n)/nrs(i,j,k)
+                       <<" "<<alpha_meb(i,j,k,n)
+                       // <<" "<<ubar(i,j,k,n)
+                       // <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==16 && j==12){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==15 && j==13){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==16 && j==13){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)
+                       <<" "<<std::endl;
+            }
         }
         else
         {
+            if ( i==15 && j==12){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)+kappa_out(i,j,k,n)*ubar(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==16 && j==12){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)+kappa_out(i,j,k,n)*ubar(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==15 && j==13){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)+kappa_out(i,j,k,n)*ubar(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+            if ( i==16 && j==13){
+                Print()<<i<<","<<j<<" U^n+1 pieces "
+                       <<kappa_out(i,j,k,n)
+                       <<" "<<ubar(i,j,k,n)
+                       <<" "<<ubar_hat(i,j,k,n)
+                       <<" "<<U_out(i,j,k,n)+kappa_out(i,j,k,n)*ubar(i,j,k,n)
+                       <<" "<<std::endl;
+            }
+
             U_out(i,j,k,n) = 1.e40;
         }
     });