! $Id$ ! ! Equation of state for an ideal gas without 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 0 ! ! PENCILS PROVIDED lnTT; glnTT(3); TT; TT1; gTT(3) ! PENCILS PROVIDED pp; del2pp; mu1; gmu1(3); glnmu(3) ! PENCILS PROVIDED rho1gpp(3); glnpp(3); del2lnTT ! ! PENCILS PROVIDED hss(3,3); hlnTT(3,3); del2ss; del6ss; del6lnTT ! PENCILS PROVIDED yH; ee; ss; delta; glnmumol(3); ppvap; csvap2; cs2 ! PENCILS PROVIDED cp1tilde; cp; gamma_m1; gamma ! PENCILS PROVIDED rho_anel; gradcp(3) ! ! !*************************************************************** module EquationOfState ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! include 'eos.h' ! integer, parameter :: ilnrho_ss=1,ilnrho_ee=2,ilnrho_pp=3 integer, parameter :: ilnrho_lnTT=4,ilnrho_cs2=5 integer, parameter :: irho_cs2=6, irho_ss=7, irho_lnTT=8, ilnrho_TT=9 integer, parameter :: ipp_ss=11, irho_TT=10, ipp_cs2=12 integer, parameter :: irho_eth=13, ilnrho_eth=14 ! integer :: iglobal_cs2, iglobal_glnTT, ics ! real :: lnTT0=impossible ! real :: mu=1. real :: cs0=1., rho0=1. real :: cs20=1., lnrho0=0. real :: gamma=5./3. real :: Rgas_cgs=0., Rgas, Rgas_unit_sys=1., error_cp=1e-6 real :: gamma_m1 !(=gamma-1) real :: gamma1 !(=1/gamma) real :: cp=impossible, cp1=impossible, cv=impossible, cv1=impossible real :: cs2top_ini=impossible, dcs2top_ini=impossible real :: cs2bot=1., cs2top=1. real :: cs2cool=0. real :: mpoly=1.5, mpoly0=1.5, mpoly1=1.5, mpoly2=1.5 real, dimension(3) :: beta_glnrho_global=0., beta_glnrho_scaled=0. integer :: isothtop=0 integer :: ieosvars=-1, ieosvar1=-1, ieosvar2=-1, ieosvar_count=0 integer :: ll1,ll2,mm1,mm2,nn1,nn2 logical :: leos_isothermal=.false., leos_isentropic=.false. logical :: leos_isochoric=.false., leos_isobaric=.false. logical :: leos_localisothermal=.false. character (len=20) :: input_file logical, SAVE :: lcheminp_eos=.false. logical :: l_gamma_m1=.false. logical :: l_gamma=.false. logical :: l_cp=.false. integer :: imass=1!, iTemp1=2,iTemp2=3,iTemp3=4 ! character (len=labellen) :: ieos_profile='nothing' real, dimension(mz) :: profz_eos=1.,dprofz_eos=0. ! real, dimension(nchemspec,18) :: species_constants real, dimension(nchemspec,7) :: tran_data real, dimension(nchemspec) :: Lewis_coef, Lewis_coef1 ! ! !NILS: Why do we spend a lot of memory allocating these variables here???? real, dimension (mx,my,mz), SAVE :: mu1_full ! namelist /eos_init_pars/ mu, cp, cs0, rho0, gamma, error_cp ! namelist /eos_run_pars/ mu, cp, cs0, rho0, gamma, error_cp ! contains !*********************************************************************** subroutine register_eos ! ! 14-jun-03/axel: adapted from register_eos ! leos_chemistry=.true. ! ilnTT = 0 ! ! Identify version number. ! if (lroot) call svn_id( & '$Id$') ! endsubroutine register_eos !*********************************************************************** subroutine units_eos ! ! This routine calculates things related to units and must be called ! before the rest of the units are being calculated. ! ! 22-jun-06/axel: adapted from initialize_eos ! 16-mar-10/Natalia ! use Mpicomm, only: stop_it ! logical :: chemin=.false.,cheminp=.false. ! ! Initialize variable selection code (needed for RELOADing) ! ieosvars=-1 ieosvar_count=0 ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs elseif (unit_system == 'SI') then Rgas_unit_sys = k_B_cgs/m_u_cgs*1.e-4 endif ! if (unit_temperature == impossible) then call stop_it('unit_temperature is not found!') else Rgas=Rgas_unit_sys*unit_temperature/unit_velocity**2 endif ! inquire(file='chem.inp',exist=cheminp) inquire(file='chem.in',exist=chemin) if(chemin .and. cheminp) call fatal_error('eos_chemistry',& 'chem.inp and chem.in found. Please decide for one') ! if (cheminp) input_file='chem.inp' if (chemin) input_file='chem.in' lcheminp_eos = cheminp .or. chemin ! inquire(FILE=input_file, EXIST=lcheminp_eos) ! if (lroot) then ! if (.not. lcheminp_eos ) then call fatal_error('initialize_eos',& 'chem.imp is not found!') else print*,'initialize_eos: chem.imp is found! Now cp, cv, gamma, mu are pencils ONLY!' endif endif ! endsubroutine units_eos !*********************************************************************** subroutine initialize_eos ! use SharedVariables, only: put_shared_variable ! ! Initialize variable selection code (needed for RELOADing) ! ieosvars=-1 ieosvar_count=0 ! ! write constants to disk. In future we may want to deal with this ! using an include file or another module. ! if (lroot) then open (1,file=trim(datadir)//'/pc_constants.pro',position="append") write (1,'(a,1pd26.16)') 'k_B=',k_B write (1,'(a,1pd26.16)') 'm_H=',m_H write (1,*) 'lnTTO=',lnTT0 write (1,*) 'cp=',cp close (1) endif ! if ((nxgrid==1) .and. (nygrid==1) .and. (nzgrid==1)) then ll1=1; ll2=mx; mm1=m1; mm2=m2; nn1=n1; nn2=n2 else if (nxgrid==1) then ll1=l1; ll2=l2 else ll1=1; ll2=mx endif ! if (nygrid==1) then mm1=m1; mm2=m2 else mm1=1; mm2=my endif ! if (nzgrid==1) then nn1=n1; nn2=n2 else nn1=1; nn2=mz endif endif if (.not.ldensity) then call put_shared_variable('rho0',rho0,caller='initialize_eos') call put_shared_variable('lnrho0',lnrho0) endif ! endsubroutine initialize_eos !*********************************************************************** subroutine select_eos_variable(variable,findex) ! ! Select eos variable ! ! 02-apr-06/tony: implemented ! use FArrayManager ! character (len=*), intent(in) :: variable integer, intent(in) :: findex integer :: this_var=0 integer, save :: ieosvar_selected=0 integer, parameter :: ieosvar_lnrho = 2**0 integer, parameter :: ieosvar_rho = 2**1 integer, parameter :: ieosvar_ss = 2**2 integer, parameter :: ieosvar_lnTT = 2**3 integer, parameter :: ieosvar_TT = 2**4 integer, parameter :: ieosvar_cs2 = 2**5 integer, parameter :: ieosvar_pp = 2**6 ! if (ieosvar_count==0) ieosvar_selected=0 ! if (ieosvar_count>=2) & call fatal_error("select_eos_variable", & "2 thermodynamic quantities have already been defined while attempting to add a 3rd: ") !//variable) ! ieosvar_count=ieosvar_count+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. ! call farray_register_global('cs2',iglobal_cs2) call farray_register_global('glnTT',iglobal_glnTT,vector=3) ! 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=='TT') then this_var=ieosvar_TT 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_count==1) then ieosvar1=findex ieosvar_selected=ieosvar_selected+this_var return endif ! ! Ensure the indexes are in the correct order. ! if (this_var0.) then do j2=mm1,mm2 do j3=nn1,nn2 mu1_full_tmp(:,j2,j3)= & mu1_full_tmp(:,j2,j3)+unit_mass*f(:,j2,j3,ichemspec(k)) & /species_constants(k,imass) enddo enddo endif enddo mu1_full=mu1_full_tmp ! endsubroutine getmu_array !*********************************************************************** subroutine rprint_eos(lreset,lwrite) ! ! Writes iyH and ilnTT to index.pro file ! ! 14-jun-03/axel: adapted from rprint_radiation ! 21-11-04/anders: moved diagnostics to entropy ! logical :: lreset logical, optional :: lwrite ! call keep_compiler_quiet(lreset) call keep_compiler_quiet(present(lwrite)) ! endsubroutine rprint_eos !*********************************************************************** subroutine get_slices_eos(f,slices) ! ! Write slices for animation of Eos variables. ! ! 26-jul-06/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f type (slice_data) :: slices ! ! Loop over slices. ! select case (trim(slices%name)) ! ! Temperature. ! case ('lnTT') slices%yz =f(ix_loc,m1:m2,n1:n2,ilnTT) slices%xz =f(l1:l2,iy_loc,n1:n2,ilnTT) slices%xy =f(l1:l2,m1:m2,iz_loc,ilnTT) slices%xy2=f(l1:l2,m1:m2,iz2_loc,ilnTT) if (lwrite_slice_xy3) slices%xy3=f(l1:l2,m1:m2,iz3_loc,ilnTT) if (lwrite_slice_xy4) slices%xy4=f(l1:l2,m1:m2,iz4_loc,ilnTT) slices%ready=.true. case ('pp') if (ldensity_nolog .or. ltemperature_nolog) & call fatal_error('get_slices_eos',& 'pp not implemented for ldensity_nolot .or. ltemperature_nolog') slices%yz =Rgas*exp(f(ix_loc,m1:m2,n1:n2,ilnTT)+f(ix_loc,m1:m2,n1:n2,ilnrho))& *mu1_full(ix_loc,m1:m2,n1:n2) slices%xz =Rgas*exp(f(l1:l2,iy_loc,n1:n2,ilnTT)+f(l1:l2,iy_loc,n1:n2,ilnrho))& *mu1_full(l1:l2,iy_loc,n1:n2) slices%xy =Rgas*exp(f(l1:l2,m1:m2,iz_loc,ilnTT)+f(l1:l2,m1:m2,iz_loc,ilnrho))& *mu1_full(l1:l2,m1:m2,iz_loc) if (lwrite_slice_xy4) slices%xy3 =Rgas*exp(f(l1:l2,m1:m2,iz3_loc,ilnTT)& +f(l1:l2,m1:m2,iz3_loc,ilnrho))*mu1_full(l1:l2,m1:m2,iz3_loc) if (lwrite_slice_xy3) slices%xy4 =Rgas*exp(f(l1:l2,m1:m2,iz4_loc,ilnTT)& +f(l1:l2,m1:m2,iz4_loc,ilnrho))*mu1_full(l1:l2,m1:m2,iz4_loc) slices%ready=.true. ! endselect ! endsubroutine get_slices_eos !*********************************************************************** subroutine pencil_criteria_eos ! ! All pencils that the EquationOfState module depends on are specified here. ! ! 02-04-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. ! lpenc_requested(i_TT)=.true. lpenc_requested(i_TT1)=.true. if (.not. ldustdensity) then lpenc_requested(i_lnTT)=.true. endif lpenc_requested(i_del2lnTT)=.true. ! if (ltemperature_nolog) then lpenc_requested(i_gTT)=.true. endif lpenc_requested(i_glnTT)=.true. lpenc_requested(i_glnrho)=.true. lpenc_requested(i_glnrho2)=.true. lpenc_requested(i_del2lnrho)=.true. ! lpenc_requested(i_glnpp)=.true. lpenc_requested(i_del2pp)=.true. lpenc_requested(i_mu1)=.true. lpenc_requested(i_gmu1)=.true. lpenc_requested(i_pp)=.true. lpenc_requested(i_rho1gpp)=.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 ! ! Modified by Natalia. Taken from eos_temperature_ionization module ! logical, dimension(npencils) :: lpencil_in ! if (lpencil_in(i_TT1)) lpencil_in(i_TT)=.true. ! if (lpencil_in(i_pp)) then lpencil_in(i_mu1)=.true. lpencil_in(i_rho)=.true. lpencil_in(i_TT)=.true. endif if (lpencil_in(i_rho1gpp)) then lpencil_in(i_mu1)=.true. lpencil_in(i_gmu1)=.true. lpencil_in(i_pp)=.true. lpencil_in(i_rho)=.true. lpencil_in(i_glnTT)=.true. endif if (lpencil_in(i_glnpp)) then lpencil_in(i_pp)=.true. lpencil_in(i_rho)=.true. lpencil_in(i_rho1gpp)=.true. endif ! if (lpencil_in(i_ee)) then lpencil_in(i_mu1)=.true. lpencil_in(i_TT)=.true. endif ! call keep_compiler_quiet(lpencil_in) ! 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(OUT) :: 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, only: grad, del2, dot2 ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p logical, dimension(npencils) :: lpenc_loc real, dimension(nx) :: glnTT2,TT1del2TT,del2lnrho real, dimension(nx) :: rho1del2rho,del2mu1,gmu12,glnpp2 ! intent(in) :: f, lpenc_loc intent(inout) :: p ! integer :: i ! ! Temperature ! if (lpenc_loc(i_lnTT)) then if (ltemperature_nolog) then p%lnTT=log(f(l1:l2,m,n,iTT)) else p%lnTT=f(l1:l2,m,n,ilnTT) endif endif ! if (lpenc_loc(i_TT)) then if (ltemperature_nolog) then p%TT=f(l1:l2,m,n,iTT) else p%TT=exp(f(l1:l2,m,n,ilnTT)) endif ! if (minval(p%TT)==0.) then call fatal_error('calc_pencils_eos','p%TT=0!') endif endif ! if (lpenc_loc(i_TT1)) then if (ltemperature_nolog) then p%TT1=1./f(l1:l2,m,n,iTT) else p%TT1=1./exp(f(l1:l2,m,n,ilnTT)) endif endif ! ! Temperature laplacian and gradient ! if (lpenc_loc(i_glnTT)) then if (ltemperature_nolog) then call grad(f,iTT,p%glnTT) p%glnTT(:,1)=p%glnTT(:,1)/p%TT(:) p%glnTT(:,2)=p%glnTT(:,2)/p%TT(:) p%glnTT(:,3)=p%glnTT(:,3)/p%TT(:) else call grad(f,ilnTT,p%glnTT) endif endif ! if (ltemperature_nolog) then if (lpenc_loc(i_gTT)) call grad(f,iTT,p%gTT) !NILS: The call below does not yield del2lnTT but rather del2TT, !NILS: this should be fixed before used. One should also look !NILS: through the chemistry module to make sure del2lnTT is used !NILS: corectly. if (lpenc_loc(i_del2lnTT)) call del2(f,iTT,p%del2lnTT) call fatal_error('calc_pencils_eos',& 'del2lnTT is not correctly implemented - this must be fixed!') else if (lpenc_loc(i_del2lnTT)) call del2(f,ilnTT,p%del2lnTT) endif ! ! Mean molecular weight ! if (lpenc_loc(i_mu1)) p%mu1=mu1_full(l1:l2,m,n) if (lpenc_loc(i_gmu1)) call grad(mu1_full,p%gmu1) if (lpenc_loc(i_glnmu)) then do i = 1,3 p%glnmu(:,i) = -p%gmu1(:,i)/p%mu1(:) enddo endif ! ! Pressure ! if (lpenc_loc(i_pp)) p%pp = Rgas*p%TT*p%mu1*p%rho ! ! Logarithmic pressure gradient ! if (lpenc_loc(i_rho1gpp)) then do i=1,3 p%rho1gpp(:,i) = p%pp/p%rho(:) & *(p%glnrho(:,i)+p%glnTT(:,i)+p%gmu1(:,i)/p%mu1(:)) enddo endif ! ! Gradient of lnpp ! if (lpenc_loc(i_glnpp)) then do i=1,3 p%glnpp(:,i)=p%rho1gpp(:,i)*p%rho(:)/p%pp(:) enddo endif ! ! Laplasian of pressure ! if (lpenc_loc(i_del2pp)) then call dot2(p%glnTT,glnTT2) TT1del2TT=p%del2lnTT+glnTT2 rho1del2rho=p%del2lnrho+p%glnrho2 call dot2(p%gmu1,gmu12) call dot2(p%glnpp,glnpp2) call del2(mu1_full,del2mu1) p%del2pp& =p%pp*glnpp2& +p%pp*(rho1del2rho+TT1del2TT+del2mu1/p%mu1)& -p%pp*(p%glnrho2+glnTT2+gmu12/p%mu1**2) endif ! ! Energy per unit mass (this has been moved to chemistry.f90 in order ! to get the correct cv. ! !if (lpenc_loc(i_ee)) p%ee = p%cv*p%TT ! ! endif ! endsubroutine calc_pencils_eos_pencpar !*********************************************************************** 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 ! call keep_compiler_quiet(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 ! call keep_compiler_quiet(f) ! endsubroutine ioncalc !*********************************************************************** subroutine getdensity(f,EE,TT,yH,rho_full) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz), intent(out) :: rho_full real, intent(in), optional :: EE,TT,yH ! if (ldensity_nolog) then rho_full=f(:,:,:,ilnrho) else rho_full=exp(f(:,:,:,ilnrho)) endif ! call keep_compiler_quiet(present(yH)) call keep_compiler_quiet(present(EE)) call keep_compiler_quiet(present(TT)) ! endsubroutine getdensity !*********************************************************************** subroutine gettemperature(f,TT_full) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz), intent(out) :: TT_full ! if (ltemperature_nolog) then TT_full=f(:,:,:,ilnTT) else TT_full=exp(f(:,:,:,ilnTT)) endif ! endsubroutine gettemperature !*********************************************************************** subroutine getpressure(pp,TT,rho,mu1) ! real, dimension (nx), intent(out) :: pp real, dimension (nx), intent(in) :: TT, rho, mu1 integer :: j2,j3 ! pp=Rgas*mu1*rho*TT ! 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_ call fatal_error('get_cp1','SHOULD NOT BE CALLED WITH eos_chemistry') 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_ call fatal_error('get_cv1','SHOULD NOT BE CALLED WITH eos_chemistry') 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 ! cs2=impossible cp1tilde=impossible call fatal_error('pressure_gradient_farray','SHOULD NOT BE CALLED WITH eos_chemistry') ! call keep_compiler_quiet(f) ! 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 ! call fatal_error('pressure_gradient_farray','SHOULD NOT BE CALLED WITH eos_chemistry') ! call keep_compiler_quiet(cs2,cp1tilde,ss,lnrho) ! endsubroutine pressure_gradient_point !*********************************************************************** subroutine temperature_gradient(f,glnrho,gss,glnTT) ! ! Calculate thermodynamical quantities ! and optionally glnPP and glnTT ! gP/rho=cs2*(glnrho+cp1*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 ! call fatal_error('temperature_gradien','SHOULD NOT BE CALLED WITH eos_chemistry') ! ! given that we just stopped, it cannot become worse by setting ! cp1tilde to impossible, which allows the compiler to compile. ! call keep_compiler_quiet(f) call keep_compiler_quiet(glnrho,gss,glnTT) ! endsubroutine temperature_gradient !*********************************************************************** subroutine temperature_laplacian(f,del2lnrho,del2ss,del2lnTT) ! ! Calculate thermodynamical quantities ! and optionally glnPP and glnTT ! gP/rho=cs2*(glnrho+cp1*gss) ! ! 17-nov-03/tobi: adapted from subroutine eoscalc ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(nx), intent(in) :: del2lnrho,del2ss real, dimension(nx), intent(out) :: del2lnTT ! call fatal_error('temperature_laplacian','SHould not be called!') ! call keep_compiler_quiet(f) call keep_compiler_quiet(del2lnrho,del2ss,del2lnTT) ! endsubroutine temperature_laplacian !*********************************************************************** subroutine temperature_hessian(f,hlnrho,hss,hlnTT) ! ! Calculate thermodynamical quantities, cs2 and cp1 ! and optionally hlnPP and hlnTT ! hP/rho=cs2*(hlnrho+cp1*hss) ! ! 17-nov-03/tobi: adapted from subroutine eoscalc ! 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 ! call fatal_error('temperature_hessian', & 'This routine is not coded for eos_chemistry') ! hlnTT(:,:,:)=0 ! call keep_compiler_quiet(f) call keep_compiler_quiet(hss) call keep_compiler_quiet(hlnrho) ! endsubroutine temperature_hessian !*********************************************************************** subroutine eosperturb(f,psize,ee,pp,ss) ! ! Set f(l1:l2,m,n,iss), depending on the valyes of ee and pp ! Adding pressure perturbations is not implemented ! real, dimension(mx,my,mz,mfarray), intent(inout) :: f integer, intent(in) :: psize real, dimension(psize), intent(in), optional :: ee, pp, ss real, dimension(psize) :: lnrho_ ! if (psize==nx) then lnrho_=f(l1:l2,m,n,ilnrho) elseif (psize==mx) then lnrho_=f(:,m,n,ilnrho) else call not_implemented("eosperturb") endif ! call fatal_error('eosperturb', & 'This routine is not coded for eos_chemistry') ! call keep_compiler_quiet(present(ee)) call keep_compiler_quiet(present(pp)) call keep_compiler_quiet(present(ss)) ! endsubroutine eosperturb !*********************************************************************** subroutine eoscalc_farray(f,psize,lnrho,yH,lnTT,ee,pp,cs2,kapparho) ! ! dummy routine to calculate thermodynamical quantities ! copied from eo_idealgas ! real, dimension(mx,my,mz,mfarray), intent(in) :: f integer, intent(in) :: psize real, dimension(psize), optional :: lnrho,lnTT real, dimension(psize), optional :: yH,ee,pp,cs2,kapparho ! call keep_compiler_quiet(f) call keep_compiler_quiet(present(lnrho),present(lnTT)) call keep_compiler_quiet(present(yH),present(ee)) call keep_compiler_quiet(present(pp),present(kapparho)) call keep_compiler_quiet(present(cs2)) ! call fatal_error('eoscalc_farray', & 'This routine is not coded for eos_chemistry') ! endsubroutine eoscalc_farray !*********************************************************************** subroutine eoscalc_point(ivars,var1,var2,lnrho,ss,yH,lnTT,ee,pp,cs2) ! ! Calculate thermodynamical quantities ! ! 22-jun-06/axel: reinstated cp,cp1,cv,cv1 in hopefully all the places. ! integer, intent(in) :: ivars real, intent(in) :: var1,var2 real, optional :: lnrho,ss real, optional :: yH,lnTT real, optional :: ee,pp,cs2 ! call fatal_error('eoscalc_point', & 'This routine is not coded for eos_chemistry') ! call keep_compiler_quiet(present(lnrho),present(lnTT)) call keep_compiler_quiet(present(pp),present(ee)) call keep_compiler_quiet(present(yH)) call keep_compiler_quiet(present(ss),present(cs2)) call keep_compiler_quiet(ivars) call keep_compiler_quiet(var1,var2) ! endsubroutine eoscalc_point !*********************************************************************** subroutine eoscalc_pencil(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 cp1 to ! subroutine pressure_gradient ! 27-mar-06/tony: Introduces cv, cv1, gamma1 to make faster ! + more explicit ! 31-mar-06/tony: I removed messy lcalc_cp stuff completely. cp=1. ! is just fine. ! 22-jun-06/axel: reinstated cp,cp1,cv,cv1 in hopefully all the places. ! integer, intent(in) :: ivars real, dimension(nx), intent(in) :: var1,var2 real, dimension(nx), optional :: lnrho,ss real, dimension(nx), optional :: yH,lnTT real, dimension(nx), optional :: ee,pp,cs2 ! call fatal_error('eoscalc_pencil', & 'This routine is not coded for eos_chemistry') ! call keep_compiler_quiet(present(lnrho),present(lnTT)) call keep_compiler_quiet(present(pp),present(ee)) call keep_compiler_quiet(present(yH)) call keep_compiler_quiet(present(ss),present(cs2)) call keep_compiler_quiet(ivars) call keep_compiler_quiet(var1,var2) ! endsubroutine eoscalc_pencil !*********************************************************************** subroutine get_soundspeed(TT,cs2) ! ! Calculate sound speed for given temperature ! ! 20-Oct-03/tobi: Coded ! real, intent(in) :: TT real, intent(out) :: cs2 ! cs2=impossible call fatal_error('get_soundspeed', & 'This routine is not coded for eos_chemistry') ! call keep_compiler_quiet(TT) ! endsubroutine get_soundspeed !*********************************************************************** 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 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 ! cs2top=cs2bot call fatal_error('isothermal_entropy', & 'This routine is not coded for eos_chemistry') ! call keep_compiler_quiet(f) call keep_compiler_quiet(T0) ! 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 fatal_error('isothermal_lnrho_ss', & 'This routine is not coded for eos_chemistry') ! 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+dhrube: coded ! 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) ! ! 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 ! ! character (len=3) :: topbot real, dimension (:,:,:,:) :: f ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! 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 (mx,my,mz,mfarray) :: 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 (mx,my,mz,mfarray) :: 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_ss_temp_y !*********************************************************************** subroutine bc_ss_temp_z(f,topbot) ! ! boundary condition for entropy: constant temperature ! ! 3-aug-2002/wolf: coded ! 26-aug-2003/tony: distributed across ionization modules ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_ss_temp_z !*********************************************************************** subroutine bc_lnrho_temp_z(f,topbot) ! ! boundary condition for lnrho *and* ss: constant temperature ! ! 27-sep-2002/axel: coded ! 19-aug-2005/tobi: distributed across ionization modules ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 lnrho: constant pressure ! ! 4-apr-2003/axel: coded ! 1-may-2003/axel: added the same for top boundary ! 19-aug-2005/tobi: distributed across ionization modules ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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_chemistry.f90') ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! ! 3-aug-2002/wolf: coded ! 26-aug-2003/tony: distributed across ionization modules ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! 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 ! ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_ss_energy !*********************************************************************** subroutine bc_stellar_surface(f,topbot) ! real, dimension (:,:,:,:) :: f character (len=3) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_stellar_surface !*********************************************************************** subroutine bc_lnrho_cfb_r_iso(f,topbot) ! ! Boundary condition for radial centrifugal balance ! ! 21-aug-2006/wlad: coded ! real, dimension (:,:,:,:), intent (inout) :: f character (len=3), intent (in) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_lnrho_cfb_r_iso !*********************************************************************** subroutine bc_lnrho_hds_z_iso(f,topbot) ! ! Boundary condition for density *and* entropy. ! ! 12-Juil-2006/dintrans: coded ! real, dimension (:,:,:,:), intent (inout) :: f character (len=3), intent (in) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_lnrho_hds_z_iso !*********************************************************************** subroutine bc_lnrho_hdss_z_iso(f,topbot) ! ! Smooth out density perturbations with respect to hydrostatic ! stratification in Fourier space. ! 05-jul-07/tobi: Adapted from bc_aa_pot3 ! real, dimension (:,:,:,:), intent (inout) :: f character (len=3), intent (in) :: topbot ! call keep_compiler_quiet(f) call keep_compiler_quiet(topbot) ! endsubroutine bc_lnrho_hdss_z_iso !*********************************************************************** subroutine find_mass(element_name,MolMass) ! ! Find mass of element ! ! 05-feb-08/nils: coded ! use Mpicomm, only: stop_it ! character (len=*), intent(in) :: element_name real, intent(out) :: MolMass ! select case (element_name) case ('H') MolMass=1.00794 case ('C') MolMass=12.0107 case ('N') MolMass=14.00674 case ('O') MolMass=15.9994 case ('Ar','AR') MolMass=39.948 case ('He','HE') MolMass=4.0026 case ('S') MolMass=32.0655 case ('CLOUD') MolMass=0. case default if (lroot) print*,'element_name=',element_name call stop_it('find_mass: Element not found!') end select ! endsubroutine find_mass !*********************************************************************** subroutine find_species_index(species_name,ind_glob,ind_chem,found_specie) ! ! Find index in the f array for specie ! ! 05-feb-08/nils: coded ! integer, intent(out) :: ind_glob integer, intent(inout) :: ind_chem character (len=*), intent(in) :: species_name integer :: k logical, intent(out) :: found_specie ! ind_glob=0 ! ind_chem=0 do k=1,nchemspec if (trim(varname(ichemspec(k)))==species_name) then ind_glob=k+ichemspec(1)-1 ind_chem=k exit endif !print*, trim(varname(ichemspec(k))),(species_name) enddo ! ! Check if the species was really found ! if ((ind_glob==0)) then found_specie=.false. ! if (lroot) print*,' no species has been found ',' species index= ', ind_glob,ind_chem,species_name ! call fatal_error('find_species_index',& ! 'no species has been found') else found_specie=.true. ! if (lroot) print*,species_name,' species index= ',ind_chem endif ! endsubroutine find_species_index !*********************************************************************** subroutine read_species(input_file) ! ! This subroutine reads all species information from chem.inp ! See the chemkin manual for more information on ! the syntax of chem.inp. ! ! 06-mar-08/nils: coded ! use Mpicomm, only: stop_it ! logical :: IsSpecie=.false., emptyfile integer :: k,file_id=123, StartInd, StopInd character (len=80) :: ChemInpLine character (len=*) :: input_file ! emptyFile=.true. k=1 open(file_id,file=input_file) dataloop: do read(file_id,'(80A)',end=1000) ChemInpLine(1:80) emptyFile=.false. ! ! Check if we are reading a line within the species section ! if (ChemInpLine(1:7)=="SPECIES") IsSpecie=.true. if (ChemInpLine(1:3)=="END" .and. IsSpecie) IsSpecie=.false. ! ! Read in species ! if (IsSpecie) then if (ChemInpLine(1:7) /= "SPECIES") then StartInd=1; StopInd =0 stringloop: do StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 if (StopInd==StartInd) then StartInd=StartInd+1 else if (k>nchemspec) then print*,'nchemspec=',nchemspec call stop_it("There were too many species, "//& "please increase nchemspec!") endif varname(ichemspec(k))=trim(ChemInpLine(StartInd:StopInd-1)) StartInd=StopInd k=k+1 endif if (StartInd==80) exit enddo stringloop endif endif enddo dataloop ! ! Stop if chem.inp is empty ! 1000 if (emptyFile) call stop_it('The input file chem.inp was empty!') ! ! Check if nchemspec where not too large ! if (k0 .and. ind_chem<=nchemspec) then ! if (existing_specie) then if (found_specie) then ! ! Find molar mass ! MolMass=0 do iElement=1,4 In1=25+(iElement-1)*5 In2=26+(iElement-1)*5 In3=27+(iElement-1)*5 In4=29+(iElement-1)*5 if (ChemInpLine(In1:In1)==' ') then MolMass(iElement)=0 else element_string=trim(ChemInpLine(In1:In2)) call find_mass(element_string,MolMass(iElement)) In5=verify(ChemInpLine(In3:In4),' ')+In3-1 NumberOfElement_string=trim(ChemInpLine(In5:In4)) read (unit=NumberOfElement_string,fmt='(I5)') & NumberOfElement_i MolMass(iElement)=MolMass(iElement)*NumberOfElement_i endif enddo species_constants(ind_chem,imass)=sum(MolMass) ! ! Find temperature-ranges for low and high temperature fitting ! do iTemperature=1,3 In1=46+(iTemperature-1)*10 In2=55+(iTemperature-1)*10 if (iTemperature==3) In2=73 In3=verify(ChemInpLine(In1:In2),' ')+In1-1 TemperatureNr_i=trim(ChemInpLine(In3:In2)) read (unit=TemperatureNr_i,fmt='(F10.1)') nne tmp_temp(iTemperature)=nne enddo species_constants(ind_chem,iTemp1)=tmp_temp(1) species_constants(ind_chem,iTemp2)=tmp_temp(3) species_constants(ind_chem,iTemp3)=tmp_temp(2) ! elseif (ChemInpLine(80:80)=="2") then ! Read iaa1(1):iaa1(5) read (unit=ChemInpLine(1:75),fmt='(5E15.8)') & species_constants(ind_chem,iaa1(1):iaa1(5)) ! elseif (ChemInpLine(80:80)=="3") then ! Read iaa1(6):iaa5(3) read (unit=ChemInpLine(1:75),fmt='(5E15.8)') & species_constants(ind_chem,iaa1(6):iaa2(3)) elseif (ChemInpLine(80:80)=="4") then ! Read iaa2(4):iaa2(7) read (unit=ChemInpLine(1:75),fmt='(4E15.8)') & species_constants(ind_chem,iaa2(4):iaa2(7)) endif ! endif endif endif !(from ind_chem>0 query) endif enddo dataloop2 1001 continue close(file_id) ! endsubroutine read_thermodyn !*********************************************************************** subroutine write_thermodyn ! ! This subroutine writes the thermodynamical data for every specie ! to ./data/chem.out. ! ! 06-mar-08/nils: coded ! use General ! character (len=fnlen) :: input_file="./data/chem.out" character (len=intlen) :: ispec integer :: file_id=123,k integer, dimension(7) :: iaa1,iaa2 integer :: iTemp1=2,iTemp3=4 ! ! Initialize some index pointers ! iaa1(1)=5;iaa1(2)=6;iaa1(3)=7;iaa1(4)=8 iaa1(5)=9;iaa1(6)=10;iaa1(7)=11 ! iaa2(1)=12;iaa2(2)=13;iaa2(3)=14;iaa2(4)=15 iaa2(5)=16;iaa2(6)=17;iaa2(7)=18 ! open(file_id,file=input_file) write(file_id,*) 'Specie' write(file_id,*) 'MolMass Temp1 Temp2 Temp3' write(file_id,*) 'a1(1) a1(2) a1(3) a1(4) a1(5) a1(6) a1(7)' write(file_id,*) 'a2(1) a2(2) a2(3) a2(4) a2(5) a2(6) a2(7)' write(file_id,*) '***********************************************' dataloop2: do k=1,nchemspec write(file_id,*) varname(ichemspec(k)) write(file_id,'(F10.2,3F10.2)') species_constants(k,imass),& species_constants(k,iTemp1:iTemp3) write(file_id,'(7E12.5)') species_constants(k,iaa1) write(file_id,'(7E12.5)') species_constants(k,iaa2) enddo dataloop2 ! close(file_id) ! if (lroot) then print*,'Write pc_constants.pro in chemistry.f90' open (143,FILE=trim(datadir)//'/pc_constants.pro',POSITION="append") write (143,*) 'specname=strarr(',nchemspec,')' write (143,*) 'specmass=fltarr(',nchemspec,')' do k=1,nchemspec ispec=itoa(k-1) write (143,*) 'specname[',trim(ispec),']=',"'",& trim(varname(ichemspec(k))),"'" write (143,*) 'specmass[',trim(ispec),']=',species_constants(k,imass) enddo close (143) endif ! endsubroutine write_thermodyn !*********************************************************************** subroutine read_transport_data ! ! Reading of the chemkin transport data ! ! 01-apr-08/natalia: coded ! use Mpicomm, only: stop_it ! logical :: emptyfile logical :: found_specie integer :: file_id=123, ind_glob, ind_chem character (len=80) :: ChemInpLine character (len=10) :: specie_string integer :: VarNumber integer :: StartInd,StopInd,StartInd_1,StopInd_1 logical :: tranin=.false. logical :: trandat=.false. ! emptyFile=.true. ! StartInd_1=1; StopInd_1 =0 inquire (file='tran.dat',exist=trandat) inquire (file='tran.in',exist=tranin) if (tranin .and. trandat) then call fatal_error('eos_chemistry',& 'tran.in and tran.dat found. Please decide which one to use.') endif if (tranin) open(file_id,file='tran.in') if (trandat) open(file_id,file='tran.dat') ! if (lroot) print*, 'the following species are found '//& 'in tran.in/dat: beginning of the list:' ! dataloop: do ! read(file_id,'(80A)',end=1000) ChemInpLine(1:80) emptyFile=.false. ! StopInd_1=index(ChemInpLine,' ') specie_string=trim(ChemInpLine(1:StopInd_1-1)) ! call find_species_index(specie_string,ind_glob,ind_chem,found_specie) ! if (found_specie) then if (lroot) print*,specie_string,' ind_glob=',ind_glob,' ind_chem=',ind_chem ! VarNumber=1; StartInd=1; StopInd =0 stringloop: do while (VarNumber<7) ! StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 ! if (StopInd==StartInd) then StartInd=StartInd+1 else if (VarNumber==1) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E1.0)') & tran_data(ind_chem,VarNumber) elseif (VarNumber==2) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & tran_data(ind_chem,VarNumber) elseif (VarNumber==3) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & tran_data(ind_chem,VarNumber) elseif (VarNumber==4) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & tran_data(ind_chem,VarNumber) elseif (VarNumber==5) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & tran_data(ind_chem,VarNumber) elseif (VarNumber==6) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & tran_data(ind_chem,VarNumber) else call stop_it("No such VarNumber!") endif ! VarNumber=VarNumber+1 StartInd=StopInd endif if (StartInd==80) exit enddo stringloop ! endif enddo dataloop ! ! Stop if tran.dat is empty ! ! 1000 if (emptyFile) call stop_it('The input file tran.dat was empty!') ! if (lroot) print*, 'the following species are found in tran.dat: end of the list:' ! close(file_id) ! endsubroutine read_transport_data !*********************************************************************** subroutine read_Lewis ! ! Reading of the species Lewis numbers in an input file ! ! 21-jun-10/julien: coded ! use Mpicomm, only: stop_it ! logical :: emptyfile logical :: found_specie integer :: file_id=123, ind_glob, ind_chem, i real :: lewisk character (len=10) :: specie_string ! emptyFile=.true. ! open(file_id,file="lewis.dat") ! if (lroot) print*, 'lewis.dat: beginning of the list:' ! i=0 dataloop: do read(file_id,*,end=1000) specie_string, lewisk emptyFile=.false. ! call find_species_index(specie_string,ind_glob,ind_chem,found_specie) ! if (found_specie) then if (lroot) print*,specie_string,' ind_glob=',ind_glob,' Lewis=', lewisk Lewis_coef(ind_chem) = lewisk Lewis_coef1(ind_chem) = 1./lewisk i=i+1 endif enddo dataloop ! ! Stop if lewis.dat is empty ! ! 1000 if (emptyFile) call stop_it('End of the file!') ! if (lroot) then print*, 'lewis.dat: end of the list:' if (i == 0) & print*, 'File lewis.dat empty ===> Lewis numbers set to unity' endif ! close(file_id) ! endsubroutine read_Lewis !*********************************************************************** 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 !*********************************************************************** endmodule EquationOfState