! $Id$ ! ! MODULE_DOC: This modules adds chemical species and reactions. ! MODULE_DOC: The units used in the chem.in files are cm3,mole,sec,kcal and K ! This was found out by comparing the mechanism found in ! samples/0D/chemistry_H2_ignition_delay ! with Flow Reactor Studies and Kinetic Modeling of the ReactionH/O22 ! of A. MUELLER, T. J. KIM, R. A. YETTER, F. L. DRYER ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this ! CPARAM logical, parameter :: lchemistry = .true. ! ! MVAR CONTRIBUTION 1 ! MAUX CONTRIBUTION 1 ! ! PENCILS PROVIDED cv; cv1; cp; cp1; glncp(3); gXXk(3,nchemspec) ! PENCILS PROVIDED nu; gradnu(3); gYYk(3,nchemspec) ! PENCILS PROVIDED DYDt_reac(nchemspec); DYDt_diff(nchemspec) ! PENCILS PROVIDED lambda; glambda(3); lambda1 ! PENCILS PROVIDED Diff_penc_add(nchemspec); H0_RT(nchemspec); hhk_full(nchemspec) ! PENCILS PROVIDED ghhk(3,nchemspec); S0_R(nchemspec); cs2 ! PENCILS PROVIDED glnpp(3); del2pp; mukmu1(nchemspec) ! PENCILS PROVIDED ccondens; ppwater ! PENCILS PROVIDED Ywater ! !*************************************************************** module Chemistry ! use Cparam use Cdata use General, only: keep_compiler_quiet use EquationOfState use Messages, only: svn_id, timing, fatal_error, inevitably_fatal_error use Mpicomm, only: stop_it ! implicit none ! include 'chemistry.h' ! real :: Rgas, Rgas_unit_sys=1. real, dimension(mx,my,mz) :: cp_full, cv_full real, dimension(:,:,:), pointer :: mu1_full real, dimension(mx,my,mz) :: lambda_full, rho_full, TT_full real, dimension(mx,my,mz,nchemspec) :: cv_R_spec_full !real, dimension (mx,my,mz) :: e_int_full,cp_R_spec real, dimension(mx,my,mz,nchemspec) :: cp_spec_glo ! parameters for simplified cases real :: lambda_const=impossible real :: visc_const=impossible real :: Diff_coef_const=impossible real :: Sc_number=0.7!, Pr_number=0.7 ! real :: Cp_const=impossible real :: Cv_const=impossible logical :: lfix_Sc=.false., lfix_Pr=.false. logical :: init_from_file, reinitialize_chemistry=.false. character(len=30) :: reac_rate_method = 'chemkin' ! parameters for initial conditions real :: init_x1=-0.2, init_x2=0.2 real :: init_y1=-0.2, init_y2=0.2 real :: init_z1=-0.2, init_z2=0.2 real :: init_TT1=298., init_TT2=2400., init_ux=0., init_uy=0., init_uz=0. real :: init_zz1=0.01, init_zz2=0.2 real :: flame_pos=0. real :: init_rho2=1. real :: init_rho=1. real :: str_thick=0.02 real :: init_pressure=10.13e5 real :: global_phi=impossible ! logical :: lone_spec=.false., lfilter_strict=.false. ! ! parameters related to chemical reactions, diffusion and advection ! logical :: lreactions=.true. logical :: ladvection=.true. logical :: ldiffusion=.true. ! logical :: lnospec_eqns=.true. ! logical :: lheatc_chemistry=.true. logical :: lspecies_cond_simplified=.false. logical :: lDiff_simple=.false. logical :: lDiff_lewis=.false. logical :: lThCond_simple=.false. logical :: lT_const=.false. logical :: lDiff_fick=.false. logical :: lFlux_simple=.false. logical :: ldiff_corr=.false. logical, save :: tran_exist=.false. logical, save :: lew_exist=.false. logical :: lSmag_heat_transport=.false. logical :: lSmag_diffusion=.false. ! logical :: lfilter=.false. logical :: lkreactions_profile=.false., lkreactions_alpha=.false. integer :: nreactions=0, nreactions1=0, nreactions2=0 integer :: ll1, ll2, mm1, mm2, nn1, nn2 real, allocatable, dimension(:) :: kreactions_profile_width, kreactions_alpha ! integer :: mreactions ! ! The stociometric factors need to be reals for arbitrary reaction orders ! real, allocatable, dimension(:,:) :: stoichio, Sijm, Sijp ! integer, allocatable, dimension(:,:) :: stoichio,Sijm,Sijp real, allocatable, dimension(:,:) :: kreactions_z real, allocatable, dimension(:) :: kreactions_m, kreactions_p logical, allocatable, dimension(:) :: back character(len=30), allocatable, dimension(:) :: reaction_name logical :: lT_tanh=.false. logical :: ldamp_zone_for_NSCBC=.false. logical :: linit_temperature=.false., linit_density=.false. logical :: lreac_as_aux=.false. ! ! 1step_test case ! possible, will be removed later ! logical :: l1step_test=.false. logical :: lflame_front=.false., ltriple_flame=.false. logical :: lFlameMaster=.false. integer :: ipr=2 real :: Tc=440., Tinf=2000., beta=1.09 real :: z_cloud=0. ! ! hydro-related parameters ! real, dimension(nchemspec) :: amplchemk=0., amplchemk2=0. real, dimension(nchemspec) :: chem_diff_prefactor=1., initial_massfractions real :: amplchem=1., kx_chem=1., ky_chem=1., kz_chem=1., widthchem=1. real :: chem_diff=0. character(len=labellen), dimension(ninit) :: initchem='nothing' character(len=labellen), allocatable, dimension(:) :: kreactions_profile character(len=60) :: prerun_directory='nothing' character(len=60) :: file_name='nothing' ! real, allocatable, dimension(:,:,:,:,:) :: Bin_Diff_coef real, allocatable, dimension(:,:,:,:) :: Diff_full, Diff_full_add real, dimension(mx,my,mz,nchemspec) :: XX_full real, dimension(mx,my,mz,nchemspec) :: species_viscosity real, dimension(mx,my,mz,nchemspec), save :: RHS_Y_full real, dimension(nchemspec) :: nu_spec=0., mobility=1. ! ! Chemkin related parameters ! logical :: lcheminp=.false., lchem_cdtc=.false. logical :: lmobility=.false. ! real, dimension(nchemspec,18) :: species_constants integer :: iTemp1=2, iTemp2=3, iTemp3=4 integer, dimension(7) :: iaa1, iaa2 real, allocatable, dimension(:) :: B_n, alpha_n, E_an real, allocatable, dimension(:,:) :: low_coeff, high_coeff, troe_coeff, a_k4 logical, allocatable, dimension(:) :: Mplus_case logical, allocatable, dimension(:) :: photochem_case real :: lamb_low, lamb_up, Pr_turb=0.7 ! ! Atmospheric physics ! logical :: latmchem=.false. logical :: lcloud=.false. integer, SAVE :: index_O2=0., index_N2=0., index_O2N2=0., index_H2O=0. ! ! Lewis coefficients ! real, dimension(nchemspec) :: Lewis_coef, Lewis_coef1 ! ! Transport data ! real, dimension(nchemspec,7) :: tran_data ! ! Diagnostics ! real, allocatable, dimension(:,:) :: net_react_m, net_react_p real, dimension(nchemspec) :: Ythresh=0. logical :: lchemistry_diag=.false. ! ! Hot spot problem ! logical :: lhotspot=.false. ! ! input parameters namelist /chemistry_init_pars/ & initchem, amplchem, kx_chem, ky_chem, kz_chem, widthchem, & amplchemk,amplchemk2, chem_diff,nu_spec,lDiff_simple, & lDiff_lewis,lFlux_simple, & lThCond_simple,lambda_const, visc_const,Cp_const,Cv_const,Diff_coef_const, & init_x1,init_x2,init_y1,init_y2,init_z1,init_z2,init_TT1,& init_TT2,init_rho, & init_ux,init_uy,init_uz,l1step_test,Sc_number,init_pressure,lfix_Sc, & str_thick,lfix_Pr, lT_tanh,lT_const, lheatc_chemistry, lspecies_cond_simplified, & ldamp_zone_for_NSCBC, latmchem, lcloud, prerun_directory, & lchemistry_diag,lfilter_strict,linit_temperature, & linit_density, init_rho2, & file_name, lreac_as_aux, init_zz1, init_zz2, flame_pos, & reac_rate_method,global_phi, lSmag_heat_transport, Pr_turb, lSmag_diffusion, z_cloud, & lhotspot ! ! ! run parameters namelist /chemistry_run_pars/ & lkreactions_profile, lkreactions_alpha, & chem_diff,chem_diff_prefactor, nu_spec, ldiffusion, ladvection, & lreactions, lchem_cdtc, lheatc_chemistry, lspecies_cond_simplified, lchemistry_diag, & lmobility,mobility, lfilter,lT_tanh,lDiff_simple,lDiff_lewis,lFlux_simple, & lThCond_simple,visc_const,cp_const,reinitialize_chemistry,init_from_file, & lfilter_strict,init_TT1,init_TT2,init_x1,init_x2, linit_temperature, & linit_density, & ldiff_corr, lDiff_fick, lreac_as_aux, reac_rate_method,global_phi, & Ythresh ! ! diagnostic variables (need to be consistent with reset list below) ! integer, dimension(nchemspec) :: idiag_Ym=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_TYm=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_dYm=0 ! DIAG_DOC: $\delta\left/\delta t$ integer, dimension(nchemspec) :: idiag_dYmax=0 ! DIAG_DOC: $max\delta\left/\delta t$ integer, dimension(nchemspec) :: idiag_Ymax=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_Ymin=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_hm=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_cpm=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_diffm=0 ! DIAG_DOC: $\left$ integer, dimension(nchemspec) :: idiag_Ymz=0 ! DIAG_DOC: $\left_{xy}(z)$ integer :: idiag_dtchem=0 ! DIAG_DOC: $dt_{chem}$ ! ! integer :: idiag_cpfull=0 integer :: idiag_cvfull=0 integer :: idiag_e_intm=0 ! integer :: idiag_lambdam=0 integer :: idiag_num=0 ! ! Auxiliaries. ! integer :: ireac=0 integer, dimension(nchemspec) :: ireaci=0 contains ! !*********************************************************************** subroutine register_chemistry ! ! Configure pre-initialised (i.e. before parameter read) variables ! which should be know to be able to evaluate ! ! 13-aug-07/steveb: coded ! 8-jan-08/axel: added modifications analogously to dustdensity ! 5-mar-08/nils: Read thermodynamical data from chem.inp ! use FArrayManager ! integer :: k, ichemspec_tmp character(len=fnlen) :: input_file logical :: chemin, cheminp ! ! 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 ! ! Set ichemistry to consecutive numbers nvar+1, nvar+2, ..., nvar+nchemspec. ! call farray_register_pde('chemspec',ichemspec_tmp,array=nchemspec) do k = 1,nchemspec ichemspec(k) = ichemspec_tmp+k-1 enddo ! ! Register viscosity ! call farray_register_auxiliary('viscosity',iviscosity,communicated=.false.) ! ! Writing files for use with IDL ! if (naux+naux_com < maux+maux_com) aux_var(aux_count) = ',viscosity $' if (naux+naux_com == maux+maux_com) aux_var(aux_count) = ',viscosity' aux_count = aux_count+1 if (lroot) write (4,*) ',visocsity $' if (lroot) write (15,*) 'viscosity = fltarr(mx,my,mz)*one' ! ! Read species to be used from chem.inp (if the file exists). ! inquire (FILE='chem.inp', EXIST=lcheminp) if (.not. lcheminp) inquire (FILE='chem.in', EXIST=lcheminp) inquire (FILE='chem.inp', EXIST=cheminp) inquire (FILE='chem.in', EXIST=chemin) if (cheminp .and. chemin) call fatal_error('chemistry', & 'chem.inp and chem.in found, please decide for one') if (cheminp) input_file='chem.inp' if (chemin) input_file='chem.in' ! if (lcheminp) call read_species(input_file) ! ! Read data on the thermodynamical properties of the different species. ! All these data are stored in the array species_constants. ! if (lcheminp) call read_thermodyn(input_file) ! ! Write all data on species and their thermodynamics to file. ! if (lcheminp) call write_thermodyn ! ! Identify version number (generated automatically by SVN). ! if (lroot) call svn_id( "$Id$") ! endsubroutine register_chemistry !*********************************************************************** subroutine initialize_chemistry(f) ! ! called by run.f90 after reading parameters, but before the time loop ! ! 13-aug-07/steveb: coded ! 19-feb-08/axel: reads in chemistry.dat file ! 21-nov-10/julien: added the reaction rates as optional auxiliary variables ! in the f array for output. ! use FArrayManager use SharedVariables, only: get_shared_variable use Messages, only: warning ! real, dimension(mx,my,mz,mfarray) :: f logical :: data_file_exit=.false. logical :: exist, exist1, exist2 character(len=15) :: file1='chemistry_m.dat', file2='chemistry_p.dat' integer :: ind_glob, ind_chem, stat integer :: i logical :: found_specie=.false. character(len=2) :: car2 ! ! initialize chemistry ! if (lcheminp) then if (unit_temperature /= 1) then call fatal_error('initialize_chemistry', & 'unit_temperature must be unity when using chemistry!') endif ! ! calculate universal gas constant based on Boltzmann constant ! and the proton mass ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif endif ! if (nchemspec == 1) then lone_spec = .true. lreactions = .false. ladvection = .false. ldiffusion = .false. endif ! ! check for the existence of chemistry input files ! inquire (file=file1,exist=exist1) inquire (file=file2,exist=exist2) inquire (file='chemistry.dat',exist=exist) ! ! Read in data file in ChemKin format ! if (lcheminp) then call chemkin_data(f) data_file_exit = .true. if (latmchem) then call find_species_index('O2',ind_glob,ind_chem,found_specie) if (found_specie) then index_O2 = ind_chem else call fatal_error('initialize_chemistry', 'no O2 has been found') endif call find_species_index('N2',ind_glob,ind_chem,found_specie) if (found_specie) then index_N2 = ind_chem else call fatal_error('initialize_chemistry', 'no N2 has been found') endif call find_species_index('O2N2',ind_glob,ind_chem,found_specie) if (found_specie) then index_O2N2 = ind_chem else call fatal_error('initialize_chemistry', 'no O2N2 has been found') endif call find_species_index('H2O',ind_glob,ind_chem,found_specie) if (found_specie) then index_H2O = ind_chem else call fatal_error('initialize_chemistry', 'no H2O has been found') endif endif if (lcloud) then call find_species_index('H2O',ind_glob,ind_chem,found_specie) if (found_specie) then index_H2O = ind_chem else call fatal_error('initialize_chemistry', 'no H2O has been found') endif endif endif ! ! Alternatively, read in stoichiometric matrices in explicit format. ! For historical reasons this is referred to as "astrobiology_data" ! if (exist1 .and. exist2) then call astrobiology_data(f) data_file_exit = .true. endif ! if (exist) then call astrobiology_data(f) data_file_exit = .true. endif ! ! check the existence of a data file ! if (.not. data_file_exit) then call stop_it('initialize_chemistry: there is no chemistry data file') 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 ! ! Reinitialize if required ! if (reinitialize_chemistry) then if (lroot) print*,'Reinitializing chemistry.' call init_chemistry(f) endif ! ! allocate memory for net_reaction diagnostics ! if (allocated(net_react_p) .and. .not. lreloading) then print*, 'this should not be here' endif if (lchemistry_diag .and. .not. lreloading) then allocate(net_react_p(nchemspec,nreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for net_react_p") net_react_p = 0. allocate(net_react_m(nchemspec,nreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for net_react_m") net_react_m = 0. endif ! ! Define the chemical reaction rates as auxiliary variables for output ! if (lreac_as_aux) then if (ireac == 0) then call farray_register_auxiliary('reac',ireac,array=nchemspec) do i = 0, nchemspec-1 ireaci(i+1) = ireac+i enddo else if (lroot) print*, 'initialize_chemistry: ireac = ', ireac call farray_index_append('ireac',ireac) do i = 1, nchemspec write (car2,'(i2)') i call farray_index_append('ireac'//trim(adjustl(car2)),ireaci(i)) enddo endif endif ! if (leos) then call get_shared_variable('mu1_full',mu1_full,caller='initialize_chemistry') else call warning('initialize_chemistry','mu1_full not provided by eos') endif ! ! write array dimension to chemistry diagnostics file ! open (1,file=trim(datadir)//'/net_reactions.dat',position='append') write (1,*) nchemspec,nreactions close (1) ! endsubroutine initialize_chemistry !*********************************************************************** subroutine init_chemistry(f) ! ! initialise chemistry initial condition; called from start.f90 ! ! 13-aug-07/steveb: coded ! jul-10/julien: Added some new initial cases ! use Initcond use InitialCondition, only: initial_condition_chemistry ! real, dimension(mx,my,mz,mfarray) :: f real :: PP integer :: j,k logical :: lnothing, air_exist ! intent(inout) :: f ! ! different initializations of nd (called from start) ! lnothing = .false. do j = 1,ninit select case (initchem(j)) ! case ('nothing') if (lroot .and. .not. lnothing) print*, 'initchem: nothing ' lnothing = .true. case ('constant') do k = 1,nchemspec f(:,:,:,ichemspec(k)) = f(:,:,:,ichemspec(k)) + amplchemk2(k) enddo case ('positive-noise') do k = 1,nchemspec call posnoise(amplchemk(k),f,ichemspec(k)) enddo case ('positive-noise-rel') do k = 1,nchemspec call posnoise_rel(amplchemk(k),amplchemk2(k),f,ichemspec(k)) enddo case ('innerbox') do k = 1,nchemspec call innerbox(amplchemk(k),amplchemk2(k),f,ichemspec(k),widthchem) enddo case ('cos2x_cos2y_cos2z') do k = 1,nchemspec call cos2x_cos2y_cos2z(amplchemk(k),f,ichemspec(k)) enddo case ('coswave-x') do k = 1,nchemspec call coswave(amplchem,f,ichemspec(k),kx=kx_chem) enddo case ('gaussian-x') do k = 1,nchemspec call gaussian(amplchem,f,ichemspec(k),kx=kx_chem) enddo case ('gaussian-pos') do k = 1,nchemspec call gaussianpos(amplchemk(k),f,ichemspec(k),widthchem,kx_chem,ky_chem,kz_chem) print*,"c(",x(l1),",",y(m1),",",z(n1),")=", f(l1,m1,n1,ichemspec(k)) enddo case ('hatwave-x') do k = 1,nchemspec call hatwave(amplchem,f,ichemspec(k),widthchem,kx=kx_chem,pos=1.) enddo case ('hatwave-y') do k = 1,nchemspec call hatwave(amplchem,f,ichemspec(k),widthchem,ky=ky_chem) enddo case ('hatwave-z') do k = 1,nchemspec call hatwave(amplchem,f,ichemspec(k),widthchem,kz=kz_chem) enddo case ('air') if (lroot ) print*, 'init_chem: air ' inquire (file='air.dat',exist=air_exist) if (.not. air_exist) then inquire (file='air.in',exist=air_exist) endif if (air_exist) then call air_field(f,PP) else call stop_it('there is no air.in/dat file') endif case ('flame_front') call flame_front(f) case ('TTD') call TTD(f) case ('triple_flame') call triple_flame(f) case ('flame') if (lroot) print*, 'initchem: flame ' call flame(f) case ('flame_blob') call flame_blob(f) case ('opposite_flames') call opposite_flames(f) case ('opposite_ignitions') call opposite_ignitions(f) case ('prerun_1D') call prerun_1D(f,prerun_directory) case ('prerun_1D_opp') call prerun_1D_opp(f,prerun_directory) case ('FlameMaster') call FlameMaster_ini(f,file_name) case ('flame_front_new') call flame_front_new(f) case default ! ! Catch unknown values ! if (lroot) print*, 'initchem: No such value for initchem: ', & trim(initchem(j)) call stop_it('') ! endselect enddo ! ! Interface for user's own initial condition ! if (linitial_condition) call initial_condition_chemistry(f) ! ! The following lines are kept temporally ! ! if (lone_spec) then ! f(:,:,:,ichemspec(1))=1. ! if (lroot) print*, 'initchem: this is one specie case' ! endif ! endsubroutine init_chemistry !*********************************************************************** subroutine pencil_criteria_chemistry ! ! All pencils that this chemistry module depends on are specified here. ! ! 13-aug-07/steveb: coded ! lpenc_requested(i_gXXk) = .true. lpenc_requested(i_gYYk) = .true. ! if (lreactions) lpenc_requested(i_ghhk) = .true. ! if (lreactions) lpenc_requested(i_DYDt_reac) = .true. lpenc_requested(i_DYDt_diff) = .true. ! if (lcheminp) then lpenc_requested(i_rho) = .true. lpenc_requested(i_cv) = .true. lpenc_requested(i_cp) = .true. lpenc_requested(i_cv1) = .true. lpenc_requested(i_cp1) = .true. ! if (lreactions) lpenc_requested(i_H0_RT) = .true. ! if (lreactions) lpenc_requested(i_S0_R) = .true. lpenc_requested(i_nu) = .true. lpenc_requested(i_gradnu) = .true. lpenc_requested(i_cs2) = .true. ! ! if (lreactions) lpenc_requested(i_hhk_full) = .true. if (lThCond_simple) lpenc_requested(i_glncp) = .true. ! if (lheatc_chemistry) then lpenc_requested(i_lambda) = .true. lpenc_requested(i_glambda) = .true. lpenc_requested(i_lambda1) = .true. if (lSmag_heat_transport) lpenc_requested(i_sij2) = .true. endif ! if (latmchem .or. lcloud) then lpenc_requested(i_ppwater) = .true. lpenc_requested(i_Ywater) = .true. endif ! if (ldiffusion .or. lparticles_chemistry) then lpenc_requested(i_Diff_penc_add) = .true. if (.not. lDiff_fick) then lpenc_requested(i_mukmu1) = .true. lpenc_requested(i_glnmu) = .true. endif endif endif ! endsubroutine pencil_criteria_chemistry !*********************************************************************** subroutine pencil_interdep_chemistry(lpencil_in) ! ! Interdependency among pencils provided by this module are specified here ! ! 02-03-08/Natalia: coded ! logical, dimension(npencils) :: lpencil_in ! if (lpencil_in(i_cv1)) lpencil_in(i_cv) = .true. if (lpencil_in(i_cp1)) lpencil_in(i_cp) = .true. if (lpencil_in(i_glambda).or. lpencil_in(i_lambda1)) & lpencil_in(i_lambda) = .true. ! if (lpencil_in(i_H0_RT).or. lpencil_in(i_S0_R)) then lpencil_in(i_TT) = .true. lpencil_in(i_TT1) = .true. endif if (lpencil_in(i_S0_R)) then lpencil_in(i_lnTT) = .true. endif if (lpencil_in(i_DYDt_reac)) then lpencil_in(i_TT) = .true. lpencil_in(i_TT1) = .true. lpencil_in(i_lnTT) = .true. lpencil_in(i_H0_RT) = .true. lpencil_in(i_mu1) = .true. lpencil_in(i_rho) = .true. endif if (lpencil_in(i_hhk_full)) then lpencil_in(i_H0_RT) = .true. lpencil_in(i_TT) = .true. endif if (lpencil_in(i_ghhk)) then lpencil_in(i_glnTT) = .true. lpencil_in(i_TT) = .true. endif if (lpencil_in(i_ppwater)) then lpencil_in(i_TT) = .true. lpencil_in(i_rho) = .true. endif ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_chemistry !*********************************************************************** subroutine calc_pencils_chemistry(f,p) ! ! Calculate chemistry pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 13-aug-07/steveb: coded ! 10-jan-11/julien: adapted for the case where chemistry is solved by LSODE ! use Sub, only: grad ! real, dimension(mx,my,mz,mfarray) :: f type (pencil_case) :: p real, dimension(nx,3) :: gXX_tmp, glncp_tmp!, ghhk_tmp ! intent(in) :: f intent(inout) :: p integer :: k,i integer :: ii1=1, ii2=2, ii3=3, ii4=4, ii5=5, ii6=6, ii7=7 real :: T_low, T_up, T_mid real, dimension(nx) :: T_loc, TT_2, TT_3, TT_4 logical :: ldiffusion2 ! ldiffusion2 = ldiffusion .and. (.not. lchemonly) ! if (lpencil(i_gYYk)) then do k = 1,nchemspec call grad(f(:,:,:,ichemspec(k)),gXX_tmp) do i = 1,3 p%gYYk(:,i,k) = gXX_tmp(:,i) enddo enddo endif ! if (lpencil(i_gXXk)) then do k = 1,nchemspec call grad(XX_full(:,:,:,k),gXX_tmp) do i = 1,3 p%gXXk(:,i,k) = gXX_tmp(:,i) enddo enddo endif ! if (lpencil(i_mukmu1)) then do k = 1,nchemspec p%mukmu1(:,k) = species_constants(k,imass)/unit_mass*p%mu1(:) enddo endif ! if (lcheminp) then ! if (lpencil(i_glncp) .and. lThCond_simple) then call grad(cp_full,glncp_tmp) do i = 1,3 p%glncp(:,i) = glncp_tmp(:,i)/cp_full(l1:l2,m,n) enddo endif ! ! Specific heat at constant volume (i.e. density) ! if (lpencil(i_cv)) p%cv = cv_full(l1:l2,m,n) if (lpencil(i_cv1)) p%cv1 = 1./p%cv if (lpencil(i_cp)) p%cp = cp_full(l1:l2,m,n) if (lpencil(i_cp1)) p%cp1 = 1./p%cp ! TT_2 = p%TT*p%TT TT_3 = TT_2*p%TT TT_4 = TT_2*TT_2 ! ! Viscosity of a mixture ! if (lpencil(i_nu).and.(.not. lchemonly)) then if (visc_const < impossible) then p%nu = visc_const else p%nu = f(l1:l2,m,n,iviscosity) endif if (lpencil(i_gradnu)) then if (visc_const < impossible) then p%gradnu = 0. else call grad(f(:,:,:,iviscosity),p%gradnu) endif endif endif ! endif ! ! Dimensionless Standard-state molar enthalpy H0/RT ! if (lpencil(i_H0_RT)) then if ((.not. lT_const).and.(ilnTT /= 0)) then do k = 1,nchemspec T_low = species_constants(k,iTemp1) T_mid = species_constants(k,iTemp2) T_up = species_constants(k,iTemp3) ! ! Natalia:pencil_check ! if lpencil_check and full compiler settings then ! the problem appears for T_loc= p%TT ! Does anybody know why it is so? ! While this problem is not resolved ! I use T_loc= exp(f(l1:l2,m,n,ilnTT)) ! NILS: This is going to fail if nologtemperature=T. The real error ! NILS: should be found instead of making a quick fix. ! NILS: I am not able to reproduce the error natalia reported. ! ! T_loc= exp(f(l1:l2,m,n,ilnTT)) T_loc = p%TT where (T_loc <= T_mid) p%H0_RT(:,k) = species_constants(k,iaa2(ii1)) & +species_constants(k,iaa2(ii2))*T_loc/2 & +species_constants(k,iaa2(ii3))*TT_2/3 & +species_constants(k,iaa2(ii4))*TT_3/4 & +species_constants(k,iaa2(ii5))*TT_4/5 & +species_constants(k,iaa2(ii6))/T_loc elsewhere p%H0_RT(:,k) = species_constants(k,iaa1(ii1)) & +species_constants(k,iaa1(ii2))*T_loc/2 & +species_constants(k,iaa1(ii3))*TT_2/3 & +species_constants(k,iaa1(ii4))*TT_3/4 & +species_constants(k,iaa1(ii5))*TT_4/5 & +species_constants(k,iaa1(ii6))/T_loc endwhere enddo ! ! Enthalpy flux ! if (lpencil(i_hhk_full) ) then do k = 1,nchemspec if (species_constants(k,imass) > 0.) then p%hhk_full(:,k) = p%H0_RT(:,k)*Rgas*T_loc & /species_constants(k,imass) endif enddo endif ! if (lpencil(i_ghhk) .and. (.not. lchemonly)) then do k = 1,nchemspec if (species_constants(k,imass) > 0.) then do i = 1,3 p%ghhk(:,i,k) = (cv_R_spec_full(l1:l2,m,n,k)+1) & /species_constants(k,imass)*Rgas*p%glnTT(:,i)*T_loc(:) enddo endif enddo endif ! endif endif ! ! Find the entropy by using fifth order temperature fitting function ! if (lpencil(i_S0_R) .and. (.not. llsode .or. lchemonly)) then !AB: Natalia, maybe we should ask earlier for lentropy? if ((.not. lT_const).and.(ilnTT /= 0)) then do k = 1,nchemspec T_low = species_constants(k,iTemp1) T_mid = species_constants(k,iTemp2) T_up = species_constants(k,iTemp3) ! ! Natalia:pencil_check ! if lpencil_check and full compiler settings then ! the problem appears for T_loc= p%TT ! Does anybody know why it is so? ! While this problem is not resolved ! I use T_loc= exp(f(l1:l2,m,n,ilnTT)) ! NILS: This is going to fail if nologtemperature=T. The real error ! NILS: should be found instead of making a quick fix. ! NILS: I am not able to reproduce the error natalia reported. ! T_loc = p%TT !T_loc= exp(f(l1:l2,m,n,ilnTT)) where (T_loc <= T_mid .and. T_low <= T_loc) p%S0_R(:,k) = species_constants(k,iaa2(ii1))*p%lnTT & +species_constants(k,iaa2(ii2))*T_loc & +species_constants(k,iaa2(ii3))*TT_2/2 & +species_constants(k,iaa2(ii4))*TT_3/3 & +species_constants(k,iaa2(ii5))*TT_4/4 & +species_constants(k,iaa2(ii7)) elsewhere (T_mid <= T_loc .and. T_loc <= T_up) p%S0_R(:,k) = species_constants(k,iaa1(ii1))*p%lnTT & +species_constants(k,iaa1(ii2))*T_loc & +species_constants(k,iaa1(ii3))*TT_2/2 & +species_constants(k,iaa1(ii4))*TT_3/3 & +species_constants(k,iaa1(ii5))*TT_4/4 & +species_constants(k,iaa1(ii7)) endwhere enddo endif endif ! ! Calculate the reaction term and the corresponding pencil ! if (lreactions .and. lpencil(i_DYDt_reac)) then if (.not. llsode .or. lchemonly) then call calc_reaction_term(f,p) else p%DYDt_reac = 0. endif else p%DYDt_reac = 0. endif ! ! Calculate the thermal diffusivity ! if (lpencil(i_lambda) .and. lheatc_chemistry) then if ((lThCond_simple) .or. (lambda_const < impossible)) then if (lThCond_simple) then ! ! 04-18-11/Julien: Modified the computation of simple heat conductivity according ! to Smooke & Giovangigli 1991. lambda_const is now equal ! to 2.58e-4 ! instead of 1e4, and represents the ratio \lambda0/cp0. ! Formula: ! \lambda = lambda_const*cp*(T/T0)**0.7, with T0 now = 298K ! if (lambda_const == impossible) lambda_const = 2.58e-4 p%lambda = lambda_const*p%cp*exp(0.7*log(p%TT(:)/298.)) if (lpencil(i_glambda)) then do i = 1,3 p%glambda(:,i) = p%lambda(:)*(0.7*p%glnTT(:,i)+p%glncp(:,i)) enddo endif elseif ((.not. lThCond_simple) .and. (lambda_const < impossible)) then p%lambda = lambda_const if (lpencil(i_glambda)) p%glambda = 0. endif elseif (lSmag_heat_transport) then ! ! Natalia ! turbulent heat transport in Smagorinsky case ! probably it should be moved to viscosity module ! p%lambda=(0.15*dxmax)**2.*sqrt(2*p%sij2)/Pr_turb*p%cv*p%rho else p%lambda = lambda_full(l1:l2,m,n) if (lpencil(i_glambda)) call grad(lambda_full,p%glambda) endif if (lpencil(i_lambda1)) p%lambda1 = 1./max(tini,p%lambda) endif ! ! Calculate the diffusion term and the corresponding pencil ! if (lcheminp) then ! ! There are 4 cases: ! 1) the case of simplifyed expression for the difusion coef. (Oran paper,) ! 2) the case of constant diffusion coefficients ! 3) the case of constant Lewis numbers with diffusion coef. depending ! on heat diffusivity ! 4) full complex diffusion coefficients ! if (lpencil(i_Diff_penc_add)) then if (lDiff_simple) then ! ! 04-18-11/Julien: Changed the value of Diff_coef_const from 10 to 2.58e-4 ! according to Smooke & Giovangigli 1991. Now Diff_coef_const ! does not represent a constant diffusion coefficient, ! but \rho0 D0. ! Diff_penc_add are still the diffusion coefficients. Formula: ! D = Diff_coef_const/\rho*(T/T0)**n0.7, with T0 now = 298K. ! if (Diff_coef_const == impossible) Diff_coef_const = 2.58e-4 do k = 1,nchemspec p%Diff_penc_add(:,k) = & Diff_coef_const*p%rho1*exp(0.7*log(p%TT(:)/298.)) if (lew_exist) p%Diff_penc_add(:,k) = & p%Diff_penc_add(:,k)*Lewis_coef1(k) enddo ! ! Constant diffusion coefficients ! elseif ((.not. lDiff_simple) .and. (Diff_coef_const < impossible)) then if (lSmag_diffusion) then if (lcloud) then if (z(n)>=z_cloud) then do k = 1,nchemspec p%Diff_penc_add(:,k) = (0.15*dxmax)**2.*sqrt(2*p%sij2)/Sc_number enddo else do k = 1,nchemspec p%Diff_penc_add(:,k) = Diff_coef_const*p%rho1 enddo endif else do k = 1,nchemspec p%Diff_penc_add(:,k) = (0.15*dxmax)**2.*sqrt(2*p%sij2)/Sc_number enddo endif else do k = 1,nchemspec p%Diff_penc_add(:,k) = Diff_coef_const*p%rho1 enddo endif ! ! Diffusion coefficient of a mixture with constant Lewis numbers and ! given heat conductivity ! elseif (lDiff_lewis .and. lew_exist) then do k = 1,nchemspec p%Diff_penc_add(:,k) = p%lambda*p%rho1*p%cp1*Lewis_coef1(k) enddo elseif (lDiff_lewis .and. l1step_test) then p%Diff_penc_add(:,k) = p%lambda*p%rho1*p%cp1 ! ! Full diffusion coefficient case ! else do k = 1,nchemspec p%Diff_penc_add(:,k) = Diff_full_add(l1:l2,m,n,k) enddo endif endif endif ! ! More initialization of pencils ! if (ldiffusion2 .and. lpencil(i_DYDt_diff)) then if (.not. lchemonly) then call calc_diffusion_term(f,p) else p%DYDt_diff = 0. endif else p%DYDt_diff = 0. endif ! if (latmchem) then RHS_Y_full(l1:l2,m,n,:) = p%DYDt_diff else RHS_Y_full(l1:l2,m,n,:) = p%DYDt_reac+p%DYDt_diff endif ! if (lpencil(i_cs2) .and. lcheminp) then if (any(p%cv == 0.0)) then else p%cs2 = p%cp*p%cv1*p%mu1*p%TT*Rgas endif endif ! if (lpencil(i_ppwater) .and. .not. lchemonly) then if (index_H2O > 0) then p%ppwater = p%rho*Rgas*p%TT/18.*f(l1:l2,m,n,ichemspec(index_H2O)) endif endif if (lpencil(i_Ywater) .and. .not. lchemonly) then if (index_H2O > 0) then p%Ywater = f(l1:l2,m,n,ichemspec(index_H2O)) endif endif ! ! Energy per unit mass ! if (lpencil(i_ee)) p%ee = p%cv*p%TT ! endsubroutine calc_pencils_chemistry !*********************************************************************** subroutine flame_front(f) ! ! 05-jun-09/Nils Erland L. Haugen: adapted from similar ! routine in special/chem_stream.f90 ! 24-jun-10/Julien Savre: Modifications for lean methane/air combustion ! This routine set up the initial profiles used in 1D flame speed measurments ! real, dimension(mx,my,mz,mfarray) :: f integer :: i, j,k ! real :: mO2=0., mH2=0., mN2=0., mH2O=0., mCH4=0., mCO2=0. real :: log_inlet_density, del, PP integer :: i_H2=0, i_O2=0, i_H2O=0, i_N2=0 integer :: ichem_H2=0, ichem_O2=0, ichem_N2=0, ichem_H2O=0 integer :: i_CH4=0, i_CO2=0, ichem_CH4=0, ichem_CO2=0 real :: initial_mu1, final_massfrac_O2, final_massfrac_CH4, & final_massfrac_H2O, final_massfrac_CO2 real :: init_H2, init_O2, init_N2, init_H2O, init_CO2, init_CH4 logical :: lH2=.false., lO2=.false., lN2=.false., lH2O=.false. logical :: lCH4=.false., lCO2=.false. ! lflame_front = .true. ! call air_field(f,PP) ! if (ltemperature_nolog) f(:,:,:,ilnTT) = log(f(:,:,:,ilnTT)) ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,lH2) if (lH2) then mH2 = species_constants(ichem_H2,imass) init_H2 = initial_massfractions(ichem_H2) endif call find_species_index('O2',i_O2,ichem_O2,lO2) if (lO2) then mO2 = species_constants(ichem_O2,imass) init_O2 = initial_massfractions(ichem_O2) endif call find_species_index('N2',i_N2,ichem_N2,lN2) if (lN2) then mN2 = species_constants(ichem_N2,imass) init_N2 = initial_massfractions(ichem_N2) else init_N2 = 0 endif call find_species_index('H2O',i_H2O,ichem_H2O,lH2O) if (lH2O) then mH2O = species_constants(ichem_H2O,imass) init_H2O = initial_massfractions(ichem_H2O) endif call find_species_index('CH4',i_CH4,ichem_CH4,lCH4) if (lCH4) then mCH4 = species_constants(ichem_CH4,imass) init_CH4 = initial_massfractions(ichem_CH4) endif call find_species_index('CO2',i_CO2,ichem_CO2,lCO2) if (lCO2) then mCO2 = species_constants(ichem_CO2,imass) init_CO2 = initial_massfractions(ichem_CO2) final_massfrac_CO2 = init_CO2 else init_CO2 = 0 final_massfrac_CO2 = init_CO2 endif ! ! Find approximate value for the mass fraction of O2 after the flame front ! Warning: These formula are only correct for lean fuel/air mixtures. They ! must be modified under rich conditions to account for the excess ! of fuel. ! final_massfrac_O2 = 0. if (lH2 .and. .not. lCH4) then final_massfrac_H2O = mH2O/mH2 * init_H2 final_massfrac_O2 = 1. - final_massfrac_H2O- init_N2 elseif (lCH4) then final_massfrac_CH4 = 0. final_massfrac_H2O = 2.*mH2O/mCH4 * init_CH4 final_massfrac_CO2 = mCO2/mCH4 * init_CH4 final_massfrac_O2 = & 1. - final_massfrac_CO2 - final_massfrac_H2O & - init_N2 endif ! if (final_massfrac_O2 < 0.) final_massfrac_O2 = 0. if (lroot) then print*, ' init final' if (lH2 .and. .not. lCH4) print*, 'H2 :', init_H2, 0. if (lCH4) print*, 'CH4 :', init_CH4, 0. if (lO2) print*, 'O2 :', init_O2, final_massfrac_O2 if (lH2O) print*, 'H2O :', init_H2O, final_massfrac_H2O if (lCO2) print*, 'CO2 :', 0., final_massfrac_CO2 endif ! ! Initialize temperature and species ! do k = 1,mx ! ! Initialize temperature ! if (lT_tanh) then del = init_x2-init_x1 f(k,:,:,ilnTT) = f(k,:,:,ilnTT)+log((init_TT2+init_TT1)*0.5 & +((init_TT2-init_TT1)*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del))) else if (x(k) <= init_x1) f(k,:,:,ilnTT) = log(init_TT1) if (x(k) >= init_x2) f(k,:,:,ilnTT) = log(init_TT2) if (x(k) > init_x1 .and. x(k) < init_x2) & f(k,:,:,ilnTT) = log((x(k)-init_x1)/(init_x2-init_x1) & *(init_TT2-init_TT1)+init_TT1) endif ! ! Initialize fuel ! if (lT_tanh) then if (lH2 .and. .not. lCH4) then del = (init_x2-init_x1) f(k,:,:,i_H2) = (0.+f(l1,:,:,i_H2))*0.5+(0.-f(l1,:,:,i_H2))*0.5 & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) endif if (lCH4) then if (lroot) print*, 'No tanh initial function available for CH4 combustion.' endif ! else if (x(k) > init_x1) then if (lH2 .and. .not. lCH4) f(k,:,:,i_H2) = init_H2* & (exp(f(k,:,:,ilnTT))-init_TT2)/(init_TT1-init_TT2) if (lCH4) f(k,:,:,i_CH4) = init_CH4*(exp(f(k,:,:,ilnTT))-init_TT2) & /(init_TT1-init_TT2) endif endif ! ! Initialize oxygen ! if (lT_tanh) then del = (init_x2-init_x1) f(k,:,:,i_O2) = (f(l2,:,:,i_O2)+f(l1,:,:,i_O2))*0.5 & +((f(l2,:,:,i_O2)-f(l1,:,:,i_O2))*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) else if (x(k) > init_x2) f(k,:,:,i_O2) = final_massfrac_O2 if (x(k) > init_x1 .and. x(k) <= init_x2) & f(k,:,:,i_O2) = (x(k)-init_x1)/(init_x2-init_x1) & *(final_massfrac_O2-init_O2)+init_O2 endif enddo ! ! Initialize products ! if (lT_tanh) then do k = 1,mx if (lH2 .and. .not. lCH4) then del = (init_x2-init_x1) f(k,:,:,i_H2O) = (f(l1,:,:,i_H2)/2.*18.+f(l1,:,:,i_H2O))*0.5 & +((f(l1,:,:,i_H2)/2.*18.-f(l1,:,:,i_H2O))*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) endif if (lCH4) then if (lroot) print*, 'No tanh initial function available for CH4 combustion.' endif enddo ! elseif (.not. l1step_test) then do k = 1,mx if (x(k) >= init_x1 .and. x(k) < init_x2) then f(k,:,:,i_H2O) = (x(k)-init_x1)/(init_x2-init_x1) & *final_massfrac_H2O if (lCO2) f(k,:,:,i_CO2) = (x(k)-init_x1)/(init_x2-init_x1) & *final_massfrac_CO2 elseif (x(k) >= init_x2) then if (lCO2) f(k,:,:,i_CO2) = final_massfrac_CO2 if (lH2O) f(k,:,:,i_H2O) = final_massfrac_H2O endif enddo endif ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif ! ! Find logaritm of density at inlet ! initial_mu1 = & initial_massfractions(ichem_O2)/(mO2) & +initial_massfractions(ichem_H2O)/(mH2O) & +initial_massfractions(ichem_N2)/(mN2) if (lH2 .and. .not. lCH4) initial_mu1 = initial_mu1+ & initial_massfractions(ichem_H2)/(mH2) if (lCO2) initial_mu1 = initial_mu1+init_CO2/(mCO2) if (lCH4) initial_mu1 = initial_mu1+init_CH4/(mCH4) log_inlet_density = & log(init_pressure)-log(Rgas)-log(init_TT1)-log(initial_mu1) ! ! Initialize density ! call getmu_array(f,mu1_full) f(l1:l2,m1:m2,n1:n2,ilnrho) = log(init_pressure)-log(Rgas) & -f(l1:l2,m1:m2,n1:n2,ilnTT)-log(mu1_full(l1:l2,m1:m2,n1:n2)) ! ! Initialize velocity ! ! f(l1:l2,m1:m2,n1:n2,iux)=exp(log_inlet_density - f(l1:l2,m1:m2,n1:n2,ilnrho)) & ! * (f(l1:l2,m1:m2,n1:n2,iux)+init_ux) f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux)+init_ux ! ! Check if we want nolog of density or nolog of temperature ! if (ldensity_nolog) & f(l1:l2,m1:m2,n1:n2,irho) = exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) if (ltemperature_nolog) & f(l1:l2,m1:m2,n1:n2,iTT) = exp(f(l1:l2,m1:m2,n1:n2,ilnTT)) ! ! Renormalize all species to be sure that the sum of all mass fractions ! are unity ! do i = 1,mx do j = 1,my do k = 1,mz f(i,j,k,ichemspec) = f(i,j,k,ichemspec)/sum(f(i,j,k,ichemspec)) enddo enddo enddo ! endsubroutine flame_front !*********************************************************************** subroutine flame_front_new(f) ! ! 10-nov-18/chengeng: adapted from flame_front, but flame front on the left, ! and added initial condition for hotspot problem. ! This version replaces experimental/test_chemistry. ! ! This routine set up the initial profiles used in 1D flame speed measurments ! real, dimension(mx,my,mz,mfarray) :: f integer :: i, j,k ! real :: mO2=0., mH2=0., mN2=0., mH2O=0., mCH4=0., mCO2=0. real :: log_inlet_density, del, PP integer :: i_H2=0, i_O2=0, i_H2O=0, i_N2=0 integer :: ichem_H2=0, ichem_O2=0, ichem_N2=0, ichem_H2O=0 integer :: i_CH4=0, i_CO2=0, ichem_CH4=0, ichem_CO2=0 real :: initial_mu1, final_massfrac_O2, final_massfrac_CH4, & final_massfrac_H2O, final_massfrac_CO2,final_massfrac_H2 real :: init_H2, init_O2, init_N2, init_H2O, init_CO2, init_CH4 logical :: lH2=.false., lO2=.false., lN2=.false., lH2O=.false. logical :: lCH4=.false., lCO2=.false. real :: theta ! lflame_front = .true. ! call air_field(f,PP) ! if (ltemperature_nolog) f(:,:,:,ilnTT) = log(f(:,:,:,ilnTT)) ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,lH2) if (lH2) then mH2 = species_constants(ichem_H2,imass) init_H2 = initial_massfractions(ichem_H2) endif call find_species_index('O2',i_O2,ichem_O2,lO2) if (lO2) then mO2 = species_constants(ichem_O2,imass) init_O2 = initial_massfractions(ichem_O2) endif call find_species_index('N2',i_N2,ichem_N2,lN2) if (lN2) then mN2 = species_constants(ichem_N2,imass) init_N2 = initial_massfractions(ichem_N2) else init_N2 = 0 endif call find_species_index('H2O',i_H2O,ichem_H2O,lH2O) if (lH2O) then mH2O = species_constants(ichem_H2O,imass) init_H2O = initial_massfractions(ichem_H2O) endif ! ! Find approximate value for the mass fraction of O2 after the flame front ! Warning: These formula are only correct for lean fuel/air mixtures. They ! must be modified under rich conditions to account for the excess ! of fuel. ! final_massfrac_O2 = 0. final_massfrac_H2 = 0. if ( lH2 ) then final_massfrac_H2O = mH2O/mH2 * init_H2 final_massfrac_O2 = 1. - final_massfrac_H2O- init_N2 endif ! if (final_massfrac_O2 < 0.) final_massfrac_O2 = 0. if (lroot) then print*, ' init final' if (lH2 ) print*, 'H2 :' , init_H2, final_massfrac_H2 if (lO2) print*, 'O2 :' , init_O2, final_massfrac_O2 if (lH2O) print*, 'H2O :', init_H2O, final_massfrac_H2O endif ! ! Initialize temperature and species ! do k = 1,mx ! ! Initialize temperature ! if( lhotspot )then if( x(k)init_x2 )then f(k,:,:,ilnTT) = log( init_TT2 ) f(k,:,:,i_H2) = init_H2 f(k,:,:,i_O2) = init_O2 f(k,:,:,i_H2O) = init_H2O f(k,:,:,iux) = init_ux*(init_TT1/init_TT2-1.0) else theta = ( x(k)-init_x1 )/( init_x2-init_x1 ) f(k,:,:,ilnTT) = log( init_TT1-theta*(init_TT1-init_TT2) ) f(k,:,:,i_H2) = init_H2*theta+final_massfrac_H2 f(k,:,:,i_O2) = init_O2*theta+final_massfrac_O2 f(k,:,:,i_H2O) = init_H2O*(1.0-theta)+final_massfrac_H2O f(k,:,:,iux) = init_ux*(init_TT1/init_TT2-1.0)*theta end if end if end do ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif ! ! Find logaritm of density at inlet ! initial_mu1 = & initial_massfractions(ichem_O2)/(mO2) & +initial_massfractions(ichem_H2O)/(mH2O) & +initial_massfractions(ichem_N2)/(mN2) if (lH2 .and. .not. lCH4) initial_mu1 = initial_mu1+ & initial_massfractions(ichem_H2)/(mH2) if (lCO2) initial_mu1 = initial_mu1+init_CO2/(mCO2) if (lCH4) initial_mu1 = initial_mu1+init_CH4/(mCH4) log_inlet_density = & log(init_pressure)-log(Rgas)-log(init_TT1)-log(initial_mu1) ! ! Initialize density ! call getmu_array(f,mu1_full) f(l1:l2,m1:m2,n1:n2,ilnrho) = log(init_pressure)-log(Rgas) & -f(l1:l2,m1:m2,n1:n2,ilnTT)-log(mu1_full(l1:l2,m1:m2,n1:n2)) ! ! Initialize velocity ! ! f(l1:l2,m1:m2,n1:n2,iux)=exp(log_inlet_density - f(l1:l2,m1:m2,n1:n2,ilnrho)) & ! * (f(l1:l2,m1:m2,n1:n2,iux)+init_ux) !f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux)+init_ux ! ! Renormalize all species to be sure that the sum of all mass fractions ! are unity ! do i = 1,mx do j = 1,my do k = 1,mz f(i,j,k,ichemspec) = f(i,j,k,ichemspec)/sum(f(i,j,k,ichemspec)) enddo enddo enddo ! endsubroutine flame_front_new !*********************************************************************** subroutine TTD(f) ! ! 15-may-03/Nils Erland L. Haugen: adapted from flame_front ! real, dimension(mx,my,mz,mfarray) :: f integer :: i, j,k ! real :: initial_mu1, ksi_TTD, dTdr_c, deltaT, PP ! call air_field(f,PP) ! ksi_TTD = 1. dTdr_c = 2000 ![K/m] ! dTdr_c=20000 ![K/m] ! ! Initialize temperature ! do k = l1,l2 if (ltemperature_nolog) then f(k,:,:,iTT) = f(k,:,:,iTT)+ksi_TTD*dTdr_c*(xyz1(1)-x(k))/100. else deltaT = ksi_TTD*dTdr_c*(xyz1(1)-x(k))/100. f(k,:,:,ilnTT) = log(exp(f(k,:,:,ilnTT))+deltaT) ! print*,'deltaT=',deltaT, exp(f(k,m1,n1,ilnTT)), f(k,m1,n1,ilnTT) endif enddo ! ! Initialize density ! call getmu_array(f,mu1_full) f(l1:l2,m1:m2,n1:n2,ilnrho) = log(PP)-log(Rgas) & -f(l1:l2,m1:m2,n1:n2,ilnTT)-log(mu1_full(l1:l2,m1:m2,n1:n2)) ! ! Initialize velocity ! f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux)+init_ux ! ! Check if we want nolog of density or nolog of temperature ! if (ldensity_nolog) & f(l1:l2,m1:m2,n1:n2,irho) = exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) if (ltemperature_nolog) & f(l1:l2,m1:m2,n1:n2,iTT) = exp(f(l1:l2,m1:m2,n1:n2,ilnTT)) ! ! Renormalize all species to be sure that the sum of all mass fractions ! are unity ! do i = 1,mx do j = 1,my do k = 1,mz f(i,j,k,ichemspec) = f(i,j,k,ichemspec)/sum(f(i,j,k,ichemspec)) enddo enddo enddo ! endsubroutine TTD !*********************************************************************** subroutine triple_flame(f) ! ! 26-jul-10/Julien Savre: Copy from the flame_front case ! real, dimension(mx,my,mz,mfarray) :: f integer :: i, j,k ! real :: mO2=0., mH2=0., mN2=0., mH2O=0., mCH4=0., mCO2=0. real :: del, PP integer :: i_H2=0, i_O2=0, i_H2O=0, i_N2=0 integer :: ichem_H2=0, ichem_O2=0, ichem_N2=0, ichem_H2O=0 integer :: i_CH4=0, i_CO2=0, ichem_CH4=0, ichem_CO2=0 real :: final_massfrac_O2, final_massfrac_CH4, & final_massfrac_H2O, final_massfrac_CO2 real :: init_H2, init_O2, init_N2, init_H2O, init_CO2, init_CH4 real :: beta real :: init_y1, init_y2 real :: init_rho, init_m1 real, dimension(ny) :: dim logical :: lH2=.false., lO2=.false., lN2=.false., lH2O=.false. logical :: lCH4=.false., lCO2=.false. ! ltriple_flame = .true. ! call air_field(f,PP) ! init_y1 = xyz0(2) + Lxyz(2)/3. init_y2 = xyz0(2) + 2.*Lxyz(2)/3. ! if (ltemperature_nolog) f(:,:,:,ilnTT) = log(f(:,:,:,ilnTT)) if (lroot) print*, 'init_chem: triple_flame ' ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,lH2) if (lH2) then mH2 = species_constants(ichem_H2,imass) init_H2 = initial_massfractions(ichem_H2) endif call find_species_index('O2',i_O2,ichem_O2,lO2) if (lO2) then mO2 = species_constants(ichem_O2,imass) init_O2 = initial_massfractions(ichem_O2) endif call find_species_index('N2',i_N2,ichem_N2,lN2) if (lN2) then mN2 = species_constants(ichem_N2,imass) init_N2 = initial_massfractions(ichem_N2) endif call find_species_index('H2O',i_H2O,ichem_H2O,lH2O) if (lH2O) then mH2O = species_constants(ichem_H2O,imass) init_H2O = initial_massfractions(ichem_H2O) endif call find_species_index('CH4',i_CH4,ichem_CH4,lCH4) if (lCH4) then mCH4 = species_constants(ichem_CH4,imass) init_CH4 = initial_massfractions(ichem_CH4) endif call find_species_index('CO2',i_CO2,ichem_CO2,lCO2) if (lCO2) then mCO2 = species_constants(ichem_CO2,imass) init_CO2 = initial_massfractions(ichem_CO2) endif ! ! Find approximate value for the mass fraction of O2 after the flame front ! final_massfrac_O2 = 0. if (lH2) then final_massfrac_H2O = mH2O/(2.*mH2) * init_H2 final_massfrac_O2 = 1. - final_massfrac_H2O - init_N2 elseif (lCH4) then final_massfrac_CH4 = 0. final_massfrac_H2O = 2.*mH2O/mCH4 * init_CH4 final_massfrac_CO2 = mCO2/mCH4 * init_CH4 final_massfrac_O2 = & 1. - final_massfrac_CO2 - final_massfrac_H2O & - init_N2 endif ! if (final_massfrac_O2 < 0.) final_massfrac_O2 = 0. if (lroot) then print*, ' init final' if (lH2) print*, 'H2 :', init_H2, 0. if (lCH4) print*, 'CH4 :', init_CH4, 0. if (lO2) print*, 'O2 :', init_O2, final_massfrac_O2 if (lH2O) print*, 'H2O :', 0., final_massfrac_H2O if (lCO2) print*, 'CO2 :', 0., final_massfrac_CO2 endif ! ! Initialize temperature and species ! if (lT_tanh) then if (lroot) print*, 'Temperature initialization: tanh function.' else if (lroot) print*, 'Temperature initialization: linear.' endif ! do k = 1,mx ! ! Initialize temperature ! if (lT_tanh) then del = init_x2-init_x1 f(k,:,:,ilnTT) = f(k,:,:,ilnTT)+log((init_TT2+init_TT1)*0.5 & +((init_TT2-init_TT1)*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del))) else if (x(k) <= init_x1) then f(k,:,:,ilnTT) = log(init_TT1) endif if (x(k) >= init_x2) then f(k,:,:,ilnTT) = log(init_TT2) endif if (x(k) > init_x1 .and. x(k) < init_x2) then f(k,:,:,ilnTT) = & log((x(k)-init_x1)/(init_x2-init_x1) & *(init_TT2-init_TT1)+init_TT1) endif endif ! ! Initialize fuel ! if (lT_tanh) then if (lH2) then del = (init_x2-init_x1) f(k,:,:,i_H2) = (0.+f(l1,:,:,i_H2))*0.5 & +(0.-f(l1,:,:,i_H2))*0.5 & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) endif if (lCH4) then if (lroot) print*, 'No tanh initial function available for CH4 combustion.' endif ! else if (x(k) > init_x1) then if (lH2) then f(k,:,:,i_H2) = init_H2*(exp(f(k,:,:,ilnTT))-init_TT2) & /(init_TT1-init_TT2) endif if (lCH4) then f(k,:,:,i_CH4) = init_CH4*(exp(f(k,:,:,ilnTT))-init_TT2) & /(init_TT1-init_TT2) endif endif endif ! ! Initialize oxygen ! if (lT_tanh) then del = (init_x2-init_x1) f(k,:,:,i_O2) = (f(l2,:,:,i_O2)+f(l1,:,:,i_O2))*0.5 & +((f(l2,:,:,i_O2)-f(l1,:,:,i_O2))*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) else ! if (x(k) > init_x2) f(k,:,:,i_O2) = final_massfrac_O2 if (x(k) > init_x1 .and. x(k) <= init_x2) & f(k,:,:,i_O2) = (x(k)-init_x1)/(init_x2-init_x1) & *(final_massfrac_O2-init_O2)+init_O2 endif enddo ! ! Initialize products ! if (lT_tanh) then do k = 1,mx if (lH2) then del = (init_x2-init_x1) f(k,:,:,i_H2O) = (f(l1,:,:,i_H2)/2.*18.+f(l1,:,:,i_H2O))*0.5 & +((f(l1,:,:,i_H2)/2.*18.-f(l1,:,:,i_H2O))*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) endif if (lCH4) then if (lroot) print*, 'No tanh initial function available for CH4 combustion.' endif enddo else do k = 1,mx if (x(k) >= init_x1 .and. x(k) < init_x2) then f(k,:,:,i_H2O) = (x(k)-init_x1)/(init_x2-init_x1) & *final_massfrac_H2O if (lCO2) f(k,:,:,i_CO2) = (x(k)-init_x1)/(init_x2-init_x1) & *final_massfrac_CO2 elseif (x(k) >= init_x2) then if (lCO2) f(k,:,:,i_CO2) = final_massfrac_CO2 if (lH2O) f(k,:,:,i_H2O) = final_massfrac_H2O endif enddo endif ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif ! ! Set the initial equivalence ratio gradient in the fresh gases ! dim(1:ny) = y(m1:m2) beta = init_N2/init_O2 do i = 1, mx if (x(i) <= init_x1) then do j = 1, ny if (dim(j) >= init_y1 .and. dim(j) <= init_y2) then if (lH2) then f(i,m1+j-1,:,i_H2) = init_zz2 - (init_zz2-init_zz1) * (init_y2 - dim(j)) / & (init_y2 - init_y1) elseif (lCH4) then f(i,m1+j-1,:,i_CH4) = init_zz2 - (init_zz2-init_zz1) * (init_y2 - dim(j)) / & (init_y2 - init_y1) endif elseif (dim(j) <= init_y1) then if (lH2) then f(i,m1+j-1,:,i_H2) = init_zz1 elseif (lCH4) then f(i,m1+j-1,:,i_CH4) = init_zz1 endif elseif (dim(j) >= init_y2) then if (lH2) then f(i,m1+j-1,:,i_H2) = init_zz2 elseif (lCH4) then f(i,m1+j-1,:,i_CH4) = init_zz2 endif endif enddo ! if (lH2) then f(i,ny:my,:,i_H2) = init_zz2 if (lO2) f(i,:,:,i_O2) = (1.-f(i,:,:,i_H2)) / (1.+beta) if (lN2) f(i,:,:,i_N2) = 1.-f(i,:,:,i_O2)-f(i,:,:,i_H2) elseif (lCH4) then f(i,ny:my,:,i_CH4) = init_zz2 if (lO2) f(i,:,:,i_O2) = (1.-f(i,:,:,i_CH4)) / (1.+beta) if (lN2) f(i,:,:,i_N2) = 1.-f(i,:,:,i_O2)-f(i,:,:,i_CH4) if (lCO2) f(i,:,:,i_CO2) = 0. if (lH2O) f(i,:,:,i_H2O) = 0. endif endif enddo ! ! Renormalize all species to be sure that the sum of all mass fractions ! are unity ! do i = 1,mx do j = 1,my do k = 1,mz f(i,j,k,ichemspec) = f(i,j,k,ichemspec)/sum(f(i,j,k,ichemspec)) enddo enddo enddo ! ! Initialize density ! call getmu_array(f,mu1_full) f(l1:l2,m1:m2,n1:n2,ilnrho) = log(init_pressure)-log(Rgas) & -f(l1:l2,m1:m2,n1:n2,ilnTT)-log(mu1_full(l1:l2,m1:m2,n1:n2)) ! ! Initialize velocity ! if (lCH4) then init_m1 = init_CH4/mCH4 + init_O2/mO2 + init_N2/mN2 init_rho = init_pressure/(init_TT1 * init_m1 * Rgas) f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux) + & init_ux * init_rho / exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) else f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux)+init_ux endif ! ! Check if we want nolog of density or nolog of temperature ! if (ldensity_nolog) & f(l1:l2,m1:m2,n1:n2,irho) = exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) if (ltemperature_nolog) & f(l1:l2,m1:m2,n1:n2,iTT) = exp(f(l1:l2,m1:m2,n1:n2,ilnTT)) ! endsubroutine triple_flame !*********************************************************************** subroutine flame(f) ! ! 05-jun-09/Nils Erland L. Haugen: adapted from similar ! routine in special/chem_stream.f90 ! This routine set up the initial profiles used in 1D flame speed measurments ! NILS: This routine is essentially the samw as flame_front, but I leave ! NILS: flame_front as it is for now for backwards compatibility. ! real, dimension(mx,my,mz,mfarray) :: f integer :: i, j,k ! real :: mO2, mH2, mN2, mH2O, mCH4, mCO2 real :: log_inlet_density, del, PP integer :: i_H2, i_O2, i_H2O, i_N2, ichem_H2, ichem_O2, ichem_N2, ichem_H2O integer :: i_CH4, i_CO2, ichem_CH4, ichem_CO2 real :: initial_mu1, final_massfrac_O2 real :: init_H2, init_O2, init_N2, init_H2O, init_CO2, init_CH4 logical :: lH2, lO2, lN2, lH2O, lCH4, lCO2 ! lflame_front = .true. ! call air_field(f,PP) ! if (ltemperature_nolog) f(:,:,:,ilnTT) = log(f(:,:,:,ilnTT)) ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,lH2) if (lH2) then mH2 = species_constants(ichem_H2,imass) init_H2 = initial_massfractions(ichem_H2) endif call find_species_index('O2',i_O2,ichem_O2,lO2) if (lO2) then mO2 = species_constants(ichem_O2,imass) init_O2 = initial_massfractions(ichem_O2) endif call find_species_index('N2',i_N2,ichem_N2,lN2) if (lN2) then mN2 = species_constants(ichem_N2,imass) init_N2 = initial_massfractions(ichem_N2) endif call find_species_index('H2O',i_H2O,ichem_H2O,lH2O) if (lH2O) then mH2O = species_constants(ichem_H2O,imass) init_H2O = initial_massfractions(ichem_H2O) endif call find_species_index('CH4',i_CH4,ichem_CH4,lCH4) if (lCH4) then mCH4 = species_constants(ichem_CH4,imass) init_CH4 = initial_massfractions(ichem_CH4) endif call find_species_index('CO2',i_CO2,ichem_CO2,lCO2) if (lCO2) then mCO2 = species_constants(ichem_CO2,imass) init_CO2 = initial_massfractions(ichem_CO2) endif ! ! Find approximate value for the mass fraction of O2 after the flame front ! final_massfrac_O2 = (init_O2/mO2 -init_H2/(2.*mH2) -init_CH4*2/(mCH4))*mO2 ! if (final_massfrac_O2 < 0.) final_massfrac_O2 = 0. ! ! Initialize temperature and species ! do k = 1,mx ! ! Initialize temperature ! if (lT_tanh) then del = init_x2-init_x1 f(k,:,:,ilnTT) = f(k,:,:,ilnTT)+log((init_TT2+init_TT1)*0.5 & +((init_TT2-init_TT1)*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del))) else if (x(k) <= init_x1) then f(k,:,:,ilnTT) = log(init_TT1) endif if (x(k) >= init_x2) then f(k,:,:,ilnTT) = log(init_TT2) endif if (x(k) > init_x1 .and. x(k) < init_x2) then f(k,:,:,ilnTT) = & log((x(k)-init_x1)/(init_x2-init_x1) & *(init_TT2-init_TT1)+init_TT1) endif endif ! ! Initialize steam and hydrogen ! if (lT_tanh) then del = (init_x2-init_x1) if (lH2) then f(k,:,:,i_H2) = (0.+f(l1,:,:,i_H2))*0.5 & +(0.-f(l1,:,:,i_H2))*0.5 & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) f(k,:,:,i_H2O) = (f(l1,:,:,i_H2)/2.*18.+f(l1,:,:,i_H2O))*0.5 & +((f(l1,:,:,i_H2)/2.*18.-f(l1,:,:,i_H2O))*0.5) & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) endif if (lCH4) then f(k,:,:,i_CH4) = init_CH4*0.5 & -init_CH4*0.5 & *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) f(k,:,:,i_H2O) = init_H2O+(init_CH4-f(k,:,:,i_CH4))*2.*mH2O/mCH4 f(k,:,:,i_CO2) = init_CO2+(init_CH4-f(k,:,:,i_CH4))*1.*mCO2/mCH4 f(k,:,:,i_O2) = init_O2 -(init_CH4-f(k,:,:,i_CH4))*2.*mO2 /mCH4 endif else if (x(k) > init_x1) then if (lH2) then f(k,:,:,i_H2) = init_H2*(exp(f(k,:,:,ilnTT))-init_TT2) & /(init_TT1-init_TT2) endif if (lCH4) then f(k,:,:,i_CH4) = init_CH4*(exp(f(k,:,:,ilnTT))-init_TT2) & /(init_TT1-init_TT2) endif endif endif ! ! Initialize oxygen ! if (lT_tanh) then ! $ del=(init_x2-init_x1) ! $ f(k,:,:,i_O2)=(f(l2,:,:,i_O2)+f(l1,:,:,i_O2))*0.5 & ! $ +((f(l2,:,:,i_O2)-f(l1,:,:,i_O2))*0.5) & ! $ *(exp(x(k)/del)-exp(-x(k)/del))/(exp(x(k)/del)+exp(-x(k)/del)) else if (x(k) > init_x2) then f(k,:,:,i_O2) = final_massfrac_O2 endif if (x(k) > init_x1 .and. x(k) < init_x2) then f(k,:,:,i_O2) = (x(k)-init_x2)/(init_x1-init_x2) & *(init_O2-final_massfrac_O2) & +final_massfrac_O2 endif endif enddo ! ! Initialize steam and CO2 ! if (.not. lT_tanh) then do k = 1,mx if (x(k) >= init_x1) then if (final_massfrac_O2 > 0.) then f(k,:,:,i_H2O) = initial_massfractions(ichem_H2)/mH2*mH2O & *(exp(f(k,:,:,ilnTT))-init_TT1) & /(init_TT2-init_TT1)+init_H2O else if (x(k) >= init_x2) then if (lCO2) f(k,:,:,i_CO2) = init_CO2+(init_CH4-f(k,:,:,i_CH4)) f(k,:,:,i_H2O) = 1.-f(k,:,:,i_N2)-f(k,:,:,i_H2)-f(k,:,:,i_O2) if (lCH4) f(k,:,:,i_H2O) = f(k,:,:,i_H2O)-f(k,:,:,i_CH4) if (lCO2) f(k,:,:,i_H2O) = f(k,:,:,i_H2O)-f(k,:,:,i_CO2) else f(k,:,:,i_H2O) = (x(k)-init_x1)/(init_x2-init_x1) & *((1.-f(l2,:,:,i_N2)-f(l2,:,:,i_H2)) & -initial_massfractions(ichem_H2O)) & +initial_massfractions(ichem_H2O) endif endif endif enddo endif ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif ! ! Find logaritm of density at inlet ! initial_mu1 & = initial_massfractions(ichem_H2)/(mH2) & +initial_massfractions(ichem_O2)/(mO2) & +initial_massfractions(ichem_H2O)/(mH2O) & +initial_massfractions(ichem_N2)/(mN2) if (lCO2) initial_mu1 = initial_mu1+init_CO2/(mCO2) if (lCH4) initial_mu1 = initial_mu1+init_CH4/(mCH4) log_inlet_density = & log(init_pressure)-log(Rgas)-log(init_TT1)-log(initial_mu1) ! call getmu_array(f,mu1_full) ! ! Initialize density ! f(l1:l2,m1:m2,n1:n2,ilnrho) = log(init_pressure)-log(Rgas) & -f(l1:l2,m1:m2,n1:n2,ilnTT)-log(mu1_full(l1:l2,m1:m2,n1:n2)) ! ! Initialize velocity ! f(l1:l2,m1:m2,n1:n2,iux) = f(l1:l2,m1:m2,n1:n2,iux)+init_ux ! ! Check if we want nolog of density or nolog of temperature ! if (ldensity_nolog) & f(l1:l2,m1:m2,n1:n2,irho) = exp(f(l1:l2,m1:m2,n1:n2,ilnrho)) if (ltemperature_nolog) & f(l1:l2,m1:m2,n1:n2,iTT) = exp(f(l1:l2,m1:m2,n1:n2,ilnTT)) ! ! Renormalize all species too be sure that the sum of all mass fractions ! are unity ! do i = 1,mx do j = 1,my do k = 1,mz f(i,j,k,ichemspec) = f(i,j,k,ichemspec)/sum(f(i,j,k,ichemspec)) enddo enddo enddo ! endsubroutine flame !*********************************************************************** subroutine flame_blob(f) ! real, dimension(mx,my,mz,mfarray) :: f integer :: j1, j2, j3 ! real :: mO2, mH2, mN2, mH2O, PP integer :: i_H2, i_O2, i_H2O, i_N2, ichem_H2, ichem_O2, ichem_N2, ichem_H2O real :: initial_mu1, final_massfrac_O2 logical :: found_specie ! real :: Rad, sz1, sz2 ! lflame_front = .true. ! call air_field(f,PP) ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,found_specie) call find_species_index('O2',i_O2,ichem_O2,found_specie) call find_species_index('N2',i_N2,ichem_N2,found_specie) call find_species_index('H2O',i_H2O,ichem_H2O,found_specie) mO2 = species_constants(ichem_O2,imass) mH2 = species_constants(ichem_H2,imass) mH2O = species_constants(ichem_H2O,imass) mN2 = species_constants(ichem_N2,imass) ! ! Find approximate value for the mass fraction of O2 after the flame front ! final_massfrac_O2 & = (initial_massfractions(ichem_O2)/mO2 & -initial_massfractions(ichem_H2)/(2*mH2))*mO2 ! ! Initialize temperature and species in air_field(f) ! if (unit_system == 'cgs') then Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy endif ! ! Find logaritm of density at inlet ! initial_mu1 & = initial_massfractions(ichem_H2)/(mH2) & +initial_massfractions(ichem_O2)/(mO2) & +initial_massfractions(ichem_H2O)/(mH2O) & +initial_massfractions(ichem_N2)/(mN2) ! call getmu_array(f,mu1_full) ! do j3 = 1,mz do j2 = 1,my do j1 = 1,mx ! Rad = 0. if (nxgrid > 1) then Rad = x(j1)*x(j1) endif if (nygrid > 1) then Rad = Rad+y(j2)*y(j2) endif if (nzgrid > 1) then Rad = Rad+z(j3)*z(j3) endif ! Rad = sqrt(Rad) ! f(j1,j2,j3,ilnTT) = log((init_TT2-init_TT1)*exp(-(Rad/init_x2)**2)+init_TT1) mu1_full(j1,j2,j3) = f(j1,j2,j3,i_H2)/(mH2)+f(j1,j2,j3,i_O2)/(mO2) & +f(j1,j2,j3,i_H2O)/(mH2O)+f(j1,j2,j3,i_N2)/(mN2) ! f(j1,j2,j3,ilnrho) = log(init_pressure)-log(Rgas)-f(j1,j2,j3,ilnTT) & -log(mu1_full(j1,j2,j3)) ! ! Initialize velocity ! f(j1,j2,j3,iux) = f(j1,j2,j3,iux) & +init_ux!*exp(log_inlet_density)/exp(f(j1,j2,j3,ilnrho)) f(j1,j2,j3,iuy) = f(j1,j2,j3,iuy)+ init_uy f(j1,j2,j3,iuz) = f(j1,j2,j3,iuz)+ init_uz ! if (nxgrid == 1) f(j1,j2,j3,iux) = 0. if (nygrid == 1) f(j1,j2,j3,iuy) = 0. if (nzgrid == 1) f(j1,j2,j3,iuz) = 0. ! if (nxgrid /= 1) then sz1 = (xyz0(1)+Lxyz(1)*0.15) sz2 = (xyz0(1)+Lxyz(1)*(1.-0.15)) endif ! if (nygrid /= 1) then sz1 = (xyz0(2)+Lxyz(2)*0.15) sz2 = (xyz0(2)+Lxyz(2)*(1.-0.15)) endif ! if (nzgrid /= 1) then sz1 = (xyz0(3)+Lxyz(3)*0.15) sz2 = (xyz0(3)+Lxyz(3)*(1.-0.15)) endif ! enddo enddo enddo ! ! Check if we want nolog of density ! if (ldensity_nolog) f(:,:,:,irho) = exp(f(:,:,:,ilnrho)) ! endsubroutine flame_blob !*********************************************************************** subroutine opposite_flames(f) ! ! 03-jan-10/nilshau: adapted from opposite_ignitions ! ! Set up two oppositely directed flame fronts in the x-direction. ! The two fronts have fresh gas between them. ! real, dimension(mx,my,mz,mfarray) :: f integer :: j1, j2, j3 ! real :: mO2, mH2, mN2, mH2O, lower, upper, PP integer :: i_H2, i_O2, i_H2O, i_N2, ichem_H2, ichem_O2, ichem_N2, ichem_H2O integer :: i_C3H8, ichem_C3H8, i_CO2, ichem_CO2 real :: final_massfrac_O2, mu1, phi, delta_O2, mC3H8, mCO2 logical :: found_specie, lH2, lCO2, lC3H8, lH2O real :: norm, flat_range ! lflame_front = .true. ! call air_field(f,PP) ! ! Initialize some indexes ! call find_species_index('H2',i_H2,ichem_H2,lH2) call find_species_index('O2',i_O2,ichem_O2,found_specie) call find_species_index('N2',i_N2,ichem_N2,found_specie) call find_species_index('H2O',i_H2O,ichem_H2O,lH2O) call find_species_index('C3H8',i_C3H8,ichem_C3H8,lC3H8) call find_species_index('CO2',i_CO2,ichem_CO2,lCO2) mO2 = species_constants(ichem_O2,imass) mH2O = species_constants(ichem_H2O,imass) mN2 = species_constants(ichem_N2,imass) if (lC3H8) mC3H8 = species_constants(ichem_C3H8,imass) if (lH2) mH2 = species_constants(ichem_H2,imass) if (lCO2) mCO2 = species_constants(ichem_CO2,imass) ! ! Find approximate value for the mass fraction of O2 after the flame front ! final_massfrac_O2 = initial_massfractions(ichem_O2) if (lH2) then delta_O2 = initial_massfractions(ichem_H2)/(2*mH2)*mO2 else delta_O2 = 0. endif final_massfrac_O2 = final_massfrac_O2-delta_O2 ! if (lC3H8) then delta_O2 = 5*initial_massfractions(ichem_C3H8)/mC3H8*mO2 else delta_O2 = 0. endif final_massfrac_O2 = final_massfrac_O2-delta_O2 ! if (ltemperature_nolog) call fatal_error('opposite_flames', & 'only implemented for ltemperature_nolog=F') ! ! Loop over all grid points ! do j3 = 1,mz do j2 = 1,my do j1 = 1,mx ! ! First define the distance from the lower and upper domain boundary. ! lower = x(j1)-xyz0(1) upper = xyz1(1)-x(j1) ! ! Find progress variable phi based on distance from boundaries. ! flat_range = init_x2*0.1 phi & = exp(-((lower-flat_range)/init_x2)**2) & +exp(-((upper-flat_range)/init_x2)**2) if (phi > 1.0) phi = 1. if (flat_range > lower) phi = 1. if (flat_range > upper) phi = 1. ! ! Find temperature, species and density based on progress variable ! f(j1,j2,j3,ilnTT) = log(init_TT1+phi*(init_TT2-init_TT1)) if (lH2) then f(j1,j2,j3,i_H2) = (1-phi)*initial_massfractions(ichem_H2) f(j1,j2,j3,i_H2O) = phi*initial_massfractions(ichem_H2)*mH2O/mH2 endif if (lC3H8) then f(j1,j2,j3,i_C3H8) = (1-phi)*initial_massfractions(ichem_C3H8) f(j1,j2,j3,i_H2O) = 4*phi*initial_massfractions(ichem_C3H8)*mH2O/mC3H8 f(j1,j2,j3,i_CO2) = 3*phi*initial_massfractions(ichem_C3H8)*mCO2/mC3H8 endif f(j1,j2,j3,i_O2) = (1-phi)*(initial_massfractions(ichem_O2) & -final_massfrac_O2)+final_massfrac_O2 ! ! Re-normalize mass fractions ! norm = f(j1,j2,j3,i_O2)+f(j1,j2,j3,i_N2) if (lC3H8) norm = norm + f(j1,j2,j3,i_C3H8) if (lCO2) norm = norm + f(j1,j2,j3,i_CO2) if (lH2O) norm = norm + f(j1,j2,j3,i_H2O) if (lH2) norm = norm + f(j1,j2,j3,i_H2) f(j1,j2,j3,minval(ichemspec):maxval(ichemspec)) & = f(j1,j2,j3,minval(ichemspec):maxval(ichemspec))/norm ! ! Find mean molecular weight and density ! mu1 & = f(j1,j2,j3,i_O2 )/mO2 & +f(j1,j2,j3,i_H2O)/mH2O & +f(j1,j2,j3,i_N2 )/mN2 if (lH2) mu1 = mu1+f(j1,j2,j3,i_H2 )/mH2 if (lCO2) mu1 = mu1+f(j1,j2,j3,i_CO2 )/mCO2 if (lC3H8) mu1 = mu1+f(j1,j2,j3,i_C3H8)/mC3H8 f(j1,j2,j3,ilnrho) = log(init_pressure)-log(Rgas)-f(j1,j2,j3,ilnTT) & -log(mu1) enddo enddo enddo ! ! Check if we want nolog of density ! if (ldensity_nolog) f(:,:,:,irho) = exp(f(:,:,:,ilnrho)) if (ltemperature_nolog) f(:,:,:,iTT) = exp(f(:,:,:,ilnTT)) ! endsubroutine opposite_flames !*********************************************************************** subroutine opposite_ignitions(f) ! ! 03-jan-10/nilshau: adapted from flame_blob ! ! Set up two oppositely directed flame fronts in the x-direction. ! The two fronts have fresh gas between them. ! real, dimension(mx,my,mz,mfarray) :: f integer :: j1, j2, j3 ! real :: lower, upper, PP real :: T0, T1, rho0, rho1 ! lflame_front = .true. ! call air_field(f,PP) ! ! Initialize some indexes ! do j3 = n1,n2 do j2 = m1,m2 do j1 = l1,l2 ! ! First define the distance from the lower and upper domain boundary. ! lower = x(j1)-xyz0(1) upper = xyz1(1)-x(j1) ! ! Find temperature and density based on distance from boundaries. ! if (ltemperature_nolog) then T0 = f(j1,j2,j3,ilnTT) else T0 = exp(f(j1,j2,j3,ilnTT)) endif if (ldensity_nolog) then rho0 = f(j1,j2,j3,ilnrho) else rho0 = exp(f(j1,j2,j3,ilnrho)) endif T1 = & (init_TT2-init_TT1)*exp(-(lower/init_x2)**2)+ & (init_TT2-init_TT1)*exp(-(upper/init_x2)**2)+ & init_TT1 if (ltemperature_nolog) then f(j1,j2,j3,ilnTT) = T1 else f(j1,j2,j3,ilnTT) = log(T1) endif rho1 = rho0*T0/T1 if (ldensity_nolog) then f(j1,j2,j3,ilnrho) = rho1 else f(j1,j2,j3,ilnrho) = log(rho1) endif enddo enddo enddo ! endsubroutine opposite_ignitions !*********************************************************************** subroutine calc_for_chem_mixture(f) ! ! Calculate quantities for a mixture ! ! 22-jun-10/julien: Added evaluation of diffusion coefficients using constant ! Lewis numers Di = lambda/(rho*Cp*Lei) ! 10-jan-11/julien: Modified for a resolution with LSODE ! 26-jui-11/julien: Replaced fatal_error by inevitably_fatal_error to allow ! proper exit when T_locT_up ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx) :: tmp_sum, tmp_sum2, nu_dyn, nuk_nuj, Phi real, dimension(mx) :: cp_R_spec, T_loc, T_loc_2, T_loc_3, T_loc_4 ! integer :: k, j, j2, j3 real :: T_up, T_mid, T_low real :: mk_mj real :: EE=0., TT=0., yH=1. ! logical, save :: lwrite=.true. ! character(len=fnlen) :: output_file="./data/mix_quant.out" character(len=15) :: writeformat integer :: file_id=123 integer :: ii1=1, ii2=2, ii3=3, ii4=4, ii5=5 ! ! Density and temperature ! ! call timing('calc_for_chem_mixture','entered') call getdensity(f,EE,TT,yH,rho_full) call gettemperature(f,TT_full) ! ! Now this routine is only for chemkin data !!! ! if (lcheminp) then ! if (unit_system == 'cgs') then ! Rgas_unit_sys = k_B_cgs/m_u_cgs Rgas = Rgas_unit_sys/unit_energy ! call getmu_array(f,mu1_full) ! if (l1step_test) then species_constants(:,imass) = 1. mu1_full = 1. endif ! ! Mole fraction XX ! do k = 1,nchemspec if (species_constants(k,imass) > 0.) then do j2 = mm1,mm2 do j3 = nn1,nn2 XX_full(:,j2,j3,k) = f(:,j2,j3,ichemspec(k))*unit_mass & /(species_constants(k,imass)*mu1_full(:,j2,j3)) enddo enddo endif enddo ! ! do m=m1,m2 ! do n=n1,n2 ! call getpressure(pp_full(l1:l2,m,n),TT_full(l1:l2,m,n),& ! rho_full(l1:l2,m,n),mu1_full(l1:l2,m,n)) ! enddo ! enddo ! ! Specific heat at constant pressure ! if ((Cp_const < impossible) .or. (Cv_const < impossible)) then ! if (Cp_const < impossible) then cp_full = Cp_const*mu1_full cv_full = (Cp_const-Rgas)*mu1_full endif ! if (Cv_const < impossible) then cv_full = Cv_const*mu1_full endif else cp_full = 0. cv_full = 0. do j3 = nn1,nn2 do j2 = mm1,mm2 T_loc = TT_full(:,j2,j3) T_loc_2 = T_loc*T_loc T_loc_3 = T_loc_2*T_loc T_loc_4 = T_loc_3*T_loc do k = 1,nchemspec if (species_constants(k,imass) > 0.) then T_low = species_constants(k,iTemp1)-1. T_mid = species_constants(k,iTemp2) T_up = species_constants(k,iTemp3) ! ! $ if (j1<=l1 .or. j2>=l2) then ! $ T_low=0. ! $ T_up=1e10 ! $ endif ! if (j2 <= m1 .or. j2 >= m2) then T_low = 0. T_up = 1e10 endif ! if (j3 <= n1 .or. j3 >= n2) then T_low = 0. T_up = 1e10 endif ! if (lcloud) T_low = 20. where (T_loc >= T_low .and. T_loc <= T_mid) cp_R_spec = species_constants(k,iaa2(ii1)) & +species_constants(k,iaa2(ii2))*T_loc & +species_constants(k,iaa2(ii3))*T_loc_2 & +species_constants(k,iaa2(ii4))*T_loc_3 & +species_constants(k,iaa2(ii5))*T_loc_4 elsewhere (T_loc >= T_mid .and. T_loc <= T_up) cp_R_spec = species_constants(k,iaa1(ii1)) & +species_constants(k,iaa1(ii2))*T_loc & +species_constants(k,iaa1(ii3))*T_loc_2 & +species_constants(k,iaa1(ii4))*T_loc_3 & +species_constants(k,iaa1(ii5))*T_loc_4 endwhere cv_R_spec_full(:,j2,j3,k) = cp_R_spec-1. ! ! Check if the temperature are within bounds ! if (maxval(T_loc) > T_up .or. minval(T_loc) < T_low) then print*,'iproc=',iproc print*,'TT_full(:,j2,j3)=',T_loc print*,'j2,j3=',j2,j3 call inevitably_fatal_error('calc_for_chem_mixture', & 'TT_full(:,j2,j3) is outside range', .true.) endif ! ! Find cp and cv for the mixture for the full domain ! cp_full(:,j2,j3) = cp_full(:,j2,j3)+f(:,j2,j3,ichemspec(k)) & *cp_R_spec/species_constants(k,imass)*Rgas cv_full(:,j2,j3) = cv_full(:,j2,j3)+f(:,j2,j3,ichemspec(k)) & *cv_R_spec_full(:,j2,j3,k)/species_constants(k,imass)*Rgas cp_spec_glo(:,j2,j3,k)=cp_R_spec/species_constants(k,imass)*Rgas endif enddo enddo enddo endif ! ! All the transport properties are calculated only if we are not using LSODE ! to solve chemistry or during the transport substep ! if (.not. lchemonly) then ! ! Viscosity of a mixture ! if (tran_exist) then call calc_diff_visc_coef(f) endif ! if (visc_const == impossible) then do j3 = nn1,nn2 do j2 = mm1,mm2 ! if (lone_spec) then f(:,j2,j3,iviscosity) = species_viscosity(:,j2,j3,1)/rho_full(:,j2,j3) else nu_dyn = 0. do k = 1,nchemspec if (species_constants(k,imass) > 0.) then tmp_sum2 = 0. do j = 1,nchemspec mk_mj = species_constants(k,imass)/species_constants(j,imass) nuk_nuj = species_viscosity(:,j2,j3,k)/species_viscosity(:,j2,j3,j) Phi = 1./sqrt(8.)*1./sqrt(1.+mk_mj)*(1.+sqrt(nuk_nuj)*mk_mj**(-0.25))**2 tmp_sum2 = tmp_sum2 + XX_full(:,j2,j3,j)*Phi enddo nu_dyn = nu_dyn + XX_full(:,j2,j3,k)*species_viscosity(:,j2,j3,k)/tmp_sum2 endif enddo f(:,j2,j3,iviscosity) = nu_dyn/rho_full(:,j2,j3) endif ! enddo enddo endif ! ! Diffusion coefficient of a mixture from tran.dat file ! if ((.not. lDiff_simple).and.(.not. lDiff_lewis)) then ! do j3 = nn1,nn2 do j2 = mm1,mm2 ! Diff_full(:,j2,j3,:) = 0. if (.not. lone_spec) then ! if (lfix_Sc) then do k = 1,nchemspec if (species_constants(k,imass) > 0.) then Diff_full(:,j2,j3,k) = species_viscosity(:,j2,j3,k) & /rho_full(:,j2,j3)/Sc_number endif enddo elseif (ldiffusion) then ! ! The mixture diffusion coefficient as described in eq. 5-45 of the Chemkin ! manual. Previously eq. 5-44 was used, but due to problems in the limit ! when the mixture becomes a pure specie we changed to the more robust eq. 5-45. ! do k = 1,nchemspec tmp_sum = 0. tmp_sum2 = 0. do j = 1,nchemspec if (species_constants(k,imass) > 0.) then if (j /= k) then tmp_sum = tmp_sum & +XX_full(:,j2,j3,j)/Bin_Diff_coef(:,j2,j3,j,k) tmp_sum2 = tmp_sum2 & +XX_full(:,j2,j3,j)*species_constants(j,imass) ! endif endif enddo Diff_full_add(:,j2,j3,k) = mu1_full(:,j2,j3)*tmp_sum2/tmp_sum enddo endif endif enddo enddo endif ! ! Thermal diffusivity ! if (lheatc_chemistry .and. (.not. lThCond_simple)) then call calc_therm_diffus_coef(f) endif endif ! else call stop_it('This case works only for cgs units system!') endif endif ! ! Write block ! if (lwrite) then open (file_id,file=output_file) write (file_id,*) 'Mixture quantities' write (file_id,*) '*******************' write (file_id,*) '' write (file_id,*) 'Mass, g/mole' write (file_id,'(7E12.4)') 1./maxval(mu1_full/unit_mass) write (file_id,*) '' write (file_id,*) 'Density, g/cm^3' write (file_id,'(7E12.4)') rho_full(l1,m1,n1)*unit_mass/unit_length**3, & rho_full(l2,m2,n2)*unit_mass/unit_length**3 write (file_id,*) '' write (file_id,*) 'Themperature, K' ! Commented the next line out because ! samples/2d-tests/chemistry_GrayScott apparently has no f(:,:,:,5) if (ilnTT > 0) write (file_id,'(7E12.4)') & exp(f(l1,m1,n1,ilnTT))*unit_temperature, & exp(f(l2,m2,n2,ilnTT))*unit_temperature write (file_id,*) '' write (file_id,*) 'Cp, erg/mole/K' write (file_id,'(7E12.4)') cp_full(l1,m1,n1)/Rgas* & Rgas_unit_sys/mu1_full(l1,m1,n1)/unit_mass,cp_full(l2,m2,n2)/Rgas* & Rgas_unit_sys/mu1_full(l2,m2,n2)/unit_mass write (file_id,*) '' write (file_id,*) 'cp, erg/g/K' write (file_id,'(7E12.4)') cp_full(l1,m1,n1)/Rgas*Rgas_unit_sys,cp_full(l2,m2,n2)/Rgas*Rgas_unit_sys write (file_id,*) '' write (file_id,*) 'gamma,max,min' write (file_id,'(7E12.4)') cp_full(l1,m1,n1)/cv_full(l1,m1,n1), & cp_full(l2,m2,n2)/cv_full(l2,m2,n2) if (.not. lchemonly) then write (file_id,*) '' write (file_id,*) 'Species viscosity, g/cm/s,' do k = 1,nchemspec write (file_id,'(7E12.4)') species_viscosity(l1,m1,n1,k), & species_viscosity(l2-1,m1,n1,k) enddo write (file_id,*) '' write (file_id,*) 'Thermal cond, erg/(cm K s),' write (file_id,'(7E12.4)') (lambda_full(l1,m1,n1)* & unit_energy/unit_time/unit_length/unit_temperature), & (lambda_full(l2,m2,n2)* & unit_energy/unit_time/unit_length/unit_temperature) write (file_id,*) '' write (file_id,*) 'Species Diffusion coefficient, cm^2/s' if (.not. lDiff_simple) then do k = 1,nchemspec write (file_id,'(7E12.4)') & Diff_full_add(l1,m1,n1,k)*unit_length**2/unit_time, & Diff_full_add(l2,m2,n2,k)*unit_length**2/unit_time enddo endif if (lparticles_chemistry) then write (file_id,*) '' writeformat = '( E12.4)' write (writeformat(2:3),'(I2)') nchemspec write (file_id,*) 'Mass fraction, -' write (file_id,writeformat) f(l1,m1,n1,ichemspec(1):ichemspec(nchemspec)) write (file_id,*) '' write (file_id,writeformat) species_constants(:,imass) endif endif write (file_id,*) '' if (lroot) print*,'calc_for_chem_mixture: writing mix_quant.out file' close (file_id) lwrite = .false. endif ! call timing('calc_for_chem_mixture','finished') ! endsubroutine calc_for_chem_mixture !*********************************************************************** subroutine chemspec_normalization(f) ! ! 20-sep-10/Natalia: coded ! renormalization of the species ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz) :: sum_Y integer :: k ! sum_Y = 0. do k = 1,nchemspec sum_Y = sum_Y+f(:,:,:,ichemspec(k)) enddo do k = 1,nchemspec f(:,:,:,ichemspec(k)) = f(:,:,:,ichemspec(k))/sum_Y enddo ! endsubroutine chemspec_normalization !*********************************************************************** subroutine astrobiology_data(f) ! ! Proceedure to read in stoichiometric matrices in explicit format for ! forward and backward reations. For historical reasons this is referred ! to as "astrobiology_data". ! ! 28-feb-08/axel: coded ! character(len=80) :: chemicals='' ! Careful, limits the absolut size of the input matrix !!! character(len=15) :: file1='chemistry_m.dat', file2='chemistry_p.dat' ! character (len=fnlen) :: input_file='chem.inp' real, dimension(mx,my,mz,mfarray) :: f real :: dummy logical :: exist, exist1, exist2 integer :: i, j, stat integer :: nchemspectemp character :: tmpchar logical :: inside ! ! Find number of reactions by reading how many lines we have in file2 ! j = 1 open (19,file=file2) read (19,*) chemicals do while (.true.) read (19,*,end=996) dummy j = j+1 enddo 996 close (19) mreactions = j-1 if (lroot) print*,'Number of reactions=',mreactions ! ! Find number of compounds by reading how many columns we have in file1 ! open (19,file=file1) read (19,fmt="(a80)") chemicals nchemspectemp = 0 inside = .true. do i = 1,len_trim(chemicals) tmpchar = chemicals(i:i) if (tmpchar == ' ') then if (.not. inside) then inside = .true. nchemspectemp = nchemspectemp+1 endif else inside = .false. endif enddo if (inside) nchemspectemp = nchemspectemp-1 close (19) if (lroot) print*,'Number of compounds=',nchemspectemp if (nchemspectemp > nchemspec) call & stop_it("Too many chemicals! Change NCHEMSPEC in src/cparam.local") ! ! Allocate reaction arrays (but not during reloading!) ! if (.not. lreloading) then allocate(stoichio(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for stoichio") allocate(Sijm(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for Sijm") allocate(Sijp(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for Sijp") allocate(kreactions_z(mz,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_z") allocate(kreactions_p(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_p") allocate(kreactions_m(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_m") allocate(reaction_name(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for reaction_name") allocate(kreactions_profile(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_profile") allocate(kreactions_profile_width(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_profile_width") allocate(kreactions_alpha(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for kreactions_alpha") allocate(back(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for back") endif ! ! Initialize data ! kreactions_z = 1. Sijp = 0 Sijm = 0 back = .true. ! ! read chemistry data ! inquire (file=file1,exist=exist1) inquire (file=file2,exist=exist2) ! if (exist1 .and. exist2) then ! ! if both chemistry1.dat and chemistry2.dat are present, ! then read Sijp and Sijm, and calculate their sum ! ! file1 ! open (19,file=file1) read (19,*) chemicals do j = 1,mreactions read (19,*,end=994) kreactions_m(j),(Sijm(i,j),i=1,nchemspectemp) enddo 994 close (19) nreactions1 = j-1 ! ! file2 ! open (19,file=file2) read (19,*) chemicals do j = 1,mreactions if (lkreactions_profile) then if (lkreactions_alpha) then read (19,*) kreactions_p(j),(Sijp(i,j),i=1,nchemspectemp),kreactions_profile(j), & kreactions_profile_width(j),kreactions_alpha(j) else read (19,*) kreactions_p(j),(Sijp(i,j),i=1,nchemspectemp),kreactions_profile(j), & kreactions_profile_width(j) endif else if (lkreactions_alpha) then read (19,*) kreactions_p(j),(Sijp(i,j),i=1,nchemspectemp),kreactions_alpha(j) else read (19,*) kreactions_p(j),(Sijp(i,j),i=1,nchemspectemp) endif endif enddo close (19) nreactions2 = j-1 ! ! calculate stoichio and nreactions ! if (nreactions1 == nreactions2) then nreactions = nreactions1 stoichio = Sijp-Sijm else call stop_it('nreactions1/=nreactions2') endif if (nreactions /= mreactions) call stop_it('nreactions/=mreactions') ! else ! ! old method: read chemistry data, if present ! inquire (file='chemistry.dat',exist=exist) if (exist) then open (19,file='chemistry.dat') read (19,*) chemicals do j = 1,mreactions read (19,*,end=990) kreactions_p(j),(stoichio(i,j),i=1,nchemspec) enddo 990 close (19) nreactions = j-1 Sijm = -min(stoichio,0.0) Sijp = +max(stoichio,0.0) else if (lroot) print*,'no chemistry.dat file to be read.' lreactions = .false. endif endif ! ! print input data for verification ! if (lroot) then ! print*,'chemicals=',chemicals print*,'kreactions_m=',kreactions_m(1:nreactions) print*,'kreactions_p=',kreactions_p(1:nreactions) print*,'Sijm:' do i = 1,nreactions print*,Sijm(:,i) enddo print*,'Sijp:' do i = 1,nreactions print*,Sijp(:,i) enddo print*,'stoichio=' do i = 1,nreactions print*,stoichio(:,i) enddo endif ! ! possibility of z-dependent kreactions_z profile ! if (lkreactions_profile) then do j = 1,nreactions if (kreactions_profile(j) == 'cosh') then do n = 1,mz kreactions_z(n,j) = 1./cosh(z(n)/kreactions_profile_width(j))**2 enddo elseif (kreactions_profile(j) == 'gauss') then do n = 1,mz kreactions_z(n,j) = exp(-((z(n)/kreactions_profile_width(j))**2)) enddo elseif (kreactions_profile(j) == 'square') then do n = 1,mz if (n < mz/2) then kreactions_z(n,j) = kreactions_profile_width(j) else kreactions_z(n,j) = 0. endif enddo elseif (kreactions_profile(j) == 'saw') then do n = 1,mz kreactions_z(n,j) = 0.51+(sin(pi*z(n)/kreactions_profile_width(j)) & +sin(2*pi*z(n)/kreactions_profile_width(j))/2 & +sin(3*pi*z(n)/kreactions_profile_width(j))/3 & +sin(4*pi*z(n)/kreactions_profile_width(j))/4)/3 enddo elseif (kreactions_profile(j) == 'sin') then do n = 1,mz kreactions_z(n,j) = 0.5*(1+cos(pi*z(n)/kreactions_profile_width(j))) enddo elseif (kreactions_profile(j) == 'sin-bg') then do n = 1,mz kreactions_z(n,j) = 0.5*(1.1+cos(pi*z(n)/kreactions_profile_width(j))) enddo elseif (kreactions_profile(j) == 'spike') then do n = 1,mz if (cos(pi*z(n)/kreactions_profile_width(j)) > 0.99) then kreactions_z(n,j) = 1 else kreactions_z(n,j) = 0 endif enddo endif enddo endif ! call keep_compiler_quiet(f) ! endsubroutine astrobiology_data !*********************************************************************** subroutine dchemistry_dt(f,df,p) ! ! calculate right hand side of ONE OR MORE extra coupled PDEs ! along the 'current' Pencil, i.e. f(l1:l2,m,n) where ! m,n are global variables looped over in equ.f90 ! ! Due to the multi-step Runge Kutta timestepping used one MUST always ! add to the present contents of the df array. NEVER reset it to zero. ! ! several precalculated Pencils of information are passed if for ! efficiency. ! ! 13-aug-07/steveb: coded ! 8-jan-08/natalia: included advection/diffusion ! 20-feb-08/axel: included reactions ! 22-jun-10/julien: modified evaluation of enthalpy fluxes with ! constant Lewis numbers ! 10-jan-11/julien: modified to solve chemistry with LSODE ! use Diagnostics use Sub, only: grad,dot_mn use Special, only: special_calc_chemistry ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz,mvar) :: df real, dimension(nx,3) :: gchemspec, dk_D, sum_diff=0. real, dimension(nx) :: ugchemspec, sum_DYDT, sum_dhhk=0. real, dimension(nx) :: sum_dk_ghk, dk_dhhk, sum_hhk_DYDt_reac type (pencil_case) :: p real, dimension(nx) :: RHS_T_full, diffus_chem ! ! indices ! integer :: j, k,i,ii integer :: i1=1, i2=2, i3=3, i4=4, i5=5, i6=6, i7=7, i8=8, i9=9, i10=10 integer :: i11=11, i12=12, i13=13, i14=14, i15=15, i16=16, i17=17, i18=18, i19=19 integer :: iz1=1, iz2=2, iz3=3, iz4=4, iz5=5, iz6=6, iz7=7, iz8=8, iz9=9, iz10=10 integer :: iz11=11, iz12=12, iz13=13, iz14=14, iz15=15, iz16=16, iz17=17 integer :: iz18=18, iz19=19 ! intent(in) :: p,f intent(inout) :: df ! logical :: ldiffusion2 ldiffusion2 = ldiffusion .and. (.not. lchemonly) ! ! identify module and boundary conditions ! call timing('dchemistry_dt','entered',mnloop=.true.) if (headtt .or. ldebug) print*,'dchemistry_dt: SOLVE dchemistry_dt' ! ! Interface for your personal subroutines calls ! if (lspecial) call special_calc_chemistry(f,df,p) ! ! loop over all chemicals ! do k = 1,nchemspec ! ! advection terms ! if (lhydro .and. ladvection .and.(.not. lchemonly)) then call grad(f,ichemspec(k),gchemspec) call dot_mn(p%uu,gchemspec,ugchemspec) if (lmobility) ugchemspec = ugchemspec*mobility(k) df(l1:l2,m,n,ichemspec(k)) = df(l1:l2,m,n,ichemspec(k))-ugchemspec endif ! ! diffusion operator ! ! Temporary we check the existence of chem.imp data, ! further one should check the existence of a file with ! binary diffusion coefficients! ! if (ldiffusion2) then df(l1:l2,m,n,ichemspec(k)) = df(l1:l2,m,n,ichemspec(k))+ & p%DYDt_diff(:,k) endif ! ! chemical reactions: ! multiply with stoichiometric matrix with reaction speed ! d/dt(x_i) = S_ij v_j ! if (lreactions) then if (lchemonly) then ! If chemistry is solved in a separate step, we want df to contain only the ! chemical contribution, that's why no sum is required df(l1:l2,m,n,ichemspec(k)) = p%DYDt_reac(:,k) else df(l1:l2,m,n,ichemspec(k)) = df(l1:l2,m,n,ichemspec(k))+ & p%DYDt_reac(:,k) endif endif ! ! Add filter for negative concentrations ! if (lfilter .and. .not. lfilter_strict) then do i = 1,mx if ((f(i,m,n,ichemspec(k))+df(i,m,n,ichemspec(k))*dt) < -1e-25 ) then df(i,m,n,ichemspec(k)) = -1e-25*dt endif if ((f(i,m,n,ichemspec(k))+df(i,m,n,ichemspec(k))*dt) > 1. ) then df(i,m,n,ichemspec(k)) = 1.*dt endif enddo endif ! ! Add strict filter for negative concentrations ! if (lfilter_strict) then do i = 1,mx if ((f(i,m,n,ichemspec(k))+df(i,m,n,ichemspec(k))*dt) < 0.0 ) then if (df(i,m,n,ichemspec(k)) < 0.) df(i,m,n,ichemspec(k)) = 0. endif if ((f(i,m,n,ichemspec(k))+df(i,m,n,ichemspec(k))*dt) > 1. ) then df(i,m,n,ichemspec(k)) = 1.*dt endif enddo endif ! enddo ! if (ldensity .and. lcheminp) then ! if (l1step_test) then sum_DYDt = 0. do i = 1,nx sum_DYDt(i) = -p%rho(1)*(p%TT(i)-Tinf)*p%TT1(i) & *Cp_const/lambda_const*beta*(beta-1.)*f(l1,m,n,iux)*f(l1,m,n,iux) enddo else sum_DYDt = 0. sum_hhk_DYDt_reac = 0. sum_dk_ghk = 0. ! do k = 1,nchemspec if (species_constants(k,imass) > 0.) then sum_DYDt = sum_DYDt+Rgas/species_constants(k,imass) & *(p%DYDt_reac(:,k)+p%DYDt_diff(:,k)) if (lreactions) then sum_hhk_DYDt_reac = sum_hhk_DYDt_reac-p%hhk_full(:,k)*p%DYDt_reac(:,k) endif ! ! Sum over all species of diffusion terms ! if (ldiffusion2) then if (lDiff_fick) then call grad(f,ichemspec(k),gchemspec) do i = 1,3 dk_D(:,i) = gchemspec(:,i)*p%Diff_penc_add(:,k) enddo elseif (lFlux_simple) then do i = 1,3 dk_D(:,i) = p%gXXk(:,i,k)*p%Diff_penc_add(:,k)*p%mukmu1(:,k) enddo else do i = 1,3 dk_D(:,i) = (p%gXXk(:,i,k) & +(XX_full(l1:l2,m,n,k)-f(l1:l2,m,n,ichemspec(k)))*p%glnpp(:,i)) & *p%Diff_penc_add(:,k)*p%mukmu1(:,k) enddo endif ! call dot_mn(dk_D,p%ghhk(:,:,k),dk_dhhk) sum_dk_ghk = sum_dk_ghk+dk_dhhk if (ldiff_corr) sum_diff(:,k) = sum_diff(:,k)+dk_D(:,k) endif ! endif enddo ! ! If the correction velocity is added ! if (ldiff_corr .and. ldiffusion2) then call fatal_error('dchemistry_dt',& 'The correction vel. is not properly implemented - pleas fix!') do k = 1,nchemspec call dot_mn(sum_diff,p%ghhk(:,:,k),sum_dhhk) sum_dk_ghk(:) = sum_dk_ghk(:)-f(l1:l2,m,n,ichemspec(k))*sum_dhhk(:) enddo endif endif ! if (l1step_test) then RHS_T_full = sum_DYDt(:) else if (ltemperature_nolog) then RHS_T_full = p%cv1*((sum_DYDt(:)-Rgas*p%mu1*p%divu)*p%TT & +sum_dk_ghk+sum_hhk_DYDt_reac) ! call stop_it('ltemperature_nolog case does not work now!') else if (lchemonly) then RHS_T_full = (sum_DYDt(:)+sum_hhk_DYDt_reac*p%TT1(:))*p%cv1 else RHS_T_full = (sum_DYDt(:)-Rgas*p%mu1*p%divu)*p%cv1 & +sum_dk_ghk*p%TT1(:)*p%cv1+sum_hhk_DYDt_reac*p%TT1(:)*p%cv1 endif endif endif ! if (.not. lT_const) then if (lchemonly) then ! If chemistry is solved in a separate step, we want df to contain only the ! chemical contribution, that's why no sum is required df(l1:l2,m,n,ilnTT) = RHS_T_full else df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + RHS_T_full endif endif ! if (lheatc_chemistry .and.(.not. lchemonly)) & call calc_heatcond_chemistry(f,df,p) endif ! if (lreactions .and. ireac /= 0 .and. ((.not. llsode).or. lchemonly)) then if (llast) call get_reac_rate(sum_hhk_DYDt_reac,f,p) endif ! ! Atmosphere case ! if (lcloud .and.(.not. lchemonly)) then ! df(l1:l2,m,n,ichemspec(index_H2O)) = df(l1:l2,m,n,ichemspec(index_H2O)) & - p%ccondens do i = 1,mx if ((f(i,m,n,ichemspec(index_H2O)) & +df(i,m,n,ichemspec(index_H2O))*dt) >= 1. ) & df(i,m,n,ichemspec(index_H2O)) = 0. if ((f(i,m,n,ichemspec(index_H2O)) & +df(i,m,n,ichemspec(index_H2O))*dt) < 0. ) & df(i,m,n,ichemspec(index_H2O)) = 0. enddo endif ! ! this damping zone is needed in a case of NSCBC ! if (ldamp_zone_for_NSCBC) call damp_zone_for_NSCBC(f,df) ! ! For the timestep calculation, need maximum diffusion ! if (lfirst .and. ldt .and. (.not. lchemonly)) then if (.not. lcheminp) then diffus_chem = chem_diff*maxval(chem_diff_prefactor)*dxyz_2 else diffus_chem=0. do j = 1,nx if (ldiffusion .and. .not. ldiff_simple) then ! !-------------------------------------- ! This expression should be discussed !-------------------------------------- ! diffus_chem(j) = diffus_chem(j)+ & maxval(Diff_full_add(l1+j-1,m,n,1:nchemspec))*dxyz_2(j) else diffus_chem(j) = 0. endif enddo endif maxdiffus=max(maxdiffus,diffus_chem) endif ! ! NB: it should be discussed ! if (lfirst .and. ldt) then if (lreactions .and.(.not. llsode .or. lchemonly)) then ! ! calculate maximum of *relative* reaction rate if decaying, ! or maximum of absolute rate, if growing. ! if (lchem_cdtc) then reac_chem = 0. do k = 1,nchemspec reac_chem = max(reac_chem, & abs(p%DYDt_reac(:,k)/max(f(l1:l2,m,n,ichemspec(k)),.001))) enddo ! elseif (lcheminp) then reac_chem = 0. !sum_reac_rate=0. do k = 1,nchemspec reac_chem = reac_chem+abs(p%DYDt_reac(:,k)/ & max(f(l1:l2,m,n,ichemspec(k)),0.001)) !sum_reac_rate=sum_reac_rate+p%DYDt_reac(:,k) enddo if (maxval(reac_chem) > 1e11) then reac_chem = 1e11 endif endif endif endif ! ! Calculate diagnostic quantities ! call timing('dchemistry_dt','before ldiagnos',mnloop=.true.) if (ldiagnos) then if (idiag_dtchem /= 0) then call max_mn_name(reac_chem/cdtc,idiag_dtchem,l_dt=.true.) endif ! ! WL: instead of hardcoding Y1-Y9, wouldn't it be possible ! to have them all in the same array? The nbody ! module, for instance, has idiag_xxq and idiag_vvq, which ! allows the user to add output the positions and velocities ! of as many particle they wants. ! RP: Totally agree... I still have to expand manually these hard-coded ! Y1-Y9 and chemspec-chemspec9 when needed, but this is just work-around... ! AB: I also agree! ! ! do ii = 1,nchemspec call sum_mn_name(f(l1:l2,m,n,ichemspec(ii)),idiag_Ym(ii)) call max_mn_name(f(l1:l2,m,n,ichemspec(ii)),idiag_Ymax(ii)) if (idiag_Ymin(ii)/= 0) & call max_mn_name(-f(l1:l2,m,n,ichemspec(ii)),idiag_Ymin(ii),lneg=.true.) if (idiag_TYm(ii)/= 0) & call sum_mn_name(max(1.-f(l1:l2,m,n,ichemspec(ii))/Ythresh(ii),0.),idiag_TYm(ii)) if (idiag_diffm(ii)/= 0) & call sum_mn_name(Diff_full_add(l1:l2,m,n,ii),idiag_diffm(ii)) enddo ! call sum_mn_name(cp_full(l1:l2,m,n),idiag_cpfull) call sum_mn_name(cv_full(l1:l2,m,n),idiag_cvfull) ! call sum_mn_name(lambda_full(l1:l2,m,n),idiag_lambdam) call sum_mn_name(f(l1:l2,m,n,iviscosity), idiag_num) ! ! Sample for hard coded diffusion diagnostics ! ! if (idiag_diff1m /= 0) call sum_mn_name(diff_full(l1:l2,m,n,i1), & ! idiag_diff1m) endif ! ! 1d-averages. Happens at every it1d timesteps, NOT at every it1 ! if (l1davgfirst) then do ii = 1,nchemspec call xysum_mn_name_z(f(l1:l2,m,n,ichemspec(ii)),idiag_Ymz(ii)) enddo endif call timing('dchemistry_dt','finished',mnloop=.true.) ! endsubroutine dchemistry_dt !*********************************************************************** subroutine read_chemistry_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read (parallel_unit, NML=chemistry_init_pars, IOSTAT=iostat) ! endsubroutine read_chemistry_init_pars !*********************************************************************** subroutine write_chemistry_init_pars(unit) ! integer, intent(in) :: unit ! write (unit, NML=chemistry_init_pars) ! endsubroutine write_chemistry_init_pars !*********************************************************************** subroutine read_chemistry_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read (parallel_unit, NML=chemistry_run_pars, IOSTAT=iostat) ! endsubroutine read_chemistry_run_pars !*********************************************************************** subroutine write_chemistry_run_pars(unit) ! integer, intent(in) :: unit ! write (unit, NML=chemistry_run_pars) ! endsubroutine write_chemistry_run_pars !*********************************************************************** subroutine rprint_chemistry(lreset,lwrite) ! ! reads and registers print parameters relevant to chemistry ! ! 13-aug-07/steveb: coded ! use Diagnostics, only: parse_name use FArrayManager, only: farray_index_append use General, only: itoa, get_species_nr ! integer :: iname, inamez,ii logical :: lreset, lwr logical, optional :: lwrite character(len=6) :: diagn_Ym, number character(len=6) :: diagn_Ymax character(len=6) :: diagn_Ymin character(len=7) :: diagn_dYmax character(len=6) :: diagn_TYm character(len=6) :: diagn_dYm character(len=6) :: diagn_hm character(len=6) :: diagn_cpm character(len=7) :: diagn_diffm character(len=fmtlen) :: sname ! lwr = .false. if (present(lwrite)) lwr = lwrite ! ! reset everything in case of reset ! (this needs to be consistent with what is defined above!) ! if (lreset) then idiag_dtchem = 0 idiag_Ym = 0 idiag_TYm = 0 idiag_dYm = 0 idiag_Ymax = 0 idiag_Ymin = 0 idiag_dYmax = 0 idiag_hm = 0 idiag_cpm = 0 idiag_diffm = 0 ! idiag_cpfull = 0 idiag_cvfull = 0 idiag_e_intm = 0 idiag_Ymz = 0 ! idiag_lambdam = 0 idiag_num = 0 ! endif ! ! check for those quantities that we want to evaluate online ! do iname = 1,nname do ii=1,nchemspec write (number,'(I2)') ii diagn_Ym = 'Y'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_Ym),idiag_Ym(ii)) diagn_Ymax = 'Y'//trim(adjustl(number))//'max' call parse_name(iname,cname(iname),cform(iname),trim(diagn_Ymax),idiag_Ymax(ii)) diagn_Ymin = 'Y'//trim(adjustl(number))//'min' call parse_name(iname,cname(iname),cform(iname),trim(diagn_Ymin),idiag_Ymin(ii)) diagn_dYm = 'dY'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_dYm),idiag_dYm(ii)) diagn_TYm = 'TY'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_TYm),idiag_TYm(ii)) diagn_dYmax = 'dY'//trim(adjustl(number))//'max' call parse_name(iname,cname(iname),cform(iname),trim(diagn_dYmax),idiag_dYmax(ii)) diagn_hm = 'h'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_hm),idiag_hm(ii)) diagn_cpm = 'cp'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_cpm),idiag_cpm(ii)) diagn_diffm = 'diff'//trim(adjustl(number))//'m' call parse_name(iname,cname(iname),cform(iname),trim(diagn_diffm),idiag_diffm(ii)) enddo call parse_name(iname,cname(iname),cform(iname),'dtchem',idiag_dtchem) call parse_name(iname,cname(iname),cform(iname),'cpfull',idiag_cpfull) call parse_name(iname,cname(iname),cform(iname),'cvfull',idiag_cvfull) ! ! Sample for hard-coded heat capacity diagnostics ! ! call parse_name(iname,cname(iname),cform(iname),'cp1m',idiag_cp1m) call parse_name(iname,cname(iname),cform(iname),'e_intm',idiag_e_intm) call parse_name(iname,cname(iname),cform(iname),'lambdam',idiag_lambdam) call parse_name(iname,cname(iname),cform(iname),'num',idiag_num) ! ! Sample for hard-coded diffusion diagnostics ! ! call parse_name(iname,cname(iname),cform(iname),'diff1m',idiag_diff1m) enddo ! ! xy-averages ! do inamez = 1,nnamez do ii=1,nchemspec write (number,'(I2)') ii diagn_Ym = 'Y'//trim(adjustl(number))//'mz' call parse_name(inamez,cnamez(inamez),cformz(inamez),trim(diagn_Ym),idiag_Ymz(ii)) enddo enddo ! ! check for those quantities for which we want video slices ! do iname=1,nnamev sname=trim(cnamev(iname)) if (sname(1:8)=='chemspec') then if (get_species_nr(sname,'chemspec',nchemspec,'rprint_chemistry')>0) & cformv(iname)='DEFINED' endif enddo ! endsubroutine rprint_chemistry !*********************************************************************** subroutine get_slices_chemistry(f,slices) ! ! Write slices for animation of Chemistry variables. ! ! 13-aug-07/steveb: dummy ! 16-may-09/raphael: added more slices ! use Slices_methods, only: assign_slices_scal real, dimension(mx,my,mz,mfarray) :: f type (slice_data) :: slices character(len=fmtlen) :: sname integer :: ispec ! ! Chemical species mass fractions. ! sname=trim(slices%name) if (sname(9:)==' ') then ! 9=len('chemspec')+1 ispec=1 else read(sname(9:),'(i3)') ispec endif ! call assign_slices_scal(slices,f,ichemspec(ispec)) ! endsubroutine get_slices_chemistry !*********************************************************************** subroutine build_stoich_matrix(StartInd,StopInd,k,ChemInpLine,product) ! ! calculation of the stoichoimetric matrix ! ! 10-mar-08/nils: coded ! integer, intent(in) :: StartInd, StopInd,k character(len=*), intent(in) :: ChemInpLine logical, intent(in) :: product integer :: StartSpecie, ind_glob, ind_chem real :: stoi logical :: found_specie, lreal=.false. integer :: stoi_int ! if ((ChemInpLine(StartInd:StopInd) /= "M" ) & .and. (ChemInpLine(StartInd:StartInd+1) /= "hv" )) then StartSpecie = verify(ChemInpLine(StartInd:StopInd), & "1234567890")+StartInd-1 ! ! Call to a routine that checks for arbitrary stoiciometric coefficents ! removes them and shifts the begin of the species index in cheminpline ! if (StopInd-StartSpecie >= 3 .and. StartSpecie > 1) & call find_remove_real_stoic(ChemInpLine(StartSpecie-1:StopInd),lreal,stoi,StartSpecie) ! call find_species_index(ChemInpLine(StartSpecie:StopInd), & ind_glob,ind_chem,found_specie) if (.not. found_specie) then print*,'ChemInpLine=',ChemInpLine print*,'StartSpecie,StopInd=',StartSpecie,StopInd print*,'ChemInpLine(StartSpecie:StopInd)=',ChemInpLine(StartSpecie:StopInd) print*,'ind_glob,ind_chem=',ind_glob,ind_chem ! if (.not. lpencil_check_small) then ! if (.not. lpencil_check) then call stop_it("build_stoich_matrix: Did not find species!") ! endif ! endif endif ! ! if (found_specie) then if (StartSpecie == StartInd) then stoi = 1.0 else if (.not. lreal) then read (unit=ChemInpLine(StartInd:StartInd),fmt='(I1)') stoi_int endif endif endif ! if (product) then Sijm(ind_chem,k) = Sijm(ind_chem,k)+stoi else Sijp(ind_chem,k) = Sijp(ind_chem,k)+stoi endif ! endif ! endsubroutine build_stoich_matrix !*********************************************************************** subroutine write_reactions ! ! write reaction coefficient in the output file ! ! 11-mar-08/nils: coded ! use General, only: itoa ! integer :: reac, spec character(len=80) :: reac_string, product_string, output_string character(len=intlen) :: Sijp_string, Sijm_string character(len=1) :: separatorp, separatorm character(len=fnlen) :: input_file="./data/chem.out" integer :: file_id=123 ! open (file_id,file=input_file,POSITION='APPEND',FORM='FORMATTED') write (file_id,*) 'REACTIONS' !open(file_id,file=input_file) ! do reac = 1,mreactions reac_string = '' product_string = '' separatorp = '' separatorm = '' ! do spec = 1,nchemspec if (Sijp(spec,reac) > 0) then Sijp_string = '' if (Sijp(spec,reac) /= 1) then write (Sijp_string,'(F3.1)') Sijp(spec,reac) endif ! if (Sijp(spec,reac)>1) Sijp_string=itoa(Sijp(spec,reac)) reac_string = trim(reac_string)//trim(separatorp)// & trim(Sijp_string)//trim(varname(ichemspec(spec))) separatorp = '+' endif if (Sijm(spec,reac) > 0) then Sijm_string = '' if (Sijm(spec,reac) /= 1) write (Sijm_string,'(F3.1)') Sijm(spec,reac) ! if (Sijm(spec,reac)>1) Sijm_string=itoa(Sijm(spec,reac)) product_string = trim(product_string)//trim(separatorm)// & trim(Sijm_string)//trim(varname(ichemspec(spec))) separatorm = '+' endif enddo ! output_string = trim(reac_string)//'='//trim(product_string) ! if (.not. photochem_case(reac)) then ! ! Note that since the B_n term is in logarithmic form within the code ! the exponential must be used for output. ! write (unit=output_string(30:45),fmt='(E14.4)') exp(B_n(reac)) write (unit=output_string(47:62),fmt='(E14.4)') alpha_n(reac) write (unit=output_string(64:79),fmt='(E14.4)') E_an(reac) endif write (file_id,*) trim(output_string) if (.not. photochem_case(reac)) then if (maxval(abs(low_coeff(:,reac))) > 0.) then write (file_id,*) 'LOW/',exp(low_coeff(1,reac)),low_coeff(2:,reac) elseif (maxval(abs(high_coeff(:,reac))) > 0.) then write (file_id,*) 'HIGH/',high_coeff(:,reac) endif if (maxval(abs(troe_coeff(:,reac))) > 0.) then write (file_id,*) 'TROE/',troe_coeff(:,reac) endif if (minval(a_k4(:,reac)) < impossible) then write (file_id,*) a_k4(:,reac) endif else write (file_id,*) ' min lambda=',lamb_low,' max lambda=',lamb_up endif ! enddo ! write (file_id,*) 'END' write (file_id,*) '(M+) case: ',Mplus_case write (file_id,*) 'photochemical case: ',photochem_case ! close (file_id) ! endsubroutine write_reactions !*********************************************************************** subroutine chemkin_data(f) ! ! if the file with chemkin data exists ! reading the Chemkin data ! ! 21-jul-10/julien: Reading lewis.dat file to collect constant Lewis numbers ! for each species ! character(len=fnlen) :: input_file real, dimension(mx,my,mz,mfarray) :: f integer :: stat, k,i character(len=fnlen) :: input_file2="./data/stoich.out" integer :: file_id=123 logical :: chemin,cheminp ! inquire(file='chem.inp',exist=cheminp) inquire(file='chem.in',exist=chemin) if (chemin .and. cheminp) call fatal_error('chemistry',& 'chem.in and chem.inp found, please decide for one') if (cheminp) input_file='chem.inp' if (chemin) input_file='chem.in' ! inquire (file='tran.dat',exist=tran_exist) if (.not. tran_exist) then inquire (file='tran.in',exist=tran_exist) endif inquire (file='lewis.dat',exist=lew_exist) ! ! Allocate binary diffusion coefficient array ! ! if (.not. lreloading) then !Natalia: !this does not work for ldiffusion=F ! if (ldiffusion .and. .not. lfix_Sc) then if (.not. lfix_Sc) then !NILS: Since Bin_diff_coeff is such a huge array we must check if it !NILS: required to define it for the full domain!!!!!! if (.not. lreloading) then allocate(Bin_Diff_coef(mx,my,mz,nchemspec,nchemspec),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory "// & "for binary diffusion coefficients") ! allocate(Diff_full(mx,my,mz,nchemspec),STAT=stat) allocate(Diff_full_add(mx,my,mz,nchemspec),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory "// & "for binary diffusion coefficients") endif ! endif ! endif ! if (tran_exist) then if (lroot) then print*,'tran.in/dat file with transport data is found.' endif call read_transport_data endif ! if (lew_exist) then if (lroot) then print*,'lewis.dat file with transport data is found.' print*,'Species diffusion coefficients calculated using constant Lewis numbers.' endif call read_Lewis endif ! if (.not. lew_exist .and. lDiff_lewis .and. lroot) then if (.not. l1step_test) then print*, 'No lewis.dat file present, switch to simplified diffusion' lDiff_lewis = .false. lDiff_simple = .true. else print*, 'Le=1' lDiff_lewis = .true. endif endif ! if (lroot .and. .not. tran_exist .and. .not. lew_exist) then if (chem_diff == 0.) & call inevitably_fatal_error('chemkin data', 'chem_diff = 0.') print*,'tran.dat file with transport data is not found.' print*,'lewis.dat file with Lewis numbers is not found.' print*,'Now diffusion coefficients is ',chem_diff print*,'Now species viscosity is ',nu_spec Bin_Diff_coef = chem_diff/(unit_length*unit_length/unit_time) do k = 1,nchemspec species_viscosity(:,:,:,k) = nu_spec(k)/ & (unit_mass/unit_length/unit_time) enddo endif ! ! Find number of ractions ! call read_reactions(input_file,NrOfReactions=mreactions) if (lroot) print*,'Number of reactions=',mreactions if (lroot) print*,'Number of species=',nchemspec nreactions = mreactions ! ! Allocate reaction arrays ! if (.not. lreloading) then allocate(stoichio(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for stoichio") allocate(Sijm(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for Sijm") allocate(Sijp(nchemspec,mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for Sijp") allocate(reaction_name(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for reaction_name") allocate(back(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for back") allocate(B_n(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for B_n") B_n = 0. allocate(alpha_n(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for alpha_n") alpha_n = 0. allocate(E_an(mreactions),STAT=stat) if (stat > 0) call stop_it("Couldn't allocate memory for E_an") E_an = 0. ! allocate(low_coeff(3,nreactions),STAT=stat) low_coeff = 0. if (stat > 0) call stop_it("Couldn't allocate memory for low_coeff") allocate(high_coeff(3,nreactions),STAT=stat) high_coeff = 0. if (stat > 0) call stop_it("Couldn't allocate memory for high_coeff") allocate(troe_coeff(3,nreactions),STAT=stat) troe_coeff = 0. if (stat > 0) call stop_it("Couldn't allocate memory for troe_coeff") allocate(a_k4(nchemspec,nreactions),STAT=stat) a_k4 = impossible if (stat > 0) call stop_it("Couldn't allocate memory for troe_coeff") allocate(Mplus_case (nreactions),STAT=stat) Mplus_case = .false. allocate(photochem_case (nreactions),STAT=stat) photochem_case = .false. if (stat > 0) call stop_it("Couldn't allocate memory for photochem_case") endif ! ! Initialize data ! Sijp = 0. Sijm = 0.0 back = .true. ! ! read chemistry data ! call read_reactions(input_file) call write_reactions ! ! calculate stoichio and nreactions ! stoichio = Sijp-Sijm ! ! print input data for verification ! if (lroot .and. nreactions > 0) then ! open (file_id,file=input_file2,POSITION='rewind',FORM='FORMATTED') write (file_id,*) 'STOICHIOMETRIC MATRIX' ! write (file_id,*) 'Species names' write (file_id,101) varname(ichemspec(:)) ! write (file_id,*) 'Sijm' do i = 1,nreactions write (file_id,100) i,Sijm(:,i) enddo write (file_id,*) 'Sijp:' do i = 1,nreactions write (file_id,100) i,Sijp(:,i) enddo write (file_id,*) 'stoichio=' do i = 1,nreactions write (file_id,100) i,stoichio(:,i) enddo close (file_id) endif ! 100 format(I1,26f6.1) 101 format(' ',26A6) ! call keep_compiler_quiet(f) ! endsubroutine chemkin_data !*********************************************************************** subroutine read_reactions(input_file,NrOfReactions) ! ! This subroutine reads all reaction information from chem.inp ! See the chemkin manual for more information on ! the syntax of chem.inp. ! ! 10-mar-08/nils: coded ! logical :: IsReaction=.false., found_new_reaction=.false. logical, save :: find_specie, found_specie integer, optional :: NrOfReactions integer, save :: ind_glob, ind_chem integer :: i, k, file_id=123, StartInd, StopInd, StartInd_add integer :: StopInd_add, StopInd_add_, StopIndName integer :: VarNumber, VarNumber_add, SeparatorInd integer :: PlusInd integer :: LastLeftCharacter, ParanthesisInd, Mplussind integer :: photochemInd, plusind_ character(len=120) :: ChemInpLine, ChemInpLine_add character(len=*) :: input_file if (present(NrOfReactions)) NrOfReactions=0 k=1 open(file_id,file=input_file) dataloop3: do if (found_new_reaction) then ChemInpLine=ChemInpLine_add else read(file_id,'(80A)',end=1012) ChemInpLine endif found_new_reaction=.false. ! ! Check if we are reading a line within the reactions section ! if (ChemInpLine(1:9)=="REACTIONS") IsReaction=.true. if (ChemInpLine(1:3)=="END" .and. IsReaction) IsReaction=.false. if (present(NrOfReactions)) then ! ! Find number of reactions ! if (IsReaction) then if (ChemInpLine(1:9) /= "REACTIONS") then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),'=')+StartInd-1 if (StopInd>0 .and. ChemInpLine(1:1) /= '!') then NrOfReactions=NrOfReactions+1 endif endif endif else ! ! Read in species ! if (IsReaction) then if (ChemInpLine(1:9) /= "REACTIONS") then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),'=')+StartInd-1 if (StopInd>0 .and. ChemInpLine(1:1) /= '!') then ! ! Fill in reaction name ! StopIndName=index(ChemInpLine(StartInd:),' ')+StartInd-1 reaction_name(k)=ChemInpLine(StartInd:StopIndName) ! ! Photochemical case ! ParanthesisInd=0 photochemInd=0 SeparatorInd=0 photochemInd=index(ChemInpLine(StartInd:),'hv') if (photochemInd>0) then photochem_case (k)=.true. ParanthesisInd=index(ChemInpLine(photochemInd:),'(') & +photochemInd-1 if ((ParanthesisInd>0) .and. (photochemInd>0)) then StopInd=index(ChemInpLine(StartInd:),'lam') SeparatorInd=index(ChemInpLine(ParanthesisInd:StopInd),'<') if (SeparatorInd>0) then SeparatorInd=SeparatorInd+ParanthesisInd-1 read (unit=ChemInpLine(ParanthesisInd+1:& SeparatorInd-1),fmt='(E15.8)') lamb_low endif ParanthesisInd=index(ChemInpLine(StopInd:),')') +StopInd-1 SeparatorInd=index(ChemInpLine(StopInd:ParanthesisInd),'<') & +StopInd-1 if (SeparatorInd>0) then read (unit=ChemInpLine(SeparatorInd+1:& ParanthesisInd-1),fmt='(E15.8)') lamb_up endif endif endif ! ! End of the photochemical case ! ! Find reactant side stoichiometric coefficients ! SeparatorInd=index(ChemInpLine(StartInd:),'<=') if (SeparatorInd==0) then SeparatorInd=index(ChemInpLine(StartInd:),'=') if (index(ChemInpLine(StartInd:),'=>')/=0) back(k)=.false. endif ParanthesisInd=0 MplussInd=0 ParanthesisInd=index(ChemInpLine(StartInd:),'(+M)') MplussInd=index(ChemInpLine(StartInd:),'+M') found_new_reaction=.false. if (ParanthesisInd>0 .or. Mplussind>0) then if (ParanthesisInd>0) then LastLeftCharacter=min(ParanthesisInd,SeparatorInd)-1 else LastLeftCharacter=min(MplussInd,SeparatorInd)-1 endif ! ! reading of the additional data for (+M) case ! 100 read(file_id,'(80A)',end=1012) ChemInpLine_add if (ChemInpLine_add(1:1) == ' ') then if (ParanthesisInd>0) then Mplus_case (k)=.true. endif i=1 do while (i<80) if (ChemInpLine_add(i:i)==' ') then i=i+1 elseif (ChemInpLine_add(i:i+2)=='LOW') then ! if (lroot) print*,ChemInpLine_add(i:i+2),& ! ' coefficients for reaction ', & ! reaction_name(k),'number ', k VarNumber_add=1; StartInd_add=i+4; StopInd_add=i+4 do while (VarNumber_add<4) StopInd_add=index(ChemInpLine_add(StartInd_add:),& ' ')+StartInd_add-2 StopInd_add_=index(ChemInpLine_add(StartInd_add:),& '/')+StartInd_add-2 StopInd_add=min(StopInd_add,StopInd_add_) if (StopInd_add==StartInd_add) then StartInd_add=StartInd_add+1 else if (VarNumber_add==1) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') low_coeff(1,k) if (low_coeff(1,k)/=0.) low_coeff(1,k)=log(low_coeff(1,k)) elseif (VarNumber_add==2) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') low_coeff(2,k) elseif (VarNumber_add==3) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') low_coeff(3,k) else call stop_it("No such VarNumber!") endif endif VarNumber_add=VarNumber_add+1 !StartInd_add=StopInd_add StartInd_add=verify(ChemInpLine_add(StopInd_add+1:),& ' ')+StopInd_add StopInd_add=StartInd_add enddo i=80 elseif (ChemInpLine_add(i:i+3)=='TROE') then ! if (lroot) print*,ChemInpLine_add(i:i+3),& ! ' coefficients for reaction ', reaction_name(k),& ! 'number ', k VarNumber_add=1; StartInd_add=i+5; StopInd_add=i+5 do while (VarNumber_add<4) StopInd_add=index(ChemInpLine_add(StartInd_add:),& ' ')+StartInd_add-2 StopInd_add_=index(ChemInpLine_add(StartInd_add:),& '/')+StartInd_add-2 StopInd_add=min(StopInd_add,StopInd_add_) if (StopInd_add==StartInd_add) then StartInd_add=StartInd_add+1 else if (VarNumber_add==1) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') troe_coeff(1,k) elseif (VarNumber_add==2) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') troe_coeff(2,k) elseif (VarNumber_add==3) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') troe_coeff(3,k) else call stop_it("No such VarNumber!") endif endif VarNumber_add=VarNumber_add+1 ! StartInd_add=StopInd_add StartInd_add=verify(ChemInpLine_add(StopInd_add+1:),& ' ')+StopInd_add StopInd_add=StartInd_add enddo i=80 elseif (ChemInpLine_add(i:i+3)=='HIGH') then ! if (lroot) print*,ChemInpLine_add(i:i+3),& ! ' coefficients for reaction ', reaction_name(k),& ! 'number ', k VarNumber_add=1; StartInd_add=i+5; StopInd_add=i+5 do while (VarNumber_add<4) StopInd_add=index(ChemInpLine_add(StartInd_add:),& ' ')+StartInd_add-2 StopInd_add_=index(ChemInpLine_add(StartInd_add:),& '/')+StartInd_add-2 StopInd_add=min(StopInd_add,StopInd_add_) if (StopInd_add==StartInd_add) then StartInd_add=StartInd_add+1 else if (VarNumber_add==1) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') high_coeff(1,k) if (high_coeff(1,k)/=0.) high_coeff(1,k)=log(high_coeff(1,k)) elseif (VarNumber_add==2) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') high_coeff(2,k) elseif (VarNumber_add==3) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') high_coeff(3,k) else call stop_it("No such VarNumber!") endif endif VarNumber_add=VarNumber_add+1 ! StartInd_add=StopInd_add StartInd_add=verify(ChemInpLine_add(StopInd_add+1:),& ' ')+StopInd_add StopInd_add=StartInd_add enddo i=80 else ! a_k4=0. StartInd_add=i; StopInd_add=0; StopInd_add_=0 do while (ChemInpLine_add(i:i+1)/=' ') find_specie=.true. do while (StartInd_add/=StopInd_add_) StopInd_add=index(ChemInpLine_add(StartInd_add:),& '/')+StartInd_add-2 StopInd_add_=index(ChemInpLine_add(StartInd_add:),& ' ')+StartInd_add-1 if (find_specie) then call find_species_index(ChemInpLine_add& (StartInd_add:StopInd_add),ind_glob,& ind_chem,found_specie) else if (found_specie) then read (unit=ChemInpLine_add(StartInd_add:& StopInd_add),fmt='(E15.8)') a_k4(ind_chem,k) else print*,'ChemInpLine=',ChemInpLine_add print*,'Specie=',& ChemInpLine_add(StartInd_add:StopInd_add) print*,'StartInd_add,StopInd_add=',& StartInd_add,StopInd_add call stop_it("read_reactions: Did not find specie!") endif endif StartInd_add=StopInd_add+2 find_specie=.false. enddo i=StopInd_add_ StartInd_add=StartInd_add+1 enddo i=80 !call find_species_index('N2',& ! ind_glob,ind_chem,found_specie) !if (found_specie) a_k4(ind_chem,k)=1. do ind_chem=1,nchemspec if (a_k4(ind_chem,k)==impossible) a_k4(ind_chem,k)=1 enddo endif enddo goto 100 ! this is an excellent example of "spaghetti code" [Bourdin.KIS] else found_new_reaction=.true. endif else LastLeftCharacter=SeparatorInd-1 endif StartInd=1 PlusInd=index(ChemInpLine(StartInd:LastLeftCharacter),'+')& +StartInd-1 do while (PlusInd0) StopInd=PlusInd-1 call build_stoich_matrix(StartInd,StopInd,k,& ChemInpLine,.false.) StartInd=StopInd+2 plusind_=index(ChemInpLine(StartInd:),'+') if (plusind_ > 0) then PlusInd=plusind_+StartInd-1 else PlusInd=10000 endif enddo StopInd=LastLeftCharacter call build_stoich_matrix(StartInd,StopInd,k,ChemInpLine,.false.) ! ! Find product side stoichiometric coefficients ! StartInd=index(ChemInpLine,'>')+1 if (StartInd==1) StartInd=index(ChemInpLine,'=')+1 SeparatorInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 ParanthesisInd=index(ChemInpLine(StartInd:),'(+M)')+StartInd-1 MplussInd=index(ChemInpLine(StartInd:),'+M')+StartInd-1 if (ParanthesisInd>StartInd) then LastLeftCharacter=min(ParanthesisInd,SeparatorInd)-1 elseif (MplussInd>StartInd) then LastLeftCharacter=min(MplussInd,SeparatorInd)-1 else LastLeftCharacter=SeparatorInd-1 endif PlusInd=index(ChemInpLine(StartInd:LastLeftCharacter),'+')& +StartInd-1 do while (PlusIndStartInd) StopInd=PlusInd-1 call build_stoich_matrix(StartInd,StopInd,k,& ChemInpLine,.true.) StartInd=StopInd+2 PlusInd=index(ChemInpLine(StartInd:),'+')+StartInd-1 enddo StopInd=LastLeftCharacter call build_stoich_matrix(StartInd,StopInd,k,ChemInpLine,.true.) ! ! Find Arrhenius coefficients ! VarNumber=1; StartInd=1; StopInd =0 stringloop: do while (VarNumber<4) 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='(E15.8)') & B_n(k) if (B_n(k)/=0.) B_n(k)=log(B_n(k)) elseif (VarNumber==2) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & alpha_n(k) elseif (VarNumber==3) then read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') & E_an(k) else call stop_it("No such VarNumber!") endif VarNumber=VarNumber+1 StartInd=StopInd endif if (StartInd==80) exit enddo stringloop ! ! Increase reaction counter by one ! k=k+1 endif endif endif endif enddo dataloop3 1012 continue close(file_id) ! endsubroutine read_reactions !*********************************************************************** subroutine get_reaction_rate(f,vreact_p,vreact_m,p) ! ! This subroutine calculates forward and reverse reaction rates, ! if chem.inp file exists. ! For more details see Chemkin Theory Manual ! ! NILS: This routine is by far the most CPU intencive routine in the code, ! NILS: so one should maybe put some effort into optimizing it more. ! ! 17-mar-08/natalia: coded ! 11-nov-10/julien: optimized, reaction rates are calculated under a ! logarithmic form ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(nx,nreactions), intent(out) :: vreact_p, vreact_m ! type (pencil_case) :: p real, dimension(nx) :: dSR=0., dHRT=0., Kp, Kc real, dimension(nx) :: prod1, prod2 real, dimension(nx) :: kf=0., kr=0. real, dimension(nx) :: rho_cgs, p_atm real, dimension(nx) :: mix_conc integer :: k, reac, i real :: sum_tmp=0., ddd real :: Rcal, Rcal1, lnRgas, l10, lnp_atm logical, save :: lwrite_first=.true. character(len=fnlen) :: input_file="./data/react.out" integer :: file_id=123 real :: B_n_0, alpha_n_0, E_an_0 real, dimension(nx) :: kf_0, Pr, sum_sp real, dimension(nx) :: Fcent, ccc, nnn, lnPr, FF, tmpF real, dimension(nx) :: TT1_loc ! ! Check which reactions rate method we will use ! if (reac_rate_method == 'chemkin') then ! TT1_loc = p%TT1 ! if (lwrite_first) open (file_id,file=input_file) ! ! p is in atm units; atm/bar=1./10.13 ! Rcal = Rgas_unit_sys/4.14*1e-7 Rcal1 = 1./Rcal lnRgas = log(Rgas) l10 = log(10.) rho_cgs = p%rho*unit_mass/unit_length**3 lnp_atm = log(1e6*unit_length**3/unit_energy) p_atm = 1e6*(unit_length**3)/unit_energy ! ! calculation of the reaction rate ! do reac = 1,nreactions ! ! Find the product of the species molar consentrations (where ! each molar consentration is taken to the power of the number) ! prod1 = 1. prod2 = 1. do k = 1,nchemspec if (abs(Sijp(k,reac)) == 1) then prod1 = prod1*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass)) elseif (abs(Sijp(k,reac)) == 2) then prod1 = prod1*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass))*(f(l1:l2,m,n,ichemspec(k)) & *rho_cgs(:)/species_constants(k,imass)) elseif (abs(Sijp(k,reac)) > 0) then prod1 = prod1*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass))**Sijp(k,reac) endif enddo do k = 1,nchemspec if (abs(Sijm(k,reac)) == 1.0) then prod2 = prod2*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass)) elseif (abs(Sijm(k,reac)) == 2.0) then prod2 = prod2*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass))*(f(l1:l2,m,n,ichemspec(k)) & *rho_cgs(:)/species_constants(k,imass)) elseif (abs(Sijm(k,reac)) > 0.0) then prod2 = prod2*(f(l1:l2,m,n,ichemspec(k))*rho_cgs(:) & /species_constants(k,imass))**Sijm(k,reac) endif enddo ! ! Find forward rate constant for reaction 'reac' ! if (latmchem) then if ((B_n(reac) == 0.) .and. (alpha_n(reac) == 0.) & .and. (E_an(reac) == 0.)) then do i = 1,nx call calc_extra_react(f,reac,kf(i),i,m,n,p) enddo kf = log(kf) else kf = B_n(reac)+alpha_n(reac)*p%lnTT-E_an(reac)*TT1_loc endif else kf = B_n(reac)+alpha_n(reac)*p%lnTT-E_an(reac)*Rcal1*TT1_loc endif ! ! Find backward rate constant for reaction 'reac' ! dSR = 0. dHRT = 0. sum_tmp = 0. do k = 1,nchemspec dSR = dSR+(Sijm(k,reac) -Sijp(k,reac))*p%S0_R(:,k) dHRT = dHRT+(Sijm(k,reac)-Sijp(k,reac))*p%H0_RT(:,k) sum_tmp = sum_tmp+(Sijm(k,reac)-Sijp(k,reac)) enddo Kp = dSR-dHRT ! if (sum_tmp == 0.) then Kc = Kp else Kc = Kp+sum_tmp*(lnp_atm-p%lnTT-lnRgas) endif ! ! Multiply by third body reaction term ! if (minval(a_k4(:,reac)) < impossible) then sum_sp = 0. do k = 1,nchemspec sum_sp = sum_sp+a_k4(k,reac)*f(l1:l2,m,n,ichemspec(k)) & *rho_cgs(:)/species_constants(k,imass) enddo mix_conc = sum_sp else sum_sp = 1. mix_conc = rho_cgs(:)*p%mu1(:)/unit_mass endif ! ! The Lindeman approach to the fall of reactions ! if (maxval(abs(low_coeff(:,reac))) > 0.) then B_n_0 = low_coeff(1,reac) alpha_n_0 = low_coeff(2,reac) E_an_0 = low_coeff(3,reac) kf_0(:) = B_n_0+alpha_n_0*p%lnTT(:)-E_an_0*Rcal1*TT1_loc(:) Pr = exp(kf_0-kf)*mix_conc kf = kf+log(Pr/(1.+Pr)) elseif (maxval(abs(high_coeff(:,reac))) > 0.) then B_n_0 = high_coeff(1,reac) alpha_n_0 = high_coeff(2,reac) E_an_0 = high_coeff(3,reac) kf_0(:) = B_n_0+alpha_n_0*p%lnTT(:)-E_an_0*Rcal1*TT1_loc(:) Pr = exp(kf_0-kf)*mix_conc kf = kf-log(1.+Pr) endif ! ! The Troe approach ! if (maxval(abs(troe_coeff(:,reac))) > 0.) then Fcent = (1.-troe_coeff(1,reac))*exp(-p%TT(:)/troe_coeff(2,reac)) & +troe_coeff(1,reac)*exp(-p%TT(:)/troe_coeff(3,reac)) ccc = -0.4-0.67*log10(Fcent) nnn = 0.75-1.27*log10(Fcent) ddd = 0.14 lnPr = log10(Pr) tmpF = ((lnPr+ccc)/(nnn-ddd*(lnPr+ccc)))**2 tmpF = 1./(1.+tmpF) FF = tmpF*log10(Fcent) FF = FF*l10 kf = kf+FF endif ! ! Find forward (vreact_p) and backward (vreact_m) rate of ! progress variable. ! (vreact_p - vreact_m) is labeled q in the chemkin manual ! if (latmchem) then kr = kf else kr = kf-Kc endif if (lpencil_check_at_work) where (kr > 32) kr = kr / exp(real(nint(alog(kr)))) ! if (Mplus_case (reac)) then where (prod1 > 0.) vreact_p(:,reac) = prod1*exp(kf) elsewhere vreact_p(:,reac) = 0. endwhere where (prod2 > 0.) vreact_m(:,reac) = prod2*exp(kr) elsewhere vreact_m(:,reac) = 0. endwhere ! else where (prod1 > 0.) vreact_p(:,reac) = prod1*exp(kf)*sum_sp elsewhere vreact_p(:,reac) = 0. endwhere where (prod2 > 0.) vreact_m(:,reac) = prod2*exp(kr)*sum_sp elsewhere vreact_m(:,reac) = 0. endwhere endif ! if (.not. back(reac)) vreact_m(:,reac) = 0. enddo ! ! This part calculates forward and reverse reaction rates ! for the test case R->P ! For more details see Doom, et al., J. Comp. Phys., 226, 2007 ! elseif (reac_rate_method == '1step_test') then do i = 1,nx if (p%TT(i) > Tc) then vreact_p(i,reac) = f(l1,m,n,iux)*f(l1,m,n,iux)*p%rho(1)*Cp_const & /lambda_const*beta*(beta-1.)*(1.-f(l1-1+i,m,n,ichemspec(ipr))) else vreact_p(i,reac) = 0. endif enddo vreact_m(:,reac) = 0. ! ! Add alternative method for finding the reaction rates based on work by ! Roux et al. (2009) ! NILS: Should split the different methods for calculating the reaction ! NILS: rates into different subroutines at some point. ! elseif (reac_rate_method == 'roux') then call roux(f,p,vreact_p,vreact_m) endif ! if (lwrite_first .and. lroot) & print*,'get_reaction_rate: writing react.out file' lwrite_first = .false. ! endsubroutine get_reaction_rate !*********************************************************************** subroutine roux(f,p,vreact_p,vreact_m) ! ! nilshau: 2011.01.11 (coded) ! ! Calculate reaction rates based on the method of Roux et al. (2009). ! This is a single step reaction mechanism for propane: ! ! C3H8 + 5O2 +XN2 -> 3CO2 + 4H2O + XN2 ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(nx,nreactions), intent(out) :: vreact_p, vreact_m type (pencil_case) :: p ! real :: mC3H8, mO2, Rcal, f_phi, E_a real, save :: init_C3H8, init_O2 integer :: i_O2, i_C3H8, ichem_O2, ichem_C3H8, j logical :: lO2, lC3H8 logical, save :: lfirsttime=.true. real, dimension(nx) :: activation_energy, pre_exp, term1, term2 ! if (nreactions /= 1) & call fatal_error('roux','nreactions should always be 1.') ! ! Check that a global equivalence ratio is given at input ! if (global_phi == impossible) call fatal_error('roux', & 'global_phi must be given as input') ! Rcal = Rgas_unit_sys/4.14*1e-7 ! ! Find indeces for oxygen and propane ! call find_species_index('O2',i_O2,ichem_O2,lO2) call find_species_index('C3H8',i_C3H8,ichem_C3H8,lC3H8) ! ! Check that oxygen and propane exist and find their molar masses ! if (lO2) then mO2 = species_constants(ichem_O2,imass) else call fatal_error('roux','O2 is not defined!') endif if (lC3H8) then mC3H8 = species_constants(ichem_C3H8,imass) else call fatal_error('roux','C3H8 is not defined!') endif ! ! Find initial mass fractions ! if (lfirsttime) then do j = 1,nchemspec initial_massfractions(j) = f(l1,m1,n1,ichemspec(j)) if (lroot) print*,'initial_massfractions=',initial_massfractions enddo init_O2 = initial_massfractions(ichem_O2) init_C3H8 = initial_massfractions(ichem_C3H8) endif ! ! Find Laminar flame speed corrector based on equivalence ratio phi ! f_phi & = 0.5*(1+tanh((0.8-global_phi)/1.5)) & +2.11/4*(1+tanh((global_phi-0.11)/0.2)) & *(1+tanh((1.355-global_phi)/0.24)) ! ! Find the classical Arrhenius terms ! E_a = 31126 activation_energy = exp(-E_a*p%TT1/Rcal) pre_exp = 3.2916e10 ! ! Find density and mass fraction dependent terms ! term1 = (f(l1:l2,m,n,i_C3H8)*p%rho & /species_constants(i_C3H8,imass))**0.856 term2 = (f(l1:l2,m,n,i_O2)*p%rho /species_constants(i_O2,imass))**0.503 ! ! Use the above to find reaction terms ! vreact_p(:,1) = f_phi*pre_exp*term1*term2*activation_energy ! ! Set reaction rate to zero when mass fractions of propane of oxygen is ! very close to zero. ! where (f(l1:l2,m,n,i_C3H8) < 1e-12) vreact_p(:,1) = 0. endwhere where (f(l1:l2,m,n,i_O2) < 1e-12) vreact_p(:,1) = 0. endwhere vreact_m(:,1) = 0. ! ! Print debugging output ! if (lfirsttime .and. lroot) then print*,'i_O2, i_C3H8, ichem_O2, ichem_C3H8=', & i_O2, i_C3H8, ichem_O2, ichem_C3H8 print*,'lO2, lC3H8=',lO2, lC3H8 print*,'init_C3H8,init_O2,mO2,mC3H8=',init_C3H8,init_O2,mO2,mC3H8 endif ! if (lfirsttime) lfirsttime = .false. ! endsubroutine roux !*********************************************************************** subroutine calc_reaction_term(f,p) ! ! Calculation of the reaction term ! use Diagnostics, only: sum_mn_name, max_mn_name ! real :: alpha, eps real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(nx,mreactions) :: vreactions, vreactions_p, vreactions_m real, dimension(nx,nchemspec) :: xdot real, dimension(nx) :: rho1 real, dimension(nx,nchemspec) :: molm type (pencil_case) :: p integer :: k,j,ii integer :: i1=1, i2=2, i3=3, i4=4, i5=5, i6=6, i7=7, i8=8, i9=9, i10=10 integer :: i11=11, i12=12, i13=13, i14=14, i15=15, i16=16, i17=17, i18=18, i19=19 ! eps = sqrt(epsilon(alpha)) ! p%DYDt_reac = 0. rho1 = 1./p%rho if (lcheminp .and. (.not. l1step_test)) then do k = 1,nchemspec molm(:,k) = rho1*species_constants(k,imass) enddo else molm = 1. endif ! ! if we do reactions, we must calculate the reaction speed vector ! outside the loop where we multiply it by the stoichiometric matrix ! if (.not. lcheminp) then ! ! Axel' case ! do j = 1,nreactions if (lkreactions_alpha) then alpha = kreactions_alpha(j) else alpha = 1. endif vreactions_p(:,j) = alpha*kreactions_p(j)*kreactions_z(n,j) vreactions_m(:,j) = alpha*kreactions_m(j)*kreactions_z(n,j) do k = 1,nchemspec vreactions_p(:,j) = vreactions_p(:,j)* & f(l1:l2,m,n,ichemspec(k))**Sijm(k,j) vreactions_m(:,j) = vreactions_m(:,j)* & f(l1:l2,m,n,ichemspec(k))**Sijp(k,j) enddo enddo vreactions_m = -vreactions_m vreactions_p = -vreactions_p else ! ! Chemkin data case ! call get_reaction_rate(f,vreactions_p,vreactions_m,p) endif ! ! Calculate rate of reactions (labeled q in the chemkin manual) ! vreactions = vreactions_p-vreactions_m ! ! Calculate production rate for all species k (called \dot(\omega)_k ! in the chemkin manual) ! xdot = 0. do k = 1,nchemspec do j = 1,nreactions xdot(:,k) = xdot(:,k)-stoichio(k,j)*vreactions(:,j)*molm(:,k) enddo enddo p%DYDt_reac = xdot*unit_time ! ! For diagnostics ! if (lchemistry_diag) then do k = 1,nchemspec do j = 1,nreactions net_react_p(k,j) = net_react_p(k,j)+stoichio(k,j) & *sum(vreactions_p(:,j)) net_react_m(k,j) = net_react_m(k,j)+stoichio(k,j) & *sum(vreactions_m(:,j)) enddo enddo endif ! ! NH: ! ! Sums for diagnostics ! !sum_omega=0. !sum_Y=0. !do k=1,nchemspec ! sum_omega=sum_omega+maxval(p%DYDt_reac(:,k)) ! sum_Y=sum_Y+maxval(f(l1:l2,m,n,ichemspec(k))) !enddo ! ! Calculate diagnostic quantities ! if (ldiagnos) then do ii=1,nchemspec if (idiag_dYm(ii) /= 0) then call sum_mn_name(p%DYDt_reac(:,ii),idiag_dYm(ii)) endif if (idiag_dYmax(ii) /= 0) then call max_mn_name(abs(p%DYDt_reac(:,ii)),idiag_dYmax(ii)) endif if (idiag_hm(ii) /= 0) then call sum_mn_name(p%H0_RT(:,ii)*Rgas* & p%TT(:)/species_constants(ii,imass),idiag_hm(ii)) endif enddo endif ! endsubroutine calc_reaction_term !*********************************************************************** subroutine write_net_reaction ! ! write net reactions to file ! open (1,file=trim(datadir)//'/net_reactions.dat',position='append') write (1,*) t write (1,'(8e10.2)') net_react_p, net_react_m close (1) ! endsubroutine write_net_reaction !*********************************************************************** subroutine calc_collision_integral(omega,lnTst,Omega_kl) ! ! Get coefficients for calculating of the collision integral ! This routine is called from calc_diff_visc_coeff, which again is called from ! calc_for_chem_mixture, which is why we work on full chunks of arrays here. ! ! 03-apr-08/natalia: coded ! character(len=*), intent(in) :: omega real, dimension(mx,my,mz), intent(in) :: lnTst real, dimension(mx,my,mz), intent(out) :: Omega_kl integer :: i real, dimension(8) :: aa ! select case (omega) case ('Omega11') aa(1) = 6.96945701E-1 aa(2) = 3.39628861E-1 aa(3) = 1.32575555E-2 aa(4) = -3.41509659E-2 aa(5) = 7.71359429E-3 aa(6) = 6.16106168E-4 aa(7) = -3.27101257E-4 aa(8) = 2.51567029E-5 case ('Omega22') aa(1) = 6.33225679E-1 aa(2) = 3.14473541E-1 aa(3) = 1.78229325E-2 aa(4) = -3.99489493E-2 aa(5) = 8.98483088E-3 aa(6) = 7.00167217E-4 aa(7) = -3.82733808E-4 aa(8) = 2.97208112E-5 case default call stop_it('Insert Omega_kl') endselect ! Omega_kl = 0. do i = 1,8 Omega_kl = Omega_kl+aa(i)*(lnTst)**(i-1) enddo Omega_kl = 1./Omega_kl ! endsubroutine calc_collision_integral !*********************************************************************** subroutine calc_diff_visc_coef(f) ! ! Calculation of the binary diffusion coefficients and the species viscosities. ! This routind is called from calc_for_chem_mixture, ! which is why we work on full chunks of arrays here. ! real, dimension(mx,my,mz,mfarray) :: f intent(in) :: f real, dimension(mx,my,mz) :: Omega_kl, prefactor real, dimension(mx,my,mz) :: lnTjk, lnTk_array integer :: k, j, j2, j3 real :: eps_jk, sigma_jk, m_jk, delta_jk, delta_st real :: Na=6.022E23, tmp_local, tmp_local2, delta_jk_star character(len=7) :: omega ! ! Find binary diffusion coefficients ! tmp_local = 3./16.*sqrt(2.*k_B_cgs**3/pi) do j3 = nn1,nn2 do j2 = mm1,mm2 prefactor(ll1:ll2,j2,j3) = tmp_local*sqrt(TT_full(ll1:ll2,j2,j3)) & *unit_length**3/(Rgas_unit_sys*rho_full(ll1:ll2,j2,j3)) enddo enddo ! omega = 'Omega11' ! ! Check if we use fixed Schmidt number for speeding up calculations ! if (.not. lfix_Sc) then ! ! Check if we use simplified version of the binary diffusion calculation ! if (ldiffusion .and. (.not. lDiff_simple) .and. (.not. lDiff_lewis)) then ! ! Do non-simplified binary diffusion coefficient ! do k = 1,nchemspec do j = k,nchemspec ! Account for the difference between eq. 5-4 and 5-31 in the Chemkin theory ! manual ! if (j /= k) then eps_jk = sqrt(tran_data(j,2)*tran_data(k,2)) sigma_jk = 0.5*(tran_data(j,3)+tran_data(k,3))*1e-8 m_jk = (species_constants(j,imass)*species_constants(k,imass)) & /(species_constants(j,imass)+species_constants(k,imass))/Na delta_jk = 0.5*tran_data(j,4)*tran_data(k,4)*1e-18*1e-18 else eps_jk = tran_data(j,2) sigma_jk = tran_data(j,3)*1e-8 m_jk = species_constants(j,imass)/(2*Na) delta_jk = 0.5*(tran_data(j,4)*1e-18)*(tran_data(j,4)*1e-18) endif ! ! Loop over all grid points ! do j3 = nn1,nn2 do j2 = mm1,mm2 if (ltemperature_nolog) then lnTjk(ll1:ll2,j2,j3) = log(f(ll1:ll2,j2,j3,ilnTT)/eps_jk) else lnTjk(ll1:ll2,j2,j3) = f(ll1:ll2,j2,j3,ilnTT)-log(eps_jk) endif ! Omega_kl(ll1:ll2,j2,j3) = & 1./(6.96945701E-1 +3.39628861E-1*lnTjk(ll1:ll2,j2,j3) & +1.32575555E-2*lnTjk(ll1:ll2,j2,j3)*lnTjk(ll1:ll2,j2,j3) & -3.41509659E-2*lnTjk(ll1:ll2,j2,j3)**3 & +7.71359429E-3*lnTjk(ll1:ll2,j2,j3)**4 & +6.16106168E-4*lnTjk(ll1:ll2,j2,j3)**5 & -3.27101257E-4*lnTjk(ll1:ll2,j2,j3)**6 & +2.51567029E-5*lnTjk(ll1:ll2,j2,j3)**7) delta_jk_star = delta_jk/(eps_jk*k_B_cgs*sigma_jk**3) ! Omega_kl(ll1:ll2,j2,j3) = Omega_kl(ll1:ll2,j2,j3) & +0.19*delta_jk_star*delta_jk_star/(TT_full(ll1:ll2,j2,j3)/eps_jk) if (j /= k) then Bin_Diff_coef(ll1:ll2,j2,j3,k,j) = prefactor(ll1:ll2,j2,j3)/mu1_full(ll1:ll2,j2,j3) & /(sqrt(m_jk)*sigma_jk*sigma_jk*Omega_kl(ll1:ll2,j2,j3)) else Bin_Diff_coef(ll1:ll2,j2,j3,k,j) = prefactor(ll1:ll2,j2,j3) & /(sqrt(m_jk)*sigma_jk*sigma_jk*Omega_kl(ll1:ll2,j2,j3))*species_constants(k,imass) ! endif enddo enddo enddo enddo ! do j3 = nn1,nn2 do j2 = mm1,mm2 do k = 1,nchemspec do j = 1,k-1 Bin_Diff_coef(ll1:ll2,j2,j3,k,j) = Bin_Diff_coef(ll1:ll2,j2,j3,j,k) enddo enddo enddo enddo ! else if (.not. lDiff_simple .and. .not. lDiff_lewis) then Bin_Diff_coef = 0. endif endif endif ! ! Calculate viscosity ! ! if (visc_const==impossible) then omega = 'Omega22' tmp_local = 5./16.*sqrt(k_B_cgs/(Na*pi)) ! do k = 1,nchemspec tmp_local2 = sqrt(species_constants(k,imass))/ & ((tran_data(k,3)*1e-8)*(tran_data(k,3)*1e-8))*tmp_local ! ! 1 Debye = 10**(-18) esu -> (1e-18*tran_data(k,4)) ! delta_st = (1e-18*tran_data(k,4))*(1e-18*tran_data(k,4))/2./ & (tran_data(k,2)*k_B_cgs*(tran_data(k,3)*1e-8)**3) ! if (ltemperature_nolog) then lnTk_array = log(f(:,:,:,ilnTT)/tran_data(k,2)) else lnTk_array = f(:,:,:,ilnTT)-log(tran_data(k,2)) endif call calc_collision_integral(omega,lnTk_array,Omega_kl) ! do j3 = nn1,nn2 do j2 = mm1,mm2 species_viscosity(ll1:ll2,j2,j3,k) = sqrt(TT_full(ll1:ll2,j2,j3)) & /(Omega_kl(ll1:ll2,j2,j3) & +0.2*delta_st*delta_st/(TT_full(ll1:ll2,j2,j3) & /tran_data(k,2)))*tmp_local2 & /(unit_mass/unit_length/unit_time) enddo enddo enddo ! endif ! endsubroutine calc_diff_visc_coef !*********************************************************************** subroutine calc_therm_diffus_coef(f) ! ! Calculate the thermal diffusion coefficient based on equation 5-17 in ! the Chemkin theory manual ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz,nchemspec) :: species_cond real, dimension(mx) :: tmp_val, ZZ, FF, tmp_sum, tmp_sum2 real, dimension(mx) :: AA, BB, f_tran, f_rot, f_vib real, dimension(mx) :: Cv_vib_R, T_st, pi_1_5, pi_2 real :: Cv_rot_R, Cv_tran_R intent(in) :: f integer :: j2, j3,k ! ! call timing('calc_therm_diffus_coef','just entered') ! pi_2 = pi*pi pi_1_5 = pi*sqrt(pi) ! ! With lspecies_cond_simplified, species conduction is calculated ! according to the book "Transport phenomena" by Bird, Warren, & Lightfoot ! page 276, Eq.(9.3-15). ! do j3 = nn1,nn2 do j2 = mm1,mm2 tmp_sum = 0. tmp_sum2 = 0. do k = 1,nchemspec if (lspecies_cond_simplified) then species_cond(:,j2,j3,k) = (species_viscosity(:,j2,j3,k))*& ( cp_spec_glo(:,j2,j3,k) + 1.25*Rgas/species_constants(k,imass) ) else ! ! Check if the molecule is a single atom (0), linear (1) or non-linear (2). ! if (tran_data(k,1) == 0.) then Cv_tran_R = 1.5 Cv_rot_R = 0. Cv_vib_R = 0. elseif (tran_data(k,1) == 1.) then Cv_tran_R = 1.5 Cv_rot_R = 1. Cv_vib_R = cv_R_spec_full(:,j2,j3,k)-2.5 elseif (tran_data(k,1) == 2.) then Cv_tran_R = 1.5 Cv_rot_R = 1.5 Cv_vib_R = cv_R_spec_full(:,j2,j3,k)-3. else Cv_tran_R = 0 Cv_rot_R = 0 Cv_vib_R = 0 call fatal_error('calc_therm_diffus_coef','No such tran_data!') endif ! ! The rotational and vibrational contributions are zero for the single ! atom molecules but not for the linear or non-linear molecules ! if (tran_data(k,1) > 0. .and. (.not. lfix_Sc)) then tmp_val = Bin_Diff_coef(:,j2,j3,k,k)*rho_full(:,j2,j3) & /species_viscosity(:,j2,j3,k) AA = 2.5-tmp_val T_st = tran_data(k,2)/298. FF = 1.+pi_1_5/2.*sqrt(T_st)+(pi_2/4.+2.) & *(T_st)+pi_1_5*(T_st)**1.5 ZZ = tran_data(k,6)*FF T_st = tran_data(k,2)/TT_full(:,j2,j3) FF = 1.+pi_1_5/2.*sqrt(T_st)+(pi_2/4.+2.) & *(T_st)+pi_1_5*(T_st)**1.5 ZZ = ZZ/FF BB = ZZ+2./pi*(5./3.*Cv_rot_R+tmp_val) f_tran = 2.5*(1.- 2./pi*Cv_rot_R/Cv_tran_R*AA/BB) f_rot = tmp_val*(1+2./pi*AA/BB) f_vib = tmp_val else f_tran = 2.5 f_rot = 0.0 f_vib = 0.0 endif species_cond(:,j2,j3,k) = (species_viscosity(:,j2,j3,k)) & /(species_constants(k,imass)/unit_mass)*Rgas* & (f_tran*Cv_tran_R+f_rot*Cv_rot_R & +f_vib*Cv_vib_R) endif ! ! tmp_sum and tmp_sum2 are used later to find the mixture averaged ! conductivity. ! tmp_sum = tmp_sum +XX_full(:,j2,j3,k)*species_cond(:,j2,j3,k) tmp_sum2 = tmp_sum2 +XX_full(:,j2,j3,k)/species_cond(:,j2,j3,k) enddo ! ! Find the mixture averaged conductivity ! where (tmp_sum2 <= 0.) lambda_full(:,j2,j3) = 0. elsewhere lambda_full(:,j2,j3) = 0.5*(tmp_sum+1./tmp_sum2) endwhere enddo enddo ! call keep_compiler_quiet(f) call timing('calc_therm_diffus_coef','just finished') ! call keep_compiler_quiet(f) ! endsubroutine calc_therm_diffus_coef !*********************************************************************** subroutine calc_diffusion_term(f,p) ! ! Calculate diffusion term, p%DYDt_diff ! ! 22-jun-10/julien: Evaluation of mass diffusion fluxes using Fick's ! law for simplified diffusion using constant Lewis numbers ! use Sub, only: del2,grad,dot_mn ! real, dimension(mx,my,mz,mfarray) :: f type (pencil_case) :: p real, dimension(nx) :: Xk_Yk real, dimension(nx,3) :: gDiff_full_add, gchemspec, gXk_Yk real, dimension(nx) :: del2chemspec real, dimension(nx) :: diff_op, diff_op1, diff_op2, diff_op3, del2XX, del2lnpp real, dimension(nx) :: glnpp_gXkYk, glnrho_glnpp, gD_glnpp, glnpp_glnpp real, dimension(nx) :: sum_gdiff=0., gY_sumdiff, glnmu_glnpp real, dimension(nx,3) :: sum_diff=0., dk_D real :: diff_k integer :: k,i ! intent(in) :: f ! p%DYDt_diff = 0. diff_k = chem_diff ! ! Loop over all chemical species. ! do k = 1,nchemspec ! ! Simple chemistry (not using ChemKin formalism). ! Eliminate need for p%glnrho if density is not included. ! if (.not. lcheminp) then if (chem_diff /= 0.) then diff_k = chem_diff*chem_diff_prefactor(k) if (headtt) print*,'dchemistry_dt: k,diff_k=',k,diff_k call del2(f,ichemspec(k),del2chemspec) call grad(f,ichemspec(k),gchemspec) if (ldensity) then call dot_mn(p%glnrho,gchemspec,diff_op) diff_op = diff_op+del2chemspec else diff_op = del2chemspec endif p%DYDt_diff(:,k) = diff_k*diff_op endif else ! ! Detailed chemistry and transport using CHEMKIN formalism. ! if (ldiffusion) then ! ! Calculate diffusion coefficient gradients gDiff_full_add in 3 cases: ! 1) Simplified diffusion coefficients ! 2) Constant Lewis numbers and heat conductivity ! 3) Detailed transport ! if (lDiff_simple) then do i = 1,3 gDiff_full_add(:,i) = p%Diff_penc_add(:,k) & *(0.7*p%glnTT(:,i)-p%glnrho(:,i)) enddo elseif (lDiff_lewis .and. lew_exist) then do i = 1,3 gDiff_full_add(:,i) = & (p%glambda(:,i)*p%lambda1(:)-p%glncp(:,i)-p%glnrho(:,i)) & *p%Diff_penc_add(:,k) enddo elseif (lDiff_lewis .and. l1step_test) then do i = 1,3 gDiff_full_add(:,i) = & (p%glambda(:,i)*p%lambda1(:)-p%glncp(:,i)-p%glnrho(:,i)) & *p%Diff_penc_add(:,k) enddo else call grad(Diff_full_add(:,:,:,k),gDiff_full_add) endif ! ! Calculate the terms needed by the diffusion fluxes in 3 cases: ! 1) Fickian diffusion law (gradient of species MASS fractions) ! 2) Simplified fluxes (gradient of species MOLAR fractions) ! 3) Detailed transport (with pressure gradients included) ! if (lDiff_fick) then call del2(f,ichemspec(k),del2chemspec) call grad(f,ichemspec(k),gchemspec) call dot_mn(p%glnrho,gchemspec,diff_op1) call dot_mn(gDiff_full_add,gchemspec,diff_op2) elseif (lFlux_simple) then call del2(XX_full(:,:,:,k),del2XX) call dot_mn(p%glnrho,p%gXXk(:,:,k),diff_op1) call dot_mn(gDiff_full_add,p%gXXk(:,:,k),diff_op2) call dot_mn(p%glnmu,p%gXXk(:,:,k),diff_op3) else call del2(XX_full(:,:,:,k),del2XX) call dot_mn(p%glnrho,p%gXXk(:,:,k),diff_op1) call dot_mn(gDiff_full_add,p%gXXk(:,:,k),diff_op2) call dot_mn(p%glnmu,p%gXXk(:,:,k),diff_op3) call dot_mn(p%glnpp,p%glnpp,glnpp_glnpp) do i = 1,3 gXk_Yk(:,i) = p%gXXk(:,i,k)-p%gYYk(:,i,k) enddo del2lnpp = p%del2pp/p%pp-glnpp_glnpp Xk_Yk = XX_full(l1:l2,m,n,k)-f(l1:l2,m,n,ichemspec(k)) call dot_mn(p%glnrho,p%glnpp,glnrho_glnpp) call dot_mn(gDiff_full_add,p%glnpp,gD_glnpp) call dot_mn(gXk_Yk,p%glnpp,glnpp_gXkYk) call dot_mn(p%glnmu,p%glnpp,glnmu_glnpp) endif ! ! Calculate the diffusion fluxes and dk_D in 3 cases: ! 1) Fickian diffusion law (gradient of species MASS fractions) ! 2) Simplified diffusion fluxes (only gradient of species MOLAR fractions) ! 3) Detailed transport (with pressure gradients included) ! Note that the ratio Wk/Wm is introduced here and not during the calculation ! of species diffusion coefficients as before. It indeed depends on the diffusion ! flux formulation and not on the diffusive properties of each species. ! if (lDiff_fick) then p%DYDt_diff(:,k) = p%Diff_penc_add(:,k)*(del2chemspec+diff_op1) & +diff_op2 do i = 1,3 dk_D(:,i) = p%Diff_penc_add(:,k)*gchemspec(:,i) enddo elseif (lFlux_simple) then p%DYDt_diff(:,k) = p%Diff_penc_add(:,k)*p%mukmu1(:,k) & *(del2XX+diff_op1-diff_op3)+ & p%mukmu1(:,k)*diff_op2 do i = 1,3 dk_D(:,i) = p%Diff_penc_add(:,k)*p%mukmu1(:,k)*p%gXXk(:,i,k) enddo else p%DYDt_diff(:,k) = p%Diff_penc_add(:,k)*p%mukmu1(:,k) & *(del2XX+diff_op1-diff_op3)+ & p%mukmu1(:,k)*diff_op2+ & p%Diff_penc_add(:,k)*p%mukmu1(:,k)*Xk_Yk(:) & *(del2lnpp+glnrho_glnpp-glnmu_glnpp)+ & Xk_Yk(:)*p%mukmu1(:,k)*gD_glnpp+p%Diff_penc_add(:,k) & *p%mukmu1(:,k)*glnpp_gXkYk do i = 1,3 dk_D(:,i) = p%Diff_penc_add(:,k)*p%mukmu1(:,k)*(p%gXXk(:,i,k)+ & Xk_Yk(:)*p%glnpp(:,i)) enddo endif ! if (ldiff_corr) then do i = 1,3 sum_diff(:,i) = sum_diff(:,i)+dk_D(:,i) enddo sum_gdiff(:) = sum_gdiff(:)+ p%DYDt_diff(:,k) endif endif endif enddo ! ! Adding correction diffusion velocity to ensure mass balance ! if (ldiffusion .and. ldiff_corr) then do k = 1,nchemspec call grad(f,ichemspec(k),gchemspec) call dot_mn(gchemspec(:,:),sum_diff,gY_sumdiff) p%DYDt_diff(:,k) = p%DYDt_diff(:,k)- & (gY_sumdiff+f(l1:l2,m,n,ichemspec(k))*sum_gdiff(:)) enddo endif ! endsubroutine calc_diffusion_term !*********************************************************************** subroutine calc_heatcond_chemistry(f,df,p) ! ! Calculate gamma*chi*(del2lnT+gradlnTT.grad(lnT+lnrho+lncp+lnchi)) ! ! 29-feb-08/natalia: coded ! use Sub, only: dot, del2 ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz,mvar) :: df type (pencil_case) :: p real, dimension(nx) :: g2TT, g2TTlambda=0., tmp1, del2TT ! if (ltemperature_nolog) then call dot(p%gTT,p%glambda,g2TTlambda) call del2(f,iTT,del2TT) else call dot(p%glnTT,p%glambda,g2TTlambda) call dot(p%glnTT,p%glnTT,g2TT) endif ! ! Add heat conduction to RHS of temperature equation ! ! if (l1step_test .or. lSmag_heat_transport) then if (l1step_test) then tmp1 = p%lambda(:)*(p%del2lnTT+g2TT)*p%cv1/p%rho(:) else if (ltemperature_nolog) then tmp1 = (p%lambda(:)*del2TT+g2TTlambda)*p%cv1/p%rho(:) ! tmp1 = (p%lambda(:)*p%del2lnTT+g2TTlambda)*p%cv1/p%rho(:) df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + tmp1 else tmp1 = (p%lambda(:)*(p%del2lnTT+g2TT)+g2TTlambda)*p%cv1/p%rho(:) df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + tmp1 endif endif ! call keep_compiler_quiet(f) ! endsubroutine calc_heatcond_chemistry !*********************************************************************** subroutine get_RHS_Y_full(RHS_Y) ! real, dimension(mx,my,mz,nchemspec) :: RHS_Y intent(out) :: RHS_Y ! RHS_Y = RHS_Y_full ! endsubroutine get_RHS_Y_full !*********************************************************************** subroutine get_cs2_full(cs2_full) ! real, dimension(mx,my,mz) :: cs2_full intent(out) :: cs2_full integer :: j,k ! do j = 1,my do k = 1,mz if (minval(cv_full(:,j,k)) <= 0) then cs2_full(:,j,k) = 0. else cs2_full(:,j,k) = cp_full(:,j,k)/cv_full(:,j,k)*mu1_full(:,j,k) & *TT_full(:,j,k)*Rgas endif enddo enddo ! endsubroutine get_cs2_full !*********************************************************************** subroutine get_cs2_slice(f,slice,dir,index) ! ! Find a slice of the speed of sound ! ! 10-dez-09/nils: coded ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(:,:), intent(out) :: slice integer, intent(in) :: index, dir ! if (dir == 1) then slice = cp_full(index,m1:m2,n1:n2)/cv_full(index,m1:m2,n1:n2) & *mu1_full(index,m1:m2,n1:n2)*TT_full(index,m1:m2,n1:n2)*Rgas elseif (dir == 2) then slice = cp_full(l1:l2,index,n1:n2)/cv_full(l1:l2,index,n1:n2) & *mu1_full(l1:l2,index,n1:n2)*TT_full(l1:l2,index,n1:n2)*Rgas elseif (dir == 3) then slice = cp_full(l1:l2,m1:m2,index)/cv_full(l1:l2,m1:m2,index) & *mu1_full(l1:l2,m1:m2,index)*TT_full(l1:l2,m1:m2,index)*Rgas else call fatal_error('get_cs2_slice','No such dir!') endif ! endsubroutine get_cs2_slice !*********************************************************************** subroutine get_gamma_full(gamma_full) ! real, dimension(mx,my,mz) :: gamma_full intent(out) :: gamma_full integer :: j,k ! do j = 1,my do k = 1,mz if (minval(cv_full(:,j,k)) <= 0) then gamma_full(:,j,k) = 0. else gamma_full(:,j,k) = cp_full(:,j,k)/cv_full(:,j,k) endif enddo enddo endsubroutine get_gamma_full !*********************************************************************** subroutine get_gamma_slice(f,slice,dir,index) ! ! Get a 2D slice of gamma ! ! 10-dez-09/Nils Erland L. Haugen: coded ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(:,:), intent(out) :: slice integer, intent(in) :: index, dir ! if (dir == 1) then slice = cp_full(index,m1:m2,n1:n2)/cv_full(index,m1:m2,n1:n2) elseif (dir == 2) then slice = cp_full(l1:l2,index,n1:n2)/cv_full(l1:l2,index,n1:n2) elseif (dir == 3) then slice = cp_full(l1:l2,m1:m2,index)/cv_full(l1:l2,m1:m2,index) else call fatal_error('get_gamma_slice','No such dir!') endif ! endsubroutine get_gamma_slice !*********************************************************************** subroutine air_field(f,PP) ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz) :: sum_Y, tmp real :: PP ! (in dynes = 1atm) ! logical :: emptyfile=.true. logical :: found_specie integer :: file_id=123, ind_glob, ind_chem character(len=80) :: ChemInpLine character(len=10) :: specie_string character(len=1) :: tmp_string integer :: i, j, k=1 real :: YY_k, air_mass, TT=300. real :: velx=0. real, dimension(nchemspec) :: stor2 integer, dimension(nchemspec) :: stor1 ! integer :: StartInd, StopInd, StartInd_1, StopInd_1 integer :: iostat logical :: airdat=.false.,airin=.false. ! air_mass = 0. StartInd_1 = 1 StopInd_1 = 0 ! inquire (file='air.dat',exist=airdat) inquire (file='air.in',exist=airin) if (airdat .and. airin) then call fatal_error('chemistry',& 'air.in and air.dat found. Please decide for one') endif if (airdat) open(file_id,file='air.dat') if (airin) open(file_id,file='air.in') ! if (lroot) print*, 'the following parameters and '//& 'species are found in air.dat (volume fraction fraction in %): ' dataloop: do read(file_id,'(80A)',IOSTAT=iostat) ChemInpLine if (iostat < 0) exit dataloop emptyFile=.false. StartInd_1=1; StopInd_1=0 StopInd_1=index(ChemInpLine,' ') specie_string=trim(ChemInpLine(1:StopInd_1-1)) tmp_string=trim(ChemInpLine(1:1)) if (tmp_string == '!' .or. tmp_string == ' ') then elseif (tmp_string == 'T') then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 read (unit=ChemInpLine(StartInd:StopInd),fmt='(E14.7)') TT if (lroot) print*, ' Temperature, K ', TT elseif (tmp_string == 'P') then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 read (unit=ChemInpLine(StartInd:StopInd),fmt='(E14.7)') PP if (lroot) print*, ' Pressure, Ba ', PP elseif (tmp_string == 'V') then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 read (unit=ChemInpLine(StartInd:StopInd),fmt='(E14.7)') velx if (lroot) print*, ' Velocity, cm/s ', velx else call find_species_index(specie_string,ind_glob,ind_chem,found_specie) if (found_specie) then StartInd=1; StopInd =0 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 StartInd=verify(ChemInpLine(StopInd:),' ')+StopInd-1 StopInd=index(ChemInpLine(StartInd:),' ')+StartInd-1 read (unit=ChemInpLine(StartInd:StopInd),fmt='(E15.8)') YY_k if (lroot) print*, ' volume fraction, %, ', YY_k, & species_constants(ind_chem,imass) if (species_constants(ind_chem,imass)>0.) then air_mass=air_mass+YY_k*0.01/species_constants(ind_chem,imass) endif if (StartInd==80) exit stor1(k)=ind_chem stor2(k)=YY_k k=k+1 endif endif enddo dataloop ! ! Stop if air.dat is empty ! if (emptyFile) call stop_it('The input file air.dat was empty!') air_mass=1./air_mass do j=1,k-1 f(:,:,:,ichemspec(stor1(j)))=stor2(j)*0.01 enddo sum_Y=0. do j=1,nchemspec sum_Y=sum_Y+f(:,:,:,ichemspec(j)) enddo do j=1,nchemspec f(:,:,:,ichemspec(j))=f(:,:,:,ichemspec(j))/sum_Y initial_massfractions(j)=f(l1,m1,n1,ichemspec(j)) enddo if (mvar < 5) then call fatal_error("air_field", "I can only set existing fields") endif if (.not. reinitialize_chemistry) then if (.not.lflame_front .and. .not.ltriple_flame .and. .not.lFlameMaster) then if (ltemperature_nolog) then f(:,:,:,iTT)=TT else f(:,:,:,ilnTT)=alog(TT)!+f(:,:,:,ilnTT) endif if (ldensity_nolog) then f(:,:,:,ilnrho)=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 else tmp=(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 f(:,:,:,ilnrho)=alog(tmp) endif if (nxgrid>1) f(:,:,:,iux)=f(:,:,:,iux)+init_ux endif endif if (init_from_file) then if (lroot) print*, 'Velocity field read from file, initialization' //& 'of density and temperature' if (ldensity_nolog) then f(:,:,:,ilnrho)=f(:,:,:,ilnrho)*(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 else tmp=f(:,:,:,ilnrho)*(PP/(k_B_cgs/m_u_cgs)*& air_mass/TT)/unit_mass*unit_length**3 f(:,:,:,ilnrho)=alog(tmp) endif if (ltemperature_nolog) then if (ldensity_nolog) then f(:,:,:,iTT)=(PP/(k_B_cgs/m_u_cgs)*& air_mass/f(:,:,:,ilnrho))/unit_mass*unit_length**3 else f(:,:,:,iTT)=(PP/(k_B_cgs/m_u_cgs)*& air_mass/exp(f(:,:,:,ilnrho)))/unit_mass*unit_length**3 endif else if (ldensity_nolog) then tmp=(PP/(k_B_cgs/m_u_cgs)*& air_mass/f(:,:,:,ilnrho))/unit_mass*unit_length**3 f(:,:,:,ilnTT)=alog(tmp)!+f(:,:,:,ilnTT) else tmp=(PP/(k_B_cgs/m_u_cgs)*& air_mass/exp(f(:,:,:,ilnrho)))/unit_mass*unit_length**3 f(:,:,:,ilnTT)=alog(tmp)!+f(:,:,:,ilnTT) endif endif if (velx/=0.) f(:,:,:,iux)=f(:,:,:,iux)+velx endif if (linit_temperature) then do i=1,mx if (x(i)<=init_x1) then f(i,:,:,ilnTT)=alog(init_TT1) endif if (x(i)>=init_x2) then f(i,:,:,ilnTT)=alog(init_TT2) endif if (x(i)>init_x1 .and. x(i)1) then f(:,:,:,iux)=f(:,:,:,iux)+velx endif if (lroot) print*, 'Air temperature, K', TT if (lroot) print*, 'Air pressure, dyn', PP if (lroot) print*, 'Air density, g/cm^3:' if (lroot) print '(E10.3)', PP/(k_B_cgs/m_u_cgs)*air_mass/TT if (lroot) print*, 'Air mean weight, g/mol', air_mass if (lroot) print*, 'R', k_B_cgs/m_u_cgs close(file_id) ! endsubroutine air_field !*********************************************************************** ! NSCBC boundary conditions !*********************************************************************** subroutine damp_zone_for_NSCBC(f,df) ! ! 16-jul-06/natalia: coded ! 24-jan-11/julien: modified ! buffer zone to damp the acustic waves!!!!!!!!!!! ! important for NSCBC ! Most of this routine is still hard-coded ! Should contain more tests to detect boundaries and set the reference conditions ! real, dimension(mx,my,mz,mfarray), intent(in) :: f real, dimension(mx,my,mz,mvar), intent(inout) :: df integer :: sz_x, sz_y, sz_z, ll1, ll2, i ! integer :: sz_r_x, sz_l_x, sz_r_y, sz_l_y, sz_r_z, sz_l_z real :: dt1, func_x !,func_y,func_z real :: ux_ref, uy_ref, uz_ref, lnTT_ref, lnrho_ref, gamma, cs real :: del ! logical :: lzone_y=.false.,lzone_z=.false. logical :: dir_damp1=.true. !, dir_damp2=.false., dir_damp3=.false. logical :: lright=.true., lleft=.true. ! ux_ref = 0. uy_ref = 0. uz_ref = 0. ! lnTT_ref=6.39693 lnTT_ref = log(init_TT1) lnrho_ref = log(init_pressure)-log(Rgas)-lnTT_ref-log(mu1_full(l1,m1,n1)) gamma = 1.4 cs = sqrt(gamma*init_pressure/exp(lnrho_ref)) ! ! The characteristic time scale is taken of the order of a CFL time scale ! dt1 = (ux_ref+cs)/(Lxyz(1)/nxgrid) del = 0.1 ! sz_x = int(del*nxgrid) sz_y = int(del*nygrid) sz_z = int(del*nzgrid) ll1 = l1 ll2 = l2 ! if (nxgrid /= 1 .and. dir_damp1) then ! ! On the right side ! if (lright) then ll1 = nxgrid-sz_x ll2 = nxgrid do i = l1,l2 if (x(i) >= xgrid(ll1)) then func_x = (x(i)-xgrid(ll1))**3/(xgrid(ll2)-xgrid(ll1))**3 df(i,m,n,iux) = df(i,m,n,iux)-func_x*(f(i,m,n,iux)-ux_ref)*dt1 df(i,m,n,iuy) = df(i,m,n,iuy)-func_x*(f(i,m,n,iuy)-uy_ref)*dt1 df(i,m,n,iuz) = df(i,m,n,iuz)-func_x*(f(i,m,n,iuz)-uz_ref)*dt1 ! df(i,m,n,ilnrho)=df(i,m,n,ilnrho)-func_x*(f(i,m,n,ilnrho)-lnrho_ref)*dt1 ! df(i,m,n,ilnTT)=df(i,m,n,ilnTT)-func_x*(f(i,m,n,ilnTT)-lnTT_ref)*dt1 endif enddo endif ! ! On the left side ! if (lleft) then ll1 = 1 ll2 = sz_x do i = l1,l2 if (x(i) <= xgrid(ll2)) then func_x = (x(i)-xgrid(ll2))**3/(xgrid(ll1)-xgrid(ll2))**3 df(i,m,n,iux) = df(i,m,n,iux)-func_x*(f(i,m,n,iux)-ux_ref)*dt1 df(i,m,n,iuy) = df(i,m,n,iuy)-func_x*(f(i,m,n,iuy)-uy_ref)*dt1 df(i,m,n,iuz) = df(i,m,n,iuz)-func_x*(f(i,m,n,iuz)-uz_ref)*dt1 ! df(i,m,n,ilnrho)=df(i,m,n,ilnrho)-func_x*(f(i,m,n,ilnrho)-lnrho_ref)*dt1 ! df(i,m,n,ilnTT)=df(i,m,n,ilnTT)-func_x*(f(i,m,n,ilnTT)-lnTT_ref)*dt1 endif enddo endif ! endif ! ! The following should be replaced to generalize the formula ! ! if (nygrid/=1 .and. dir_damp2) then ! ! if (sz_r_y<=m1) call fatal_error('to use ldamp_zone_NSCBC',& ! 'you should increase nygrid!') ! ! if ((m<=sz_l_y) .and. (m>=m1)) then ! func_y=(y(m)-y(sz_l_y))**3/(y(m1)-y(sz_l_y))**3 ! lzone_y=.true. ! elseif ((m>=sz_r_y) .and. (m<=m2)) then ! func_y= (y(m)-y(sz_r_y))**3/(y(m2)-y(sz_r_y))**3 ! lzone_y=.true. ! endif ! ! if (lzone_y) then ! df(sz_l_x:sz_r_x,m,n,iux)=df(sz_l_x:sz_r_x,m,n,iux)& ! -func_y*(f(sz_l_x:sz_r_x,m,n,iux)-ux_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,iuy)=df(sz_l_x:sz_r_x,m,n,iuy)& ! -func_y*(f(sz_l_x:sz_r_x,m,n,iuy)-uy_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,iuz)=df(sz_l_x:sz_r_x,m,n,iuz)& ! -func_y*(f(sz_l_x:sz_r_x,m,n,iuz)-uz_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,ilnrho)=df(sz_l_x:sz_r_x,m,n,ilnrho)& ! -func_y*(f(sz_l_x:sz_r_x,m,n,ilnrho)-lnrho_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,ilnTT)=df(sz_l_x:sz_r_x,m,n,ilnTT)& ! -func_y*(f(sz_l_x:sz_r_x,m,n,ilnTT)-lnTT_ref)*dt1 ! lzone_y=.false. ! endif ! endif ! ! if (nzgrid/=1 .and. dir_damp3) then ! if (sz_r_z<=n1) call fatal_error('to use ldamp_zone_NSCBC',& ! 'you should increase nzgrid!') ! ! if ((n<=sz_l_z) .and. (n>=n1)) then ! func_z=(z(n)-z(sz_l_z))**3/(z(n1)-z(sz_l_z))**3 ! lzone_z=.true. ! elseif ((n>=sz_r_z) .and. (n<=n2)) then ! func_z= (z(n)-z(sz_r_z))**3/(z(n2)-z(sz_r_z))**3 ! lzone_z=.true. ! endif ! ! if (lzone_z) then ! df(sz_l_x:sz_r_x,m,n,iux)=df(sz_l_x:sz_r_x,m,n,iux)& ! -func_z*(f(sz_l_x:sz_r_x,m,n,iux)-ux_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,iuy)=df(sz_l_x:sz_r_x,m,n,iuy)& ! -func_z*(f(sz_l_x:sz_r_x,m,n,iuy)-uy_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,iuz)=df(sz_l_x:sz_r_x,m,n,iuz)& ! -func_z*(f(sz_l_x:sz_r_x,m,n,iuz)-uz_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,ilnrho)=df(sz_l_x:sz_r_x,m,n,ilnrho)& ! -func_z*(f(sz_l_x:sz_r_x,m,n,ilnrho)-lnrho_ref)*dt1 ! df(sz_l_x:sz_r_x,m,n,ilnTT)=df(sz_l_x:sz_r_x,m,n,ilnTT)& ! -func_z*(f(sz_l_x:sz_r_x,m,n,ilnTT)-lnTT_ref)*dt1 ! lzone_z=.false. ! endif ! endif ! endsubroutine damp_zone_for_NSCBC !*********************************************************************** subroutine jacobn(f,jacob) ! ! Compute the jacobian, i.e. the matrix jacob(nchemspec x nchemspec) ! where jacob(i,j)=dv_i/dc_j ! v is the vector of dimension nchemspec of the rates dc_j/dt ! (the c_j being concentrations, stocked in f among other) ! ! 28-may-09/rplasson: coded ! ! exchange data real, dimension(mx,my,mz,mfarray) :: f real, dimension(mx,my,mz,nchemspec,nchemspec) :: jacob ! intent(in) :: f intent(out) :: jacob ! internal data ! ! indices integer :: i, j, k, l, ii ! ! temporary real :: tmp_p, tmp_m ! Code ! jacob = 0. ! ! identify module ! if (headtt .or. ldebug) print*,'jacobn: compute the jacobian matrix' ! if (lreactions) then do n = n1,n2 do m = m1,m2 do l = l1,l2 do i = 1,nchemspec do j = 1,nchemspec ! Compute dv_i/dc_j ! print*,"dv_",i,"/dc_",j,"=" do k = 1,nreactions ! Check if compound i participate in reaction k if (Sijp(i,k) /= Sijm(i,k)) then ! Compute the contribution of reaction k to dv_i/dc_j ! print*,"+(",(Sijp(i,k)-Sijm(i,k)),"k_",k,"(+/-)" tmp_p = (Sijp(i,k)-Sijm(i,k))*kreactions_p(k)*kreactions_z(n,k) tmp_m = (Sijm(i,k)-Sijp(i,k))*kreactions_m(k)*kreactions_z(n,k) do ii = 1,nchemspec ! Compute the contribution of compound ii in reaction k ! print*,"**-",Sijm(ii,k),Sijp(ii,k),"-**" if (ii /= j) then ! print*,"c_",ii,"^",Sijm(ii,k)," (/) ","c_",ii,"^",Sijp(ii,k) tmp_p = tmp_p*f(l,m,n,ichemspec(ii))**Sijm(ii,k) tmp_m = tmp_m*f(l,m,n,ichemspec(ii))**Sijp(ii,k) else if (Sijm(ii,k) == 0) then ! print*,"0*c_",ii tmp_p = 0. elseif (Sijm(ii,k) > 1) then ! print*,Sijm(ii,k),"*c_",ii,"^",(Sijm(ii,k)-1) tmp_p = Sijm(ii,k)*tmp_p*f(l,m,n,ichemspec(ii))**(Sijm(ii,k)-1) ! else ! print*,"c_",ii,"^0" endif ! print*," (/) " if (Sijp(ii,k) == 0) then ! print*,"0*c_",ii tmp_m = 0. elseif (Sijp(ii,k) > 1) then ! print*,Sijp(ii,k),"*c_",ii,"^",(Sijp(ii,k)-1) tmp_m = Sijp(ii,k)*tmp_m*f(l,m,n,ichemspec(ii))**(Sijp(ii,k)-1) ! else ! print*,"c_",ii,"^0" endif endif enddo ! print*,")" ! Add the contribution of reaction k to dv_i/dc_j jacob(l,m,n,i,j) = jacob(l,m,n,i,j)+tmp_p+tmp_m ! print*,"(=",tmp_p," (-) ",tmp_m,")" endif enddo enddo enddo enddo enddo enddo endif ! endsubroutine jacobn !*********************************************************************** subroutine get_mu1_slice(f,slice,grad_slice,index,sgn,direction) ! ! For the NSCBC boudary conditions the slice of mu1 at the boundary, and ! its gradient, is required. ! ! 10-dez-09/Nils Erland L. Haugen: coded ! use Deriv, only: der_onesided_4_slice_other ! real, dimension(:,:), intent(out) :: slice real, dimension(:,:), intent(out) :: grad_slice integer, intent(in) :: index, sgn, direction real, dimension(mx,my,mz,mfarray) :: f ! if (direction == 1) then slice = mu1_full(index,m1:m2,n1:n2) call der_onesided_4_slice_other(mu1_full,sgn,grad_slice,index,direction) elseif (direction == 2) then slice = mu1_full(l1:l2,index,n1:n2) call der_onesided_4_slice_other(mu1_full,sgn,grad_slice,index,direction) else slice = mu1_full(l1:l2,m1:m2,index) call der_onesided_4_slice_other(mu1_full,sgn,grad_slice,index,direction) endif ! endsubroutine get_mu1_slice !*********************************************************************** subroutine calc_extra_react(f,reac,kf_loc,i,mm,nn,p) ! ! character (len=*), intent(in) :: element_name real, dimension(mx,my,mz,mfarray) :: f type (pencil_case), intent(in) :: p integer, intent(in) :: reac, i, mm, nn real, intent(out) :: kf_loc real :: X_O2, X_N2, X_O2N2, X_M, X_H2O real :: K1, K2, K3, K4, KMT06 ! ! select case (reaction_name(reac)) case ('O=O3') X_O2 = f(l1+i-1,mm,nn,ichemspec(index_O2))*unit_mass & /species_constants(index_O2,imass)*p%rho(i) X_N2 = f(l1+i-1,mm,nn,ichemspec(index_N2))*unit_mass & /species_constants(index_N2,imass)*p%rho(i) kf_loc = 5.60D-34*X_O2*X_N2*((1./300.*p%TT(i))**(-2.6)) & +6.00D-34*X_O2**2*((1./300.*p%TT(i))**(-2.6)) case ('O1D=O') X_O2 = f(l1+i-1,mm,nn,ichemspec(index_O2))*unit_mass & /species_constants(index_O2,imass)*p%rho(i) X_N2 = f(l1+i-1,mm,nn,ichemspec(index_N2))*unit_mass & /species_constants(index_N2,imass)*p%rho(i) kf_loc = 3.20D-11*X_O2*exp(67.*p%TT1(i))+1.80D-11*X_N2*exp(107.*p%TT1(i)) case ('OH+CO=HO2') X_O2N2 = f(l1+i-1,mm,nn,ichemspec(index_O2N2))*unit_mass & /species_constants(index_O2N2,imass)*p%rho(i) kf_loc = 1.30D-13*(1+((0.6*index_O2N2)/(2.652E+19*(300.*p%TT1(i))))) case ('2HO2=H2O2') X_M = p%rho(i)*mu1_full(l1+i-1,mm,nn) X_H2O = f(l1+i-1,mm,nn,ichemspec(index_H2O))*unit_mass & /species_constants(index_H2O,imass)*p%rho(i) KMT06 = 1 + (1.4E-21 * EXP(2200.*p%TT1(i)) * X_H2O) kf_loc = 2.20D-13*KMT06*EXP(600.*p%TT1(i)) & + 1.90D-33*X_M*KMT06*EXP(980.*p%TT1(i)) case ('OH+HNO3=NO3') X_O2N2 = f(l1+i-1,mm,nn,ichemspec(index_O2N2))*unit_mass & /species_constants(index_O2N2,imass)*p%rho(i) K1 = 2.4E-14 * EXP(460.*p%TT1(i)) K3 = 6.5E-34 * EXP(1335.*p%TT1(i)) K4 = 2.7E-17 * EXP(2199.*p%TT1(i)) K2 = (K3 * X_O2N2) / (1 + (K3*X_O2N2/K4)) kf_loc = K1 + K2 case default if (lroot) print*,'reaction_name=', reaction_name(reac) call stop_it('calc_extra_react: Element not found!') endselect ! endsubroutine calc_extra_react !*********************************************************************** subroutine prerun_1D(f,prerun_dir) ! ! read snapshot file, possibly with mesh and time (if mode=1) ! 11-apr-97/axel: coded ! character(len=*) :: prerun_dir character(len=100) :: file character(len=10) :: processor real, dimension(mx,my,mz,mfarray) :: f real, dimension(:,:,:,:), allocatable :: a real, dimension(:,:,:), allocatable :: grid real, dimension(:), allocatable :: cc real :: t, sum integer :: j, k, i, ii, imid, ipos integer :: mx2, my2, mz2, mv2 integer :: l22, m22 ! ! NILS: For the time beeing the prerun simulation can not have been run in parallell ! NILS: Should fix this soon. ! ! Read dimension of the array stored during the 1D prerun ! processor = 'proc0' file = trim(prerun_dir)//'/data/'//trim(processor)//'/dim.dat' print*,'Reading prerun dim from ',file open (1,FILE=file) read (1,*) mx2, my2, mz2, mv2 close (1) l22 = mx2-l1+1 m22 = mx2-2*(l1-1) allocate(a(mx2,my2,mz2,mv2), grid(m22,my2-6,mz2-6), cc(m22)) ! ! Read the grid used during the prerun ! file = trim(prerun_dir)//'/data/'//trim(processor)//'/grid.dat' print*,'Reading prerun grid from ',file open (1,FILE=file,FORM='unformatted') read (1) t,grid(:,1,1),grid(1,:,1),grid(1,1,:) close (1) ! ! Read the data stored during the prerun ! file = trim(prerun_dir)//'/data/'//trim(processor)//'/var.dat' print*,'Reading inlet data from ',file open (1,FILE=file,FORM='unformatted') read (1) a close (1) ! ! Definition of a progress variable ! cc(1) = 0. ipos = 0 imid = 0 do ii = 2, m22-1 cc(ii) = (exp(a(ii,m1,n1,iuz+2)) - exp(a(l1,m1,n1,iuz+2)))/ & (exp(a(m22-1,m1,n1,iuz+2)) - exp(a(l1,m1,n1,iuz+2))) if (cc(ii) > 0.7 .and. cc(ii-1) <= 0.7) imid = ii if (grid(ii,1,1) > flame_pos .and. grid(ii-1,1,1) <= flame_pos) ipos = ii if (ipos > 0 .and. imid > 0) exit enddo ! ! The center of the flame, arbitrarily identified by cc=0.7, is ! located at the position flame_pos (chemistry_init_pars) ! The flame is thus moved in its domain ! ! Spread the data on the f-array ! do j = 1,my do k = 1,mz do i = 1,mx ! do ii = 2, m22-1 if (x(i) > grid(ii,1,1)-(grid(imid,1,1)-grid(ipos,1,1)) .and. x(i) & <= grid(ii+1,1,1)-(grid(imid,1,1)-grid(ipos,1,1))) then if (.not. lperi(1)) then f(i,j,k,iux) = f(i,j,k,iux)+a(ii,m1,n1,iux)+(x(i)-grid(ii,1,1) & +(grid(imid,1,1)-grid(ipos,1,1)))*(a(ii+1,m1,n1,iux) & -a(ii,m1,n1,iux))/(grid(ii+1,1,1)-grid(ii,1,1)) endif f(i,j,k,iuz+1:mvar) = a(ii,m1,n1,iuz+1:mvar)+(x(i)-grid(ii,1,1)+ & (grid(imid,1,1)-grid(ipos,1,1)))*(a(ii+1,m1,n1,iuz+1:mvar) & -a(ii,m1,n1,iuz+1:mvar))/(grid(ii+1,1,1)-grid(ii,1,1)) elseif (x(i) <= grid(2,1,1)-(grid(imid,1,1)-grid(ipos,1,1))) then if (.not. lperi(1)) f(i,j,k,iux) = f(i,j,k,iux)+a(l1,m1,n1,iux) f(i,j,k,iuz+1:mvar) = a(l1,m1,n1,iuz+1:mvar) exit elseif (x(i) >= grid(m22-1,1,1)-(grid(imid,1,1)-grid(ipos,1,1))) then if (.not. lperi(1)) f(i,j,k,iux) = f(i,j,k,iux)+a(l22,m1,n1,iux) f(i,j,k,iuz+1:mvar) = a(l22,m1,n1,iuz+1:mvar) exit endif enddo ! ! Renormalize the species mass fractions ! sum = 0. do ii = 1, nchemspec sum = sum + f(i,j,k,ichemspec(ii)) enddo f(i,j,k,ichemspec(1):ichemspec(nchemspec)) = & f(i,j,k,ichemspec(1):ichemspec(nchemspec))/sum ! enddo enddo enddo ! ! Set the y and z velocities to zero in order to avoid random noise ! ! if (.not. reinitialize_chemistry) then ! f(:,:,:,iuy:iuz)=0 ! endif ! deallocate(a,grid,cc) ! endsubroutine prerun_1D !*********************************************************************** subroutine prerun_1D_opp(f,prerun_dir) ! ! read snapshot file, possibly with mesh and time (if mode=1) ! 04-aug-10/julien: coded ! character(len=*) :: prerun_dir character(len=100) :: file character(len=10) :: processor real, dimension(mx,my,mz,mfarray) :: f real, dimension(:,:,:,:), allocatable :: a real, dimension(:,:,:), allocatable :: grid real :: t real :: x1, x2, xm, xfl integer :: j, k, i, ii, ifl integer :: mx2, my2, mz2, mv2 integer :: l22, m22 ! x1 = xyz0(1)+Lxyz(1) / 4. x2 = x1+Lxyz(1) / 2. xm = xyz0(1)+Lxyz(1) / 2. ! ! Read dimension of the array stored during the 1D prerun ! processor = 'proc0' file = trim(prerun_dir)//'/data/'//trim(processor)//'/dim.dat' print*,'Reading prerun dim from ',file open (1,FILE=file) read (1,*) mx2, my2, mz2, mv2 close (1) l22 = mx2-l1+1 m22 = mx2-2*(l1-1) allocate(a(mx2,my2,mz2,mv2), grid(m22,my2-6,mz2-6)) ! ! Read the grid used during the prerun ! file = trim(prerun_dir)//'/data/'//trim(processor)//'/grid.dat' print*,'Reading prerun grid from ',file open (1,FILE=file,FORM='unformatted') read (1) t,grid(:,1,1),grid(1,:,1),grid(1,1,:) close (1) ! ! Read the data stored during the prerun ! file = trim(prerun_dir)//'/data/'//trim(processor)//'/var.dat' print*,'Reading inlet data from ',file open (1,FILE=file,FORM='unformatted') read (1) a close (1) ! do ii = 1, m22-1 if (exp(a(ii,m1,n1,iuz+2)) < 1200. .and. exp(a(ii+1,m1,n1,iuz+2)) >= 1200.) then xfl = grid(ii,1,1) ifl = ii exit endif enddo ! ! Spread the data on the f-array : interpolations ! do j = 1,my do k = 1,mz do i = 1,mx ! if (x(i) < xm) then do ii = 2, m22-1 if (x(i)-x1 > grid(ii,1,1)-xfl .and. x(i)-x1 <= grid(ii+1,1,1)-xfl) then f(i,j,k,iuz+1:mvar) = a(ii,m1,n1,iuz+1:mvar)+((x(i)-x1)-(grid(ii,1,1)-xfl))* & (a(ii+1,m1,n1,iuz+1:mvar)-a(ii,m1,n1,iuz+1:mvar))/ & (grid(ii+1,1,1)-grid(ii,1,1)) exit elseif (x(i)-x1 <= grid(l1,1,1)-xfl) then f(i,j,k,iuz+1:mvar) = a(l1,m1,n1,iuz+1:mvar) exit elseif (x(i)-x1 >= grid(m22,1,1)-xfl) then f(i,j,k,iuz+1:mvar) = a(l22,m1,n1,iuz+1:mvar) exit endif enddo ! elseif (x(i) >= xm) then do ii = 2, m22-1 if (x(i)-x2 > grid(ii,1,1)-xfl .and. x(i)-x2 <= grid(ii+1,1,1)-xfl & .and. 2*ifl > ii+1) then f(i,j,k,iuz+1:mvar) = a(2*ifl-ii,m1,n1,iuz+1:mvar)+((x(i)-x2)-(grid(ii,1,1)-xfl)) & *(a(2*ifl-ii-1,m1,n1,iuz+1:mvar)-a(2*ifl-ii,m1,n1,iuz+1:mvar))/(grid(ii+1,1,1)-grid(ii,1,1)) exit elseif (x(i)-x2 <= grid(l1,1,1)-xfl) then f(i,j,k,iuz+1:mvar) = a(l22,m1,n1,iuz+1:mvar) exit elseif (x(i)-x2 >= grid(m22-1,1,1)-xfl .or. 2*ifl <= ii+1) then f(i,j,k,iuz+1:mvar) = a(l1,m1,n1,iuz+1:mvar) exit endif enddo endif enddo enddo enddo ! ! Set the y and z velocities to zero in order to avoid random noise ! ! if (.not. reinitialize_chemistry) then ! f(:,:,:,iuy:iuz)=0 ! endif ! deallocate(a,grid) ! endsubroutine prerun_1D_opp !*********************************************************************** subroutine FlameMaster_ini(f,file_name) ! ! read FlameMaster file for flame initialization ! 11-nov-10/julien: coded ! character(len=*) :: file_name real, dimension(mx,my,mz,mfarray) :: f real, dimension(:,:), allocatable :: a real, dimension(:), allocatable :: grid real, dimension(:), allocatable :: cc integer :: j, k, i, ii, imid, ipos integer :: is, js integer :: nsp, npts character(len=8), dimension(:), allocatable :: name_sp character(len=10) :: car10='nothing' character(len=12) :: car12 character(len=20) :: car20='nothing' character(len=1) :: car1 ! lFlameMaster = .true. ! ! Read dimension of the array stored in the FlameMaster initial file ! open (1,FILE=trim(file_name)) if (lroot) print*, 'Reading initial conditions in file ', trim(file_name) do while (car10 /= 'FlameThick') read (1,*) car10 enddo read (1,*) car12, car1, nsp read (1,*) car12, car1, npts allocate(a(npts,nsp+3), grid(npts), cc(npts)) allocate(name_sp(nsp)) ! ! Read the data stored during the prerun ! do while (car10 /= 'body') read (1,*) car10 enddo read (1,*) read (1,*) grid grid = grid*100. ! i = 0 do while (trim(car20) /= 'trailer') read (1,*) car20 select case (car20(1:12)) case ('massflowrate') read (1,*) a(:,1) case ('temperature') read (1,*) a(:,3) case ('massfraction') i = i+1 read (1,*) a(:,3+i) name_sp(i) = car20(14:20) case ('density') read (1,*) a(:,2) endselect enddo close (1) a(:,1) = a(:,1) / a(:,2) a(:,1) = a(:,1)*100. a(:,2) = a(:,2)/1000. ! ! Definition of a progress variable ! cc(1) = 0. ipos = 0 imid = 0 do ii = 2, npts cc(ii) = (a(ii,3) - a(1,3)) / (a(npts,3) - a(1,3)) if (cc(ii) > 0.7 .and. cc(ii-1) <= 0.7) imid = ii if (grid(ii) > flame_pos .and. grid(ii-1) <= flame_pos) ipos = ii if (ipos > 0 .and. imid > 0) exit enddo ! ! The center of the flame, arbitrarily identified by cc=0.7, is ! located at the position flame_pos (chemistry_init_pars) ! The flame is thus moved in its domain ! ! Spread the data on the f-array : interpolations ! do j = 1,my do k = 1,mz do i = 1,mx ! do ii = 2, npts-1 if (x(i) > grid(ii)-(grid(imid)-grid(ipos)) .and. x(i) & <= grid(ii+1)-(grid(imid)-grid(ipos))) then if (.not. init_from_file) & f(i,j,k,iux) = f(i,j,k,iux)+a(ii,1)+(x(i)-grid(ii)+(grid(imid)-grid(ipos)))* & (a(ii+1,1)-a(ii,1)) / (grid(ii+1)-grid(ii)) f(i,j,k,iuz+2) = a(ii,3)+(x(i)-grid(ii)+(grid(imid)-grid(ipos)))* & (a(ii+1,3)-a(ii,3)) / (grid(ii+1)-grid(ii)) f(i,j,k,iuz+1) = a(ii,2)+(x(i)-grid(ii)+(grid(imid)-grid(ipos)))* & (a(ii+1,2)-a(ii,2)) / (grid(ii+1)-grid(ii)) do is = 1, nsp do js = 1, nchemspec if (trim(name_sp(is)) == trim(varname(ichemspec(js)))) then f(i,j,k,iuz+2+js) = a(ii,is+3)+(x(i)-grid(ii)+(grid(imid)-grid(ipos)))* & (a(ii+1,is+3)-a(ii,is+3)) / (grid(ii+1)-grid(ii)) exit endif enddo enddo ! elseif (x(i) <= grid(2)-(grid(imid)-grid(ipos))) then if (.not. init_from_file) f(i,j,k,iux) = f(i,j,k,iux)+a(1,1) f(i,j,k,iuz+2) = a(1,3) f(i,j,k,iuz+1) = a(1,2) do is = 1, nsp do js = 1, nchemspec if (trim(name_sp(is)) == trim(varname(ichemspec(js)))) then f(i,j,k,iuz+2+js) = a(1,is+3) exit endif enddo enddo exit ! elseif (x(i) >= grid(npts)-(grid(imid)-grid(ipos))) then if (.not. init_from_file) f(i,j,k,iux) = f(i,j,k,iux)+a(npts,1) f(i,j,k,iuz+2) = a(npts,3) f(i,j,k,iuz+1) = a(npts,2) do is = 1, nsp do js = 1, nchemspec if (trim(name_sp(is)) == trim(varname(ichemspec(js)))) then f(i,j,k,iuz+2+js) = a(npts,is+3) exit endif enddo enddo exit endif enddo ! ! Renormalize the species mass fractions ! f(i,j,k,iuz+3:iuz+2+nchemspec) = f(i,j,k,iuz+3:iuz+2+nchemspec) & / sum(f(i,j,k,iuz+3:iuz+2+nchemspec)) ! enddo enddo enddo if (.not. ldensity_nolog) f(:,:,:,iuz+1) = log(f(:,:,:,iuz+1)) if (.not. ltemperature_nolog) f(:,:,:,iuz+2) = log(f(:,:,:,iuz+2)) ! ! Set the y and z velocities to zero in order to avoid random noise ! ! if (.not. reinitialize_chemistry) then ! f(:,:,:,iuy:iuz)=0 ! endif ! deallocate(a,grid,name_sp) ! endsubroutine FlameMaster_ini !*********************************************************************** subroutine get_reac_rate(wt,f,p) ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(nx) :: wt real, dimension(nx,nchemspec) :: ydot type (pencil_case) :: p ! ydot = p%DYDt_reac f(l1:l2,m,n,ireaci(1):ireaci(nchemspec)) = ydot ! f(l1:l2,m,n,ireaci(1):ireaci(nchemspec))+ydot ! if (maux == nchemspec+1) f(l1:l2,m,n,ireaci(nchemspec)+1) = wt ! f(l1:l2,m,n,ireaci(nchemspec)+1)+wt ! endsubroutine get_reac_rate !*********************************************************************** subroutine chemistry_clean_up ! if (allocated(Bin_diff_coef)) deallocate(Bin_diff_coef) if (allocated(stoichio)) deallocate(stoichio) if (allocated(Sijm)) deallocate(Sijm) if (allocated(Sijp)) deallocate(Sijp) if (allocated(kreactions_z)) deallocate(kreactions_z) if (allocated(kreactions_m)) deallocate(kreactions_m) if (allocated(kreactions_p)) deallocate(kreactions_p) if (allocated(reaction_name)) deallocate(reaction_name) if (allocated(B_n)) deallocate(B_n) if (allocated(alpha_n)) deallocate(alpha_n) if (allocated(E_an)) deallocate(E_an) if (allocated(low_coeff)) deallocate(low_coeff) if (allocated(high_coeff)) deallocate(high_coeff) if (allocated(troe_coeff)) deallocate(troe_coeff) if (allocated(a_k4)) deallocate(a_k4) if (allocated(Mplus_case)) deallocate(Mplus_case) if (allocated(photochem_case)) deallocate(photochem_case) if (allocated(net_react_m)) deallocate(net_react_m) if (allocated(net_react_p)) deallocate(net_react_p) if (allocated(back)) deallocate(back) if (allocated(Diff_full)) deallocate(Diff_full) if (allocated(Diff_full_add)) deallocate(Diff_full_add) ! endsubroutine chemistry_clean_up !*********************************************************************** subroutine find_remove_real_stoic(Speciesstring,lreal,stoi,startindex) ! character(len=*), intent(in) :: Speciesstring logical, intent(inout) :: lreal real, intent(inout) :: stoi integer, intent(inout) :: startindex integer :: lreadable ! read (Speciesstring(1:3),'(F3.1)',iostat=lreadable) stoi if (lreadable == 0) then startindex = startindex+2 lreal = .True. endif ! endsubroutine find_remove_real_stoic !!*********************************************************************** subroutine chemspec_normalization_N2(f) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: sum_Y !, sum_Y2 integer :: k ,isN2, ichemsN2 logical :: lsN2 ! call find_species_index('N2', isN2, ichemsN2, lsN2 ) sum_Y=0.0 !; sum_Y2=0.0 do k=1,nchemspec if (k/=ichemsN2) then sum_Y=sum_Y+f(:,:,:,ichemspec(k)) endif enddo f(:,:,:,isN2)=1.0-sum_Y ! endsubroutine chemspec_normalization_N2 !*********************************************************************** subroutine getmu_array(f,mu1_full) ! ! Calculate mean molecular weight ! ! 16-mar-10/natalia ! 30-jun-17/MR: moved here from eos_chemistry. ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: mu1_full integer :: k,j2,j3 ! ! Mean molecular weight ! mu1_full=0. do k=1,nchemspec if (species_constants(k,imass)>0.) then do j2=mm1,mm2 do j3=nn1,nn2 mu1_full(:,j2,j3)= & mu1_full(:,j2,j3)+unit_mass*f(:,j2,j3,ichemspec(k)) & /species_constants(k,imass) enddo enddo endif enddo ! endsubroutine getmu_array !*********************************************************************** subroutine read_Lewis ! ! Reading of the species Lewis numbers in an input file ! ! 21-jun-10/julien: coded ! 30-jun-17/MR: moved here from eos_chemistry. ! 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 read_transport_data ! ! Reading of the chemkin transport data ! ! 01-apr-08/natalia: coded ! 30-jun-17/MR: moved here from eos_chemistry. ! 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 !*********************************************************************** endmodule Chemistry