From: Ken S. <ke...@as...> - 2011-09-23 18:59:20
|
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 > |