diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index 74fc7d8a3..b0dd49652 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -30,6 +30,14 @@ New Features
 
 MESA no longer stops when reactions for which special rates are set are not in the nuclear network, only a warning is printed. This is intended to make it easier to test various network sizes without having to also change the list of special reactions.
 
+Two new profile columns, ``superad_reduction_Lrad_div_Ledd`` and ``superad_reduction_trigger``, and one new history column, ``num_cells_with_superad_reduction``, expose what the ``use_superad_reduction`` throttle is doing per cell (which threshold fired and at what proximity to Eddington). The existing ``superad_reduction_factor`` profile column and ``max_superad_reduction_factor`` history column are unchanged. The new ``superad_reduction_trigger`` is an integer flag: 1 = Eddington proximity (``Lrad/Ledd > Gamma_limit``); 2 = density-inversion criterion (Joss et al. 1973, Paxton, Cantiello et al. 2013 eq. 17); 3 = both.
+
+Opt-in convective-turnover-time limiter for ``use_superad_reduction``: when the new control ``superad_reduction_use_turnover_limit`` is set to ``.true.``, the throttle excess ``(Gamma_factor - 1)`` is multiplied by ``f_tau = 1 - exp(-dt / tau_conv)``, where ``tau_conv = scale_height(k) / mlt_vc_old(k)`` is the local convective turnover time built from the previous-step converged convective velocity. The limiter recovers the current behavior when ``dt >> tau_conv`` and smoothly suppresses the throttle when ``dt << tau_conv`` (where the implicit-step throttle would be acting faster than convection itself can respond). Default ``.false.``; bit-for-bit unchanged on every shipped test_suite case.
+
+Companion control to the turnover-time limiter: ``superad_reduction_turnover_vc_floor_frac`` (real, default ``0d0``). When set to a small positive value, floors ``mlt_vc_old(k)`` at ``frac * csound_face(k)`` before computing ``tau_conv``. This prevents ``f_tau`` from collapsing to ~0 in slow-convection iron-bump cells where the bare limiter would fully suppress the throttle and re-expose the radiation-pressure inversions to the Newton solver. With ``frac = 1d-3`` the floor restores baseline-quality numerics through the post-He burning crunch (60 M_sun benchmark: 15 retries through fe_core_infall, matching the limiter-off baseline's 14 retries; the bare limiter accumulates ~780 retries in the same window). Default ``0d0`` disables the floor and recovers the bare-limiter behavior exactly.
+
+The turnover-time limiter formula now anchors the relaxation at the previous step's converged ``Gamma_factor`` rather than at 1 (no throttle): ``Gamma_factor_new = Gamma_factor_old + f_tau * (Gamma_factor_inst - Gamma_factor_old)`` instead of the earlier ``Gamma_factor_new = 1 + f_tau * (Gamma_factor_inst - 1)``. The two limits remain physically correct (``dt >> tau_conv`` recovers the instantaneous value, ``dt << tau_conv`` holds the throttle at its previous value) but the new formula no longer asks convection to *release* the throttle on timescales shorter than ``tau_conv`` -- which the old formula did, and which caused a small but visible HRD jump at the end of the to_cc phase when ``dt`` collapsed for iron-core formation. The change required a new paired ``superad_reduction_factor_old(:)`` array (mirroring the ``mlt_vc/mlt_vc_old`` pattern). On the very first step ``Gamma_factor_old`` is unset and the formula falls through to the old behavior, preserving first-step output bit-for-bit. No new control: the new formula is strictly more correct on every step after the first and is enabled whenever ``superad_reduction_use_turnover_limit = .true.``.
+
 .. _Bug Fixes main:
 
 Bug Fixes
@@ -45,6 +53,10 @@ The parameter ``report_max_infall_inside_fe_core`` was ignored in versions r25.1
 
 ``fe_core_infall_limit`` now obeys ``when_to_stop_rtol`` and ``when_to_stop_atol`` again (broken since r11532).
 
+In ``set_superad_reduction`` (``star/private/turb_support.f90``) the density-inversion contribution to ``Gamma_term`` was being scaled by ``superad_reduction_Gamma_limit_scale`` instead of the intended ``superad_reduction_Gamma_inv_scale``, making the latter control silently inert. Restored to the as-published algorithm (Jermyn et al. 2023 eq. 64, ``alpha_2`` term). Existing test_suite cases that use ``superad_reduction`` all set ``Gamma_limit_scale = Gamma_inv_scale``, so they are bit-for-bit unaffected; the fix is observable only when the two scales are set to different values.
+
+The default value of ``superad_reduction_limit`` is changed from ``-1d0`` (uncapped) to ``100d0`` (logistic saturation at ``Gamma_factor ~= 100``). The cap was already documented as the recommended setting, but the shipped default disabled it. With the cap inactive, the raw ``1/sqrt(beta)`` boost in radiation-pressure-dominated cells can drive ``Gamma_factor`` arbitrarily large in a single Newton iteration. Existing test_suite cases that exercise ``superad_reduction`` all explicitly set this control to ``-1d0``, so they are bit-for-bit unaffected.
+
 
 .. note:: Before releasing a new version of MESA, move `Changes in main` to a new section below with the version number as the title, and add a new `Changes in main` section at the top of the file (see ``changelog_template.rst``).
 
diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults
index 357fe8631..950aec433 100644
--- a/star/defaults/controls.defaults
+++ b/star/defaults/controls.defaults
@@ -2655,27 +2655,168 @@
 
       ! use_superad_reduction
       ! ~~~~~~~~~~~~~~~~~~~~~
+
+      ! Implicit alternative to ``okay_to_reduce_gradT_excess`` (MLT++).
+      ! When ``.true.``, throttles the radiative gradient ``gradr`` in
+      ! convective cells that either (a) carry a radiative luminosity
+      ! approaching the Eddington limit, or (b) sit in the
+      ! density-inversion regime (Joss--Salpeter--Ostriker). MLT/TDC is then
+      ! re-solved on the throttled ``gradr``, so the final ``gradT`` is
+      ! self-consistent with MLT for the reduced radiative input.
+      ! Acts on the *input* to MLT, not the output --- unlike MLT++.
+
+      ! ::
+
+      use_superad_reduction = .false.
+
+
       ! superad_reduction_Gamma_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Radiative-Eddington trigger threshold:
+      ! ``Lrad/Ledd = (4 a T^4 / 3 P) * gradT``.
+      ! When this ratio exceeds ``Gamma_limit`` in a convective cell,
+      ! the Eddington-proximity contribution to the throttle activates.
+
+      ! ::
+
+      superad_reduction_Gamma_limit = 0.5d0
+
+
       ! superad_reduction_Gamma_limit_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Sensitivity of the throttle to ``Lrad/Ledd`` excess above
+      ! ``Gamma_limit``. The Eddington contribution to ``Gamma_factor``
+      ! grows quadratically (then linearly) with the fractional excess;
+      ! ``Gamma_limit_scale`` is the proportionality constant. Larger
+      ! values produce stronger throttling for a given excess.
+
+      ! ::
+
+      superad_reduction_Gamma_limit_scale = 5d0
+
+
       ! superad_reduction_Gamma_inv_scale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Sensitivity of the throttle to the density-inversion criterion
+      ! ``Lrad/Ledd > 4(1-beta)/(4-3 beta)``, where ``beta = Pgas/Ptot``.
+      ! Analogous to ``Gamma_limit_scale`` but for the inversion term.
+      ! Tune independently when inversion proximity dominates over
+      ! Eddington proximity (or vice versa).
+
+      ! ::
+
+      superad_reduction_Gamma_inv_scale = 5d0
+
+
       ! superad_reduction_diff_grads_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Superadiabaticity floor below which the throttle is silenced.
+      ! When ``gradT - gradL < diff_grads_limit``, the cell is treated
+      ! as essentially adiabatic and the throttle contribution drops
+      ! smoothly to zero (quintic smoothstep). Prevents the throttle
+      ! from amplifying numerical noise in nearly-adiabatic cells.
+
+      ! ::
+
+      superad_reduction_diff_grads_limit = 1d-3
+
+
       ! superad_reduction_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! Implicit alternative to ``okay_to_reduce_gradT_excess``
+      ! Upper saturation cap on the throttle multiplier ``Gamma_factor``.
+      ! The raw ``Gamma_factor`` includes a ``1/sqrt(beta)`` boost that
+      ! can blow up in radiation-dominated cells; this control imposes
+      ! a logistic squash so ``Gamma_factor -> reduction_limit`` as the
+      ! raw value diverges. Set ``-1d0`` to disable the cap (legacy
+      ! behavior); finite values guard against the second MLT call
+      ! seeing ``gradr_scaled`` jumping all the way to ``gradL`` in a
+      ! single Newton iteration.
 
       ! ::
 
-      use_superad_reduction = .false.
-      superad_reduction_Gamma_limit = 0.5d0
-      superad_reduction_Gamma_limit_scale = 5d0
-      superad_reduction_Gamma_inv_scale = 5d0
-      superad_reduction_diff_grads_limit = 1d-3
-      superad_reduction_limit = -1d0
+      superad_reduction_limit = 100d0
+
+
+      ! superad_reduction_use_turnover_limit
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Optional convective-turnover-time limiter for the throttle.
+      ! When ``.true.``, relaxes the throttle using a fraction
+      ! ``f_tau`` set by ``superad_reduction_turnover_limit_function``,
+      ! where ``tau_conv = scale_height(k) / mlt_vc_old(k)`` is the
+      ! local convective turnover time computed from the previous
+      ! step's converged convective velocity. This reflects the
+      ! physical constraint that convection cannot redistribute the
+      ! radiative input on timescales shorter than its own turnover
+      ! time --- so when the global time step ``dt`` is much shorter
+      ! than ``tau_conv`` the throttle should not act at full strength.
+      !
+      ! - ``dt >> tau_conv``: ``f_tau -> 1``, current behavior recovered.
+      ! - ``dt << tau_conv``: ``f_tau -> 0``, throttle smoothly off.
+      !
+      ! ``mlt_vc_old`` is used (rather than the current iterate) so
+      ! that ``tau_conv`` is a frozen parameter w.r.t. the inner
+      ! Newton solve: the auto-diff Jacobian remains clean. When the
+      ! cell was radiative in the previous step
+      ! (``mlt_vc_old = 0``), ``f_tau -> 0`` and the throttle is off
+      ! for this step.
+      !
+      ! Default ``.false.``: the limiter block is skipped and the
+      ! routine behaves exactly as it does without this control.
+
+      ! ::
+
+      superad_reduction_use_turnover_limit = .false.
+
+      ! superad_reduction_turnover_limit_function
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Functional form for the turnover-time relaxation fraction. Allowed
+      ! values are:
+      !
+      ! - ``'exponential'``: ``f_tau = 1 - exp(-dt/tau_conv)``
+      ! - ``'linear'``: ``f_tau = min(dt/tau_conv, 1)``
+
+      ! ::
+
+      superad_reduction_turnover_limit_function = 'exponential'
+
+      ! superad_reduction_turnover_vc_floor_frac
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Fractional floor on mlt_vc_old, expressed as a fraction of the
+      ! local face sound speed csound_face(k), used only when
+      ! superad_reduction_use_turnover_limit = .true.
+      !
+      ! The turnover-time limiter uses
+      !   tau_conv(k) = scale_height(k) / mlt_vc_old(k)
+      ! to gauge whether convection has had time to respond to a step.
+      ! Slow-convection cells (e.g., the surface iron-bump zone) have
+      ! tiny mlt_vc_old, hence huge tau_conv, hence f_tau -> 0 and the
+      ! throttle is suppressed there even when the rest of the star is
+      ! evolving normally. That can leave radiation-pressure inversions
+      ! unsmoothed and make the Newton solve hard.
+      !
+      ! Setting this control to a small positive value floors
+      ! mlt_vc_old at frac * csound_face(k) before computing tau_conv.
+      ! Physical interpretation: even in nominally radiative cells
+      ! perturbations cannot propagate slower than ~ acoustic speed, so
+      ! tau_conv cannot meaningfully exceed scale_height / (frac * cs).
+      ! Numerically this prevents f_tau from dropping arbitrarily close
+      ! to zero in the slow-convection cells most prone to retries.
+      !
+      ! Default 0d0 disables the floor (current behavior). Reasonable
+      ! starting value: 1d-3 (cap tau_conv at ~ scale_height / (1e-3 cs)
+      ! = ~ 1000 sound-crossing times of the scale height).
+
+      ! ::
+
+      superad_reduction_turnover_vc_floor_frac = 0d0
 
 
    ! overshooting
diff --git a/star/defaults/history_columns.list b/star/defaults/history_columns.list
index 858f88e3d..06fad8c28 100644
--- a/star/defaults/history_columns.list
+++ b/star/defaults/history_columns.list
@@ -935,6 +935,10 @@
       !max_L_rad_div_Ledd
       !max_L_rad_div_Ledd_div_phi_Joss
 
+   !## superad_reduction
+      !max_superad_reduction_factor ! peak Gamma_factor over the star (1 = no reduction)
+      !num_cells_with_superad_reduction ! count of cells where Gamma_factor > 1
+
 
    !## RTI
       !rti_regions <num>
diff --git a/star/defaults/profile_columns.list b/star/defaults/profile_columns.list
index 2bf7338c2..4b683c4c8 100644
--- a/star/defaults/profile_columns.list
+++ b/star/defaults/profile_columns.list
@@ -511,6 +511,8 @@
    !dvc_dt_TDC_div_g ! dimensionless ratio of convective velocity to g
 
    !superad_reduction_factor ! gamma_factor from superad_reduction
+   !superad_reduction_Lrad_div_Ledd ! local Lrad/Ledd used by superad_reduction (0 if inactive)
+   !superad_reduction_trigger ! 0=untouched, 1=Eddington, 2=density-inversion, 3=both
    !conv_vel_div_mlt_vc ! conv_vel from any convection model (including RSP) / mlt_vc, for comparison 
 
    !log_Lconv
diff --git a/star/private/adjust_mesh.f90 b/star/private/adjust_mesh.f90
index a243f2442..1d46dbefd 100644
--- a/star/private/adjust_mesh.f90
+++ b/star/private/adjust_mesh.f90
@@ -410,6 +410,7 @@ integer function remesh(s)
             prv% j_rot, prv% i_rot, prv% omega, prv% D_omega, &
             prv% mlt_vc, prv% lnT, prv% w, specific_PE, specific_KE, &
             prv% m, prv% r, prv% rho, prv% dPdr_dRhodr_info, prv% D_mix, &
+            prv% superad_reduction_factor, &
             cell_type, comes_from, prv% dq, xq_old, s% xh, s% xa, s% dq, xq_new, ierr)
          if (ierr /= 0) then
             s% retry_message = 'do_mesh_adjust failed in mesh_adjust'
@@ -429,6 +430,7 @@ integer function remesh(s)
             s% prev_mesh_omega(k) = prv% prev_mesh_omega(k)
             s% prev_mesh_dq(k) = prv% prev_mesh_dq(k)
             s% prev_mesh_mlt_vc(k) = prv% prev_mesh_mlt_vc(k)
+            s% prev_mesh_superad_reduction_factor(k) = prv% prev_mesh_superad_reduction_factor(k)
          end do
 
          ! restore ST info (for time smoothing)
diff --git a/star/private/alloc.f90 b/star/private/alloc.f90
index 675345b21..9c3c58e8c 100644
--- a/star/private/alloc.f90
+++ b/star/private/alloc.f90
@@ -383,6 +383,8 @@ subroutine star_info_old_arrays(s, action, ierr)
          if (failed('dq_old')) return
          call do1D(s, s% mlt_vc_old, nz, action, ierr)
          if (failed('mlt_vc_old')) return
+         call do1D(s, s% superad_reduction_factor_old, nz, action, ierr)
+         if (failed('superad_reduction_factor_old')) return
          call do1D(s, s% omega_old, nz, action, ierr)
          if (failed('omega_old')) return
          call do1D(s, s% j_rot_old, nz, action, ierr)
@@ -776,6 +778,10 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('gradT_excess_effect')) exit
             call do1(s% superad_reduction_factor, c% superad_reduction_factor)
             if (failed('superad_reduction_factor')) exit
+            call do1(s% superad_reduction_Lrad_div_Ledd, c% superad_reduction_Lrad_div_Ledd)
+            if (failed('superad_reduction_Lrad_div_Ledd')) exit
+            call do1_integer(s% superad_reduction_trigger, c% superad_reduction_trigger)
+            if (failed('superad_reduction_trigger')) exit
 
             call do1(s% domega_dlnR, c% domega_dlnR)
             if (failed('domega_dlnR')) exit
@@ -1397,6 +1403,8 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('prev_mesh_omega')) exit
             call do1(s% prev_mesh_mlt_vc, c% prev_mesh_mlt_vc)
             if (failed('prev_mesh_mlt_vc')) exit
+            call do1(s% prev_mesh_superad_reduction_factor, c% prev_mesh_superad_reduction_factor)
+            if (failed('prev_mesh_superad_reduction_factor')) exit
             call do1(s% prev_mesh_dq, c% prev_mesh_dq)
             if (failed('prev_mesh_dq')) exit
             ! These are needed for time-smoothing of ST mixing
diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90
index 3973e3970..507a4bd0e 100644
--- a/star/private/ctrls_io.f90
+++ b/star/private/ctrls_io.f90
@@ -136,6 +136,8 @@ module ctrls_io
    gradT_excess_min_center_he4, gradT_excess_max_logT, gradT_excess_min_log_tau_full_on, gradT_excess_max_log_tau_full_off, &
     use_superad_reduction, superad_reduction_gamma_limit, superad_reduction_gamma_limit_scale, D_mix_zero_region_top_q, &
     superad_reduction_gamma_inv_scale, superad_reduction_diff_grads_limit, superad_reduction_limit, &
+    superad_reduction_use_turnover_limit, superad_reduction_turnover_limit_function, &
+    superad_reduction_turnover_vc_floor_frac, &
     make_gradr_sticky_in_solver_iters, min_logT_for_make_gradr_sticky_in_solver_iters, &
     max_logT_for_mlt, thermohaline_coeff, thermohaline_option, mixing_length_alpha, remove_small_D_limit, &
     alt_scale_height_flag, Henyey_MLT_y_param, Henyey_MLT_nu_param, no_MLT_below_shock, mlt_make_surface_no_mixing, &
@@ -1078,6 +1080,9 @@ subroutine store_controls(s, ierr)
  s% superad_reduction_gamma_inv_scale = superad_reduction_gamma_inv_scale
  s% superad_reduction_diff_grads_limit = superad_reduction_diff_grads_limit
  s% superad_reduction_limit = superad_reduction_limit
+ s% superad_reduction_use_turnover_limit = superad_reduction_use_turnover_limit
+ s% superad_reduction_turnover_limit_function = superad_reduction_turnover_limit_function
+ s% superad_reduction_turnover_vc_floor_frac = superad_reduction_turnover_vc_floor_frac
 
  s% max_logT_for_mlt = max_logT_for_mlt
  s% mlt_make_surface_no_mixing = mlt_make_surface_no_mixing
@@ -2801,6 +2806,9 @@ subroutine set_controls_for_writing(s, ierr)
  superad_reduction_gamma_inv_scale = s% superad_reduction_gamma_inv_scale
  superad_reduction_diff_grads_limit = s% superad_reduction_diff_grads_limit
  superad_reduction_limit = s% superad_reduction_limit
+ superad_reduction_use_turnover_limit = s% superad_reduction_use_turnover_limit
+ superad_reduction_turnover_limit_function = s% superad_reduction_turnover_limit_function
+ superad_reduction_turnover_vc_floor_frac = s% superad_reduction_turnover_vc_floor_frac
 
  max_logT_for_mlt = s% max_logT_for_mlt
  mlt_make_surface_no_mixing = s% mlt_make_surface_no_mixing
diff --git a/star/private/evolve.f90 b/star/private/evolve.f90
index dece07d13..f936707ad 100644
--- a/star/private/evolve.f90
+++ b/star/private/evolve.f90
@@ -1864,6 +1864,7 @@ integer function prepare_for_new_step(s)
                s% prev_mesh_omega(k) = s% omega(k)
                s% prev_mesh_dq(k) = s% dq(k)
                s% prev_mesh_mlt_vc(k) = s% mlt_vc(k)
+               s% prev_mesh_superad_reduction_factor(k) = s% superad_reduction_factor(k)
                s% prev_mesh_species_or_nvar_hydro_changed = .false.
             end do
             s% prev_mesh_nz = s% nz
@@ -1888,6 +1889,21 @@ integer function prepare_for_new_step(s)
             end if
          end if
 
+         ! Snapshot s%superad_reduction_factor into _old here -- AFTER mesh
+         ! adjust has remapped it onto the new mesh, but BEFORE
+         ! set_vars_if_needed below would overwrite it with the new step's
+         ! first-pass value. Used by the smooth-relaxation turnover-time
+         ! limiter as the Gamma_factor_old anchor.
+         if (s% use_superad_reduction .and. .not. s% RSP_flag) then
+            if (associated(s% superad_reduction_factor) .and. &
+                associated(s% superad_reduction_factor_old)) then
+               do k = 1, s% nz
+                  s% superad_reduction_factor_old(k) = s% superad_reduction_factor(k)
+               end do
+               s% have_superad_reduction_factor = .true.
+            end if
+         end if
+
          call set_vars_if_needed(s, s% dt_next, 'prepare_for_new_step', ierr)
          if (failed('set_vars_if_needed ierr')) return
          call set_phot_info(s)  ! this sets Teff and other stuff
@@ -2198,6 +2214,7 @@ integer function prepare_to_retry(id)
                   s% omega_old(k) = s% prev_mesh_omega(k)
                   s% j_rot_old(k) = s% prev_mesh_j_rot(k)
                   s% mlt_vc_old(k) = s% prev_mesh_mlt_vc(k)
+                  s% superad_reduction_factor_old(k) = s% prev_mesh_superad_reduction_factor(k)
                end do
                call normalize_dqs(s, s% prev_mesh_nz, s% dq_old, ierr)
                if (ierr /= 0) then
diff --git a/star/private/evolve_support.f90 b/star/private/evolve_support.f90
index 46c692e92..68ffec86f 100644
--- a/star/private/evolve_support.f90
+++ b/star/private/evolve_support.f90
@@ -58,6 +58,18 @@ subroutine new_generation(s, ierr)
             call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr)
             if (ierr /= 0) return
 
+            ! Always ensure superad_reduction_factor_old is allocated, even
+            ! when use_superad_reduction is .false., so check_sizes does not
+            ! trip on a null pointer. The actual snapshot used by the
+            ! smooth-relaxation limiter is done in evolve.f90 between the
+            ! mesh adjust and set_vars_if_needed (and the
+            ! have_superad_reduction_factor flag is set there).
+            if (associated(s% superad_reduction_factor)) then
+               call copy_to_old(s% superad_reduction_factor, &
+                                s% superad_reduction_factor_old, ierr)
+               if (ierr /= 0) return
+            end if
+
             call enlarge_if_needed_2(s% xh_old,s% nvar_hydro,nz,nz_alloc_extra,ierr)
             if (ierr /= 0) return
             if (s% fill_arrays_with_NaNs) call fill_with_NaNs_2d(s% xh_old)
diff --git a/star/private/history.f90 b/star/private/history.f90
index 66540268e..3731311e7 100644
--- a/star/private/history.f90
+++ b/star/private/history.f90
@@ -1993,6 +1993,22 @@ subroutine history_getval(&
          case(h_gradT_excess_max_lambda)
             val = s% gradT_excess_max_lambda
 
+         case(h_max_superad_reduction_factor)
+            val = 1d0
+            if (s% use_superad_reduction) then
+               do k = 1, nz
+                  if (s% superad_reduction_factor(k) > val) val = s% superad_reduction_factor(k)
+               end do
+            end if
+         case(h_num_cells_with_superad_reduction)
+            int_val = 0
+            if (s% use_superad_reduction) then
+               do k = 1, nz
+                  if (s% superad_reduction_factor(k) > 1d0) int_val = int_val + 1
+               end do
+            end if
+            is_int_val = .true.
+
          case(h_max_L_rad_div_Ledd)
             do k = 1, nz
                tmp = get_Lrad_div_Ledd(s, k)
diff --git a/star/private/init.f90 b/star/private/init.f90
index 72b2605fe..421863622 100644
--- a/star/private/init.f90
+++ b/star/private/init.f90
@@ -538,6 +538,7 @@ subroutine set_starting_star_data(s, ierr)
          s% okay_to_set_mixing_info = .true.
          s% okay_to_set_mlt_vc = .false.  ! not until have set mlt_cv_old
          s% have_mlt_vc = .false.
+         s% have_superad_reduction_factor = .false.
 
          s% have_ST_start_info = .false.
          s% prev_mesh_have_ST_start_info = .false.
diff --git a/star/private/mesh_adjust.f90 b/star/private/mesh_adjust.f90
index 3e939771a..c9e0b9ea6 100644
--- a/star/private/mesh_adjust.f90
+++ b/star/private/mesh_adjust.f90
@@ -48,6 +48,7 @@ subroutine do_mesh_adjust( &
             j_rot_old, i_rot_old, omega_old, D_omega_old, &
             mlt_vc_old, lnT_old, w_old, specific_PE_old, specific_KE_old, &
             old_m, old_r, old_rho, dPdr_dRhodr_info_old, D_mix_old, &
+            superad_reduction_factor_old, &
             cell_type, comes_from, dq_old, xq_old, xh, xa, dq, xq, ierr)
          use chem_lib, only: basic_composition_info
          use interp_1d_def
@@ -62,7 +63,8 @@ subroutine do_mesh_adjust( &
             lnd_old, lnPgas_old, mlt_vc_old, lnT_old, w_old, &
             specific_PE_old, specific_KE_old, &
             old_m, old_r, old_rho, dPdr_dRhodr_info_old, &
-            j_rot_old, omega_old, D_omega_old, D_mix_old
+            j_rot_old, omega_old, D_omega_old, D_mix_old, &
+            superad_reduction_factor_old
          type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old
          real(dp), dimension(:,:), pointer :: xh_old, xa_old
          real(dp), dimension(:,:), pointer :: xh, xa
@@ -247,6 +249,16 @@ subroutine do_mesh_adjust( &
             0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
          if (failed('mlt_cv')) return
 
+         ! Remap the previous-step superadiabatic-reduction factor onto the
+         ! new mesh so the smooth-relaxation limiter reads the correct local
+         ! Gamma_factor_old(k) after mesh adjustments. Neutral fill = 1d0
+         ! (no throttle) for cells that have no overlap with the old mesh.
+         call do_interp_pt_val( &
+            s, nz, nz_old, nzlo, nzhi, &
+            s% superad_reduction_factor, superad_reduction_factor_old, &
+            1d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
+         if (failed('superad_reduction_factor')) return
+
          call do_interp_pt_val( &
             s, nz, nz_old, nzlo, nzhi, s% D_mix, D_mix_old, &
             0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
@@ -556,6 +568,7 @@ subroutine do_prune_mesh_surface( &
             mlt_vc_old, lnT_old, &
             dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, &
             D_SSI_old, D_ES_old, D_GSF_old, D_mix_old, &
+            superad_reduction_factor_old, &
             xh, xa, ierr)
          use auto_diff_support
          type (star_info), pointer :: s
@@ -564,7 +577,8 @@ subroutine do_prune_mesh_surface( &
             j_rot_old, omega_old, &
             D_omega_old, am_nu_rot_old, mlt_vc_old, lnT_old, &
             dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, &
-            D_SSI_old, D_ES_old, D_GSF_old, D_mix_old
+            D_SSI_old, D_ES_old, D_GSF_old, D_mix_old, &
+            superad_reduction_factor_old
          type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old
          real(dp), dimension(:,:), pointer :: xh_old, xa_old
          real(dp), dimension(:,:), pointer :: xh, xa
@@ -595,6 +609,7 @@ subroutine do_prune_mesh_surface( &
          call prune1(s% lnT, lnT_old, skip)
          call prune1(s% D_mix, D_mix_old, skip)
          call prune1(s% mlt_vc, mlt_vc_old, skip)
+         call prune1(s% superad_reduction_factor, superad_reduction_factor_old, skip)
 
          if (s% rotation_flag) then
             call prune1(s% j_rot, j_rot_old, skip)
diff --git a/star/private/profile_getval.f90 b/star/private/profile_getval.f90
index 471399f94..9e2cf4bfe 100644
--- a/star/private/profile_getval.f90
+++ b/star/private/profile_getval.f90
@@ -499,6 +499,11 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val)
 
             case (p_superad_reduction_factor)
                val = s% superad_reduction_factor(k)
+            case (p_superad_reduction_Lrad_div_Ledd)
+               val = s% superad_reduction_Lrad_div_Ledd(k)
+            case (p_superad_reduction_trigger)
+               int_val = s% superad_reduction_trigger(k)
+               int_flag = .true.
             case (p_gradT_excess_effect)
                val = s% gradT_excess_effect(k)
             case (p_diff_grads)
diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90
index 689974af6..4f404ea28 100644
--- a/star/private/remove_shells.f90
+++ b/star/private/remove_shells.f90
@@ -1111,6 +1111,7 @@ subroutine do_remove_surface(id, surface_k, ierr)
             prv% mlt_vc, prv% lnT, &
             prv% dPdr_dRhodr_info, prv% nu_ST, prv% D_ST, prv% D_DSI, prv% D_SH, &
             prv% D_SSI, prv% D_ES, prv% D_GSF, prv% D_mix, &
+            prv% superad_reduction_factor, &
             s% xh, s% xa, ierr)
          if (ierr /= 0) then
             return
@@ -1234,6 +1235,7 @@ subroutine do_relax_to_star_cut( &
          ! save have_mlt_vc and set to false (to load ZAMS model)
          save_have_mlt_vc = s% have_mlt_vc
          s% have_mlt_vc = .false.
+         s% have_superad_reduction_factor = .false.
 
          !save composition and entropy profiles
          xa(:,:) = s% xa(:,k_remove:s% nz)
diff --git a/star/private/star_history_def.f90 b/star/private/star_history_def.f90
index e569989d9..b87f2426d 100644
--- a/star/private/star_history_def.f90
+++ b/star/private/star_history_def.f90
@@ -402,7 +402,10 @@ module star_history_def
       integer, parameter :: h_gradT_excess_max_lambda = h_gradT_excess_min_beta + 1
       integer, parameter :: h_gradT_excess_alpha = h_gradT_excess_max_lambda + 1
 
-      integer, parameter :: h_log_Ledd = h_gradT_excess_alpha + 1
+      integer, parameter :: h_max_superad_reduction_factor = h_gradT_excess_alpha + 1
+      integer, parameter :: h_num_cells_with_superad_reduction = h_max_superad_reduction_factor + 1
+
+      integer, parameter :: h_log_Ledd = h_num_cells_with_superad_reduction + 1
       integer, parameter :: h_compactness = h_log_Ledd + 1
       integer, parameter :: h_compactness_parameter = h_compactness + 1
       integer, parameter :: h_mu4 = h_compactness_parameter + 1
@@ -915,6 +918,9 @@ subroutine history_column_names_init(ierr)
          history_column_name(h_gradT_excess_min_beta) = 'gradT_excess_min_beta'
          history_column_name(h_gradT_excess_max_lambda) = 'gradT_excess_max_lambda'
 
+         history_column_name(h_max_superad_reduction_factor) = 'max_superad_reduction_factor'
+         history_column_name(h_num_cells_with_superad_reduction) = 'num_cells_with_superad_reduction'
+
          history_column_name(h_gamma1_min) = 'gamma1_min'
          history_column_name(h_logT_max) = 'logT_max'
          history_column_name(h_logQ_max) = 'logQ_max'
diff --git a/star/private/star_profile_def.f90 b/star/private/star_profile_def.f90
index a51dd25c6..bfd722683 100644
--- a/star/private/star_profile_def.f90
+++ b/star/private/star_profile_def.f90
@@ -72,7 +72,9 @@ module star_profile_def
       integer, parameter :: p_vel_km_per_s = p_v_kms + 1
       integer, parameter :: p_log_abs_v = p_vel_km_per_s + 1
       integer, parameter :: p_superad_reduction_factor = p_log_abs_v + 1
-      integer, parameter :: p_gradT_excess_effect = p_superad_reduction_factor + 1
+      integer, parameter :: p_superad_reduction_Lrad_div_Ledd = p_superad_reduction_factor + 1
+      integer, parameter :: p_superad_reduction_trigger = p_superad_reduction_Lrad_div_Ledd + 1
+      integer, parameter :: p_gradT_excess_effect = p_superad_reduction_trigger + 1
       integer, parameter :: p_log_diff_grads = p_gradT_excess_effect + 1
       integer, parameter :: p_diff_grads = p_log_diff_grads + 1
 
@@ -773,6 +775,8 @@ subroutine profile_column_names_init(ierr)
          profile_column_name(p_diff_grads) = 'diff_grads'
          profile_column_name(p_gradT_excess_effect) = 'gradT_excess_effect'
          profile_column_name(p_superad_reduction_factor) = 'superad_reduction_factor'
+         profile_column_name(p_superad_reduction_Lrad_div_Ledd) = 'superad_reduction_Lrad_div_Ledd'
+         profile_column_name(p_superad_reduction_trigger) = 'superad_reduction_trigger'
 
          profile_column_name(p_v) = 'v'
          profile_column_name(p_velocity) = 'velocity'
diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90
index 67eda0362..834d5ceec 100644
--- a/star/private/turb_support.f90
+++ b/star/private/turb_support.f90
@@ -221,6 +221,8 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
 
       ! these are used by use_superad_reduction
       real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit
+      real(dp) :: vc_old_local, vc_old_floor, Gamma_factor_old_local
+      type(auto_diff_real_star_order1) :: tau_conv, f_turnover
       type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, &
          diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq, mlt_Pturb, Ptot
       logical ::  test_partials, using_TDC, have_Y_face_guess
@@ -287,7 +289,11 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
       conv_vel = 0d0
       D = 0d0
       Gamma = 0d0
-      if (k /= 0) s% superad_reduction_factor(k) = 1d0
+      if (k /= 0) then
+         s% superad_reduction_factor(k) = 1d0
+         s% superad_reduction_Lrad_div_Ledd(k) = 0d0
+         s% superad_reduction_trigger(k) = 0
+      end if
 
       ! Bail if we asked for no mixing, or if parameters are bad.
       if (MLT_option == 'none' .or. beta < 1d-10 .or. mixing_length_alpha <= 0d0 .or. &
@@ -465,6 +471,11 @@ subroutine set_superad_reduction()
          Lrad_div_Ledd = 4d0*crad/3d0*pow4(T)/P*gradT
          Gamma_inv_threshold = 4d0*(1d0-beta)/(4d0-3*beta)
 
+         if (k /= 0) then
+            s% superad_reduction_Lrad_div_Ledd(k) = Lrad_div_Ledd% val
+            s% superad_reduction_trigger(k) = 0
+         end if
+
          Gamma_factor = 1d0
          if (gradT > gradL) then
             if (Lrad_div_Ledd > Gamma_limit .or. Lrad_div_Ledd > Gamma_inv_threshold) then
@@ -483,6 +494,7 @@ subroutine set_superad_reduction()
                !   Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
                !end if
                if (Lrad_div_Ledd > Gamma_limit) then
+                  if (k /= 0) s% superad_reduction_trigger(k) = s% superad_reduction_trigger(k) + 1
                   alfa0 = Lrad_div_Ledd/Gamma_limit-1d0
                   if (alfa0 < 1d0) then
                      Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
@@ -492,11 +504,12 @@ subroutine set_superad_reduction()
                   !Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0)
                end if
                if (Lrad_div_Ledd% val > Gamma_inv_threshold) then
+                  if (k /= 0) s% superad_reduction_trigger(k) = s% superad_reduction_trigger(k) + 2
                   alfa0 = Lrad_div_Ledd/Gamma_inv_threshold-1d0
                   if (alfa0 < 1d0) then
-                     Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
+                     Gamma_term = Gamma_term + scale_value2*(0.5d0*alfa0*alfa0)
                   else
-                     Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0)
+                     Gamma_term = Gamma_term + scale_value2*(alfa0-0.5d0)
                   end if
                   !Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
                end if
@@ -512,6 +525,53 @@ subroutine set_superad_reduction()
                end if
             end if
          end if
+
+         ! Convective-turnover-time limiter (placed OUTSIDE the
+         ! Gamma_term>0d0 block so it also handles the case where the
+         ! instantaneous Gamma_factor is 1 but the previous step's value
+         ! was > 1 -- i.e., we must relax the throttle DOWN, not snap it
+         ! to 1 in a single step). Relaxes Gamma_factor toward its
+         ! instantaneous (capped) value over the local turnover time
+         ! tau_conv = scale_height / mlt_vc_old, anchored at the
+         ! previous-step's converged Gamma_factor:
+         !
+         !    d Gamma_fac / dt = (Gamma_fac_inst - Gamma_fac) / tau_conv
+         !    Gamma_fac_new = Gamma_fac_old + f_turnover * (Gamma_fac_inst - Gamma_fac_old)
+         !
+         ! dt >> tau_conv: f_turnover -> 1, Gamma_fac -> Gamma_fac_inst (no smoothing).
+         ! dt << tau_conv: f_turnover -> 0, Gamma_fac -> Gamma_fac_old (convection
+         !   had no time to adapt, throttle held at previous-step value
+         !   independently of whether the new step needs more or less throttling).
+         !
+         ! Skip until the previous-step mlt_vc and Gamma_factor_old have both
+         ! been populated. Fires only when there is throttling to track --
+         ! either an instantaneous Gamma_factor>1 OR a non-trivial old value
+         ! that has to relax back to 1.
+         if (s% superad_reduction_use_turnover_limit .and. k > 0 .and. &
+             s% have_mlt_vc .and. associated(s% mlt_vc_old) .and. &
+             s% have_superad_reduction_factor .and. &
+             associated(s% superad_reduction_factor_old) .and. &
+             s% dt > 0d0) then
+            ! Clamp at >= 1 so we never relax from a non-physical sub-1 anchor.
+            Gamma_factor_old_local = max(s% superad_reduction_factor_old(k), 1d0)
+            if (Gamma_factor > 1d0 .or. Gamma_factor_old_local > 1d0) then
+               ! Optional floor on mlt_vc_old at a fraction of the local face
+               ! sound speed; see controls.defaults for the physical motivation.
+               vc_old_floor = s% superad_reduction_turnover_vc_floor_frac &
+                              * s% csound_face(k)
+               vc_old_local = max(s% mlt_vc_old(k), vc_old_floor, 1d-30)
+               tau_conv = scale_height / vc_old_local
+               select case (trim(s% superad_reduction_turnover_limit_function))
+               case ('linear')
+                  f_turnover = min(s% dt / tau_conv, 1d0)
+               case default
+                  f_turnover = 1d0 - exp(-s% dt / tau_conv)
+               end select
+               Gamma_factor = Gamma_factor_old_local + &
+                              f_turnover * (Gamma_factor - Gamma_factor_old_local)
+            end if
+         end if
+
          if (k /= 0) s% superad_reduction_factor(k) = Gamma_factor% val
          if (Gamma_factor > 1d0) then
             grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr
diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc
index 2d0c29d8d..56e8b19bd 100644
--- a/star_data/private/star_controls.inc
+++ b/star_data/private/star_controls.inc
@@ -306,6 +306,9 @@
             superad_reduction_Gamma_inv_scale, &
             superad_reduction_diff_grads_limit, &
             superad_reduction_limit
+         logical :: superad_reduction_use_turnover_limit
+         character(len=strlen) :: superad_reduction_turnover_limit_function
+         real(dp) :: superad_reduction_turnover_vc_floor_frac
 
       ! mixing parameters
 
diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc
index 1991b965d..f355021a8 100644
--- a/star_data/public/star_data_step_work.inc
+++ b/star_data/public/star_data_step_work.inc
@@ -333,7 +333,8 @@
       real(dp), pointer, dimension(:,:) :: prev_mesh_xh ! (nvar_hydro,prev_mesh_nz)
       real(dp), pointer, dimension(:,:) :: prev_mesh_xa ! (species,prev_mesh_nz)
       real(dp), pointer, dimension(:) :: & ! (prev_mesh_nz)
-         prev_mesh_j_rot, prev_mesh_omega, prev_mesh_dq, prev_mesh_mlt_vc
+         prev_mesh_j_rot, prev_mesh_omega, prev_mesh_dq, prev_mesh_mlt_vc, &
+         prev_mesh_superad_reduction_factor
       logical :: prev_mesh_species_or_nvar_hydro_changed ! specifies if either species or number of hydro variables
                                                          ! have been altered since prev_mesh info was stored
       real(dp), pointer, dimension(:) :: prev_mesh_D_ST_start, prev_mesh_nu_ST_start
@@ -355,6 +356,19 @@
       real(dp), pointer :: superad_reduction_factor(:) ! only set when use_superad_reduction = .true.
          ! reports effect of superad_reduction.  value = 1d0 for no effect.
 
+      real(dp), pointer :: superad_reduction_Lrad_div_Ledd(:)
+         ! only set when use_superad_reduction = .true.
+         ! reports Lrad/Ledd = (4 a T^4 / 3 P) * gradT, computed from the
+         ! first-pass (unthrottled) gradT inside set_superad_reduction.
+         ! Diagnostic only; not used by the solver.
+
+      integer, pointer :: superad_reduction_trigger(:)
+         ! only set when use_superad_reduction = .true.
+         ! 0 = throttle inactive in this cell.
+         ! 1 = activated by Eddington threshold (Lrad/Ledd > Gamma_limit).
+         ! 2 = activated by density-inversion threshold only.
+         ! 3 = both triggers active.
+
       ! mixing length alpha for MLT (can be set by user and can vary within model)
       real(dp), pointer :: alpha_mlt(:)
 
@@ -557,7 +571,14 @@
 
          real(dp), pointer, dimension(:,:) :: xh_old, xa_old
          real(dp), pointer, dimension(:) :: &
-            dq_old, q_old, j_rot_old, omega_old, mlt_vc_old
+            dq_old, q_old, j_rot_old, omega_old, mlt_vc_old, &
+            superad_reduction_factor_old
+            ! superad_reduction_factor_old: previous step's converged
+            ! Gamma_factor (= 1 when no throttle was active). Used by the
+            ! turnover-time limiter to relax Gamma_factor smoothly toward
+            ! the new instantaneous target rather than snapping it back to 1
+            ! when dt collapses below tau_conv (which would cause an
+            ! unphysical HRD jump). Companion to mlt_vc_old.
          real(dp) :: time_old, dt_old, mstar_old, xmstar_old, mstar_dot_old, &
             M_center_old, v_center_old, R_center_old, L_center_old
          integer :: nz_old, model_number_old
@@ -893,7 +914,12 @@
             RSP_just_set_velocities, &
             have_new_cz_bdy_info, using_gold_tolerances, need_to_setvars, have_new_generation, &
             using_velocity_time_centering, okay_to_set_mixing_info, &
-            okay_to_set_mlt_vc, have_mlt_vc
+            okay_to_set_mlt_vc, have_mlt_vc, &
+            have_superad_reduction_factor
+                                            ! true once a step has set
+                                            ! superad_reduction_factor and we've
+                                            ! copied it into _old (consumed by
+                                            ! the smooth-relaxation limiter)
 
          real(dp), pointer :: xa_removed(:) ! mass fractions for removed mass (1:species)
          real(dp) :: h1_czb_mass