[mesa-users] Other kap

[Ken S. <ke...@as...>]

Hi all.  I wanted to modify the opacity used in MESA/star in a way
that the current controls can't do, so I tried to edit my
run_star_extras.f as outlined in other_kap.f.  Namely, I added these
lines to the extra_controls subroutine in run_star_extras.f:


----------------
         s% use_other_kap = .true.
         s% other_kap_get_Type1 => my_other_kap_get_Type1
----------------


As a simple test, my_other_kap_get_Type1 looks like this:


----------------
      subroutine my_other_kap_get_Type1( &
            id, handle, zbar, X, Zbase, log10_rho, log10_T,  &
            species, chem_id, net_iso, xa, &
            lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
            kap, dln_kap_dlnRho, dln_kap_dlnT, ierr)
         use kap_lib, only : kap_get_Type1

         ! INPUT
         integer, intent(in) :: id ! star id if available; 0 otherwise
         integer, intent(in) :: handle ! from alloc_other_kap_handle
         real*8, intent(in) :: zbar ! average ion charge
         real*8, intent(in) :: X ! the hydrogen mass fraction
         real*8, intent(in) :: Zbase ! the metallicity
         real*8, intent(in) :: log10_rho ! the density
         real*8, intent(in) :: log10_T ! the temperature
         integer, intent(in) :: species
         integer, pointer :: chem_id(:) ! maps species to chem id
            ! index from 1 to species
            ! value is between 1 and num_chem_isos
         integer, pointer :: net_iso(:) ! maps chem id to species number
            ! index from 1 to num_chem_isos (defined in chem_def)
            ! value is 0 if the iso is not in the current net
            ! else is value between 1 and number of species in current net
         real*8, intent(in) :: xa(:) ! mass fractions
         double precision, intent(in) :: lnfree_e, d_lnfree_e_dlnRho,
d_lnfree_e_dlnT
            ! free_e := total combined number per nucleon of free
electrons and positrons

         ! OUTPUT
         real*8, intent(out) :: kap ! opacity
         real*8, intent(out) :: dln_kap_dlnRho ! partial derivative at
constant T
         real*8, intent(out) :: dln_kap_dlnT   ! partial derivative at
constant Rho
         integer, intent(out) :: ierr ! 0 means AOK.


         call kap_get_Type1( &
            handle, zbar, X, Zbase, log10_rho, log10_T,  &
            lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
            kap, dln_kap_dlnRho, dln_kap_dlnT, ierr)

      end subroutine my_other_kap_get_Type1
----------------


However, these changes cause the run to exit immediately, without
producing any error messages:


----------------
DATE: 2011-09-19
TIME: 12:27:28
 read inlist_wd
 set_eos_PC_parameters
              mass_fraction_limit_for_PC    1.0000000000000000E-02
                        logRho1_PC_limit    2.0000000000000000E+01
                        logRho2_PC_limit    2.0000000000000000E+01
                      log_Gamma_all_HELM    1.6020600000000000E+00
                        log_Gamma_all_PC    1.9030899999999999E+00
                                PC_min_Z    9.9900000000000000E-01
load saved model 1_0_x12_0_50.mod

DATE: 2011-09-19
TIME: 12:27:32
----------------


>From inserting debugging messages, I know the lines "s% use_other_kap
= .true." and "s% other_kap_get_Type1 => my_other_kap_get_Type1" are
reached, but my_other_kap_get_Type1 doesn't ever seem to be called.
Any thoughts out there?


Thanks,
Ken