Wiki

2.6 Modelling the dust using size-bins

KROME can handle the dust by considering it as an additional set of species (and relative ODE) subdivided into different bins that correspond to different grain sizes.
To activate this option in KROME you should add the option -dust=N,types where N is the number of bins, and types are the list of dust types employed, e.g. -dust=10,C,Si will activate 10 bins of carbon-based AND 10 bins of silicon-based dust, for a total of 10x2=20 bins. Note that in this case adding 20 bins means to add 20 more ODE equations to the initial set. The influence on the global computational performance depends on the type of chemistry and microphysics where the dust is involved.
To initialize the bins you just need to write into KROME the following line of code

  call krome_init_dust_distribution(x(:), dust2gas)

that is a function included in the krome_user module. In this way you are telling to KROME that you want a MRN distribution (-3.5 power law) in the range 5e-7 cm and 2.5e-5 cm. The bins will be logarithmically spaced.
x(:) is the array of the abundances in number densities of the chemical species, and dust2gas is the dust-to-gas ratio: these are used by KROME to retrieve the amount of dust.
If you want to change this you may add three additional arguments

  call krome_init_dust_distribution(x(:), dust2gas, smallest, largest, expMRN)

where the first two arguments are the smallest and the largest dust size in cm, and expMRN is the exponential (including its sign) of the power law.

The subroutine krome_set_dust_distribution sets the size distribution in order to have all the bins with the same amount of dust. The default is that the sum of the bins is 1.

To obtain the distribution of dust from KROME use the function

krome_get_dust_distribution()

that returns an array (real8) of size krome_ndust* containing the number densities of all the dust bins.
So for example, in the case -dust=10,C,Si the size of the array will be 10x2=20, since there are 2 dust types, namely C and Si, with 10 bins each.

2.6.1 Freeze-out/evaporation (gas-like rates without bins!)

In KROME it is possible to use size-independent gas-dust reactions as gas-like rate reactions.
If your dust distribution is constant during the time evolution you could use a gas-gas reactions that pretend to be a gas-dust one.
More in detail. Consider e.g. the CO freeze-out reaction

 CO + dust -> CO_dust

where CO_dust is the species frozen onto the dust grain.
KROME can handle CO_dust as it was a gas species, this without using the dust machinery, so without -dust option.
This is possible because there are rates that are specifically designed to do this, check for example
http://kida.obs.u-bordeaux1.fr/networks.html
If you take e.g. Reboussin+2014 network, in the nls_grain.dat there are some species starting with the letter J.
These are species onto the dust phase, e.g. JCO is CO on dust surface.
This workaround avoids to explicitly indicate the dust distribution, since they assume a fixed one (check the reference).
In KROME instead of using J you have to use CO_dust. Your chemical network will look like

 @format:idx,R,P,rate
 1,CO,CO_dust, 1e-8

where 1e-8 is an example!
Analogously for evaporation, we have e.g.

 @format:idx,R,P,rate
 2,CO_dust,CO, 1e12*exp(-1100./Tgas)

where again the rate is an example.

2.6.2 Solver-friendly freeze-out/evaporation (also without bins!)

Including evaporation could cause solver's instability, especially when dust temperature is time-dependent.
To cope with this, we introduced a feature that exploits the fact that n(CO_tot) = n(CO_ice)+n(CO_gas), where this example is for CO, but applies to any ice species (W-F. Thi, private communication).
To activate this feature (for e.g. CO) you have to add this to your network file

@ice:CO,freezeout,krate_stickSi(n(:),idx_CO,Tdust)
@ice:CO,evaporation,krate_evaporation(n(:),idx_CO,Tdust)

where Tdust is the dust temperature (you could assume that Tdust=Tgas, but it clearly depends on the environment).
After doing this you have access to CO in the gas phase, and total CO.
Then, to retrieve the ice content in your main program

myIce = x(krome_idx_CO_total)-x(krome_idx_CO)

In the example above we used KROME default rate functions (i.e.krate_stickSi and krate_evaporation), but in principle you can use whatever you want

@ice:CO,freezeout,YOUR FAVOURITE FREEZE-OUT RATE
@ice:CO,evaporation, YOUR FAVOURITE EVAPORATION RATE

How it works?
The ice and the gas differentials are

dgas = dchem +ke*ice -kf*gas
dice = -ke*ice +kf*gas

where dchem is the gas chemistry, and ke and kf are the evaporation and freeze-out rate coefficients.
If we define tot=gas+ice, we rewrite the previous set as

dtot = dgas+dice = dchem
dgas = dchem - gas*(ke+kf) + tot*ke

that is the same as before, but with a variable change. This form is more stable.

2.6.3 Charged grains

KROME handles also charged grains by using the gas-like method described above.
In particular we have GRAIN0 (neutral), GRAIN+, and GRAIN-, as well as PAH (neutral, note 0 missing), PAH+, PAH-.