diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
index 9e749874c3..332321b209 100644
--- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90
+++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
@@ -34,6 +34,9 @@ module MOM_entrain_diffusive
                              !! calculate the diapycnal entrainment.
   real    :: Tolerance_Ent   !< The tolerance with which to solve for entrainment values
                              !! [H ~> m or kg m-2].
+  real    :: max_Ent         !< A large ceiling on the maximum permitted amount of entrainment
+                             !! across each interface between the mixed and buffer layers within
+                             !! a timestep [H ~> m or kg m-2].
   real    :: Rho_sig_off     !< The offset between potential density and a sigma value [R ~> kg m-3]
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                              !! regulate the timing of diagnostic output.
@@ -1053,7 +1056,6 @@ subroutine set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl,
   real, dimension(SZI_(G), SZK_(GV)) :: &
     S_est ! An estimate of the coordinate potential density - 1000 after
           ! entrainment for each layer [R ~> kg m-3].
-  real :: max_ent  ! The maximum possible entrainment [H ~> m or kg m-2].
   real :: dh       ! An available thickness [H ~> m or kg m-2].
   real :: Kd_x_dt  ! The diffusion that remains after thin layers are
                    ! entrained [H2 ~> m2 or kg2 m-4].
@@ -1063,8 +1065,6 @@ subroutine set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl,
   integer :: i, k, is, ie, nz
   is = G%isc ; ie = G%iec ; nz = GV%ke
 
-!  max_ent = 1.0e14*GV%Angstrom_H ! This is set to avoid roundoff problems.
-  max_ent = 1.0e4*GV%m_to_H
   h_neglect = GV%H_subroundoff
 
   do i=is,ie ; pres(i) = tv%P_Ref ; enddo
@@ -1084,8 +1084,7 @@ subroutine set_Ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl,
 
   do k=2,kmb ; do i=is,ie
     if (do_i(i)) then
-      Ent_bl(i,K) = min(2.0 * dtKd_int(i,K) / (h(i,j,k-1) + h(i,j,k) + h_neglect), &
-                        max_ent)
+      Ent_bl(i,K) = min(2.0 * dtKd_int(i,K) / (h(i,j,k-1) + h(i,j,k) + h_neglect), CS%max_Ent)
     else ; Ent_bl(i,K) = 0.0 ; endif
   enddo ; enddo
 
@@ -2116,6 +2115,10 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re
                  "The tolerance with which to solve for entrainment values.", &
                  units="m", default=US%Z_to_m*MAX(100.0*GV%Angstrom_Z,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H, &
                  do_not_log=just_read_params)
+  call get_param(param_file, mdl, "ENTRAIN_DIFFUSIVE_MAX_ENT", CS%max_Ent, &
+                 "A large ceiling on the maximum permitted amount of entrainment across each "//&
+                 "interface between the mixed and buffer layers within a timestep.", &
+                 units="m", default=1.0e4, scale=GV%m_to_H, do_not_log=.not.CS%bulkmixedlayer)
 
   CS%Rho_sig_off = 1000.0*US%kg_m3_to_R