diff --git a/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag
index 83175b8a8..f31d16159 100644
--- a/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag
+++ b/Exec/RegTests/GABLS1/inputs_gabls1_mynn25_smag
@@ -15,9 +15,8 @@ geometry.is_periodic = 1 1 0
 zlo.type                    = "Most"
 erf.most.z0                 = 0.1  # from Cuxart et al. 2006
 erf.most.zref               = 3.125  # == dz/2
-erf.most.surf_temp_flux     = 0.0
-#erf.most.surf_temp          = 265.0 # initial value, should match input_sounding
-#erf.most.surf_heating_rate  = -0.25 # [K/h] from Cuxart et al. 2006
+erf.most.surf_temp          = 265.0 # initial value, should match input_sounding
+erf.most.surf_heating_rate  = -0.25 # [K/h] from Cuxart et al. 2006
 
 zhi.type        = "SlipWall"
 zhi.theta_grad  = 0.01  # [K/m] to match the input sounding
@@ -28,8 +27,10 @@ erf.init_sounding_ideal = 1
 erf.input_sounding_file = "input_sounding_GABLS1"
 
 # TIME STEP CONTROL
-erf.fixed_dt        = 10.0  # fixed time step (Cuxart et al. 2006)
-erf.no_substepping  = 1
+#erf.fixed_dt        = 10.0  # fixed time step (Cuxart et al. 2006)
+#erf.no_substepping  = 1
+erf.fixed_dt        = 1.0  # largest stable low Mach dt
+erf.fixed_mri_dt_ratio = 6
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval    = 1       # timesteps between computing mass
@@ -41,36 +42,36 @@ amr.max_level       = 0       # maximum level number allowed
 
 # CHECKPOINT FILES
 erf.check_file      = chk        # root name of checkpoint file
-erf.check_int       = 100        # number of timesteps between checkpoints
+erf.check_int       = 600        # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1     = plt       # prefix of plotfile name
-erf.plot_int_1      = 1        # number of timesteps between plotfiles
-erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure temp theta
+erf.plot_int_1      = 600        # number of timesteps between plotfiles
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv
 
 
 # SOLVER CHOICE
-erf.alpha_T = 0.0
-erf.alpha_C = 1.0
+erf.dycore_vert_adv_type = "Upwind_3rd"
+
+erf.molec_diff_type = "None"
+#erf.alpha_T = 0.0
+#erf.alpha_C = 1.0
+
 erf.use_gravity = true
 
-# Geostrophic wind (Cuxart et al. 2006)
-# f = 1.39e-4 s^-1
+# Coriolis parameter f = 1.39e-4 s^-1 (Cuxart et al. 2006)
 erf.use_coriolis = true
 erf.latitude = 73.0
 erf.rotational_time_period = 86455.2516813368
 
-#erf.abl_driver_type = "GeostrophicWind"
-#erf.abl_geo_wind = 8.0 0.0 0.0
-# manually calculate:
-erf.abl_geo_forcing = 0.0 0.001112 0.0  # neglect vertical component
+# Geostrophic wind (Cuxart et al. 2006)
+erf.abl_driver_type = "GeostrophicWind"
+erf.abl_geo_wind = 8.0 0.0 0.0
 
 # Turbulence closure
-erf.molec_diff_type = "None"
 erf.les_type        = "Smagorinsky"
 erf.Cs              = 0.1
 erf.rho0_trans      = 1.3223 # from Cuxart et al. 2006
 erf.theta_ref       = 263.5 # from Cuxart et al. 2006
 
 erf.pbl_type    = "MYNN2.5"
-#erf.QKE_0       = 0.8  # == 2*KE_0; note: Cuxart et al. 2006 has an initial TKE profile
diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp
index 25d70f4f0..c9bffb30f 100644
--- a/Exec/RegTests/GABLS1/prob.cpp
+++ b/Exec/RegTests/GABLS1/prob.cpp
@@ -97,6 +97,7 @@ Problem::init_custom_pert(
 
     // Set an initial value for QKE = 2*e
     if (z < 250) {
+        // Following Cuxart et al. 2006
         state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * 0.8 * std::pow(1 - z/250, 3);
     } else {
         state(i, j, k, RhoQKE_comp) = r_hse(i,j,k) * parms.QKE_0;
diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp
index 1e7ac7b6a..c890a59d4 100644
--- a/Source/Diffusion/PBLModels.cpp
+++ b/Source/Diffusion/PBLModels.cpp
@@ -227,7 +227,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
 
           amrex::Real SM = (R2*E2 - R1*E4)/(E2*E3 - E1*E4);
           amrex::Real SH = (R1*E3 - R2*E1)/(E2*E3 - E1*E4);
-          amrex::Real SQ = 3.0 * SM;
+          amrex::Real SQ = 3.0 * SM; // Nakanishi & Niino 2009
 
           // Finally, compute the eddy viscosity/diffusivities
           const amrex::Real rho = cell_data(i,j,k,Rho_comp);