diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp
index 16617159662..7a4f032254b 100644
--- a/Source/Diagnostics/BTDiagnostics.cpp
+++ b/Source/Diagnostics/BTDiagnostics.cpp
@@ -308,7 +308,7 @@ BTDiagnostics::DoDump (int step, int i_buffer, bool force_flush)
         // was already fully written and buffer was reset to zero size or that
         // lab snapshot was not even started to be backtransformed yet
         const auto do_forced_flush = (force_flush && !buffer_empty(i_buffer));
-
+        if (do_forced_flush) amrex::Print() << "forcing flush \n";
         return is_buffer_full || last_z_slice_filled || do_forced_flush;
     }
 
@@ -1015,6 +1015,66 @@ BTDiagnostics::Flush (int i_buffer, bool force_flush)
     bool const isBTD = true;
     double const labtime = m_t_lab[i_buffer];
 
+    amrex::Vector< amrex::MultiFab > out;
+    out.resize(nlev_output);
+    for(int lev = 0; lev < nlev_output; ++lev) {
+        if (!force_flush) {
+        out[lev] = amrex::MultiFab(m_mf_output[i_buffer][lev].boxArray(),
+                 m_mf_output[i_buffer][lev].DistributionMap(),
+                 m_mf_output[i_buffer][lev].nComp(), 0);
+        out[lev].ParallelCopy( m_mf_output[i_buffer][lev], 0, 0, m_mf_output[i_buffer][lev].nComp() );
+        } else {
+            auto const& ba = m_mf_output[i_buffer][lev].boxArray();
+            auto const& dm = m_mf_output[i_buffer][lev].DistributionMap();
+            amrex::BoxList bl(ba.ixType());
+            amrex::Vector<int> proc;
+            amrex::Vector<int> imap;
+            for (int i = 0; i < int(ba.size()); ++i) {
+                auto b = ba[i];
+                int new_small = b.bigEnd(m_moving_window_dir)
+                                - (m_buffer_counter[i_buffer]-1);
+                b.setSmall(m_moving_window_dir,
+                         std::max(b.smallEnd(m_moving_window_dir),new_small));
+                if (b.ok()) {
+                    bl.push_back(b);
+                    proc.push_back(dm[i]);
+                    if (dm[i] == amrex::ParallelDescriptor::MyProc()) {
+                        imap.push_back(i);
+                    }
+                }
+            }
+            out[lev] = amrex::MultiFab(amrex::BoxArray(std::move(bl)),
+                       amrex::DistributionMapping(std::move(proc)),
+                       m_mf_output[i_buffer][lev].nComp(),0);
+            for (amrex::MFIter mfi(out[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
+            {
+                auto const& bx = mfi.tilebox();
+                (out[lev])[mfi].template copy<amrex::RunOn::Device>
+                    ( (m_mf_output[i_buffer][lev])[imap[mfi.LocalIndex()]],
+                       bx, 0, bx, 0, m_mf_output[i_buffer][lev].nComp() );
+            }
+            amrex::BoxArray newba = out[lev].boxArray();
+            amrex::Box newbbox = newba.minimalBox();
+            //int new_small = std::max(m_snapshot_box[i_buffer].smallEnd(m_moving_window_dir),
+            //                newbbox.smallEnd(m_moving_window_dir));
+            int ncellsflushed = m_buffer_flush_counter[i_buffer]*m_buffer_size;
+            int currentncells = newbbox.bigEnd(m_moving_window_dir) - newbbox.smallEnd(m_moving_window_dir);
+            int new_small = newbbox.smallEnd(m_moving_window_dir);
+            m_snapshot_box[i_buffer].setSmall(m_moving_window_dir, new_small);
+            auto ref_ratio = amrex::IntVect(1);
+            if (lev > 0 ) { ref_ratio = WarpX::RefRatio(lev-1); }
+            const amrex::Real new_zlo = m_snapshot_domain_lab[i_buffer].hi(m_moving_window_dir) -
+                                 (currentncells + ncellsflushed) *
+                                 dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]);
+            m_snapshot_domain_lab[i_buffer].setLo(m_moving_window_dir, new_zlo);
+            amrex::Vector<int> BTdiag_periodicity(AMREX_SPACEDIM, 0);
+            m_geom_snapshot[i_buffer][lev].define( m_snapshot_box[i_buffer],
+                                               &m_snapshot_domain_lab[i_buffer],
+                                               amrex::CoordSys::cartesian,
+                                               BTdiag_periodicity.data() );
+        }
+    }
+
     amrex::Vector<amrex::BoxArray> vba;
     amrex::Vector<amrex::DistributionMapping> vdmap;
     amrex::Vector<amrex::Geometry> vgeom;
@@ -1064,7 +1124,7 @@ BTDiagnostics::Flush (int i_buffer, bool force_flush)
         }
     }
     m_flush_format->WriteToFile(
-        m_varnames, m_mf_output.at(i_buffer), m_geom_output.at(i_buffer), warpx.getistep(),
+        m_varnames, out, m_geom_output.at(i_buffer), warpx.getistep(),
         labtime,
         m_output_species.at(i_buffer), nlev_output, file_name, m_file_min_digits,
         m_plot_raw_fields, m_plot_raw_fields_guards,