Merge pull request #525 from danieljprice/radstar

src/main/readwrite_infile.F90

@@ -111,7 +111,7 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
  use dust_formation,  only:write_options_dust_formation
  use nicil_sup,       only:write_options_nicil
  use metric,          only:write_options_metric
- use eos,             only:write_options_eos,ieos
+ use eos,             only:write_options_eos,ieos,X_in,Z_in
  use ptmass,          only:write_options_ptmass
  use ptmass_radiation,only:write_options_ptmass_radiation
  use cooling,         only:write_options_cooling
@@ -283,7 +283,12 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
     call write_inopt(exchange_radiation_energy,'gas-rad_exchange','exchange energy between gas and radiation',iwritein)
     call write_inopt(limit_radiation_flux,'flux_limiter','limit radiation flux',iwritein)
     call write_inopt(iopacity_type,'iopacity_type','opacity method (0=inf,1=mesa,2=constant,-1=preserve)',iwritein)
-    if (iopacity_type == 2) call write_inopt(kappa_cgs,'kappa_cgs','constant opacity value in cm2/g',iwritein)
+    if (iopacity_type == 1) then
+       call write_inopt(X_in,'X','hydrogen mass fraction for MESA opacity table',iwritein)
+       call write_inopt(Z_in,'Z','metallicity for MESA opacity table',iwritein)
+    elseif (iopacity_type == 2) then
+       call write_inopt(kappa_cgs,'kappa_cgs','constant opacity value in cm2/g',iwritein)
+    endif
     if (implicit_radiation) then
        call write_inopt(tol_rad,'tol_rad','tolerance on backwards Euler implicit solve of dxi/dt',iwritein)
        call write_inopt(itsmax_rad,'itsmax_rad','max number of iterations allowed in implicit solver',iwritein)

src/utils/moddump_LTE_to_rad.f90

@@ -39,7 +39,10 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
 
  if (.not. do_radiation) call fatal("moddump_LTE_to_rad","Not compiled with radiation")
 
- mu = gmw
+ X_in=0.687
+ Z_in=0.0142
+ mu = 0.61821
+ gmw = mu
  gamma_fixed = 5/3.  ! gamma should be exactly 5/3, because that is what ieos=12 assumes
  gamma = gamma_fixed
  print*,'Assuming gmw = ',mu,' and gamma=',gamma,'X = ',X_in,'Z = ',Z_in  ! X and Z are only used for calculating opacity
@@ -52,7 +55,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
     call get_idealplusrad_temp(rho_cgs,vxyzu(4,i)*unit_ergg,mu,tempi,ierr)
 
     ! calculate u and xi
-    ugasi = ugas_from_Tgas(tempi,gamma,gmw)
+    ugasi = ugas_from_Tgas(tempi,gamma,mu)
     vxyzu(4,i) = ugasi
     rad(iradxi,i) = radiation_and_gas_temperature_equal(rhoi,ugasi,gamma,mu)