Re: [mesa-users] Other kap

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

Not sure if anyone else had this problem, but I figured out the
solution with help from Aaron.  In my star/private/opacities.f file,
the call to s% other_kap_get_Type2() was missing the line "lnfree_e,
d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &".




On Mon, Sep 19, 2011 at 12:38 PM, Ken Shen <ke...@as...> wrote:
> 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
>