! $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 ! 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(ndustrad) :: ap0=0.0 real, dimension(ndustrad) :: 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, apmid=1.5, aphigh=2.0, mbar=1.0 real :: ap1=1.0, qplaw=0.0,vapor_mixing_ratio_qvs=0., GS_condensation=0.0 real :: modified_vapor_diffusivity=6.74e-6, n0_mean = 1.e8 real, pointer :: G_condensation, ssat0 real :: sigma_initdist=0.2, a0_initdist=5e-6, rpbeta0=0.0 real :: xi_accretion = 0., lambda = 1. real :: n0 = 1., alpha = 1., b=1. integer :: nbin_initdist=20, ip1=npar/2 logical :: lsweepup_par=.false., lcondensation_par=.false. logical :: llatent_heat=.true., lborder_driving_ocean=.false. logical :: lcondensation_simplified=.false. logical :: lcondensation_rate=.false., ldt_evaporation=.false. logical :: ldt_condensation=.false., ldt_condensation_off=.false. logical :: lconstant_radius_w_chem=.false. logical :: lfixed_particles_radius=.false. logical :: ltauascalar = .false., ldust_condensation=.false. logical :: ldust_accretion = .false., lreinitialize_ap=.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, apmid, aphigh, mbar, ap1, ip1, qplaw, eps_dtog, nbin_initdist, & sigma_initdist, a0_initdist, rpbeta0, lparticles_radius_rpbeta, & lfixed_particles_radius, lambda, & n0, b, alpha ! 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, rpbeta0, & lfixed_particles_radius, & lconstant_radius_w_chem, & initap, lreinitialize_ap, ap0, & lcondensation_rate, vapor_mixing_ratio_qvs, & ltauascalar, modified_vapor_diffusivity, ldt_evaporation, & ldt_condensation, ldt_condensation_off, & ldust_condensation, xi_accretion, ldust_accretion, & tstart_condensation_par ! integer :: idiag_apm=0, idiag_ap2m=0, idiag_ap3m=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 ! if (lparticles_chemistry) call append_npaux('ieffp',ieffp) ! if (lparticles_radius_rpbeta) call append_npvar('irpbeta',irpbeta) ! endsubroutine register_particles_radius !*********************************************************************** subroutine initialize_particles_radius(f,fp) ! ! Perform any post-parameter-read initialization i.e. calculate derived ! parameters. ! ! 22-aug-05/anders: coded ! use SharedVariables, only: put_shared_variable, get_shared_variable use General, only: random_number_wrapper ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(mpar_loc,mparray) :: fp real :: radius_fraction integer :: j, k, ind integer :: pos ! ! reinitialize ! if (lreinitialize_ap) then ! ! Copied the following 6 code lines from subroutine read_particles_rad_init_pars, ! which is not called on lreinitialize_ap, so therefore the following lines here. ! Find how many different particle radii we are using. This must be done ! because not all parts of the code are adapted to work with more than one ! particle radius. ! npart_radii=0 do pos = 1,ndustrad if (ap0(pos) /= 0) then npart_radii = npart_radii+1 endif enddo ! do j = 1,ninit select case (initap(j)) ! case ('constant') if (lroot) print*, 'set_particles_radius: constant radius' do k=1,npar_loc 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 ('nothing') if (lroot .and. j == 1) print*, 'set_particles_radius: nothing' ! case default if (lroot) print*, 'initialize_particles_radius: '// & 'No such such value for initap: ', trim(initap(j)) call fatal_error('initialize_particles_radius','') endselect enddo endif ! ! Calculate the number density of bodies within a superparticle. ! if (npart_radii > 1 .and. & (.not. lcartesian_coords .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. lascalar & .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') ! ! If we have decided to hold the radius of the particles to be fixed then ! we should not have any process that changes the radius. ! if (lfixed_particles_radius) then if (lsweepup_par) & call fatal_error('initialize_particles_radius', & 'incosistency: lfixed_particles_radius and lsweepup_par cannot both be true') if (lcondensation_par) & call fatal_error('initialize_particles_radius', & 'incosistency: lfixed_particles_radius and lcondensation_par cannot both be true') if (lparticles_chemistry) & call fatal_error('initialize_particles_radius', & 'incosistency: lfixed_particles_radius and lparticles_chemistry cannot both be true') endif ! call keep_compiler_quiet(f) ! if (lascalar) then call get_shared_variable('G_condensation', G_condensation) if (lcondensation_rate) call get_shared_variable('ssat0', ssat0) endif ! 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(ndustrad) :: radii_cumulative real, dimension(nbin_initdist) :: n_initdist, a_initdist integer, dimension(nbin_initdist) :: nn_initdist real :: radius_fraction, mcen, mmin, mmax, fcen, p real :: lna0, lna1, lna, lna0_initdist, fp_sum_temp 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 ! ! Set lucky droplet radius, ap1 ! 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 ! ! Set lucky droplet radius, one per processor ! case ('constant-proc') if (initial .and. lroot) print*, 'set_particles_radius: set lucky droplet radius per processor' fp(npar_low,iap) = ap1 ! case ('2-size') if (initial .and. lroot) print*, 'set_particles_radius: give particles two radii' do k = npar_low,npar_high if (ipar(k)<=ip1) then fp(k,iap)=aplow else fp(k,iap)=aphigh endif enddo ! case ('2-size-alternate') if (initial .and. lroot) print*, 'set_particles_radius: give particles alternating two radii' do k = npar_low,npar_high if (mod(k,2)==0) then fp(k,iap)=aplow else fp(k,iap)=aphigh endif enddo ! case ('3-size-alternate') if (initial .and. lroot) print*, 'set_particles_radius: give particles alternating three radii' do k = npar_low,npar_high if (mod(k,3)==0) then fp(k,iap)=aplow else if (mod(k,3)==1) then fp(k,iap)=apmid else fp(k,iap)=aphigh endif enddo ! 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 ! lspd: if (lparticles_density) then lsqp: if (qplaw /= 4) then lsk: do k = npar_low, npar_high aplow = 10**(alog10(fp(k,iap)) - 0.5 * (alog10(ap1) - alog10(ap0(j))) / npar) aphigh = 10**(alog10(fp(k,iap)) + 0.5 * (alog10(ap1) - alog10(ap0(j))) / npar) fp(k,irhopswarm) = (aphigh**(4-qplaw) - aplow**(4-qplaw)) / & (ap1**(4-qplaw) - ap0(j)**(4-qplaw)) * rhop_swarm * real(npar) enddo lsk else lsqp fp(npar_low:npar_high,irhopswarm) = rhop_swarm endif lsqp endif lspd ! ! exponential distribution, which could be extended to the Gamma distribution case ('exponential') ! if (initial .and. lroot) print*, 'set_particles_radius: '// & 'exponential=', a0_initdist call random_number_wrapper(fp(npar_low:npar_high,iap)) fp(npar_low:npar_high,iap) = lambda*exp(-lambda*(fp(npar_low:npar_high,iap)-a0_initdist)) ! ! exponential distribution, which could be extended to the Gamma distribution case ('weibull') ! if (initial .and. lroot) print*, 'set_particles_radius: '// & 'weibull=', a0_initdist call random_number_wrapper(fp(npar_low:npar_high,iap)) fp_sum_temp = sum(fp(npar_low:npar_high,iap)) fp(npar_low:npar_high,iap) = a0_initdist*(1./alpha)*(-log(fp(npar_low:npar_high,iap)))**(1./b)/fp_sum_temp ! ! 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)) < p) then mmax = mmax+1.0 else exit endif enddo ! mcen = 0.5*(mmin+mmax) fcen = 1.0-exp(-mcen)*(1+mcen) ! do while (abs(p-fcen) > 1.0e-6) if (fcen < p) then mmin = mcen else mmax = mcen endif mcen = 0.5*(mmin+mmax) fcen = 1.0-exp(-mcen)*(1+mcen) enddo ! fp(k,iap) = (mcen*mbar/four_pi_rhopmat_over_three)**(1.0/3.0) ! 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 default if (lroot) print*, 'init_particles_radius: '// & 'No such such value for initap: ', trim(initap(j)) call fatal_error('init_particles_radius','') endselect enddo ! ! Set initial particle radius if lparticles_mass=T ! if (lparticles_mass) fp(:,iapinit) = fp(:,iap) ! if (lparticles_radius_rpbeta) & fp(npar_low:npar_high,irpbeta) = rpbeta0/(fp(npar_low:npar_high,iap)*rhopmat) ! call keep_compiler_quiet(f) ! endsubroutine set_particle_radius !*********************************************************************** subroutine pencil_criteria_par_radius() ! ! All pencils that the Particles_radius module depends on are specified here. ! ! 21-nov-06/anders: coded ! if (lsweepup_par) then lpenc_requested(i_uu) = .true. lpenc_requested(i_rho) = .true. lpenc_requested(i_cc) = .true. endif if (lcondensation_par) then lpenc_requested(i_csvap2) = .true. lpenc_requested(i_TT1) = .true. lpenc_requested(i_ppvap) = .true. lpenc_requested(i_rho) = .true. lpenc_requested(i_rho1) = .true. lpenc_requested(i_cc) = .true. lpenc_requested(i_cc1) = .true. lpenc_requested(i_np) = .true. lpenc_requested(i_rhop) = .true. if (ltemperature) then lpenc_requested(i_cv1) = .true. lpenc_requested(i_TT1) = .true. lpenc_requested(i_TT) = .true. endif if (lascalar) then lpenc_requested(i_acc) = .true. lpenc_requested(i_ssat) = .true. endif endif ! endsubroutine pencil_criteria_par_radius !*********************************************************************** subroutine dap_dt_pencil(f,df,fp,dfp,p,ineargrid) ! ! Evolution of particle radius. ! ! 22-aug-05/anders: coded ! 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 logical :: lfirstcall=.true., lheader integer :: k, k1, k2 real :: mass_per_radius, rho real, dimension(:), allocatable :: effectiveness_factor, mass_loss ! intent(in) :: f intent(out) :: dfp intent(inout) :: fp ! ! Print out header information in first time step. ! lheader = lfirstcall .and. lroot .and.(.not. lpencil_check_at_work) ! ! Identify module. ! if (lheader) print*,'dap_dt_pencil: Calculate dap/dt' ! if (lsweepup_par) call dap_dt_sweepup_pencil(f,df,fp,dfp,p,ineargrid) if (lcondensation_par) & call dap_dt_condensation_pencil(f,df,fp,dfp,p,ineargrid) ! ! lfirstcall = .false. ! ! Decrease particle radius with dmass if the density in the outer shell ! is wholly consumed (equation 51 in Equations to be solved for DNS ! reactive particles in a turbulent flow) ! if (lparticles_chemistry .and. npar_imn(imn) /= 0) then k1 = k1_imn(imn) k2 = k2_imn(imn) ! ! Particles can lose radius under consumption ! if (.not. lconstant_radius_w_chem) then ! ! mass loss has a positive value -> 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) else ! ! Constant particle radius with activated chemistry ! dfp(k1:k2,iap) = 0.0 ! endif 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 (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, rho0 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, dimension(nx) :: tau_phase1, tau_evaporation1, tau_evaporation 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(:,1)*p%rho !DMDM 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 tau_phase1 = 0.0 tau_evaporation1 = 0.0 tau_evaporation = 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 ap rhosat(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 elseif (lfirst .and. ldt .and. lascalar .and. ldt_condensation) then do ix = 1,nx dt1_condensation(ix) = tau_phase1(ix) if (ldt_evaporation) dt1_condensation(ix) = max(tau_phase1(ix), tau_evaporation1(ix)) enddo elseif (lfirst .and. ldt .and. lcondensation_simplified .and. ldt_condensation) then do ix = 1,nx if (ldt_evaporation) then tau_evaporation1(ix) = -(2*GS_condensation)/ap0(1)**2 dt1_condensation(ix) = max(tau_phase1(ix), tau_evaporation1(ix)) else dt1_condensation(ix) = tau_phase1(ix) endif enddo elseif (lfirst .and. ldt .and. ldt_condensation_off) then dt1_condensation = 0. endif endif ! if (lborder_driving_ocean) then if (z(n) < ztop_ocean) then df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) - & (1.0-TTocean*p%TT1)*tau_ocean_driving1 endif endif ! ! Time-step contribution of condensation. ! if (lfirst .and. ldt) then dt1_condensation = dt1_condensation/cdtpc dt1_max = max(dt1_max,dt1_condensation) endif ! endif ! call keep_compiler_quiet(f) ! endsubroutine dap_dt_condensation_pencil !*********************************************************************** subroutine dap_dt(f,df,fp,dfp,ineargrid) ! ! Evolution of particle radius. ! ! 21-nov-06/anders: coded use Sub ! 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 integer, dimension(mpar_loc,3) :: ineargrid ! ! Diagnostic output. ! if (ldiagnos) then if (idiag_apm /= 0) call sum_par_name(fp(1:npar_loc,iap),idiag_apm) if (idiag_ap2m /= 0) call sum_par_name(fp(1:npar_loc,iap)**2,idiag_ap2m) if (idiag_ap3m /= 0) call sum_par_name(fp(1:npar_loc,iap)**3,idiag_ap3m) if (idiag_apmin /= 0) & call max_par_name(-fp(1:npar_loc,iap),idiag_apmin,lneg=.true.) if (idiag_apmax /= 0) call max_par_name(fp(1:npar_loc,iap),idiag_apmax) if (idiag_npswarmm /= 0) & call sum_par_name(rhop_swarm/ & (four_pi_rhopmat_over_three*fp(1:npar_loc,iap)**3),idiag_npswarmm) if (idiag_ieffp /= 0) & call sum_par_name(fp(1:npar_loc,ieffp),idiag_ieffp) endif ! call keep_compiler_quiet(f,df) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine dap_dt !*********************************************************************** subroutine read_particles_rad_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat integer :: pos ! read (parallel_unit, NML=particles_radius_init_pars, IOSTAT=iostat) ! ! Find how many different particle radii we are using. This must be done ! because not all parts of the code are adapted to work with more than one ! particle radius. ! do pos = 1,ndustrad if (ap0(pos) /= 0) then npart_radii = npart_radii+1 endif enddo ! endsubroutine read_particles_rad_init_pars !*********************************************************************** subroutine write_particles_rad_init_pars(unit) ! integer, intent(in) :: unit ! write (unit, NML=particles_radius_init_pars) ! endsubroutine write_particles_rad_init_pars !*********************************************************************** subroutine read_particles_rad_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read (parallel_unit, NML=particles_radius_run_pars, IOSTAT=iostat) ! endsubroutine read_particles_rad_run_pars !*********************************************************************** subroutine write_particles_rad_run_pars(unit) ! integer, intent(in) :: unit ! write (unit, NML=particles_radius_run_pars) ! endsubroutine write_particles_rad_run_pars !*********************************************************************** subroutine rprint_particles_radius(lreset,lwrite) ! ! Read and register print parameters relevant for particles radius. ! ! 22-aug-05/anders: coded ! use Diagnostics, only: parse_name ! logical :: lreset logical, optional :: lwrite ! integer :: iname ! ! Reset everything in case of reset. ! if (lreset) then idiag_apm = 0 idiag_ap2m = 0 idiag_ap3m = 0 idiag_apmin = 0 idiag_apmax = 0 idiag_dvp12m = 0 idiag_dtsweepp = 0 idiag_npswarmm = 0 idiag_ieffp = 0 endif ! ! Run through all possible names that may be listed in print.in. ! if (lroot .and. ip < 14) & print*, 'rprint_particles_radius: run through parse list' do iname = 1,nname call parse_name(iname,cname(iname),cform(iname),'apm',idiag_apm) call parse_name(iname,cname(iname),cform(iname),'ap2m',idiag_ap2m) call parse_name(iname,cname(iname),cform(iname),'ap3m',idiag_ap3m) call parse_name(iname,cname(iname),cform(iname),'apmin',idiag_apmin) call parse_name(iname,cname(iname),cform(iname),'apmax',idiag_apmax) call parse_name(iname,cname(iname),cform(iname),'dvp12m',idiag_dvp12m) call parse_name(iname,cname(iname),cform(iname),'dtsweepp', & idiag_dtsweepp) call parse_name(iname,cname(iname),cform(iname),'ieffp',idiag_ieffp) if (.not. lparticles_number) call parse_name(iname,cname(iname), & cform(iname),'npswarmm',idiag_npswarmm) enddo ! endsubroutine rprint_particles_radius !*********************************************************************** subroutine get_stbin(api,iStbin) integer, intent(out) :: iStbin integer :: k=0 real :: api k = 1 if (lfixed_particles_radius) then do while((api >= ap0(k)).and.(k <= ndustrad)) iStbin = k k = k+1 enddo endif endsubroutine get_stbin !*********************************************************************** subroutine get_mass_from_radius(mpi,fpwn,ip) real, dimension (lpar_max,mparray) :: fpwn integer,intent(in) :: ip real,intent(out) :: mpi real :: api api = fpwn(ip,iap) mpi=(4./3.)*pi*rhopmat*(api**3) endsubroutine get_mass_from_radius !*********************************************************************** subroutine get_maxrad(apmax) real :: apmax apmax=maxval(ap0) endsubroutine get_maxrad !*********************************************************************** endmodule Particles_radius