diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp
index c251cda5b..82b7a346b 100644
--- a/Exec/DevTests/MovingTerrain/prob.cpp
+++ b/Exec/DevTests/MovingTerrain/prob.cpp
@@ -155,47 +155,57 @@ Problem::init_custom_terrain (const Geometry& geom,
                               MultiFab& z_phys_nd,
                               const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    Real Ampl        = parms.Ampl;
-    Real wavelength  = 100.;
-    Real kp          = 2.0 * PI / wavelength;
-    Real g           = CONST_GRAV;
-    Real omega       = std::sqrt(g * kp);
-
-    for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-
-            // Location of nodes
-            Real x = ii  * dx[0];
-
-            // Wave height
-            Real height = Ampl * std::sin(kp * x - omega * time);
-
-            // Populate terrain height
-            z_arr(i,j,k0) = height;
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        Real Ampl        = parms.Ampl;
+        Real wavelength  = 100.;
+        Real kp          = 2.0 * PI / wavelength;
+        Real g           = CONST_GRAV;
+        Real omega       = std::sqrt(g * kp);
+
+        for (MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+
+                // Location of nodes
+                Real x = ii  * dx[0];
+
+                // Wave height
+                Real height = Ampl * std::sin(kp * x - omega * time);
+
+                // Populate terrain height
+                z_arr(i,j,k0) = height;
+            });
+        }
     }
 }
diff --git a/Exec/RegTests/ParticlesOverWoA/prob.cpp b/Exec/RegTests/ParticlesOverWoA/prob.cpp
index 422823dde..b65838449 100644
--- a/Exec/RegTests/ParticlesOverWoA/prob.cpp
+++ b/Exec/RegTests/ParticlesOverWoA/prob.cpp
@@ -126,54 +126,64 @@ void
 Problem::init_custom_terrain(
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real num  = 8 * a * a * a;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-            // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
-
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-            // Real y = (jj  * dx[1] - ycen);
-
-            // WoA Hill in x-direction
-            Real height = num / (x*x + 4 * a * a);
-
-            // Populate terrain height
-            z_arr(i,j,k0) = height;
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real num  = 8 * a * a * a;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+                // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
+
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+                // Real y = (jj  * dx[1] - ycen);
+
+                // WoA Hill in x-direction
+                Real height = num / (x*x + 4 * a * a);
+
+                // Populate terrain height
+                z_arr(i,j,k0) = height;
+            });
+        }
     }
 }
diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp
index 948d8c9b2..f99f0e224 100644
--- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp
+++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp
@@ -87,54 +87,64 @@ void
 Problem::init_custom_terrain(
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real num  = 8 * a * a * a;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-
-            if(fabs(x)<a){
-                z_arr(i,j,k0) = pow(a*a - x*x, 0.5);
-            }
-            else{
-                z_arr(i,j,k0) = 0.0;
-            }
-
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real num  = 8 * a * a * a;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+
+                if(fabs(x)<a){
+                    z_arr(i,j,k0) = pow(a*a - x*x, 0.5);
+                }
+                else{
+                    z_arr(i,j,k0) = 0.0;
+                }
+
+            });
+        }
     }
 }
diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp
index af3fae86d..2c90b53ec 100644
--- a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp
+++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp
@@ -146,55 +146,65 @@ void
 Problem::init_custom_terrain (
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
 
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
 
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
 
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-            int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+                int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
 
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-            Real y = (jj  * dx[1] - ycen);
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+                Real y = (jj  * dx[1] - ycen);
 
-            if(std::pow(x*x + y*y, 0.5) < a){
-                z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5);
-            }
-            else{
-                z_arr(i,j,k0) = 0.0;
-            }
+                if(std::pow(x*x + y*y, 0.5) < a){
+                    z_arr(i,j,k0) = std::pow(a*a - x*x - y*y , 0.5);
+                }
+                else{
+                    z_arr(i,j,k0) = 0.0;
+                }
 
-        });
+            });
+        }
     }
 }
diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp
index 45c64c93d..82e89f3a1 100644
--- a/Exec/RegTests/WitchOfAgnesi/prob.cpp
+++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp
@@ -4,14 +4,14 @@
 using namespace amrex;
 
 std::unique_ptr<ProblemBase>
-amrex_probinit(
+amrex_probinit (
     const amrex_real* /*problo*/,
     const amrex_real* /*probhi*/)
 {
     return std::make_unique<Problem>();
 }
 
-Problem::Problem()
+Problem::Problem ()
 {
   // Parse params
   amrex::ParmParse pp("prob");
@@ -25,7 +25,7 @@ Problem::Problem()
 }
 
 void
-Problem::init_custom_pert(
+Problem::init_custom_pert (
     const Box& bx,
     const Box& xbx,
     const Box& ybx,
@@ -35,119 +35,121 @@ Problem::init_custom_pert(
     Array4<Real      > const& x_vel_pert,
     Array4<Real      > const& y_vel_pert,
     Array4<Real      > const& z_vel_pert,
-    Array4<Real      > const& r_hse,
-    Array4<Real      > const& p_hse,
+    Array4<Real      > const& /*r_hse*/,
+    Array4<Real      > const& /*p_hse*/,
     Array4<Real const> const& z_nd,
-    Array4<Real const> const& z_cc,
+    Array4<Real const> const& /*z_cc*/,
     GeometryData const& geomdata,
     Array4<Real const> const& /*mf_m*/,
     Array4<Real const> const& /*mf_u*/,
     Array4<Real const> const& /*mf_v*/,
     const SolverChoice& sc)
 {
-  const int khi = geomdata.Domain().bigEnd()[2];
+    const int khi = geomdata.Domain().bigEnd()[2];
 
     const bool use_moisture = (sc.moisture_type != MoistureType::None);
 
-  AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
+    AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
 
-  amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-    // Geometry (note we must include these here to get the data on device)
-    const auto prob_lo  = geomdata.ProbLo();
-    const auto dx       = geomdata.CellSize();
-    const Real x = prob_lo[0] + (i + 0.5) * dx[0];
-    const Real z = z_cc(i,j,k);
+    amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        // Set scalar = 0 everywhere
+        state_pert(i, j, k, RhoScalar_comp) = 0.0;
 
-    // Set scalar = 0 everywhere
-    state_pert(i, j, k, RhoScalar_comp) = 0.0;
+        if (use_moisture) {
+            state_pert(i, j, k, RhoQ1_comp) = 0.0;
+            state_pert(i, j, k, RhoQ2_comp) = 0.0;
+        }
+    });
 
-    if (use_moisture) {
-        state_pert(i, j, k, RhoQ1_comp) = 0.0;
-        state_pert(i, j, k, RhoQ2_comp) = 0.0;
-    }
-  });
-
-  // Set the x-velocity
-  amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      x_vel_pert(i, j, k) = parms.U_0;
-  });
-
-  // Set the y-velocity
-  amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      y_vel_pert(i, j, k) = 0.0;
-  });
-
-  const auto dx = geomdata.CellSize();
-  amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
-  dxInv[0] = 1. / dx[0];
-  dxInv[1] = 1. / dx[1];
-  dxInv[2] = 1. / dx[2];
-
-  // Set the z-velocity from impenetrable condition
-  amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-  {
-      z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv);
-  });
-
-  amrex::Gpu::streamSynchronize();
+    // Set the x-velocity
+    amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        x_vel_pert(i, j, k) = parms.U_0;
+    });
+
+    // Set the y-velocity
+    amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        y_vel_pert(i, j, k) = 0.0;
+    });
+
+    const auto dx = geomdata.CellSize();
+    amrex::GpuArray<Real, AMREX_SPACEDIM> dxInv;
+    dxInv[0] = 1. / dx[0];
+    dxInv[1] = 1. / dx[1];
+    dxInv[2] = 1. / dx[2];
+
+    // Set the z-velocity from impenetrable condition
+    amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        z_vel_pert(i, j, k) = WFromOmega(i, j, k, 0.0, x_vel_pert, y_vel_pert, z_nd, dxInv);
+    });
 
+    amrex::Gpu::streamSynchronize();
 }
 
 void
-Problem::init_custom_terrain(
+Problem::init_custom_terrain (
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real num  = 8 * a * a * a;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-            // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
-
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-            // Real y = (jj  * dx[1] - ycen);
-
-            // WoA Hill in x-direction
-            Real height = num / (x*x + 4 * a * a);
 
-            // Populate terrain height
-            z_arr(i,j,k0) = height;
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real num  = 8 * a * a * a;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        // Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+                // int jj = amrex::min(amrex::max(j,domlo_y),domhi_y);
+
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+                // Real y = (jj  * dx[1] - ycen);
+
+                // WoA Hill in x-direction
+                Real height = num / (x*x + 4 * a * a);
+
+                // Populate terrain height
+                z_arr(i,j,k0) = height;
+            });
+        }
     }
 }
diff --git a/Exec/SpongeTest/prob.cpp b/Exec/SpongeTest/prob.cpp
index 469ba545a..2cb3b36e2 100644
--- a/Exec/SpongeTest/prob.cpp
+++ b/Exec/SpongeTest/prob.cpp
@@ -87,54 +87,64 @@ void
 Problem::init_custom_terrain(
     const Geometry& geom,
     MultiFab& z_phys_nd,
-    const Real& /*time*/)
+    const Real& time)
 {
-    // Domain cell size and real bounds
-    auto dx = geom.CellSizeArray();
-    auto ProbLoArr = geom.ProbLoArray();
-    auto ProbHiArr = geom.ProbHiArray();
-
-    // Domain valid box (z_nd is nodal)
-    const amrex::Box& domain = geom.Domain();
-    int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
-    // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
-    int domlo_z = domain.smallEnd(2);
-
-    // User function parameters
-    Real a    = 0.5;
-    Real num  = 8 * a * a * a;
-    Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
-    Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
-
-    // Number of ghost cells
-    int ngrow = z_phys_nd.nGrow();
-
-    // Populate bottom plane
-    int k0 = domlo_z;
-
-    for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
-    {
-        // Grown box with no z range
-        amrex::Box xybx = mfi.growntilebox(ngrow);
-        xybx.setRange(2,0);
-
-        amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
-
-        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
-
-            // Clip indices for ghost-cells
-            int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
-
-            // Location of nodes
-            Real x = (ii  * dx[0] - xcen);
-
-            if(fabs(x)<a){
-                z_arr(i,j,k0) = pow(a*a - x*x, 0.5);
-            }
-            else{
-                z_arr(i,j,k0) = 0.0;
-            }
-
-        });
+    // Check if a valid csv file exists for the terrain
+    std::string fname;
+    amrex::ParmParse pp("erf");
+    auto valid_fname = pp.query("terrain_file_name",fname);
+    if (valid_fname) {
+        this->read_custom_terrain(fname,geom,z_phys_nd,time);
+    } else {
+
+        // Domain cell size and real bounds
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+
+        // Domain valid box (z_nd is nodal)
+        const amrex::Box& domain = geom.Domain();
+        int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
+        // int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
+        int domlo_z = domain.smallEnd(2);
+
+        // User function parameters
+        Real a    = 0.5;
+        Real num  = 8 * a * a * a;
+        Real xcen = 0.5 * (ProbLoArr[0] + ProbHiArr[0]);
+        Real ycen = 0.5 * (ProbLoArr[1] + ProbHiArr[1]);
+
+        // Number of ghost cells
+        int ngrow = z_phys_nd.nGrow();
+
+        // Populate bottom plane
+        int k0 = domlo_z;
+
+        for ( amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<Real> const& z_arr = z_phys_nd.array(mfi);
+
+            ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int)
+            {
+
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,domlo_x),domhi_x);
+
+                // Location of nodes
+                Real x = (ii  * dx[0] - xcen);
+
+                if(fabs(x)<a){
+                    z_arr(i,j,k0) = pow(a*a - x*x, 0.5);
+                }
+                else{
+                    z_arr(i,j,k0) = 0.0;
+                }
+
+            });
+        }
     }
 }
diff --git a/Source/prob_common.H b/Source/prob_common.H
index acba961b7..56190952e 100644
--- a/Source/prob_common.H
+++ b/Source/prob_common.H
@@ -323,6 +323,86 @@ public:
         });
     }
 
+    /**
+     * Function to perform custom initialization of terrain from an input file
+     *
+     * Note: Terrain functionality can also be used to provide grid stretching.
+     *
+     * @param[in] geom container for geometric information
+     * @param[out] z_phys_nd height coordinate at nodes
+     * @param[in] time current time
+    */
+    virtual void
+    read_custom_terrain (const std::string& fname,
+                         const amrex::Geometry& geom,
+                         amrex::MultiFab& z_phys_nd,
+                         const amrex::Real& /*time*/)
+    {
+        // Read terrain file on the host
+        amrex::Print()<<"Reading terrain file: "<< fname << std::endl;
+        std::ifstream file(fname);
+        amrex::Gpu::HostVector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
+        amrex::Real value1,value2,value3;
+        while(file>>value1>>value2>>value3){
+            m_xterrain.push_back(value1);
+            m_yterrain.push_back(value2);
+            m_zterrain.push_back(value3);
+        }
+        file.close();
+
+        // Copy data to the GPU
+        int ncell = m_xterrain.size();
+        amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(ncell),d_yterrain(ncell),d_zterrain(ncell);
+        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
+        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
+        amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
+        amrex::Real* d_xt = d_xterrain.data();
+        amrex::Real* d_yt = d_yterrain.data();
+        amrex::Real* d_zt = d_zterrain.data();
+
+        // Populate z_phys data
+        amrex::Real tol = 1.0e-4;
+        int ngrow = z_phys_nd.nGrow();
+        auto dx = geom.CellSizeArray();
+        auto ProbLoArr = geom.ProbLoArray();
+        int ilo = geom.Domain().smallEnd(0);
+        int jlo = geom.Domain().smallEnd(1);
+        int klo = geom.Domain().smallEnd(2);
+        int ihi = geom.Domain().bigEnd(0);
+        int jhi = geom.Domain().bigEnd(1);
+        for (amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            // Grown box with no z range
+            amrex::Box xybx = mfi.growntilebox(ngrow);
+            xybx.setRange(2,0);
+
+            amrex::Array4<amrex::Real> const& z_arr = z_phys_nd.array(mfi);
+            amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
+            {
+                // Clip indices for ghost-cells
+                int ii = amrex::min(amrex::max(i,ilo),ihi);
+                int jj = amrex::min(amrex::max(j,jlo),jhi);
+
+                // Location of nodes
+                bool found = false;
+                amrex::Real x = ProbLoArr[0]  + ii  * dx[0];
+                amrex::Real y = ProbLoArr[1]  + jj  * dx[1];
+                amrex::Real zloc = 0.0;
+                for(int n=0; n<ncell; ++n)
+                {
+                    amrex::Real delta=std::sqrt(std::pow(x-d_xt[n],2)+std::pow(y-d_yt[n],2));
+                    if(delta<tol){
+                        found=true;
+                        zloc=d_zt[n];
+                        continue;
+                    }
+                }
+                AMREX_ASSERT_WITH_MESSAGE(found, "Location read from terrain file does not match the grid!");
+                z_arr(i,j,klo) = zloc;
+            });
+        }
+    }
+
 #ifdef ERF_USE_TERRAIN_VELOCITY
     virtual amrex::Real compute_terrain_velocity(const amrex::Real /*time*/)
     {