Wiki

2.5 Getting started with the bin-based photochemistry.

KROME is capable of using photo chemistry based on frequency/energy bins. This is a standard approach useful for 3D ray-tracing codes, where the number of frequencies is very limited. In order to do that you have to set the frequency intervals and the corresponding intensity intervals. To activate this mode you should run KROME pre-processor (./krome) with the option -photoBins=N, where N is the number of bins that you want to use (e.g. -photoBins=10). Now in your framework code you need to initialize these bins (both frequency/energy and intensity). There are several functions/subroutine to do that. (Here you find a list of the photochemistry-related functions)

Probably, the most straightforward is to employ

call krome_set_photoBin_J21log(lowest, highest)

where lowest and highest are the smallest and the largest energy (eV) values of your distribution, e.g.

call krome_set_photoBin_J21log(13.6d0, 1d3)

This subroutine does all the job. It fills the internal energy bins array with the values (using a log spacing) and the corresponding intensity bin array with J = 6.2415d-10 * (13.6d0/E)**1.5 ev/s/cm2/sr/Hz. If you want a linear spacing energy distribution simply use the similar subroutine

call krome_set_photoBin_J21lin(13.6d0, 1d3)

where note "lin" instead of "log". Now, if you desire to display what KROME has initialized write

phl(:) = krome_get_photoBinE_left() !returns left bin points
phr(:) = krome_get_photoBinE_right() !returns right bin points
phm(:) = krome_get_photoBinE_mid() !returns left middle points
phJ(:) = krome_get_photoBinJ() !returns bin intensities
do i=1,krome_nPhotoBins
   print '(I5,99E17.8)',i,phl(i),phm(i),phr(i),phJ(i)
end do

where the ph* arrays are real*8 arrays of size krome_nPhotoBins (a common variable in the krome_user module). krome_get_photoBinE_left and krome_get_photoBinE_right() return two arrays containing the left and the right energy values for each bin stored inside KROME, while krome_get_photoBinE_mid() returns an array of the middle points (right+left)/2. Energies are in eV. Analogously, krome_get_photoBinJ() returns an array with the corresponding intensities in eV/s/cm2/sr/Hz.

If you want to set your own distribution you can use

krome_set_photobinE_lin(lower, upper)

that initializes a linear energy distribution of bins in the range [lower,upper] or

krome_set_photobinE_log(lower, upper)

which is the same as above, but with a logarithmic distribution. Than you can decide your own custom corresponding intensity bins distribution, e.g.

phm(:) = krome_get_photoBinE_mid() !returns left middle points
myJ(:) = 6.2415d-10 * (13.6d0/phm(:))**1.1 !computes custom intensities
call krome_set_photobinJ(myJ(:)) !initialize custom intensities

where the intensity must be in eV/s/cm2/sr/Hz. When you have set both the energy/frequency bin array, and the corresponding intensities, the initialization is done!

Note that every time you change the intensity the rate coefficients are automatically recalculated.

2.5.1 Change intensity at run-time

Once you have set the intensity by using the subroutine krome_set_photobinJ, you can rescale the function in order to change the flux during the execution of the program. If you want to rescale (i.e. multiply) all the bins by a given double value (e.g. xscale) just use

call krome_photoBin_scale(xscale)

You can also multiply each frequency bin by a different value. To do this use

call krome_photoBin_scale_array(arrayscale(:))

where arrayscale(:) is an array of double of size krome_nPhotoBins. This is very useful to rescale the intensity according to the opacity found.

Note that the two functions indicated above affect the intensity: to restore the intensity initialized with the subroutine krome_get_photoBinJ or any other "automatic" subroutine (e.g. krome_set_photoBin_J21lin) use

call krome_photoBin_restore()

Intensity is stored automatically when the any of the intesity function is created (e.g. krome_set_photoBin_J21lin), but you can store the values to be restored at any time in your code by calling

call krome_photoBin_store()

An example is

!initialize a Draine flux in the range 5-13.6 eV
call krome_set_photoBin_draineLin(5d0, 13.6d0)

!scale flux by 10 times
call krome_photoBin_scale(1d1)

!store current flux (i.e. Drainex10)
call krome_photoBin_store()

!scale flux 50 times MORE (i.e. Drainex10x50)
call krome_photoBin_scale(5d1)

!return to stored flux (i.e. Drainex10)
call krome_photoBin_restore()

2.5.2 Opacity

KROME computes self-consistent opacity using effective chemistry, photochemical cross sections and dust data.

Note that:

when -dust option is enabled all these functions return both gas and dust opacity.
these functions are mutually exclusive

There are several user functions:

To return gas opacity from cross sections and species abundances use

krome_get_opacity(n(:),Tgas)

where n(:) is the real*8 array with the species abundances (same as KROME input), and Tgas is the gas temperature. It returns the opacity per each photobin (i.e. an array of size krome_nPhotoBins). This function assumes the num2col method employed to estimate the column density. More methods with non-assumed column density are explained here below.

If you have more information on the geometry you can use

krome_get_opacity_size(n(:), Tgas, dr)

which is as above, but it computes the column density assuming a gas depth of dr (i.e. surface_density = number_density*dr). It returns the opacity per each photobin (i.e. an array of size krome_nPhotoBins). If you are using dust bins, consistent opacity from dust is automatically included.

If you want to include dust using opacity tables, first load the tables using

call krome_load_opacity_table(fileName)

the file (with path fileName) should be column 1 wavelength (in micron), column 2 opacity (cm^2/g), e.g. draine opacity.

The opacity function works as above

krome_get_opacity_size_d2g(n(:), Tgas, dr, dust2gas)

the only difference is given by the dust-to-gas mass ratio, that should be provided.

2.5.3 Photoheating

To compute the photoheating associated to the rate just add PHOTO in the list of the heating options, e.g. -heating=PHOTO as usual.

2.5.4 Cross section from fortran expression

There are several ways to tell to KROME what cross section to use. In this case we use a fortran expression in the network file:

@photo_start
 1,H,H+,E,6.3d-18/(energy_eV/13.6)**3
@photo_end

As you can see, it is typed exactly like a rate, but the difference here is that the independent variable is not Tgas but energy_eV (clearly in eV). Note that the cross section is in units of cm^2. Don't forget to wrap the rate(s) around the @photo_start @photo_stop tokens.

2.5.5 Cross section from file

KROME allows to use cross sections from file. They should be provided in a two column file: (i) energy in eV, and (ii) cross section in cm2 as usual. The rate is then indicated in the reaction file by using the token @xsecFile as in the following example:

@photo_start
 1,C,C+,E,@xsecFile=photoC.dat
 2,O,O+,E,@xsecFile=photoO.dat
@photo_end

Do not forget to wrap the reactions in the photo_start/stop token (as in the example above).

2.5.6 Cross section from SWRI file

KROME can use the cross section from the SWRI photochemistry database at http://phidrates.space.swri.edu. For example, to include water photodissociation/ionization write

@photo_start
 1,H2O,H,OH,@xsecFile=SWRI
 2,H2O,H2O+,E,@xsecFile=SWRI
@photo_end

The corresponding data file should be located in the directory /data/database/swri_xsecs/ with the name of the reactant, e.g. H2O.dat for water. Additional files can be downloaded from here.

2.5.7 Cross section from Leiden file

KROME can use the cross section from the Leiden photochemistry database at http://home.strw.leidenuniv.nl/~ewine/photo/ For example, to include water photodissociation/ionization write

@photo_start
 1,H2O,H,OH,@xsecFile=LEIDEN
 2,H2O,H2O+,E,@xsecFile=LEIDEN
@photo_end

The corresponding data file should be located in the directory /data/database/leiden_xsecs/ with the name of the reactant and the products, e.g. H2O__H_OH.dat for reaction 1 above, and H2O__H2O+_E.dat for the second. Additional files can be downloaded from here.

Note that even photoionization and photodissociation data are in the same file, you have to provide one file per branching ratio.

2.5.8 Cross sections not present in the database

For the cross sections that are not present in the database a rescaling could be used. The photo-rates in KIDA database are computed assuming a Draine flux in the range 5-13.6 eV. This means that, if you know the attenuation of your bin-based flux, you can obtain the effective rate (k') by rescaling the unattenuated rate (k) by

\begin{equation*} k' = k\, \times\,\frac{\int_{0}^\infty I(E)/E\,dE}{\int_{0}^\infty F_D(E)/E\,dE} \end{equation*}

where I(E) is your flux and F_D(E) the Draine flux. The unattenuated rate (k) is the rate when Av=0, i.e. k=alpha of KIDA rates alpha*exp(-gamma*Av). To do this modify your network file as e.g.

@var: fscale = get_ratioFluxDraine()
226,H2+,H,H+,1.1d-9*fscale
227,H3+,H2,H+,4.9d-13*fscale
228,H3+,H2+,H,4.9d-13*fscale
...

2.5.9 Impinging radiation from file

To load a given spectral energy distribution from file (e.g. named "fname") use the subroutine

call krome_load_photoBin_file("fname")

It automatically sets the interval size and the intensity for each bin. The file should be a 3-column file as

left_interval     right_interval    intensity

where left/right_interval are the intervals in eV of each bin, while intensity (eV/cm2/sr) is the radiation intensity. So for example, if you have 3 energy bins with three energy ranges your file should be something like

8.0   8.5   1.21e21
8.5   9.0   3.33e21
9.0   9.9   8.11e21

In this example the last bin is larger than the previous two (0.9 eV instead of 0.5). Note that the number of data in the file MUST be equal to the number of the photo-bins set with -photoBins. Finally, the file should be in build/ directory, and you can use the relative path to this folder.