diff --git a/Docs/sphinx_doc/Verification.rst b/Docs/sphinx_doc/Verification.rst
index 589e86ebd..43bd94dd1 100644
--- a/Docs/sphinx_doc/Verification.rst
+++ b/Docs/sphinx_doc/Verification.rst
@@ -137,3 +137,4 @@ with the exact solution at two different horizontal locations are shown below.
 
 .. image:: figures/Terrain3d_Hemisphere.png
   :width: 600
+
diff --git a/Exec/RegTests/DynamicRefinement/inputs b/Exec/RegTests/DynamicRefinement/inputs
index ead24d4e8..06ed35c52 100644
--- a/Exec/RegTests/DynamicRefinement/inputs
+++ b/Exec/RegTests/DynamicRefinement/inputs
@@ -80,6 +80,6 @@ erf.hi_scal2.value_greater = 10.
 
 amr.n_error_buf  = 4
 erf.regrid_int   = 2 
-erf.cf_width     = 4
-erf.cf_set_width = 1
+erf.cf_width     = 0 # Internal relaxation
+erf.cf_set_width = 0 # Internal Dirichlet
 
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
index 282868c56..66864fb37 100644
--- a/Source/IO/Plotfile.cpp
+++ b/Source/IO/Plotfile.cpp
@@ -759,6 +759,20 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
 #endif
     }
 
+    // Fill terrain distortion MF
+    if (solverChoice.use_terrain) {
+        for (int lev(0); lev <= finest_level; ++lev) {
+            MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0);
+            Real dz = Geom()[lev].CellSizeArray()[2];
+            for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+                const Box& bx = mfi.tilebox();
+                Array4<      Real> mf_arr = mf_nd[lev].array(mfi);
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    mf_arr(i,j,k,2) -= k * dz;
+                });
+            }
+        }
+    }
 
     std::string plotfilename;
     if (which == 1)
@@ -771,17 +785,6 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         if (plotfile_type == "amrex") {
             amrex::Print() << "Writing plotfile " << plotfilename << "\n";
             if (solverChoice.use_terrain) {
-                // We started with mf_nd holding 0 in every component; here we fill only the offset in z
-                int lev = 0;
-                MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0);
-                Real dz = Geom()[lev].CellSizeArray()[2];
-                for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-                    const Box& bx = mfi.tilebox();
-                    Array4<      Real> mf_arr = mf_nd[lev].array(mfi);
-                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                        mf_arr(i,j,k,2) -= k * dz;
-                    });
-                }
                 WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
                                                    GetVecOfConstPtrs(mf),
                                                    GetVecOfConstPtrs(mf_nd),
@@ -789,9 +792,9 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                                                    t_new[0], istep);
             } else {
                 WriteMultiLevelPlotfile(plotfilename, finest_level+1,
-                                               GetVecOfConstPtrs(mf),
-                                               varnames,
-                                               Geom(), t_new[0], istep, refRatio());
+                                        GetVecOfConstPtrs(mf),
+                                        varnames,
+                                        Geom(), t_new[0], istep, refRatio());
             }
             writeJobInfo(plotfilename);
 #ifdef ERF_USE_HDF5
@@ -864,8 +867,19 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
             }
 
-            WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames,
-                                           g2, t_new[0], istep, rr);
+            amrex::Print() << "Writing plotfile " << plotfilename << "\n";
+            if (solverChoice.use_terrain) {
+                WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
+                                                   GetVecOfConstPtrs(mf),
+                                                   GetVecOfConstPtrs(mf_nd),
+                                                   varnames,
+                                                   t_new[0], istep);
+            } else {
+                WriteMultiLevelPlotfile(plotfilename, finest_level+1,
+                                        GetVecOfConstPtrs(mf2), varnames,
+                                        g2, t_new[0], istep, rr);
+            }
+
             writeJobInfo(plotfilename);
 #ifdef ERF_USE_NETCDF
         } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") {
@@ -958,7 +972,7 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl
             } else {
                 data = mf[level];
             }
-            VisMF::Write(*data       , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
+            VisMF::Write(*data        , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
             VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix));
         }
     }