Coding Guides

Adding a new line/coolant

Here we describe the steps to add a new coolant to ProDiMo if the LAMDA-format file is available. Some of this describption is also relevant for adding a different type of coolant, however, the details will be different.

Please also take a look at Custom Molecules/Lines, which allows to add certain types of coolants to a model without changing the code. This can be especially useful if you just want to test a new line data file.

There are two options. One can simply add a new line/molecule but without making it a cooling and/or heating process. This option will be described first, what needs to be done to have the line also as coolant, will be show later.

Adding a new line (SYS)

dataNLTE.f

you need to increase the value of N_bb_species. 37 (likely different now) is the number of the coolant HNC in the variable SYS. For a species with ortho/para, one has to increase by 2.

!----------------------------------------------------------------------
      MODULE NLTE
!----------------------------------------------------------------------
! ***   data structures for non-LTE N-Level systems                 ***
!----------------------------------------------------------------------
      use SPECIES,ONLY: SPlen
      use PARAMETERS,ONLY: H2O_rovibration
      integer,parameter :: N_bb_species=37  ! number of bound-bound coolants

init_heatcool.f

Now we can add the new line. How this is done depends on the type of the line (e.g. LAMDA format, HITRAN, ...). But in any case one has to add a new line in the subroutine setup_SYSconfs in init_heatcool.f.

LAMDA format

Here is an example for a line in the LAMDA format, see also the existing entries in the code for further examples.

    ! H13CN is the new SYS, acivated by default, HCN is the chemical 
    ! species, -1,-1 means not included as a coolant for the heating/
    ! cooling balance, the last argument is the name of the data  file 
    ! in the data directory
    call init_SYS_LAMDA(109,"H13CN",.true.,"HCN",-1,-1,
   >                      "H13CN_Lambda.dat")

Here we initialise a new SYS entry, by setting the index (109 in the example) accordingly. As a minimum one needs to provide a name, which is usually the molecule (, here H13CN or e.g o-H2 for an ortho species), and the cspname which is the according chemical species name in ProDiMo. Via the switch active one can switch on/off the line. Please note only if the chemical species given as cspname is actually part of the chemical network used in the model, the line can be active (i.e. will be automatically deactivated if the species is not in the network). Please note, in case of ortho/para lines you have to have two such entries (see the examples in the code).

The two -1 are the coolnum and heatnum, which means that this line is not included in the heating/cooling balance. If you want to include it, you need to set these numbers to the according numbers of the cooling and heating process (see below). The last argument is the name of the data file in the data directory.

If you chose to not include the line as a coolant, you are almost done, but don't forget to add the data file in the data directory and to test the code changes.

If you add a species that is an isotopologue or an ortho/para species one also needs to adapt the routine fac_nmol. Here the calculation for the actual number species is done for those special cases. Here is an example for an ortho/para species:

      if (name.eq."o-H2") then
        call H2_ORTHO_PARA(Tg,Td,fortho,fpara)
        fac_nmol = fortho
      else if (name.eq."p-H2") then
        call H2_ORTHO_PARA(Tg,Td,fortho,fpara)
        fac_nmol = fpara
      ...
Special molecules

If you read data from multiple files or different formats, please look at the code. For example water is such a case. It might then also be necessary to set the dimensions for the SYS data structure by hand. Examples are at the beginning of the init_heatcool routine.

      ! Setting of NDIMS for special coolants

      !-------------------------------------------------------------------------
      ! Atoms, because of UVpumping ... if that is false dynamic alloc might work
      SYS(1)%Ndims(:) = (/ 18,   57,    4,    7,   12/)  ! CII, C+ can be that READ_NLTE_KAMP is used
      SYS(4)%Ndims(:) = (/ 91,  647,    5,   10,   28/)  ! OI
      SYS(5)%Ndims(:) = (/ 59,  117,    6,   10,   21/)  ! CI

Make it a coolant

This is only possible if you have already added the new line/molecule to SYS as described above. To also make it a cooling and heating process a few more changes are required.

datamod.f

you need to increase the value of COOLMAX if it is not high enough already.

      module HEATCOOL
      !--------------------------------------------------------------
      ! ***  definition of variables for gas heating and cooling  ***
      !--------------------------------------------------------------
      use REACTIONS,ONLY: REACMAX
      integer,parameter :: COOLMAX=47        ! max no. of heat/cool processes

init_heatcool.f

Increase the values

      NHEAT=47
      NCOOL=42

Add two if you have ortho and para since in that case you have two SYS entries.

Add a name for the heating in init_heatcoolnames

      heatname(47) = "IR background heating by HNC"

and a name for the cooling by HNC

      coolname(42) = "HNC cooling"

Now set the according coolnum and heatnum for the SYS we have added before (in subroutine setup_SYSconf). This will automatically enable the new heating/cooling process and it will be considered in the heating/cooling balance.

write_MCFOST.f90

TODO: this is still not updated ... but likely will also not work any more

In case the transition name is not the same as the species name and if the number density of the species needs to be tweaked (e.g. isotopologue).

Adding a HITRAN molecule

TODO: one can use custom molecule or add it in the code. Similar to the others but also different (selection rules etc.).

When you are done

Test the code changes (check also if the line transfer is working, and eventually FLiTs) and commit the changed sources files and the corresponding data file (e.g. data/hnc_lambda.dat)

Electrons as collisional partner

Collision rates with electrons can be approximated within ProDiMo for linear molecules (CO, HCN, SiO, ...) or molecular ions (HCO+, ...) using the approximation proposed by Dickinson et al. (1977, Astr. Ap., 54, 645).

In init_heatcool.f, add the following lines (here for CO):

      SYS(3)%Brot = 1.922529 ! rotational constant cm^-1
      SYS(3)%Drot = 6.12107794e-6 ! centrifugal distortion cm^-1
      SYS(3)%dipole = 0.112 ! Debye
      CALL coll_rate_electrons(3,5,.false.)

Adding photo cross-sections

in data/ChemicalNetwork/pd or data/ChemicalNetwork/pd2017_version2, depending on the format of the files you have, edit the pp_list.txt file. The first line is the number of entries (Reactions) in the file. Towards the end of the file add the reaction and the file name of the cross-section. See also the README files in the corresponding directories.

pp_list.txt contains a list of photoprocesses and the associated datafile names. The data are in the Leiden photorate format.

Line 1: number of photoprocesses
Column
   1         2         3         4            5
Reactant  Product1  Product2  Filename  Branching ratio