From ae4228d12a00750f435f6d3a310db6a1de0a6558 Mon Sep 17 00:00:00 2001 From: Philip Mocz Date: Fri, 9 Aug 2024 18:54:23 -0400 Subject: [PATCH] starspots Part 2 (#713) * [ci skip] adding some stuff, not yet done * [ci skip] updates to starspots * [ci skip] run_stars_extras * [ci skip] fix controls formatting * [ci skip] remove redundant library * [ci skip] modify turb_support * [ci skip] switch run_stars_extra to plain again * [ci skip] move need_atm_Psurf to right place * [ci skip] mv starspot stuff in hydro_vars.f90 * oh I think I found it --------- Co-authored-by: Meridith Joyce --- star/defaults/controls.defaults | 26 +- star/make/makefile_base | 1 + star/private/ctrls_io.f90 | 4 +- star/private/hydro_vars.f90 | 20 +- star/private/starspots.f90 | 111 +++++++ star/private/turb_support.f90 | 17 +- star/test_suite/starspots/README.rst | 2 +- star/test_suite/starspots/inlist_starspots | 6 +- .../starspots/src/run_star_extras.f90 | 299 ++---------------- star_data/private/star_controls.inc | 1 + 10 files changed, 199 insertions(+), 288 deletions(-) create mode 100644 star/private/starspots.f90 diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index eb00f74a0..68077fa56 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -4640,20 +4640,34 @@ ! ========= - ! fspot - ! ~~~~~~~~~~~ + ! do_starspots + ! ~~~~~~~~~~~~ + + ! If ``.true.``, switch on impedence of the surface flux due to magnetic pressure from starspots, + ! parameterized in the style of an atmospheric boundary modification. First described by + ! `Somers et al. (2015; ApJ) `__. + ! Detailed discussion of this functionality can be found in |MESA VI|. + + ! :: - ! Filling Factor of starspots (in [0,1]) + do_starspots = .false. + + ! fspot + ! ~~~~~ + + ! Filling Factor of starspots. Valid values between 0.0 and 1.0 (no spots to 100% coverage) ! :: fspot = 0d0 ! xspot - ! ~~~~~~~~~~~ + ! ~~~~~ - ! Temperature contrast between the spotted and unspotted regions (in [0,1]) + ! Temperature contrast between the spotted and unspotted regions. + ! Valid values are between 1.0 (no contribution from magnetic pressure) and 0.5 + ! (half of the total pressure is due to magnetic pressure) ! :: @@ -8291,7 +8305,7 @@ ! Note that these parameters are not exactly the same ! as used by Paul Duffell. - ! His calibrated D is 2, where mesa has default D = 3 (see mesaIV paper). + ! His calibrated D is 2, where mesa has default D = 3 (see |MESA IV|). ! Users should try various values since the choice is not clear cut. ! :: diff --git a/star/make/makefile_base b/star/make/makefile_base index 9568ed966..4afdafa85 100644 --- a/star/make/makefile_base +++ b/star/make/makefile_base @@ -205,6 +205,7 @@ SRCS = \ star_solver.f90 \ struct_burn_mix.f90 \ winds.f90 \ + starspots.f90 \ gravity_darkening.f90 \ mass_utils.f90 \ eps_mdot.f90 \ diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 251379874..6a277c844 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -492,7 +492,7 @@ module ctrls_io use_T_tau_gradr_factor, & ! starspots - fspot, xspot, & + do_starspots, fspot, xspot, & ! extra heat near surface to model irradiation irradiation_flux, column_depth_for_irradiation, & @@ -1289,6 +1289,7 @@ subroutine store_controls(s, ierr) s% use_T_tau_gradr_factor = use_T_tau_gradr_factor ! starspots + s% do_starspots = do_starspots s% fspot = fspot s% xspot = xspot @@ -2974,6 +2975,7 @@ subroutine set_controls_for_writing(s, ierr) use_T_tau_gradr_factor = s% use_T_tau_gradr_factor ! starspots + do_starspots = s% do_starspots fspot = s% fspot xspot = s% xspot diff --git a/star/private/hydro_vars.f90 b/star/private/hydro_vars.f90 index bb1a2fab7..e190885cb 100644 --- a/star/private/hydro_vars.f90 +++ b/star/private/hydro_vars.f90 @@ -387,22 +387,27 @@ end subroutine set_Teff subroutine set_Teff_info_for_eqns(s, skip_partials, & - need_atm_Psurf, need_atm_Tsurf, r_surf, L_surf, Teff, & + need_atm_Psurf_in, need_atm_Tsurf_in, r_surf, L_surf, Teff, & lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, & ierr) use star_utils, only: set_phot_info use atm_lib, only: atm_Teff + use starspots, only: starspot_tweak_PT, starspot_restore_PT type (star_info), pointer :: s logical, intent(in) :: skip_partials, & - need_atm_Psurf, need_atm_Tsurf + need_atm_Psurf_in, need_atm_Tsurf_in real(dp), intent(out) :: r_surf, L_surf, Teff, & lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap integer, intent(out) :: ierr + logical :: need_atm_Psurf, need_atm_Tsurf include 'formats' + need_atm_Psurf = need_atm_Psurf_in + need_atm_Tsurf = need_atm_Tsurf_in + ierr = 0 r_surf = s% r(1) @@ -435,11 +440,21 @@ subroutine set_Teff_info_for_eqns(s, skip_partials, & lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, & ierr) else + ! starspot YREC routine + if (s% do_starspots) then + need_atm_Psurf = .true. + need_atm_Tsurf = .true. + call starspot_tweak_PT(s) + end if call get_surf_PT( & s, skip_partials, need_atm_Psurf, need_atm_Tsurf, & lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, & ierr) + ! starspot YREC routine + if (s% do_starspots) then + call starspot_restore_PT(s) + end if end if if (ierr /= 0) then if (s% report_ierr) then @@ -785,7 +800,6 @@ subroutine get_surf_PT( & real(dp) :: kap_surf real(dp) :: M_surf - include 'formats' ! Set up stellar surface parameters diff --git a/star/private/starspots.f90 b/star/private/starspots.f90 new file mode 100644 index 000000000..da13b27b9 --- /dev/null +++ b/star/private/starspots.f90 @@ -0,0 +1,111 @@ +! *********************************************************************** +! +! Copyright (C) 2010-2019 Meridith Joyce & The MESA Team +! +! MESA is free software; you can use it and/or modify +! it under the combined terms and restrictions of the MESA MANIFESTO +! and the GNU General Library Public License as published +! by the Free Software Foundation; either version 2 of the License, +! or (at your option) any later version. +! +! You should have received a copy of the MESA MANIFESTO along with +! this software; if not, it is available at the mesa website: +! http://mesa.sourceforge.net/ +! +! MESA is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +! See the GNU Library General Public License for more details. +! +! You should have received a copy of the GNU Library General Public License +! along with this software; if not, write to the Free Software +! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +! +! *********************************************************************** + +module starspots + + use star_private_def + use const_def + use utils_lib + + implicit none + + private + + public :: starspot_tweak_gradr + public :: starspot_tweak_PT + public :: starspot_restore_PT + + real(dp) :: L_init + +contains + +! ------------- +! parameterized YREC routines +! MESA models a pressure contrast rather than temperature contrast +! ------------- + + subroutine starspot_tweak_gradr(s, P, gradr, gradr_spot) + ! adjusts the gradient of the radius to account for starspots. + ! This subroutine is called at the beginning of Get_results() + ! in turb_support.f90 + ! ------------------------------------------------------------ + use auto_diff_support + type(star_info), pointer :: s + type(auto_diff_real_star_order1), intent(in) :: P + type(auto_diff_real_star_order1), intent(in) :: gradr + type(auto_diff_real_star_order1), intent(out) :: gradr_spot + + real(dp) :: mu_ideal_gas, R2, Teff_local, PB_i + type(auto_diff_real_star_order1) :: xspot_of_r ! xspot4 + + if (.not. s%doing_relax .and. .not. s%doing_first_model_of_run) then + mu_ideal_gas = s%mu(1) !1.00794d0 ! for hydrogen, 1 gram per mole + R2 = pow2(s%R(1)) + Teff_local = pow(s%L(1)/(pi4*boltz_sigma*R2), 0.25d0) + PB_i = (cgas*s%rho(1)/mu_ideal_gas)*(1.0 - s%xspot)*Teff_local + xspot_of_r = (P - PB_i)/P + gradr_spot = gradr/(s%fspot*pow4(xspot_of_r) + 1d0 - s%fspot) + else + gradr_spot = gradr + end if + end subroutine starspot_tweak_gradr + + subroutine starspot_tweak_PT(s) + ! saves the surface luminosity and adjusts it and the effective + ! temperature to account for starspots. + ! This subroutine is called at the beginning of get_surf_PT() + ! in hydro_vars.f90 + ! ------------------------------------------------------------ + + type(star_info), pointer :: s + + real(dp) :: alp + + alp = 1d0 - s%fspot + s%fspot*pow4(s%xspot) + + ! This is the surface-average value for luminosity + L_init = s%L(1) + + ! Set the surface L to the unspotted, ambient L + s%L(1) = s%L(1)/alp + + ! Now, set the Teff. Used in atm table lookup to set boundary conditions + s%Teff = pow(s%L(1)/(pi4*pow2(s%r(1))*boltz_sigma), 0.25_dp) + + end subroutine starspot_tweak_PT + + subroutine starspot_restore_PT(s) + ! restores the surface luminosity effective temeperature. + ! Called at the end of get_surf_PT() + ! ------------------------------------------------------------ + + type(star_info), pointer :: s + + s%Teff = pow(L_init/(pi4*pow2(s%r(1))*boltz_sigma), 0.25_dp) + s%L(1) = L_init + + end subroutine starspot_restore_PT + +end module starspots diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 08ec63107..6fa8326d8 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -105,13 +105,14 @@ end subroutine get_gradT subroutine do1_mlt_eval( & s, k, MLT_option, gradL_composition_term, & - gradr, grada, scale_height, mixing_length_alpha, & + gradr_in, grada, scale_height, mixing_length_alpha, & mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) use chem_def, only: ih1 + use starspots, only: starspot_tweak_gradr type (star_info), pointer :: s integer, intent(in) :: k character (len=*), intent(in) :: MLT_option - type(auto_diff_real_star_order1), intent(in) :: gradr, grada, scale_height + type(auto_diff_real_star_order1), intent(in) :: gradr_in, grada, scale_height real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha integer, intent(out) :: mixing_type type(auto_diff_real_star_order1), intent(out) :: & @@ -120,9 +121,11 @@ subroutine do1_mlt_eval( & real(dp) :: cgrav, m, XH1, gradL_old, grada_face_old integer :: iso, old_mix_type - type(auto_diff_real_star_order1) :: r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp + type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp include 'formats' ierr = 0 + + gradr = gradr_in cgrav = s% cgrav(k) m = s% m_grav(k) @@ -145,7 +148,12 @@ subroutine do1_mlt_eval( & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) - else + else + ! starspot YREC routine + if (s% do_starspots) then + dV = 0d0 ! dV = 1/rho - 1/rho_start and we assume rho = rho_start. + call starspot_tweak_gradr(s, P, gradr_in, gradr) + end if call Get_results(s, k, MLT_option, & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & @@ -228,7 +236,6 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg end if - ! check if this particular k can be done with TDC using_TDC = .false. if (s% MLT_option == 'TDC') using_TDC = .true. diff --git a/star/test_suite/starspots/README.rst b/star/test_suite/starspots/README.rst index f7b60e9c1..f80abafd3 100644 --- a/star/test_suite/starspots/README.rst +++ b/star/test_suite/starspots/README.rst @@ -34,6 +34,6 @@ A plot of the HR diagram for the problem is included below: :width: 100% -Last-Updated: 06Aug2024 by Meridith Joyce +Last-Updated: 09Aug2024 by Meridith Joyce Last-Run: 06Aug2024 (MESA 7890d305) by pmocz on C916PXT6XW in 629 seconds using 8 threads. diff --git a/star/test_suite/starspots/inlist_starspots b/star/test_suite/starspots/inlist_starspots index ef54ff85d..df1f3213b 100644 --- a/star/test_suite/starspots/inlist_starspots +++ b/star/test_suite/starspots/inlist_starspots @@ -61,10 +61,8 @@ star_history_name = 'history.data' max_age = 24.63d9 !26d9 !0.8d9 - - use_other_mlt_results = .true. - use_other_surface_PT = .true. - + + do_starspots = .true. fspot = 0.34d0 xspot = 0.85d0 diff --git a/star/test_suite/starspots/src/run_star_extras.f90 b/star/test_suite/starspots/src/run_star_extras.f90 index 9145d2541..89743ce76 100644 --- a/star/test_suite/starspots/src/run_star_extras.f90 +++ b/star/test_suite/starspots/src/run_star_extras.f90 @@ -1,6 +1,6 @@ ! *********************************************************************** ! -! Copyright (C) 2010-2019 Bill Paxton & The MESA Team +! Copyright (C) 2011 The MESA Team ! ! this file is part of mesa. ! @@ -26,21 +26,18 @@ module run_star_extras use star_def use const_def use math_lib - use chem_def !! maybe taking up uncessary space but whatever - - + use auto_diff + implicit none - real(dp) :: fspot, xspot, PB_i, Teff_local include "test_suite_extras_def.inc" - + + ! these routines are called by the standard run_star check_model contains include "test_suite_extras.inc" - - - ! these routines are called by the standard run_star check_model - + + subroutine extras_controls(id, ierr) integer, intent(in) :: id integer, intent(out) :: ierr @@ -48,40 +45,17 @@ subroutine extras_controls(id, ierr) ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - - ! this is the place to set any procedure pointers you want to change - ! e.g., other_wind, other_mixing, other_energy (see star_data.inc) - - ! the extras functions in this file will not be called - ! unless you set their function pointers as done below. - ! otherwise we use a null_ version which does nothing (except warn). - - - fspot = s% fspot - xspot = s% xspot - - s% extras_startup => extras_startup - s% extras_start_step => extras_start_step s% extras_check_model => extras_check_model s% extras_finish_step => extras_finish_step s% extras_after_evolve => extras_after_evolve s% how_many_extra_history_columns => how_many_extra_history_columns s% data_for_extra_history_columns => data_for_extra_history_columns s% how_many_extra_profile_columns => how_many_extra_profile_columns - s% data_for_extra_profile_columns => data_for_extra_profile_columns - - s% how_many_extra_history_header_items => how_many_extra_history_header_items - s% data_for_extra_history_header_items => data_for_extra_history_header_items - s% how_many_extra_profile_header_items => how_many_extra_profile_header_items - s% data_for_extra_profile_header_items => data_for_extra_profile_header_items - - s% other_surface_PT => starspot_tweak_PT - s% other_mlt_results => YREC_spots_other_mlt_results - + s% data_for_extra_profile_columns => data_for_extra_profile_columns end subroutine extras_controls - - + + subroutine extras_startup(id, restart, ierr) integer, intent(in) :: id logical, intent(in) :: restart @@ -92,153 +66,19 @@ subroutine extras_startup(id, restart, ierr) if (ierr /= 0) return call test_suite_startup(s, restart, ierr) end subroutine extras_startup - - - integer function extras_start_step(id) + + + subroutine extras_after_evolve(id, ierr) integer, intent(in) :: id - integer :: ierr - real(dp) :: mu_ideal_gas, R2 + integer, intent(out) :: ierr type (star_info), pointer :: s - real(dp) :: power_he_burn, power_c_burn, power_neutrinos, & - center_h1, center_he4, ocz_top_mass, ocz_bot_mass, & - ocz_top_radius, ocz_bot_radius!, mass_difference_prev !! no!! - integer :: nz, j, i, k, k_ocz_top, k_ocz_bot, n_conv_bdy -! include 'formats' + real(dp) :: dt ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - extras_start_step = 0 - - !! to use constants as defined in MESA, import const_def and call - !! variables by their names there - - ! set PB_i here and then other_mlt will know about it, because - ! quantities set here are carried over the course of a timestep, - ! AKA all newton iterations - - mu_ideal_gas = s% mu(1) !1.00794d0 ! for hydrogen, 1 gram per mole - !write(*,*) 'MESA def, my def: ', s% mu(1), mu_ideal_gas - - R2 = pow2(s%R(1)) - Teff_local = pow(s%L(1)/(pi4*boltz_sigma*R2), 0.25d0) - PB_i = (cgas* s%rho(1)/mu_ideal_gas) * (1.0 - xspot) * Teff_local - - end function extras_start_step - - -!----------------------------------------------------------------------------------------- -! -! YREC routines -! -!---------------------------------------------------------------------------------------- - subroutine YREC_spots_other_mlt_results(id, k, MLT_option, & ! NOTE: k=0 is a valid arg - r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, & - iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & - alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) - use const_def, only: dp - use auto_diff - integer, intent(in) :: id - integer, intent(in) :: k - ! integer, intent(out) :: ierr - ! type (star_info), pointer :: s - ! integer :: nz, j - character (len=*), intent(in) :: MLT_option - type(auto_diff_real_star_order1), intent(in) :: & - r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height - integer, intent(in) :: iso - real(dp), intent(in) :: & - XH1, cgrav, m, gradL_composition_term, & - mixing_length_alpha, alpha_semiconvection, thermohaline_coeff - integer, intent(out) :: mixing_type - type(auto_diff_real_star_order1), intent(out) :: & - gradT, Y_face, conv_vel, D, Gamma - integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: xspot_of_r !, xspot4 - type(auto_diff_real_star_order1) :: gradr_spot - !ierr = 0 - ! ------------------------------ 10/26/21 - type (star_info), pointer :: s - ierr = 0 - call star_ptr(id, s, ierr) - if (ierr /= 0) return - - !------------------------------ - !if (s% star_age >= 10d0) then - if (.not. s% doing_relax .and. .not. s% doing_first_model_of_run) then - xspot_of_r = (P - PB_i)/P - gradr_spot = gradr/( fspot*pow4(xspot_of_r) + 1d0 - fspot) - else - gradr_spot = gradr - end if - - call star_mlt_results(id, k, MLT_option, & - r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr_spot, grada, scale_height, & - iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & - alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) - end subroutine YREC_spots_other_mlt_results - - - subroutine starspot_tweak_PT(id, skip_partials, & - lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & - lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, ierr) - - use const_def, only: dp - - integer, intent(in) :: id - logical, intent(in) :: skip_partials - - logical :: need_atm_Psurf, need_atm_Tsurf - - real(dp), intent(out) :: & - lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & - lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap - integer, intent(out) :: ierr - - ! For my tweaks - real(dp) :: alp - - ! Call the stock get_surf_PT - type (star_info), pointer :: s - real(dp) :: L_init - - include 'formats' - - ierr = 0 - call star_ptr(id, s, ierr) - if (ierr /= 0) return - - need_atm_Psurf = .true. - need_atm_Tsurf = .true. - - alp = 1d0 - fspot + fspot*pow4(xspot) - - ! This is the surface-average value for luminosity - L_init = s% L(1) - - ! Set the surface L to the unspotted, ambient L - s% L(1) = s% L(1) / alp - - ! Now, set the Teff. Used in atm table lookup to set boundary conditions - s% Teff = pow(s% L(1)/(pi4*pow2(s% r(1))*boltz_sigma), 0.25_dp) - - ! Set everything with Lamb. - call star_get_surf_PT(id, skip_partials, need_atm_Psurf, need_atm_Tsurf, & - lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, & - lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, & - ierr) - - s% Teff = pow(L_init/(pi4*pow2(s% r(1))*boltz_sigma), 0.25_dp) - s% L(1) = L_init - - end subroutine starspot_tweak_PT - -!------------------------------------------------------------------------------ -! -! end YREC routines -! -!------------------------------------------------------------------------------ + call test_suite_after_evolve(s, ierr) + end subroutine extras_after_evolve + ! returns either keep_going, retry, or terminate. integer function extras_check_model(id) @@ -248,17 +88,7 @@ integer function extras_check_model(id) ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - extras_check_model = keep_going - - if (.false. .and. s% star_mass_h1 < 0.35d0) then - ! stop when star hydrogen mass drops to specified level - extras_check_model = terminate - write(*, *) 'have reached desired hydrogen mass' - return - end if - - ! by default, indicate where (in the code) MESA terminated - if (extras_check_model == terminate) s% termination_code = t_extras_check_model + extras_check_model = keep_going end function extras_check_model @@ -269,11 +99,10 @@ integer function how_many_extra_history_columns(id) ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - how_many_extra_history_columns = 0 end function how_many_extra_history_columns - - + + subroutine data_for_extra_history_columns(id, n, names, vals, ierr) integer, intent(in) :: id, n character (len=maxlen_history_column_name) :: names(n) @@ -283,15 +112,11 @@ subroutine data_for_extra_history_columns(id, n, names, vals, ierr) ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - - ! note: do NOT add the extras names to history_columns.list - ! the history_columns.list is only for the built-in history column options. - ! it must not include the new column names you are adding here. - end subroutine data_for_extra_history_columns - + integer function how_many_extra_profile_columns(id) + use star_def, only: star_info integer, intent(in) :: id integer :: ierr type (star_info), pointer :: s @@ -300,9 +125,11 @@ integer function how_many_extra_profile_columns(id) if (ierr /= 0) return how_many_extra_profile_columns = 0 end function how_many_extra_profile_columns - - + + subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr) + use star_def, only: star_info, maxlen_profile_column_name + use const_def, only: dp integer, intent(in) :: id, n, nz character (len=maxlen_profile_column_name) :: names(n) real(dp) :: vals(nz,n) @@ -313,56 +140,9 @@ subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr) call star_ptr(id, s, ierr) if (ierr /= 0) return end subroutine data_for_extra_profile_columns - - - integer function how_many_extra_history_header_items(id) - integer, intent(in) :: id - integer :: ierr - type (star_info), pointer :: s - ierr = 0 - call star_ptr(id, s, ierr) - if (ierr /= 0) return - how_many_extra_history_header_items = 0 - end function how_many_extra_history_header_items - - - subroutine data_for_extra_history_header_items(id, n, names, vals, ierr) - integer, intent(in) :: id, n - character (len=maxlen_history_column_name) :: names(n) - real(dp) :: vals(n) - type(star_info), pointer :: s - integer, intent(out) :: ierr - ierr = 0 - call star_ptr(id,s,ierr) - if(ierr/=0) return - end subroutine data_for_extra_history_header_items - - - integer function how_many_extra_profile_header_items(id) - integer, intent(in) :: id - integer :: ierr - type (star_info), pointer :: s - ierr = 0 - call star_ptr(id, s, ierr) - if (ierr /= 0) return - how_many_extra_profile_header_items = 0 - end function how_many_extra_profile_header_items - - - subroutine data_for_extra_profile_header_items(id, n, names, vals, ierr) - integer, intent(in) :: id, n - character (len=maxlen_profile_column_name) :: names(n) - real(dp) :: vals(n) - type(star_info), pointer :: s - integer, intent(out) :: ierr - ierr = 0 - call star_ptr(id,s,ierr) - if(ierr/=0) return - end subroutine data_for_extra_profile_header_items - + ! returns either keep_going or terminate. - ! note: cannot request retry; extras_check_model can do that. integer function extras_finish_step(id) integer, intent(in) :: id integer :: ierr @@ -371,26 +151,9 @@ integer function extras_finish_step(id) call star_ptr(id, s, ierr) if (ierr /= 0) return extras_finish_step = keep_going - - ! to save a profile, - ! s% need_to_save_profiles_now = .true. - ! to update the star log, - ! s% need_to_update_history_now = .true. - - ! see extras_check_model for information about custom termination codes - ! by default, indicate where (in the code) MESA terminated - if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step end function extras_finish_step - - - subroutine extras_after_evolve(id, ierr) - integer, intent(in) :: id - integer, intent(out) :: ierr - type (star_info), pointer :: s - ierr = 0 - call star_ptr(id, s, ierr) - if (ierr /= 0) return - call test_suite_after_evolve(s, ierr) - end subroutine extras_after_evolve + + end module run_star_extras + diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index 64233afca..855cc6a5f 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -433,6 +433,7 @@ logical :: use_T_tau_gradr_factor ! starspots + logical :: do_starspots real(dp) :: fspot, xspot ! mass gain or loss