! $Id$ ! ! This module takes care of everything related to particle radius. ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! MPVAR CONTRIBUTION 1 ! MPAUX CONTRIBUTION 1 ! CPARAM logical, parameter :: lparticles_radius=.true. ! !*************************************************************** module Particles_radius ! use Cdata use General, only: keep_compiler_quiet use Messages use Particles_cdata use Particles_sub use Particles_chemistry ! implicit none ! include '../particles_radius.h' ! real :: vthresh_sweepup=-1.0, deltavp12_floor=0.0 real, dimension (ninit) :: ap0=0.0 real, dimension (ninit) :: radii_distribution=0.0 real :: tstart_sweepup_par=0.0, cdtps=0.2, cdtpc=0.2 real :: tstart_condensation_par=0.0 real :: apmin=0.0, latent_heat_SI=2.257e6, alpha_cond=1.0, alpha_cond1=1.0 real :: diffusion_coefficient=1.0, diffusion_coefficient1=1.0 real :: tau_damp_evap=0.0, tau_damp_evap1=0.0 real :: tau_ocean_driving=0.0, tau_ocean_driving1=0.0 real :: ztop_ocean=0.0, TTocean=300.0 real :: aplow=1.0, aphigh=2.0, mbar=1.0 real :: ap1=1.0, qplaw=0.0, GS_condensation=0. real :: sigma_initdist=0.2, a0_initdist=5e-6, rpbeta0=0.0 integer :: nbin_initdist=20 logical :: lsweepup_par=.false., lcondensation_par=.false. logical :: llatent_heat=.true., lborder_driving_ocean=.false. logical :: lcondensation_simplified=.false. character (len=labellen), dimension(ninit) :: initap='nothing' character (len=labellen) :: condensation_coefficient_type='constant' ! namelist /particles_radius_init_pars/ & initap, ap0, rhopmat, vthresh_sweepup, deltavp12_floor, & lsweepup_par, lcondensation_par, tstart_sweepup_par, cdtps, apmin, & condensation_coefficient_type, alpha_cond, diffusion_coefficient, & tau_damp_evap, llatent_heat, cdtpc, tau_ocean_driving, & lborder_driving_ocean, ztop_ocean, radii_distribution, TTocean, & aplow, aphigh, mbar, ap1, qplaw, eps_dtog, nbin_initdist, & sigma_initdist, a0_initdist, lparticles_radius_rpbeta, rpbeta0 ! namelist /particles_radius_run_pars/ & rhopmat, vthresh_sweepup, deltavp12_floor, & lsweepup_par, lcondensation_par, tstart_sweepup_par, cdtps, apmin, & condensation_coefficient_type, alpha_cond, diffusion_coefficient, & tau_damp_evap, llatent_heat, cdtpc, tau_ocean_driving, & lborder_driving_ocean, ztop_ocean, TTocean, & lcondensation_simplified, GS_condensation, lparticles_radius_rpbeta, rpbeta0 ! integer :: idiag_apm=0, idiag_ap2m=0, idiag_apmin=0, idiag_apmax=0 integer :: idiag_dvp12m=0, idiag_dtsweepp=0, idiag_npswarmm=0 integer :: idiag_ieffp=0 ! contains !*********************************************************************** subroutine register_particles_radius() ! ! Set up indices for access to the fp and dfp arrays. ! ! 22-aug-05/anders: coded ! if (lroot) call svn_id( & "$Id$") ! ! Index for particle radius. ! call append_npvar('iap',iap) ! ! Index for the effectiveness factor of surface reactions ! call append_npaux('ieffp',ieffp) ! if (lparticles_radius_rpbeta) call append_npvar('irpbeta',irpbeta) ! endsubroutine register_particles_radius !*********************************************************************** subroutine initialize_particles_radius(f) ! ! Perform any post-parameter-read initialization i.e. calculate derived ! parameters. ! ! 22-aug-05/anders: coded ! use SharedVariables, only: put_shared_variable ! real, dimension (mx,my,mz,mfarray) :: f ! ! Calculate the number density of bodies within a superparticle. ! if (npart_radii>1 .and. & (.not. lcartesian_coords .or. lparticles_nbody .or. & lparticles_number .or. lparticles_spin)) then call fatal_error('initialize_particles_radius: npart_radii > 1','') else mpmat=4/3.0*pi*rhopmat*ap0(1)**3 if (lroot) print*, 'initialize_particles_radius: '// & 'mass per dust grain mpmat=', mpmat endif ! if ((lsweepup_par.or.lcondensation_par).and..not.lpscalar & .and..not.lcondensation_simplified) then call fatal_error('initialize_particles_radius', & 'must have passive scalar module for sweep-up and condensation') endif ! ! Short hand for spherical particle prefactor. ! four_pi_rhopmat_over_three=four_pi_over_three*rhopmat ! ! Inverse coefficients. ! alpha_cond1=1/alpha_cond diffusion_coefficient1=1/diffusion_coefficient if (tau_damp_evap/=0.0) tau_damp_evap1=1/tau_damp_evap if (tau_ocean_driving/=0.0) tau_ocean_driving1=1/tau_ocean_driving ! call put_shared_variable('ap0',ap0,caller='initialize_particles_radius') ! call keep_compiler_quiet(f) ! endsubroutine initialize_particles_radius !*********************************************************************** subroutine set_particle_radius(f,fp,npar_low,npar_high,init) ! ! Set radius of new particles. ! ! 18-sep-09/nils: adapted from init_particles_radius ! use General, only: random_number_wrapper ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mpar_loc,mparray) :: fp integer :: npar_low,npar_high logical, optional :: init ! real, dimension (mpar_loc) :: r_mpar_loc, p_mpar_loc, tmp_mpar_loc real, dimension (ninit) :: radii_cumulative real, dimension (nbin_initdist) :: n_initdist, a_initdist integer, dimension (nbin_initdist) :: nn_initdist real :: radius_fraction, mcen, mmin, mmax, fcen, p, rhopm real :: lna0, lna1, lna, lna0_initdist integer :: i, j, k, kend, ind, ibin logical :: initial ! initial=.false. if (present(init)) then if (init) initial=.true. endif ! do j=1,ninit ! select case (initap(j)) ! case ('nothing') if (initial.and.lroot.and.j==1) print*, 'set_particles_radius: nothing' ! case ('constant') if (initial.and.lroot) print*, 'set_particles_radius: constant radius' ind=1 do k=npar_low,npar_high if (npart_radii>1) then call random_number_wrapper(radius_fraction) ind=ceiling(npart_radii*radius_fraction) endif fp(k,iap)=ap0(ind) enddo ! case ('constant-1') if (initial.and.lroot) print*, 'set_particles_radius: set particle 1 radius' do k=npar_low,npar_high if (ipar(k)==1) fp(k,iap)=ap1 enddo ! case ('power-law') call random_number_wrapper(fp(npar_low:npar_high,iap)) fp(npar_low:npar_high,iap) = ((aphigh**(qplaw+1.0)-aplow**(qplaw+1.0)) & *fp(npar_low:npar_high,iap)+aplow**(qplaw+1.0))**(1.0/(qplaw+1.0)) ! case ('random') if (initial.and.lroot) print*, 'set_particles_radius: random radius' do k=npar_low,npar_high call random_number_wrapper(fp(k,iap)) fp(k,iap)=fp(k,iap)*(aphigh-aplow)+aplow enddo ! case ('logarithmically-spaced') if (initial.and.lroot) print*, 'set_particles_radius: '// & 'logarithmically spaced with ap0, ap1=', ap0(j), ap1 do k=npar_low,npar_high fp(k,iap)=10**(alog10(ap0(j))+(ipar(k)-0.5)* & (alog10(ap1)-alog10(ap0(j)))/npar) enddo if (lparticles_density) then if (lgravz .and. lgravz_gas) then rhopm = eps_dtog*sqrt(2*pi)*1.0/Lz else rhopm = eps_dtog*1.0 endif do k=npar_low,npar_high aplow =10**(alog10(fp(k,iap))-(alog10(ap1)-alog10(ap0(j)))/npar/2) aphigh=10**(alog10(fp(k,iap))+(alog10(ap1)-alog10(ap0(j)))/npar/2) fp(k,irhopswarm)=(aphigh**(4-qplaw)-aplow**(4-qplaw))/ & (ap1**(4-qplaw)-ap0(j)**(4-qplaw))*rhopm*nwgrid enddo endif ! ! Lognormal distribution. Here, ap1 is the largest value in the distribution. ! Initialize particle radii by a direct probabilistic calculation using ! gaussian noise for ln(a/a0)/sigma. ! case ('lognormal') if (initial.and.lroot) print*, 'set_particles_radius: '// & 'lognormal=', a0_initdist call random_number_wrapper(r_mpar_loc) call random_number_wrapper(p_mpar_loc) tmp_mpar_loc=sqrt(-2*log(r_mpar_loc))*sin(2*pi*p_mpar_loc) fp(:,iap)=a0_initdist*exp(sigma_initdist*tmp_mpar_loc) ! ! Lognormal distribution. Here, ap1 is the largest value in the distribution ! and ap0 is the smallest radius initially. ! case ('old_lognormal') if (initial.and.lroot) print*, 'set_particles_radius: '// & 'lognormal=', ap0(1), ap1 lna0=log(ap0(1)) lna1=log(ap1) lna0_initdist=log(a0_initdist) do ibin=1,nbin_initdist lna=lna0+(lna1-lna0)*(ibin-1)/(nbin_initdist-1) a_initdist(ibin)=exp(lna) n_initdist(ibin)=exp(-0.5*(lna-lna0_initdist)**2/sigma_initdist**2) & /(sqrt(twopi)*sigma_initdist*a_initdist(ibin)) enddo nn_initdist=nint(n_initdist) nn_initdist=(npar_high-npar_low+1)*nn_initdist/(sum(nn_initdist)-1) ! ! is now normalized to the number of particles, ! so set the corresponding distribution. ! Normally, the number of bins is less than the number of available ones, ! so then we patch the rest with fp(kend+1:,iap)=a0_initdist. ! k=npar_low do ibin=1,nbin_initdist kend=k+nn_initdist(ibin) if (kend > k) then fp(k:kend,iap)=a_initdist(ibin) k=kend endif enddo ! ! put all the remaining particles (from kend+1 to the end of the array) ! in the bin corresponding to the middle of the distribution. ! fp(kend+1:,iap)=a0_initdist ! case ('specify') ! ! User specified particle size distribution with constant radii. ! if (initial.and.lroot) & print*, 'set_particles_radius: constant radius, user specified distribution' radii_cumulative=0.0 radii_cumulative(1)=radii_distribution(1) do i=2,npart_radii radii_cumulative(i) = radii_cumulative(i-1) + radii_distribution(i) enddo if (radii_cumulative(npart_radii) /= 1.0) then ! ! Renormalize. ! do i=1,npart_radii radii_cumulative(i)=radii_cumulative(i)/radii_cumulative(npart_radii) enddo endif ! do k=npar_low,npar_high call random_number_wrapper(radius_fraction) do i=1,npart_radii if (radius_fraction <= radii_cumulative(i)) then fp(k,iap)=ap0(i) exit endif enddo enddo ! ! Coagulation test with linear kernel. We initially put particles according ! to the distribution function ! ! fk = dn/dm = n_0/mbar_0*exp(-m/mbar_0) ! => rhok = m_k*n_0/mbar_0*exp(-m/mbar_0) ! ! Integrating rhok over all mass and normalising by n_0*mbar_0 gives ! ! I = int_0^m[rhok]/int_0^oo[rhok] = 1 - exp(-x) - x*exp(-x) ! ! where x=mk/mbar_0. We place particles equidistantly along the integral ! to obtain the right distribution function. ! case ('kernel-lin') if (initial.and.lroot) print*, 'set_particles_radius: '// & 'initial condition for linear kernel test' do k=npar_low,npar_high p=(ipar(k)-0.5)/float(npar) mmin=0.0 mmax=1.0 do while (.true.) if ((1.0-exp(-mmax)*(1+mmax))1.0e-6) if (fcen particle is losing mass ! (the mass vector is pointing out of the particle) ! allocate(mass_loss(k1:k2)) allocate(effectiveness_factor(k1:k2)) ! call get_radius_chemistry(mass_loss,effectiveness_factor) ! if (.not. lsurface_nopores) then do k = k1,k2 if (fp(k,irhosurf) < 0) then rho = fp(k,imp) / (fp(k,iap)**3 * 4./3. * pi ) mass_per_radius = 4. * pi * rho * fp(k,iap)**2 dfp(k,iap) = dfp(k,iap) - mass_loss(k) *(1-effectiveness_factor(k))/mass_per_radius endif fp(k,ieffp) = effectiveness_factor(k) enddo else do k = k1,k2 rho = fp(k,imp) / (fp(k,iap)**3 * 4./3. * pi ) mass_per_radius = 4. * pi * rho * fp(k,iap)**2 dfp(k,iap) = dfp(k,iap) - mass_loss(k)/mass_per_radius enddo endif deallocate(mass_loss) deallocate(effectiveness_factor) endif ! endsubroutine dap_dt_pencil !*********************************************************************** subroutine dap_dt_sweepup_pencil(f,df,fp,dfp,p,ineargrid) ! ! Increase in particle radius due to sweep-up of small grains in the gas. ! ! 22-aug-05/anders: coded ! use Diagnostics, only: max_mn_name use Particles_number ! 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 type (pencil_case) :: p integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension (nx) :: dt1_sweepup real :: deltavp integer :: k, ix0, ix ! intent (in) :: f, fp intent (inout) :: dfp ! ! Increase in particle radius due to sweep-up of small grains in the gas. ! if (t>=tstart_sweepup_par) then ! if (lfirst.and.ldt) dt1_sweepup=0.0 ! if (npar_imn(imn)/=0) then do k=k1_imn(imn),k2_imn(imn) ix0=ineargrid(k,1) ix=ix0-nghost ! No interpolation needed here. ! Relative speed. deltavp=sqrt( & (fp(k,ivpx)-p%uu(ix,1))**2 + & (fp(k,ivpy)-p%uu(ix,2))**2 + & (fp(k,ivpz)-p%uu(ix,3))**2 ) if (deltavp12_floor/=0.0) & deltavp=sqrt(deltavp**2+deltavp12_floor**2) ! Allow boulders to sweep up small grains if relative velocity not too high. if (deltavp<=vthresh_sweepup .or. vthresh_sweepup<0.0) then ! Radius increase due to sweep-up. dfp(k,iap) = dfp(k,iap) + & 0.25*deltavp*p%cc(ix)*p%rho(ix)*rhopmat1 ! ! Deplete gas of small grains. ! if (lparticles_number) np_swarm=fp(k,inpswarm) if (lpscalar_nolog) then df(ix0,m,n,icc) = df(ix0,m,n,icc) - & np_swarm*pi*fp(k,iap)**2*deltavp*p%cc(ix) else df(ix0,m,n,ilncc) = df(ix0,m,n,ilncc) - & np_swarm*pi*fp(k,iap)**2*deltavp endif ! ! Time-step contribution of sweep-up. ! if (lfirst.and.ldt) then dt1_sweepup(ix) = dt1_sweepup(ix) + & np_swarm*pi*fp(k,iap)**2*deltavp endif ! endif ! if (ldiagnos) then if (idiag_dvp12m/=0) call sum_par_name((/deltavp/),idiag_dvp12m) endif enddo endif ! ! Time-step contribution of sweep-up. ! if (lfirst.and.ldt) then dt1_sweepup=dt1_sweepup/cdtps dt1_max=max(dt1_max,dt1_sweepup) if (ldiagnos.and.idiag_dtsweepp/=0) & call max_mn_name(dt1_sweepup,idiag_dtsweepp,l_dt=.true.) endif ! endif ! call keep_compiler_quiet(f) ! endsubroutine dap_dt_sweepup_pencil !*********************************************************************** subroutine dap_dt_condensation_pencil(f,df,fp,dfp,p,ineargrid) ! ! Change in particle radius due to condensation of monomers from the gas and ! evaporation of solid/liquid droplets. ! ! 15-jan-10/anders: coded ! use EquationOfState, only: gamma use Particles_number ! 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 type (pencil_case) :: p integer, dimension (mpar_loc,3) :: ineargrid ! real, dimension (nx) :: ap_equi, vth, dt1_condensation, rhovap real, dimension (nx) :: total_surface_area, ppsat real, dimension (nx) :: rhocond_tot, rhosat, np_total real :: dapdt, drhocdt, alpha_cond_par integer :: k, ix, ix0 ! intent (in) :: f, fp intent (inout) :: dfp ! ! Change in particle radius due to condensation and evaporation. ! if (t>=tstart_condensation_par) then ! if (npar_imn(imn)/=0) then rhovap=p%cc*p%rho ppsat=6.035e11*exp(-5938*p%TT1) ! Valid for water vth=sqrt(p%csvap2) rhosat=gamma*ppsat/p%csvap2 rhocond_tot=p%rhop+rhovap if (lfirst.and.ldt) then np_total=0.0 total_surface_area=0.0 dt1_condensation=0.0 endif do k=k1_imn(imn),k2_imn(imn) ix0=ineargrid(k,1) ix=ix0-nghost ! ! The condensation/evaporation mass flux is ! ! F = vth*rhovap*alpha ! ! where vth is the thermal speed, rhopvap is the mass density of vapor, ! and alpha \in [0,1] is the condensation coefficient. Various authors argue ! for different choices of alpha. Following Barkstrom 1978 we take ! ! alpha = 1/(vth*ap/D + 1/alpha0) ! ! where D is the diffusion coefficient of vapor and ap the particle radius. ! Small particles, with aprhosat(ix)) then dt1_condensation(ix) = max(total_surface_area(ix)*vth(ix), & pi*vth(ix)*np_total(ix)*ap_equi(ix)**2*alpha_cond) endif enddo endif endif ! if (lborder_driving_ocean) then if (z(n)