! $Id$
!
! This modules contains the routines for simulation with
! simple hydrogen ionization.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: leos = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 2
!
! PENCILS PROVIDED ss; gss(3); ee; pp; lnTT; cs2; cp; cp1; cp1tilde
! PENCILS PROVIDED glnTT(3); TT; TT1; gTT(3); yH; hss(3,3); hlnTT(3,3); del2TT; del6TT; del6lnTT
! PENCILS PROVIDED del2ss; del6ss; del2lnTT; cv; cv1; glnmumol(3); ppvap; csvap2
! PENCILS PROVIDED rho_anel; rho1gpp(3)
!
!***************************************************************
module EquationOfState
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
!
implicit none
!
include 'eos.h'
! integers specifying which independent variables to use in eoscalc
integer, parameter :: ilnrho_ss=1,ilnrho_ee=2,ilnrho_pp=3,ilnrho_lnTT=4
integer, parameter :: irho_ss=7, ilnrho_TT=9, irho_TT=10, ipp_ss=11
integer, parameter :: ipp_cs2=12
integer, parameter :: irho_eth=13, ilnrho_eth=14
! secondary parameters calculated in initialize
real :: TT_ion,lnTT_ion,TT_ion_,lnTT_ion_
real :: ss_ion,ee_ion,kappa0,xHe_term,ss_ion1,Srad0
real :: lnrho_e,lnrho_e_,lnrho_H,lnrho_He
! integer :: icp,ics
integer :: ics
! namelist parameters
real, parameter :: yHmin=tiny(TT_ion), yHmax=1-epsilon(TT_ion)
real :: xHe=0.1
real :: yMetals=0
real :: yHacc=1e-5
!ajwm Moved here from Density.f90
!ajwm Completely irrelevant to eos_ionization but density and entropy need
!ajwm reworking to be independent of these things first
!ajwm can't use impossible else it breaks reading param.nml
!ajwm SHOULDN'T BE HERE... But can wait till fully unwrapped
real :: cs0=impossible, rho0=impossible, cp=impossible, cv=impossible
real :: cs20=impossible, lnrho0=impossible
logical :: lpp_as_aux=.false., lcp_as_aux=.false.
real :: gamma=impossible, gamma_m1=impossible,gamma1=impossible
real :: cs2bot=impossible, cs2top=impossible
! input parameters
namelist /eos_init_pars/ xHe,yMetals,yHacc,lpp_as_aux,lcp_as_aux
! run parameters
namelist /eos_run_pars/ xHe,yMetals,yHacc,lpp_as_aux,lcp_as_aux
!ajwm Not sure this should exist either...
integer :: imass=1
!
real, dimension(nchemspec,18) :: species_constants
!
real :: Cp_const=impossible
real :: Pr_number=0.7
logical :: lpres_grad=.false.
!
contains
!***********************************************************************
subroutine register_eos
!
! 2-feb-03/axel: adapted from Interstellar module
! 13-jun-03/tobi: re-adapted from visc_shock module
!
use FArrayManager
use Sub
!
leos_ionization=.true.
!
! Set indices for auxiliary variables.
!
call farray_register_auxiliary('yH',iyH)
if (ilnTT==0) call farray_register_auxiliary('lnTT',ilnTT)
!
! Identify version number (generated automatically by SVN).
!
if (lroot) call svn_id( &
"$Id$")
!
! Writing files for use with IDL.
!
aux_var(aux_count)=',yh $'
aux_count=aux_count+1
if (naux < maux) aux_var(aux_count)=',lnTT $'
if (naux == maux) aux_var(aux_count)=',lnTT'
aux_count=aux_count+1
if (lroot) then
write(15,*) 'yH = fltarr(mx,my,mz)*one'
write(15,*) 'lnTT = fltarr(mx,my,mz)*one'
endif
!
endsubroutine register_eos
!***********************************************************************
subroutine units_eos
!
! If unit_temperature hasn't been specified explictly in start.in,
! set it to 1 (Kelvin).
!
! 24-jun-06/tobi: coded
!
if (unit_temperature==impossible) then
if (lfix_unit_std) then
unit_temperature=unit_density*unit_velocity**2/k_B_cgs
else
unit_temperature=1.
endif
endif
!
endsubroutine units_eos
!***********************************************************************
subroutine initialize_eos
!
! Perform any post-parameter-read initialization, e.g. set derived
! parameters.
!
! 2-feb-03/axel: adapted from Interstellar module
!
use General
use Mpicomm, only: stop_it
use Sub, only: register_report_aux
use SharedVariables,only: put_shared_variable
!
real :: mu1yHxHe
!
if (lroot) print*,'initialize_eos: ENTER'
!
! ionization parameters
! since m_e and chiH, as well as hbar are all very small
! it is better to divide m_e and chiH separately by hbar.
!
mu1yHxHe=1+3.97153*xHe
TT_ion=chiH/k_B
lnTT_ion=log(TT_ion)
TT_ion_=chiH_/k_B
lnTT_ion_=log(chiH_/k_B)
lnrho_e=1.5*log((m_e/hbar)*(chiH/hbar)/2./pi)+log(m_H)+log(mu1yHxHe)
lnrho_H=1.5*log((m_H/hbar)*(chiH/hbar)/2./pi)+log(m_H)+log(mu1yHxHe)
lnrho_He=1.5*log((m_He/hbar)*(chiH/hbar)/2./pi)+log(m_H)+log(mu1yHxHe)
lnrho_e_=1.5*log((m_e/hbar)*(chiH_/hbar)/2./pi)+log(m_H)+log(mu1yHxHe)
ss_ion=k_B/m_H/mu1yHxHe
ss_ion1=1/ss_ion
ee_ion=ss_ion*TT_ion
kappa0=sigmaH_/m_H/mu1yHxHe/4.0
Srad0=sigmaSB*TT_ion**4.0D0/pi
!
if (xHe>0) then
xHe_term=xHe*(log(xHe)-lnrho_He)
elseif (xHe<0) then
call stop_it('initialize_eos: xHe lower than zero makes no sense')
else
xHe_term=0
endif
!
! pressure and cp as optional auxiliary variable
!
if (lpp_as_aux) call register_report_aux('pp',ipp)
if (lcp_as_aux) call register_report_aux('cp',icp)
!
if (lroot.and.ip<14) then
print*,'initialize_eos: reference values for ionization'
print*,'initialize_eos: yHmin,yHmax,yMetals=',yHmin,yHmax,yMetals
print*,'initialize_eos: TT_ion,ss_ion,kappa0=', &
TT_ion,ss_ion,kappa0
print*,'initialize_eos: lnrho_e,lnrho_H,lnrho_He,lnrho_e_=', &
lnrho_e,lnrho_H,lnrho_He,lnrho_e_
endif
if (.not.ldensity) then
call put_shared_variable('rho0',rho0,caller='initialize_eos')
call put_shared_variable('lnrho0',lnrho0)
endif
!
! write scale non-free constants to file; to be read by idl
!
if (lroot) then
open (1,file=trim(datadir)//'/pc_constants.pro',position="append")
write (1,*) 'TT_ion=',TT_ion
write (1,*) 'lnTT_ion=',lnTT_ion
write (1,*) 'TT_ion_=',TT_ion_
write (1,*) 'lnTT_ion_=',lnTT_ion_
write (1,*) 'lnrho_e=',lnrho_e
write (1,*) 'lnrho_H=',lnrho_H
write (1,*) 'lnrho_p=',lnrho_H
write (1,*) 'lnrho_He=',lnrho_He
write (1,*) 'lnrho_e_=',lnrho_e_
write (1,*) 'ss_ion=',ss_ion
write (1,*) 'ee_ion=',ee_ion
write (1,*) 'kappa0=',kappa0
write (1,*) 'Srad0=',Srad0
write (1,*) 'k_B=',k_B
write (1,*) 'm_H=',m_H
close (1)
endif
!
endsubroutine initialize_eos
!***********************************************************************
subroutine select_eos_variable(variable,findex)
!
! Calculate average particle mass in the gas relative to
!
! 02-apr-06/tony: implemented
!
character (len=*), intent(in) :: variable
integer, intent(in) :: findex
!
call keep_compiler_quiet(variable)
call keep_compiler_quiet(findex)
! DUMMY ideagas version below
!! integer :: this_var=0
!! integer, save :: ieosvar=0
!! integer, save :: ieosvar_selected=0
!! integer, parameter :: ieosvar_lnrho = 1
!! integer, parameter :: ieosvar_rho = 2
!! integer, parameter :: ieosvar_ss = 4
!! integer, parameter :: ieosvar_lnTT = 8
!! integer, parameter :: ieosvar_TT = 16
!! integer, parameter :: ieosvar_cs2 = 32
!! integer, parameter :: ieosvar_pp = 64
!!!
!! if (ieosvar>=2) &
!! call fatal_error("select_eos_variable", &
!! "2 thermodynamic quantities have already been defined while attempting to add a 3rd: ") !//variable)
!!
!! ieosvar=ieosvar+1
!!
!!! select case (variable)
!! if (variable=='ss') then
!! this_var=ieosvar_ss
!! if (findex<0) then
!! leos_isentropic=.true.
!! endif
!! elseif (variable=='cs2') then
!! this_var=ieosvar_cs2
!! if (findex==-2) then
!! leos_localisothermal=.true.
!! elseif (findex<0) then
!! leos_isothermal=.true.
!! endif
!! elseif (variable=='lnTT') then
!! this_var=ieosvar_lnTT
!! if (findex<0) then
!! leos_isothermal=.true.
!! endif
!! elseif (variable=='lnrho') then
!! this_var=ieosvar_lnrho
!! if (findex<0) then
!! leos_isochoric=.true.
!! endif
!! elseif (variable=='rho') then
!! this_var=ieosvar_rho
!! if (findex<0) then
!! leos_isochoric=.true.
!! endif
!! elseif (variable=='pp') then
!! this_var=ieosvar_pp
!! if (findex<0) then
!! leos_isobaric=.true.
!! endif
!! else
!! call fatal_error("select_eos_variable", &
!! "unknown thermodynamic variable")
!! endif
!! if (ieosvar==1) then
!! ieosvar1=findex
!! ieosvar_selected=ieosvar_selected+this_var
!! return
!! endif
!!!
!!! Ensure the indexes are in the correct order.
!!!
!! if (this_var<ieosvar_selected) then
!! ieosvar2=ieosvar1
!! ieosvar1=findex
!! else
!! ieosvar2=findex
!! endif
!! ieosvar_selected=ieosvar_selected+this_var
!! select case (ieosvar_selected)
!! case (ieosvar_lnrho+ieosvar_ss)
!! ieosvars=ilnrho_ss
!! case (ieosvar_lnrho+ieosvar_lnTT)
!! ieosvars=ilnrho_lnTT
!! case (ieosvar_lnrho+ieosvar_cs2)
!! ieosvars=ilnrho_cs2
!! case default
!! print*,"select_eos_variable: Thermodynamic variable combination, ieosvar_selected= ",ieosvar_selected
!! call fatal_error("select_eos_variable", &
!! "This thermodynamic variable combination is not implemented: ")
!! endselect
!
endsubroutine select_eos_variable
!***********************************************************************
subroutine rprint_eos(lreset,lwrite)
!
logical :: lreset
logical, optional :: lwrite
!
call keep_compiler_quiet(lreset)
call keep_compiler_quiet(present(lwrite))
!
! check for those quantities for which we want video slices
!
if (lwrite_slices) then
where(cnamev=='lnTT'.or.cnamev=='yH') cformv='DEFINED'
endif
!
endsubroutine rprint_eos
!***********************************************************************
subroutine get_slices_eos(f,slices)
!
! Write slices for animation of Eos variables.
!
! 26-jul-06/tony: coded
!
use Slices_methods, only: assign_slices_scal
!
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
!
! Loop over slices.
!
select case (trim(slices%name))
!
! Temperature.
!
case ('lnTT'); call assign_slices_scal(slices,f,ilnTT)
!
! Degree of ionization.
!
case ('yH'); call assign_slices_scal(slices,f,iyH)
!
endselect
!
endsubroutine get_slices_eos
!***********************************************************************
subroutine pencil_criteria_eos
!
! All pencils that the EquationOfState module depends on are specified here.
!
! 02-apr-06/tony: coded
!
! EOS is a pencil provider but evolves nothing so it is unlokely that
! it will require any pencils for it's own use.
!
! pp pencil if lpp_as_aux
!
if (lpp_as_aux) lpenc_requested(i_pp)=.true.
if (lcp_as_aux) lpenc_requested(i_cp1tilde)=.true.
!
endsubroutine pencil_criteria_eos
!***********************************************************************
subroutine pencil_interdep_eos(lpencil_in)
!
! Interdependency among pencils from the Entropy module is specified here.
!
! 20-11-04/anders: coded
!
logical, dimension(npencils) :: lpencil_in
!
if (lpencil_in(i_gTT)) then
lpencil_in(i_glnTT)=.true.
lpencil_in(i_TT)=.true.
endif
if (lpencil_in(i_del2lnTT)) then
lpencil_in(i_del2lnrho)=.true.
lpencil_in(i_del2ss)=.true.
endif
if (lpencil_in(i_glnTT)) then
lpencil_in(i_glnrho)=.true.
lpencil_in(i_gss)=.true.
endif
if (lpencil_in(i_TT)) lpencil_in(i_lnTT)=.true.
if (lpencil_in(i_TT1)) lpencil_in(i_lnTT)=.true.
!
if (lpencil_in(i_hlnTT)) then
lpencil_in(i_hss)=.true.
if (.not.pretend_lnTT) lpencil_in(i_hlnrho)=.true.
endif
!
endsubroutine pencil_interdep_eos
!***********************************************************************
subroutine calc_pencils_eos_std(f,p)
!
! Envelope adjusting calc_pencils_eos_pencpar to the standard use with
! lpenc_loc=lpencil
!
! 9-oct-15/MR: coded
!
real, dimension (mx,my,mz,mfarray),intent(INOUT):: f
type (pencil_case), intent(INOUT):: p
!
call calc_pencils_eos_pencpar(f,p,lpencil)
!
endsubroutine calc_pencils_eos_std
!***********************************************************************
subroutine calc_pencils_eos_pencpar(f,p,lpenc_loc)
!
! Calculate Entropy pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 02-apr-06/tony: coded
! 09-oct-15/MR: added mask parameter lpenc_loc
!
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
logical, dimension(npencils) :: lpenc_loc
!
intent(inout):: f,p
intent(in) :: lpenc_loc
!
integer :: i
!
! THE FOLLOWING 2 ARE CONCEPTUALLY WRONG
! FOR pretend_lnTT since iss actually contain lnTT NOT entropy!
! The code is not wrong however since this is correctly
! handled by the eos module.
! ss
if (lpenc_loc(i_ss)) p%ss=f(l1:l2,m,n,iss)
!
! gss
if (lpenc_loc(i_gss)) call grad(f,iss,p%gss)
! pp
if (lpenc_loc(i_pp)) call eoscalc(f,nx,pp=p%pp)
! ee
if (lpenc_loc(i_ee)) call eoscalc(f,nx,ee=p%ee)
! lnTT
if (lpenc_loc(i_lnTT)) call eoscalc(f,nx,lnTT=p%lnTT)
! yH
if (lpenc_loc(i_yH)) call eoscalc(f,nx,yH=p%yH)
! TT
if (lpenc_loc(i_TT)) p%TT=exp(p%lnTT)
! TT1
if (lpenc_loc(i_TT1)) p%TT1=exp(-p%lnTT)
! cs2 and cp1tilde
if (lpenc_loc(i_cs2) .or. lpenc_loc(i_cp1tilde)) &
call pressure_gradient(f,p%cs2,p%cp1tilde)
! glnTT
if (lpenc_loc(i_glnTT)) then
call temperature_gradient(f,p%glnrho,p%gss,p%glnTT)
endif
! gTT
if (lpenc_loc(i_gTT)) then
do i=1,3; p%gTT(:,i)=p%glnTT(:,i)*p%TT; enddo
endif
! hss
if (lpenc_loc(i_hss)) then
call g2ij(f,iss,p%hss)
endif
! del2ss
if (lpenc_loc(i_del2ss)) then
call del2(f,iss,p%del2ss)
endif
! del2lnTT
if (lpenc_loc(i_del2lnTT)) then
call temperature_laplacian(f,p)
endif
! del2TT
if (lpenc_loc(i_del2TT)) then
call fatal_error('calc_pencils_eos_pencpar','del2TT not implemented')
endif
! del6lnTT
if (lpenc_loc(i_del6lnTT)) then
call fatal_error('calc_pencils_eos_pencpar','del6lnTT not implemented')
endif
! del6TT
if (lpenc_loc(i_del6TT)) then
call fatal_error('calc_pencils_eos_pencpar','del6TT not implemented')
endif
! del6ss
if (lpenc_loc(i_del6ss)) then
call del6(f,iss,p%del6ss)
endif
! hlnTT
if (lpenc_loc(i_hlnTT)) then
call temperature_hessian(f,p%hlnrho,p%hss,p%hlnTT)
endif
!
if (lpenc_loc(i_glnmumol)) p%glnmumol(:,:)=0.
!
! pressure and cp as optional auxiliary pencils
!
if (lpp_as_aux) f(l1:l2,m,n,ipp)=p%pp
if (lcp_as_aux) f(l1:l2,m,n,icp)=p%cp1tilde
!
! This routine does not yet compute cv or cv1, but since those pencils
! are supposed to be provided here, we better set them to impossible.
!
if (lpenc_loc(i_cv1)) p%cv1=impossible
if (lpenc_loc(i_cp1)) p%cp1=impossible
if (lpenc_loc(i_cv)) p%cv=impossible
if (lpenc_loc(i_cp)) p%cp=impossible
!
endsubroutine calc_pencils_eos_pencpar
!***********************************************************************
subroutine getmu(f,mu_tmp)
!
! Calculate average particle mass.
! Note that the particles density is N = nHI + nHII + ne + nHe
! = (1-y)*nH + y*nH + y*nH + xHe*nH = (1 + yH + xHe) * nH, where
! nH is the number of protons per cubic centimeter.
! The number of particles per mole is therefore 1 + yH + xHe.
! The mass per mole is M=1.+3.97153*xHe, so the mean molecular weight
! per particle is M/N = (1.+3.97153*xHe)/(1 + yH + xHe).
!
! 12-aug-03/tony: implemented
!
real, dimension (mx,my,mz,mfarray), optional :: f
real, optional, intent(out) :: mu_tmp
!
mu_tmp=1.+3.97153*xHe
!
! tobi: the real mean molecular weight would be:
!
! mu_tmp=(1.+3.97153*xHe)/(1+yH+xHe)
!
call keep_compiler_quiet(present(f))
!
endsubroutine getmu
!***********************************************************************
subroutine ioninit(f)
!
! the ionization fraction has to be set to a value yH0 < yH < yHmax before
! rtsafe is called for the first time
!
! 12-jul-03/tobi: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
f(:,:,:,iyH) = 0.5*(yHmax-yHmin)
call ioncalc(f)
!
endsubroutine ioninit
!***********************************************************************
subroutine ioncalc(f)
!
! calculate degree of ionization and temperature
! This routine is called from equ.f90 and operates on the full 3-D array.
!
! 13-jun-03/tobi: coded
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx) :: lnrho,ss,yH,lnTT
!
do n=1,mz
do m=1,my
if (ldensity_nolog) then
lnrho=log(f(:,m,n,ilnrho))
else
lnrho=f(:,m,n,ilnrho)
endif
ss=f(:,m,n,iss)
yH=f(:,m,n,iyH)
call rtsafe_pencil(lnrho,ss,yH)
f(:,m,n,iyH)=yH
lnTT=(ss/ss_ion+(1-yH)*(log(1-yH+epsi)-lnrho_H) &
+yH*(2*log(yH)-lnrho_e-lnrho_H)+xHe_term)/(1+yH+xHe)
lnTT=(2.0/3.0)*(lnTT+lnrho-2.5)+lnTT_ion
f(:,m,n,ilnTT)=lnTT
enddo
enddo
!
endsubroutine ioncalc
!***********************************************************************
subroutine getdensity(EE,TT,yH,rho)
!
! calculate density. Is currently only being used by the interstellar
! module. I guess we can/should replace this now by a call to eoscalc.
!
real, intent(in) :: EE,TT,yH
real, intent(out) :: rho
!
rho=EE/(1.5*(1.+yH+xHe)*ss_ion*TT+yH*ee_ion)
!
endsubroutine getdensity
!***********************************************************************
subroutine gettemperature(f,TT_tmp)
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz), intent(out) :: TT_tmp
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(TT_tmp)
!
endsubroutine gettemperature
!***********************************************************************
subroutine getpressure(pp_tmp,TT_tmp,rho_tmp,mu1_tmp)
!
real, dimension (nx), intent(out) :: pp_tmp
real, dimension (nx), intent(in) :: TT_tmp,rho_tmp,mu1_tmp
!
call fatal_error('getpressure','Should not be called with eos_ionization.')
!
call keep_compiler_quiet(pp_tmp)
call keep_compiler_quiet(TT_tmp)
call keep_compiler_quiet(rho_tmp)
call keep_compiler_quiet(mu1_tmp)
!
endsubroutine getpressure
!***********************************************************************
subroutine get_cp1(cp1_)
!
! 04-nov-06/axel: added to alleviate spurious use of pressure_gradient
!
! return the value of cp1 to outside modules
!
real, intent(out) :: cp1_
!
! for variable ionization, it doesn't make sense to calculate
! just a single value of cp1, because it must depend on position.
! Therefore, return impossible, so one can reconsider this case.
!
call fatal_error('get_cp1','SHOULD NOT BE CALLED WITH eos_ionization')
cp1_=impossible
!
endsubroutine get_cp1
!***********************************************************************
subroutine get_cv1(cv1_)
!
! 22-dec-10/PJK: adapted from get_cp1
!
! return the value of cv1 to outside modules
!
real, intent(out) :: cv1_
!
! for variable ionization, it doesn't make sense to calculate
! just a single value of cv1, because it must depend on position.
! Therefore, return impossible, so one can reconsider this case.
!
call fatal_error('get_cv1','SHOULD NOT BE CALLED WITH eos_ionization')
cv1_=impossible
!
endsubroutine get_cv1
!***********************************************************************
subroutine pressure_gradient_farray(f,cs2,cp1tilde)
!
! Calculate thermodynamical quantities, cs2 and cp1tilde
! and optionally glnPP and glnTT
! gP/rho=cs2*(glnrho+cp1tilde*gss)
!
! 17-nov-03/tobi: adapted from subroutine eoscalc
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
real, dimension(nx), intent(out) :: cs2,cp1tilde
real, dimension(nx) :: lnrho,yH,lnTT
real, dimension(nx) :: R,dlnTTdy,dRdy,temp
real, dimension(nx) :: dlnPPdlnrho,fractions,fractions1
real, dimension(nx) :: dlnPPdss,TT1
!
lnrho=f(l1:l2,m,n,ilnrho)
yH=f(l1:l2,m,n,iyH)
lnTT=f(l1:l2,m,n,ilnTT)
TT1=exp(-lnTT)
fractions=(1+yH+xHe)
fractions1=1/fractions
!
R=lnrho_e-lnrho+1.5*(lnTT-lnTT_ion)-TT_ion*TT1+log(1-yH+epsi)-2*log(yH)
dlnTTdy=(2*(-R-TT_ion*TT1)-3)/3*fractions1
dRdy=dlnTTdy*(1.5+TT_ion*TT1)-1/(1-yH+epsi)-2/yH
temp=(dlnTTdy+fractions1)/dRdy
dlnPPdlnrho=(5-2*TT_ion*TT1*temp)/3
dlnPPdss=ss_ion1*fractions1*(dlnPPdlnrho-temp-1)
cs2=fractions*ss_ion*dlnPPdlnrho/TT1
cp1tilde=dlnPPdss/dlnPPdlnrho
!
endsubroutine pressure_gradient_farray
!***********************************************************************
subroutine pressure_gradient_point(lnrho,ss,cs2,cp1tilde)
!
! Calculate thermodynamical quantities, cs2 and cp1tilde
! and optionally glnPP and glnTT
! gP/rho=cs2*(glnrho+cp1tilde*gss)
!
! 17-nov-03/tobi: adapted from subroutine eoscalc
!
real, intent(in) :: lnrho,ss
real, intent(out) :: cs2,cp1tilde
real :: yH,lnTT
real :: R,dlnTTdy,dRdy,temp
real :: dlnPPdlnrho,fractions,fractions1
real :: dlnPPdss,TT1
!
yH=0.5
call rtsafe(ilnrho_ss,lnrho,ss,yHmin,yHmax,yH)
fractions=(1+yH+xHe)
fractions1=1/fractions
lnTT=(2.0/3.0)*((ss/ss_ion+(1-yH)*(log(1-yH+epsi)-lnrho_H) &
+yH*(2*log(yH)-lnrho_e-lnrho_H) &
+xHe_term)/fractions+lnrho-2.5)+lnTT_ion
TT1=exp(-lnTT)
!
R=lnrho_e-lnrho+1.5*(lnTT-lnTT_ion)-TT_ion*TT1+log(1-yH+epsi)-2*log(yH)
dlnTTdy=(2*(-R-TT_ion*TT1)-3)/3*fractions1
dRdy=dlnTTdy*(1.5+TT_ion*TT1)-1/(1-yH+epsi)-2/yH
temp=(dlnTTdy+fractions1)/dRdy
dlnPPdlnrho=(5-2*TT_ion*TT1*temp)/3
dlnPPdss=ss_ion1*fractions1*(dlnPPdlnrho-temp-1)
cs2=fractions*ss_ion*dlnPPdlnrho/TT1
cp1tilde=dlnPPdss/dlnPPdlnrho
!
endsubroutine pressure_gradient_point
!***********************************************************************
subroutine temperature_gradient(f,glnrho,gss,glnTT)
!
! Calculate thermodynamical quantities, cs2 and cp1tilde
! and optionally glnPP and glnTT
! gP/rho=cs2*(glnrho+cp1tilde*gss)
!
! 17-nov-03/tobi: adapted from subroutine eoscalc
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
real, dimension(nx,3), intent(in) :: glnrho,gss
real, dimension(nx,3), intent(out) :: glnTT
real, dimension(nx) :: lnrho,yH,lnTT,TT1,fractions1
real, dimension(nx) :: R,dlnTTdy,dRdy
real, dimension(nx) :: dlnTTdydRdy,dlnTTdlnrho,dlnTTdss
integer :: j
!
lnrho=f(l1:l2,m,n,ilnrho)
yH=f(l1:l2,m,n,iyH)
lnTT=f(l1:l2,m,n,ilnTT)
TT1=exp(-lnTT)
fractions1=1/(1+yH+xHe)
!
R=lnrho_e-lnrho+1.5*(lnTT-lnTT_ion)-TT_ion*TT1+log(1-yH+epsi)-2*log(yH)
dlnTTdy=((2.0/3.0)*(-R-TT_ion*TT1)-1)*fractions1
dRdy=dlnTTdy*(1.5+TT_ion*TT1)-1/(1-yH+epsi)-2/yH
dlnTTdydRdy=dlnTTdy/dRdy
dlnTTdlnrho=(2.0/3.0)*(1-TT_ion*TT1*dlnTTdydRdy)
dlnTTdss=(dlnTTdlnrho-dlnTTdydRdy)*fractions1*ss_ion1
do j=1,3
glnTT(:,j)=dlnTTdlnrho*glnrho(:,j)+dlnTTdss*gss(:,j)
enddo
!
endsubroutine temperature_gradient
!***********************************************************************
subroutine temperature_laplacian(f,p)
!
! Calculate thermodynamical quantities, cs2 and cp1tilde
! and optionally glnPP and glnTT
! gP/rho=cs2*(glnrho+cp1tilde*gss)
!
! 12-dec-05/tony: adapted from subroutine temperature_gradient
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
type (pencil_case) :: p
!
call not_implemented('temperature_laplacian')
!
p%del2lnTT=0.0
call keep_compiler_quiet(f)
!
endsubroutine temperature_laplacian
!***********************************************************************
subroutine temperature_hessian(f,hlnrho,hss,hlnTT)
!
! Calculate thermodynamical quantities, cs2 and cp1tilde
! and optionally hlnPP and hlnTT
! hP/rho=cs2*(hlnrho+cp1tilde*hss)
!
! 10-apr-04/axel: adapted from temperature_gradient
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
real, dimension(nx,3,3), intent(in) :: hlnrho,hss
real, dimension(nx,3,3), intent(out) :: hlnTT
real, dimension(nx) :: lnrho,yH,lnTT,TT1,fractions1
real, dimension(nx) :: R,dlnTTdy,dRdy
real, dimension(nx) :: dlnTTdydRdy,dlnTTdlnrho,dlnTTdss
integer :: i,j
!
lnrho=f(l1:l2,m,n,ilnrho)
yH=f(l1:l2,m,n,iyH)
lnTT=f(l1:l2,m,n,ilnTT)
TT1=exp(-lnTT)
fractions1=1/(1+yH+xHe)
!
R=lnrho_e-lnrho+1.5*(lnTT-lnTT_ion)-TT_ion*TT1+log(1-yH+epsi)-2*log(yH)
dlnTTdy=((2.0/3.0)*(-R-TT_ion*TT1)-1)*fractions1
dRdy=dlnTTdy*(1.5+TT_ion*TT1)-1/(1-yH+epsi)-2/yH
dlnTTdydRdy=dlnTTdy/dRdy
dlnTTdlnrho=(2.0/3.0)*(1-TT_ion*TT1*dlnTTdydRdy)
dlnTTdss=(dlnTTdlnrho-dlnTTdydRdy)*fractions1*ss_ion1
do j=1,3
do i=1,3
hlnTT(:,i,j)=dlnTTdlnrho*hlnrho(:,i,j)+dlnTTdss*hss(:,i,j)
enddo
enddo
!
endsubroutine temperature_hessian
!***********************************************************************
subroutine eosperturb(f,psize,ee,pp,ss)
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
integer, intent(in) :: psize
real, dimension(psize), intent(in), optional :: ee,pp,ss
!
call not_implemented("eosperturb")
call keep_compiler_quiet(f)
call keep_compiler_quiet(present(ee),present(pp),present(ss))
!
endsubroutine eosperturb
!***********************************************************************
subroutine eoscalc_farray(f,psize,lnrho,ss,yH,lnTT,ee,pp,cs2,kapparho)
!
! Calculate thermodynamical quantities
!
! 2-feb-03/axel: simple example coded
! 13-jun-03/tobi: the ionization fraction as part of the f-array
! now needs to be given as an argument as input
! 17-nov-03/tobi: moved calculation of cs2 and cp1tilde to
! subroutine pressure_gradient
!
use Sub
use Mpicomm, only: stop_it
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
integer, intent(in) :: psize
real, dimension(psize), intent(out), optional :: lnrho,ss
real, dimension(psize), intent(out), optional :: yH,lnTT
real, dimension(psize), intent(out), optional :: ee,pp,kapparho
real, dimension(psize), optional :: cs2
real, dimension(psize) :: lnrho_,ss_,yH_,lnTT_,TT_,fractions,exponent
!
select case (psize)
!
case (nx)
lnrho_=f(l1:l2,m,n,ilnrho)
ss_=f(l1:l2,m,n,iss)
yH_=f(l1:l2,m,n,iyH)
lnTT_=f(l1:l2,m,n,ilnTT)
!
case (mx)
lnrho_=f(:,m,n,ilnrho)
ss_=f(:,m,n,iss)
yH_=f(:,m,n,iyH)
lnTT_=f(:,m,n,ilnTT)
!
case default
call stop_it("eoscalc: no such pencil size")
!
end select
!
TT_=exp(lnTT_)
fractions=(1+yH_+xHe)
!
if (present(lnrho)) lnrho=lnrho_
if (present(ss)) ss=ss_
if (present(yH)) yH=yH_
if (present(lnTT)) lnTT=lnTT_
if (present(ee)) ee=1.5*fractions*ss_ion*TT_+yH_*ee_ion
if (present(pp)) pp=fractions*exp(lnrho_)*TT_*ss_ion
if (present(cs2)) &
call fatal_error('eoscalc_farray','calculation of cs2 not implemented')
!
! Hminus opacity
!
if (present(kapparho)) then
! lnchi=2*lnrho_-lnrho_e_+1.5*(lnTT_ion_-lnTT_) &
! +TT_ion_/TT_+log(yH_+yMetals)+log(1-yH_+epsi)+lnchi0
exponent = (2*lnrho_-lnrho_e_+1.5*(lnTT_ion_-lnTT_)+TT_ion_/TT_)
!
! Ensure exponentiation and successive multiplication with
! numbers up to ~ 1e4 avoids overflow:
!
exponent = min(exponent,log(huge1)-5.)
kapparho = exp(exponent + alog(yH_+yMetals))*(1-yH_)*kappa0
endif
!
endsubroutine eoscalc_farray
!***********************************************************************
subroutine eoscalc_pencil(ivars,var1,var2,lnrho,ss,yH,lnTT,ee,pp)
!
! Calculate thermodynamical quantities
!
! i13-mar-04/tony: modified
!
use Mpicomm, only: stop_it
!
integer, intent(in) :: ivars
real, dimension(nx), intent(in) :: var1,var2
real, dimension(nx), intent(out), optional :: lnrho,ss
real, dimension(nx), intent(out), optional :: yH,lnTT
real, dimension(nx), intent(out), optional :: ee,pp
real, dimension(nx) :: lnrho_,ss_,yH_,lnTT_,TT_,TT1_,rho_,ee_,pp_
real, dimension(nx) :: fractions,rhs,sqrtrhs
integer :: i
!
select case (ivars)
!
case (ilnrho_ss,irho_ss)
if (ivars==ilnrho_ss) then
lnrho_=var1
else
lnrho_=alog(var1)
endif
ss_=var2
yH_=0.5*yHmax
do i=1,nx
call rtsafe(ilnrho_ss,lnrho_(i),ss_(i),yHmin,yHmax,yH_(i))
enddo
fractions=(1+yH_+xHe)
lnTT_=(2.0/3.0)*((ss_/ss_ion+(1-yH_)*(log(1-yH_+epsi)-lnrho_H) &
+yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
+xHe_term)/fractions+lnrho_-2.5)+lnTT_ion
TT_=exp(lnTT_)
rho_=exp(lnrho_)
ee_=1.5*fractions*ss_ion*TT_+yH_*ee_ion
pp_=fractions*rho_*TT_*ss_ion
!
case (ilnrho_lnTT)
lnrho_=var1
lnTT_=var2
TT_=exp(lnTT_)
TT1_=exp(-lnTT_)
rhs=exp(lnrho_e-lnrho_+1.5*(lnTT_-lnTT_ion)-TT_ion*TT1_)
rhs=max(rhs,tini) ! avoid log(0.) below
sqrtrhs=sqrt(rhs)
yH_=2*sqrtrhs/(sqrtrhs+sqrt(4+rhs))
fractions=(1+yH_+xHe)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
ee_=1.5*fractions*ss_ion*TT_+yH_*ee_ion
pp_=(1+yH_+xHe)*exp(lnrho_)*TT_*ss_ion
!
case (ilnrho_ee)
lnrho_=var1
ee_=var2
yH_=yHmax
yH_=0.5*min(ee_/ee_ion,yH_)
do i=1,nx
call rtsafe(ilnrho_ee,lnrho_(i),ee_(i),yHmin,yHmax*min(ee_(i)/ee_ion,1.0),yH_(i))
enddo
fractions=(1+yH_+xHe)
TT_=(ee_-yH_*ee_ion)/(1.5*fractions*ss_ion)
lnTT_=log(TT_)
rho_=exp(lnrho_)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
pp_=fractions*rho_*TT_*ss_ion
!
case (ilnrho_pp)
lnrho_=var1
pp_=var2
yH_=0.5*yHmax
do i=1,nx
call rtsafe(ilnrho_pp,lnrho_(i),pp_(i),yHmin,yHmax,yH_(i))
enddo
fractions=(1+yH_+xHe)
rho_=exp(lnrho_)
TT_=pp_/(fractions*ss_ion*rho_)
lnTT_=log(TT_)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
ee_=1.5*fractions*ss_ion*TT_+yH_*ee_ion
!
case default
call stop_it("eoscalc_pencil: I don't get what the independent variables are.")
!
end select
!
if (present(lnrho)) lnrho=lnrho_
if (present(ss)) ss=ss_
if (present(yH)) yH=yH_
if (present(lnTT)) lnTT=lnTT_
if (present(ee)) ee=ee_
if (present(pp)) pp=pp_
!
endsubroutine eoscalc_pencil
!***********************************************************************
subroutine eoscalc_point(ivars,var1,var2,lnrho,ss,yH,lnTT,ee,pp,cs2)
!
! Calculate thermodynamical quantities
!
! 2-feb-03/axel: simple example coded
! 13-jun-03/tobi: the ionization fraction as part of the f-array
! now needs to be given as an argument as input
! 17-nov-03/tobi: moved calculation of cs2 and cp1tilde to
! subroutine pressure_gradient
!
use Mpicomm, only: stop_it
!
integer, intent(in) :: ivars
real, intent(in) :: var1,var2
real, intent(out), optional :: lnrho,ss
real, intent(out), optional :: yH,lnTT
real, intent(out), optional :: ee,pp,cs2
real :: lnrho_,ss_,yH_,lnTT_,TT_,TT1_,rho_,ee_,pp_,cs2_
real :: fractions,rhs,sqrtrhs
!
select case (ivars)
!
case (ilnrho_ss)
lnrho_=var1
ss_=var2
yH_=0.5*yHmax
call rtsafe(ilnrho_ss,lnrho_,ss_,yHmin,yHmax,yH_)
lnTT_=(ss_/ss_ion+(1-yH_)*(log(1-yH_+epsi)-lnrho_H) &
+yH_*(2*log(yH_)-lnrho_e-lnrho_H)+xHe_term)/(1+yH_+xHe)
lnTT_=(2.0/3.0)*(lnTT_+lnrho_-2.5)+lnTT_ion
!
TT_=exp(lnTT_)
rho_=exp(lnrho_)
ee_=1.5*(1+yH_+xHe)*ss_ion*TT_+yH_*ee_ion
pp_=(1+yH_+xHe)*rho_*TT_*ss_ion
cs2_=impossible
!
case (ilnrho_lnTT)
lnrho_=var1
lnTT_=var2
TT_=exp(lnTT_)
TT1_=1/TT_
rhs=exp(lnrho_e-lnrho_+1.5*(lnTT_-lnTT_ion)-TT_ion*TT1_)
rhs = max(rhs,tini) ! avoid log(0.) below
sqrtrhs=sqrt(rhs)
yH_=2*sqrtrhs/(sqrtrhs+sqrt(4+rhs))
fractions=(1+yH_+xHe)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
ee_=1.5*fractions*ss_ion*TT_+yH_*ee_ion
pp_=(1+yH_+xHe)*exp(lnrho_)*TT_*ss_ion
cs2_=impossible
!
case (ilnrho_ee)
lnrho_=var1
ee_=var2
yH_=0.5*min(ee_/ee_ion,yHmax)
call rtsafe(ilnrho_ee,lnrho_,ee_,yHmin,yHmax*min(ee_/ee_ion,1.0),yH_)
fractions=(1+yH_+xHe)
TT_=(ee_-yH_*ee_ion)/(1.5*fractions*ss_ion)
lnTT_=log(TT_)
rho_=exp(lnrho_)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
pp_=fractions*rho_*TT_*ss_ion
cs2_=impossible
!
case (ilnrho_pp)
lnrho_=var1
pp_=var2
yH_=0.5*yHmax
call rtsafe(ilnrho_pp,lnrho_,pp_,yHmin,yHmax,yH_)
fractions=(1+yH_+xHe)
rho_=exp(lnrho_)
TT_=pp_/(fractions*ss_ion*rho_)
lnTT_=log(TT_)
ss_=ss_ion*(fractions*(1.5*(lnTT_-lnTT_ion)-lnrho_+2.5) &
-yH_*(2*log(yH_)-lnrho_e-lnrho_H) &
-(1-yH_)*(log(1-yH_+epsi)-lnrho_H)-xHe_term)
ee_=1.5*fractions*ss_ion*TT_+yH_*ee_ion
cs2_=impossible
!
case default
call stop_it("eoscalc_point: I don't get what the independent variables are.")
!
end select
!
if (present(lnrho)) lnrho=lnrho_
if (present(ss)) ss=ss_
if (present(yH)) yH=yH_
if (present(lnTT)) lnTT=lnTT_
if (present(ee)) ee=ee_
if (present(pp)) pp=pp_
if (present(cs2)) cs2=cs2_
!
endsubroutine eoscalc_point
!***********************************************************************
subroutine read_eos_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=eos_init_pars, IOSTAT=iostat)
!
endsubroutine read_eos_init_pars
!***********************************************************************
subroutine write_eos_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=eos_init_pars)
!
endsubroutine write_eos_init_pars
!***********************************************************************
subroutine read_eos_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=eos_run_pars, IOSTAT=iostat)
!
endsubroutine read_eos_run_pars
!***********************************************************************
subroutine write_eos_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=eos_run_pars)
!
endsubroutine write_eos_run_pars
!***********************************************************************
subroutine rtsafe_pencil(lnrho,ss,yH)
!
! safe newton raphson algorithm (adapted from NR) !
! 09-apr-03/tobi: changed to subroutine
!
real, dimension(mx), intent(in) :: lnrho,ss
real, dimension(mx), intent(inout) :: yH
!
real, dimension(mx) :: dyHold,dyH,yHlow,yHhigh,f,df
real, dimension(mx) :: lnTT_,dlnTT_,TT1_,fractions1
logical, dimension(mx) :: found
integer :: i
integer, parameter :: maxit=1000
!
yHlow=yHmin
yHhigh=yHmax
dyH=yHhigh-yHlow
dyHold=dyH
!
found=.false.
!
fractions1=1/(1+yH+xHe)
lnTT_=(2.0/3.0)*((ss/ss_ion+(1-yH)*(log(1-yH+epsi)-lnrho_H) &
+yH*(2*log(yH)-lnrho_e-lnrho_H) &
+xHe_term)*fractions1+lnrho-2.5)
TT1_=exp(-lnTT_)
f=lnrho_e-lnrho+1.5*lnTT_-TT1_+log(1-yH+epsi)-2*log(yH)
dlnTT_=((2.0/3.0)*(-f-TT1_)-1)*fractions1
! wd: Need to add epsi at the end, as otherwise (yH-yHlow)*df will
! wd: eventually yield an overflow.
! wd: Something like sqrt(tini) would probably also do instead of epsi,
! wd: but even epsi does not affect the auto-tests so far.
! df=dlnTT_*(1.5+TT1_)-1/(1-yH+epsi)-2/yH
df=dlnTT_*(1.5+TT1_)-1/(1-yH+epsi)-2/(yH+epsi)
!
do i=1,maxit
where (.not.found)
where ( sign(1.,((yH-yHlow)*df-f)) &
== sign(1.,((yH-yHhigh)*df-f)) &
.or. abs(2*f) > abs(dyHold*df) )
!
! Bisection
!
dyHold=dyH
dyH=0.5*(yHhigh-yHlow)
yH=yHhigh-dyH
elsewhere
!
! Newton-Raphson
!
dyHold=dyH
dyH=f/df
! Apply floor to dyH (necessary to avoid negative yH in samples
! /0d-tests/heating_ionize)
dyH=min(dyH,yH-yHmin)
dyH=max(dyH,yH-yHmax) ! plausibly needed as well
yH=yH-dyH
endwhere
endwhere
where (abs(dyH)>max(yHacc,1e-31)*max(yH,1e-31)) ! use max to avoid underflow
fractions1=1/(1+yH+xHe)
lnTT_=(2.0/3.0)*((ss/ss_ion+(1-yH)*(log(1-yH+epsi)-lnrho_H) &
+yH*(2*log(yH)-lnrho_e-lnrho_H) &
+xHe_term)*fractions1+lnrho-2.5)
TT1_=exp(-lnTT_)
f=lnrho_e-lnrho+1.5*lnTT_-TT1_+log(1-yH+epsi)-2*log(yH)
dlnTT_=((2.0/3.0)*(-f-TT1_)-1)*fractions1
df=dlnTT_*(1.5+TT1_)-1/(1-yH+epsi)-2/yH
where (f<0)
yHhigh=yH
elsewhere
yHlow=yH
endwhere
elsewhere
found=.true.
endwhere
if (all(found)) return
enddo
!
endsubroutine rtsafe_pencil
!***********************************************************************
subroutine rtsafe(ivars,var1,var2,yHlb,yHub,yH,rterror,rtdebug)
!
! safe newton raphson algorithm (adapted from NR) !
! 09-apr-03/tobi: changed to subroutine
!
integer, intent(in) :: ivars
real, intent(in) :: var1,var2
real, intent(in) :: yHlb,yHub
real, intent(inout) :: yH
logical, intent(out), optional :: rterror
logical, intent(in), optional :: rtdebug
!
real :: dyHold,dyH,yHl,yHh,f,df,temp
integer :: i
integer, parameter :: maxit=1000
!
if (present(rterror)) rterror=.false.
if (present(rtdebug)) then
if (rtdebug) print*,'rtsafe: i,yH=',0,yH
endif
!
yHl=yHlb
yHh=yHub
dyH=1
dyHold=dyH
!
call saha(ivars,var1,var2,yH,f,df)
!
do i=1,maxit
if (present(rtdebug)) then
if (rtdebug) print*,'rtsafe: i,yH=',i,yH
endif
if ( sign(1.,((yH-yHl)*df-f)) &
== sign(1.,((yH-yHh)*df-f)) &
.or. abs(2*f) > abs(dyHold*df) ) then
dyHold=dyH
dyH=0.5*(yHl-yHh)
yH=yHh+dyH
if (yHh==yH) return
else
dyHold=dyH
dyH=f/df
temp=yH
yH=yH-dyH
if (temp==yH) return
endif
if (abs(dyH)<yHacc*yH) return
call saha(ivars,var1,var2,yH,f,df)
if (f<0) then
yHh=yH
else
yHl=yH
endif
enddo
!
if (present(rterror)) rterror=.true.
!
endsubroutine rtsafe
!***********************************************************************
subroutine saha(ivars,var1,var2,yH,f,df)
!
! we want to find the root of f
!
! 23-feb-03/tobi: errors fixed
!
use Mpicomm, only: stop_it
!
integer, intent(in) :: ivars
real, intent(in) :: var1,var2,yH
real, intent(out) :: f,df
!
real :: lnrho,ss,ee,pp
real :: lnTT_,dlnTT_,TT1_,fractions1
!
fractions1=1/(1+yH+xHe)
!
select case (ivars)
case (ilnrho_ss)
lnrho=var1
ss=var2
lnTT_=(2.0/3.0)*((ss/ss_ion+(1-yH)*(log(1-yH+epsi)-lnrho_H) &
+yH*(2*log(yH)-lnrho_e-lnrho_H) &
+xHe_term)*fractions1+lnrho-2.5)
case (ilnrho_ee)
lnrho=var1
ee=var2
!print*,'saha: TT_',2.0/3.0*(ee-yH*ee_ion)*fractions1
!lnTT_=log(2.0/3.0*(ee/ee_ion-yH)*fractions1)
! if (ee<yH*ee_ion) then
! lnTT_=-25.
! else
lnTT_=log(2.0/3.0*(ee-yH*ee_ion)*fractions1)
! endif
case (ilnrho_pp)
lnrho=var1
pp=var2
lnTT_=log(pp/ss_ion*fractions1)-lnrho
case default
call stop_it("saha: I don't get what the independent variables are.")
lnTT_=0.
end select
!
TT1_=exp(-lnTT_)
f=lnrho_e-lnrho+1.5*lnTT_-TT1_+log(1-yH+epsi)-2*log(yH)
dlnTT_=((2.0/3.0)*(-f-TT1_)-1)*fractions1
df=dlnTT_*(1.5+TT1_)-1/(1-yH+epsi)-2/yH
!
endsubroutine saha
!***********************************************************************
subroutine get_soundspeed(TT,cs2)
!
! Calculate sound speed for given temperature
!
! 20-Oct-03/tobi: coded
!
use Mpicomm
!
real, intent(in) :: TT
real, intent(out) :: cs2
!
call stop_it("get_soundspeed: with ionization, lnrho needs to be known here")
!
call keep_compiler_quiet(TT,cs2)
!
endsubroutine get_soundspeed
!***********************************************************************
subroutine isothermal_entropy(f,T0)
!
! Isothermal stratification (for lnrho and ss)
! This routine should be independent of the gravity module used.
! When entropy is present, this module also initializes entropy.
!
! Sound speed (and hence Temperature), is
! initialised to the reference value:
! sound speed: cs^2_0 from start.in
! density: rho0 = exp(lnrho0)
!
! 11-jun-03/tony: extracted from isothermal routine in Density module
! to allow isothermal condition for arbitrary density
! 17-oct-03/nils: works also with leos_ionization=T
! 18-oct-03/tobi: distributed across ionization modules
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
real, intent(in) :: T0
real, dimension(nx) :: lnrho,ss,yH,K,sqrtK,yH_term,one_yH_term
!
do n=n1,n2
do m=m1,m2
!
lnrho=f(l1:l2,m,n,ilnrho)
!
K=exp(lnrho_e-lnrho-TT_ion/T0)*(T0/TT_ion)**1.5
sqrtK=sqrt(K)
yH=2*sqrtK/(sqrtK+sqrt(4+K))
!
where (yH>0)
yH_term=yH*(2*log(yH)-lnrho_e-lnrho_H)
elsewhere
yH_term=0
endwhere
!
where (yH<1)
one_yH_term=(1-yH)*(log(1-yH+epsi)-lnrho_H)
elsewhere
one_yH_term=0
endwhere
!
ss=ss_ion*((1+yH+xHe)*(1.5*log(T0/TT_ion)-lnrho+2.5) &
-yH_term-one_yH_term-xHe_term)
!
f(l1:l2,m,n,iss)=ss
!
enddo
enddo
!
endsubroutine isothermal_entropy
!***********************************************************************
subroutine isothermal_lnrho_ss(f,T0,rho0)
!
! Isothermal stratification for lnrho and ss (for yH=0!)
!
! Currently only implemented for ionization_fixed.
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
real, intent(in) :: T0,rho0
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(T0)
call keep_compiler_quiet(rho0)
!
endsubroutine isothermal_lnrho_ss
!***********************************************************************
subroutine get_average_pressure(average_density,average_pressure)
!
! 01-dec-2009/piyali+dhruba: coded
!
use Cdata
!
real, intent(in):: average_density
real, intent(out):: average_pressure
!
call keep_compiler_quiet(average_density)
call keep_compiler_quiet(average_pressure)
!
endsubroutine get_average_pressure
!***********************************************************************
subroutine bc_ss_flux(f,topbot,lone_sided)
!
! constant flux boundary condition for entropy (called when bcz='c1')
!
! 23-jan-2002/wolf: coded
! 11-jun-2002/axel: moved into the entropy module
! 8-jul-2002/axel: split old bc_ss into two
! 26-aug-2003/tony: distributed across ionization modules
! 3-oct-16/MR: added new optional switch lone_sided
!
use Mpicomm, only: stop_it
use Gravity
use SharedVariables,only:get_shared_variable
!
real, pointer :: Fbot,Ftop,FtopKtop,FbotKbot,hcond0,hcond1,chi
logical, pointer :: lmultilayer, lheatc_chiconst
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
logical, optional :: lone_sided
real, dimension (size(f,1),size(f,2)) :: tmp_xy,TT_xy,rho_xy,yH_xy
integer :: i
!
! Do the `c1' boundary condition (constant heat flux) for entropy.
! check whether we want to do top or bottom (this is precessor dependent)
!
! Get the shared variables
!
call get_shared_variable('hcond0',hcond0,caller='bc_ss_flux')
call get_shared_variable('hcond1',hcond1)
call get_shared_variable('Fbot',Fbot)
call get_shared_variable('Ftop',Ftop)
call get_shared_variable('FbotKbot',FbotKbot)
call get_shared_variable('FtopKtop',FtopKtop)
call get_shared_variable('chi',chi)
call get_shared_variable('lmultilayer',lmultilayer)
call get_shared_variable('lheatc_chiconst',lheatc_chiconst)
!
select case (topbot)
!
! bottom boundary
! ---------------
!
case ('bot')
if (lmultilayer) then
if (headtt) print*,'bc_ss_flux: Fbot,hcond=',Fbot,hcond0*hcond1
else
if (headtt) print*,'bc_ss_flux: Fbot,hcond=',Fbot,hcond0
endif
!
! calculate Fbot/(K*cs2)
!
rho_xy=exp(f(:,:,n1,ilnrho))
TT_xy=exp(f(:,:,n1,ilnTT))
!
! check whether we have chi=constant at bottom, in which case
! we have the nonconstant rho_xy*chi in tmp_xy.
!
if (lheatc_chiconst) then
tmp_xy=Fbot/(rho_xy*chi*TT_xy)
else
tmp_xy=FbotKbot/TT_xy
endif
!
! get ionization fraction at bottom boundary
!
yH_xy=f(:,:,n1,iyH)
!
! enforce ds/dz + gamma_m1/gamma*dlnrho/dz = - gamma_m1/gamma*Fbot/(K*cs2)
!
do i=1,nghost
f(:,:,n1-i,iss)=f(:,:,n1+i,iss)+ss_ion*(1+yH_xy+xHe)* &
(f(:,:,n1+i,ilnrho)-f(:,:,n1-i,ilnrho)+3*i*dz*tmp_xy)
enddo
!
! top boundary
! ------------
!
case ('top')
if (lmultilayer) then
if (headtt) print*,'bc_ss_flux: Ftop,hcond=',Ftop,hcond0*hcond1
else
if (headtt) print*,'bc_ss_flux: Ftop,hcond=',Ftop,hcond0
endif
!
! calculate Fbot/(K*cs2)
!
rho_xy=exp(f(:,:,n2,ilnrho))
TT_xy=exp(f(:,:,n2,ilnTT))
!
! check whether we have chi=constant at bottom, in which case
! we have the nonconstant rho_xy*chi in tmp_xy.
!
if (lheatc_chiconst) then
tmp_xy=Ftop/(rho_xy*chi*TT_xy)
else
tmp_xy=FtopKtop/TT_xy
endif
!
! get ionization fraction at top boundary
!
yH_xy=f(:,:,n2,iyH)
!
! enforce ds/dz + gamma_m1/gamma*dlnrho/dz = - gamma_m1/gamma*Fbot/(K*cs2)
!
do i=1,nghost
f(:,:,n2+i,iss)=f(:,:,n2-i,iss)+ss_ion*(1+yH_xy+xHe)* &
(f(:,:,n2-i,ilnrho)-f(:,:,n2+i,ilnrho)-3*i*dz*tmp_xy)
enddo
!
case default
print*,"bc_ss_flux: invalid argument"
call stop_it("")
endselect
call keep_compiler_quiet(present(lone_sided))
!
endsubroutine bc_ss_flux
!***********************************************************************
subroutine bc_ss_flux_turb(f,topbot)
!
! dummy routine
!
! 4-may-2009/axel: dummy routine
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_flux_turb
!***********************************************************************
subroutine bc_ss_flux_turb_x(f,topbot)
!
! dummy routine
!
! 31-may-2010/pete: dummy routine
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_flux_turb_x
!***********************************************************************
subroutine bc_ss_flux_condturb_x(f,topbot)
!
! 23-apr-2014/pete: dummy
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_flux_condturb_x
!***********************************************************************
subroutine bc_ss_flux_condturb_mean_x(f,topbot)
!
! 07-jan-2015/pete: dummy
!
character (len=3) :: topbot
real, dimension (mx,my,mz,mfarray) :: f
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_flux_condturb_mean_x
!***********************************************************************
subroutine bc_ss_flux_condturb_z(f,topbot)
!
! 15-jul-2014/pete: dummy
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_flux_condturb_z
!***********************************************************************
subroutine bc_ss_temp_old(f,topbot)
!
! boundary condition for entropy: constant temperature
!
! 23-jan-2002/wolf: coded
! 11-jun-2002/axel: moved into the entropy module
! 8-jul-2002/axel: split old bc_ss into two
! 23-jun-2003/tony: implemented for leos_fixed_ionization
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_temp_old: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_temp_old
!***********************************************************************
subroutine bc_ss_temp_x(f,topbot)
!
! boundary condition for entropy: constant temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_temp_x: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_temp_x
!***********************************************************************
subroutine bc_ss_temp_y(f,topbot)
!
! boundary condition for entropy: constant temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_temp_y: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_temp_y
!***********************************************************************
subroutine bc_ss_temp_z(f,topbot,lone_sided)
!
! boundary condition for entropy: constant temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
logical, optional :: lone_sided
!
call stop_it("bc_ss_temp_z: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
call keep_compiler_quiet(present(lone_sided))
!
endsubroutine bc_ss_temp_z
!***********************************************************************
subroutine bc_lnrho_temp_z(f,topbot)
!
! boundary condition for density: constant temperature
!
! 19-aug-2005/tobi: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_lnrho_temp_z: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_lnrho_temp_z
!***********************************************************************
subroutine bc_lnrho_pressure_z(f,topbot)
!
! boundary condition for density: constant pressure
!
! 19-aug-2005/tobi: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_lnrho_pressure_z: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_lnrho_pressure_z
!***********************************************************************
subroutine bc_ss_temp2_z(f,topbot)
!
! boundary condition for entropy: constant temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_temp2_z: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_temp2_z
!***********************************************************************
subroutine bc_ss_temp3_z(f,topbot)
!
! 31-jan-2013/axel: coded to impose cs2bot and dcs2bot at bottom
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call fatal_error('bc_ss_temp3_z', &
'not implemented in eos_ionization')
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_temp3_z
!***********************************************************************
subroutine bc_ss_stemp_x(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_stemp_x: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_stemp_x
!***********************************************************************
subroutine bc_ss_stemp_y(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_stemp_y: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_stemp_y
!***********************************************************************
subroutine bc_ss_stemp_z(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 26-sep-2003/tony: coded
!
use Mpicomm, only: stop_it
use Gravity
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
real, dimension (size(f,1),size(f,2),nghost) :: lnrho,ss,yH,lnTT,TT,K,sqrtK,yH_term,one_yH_term
integer :: i
!
! Symmetric temperature/sound speed for entropy.
! This assumes that the density is already set (ie density _must_ register
! first!)
!
! check whether we want to do top or bottom (this is processor dependent)
!
select case (topbot)
!
! bottom boundary
!
case ('bot')
do i=1,nghost
f(:,:,n1-i,ilnTT) = f(:,:,n1+i,ilnTT)
enddo
!
lnrho=f(:,:,1:n1-1,ilnrho)
lnTT=f(:,:,1:n1-1,ilnTT)
TT=exp(lnTT)
!
K=exp(lnrho_e-lnrho-TT_ion/TT)*(TT/TT_ion)**1.5
sqrtK=sqrt(K)
yH=2*sqrtK/(sqrtK+sqrt(4+K))
!
where (yH>0)
yH_term=yH*(2*log(yH)-lnrho_e-lnrho_H)
elsewhere
yH_term=0
endwhere
!
where (yH<1)
one_yH_term=(1-yH)*(log(1-yH+epsi)-lnrho_H)
elsewhere
one_yH_term=0
endwhere
!
ss=ss_ion*((1+yH+xHe)*(1.5*log(TT/TT_ion)-lnrho+2.5) &
-yH_term-one_yH_term-xHe_term)
!
f(:,:,1:n1-1,iyH)=yH
f(:,:,1:n1-1,iss)=ss
!
! top boundary
!
case ('top')
do i=1,nghost
f(:,:,n2+i,ilnTT) = f(:,:,n2-i,ilnTT)
enddo
!
lnrho=f(:,:,n2+1:,ilnrho)
lnTT =f(:,:,n2+1:,ilnTT)
TT=exp(lnTT)
!
K=exp(lnrho_e-lnrho-TT_ion/TT)*(TT/TT_ion)**1.5
sqrtK=sqrt(K)
yH=2*sqrtK/(sqrtK+sqrt(4+K))
!
where (yH>0)
yH_term=yH*(2*log(yH)-lnrho_e-lnrho_H)
elsewhere
yH_term=0
endwhere
!
where (yH<1)
one_yH_term=(1-yH)*(log(1-yH+epsi)-lnrho_H)
elsewhere
one_yH_term=0
endwhere
!
ss=ss_ion*((1+yH+xHe)*(1.5*log(TT/TT_ion)-lnrho+2.5) &
-yH_term-one_yH_term-xHe_term)
!
f(:,:,n2+1:,iyH)=yH
f(:,:,n2+1:,iss)=ss
!
case default
print*,"bc_ss_stemp_z: invalid argument"
call stop_it("")
endselect
!
endsubroutine bc_ss_stemp_z
!***********************************************************************
subroutine bc_ss_a2stemp_x(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_a2stemp_x: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_a2stemp_x
!***********************************************************************
subroutine bc_ss_a2stemp_y(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_a2stemp_y: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_a2stemp_y
!***********************************************************************
subroutine bc_ss_a2stemp_z(f,topbot)
!
! boundary condition for entropy: symmetric temperature
!
! 3-aug-2002/wolf: coded
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_ss_a2stemp_z: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_a2stemp_z
!***********************************************************************
subroutine bc_ss_energy(f,topbot)
!
! boundary condition for entropy
!
! may-2002/nils: coded
! 11-jul-2002/nils: moved into the entropy module
! 26-aug-2003/tony: distributed across ionization modules
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
! The 'ce' boundary condition for entropy makes the energy constant at
! the boundaries.
! This assumes that the density is already set (ie density must register
! first!)
!
call stop_it("bc_ss_stemp_y: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_ss_energy
!***********************************************************************
subroutine bc_stellar_surface(f,topbot)
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_stellar_surface: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_stellar_surface
!***********************************************************************
subroutine bc_lnrho_cfb_r_iso(f,topbot)
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_lnrho_cfb_r_iso: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_lnrho_cfb_r_iso
!***********************************************************************
subroutine bc_lnrho_hds_z_iso(f,topbot)
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_lnrho_hds_z_iso: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_lnrho_hds_z_iso
!***********************************************************************
subroutine bc_lnrho_hdss_z_iso(f,topbot)
!
use Mpicomm, only: stop_it
!
character (len=3) :: topbot
real, dimension (:,:,:,:) :: f
!
call stop_it("bc_lnrho_hdss_z_iso: NOT IMPLEMENTED IN EOS_IONIZATION")
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(topbot)
!
endsubroutine bc_lnrho_hdss_z_iso
!***********************************************************************
subroutine write_thermodyn
!
endsubroutine write_thermodyn
!***********************************************************************
subroutine read_thermodyn(input_file)
!
character (len=*), intent(in) :: input_file
!
call keep_compiler_quiet(input_file)
!
endsubroutine read_thermodyn
!***********************************************************************
subroutine read_species(input_file)
!
character (len=*) :: input_file
!
call keep_compiler_quiet(input_file)
!
endsubroutine read_species
!***********************************************************************
subroutine find_species_index(species_name,ind_glob,ind_chem,found_specie)
!
integer, intent(out) :: ind_glob
integer, intent(inout) :: ind_chem
character (len=*), intent(in) :: species_name
logical, intent(out) :: found_specie
!
call keep_compiler_quiet(ind_glob)
call keep_compiler_quiet(ind_chem)
call keep_compiler_quiet(species_name)
call keep_compiler_quiet(found_specie)
!
endsubroutine find_species_index
!***********************************************************************
subroutine find_mass(element_name,MolMass)
!
character (len=*), intent(in) :: element_name
real, intent(out) :: MolMass
!
call keep_compiler_quiet(element_name)
call keep_compiler_quiet(MolMass)
!
endsubroutine find_mass
!***********************************************************************
subroutine get_stratz(z, rho0z, dlnrho0dz, eth0z)
!
! Get background stratification in z direction.
!
! 13-oct-14/ccyang: dummy
!
real, dimension(:), intent(in) :: z
real, dimension(:), intent(out), optional :: rho0z, dlnrho0dz, eth0z
!
call fatal_error('get_stratz', 'Stratification for this EOS is not implemented. ')
!
call keep_compiler_quiet(z)
if (present(rho0z)) call keep_compiler_quiet(rho0z)
if (present(dlnrho0dz)) call keep_compiler_quiet(dlnrho0dz)
if (present(eth0z)) call keep_compiler_quiet(eth0z)
!
endsubroutine get_stratz
!***********************************************************************
subroutine pushdiags2c(p_diag)
integer, parameter :: n_diags=0
integer(KIND=ikind8), dimension(:) :: p_diag
call keep_compiler_quiet(p_diag)
endsubroutine pushdiags2c
!***********************************************************************
subroutine pushpars2c(p_par)
use Syscalls, only: copy_addr
integer, parameter :: n_pars=1
integer(KIND=ikind8), dimension(n_pars) :: p_par
call copy_addr(cs20,p_par(1))
endsubroutine pushpars2c
!***********************************************************************
endmodule EquationOfState