diff --git a/amr-wind/core/SimTime.H b/amr-wind/core/SimTime.H
index c939560ce2..070b097648 100644
--- a/amr-wind/core/SimTime.H
+++ b/amr-wind/core/SimTime.H
@@ -40,6 +40,8 @@ public:
      */
     bool new_timestep();
 
+    void increment_timestep() {++m_time_index;}
+
     /** Return true if simulation should continue
      */
     bool continue_simulation() const;
diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp
index 25348f3066..6127034369 100644
--- a/amr-wind/equation_systems/icns/icns_advection.cpp
+++ b/amr-wind/equation_systems/icns/icns_advection.cpp
@@ -50,6 +50,11 @@ MacProjOp::MacProjOp(
 {
     amrex::ParmParse pp("incflo");
     pp.query("density", m_rho_0);
+    bool disable_omac{true};
+    pp.query("disable_overset_mac", disable_omac);
+    if (m_has_overset && disable_omac) {
+        m_has_overset = false;
+    }
 }
 
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H
index 34068e43a1..1dbf433c61 100644
--- a/amr-wind/equation_systems/vof/vof_advection.H
+++ b/amr-wind/equation_systems/vof/vof_advection.H
@@ -23,6 +23,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
     {
         amrex::ParmParse pp_multiphase("VOF");
         pp_multiphase.query("remove_debris", m_rm_debris);
+        pp_multiphase.query("sharpen_acq_field", m_sharpen_acq);
+        pp_multiphase.query("sharpen_dest_field", m_sharpen_dest);
+        pp_multiphase.query("replace_masked", m_replace_mask);
 
         // Setup density factor arrays for multiplying velocity flux
         fields_in.repo.declare_face_normal_field(
@@ -34,7 +37,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
         const FieldState /*unused*/,
         const amrex::Real /*unused*/,
         const amrex::Real /*unused*/)
-    {}
+    {
+        // Would be best to sharpen old vof and update old density here
+    }
 
     void operator()(const FieldState /*unused*/, const amrex::Real dt)
     {
@@ -43,13 +48,38 @@ struct AdvectionOp<VOF, fvm::Godunov>
 
         auto& repo = fields.repo;
         const auto& geom = repo.mesh().Geom();
+        const int nlevels = repo.num_active_levels();
 
         auto& aa_x = repo.get_field("advalpha_x");
         auto& aa_y = repo.get_field("advalpha_y");
         auto& aa_z = repo.get_field("advalpha_z");
 
-        // cppcheck-suppress constVariable
-        auto& dof_field = fields.field;
+        // Old and new states
+        auto& dof_old = fields.field.state(amr_wind::FieldState::Old);
+        auto& dof_new = fields.field;
+        // Working state of vof is nph, to keep others untouched during step
+        auto& dof_field = fields.field.state(amr_wind::FieldState::NPH);
+
+        // Sharpen starting vof field if hybrid solver is being used
+        if (repo.int_field_exists("iblank_cell") && m_sharpen_acq) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::sharpen_acquired_vof(
+                nlevels, f_iblank, dof_old);
+        }
+        // Sharpen future vof field if hybrid solver is being used
+        if (repo.int_field_exists("iblank_cell") && m_sharpen_dest) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::sharpen_acquired_vof(
+                nlevels, f_iblank, dof_new);
+        }
+
+        // Initialize as old state values
+        for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
+            amrex::MultiFab::Copy(
+                dof_field(lev), dof_old(lev), 0, 0, dof_field.num_comp(),
+                dof_field.num_grow());
+        }
+
         //
         // Advect volume using Implicit Eulerian Sweeping method with PLIC
         // reconstruction
@@ -71,8 +101,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
             isweep = 1;
         }
 
-        const int nlevels = repo.num_active_levels();
-
         amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
             nlevels);
         amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> advas(
@@ -86,20 +114,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
             advas[lev][2] = &aa_z(lev);
         }
 
-        // Sharpen acquired vof field if hybrid solver is being used
-        if (repo.int_field_exists("iblank_cell")) {
-            auto& f_iblank = repo.get_int_field("iblank_cell");
-            amr_wind::multiphase::sharpen_acquired_vof(
-                nlevels, f_iblank, dof_field);
-            // Update old vof so that old density can be updated
-            // cppcheck-suppress constVariableReference
-            auto& old_dof_field = dof_field.state(amr_wind::FieldState::Old);
-            for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
-                amrex::MultiFab::Copy(
-                    old_dof_field(lev), dof_field(lev), 0, 0,
-                    dof_field.num_comp(), dof_field.num_grow());
-            }
-        }
         // Split advection step 1, with cmask calculation
         multiphase::split_advection_step(
             isweep, 0, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
@@ -112,6 +126,20 @@ struct AdvectionOp<VOF, fvm::Godunov>
         multiphase::split_advection_step(
             isweep, 2, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
             v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);
+
+        // Replace masked cells using overset
+        if (repo.int_field_exists("iblank_cell") && m_replace_mask) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::replace_masked_vof(
+                nlevels, f_iblank, dof_field, dof_new);
+        }
+
+        // Copy working version of vof to new state
+        for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
+            amrex::MultiFab::Copy(
+                dof_new(lev), dof_field(lev), 0, 0, dof_field.num_comp(),
+                dof_field.num_grow());
+        }
     }
 
     PDEFields& fields;
@@ -120,6 +148,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
     Field& w_mac;
     int isweep = 0;
     bool m_rm_debris{true};
+    bool m_sharpen_acq{false};
+    bool m_sharpen_dest{false};
+    bool m_replace_mask{true};
     // Lagrangian transport is deprecated, only Eulerian is supported
 };
 
diff --git a/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H b/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
index fefa83f454..f65ecaca40 100644
--- a/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
+++ b/amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
@@ -45,6 +45,36 @@ static void sharpen_acquired_vof(
     }
 }
 
+static void replace_masked_vof(
+    const int nlevels,
+    amr_wind::IntField& f_iblank,
+    amr_wind::Field& f_vof,
+    amr_wind::Field& f_vof_new)
+{
+    // Sharpen data from nalu-wind (in iblank regions)
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& iblank = f_iblank(lev);
+        auto& vof = f_vof(lev);
+        const auto& vof_new = f_vof_new(lev);
+
+        for (amrex::MFIter mfi(iblank); mfi.isValid(); ++mfi) {
+            const auto& gbx = mfi.growntilebox();
+            const amrex::Array4<const int>& native_flag =
+                iblank.const_array(mfi);
+            const amrex::Array4<amrex::Real>& volfrac = vof.array(mfi);
+            const amrex::Array4<const amrex::Real>& vfmasked =
+                vof_new.const_array(mfi);
+            amrex::ParallelFor(
+                gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    // In iblanked regions, sharpen VOF and limit it
+                    volfrac(i, j, k) = (native_flag(i, j, k) > 0)
+                                           ? volfrac(i, j, k)
+                                           : vfmasked(i, j, k);
+                });
+        }
+    }
+}
+
 } // namespace amr_wind::multiphase
 
 #endif // VOF_HYBRIDSOLVER_OPS.H
diff --git a/amr-wind/equation_systems/vof/volume_fractions.H b/amr-wind/equation_systems/vof/volume_fractions.H
index 89284bd535..0503fb411b 100644
--- a/amr-wind/equation_systems/vof/volume_fractions.H
+++ b/amr-wind/equation_systems/vof/volume_fractions.H
@@ -64,6 +64,65 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
     mz = mm1 - mm2;
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal_neumann(
+    int i,
+    int j,
+    int k,
+    int ibdy,
+    int jbdy,
+    int kbdy,
+    amrex::Array4<amrex::Real const> const& volfrac,
+    amrex::Real& mx,
+    amrex::Real& my,
+    amrex::Real& mz) noexcept
+{
+    amrex::Real mm1, mm2;
+
+    // Do neumann condition via indices
+    int im1 = ibdy == -1 ? i : i - 1;
+    int jm1 = jbdy == -1 ? j : j - 1;
+    int km1 = kbdy == -1 ? k : k - 1;
+    int ip1 = ibdy == +1 ? i : i + 1;
+    int jp1 = jbdy == +1 ? j : j + 1;
+    int kp1 = kbdy == +1 ? k : k + 1;
+
+    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
+          volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
+          2.0 * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) +
+                 volfrac(im1, j, km1) + volfrac(im1, j, kp1)) +
+          4.0 * volfrac(im1, j, k);
+    mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
+          volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
+          2.0 * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) +
+                 volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) +
+          4.0 * volfrac(ip1, j, k);
+    mx = mm1 - mm2;
+
+    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
+          volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
+          2.0 * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) +
+                 volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) +
+          4.0 * volfrac(i, jm1, k);
+    mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
+          volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
+          2.0 * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) +
+                 volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) +
+          4.0 * volfrac(i, jp1, k);
+    my = mm1 - mm2;
+
+    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) +
+          volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) +
+          2.0 * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) +
+                 volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) +
+          4.0 * volfrac(i, j, km1);
+    mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) +
+          volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) +
+          2.0 * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) +
+                 volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) +
+          4.0 * volfrac(i, j, kp1);
+    mz = mm1 - mm2;
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
     int i,
     int j,
diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H
index 3eca0d7dff..ba97cc6e20 100644
--- a/amr-wind/incflo.H
+++ b/amr-wind/incflo.H
@@ -173,10 +173,25 @@ private:
     bool m_do_initial_proj = true;
     int m_initial_iterations = 3;
 
+    // Sharpening multiphase for overset
+    int m_sharpen_iterations = 10;
+    amrex::Real m_sharpen_tolerance = 1e-12;
+    int m_sharpen_calctolniter = 1;
+    amrex::Real m_sharpen_rlscale = 1.5;
+    amrex::Real m_sharpen_margin = 0.1;
+    amrex::Real m_sharpen_proctg_tol = 0.0;
+    bool m_sharpen_hs_pressure = false;
+    bool m_sharpen_hsp_guess = true;
+    bool m_sharpen_hsp_replace = true;
+    bool m_sharpen_gradp = false;
+
     bool m_constant_density = true;
 
     bool m_use_godunov = true;
 
+    // Other choices for overset
+    bool m_disable_onodal = true;
+
     // Prescribe advection velocity
     bool m_prescribe_vel = false;
 
diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp
index d93ec12467..eacc5d40b6 100644
--- a/amr-wind/incflo.cpp
+++ b/amr-wind/incflo.cpp
@@ -388,6 +388,21 @@ void incflo::init_physics_and_pde()
         // Check for if velocity is prescribed
         amrex::ParmParse pp("incflo");
         pp.query("prescribe_velocity", m_prescribe_vel);
+
+        // Check for sharpening iterations (overset feature)
+        pp.query("sharpen_iterations", m_sharpen_iterations);
+        pp.query("sharpen_tolerance", m_sharpen_tolerance);
+        pp.query("sharpen_calctol_niter", m_sharpen_calctolniter);
+        pp.query("sharpen_rlscale", m_sharpen_rlscale);
+        pp.query("sharpen_margin", m_sharpen_margin);
+        pp.query("sharpen_targetvof_tol", m_sharpen_proctg_tol);
+        pp.query("sharpen_hs_pressure", m_sharpen_hs_pressure);
+        pp.query("sharpen_guess_hsp", m_sharpen_hsp_guess);
+        pp.query("sharpen_replace_hsp", m_sharpen_hsp_replace);
+        pp.query("sharpen_pressure_grad", m_sharpen_gradp);
+
+        // Determine if overset values should be forced into projection
+        pp.query("disable_overset_nodal", m_disable_onodal);
     }
     m_sim.create_turbulence_model();
 
diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp
index b4ef5d95ac..6242837390 100644
--- a/amr-wind/incflo_advance.cpp
+++ b/amr-wind/incflo_advance.cpp
@@ -8,6 +8,7 @@
 #include "amr-wind/utilities/console_io.H"
 #include "amr-wind/utilities/PostProcessing.H"
 #include "AMReX_MultiFabUtil.H"
+#include "amr-wind/overset/sharpen_nalu_data.H"
 
 using namespace amrex;
 
@@ -190,11 +191,30 @@ void incflo::ApplyPredictor(bool incremental_projection)
     const auto& density_old = density_new.state(amr_wind::FieldState::Old);
     auto& density_nph = density_new.state(amr_wind::FieldState::NPH);
 
-    // Recalculate incoming pressure gradient field if overset
+    auto gp_copy = density().repo().create_scratch_field(3);
+    auto& gp = density().repo().get_field("gp");
+
+    // Process data for overset multiphase
     if (sim().has_overset()) {
-        UpdateGradP(
-            (density_old).vec_const_ptrs(), m_time.current_time(),
-            m_time.deltaT());
+        // Sharpen nalu fields
+        amr_wind::overset::SharpenNaluDataDiscrete(
+            sim(), m_sharpen_iterations, m_sharpen_tolerance,
+            m_sharpen_calctolniter, m_sharpen_rlscale, m_sharpen_margin,
+            m_sharpen_proctg_tol, m_sharpen_hs_pressure, m_sharpen_gradp);
+        // Recalculate pressure gradient with incoming sharpened field
+        if (!m_sharpen_gradp || sim().time().current_time() == 0.0) {
+            UpdateGradP(
+                (density_old).vec_const_ptrs(), m_time.current_time(),
+                m_time.deltaT());
+        }
+        if (m_sharpen_hsp_guess ||
+            (m_sharpen_gradp && sim().time().current_time() == 0.0)) {
+            amr_wind::overset::ReplaceMaskedGradP(sim());
+        }
+        if (m_sharpen_gradp) {
+            amr_wind::overset::CopyGradP(
+                *gp_copy, gp, sim().repo().num_active_levels());
+        }
     }
 
     // *************************************************************************************
@@ -269,6 +289,10 @@ void incflo::ApplyPredictor(bool incremental_projection)
         }
     }
 
+    if (m_verbose > 2) {
+        PrintMaxVelLocations("before pre advection");
+    }
+
     // Extrapolate and apply MAC projection for advection velocities
     icns().pre_advection_actions(amr_wind::FieldState::Old);
 
@@ -396,6 +420,13 @@ void incflo::ApplyPredictor(bool incremental_projection)
         (density_new).vec_const_ptrs(), new_time, m_time.deltaT(),
         incremental_projection);
 
+    if (m_sharpen_hsp_replace) {
+        amr_wind::overset::ReapplyModifiedGradP(sim());
+    }
+    if (m_sharpen_gradp) {
+        amr_wind::overset::ReapplyOversetGradP(*gp_copy, sim());
+    }
+
     if (m_verbose > 2) {
         PrintMaxVelLocations("after nodal projection");
     }
diff --git a/amr-wind/overset/TiogaInterface.H b/amr-wind/overset/TiogaInterface.H
index 7a2c513f0d..ba77f43d66 100644
--- a/amr-wind/overset/TiogaInterface.H
+++ b/amr-wind/overset/TiogaInterface.H
@@ -104,6 +104,8 @@ private:
 
     CFDSim& m_sim;
 
+    bool m_disable_pressure_from_nalu{true};
+
     //! IBLANK on cell centered fields
     IntField& m_iblank_cell;
 
diff --git a/amr-wind/overset/TiogaInterface.cpp b/amr-wind/overset/TiogaInterface.cpp
index 59e42d414d..cf96ae389d 100644
--- a/amr-wind/overset/TiogaInterface.cpp
+++ b/amr-wind/overset/TiogaInterface.cpp
@@ -4,6 +4,7 @@
 #include "amr-wind/equation_systems/PDEBase.H"
 #include "amr-wind/core/field_ops.H"
 #include "amr-wind/utilities/IOManager.H"
+#include "AMReX_ParmParse.H"
 
 #include <memory>
 #include <numeric>
@@ -77,6 +78,9 @@ TiogaInterface::TiogaInterface(CFDSim& sim)
           FieldLoc::NODE))
 {
     m_sim.io_manager().register_output_int_var(m_iblank_cell.name());
+
+    amrex::ParmParse pp("Overset");
+    pp.query("ignore_nalu_pressure_field", m_disable_pressure_from_nalu);
 }
 
 // clang-format on
@@ -127,6 +131,9 @@ void TiogaInterface::post_overset_conn_work()
     for (int lev = 0; lev < nlevels; ++lev) {
         htod_memcpy(m_iblank_cell(lev), (*m_iblank_cell_host)(lev), 0, 0, 1);
         htod_memcpy(m_iblank_node(lev), (*m_iblank_node_host)(lev), 0, 0, 1);
+
+        m_iblank_cell(lev).FillBoundary(m_sim.mesh().Geom()[lev].periodicity());
+        m_iblank_node(lev).FillBoundary(m_sim.mesh().Geom()[lev].periodicity());
     }
 
     iblank_to_mask(m_iblank_cell, m_mask_cell);
@@ -282,7 +289,8 @@ void TiogaInterface::update_solution()
     }
 
     // Update nodal variables on device
-    {
+    if (!m_disable_pressure_from_nalu) {
+        // Pressure is the only nodal variable - skip copy if requested
         int icomp = 0;
         for (const auto& cvar : m_node_vars) {
             auto& fld = repo.get_field(cvar);
diff --git a/amr-wind/overset/sharpen_nalu_data.H b/amr-wind/overset/sharpen_nalu_data.H
new file mode 100644
index 0000000000..89942a24b9
--- /dev/null
+++ b/amr-wind/overset/sharpen_nalu_data.H
@@ -0,0 +1,1020 @@
+#ifndef SHARPEN_NALU_DATA_H_
+#define SHARPEN_NALU_DATA_H_
+
+#include "amr-wind/CFDSim.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
+#include "amr-wind/equation_systems/vof/volume_fractions.H"
+#include "amr-wind/utilities/IOManager.H"
+
+namespace amr_wind {
+namespace overset {
+
+// Approximate signed distance function
+amrex::Real asdf(const amrex::Real a_vof, const amrex::Real i_th)
+{
+    // function of local vof value and interface thickness
+    return (i_th * log((a_vof + 1e-12) / (1. - a_vof + 1e-12)));
+}
+
+void process_vof(amrex::MultiFab& mf_vof, const amrex::Real vof_tol)
+{
+    for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) {
+        const auto& gbx = mfi.growntilebox();
+        const amrex::Array4<amrex::Real>& vof = mf_vof.array(mfi);
+        // Populate approximate signed distance function
+        amrex::ParallelFor(
+            gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // This conditional is just for the sake of debugging
+                if (vof(i, j, k) < 1.0 && vof(i, j, k) > 0.0) {
+                    vof(i, j, k) =
+                        vof(i, j, k) < vof_tol
+                            ? 0.0
+                            : (vof(i, j, k) > 1. - vof_tol ? 1. : vof(i, j, k));
+                }
+            });
+    }
+}
+
+void populate_psi(
+    amrex::MultiFab& mf_psi, amrex::MultiFab& mf_vof, const amrex::Real i_th)
+{
+    for (amrex::MFIter mfi(mf_psi); mfi.isValid(); ++mfi) {
+        const auto& gbx = mfi.growntilebox();
+        const amrex::Array4<amrex::Real>& psi = mf_psi.array(mfi);
+        const amrex::Array4<const amrex::Real>& vof = mf_vof.const_array(mfi);
+        // Populate approximate signed distance function
+        amrex::ParallelFor(
+            gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                psi(i, j, k) = asdf(vof(i, j, k), i_th);
+            });
+    }
+}
+
+void populate_normal_vector(
+    amrex::MultiFab& mf_normvec,
+    amrex::MultiFab& mf_vof,
+    amrex::iMultiFab& mf_iblank)
+{
+    for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) {
+        const auto& gbxm1 = grow(mfi.growntilebox(), -1);
+        const amrex::Array4<amrex::Real>& normvec = mf_normvec.array(mfi);
+        const amrex::Array4<const amrex::Real>& vof = mf_vof.const_array(mfi);
+        const amrex::Array4<const int>& iblank = mf_iblank.const_array(mfi);
+        // Calculate gradients in each direction with centered diff
+        // (Should this be replaced with youngs_fd_normal?)
+        amrex::ParallelFor(
+            gbxm1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // Neumann condition across nalu bdy
+                int ibdy = 0;
+                int jbdy = 0;
+                int kbdy = 0;
+                if (iblank(i, j, k) != iblank(i - 1, j, k)) {
+                    ibdy = -1;
+                }
+                if (iblank(i, j, k) != iblank(i, j - 1, k)) {
+                    jbdy = -1;
+                }
+                if (iblank(i, j, k) != iblank(i, j, k - 1)) {
+                    kbdy = -1;
+                }
+                // no cell should be isolated such that -1 and 1 are needed
+                if (iblank(i, j, k) != iblank(i + 1, j, k)) {
+                    ibdy = +1;
+                }
+                if (iblank(i, j, k) != iblank(i, j + 1, k)) {
+                    jbdy = +1;
+                }
+                if (iblank(i, j, k) != iblank(i, j, k + 1)) {
+                    kbdy = +1;
+                }
+                // Calculate normal
+                amrex::Real mx, my, mz, mmag;
+                multiphase::youngs_fd_normal_neumann(
+                    i, j, k, ibdy, jbdy, kbdy, vof, mx, my, mz);
+                // Normalize normal
+                mmag = std::sqrt(mx * mx + my * my + mz * mz + 1e-20);
+                // Save normal
+                normvec(i, j, k, 0) = mx / mmag;
+                normvec(i, j, k, 1) = my / mmag;
+                normvec(i, j, k, 2) = mz / mmag;
+            });
+    }
+}
+
+amrex::Real alpha_discrete_flux(
+    int i,
+    int j,
+    int k,
+    int dir,
+    const amrex::Real margin,
+    amrex::Array4<amrex::Real const> const& vof,
+    amrex::Array4<amrex::Real const> const& tg_vof,
+    amrex::Array4<amrex::Real const> const& normal)
+{
+    // Set up neighbor indices
+    int ii = i;
+    int jj = j;
+    int kk = k;
+    ii += (dir == 0) ? -1 : 0;
+    jj += (dir == 1) ? -1 : 0;
+    kk += (dir == 2) ? -1 : 0;
+
+    // Gradient of phi normal to interface
+    const amrex::Real gphi = (vof(i, j, k) - vof(ii, jj, kk));
+    // Normal vector in each cell (already normalized)
+    const amrex::Real norm_ = normal(i, j, k, dir);
+    const amrex::Real norm_nb = normal(ii, jj, kk, dir);
+
+    // Determine which delta_phi (and multiply by normal)
+    // The sign depends on side of flux face (like upwinding)
+    const amrex::Real dphi_ = (tg_vof(i, j, k) - vof(i, j, k)) * (-norm_);
+    const amrex::Real dphi_nb =
+        (tg_vof(ii, jj, kk) - vof(ii, jj, kk)) * norm_nb;
+    // Average value used across the interface
+    amrex::Real dphi_eval = 0.5 * (dphi_ + dphi_nb);
+    // Upwinding when on the gas side, downwinding on the liquid
+    // Across the interface defined as crossing 0.5 or within margin of 0.5
+    if ((std::abs(vof(i, j, k) - 0.5) > margin ||
+         std::abs(vof(ii, jj, kk) - 0.5) > margin)) {
+        if (gphi > 0.0) {
+            dphi_eval = (vof(ii, jj, kk) < 0.5 && vof(i, j, k) <= 0.5 + margin)
+                            ? dphi_nb
+                            : dphi_eval;
+            dphi_eval = (vof(ii, jj, kk) >= 0.5 - margin && vof(i, j, k) > 0.5)
+                            ? dphi_
+                            : dphi_eval;
+        }
+        if (gphi < 0.0) {
+            dphi_eval = (vof(i, j, k) < 0.5 && vof(ii, jj, kk) <= 0.5 + margin)
+                            ? dphi_
+                            : dphi_eval;
+            dphi_eval = (vof(i, j, k) >= 0.5 - margin && vof(ii, jj, kk) > 0.5)
+                            ? dphi_nb
+                            : dphi_eval;
+        }
+    }
+    return dphi_eval;
+}
+
+void velocity_discrete_face(
+    int i,
+    int j,
+    int k,
+    int dir,
+    amrex::Array4<amrex::Real const> const& vof,
+    amrex::Array4<amrex::Real const> const& velocity,
+    amrex::Real& uface,
+    amrex::Real& vface,
+    amrex::Real& wface)
+{
+    // Set up neighbor indices
+    int ii = i;
+    int jj = j;
+    int kk = k;
+    ii += (dir == 0) ? -1 : 0;
+    jj += (dir == 1) ? -1 : 0;
+    kk += (dir == 2) ? -1 : 0;
+
+    // Gradient of phi normal to interface
+    const amrex::Real gphi = (vof(i, j, k) - vof(ii, jj, kk));
+
+    // Get velocities on both sides
+    const amrex::Real u_ = velocity(i, j, k, 0);
+    const amrex::Real v_ = velocity(i, j, k, 1);
+    const amrex::Real w_ = velocity(i, j, k, 2);
+    const amrex::Real u_nb = velocity(ii, jj, kk, 0);
+    const amrex::Real v_nb = velocity(ii, jj, kk, 1);
+    const amrex::Real w_nb = velocity(ii, jj, kk, 2);
+    // Average value when gphi = 0
+    uface = 0.5 * (u_ + u_nb);
+    vface = 0.5 * (v_ + v_nb);
+    wface = 0.5 * (w_ + w_nb);
+    // Use simple upwinding
+    if (gphi > 0.0) {
+        uface = u_nb;
+        vface = v_nb;
+        wface = w_nb;
+    }
+    if (gphi < 0.0) {
+        uface = u_;
+        vface = v_;
+        wface = w_;
+    }
+}
+
+void gp_rho_discrete_face(
+    int i,
+    int j,
+    int k,
+    int dir,
+    amrex::Array4<amrex::Real const> const& vof,
+    amrex::Array4<amrex::Real const> const& gp,
+    amrex::Array4<amrex::Real const> const& rho,
+    amrex::Real& uface,
+    amrex::Real& vface,
+    amrex::Real& wface)
+{
+    // Set up neighbor indices
+    int ii = i;
+    int jj = j;
+    int kk = k;
+    ii += (dir == 0) ? -1 : 0;
+    jj += (dir == 1) ? -1 : 0;
+    kk += (dir == 2) ? -1 : 0;
+
+    // Gradient of phi normal to interface
+    const amrex::Real gphi = (vof(i, j, k) - vof(ii, jj, kk));
+
+    // Get velocities on both sides
+    const amrex::Real u_ = gp(i, j, k, 0) / rho(i, j, k);
+    const amrex::Real v_ = gp(i, j, k, 1) / rho(i, j, k);
+    const amrex::Real w_ = gp(i, j, k, 2) / rho(i, j, k);
+    const amrex::Real u_nb = gp(ii, jj, kk, 0) / rho(ii, jj, kk);
+    const amrex::Real v_nb = gp(ii, jj, kk, 1) / rho(ii, jj, kk);
+    const amrex::Real w_nb = gp(ii, jj, kk, 2) / rho(ii, jj, kk);
+    // Average value when gphi = 0
+    uface = 0.5 * (u_ + u_nb);
+    vface = 0.5 * (v_ + v_nb);
+    wface = 0.5 * (w_ + w_nb);
+    // Use simple upwinding
+    if (gphi > 0.0) {
+        uface = u_nb;
+        vface = v_nb;
+        wface = w_nb;
+    }
+    if (gphi < 0.0) {
+        uface = u_;
+        vface = v_;
+        wface = w_;
+    }
+}
+
+void populate_sharpen_discrete_fluxes(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_vof,
+    amrex::MultiFab& mf_target_vof,
+    amrex::MultiFab& mf_norm,
+    amrex::MultiFab& mf_velocity,
+    amrex::MultiFab& mf_gp,
+    amrex::MultiFab& mf_density,
+    const amrex::Real margin,
+    const amrex::Real rho1,
+    const amrex::Real rho2)
+{
+    for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        const auto& xbx = amrex::surroundingNodes(vbx, 0);
+        const auto& ybx = amrex::surroundingNodes(vbx, 1);
+        // Extra points for zbox are for pressure source
+        const auto& zbxp1 = grow(amrex::surroundingNodes(vbx, 2), 1);
+        const amrex::Array4<amrex::Real>& fx = mf_fx.array(mfi);
+        const amrex::Array4<amrex::Real>& fy = mf_fy.array(mfi);
+        const amrex::Array4<amrex::Real>& fz = mf_fz.array(mfi);
+        const amrex::Array4<const amrex::Real>& vof = mf_vof.const_array(mfi);
+        const amrex::Array4<const amrex::Real>& tg_vof =
+            mf_target_vof.const_array(mfi);
+        const amrex::Array4<const amrex::Real>& norm = mf_norm.const_array(mfi);
+        const amrex::Array4<const amrex::Real>& vel =
+            mf_velocity.const_array(mfi);
+        /*const amrex::Array4<const amrex::Real>& gp_rho =
+            mf_gp_rho.const_array(mfi);*/
+        const amrex::Array4<const amrex::Real>& gp = mf_gp.const_array(mfi);
+        const amrex::Array4<const amrex::Real>& rho =
+            mf_density.const_array(mfi);
+        // Populate vof and density fluxes for each direction
+        amrex::ParallelFor(
+            xbx, ybx, zbxp1,
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // vof flux
+                amrex::Real flux =
+                    alpha_discrete_flux(i, j, k, 0, margin, vof, tg_vof, norm);
+                fx(i, j, k, 0) = flux;
+                // density flux
+                flux *= (rho1 - rho2);
+                fx(i, j, k, 1) = flux;
+                // momentum fluxes (dens flux * face vel)
+                amrex::Real uf, vf, wf;
+                velocity_discrete_face(i, j, k, 0, vof, vel, uf, vf, wf);
+                fx(i, j, k, 2) = flux * uf;
+                fx(i, j, k, 3) = flux * vf;
+                fx(i, j, k, 4) = flux * wf;
+                gp_rho_discrete_face(i, j, k, 0, vof, gp, rho, uf, vf, wf);
+                fx(i, j, k, 5) = flux * uf;
+                fx(i, j, k, 6) = flux * vf;
+                fx(i, j, k, 7) = flux * wf;
+            },
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                amrex::Real flux =
+                    alpha_discrete_flux(i, j, k, 1, margin, vof, tg_vof, norm);
+                fy(i, j, k, 0) = flux;
+                flux *= (rho1 - rho2);
+                fy(i, j, k, 1) = flux;
+                amrex::Real uf, vf, wf;
+                velocity_discrete_face(i, j, k, 1, vof, vel, uf, vf, wf);
+                fy(i, j, k, 2) = flux * uf;
+                fy(i, j, k, 3) = flux * vf;
+                fy(i, j, k, 4) = flux * wf;
+                gp_rho_discrete_face(i, j, k, 1, vof, gp, rho, uf, vf, wf);
+                fy(i, j, k, 5) = flux * uf;
+                fy(i, j, k, 6) = flux * vf;
+                fy(i, j, k, 7) = flux * wf;
+            },
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                amrex::Real flux =
+                    alpha_discrete_flux(i, j, k, 2, margin, vof, tg_vof, norm);
+                fz(i, j, k, 0) = flux;
+                flux *= (rho1 - rho2);
+                fz(i, j, k, 1) = flux;
+                amrex::Real uf, vf, wf;
+                velocity_discrete_face(i, j, k, 2, vof, vel, uf, vf, wf);
+                fz(i, j, k, 2) = flux * uf;
+                fz(i, j, k, 3) = flux * vf;
+                fz(i, j, k, 4) = flux * wf;
+                gp_rho_discrete_face(i, j, k, 2, vof, gp, rho, uf, vf, wf);
+                fz(i, j, k, 5) = flux * uf;
+                fz(i, j, k, 6) = flux * vf;
+                fz(i, j, k, 7) = flux * wf;
+                // Turn "on" all z faces, later modified in process_fluxes
+                fz(i, j, k, 8) = 1.0;
+            });
+    }
+}
+
+void process_fluxes_calc_src(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_psource,
+    amrex::iMultiFab& mf_iblank,
+    const amrex::Real grav_z)
+{
+    int ncompx = mf_fx.nComp();
+    int ncompy = mf_fy.nComp();
+    int ncompz = mf_fz.nComp();
+    for (amrex::MFIter mfi(mf_iblank); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        const auto& xbx = amrex::surroundingNodes(vbx, 0);
+        const auto& ybx = amrex::surroundingNodes(vbx, 1);
+        const auto& zbxp1 = grow(amrex::surroundingNodes(vbx, 2), 1);
+        const amrex::Array4<amrex::Real>& fx = mf_fx.array(mfi);
+        const amrex::Array4<amrex::Real>& fy = mf_fy.array(mfi);
+        const amrex::Array4<amrex::Real>& fz = mf_fz.array(mfi);
+        const amrex::Array4<amrex::Real>& sp = mf_psource.array(mfi);
+        const amrex::Array4<const int>& iblank = mf_iblank.const_array(mfi);
+        // Only faces with iblank = -1 on both sides can have nonzero flux
+        amrex::ParallelFor(
+            xbx, ybx, zbxp1,
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                bool zero_all = (iblank(i - 1, j, k) + iblank(i, j, k) > -2);
+                for (int n = 0; n < ncompx; ++n) {
+                    fx(i, j, k, n) *= zero_all ? 0. : 1.;
+                }
+            },
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                bool zero_all = (iblank(i, j - 1, k) + iblank(i, j, k) > -2);
+                for (int n = 0; n < ncompy; ++n) {
+                    fy(i, j, k, n) *= zero_all ? 0. : 1.;
+                }
+            },
+            [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                bool zero_all = (iblank(i, j, k - 1) + iblank(i, j, k) > -2);
+                for (int n = 0; n < ncompz; ++n) {
+                    fz(i, j, k, n) *= zero_all ? 0. : 1.;
+                }
+            });
+        // With knowledge of fluxes, compute pressure source term
+        amrex::Box const& nbx = mfi.nodaltilebox();
+        amrex::ParallelFor(
+            nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // To get the local density flux corresponding to the node,
+                // average the surrounding fluxes. Fluxes in AMR-Wind domain are
+                // zeroed, and they are flagged as 0 to avoid dividing by the
+                // wrong number in the average. For AMR-Wind nodes (not
+                // bordering nalu-wind cells), psource = 0 / (0 + tiny)
+                const amrex::Real rho_flux =
+                    (fz(i, j, k, 1) + fz(i - 1, j, k, 1) + fz(i, j - 1, k, 1) +
+                     fz(i - 1, j - 1, k, 1)) /
+                    (fz(i, j, k, 8) + fz(i - 1, j, k, 8) + fz(i, j - 1, k, 8) +
+                     fz(i - 1, j - 1, k, 8) + 1e-100);
+                // Pressure source is proportional to local density source
+                // (because more or less density is above node)
+                sp(i, j, k) = rho_flux * grav_z;
+            });
+    }
+}
+
+void apply_fluxes(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_psource,
+    amrex::MultiFab& mf_vof,
+    amrex::MultiFab& mf_dens,
+    amrex::MultiFab& mf_vel,
+    amrex::MultiFab& mf_gp,
+    amrex::MultiFab& mf_pressure,
+    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx,
+    amrex::Real ptfac,
+    const bool mod_gp)
+{
+    constexpr amrex::Real tiny = 1e-12;
+    // Pseudo-time factor
+    // const amrex::Real ptfac = 1.0;
+    for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        amrex::Box const& nbx = mfi.nodaltilebox();
+        const amrex::Array4<amrex::Real>& fx = mf_fx.array(mfi);
+        const amrex::Array4<amrex::Real>& fy = mf_fy.array(mfi);
+        const amrex::Array4<amrex::Real>& fz = mf_fz.array(mfi);
+        const amrex::Array4<amrex::Real>& vof = mf_vof.array(mfi);
+        const amrex::Array4<amrex::Real>& dens = mf_dens.array(mfi);
+        const amrex::Array4<amrex::Real>& vel = mf_vel.array(mfi);
+        const amrex::Array4<amrex::Real>& gp = mf_gp.array(mfi);
+        const amrex::Array4<amrex::Real>& p = mf_pressure.array(mfi);
+        const amrex::Array4<amrex::Real>& sp = mf_psource.array(mfi);
+        amrex::ParallelFor(
+            vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real oldvof = vof(i, j, k);
+                const amrex::Real olddens = dens(i, j, k);
+                vof(i, j, k) += ptfac * ((fx(i + 1, j, k, 0) - fx(i, j, k, 0)) +
+                                         (fy(i, j + 1, k, 0) - fy(i, j, k, 0)) +
+                                         (fz(i, j, k + 1, 0) - fz(i, j, k, 0)));
+                dens(i, j, k) +=
+                    ptfac * ((fx(i + 1, j, k, 1) - fx(i, j, k, 1)) +
+                             (fy(i, j + 1, k, 1) - fy(i, j, k, 1)) +
+                             (fz(i, j, k + 1, 1) - fz(i, j, k, 1)));
+                vel(i, j, k, 0) =
+                    1.0 / dens(i, j, k) *
+                    (olddens * vel(i, j, k, 0) +
+                     ptfac * ((fx(i + 1, j, k, 2) - fx(i, j, k, 2)) +
+                              (fy(i, j + 1, k, 2) - fy(i, j, k, 2)) +
+                              (fz(i, j, k + 1, 2) - fz(i, j, k, 2))));
+                vel(i, j, k, 1) =
+                    1.0 / dens(i, j, k) *
+                    (olddens * vel(i, j, k, 1) +
+                     ptfac * ((fx(i + 1, j, k, 3) - fx(i, j, k, 3)) +
+                              (fy(i, j + 1, k, 3) - fy(i, j, k, 3)) +
+                              (fz(i, j, k + 1, 3) - fz(i, j, k, 3))));
+                vel(i, j, k, 2) =
+                    1.0 / dens(i, j, k) *
+                    (olddens * vel(i, j, k, 2) +
+                     ptfac * ((fx(i + 1, j, k, 4) - fx(i, j, k, 4)) +
+                              (fy(i, j + 1, k, 4) - fy(i, j, k, 4)) +
+                              (fz(i, j, k + 1, 4) - fz(i, j, k, 4))));
+
+                if (mod_gp) {
+
+                    gp(i, j, k, 0) +=
+                        ptfac * ((fx(i + 1, j, k, 5) - fx(i, j, k, 5)) +
+                                 (fy(i, j + 1, k, 5) - fy(i, j, k, 5)) +
+                                 (fz(i, j, k + 1, 5) - fz(i, j, k, 5)));
+                    gp(i, j, k, 1) +=
+                        ptfac * ((fx(i + 1, j, k, 6) - fx(i, j, k, 6)) +
+                                 (fy(i, j + 1, k, 6) - fy(i, j, k, 6)) +
+                                 (fz(i, j, k + 1, 6) - fz(i, j, k, 6)));
+                    gp(i, j, k, 2) +=
+                        ptfac * ((fx(i + 1, j, k, 7) - fx(i, j, k, 7)) +
+                                 (fy(i, j + 1, k, 7) - fy(i, j, k, 7)) +
+                                 (fz(i, j, k + 1, 7) - fz(i, j, k, 7)));
+                }
+
+                /*if (gp(i, j, k, 2) > 100.0)
+                    std::cout
+                        << i << " " << j << " " << k << " gp fluxes "
+                        << fx(i + 1, j, k, 7) << " " << -fx(i, j, k, 7) << " "
+                        << fy(i, j + 1, k, 7) << " " << -fy(i, j, k, 7) << " "
+                        << fz(i, j, k + 1, 7) << " " << -fz(i, j, k, 7)
+                        << " density fluxes " << fx(i + 1, j, k, 1) << " "
+                        << -fx(i, j, k, 1) << " " << fy(i, j + 1, k, 1) << " "
+                        << -fy(i, j, k, 1) << " " << fz(i, j, k + 1, 1) << " "
+                        << -fz(i, j, k, 1) << std::endl;*/
+
+                // Ensure vof is bounded
+                vof(i, j, k) =
+                    vof(i, j, k) < tiny
+                        ? 0.0
+                        : (vof(i, j, k) > 1. - tiny ? 1. : vof(i, j, k));
+                // Density is corrected later
+            });
+
+        // Apply pressure source/sink
+        amrex::ParallelFor(
+            nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // Move hydrostatic pressure around due to density flux
+                p(i, j, k) += ptfac * sp(i, j, k) * dx[2];
+            });
+    }
+}
+
+amrex::Real measure_convergence(
+    amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz)
+{
+    // Get the maximum flux magnitude, but just for vof fluxes
+    const amrex::Real err_fx = amrex::ReduceMax(
+        mf_fx, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fx) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fx(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err_fy = amrex::ReduceMax(
+        mf_fy, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fy) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fy(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err_fz = amrex::ReduceMax(
+        mf_fz, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx,
+            amrex::Array4<amrex::Real const> const& fz) -> amrex::Real {
+            amrex::Real err_fab = -1.0;
+            amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept {
+                err_fab = amrex::max(err_fab, std::abs(fz(i, j, k, 0)));
+            });
+            return err_fab;
+        });
+    const amrex::Real err = amrex::max(err_fx, amrex::max(err_fy, err_fz));
+    return err;
+}
+
+amrex::Real calculate_pseudo_dt_flux(
+    amrex::MultiFab& mf_fx,
+    amrex::MultiFab& mf_fy,
+    amrex::MultiFab& mf_fz,
+    amrex::MultiFab& mf_vof,
+    amrex::Real tol)
+{
+    // Get the maximum flux magnitude, but just for vof fluxes
+    const amrex::Real pdt_fx = amrex::ReduceMin(
+        mf_fx, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fx,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fx(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell i
+                    pdt_lim = vof(i, j, k) / fx(i, j, k, 0);
+                } else if (fx(i, j, k, 0) < -tol && vof(i - 1, j, k) > tol) {
+                    // VOF is removed from cell i-1
+                    pdt_lim = vof(i - 1, j, k) / -fx(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt_fy = amrex::ReduceMin(
+        mf_fy, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fy,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fy(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell j
+                    pdt_lim = vof(i, j, k) / fy(i, j, k, 0);
+                } else if (fy(i, j, k, 0) < -tol && vof(i, j - 1, k) > tol) {
+                    // VOF is removed from cell j-1
+                    pdt_lim = vof(i, j - 1, k) / -fy(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt_fz = amrex::ReduceMin(
+        mf_fz, mf_vof, 0,
+        [=] AMREX_GPU_HOST_DEVICE(
+            amrex::Box const& bx, amrex::Array4<amrex::Real const> const& fz,
+            amrex::Array4<amrex::Real const> const& vof) -> amrex::Real {
+            amrex::Real pdt_fab = 1.0;
+            amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept {
+                amrex::Real pdt_lim = 1.0;
+                if (fz(i, j, k, 0) > tol && vof(i, j, k) > tol) {
+                    // VOF is removed from cell k
+                    pdt_lim = vof(i, j, k) / fz(i, j, k, 0);
+                } else if (fz(i, j, k, 0) < -tol && vof(i, j, k - 1) > tol) {
+                    // VOF is removed from cell k-1
+                    pdt_lim = vof(i, j, k - 1) / -fz(i, j, k, 0);
+                }
+                pdt_fab = amrex::min(pdt_fab, pdt_lim);
+            });
+            return pdt_fab;
+        });
+    const amrex::Real pdt = amrex::min(pdt_fx, amrex::min(pdt_fy, pdt_fz));
+    return pdt;
+}
+
+// for debugging
+void equate_field(amrex::MultiFab& mf_dest, amrex::MultiFab& mf_src)
+{
+    for (amrex::MFIter mfi(mf_dest); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        const amrex::Array4<const amrex::Real>& src = mf_src.const_array(mfi);
+        const amrex::Array4<amrex::Real>& dest = mf_dest.array(mfi);
+        amrex::ParallelFor(
+            vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                dest(i, j, k) = std::sqrt(
+                    src(i, j, k, 0) * src(i, j, k, 0) +
+                    src(i, j, k, 1) * src(i, j, k, 1) +
+                    src(i, j, k, 2) * src(i, j, k, 2));
+            });
+    }
+}
+
+void harmonize_vof(
+    amrex::MultiFab& mf_vof_target,
+    amrex::MultiFab& mf_vof_original,
+    amrex::iMultiFab& mf_iblank)
+{
+    for (amrex::MFIter mfi(mf_vof_target); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        const amrex::Array4<amrex::Real>& tg_vof = mf_vof_target.array(mfi);
+        const amrex::Array4<const amrex::Real>& og_vof =
+            mf_vof_original.const_array(mfi);
+        const amrex::Array4<const int>& iblank = mf_iblank.const_array(mfi);
+        amrex::ParallelFor(
+            vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // Replace amr-wind vof values with originals
+                if (iblank(i, j, k) != -1) {
+                    tg_vof(i, j, k) = og_vof(i, j, k);
+                }
+            });
+    }
+}
+
+void replace_gradp_hs(
+    amrex::MultiFab& mf_gp,
+    amrex::MultiFab& mf_density,
+    amrex::iMultiFab& mf_iblank,
+    const amrex::Real grav_z)
+{
+    for (amrex::MFIter mfi(mf_gp); mfi.isValid(); ++mfi) {
+        const auto& gbx = mfi.growntilebox();
+        const amrex::Array4<amrex::Real>& gp = mf_gp.array(mfi);
+        const amrex::Array4<const amrex::Real>& rho =
+            mf_density.const_array(mfi);
+        const amrex::Array4<const int>& iblank = mf_iblank.const_array(mfi);
+        amrex::ParallelFor(
+            gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // Replace amr-wind vof values with originals
+                if (iblank(i, j, k) == -1) {
+                    gp(i, j, k, 0) = 0.;
+                    gp(i, j, k, 1) = 0.;
+                    gp(i, j, k, 2) = rho(i, j, k) * grav_z;
+                }
+            });
+    }
+}
+
+void replace_gradp(
+    amrex::MultiFab& mf_gp,
+    amrex::MultiFab& mf_gp0,
+    amrex::iMultiFab& mf_iblank)
+{
+    for (amrex::MFIter mfi(mf_gp); mfi.isValid(); ++mfi) {
+        const auto& gbx = mfi.growntilebox();
+        const amrex::Array4<amrex::Real>& gp = mf_gp.array(mfi);
+        const amrex::Array4<const amrex::Real>& gp0 = mf_gp0.const_array(mfi);
+        const amrex::Array4<const int>& iblank = mf_iblank.const_array(mfi);
+        amrex::ParallelFor(
+            gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // Replace amr-wind vof values with originals
+                if (iblank(i, j, k) == -1) {
+                    gp(i, j, k, 0) = gp0(i, j, k, 0);
+                    gp(i, j, k, 1) = gp0(i, j, k, 1);
+                    gp(i, j, k, 2) = gp0(i, j, k, 2);
+                }
+            });
+    }
+}
+
+void apply_pressure_gradient(
+    amrex::MultiFab& mf_vel,
+    amrex::MultiFab& mf_density,
+    amrex::MultiFab& mf_gp,
+    const amrex::Real scaling_factor)
+{
+    for (amrex::MFIter mfi(mf_gp); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        const amrex::Array4<amrex::Real>& vel = mf_vel.array(mfi);
+        const amrex::Array4<const amrex::Real>& rho =
+            mf_density.const_array(mfi);
+        const amrex::Array4<const amrex::Real>& gp = mf_gp.const_array(mfi);
+        amrex::ParallelFor(
+            vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real soverrho = scaling_factor / rho(i, j, k);
+                vel(i, j, k, 0) -= gp(i, j, k, 0) * soverrho;
+                vel(i, j, k, 1) -= gp(i, j, k, 1) * soverrho;
+                vel(i, j, k, 2) -= gp(i, j, k, 2) * soverrho;
+            });
+    }
+}
+
+void SharpenNaluDataDiscrete(
+    CFDSim& sim,
+    const int nit,
+    const amrex::Real tol,
+    const int ncconv,
+    const amrex::Real rlscale,
+    const amrex::Real margin,
+    const amrex::Real proc_tgvof_tol,
+    const bool sh_pressure,
+    const bool sh_gradp)
+{
+    auto& repo = sim.repo();
+    auto nlevels = repo.num_active_levels();
+    auto geom = sim.mesh().Geom();
+
+    // Check for vof variable -- sharpening only required with VOF
+    bool yesVOF = repo.field_exists("vof");
+    if (!yesVOF) return;
+
+    // Get phase densities
+    auto& mphase = sim.physics_manager().get<MultiPhase>();
+    const amrex::Real rhol = mphase.rho1();
+    const amrex::Real rhog = mphase.rho2();
+    // If pressure is not to be modified, set gravity to 0
+    const amrex::Real grav_z = sh_pressure ? mphase.gravity()[2] : 0.0;
+
+    // Get blanking for cells
+    auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified
+    auto& vof = repo.get_field("vof");
+    auto& levelset = repo.get_field("levelset");
+    auto& rho = repo.get_field("density");
+    auto& velocity = repo.get_field("velocity");
+    auto& gp = repo.get_field("gp");
+    auto& p = repo.get_field("p");
+
+    // Create scratch fields for fluxes - 5 components are vof, density, and
+    // 3 of velocity (flux_z has 6th component for flag in p_src calc)
+    auto flux_x = repo.create_scratch_field(8, 0, amr_wind::FieldLoc::XFACE);
+    auto flux_y = repo.create_scratch_field(8, 0, amr_wind::FieldLoc::YFACE);
+    auto flux_z = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::ZFACE);
+    // Create scratch field for pressure source term
+    auto p_src = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::NODE);
+    // Create scratch field for approximate signed distance function and grad
+    // (components 0-2 are gradient, 3 is asdf)
+    auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
+
+    auto target_vof = repo.create_scratch_field(1, vof.num_grow()[0]);
+
+    // Create levelset field
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // Thickness used here is user parameter, whatever works best
+        auto dx = (geom[lev]).CellSizeArray();
+        const amrex::Real i_th = rlscale * std::cbrt(dx[0] * dx[1] * dx[2]);
+
+        // Populate approximate signed distance function
+        populate_psi(levelset(lev), vof(lev), i_th);
+    }
+
+    // Convert levelset to vof to get target_vof
+    mphase.levelset2vof(iblank_cell, *target_vof);
+
+    // Process target vof for tiny margins from single-phase
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // A tolerance of 0 should do nothing
+        process_vof((*target_vof)(lev), proc_tgvof_tol);
+    }
+
+    // Replace vof with original values in amr domain
+    for (int lev = 0; lev < nlevels; ++lev) {
+        harmonize_vof((*target_vof)(lev), vof(lev), iblank_cell(lev));
+    }
+
+    // Put fluxes in vector for averaging down during iterations
+    amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
+        repo.num_active_levels());
+    for (int lev = 0; lev < nlevels; ++lev) {
+        fluxes[lev][0] = &(*flux_x)(lev);
+        fluxes[lev][1] = &(*flux_y)(lev);
+        fluxes[lev][2] = &(*flux_z)(lev);
+    }
+
+    // Pseudo-time loop
+    amrex::Real err = 100.0 * tol;
+    int n = 0;
+    while (n < nit && err > tol) {
+        // Increment step counter
+        ++n;
+        // Determine if convergence error is calculated this step
+        bool cconv = n % ncconv == 0;
+        // Zero error if being calculated this step
+        err = cconv ? 0.0 : err;
+
+        // Maximum possible value of pseudo time factor (dtau)
+        amrex::Real ptfac = 1.0;
+
+        // Maximum pseudoCFL, 0.5 seems to work well
+        const amrex::Real pCFL = 0.5;
+
+        for (int lev = 0; lev < nlevels; ++lev) {
+            // Populate normal vector
+            populate_normal_vector(
+                (*normal_vec)(lev), vof(lev), iblank_cell(lev));
+
+            // Sharpening fluxes for vof, density, and momentum
+            populate_sharpen_discrete_fluxes(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
+                (*target_vof)(lev), (*normal_vec)(lev), velocity(lev), gp(lev),
+                rho(lev), margin, rhol, rhog);
+
+            // Process fluxes and get pressure source term
+            process_fluxes_calc_src(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
+                iblank_cell(lev), grav_z);
+
+            // Measure convergence to determine if loop can stop
+            if (cconv) {
+                // Update error at specified interval of steps
+                const amrex::Real err_lev = measure_convergence(
+                    (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev));
+                err = amrex::max(err, err_lev);
+            }
+        }
+
+        // Average down fluxes across levels for consistency
+        for (int lev = nlevels - 1; lev > 0; --lev) {
+            amrex::IntVect rr =
+                geom[lev].Domain().size() / geom[lev - 1].Domain().size();
+            amrex::average_down_faces(
+                GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr,
+                geom[lev - 1]);
+        }
+
+        // Get pseudo dt (dtau)
+        for (int lev = 0; lev < nlevels; ++lev) {
+            // Compare vof fluxes to vof in source cells
+            // Convergence tolerance determines what size of fluxes matter
+            const amrex::Real ptfac_lev = calculate_pseudo_dt_flux(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev), tol);
+            ptfac = amrex::min(ptfac, ptfac_lev);
+        }
+        amrex::ParallelDescriptor::ReduceRealMin(ptfac);
+
+        // Conform pseudo dt (dtau) to pseudo CFL
+        ptfac = pCFL * ptfac;
+
+        // Apply fluxes
+        for (int lev = 0; lev < nlevels; ++lev) {
+            // Height of cell is needed for pressure source
+            auto dx = (geom[lev]).CellSizeArray();
+
+            apply_fluxes(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
+                vof(lev), rho(lev), velocity(lev), gp(lev), p(lev), dx, ptfac,
+                sh_gradp);
+        }
+
+        // sim.io_manager().write_plot_file();
+        // sim.time().increment_timestep();
+
+        // Fillpatch for ghost cells
+        vof.fillpatch(sim.time().current_time());
+        velocity.fillpatch(sim.time().current_time());
+
+        // Update density (fillpatch built in)
+        mphase.set_density_via_vof();
+
+        // Ensure that err is same across processors
+        if (cconv) {
+            amrex::ParallelDescriptor::ReduceRealMax(err);
+        }
+
+        amrex::Print() << "sharpen step " << n << " " << err << " " << tol
+                       << std::endl;
+    }
+
+    // Purely for debugging via visualization, should be removed later
+    // Currently set up to overwrite the levelset field (not used as time
+    // evolves) with the post-sharpening velocity magnitude
+    for (int lev = 0; lev < nlevels; ++lev) {
+        equate_field(levelset(lev), velocity(lev));
+    }
+
+    // Copy to old
+    amr_wind::field_ops::copy(
+        vof.state(amr_wind::FieldState::Old), vof, 0, 0, vof.num_comp(),
+        vof.num_grow());
+    amr_wind::field_ops::copy(
+        rho.state(amr_wind::FieldState::Old), rho, 0, 0, rho.num_comp(),
+        rho.num_grow());
+    amr_wind::field_ops::copy(
+        velocity.state(amr_wind::FieldState::Old), velocity, 0, 0,
+        velocity.num_comp(), velocity.num_grow());
+
+    // Fillpatch for pressure to make sure pressure stencil has all points
+    p.fillpatch(sim.time().current_time());
+}
+
+void ReplaceMaskedGradP(CFDSim& sim)
+{
+    auto& repo = sim.repo();
+    auto nlevels = repo.num_active_levels();
+    auto geom = sim.mesh().Geom();
+
+    // Get gravity
+    auto& mphase = sim.physics_manager().get<MultiPhase>();
+    const amrex::Real grav_z = mphase.gravity()[2];
+
+    // Get blanking for cells
+    auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified or used
+    auto& rho = repo.get_field("density");
+    auto& gp = repo.get_field("gp");
+
+    // Replace initial gp with best guess (hydrostatic)
+    for (int lev = 0; lev < nlevels; ++lev) {
+        replace_gradp_hs(gp(lev), rho(lev), iblank_cell(lev), grav_z);
+    }
+}
+
+void ReapplyModifiedGradP(CFDSim& sim)
+{
+    auto& repo = sim.repo();
+    auto nlevels = repo.num_active_levels();
+    auto geom = sim.mesh().Geom();
+
+    // Get gravity
+    auto& mphase = sim.physics_manager().get<MultiPhase>();
+    const amrex::Real grav_z = mphase.gravity()[2];
+
+    // Get timestep
+    const amrex::Real dt = sim.time().deltaT();
+
+    // Get blanking for cells
+    auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified or used
+    auto& vel = repo.get_field("velocity");
+    auto& rho = repo.get_field("density");
+    auto& gp = repo.get_field("gp");
+
+    // For iblanked cells, replace gp with original gp, to get original vel
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // Remove pressure gradient term
+        apply_pressure_gradient(vel(lev), rho(lev), gp(lev), -dt);
+        // Modify pressure gradient
+        replace_gradp_hs(gp(lev), rho(lev), iblank_cell(lev), grav_z);
+        // Reapply pressure gradient term
+        apply_pressure_gradient(vel(lev), rho(lev), gp(lev), dt);
+    }
+}
+
+void CopyGradP(ScratchField& gp_copy, Field& gp, const int nlevels)
+{
+    // For iblanked cells, replace gp with original gp, to get original vel
+    for (int lev = 0; lev < nlevels; ++lev) {
+        amrex::MultiFab::Copy(
+            gp_copy(lev), gp(lev), 0, 0, gp(lev).nComp(), gp(lev).nGrow());
+    }
+}
+
+void ReapplyOversetGradP(ScratchField& gp_copy, CFDSim& sim)
+{
+    auto& repo = sim.repo();
+    auto nlevels = repo.num_active_levels();
+
+    // Get timestep
+    const amrex::Real dt = sim.time().deltaT();
+
+    // Get blanking for cells
+    auto& iblank_cell = repo.get_int_field("iblank_cell");
+
+    // Get fields that will be modified or used
+    auto& vel = repo.get_field("velocity");
+    auto& rho = repo.get_field("density");
+    auto& gp = repo.get_field("gp");
+
+    // For iblanked cells, replace gp with original gp, to get original vel
+    for (int lev = 0; lev < nlevels; ++lev) {
+        // Remove pressure gradient term
+        apply_pressure_gradient(vel(lev), rho(lev), gp(lev), -dt);
+        // Modify pressure gradient
+        replace_gradp(gp(lev), gp_copy(lev), iblank_cell(lev));
+        // Reapply pressure gradient term
+        apply_pressure_gradient(vel(lev), rho(lev), gp(lev), dt);
+    }
+}
+
+} // namespace overset
+} // namespace amr_wind
+
+#endif
\ No newline at end of file
diff --git a/amr-wind/physics/multiphase/CMakeLists.txt b/amr-wind/physics/multiphase/CMakeLists.txt
index e6fe108096..0c2f5b7406 100644
--- a/amr-wind/physics/multiphase/CMakeLists.txt
+++ b/amr-wind/physics/multiphase/CMakeLists.txt
@@ -9,4 +9,5 @@ target_sources(${amr_wind_lib_name}
   SloshingTank.cpp
   RainDrop.cpp
   BreakingWaves.cpp
+  FuzzyInterface.cpp
   )
\ No newline at end of file
diff --git a/amr-wind/physics/multiphase/FuzzyInterface.H b/amr-wind/physics/multiphase/FuzzyInterface.H
new file mode 100644
index 0000000000..53d1420700
--- /dev/null
+++ b/amr-wind/physics/multiphase/FuzzyInterface.H
@@ -0,0 +1,55 @@
+#ifndef FuzzyInterface_H
+#define FuzzyInterface_H
+
+#include "amr-wind/core/Physics.H"
+#include "amr-wind/core/Field.H"
+
+/** Multiphase 3D physics
+ *  \ingroup multiphase_physics
+ *
+ */
+
+namespace amr_wind {
+
+class FuzzyInterface : public Physics::Register<FuzzyInterface>
+{
+public:
+    static std::string identifier() { return "FuzzyInterface"; }
+
+    explicit FuzzyInterface(CFDSim& sim);
+
+    ~FuzzyInterface() override = default;
+
+    //! Initialize the levelset and velocity fields for Sloshing Tank
+    //! simulations
+    void initialize_fields(int level, const amrex::Geometry& geom) override;
+
+    void post_init_actions() override {}
+
+    void post_regrid_actions() override {}
+
+    void pre_advance_work() override {}
+
+    void post_advance_work() override {}
+
+private:
+    Field& m_velocity;
+    Field& m_vof;
+    Field& m_pressure;
+
+    //! Initial zero-level free-surface water depth
+    amrex::Real m_waterlevel{0.0};
+    //! Divide vof between smooth and discrete at x
+    amrex::Real m_lx_vofjump{-0.25};
+    amrex::Real m_hx_vofjump{0.25};
+    // Thickness of interface in smooth section
+    amrex::Real m_intf_th{0.0};
+
+    //! Stuff to get from MultiPhase physics
+    amrex::Real m_rho1{1000.};
+    amrex::Real m_rho2{1.};
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
+};
+
+} // namespace amr_wind
+#endif /* FuzzyInterface_H */
diff --git a/amr-wind/physics/multiphase/FuzzyInterface.cpp b/amr-wind/physics/multiphase/FuzzyInterface.cpp
new file mode 100644
index 0000000000..3f6b8a37b6
--- /dev/null
+++ b/amr-wind/physics/multiphase/FuzzyInterface.cpp
@@ -0,0 +1,121 @@
+#include "amr-wind/physics/multiphase/FuzzyInterface.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
+#include "amr-wind/CFDSim.H"
+#include "AMReX_ParmParse.H"
+#include "amr-wind/fvm/gradient.H"
+#include "amr-wind/core/field_ops.H"
+
+namespace amr_wind {
+
+FuzzyInterface::FuzzyInterface(CFDSim& sim)
+    : m_velocity(sim.repo().get_field("velocity"))
+    , m_vof(sim.repo().get_field("vof"))
+    , m_pressure(sim.repo().get_field("p"))
+{
+    amrex::ParmParse pp(identifier());
+    pp.query("water_level", m_waterlevel);
+    pp.get("lo_x_vofjump", m_lx_vofjump);
+    pp.get("hi_x_vofjump", m_hx_vofjump);
+    pp.query("interface_thickness", m_intf_th);
+
+    amrex::ParmParse pp_multiphase("MultiPhase");
+    pp_multiphase.add("water_level", m_waterlevel);
+
+    auto& mphase = sim.physics_manager().get<MultiPhase>();
+    m_rho1 = mphase.rho1();
+    m_rho2 = mphase.rho2();
+    m_gravity = mphase.gravity();
+
+    // This physics case specifies vof directly
+    mphase.set_intf_init_method_to_vof();
+}
+
+/** Initialize the velocity and levelset fields at the beginning of the
+ *  simulation.
+ *
+ */
+void FuzzyInterface::initialize_fields(int level, const amrex::Geometry& geom)
+{
+    auto& velocity = m_velocity(level);
+    velocity.setVal(0.0);
+
+    auto& volfrac = m_vof(level);
+    auto& pressure = m_pressure(level);
+    const auto& dx = geom.CellSizeArray();
+    const auto& problo = geom.ProbLoArray();
+    const auto& probhi = geom.ProbHiArray();
+    const amrex::Real water_level = m_waterlevel;
+    const amrex::Real i_th = m_intf_th;
+    const amrex::Real lx_vj = m_lx_vofjump;
+    const amrex::Real hx_vj = m_hx_vofjump;
+    const amrex::Real rho1 = m_rho1;
+    const amrex::Real rho2 = m_rho2;
+    const amrex::Real grav_z = m_gravity[2];
+
+    for (amrex::MFIter mfi(volfrac); mfi.isValid(); ++mfi) {
+        const auto& vbx = mfi.validbox();
+        auto vof = volfrac.array(mfi);
+        auto p = pressure.array(mfi);
+
+        amrex::ParallelFor(
+            vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real x = problo[0] + (i + 0.5) * dx[0];
+                const amrex::Real z = problo[2] + (k + 0.5) * dx[2];
+                const amrex::Real zbtm = problo[2] + (k + 0.0) * dx[2];
+                const amrex::Real vof_sharp =
+                    std::max(0.0, std::min(1.0, (water_level - zbtm) / dx[2]));
+                const amrex::Real vof_smooth =
+                    -0.5 * (std::erf((z - water_level) / m_intf_th) + 1.0) +
+                    1.0;
+                if (x < lx_vj || x > hx_vj) {
+                    // Sharp vof
+                    vof(i, j, k) = vof_sharp;
+                } else {
+                    vof(i, j, k) = vof_smooth;
+                }
+            });
+
+        amrex::Box const& nbx = mfi.grownnodaltilebox();
+        amrex::ParallelFor(
+            nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // For pressure nodes, no offset
+                const amrex::Real x = problo[0] + i * dx[0];
+                const amrex::Real y = problo[1] + j * dx[1];
+                const amrex::Real z = problo[2] + k * dx[2];
+
+                // Sharp interpretation
+                amrex::Real ih_g =
+                    amrex::max(0.0, amrex::min(probhi[2] - water_level, probhi[2] - z));
+                amrex::Real ih_l =
+                    amrex::max(0.0, amrex::min(water_level - z, water_level - problo[2]));
+                const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
+                const amrex::Real psharp = - irho * grav_z;
+
+                // Smooth interpretation
+                const double x0 = (z - water_level) / i_th;
+                const double x1 = (probhi[2] - water_level) / i_th;
+                const double int_erf0 =
+                    (x0 * std::erf(x0) +
+                     std::exp(-std::pow(x0, 2)) / std::sqrt(M_PI)) *
+                    i_th;
+                const double int_erf1 =
+                    (x1 * std::erf(x1) +
+                     std::exp(-std::pow(x1, 2)) / std::sqrt(M_PI)) *
+                    i_th;
+                const double int_vof0 = 0.5 * z - 0.5 * int_erf0;
+                const double int_vof1 = 0.5 * probhi[2] - 0.5 * int_erf1;
+                const double int_rho0 =
+                    rho1 * int_vof0 + rho2 * (z - int_vof0);
+                const double int_rho1 =
+                    rho1 * int_vof1 + rho2 * (probhi[2] - int_vof1);
+
+                // g * integral(rho)dz
+                const amrex::Real psmooth = -(int_rho1 - int_rho0) * grav_z;
+
+                // Populate pressure
+                p(i, j, k) = (x < lx_vj || x > hx_vj)? psharp : psmooth;
+            });
+    }
+}
+
+} // namespace amr_wind
diff --git a/amr-wind/physics/multiphase/MultiPhase.H b/amr-wind/physics/multiphase/MultiPhase.H
index d31affbf18..f7b80fe451 100644
--- a/amr-wind/physics/multiphase/MultiPhase.H
+++ b/amr-wind/physics/multiphase/MultiPhase.H
@@ -3,6 +3,8 @@
 
 #include "amr-wind/core/Physics.H"
 #include "amr-wind/core/Field.H"
+#include "amr-wind/core/IntField.H"
+#include "amr-wind/core/ScratchField.H"
 
 /** Multiphase physics
  *
@@ -51,8 +53,20 @@ public:
 
     void levelset2vof();
 
+    void levelset2vof(IntField& iblank_cell, ScratchField& vof_scr);
+
     void favre_filtering();
 
+    void set_intf_init_method_to_vof()
+    {
+        m_interface_init_method = InterfaceCapturingMethod::VOF;
+    }
+
+    void set_intf_init_method_to_ls()
+    {
+        m_interface_init_method = InterfaceCapturingMethod::LS;
+    };
+
     amrex::Real volume_fraction_sum();
 
     amrex::Real momentum_sum(int n);
@@ -63,6 +77,8 @@ public:
 
     amrex::Real rho2() const { return m_rho2; }
 
+    amrex::Vector<amrex::Real> gravity() const { return m_gravity; }
+
 private:
     const CFDSim& m_sim;
 
@@ -87,6 +103,8 @@ private:
     // Reconstructing true pressure field at end of timestep
     bool is_ptrue{false};
 
+    bool m_init_p{false};
+
     // Info to create rho0
     amrex::Real water_level0{0.0};
     // Info to reconstruct true pressure
@@ -101,6 +119,9 @@ private:
     InterfaceCapturingMethod m_interface_capturing_method =
         InterfaceCapturingMethod::VOF;
 
+    InterfaceCapturingMethod m_interface_init_method =
+        InterfaceCapturingMethod::LS;
+
     // Verbose flag for multiphase
     int m_verbose{0};
 
diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp
index 7b85bb0008..e1a6b5daf5 100644
--- a/amr-wind/physics/multiphase/MultiPhase.cpp
+++ b/amr-wind/physics/multiphase/MultiPhase.cpp
@@ -66,6 +66,8 @@ MultiPhase::MultiPhase(CFDSim& sim)
                 "reference_pressure", 1, (*m_vof).num_grow()[0], 1);
         }
     }
+
+    pp_multiphase.query("initialize_pressure", m_init_p);
 }
 
 InterfaceCapturingMethod MultiPhase::interface_capturing_method()
@@ -78,13 +80,32 @@ void MultiPhase::post_init_actions()
 
     const auto& io_mgr = m_sim.io_manager();
     if (!io_mgr.is_restart()) {
-        switch (m_interface_capturing_method) {
-        case InterfaceCapturingMethod::VOF:
-            levelset2vof();
-            set_density_via_vof();
-            break;
+        // Different steps depending on how case is initialized, solved
+        switch (m_interface_init_method) {
         case InterfaceCapturingMethod::LS:
-            set_density_via_levelset();
+            switch (m_interface_capturing_method) {
+            case InterfaceCapturingMethod::VOF:
+                levelset2vof();
+                set_density_via_vof();
+                break;
+            case InterfaceCapturingMethod::LS:
+                set_density_via_levelset();
+                break;
+            };
+            break;
+
+        case InterfaceCapturingMethod::VOF:
+            switch (m_interface_capturing_method) {
+            case InterfaceCapturingMethod::VOF:
+                // VOF has been specified directly, just update density
+                set_density_via_vof();
+                break;
+            case InterfaceCapturingMethod::LS:
+                amrex::Abort(
+                    "Initialization failed. Case cannot be initialized with "
+                    "VOF and then simulated with Levelset");
+                break;
+            };
             break;
         };
     }
@@ -122,6 +143,20 @@ void MultiPhase::post_init_actions()
                 m_sim.mesh().Geom());
         }
     }
+
+    if (m_init_p && !is_wlev) {
+        amrex::Abort(
+            "Initialize pressure requested, but physics case does not "
+            "specify water level.");
+    }
+    // Make p field if both are specified
+    if (m_init_p && is_wlev) {
+        pp_multiphase.get("water_level", water_level0);
+        // Initialize rho0 field for perturbational density, pressure
+        auto& p = m_sim.repo().get_field("p");
+        hydrostatic::define_p0(
+            p, m_rho1, m_rho2, water_level0, m_gravity[2], m_sim.mesh().Geom());
+    }
 }
 
 void MultiPhase::post_regrid_actions()
@@ -163,7 +198,8 @@ void MultiPhase::post_advance_work()
 {
     switch (m_interface_capturing_method) {
     case InterfaceCapturingMethod::VOF:
-        // Compute and print the total volume fraction, momenta, and differences
+        // Compute and print the total volume fraction, momenta, and
+        // differences
         if (m_verbose > 0) {
             m_total_volfrac = volume_fraction_sum();
             amrex::Real mom_x = momentum_sum(0) - q0;
@@ -511,4 +547,82 @@ void MultiPhase::levelset2vof()
     (*m_vof).fillpatch(0.0);
 }
 
+void MultiPhase::levelset2vof(IntField& iblank_cell, ScratchField& vof_scr)
+{
+    const int nlevels = m_sim.repo().num_active_levels();
+    (*m_levelset).fillpatch(m_sim.time().current_time());
+    const auto& geom = m_sim.mesh().Geom();
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& levelset = (*m_levelset)(lev);
+        auto& vof = vof_scr(lev);
+        const auto& dx = geom[lev].CellSizeArray();
+
+        for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
+            const auto& vbx = mfi.validbox();
+            const amrex::Array4<amrex::Real>& phi = levelset.array(mfi);
+            const amrex::Array4<amrex::Real>& volfrac = vof.array(mfi);
+            const amrex::Array4<const int>& iblank =
+                iblank_cell(lev).const_array(mfi);
+            const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]);
+            amrex::ParallelFor(
+                vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    // Neumann of levelset across iblank boundaries
+                    int ibdy = 0;
+                    int jbdy = 0;
+                    int kbdy = 0;
+                    if (iblank(i, j, k) != iblank(i - 1, j, k)) {
+                        ibdy = -1;
+                    }
+                    if (iblank(i, j, k) != iblank(i, j - 1, k)) {
+                        jbdy = -1;
+                    }
+                    if (iblank(i, j, k) != iblank(i, j, k - 1)) {
+                        kbdy = -1;
+                    }
+                    // no cell should be isolated such that -1 and 1 are
+                    // needed
+                    if (iblank(i, j, k) != iblank(i + 1, j, k)) {
+                        ibdy = +1;
+                    }
+                    if (iblank(i, j, k) != iblank(i, j + 1, k)) {
+                        jbdy = +1;
+                    }
+                    if (iblank(i, j, k) != iblank(i, j, k + 1)) {
+                        kbdy = +1;
+                    }
+                    amrex::Real mx, my, mz;
+                    multiphase::youngs_fd_normal_neumann(
+                        i, j, k, ibdy, jbdy, kbdy, phi, mx, my, mz);
+                    mx = std::abs(mx / 32.);
+                    my = std::abs(my / 32.);
+                    mz = std::abs(mz / 32.);
+                    amrex::Real normL1 = (mx + my + mz);
+                    mx = mx / normL1;
+                    my = my / normL1;
+                    mz = mz / normL1;
+                    // Make sure that alpha is negative far away from the
+                    // interface
+                    amrex::Real alpha;
+                    if (phi(i, j, k) < -eps) {
+                        alpha = -1.0;
+                    } else {
+                        alpha = phi(i, j, k) / normL1;
+                        alpha = alpha + 0.5;
+                    }
+                    if (alpha >= 1.0) {
+                        volfrac(i, j, k) = 1.0;
+                    } else if (alpha <= 0.0) {
+                        volfrac(i, j, k) = 0.0;
+                    } else {
+                        volfrac(i, j, k) =
+                            multiphase::cut_volume(mx, my, mz, alpha, 0.0, 1.0);
+                    }
+                });
+        }
+    }
+    // Fill ghost and boundary cells before simulation begins
+    vof_scr.fillpatch(0.0);
+}
+
 } // namespace amr_wind
diff --git a/amr-wind/physics/multiphase/SloshingTank.H b/amr-wind/physics/multiphase/SloshingTank.H
index bf7b24b763..c6097798d9 100644
--- a/amr-wind/physics/multiphase/SloshingTank.H
+++ b/amr-wind/physics/multiphase/SloshingTank.H
@@ -35,6 +35,7 @@ public:
 private:
     Field& m_velocity;
     Field& m_levelset;
+    Field& m_pressure;
 
     //! Initial free surface amplitude magnitude
     amrex::Real m_amplitude{0.1};
@@ -44,6 +45,11 @@ private:
 
     //! Initial zero-level free-surface water depth
     amrex::Real m_waterlevel{0.0};
+
+    //! Stuff to get from MultiPhase physics
+    amrex::Real m_rho1{1000.};
+    amrex::Real m_rho2{1.};
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 };
 
 } // namespace amr_wind
diff --git a/amr-wind/physics/multiphase/SloshingTank.cpp b/amr-wind/physics/multiphase/SloshingTank.cpp
index 7558141c7c..43dfb02a27 100644
--- a/amr-wind/physics/multiphase/SloshingTank.cpp
+++ b/amr-wind/physics/multiphase/SloshingTank.cpp
@@ -1,4 +1,5 @@
 #include "amr-wind/physics/multiphase/SloshingTank.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
 #include "amr-wind/CFDSim.H"
 #include "AMReX_ParmParse.H"
 #include "amr-wind/fvm/gradient.H"
@@ -9,11 +10,17 @@ namespace amr_wind {
 SloshingTank::SloshingTank(CFDSim& sim)
     : m_velocity(sim.repo().get_field("velocity"))
     , m_levelset(sim.repo().get_field("levelset"))
+    , m_pressure(sim.repo().get_field("p"))
 {
     amrex::ParmParse pp(identifier());
     pp.query("amplitude", m_amplitude);
     pp.query("peak_enhance", m_kappa);
     pp.query("water_level", m_waterlevel);
+
+    const auto& mphase = sim.physics_manager().get<MultiPhase>();
+    m_rho1 = mphase.rho1();
+    m_rho2 = mphase.rho2();
+    m_gravity = mphase.gravity();
 }
 
 /** Initialize the velocity and levelset fields at the beginning of the
@@ -26,6 +33,7 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     velocity.setVal(0.0);
 
     auto& levelset = m_levelset(level);
+    auto& pressure = m_pressure(level);
     const auto& dx = geom.CellSizeArray();
     const auto& problo = geom.ProbLoArray();
     const auto& probhi = geom.ProbHiArray();
@@ -34,10 +42,14 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real water_level = m_waterlevel;
     const amrex::Real Lx = probhi[0] - problo[0];
     const amrex::Real Ly = probhi[1] - problo[1];
+    const amrex::Real rho1 = m_rho1;
+    const amrex::Real rho2 = m_rho2;
+    const amrex::Real grav_z = m_gravity[2];
 
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
+        auto p = pressure.array(mfi);
 
         amrex::ParallelFor(
             vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -51,6 +63,34 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
                                         std::pow(y - problo[1] - 0.5 * Ly, 2)));
                 phi(i, j, k) = z0 - z;
             });
+
+        amrex::Box const& nbx = mfi.grownnodaltilebox();
+        amrex::ParallelFor(
+            nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                // For pressure nodes, no offset
+                const amrex::Real x = problo[0] + i * dx[0];
+                const amrex::Real y = problo[1] + j * dx[1];
+                const amrex::Real z = problo[2] + k * dx[2];
+                const amrex::Real z0 =
+                    water_level +
+                    Amp * std::exp(
+                              -kappa * (std::pow(x - problo[0] - 0.5 * Lx, 2) +
+                                        std::pow(y - problo[1] - 0.5 * Ly, 2)));
+                // Liquid height
+                const amrex::Real hliq = z0 - problo[2];
+                // Hight of current position from base
+                const amrex::Real hz = z - problo[2];
+                // Integrated (top-down in z) phase heights to pressure node
+                amrex::Real ih_g = amrex::max(
+                    0.0, amrex::min(probhi[2] - z0, probhi[2] - z));
+                amrex::Real ih_l =
+                    amrex::max(0.0, amrex::min(z0 - z, z0 - problo[2]));
+                // Integrated rho at pressure node
+                const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
+
+                // Add term to reference pressure
+                p(i, j, k) = -irho * grav_z;
+            });
     }
 }
 
diff --git a/amr-wind/physics/multiphase/hydrostatic_ops.H b/amr-wind/physics/multiphase/hydrostatic_ops.H
index e0b74a6d9b..de7491e33f 100644
--- a/amr-wind/physics/multiphase/hydrostatic_ops.H
+++ b/amr-wind/physics/multiphase/hydrostatic_ops.H
@@ -46,15 +46,12 @@ static void define_p0(
             auto p0_arr = p0(lev).array(mfi);
             amrex::ParallelFor(
                 nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
-                    // Height of pressure node
-                    const amrex::Real hnode = k * dx[2];
-                    // Liquid height
-                    const amrex::Real hliq = wlev - problo[2];
+                    const amrex::Real z = problo[2] + k * dx[2];
                     // Integrated (top-down in z) phase heights to pressure node
                     amrex::Real ih_g = amrex::max(
-                        0.0, amrex::min(probhi[2] - hliq, probhi[2] - hnode));
+                        0.0, amrex::min(probhi[2] - wlev, probhi[2] - z));
                     amrex::Real ih_l = amrex::max(
-                        0.0, amrex::min(hliq - hnode, hliq - problo[2]));
+                        0.0, amrex::min(wlev - z, wlev - problo[2]));
                     // Integrated rho at pressure node
                     const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
 
diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp
index f4477547c6..d724f32de8 100644
--- a/amr-wind/projection/incflo_apply_nodal_projection.cpp
+++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp
@@ -333,7 +333,7 @@ void incflo::ApplyProjection(
     }
 
     // Setup masking for overset simulations
-    if (sim().has_overset()) {
+    if (sim().has_overset() && !m_disable_onodal) {
         auto& linop = nodal_projector->getLinOp();
         const auto& imask_node = repo().get_int_field("mask_node");
         for (int lev = 0; lev <= finest_level; ++lev) {
@@ -341,7 +341,7 @@ void incflo::ApplyProjection(
         }
     }
 
-    if (m_sim.has_overset()) {
+    if (m_sim.has_overset() && !m_disable_onodal) {
         auto phif = m_repo.create_scratch_field(1, 1, amr_wind::FieldLoc::NODE);
         if (incremental) {
             for (int lev = 0; lev <= finestLevel(); ++lev) {