colors: use Fe/H instead of zbase (#939)

Author: Debraheem

colors/defaults/colors.defaults

@@ -69,19 +69,65 @@
       sed_per_model = .false.
 
 
+   ! ``stellar_atm``
+   ! ~~~~~~~~~~~~~~~
+
+   ! Path to the directory containing the stellar atmosphere lookup table and
+   ! spectra used by the colors module. This directory should contain a
+   ! ``lookup_table.csv`` file and the corresponding SED files. The default
+   ! uses the ``Kurucz2003all`` ATLAS9 atmosphere grid shipped with MESA.
+   
+   ! ::
+
+      stellar_atm = '/data/colors_data/stellar_models/Kurucz2003all/'
+
+
+   ! ``z_over_x_ref``
+   ! ~~~~~~~~~~~~~~~~
+
+   ! ``z_over_x_ref`` is the reference metal-to-hydrogen ratio used to map
+   ! the photospheric composition onto the atmosphere-table metallicity axis:
+   ! [Fe/H] = log10((Z/X)/z_over_x_ref).
+   
+   ! The default matches the GS98 solar mixture used by the default
+   ! Kurucz2003all ATLAS9 atmosphere grid and equals
+   ! 0.0169/(1 - 0.0169 - 0.2485) = 2.30057173972d-2.
+   
+   ! If the photospheric metallicity, expressed as Z/X, falls outside the
+   ! metallicity range available in the tabulated atmosphere lookup table,
+   ! then MESA uses the nearest metallicity in the table for the atmosphere
+   ! interpolation. If the photospheric hydrogen mass fraction or metal mass
+   ! fraction is not positive, then MESA cannot form
+   ! log10((Z/X)/z_over_x_ref), so it uses the lowest metallicity in the
+   ! table instead. This is a fallback safeguard and not a dedicated treatment
+   ! for H-free atmospheres.
+   
+   ! For the shipped Kurucz2003all table in this checkout, the
+   ! metallicity range is [Fe/H] = -2.5 to [Fe/H] = 4.0. With the
+   ! default ``z_over_x_ref``, the lower edge corresponds to
+   ! Z/X ~ 7.275d-5.
+   
+   ! ::
+
+      z_over_x_ref = 2.30057173972d-2
+
+
 ! Extra inlist controls
 ! ---------------------
 
    ! One can split a colors inlist into pieces using the following parameters.
    ! It works recursively, so the extras can read extras too.
 
-         ! ``read_extra_colors_inlist(1..5)``
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ! ``extra_colors_inlist_name(1..5)``
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   ! ``read_extra_colors_inlist(1..5)``
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   ! ``extra_colors_inlist_name(1..5)``
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-         ! If ``read_extra_colors_inlist(i)`` is true, then read ``&colors`` from the file ``extra_colors_inlist_name(i)``.
-         ! ::
+   ! If ``read_extra_colors_inlist(i)`` is true, then read ``&colors``
+   ! from the file ``extra_colors_inlist_name(i)``.
+   !
+   
+   ! ::
 
       read_extra_colors_inlist(:) = .false.
       extra_colors_inlist_name(:) = 'undefined'

colors/private/colors_ctrls_io.f90

@@ -38,6 +38,7 @@ module colors_ctrls_io
    character(len=256) :: mag_system
 
    real(dp) :: distance
+   real(dp) :: z_over_x_ref
    logical :: make_csv
    logical :: sed_per_model
    logical :: use_colors
@@ -47,6 +48,7 @@ module colors_ctrls_io
       vega_sed, &
       stellar_atm, &
       distance, &
+      z_over_x_ref, &
       make_csv, &
       sed_per_model, &
       mag_system, &
@@ -155,6 +157,7 @@ subroutine store_controls(rq, ierr)
       rq%vega_sed = vega_sed
       rq%stellar_atm = stellar_atm
       rq%distance = distance
+      rq%z_over_x_ref = z_over_x_ref
       rq%make_csv = make_csv
       rq%sed_per_model = sed_per_model
       rq%colors_results_directory = colors_results_directory
@@ -192,6 +195,7 @@ subroutine set_controls_for_writing(rq)
       vega_sed = rq%vega_sed
       stellar_atm = rq%stellar_atm
       distance = rq%distance
+      z_over_x_ref = rq%z_over_x_ref
       make_csv = rq%make_csv
       sed_per_model = rq%sed_per_model
       colors_results_directory = rq%colors_results_directory

colors/private/colors_history.f90

@@ -112,8 +112,10 @@ subroutine data_for_colors_history_columns( &
             filter_name = trim(remove_dat(color_filter_names(i)))
             names(i + filter_offset) = filter_name
 
-            if (t_eff >= 0 .and. metallicity >= 0) then
-               ! pick the precomputed zero-point for the requested mag system
+            ! Negative [M/H] values are valid for metal-poor atmosphere grids.
+            ! so we do not apply a limit on the "metallicity" parameter.
+            if (t_eff >= 0) then
+               ! Select precomputed zero-point based on magnitude system
                select case (trim(cs%mag_system))
                case ('VEGA', 'Vega', 'vega')
                   zero_point = cs%filters(i)%vega_zero_point
@@ -152,4 +154,4 @@ subroutine data_for_colors_history_columns( &
 
    end subroutine data_for_colors_history_columns
 
-end module colors_history
+end module colors_history

colors/public/colors_def.f90

@@ -49,6 +49,7 @@ module colors_def
       character(len=256) :: mag_system
       real(dp) :: metallicity
       real(dp) :: distance
+      real(dp) :: z_over_x_ref
       logical :: make_csv
       logical :: sed_per_model
       logical :: use_colors
@@ -315,4 +316,4 @@ subroutine do_free_colors_tables
 
    end subroutine do_free_colors_tables
 
-end module colors_def
+end module colors_def

star/private/history.f90

@@ -70,7 +70,8 @@ subroutine do_history_info(s, write_flag, ierr)
          num_extra_cols, num_binary_cols, num_extra_binary_cols, num_colors_cols,num_extra_header_items, n
       integer, parameter :: num_epsnuc_out = 12
       real(dp) :: &
-         epsnuc_out(num_epsnuc_out), csound_surf, v_surf, envelope_fraction_left
+         epsnuc_out(num_epsnuc_out), csound_surf, v_surf, envelope_fraction_left, m_div_h, &
+         min_m_div_h, max_m_div_h
       integer :: mixing_regions, mix_relr_regions, burning_regions, burn_relr_regions
       integer, pointer :: mixing_type(:), mix_relr_type(:), burning_type(:), burn_relr_type(:)
       character (len = maxlen_history_column_name), pointer, dimension(:) :: &
@@ -268,9 +269,28 @@ subroutine do_history_info(s, write_flag, ierr)
          colors_col_names(1:num_colors_cols) = 'unknown'
          colors_col_vals(1:num_colors_cols) = -1d99
 
-         call data_for_colors_history_columns(s%T(1), safe_log10(s%grav(1)), s%R(1), s%kap_rq%Zbase, &
-            s% model_number, s% colors_handle, num_colors_cols, colors_col_names, colors_col_vals, ierr)
+         ! Here we compute [Fe/H], and thencall colors to compute history columns.
+         if (colors_settings% z_over_x_ref <= 0d0) then
+            write(*, *) 'colors error: z_over_x_ref must be positive'
+            ierr = -1
+            call dealloc
+            return
+         end if
+
+         min_m_div_h = minval(colors_settings% lu_meta)
+         max_m_div_h = maxval(colors_settings% lu_meta)
+
+         ! Map the current photospheric Z/X onto the atmosphere table metallicity axis.
+         if (s% X(s% photosphere_cell_k) > 0d0 .and. s% Z(s% photosphere_cell_k) > 0d0) then
+            m_div_h = log10((s% Z(s% photosphere_cell_k)/s% X(s% photosphere_cell_k)) / &
+               colors_settings% z_over_x_ref)
+            m_div_h = max(min_m_div_h, min(max_m_div_h, m_div_h))
+         else
+            m_div_h = min_m_div_h
+         end if
 
+         call data_for_colors_history_columns(s%T(1), log10(s%grav(1)), s%R(1), m_div_h, &
+            s% colors_handle, num_colors_cols, colors_col_names, colors_col_vals, ierr)
          if (ierr /= 0) then
             call dealloc
             return