diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp
index 5aed34a33..436da69ce 100644
--- a/Exec/RegTests/GABLS1/prob.cpp
+++ b/Exec/RegTests/GABLS1/prob.cpp
@@ -55,7 +55,7 @@ Problem::init_custom_pert(
     amrex::Array4<amrex::Real      > const& r_hse,
     amrex::Array4<amrex::Real      > const& /*p_hse*/,
     amrex::Array4<amrex::Real const> const& /*z_nd*/,
-    amrex::Array4<amrex::Real const> const& /*z_cc*/,
+    amrex::Array4<amrex::Real const> const& z_cc,
     amrex::GeometryData const& geomdata,
     amrex::Array4<amrex::Real const> const& /*mf_m*/,
     amrex::Array4<amrex::Real const> const& /*mf_u*/,
@@ -71,21 +71,35 @@ Problem::init_custom_pert(
     amrex::Print() << "Initializing QKE" << std::endl;
   }
 
+  const bool use_terrain = sc.use_terrain;
+
   if (parms.pert_ref_height > 0) {
     if (parms.pert_deltaU != 0.0) {
       amrex::Print() << "Adding divergence-free x-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("u perturbations not supported for terrain");
+      }
     }
     if (parms.pert_deltaV != 0.0) {
       amrex::Print() << "Adding divergence-free y-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("v perturbations not supported for terrain");
+      }
     }
     if (parms.T_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random theta perturbations" << std::endl;
     }
     if (parms.U_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random x-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("u perturbations not supported for terrain");
+      }
     }
     if (parms.V_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random y-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("v perturbations not supported for terrain");
+      }
     }
   }
 
@@ -96,7 +110,7 @@ Problem::init_custom_pert(
     const Real* dx = geomdata.CellSize();
     const Real x = prob_lo[0] + (i + 0.5) * dx[0];
     const Real y = prob_lo[1] + (j + 0.5) * dx[1];
-    const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+    const Real z = use_terrain ? z_cc(i,j,k) : prob_lo[2] + (k + 0.5) * dx[2];
 
     // Define a point (xc,yc,zc) at the center of the domain
     const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]);
@@ -196,4 +210,3 @@ Problem::init_custom_pert(
     }
   });
 }
-
diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp
index 53b7d3a21..81740311c 100644
--- a/Source/Diffusion/PBLModels.cpp
+++ b/Source/Diffusion/PBLModels.cpp
@@ -168,7 +168,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
                 // First Length Scale
                 AMREX_ASSERT(l_obukhov != 0);
                 int lk = amrex::max(k,0);
-                const Real zval = gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
+                const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,lk,z_nd_arr) : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
                 const Real zeta = zval/l_obukhov;
                 Real l_S;
                 if (zeta >= 1.0) {
diff --git a/Source/ERF.H b/Source/ERF.H
index 4d2c65827..d3be2249a 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -656,7 +656,9 @@ private:
                                                     // eddy viscosity
                                                     "Kmv","Kmh",
                                                     // eddy diffusivity of heat
-                                                    "Khv","Khh"
+                                                    "Khv","Khh",
+                                                    // mynn pbl lengthscale
+                                                    "Lpbl"
                                                      ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup"
 #ifdef ERF_COMPUTE_ERROR
                                                     ,"xvel_err", "yvel_err", "zvel_err", "pp_err"
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
index a688fa5d2..1f75e5742 100644
--- a/Source/IO/Plotfile.cpp
+++ b/Source/IO/Plotfile.cpp
@@ -631,6 +631,10 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Theta_h,mf_comp,1,0);
             mf_comp ++;
         }
+        if (containerHasElement(plot_var_names, "Lpbl")) {
+            MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::PBL_lengthscale,mf_comp,1,0);
+            mf_comp ++;
+        }
 
         // NOTE: Protect against accessing non-existent data
         if (use_moisture) {