## Learning outcomes

• Learn how to write your own run_stars_extra.f to perform calculations not already in MESA
• Pass arguments through the inlist to custom routines
• Add custom defined variables to history log file

NOTE: Global variables can be found in $MESA_DIR/const/public/const_def.f90. The default run_stars_extras.f file, which can be found in /src from a MESA work directory, looks as follows:  module run_star_extras use star_lib use star_def use const_def use crlibm_lib implicit none ! these routines are called by the standard run_star check_model contains include 'standard_run_star_extras.inc' end module run_star_extras  As can be seen, all of the heavy lifting is found in standard_run_star_extras.inc. This file can be found in $MESA_DIR/include/standard_run_star_extras.inc. To create your own run_stars_extra.f routines, you will need to first paste the contents of that file here.

The execution of MESA proceeds as follows:

(Courtesy Josiah Schwab, MESA Summer School 2016)

All of the boxes to the right are places in run_stars_extra.f where control can be taken. The gray take_step box is where all of the heavy machinery of MESA takes place.

## An example extension

We will now modify MESA to calculate the blackbody temperature of the Earth. We will further make MESA stop if it exceeds a certain value. That temperature can be calculated with

$$\frac{T_\oplus}{T_\odot} = \sqrt{\frac{R_\odot}{2\text{AU}}}.$$

In order to insert a stopping condition, we must modify extras_finish_step. We need to declare the variable we want to use (Tearth) and then return a command to terminate if our condition is met. We furthermore want the user to be able to specify the stopping value. We can do that as follows:

      ! returns either keep_going or terminate.
! note: cannot request retry or backup; extras_check_model can do that.
integer function extras_finish_step(id, id_extra)
integer, intent(in) :: id, id_extra
integer :: ierr
type (star_info), pointer :: s

real(dp) :: Tearth

ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
extras_finish_step = keep_going
call store_extra_info(s)

! calculate blackbody temperature of earth
Tearth = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU))
write(*,*) "Tearth =", Tearth

! stop if it exceeds 300 K
if (Tearth > s% x_ctrl(1)) extras_finish_step = terminate

! 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


NOTE 1: Structure variables in Fortran can be accessed with %, e.g. s% Teff accesses the effective temperature of variable s. This is analogous to s.Teff in C-like languages.

NOTE 2: Constants (i.e. values defined in $MESA_DIR/const) do not need s% to be accessed; e.g. Rsun can be accessed directly. NOTE 3: The comments before the header of the function call says the possible return values: keep_going and terminate. NOTE 4: x_ctrl is a vector of real numbers that the user can set in the inlist. In this case, we have designated x_ctrl(1) to denote the maximum blackbody temperature of the Earth. One can use as many as desired. Also available are x_integer_ctrl and x_logical_ctrl which are vectors for integers and booleans, respectively. Now we can modify our inlist to have, for example, x_ctrl(1) = 300 ! K; Tearth cut-off value  inside of the &controls section. That's pretty good, but if we were to ./mk and ./rn this, we would obtain a model whose temperature exceeds the value we set. Instead, it might be better to iterate until we get one exactly on the value that we want. ## Iteratively finding a specified value For this, we can modify extras_check_model and retry until we get the exact value -- or at least one within a user-specified tolerance:  ! returns either keep_going, retry, backup, or terminate. integer function extras_check_model(id, id_extra) integer, intent(in) :: id, id_extra integer :: ierr type (star_info), pointer :: s real(dp) :: Tearth, Tearth_cutoff, Tearth_tol 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 Tearth_cutoff = s% x_ctrl(1) Tearth_tol = s% x_ctrl(2) Tearth = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU)) if (Tearth < Tearth_cutoff - Tearth_tol) then extras_check_model = keep_going write(*, *) 'Tearth too low, keep going' end if if (Tearth > Tearth_cutoff + Tearth_tol) then extras_check_model = retry s% dt = s% dt / 10 write(*, *) 'Tearth exceeds limit, backing up' end if if (Tearth >= Tearth_cutoff - Tearth_tol .and. & Tearth <= Tearth_cutoff + Tearth_tol) then extras_check_model = terminate write(*, *) 'Bullseye!' end if ! if you want to check multiple conditions, it can be useful ! to set a different termination code depending on which ! condition was triggered. MESA provides 9 customizeable ! termination codes, named t_xtra1 .. t_xtra9. You can ! customize the messages that will be printed upon exit by ! setting the corresponding termination_code_str value. ! termination_code_str(t_xtra1) = 'my termination condition' ! by default, indicate where (in the code) MESA terminated if (extras_check_model == terminate) s% termination_code = t_extras_check_model end function extras_check_model  Now we can modify our inlist to have, for example, x_ctrl(1) = 300 ! K; Tearth cut-off value x_ctrl(2) = 1d-6 ! K; Tearth tolerance  inside of the &controls section. When we ./mk and ./rn this, we will get a model having $T_\oplus\$ within a millionth of the specified value of 300 K.

Finally, we want these values to be reported in the LOGS/history.data file. We must first tell MESA that we want another column:

      integer function how_many_extra_history_columns(id, id_extra)
integer, intent(in) :: id, id_extra
integer :: ierr
type (star_info), pointer :: s
ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
how_many_extra_history_columns = 1
end function how_many_extra_history_columns


Now we need only to add the calculation to the names and vals vectors and it will appear in LOGS/history.data:

      subroutine data_for_extra_history_columns(id, id_extra, n, names, vals, ierr)
integer, intent(in) :: id, id_extra, n
character (len=maxlen_history_column_name) :: names(n)
real(dp) :: vals(n)
integer, intent(out) :: ierr
type (star_info), pointer :: s
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 log column options.
! it must not include the new column names you are adding here.

names(1) = "Tearth"
vals(1) = s% Teff * sqrt(s% photosphere_r * Rsun / (2.0 * AU))

end subroutine data_for_extra_history_columns


Of course we must ./mk again, and then if we ./rn we will obtain these numbers at each time step.

Figure: Blackbody temperature evolution over the history of planet Earth