From d7388acd38d86a40b7c12955876b30b81219fe70 Mon Sep 17 00:00:00 2001 From: AlagueroA Date: Thu, 11 Jun 2026 15:07:18 +0200 Subject: [PATCH 1/6] attempt to moddump_nucleationtogrowth --- src/main/dust_formation.f90 | 4 +- src/main/growth.f90 | 48 ++++-- src/main/physcon.f90 | 4 + src/utils/moddump_nucleationtodustgrowth.f90 | 145 +++++++++++++++++++ 4 files changed, 185 insertions(+), 16 deletions(-) create mode 100644 src/utils/moddump_nucleationtodustgrowth.f90 diff --git a/src/main/dust_formation.f90 b/src/main/dust_formation.f90 index b9e154999..45edc220d 100644 --- a/src/main/dust_formation.f90 +++ b/src/main/dust_formation.f90 @@ -25,11 +25,12 @@ module dust_formation ! use part, only:idJstar,idK0,idK1,idK2,idK3,idmu,idgamma,idsat,idkappa - use physcon, only:kboltz,pi,atomic_mass_unit,mass_proton_cgs,patm + use physcon, only:kboltz,pi,atomic_mass_unit,mass_proton_cgs,patm,rho_Cdust,a0,mc use dim, only:nElements implicit none integer, public :: idust_opacity = 0 +! real, public :: rho_Cdust = 2.62, a0 = 1.28e-8, mc = 12.*atomic_mass_unit !cgs public :: set_abundances,evolve_dust,evolve_chem,calc_kappa_dust,& calc_kappa_bowen,chemical_equilibrium_light,psat_C,calc_nucleation,& @@ -233,7 +234,6 @@ pure real function calc_kappa_dust(K3, Tdust, rho_cgs) real, intent(in) :: K3, Tdust, rho_cgs real :: kappa_cgs, fac - real, parameter :: rho_Cdust = 2.62, mc = 12.*atomic_mass_unit real, parameter :: Qplanck_abs = 1.6846124267740528e+04 real, parameter :: Qross_ext = 9473.2722815583311 diff --git a/src/main/growth.f90 b/src/main/growth.f90 index 09f97f867..5cc0a71dd 100644 --- a/src/main/growth.f90 +++ b/src/main/growth.f90 @@ -217,10 +217,12 @@ end subroutine print_growthinfo !----------------------------------------------------------------------- subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,dmdt) use part, only:rhoh,idust,igas,iamtype,iphase,isdead_or_accreted,& - massoftype,Omega_k,dustfrac,tstop,deltav + massoftype,Omega_k,dustfrac,tstop,deltav,& + nucleation use options, only:use_dustfrac,use_porosity use physcon, only:fourpi use eos, only:ieos,get_spsound + use dim, only:do_nucleation real, intent(in) :: dustprop(:,:) real, intent(inout) :: dustgasprop(:,:) real, intent(in) :: xyzh(:,:) @@ -244,6 +246,7 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d !$omp shared(npart,iphase,ieos,massoftype,use_dustfrac,dustfrac,use_porosity,grainsizemin) & !$omp shared(ifrag,ieros,utime,umass,dsize,cohacc) & !$omp shared(xyzh,vxyzu,dustprop,dustgasprop,dmdt,filfac,VrelVf,tstop,deltav) & + !$omp shared(do_nucleation,nucleation) & !$omp private(i,iam,rho,rhog,rhod,vrel,sdust) & !$omp private(mass_min,massgrain,rhograin) & !$omp firstprivate(filfaci) @@ -258,7 +261,11 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d rho = rhoh(xyzh(4,i),massoftype(igas)) rhog = rho*(1-dustfrac(1,i)) rhod = rho*dustfrac(1,i) - dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i)) + if (do_nucleation) then + dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i),nucleation(7,i)) + else + dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i)) + endif dustgasprop(2,i) = rhog dustgasprop(3,i) = tstop(1,i) * Omega_k(i) dustgasprop(4,i) = sqrt(deltav(1,1,i)**2 + deltav(2,1,i)**2 + deltav(3,1,i)**2) @@ -596,9 +603,11 @@ end subroutine check_dustprop !----------------------------------------------------------------------- subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_index) use dust, only:grainsizecgs,graindenscgs - use part, only:iamtype,iphase,idust,igas,dustprop,filfac,probastick + use part, only:iamtype,iphase,idust,igas,dustprop,filfac,probastick,nucleation use physcon, only:fourpi use options, only:use_dustfrac + use dim, only:do_nucleation + use physcon, only:rho_Cdust,a0 integer, intent(in) :: npart real, intent(in) :: xyzh(:,:) logical, intent(in), optional :: sizedistrib @@ -609,18 +618,29 @@ subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_i do i=1,npart iam = iamtype(iphase(i)) if (iam == idust .or. (iam == igas .and. use_dustfrac)) then - dustprop(2,i) = graindenscgs / unit_density - if (sizedistrib) then - r = sqrt(xyzh(1,i)**2 + xyzh(2,i)**2) - h = H_R_ref * R_ref * au / udist * (r * udist / au / R_ref)**(1.5-q_index) - dustprop(1,i) = grainsizecgs/udist * (r * udist / au / R_ref)**pwl_sizedistrib * exp(-0.5*xyzh(3,i)**2/h**2) - if (dustprop(1,i) < 2.e-5/udist) then - dustprop(1,i) = 2.e-5/udist - endif - dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3 + + if (do_nucleation) then !used by moddump_nucleationtodustgrowth + dustprop(2,i) = rho_Cdust / unit_density + if (nucleation(2,i) < 1e-30) then !arbitrary threshold + dustprop(1,i) = tiny(dustprop(1,i)) + else + dustprop(1,i) = a0*nucleation(3,i)/nucleation(2,i) / udist !avg grain size = a0 K1/K0 + dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3 ! convert into dust mass + endif else - dustprop(1,i) = fourpi/3. * dustprop(2,i) * (grainsizecgs / udist)**3 - endif + dustprop(2,i) = graindenscgs / unit_density + if (sizedistrib) then + r = sqrt(xyzh(1,i)**2 + xyzh(2,i)**2) + h = H_R_ref * R_ref * au / udist * (r * udist / au / R_ref)**(1.5-q_index) + dustprop(1,i) = grainsizecgs/udist * (r * udist / au / R_ref)**pwl_sizedistrib * exp(-0.5*xyzh(3,i)**2/h**2) + if (dustprop(1,i) < 2.e-5/udist) then + dustprop(1,i) = 2.e-5/udist + endif + dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3 + else + dustprop(1,i) = fourpi/3. * dustprop(2,i) * (grainsizecgs / udist)**3 + endif + endif else dustprop(:,i) = 0. endif diff --git a/src/main/physcon.f90 b/src/main/physcon.f90 index 537b451e0..b6ac23bdc 100644 --- a/src/main/physcon.f90 +++ b/src/main/physcon.f90 @@ -101,4 +101,8 @@ module physcon ! real(kind=8), parameter :: eVtoK = 1.1604519d4 !Degrees kelvin per eV K/eV +! +!--Dust nucleation +! +real(kind=8), parameter :: rho_Cdust = 2.62, a0 = 1.28e-8, mc = 12.*atomic_mass_unit end module physcon diff --git a/src/utils/moddump_nucleationtodustgrowth.f90 b/src/utils/moddump_nucleationtodustgrowth.f90 new file mode 100644 index 000000000..d1061c50c --- /dev/null +++ b/src/utils/moddump_nucleationtodustgrowth.f90 @@ -0,0 +1,145 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2026 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module moddump +! +! default moddump routine: does not make any modifications +! +! :References: None +! +! :Owner: Daniel Price +! +! :Runtime parameters: None +! +! :Dependencies: None +! + implicit none + character(len=*), parameter, public :: moddump_flags = '' + +contains + +subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) +use dim, only:use_dust,maxdusttypes,maxdustlarge,maxdustsmall,use_dustgrowth,& + update_max_sizes,isothermal +use part, only:igas,idust,ndusttypes,ndustsmall,ndustlarge,& + grainsize,graindens,dustfrac,& + nucleation +use dust_formation, only:mass_per_H +use set_dust, only:set_dustfrac,set_dustbinfrac +use options, only:use_dustfrac,use_porosity +use growth, only:set_dustprop,convert_to_twofluid,iporosity,init_growth,ifrag,gsizemincgs +use prompting, only:prompt +use dust, only:grainsizecgs,graindenscgs +use table_utils, only:logspace +use units, only:umass,udist +use physcon, only:rho_Cdust,a0,mc +use cooling_solver, only:excitation_HI + integer, intent(inout) :: npart + integer, intent(inout) :: npartoftype(:) + real, intent(inout) :: massoftype(:) + real, intent(inout) :: xyzh(:,:),vxyzu(:,:) + integer :: i,ierr,j,itype,ipart,iloc,dust_method,np_ratio,np_gas,np_dust,maxdust + real :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens + real :: inradius,outradius,pwl_sizedistrib,R_ref,H_R_ref,q_index + logical :: sizedistrib + + +! cooling values + excitation_HI = 1 + +! dust default values + udens = umass/(udist**3) + dust_method = 1 + np_ratio = 5 + dust_to_gas = 0.01 + ndusttypes = 1 + smincgs = 1.e-5 + smaxcgs = 1. + sindex = 3.5 + dustbinfrac = 0. + pwl_sizedistrib = -2 + R_ref = 100 + H_R_ref = 0.0895 + q_index = 0.25 + + sizedistrib = .false. + + grainsizecgs = 1. !dummy + graindenscgs = 1. + + gsizemincgs = 1.e-8 + + +! Interactive prompts +!call prompt('Which dust method do you want? (1=one fluid,2=two fluid)',dust_method,1,2) +dust_method = 1 +if (dust_method==1) then + maxdust = maxdustsmall +elseif (dust_method==2) then + maxdust = maxdustlarge + call prompt('Enter ratio between number of gas particles and dust particles',np_ratio,1) + !--We do not care if modulo(npart,np_ratio) is stricly zero, since npart can + ! be a weird value depdending on the simulation it comes from. + endif + +! call prompt('Enter total dust to gas ratio',dust_to_gas,0.) computed from nuc. + +!call prompt('How many grain sizes do you want?',ndusttypes,1,maxdust) +ndusttypes = 1 +! if (ndusttypes > 1) then !uncomment if several dust populations +! !--grainsizes +! call prompt('Enter minimum grain size in cm',smincgs,0.) +! call prompt('Enter maximum grain size in cm',smaxcgs,0.) +! !--mass distribution +! call prompt('Enter power-law index, e.g. MRN',sindex) +! call set_dustbinfrac(smincgs/udist,smaxcgs/udist,sindex,dustbinfrac(1:ndusttypes),grainsize(1:ndusttypes)) +! !--grain density +! call prompt('Enter grain density in g/cm^3',graindens(1),0.) +! graindens = graindens(1)/udens +! else +if (use_dustgrowth) then + ifrag = 1 + iporosity = 0 !turn on use_porosity if you change that + grainsize(1) = grainsizecgs/udist + graindens(1) = graindenscgs/udens +endif + +np_gas = npartoftype(igas) + +if (dust_method == 1) then + + use_dustfrac = .true. + ndustsmall = ndusttypes + + do i=1,np_gas + if (ndusttypes > 1) then + dustfrac(1:ndusttypes,i) = dust_to_gas*dustbinfrac(1:ndusttypes) + else + if (nucleation(5,i) Date: Thu, 11 Jun 2026 15:33:07 +0200 Subject: [PATCH 2/6] edits in growth --- src/main/growth.f90 | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/main/growth.f90 b/src/main/growth.f90 index f68aacaa3..f4e9d724b 100644 --- a/src/main/growth.f90 +++ b/src/main/growth.f90 @@ -232,7 +232,7 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d real, intent(inout) :: dustgasprop(:,:) real, intent(in) :: xyzh(:,:) real, intent(in) :: filfac(:) - real, intent(inout) :: VrelVf(:),vxyzu(:,:) + real, intent(inout) :: VrelVf(:,:),vxyzu(:,:) real, intent(out) :: dmdt(:) integer, intent(in) :: npart ! @@ -651,7 +651,8 @@ subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_i logical, intent(in), optional :: sizedistrib real, intent(in), optional :: pwl_sizedistrib,R_ref,H_R_ref,q_index integer :: i,iam - real :: r,h + real :: r,h,grainmassmin + do i=1,npart iam = iamtype(iphase(i)) @@ -659,11 +660,13 @@ subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_i if (do_nucleation) then !used by moddump_nucleationtodustgrowth dustprop(2,i) = rho_Cdust / unit_density - if (nucleation(2,i) < 1e-30) then !arbitrary threshold - dustprop(1,i) = tiny(dustprop(1,i)) + grainmassmin = fourpi/3 * grainsizemin**3 * dustprop(2,i) + if (nucleation(2,i) < 1e-30) then !arbitrary threshold in number of dust grain + dustprop(1,i) = grainmassmin !tiny(dustprop(1,i)) else dustprop(1,i) = a0*nucleation(3,i)/nucleation(2,i) / udist !avg grain size = a0 K1/K0 dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3 ! convert into dust mass + if (dustprop(1,i) Date: Thu, 11 Jun 2026 17:54:02 +0200 Subject: [PATCH 3/6] apparently working version of moddump_nucleationtodustgrowth.f90 --- src/main/growth.f90 | 2 +- src/utils/moddump_nucleationtodustgrowth.f90 | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/main/growth.f90 b/src/main/growth.f90 index f4e9d724b..0346441b3 100644 --- a/src/main/growth.f90 +++ b/src/main/growth.f90 @@ -660,7 +660,7 @@ subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_i if (do_nucleation) then !used by moddump_nucleationtodustgrowth dustprop(2,i) = rho_Cdust / unit_density - grainmassmin = fourpi/3 * grainsizemin**3 * dustprop(2,i) + grainmassmin = fourpi/3 * (1.e-8/udist)**3 * dustprop(2,i) ! arbitrary threshold if (nucleation(2,i) < 1e-30) then !arbitrary threshold in number of dust grain dustprop(1,i) = grainmassmin !tiny(dustprop(1,i)) else diff --git a/src/utils/moddump_nucleationtodustgrowth.f90 b/src/utils/moddump_nucleationtodustgrowth.f90 index d1061c50c..694db7652 100644 --- a/src/utils/moddump_nucleationtodustgrowth.f90 +++ b/src/utils/moddump_nucleationtodustgrowth.f90 @@ -26,11 +26,12 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) update_max_sizes,isothermal use part, only:igas,idust,ndusttypes,ndustsmall,ndustlarge,& grainsize,graindens,dustfrac,& - nucleation + nucleation,dustprop use dust_formation, only:mass_per_H use set_dust, only:set_dustfrac,set_dustbinfrac use options, only:use_dustfrac,use_porosity -use growth, only:set_dustprop,convert_to_twofluid,iporosity,init_growth,ifrag,gsizemincgs +use growth, only:set_dustprop,convert_to_twofluid,iporosity,init_growth,ifrag,gsizemincgs,& + ivrelkin use prompting, only:prompt use dust, only:grainsizecgs,graindenscgs use table_utils, only:logspace @@ -71,6 +72,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) graindenscgs = 1. gsizemincgs = 1.e-8 + ivrelkin = 0 ! relative dust motions from 1-fluid dust velocities not yet available ! Interactive prompts From 276f40d63bf5a05d255ba8df22d66cd44167cf4c Mon Sep 17 00:00:00 2001 From: AlagueroA Date: Thu, 11 Jun 2026 17:55:07 +0200 Subject: [PATCH 4/6] renamed moddump_nucleationtodustgrowth.f90 into moddump_nucleationtogrowth.f90 --- ..._nucleationtodustgrowth.f90 => moddump_nucleationtogrowth.f90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/utils/{moddump_nucleationtodustgrowth.f90 => moddump_nucleationtogrowth.f90} (100%) diff --git a/src/utils/moddump_nucleationtodustgrowth.f90 b/src/utils/moddump_nucleationtogrowth.f90 similarity index 100% rename from src/utils/moddump_nucleationtodustgrowth.f90 rename to src/utils/moddump_nucleationtogrowth.f90 From dd252ee25e37736c33c9519607ad4c487b837264 Mon Sep 17 00:00:00 2001 From: AlagueroA Date: Fri, 12 Jun 2026 14:21:59 +0200 Subject: [PATCH 5/6] added utility function to get dynamical timescale in order to compute St in a more modular way --- src/main/growth.f90 | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/main/growth.f90 b/src/main/growth.f90 index 0346441b3..4e4d2c533 100644 --- a/src/main/growth.f90 +++ b/src/main/growth.f90 @@ -238,7 +238,7 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d ! real :: rhog,rhod,vrel,rho,sdust real :: massgrain,rhograin,filfaci - real :: dtarb,frac_masschange,att_factor + real :: dtarb,frac_masschange,att_factor,tdyn integer :: i,iam @@ -254,13 +254,15 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d !$omp shared(ifrag,ieros,utime,umass,dsize,cohacc,dtmax) & !$omp shared(xyzh,vxyzu,dustprop,dustgasprop,dmdt,filfac,VrelVf,tstop,deltav,Vrel_disp) & !$omp shared(do_nucleation,nucleation) & - !$omp private(i,iam,rho,rhog,rhod,vrel,sdust) & + !$omp private(i,iam,rho,rhog,rhod,vrel,sdust,tdyn) & !$omp private(massgrain,rhograin,dtarb,frac_masschange,att_factor) & !$omp firstprivate(filfaci) do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then iam = iamtype(iphase(i)) + call get_dynamical_timescale(i,tdyn) + if (iam == idust .or. (iam == igas .and. use_dustfrac)) then if (use_dustfrac .and. iam == igas) then @@ -274,7 +276,7 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i)) endif dustgasprop(2,i) = rhog - dustgasprop(3,i) = tstop(1,i) * Omega_k(i) + dustgasprop(3,i) = tstop(1,i) / tdyn dustgasprop(4,i) = sqrt(deltav(1,1,i)**2 + deltav(2,1,i)**2 + deltav(3,1,i)**2) else rhod = rhoh(xyzh(4,i),massoftype(idust)) @@ -1094,4 +1096,24 @@ real function get_size(mass,dens,filfac) end function get_size + +!----------------------------------------------------------------------- +!+ +! Compute a characteristic dynamical timescale. Used in growth.f90 to get the St. +!+ +!----------------------------------------------------------------------- + subroutine get_dynamical_timescale(i,tdyn) + use part, only:xyzmh_ptmass,ivwind,iReff,Omega_k + use dim, only:do_nucleation + integer, intent(in) :: i + real, intent(inout) :: tdyn + + if (do_nucleation) then + tdyn = 1/(xyzmh_ptmass(ivwind,1)/xyzmh_ptmass(iReff,1)) + else + tdyn = 1/Omega_k(i) + endif + + end subroutine get_dynamical_timescale + end module growth From 9f9150eaf750c3990fe28cb177f828b1ec2db1da Mon Sep 17 00:00:00 2001 From: AlagueroA Date: Fri, 12 Jun 2026 14:57:11 +0200 Subject: [PATCH 6/6] cleaned moddump_nucleationtogrowth.f90 --- src/utils/moddump_nucleationtogrowth.f90 | 76 +++++++++++------------- 1 file changed, 35 insertions(+), 41 deletions(-) diff --git a/src/utils/moddump_nucleationtogrowth.f90 b/src/utils/moddump_nucleationtogrowth.f90 index 694db7652..cdb0fb836 100644 --- a/src/utils/moddump_nucleationtogrowth.f90 +++ b/src/utils/moddump_nucleationtogrowth.f90 @@ -6,7 +6,7 @@ !--------------------------------------------------------------------------! module moddump ! -! default moddump routine: does not make any modifications +! Convert a dust-nucleation file into a dust-growth file. Initial conditions of the dust-growth simulation are provided by the dust-nucleation simulation. ! ! :References: None ! @@ -23,13 +23,13 @@ module moddump subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) use dim, only:use_dust,maxdusttypes,maxdustlarge,maxdustsmall,use_dustgrowth,& - update_max_sizes,isothermal -use part, only:igas,idust,ndusttypes,ndustsmall,ndustlarge,& + update_max_sizes,isothermal,do_nucleation +use part, only:igas,idust,ndusttypes,ndustsmall,& grainsize,graindens,dustfrac,& - nucleation,dustprop + nucleation use dust_formation, only:mass_per_H use set_dust, only:set_dustfrac,set_dustbinfrac -use options, only:use_dustfrac,use_porosity +use options, only:use_dustfrac use growth, only:set_dustprop,convert_to_twofluid,iporosity,init_growth,ifrag,gsizemincgs,& ivrelkin use prompting, only:prompt @@ -42,66 +42,60 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) integer, intent(inout) :: npartoftype(:) real, intent(inout) :: massoftype(:) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - integer :: i,ierr,j,itype,ipart,iloc,dust_method,np_ratio,np_gas,np_dust,maxdust - real :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens - real :: inradius,outradius,pwl_sizedistrib,R_ref,H_R_ref,q_index + integer :: i,dust_method,np_ratio,np_gas,maxdust + real :: dust_to_gas,smincgs,smaxcgs,sindex,udens + real :: pwl_sizedistrib,R_ref,H_R_ref,q_index logical :: sizedistrib -! cooling values - excitation_HI = 1 +! Add flag for dust nucleation + if (.not. (do_nucleation .and. use_dustgrowth .and. (.not. isothermal))) then + print*,' DOING NOTHING: COMPILE WITH DUSTGROWTH=yes DUST_NUCLEATION=yes ISOTHERMAL=no' + stop + endif -! dust default values +! units udens = umass/(udist**3) - dust_method = 1 - np_ratio = 5 - dust_to_gas = 0.01 - ndusttypes = 1 + +! parameters + ndusttypes = 1 ! growth can only work with 1 dust population + +! cooling + excitation_HI = 1 + +! dust + dust_method = 1 ! 1:dust as a mxiture ; 2:dust as particles (not working yet) + np_ratio = 5 ! np_dust/np_gas if dust_method=2 + dust_to_gas = 0.01 ! only to inialise the variable, changed later + + sizedistrib = .false. ! initial grain size distribution smincgs = 1.e-5 smaxcgs = 1. sindex = 3.5 - dustbinfrac = 0. - pwl_sizedistrib = -2 + + pwl_sizedistrib = -2 ! spatial distribution of dust, if sizedistrib = .true. R_ref = 100 H_R_ref = 0.0895 q_index = 0.25 - sizedistrib = .false. + gsizemincgs = 1.e-8 ! minimum allowed grainsize in cm + ivrelkin = 0 ! keep 0 until relative dust motions from 1-fluid are available - grainsizecgs = 1. !dummy + grainsizecgs = 1. !dummies but need definition to be used by dust.f90 graindenscgs = 1. - gsizemincgs = 1.e-8 - ivrelkin = 0 ! relative dust motions from 1-fluid dust velocities not yet available - ! Interactive prompts -!call prompt('Which dust method do you want? (1=one fluid,2=two fluid)',dust_method,1,2) -dust_method = 1 if (dust_method==1) then maxdust = maxdustsmall elseif (dust_method==2) then - maxdust = maxdustlarge - call prompt('Enter ratio between number of gas particles and dust particles',np_ratio,1) + print*,' DOING NOTHING: NOT YET READY WITH DUST AS PARTICLES' + stop !--We do not care if modulo(npart,np_ratio) is stricly zero, since npart can ! be a weird value depdending on the simulation it comes from. endif -! call prompt('Enter total dust to gas ratio',dust_to_gas,0.) computed from nuc. - -!call prompt('How many grain sizes do you want?',ndusttypes,1,maxdust) ndusttypes = 1 -! if (ndusttypes > 1) then !uncomment if several dust populations -! !--grainsizes -! call prompt('Enter minimum grain size in cm',smincgs,0.) -! call prompt('Enter maximum grain size in cm',smaxcgs,0.) -! !--mass distribution -! call prompt('Enter power-law index, e.g. MRN',sindex) -! call set_dustbinfrac(smincgs/udist,smaxcgs/udist,sindex,dustbinfrac(1:ndusttypes),grainsize(1:ndusttypes)) -! !--grain density -! call prompt('Enter grain density in g/cm^3',graindens(1),0.) -! graindens = graindens(1)/udens -! else if (use_dustgrowth) then ifrag = 1 iporosity = 0 !turn on use_porosity if you change that @@ -118,7 +112,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) do i=1,np_gas if (ndusttypes > 1) then - dustfrac(1:ndusttypes,i) = dust_to_gas*dustbinfrac(1:ndusttypes) + print*,' DOING NOTHING: GROWTH NOT WORKING FOR MORE THAN 1 DUST POPULATION' else if (nucleation(5,i)