diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H
index 272213257..b15bd2410 100644
--- a/Source/BoundaryConditions/MOSTStress.H
+++ b/Source/BoundaryConditions/MOSTStress.H
@@ -256,8 +256,8 @@ struct adiabatic_wave_coupled
         u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k));
         do {
             ustar = u_star_arr(i,j,k);
-            z0    = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
-                             + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps);
+            z0    = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
+                                      + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
             u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0);
             ++iter;
         } while ((std::abs(u_star_arr(i,j,k) - ustar) > tol) && iter <= max_iters);
@@ -273,6 +273,7 @@ private:
     const amrex::Real tol = 1.0e-5;
     const amrex::Real eps = 1e-15;
     const amrex::Real z0_eps = 1.0e-6;
+    const amrex::Real z0_max = 0.1;
 };
 
 
@@ -515,8 +516,8 @@ struct surface_flux_wave_coupled
         u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k));
         do {
             ustar = u_star_arr(i,j,k);
-            z0    = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
-                             + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps);
+            z0    = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
+                                      + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
             Olen = -ustar * ustar * ustar * tm_arr(i,j,k) /
                    (mdata.kappa * mdata.gravity * mdata.surf_temp_flux);
             zeta  = mdata.zref / Olen;
@@ -539,6 +540,7 @@ private:
     const amrex::Real tol = 1.0e-5;
     const amrex::Real eps = 1e-15;
     const amrex::Real z0_eps = 1.0e-6;
+    const amrex::Real z0_max = 0.1;
 };
 
 
@@ -788,8 +790,8 @@ struct surface_temp_wave_coupled
         u_star_arr(i,j,k) = mdata.kappa * umm_arr(i,j,k) / std::log(mdata.zref / z0_arr(i,j,k));
         do {
             ustar = u_star_arr(i,j,k);
-            z0    = std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
-                             + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps);
+            z0    = std::min( std::max(1200.0 * Hwave_arr(i,j,k) * std::pow( Hwave_arr(i,j,k)/(Lwave_arr(i,j,k)+eps), 4.5 )
+                                      + 0.11 * eta_arr(ie,je,k,EddyDiff::Mom_v) / ustar, z0_eps), z0_max );
             tflux = -(tm_arr(i,j,k) - t_surf_arr(i,j,k)) * ustar * mdata.kappa /
                      (std::log(mdata.zref / z0) - psi_h);
             Olen = -ustar * ustar * ustar * tm_arr(i,j,k) /
@@ -813,6 +815,7 @@ private:
     const amrex::Real tol = 1.0e-5;
     const amrex::Real eps = 1e-15;
     const amrex::Real z0_eps = 1.0e-6;
+    const amrex::Real z0_max = 0.1;
 };