diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 3ed6b8e2ee..e6f5d4abb5 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -41,6 +41,8 @@ GYRE has been upgraded to 9.1.1, the most recent stable release. Changes since t The `MESA SDK `__ recommended for compiling MESA has been updated to 26.6.1. Although this newer SDK is not required to successfully build MESA, it brings the benefit of restoring the cross-platform bit-for-bit compatibility that MESA once enjoyed (meaning that runs on Linux/Intel, MacOS/Intel and MacOS/ARM give identical results). +The new control ``superad_reduction_use_turnover_limit`` relaxes the applied superadiabatic reduction from its previous accepted value toward the instantaneous value. ``superad_reduction_turnover_limit_function`` selects either the exponential response ``1-exp(-dt/tau_conv)`` or the linear response ``min(dt/tau_conv,1)``. The limiter acts on the applied reduction ``1/Gamma_factor``. Zones with a lagged convective velocity use ``scale_height/max(mlt_vc,1d-10 cm/s)``; other zones use the Brunt frequency. The timescale is set at the start of the step and held fixed during solver iterations. With ``use_face_reconstruction``, this calculation uses the reconstructed face thermodynamic state. The scale height is the interpolated or reconstructed face value used by MLT and TDC. The previous reduction is preserved across retries, remeshing, and photo restarts. For ``k > 0``, ``superad_reduction_max_logT`` restricts the reduction to faces whose start-of-step temperature is below the selected ``logT``. Its default is ``7d0``, corresponding to :math:`10^7\,\mathrm{K}`. The ``k=0`` model-construction path is unchanged. + Additional controls are available for TDC envelope remeshing: - ``remesh_for_TDC_pulsations_when_load`` remeshes a model after loading it from a ``.mod`` file. @@ -68,6 +70,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`` used ``superad_reduction_Gamma_limit_scale`` instead of ``superad_reduction_Gamma_inv_scale``. This made the latter control inactive. The published algorithm is restored (Jermyn et al. 2023, equation 64, ``alpha_2`` term). Existing test suite cases set the two scales to the same value and are bit-for-bit unaffected. The fix changes results only when the scales differ. + +Photo restarts now preserve the stored superadiabatic turnover-limiter state while ``finish_load_model`` rebuilds derived quantities. Previously, the stored reduction could be reset before ``extras_startup`` restored controls selected dynamically by ``run_star_extras``. Runs that leave superadiabatic reduction disabled still reset the factor on the first evolution step. + .. 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``). Other Changes diff --git a/linters/check_photos.py b/linters/check_photos.py index 5ae5d24664..dabee54228 100755 --- a/linters/check_photos.py +++ b/linters/check_photos.py @@ -54,6 +54,9 @@ def discard(self, value): "net_iso", "model_profile_filename", "include_binary_history_in_log_file", + "superad_reduction_factor_old", # Initialized from superad_reduction_factor + "have_superad_reduction_factor", # Set when superad_reduction_factor is read + "okay_to_set_superad_reduction_factor", # Reset during photo loading "i_dj_rot_dt", # Shadows s% i_j_rot "i_equ_w_div_wc", # Shadows s% i_w_div_wc "have_j_rot", # Set from s% rotation_flag diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index ebd7b06dd0..8c77b04e01 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2658,29 +2658,113 @@ ! use_superad_reduction ! ~~~~~~~~~~~~~~~~~~~~~ + + ! Implicit alternative to ``okay_to_reduce_gradT_excess``. MLT or TDC + ! is repeated with ``gradr`` moved toward ``gradL`` in cells near the + ! Eddington or density-inversion limits. The two contributions are + ! evaluated separately and added to ``Gamma_factor``, with + ! ``gradr_scaled = gradL + (gradr-gradL)/Gamma_factor``. + + ! :: + + use_superad_reduction = .false. + + + ! superad_reduction_max_logT + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! For k > 0, apply superadiabatic reduction only where the start-of-step + ! ``logT`` is less than this limit. The start-of-step value keeps the + ! selection fixed during solver iterations. + + ! :: + + superad_reduction_max_logT = 7d0 + + ! superad_reduction_Gamma_limit ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Threshold in ``Lrad/Ledd`` for the Eddington contribution. + + ! :: + + superad_reduction_Gamma_limit = 0.5d0 + + ! superad_reduction_Gamma_limit_scale ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Scale factor for the Eddington contribution. Set to zero to disable + ! this contribution. + + ! :: + + superad_reduction_Gamma_limit_scale = 5d0 + + ! superad_reduction_Gamma_inv_scale ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Scale factor for the density-inversion contribution, which activates + ! when ``Lrad/Ledd > 4*(1-beta)/(4-3*beta)``. Set to zero to disable + ! this contribution. + + ! :: + + superad_reduction_Gamma_inv_scale = 5d0 + + ! superad_reduction_diff_grads_limit ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Smoothly turns off the reduction as ``gradT - gradL`` approaches + ! this value. + + ! :: + + superad_reduction_diff_grads_limit = 1d-3 + + ! superad_reduction_limit ! ~~~~~~~~~~~~~~~~~~~~~~~ - ! Implicit alternative to ``okay_to_reduce_gradT_excess`` + ! Upper limit for ``Gamma_factor``. A value less than or equal to one + ! disables the limit. ! :: - 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_use_turnover_limit + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Relaxes the applied reduction from its previous accepted value toward + ! the instantaneous value over the local response time. The response + ! fraction is set by ``superad_reduction_turnover_limit_function``. + ! In zones with a lagged convective velocity, ``tau_conv`` is + ! ``scale_height/max(mlt_vc,1d-10 cm/s)``. Other zones use the Brunt + ! frequency. The timescale is set at the start of the step and held fixed + ! during solver iterations. The limiter acts on ``1/Gamma_factor``, which + ! multiplies ``gradr-gradL``, and uses the same interpolated or + ! reconstructed face state as MLT or TDC. + + ! :: + + superad_reduction_use_turnover_limit = .false. + + ! superad_reduction_turnover_limit_function + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Sets the response fraction used by the turnover limiter. The + ! ``'exponential'`` option uses ``1-exp(-dt/tau_conv)`` and the + ! ``'linear'`` option uses ``min(dt/tau_conv,1)``. + + ! :: + + superad_reduction_turnover_limit_function = 'linear' + ! overshooting ! ____________ diff --git a/star/private/adjust_mesh.f90 b/star/private/adjust_mesh.f90 index a243f2442f..64bd39c99f 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,8 @@ 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/adjust_mesh_split_merge.f90 b/star/private/adjust_mesh_split_merge.f90 index 906e3aa85e..a6d951b5b4 100644 --- a/star/private/adjust_mesh_split_merge.f90 +++ b/star/private/adjust_mesh_split_merge.f90 @@ -657,6 +657,7 @@ subroutine do_merge(s, i_merge, species, new_xa, ierr) s% lnT(im) = s% lnT(i0) s% D_mix(im) = s% D_mix(i0) s% mlt_vc(im) = s% mlt_vc(i0) + s% superad_reduction_factor(im) = s% superad_reduction_factor(i0) s% csound(im) = s% csound(i0) s% tau(im) = s% tau(i0) s% opacity(im) = s% opacity(i0) @@ -831,7 +832,8 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr) sumx, sumxp, new_xaL, new_xaR, star_PE0, star_PE1, & grad_alpha, f, new_alphaL, new_alphaR, v_R, v_C, v_L, min_dm, & mlt_vcL, mlt_vcR, tauL, tauR, etrb, etrb_L, etrb_C, etrb_R, grad_etrb, & - j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old + j_rot_new, dmbar_old, dmbar_p1_old, dmbar_new, dmbar_p1_new, dmbar_p2_new, J_old, & + superad_reduction_factorL, superad_reduction_factorR logical :: done, use_new_grad_rho include 'formats' @@ -868,13 +870,16 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr) rR = s% r(i) mlt_vcR = s% mlt_vc(i) + superad_reduction_factorR = s% superad_reduction_factor(i) if (i == nz) then rL = s% R_center mlt_vcL = 0d0 + superad_reduction_factorL = 1d0 tauL = tau_center else rL = s% r(ip) mlt_vcL = s% mlt_vc(ip) + superad_reduction_factorL = s% superad_reduction_factor(ip) tauL = s% tau(ip) end if @@ -1038,6 +1043,7 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr) s% lnT(jp) = s% lnT(j) s% D_mix(jp) = s% D_mix(j) s% mlt_vc(jp) = s% mlt_vc(j) + s% superad_reduction_factor(jp) = s% superad_reduction_factor(j) s% csound(jp) = s% csound(j) s% tau(jp) = s% tau(j) s% opacity(jp) = s% opacity(j) @@ -1175,9 +1181,11 @@ subroutine do_split(s, i_split, species, tau_center, grad_xa, new_xa, ierr) s% dPdr_dRhodr_info(ip) = s% dPdr_dRhodr_info(i) end if - ! mlt_vc is face-based, so a split creates a new interior face value here. + ! These are face-based, so a split creates new interior face values here. ! Interpolate using the same left/right orientation as tau(ip). s% mlt_vc(ip) = mlt_vcR + (mlt_vcL - mlt_vcR)*dMR/dM + s% superad_reduction_factor(ip) = superad_reduction_factorR + & + (superad_reduction_factorL - superad_reduction_factorR)*dMR/dM s% tau(ip) = tauR + (tauL - tauR)*dMR/dM if (is_bad(s% tau(ip))) then diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 53f5ac571c..088503f642 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,7 @@ 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 + if (action == do_allocate) s% superad_reduction_factor = 1d0 call do1(s% domega_dlnR, c% domega_dlnR) if (failed('domega_dlnR')) exit @@ -1163,6 +1166,9 @@ subroutine star_info_arrays(s, c_in, action_in, ierr) if (failed('unsmoothed_brunt_B')) exit call do1(s% smoothed_brunt_B, c% smoothed_brunt_B) if (failed('smoothed_brunt_B')) exit + call do1(s% tau_conv_start, c% tau_conv_start) + if (failed('tau_conv_start')) exit + if (action == do_allocate) s% tau_conv_start = 0d0 call do1(s% RTI_du_diffusion_kick, c% RTI_du_diffusion_kick) if (failed('RTI_du_diffusion_kick')) exit @@ -1427,6 +1433,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 921befa107..f5c06c51d4 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -138,8 +138,10 @@ module ctrls_io gradT_excess_beta1, gradT_excess_lambda2, gradT_excess_beta2, & gradT_excess_dlambda, gradT_excess_dbeta, gradT_excess_max_center_h1, & 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, & + use_superad_reduction, superad_reduction_max_logT, & + 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, & 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, & @@ -656,6 +658,16 @@ subroutine check_controls(s, ierr) return end if + if (.not. (trim(s% superad_reduction_turnover_limit_function) == 'linear' .or. & + trim(s% superad_reduction_turnover_limit_function) == 'exponential')) then + write(*,'(A)') + write(*,*) "Invalid choice for superad_reduction_turnover_limit_function" + write(*,*) "Available options are 'linear' or 'exponential'" + write(*,'(A)') + ierr = -1 + return + end if + end subroutine check_controls @@ -1041,11 +1053,14 @@ subroutine store_controls(s) s% dq_D_mix_zero_at_H_C_crossover = dq_D_mix_zero_at_H_C_crossover s% use_superad_reduction = use_superad_reduction + s% superad_reduction_max_logT = superad_reduction_max_logT s% superad_reduction_gamma_limit = superad_reduction_gamma_limit s% superad_reduction_gamma_limit_scale = superad_reduction_gamma_limit_scale 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% max_logT_for_mlt = max_logT_for_mlt s% mlt_make_surface_no_mixing = mlt_make_surface_no_mixing @@ -2769,11 +2784,14 @@ subroutine set_controls_for_writing(s, ierr) dq_D_mix_zero_at_H_C_crossover = s% dq_D_mix_zero_at_H_C_crossover use_superad_reduction = s% use_superad_reduction + superad_reduction_max_logT = s% superad_reduction_max_logT superad_reduction_gamma_limit = s% superad_reduction_gamma_limit superad_reduction_gamma_limit_scale = s% superad_reduction_gamma_limit_scale 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 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 dece07d13d..0f4110886a 100644 --- a/star/private/evolve.f90 +++ b/star/private/evolve.f90 @@ -319,6 +319,7 @@ integer function do_step_part1(id, first_try) s% need_to_setvars = .true. ! always start fresh s% okay_to_set_mixing_info = .true. ! set false by element diffusion s% okay_to_set_mlt_vc = .false. ! don't change mlt_vc until have set mlt_vc_old + s% okay_to_set_superad_reduction_factor = .false. if (s% timestep_hold > s% model_number + 10000) then write(*,3) 'ERROR: s% timestep_hold', s% timestep_hold, s% model_number @@ -698,6 +699,7 @@ integer function do_step_part2(id, first_try) s% have_mlt_vc = .true. end if s% okay_to_set_mlt_vc = .false. + s% okay_to_set_superad_reduction_factor = .false. end if if (.not. okay_energy_conservation()) return @@ -1864,6 +1866,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 @@ -2198,6 +2201,8 @@ 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 46c692e921..e90a91b2fe 100644 --- a/star/private/evolve_support.f90 +++ b/star/private/evolve_support.f90 @@ -58,6 +58,12 @@ subroutine new_generation(s, ierr) call copy_to_old(s% mlt_vc, s% mlt_vc_old, ierr) if (ierr /= 0) return + call copy_to_old(s% superad_reduction_factor, & + s% superad_reduction_factor_old, ierr) + if (ierr /= 0) return + s% have_superad_reduction_factor = .true. + s% okay_to_set_superad_reduction_factor = .true. + 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) @@ -191,8 +197,10 @@ subroutine set_current_to_old(s) end do s% dq(k) = s% dq_old(k) s% mlt_vc(k) = s% mlt_vc_old(k) + s% superad_reduction_factor(k) = s% superad_reduction_factor_old(k) end do s% okay_to_set_mlt_vc = .true. + s% okay_to_set_superad_reduction_factor = .true. call set_qs(s, s% nz, s% q, s% dq, ierr) if (ierr /= 0) then diff --git a/star/private/hydro_vars.f90 b/star/private/hydro_vars.f90 index 9c93035bea..af9865614b 100644 --- a/star/private/hydro_vars.f90 +++ b/star/private/hydro_vars.f90 @@ -542,6 +542,10 @@ subroutine set_hydro_vars( & if (failed('set_rotation_info')) return end if + ! Invalidate before MLT rebuilds the face cache for the current state. + if (.not. skip_mlt .and. .not. s% RSP_flag) & + s% reconstructed_face_state_valid(1:s%nz) = .false. + if (.not. skip_grads) then if (dbg) write(*,*) 'call do_brunt_B' call do_brunt_B(s, nzlo, nzhi, ierr) ! for unsmoothed_brunt_B @@ -563,9 +567,6 @@ subroutine set_hydro_vars( & end if if (.not. skip_mlt .and. .not. s% RSP_flag) then - - s% reconstructed_face_state_valid(1:s%nz) = .false. - if (.not. skip_mixing_info) then if (s% make_gradr_sticky_in_solver_iters) then s% fixed_gradr_for_rest_of_solver_iters(nzlo:nzhi) = .false. @@ -611,7 +612,8 @@ subroutine set_hydro_vars( & end if end if - call set_mlt_vars(s, nzlo, nzhi, ierr) + call set_mlt_vars(s, nzlo, nzhi, ierr, & + set_tau_conv = .not. skip_grads) if (failed('set_mlt_vars')) return if (dbg) write(*,*) 'call check_for_redo_MLT' diff --git a/star/private/init.f90 b/star/private/init.f90 index 72b2605fe0..59e09a03a4 100644 --- a/star/private/init.f90 +++ b/star/private/init.f90 @@ -538,6 +538,8 @@ 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% okay_to_set_superad_reduction_factor = .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 78c5d47bec..121b579970 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,12 @@ subroutine do_mesh_adjust( & 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr) if (failed('mlt_cv')) return + 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 +564,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 +573,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 +605,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/photo_in.f90 b/star/private/photo_in.f90 index aa8879bf5c..c4efa8401f 100644 --- a/star/private/photo_in.f90 +++ b/star/private/photo_in.f90 @@ -134,8 +134,13 @@ subroutine read_star_photo(s, fname, ierr) read(iounit, iostat=ierr) & s% dq(1:nz), s% xa(:,1:nz), s% xh(:,1:nz), & s% omega(1:nz), s% j_rot(1:nz), s% mlt_vc(1:nz), s% conv_vel(1:nz), & + s% superad_reduction_factor(1:nz), & s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing s% have_ST_start_info + if (failed('superad_reduction_factor')) return + s% superad_reduction_factor_old(1:nz) = s% superad_reduction_factor(1:nz) + s% have_superad_reduction_factor = .true. + s% okay_to_set_superad_reduction_factor = .false. call read_part_number(iounit) if (failed('rsp_num_periods')) return diff --git a/star/private/photo_out.f90 b/star/private/photo_out.f90 index 5967e184e2..6222abf8c1 100644 --- a/star/private/photo_out.f90 +++ b/star/private/photo_out.f90 @@ -77,6 +77,7 @@ subroutine output_star_photo(s,iounit,ierr) write(iounit) & s% dq(1:nz), s% xa(:,1:nz), s% xh(:,1:nz), & s% omega(1:nz), s% j_rot(1:nz), s% mlt_vc(1:nz), s% conv_vel(1:nz), & + s% superad_reduction_factor(1:nz), & s% D_ST_start(1:nz), s% nu_ST_start(1:nz), & ! needed for ST time smoothing s% have_ST_start_info diff --git a/star/private/remove_shells.f90 b/star/private/remove_shells.f90 index 689974af6b..e702c4b276 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,8 @@ 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. + s% okay_to_set_superad_reduction_factor = .false. !save composition and entropy profiles xa(:,:) = s% xa(:,k_remove:s% nz) diff --git a/star/private/star_utils.f90 b/star/private/star_utils.f90 index 051678cdf6..4ad8f8ef57 100644 --- a/star/private/star_utils.f90 +++ b/star/private/star_utils.f90 @@ -3709,18 +3709,26 @@ real(dp) function conv_time_scale(s,k_in) result(tau_conv) end if k = max(2,k_in) brunt_B = s% brunt_B(k) - call get_face_weights(s, k, alfa, beta) - rho_face = alfa*s% rho(k) + beta*s% rho(k-1) - Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1) - chiT_face = alfa*s% chiT(k) + beta*s% chiT(k-1) - chiRho_face = alfa*s% chiRho(k) + beta*s% chiRho(k-1) + if (s% use_face_reconstruction .and. s% reconstructed_face_state_valid(k)) then + rho_face = s% reconstructed_rho_face_ad(k)% val + Peos_face = s% reconstructed_P_face_ad(k)% val + chiT_face = s% reconstructed_ChiT_face_ad(k)% val + chiRho_face = s% reconstructed_ChiRho_face_ad(k)% val + grada_face = s% reconstructed_grada_face_ad(k)% val + else + call get_face_weights(s, k, alfa, beta) + rho_face = alfa*s% rho(k) + beta*s% rho(k-1) + Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1) + chiT_face = alfa*s% chiT(k) + beta*s% chiT(k-1) + chiRho_face = alfa*s% chiRho(k) + beta*s% chiRho(k-1) + grada_face = alfa*s% grada(k) + beta*s% grada(k-1) + end if f = pow2(s% grav(k))*rho_face/Peos_face*chiT_face/chiRho_face dlnP = s% lnPeos(k-1) - s% lnPeos(k) dlnT = s% lnT(k-1) - s% lnT(k) - grada_face = alfa*s% grada(k) + beta*s% grada(k-1) gradT_actual = safe_div_val(s, dlnT, dlnP) ! mlt has not been called yet when doing this brunt_N2 = f*(brunt_B - (gradT_actual - grada_face)) - if(abs(brunt_B) > 0d0) then + if (abs(brunt_N2) > 0d0) then tau_conv = 1d0/sqrt(abs(brunt_N2)) else tau_conv = 0d0 @@ -3732,15 +3740,21 @@ subroutine set_conv_time_scales(s) type (star_info), pointer :: s integer :: k real(dp) :: tau_conv + logical :: set_tau_conv include 'formats' s% min_conv_time_scale = 1d99 s% max_conv_time_scale = 0d0 + set_tau_conv = s% superad_reduction_use_turnover_limit do k=1,s%nz + if (set_tau_conv) then + tau_conv = conv_time_scale(s,k) + s% tau_conv_start(k) = tau_conv + end if if (s% X(k) > s% max_X_for_conv_timescale) cycle if (s% X(k) < s% min_X_for_conv_timescale) cycle if (s% q(k) > s% max_q_for_conv_timescale) cycle if (s% q(k) < s% min_q_for_conv_timescale) exit - tau_conv = conv_time_scale(s,k) + if (.not. set_tau_conv) tau_conv = conv_time_scale(s,k) if (tau_conv < s% min_conv_time_scale) & s% min_conv_time_scale = tau_conv if (tau_conv > s% max_conv_time_scale) & diff --git a/star/private/turb_info.f90 b/star/private/turb_info.f90 index b88c314cf6..4c98d2684f 100644 --- a/star/private/turb_info.f90 +++ b/star/private/turb_info.f90 @@ -36,24 +36,37 @@ module turb_info public :: check_for_redo_MLT ! for hydro_vars public :: set_gradT_excess_alpha ! for evolve + real(dp), parameter :: min_mlt_vc_for_tau_conv = 1d-10 ! cm/s + contains - subroutine set_mlt_vars(s, nzlo, nzhi, ierr) + subroutine set_mlt_vars(s, nzlo, nzhi, ierr, set_tau_conv) use star_utils, only: start_time, update_time type (star_info), pointer :: s integer, intent(in) :: nzlo, nzhi integer, intent(out) :: ierr + logical, intent(in), optional :: set_tau_conv integer :: k, op_err integer(i8) :: time0 real(dp) :: total - logical :: make_gradr_sticky_in_solver_iters + logical :: make_gradr_sticky_in_solver_iters, set_tau, use_mlt_vc_tau_surface include 'formats' ierr = 0 + ! Keep the start-state timescale fixed during solver iterations. + set_tau = .false. + if (present(set_tau_conv)) set_tau = set_tau_conv + set_tau = set_tau .and. s% superad_reduction_use_turnover_limit + use_mlt_vc_tau_surface = .false. + if (set_tau .and. s% have_mlt_vc .and. nzlo == 1) then + if (associated(s% mlt_vc)) & + use_mlt_vc_tau_surface = s% mlt_vc(1) > 0d0 + end if if (s% doing_timing) call start_time(s, time0, total) !$OMP PARALLEL DO PRIVATE(k,op_err,make_gradr_sticky_in_solver_iters) SCHEDULE(dynamic,2) do k = nzlo, nzhi op_err = 0 - call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, op_err) + call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, op_err, & + set_tau_conv_in = set_tau) if (make_gradr_sticky_in_solver_iters .and. s% solver_iter > 3) then if (.not. s% fixed_gradr_for_rest_of_solver_iters(k)) then s% fixed_gradr_for_rest_of_solver_iters(k) = & @@ -63,6 +76,10 @@ subroutine set_mlt_vars(s, nzlo, nzhi, ierr) if (op_err /= 0) ierr = op_err end do !$OMP END PARALLEL DO + ! The surface Brunt timescale uses the first interior face. + if (ierr == 0 .and. set_tau .and. s% use_face_reconstruction .and. & + nzlo == 1 .and. nzhi >= 2 .and. .not. use_mlt_vc_tau_surface) & + s% tau_conv_start(1) = s% tau_conv_start(2) if (s% doing_timing) call update_time(s, time0, total, s% time_mlt) end subroutine set_mlt_vars @@ -70,7 +87,7 @@ end subroutine set_mlt_vars subroutine do1_mlt_2(s, k, & make_gradr_sticky_in_solver_iters, ierr, & - mixing_length_alpha_in, gradL_composition_term_in) + mixing_length_alpha_in, gradL_composition_term_in, set_tau_conv_in) ! get convection info for point k use star_utils use turb_support, only: do1_mlt_eval @@ -82,14 +99,15 @@ subroutine do1_mlt_2(s, k, & integer, intent(out) :: ierr real(dp), intent(in), optional :: & mixing_length_alpha_in, gradL_composition_term_in + logical, intent(in), optional :: set_tau_conv_in type(auto_diff_real_star_order1) :: gradr_factor - real(dp) :: f, gradL_composition_term, abs_du_div_cs, cs, mixing_length_alpha + real(dp) :: f, gradL_composition_term, abs_du_div_cs, cs, mixing_length_alpha, tau_conv real(dp), pointer :: vel(:) integer :: i, mixing_type, nz, k_T_max real(dp), parameter :: conv_vel_mach_limit = 0.9d0 real(dp) :: crystal_pad - logical :: no_mix + logical :: no_mix, set_tau_conv type(auto_diff_real_star_order1) :: & T_face_ad, P_face_ad, energy_face_ad, opacity_face_ad, rho_face_ad, chiRho_face_ad, chiT_face_ad, Cp_face_ad, & grada_face_ad, scale_height_ad, gradr_ad, & @@ -98,6 +116,8 @@ subroutine do1_mlt_2(s, k, & ierr = 0 nz = s% nz + set_tau_conv = .false. + if (present(set_tau_conv_in)) set_tau_conv = set_tau_conv_in if (k < 1 .or. k > nz) then write(*,3) 'bad k for do1_mlt', k, nz @@ -129,6 +149,21 @@ subroutine do1_mlt_2(s, k, & grada_face_ad, opacity_face_ad, scale_height_ad, gradr_ad, ierr) if (ierr /= 0) return + ! Save the timescale after face reconstruction and before MLT updates mlt_vc. + if (set_tau_conv) then + tau_conv = s% tau_conv_start(k) + if (s% use_face_reconstruction .and. k > 1) & + tau_conv = conv_time_scale(s,k) + if (s% have_mlt_vc) then + if (associated(s% mlt_vc)) then + if (s% mlt_vc(k) > 0d0) & + tau_conv = scale_height_ad% val/ & + max(s% mlt_vc(k), min_mlt_vc_for_tau_conv) + end if + end if + s% tau_conv_start(k) = tau_conv + end if + if (s% rotation_flag .and. s% mlt_use_rotation_correction) then gradr_factor = s% ft_rot(k)/s% fp_rot(k)*s% gradr_factor(k) else diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 3e6b19357e..6f0c3aab92 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -20,7 +20,7 @@ module turb_support use star_private_def -use const_def, only: dp, crad, no_mixing +use const_def, only: dp, crad, ln10, no_mixing use num_lib use utils_lib use auto_diff_support @@ -228,9 +228,10 @@ 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 - type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, & + real(dp) :: f_turnover, Gamma_factor_old_local, eta_old + type(auto_diff_real_star_order1) :: eta_inst, 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 + logical :: test_partials, using_TDC, have_Y_face_guess, hold_superad_reduction_factor logical, parameter :: report = .false. include 'formats' @@ -298,7 +299,16 @@ 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 + hold_superad_reduction_factor = .false. + ! Preserve stored photo state while finish_load_model rebuilds derived quantities. + if (k > 0 .and. s% have_superad_reduction_factor) & + hold_superad_reduction_factor = & + s% doing_finish_load_model .or. & + (s% use_superad_reduction .and. & + s% superad_reduction_use_turnover_limit .and. & + .not. s% okay_to_set_superad_reduction_factor) + if (k /= 0 .and. .not. hold_superad_reduction_factor) & + s% superad_reduction_factor(k) = 1d0 ! 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. & @@ -468,6 +478,11 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg contains subroutine set_superad_reduction() + Gamma_factor = 1d0 + if (k > 0) then + if (s% lnT_start(k)/ln10 >= s% superad_reduction_max_logT) return + end if + Gamma_limit = s% superad_reduction_Gamma_limit scale_value1 = s% superad_reduction_Gamma_limit_scale scale_value2 = s% superad_reduction_Gamma_inv_scale @@ -476,7 +491,6 @@ subroutine set_superad_reduction() Lrad_div_Ledd = 4d0*crad/3d0*pow4(T)/P*gradT Gamma_inv_threshold = 4d0*(1d0-beta)/(4d0-3*beta) - Gamma_factor = 1d0 if (gradT > gradL) then if (Lrad_div_Ledd > Gamma_limit .or. Lrad_div_Ledd > Gamma_inv_threshold) then alfa0 = (gradT-gradL)/diff_grads_limit @@ -505,9 +519,9 @@ subroutine set_superad_reduction() if (Lrad_div_Ledd% val > Gamma_inv_threshold) then 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 @@ -523,6 +537,30 @@ subroutine set_superad_reduction() end if end if end if + + if (hold_superad_reduction_factor) then + ! Start-of-step setup must not advance the stored response. + Gamma_factor = max(s% superad_reduction_factor(k), 1d0) + else if (s% superad_reduction_use_turnover_limit .and. k > 0) then + if (s% dt > 0d0 .and. s% tau_conv_start(k) > 0d0 .and. & + s% have_superad_reduction_factor) then + if (associated(s% superad_reduction_factor_old)) then + Gamma_factor_old_local = max(s% superad_reduction_factor_old(k), 1d0) + if (Gamma_factor > 1d0 .or. Gamma_factor_old_local > 1d0) then + select case (trim(s% superad_reduction_turnover_limit_function)) + case ('linear') + f_turnover = min(s% dt/s% tau_conv_start(k), 1d0) + case ('exponential') + f_turnover = -expm1(-s% dt/s% tau_conv_start(k)) + end select + eta_old = 1d0/Gamma_factor_old_local + eta_inst = 1d0/Gamma_factor + Gamma_factor = 1d0/(eta_old + f_turnover*(eta_inst - eta_old)) + end if + end if + 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/test_suite/zams_to_cc_80/inlist_common b/star/test_suite/zams_to_cc_80/inlist_common index d2056e06f9..681699e07e 100644 --- a/star/test_suite/zams_to_cc_80/inlist_common +++ b/star/test_suite/zams_to_cc_80/inlist_common @@ -120,23 +120,18 @@ -! superadiabatic convection routines, it's a choice: -! typically only necessary at higher metallicities 1/5 solar + +! superadiabatic reduction - ! superadiabatic reduction, implicit, new - use_superad_reduction = .false. + use_superad_reduction = .true. superad_reduction_Gamma_limit = 0.5d0 !default is 0.5d0 superad_reduction_Gamma_limit_scale = 5d0 superad_reduction_Gamma_inv_scale = 5d0 superad_reduction_diff_grads_limit = 1d-2 ! default is 1d-3 superad_reduction_limit = -1d0 + superad_reduction_use_turnover_limit = .true. + superad_reduction_turnover_limit_function = 'linear' - - ! MLT ++, explicit, well tested - okay_to_reduce_gradT_excess = .true. - gradT_excess_f1 = 1d-4 - gradT_excess_f2 = 1d-2 - !gradT_excess_lambda1 = -1d0 ! full on + okay_to_reduce_gradT_excess = .false. ! mixing D_omega_mixing_rate = 1d0 diff --git a/star/test_suite/zams_to_cc_80/inlist_to_zams b/star/test_suite/zams_to_cc_80/inlist_to_zams index 7960b4d2dc..26ed9f20d6 100644 --- a/star/test_suite/zams_to_cc_80/inlist_to_zams +++ b/star/test_suite/zams_to_cc_80/inlist_to_zams @@ -51,23 +51,6 @@ ! mlt - ! turn off mlt++ on the pre-zams. - - ! superadiabatic reduction, implicit, new - use_superad_reduction = .false. - superad_reduction_Gamma_limit = 0.5d0 !default is 0.5d0 - superad_reduction_Gamma_limit_scale = 5d0 - superad_reduction_Gamma_inv_scale = 5d0 - superad_reduction_diff_grads_limit = 1d-2 ! default is 1d-3 - superad_reduction_limit = -1d0 - - - ! MLT ++, explicit, well tested - okay_to_reduce_gradT_excess = .false. - gradT_excess_f1 = 1d-4 - gradT_excess_f2 = 1d-2 - !gradT_excess_lambda1 = -1d0 ! full on - ! mixing ! timesteps diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index 2d93c0483c..de0ded144a 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -301,11 +301,14 @@ logical :: use_superad_reduction real(dp) :: & + superad_reduction_max_logT, & superad_reduction_Gamma_limit, & superad_reduction_Gamma_limit_scale, & 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 ! mixing parameters diff --git a/star_data/public/star_data_def.inc b/star_data/public/star_data_def.inc index 90c31b5bf4..4bbfdce821 100644 --- a/star_data/public/star_data_def.inc +++ b/star_data/public/star_data_def.inc @@ -1,7 +1,7 @@ character(len=24) :: version_number ! mesa version from file $MESA_DIR/data/version_number - integer, parameter :: star_def_version = 18 + integer, parameter :: star_def_version = 19 integer, parameter :: nz_alloc_extra = 200 diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index fa97ddb85f..4f9f176c03 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 @@ -507,6 +508,7 @@ ! this is the Brassard et al B := -(chiY/chiT)*(dlnY/dlnP) real(dp), pointer :: unsmoothed_brunt_B(:) ! pre-smoothing real(dp), pointer :: smoothed_brunt_B(:) ! with smoothing applied for gradL + real(dp), pointer :: tau_conv_start(:) ! seconds; fixed during solver iterations ! asteroseismology info @@ -564,7 +566,9 @@ 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 + ! previous accepted Gamma_factor for the turnover limiter 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 @@ -901,7 +905,9 @@ 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, & + okay_to_set_superad_reduction_factor, have_superad_reduction_factor + ! current factor is frozen until it has been copied to the old state real(dp), pointer :: xa_removed(:) ! mass fractions for removed mass (1:species) real(dp) :: h1_czb_mass