diff --git a/Tests/Particles/Mapped/InitHelix.cpp b/Tests/Particles/Mapped/InitHelix.cpp
new file mode 100644
index 00000000000..1f2d28f7a3b
--- /dev/null
+++ b/Tests/Particles/Mapped/InitHelix.cpp
@@ -0,0 +1,55 @@
+#include <AMReX.H>
+#include <AMReX_MultiFab.H>
+
+using namespace amrex;
+
+void
+InitHelix (amrex::MultiFab& a_xyz_loc, amrex::Geometry& geom)
+{
+    const Real tpi = 2.0* amrex::Math::pi<Real>();
+
+    auto domain = geom.Domain();
+    auto probhi = geom.ProbHiArray();
+    auto problo = geom.ProbLoArray();
+
+    Real ilen = static_cast<Real>(domain.length(0));
+    Real jlen = static_cast<Real>(domain.length(1));
+    Real klen = static_cast<Real>(domain.length(2));
+
+    // Center of the toroidal quadrilateral
+    Real cx = 0.5 * (problo[0]+probhi[0]);
+    Real cy = 0.5 * (problo[1]+probhi[1]);
+    Real cz = 0.5 * (problo[2]+probhi[2]);
+
+    // This is "delta eta" which is in the j direction
+    Real d_eta = 1. / jlen;
+
+    // This is "delta xi" which is in the i direction
+    Real d_xi = 1. / ilen;
+
+    // This is "delta zeta" which is in the k direction
+    Real d_zeta = 1. / klen;
+    // loc_arr is nodal so no offset
+    for(MFIter mfi(a_xyz_loc); mfi.isValid(); ++mfi)
+    {
+        const Box& gtbx = mfi.growntilebox();
+        amrex::Print()  << " grown tile box : " << gtbx << "\n";
+        auto loc_arr = a_xyz_loc.array(mfi);
+
+        ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            Real xi     = (static_cast<int>(i)) * d_xi;
+            Real eta    = (static_cast<int>(j)) * d_eta;
+            Real zeta   = (static_cast<int>(k)) * d_zeta;
+           
+            loc_arr(i,j,k,0) = cx + (0.1 * xi + 0.1) * cos(tpi * eta);
+            loc_arr(i,j,k,1) = cy + (0.1 * xi + 0.1) * sin(tpi * eta);
+            //loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); 3D torus
+            loc_arr(i,j,k,2) = zeta*0.02 + 0.02*tpi*eta;// ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi) * tpi*eta;
+            amrex::Real rad = std::sqrt((loc_arr(i,j,k,0)-0.5)*(loc_arr(i,j,k,0)-0.5) + (loc_arr(i,j,k,1)-0.5)*(loc_arr(i,j,k,1)-0.5)) ;
+        });
+
+    }
+    a_xyz_loc.FillBoundary(geom.periodicity());
+
+}
diff --git a/Tests/Particles/Mapped/InitTorus.cpp b/Tests/Particles/Mapped/InitTorus.cpp
index d4b032516ad..fdbc147440c 100644
--- a/Tests/Particles/Mapped/InitTorus.cpp
+++ b/Tests/Particles/Mapped/InitTorus.cpp
@@ -44,8 +44,7 @@ InitTorus (amrex::MultiFab& a_xyz_loc, amrex::Geometry& geom)
            
             loc_arr(i,j,k,0) = cx + (0.1 * xi + 0.1) * cos(tpi * eta);
             loc_arr(i,j,k,1) = cy + (0.1 * xi + 0.1) * sin(tpi * eta);
-            //loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); 3D torus
-            loc_arr(i,j,k,2) = cz + zeta*0.2 - 0.1 + 0.3;
+            loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); //3D torus
             amrex::Real rad = std::sqrt((loc_arr(i,j,k,0)-0.5)*(loc_arr(i,j,k,0)-0.5) + (loc_arr(i,j,k,1)-0.5)*(loc_arr(i,j,k,1)-0.5)) ;
         });
 
diff --git a/Tests/Particles/Mapped/InitUND.cpp b/Tests/Particles/Mapped/InitUND.cpp
index 8aee64eb27a..7bbbbedd66b 100644
--- a/Tests/Particles/Mapped/InitUND.cpp
+++ b/Tests/Particles/Mapped/InitUND.cpp
@@ -4,7 +4,7 @@
 using namespace amrex;
 
 enum struct ProbType {
-    Torus, Annulus, Stretched, Hill
+    Helix, Torus, Annulus, Stretched, Hill
 };
 
 void
@@ -18,11 +18,50 @@ InitUND_map (MultiFab& und, const MultiFab& a_xyz_loc, Geometry& geom, int flow_
     // Center of the annulus
     Real cx = 0.5 * (problo[0]+probhi[0]);
     Real cy = 0.5 * (problo[1]+probhi[1]);
-
+    const Real tpi = 2.0* amrex::Math::pi<Real>();
     //
     // ANNULUS
     //
-    if (prob_type == ProbType::Torus) {
+    if (prob_type == ProbType::Helix) {
+        // We only do this problem with fully mapped coordinates
+        AMREX_ALWAYS_ASSERT(a_xyz_loc.nComp() == AMREX_SPACEDIM);
+        for (MFIter mfi(und); mfi.isValid(); ++mfi)
+        {
+            const Box& tile_box = mfi.growntilebox();
+            auto loc_arr = a_xyz_loc.const_array(mfi);
+            auto und_arr = und.array(mfi);
+
+            ParallelFor( tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                // Physical location of cell center
+                Real x = loc_arr(i,j,k,0);
+                Real y = loc_arr(i,j,k,1);
+                Real theta;
+                //if (x == cx) {
+                //   theta = 3.14/2.;
+                //} else {
+                //   theta = atan((y-cy)/(x-cx));
+                //}
+                Real    rad = sqrt( (x-cx)*( x-cx) + (y-cy)*(y-cy));
+                theta = std::atan2((y-cy),(x-cx));
+                und_arr(i,j,k,0) = -tpi* rad * std::sin(theta); // rad*sin(theta);
+                und_arr(i,j,k,1) = tpi * rad * std::cos(theta); // -1.*rad*cos(theta);
+
+#if (AMREX_SPACEDIM == 3)
+                // Real z = loc_arr(i,j,k,2);
+                und_arr(i,j,k,2) =  0.02*tpi;
+#endif
+                if (i == 20 && j==5 && k==5) amrex::Print() << "UND AT " <<  IntVect(AMREX_D_DECL(i,j,k)) << " "
+                                                        << RealVect(AMREX_D_DECL(und_arr(i,j,k,0),und_arr(i,j,k,1),und_arr(i,j,k,2)))
+                                                        << std::endl;
+                if (i == 20 && j==23 && k==5) amrex::Print() << "UND AT " <<  IntVect(AMREX_D_DECL(i,j,k)) << " "
+                                                        << RealVect(AMREX_D_DECL(und_arr(i,j,k,0),und_arr(i,j,k,1),und_arr(i,j,k,2)))
+                                                        << std::endl;
+
+
+            });
+        }
+    } else if (prob_type == ProbType::Torus) {
         // We only do this problem with fully mapped coordinates
         AMREX_ALWAYS_ASSERT(a_xyz_loc.nComp() == AMREX_SPACEDIM);
         for (MFIter mfi(und); mfi.isValid(); ++mfi)
diff --git a/Tests/Particles/Mapped/Make.package b/Tests/Particles/Mapped/Make.package
index 5ec86d618c3..14c2565aa64 100644
--- a/Tests/Particles/Mapped/Make.package
+++ b/Tests/Particles/Mapped/Make.package
@@ -1,6 +1,7 @@
 CEXE_sources += main.cpp
 
 CEXE_sources += InitTorus.cpp
+CEXE_sources += InitHelix.cpp
 CEXE_sources += InitAnnulus.cpp
 CEXE_sources += InitHill.cpp
 CEXE_sources += InitStretched.cpp
diff --git a/Tests/Particles/Mapped/MappedPC.cpp b/Tests/Particles/Mapped/MappedPC.cpp
index b260fd3ac51..5ce0a0677d8 100644
--- a/Tests/Particles/Mapped/MappedPC.cpp
+++ b/Tests/Particles/Mapped/MappedPC.cpp
@@ -43,7 +43,7 @@ InitParticles (MultiFab& a_xyz_loc)
 #elif (AMREX_SPACEDIM == 3)
             int k = iv[2];
             //if (iv[0] == 0 && iv[1] == 0 && iv[2] >= dom_lo.z && iv[2] <= dom_hi.z/2) { // FOR EVERYTHING BUT TORUS
-            if (  ( (iv[0] > 1) && (iv[0] < 30) ) && (iv[1] == 0 ) && (iv[2]>1 && iv[2] < 10 )) {  // FOR TORUS
+            if (  ( (iv[0] > 4) && (iv[0] < 28) ) && (iv[1] == 1 ) && (iv[2]>1 && iv[2] < 50 )) {  // FOR TORUS
 #endif
                 int i = iv[0];
                 int j = iv[1];
diff --git a/Tests/Particles/Mapped/inputs b/Tests/Particles/Mapped/inputs
index af823d1936c..43ca1dc0e27 100644
--- a/Tests/Particles/Mapped/inputs
+++ b/Tests/Particles/Mapped/inputs
@@ -2,11 +2,11 @@ size = (32,32,64 )
 max_grid_size = 128
 is_periodic = (1,1,1)
 num_ppc = 1
-nsteps = 750
+nsteps = 1
 nlevs = 1
 
 grid_type = mapped
 vel_type  = nd
-prob_type = torus
+prob_type = helix
 
 do_plotfile = 1
diff --git a/Tests/Particles/Mapped/main.cpp b/Tests/Particles/Mapped/main.cpp
index d7a17f5c126..86efb22c62f 100644
--- a/Tests/Particles/Mapped/main.cpp
+++ b/Tests/Particles/Mapped/main.cpp
@@ -11,12 +11,37 @@
 
 using namespace amrex;
 
+void
+WriteMultiLevelPlotfileWithMapped (const std::string& plotfilename, int nlevels,
+                                         const Vector<const MultiFab*>& mf,
+                                         const Vector<const MultiFab*>& mf_nd,
+                                         const Vector<std::string>& varnames,
+                                         Real time,
+                                         const Vector<int>& level_steps,
+                                         const std::string &versionName,
+                                         const std::string &levelPrefix,
+                                         const std::string &mfPrefix,
+                                         const Vector<std::string>& extra_dirs,
+                                         const Vector<Geometry> geom);
+
+
+void
+WriteGenericPlotfileHeaderWithMapped (std::ostream &HeaderFile,
+                                            int nlevels,
+                                            const Vector<BoxArray> &bArray,
+                                            const Vector<std::string> &varnames,
+                                            Real time,
+                                            const Vector<int> &level_steps,
+                                            const std::string &versionName,
+                                            const std::string &levelPrefix,
+                                            const std::string &mfPrefix,
+                                            const Vector<Geometry> geom);
 enum struct GridType {
     Regular, Terrain, Mapped
 };
 
 enum struct ProbType {
-    Torus, Annulus, Stretched, Unstretched, Hill
+    Helix, Torus, Annulus, Stretched, Unstretched, Hill
 };
 
 enum struct VelType {
@@ -75,6 +100,7 @@ struct TestParams
 
 void Test ();
 
+void InitHelix       (MultiFab& a_xyz_loc, Geometry& geom);
 void InitTorus       (MultiFab& a_xyz_loc, Geometry& geom);
 void InitAnnulus     (MultiFab& a_xyz_loc, Geometry& geom);
 void InitUnstretched (MultiFab& a_xyz_loc  , Geometry& geom);
@@ -123,11 +149,13 @@ void get_test_params(TestParams& params)
 
     std::string prob_type_string;
     pp.get("prob_type", prob_type_string);
-    AMREX_ALWAYS_ASSERT(prob_type_string == "torus"       ||
+    AMREX_ALWAYS_ASSERT(prob_type_string == "helix"       ||
+                        prob_type_string == "torus"     ||
                         prob_type_string == "annulus"     ||
                         prob_type_string == "stretched"   ||
                         prob_type_string == "unstretched" ||
                         prob_type_string == "hill");
+    if (prob_type_string == "helix"      ) params.prob_type = ProbType::Helix;
     if (prob_type_string == "torus"      ) params.prob_type = ProbType::Torus;
     if (prob_type_string == "annulus"    ) params.prob_type = ProbType::Annulus;
     if (prob_type_string == "unstretched") params.prob_type = ProbType::Unstretched;
@@ -159,6 +187,189 @@ int main (int argc, char* argv[])
     amrex::Finalize();
 }
 
+void
+WriteMultiLevelPlotfileWithMapped (const std::string& plotfilename, int nlevels,
+                                         const Vector<const MultiFab*>& mf,
+                                         const Vector<const MultiFab*>& mf_nd,
+                                         const Vector<std::string>& varnames,
+                                         Real time,
+                                         const Vector<int>& level_steps,
+                                         const std::string &versionName,
+                                         const std::string &levelPrefix,
+                                         const std::string &mfPrefix,
+                                         const Vector<std::string>& extra_dirs,
+                                         const Vector<Geometry> geom)
+{
+    BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()");
+    amrex::Vector<IntVect> ref_ratio(nlevels);
+    for (int i = 0 ; i <nlevels; ++i) {
+        ref_ratio[i] = amrex::IntVect::TheUnitVector();
+    }
+    BL_ASSERT(nlevels <= mf.size());
+    BL_ASSERT(nlevels <= ref_ratio.size()+1);
+    BL_ASSERT(nlevels <= level_steps.size());
+    BL_ASSERT(mf[0]->nComp() == varnames.size());
+    int finest_level = 0;
+    bool callBarrier(false);
+    PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier);
+    if (!extra_dirs.empty()) {
+        for (const auto& d : extra_dirs) {
+            const std::string ed = plotfilename+"/"+d;
+            PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier);
+        }
+    }
+    ParallelDescriptor::Barrier();
+
+    if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) {
+        Vector<BoxArray> boxArrays(nlevels);
+        for(int level(0); level < boxArrays.size(); ++level) {
+            boxArrays[level] = mf[level]->boxArray();
+        }
+
+        auto f = [=]() {
+            VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
+            std::string HeaderFileName(plotfilename + "/Header");
+            std::ofstream HeaderFile;
+            HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
+            HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out   |
+                                                    std::ofstream::trunc |
+                                                    std::ofstream::binary);
+            if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName);
+            WriteGenericPlotfileHeaderWithMapped(HeaderFile, nlevels, boxArrays, varnames,
+                                                  time, level_steps, versionName,
+                                                  levelPrefix, mfPrefix, geom);
+        };
+
+        if (AsyncOut::UseAsyncOut()) {
+            AsyncOut::Submit(std::move(f));
+        } else {
+            f();
+        }
+    }
+
+    std::string mf_nodal_prefix = "Nu_nd";
+    for (int level = 0; level <= finest_level; ++level)
+    {
+        if (AsyncOut::UseAsyncOut()) {
+            VisMF::AsyncWrite(*mf[level],
+                              MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix),
+                              true);
+            VisMF::AsyncWrite(*mf_nd[level],
+                              MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix),
+                              true);
+        } else {
+            const MultiFab* data;
+            std::unique_ptr<MultiFab> mf_tmp;
+            if (mf[level]->nGrowVect() != 0) {
+                mf_tmp = std::make_unique<MultiFab>(mf[level]->boxArray(),
+                                                    mf[level]->DistributionMap(),
+                                                    mf[level]->nComp(), 0, MFInfo(),
+                                                    mf[level]->Factory());
+                MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0);
+                data = mf_tmp.get();
+            } else {
+                data = mf[level];
+            }
+            VisMF::Write(*data        , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
+            VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix));
+        }
+    }
+}
+
+void
+WriteGenericPlotfileHeaderWithMapped (std::ostream &HeaderFile,
+                                            int nlevels,
+                                            const Vector<BoxArray> &bArray,
+                                            const Vector<std::string> &varnames,
+                                            Real time,
+                                            const Vector<int> &level_steps,
+                                            const std::string &versionName,
+                                            const std::string &levelPrefix,
+                                            const std::string &mfPrefix,
+                                            const Vector<Geometry> geom)
+{
+        amrex::Vector<IntVect> ref_ratio(nlevels);
+        for (int i = 0 ; i <nlevels; ++i) {
+            ref_ratio[i] = amrex::IntVect::TheUnitVector();
+        } 
+        BL_ASSERT(nlevels <= bArray.size());
+        BL_ASSERT(nlevels <= ref_ratio.size()+1);
+        BL_ASSERT(nlevels <= level_steps.size());
+        int finest_level = 0;
+        HeaderFile.precision(17);
+
+        // ---- this is the generic plot file type name
+        HeaderFile << versionName << '\n';
+
+        HeaderFile << varnames.size() << '\n';
+
+        for (int ivar = 0; ivar < varnames.size(); ++ivar) {
+            HeaderFile << varnames[ivar] << "\n";
+        }
+        HeaderFile << AMREX_SPACEDIM << '\n';
+        HeaderFile << time << '\n';
+        HeaderFile << finest_level << '\n';
+        for (int i = 0; i < AMREX_SPACEDIM; ++i) {
+            HeaderFile << geom[0].ProbLo(i) << ' ';
+        }
+        HeaderFile << '\n';
+        for (int i = 0; i < AMREX_SPACEDIM; ++i) {
+            HeaderFile << geom[0].ProbHi(i) << ' ';
+        }
+        HeaderFile << '\n';
+        for (int i = 0; i < finest_level; ++i) {
+            HeaderFile << ref_ratio[i][0] << ' ';
+        }
+        HeaderFile << '\n';
+        for (int i = 0; i <= finest_level; ++i) {
+            HeaderFile << geom[i].Domain() << ' ';
+        }
+        HeaderFile << '\n';
+        for (int i = 0; i <= finest_level; ++i) {
+            HeaderFile << level_steps[i] << ' ';
+        }
+        HeaderFile << '\n';
+        for (int i = 0; i <= finest_level; ++i) {
+            for (int k = 0; k < AMREX_SPACEDIM; ++k) {
+                HeaderFile << geom[i].CellSize()[k] << ' ';
+            }
+            HeaderFile << '\n';
+        }
+        HeaderFile << (int) geom[0].Coord() << '\n';
+        HeaderFile << "0\n";
+
+        for (int level = 0; level <= finest_level; ++level) {
+            HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n';
+            HeaderFile << level_steps[level] << '\n';
+
+            const IntVect& domain_lo = geom[level].Domain().smallEnd();
+            for (int i = 0; i < bArray[level].size(); ++i)
+            {
+                // Need to shift because the RealBox ctor we call takes the
+                // physical location of index (0,0,0).  This does not affect
+                // the usual cases where the domain index starts with 0.
+                const Box& b = shift(bArray[level][i], -domain_lo);
+                RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
+                for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+                    HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
+                }
+            }
+
+            HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n';
+        }
+        HeaderFile << "1" << "\n";
+        HeaderFile << "3" << "\n";
+        HeaderFile << "amrexvec_nu_x" << "\n";
+        HeaderFile << "amrexvec_nu_y" << "\n";
+        HeaderFile << "amrexvec_nu_z" << "\n";
+        std::string mf_nodal_prefix = "Nu_nd";
+        for (int level = 0; level <= finest_level; ++level) {
+            HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n';
+        }
+}
+
+
+
 void Test()
 {
     BL_PROFILE("Test");
@@ -223,7 +434,11 @@ void Test()
     MultiFab a_xyz_loc(ba_nd,dm[lev],AMREX_SPACEDIM,2);
 
     // Annulus
-    if (params.prob_type == ProbType::Torus) {
+    if (params.prob_type == ProbType::Helix) {
+        AMREX_ALWAYS_ASSERT(params.grid_type == GridType::Mapped);
+        AMREX_ALWAYS_ASSERT(AMREX_SPACEDIM==3);
+        InitHelix(a_xyz_loc, geom[lev]);
+    } else if (params.prob_type == ProbType::Torus) {
         AMREX_ALWAYS_ASSERT(params.grid_type == GridType::Mapped);
         AMREX_ALWAYS_ASSERT(AMREX_SPACEDIM==3);
         InitTorus(a_xyz_loc, geom[lev]);
@@ -416,11 +631,12 @@ void Test()
         }
 
 
-    if (nt%20 ==0){
+    if (nt%2 ==0){
+        amrex::Print() << " writing plotfile \n";
         plotfilename = Concatenate("plt", nt, 5);
-        Vector<std::string> varname = {"ux", "uy"};
+        Vector<std::string> varname = {"ux", "uy", "uz"};
         amrex::MultiFab plotmf(ba[0], dm[0], varname.size(), 0 );
-        amrex::average_node_to_cellcenter (plotmf, 0, und, 0, 2, 0);
+        amrex::average_node_to_cellcenter (plotmf, 0, und, 0, 3, 0);
         // if (params.grid_type == GridType::Terrain) {
         //     terrain_pc.WritePlotFile(plotfilename, "particles");
         // } else if (params.grid_type == GridType::Mapped) {
@@ -428,9 +644,48 @@ void Test()
         // } else if (params.grid_type == GridType::Regular) {
         //     regular_pc.WritePlotFile(plotfilename, "particles");
         // }
-        WriteSingleLevelPlotfile(plotfilename, plotmf, varname, geom[0],0.0,0);
-        //WriteSingleLevelPlotfile("plt_grid",a_xyz_loc,{"gridmap",AMREX_D_DECL("x1","y1","z1")},geom,0.0,0);
+        // Make a nodal multifab that stores the location of the grid as x,y,z (3 components)
+        amrex::Vector<amrex::MultiFab> mf_nd(1);
+        if (params.grid_type == GridType::Mapped) {
+            amrex::BoxArray nodal_ba(ba[0]);
+            nodal_ba.surroundingNodes();
+            mf_nd[0].define(nodal_ba,dm[0],3,0);
+            mf_nd[0].setVal(0.);
+            amrex::Print() << nodal_ba[0] << "\n";
+            //Copy data from xyz location mf (copy 3 components)
+            amrex::MultiFab::Copy(mf_nd[0],a_xyz_loc,0,0,3,0);
+            amrex::Real dx = geom[0].CellSizeArray()[0];
+            amrex::Real dy = geom[0].CellSizeArray()[1];
+            amrex::Real dz = geom[0].CellSizeArray()[2];
+            for (amrex::MFIter mfi(mf_nd[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+                const amrex::Box & bx = mfi.tilebox();
+                amrex::Array4<amrex::Real> mf_arr = mf_nd[0].array(mfi);
+                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    mf_arr(i,j,k,0) -= i*dx;
+                    mf_arr(i,j,k,1) -= j*dy;
+                    mf_arr(i,j,k,2) -= k*dz;
+                });
+            }
+        }
+
+
+        const std::string versionName = "HyperCLaw-V1.1";
+        const std::string levelPrefix = "Level_";
+        const std::string mfPrefix = "Cell";
+        const amrex::Vector<std::string> extra_dirs = amrex::Vector<std::string>();
+        Vector<const MultiFab*> mfarr(1,&plotmf);
+        Vector<Geometry> geomarr(1,geom[0]);
+        Vector<int> level_steps(1,nt);
+        Vector<IntVect> ref_ratio;
+        ///WriteSingleLevelPlotfile(plotfilename, plotmf, varname, geom[0],0.0,0);
+        WriteMultiLevelPlotfileWithMapped(plotfilename, 1,
+                                           GetVecOfConstPtrs(mfarr),
+                                           GetVecOfConstPtrs(mf_nd),
+                                           varname,
+                                           nt, level_steps, versionName, levelPrefix, mfPrefix, extra_dirs, geomarr); 
         mapped_pc.WritePlotFile(plotfilename,"particles");
+
+        //WriteSingleLevelPlotfile("plt_grid",a_xyz_loc,{"gridmap",AMREX_D_DECL("x1","y1","z1")},geom,0.0,0);
         }
     } // nt
 //    plotfilename = Concatenate("plt", nt, 5);