diff --git a/amr-wind/utilities/sampling/CMakeLists.txt b/amr-wind/utilities/sampling/CMakeLists.txt
index 8c2ac5366b..f5de390008 100644
--- a/amr-wind/utilities/sampling/CMakeLists.txt
+++ b/amr-wind/utilities/sampling/CMakeLists.txt
@@ -13,7 +13,7 @@ target_sources(${amr_wind_lib_name}
   FieldNorms.cpp
   KineticEnergy.cpp
   Enstrophy.cpp
-  FreeSurface.cpp
+  FreeSurfaceSampler.cpp
   VolumeSampler.cpp
   WaveEnergy.cpp
   )
diff --git a/amr-wind/utilities/sampling/FreeSurface.H b/amr-wind/utilities/sampling/FreeSurface.H
deleted file mode 100644
index 27106b50e8..0000000000
--- a/amr-wind/utilities/sampling/FreeSurface.H
+++ /dev/null
@@ -1,133 +0,0 @@
-#ifndef FREESURFACE_H
-#define FREESURFACE_H
-
-#include <memory>
-
-#include "amr-wind/CFDSim.H"
-#include "amr-wind/utilities/PostProcessing.H"
-
-/**
- *  \defgroup findinterface Multiphase-sampling utilities
- *  Multiphase-sampling utilities
- *
- *  findinterface uses level_mask to find the location of the free surface
- *  between liquid and gas phases, given a user-defined 2D grid. It supports
- *  output in ascii as well as NetCDF format.
- *
- *  \ingroup utilities
- */
-
-namespace amr_wind::free_surface {
-
-/** Collection of data sampling objects
- *  \ingroup findinterface
- *
- *  A concrete implementation of the post-processing interface that deals with
- *  data probes. Sampling positions are defined with a plane and a mesh-aligned
- *  orientation.
- */
-class FreeSurface : public PostProcessBase::Register<FreeSurface>
-{
-public:
-    static std::string identifier() { return "FreeSurface"; }
-
-    FreeSurface(CFDSim& /*sim*/, std::string /*label*/);
-
-    ~FreeSurface() override;
-
-    //! Perform actions before mesh is created
-    void pre_init_actions() override {}
-
-    //! Read user inputs and create the different data probe instances
-    void initialize() override;
-
-    //! Interpolate fields at a given timestep and output to disk
-    void post_advance_work() override;
-
-    //! Redo some of the initialization work when the grid changes
-    void post_regrid_actions() override;
-
-    //! Output functions for private variables
-    //! Number of points
-    int num_gridpoints() const { return m_npts; }
-    //! Number of instances
-    int num_instances() const { return m_ninst; }
-    //! Locations
-    amrex::Vector<amrex::Array<amrex::Real, 2>> locations() const
-    {
-        return m_locs;
-    }
-    //! Outputs (heights)
-    amrex::Vector<amrex::Real> heights() const { return m_out; }
-
-protected:
-    //! Output data based on user-defined format
-    virtual void process_output();
-
-    //! Prepare NetCDF metadata
-    virtual void prepare_netcdf_file();
-
-    //! Write sampled data into a NetCDF file
-    void write_netcdf();
-
-    /** Output sampled data in ASCII format
-     *
-     *  Note that this should be used for debugging only and not in production
-     *  runs as it can have significant impacts on code performance.
-     */
-    void write_ascii();
-
-private:
-    CFDSim& m_sim;
-
-    /** Name of this sampling object.
-     *
-     *  The label is used to read user inputs from file and is also used for
-     *  naming files directories depending on the output format.
-     */
-    const std::string m_label;
-
-    //! reference to VOF
-    const Field& m_vof;
-
-    //! Format of the data output (ascii, netcdf, etc.)
-#ifdef AMR_WIND_USE_NETCDF
-    std::string m_out_fmt{"netcdf"};
-    std::string m_ncfile_name;
-#else
-    std::string m_out_fmt{"ascii"};
-#endif
-
-    //! Number or points on 2D grid in each direction
-    amrex::Vector<int> m_npts_dir;
-    int m_npts{0};
-    //! Number of instances (possible sampling points per location, like in the
-    //! case of a breaking wave)
-    int m_ninst{1};
-
-    //! Coordinate direction to search along, default is z
-    //! (this also determines the meaning of start and end points)
-    int m_coorddir{2};
-    //! Grid coordinates, determined as a function of m_coorddir
-    int m_gc0 = 0;
-    int m_gc1 = 1;
-
-    //! Parameters to set up plane
-    amrex::Vector<amrex::Real> m_start, m_end;
-    //! Locations of points in 2D grid
-    amrex::Vector<amrex::Array<amrex::Real, 2>> m_locs;
-    //! Output coordinate
-    amrex::Vector<amrex::Real> m_out;
-
-    //! Frequency of data sampling and output
-    int m_out_freq{100};
-
-    //! Max number of sample points found in a single cell
-    int m_ncomp{1};
-    //! Max number of sample points allowed in a single cell
-    int m_ncmax{8};
-};
-
-} // namespace amr_wind::free_surface
-
-#endif /* FREESURFACE_H */
diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.H b/amr-wind/utilities/sampling/FreeSurfaceSampler.H
new file mode 100644
index 0000000000..63fe19efd1
--- /dev/null
+++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.H
@@ -0,0 +1,114 @@
+#ifndef FREESURFACESAMPLER_H
+#define FREESURFACESAMPLER_H
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/utilities/sampling/SamplerBase.H"
+
+namespace amr_wind::sampling {
+
+class FreeSurfaceSampler : public SamplerBase::Register<FreeSurfaceSampler>
+{
+public:
+    static std::string identifier() { return "FreeSurfaceSampler"; }
+
+    FreeSurfaceSampler(CFDSim& /*sim*/);
+
+    ~FreeSurfaceSampler() override;
+
+    //! Read user inputs and create mesh-based information for locating
+    //! interface
+    void initialize(const std::string& key) override;
+
+    //! Populate and return a vector of probe locations to be sampled
+    void sampling_locations(SampleLocType& /*locs*/) const override;
+    void output_locations(SampleLocType& locs) const override
+    {
+        return sampling_locations(locs);
+    }
+
+    //! Find heights associated with 2D sample locations
+    void update_sampling_locations() override;
+
+    //! Redo some of the initialization work when the grid changes
+    void post_regrid_actions() override;
+
+    void
+    define_netcdf_metadata(const ncutils::NCGroup& /*unused*/) const override;
+    void
+    populate_netcdf_metadata(const ncutils::NCGroup& /*unused*/) const override;
+    void output_netcdf_data(
+        const ncutils::NCGroup& /*unused*/,
+        const size_t /*unused*/) const override;
+
+    //! Name of this sampling object
+    std::string label() const override { return m_label; }
+    std::string& label() override { return m_label; }
+
+    //! Type of this sampling object
+    std::string sampletype() const override { return identifier(); }
+
+    //! Unique identifier for this set of probe locations
+    int id() const override { return m_id; }
+    int& id() override { return m_id; }
+
+    //! Output functions for private variables
+    //! Number of 2D grid points
+    int num_gridpoints() const { return m_npts; }
+    //! Number of instances
+    int num_instances() const { return m_ninst; }
+    //! Number of points, total
+    long num_points() const override
+    {
+        return static_cast<long>(m_npts) * m_ninst;
+    }
+    long num_output_points() const override
+    {
+        return static_cast<long>(m_npts) * m_ninst;
+    }
+    //! Locations
+    amrex::Vector<amrex::Array<amrex::Real, 2>> grid_locations() const
+    {
+        return m_grid_locs;
+    }
+    //! Outputs (heights)
+    amrex::Vector<amrex::Real> heights() const { return m_out; }
+
+private:
+    CFDSim& m_sim;
+
+    //! reference to VOF
+    const Field& m_vof;
+
+    //! Number or points on 2D grid in each direction
+    amrex::Vector<int> m_npts_dir;
+    int m_npts{0};
+    //! Number of instances (possible sampling points per location, like in the
+    //! case of a breaking wave)
+    int m_ninst{1};
+
+    //! Coordinate direction to search along, default is z
+    //! (this also determines the meaning of start and end points)
+    int m_coorddir{2};
+    //! Grid coordinates, determined as a function of m_coorddir
+    int m_gc0 = 0;
+    int m_gc1 = 1;
+
+    //! Parameters to set up plane
+    amrex::Vector<amrex::Real> m_start, m_end;
+    //! Locations of points in 2D grid
+    amrex::Vector<amrex::Array<amrex::Real, 2>> m_grid_locs;
+    //! Output coordinate
+    amrex::Vector<amrex::Real> m_out;
+
+    //! Max number of sample points found in a single cell
+    int m_ncomp{1};
+    //! Max number of sample points allowed in a single cell
+    int m_ncmax{8};
+
+    std::string m_label;
+    int m_id{-1};
+};
+
+} // namespace amr_wind::sampling
+
+#endif /* FREESURFACESAMPLER_H */
diff --git a/amr-wind/utilities/sampling/FreeSurface.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
similarity index 82%
rename from amr-wind/utilities/sampling/FreeSurface.cpp
rename to amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
index 226a5d1f93..a2b1af5aca 100644
--- a/amr-wind/utilities/sampling/FreeSurface.cpp
+++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp
@@ -1,35 +1,30 @@
-#include "amr-wind/utilities/sampling/FreeSurface.H"
+#include "amr-wind/utilities/sampling/FreeSurfaceSampler.H"
 #include "amr-wind/utilities/io_utils.H"
 #include <AMReX_MultiFabUtil.H>
 #include <utility>
-#include "amr-wind/utilities/ncutils/nc_interface.H"
 #include "amr-wind/equation_systems/vof/volume_fractions.H"
-#include "amr-wind/utilities/IOManager.H"
 
 #include "AMReX_ParmParse.H"
 
-namespace amr_wind::free_surface {
+namespace amr_wind::sampling {
 
-FreeSurface::FreeSurface(CFDSim& sim, std::string label)
-    : m_sim(sim), m_label(std::move(label)), m_vof(sim.repo().get_field("vof"))
+FreeSurfaceSampler::FreeSurfaceSampler(CFDSim& sim)
+    : m_sim(sim), m_vof(sim.repo().get_field("vof"))
 {}
 
-FreeSurface::~FreeSurface() = default;
+FreeSurfaceSampler::~FreeSurfaceSampler() = default;
 
-void FreeSurface::initialize()
+void FreeSurfaceSampler::initialize(const std::string& key)
 {
-    BL_PROFILE("amr-wind::FreeSurface::initialize");
+    BL_PROFILE("amr-wind::FreeSurfaceSampler::initialize");
 
     {
-        amrex::ParmParse pp(m_label);
-        pp.query("output_frequency", m_out_freq);
-        pp.query("output_format", m_out_fmt);
-        // Load parameters of freesurface sampling
-        pp.query("num_instances", m_ninst);
+        amrex::ParmParse pp(key);
+        pp.getarr("plane_start", m_start);
+        pp.getarr("plane_end", m_end);
+        pp.getarr("plane_num_points", m_npts_dir);
         pp.query("search_direction", m_coorddir);
-        pp.getarr("num_points", m_npts_dir);
-        pp.getarr("start", m_start);
-        pp.getarr("end", m_end);
+        pp.query("num_instances", m_ninst);
         pp.query("max_sample_points_per_cell", m_ncmax);
         AMREX_ALWAYS_ASSERT(static_cast<int>(m_start.size()) == AMREX_SPACEDIM);
         AMREX_ALWAYS_ASSERT(static_cast<int>(m_end.size()) == AMREX_SPACEDIM);
@@ -53,7 +48,7 @@ void FreeSurface::initialize()
         }
         default: {
             amrex::Abort(
-                "FreeSurface: Invalid coordinate search direction "
+                "FreeSurfaceSampler: Invalid coordinate search direction "
                 "encountered");
             break;
         }
@@ -66,7 +61,7 @@ void FreeSurface::initialize()
     m_npts = m_npts_dir[0] * m_npts_dir[1];
 
     // Turn parameters into 2D grid
-    m_locs.resize(m_npts);
+    m_grid_locs.resize(m_npts);
     m_out.resize(static_cast<long>(m_npts) * m_ninst);
 
     // Get size of sample grid spacing
@@ -84,9 +79,9 @@ void FreeSurface::initialize()
                 m_out[idx * m_ninst + ni] = m_start[m_coorddir];
             }
             // Grid direction 1
-            m_locs[idx][0] = m_start[m_gc0] + dxs0 * i;
+            m_grid_locs[idx][0] = m_start[m_gc0] + dxs0 * i;
             // Grid direction 2
-            m_locs[idx][1] = m_start[m_gc1] + dxs1 * j;
+            m_grid_locs[idx][1] = m_start[m_gc1] + dxs1 * j;
 
             ++idx;
         }
@@ -335,22 +330,34 @@ void FreeSurface::initialize()
                 });
         }
     }
+}
 
-    if (m_out_fmt == "netcdf") {
-        prepare_netcdf_file();
+void FreeSurfaceSampler::sampling_locations(SampleLocType& locs) const
+{
+    locs.resize(num_output_points());
+
+    int idx = 0;
+    for (int j = 0; j < m_npts_dir[1]; ++j) {
+        for (int i = 0; i < m_npts_dir[0]; ++i) {
+            // Initialize output values to 0.0
+            for (int ni = 0; ni < m_ninst; ++ni) {
+                // Grid direction 1
+                locs[idx * m_ninst + ni][m_gc0] = m_grid_locs[idx][0];
+                // Grid direction 2
+                locs[idx * m_ninst + ni][m_gc1] = m_grid_locs[idx][1];
+                // Output direction
+                locs[idx * m_ninst + ni][m_coorddir] =
+                    m_out[idx * m_ninst + ni];
+            }
+            ++idx;
+        }
     }
 }
 
-void FreeSurface::post_advance_work()
+void FreeSurfaceSampler::update_sampling_locations()
 {
 
-    BL_PROFILE("amr-wind::FreeSurface::post_advance_work");
-    const auto& time = m_sim.time();
-    const int tidx = time.time_index();
-    // Skip processing if it is not an output timestep
-    if (!(tidx % m_out_freq == 0)) {
-        return;
-    }
+    BL_PROFILE("amr-wind::FreeSurfaceSampler::update_sampling_locations");
 
     // Zero data in output array
     for (int n = 0; n < m_npts * m_ninst; n++) {
@@ -518,13 +525,11 @@ void FreeSurface::post_advance_work()
             dout_ptr[n] = plo0[coorddir];
         });
     }
-
-    process_output();
 }
 
-void FreeSurface::post_regrid_actions()
+void FreeSurfaceSampler::post_regrid_actions()
 {
-    BL_PROFILE("amr-wind::FreeSurface::post_regrid_actions");
+    BL_PROFILE("amr-wind::FreeSurfaceSampler::post_regrid_actions");
     // Small number for floating-point comparisons
     constexpr amrex::Real eps = 1.0e-16;
     // Get working fields
@@ -665,175 +670,45 @@ void FreeSurface::post_regrid_actions()
     }
 }
 
-void FreeSurface::process_output()
-{
-    if (m_out_fmt == "ascii") {
-        write_ascii();
-    } else if (m_out_fmt == "netcdf") {
-        write_netcdf();
-    } else {
-        amrex::Abort("FreeSurface: Invalid output format encountered");
-    }
-}
-
-void FreeSurface::write_ascii()
-{
-    BL_PROFILE("amr-wind::FreeSurface::write_ascii");
-    amrex::Print()
-        << "WARNING: FreeSurface: ASCII output will impact performance"
-        << std::endl;
-
-    const std::string post_dir = m_sim.io_manager().post_processing_directory();
-    const std::string sname =
-        amrex::Concatenate(m_label, m_sim.time().time_index());
-
-    const std::string fname = post_dir + "/" + sname + ".txt";
-
-    if (amrex::ParallelDescriptor::IOProcessor()) {
-        //
-        // Have I/O processor open file and write everything.
-        //
-        std::ofstream File;
-
-        File.open(fname.c_str(), std::ios::out | std::ios::trunc);
-
-        if (!File.good()) {
-            amrex::FileOpenFailed(fname);
-        }
-
-        std::string str1 = "x";
-        std::string str2 = "y";
-        switch (m_coorddir) {
-        case 0: {
-            str1 = "y";
-            str2 = "z";
-            break;
-        }
-        case 1: {
-            str1 = "x";
-            str2 = "z";
-            break;
-        }
-        default: {
-            // Do nothing, initial settings are good
-            break;
-        }
-        }
-
-        // Metadata
-        File << m_sim.time().new_time() << '\n';
-        File << m_npts << '\n';
-        File << str1 << ' ' << m_npts_dir[0] << ' ' << str2 << ' '
-             << m_npts_dir[1] << '\n';
-
-        // Points in grid (x, y, z0, z1, ...)
-        for (int n = 0; n < m_npts; ++n) {
-            File << m_locs[n][0] << ' ' << m_locs[n][1];
-            for (int ni = 0; ni < m_ninst; ++ni) {
-                File << ' ' << m_out[ni * m_npts + n];
-            }
-            File << '\n';
-        }
-
-        File.flush();
-
-        File.close();
-
-        if (!File.good()) {
-            amrex::Abort("FreeSurface::write_ascii(): problem writing file");
-        }
-    }
-}
-
-void FreeSurface::prepare_netcdf_file()
-{
 #ifdef AMR_WIND_USE_NETCDF
-
-    const std::string post_dir = m_sim.io_manager().post_processing_directory();
-    const std::string sname =
-        amrex::Concatenate(m_label, m_sim.time().time_index());
-
-    m_ncfile_name = post_dir + "/" + sname + ".nc";
-
-    // Only I/O processor handles NetCDF generation
-    if (!amrex::ParallelDescriptor::IOProcessor()) {
-        return;
-    }
-
-    auto ncf = ncutils::NCFile::create(m_ncfile_name, NC_CLOBBER | NC_NETCDF4);
-    const std::string nt_name = "num_time_steps";
-    const std::string ngp_name = "num_grid_points";
-    const std::string ninst_name = "num_instances";
-    ncf.enter_def_mode();
-    ncf.put_attr("title", "AMR-Wind data sampling output");
-    ncf.put_attr("version", ioutils::amr_wind_version());
-    ncf.put_attr("created_on", ioutils::timestamp());
-    ncf.def_dim(nt_name, NC_UNLIMITED);
-    ncf.def_dim(ngp_name, m_npts);
-    ncf.def_dim(ninst_name, m_ninst);
-    ncf.def_dim("ndim2", 2);
-    ncf.def_var("time", NC_DOUBLE, {nt_name});
+void FreeSurfaceSampler::define_netcdf_metadata(
+    const ncutils::NCGroup& grp) const
+{
 
     // Metadata related to the 2D grid used to sample
-    const std::vector<int> ijk{m_npts_dir[0], m_npts_dir[1]};
-    ncf.put_attr("ijk_dims", ijk);
-    ncf.put_attr("start", m_start);
-    ncf.put_attr("end", m_end);
-
-    // Set up array of data for locations in 2D grid
-    ncf.def_var("coordinates2D", NC_DOUBLE, {ngp_name, "ndim2"});
-
-    // Set up array for height outputs
-    ncf.def_var("heights", NC_DOUBLE, {nt_name, ninst_name, ngp_name});
-
-    ncf.exit_def_mode();
-
-    // Populate data that doesn't change
-    const std::vector<size_t> start{0, 0};
-    std::vector<size_t> count{0, 2};
-    count[0] = m_npts;
-    auto xy = ncf.var("coordinates2D");
-    xy.put(m_locs[0].data(), start, count);
-
-#else
-    amrex::Abort(
-        "NetCDF support was not enabled during build time. Please "
-        "recompile or "
-        "use native format");
-#endif
+    const std::vector<int> ijk{m_npts_dir[0], m_npts_dir[1], m_ninst};
+    grp.put_attr("sampling_type", identifier());
+    grp.put_attr("ijk_dims", ijk);
+    grp.put_attr("plane_start", m_start);
+    grp.put_attr("plane_end", m_end);
+    grp.def_var("points", NC_DOUBLE, {"num_time_steps", "num_points", "ndim"});
 }
 
-void FreeSurface::write_netcdf()
+void FreeSurfaceSampler::populate_netcdf_metadata(
+    const ncutils::NCGroup& /*unused*/) const
+{}
+void FreeSurfaceSampler::output_netcdf_data(
+    const ncutils::NCGroup& grp, const size_t nt) const
 {
-#ifdef AMR_WIND_USE_NETCDF
-
-    if (!amrex::ParallelDescriptor::IOProcessor()) {
-        return;
-    }
-    auto ncf = ncutils::NCFile::open(m_ncfile_name, NC_WRITE);
-    const std::string nt_name = "num_time_steps";
-    // Index of the next timestep
-    const size_t nt = ncf.dim(nt_name).len();
-    {
-        auto time = m_sim.time().new_time();
-        ncf.var("time").put(&time, {nt}, {1});
-    }
-
+    // Write the coordinates every time
     std::vector<size_t> start{nt, 0, 0};
-    std::vector<size_t> count{1, 0, 0};
-
-    count[1] = 1;
-    count[2] = m_npts;
-    auto var = ncf.var("heights");
-    for (int ni = 0; ni < m_ninst; ++ni) {
-        var.put(
-            &m_out[static_cast<size_t>(ni) * static_cast<size_t>(m_npts)],
-            start, count);
-        ++start[1];
-    }
-
-    ncf.close();
-#endif
+    std::vector<size_t> count{1, 0, AMREX_SPACEDIM};
+    SamplerBase::SampleLocType locs;
+    sampling_locations(locs);
+    auto xyz = grp.var("points");
+    count[1] = num_output_points();
+    xyz.put(locs[0].data(), start, count);
 }
+#else
+void FreeSurfaceSampler::define_netcdf_metadata(
+    const ncutils::NCGroup& /*unused*/) const
+{}
+void FreeSurfaceSampler::populate_netcdf_metadata(
+    const ncutils::NCGroup& /*unused*/) const
+{}
+void FreeSurfaceSampler::output_netcdf_data(
+    const ncutils::NCGroup& /*unused*/, const size_t /*unused*/) const
+{}
+#endif
 
-} // namespace amr_wind::free_surface
+} // namespace amr_wind::sampling
diff --git a/amr-wind/utilities/sampling/SamplerBase.H b/amr-wind/utilities/sampling/SamplerBase.H
index 440b031315..9151693462 100644
--- a/amr-wind/utilities/sampling/SamplerBase.H
+++ b/amr-wind/utilities/sampling/SamplerBase.H
@@ -61,6 +61,9 @@ public:
     //! Run actions after sample (useful in interpolated subsampling)
     virtual void post_sample_actions() {}
 
+    //! Run actions after regrid (important when using field-based quantities)
+    virtual void post_regrid_actions() {}
+
     //! Run specific output for the sampler
     virtual bool output_netcdf_field(
         const std::vector<double>& /*unused*/,
diff --git a/amr-wind/utilities/sampling/Sampling.cpp b/amr-wind/utilities/sampling/Sampling.cpp
index c7d7e70219..eaf23cee4a 100644
--- a/amr-wind/utilities/sampling/Sampling.cpp
+++ b/amr-wind/utilities/sampling/Sampling.cpp
@@ -193,6 +193,11 @@ void Sampling::sampling_post()
 void Sampling::post_regrid_actions()
 {
     BL_PROFILE("amr-wind::Sampling::post_regrid_actions");
+
+    for (const auto& obj : m_samplers) {
+        obj->post_regrid_actions();
+    }
+
     m_scontainer->Redistribute();
 }
 
diff --git a/docs/sphinx/user/inputs_Sampling.rst b/docs/sphinx/user/inputs_Sampling.rst
index 7dbd0f441e..6cf7295a9a 100644
--- a/docs/sphinx/user/inputs_Sampling.rst
+++ b/docs/sphinx/user/inputs_Sampling.rst
@@ -158,3 +158,29 @@ Example::
   sampling.volume1.hi        = 3.0 1.0 0.5
   sampling.volume1.lo      = 0.0 0.0 -0.5
   sampling.volume1.num_points  = 30 10 10
+
+Sampling on the air-water interface
+```````````````````````````````````
+
+The ``FreeSurfaceSampler`` samples on the air-water interface, and it requires the 
+vof (volume-of-fluid) field to be present in order to function. The sample locations
+are specified using a grid that starts at ``plane_start`` and
+extends to ``plane_end``. The resolution in each direction is specified by
+``plane_num_points``. The coordinates of the sampling
+locations are determined by the location of the air-water interface in the search
+direction, specified by ``search_direction``, and the other coordinates are 
+determined by the ``plane_`` parameters. The default search direction parameter
+is 2, indicating the samplers will search for the interface along the z-direction. 
+Due to this design, it is best to specify a plane that is normal to the intended 
+search direction. Another optional parameter is ``num_instances``, which is available
+for cases where the interface location is multivalued along the search direction,
+such as during wave breaking. This parameter defaults to 1, and the sampler will
+automatically select the highest position along the search direction when the interface
+location is multivalued.
+
+Example::
+
+  sampling.fs1.type             = FreeSurfaceSampler
+  sampling.fs1.plane_start      = 4.0 -1.0 0.0
+  sampling.fs1.plane_end        = 0.0 1.0  0.0
+  sampling.fs1.plane_num_points = 20 10
diff --git a/docs/sphinx/user/plot_sampling_freesurface.pdf b/docs/sphinx/user/plot_sampling_freesurface.pdf
new file mode 100644
index 0000000000..6079cde939
Binary files /dev/null and b/docs/sphinx/user/plot_sampling_freesurface.pdf differ
diff --git a/docs/sphinx/user/post_processing_examples.rst b/docs/sphinx/user/post_processing_examples.rst
new file mode 100644
index 0000000000..82145b107b
--- /dev/null
+++ b/docs/sphinx/user/post_processing_examples.rst
@@ -0,0 +1,24 @@
+Post-processing Examples
+========================
+
+This section provides examples of manipulating and plotting the
+post-processing data output by AMR-Wind (i.e., the files written
+to the post_processing directory, not checkpoint files or plotfiles).
+
+.. toctree::
+   :glob:
+   :maxdepth: 3
+
+   sampling.rst
+
+Simple Outputs
+--------------
+   
+The following types of post_processing outputs produce text files
+with each column labeled in a 1-line header, and due to their
+simplicity, no example scripts are provided:
+
+* FieldNorms
+* Enstrophy
+* KineticEnergy
+* WaveEnergy
diff --git a/docs/sphinx/user/sampling.rst b/docs/sphinx/user/sampling.rst
new file mode 100644
index 0000000000..e30b61f043
--- /dev/null
+++ b/docs/sphinx/user/sampling.rst
@@ -0,0 +1,20 @@
+.. _sampling:
+
+Sampling
+========
+
+AMR-Wind provides a sampling utility that measures field
+variables at specific locations within the flow. A variety of
+sampler types are available, which determine the locations to 
+be sampled and, depending on the type, update the sampling 
+locations as time progresses. The sampling data can be output
+in AMReX particle native format (``native``), NetCDF format 
+(``netcdf``), or AMReX particle ASCII format (``ascii``).
+The native and ascii formats are identical for every sampler,
+but the netcdf format can include additional output variables
+depending on the sampler type.
+
+.. toctree::
+   :maxdepth: 2
+
+   sampling_freesurfacesampler.rst
diff --git a/docs/sphinx/user/sampling_freesurfacesampler.rst b/docs/sphinx/user/sampling_freesurfacesampler.rst
new file mode 100644
index 0000000000..502bd04a8b
--- /dev/null
+++ b/docs/sphinx/user/sampling_freesurfacesampler.rst
@@ -0,0 +1,54 @@
+.. _sampling_freesurface_sampler.rst:
+
+FreeSurfaceSampler
+==================
+
+These examples plot the evolution of the liquid-gas interface as a liquid column 
+collapses and flows into the rest of the domain, and they show the velocity 
+magnitude (in color) as a function of interface position. At each instant plotted,
+the shading of the line color represents the progression of time. Below, the scripts
+for processing the data in each available format (``native``, ``netcdf``, and ``ascii``)
+are shown, and these scripts are also available in the ``tools`` directory.
+
+Note that the data accesses in the native and ascii examples can be applied to 
+other sampler types because the format is identical. However, the NetCDF format
+can vary depending on the sampler type, and further detail is provided in that example.
+
+The three examples will generate the same plot, shown here:
+
+.. image:: ./plot_sampling_freesurface.pdf
+   :width: 300pt
+
+|
+
+Native format: dam break example
+--------------------------------
+To generate the data required to replicate this example, run the simulation contained in 
+test/test_files/dam_break_godunov with the sampling format set to "native".
+
+.. literalinclude:: ../../../tools/sampling_dam_break_godunov_native.py
+   :language: python
+
+NetCDF format: dam break example
+--------------------------------
+To generate the data required to replicate this example, run the simulation contained in 
+test/test_files/dam_break_godunov with the sampling format set to "netcdf" (AMR-Wind must 
+also be compiled with NetCDF).
+
+Because the FreeSurfaceSampler tracks the interface location, the sample points change position 
+with time, and these changing positions are output as the "points" field. However, samplers that 
+do not have moving sample locations do not have the "points" field in the NetCDF dataset. For 
+static samplers, the locations are provided solely through the variable "coordinates", which represents
+the initial positions of the sampler points. Otherwise, the data accesses in this example can be 
+applied to other sampler types.
+
+.. literalinclude:: ../../../tools/sampling_dam_break_godunov_netcdf.py
+   :language: python
+
+ASCII format: dam break example
+--------------------------------
+To generate the data required to replicate this example, run the simulation contained in 
+test/test_files/dam_break_godunov with the sampling format set to "ascii".
+
+.. literalinclude:: ../../../tools/sampling_dam_break_godunov_ascii.py
+   :language: python
\ No newline at end of file
diff --git a/docs/sphinx/user/user.rst b/docs/sphinx/user/user.rst
index a874eeaeb8..af3c6c15a6 100644
--- a/docs/sphinx/user/user.rst
+++ b/docs/sphinx/user/user.rst
@@ -9,3 +9,4 @@ User Manual
    inputs.rst
    abl_bndry_io.rst
    run.rst
+   post_processing_examples.rst
diff --git a/test/test_files/dam_break_godunov/dam_break_godunov.inp b/test/test_files/dam_break_godunov/dam_break_godunov.inp
index ca6bae8937..a7358654f4 100644
--- a/test/test_files/dam_break_godunov/dam_break_godunov.inp
+++ b/test/test_files/dam_break_godunov/dam_break_godunov.inp
@@ -1,5 +1,5 @@
 time.stop_time               =   6     # Max (simulated) time to evolve
-time.max_step                =   20          # Max number of time steps
+time.max_step                =   150          # Max number of time steps
 
 time.fixed_dt         =   0.001        # Use this constant dt if > 0
 time.cfl              =   0.95         # CFL factor
@@ -7,6 +7,16 @@ time.cfl              =   0.95         # CFL factor
 time.plot_interval            =  10       # Steps between plot files
 time.checkpoint_interval      =  -100       # Steps between checkpoint files
 
+incflo.post_processing        = sampling
+sampling.output_frequency     = 30
+sampling.output_format        = native
+sampling.labels               = fs
+sampling.fields               = velocity
+sampling.fs.type              = FreeSurfaceSampler
+sampling.fs.plane_num_points  = 40 1
+sampling.fs.plane_start       = 0.0 0.0625 0.0
+sampling.fs.plane_end         = 0.5 0.0625 0.0
+
 io.output_default_variables = 0
 io.outputs = density vof
 io.derived_outputs = "components(velocity,0,2)" "components(gp,0,2)"
diff --git a/tools/amrex_particle.py b/tools/amrex_particle.py
index 0e5d7ee840..4bdb54a92f 100644
--- a/tools/amrex_particle.py
+++ b/tools/amrex_particle.py
@@ -20,14 +20,14 @@ class AmrexParticleFile:
     """
 
     @classmethod
-    def load(cls, time_index, label="sampling", root_dir="post_processing"):
+    def load(cls, time_index, label="sampling", root_dir="post_processing", suffix="/particles"):
         """
         Args:
             time_index (int): Time index to load
             label (str): Label used for this set of probes
             root_dir (path): Path to the post_processing directory
         """
-        fpath = Path(root_dir) / ("%s%05d"%(label, time_index))
+        fpath = Path(root_dir) / ("%s%05d%s"%(label, time_index, suffix))
         return cls(fpath)
 
     def __init__(self, pdir):
@@ -67,7 +67,7 @@ def parse_header(self):
             self.next_id = int(fh.readline().strip())
             self.finest_level = int(fh.readline().strip())
 
-            self.num_grids = np.empty((self.finest_level+1,), dtype=np.int)
+            self.num_grids = np.empty((self.finest_level+1,), dtype=int)
             self.grid_info = []
             for lev in range(self.finest_level+1):
                 ginfo = []
@@ -80,9 +80,9 @@ def parse_header(self):
     def load_binary_data(self):
         """Read binary data into memory"""
         self.real_data = np.empty((self.num_particles,
-                                   self.ndim+self.num_reals), dtype=np.float)
+                                   self.ndim+self.num_reals), dtype=float)
         self.int_data = np.empty((self.num_particles,
-                                  self.num_ints), dtype=np.int)
+                                  self.num_ints), dtype=int)
 
         idata = self.int_data
         rdata = self.real_data
@@ -100,7 +100,7 @@ def load_binary_data(self):
                 fh = bfiles[fstr]
                 fh.seek(offset)
                 ivals = np.fromfile(fh, dtype=np.int32, count=nints*npts)
-                rvals = np.fromfile(fh, dtype=np.float, count=nreals*npts)
+                rvals = np.fromfile(fh, dtype=float, count=nreals*npts)
 
                 for i, ii in enumerate(range(0, nints * npts, nints)):
                     pidx = ivals[ii + 2]
diff --git a/tools/freesurface_plots.py b/tools/freesurface_plots.py
deleted file mode 100644
index 41881b286d..0000000000
--- a/tools/freesurface_plots.py
+++ /dev/null
@@ -1,85 +0,0 @@
-# -*- coding: utf-8 -*-
-# pylint: disable=too-many-locals
-
-"""
-ABL Statistics plotting utility
-"""
-
-import numpy as np
-from matplotlib.backends.backend_pdf import PdfPages
-import matplotlib.pyplot as plt
-import netCDF4 as ncdf
-import pandas as pd
-import sys
-
-
-class FreeSurfaceStatsFile(object):
-    """Interface to ABL Statistics NetCDF file"""
-
-    def __init__(self, stats_file="FreeSurface.nc"):
-        """
-        Args:
-            stats_file (path): Absolute path to the NetCDF file
-        """
-        self.stats_file = stats_file
-        self.fs_stats = ncdf.Dataset(stats_file)
-        self._time = self.fs_stats["time"][:]
-        self._coords = self.fs_stats["coordinates2D"][:]
-        self._heights = self.fs_stats["heights"][:]
-        self.start = self.fs_stats.getncattr("start")
-        self.end = self.fs_stats.getncattr("end")
-        self.ninst = len(self.heights[0,:,0])
-        ijk_dims = self.fs_stats.getncattr("ijk_dims")
-        self.ni = ijk_dims[0]
-        self.nj = ijk_dims[1]
-
-    @property
-    def time(self):
-        """The time array from the FreeSurface file"""
-        return self._time
-
-    @property
-    def coords(self):
-        """The coordinates array where data is available"""
-        return self._coords
-
-    @property
-    def heights(self):
-        """Heights where data is available"""
-        return self._heights
-
-def plot_interface(stats, pdffile, num_step):
-    """Generate interface plot"""
-    ht_full = stats.heights
-    coords = stats.coords
-    time = stats.time
-    start = stats.start
-    end = stats.end
-    nx = stats.ni
-    ny = stats.nj
-    ninst = stats.ninst
-    print("Instant displayed = %.3f s"%(time[num_step]))
-    ht_now = ht_full[num_step, :, :]
-
-    # Plot 2D color plot for each instance that is available
-    for n in range(ninst):
-        plt.figure()
-        arr = np.reshape(ht_now[n,:],(nx, ny))
-        plt.imshow(arr, extent=[start[0],end[0],start[1],end[1]],
-            origin="lower",interpolation='bilinear')
-        plt.xlabel("x (m)")
-        plt.ylabel("y (m)")
-        plt.colorbar()
-        pdffile.savefig(dpi=300, bbox_inches="tight")
-        plt.close()
-
-if __name__ == "__main__":
-    clarg = sys.argv
-    if (len(clarg) != 2):
-        print("Program cancelled. Please input desired timestep for plotting.")
-        quit()
-    num_step = int(sys.argv[1])
-    print("Timestep: %d"%(num_step))
-    statsfile = FreeSurfaceStatsFile()
-    with PdfPages("interface_test.pdf") as pdf:
-        plot_interface(statsfile, pdf, num_step)
diff --git a/tools/sampling_dam_break_godunov_ascii.py b/tools/sampling_dam_break_godunov_ascii.py
new file mode 100644
index 0000000000..6cd22b9d42
--- /dev/null
+++ b/tools/sampling_dam_break_godunov_ascii.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+loc_dir = "."
+pp_dir = loc_dir + "/post_processing"
+
+nt = 5 + 1
+out_int = 30
+
+plt.figure()
+base_color = 0.8
+
+for n in range(nt):
+    x_oo=np.genfromtxt(pp_dir+"/sampling"+str(n*out_int).zfill(5)+".txt",delimiter=' ',skip_header=5)
+    # Reorder array
+    ind_0 = np.argsort(x_oo[:,0]) 
+    x = np.zeros(np.shape(x_oo))
+    flag = np.zeros(np.shape(x_oo)[0])
+    for i in range(0, len(ind_0)): 
+        x[i,:]= x_oo[ind_0[i],:]
+        # When interface is not present, z location is set to zlo
+        flag[i] = 1.0 if x[i,2] > 1e-8 else 0.
+  
+    # Shorten arrays (exclude points where interface not detected)
+    nvalid = int(flag.sum())
+    # Columns are position (3), id, cpu, ints of struct (uid, sid, nid), field components
+    xshort = x[0:nvalid,0]
+    zshort = x[0:nvalid,2]
+    vshort = np.sqrt(x[0:nvalid,8]**2+x[0:nvalid,9]**2+x[0:nvalid,10]**2)
+
+    color = base_color - (base_color) * (n/nt)
+    cstr = str(color)
+
+    plt.plot(xshort,zshort,color=cstr)
+    plt.scatter(xshort,zshort,c=vshort,cmap="jet",vmin=0.,vmax=2.)
+
+plt.ylabel('$z$',fontsize=16)
+plt.xlabel('$x$',fontsize=16)
+plt.colorbar()
+plt.savefig('plot_sampling_ascii.pdf',format='pdf',dpi=300,bbox_inches="tight")
\ No newline at end of file
diff --git a/tools/sampling_dam_break_godunov_native.py b/tools/sampling_dam_break_godunov_native.py
new file mode 100644
index 0000000000..7f980648cd
--- /dev/null
+++ b/tools/sampling_dam_break_godunov_native.py
@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+AMR_WIND_PATH = '.'
+sys.path.append(AMR_WIND_PATH+'/tools/')
+import amrex_particle
+from amrex_particle import AmrexParticleFile
+
+loc_dir = "."
+pp_dir = loc_dir + "/post_processing"
+
+pfile = AmrexParticleFile(loc_dir)
+
+nt = 5 + 1
+out_int = 30
+
+plt.figure()
+base_color = 0.8
+
+for n in range(nt):
+    pt = pfile.load(n * out_int, root_dir = pp_dir)
+    pt.parse_header()
+    pt.load_binary_data()
+    data = pt.df
+    x_oo = data.xco
+    z_oo = data.zco
+    u_oo = data.velocityx
+    v_oo = data.velocityy
+    w_oo = data.velocityz
+    # Reorder arrays
+    ind_0 = np.argsort(x_oo) 
+    x = np.zeros(np.shape(x_oo))
+    z = np.zeros(np.shape(x_oo))
+    Vmag = np.zeros(np.shape(x_oo))
+    flag = np.zeros(np.shape(x_oo))
+    for i in range(0, len(ind_0)): 
+        x[i]= x_oo[ind_0[i]]
+        z[i]= z_oo[ind_0[i]]
+        Vmag[i] = np.sqrt(u_oo[ind_0[i]]**2+v_oo[ind_0[i]]**2+w_oo[ind_0[i]]**2)
+        # When interface is not present, z location is set to zlo
+        flag[i] = 1.0 if z[i] > 1e-8 else 0.
+    # Shorten arrays (exclude points where interface not detected)
+    nvalid = int(flag.sum())
+    xshort = x[0:nvalid]
+    zshort = z[0:nvalid]
+    vshort = Vmag[0:nvalid]
+
+    color = base_color - (base_color) * (n/nt)
+    cstr = str(color)
+
+    plt.plot(xshort,zshort,color=cstr)
+    plt.scatter(xshort,zshort,c=vshort,cmap="jet",vmin=0.,vmax=2.)
+
+plt.ylabel(r'$z$',fontsize=16)
+plt.xlabel(r'$x$',fontsize=16)
+plt.colorbar()
+plt.savefig('plot_sampling_native.pdf',format='pdf',dpi=300,bbox_inches="tight")
\ No newline at end of file
diff --git a/tools/sampling_dam_break_godunov_netcdf.py b/tools/sampling_dam_break_godunov_netcdf.py
new file mode 100644
index 0000000000..ef3b6ff261
--- /dev/null
+++ b/tools/sampling_dam_break_godunov_netcdf.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import netCDF4 as ncdf
+import matplotlib.pyplot as plt
+
+loc_dir = "."
+pp_dir = loc_dir + "/post_processing"
+
+nt = 5 + 1
+out_int = 30
+
+class SamplingFile(object):
+    """Interface to Sampling NetCDF file"""
+
+    def __init__(self, sampling_file = "sampling.nc", group_name = "s1"):
+        """
+        Args:
+            stats_file (path): Absolute path to the NetCDF file
+        """
+        self.sampling_file = sampling_file
+        self.sampling = ncdf.Dataset(self.sampling_file)
+        self.fs = self.sampling["/"+group_name]
+        self._points = self.fs.variables["points"][:,:,:]
+        self._velocityx = self.fs.variables["velocityx"][:,:]
+        self._velocityy = self.fs.variables["velocityy"][:,:]
+        self._velocityz = self.fs.variables["velocityz"][:,:]
+
+    @property
+    def locations(self):
+        return self._points
+
+    @property
+    def vmag(self):
+        return np.sqrt(self._velocityx**2+self._velocityy**2+self._velocityz**2)
+
+plt.figure()
+base_color = 0.8
+
+sampl = SamplingFile(sampling_file=pp_dir+"/sampling00000.nc", \
+                    group_name = "fs")
+
+for n in range(nt):
+    x_oo=sampl.locations[n,:,:]
+    Vmag_oo = sampl.vmag[n,:]
+
+    # Reorder arrays
+    ind_0 = np.argsort(x_oo[:,0]) 
+    x = np.zeros(np.shape(x_oo)[0])
+    z = np.zeros(np.shape(x_oo)[0])
+    Vmag = np.zeros(np.shape(x_oo)[0])
+    flag = np.zeros(np.shape(x_oo)[0])
+    for i in range(0, len(ind_0)): 
+        x[i]= x_oo[ind_0[i],0]
+        z[i]= x_oo[ind_0[i],2]
+        Vmag[i] = Vmag_oo[ind_0[i]]
+        # When interface is not present, z location is set to 0
+        flag[i] = 1.0 if z[i] > 1e-8 else 0.
+    # Shorten arrays (exclude points where interface not detected)
+    nvalid = int(flag.sum())
+    xshort = x[0:nvalid]
+    zshort = z[0:nvalid]
+    vshort = Vmag[0:nvalid]
+
+    color = base_color - (base_color) * (n/nt)
+    cstr = str(color)
+
+    plt.plot(xshort,zshort,color=cstr)
+    plt.scatter(xshort,zshort,c=vshort,cmap="jet",vmin=0.,vmax=2.)
+
+plt.ylabel('$z$',fontsize=16)
+plt.xlabel('$x$',fontsize=16)
+plt.colorbar()
+plt.savefig('plot_sampling_netcdf.pdf',format='pdf',dpi=300,bbox_inches="tight")
\ No newline at end of file
diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp
index 847396b3a1..d8b7cbc680 100644
--- a/unit_tests/utilities/test_free_surface.cpp
+++ b/unit_tests/utilities/test_free_surface.cpp
@@ -1,5 +1,5 @@
 #include "aw_test_utils/MeshTest.H"
-#include "amr-wind/utilities/sampling/FreeSurface.H"
+#include "amr-wind/utilities/sampling/FreeSurfaceSampler.H"
 #include "amr-wind/utilities/tagging/FieldRefinement.H"
 
 namespace amr_wind_tests {
@@ -191,22 +191,20 @@ class FSRefineMesh : public AmrTestMesh
     std::unique_ptr<amr_wind::RefineCriteriaManager> m_mesh_refiner;
 };
 
-class FreeSurfaceImpl : public amr_wind::free_surface::FreeSurface
+class FreeSurfaceImpl : public amr_wind::sampling::FreeSurfaceSampler
 {
 public:
-    FreeSurfaceImpl(amr_wind::CFDSim& sim, const std::string& label)
-        : amr_wind::free_surface::FreeSurface(sim, label)
+    FreeSurfaceImpl(amr_wind::CFDSim& sim)
+        : amr_wind::sampling::FreeSurfaceSampler(sim)
     {}
     int check_output(const std::string& op, amrex::Real check_val);
     int check_output(int cidx, const std::string& op, amrex::Real check_val);
     int check_output_vec(
         const std::string& op, amrex::Vector<amrex::Real> check_val);
     int check_pos(int cidx, const std::string& op, amrex::Real check_val);
+    int check_sloc(const std::string& op);
 
 protected:
-    // No file output during test
-    void prepare_netcdf_file() override {}
-    void process_output() override {}
     const amrex::Real m_tol = 1e-8;
 };
 
@@ -274,7 +272,7 @@ int FreeSurfaceImpl::check_pos(
 {
     // Get number of points and position array
     auto npts_tot = num_gridpoints();
-    auto locs = locations();
+    auto locs = grid_locations();
     // Loop through grid points and check output
     int icheck = 0;
     for (int n = 0; n < npts_tot; ++n) {
@@ -296,6 +294,42 @@ int FreeSurfaceImpl::check_pos(
     return icheck;
 }
 
+int FreeSurfaceImpl::check_sloc(const std::string& op)
+{
+    // Get number of points and sampling locations array
+    auto npts_tot = num_points();
+    amrex::Vector<amrex::Array<amrex::Real, AMREX_SPACEDIM>> locs;
+    output_locations(locs);
+    // Get locations from other functions
+    auto gridlocs = grid_locations();
+    auto out = heights();
+    // Loop through grid points and check output
+    int icheck = 0;
+    for (int n = 0; n < npts_tot; ++n) {
+        if (op == "=") {
+            EXPECT_EQ(locs[n][0], gridlocs[n][0]);
+            EXPECT_EQ(locs[n][1], gridlocs[n][1]);
+            EXPECT_EQ(locs[n][2], out[n]);
+            ++icheck;
+        } else {
+            if (op == "<") {
+                EXPECT_LT(locs[n][0], gridlocs[n][0]);
+                EXPECT_LT(locs[n][1], gridlocs[n][1]);
+                EXPECT_LT(locs[n][2], out[n]);
+                ++icheck;
+            } else {
+                if (op == "~") {
+                    EXPECT_NEAR(locs[n][0], gridlocs[n][0], m_tol);
+                    EXPECT_NEAR(locs[n][1], gridlocs[n][1], m_tol);
+                    EXPECT_NEAR(locs[n][2], out[n], m_tol);
+                    ++icheck;
+                }
+            }
+        }
+    }
+    return icheck;
+}
+
 } // namespace
 
 class FreeSurfaceTest : public MeshTest
@@ -324,9 +358,9 @@ class FreeSurfaceTest : public MeshTest
         amrex::ParmParse pp(fsname);
         pp.add("output_frequency", 1);
         pp.add("num_instances", ninst);
-        pp.addarr("num_points", amrex::Vector<int>{1, 1});
-        pp.addarr("start", m_pt_coord);
-        pp.addarr("end", m_pt_coord);
+        pp.addarr("plane_num_points", amrex::Vector<int>{1, 1});
+        pp.addarr("plane_start", m_pt_coord);
+        pp.addarr("plane_end", m_pt_coord);
     }
     void setup_grid_0d(int ninst) { setup_grid_0d(ninst, "freesurface"); }
     void setup_grid_2d(int ninst)
@@ -334,19 +368,20 @@ class FreeSurfaceTest : public MeshTest
         amrex::ParmParse pp("freesurface");
         pp.add("output_frequency", 1);
         pp.add("num_instances", ninst);
-        pp.addarr("num_points", amrex::Vector<int>{npts, npts});
-        pp.addarr("start", m_pl_start);
-        pp.addarr("end", m_pl_end);
+        pp.addarr("plane_num_points", amrex::Vector<int>{npts, npts});
+        pp.addarr("plane_start", m_pl_start);
+        pp.addarr("plane_end", m_pl_end);
     }
-    void setup_grid_2d_narrow()
+    void setup_grid_2d_narrow(const std::string& fsname)
     {
-        amrex::ParmParse pp("freesurface");
+        amrex::ParmParse pp(fsname);
         pp.add("output_frequency", 1);
         pp.add("num_instances", 1);
-        pp.addarr("num_points", amrex::Vector<int>{npts, npts});
-        pp.addarr("start", m_plnarrow_s);
-        pp.addarr("end", m_plnarrow_e);
+        pp.addarr("plane_num_points", amrex::Vector<int>{npts, npts});
+        pp.addarr("plane_start", m_plnarrow_s);
+        pp.addarr("plane_end", m_plnarrow_e);
     }
+    void setup_grid_2d_narrow() { setup_grid_2d_narrow("freesurface"); }
     void setup_fieldrefinement()
     {
         amrex::ParmParse pp("tagging");
@@ -382,9 +417,9 @@ TEST_F(FreeSurfaceTest, point)
 
     init_vof(vof, m_water_level0);
     auto& m_sim = sim();
-    FreeSurfaceImpl tool(m_sim, "freesurface");
-    tool.initialize();
-    tool.post_advance_work();
+    FreeSurfaceImpl tool(m_sim);
+    tool.initialize("freesurface");
+    tool.update_sampling_locations();
 
     // Check number of points
     auto ngp = tool.num_gridpoints();
@@ -397,6 +432,9 @@ TEST_F(FreeSurfaceTest, point)
     // Check output value
     int nout = tool.check_output("~", m_water_level0);
     ASSERT_EQ(nout, 1);
+    // Check sampling locations
+    int nsloc = tool.check_sloc("~");
+    ASSERT_EQ(nsloc, 1);
 }
 
 TEST_F(FreeSurfaceTest, plane)
@@ -408,9 +446,9 @@ TEST_F(FreeSurfaceTest, plane)
 
     init_vof(vof, m_water_level1);
     auto& m_sim = sim();
-    FreeSurfaceImpl tool(m_sim, "freesurface");
-    tool.initialize();
-    tool.post_advance_work();
+    FreeSurfaceImpl tool(m_sim);
+    tool.initialize("freesurface");
+    tool.update_sampling_locations();
 
     // Check number of points
     auto ngp = tool.num_gridpoints();
@@ -434,9 +472,9 @@ TEST_F(FreeSurfaceTest, multivalued)
 
     init_vof_multival(vof, wl0, wl1, wl2);
     auto& m_sim = sim();
-    FreeSurfaceImpl tool(m_sim, "freesurface");
-    tool.initialize();
-    tool.post_advance_work();
+    FreeSurfaceImpl tool(m_sim);
+    tool.initialize("freesurface");
+    tool.update_sampling_locations();
 
     // Check number of outputs
     auto heights = tool.heights();
@@ -462,9 +500,9 @@ TEST_F(FreeSurfaceTest, sloped)
     amrex::Real domain_l = m_probhi[0];
     init_vof_slope(vof, m_water_level2, slope, domain_l);
     auto& m_sim = sim();
-    FreeSurfaceImpl tool(m_sim, "freesurface");
-    tool.initialize();
-    tool.post_advance_work();
+    FreeSurfaceImpl tool(m_sim);
+    tool.initialize("freesurface");
+    tool.update_sampling_locations();
 
     // Calculate expected output values
     amrex::Vector<amrex::Real> out_vec(static_cast<long>(npts * npts), 0.0);
@@ -493,19 +531,22 @@ TEST_F(FreeSurfaceTest, multisampler)
     setup_grid_0d(1, "freesurface0");
 
     // Set up parameters for another sampler
-    setup_grid_2d_narrow();
+    setup_grid_2d_narrow("freesurface1");
 
     // Initialize VOF distribution and access sim
     init_vof(vof, m_water_level1);
     auto& m_sim = sim();
 
     // Initialize first sampler
-    FreeSurfaceImpl tool1(m_sim, "freesurface0");
-    tool1.initialize();
+    FreeSurfaceImpl tool1(m_sim);
+    // Populate label, would be done by Sampling
+    tool1.label() = "0";
+    tool1.initialize("freesurface0");
 
     // Initialize second sampler
-    FreeSurfaceImpl tool2(m_sim, "freesurface");
-    tool2.initialize();
+    FreeSurfaceImpl tool2(m_sim);
+    tool1.label() = "1";
+    tool2.initialize("freesurface1");
 }
 
 TEST_F(FreeSurfaceTest, regrid)
@@ -540,9 +581,9 @@ TEST_F(FreeSurfaceTest, regrid)
     auto& rsim = rmesh.sim();
 
     // Initialize sampler and check result on initial mesh
-    FreeSurfaceImpl tool(rsim, "freesurface");
-    tool.initialize();
-    tool.post_advance_work();
+    FreeSurfaceImpl tool(rsim);
+    tool.initialize("freesurface");
+    tool.update_sampling_locations();
     tool.check_output("~", m_water_level1);
 
     // Change scalar for determining refinement - no fine level
@@ -553,7 +594,7 @@ TEST_F(FreeSurfaceTest, regrid)
     tool.post_regrid_actions();
 
     // Check that result is unchanged on new mesh
-    tool.post_advance_work();
+    tool.update_sampling_locations();
     tool.check_output("~", m_water_level1);
 }