!$Id: particles_surfspec.f90 21950 2014-07-08 08:53:00Z jonas.kruger $
!
! MOUDLE_DOC: This module takes care the gas phase species in the
! MODULE_DOC: immediate vicinity of reactive particles.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lparticles_surfspec=.true.
!
!***************************************************************
module Particles_surfspec
!
use Cdata
use Cparam
use General, only: keep_compiler_quiet
use Messages
use Particles_cdata
use Particles_mpicomm
use Particles_sub
use Particles_chemistry
use Chemistry
use EquationOfState
use SharedVariables
!
implicit none
!
include 'particles_surfspec.h'
!
!*********************************************************************!
! Particle independent variables below here !
!*********************************************************************!
!
character(len=labellen), dimension(ninit) :: init_surf='nothing'
character(len=10), dimension(:), allocatable :: solid_species
real, dimension(10) :: init_surf_mol_frac ! INIT_DOC: Initial surface fraction
real, dimension(:,:), allocatable :: nu_power
integer, dimension(nchemspec) :: ispecaux=0
integer :: ispecenth
integer :: ndiffsteps=3 ! RUN_DOC: Number of mass transfer diffusion steps
integer :: Ysurf
logical :: lspecies_transfer=.true. ! RUN_DOC: species transfer between particle and gas
logical :: linfinite_diffusion=.true. ! RUN_DOC: infinitely fast diffusion between particle and gas
logical :: lboundary_explicit=.true. ! RUN_DOC: explicit evolution of particle surface species
logical :: lpchem_cdtc=.false.
logical :: lpchem_mass_enth=.true.
logical :: lpfilter =.true.
logical :: lwrite=.true.
real :: rdiffconsts=0.1178
logical :: ldiffuse_backspec=.false., ldiffs=.false.
logical :: ldiffuse_backenth=.false., ldiffenth=.false.
!
integer :: jH2, jO2, jCO2, jCO, jCH4, jN2, jH2O
integer :: jOH, jAR, jO, jH, jCH, jCH2, jHCO, jCH3
integer, dimension(:), allocatable :: idiag_surf
!
!*********************************************************************!
! Particle dependent variables below here !
!*********************************************************************!
!
real, dimension(:,:), allocatable :: X_infty_reactants
real, dimension(:,:), allocatable :: mass_trans_coeff_reactants
real, dimension(:,:), allocatable :: mass_trans_coeff_species
real, dimension(:,:,:), allocatable :: weight_array, dmass_frac_dt
!
namelist /particles_surf_init_pars/ init_surf, init_surf_mol_frac
!
namelist /particles_surf_run_pars/ &
lboundary_explicit, linfinite_diffusion, lspecies_transfer, &
lpchem_mass_enth, ldiffuse_backspec, ldiffs,rdiffconsts, &
ndiffsteps,ldiffuse_backenth,ldiffenth,lpfilter
!
integer :: idiag_dtpchem=0 ! DIAG_DOC: $dt_{particle,chemistry}$
!
contains
! ******************************************************************************
! This is a wrapper routine for particle dependent and particle
! independent variables
! JONAS: Back to standalone via mpar_loc=1?
!
subroutine register_particles_surfspec()
!
use FArrayManager, only: farray_register_auxiliary
!
!
character(len=11) :: chemspecaux
integer :: i
!
! if (lroot) call svn_id( &
! "$Id: particles_surfspec.f90 20849 2014-10-06 18:45:43Z jonas.kruger $")
!
call register_indep_psurfspec()
call register_dep_psurfspec()
!
! We need to register an auxiliary array to diffuse the species transfer
!
if (ldiffuse_backspec .and. lspecies_transfer) then
chemspecaux = 'ichemspec '
do i = 1,nchemspec
write (chemspecaux(10:11),'(I2)') i
call farray_register_auxiliary(chemspecaux,ispecaux(i),communicated=.true.)
enddo
elseif (ldiffuse_backspec .and. .not. lspecies_transfer) then
call fatal_error('particles_surfspec:', &
'diffusion of the species transfer needs lspecies_transfer')
endif
!
! Diffusion of mass bound enthalpy
!
if (ldiffuse_backenth .and. lspecies_transfer) then
call farray_register_auxiliary('ispecenth',ispecenth,communicated=.true.)
elseif (ldiffuse_backenth .and. .not. lspecies_transfer) then
call fatal_error('particles_surfspec:', &
'diffusion of the mass bound enthalpy needs lspecies_transfer')
endif
!
endsubroutine register_particles_surfspec
! ******************************************************************************
!
subroutine register_indep_psurfspec()
integer :: i, k, stat, k1, k2
character(len=10), dimension(40) :: reactants
!
if (nsurfreacspec /= N_surface_species) then
print*,'N_surface_species: ', N_surface_species
print*,'NSURFREACSPEC :', nsurfreacspec
call fatal_error('register_particles_surf', &
'wrong size of storage for surface species allocated')
endif
!
! Increase of npvar according to N_surface_species, which is
! the concentration of gas phase species at the particle surface
!
if (N_surface_species <= 1) then
call fatal_error('register_particles_', 'N_surface_species must be > 1')
endif
!
! Allocate memory for a number of arrays
allocate(dependent_reactant(N_surface_reactions),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for dependent_reactant')
allocate(nu(N_surface_species,N_surface_reactions),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for nu')
allocate(nu_prime(N_surface_species,N_surface_reactions),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for nu_prime')
allocate(ac(N_surface_species),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for ac')
allocate(jmap(N_surface_species),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for jmap')
allocate(solid_species(N_surface_species),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for solid_species')
allocate(nu_power(N_surface_species,N_surface_reactions),STAT=stat)
if (stat > 0) call fatal_error('register_indep_chem', &
'Could not allocate memory for nu_power')
allocate(idiag_surf(N_surface_species),STAT=stat)
if (stat > 0) call fatal_error('register_indep_chem', &
'Could not allocate memory for idiag_surf')
idiag_surf = 0
!
call create_ad_sol_lists(solid_species,'sol')
call get_reactants(reactants)
call sort_compounds(reactants,solid_species,n_surface_species)
!
do i = 1,N_surface_species
call append_npvar('i'//solid_species(i),isurf)
enddo
isurf_end = isurf
!
! create binding between gas phase and near field gas species
! chemistry and particles_chemistry module
!
call create_jmap()
!
! print*, jH2O,JCO2,jH2,jO2,jCO,jCH,jHCO,jCH2,jCH3
!
! Set number of carbon atoms for each surface species
!
call get_ac(ac,solid_species,N_surface_species)
!
! Set the stoichiometric matrixes
! print*, 'Number of surface species: ',N_surface_species
!
call create_stoc(solid_species,nu,.true., N_surface_species,nu_power)
call create_stoc(solid_species,nu_prime,.false., N_surface_species)
!
! Define which gas phase reactants the given reaction depends on
call create_dependency(nu,dependent_reactant, &
n_surface_reactions,n_surface_reactants)
!
! Find the mole production of the forward reaction
call create_dngas()
!
if (lwrite) then
call write_outputfile()
endif
lwrite = .false.
!
endsubroutine register_indep_psurfspec
! ******************************************************************************
subroutine register_dep_psurfspec()
endsubroutine register_dep_psurfspec
! ******************************************************************************
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 29-sep-14/jonas coded
!
subroutine initialize_particles_surf(f)
real, dimension(mx,my,mz,mfarray) :: f
integer :: dimx, dimy, dimz, ncells=1
!
! print*, weight_array
!
if (lparticlemesh_gab) ncells = 7
if (lparticlemesh_tsc) ncells = 3
if (lparticlemesh_cic) ncells = 2
!
if (nxgrid /= 1) then
dimx = ncells
else
dimx = 1
endif
if (nygrid /= 1) then
dimy = ncells
else
dimy = 1
endif
if (nzgrid /= 1) then
dimz = ncells
else
dimz = 1
endif
!
if (allocated(weight_array)) deallocate(weight_array)
if (lparticlemesh_gab .or. lparticlemesh_tsc .or. lparticlemesh_cic) then
allocate(weight_array(dimx,dimy,dimz))
endif
if (.not. allocated(weight_array)) then
allocate(weight_array(1,1,1))
endif
call precalc_weights(weight_array)
!
call keep_compiler_quiet(f)
endsubroutine initialize_particles_surf
! ******************************************************************************
subroutine read_particles_surf_init_pars(iostat)
use File_io, only: parallel_unit
integer, intent(out) :: iostat
!
read (parallel_unit, NML=particles_surf_init_pars, IOSTAT=iostat)
endsubroutine read_particles_surf_init_pars
! ******************************************************************************
subroutine write_particles_surf_init_pars(unit)
integer, intent(in) :: unit
!
write (unit, NML=particles_surf_init_pars)
endsubroutine write_particles_surf_init_pars
! ******************************************************************************
subroutine read_particles_surf_run_pars(iostat)
use File_io, only: parallel_unit
integer, intent(out) :: iostat
!
read (parallel_unit, NML=particles_surf_run_pars, IOSTAT=iostat)
endsubroutine read_particles_surf_run_pars
! ******************************************************************************
subroutine write_particles_surf_run_pars(unit)
integer, intent(in) :: unit
!
write (unit, NML=particles_surf_run_pars)
endsubroutine write_particles_surf_run_pars
! ******************************************************************************
! Initial particle surface fractions
!
! 01-sep-14/jonas: coded
!
subroutine init_particles_surf(f,fp,ineargrid)
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mpar_loc,mparray) :: fp
integer, dimension(mpar_loc,3) :: ineargrid
real :: sum_surf_spec
real :: mean_molar_mass
integer :: i, j, k, ix0, igas, iy0, iz0
!
intent(in) :: f
intent(out) :: fp
!
fp(:,isurf:isurf_end) = 0.
do j = 1,ninit
!
! Writing the initial surface species fractions
select case (init_surf(j))
!
case ('nothing')
if (lroot .and. j == 1) print*, 'init_particles_surf,gas phase: nothing'
!
case ('constant')
if (lroot) print*, 'init_particles_surf: Initial Surface Fractions'
!
! This ensures that we don't have unphysical values as init
sum_surf_spec = sum(init_surf_mol_frac(1:N_surface_species))
if (sum_surf_spec > 1) then
if (lroot) print*, 'Sum of all surface fractions >1, normalizing...'
init_surf_mol_frac(1:N_surface_species) = &
init_surf_mol_frac(1:N_surface_species) / sum_surf_spec
endif
!
if (sum_surf_spec == 0.0) then
call fatal_error('particles_surfspec', &
'initial molefraction was given without value. '// &
'please specify the initial gas fraction')
endif
!
do k = 1,mpar_loc
fp(k,isurf:isurf_end) = init_surf_mol_frac(1:N_surface_species)
enddo
case ('gas')
!
! The starting particle surface mole fraction is equal to the
! gas phase composition at the particles position.
! This is not functional, and this would need the initialization of
! The gas field before
do k = 1,mpar_loc
mean_molar_mass = 0.0
!
do i = 1,nchemspec
mean_molar_mass = mean_molar_mass + &
species_constants(i,imass) * f(4,4,4,ichemspec(i))
enddo
!
!
do i = 1, N_surface_species
igas = ichemspec(jmap(i))
fp(k,isurf+i-1) = f(4,4,4,igas) / &
species_constants(jmap(i),imass)*mean_molar_mass
enddo
enddo
!
case default
if (lroot) &
print*, 'init_particles_ads: No such such value for init_surf: ', &
trim(init_surf(j))
call fatal_error('init_particles_surf','')
endselect
enddo
endsubroutine init_particles_surf
! ******************************************************************************
! evolution of particle surface fractions
! (all particles on one node)
!
! 1-oct-14/Jonas: coded
!
subroutine dpsurf_dt(f,df,fp,dfp,ineargrid)
!
use Boundcond
use Mpicomm
!
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,my,mz,mvar) :: df
real, dimension(mpar_loc,mparray), intent(in) :: fp
real, dimension(mpar_loc,mpvar) :: dfp
real :: dt1
integer, dimension(mpar_loc,3) :: ineargrid
integer :: i, j,g
integer :: li,mi,ni
!
! equations.tex eq 37
! call keep_compiler_quiet(fp)
call keep_compiler_quiet(dfp)
call keep_compiler_quiet(ineargrid)
dt1 = 1./dt
!
! Diffuse the transfer of species from particle to fluid, and adapt
! the bulk species (ichemspec(nchemspec)) to ensure species conservation
!
if (lspecies_transfer .and. ldiffuse_backspec) then
if (ldensity_nolog) call fatal_error('particles_surf', &
'not implemented for ldensity_nolog')
do j = 1,nchemspec-1
do i = 1,ndiffsteps
call boundconds_x(f,ispecaux(j),ispecaux(j))
call initiate_isendrcv_bdry(f,ispecaux(j),ispecaux(j))
call finalize_isendrcv_bdry(f,ispecaux(j),ispecaux(j))
call boundconds_y(f,ispecaux(j),ispecaux(j))
call boundconds_z(f,ispecaux(j),ispecaux(j))
!
call diffuse_interaction(f(:,:,:,ispecaux(j)),ldiffs,.False.,rdiffconsts)
enddo
!
! Control that we don't influence points that would bring the mass fraction into the negative
!
if (.not. lpfilter) then
df(l1:l2,m1:m2,n1:n2,ichemspec(j)) = df(l1:l2,m1:m2,n1:n2,ichemspec(j)) + &
f(l1:l2,m1:m2,n1:n2,ispecaux(j))
else
where ((f(l1:l2,m1:m2,n1:n2,ichemspec(j))+ &
f(l1:l2,m1:m2,n1:n2,ispecaux(j))*dt)< 0.0)
f(l1:l2,m1:m2,n1:n2,ispecaux(j)) = f(l1:l2,m1:m2,n1:n2,ichemspec(j))*dt1*(-0.9)
where ((f(l1:l2,m1:m2,n1:n2,ichemspec(j)) <= 1E-25))
f(l1:l2,m1:m2,n1:n2,ispecaux(j)) = 0.0
end where
end where
df(l1:l2,m1:m2,n1:n2,ichemspec(j)) = df(l1:l2,m1:m2,n1:n2,ichemspec(j)) + &
f(l1:l2,m1:m2,n1:n2,ispecaux(j))
endif
enddo
df(l1:l2,m1:m2,n1:n2,ichemspec(nchemspec)) = &
df(l1:l2,m1:m2,n1:n2,ichemspec(nchemspec))- &
sum(df(l1:l2,m1:m2,n1:n2,ichemspec(:)),DIM=4)
endif
!
! Diffusion of the mass bound enthalpy
!
if (lspecies_transfer .and. ldiffuse_backenth) then
if (ldensity_nolog .and. ltemperature_nolog) &
call fatal_error('particles_surf', 'diffusion of mass bound enthalpy &
& not implemented for ldensity_nolog')
do i = 1,ndiffsteps
call boundconds_x(f,ispecenth,ispecenth)
call initiate_isendrcv_bdry(f,ispecenth,ispecenth)
call finalize_isendrcv_bdry(f,ispecenth,ispecenth)
call boundconds_y(f,ispecenth,ispecenth)
call boundconds_z(f,ispecenth,ispecenth)
call diffuse_interaction(f(:,:,:,ispecenth),ldiffenth,.False.,rdiffconsts)
!
enddo
df(l1:l2,m1:m2,n1:n2,ilnTT) = df(l1:l2,m1:m2,n1:n2,ilnTT) + &
f(l1:l2,m1:m2,n1:n2,ispecenth)
endif
!
if (ldiagnos) then
do i = 1,N_surface_species
if (idiag_surf(i) /= 0) then
call sum_par_name(fp(1:npar_loc,isurf+i-1),idiag_surf(i))
endif
enddo
endif
!
endsubroutine dpsurf_dt
! ******************************************************************************
subroutine dpsurf_dt_pencil(f,df,fp,dfp,p,ineargrid)
!
! Evolution of particle surface fraction.
! (all particles on one pencil)
!
! 23-sep-14/Nils: coded
!
use Diagnostics
!
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,my,mz,mvar) :: df
real, dimension(mpar_loc,mparray) :: fp
real, dimension(mpar_loc,mpvar) :: dfp
real, dimension(nx,nchemspec) :: chem_reac
real, dimension(:,:), allocatable :: term, ndot
real, dimension(:), allocatable :: Cg_surf, mass_loss
real :: porosity
type (pencil_case) :: p
integer, dimension(mpar_loc,3) :: ineargrid
integer :: k, k1, k2, i, ix0, iy0, iz0
real :: weight, volume_cell, rho1_point
real :: mean_molar_mass, dmass, A_p, denth
real :: reac_pchem_weight, max_reac_pchem
real :: summan, m1_cell
real :: dmass_ndot
integer :: ix1, iy1, iz1, ierr
integer :: ixx, iyy, izz
integer :: ixx0, iyy0, izz0
integer :: ixx1, iyy1, izz1
integer :: index1, index2
real, pointer :: true_density_carbon
integer :: density_index
real :: dmass_frac_dt=0.0
real :: diffusion_transfer=0.0
!
intent(in) :: ineargrid
intent(inout) :: dfp, df, fp,f
!
call keep_compiler_quiet(df)
call keep_compiler_quiet(p)
call keep_compiler_quiet(ineargrid)
!
! Set the particle reaction time to 0 if no particles are present
!
reac_pchem = 0.0
!
if (lpreactions) then
! initializing the auxiliary pencils for the mass and mass bound enthalpy diffusi
!
volume_cell = (lxyz(1)*lxyz(2)*lxyz(3))/(nxgrid*nygrid*nzgrid)
if (ldiffuse_backenth) f(l1:l2,m,n,ispecenth) = 0.0
if (ldiffuse_backspec) then
do i = 1,nchemspec
f(l1:l2,m,n,ispecaux(i)) = 0.0
enddo
endif
!
! Do only if particles are present on the current pencil
!
if (npar_imn(imn) /= 0) then
!
k1 = k1_imn(imn)
k2 = k2_imn(imn)
!
allocate(Cg_surf(k1:k2))
allocate(mass_loss(k1:k2))
allocate(ndot(k1:k2,N_surface_species))
allocate(term(k1:k2,N_surface_reactants))
!
call get_shared_variable('true_density_carbon',true_density_carbon,ierr)
!
call calc_mass_trans_reactants()
call get_surface_chemistry(Cg_surf,ndot,mass_loss)
!
! set surface gas species composition
! (infinite diffusion, set as far field and convert from
! mass fractions to mole fractions)
!
do k = k1,k2
!
! particle on pencil loop
!
porosity = 1.0 - (fp(k,imp)/(fp(k,iap)**3*4./3.*pi))/true_density_carbon
A_p = 4.*pi*(fp(k,iap)**2)
mean_molar_mass = (interp_rho(k)*Rgas*interp_TT(k)/ interp_pp(k))
!
if (lboundary_explicit) then
if (lparticles_adsorbed) print*, 'values in surfspec begin',fp(k,isurf:isurf_end)
do i = 1,N_surface_reactants
diffusion_transfer = mass_trans_coeff_reactants(k,i) * (interp_species(k,jmap(i)) / &
species_constants(jmap(i),imass) * mean_molar_mass-fp(k,isurf+i-1))
term(k,i) = ndot(k,i) - fp(k,isurf+i-1) * sum(ndot(k,:)) + diffusion_transfer
if (lparticles_adsorbed) then
print*, '---------------------------'
print*, ' mass_trans_', mass_trans_coeff_reactants(k,i)
print*, ' ndot(k,i) ', ndot(k,i)
print*, 'fp(k,isurf+)',fp(k,isurf+i-1)
print*, 'sum(ndot(k,:',sum(ndot(k,:))
print*, 'x_spec ',interp_species(k,jmap(i)) / species_constants(jmap(i),imass) * &
mean_molar_mass
print*, '---------------------------'
endif
!
! the term 3/fp(k,iap) is ratio of the surface of a sphere to its volume
!
dfp(k,isurf+i-1) = dfp(k,isurf+i-1) + 3*term(k,i)/(porosity*Cg_surf(k)*fp(k,iap))
if (lparticles_adsorbed) then
print*, 'term(k,i)',term(k,i)
print*, 'dfp(k,isurf+i-1): ',k,i,dfp(k,isurf+i-1)
endif
enddo
! if (lparticles_adsorbed) stop
else
if (linfinite_diffusion .or. lbaum_and_street) then
do i = 1,N_surface_reactants
fp(k,isurf+i-1) = interp_species(k,jmap(i)) / &
species_constants(jmap(i),imass) * mean_molar_mass
dfp(k,isurf+i-1) = 0.
enddo
else
print*,'Must set linfinite_diffusion=T if lboundary_explicit=F.'
call fatal_error('dpsurf_dt_pencil', &
'Implicit solver for surface consentrations is not implemented.')
endif
endif
!
! the following block is thoroughly commented in particles_temperature
! find values for transfer of variables from particle to fluid
if (lspecies_transfer) then
ix0 = ineargrid(k,1)
iy0 = ineargrid(k,2)
iz0 = ineargrid(k,3)
call find_interpolation_indeces(ixx0,ixx1,iyy0,iyy1,izz0,izz1, &
fp,k,ix0,iy0,iz0)
!
! positive dmass means particle is losing mass
! jmap: gives the ichemspec of a surface specie
dmass = 0.
denth = 0.
!
do i = 1,N_surface_species
index1 = jmap(i)
dmass = dmass-ndot(k,i)*species_constants(index1,imass)*A_p
if (lpchem_mass_enth) then
if (ndot(k,i) > 0.0) then
!
! Species enthalpy at the particle phase temperature
! Qenth = ndot(mol/s/m2)*A(m2)*enth(J/mol)=J/s
! to convert to erg, the value has to be multiplied with 1e7 for values
! fetched from particles_chemistry
!
! Enthalpy from Pencil code at particles temperature
!
if (lpencil(i_H0_RT)) then
denth = denth+ndot(k,i)*A_p*p%cv(ix0-nghost)*(fp(k,iTp)-298.15)
else
call fatal_error('particles_surfspec','mass bound enthalpy transfer needs p%H0_RT')
endif
else
!
! Species enthalpy at the gas phase temperature
!
if (lpencil(i_H0_RT)) then
denth = denth+ndot(k,i)*A_p*p%cv(ix0-nghost)*(p%TT(ix0-nghost)-298.15)
else
call fatal_error('particles_surfspec', &
'mass bound enthalpy transfer needs p%H0_RT')
endif
endif
endif
enddo
!
! Prepare the max_reac_pchem
!
if (lfirst .and. ldt) max_reac_pchem = 0.0
!
! Sum for testing
!
summan = 0.0
!
! The whole following with find_index is a workaround
! to enable the accessing of nonreacting gas phase species by this routine
! since they are as well affected when the particle is adding species to the
! gas phase.
! find_index maps the i from the gas phase chemistry to the surface_species
! from particles_chemistry which is needed to access ndot, since this
! only contains surface species
!
if (ldensity_nolog) then
if (.not. ldiffuse_backspec) then
do i = 1,nchemspec-1
reac_pchem_weight = 0.0
if (find_index(i,jmap,N_surface_species) > 0 ) then
!NILS: Isn't index1=i?
! index1 = jmap(find_index(i,jmap,N_surface_species))
index2 = ichemspec(i)
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2) = &
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2)+(A_p* &
ndot(k,find_index(i,jmap,N_surface_species))* &
species_constants(i,imass)+ &
dmass*interp_species(k,i))*weight_array/ &
(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,irho)*volume_cell)
if (lfirst .and. ldt) then
reac_pchem_weight = max(reac_pchem_weight, &
abs(maxval(df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2))/ &
max(maxval(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2)),1e-10)))
endif
else
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(i)) = &
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(i)) + dmass * &
interp_species(k,i)*weight_array/ &
(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,irho)*volume_cell)
endif
enddo
else
call fatal_error('particles_surf','backdiffusion not yet implemented &
& for ldensity_nolog')
endif
else
!
! For density_log, there is the possibility to diffuse the species influence
! this happens in the else block of the next if clause
!
if (.not. ldiffuse_backspec) then
do i = 1,nchemspec-1
reac_pchem_weight = 0.0
if (find_index(i,jmap,N_surface_species) > 0 ) then
!NILS: Isn't index1=i?
! index1 = jmap(find_index(i,jmap,N_surface_species))
index2 = ichemspec(i)
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2) = &
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2)+(A_p* &
ndot(k,find_index(i,jmap,N_surface_species))* &
species_constants(i,imass)+ &
dmass*interp_species(k,i))*weight_array/ &
(exp(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnrho))*volume_cell)
if (lfirst .and. ldt) then
reac_pchem_weight = max(reac_pchem_weight, &
abs(max(maxval(df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2)/ &
f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,index2)),1e-10)))
endif
else
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(i)) = &
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(i)) + &
dmass*interp_species(k,i)*weight_array/ &
(exp(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnrho))*volume_cell)
endif
! summan = summan+dmass_frac_dt
enddo
else
do i = 1,nchemspec-1
if (find_index(i,jmap,N_surface_species) > 0 ) then
! NILS: Isn't index1=i?
! index1 = jmap(find_index(i,jmap,N_surface_species))
f(ix0,iy0,iz0,ispecaux(i)) = &
f(ix0,iy0,iz0,ispecaux(i))+(A_p* &
ndot(k,find_index(i,jmap,N_surface_species))* &
species_constants(i,imass)+ &
dmass*interp_species(k,i))/ &
(exp(f(ix0,iy0,iz0,ilnrho))*volume_cell)
else
f(ix0,iy0,iz0,ispecaux(i)) = &
f(ix0,iy0,iz0,ispecaux(i)) + &
dmass*interp_species(k,i)/ &
(exp(f(ix0,iy0,iz0,ilnrho))*volume_cell)
endif
enddo
endif
endif
!
! Solving for all but the other values, setting the last one to the
! negative values of all others.
!
if (.not. ldiffuse_backspec) then
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(nchemspec)) = &
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(nchemspec))- &
sum(df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ichemspec(:)),DIM=4)
endif
!
! Compare the current maximum reaction rate to the previous one
!
if (lfirst .and. ldt) max_reac_pchem = &
max(max_reac_pchem, reac_pchem_weight)
!
! Enthalpy transfer via mass transfer!
!
if (lpchem_mass_enth .and. .not. ldiffuse_backenth) then
if (ldensity_nolog) then
if (ltemperature_nolog) then
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,iTT) = df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,iTT) &
+denth*p%cv1(ix0-nghost)*weight_array/ &
(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,irho)*volume_cell)
else
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnTT) = df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnTT) &
+denth*p%cv1(ix0-nghost)*p%TT1(ix0-nghost)*weight_array/ &
(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,irho)*volume_cell)
endif
else
if (ltemperature_nolog) then
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,iTT) = df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,iTT) &
+denth*p%cv1(ix0-nghost)*weight_array/ &
(exp(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnrho))*volume_cell)
else
df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnTT) = df(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnTT) &
+denth*p%cv1(ix0-nghost)*p%TT1(ix0-nghost)*weight_array/ &
(exp(f(ixx0:ixx1,iyy0:iyy1,izz0:izz1,ilnrho))*volume_cell)
endif
endif
!
! Diffusion of mass bound enthalpy transfer
!
elseif (lpchem_mass_enth .and. ldiffuse_backenth) then
if (ldensity_nolog .and. ltemperature_nolog) then
call fatal_error('particles_surf','diffusion not implemented for &
& ldensity_nolog')
elseif (ldensity_nolog .and. .not. ltemperature_nolog) then
call fatal_error('particles_surf','diffusion not implemented for &
& ldensity_nolog')
elseif (.not. ldensity_nolog .and. ltemperature_nolog) then
call fatal_error('particles_surf','diffusion not implemented for &
& ltemperature_nolog')
else
f(ix0,iy0,iz0,ispecenth) = f(ix0,iy0,iz0,ispecenth) &
+denth*p%cv1(ix0-nghost)*p%TT1(ix0-nghost) / &
(exp(f(ix0,iy0,iz0,ilnrho))*volume_cell)
endif
endif
!
! Compare the current maximum reaction rate to the current maximum
! reaction rate in the current pencil
!
if (lfirst .and. ldt) reac_pchem = max(reac_pchem,max_reac_pchem)
!
endif
!
! end of particle on pencil loop
!
! if (lparticles_adsorbed) print*, 'values in surfspec end',fp(k,isurf:isurf_end)
enddo
!
!
if (ldiagnos) then
if (idiag_dtpchem /= 0 ) call max_name(reac_pchem/cdtc,idiag_dtpchem,l_dt=.true.)
endif
!
if (allocated(term)) deallocate(term)
if (allocated(ndot)) deallocate(ndot)
if (allocated(Cg_surf)) deallocate(Cg_surf)
if (allocated(mass_loss)) deallocate(mass_loss)
endif
!
endif
!
endsubroutine dpsurf_dt_pencil
! ******************************************************************************
! Read and register print parameters relevant for
! particles near field gas composition
!
! 06-oct-14/jonas: adapted
!
subroutine rprint_particles_surf(lreset,lwrite)
use Diagnostics
!
logical :: lreset
logical, optional :: lwrite
!
logical :: lwr
integer :: iname,i
character(len=7) :: diagn_surf, number
!
! Write information to index.pro
lwr = .false.
if (present(lwrite)) lwr = lwrite
if (lreset) then
idiag_dtpchem = 0
idiag_surf = 0
endif
!
if (lroot .and. ip < 14) print*,'rprint_particles_surf: run through parse list'
!
do iname = 1,nname
do i = 1,N_surface_species
write (number,'(I2)') i
diagn_surf = 'Ysurf'//trim(adjustl(number))
call parse_name(iname,cname(iname),cform(iname),trim(diagn_surf),idiag_surf(i))
enddo
call parse_name(iname,cname(iname),cform(iname),'dtpchem',idiag_dtpchem)
enddo
!
endsubroutine rprint_particles_surf
! ******************************************************************************
! 07-oct-2014/jonas: coded
! 29-oct-2014/jonas: moved parts of the code into routine, implemented
! error when gas phase species is in
! particles_chemistry, but not in chemistry
!
! find the indexes used in chemistry.f90 of species present
! in the near field of the particle
!
subroutine create_jmap()
use EquationOfState
!
integer :: index_glob, index_chem
logical :: found_species=.false.
!
inuH2O = find_species('H2O',solid_species,n_surface_species)
inuCO2 = find_species('CO2',solid_species,n_surface_species)
inuH2 = find_species('H2',solid_species,n_surface_species)
inuO2 = find_species('O2',solid_species,n_surface_species)
inuCO = find_species('CO',solid_species,n_surface_species)
inuCH = find_species('CH',solid_species,n_surface_species)
inuHCO = find_species('HCO',solid_species,n_surface_species)
inuCH2 = find_species('CH2',solid_species,n_surface_species)
inuCH3 = find_species('CH3',solid_species,n_surface_species)
!
call find_species_index('H2',index_glob,index_chem,found_species)
if (found_species) then
jH2 = index_chem
else
if (inuH2 > 0) call fatal_error('create_jmap','no H2 found')
endif
!
call find_species_index('O2',index_glob,index_chem,found_species)
if (found_species) then
jO2 = index_chem
else
if (inuO2 > 0) call fatal_error('create_jmap','no O2 found')
endif
!
call find_species_index('CO2',index_glob,index_chem,found_species)
if (found_species) then
jCO2 = index_chem
else
if (inuCO2 > 0) call fatal_error('create_jmap','no CO2 found')
endif
!
call find_species_index('CO',index_glob,index_chem,found_species)
if (found_species) then
jCO = index_chem
else
if (inuCO > 0) call fatal_error('create_jmap','no CO found')
endif
!
call find_species_index('CH4',index_glob,index_chem,found_species)
if (found_species) then
jCH4 = index_chem
else
if (inuCH4 > 0) call fatal_error('create_jmap','no CH4 found')
endif
!
call find_species_index('H2O',index_glob,index_chem,found_species)
if (found_species) then
jH2O = index_chem
else
if (inuH2O > 0) call fatal_error('create_jmap','no H2O found')
endif
!
call find_species_index('CH3',index_glob,index_chem,found_species)
if (found_species) then
jCH3 = index_chem
else
if (inuCH3 > 0) call fatal_error('create_jmap','no CH3 found')
endif
!
call find_species_index('CH',index_glob,index_chem,found_species)
if (found_species) then
jCH = index_chem
else
if (inuCH > 0) call fatal_error('create_jmap','no CH found')
endif
!
call find_species_index('CH2',index_glob,index_chem,found_species)
if (found_species) then
jCH2 = index_chem
else
if (inuCH2 > 0) call fatal_error('create_jmap','no CH2 found')
endif
! call find_species_index('HCO',index_glob,index_chem,found_species)
if (found_species) then
jHCO = index_chem
else
if (inuHCO > 0) call fatal_error('create_jmap','no HCO found')
endif
!
if (inuH2O > 0) jmap(inuH2O) = jH2O
if (inuCO2 > 0) jmap(inuCO2) = jCO2
if (inuH2 > 0) jmap(inuH2) = jH2
if (inuO2 > 0) jmap(inuO2) = jO2
if (inuCO > 0) jmap(inuCO) = jCO
if (inuCH > 0) jmap(inuCH) = jCH
if (inuHCO > 0) jmap(inuHCO) = jHCO
if (inuCH2 > 0) jmap(inuCH2) = jCH2
if (inuCH3 > 0) jmap(inuCH3) = jCH3
endsubroutine create_jmap
! ******************************************************************************
! restrict mass trans coefficients to surface reactants only
!
subroutine calc_mass_trans_reactants()
integer :: k, i, k1, k2
!
k1 = k1_imn(imn)
k2 = k2_imn(imn)
!
do k = k1,k2
do i = 1, N_surface_reactants
mass_trans_coeff_reactants(k,i) = mass_trans_coeff_species(k,i)
enddo
enddo
endsubroutine calc_mass_trans_reactants
! ******************************************************************************
! Allocate pencils for mass transport coefficiens
!
! nov-14/jonas: coded
!
subroutine allocate_surface_pencils()
integer :: k1, k2, stat
!
k1 = k1_imn(imn)
k2 = k2_imn(imn)
!
allocate(mass_trans_coeff_species(k1:k2,N_species), STAT=stat)
if (stat > 0) call fatal_error('allocate_variable_pencils', &
'Could not allocate memory for mass_trans_coeff_species')
allocate(mass_trans_coeff_reactants(k1:k2,N_surface_reactants),STAT=stat)
if (stat > 0) call fatal_error('register_indep_psurfchem', &
'Could not allocate memory for mass_trans_coeff_reactants')
!
endsubroutine allocate_surface_pencils
! ******************************************************************************
!
! Clean up surface pencils
!
! nov-14/jonas: coded
!
subroutine cleanup_surf_pencils()
!
deallocate(mass_trans_coeff_species)
deallocate(mass_trans_coeff_reactants)
!
endsubroutine cleanup_surf_pencils
! ******************************************************************************
!
! Calculate pencils used for surface gas phase fraction evolution
!
! nov-14/jonas: coded
!
subroutine calc_psurf_pencils(f,fp,p,ineargrid)
real, dimension(mpar_loc,mparray), intent(in) :: fp
real, dimension(mx,my,mz,mfarray), intent(in) :: f
integer, dimension(mpar_loc,3), intent(in) :: ineargrid
type (pencil_case) :: p
!
call allocate_surface_pencils()
call calc_mass_trans_coeff(f,fp,p,ineargrid)
!
endsubroutine calc_psurf_pencils
! ******************************************************************************
! diffusion coefficients of gas species at nearest grid point
!
! oct-14/Jonas: coded
!
subroutine calc_mass_trans_coeff(f,fp,p,ineargrid)
real, dimension(mx,my,mz,mparray), intent(in) :: f
real, dimension(mpar_loc,mparray), intent(in) :: fp
type (pencil_case) :: p
integer, dimension(mpar_loc,3), intent(in) :: ineargrid
integer :: k, k1, k2,i
integer :: ix0
integer :: spec_glob, spec_chem
real :: diff_coeff_species
real, dimension(:), allocatable :: Cg_diff
!
if (npar_imn(imn) /= 0) then
k1 = k1_imn(imn)
k2 = k2_imn(imn)
!
! allocate(diff_coeff_species(k1:k2,N_species))
allocate(Cg_diff(k1:k2))
!
! do k = k1,k2
! ix0 = ineargrid(k,1)
! do i = 1,N_surface_species
! diff_coeff_species(k,i) = p%Diff_penc_add(ix0-nghost,jmap(i))
! enddo
! enddo
!
do k = k1,k2
Cg_diff(k) = k
enddo
call get_surface_chemistry(Cg_diff(:))
do k = k1,k2
ix0 = ineargrid(k,1)
do i = 1,N_surface_species
diff_coeff_species = p%Diff_penc_add(ix0-nghost,jmap(i))
mass_trans_coeff_species(k,i) = Cg_diff(k)*diff_coeff_species/ fp(k,iap)
enddo
enddo
!
if (allocated(Cg_diff)) deallocate(Cg_diff)
endif
!
endsubroutine calc_mass_trans_coeff
! ******************************************************************************
subroutine particles_surfspec_clean_up()
!
if (allocated(dependent_reactant)) deallocate(dependent_reactant)
if (allocated(nu)) deallocate(nu)
if (allocated(nu_prime)) deallocate(nu_prime)
if (allocated(ac)) deallocate(ac)
if (allocated(jmap)) deallocate(jmap)
if (allocated(solid_species)) deallocate(solid_species)
if (allocated(nu_power)) deallocate(nu_power)
if (allocated(idiag_surf)) deallocate(idiag_surf)
!
endsubroutine particles_surfspec_clean_up
! ******************************************************************************
function find_index(element, list,lengthlist)
!
implicit none
!
integer :: find_index
integer :: element,i
integer, dimension(:) :: list
integer :: lengthlist
!
find_index = 0
!
do i = 1,lengthlist
if (element == list(i)) find_index = i
enddo
!
endfunction find_index
!***********************************************************************
subroutine write_outputfile()
!
! Write particle chemistry info to ./data/particle_chemistry.out
!
integer :: file_id=123
integer :: k
character(len=20) :: writeformat
character(len=30) :: output='./data/particle_chemistry.out'
!
open (file_id,file=output,position='append')
!
writeformat = '( A8)'
write (writeformat(2:3),'(I2)') N_surface_species
write (file_id,*) 'Gas phase surface species'
write (file_id,writeformat) solid_species
write (file_id,*) ''
!
writeformat = '(I2, F5.2)'
write (writeformat(5:6),'(I2)') N_surface_species
write (file_id,*) 'Nu'
do k = 1,N_surface_reactions
write (file_id,writeformat) k,nu(:,k)
enddo
write (file_id,*) ''
!
write (file_id,*) 'Nu*'
do k = 1,N_surface_reactions
write (file_id,writeformat) k,nu_prime(:,k)
enddo
write (file_id,*) ''
!
writeformat = '( I3)'
write (writeformat(2:3),'(I2)') N_surface_species
write (file_id,*) 'J_map, mapping to nchemspec'
write (file_id,writeformat) jmap
write (file_id,*) ''
!
endsubroutine write_outputfile
!***********************************************************************
endmodule Particles_surfspec