diff --git a/doc/source/run/parameters.rst b/doc/source/run/parameters.rst
index ad08b7ca03..48e74492ce 100644
--- a/doc/source/run/parameters.rst
+++ b/doc/source/run/parameters.rst
@@ -39,14 +39,6 @@ General parameters
     | Output period. No output is given for `hipace.output_period = -1`.
     | **Warning:** `hipace.output_period = 0` will make the simulation crash.
 
-* ``hipace.output_slice`` (`bool`) optional (default `0`)
-    | Gives only a 2D slice output in the XZ-plane. The output is averaged over
-      the two central grid points of the y-axis.
-
-* ``hipace.3d_on_host`` (`bool`) optional (default `0`)
-    | Allocates the 3D arrays of the fields for I/O on the host, and not the device.
-      This enables to run simulations, where the 3D array exceeds the GPU memory.
-
 * ``hipace.beam_injection_cr`` (`integer`) optional (default `1`)
     | Using a temporary coarsed grid for beam particle injection for a fixed particle-per-cell beam.
       For very high-resolution simulations, where the number of grid points (`nx*ny*nz`)
diff --git a/examples/beam_in_vacuum/inputs_SI b/examples/beam_in_vacuum/inputs_SI
index f70bd6154c..fed13b287c 100644
--- a/examples/beam_in_vacuum/inputs_SI
+++ b/examples/beam_in_vacuum/inputs_SI
@@ -36,3 +36,5 @@ beam2.density = 1.4119793504295784e23 # at this density, 1/kp = 10um
 beam2.u_mean = 0. 0. 1.e3
 beam2.u_std = 0. 0. 0.
 beam2.ppc = 2 2 1
+
+diagnostic.diag_type = xyz
diff --git a/examples/beam_in_vacuum/inputs_normalized b/examples/beam_in_vacuum/inputs_normalized
index 467bf4f8b6..81c909c732 100644
--- a/examples/beam_in_vacuum/inputs_normalized
+++ b/examples/beam_in_vacuum/inputs_normalized
@@ -27,3 +27,5 @@ beam.density = 1.0
 beam.u_mean = 0. 0. 1.e3
 beam.u_std = 0. 0. 0.
 beam.ppc = 2 2 1
+
+diagnostic.diag_type = xyz
diff --git a/examples/blowout_wake/analysis_slice_IO.py b/examples/blowout_wake/analysis_slice_IO.py
index bc2df45f86..c34ce4f85f 100755
--- a/examples/blowout_wake/analysis_slice_IO.py
+++ b/examples/blowout_wake/analysis_slice_IO.py
@@ -22,33 +22,57 @@
                                      left_edge=ds1.domain_left_edge,
                                      dims=ds1.domain_dimensions)
 F_full = all_data_level_0[field].v.squeeze()
-F_full_sliced = (F_full[:,F_full.shape[1]//2,:].squeeze() +
-                  F_full[:,F_full.shape[1]//2-1,:].squeeze())/2.
-ds2 = AMReXDataset('slice_io')
+F_full_xz = (F_full[:,F_full.shape[1]//2,:].squeeze() +
+             F_full[:,F_full.shape[1]//2-1,:].squeeze())/2.
+F_full_yz = (F_full[F_full.shape[0]//2,:,:].squeeze() +
+             F_full[F_full.shape[0]//2-1,:,:].squeeze())/2.
+
+ds2 = AMReXDataset('slice_io_xz')
 ad = ds2.all_data()
-F_slice = ad[field].reshape(ds2.domain_dimensions).v.squeeze()
+F_slice_xz = ad[field].reshape(ds2.domain_dimensions).v.squeeze()
+
+ds3 = AMReXDataset('slice_io_yz')
+ad = ds3.all_data()
+F_slice_yz = ad[field].reshape(ds3.domain_dimensions).v.squeeze()
 
 if do_plot:
-    plt.figure(figsize=(12,4))
-    plt.subplot(131)
-    plt.title('full')
-    plt.imshow(F_full_sliced)
+    plt.figure(figsize=(12,8))
+    plt.subplot(231)
+    plt.title('full_xz')
+    plt.imshow(F_full_xz)
+    plt.colorbar()
+    plt.subplot(232)
+    plt.title('slice xz')
+    plt.imshow(F_slice_xz)
+    plt.colorbar()
+    plt.subplot(233)
+    plt.title('difference')
+    plt.imshow(F_slice_xz-F_full_xz)
+    plt.colorbar()
+
+    plt.subplot(234)
+    plt.title('full_yz')
+    plt.imshow(F_full_yz)
     plt.colorbar()
-    plt.subplot(132)
-    plt.title('slice')
-    plt.imshow(F_slice)
+    plt.subplot(235)
+    plt.title('slice yz')
+    plt.imshow(F_slice_yz)
     plt.colorbar()
-    plt.subplot(133)
+    plt.subplot(236)
     plt.title('difference')
-    plt.imshow(F_slice-F_full_sliced)
+    plt.imshow(F_slice_yz-F_full_yz)
     plt.colorbar()
     plt.tight_layout()
     plt.savefig("image.pdf", bbox_inches='tight')
 
-error = np.max(np.abs(F_slice-F_full_sliced)) / np.max(np.abs(F_full_sliced))
+error_xz = np.max(np.abs(F_slice_xz-F_full_xz)) / np.max(np.abs(F_full_xz))
+error_yz = np.max(np.abs(F_slice_yz-F_full_yz)) / np.max(np.abs(F_full_yz))
 
 print("F_full.shape", F_full.shape)
-print("F_slice.shape", F_slice.shape)
-print("error", error)
+print("F_slice_xz.shape", F_slice_xz.shape)
+print("F_slice_yz.shape", F_slice_yz.shape)
+print("error_xz", error_xz)
+print("error_yz", error_yz)
 
-assert(error < 1.e-14)
+assert(error_xz < 1.e-14)
+assert(error_yz < 1.e-14)
diff --git a/examples/blowout_wake/inputs_SI b/examples/blowout_wake/inputs_SI
index 2e25755ce3..94533c27ee 100644
--- a/examples/blowout_wake/inputs_SI
+++ b/examples/blowout_wake/inputs_SI
@@ -36,3 +36,5 @@ beam.ppc = 1 1 1
 plasma.density =  2.8239587008591567e23 # at this density, 1/kp = 10um
 plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
+
+diagnostic.diag_type = xyz
diff --git a/examples/blowout_wake/inputs_normalized b/examples/blowout_wake/inputs_normalized
index b8259113b9..4215487452 100644
--- a/examples/blowout_wake/inputs_normalized
+++ b/examples/blowout_wake/inputs_normalized
@@ -37,3 +37,5 @@ beam.ppc = 1 1 1
 plasma.density = 1.
 plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
+
+diagnostic.diag_type = xyz
diff --git a/examples/gaussian_weight/inputs_SI b/examples/gaussian_weight/inputs_SI
index 66e127edc9..26e3f5e457 100644
--- a/examples/gaussian_weight/inputs_SI
+++ b/examples/gaussian_weight/inputs_SI
@@ -26,3 +26,5 @@ beam.radius = 1.
 beam.u_mean = 0. 0. 1.e3
 beam.u_std = 0. 0. 0.
 beam.ppc = 2 2 1
+
+diagnostic.diag_type = xyz
diff --git a/examples/gaussian_weight/inputs_normalized b/examples/gaussian_weight/inputs_normalized
index 4508ec9b75..0c1ca14cf3 100644
--- a/examples/gaussian_weight/inputs_normalized
+++ b/examples/gaussian_weight/inputs_normalized
@@ -26,3 +26,5 @@ beam.zmax =  20.
 beam.radius = 40.
 beam.u_mean = 0. 0. 1.e3
 beam.u_std = 3. 4. 0.
+
+diagnostic.diag_type = xyz
diff --git a/examples/linear_wake/inputs_SI b/examples/linear_wake/inputs_SI
index 6b2cd7d844..3ceb84ad35 100644
--- a/examples/linear_wake/inputs_SI
+++ b/examples/linear_wake/inputs_SI
@@ -1,6 +1,5 @@
 amr.n_cell = 32 32 200
 
-hipace.3d_on_host = 1
 hipace.output_plasma = 1
 
 amr.blocking_factor = 2
@@ -33,3 +32,5 @@ beam.ppc = 1 1 1
 plasma.density =  2.8239587008591567e23 # at this density, 1/kp = 10um
 plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
+
+diagnostic.diag_type = xyz
diff --git a/examples/linear_wake/inputs_normalized b/examples/linear_wake/inputs_normalized
index 8a558fb725..1c32a77563 100644
--- a/examples/linear_wake/inputs_normalized
+++ b/examples/linear_wake/inputs_normalized
@@ -1,7 +1,6 @@
 amr.n_cell = 32 32 200
 
 hipace.normalized_units=1
-hipace.3d_on_host = 1
 
 hipace.output_period = 1
 
@@ -34,3 +33,5 @@ beam.ppc = 1 1 1
 plasma.density = 1.
 plasma.ppc = 1 1
 plasma.u_mean = 0.0 0.0 0.
+
+diagnostic.diag_type = xyz
diff --git a/src/Hipace.H b/src/Hipace.H
index ac1817cf8f..8a324ed90c 100644
--- a/src/Hipace.H
+++ b/src/Hipace.H
@@ -193,13 +193,8 @@ public:
     /** Mixing factor between the transverse B field iterations in the predictor corrector loop
      */
     static amrex::Real m_predcorr_B_mixing_factor;
-    /** whether the 3D array stays in host memory (Pinned memory) or is allowed
-     * to be in device memory (allocated in Managed memory). */
-    static bool m_3d_on_host;
     /** Whether to call amrex::Gpu::synchronize() around all profiler region */
     static bool m_do_device_synchronize;
-    /** whether to output the central y=0 slice (if true), or the whole 3D array (if false) */
-    static bool m_slice_F_xz;
     /** Whether to dump plasma particles */
     static bool m_output_plasma;
     /** How much the box is coarsened for beam injection, to avoid exceeding max int in cell count.
diff --git a/src/Hipace.cpp b/src/Hipace.cpp
index 4ea3aaa2c3..98baf7efef 100644
--- a/src/Hipace.cpp
+++ b/src/Hipace.cpp
@@ -31,9 +31,7 @@ int Hipace::m_depos_order_z = 0;
 amrex::Real Hipace::m_predcorr_B_error_tolerance = 4e-2;
 int Hipace::m_predcorr_max_iterations = 30;
 amrex::Real Hipace::m_predcorr_B_mixing_factor = 0.05;
-bool Hipace::m_3d_on_host = false;
 bool Hipace::m_do_device_synchronize = false;
-bool Hipace::m_slice_F_xz = false;
 bool Hipace::m_output_plasma = false;
 int Hipace::m_beam_injection_cr = 1;
 amrex::Real Hipace::m_external_focusing_field_strength = 0.;
@@ -76,9 +74,7 @@ Hipace::Hipace () :
     pph.query("predcorr_max_iterations", m_predcorr_max_iterations);
     pph.query("predcorr_B_mixing_factor", m_predcorr_B_mixing_factor);
     pph.query("output_period", m_output_period);
-    pph.query("output_slice", m_slice_F_xz);
     pph.query("beam_injection_cr", m_beam_injection_cr);
-    pph.query("3d_on_host", m_3d_on_host);
     m_numprocs_z = amrex::ParallelDescriptor::NProcs() / (m_numprocs_x*m_numprocs_y);
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_numprocs_x*m_numprocs_y*m_numprocs_z
                                      == amrex::ParallelDescriptor::NProcs(),
@@ -307,7 +303,7 @@ Hipace::Evolve ()
         const amrex::Vector<int> index_array = fields.IndexArray();
         for (auto it = index_array.rbegin(); it != index_array.rend(); ++it)
         {
-            const amrex::Box& bx = fields.box(*it);
+            const amrex::Box& bx = boxArray(lev)[*it];
             amrex::Vector<amrex::DenseBins<BeamParticleContainer::ParticleType>> bins;
             bins = m_multi_beam.findParticlesInEachSlice(lev, *it, bx, geom[lev]);
 
@@ -391,7 +387,7 @@ Hipace::SolveOneSlice (int islice, int lev, amrex::Vector<amrex::DenseBins<BeamP
     // Push beam particles
     m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice, bins);
 
-    m_fields.Copy(lev, islice, FieldCopyType::StoF, 0, 0, FieldComps::nfields);
+    m_fields.FillDiagnostics(lev, islice);
 
     m_fields.ShiftSlices(lev);
 
@@ -747,6 +743,7 @@ void
 Hipace::WriteDiagnostics (int output_step, bool force_output)
 {
     HIPACE_PROFILE("Hipace::WriteDiagnostics()");
+    constexpr int lev = 0;
 
     // Dump before first and after last step, and every m_output_period steps in-between
     if (m_output_period < 0 ||
@@ -763,30 +760,17 @@ Hipace::WriteDiagnostics (int output_step, bool force_output)
         {"ExmBy", "EypBx", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz", "rho", "Psi"};
 
     amrex::Vector<std::string> rfs;
-    amrex::Vector<amrex::Geometry> geom_io = Geom();
-
-    constexpr int lev = 0; // we only have one level
-    constexpr int idim = 1;
-    amrex::RealBox prob_domain = Geom(lev).ProbDomain();
-    amrex::Box domain = Geom(lev).Domain();
-    // Define slice box
-    int const icenter = domain.length(idim)/2;
-    domain.setSmall(idim, icenter);
-    domain.setBig(idim, icenter);
-    if (m_slice_F_xz){
-        geom_io[lev] = amrex::Geometry(domain, &prob_domain, Geom(lev).Coord());
-    }
 
 #ifdef HIPACE_USE_OPENPMD
-    m_openpmd_writer.WriteDiagnostics(m_fields, m_multi_beam, geom_io[lev], m_physical_time,
-                                      output_step,  lev, m_slice_F_xz, varnames);
+    m_openpmd_writer.WriteDiagnostics(m_fields.getDiagF(), m_multi_beam, m_fields.getDiagGeom(),
+                                      m_physical_time, output_step,  lev, m_fields.getDiagSliceDir(), varnames);
 #else
     constexpr int nlev = 1;
     const amrex::IntVect local_ref_ratio {1, 1, 1};
 
     amrex::WriteMultiLevelPlotfile(
-        filename, nlev, amrex::GetVecOfConstPtrs(m_fields.getF()), varnames,
-        geom_io, m_physical_time, {output_step}, {local_ref_ratio},
+        filename, nlev, amrex::GetVecOfConstPtrs(m_fields.getDiagF()), varnames,
+        m_fields.getDiagGeom(), m_physical_time, {output_step}, {local_ref_ratio},
         "HyperCLaw-V1.1", "Level_", "Cell", rfs);
 
     // Write beam particles
diff --git a/src/diagnostics/CMakeLists.txt b/src/diagnostics/CMakeLists.txt
index f47ffd065f..4f03b0ac60 100644
--- a/src/diagnostics/CMakeLists.txt
+++ b/src/diagnostics/CMakeLists.txt
@@ -1,4 +1,5 @@
 target_sources(HiPACE
   PRIVATE
     OpenPMDWriter.cpp
+    FieldDiagnostic.cpp
 )
diff --git a/src/diagnostics/FieldDiagnostic.H b/src/diagnostics/FieldDiagnostic.H
new file mode 100644
index 0000000000..6267e98c41
--- /dev/null
+++ b/src/diagnostics/FieldDiagnostic.H
@@ -0,0 +1,53 @@
+#ifndef FIELDDIAGNOSTIC_H_
+#define FIELDDIAGNOSTIC_H_
+
+#include <AMReX_MultiFab.H>
+
+/** type of diagnostics: full xyz array or xz slice or yz slice */
+enum struct DiagType{xyz, xz, yz};
+
+/** \brief This class holds data for 1 diagnostics (full or slice) */
+class FieldDiagnostic
+{
+
+public:
+
+    /** \brief Constructor */
+    explicit FieldDiagnostic (int nlev);
+
+    /** \brief allocate arrays of this MF
+     *
+     * \param[in] lev MR level
+     * \param[in] ba BoxArray of the full simulation domain
+     * \param[in] nfields number of field components
+     * \param[in] dm DistributionMapping of the full simulation domain
+     * \param[in] geom geometry of the full simulation domain
+     */
+    void AllocData (int lev, const amrex::BoxArray& ba, int nfields, const amrex::DistributionMapping& dm, amrex::Geometry const& geom);
+
+    /** \brief return the main diagnostics multifab */
+    amrex::Vector<amrex::MultiFab>& getF () { return m_F; };
+
+    /** \brief return the main diagnostics multifab
+     *
+     * \param[in] lev MR level
+     */
+    amrex::MultiFab& getF (int lev) {return m_F[lev];};
+
+    /** return the diagnostics geometry */
+    amrex::Vector<amrex::Geometry>& getGeom () { return m_geom_io; };
+
+    /** return slice direction of the diagnostics */
+    int sliceDir () {return m_slice_dir;};
+
+private:
+
+    /** Vector over levels, all fields */
+    amrex::Vector<amrex::MultiFab> m_F;
+    DiagType m_diag_type; /**< Type of diagnostics (xyz xz yz) */
+    int m_slice_dir; /**< Slicing direction */
+    int m_nfields; /**< Number of physical fields to write */
+    amrex::Vector<amrex::Geometry> m_geom_io; /**< Diagnostics geometry */
+};
+
+#endif // FIELDDIAGNOSTIC_H_
diff --git a/src/diagnostics/FieldDiagnostic.cpp b/src/diagnostics/FieldDiagnostic.cpp
new file mode 100644
index 0000000000..e382bc6afc
--- /dev/null
+++ b/src/diagnostics/FieldDiagnostic.cpp
@@ -0,0 +1,61 @@
+#include "FieldDiagnostic.H"
+
+#include <AMReX_ParmParse.H>
+
+FieldDiagnostic::FieldDiagnostic (int nlev)
+    : m_F(nlev),
+      m_geom_io(nlev)
+{
+    amrex::ParmParse ppd("diagnostic");
+    std::string str_type;
+    ppd.get("diag_type", str_type);
+    if        (str_type == "xyz"){
+        m_diag_type = DiagType::xyz;
+        m_slice_dir = -1;
+    } else if (str_type == "xz") {
+        m_diag_type = DiagType::xz;
+        m_slice_dir = 1;
+    } else if (str_type == "yz") {
+        m_diag_type = DiagType::yz;
+        m_slice_dir = 0;
+    } else {
+        amrex::Abort("Unknown diagnostics type: must be xyz, xz or yz.");
+    }
+}
+
+void
+FieldDiagnostic::AllocData (int lev, const amrex::BoxArray& ba, int nfields, const amrex::DistributionMapping& dm, amrex::Geometry const& geom)
+{
+    m_nfields = nfields;
+    // Create a xz slice BoxArray
+    amrex::BoxList F_boxes;
+    if (m_slice_dir >= 0){
+        for (int i = 0; i < ba.size(); ++i){
+            amrex::Box bx = ba[i];
+            // Flatten the box down to 1 cell in the approprate direction.
+            bx.setSmall(m_slice_dir, ba[i].length(m_slice_dir)/2);
+            bx.setBig  (m_slice_dir, ba[i].length(m_slice_dir)/2);
+            // Note: the MR is still cell-centered, although the data will be averaged to nodal.
+            F_boxes.push_back(bx);
+        }
+    }
+    amrex::BoxArray F_slice_ba(std::move(F_boxes));
+    // m_F is defined on F_ba, the full or the slice BoxArray
+    amrex::BoxArray F_ba = m_slice_dir >= 0 ? F_slice_ba : ba;
+    // Only xy slices need guard cells, there is no deposition to/gather from the output array F.
+    amrex::IntVect nguards_F = amrex::IntVect(0,0,0);
+    // The Arena uses pinned memory.
+    m_F[lev].define(F_ba, dm, m_nfields, nguards_F,
+                    amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()));
+
+    m_geom_io[lev] = geom;
+    amrex::RealBox prob_domain = geom.ProbDomain();
+    amrex::Box domain = geom.Domain();
+    // Define slice box
+    if (m_slice_dir >= 0){
+        int const icenter = domain.length(m_slice_dir)/2;
+        domain.setSmall(m_slice_dir, icenter);
+        domain.setBig(m_slice_dir, icenter);
+        m_geom_io[lev] = amrex::Geometry(domain, &prob_domain, geom.Coord());
+    }
+}
diff --git a/src/diagnostics/OpenPMDWriter.H b/src/diagnostics/OpenPMDWriter.H
index 3ec26bbfe0..1efd6d5cc2 100644
--- a/src/diagnostics/OpenPMDWriter.H
+++ b/src/diagnostics/OpenPMDWriter.H
@@ -58,16 +58,16 @@ private:
 
     /** \brief writing openPMD field data
      *
-     * \param[in] a_fields the field class
+     * \param[in] mf the MultiFab to dump
      * \param[in] geom Geometry of the simulation, to get the cell size etc.
-     * \param[in] lev MR level
-     * \param[in] slice_F_xz whether to use slice IO or not
+     * \param[in] slice_dir direction of slicing. 0=x, 1=y, -1=no slicing (3D array)
      * \param[in] varnames list of variable names for the fields (ExmBy, EypBx, Ey, ...)
      * \param[in,out] iteration openPMD iteration to which the data is written
      */
-    void WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
-                         const int lev, const bool slice_F_xz,
-                         const amrex::Vector< std::string > varnames, openPMD::Iteration iteration);
+    void WriteFieldData (
+        amrex::MultiFab const& mf, amrex::Geometry const& geom,
+        const int slice_dir, const amrex::Vector< std::string > varnames,
+        openPMD::Iteration iteration);
 
     /** Named Beam SoA attributes per particle as defined in BeamIdx
      */
@@ -85,18 +85,20 @@ public:
 
     /** \brief writing openPMD data
      *
-     * \param[in] a_fields the field class
-     * \param[in] beams multi beam container which is written to openPMD file
+     * \param[in] a_mf Vector (levels) of multifabs
+     * \param[in] a_multi_beams multi beam container which is written to openPMD file
      * \param[in] geom Geometry of the simulation, to get the cell size etc.
      * \param[in] physical_time Physical time of the currenerationt it.
      * \param[in] output_step current iteration to be written to file
      * \param[in] lev MR level
-     * \param[in] slice_F_xz whether to use slice IO or not
+     * \param[in] slice_dir direction of slicing. 0=x, 1=y, -1=no slicing (3D array)
      * \param[in] varnames list of variable names for the fields (ExmBy, EypBx, Ey, ...)
      */
-    void WriteDiagnostics(Fields& a_fields, MultiBeam& a_multi_beam, amrex::Geometry const& geom,
-                          const amrex::Real physical_time, const int output_step, const int lev,
-                          const bool slice_F_xz, const amrex::Vector< std::string > varnames);
+    void WriteDiagnostics(
+        amrex::Vector<amrex::MultiFab> const& a_mf, MultiBeam& a_multi_beam,
+        amrex::Vector<amrex::Geometry> const& geom,
+        const amrex::Real physical_time, const int output_step, const int lev,
+        const int slice_dir, const amrex::Vector< std::string > varnames);
 
     void reset () {m_outputSeries.reset();};
 };
diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp
index 6495c6abce..312a24aad8 100644
--- a/src/diagnostics/OpenPMDWriter.cpp
+++ b/src/diagnostics/OpenPMDWriter.cpp
@@ -35,15 +35,16 @@ OpenPMDWriter::InitDiagnostics ()
 }
 
 void
-OpenPMDWriter::WriteDiagnostics (Fields& a_fields, MultiBeam& a_multi_beam,
-                                 amrex::Geometry const& geom, amrex::Real physical_time,
-                                 const int output_step, const int lev, const bool slice_F_xz,
-                                 const amrex::Vector< std::string > varnames)
+OpenPMDWriter::WriteDiagnostics (
+    amrex::Vector<amrex::MultiFab> const& a_mf, MultiBeam& a_multi_beam,
+    amrex::Vector<amrex::Geometry> const& geom,
+    const amrex::Real physical_time, const int output_step, const int lev,
+    const int slice_dir, const amrex::Vector< std::string > varnames)
 {
     io::Iteration iteration = m_outputSeries->iterations[output_step];
     iteration.setTime(physical_time);
 
-    WriteFieldData(a_fields, geom, lev, slice_F_xz, varnames, iteration);
+    WriteFieldData(a_mf[lev], geom[lev], slice_dir, varnames, iteration);
 
     a_multi_beam.ConvertUnits(ConvertDirection::HIPACE_to_SI);
     WriteBeamParticleData(a_multi_beam, iteration);
@@ -55,10 +56,10 @@ OpenPMDWriter::WriteDiagnostics (Fields& a_fields, MultiBeam& a_multi_beam,
 }
 
 void
-OpenPMDWriter::WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
-                               const int lev, const bool slice_F_xz,
-                               const amrex::Vector< std::string > varnames,
-                               openPMD::Iteration iteration)
+OpenPMDWriter::WriteFieldData (
+    amrex::MultiFab const& mf, amrex::Geometry const& geom,
+    const int slice_dir, const amrex::Vector< std::string > varnames,
+    openPMD::Iteration iteration)
 {
     // todo: periodicity/boundary, field solver, particle pusher, etc.
     auto meshes = iteration.meshes;
@@ -66,8 +67,6 @@ OpenPMDWriter::WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
     // loop over field components
     for (int icomp = 0; icomp < FieldComps::nfields; ++icomp)
     {
-        auto const& mf = a_fields.getF(lev);
-
         std::string fieldname = varnames[icomp];
         //                      "B"                "x" (todo)
         //                      "Bx"               ""  (just for now)
@@ -83,11 +82,13 @@ OpenPMDWriter::WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
         std::vector< std::string > axisLabels {"z", "y", "x"};
         auto dCells = utils::getReversedVec(geom.CellSize()); // dx, dy, dz
         auto offWindow = utils::getReversedVec(geom.ProbLo()); // start of moving window
-        if (slice_F_xz) {
-            relative_cell_pos.erase(relative_cell_pos.begin() + 1);  // remove for y
-            axisLabels.erase(axisLabels.begin() + 1); // remove y
-            dCells.erase(dCells.begin() + 1); // remove dy
-            offWindow.erase(offWindow.begin() + 1); // remove offset in y
+        if (slice_dir >= 0) {
+            // User requested slice IO
+            // remove the slicing direction in position, label, resolution, offset
+            relative_cell_pos.erase(relative_cell_pos.begin() + 2-slice_dir);
+            axisLabels.erase(axisLabels.begin() + 2-slice_dir);
+            dCells.erase(dCells.begin() + 2-slice_dir);
+            offWindow.erase(offWindow.begin() + 2-slice_dir);
         }
         field_comp.setPosition(relative_cell_pos);
         field.setAxisLabels(axisLabels);
@@ -97,7 +98,8 @@ OpenPMDWriter::WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
         // data type and global size of the simulation
         io::Datatype datatype = io::determineDatatype< amrex::Real >();
         io::Extent global_size = utils::getReversedVec(geom.Domain().size());
-        if (slice_F_xz) global_size.erase(global_size.begin() + 1);  // remove Ny
+        // If slicing requested, remove number of points for the slicing direction
+        if (slice_dir >= 0) global_size.erase(global_size.begin() + 2-slice_dir);
 
         io::Dataset dataset(datatype, global_size);
         field_comp.resetDataset(dataset);
@@ -123,9 +125,9 @@ OpenPMDWriter::WriteFieldData (Fields& a_fields, amrex::Geometry const& geom,
             amrex::IntVect const box_offset = data_box.smallEnd();
             io::Offset chunk_offset = utils::getReversedVec(box_offset);
             io::Extent chunk_size = utils::getReversedVec(data_box.size());
-            if (slice_F_xz) { // remove Ny components
-                chunk_offset.erase(chunk_offset.begin() + 1);
-                chunk_size.erase(chunk_size.begin() + 1);
+            if (slice_dir >= 0) { // remove Ny components
+                chunk_offset.erase(chunk_offset.begin() + 2-slice_dir);
+                chunk_size.erase(chunk_size.begin() + 2-slice_dir);
             }
 
             field_comp.storeChunk(data, chunk_offset, chunk_size);
diff --git a/src/fields/Fields.H b/src/fields/Fields.H
index d75fbaa27f..970a32457d 100644
--- a/src/fields/Fields.H
+++ b/src/fields/Fields.H
@@ -2,6 +2,7 @@
 #define FIELDS_H_
 
 #include "fft_poisson_solver/FFTPoissonSolver.H"
+#include "diagnostics/FieldDiagnostic.H"
 
 #include <AMReX_MultiFab.H>
 #include <AMReX_Vector.H>
@@ -67,11 +68,6 @@ public:
     /** get function for the main 3D array F
      * \param lev MR level */
     amrex::MultiFab& getF (int lev) { return m_F[lev]; }
-    /** get function for the main 3D array F
-     * \param lev MR level
-     * \param icomp component of the field */
-    amrex::MultiFab getF (int lev, int icomp );
-
     /** get function for the 2D slices */
     amrex::Vector<std::array<amrex::MultiFab, m_nslices>>& getSlices () {return m_slices; }
     /** get function for the 2D slices
@@ -88,10 +84,31 @@ public:
      * \param[in] islice slice index
      */
     amrex::MultiFab& getSlices (int lev, int islice) {return m_slices[lev][islice]; }
-
-    /** Copy between the full MultiFab and slice MultiFab */
+    /** \brief get diagnostics multifab */
+    amrex::Vector<amrex::MultiFab>& getDiagF () { return m_diags.getF(); };
+    /** \brief get diagnostics geometry */
+    amrex::Vector<amrex::Geometry>& getDiagGeom () { return m_diags.getGeom(); };
+    /** \brief get slicing direction of the diagnostics */
+    int getDiagSliceDir () { return m_diags.sliceDir(); };
+    /** \brief Copy the data from xy slices to the field diagnostics.
+     *
+     * \param[in] lev MR leve
+     * \param[in] i_slice z slice in which to write the data
+     */
+    void FillDiagnostics(int lev, int i_slice);
+    /** \brief Copy between the full MultiFab and slice MultiFab.
+     *
+     * \param[in] lev MR level
+     * \param[in] i_slice z slice in which to write the data
+     * \param[in] copy_type whether to copy from the xy slice to the main array or opposite
+     * \param[in] slice_comp first component of the xy slice to copy from/to
+     * \param[in] full_comp first component of the full array to copy from/to
+     * \param[in] ncomp number of components to copy
+     * \param[in,out] full_mf full MultiFab
+     * \param[in] slice_dir slicing direction. 0=x, 1=y, -1=no slicing (full 3D)
+     */
     void Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp,
-               int ncomp);
+               int ncomp, amrex::MultiFab& full_mf, int slice_dir);
 
     /** \brief Shift slices by 1 element: slices (1,2) are then stored in (2,3).
      *
@@ -208,8 +225,6 @@ public:
                                        const int lev);
 
 private:
-    /** Pointer to const singleton of class Hipace */
-    Hipace const* m_hipace;
     /** Vector over levels, all fields */
     amrex::Vector<amrex::MultiFab> m_F;
     /** Vector over levels, array of 4 slices required to compute current slice */
@@ -218,6 +233,8 @@ private:
     amrex::IntVect m_slices_nguards {-1, -1, -1};
     /** Whether to use Dirichlet BC for the Poisson solver. Otherwise, periodic */
     bool m_do_dirichlet_poisson = true;
+    /** Diagnostics */
+    FieldDiagnostic m_diags;
 };
 
 #endif
diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp
index e2210de0ee..0170255d8c 100644
--- a/src/fields/Fields.cpp
+++ b/src/fields/Fields.cpp
@@ -6,9 +6,9 @@
 #include "utils/Constants.H"
 
 Fields::Fields (Hipace const* a_hipace)
-    : m_hipace(a_hipace),
-      m_F(a_hipace->maxLevel()+1),
-      m_slices(a_hipace->maxLevel()+1)
+    : m_F(a_hipace->maxLevel()+1),
+      m_slices(a_hipace->maxLevel()+1),
+      m_diags(a_hipace->maxLevel()+1)
 {
     amrex::ParmParse ppf("fields");
     ppf.query("do_dirichlet_poisson", m_do_dirichlet_poisson);
@@ -24,37 +24,12 @@ Fields::AllocData (
     int nguards_xy = std::max(1, Hipace::m_depos_order_xy);
     m_slices_nguards = {nguards_xy, nguards_xy, 0};
 
-    // Create a xz slice BoxArray
-    amrex::BoxList F_boxes;
-    if (Hipace::m_slice_F_xz){
-        for (int i = 0; i < ba.size(); ++i){
-            amrex::Box bx = ba[i];
-            // Flatten the box down to 1 cell in the y direction.
-            constexpr int idim = 1;
-            bx.setSmall(idim, ba[i].length(idim)/2);
-            bx.setBig(idim, ba[i].length(idim)/2);
-            // Make this MF node-centered so it is exactly at the center of the box.
-            // bx.setType(amrex::IndexType({0,1,0}));
-            F_boxes.push_back(bx);
-        }
-    }
-    amrex::BoxArray F_slice_ba(std::move(F_boxes));
-
-    // m_F is defined on F_ba, the full or the xz slice BoxArray
-    amrex::BoxArray F_ba = Hipace::m_slice_F_xz ? F_slice_ba : ba;
-
     // Only xy slices need guard cells, there is no deposition to/gather from the output array F.
     amrex::IntVect nguards_F = amrex::IntVect(0,0,0);
 
-    if (Hipace::m_3d_on_host){
-        // The Arena uses pinned memory.
-        m_F[lev].define(F_ba, dm, FieldComps::nfields, nguards_F,
-                        amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()));
-    } else {
-        // The Arena uses managed memory.
-        m_F[lev].define(F_ba, dm, FieldComps::nfields, nguards_F,
-                        amrex::MFInfo().SetArena(amrex::The_Arena()));
-    }
+    // The Arena uses pinned memory.
+    m_F[lev].define(ba, dm, FieldComps::nfields, nguards_F, amrex::MFInfo().SetAlloc(false));
+    m_diags.AllocData(lev, ba, FieldComps::nfields, dm, geom);
 
     for (int islice=0; islice<(int) WhichSlice::N; islice++) {
         m_slices[lev][islice].define(slice_ba, slice_dm, FieldComps::nfields, m_slices_nguards,
@@ -163,11 +138,10 @@ Fields::LongitudinalDerivative (const amrex::MultiFab& src1, const amrex::MultiF
 
 void
 Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp,
-              int ncomp)
+              int ncomp, amrex::MultiFab& full_mf, int slice_dir)
 {
     using namespace amrex::literals;
     HIPACE_PROFILE("Fields::Copy()");
-    const bool do_node_center = Hipace::m_slice_F_xz;
     auto& slice_mf = m_slices[lev][(int) WhichSlice::This]; // copy from/to the current slice
     amrex::Array4<amrex::Real> slice_array; // There is only one Box.
     for (amrex::MFIter mfi(slice_mf); mfi.isValid(); ++mfi) {
@@ -179,7 +153,6 @@ Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int
         // slice_array's longitude index is i_slice.
     }
 
-    auto& full_mf = m_F[lev];
     for (amrex::MFIter mfi(full_mf); mfi.isValid(); ++mfi) {
         amrex::Box const& vbx = mfi.validbox();
         if (vbx.smallEnd(Direction::z) <= i_slice and
@@ -199,11 +172,14 @@ Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int
                 amrex::ParallelFor(copy_box, ncomp,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
-                    if (do_node_center){
+                    if        (slice_dir ==-1 /* 3D data */){
+                        full_array(i,j,k,n+full_comp) = slice_array(i,j,k,n+slice_comp);
+                    } else if (slice_dir == 0 /* yz slice */){
+                        full_array(i,j,k,n+full_comp) = 0.5_rt *
+                            (slice_array(i-1,j,k,n+slice_comp)+slice_array(i,j,k,n+slice_comp));
+                    } else /* slice_dir == 1, xz slice */{
                         full_array(i,j,k,n+full_comp) = 0.5_rt *
                             (slice_array(i,j-1,k,n+slice_comp)+slice_array(i,j,k,n+slice_comp));
-                    } else {
-                        full_array(i,j,k,n+full_comp) = slice_array(i,j,k,n+slice_comp);
                     }
                 });
             }
@@ -211,6 +187,13 @@ Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int
     }
 }
 
+void
+Fields::FillDiagnostics (int lev, int i_slice)
+{
+    Copy(lev, i_slice, FieldCopyType::StoF, 0, 0, FieldComps::nfields,
+         m_diags.getF(lev), m_diags.sliceDir());
+}
+
 void
 Fields::ShiftSlices (int lev)
 {
@@ -221,13 +204,6 @@ Fields::ShiftSlices (int lev)
               m_slices[lev][(int) WhichSlice::Previous1]);
 }
 
-amrex::MultiFab
-Fields::getF (int lev, int icomp )
-{
-    amrex::MultiFab F_comp(m_F[lev], amrex::make_alias, icomp, 1);
-    return F_comp;
-}
-
 void
 Fields::AddRhoIons (const int lev)
 {
diff --git a/tests/adaptive_time_step.1Rank.sh b/tests/adaptive_time_step.1Rank.sh
index fe47951089..217782e6f3 100755
--- a/tests/adaptive_time_step.1Rank.sh
+++ b/tests/adaptive_time_step.1Rank.sh
@@ -23,7 +23,6 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
         max_step = 20 \
         geometry.prob_lo = -2. -2. -2. \
         geometry.prob_hi =  2.  2.  2. \
-        hipace.3d_on_host = 1 \
         hipace.dt = 0. \
         hipace.output_period = 20 \
         beam.density = 1.e-8 \
@@ -42,7 +41,6 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
         max_step = 20 \
         geometry.prob_lo = -2. -2. -2. \
         geometry.prob_hi =  2.  2.  2. \
-        hipace.3d_on_host = 1 \
         hipace.dt = 0. \
         hipace.output_period = 20 \
         beam.density = 1.e-8 \
diff --git a/tests/beam_evolution.1Rank.sh b/tests/beam_evolution.1Rank.sh
index e0a8ae8277..013387f78f 100755
--- a/tests/beam_evolution.1Rank.sh
+++ b/tests/beam_evolution.1Rank.sh
@@ -20,7 +20,6 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
         max_step = 20 \
         geometry.prob_lo = -2. -2. -2. \
         geometry.prob_hi =  2.  2.  2. \
-        hipace.3d_on_host = 1 \
         hipace.dt = 3. \
         hipace.output_period = 20 \
         beam.density = 1.e-8 \
diff --git a/tests/slice_IO.1Rank.sh b/tests/slice_IO.1Rank.sh
index 17872e7a4c..cc4bc834fe 100755
--- a/tests/slice_IO.1Rank.sh
+++ b/tests/slice_IO.1Rank.sh
@@ -15,21 +15,30 @@ HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests
 
 rm -rf plt00001
 # Run the simulation
-mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized hipace.output_slice=0
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
+        diagnostic.diag_type=xyz \
+        amr.n_cell = 64 86 100
 rm -rf full_io
 mv plt00001 full_io
-mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized hipace.output_slice=1
-rm -rf slice_io
-mv plt00001 slice_io
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
+        diagnostic.diag_type=xz \
+        amr.n_cell = 64 86 100
+rm -rf slice_io_xz
+mv plt00001 slice_io_xz
+mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
+        diagnostic.diag_type=yz \
+        amr.n_cell = 64 86 100
+rm -rf slice_io_yz
+mv plt00001 slice_io_yz
 
 # assert whether the two IO types match
 $HIPACE_EXAMPLE_DIR/analysis_slice_IO.py
 
 # Make sure that the slice is much smaller than the full diagnostics
 size_full=$(du -s full_io/Level_0 | cut -f1)
-size_slice=$(du -s slice_io/Level_0 | cut -f1)
+size_slice=$(du -s slice_io_yz/Level_0 | cut -f1)
 
-if [[ $((size_full/size_slice<120)) == 1 ]]; then
+if [[ $((size_full/size_slice<60)) == 1 ]]; then
     echo $size_full
     echo $size_slice
     echo "ERROR: field data should be ~128x smaller for slice than for full diagnostics"