!--------------------------------------------------------------------------------------------------
!     FILE customkinetics.f90
!>    Module for calculating the analytical source terms for the C3H8_22_173_12_FC scheme.
!!    @file    C3H8_22_173_12_FC.f90
!!    @authors F. Collin
!!    @date    22/03/2017
!!    @since   V7.0.2
!!    @note    Réduction analytique avec YARC du schéma de Jerzembeck-Pepiot "High temperature",
!!             lui-même issu d'une réduction analytique des schémas détaillés de l'iso-octane et
!!             du n-heptane du LLNL.
!!             Validation:
!!             - Validé en Sl pour T=288K, phi de 0.5 à 1.6.
!!             - Pas encore de validation à T plus élevée.
!!             - Pas encore de validation en pression.
!!             - a priori faux en auto-allumage car le schema d'origine n'est pas construit pour
!--------------------------------------------------------------------------------------------------

MODULE mC3H82217312FC

  IMPLICIT NONE

  INTEGER, PARAMETER     :: pr = SELECTED_REAL_KIND(14,300)
  REAL(KIND=pr)          :: working_precision_variables

  INTEGER, PARAMETER :: nspec = 22                                      !< Number of species (transported)
  INTEGER, PARAMETER :: nreac = 173                                     !< Number of reactions
  INTEGER, PARAMETER :: isc_T = 1
  INTEGER, PARAMETER :: neq = nspec + 1

! QSS variables
! Number of QSS species
  INTEGER, PARAMETER :: nqss = 12                                        !< Number of QSS species
! QSS species
  INTEGER, DIMENSION(nspec + nqss) :: iqss

  REAL(pr), DIMENSION(nspec + nqss) :: W_sp,Cp_sp,h_sp,dh_sp
  CHARACTER(LEN=15), DIMENSION(nspec + nqss) :: gname
  CHARACTER(LEN=5), DIMENSION(nreac) :: rname

! Post processing variables
! Number of groups
  INTEGER, PARAMETER :: ng = 35
! Max number of species in groups
  INTEGER, PARAMETER :: maxppn = 1
! Number of species in each group
  INTEGER, DIMENSION(ng) :: ppn
! Species in each group
  INTEGER, DIMENSION(ng,maxppn) :: pp
! Name of species in each group
  CHARACTER(LEN=30), DIMENSION(ng) :: ppname
! Actual expression of each reaction
  CHARACTER(LEN=65), DIMENSION(nreac) :: reacexp
! Link between backward and forward rates
  INTEGER, DIMENSION(nreac) :: fofb

! Index of species
  INTEGER, PARAMETER :: sN2 = 1
  INTEGER, PARAMETER :: sO = 2
  INTEGER, PARAMETER :: sO2 = 3
  INTEGER, PARAMETER :: sH = 4
  INTEGER, PARAMETER :: sOH = 5
  INTEGER, PARAMETER :: sH2 = 6
  INTEGER, PARAMETER :: sH2O = 7
  INTEGER, PARAMETER :: sH2O2 = 8
  INTEGER, PARAMETER :: sHO2 = 9
  INTEGER, PARAMETER :: sCO = 10
  INTEGER, PARAMETER :: sCH2O = 11
  INTEGER, PARAMETER :: sCH3 = 12
  INTEGER, PARAMETER :: sC2H6 = 13
  INTEGER, PARAMETER :: sCH4 = 14
  INTEGER, PARAMETER :: sC2H4 = 15
  INTEGER, PARAMETER :: sCO2 = 16
  INTEGER, PARAMETER :: sCH3O2 = 17
  INTEGER, PARAMETER :: sCH3O2H = 18
  INTEGER, PARAMETER :: sC2H2 = 19
  INTEGER, PARAMETER :: sC3H6 = 20
  INTEGER, PARAMETER :: sC3H5O = 21
  INTEGER, PARAMETER :: sC3H8 = 22
  INTEGER, PARAMETER :: sCH2GSGXCH2 = 23
  INTEGER, PARAMETER :: sCH3O = 24
  INTEGER, PARAMETER :: sC2H5 = 25
  INTEGER, PARAMETER :: sHCO = 26
  INTEGER, PARAMETER :: sHCCO = 27
  INTEGER, PARAMETER :: sC2H3 = 28
  INTEGER, PARAMETER :: sCH2CHO = 29
  INTEGER, PARAMETER :: sC3H5XAXC3H5 = 30
  INTEGER, PARAMETER :: sIXC3H7 = 31
  INTEGER, PARAMETER :: sNXC3H7 = 32
  INTEGER, PARAMETER :: sIXC3H7O2 = 33
  INTEGER, PARAMETER :: sNXC3H7O2 = 34

! Index of reactions
  INTEGER, PARAMETER :: rH1 = 1
  INTEGER, PARAMETER :: rH2 = 2
  INTEGER, PARAMETER :: rH3f = 3
  INTEGER, PARAMETER :: rH3b = 4
  INTEGER, PARAMETER :: rH4 = 5
  INTEGER, PARAMETER :: rH5f = 6
  INTEGER, PARAMETER :: rH6f = 7
  INTEGER, PARAMETER :: rH6b = 8
  INTEGER, PARAMETER :: rH7f = 9
  INTEGER, PARAMETER :: rH7b = 10
  INTEGER, PARAMETER :: rH8f = 11
  INTEGER, PARAMETER :: rH9f = 12
  INTEGER, PARAMETER :: rH9b = 13
  INTEGER, PARAMETER :: rH10f = 14
  INTEGER, PARAMETER :: rH10b = 15
  INTEGER, PARAMETER :: rH11f = 16
  INTEGER, PARAMETER :: rH11b = 17
  INTEGER, PARAMETER :: rH12 = 18
  INTEGER, PARAMETER :: rH13 = 19
  INTEGER, PARAMETER :: rH14 = 20
  INTEGER, PARAMETER :: rH15 = 21
  INTEGER, PARAMETER :: rH16f = 22
  INTEGER, PARAMETER :: rH16b = 23
  INTEGER, PARAMETER :: rH17 = 24
  INTEGER, PARAMETER :: rH18 = 25
  INTEGER, PARAMETER :: rH19f = 26
  INTEGER, PARAMETER :: rH19b = 27
  INTEGER, PARAMETER :: rH20 = 28
  INTEGER, PARAMETER :: rH21 = 29
  INTEGER, PARAMETER :: rH22 = 30
  INTEGER, PARAMETER :: rH23 = 31
  INTEGER, PARAMETER :: rH24f = 32
  INTEGER, PARAMETER :: rH24b = 33
  INTEGER, PARAMETER :: rH25f = 34
  INTEGER, PARAMETER :: rH26f = 35
  INTEGER, PARAMETER :: rH26b = 36
  INTEGER, PARAMETER :: rH27 = 37
  INTEGER, PARAMETER :: rH28f = 38
  INTEGER, PARAMETER :: rH28b = 39
  INTEGER, PARAMETER :: rH29 = 40
  INTEGER, PARAMETER :: rH30f = 41
  INTEGER, PARAMETER :: rH31 = 42
  INTEGER, PARAMETER :: rH32 = 43
  INTEGER, PARAMETER :: rH33 = 44
  INTEGER, PARAMETER :: rH34 = 45
  INTEGER, PARAMETER :: rH35 = 46
  INTEGER, PARAMETER :: rH36 = 47
  INTEGER, PARAMETER :: rH37f = 48
  INTEGER, PARAMETER :: rH37b = 49
  INTEGER, PARAMETER :: rH38f = 50
  INTEGER, PARAMETER :: rH38b = 51
  INTEGER, PARAMETER :: rH39f = 52
  INTEGER, PARAMETER :: rH39b = 53
  INTEGER, PARAMETER :: rH40f = 54
  INTEGER, PARAMETER :: rH40b = 55
  INTEGER, PARAMETER :: rH41 = 56
  INTEGER, PARAMETER :: rH42 = 57
  INTEGER, PARAMETER :: rH43f = 58
  INTEGER, PARAMETER :: rH43b = 59
  INTEGER, PARAMETER :: rH44f = 60
  INTEGER, PARAMETER :: rH44b = 61
  INTEGER, PARAMETER :: rH45 = 62
  INTEGER, PARAMETER :: rH46 = 63
  INTEGER, PARAMETER :: rH47 = 64
  INTEGER, PARAMETER :: rH48f = 65
  INTEGER, PARAMETER :: rH48b = 66
  INTEGER, PARAMETER :: rH49 = 67
  INTEGER, PARAMETER :: rH50 = 68
  INTEGER, PARAMETER :: rH51 = 69
  INTEGER, PARAMETER :: rH52 = 70
  INTEGER, PARAMETER :: rH53 = 71
  INTEGER, PARAMETER :: rH54 = 72
  INTEGER, PARAMETER :: rH55 = 73
  INTEGER, PARAMETER :: rH56 = 74
  INTEGER, PARAMETER :: rH57 = 75
  INTEGER, PARAMETER :: rH58 = 76
  INTEGER, PARAMETER :: rH59 = 77
  INTEGER, PARAMETER :: rH60 = 78
  INTEGER, PARAMETER :: rH61 = 79
  INTEGER, PARAMETER :: rH62 = 80
  INTEGER, PARAMETER :: rH63f = 81
  INTEGER, PARAMETER :: rH63b = 82
  INTEGER, PARAMETER :: rH64 = 83
  INTEGER, PARAMETER :: rH65 = 84
  INTEGER, PARAMETER :: rH66f = 85
  INTEGER, PARAMETER :: rH66b = 86
  INTEGER, PARAMETER :: rH67 = 87
  INTEGER, PARAMETER :: rH68 = 88
  INTEGER, PARAMETER :: rH69 = 89
  INTEGER, PARAMETER :: rH70 = 90
  INTEGER, PARAMETER :: rH71 = 91
  INTEGER, PARAMETER :: rH72 = 92
  INTEGER, PARAMETER :: rH73 = 93
  INTEGER, PARAMETER :: rH74 = 94
  INTEGER, PARAMETER :: rH75f = 95
  INTEGER, PARAMETER :: rH75b = 96
  INTEGER, PARAMETER :: rH76 = 97
  INTEGER, PARAMETER :: rH77f = 98
  INTEGER, PARAMETER :: rH77b = 99
  INTEGER, PARAMETER :: rH78f = 100
  INTEGER, PARAMETER :: rH78b = 101
  INTEGER, PARAMETER :: rH79 = 102
  INTEGER, PARAMETER :: rH80f = 103
  INTEGER, PARAMETER :: rH81 = 104
  INTEGER, PARAMETER :: rH82 = 105
  INTEGER, PARAMETER :: rH83 = 106
  INTEGER, PARAMETER :: rH84 = 107
  INTEGER, PARAMETER :: rH85f = 108
  INTEGER, PARAMETER :: rH85b = 109
  INTEGER, PARAMETER :: rH86f = 110
  INTEGER, PARAMETER :: rH86b = 111
  INTEGER, PARAMETER :: rH87 = 112
  INTEGER, PARAMETER :: rH88 = 113
  INTEGER, PARAMETER :: rH89 = 114
  INTEGER, PARAMETER :: rH90 = 115
  INTEGER, PARAMETER :: rH91 = 116
  INTEGER, PARAMETER :: rH92 = 117
  INTEGER, PARAMETER :: rH93f = 118
  INTEGER, PARAMETER :: rH93b = 119
  INTEGER, PARAMETER :: rH94 = 120
  INTEGER, PARAMETER :: rH95 = 121
  INTEGER, PARAMETER :: rH96 = 122
  INTEGER, PARAMETER :: rH97 = 123
  INTEGER, PARAMETER :: rH98 = 124
  INTEGER, PARAMETER :: rH99 = 125
  INTEGER, PARAMETER :: rH100 = 126
  INTEGER, PARAMETER :: rH101 = 127
  INTEGER, PARAMETER :: rH102 = 128
  INTEGER, PARAMETER :: rH103f = 129
  INTEGER, PARAMETER :: rH103b = 130
  INTEGER, PARAMETER :: rH104 = 131
  INTEGER, PARAMETER :: rH105 = 132
  INTEGER, PARAMETER :: rH106f = 133
  INTEGER, PARAMETER :: rH106b = 134
  INTEGER, PARAMETER :: rH107f = 135
  INTEGER, PARAMETER :: rH107b = 136
  INTEGER, PARAMETER :: rH108 = 137
  INTEGER, PARAMETER :: rH109f = 138
  INTEGER, PARAMETER :: rH109b = 139
  INTEGER, PARAMETER :: rH110 = 140
  INTEGER, PARAMETER :: rH111 = 141
  INTEGER, PARAMETER :: rH112f = 142
  INTEGER, PARAMETER :: rH112b = 143
  INTEGER, PARAMETER :: rH113 = 144
  INTEGER, PARAMETER :: rH114f = 145
  INTEGER, PARAMETER :: rH114b = 146
  INTEGER, PARAMETER :: rH115f = 147
  INTEGER, PARAMETER :: rH115b = 148
  INTEGER, PARAMETER :: rH116 = 149
  INTEGER, PARAMETER :: rH117 = 150
  INTEGER, PARAMETER :: rH118 = 151
  INTEGER, PARAMETER :: rH119 = 152
  INTEGER, PARAMETER :: rH120 = 153
  INTEGER, PARAMETER :: rH121 = 154
  INTEGER, PARAMETER :: rH122 = 155
  INTEGER, PARAMETER :: rH123 = 156
  INTEGER, PARAMETER :: rH124f = 157
  INTEGER, PARAMETER :: rH124b = 158
  INTEGER, PARAMETER :: rH125 = 159
  INTEGER, PARAMETER :: rH126f = 160
  INTEGER, PARAMETER :: rH127f = 161
  INTEGER, PARAMETER :: rH127b = 162
  INTEGER, PARAMETER :: rH128 = 163
  INTEGER, PARAMETER :: rH129f = 164
  INTEGER, PARAMETER :: rH129b = 165
  INTEGER, PARAMETER :: rH130f = 166
  INTEGER, PARAMETER :: rH130b = 167
  INTEGER, PARAMETER :: rH5b = 168
  INTEGER, PARAMETER :: rH8b = 169
  INTEGER, PARAMETER :: rH25b = 170
  INTEGER, PARAMETER :: rH30b = 171
  INTEGER, PARAMETER :: rH80b = 172
  INTEGER, PARAMETER :: rH126b = 173

! Index of third bodies
  INTEGER, PARAMETER :: mM7 = 1
  INTEGER, PARAMETER :: mM5 = 2
  INTEGER, PARAMETER :: mM2 = 3
  INTEGER, PARAMETER :: mM9 = 4
  INTEGER, PARAMETER :: mM4 = 5
  INTEGER, PARAMETER :: mM3 = 6
  INTEGER, PARAMETER :: mM1 = 7
  INTEGER, PARAMETER :: mM8 = 8
  INTEGER, PARAMETER :: mM12 = 9
  INTEGER, PARAMETER :: mM6 = 10
  INTEGER, PARAMETER :: mM11 = 11
  INTEGER, PARAMETER :: mM0 = 12

  CONTAINS

! Name of mechanism
! subroutine which_mechanism
! use parallel
!implicit none

! if (irank==iroot) print*, 'Using mechanism Spec_34_Reac_173_QSS_12.mech'

! return
! end subroutine which_mechanism

! Subroutine to define groups for post-processing
  SUBROUTINE pp_data
    IMPLICIT NONE

! Number of species in each group
    ppn(1) = 1
    ppn(2) = 1
    ppn(3) = 1
    ppn(4) = 1
    ppn(5) = 1
    ppn(6) = 1
    ppn(7) = 1
    ppn(8) = 1
    ppn(9) = 1
    ppn(10) = 1
    ppn(11) = 1
    ppn(12) = 1
    ppn(13) = 1
    ppn(14) = 1
    ppn(15) = 1
    ppn(16) = 1
    ppn(17) = 1
    ppn(18) = 1
    ppn(19) = 1
    ppn(20) = 1
    ppn(21) = 1
    ppn(22) = 1
    ppn(23) = 1
    ppn(24) = 1
    ppn(25) = 1
    ppn(26) = 1
    ppn(27) = 1
    ppn(28) = 1
    ppn(29) = 1
    ppn(30) = 1
    ppn(31) = 1
    ppn(32) = 1
    ppn(33) = 1
    ppn(34) = 1
    ppn(35) = 1

! Indices of species in each group
    pp(1,1) = 1
    pp(2,1) = 2
    pp(3,1) = 3
    pp(4,1) = 4
    pp(5,1) = 5
    pp(6,1) = 6
    pp(7,1) = 7
    pp(8,1) = 8
    pp(9,1) = 9
    pp(10,1) = 10
    pp(11,1) = 11
    pp(12,1) = 12
    pp(13,1) = 13
    pp(14,1) = 14
    pp(15,1) = 15
    pp(16,1) = 16
    pp(17,1) = 17
    pp(18,1) = 18
    pp(19,1) = 19
    pp(20,1) = 20
    pp(21,1) = 21
    pp(22,1) = 22
    pp(23,1) = 23
    pp(24,1) = 24
    pp(25,1) = 25
    pp(26,1) = 26
    pp(27,1) = 27
    pp(28,1) = 28
    pp(29,1) = 29
    pp(30,1) = 30
    pp(31,1) = 31
    pp(32,1) = 32
    pp(33,1) = 33
    pp(34,1) = 34
    pp(35,1) = 1

! Name of group of species
    ppname(1) = TRIM(gname(sN2))
    ppname(2) = TRIM(gname(sO))
    ppname(3) = TRIM(gname(sO2))
    ppname(4) = TRIM(gname(sH))
    ppname(5) = TRIM(gname(sOH))
    ppname(6) = TRIM(gname(sH2))
    ppname(7) = TRIM(gname(sH2O))
    ppname(8) = TRIM(gname(sH2O2))
    ppname(9) = TRIM(gname(sHO2))
    ppname(10) = TRIM(gname(sCH2GSGXCH2))
    ppname(11) = TRIM(gname(sCO))
    ppname(12) = TRIM(gname(sCH2O))
    ppname(13) = TRIM(gname(sCH3))
    ppname(14) = TRIM(gname(sC2H6))
    ppname(15) = TRIM(gname(sCH3O))
    ppname(16) = TRIM(gname(sCH4))
    ppname(17) = TRIM(gname(sC2H5))
    ppname(18) = TRIM(gname(sC2H4))
    ppname(19) = TRIM(gname(sCO2))
    ppname(20) = TRIM(gname(sHCO))
    ppname(21) = TRIM(gname(sCH3O2))
    ppname(22) = TRIM(gname(sCH3O2H))
    ppname(23) = TRIM(gname(sC2H2))
    ppname(24) = TRIM(gname(sHCCO))
    ppname(25) = TRIM(gname(sC2H3))
    ppname(26) = TRIM(gname(sCH2CHO))
    ppname(27) = TRIM(gname(sC3H5XAXC3H5))
    ppname(28) = TRIM(gname(sC3H6))
    ppname(29) = TRIM(gname(sC3H5O))
    ppname(30) = TRIM(gname(sIXC3H7))
    ppname(31) = TRIM(gname(sNXC3H7))
    ppname(32) = TRIM(gname(sC3H8))
    ppname(33) = TRIM(gname(sIXC3H7O2))
    ppname(34) = TRIM(gname(sNXC3H7O2))
    ppname(35) = 'N2X'

  END SUBROUTINE pp_data

! Molar mass
  SUBROUTINE molar_mass
    IMPLICIT NONE

    W_sp(sN2) = 0.02801348_pr
    W_sp(sO) = 0.0159994_pr
    W_sp(sO2) = 0.0319988_pr
    W_sp(sH) = 0.00100794_pr
    W_sp(sOH) = 0.01700734_pr
    W_sp(sH2) = 0.00201588_pr
    W_sp(sH2O) = 0.01801528_pr
    W_sp(sH2O2) = 0.03401468_pr
    W_sp(sHO2) = 0.03300674_pr
    W_sp(sCH2GSGXCH2) = 0.01402688_pr
    W_sp(sCO) = 0.0280104_pr
    W_sp(sCH2O) = 0.03002628_pr
    W_sp(sCH3) = 0.01503482_pr
    W_sp(sC2H6) = 0.03006964_pr
    W_sp(sCH3O) = 0.03103422_pr
    W_sp(sCH4) = 0.01604276_pr
    W_sp(sC2H5) = 0.0290617_pr
    W_sp(sC2H4) = 0.02805376_pr
    W_sp(sCO2) = 0.0440098_pr
    W_sp(sHCO) = 0.02901834_pr
    W_sp(sCH3O2) = 0.04703362_pr
    W_sp(sCH3O2H) = 0.04804156_pr
    W_sp(sC2H2) = 0.02603788_pr
    W_sp(sHCCO) = 0.04102934_pr
    W_sp(sC2H3) = 0.02704582_pr
    W_sp(sCH2CHO) = 0.04304522_pr
    W_sp(sC3H5XAXC3H5) = 0.0410727_pr
    W_sp(sC3H6) = 0.04208064_pr
    W_sp(sC3H5O) = 0.0570721_pr
    W_sp(sIXC3H7) = 0.04308858_pr
    W_sp(sNXC3H7) = 0.04308858_pr
    W_sp(sC3H8) = 0.04409652_pr
    W_sp(sIXC3H7O2) = 0.07508738_pr
    W_sp(sNXC3H7O2) = 0.07508738_pr

  END SUBROUTINE molar_mass

! Species names
  SUBROUTINE species_name
    IMPLICIT NONE

    gname(sN2) = 'N2'
    gname(sO) = 'O'
    gname(sO2) = 'O2'
    gname(sH) = 'H'
    gname(sOH) = 'OH'
    gname(sH2) = 'H2'
    gname(sH2O) = 'H2O'
    gname(sH2O2) = 'H2O2'
    gname(sHO2) = 'HO2'
    gname(sCH2GSGXCH2) = 'CH2GSG-CH2'
    gname(sCO) = 'CO'
    gname(sCH2O) = 'CH2O'
    gname(sCH3) = 'CH3'
    gname(sC2H6) = 'C2H6'
    gname(sCH3O) = 'CH3O'
    gname(sCH4) = 'CH4'
    gname(sC2H5) = 'C2H5'
    gname(sC2H4) = 'C2H4'
    gname(sCO2) = 'CO2'
    gname(sHCO) = 'HCO'
    gname(sCH3O2) = 'CH3O2'
    gname(sCH3O2H) = 'CH3O2H'
    gname(sC2H2) = 'C2H2'
    gname(sHCCO) = 'HCCO'
    gname(sC2H3) = 'C2H3'
    gname(sCH2CHO) = 'CH2CHO'
    gname(sC3H5XAXC3H5) = 'C3H5-A-C3H5'
    gname(sC3H6) = 'C3H6'
    gname(sC3H5O) = 'C3H5O'
    gname(sIXC3H7) = 'I-C3H7'
    gname(sNXC3H7) = 'N-C3H7'
    gname(sC3H8) = 'C3H8'
    gname(sIXC3H7O2) = 'I-C3H7O2'
    gname(sNXC3H7O2) = 'N-C3H7O2'

  END SUBROUTINE species_name

! Reaction names
  SUBROUTINE reaction_name
    IMPLICIT NONE

    rname(rH1) = 'H1'
    rname(rH2) = 'H2'
    rname(rH3f) = 'H3f'
    rname(rH3b) = 'H3b'
    rname(rH4) = 'H4'
    rname(rH5f) = 'H5f'
    rname(rH6f) = 'H6f'
    rname(rH6b) = 'H6b'
    rname(rH7f) = 'H7f'
    rname(rH7b) = 'H7b'
    rname(rH8f) = 'H8f'
    rname(rH9f) = 'H9f'
    rname(rH9b) = 'H9b'
    rname(rH10f) = 'H10f'
    rname(rH10b) = 'H10b'
    rname(rH11f) = 'H11f'
    rname(rH11b) = 'H11b'
    rname(rH12) = 'H12'
    rname(rH13) = 'H13'
    rname(rH14) = 'H14'
    rname(rH15) = 'H15'
    rname(rH16f) = 'H16f'
    rname(rH16b) = 'H16b'
    rname(rH17) = 'H17'
    rname(rH18) = 'H18'
    rname(rH19f) = 'H19f'
    rname(rH19b) = 'H19b'
    rname(rH20) = 'H20'
    rname(rH21) = 'H21'
    rname(rH22) = 'H22'
    rname(rH23) = 'H23'
    rname(rH24f) = 'H24f'
    rname(rH24b) = 'H24b'
    rname(rH25f) = 'H25f'
    rname(rH26f) = 'H26f'
    rname(rH26b) = 'H26b'
    rname(rH27) = 'H27'
    rname(rH28f) = 'H28f'
    rname(rH28b) = 'H28b'
    rname(rH29) = 'H29'
    rname(rH30f) = 'H30f'
    rname(rH31) = 'H31'
    rname(rH32) = 'H32'
    rname(rH33) = 'H33'
    rname(rH34) = 'H34'
    rname(rH35) = 'H35'
    rname(rH36) = 'H36'
    rname(rH37f) = 'H37f'
    rname(rH37b) = 'H37b'
    rname(rH38f) = 'H38f'
    rname(rH38b) = 'H38b'
    rname(rH39f) = 'H39f'
    rname(rH39b) = 'H39b'
    rname(rH40f) = 'H40f'
    rname(rH40b) = 'H40b'
    rname(rH41) = 'H41'
    rname(rH42) = 'H42'
    rname(rH43f) = 'H43f'
    rname(rH43b) = 'H43b'
    rname(rH44f) = 'H44f'
    rname(rH44b) = 'H44b'
    rname(rH45) = 'H45'
    rname(rH46) = 'H46'
    rname(rH47) = 'H47'
    rname(rH48f) = 'H48f'
    rname(rH48b) = 'H48b'
    rname(rH49) = 'H49'
    rname(rH50) = 'H50'
    rname(rH51) = 'H51'
    rname(rH52) = 'H52'
    rname(rH53) = 'H53'
    rname(rH54) = 'H54'
    rname(rH55) = 'H55'
    rname(rH56) = 'H56'
    rname(rH57) = 'H57'
    rname(rH58) = 'H58'
    rname(rH59) = 'H59'
    rname(rH60) = 'H60'
    rname(rH61) = 'H61'
    rname(rH62) = 'H62'
    rname(rH63f) = 'H63f'
    rname(rH63b) = 'H63b'
    rname(rH64) = 'H64'
    rname(rH65) = 'H65'
    rname(rH66f) = 'H66f'
    rname(rH66b) = 'H66b'
    rname(rH67) = 'H67'
    rname(rH68) = 'H68'
    rname(rH69) = 'H69'
    rname(rH70) = 'H70'
    rname(rH71) = 'H71'
    rname(rH72) = 'H72'
    rname(rH73) = 'H73'
    rname(rH74) = 'H74'
    rname(rH75f) = 'H75f'
    rname(rH75b) = 'H75b'
    rname(rH76) = 'H76'
    rname(rH77f) = 'H77f'
    rname(rH77b) = 'H77b'
    rname(rH78f) = 'H78f'
    rname(rH78b) = 'H78b'
    rname(rH79) = 'H79'
    rname(rH80f) = 'H80f'
    rname(rH81) = 'H81'
    rname(rH82) = 'H82'
    rname(rH83) = 'H83'
    rname(rH84) = 'H84'
    rname(rH85f) = 'H85f'
    rname(rH85b) = 'H85b'
    rname(rH86f) = 'H86f'
    rname(rH86b) = 'H86b'
    rname(rH87) = 'H87'
    rname(rH88) = 'H88'
    rname(rH89) = 'H89'
    rname(rH90) = 'H90'
    rname(rH91) = 'H91'
    rname(rH92) = 'H92'
    rname(rH93f) = 'H93f'
    rname(rH93b) = 'H93b'
    rname(rH94) = 'H94'
    rname(rH95) = 'H95'
    rname(rH96) = 'H96'
    rname(rH97) = 'H97'
    rname(rH98) = 'H98'
    rname(rH99) = 'H99'
    rname(rH100) = 'H100'
    rname(rH101) = 'H101'
    rname(rH102) = 'H102'
    rname(rH103f) = 'H103f'
    rname(rH103b) = 'H103b'
    rname(rH104) = 'H104'
    rname(rH105) = 'H105'
    rname(rH106f) = 'H106f'
    rname(rH106b) = 'H106b'
    rname(rH107f) = 'H107f'
    rname(rH107b) = 'H107b'
    rname(rH108) = 'H108'
    rname(rH109f) = 'H109f'
    rname(rH109b) = 'H109b'
    rname(rH110) = 'H110'
    rname(rH111) = 'H111'
    rname(rH112f) = 'H112f'
    rname(rH112b) = 'H112b'
    rname(rH113) = 'H113'
    rname(rH114f) = 'H114f'
    rname(rH114b) = 'H114b'
    rname(rH115f) = 'H115f'
    rname(rH115b) = 'H115b'
    rname(rH116) = 'H116'
    rname(rH117) = 'H117'
    rname(rH118) = 'H118'
    rname(rH119) = 'H119'
    rname(rH120) = 'H120'
    rname(rH121) = 'H121'
    rname(rH122) = 'H122'
    rname(rH123) = 'H123'
    rname(rH124f) = 'H124f'
    rname(rH124b) = 'H124b'
    rname(rH125) = 'H125'
    rname(rH126f) = 'H126f'
    rname(rH127f) = 'H127f'
    rname(rH127b) = 'H127b'
    rname(rH128) = 'H128'
    rname(rH129f) = 'H129f'
    rname(rH129b) = 'H129b'
    rname(rH130f) = 'H130f'
    rname(rH130b) = 'H130b'
    rname(rH5b) = 'H5b'
    rname(rH8b) = 'H8b'
    rname(rH25b) = 'H25b'
    rname(rH30b) = 'H30b'
    rname(rH80b) = 'H80b'
    rname(rH126b) = 'H126b'

  END SUBROUTINE reaction_name

! List of QSS species
  SUBROUTINE QSS_list
    IMPLICIT NONE

    iqss(sN2) = 0
    iqss(sO) = 0
    iqss(sO2) = 0
    iqss(sH) = 0
    iqss(sOH) = 0
    iqss(sH2) = 0
    iqss(sH2O) = 0
    iqss(sH2O2) = 0
    iqss(sHO2) = 0
    iqss(sCH2GSGXCH2) = 1
    iqss(sCO) = 0
    iqss(sCH2O) = 0
    iqss(sCH3) = 0
    iqss(sC2H6) = 0
    iqss(sCH3O) = 1
    iqss(sCH4) = 0
    iqss(sC2H5) = 1
    iqss(sC2H4) = 0
    iqss(sCO2) = 0
    iqss(sHCO) = 1
    iqss(sCH3O2) = 0
    iqss(sCH3O2H) = 0
    iqss(sC2H2) = 0
    iqss(sHCCO) = 1
    iqss(sC2H3) = 1
    iqss(sCH2CHO) = 1
    iqss(sC3H5XAXC3H5) = 1
    iqss(sC3H6) = 0
    iqss(sC3H5O) = 0
    iqss(sIXC3H7) = 1
    iqss(sNXC3H7) = 1
    iqss(sC3H8) = 0
    iqss(sIXC3H7O2) = 1
    iqss(sNXC3H7O2) = 1

  END SUBROUTINE QSS_list

! Subroutine for pressure dependent coefficients
  REAL(pr) FUNCTION getlindratecoeff (Tloc,k0,kinf,fc,concin,Ploc)
    IMPLICIT NONE

    REAL(pr) ::  Tloc,k0,kinf,fc,Ploc
    REAL(pr), PARAMETER :: R = 8.31434_pr
    REAL(pr) :: ntmp,ccoeff,dcoeff,lgknull
    REAL(pr) :: f
    REAL(pr) :: conc, concin

    IF (concin>0.0_pr) THEN
      conc = concin
    ELSE
      conc = Ploc / ( R * Tloc )
    END IF
    ntmp = 0.75_pr - 1.27_pr * dlog10( fc )

    ccoeff = - 0.4_pr - 0.67_pr * dlog10( fc )
    dcoeff = 0.14_pr
    k0 = k0 * conc / MAX(kinf, 1.0e-60_pr)
    lgknull = dlog10(k0)
    f = (lgknull+ccoeff)/(ntmp-dcoeff*(lgknull+ccoeff))
    f = fc**(1.0_pr / ( f * f + 1.0_pr ))
    getlindratecoeff = kinf * f * k0 / ( 1.0_pr + k0 )

  END FUNCTION getlindratecoeff

  SUBROUTINE YtoC(c,P,T,Y)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nspec) :: c,Y
    REAL(pr) :: bla, P, T
    INTEGER :: K

    CALL molar_mass

    c(sN2) = Y(sN2) / W_sp(sN2)
    c(sO) = Y(sO) / W_sp(sO)
    c(sO2) = Y(sO2) / W_sp(sO2)
    c(sH) = Y(sH) / W_sp(sH)
    c(sOH) = Y(sOH) / W_sp(sOH)
    c(sH2) = Y(sH2) / W_sp(sH2)
    c(sH2O) = Y(sH2O) / W_sp(sH2O)
    c(sH2O2) = Y(sH2O2) / W_sp(sH2O2)
    c(sHO2) = Y(sHO2) / W_sp(sHO2)
!c(sCH2GSGXCH2) = Y(sCH2GSGXCH2) / W_sp(sCH2GSGXCH2)
    c(sCO) = Y(sCO) / W_sp(sCO)
    c(sCH2O) = Y(sCH2O) / W_sp(sCH2O)
    c(sCH3) = Y(sCH3) / W_sp(sCH3)
    c(sC2H6) = Y(sC2H6) / W_sp(sC2H6)
!c(sCH3O) = Y(sCH3O) / W_sp(sCH3O)
    c(sCH4) = Y(sCH4) / W_sp(sCH4)
!c(sC2H5) = Y(sC2H5) / W_sp(sC2H5)
    c(sC2H4) = Y(sC2H4) / W_sp(sC2H4)
    c(sCO2) = Y(sCO2) / W_sp(sCO2)
!c(sHCO) = Y(sHCO) / W_sp(sHCO)
    c(sCH3O2) = Y(sCH3O2) / W_sp(sCH3O2)
    c(sCH3O2H) = Y(sCH3O2H) / W_sp(sCH3O2H)
    c(sC2H2) = Y(sC2H2) / W_sp(sC2H2)
!c(sHCCO) = Y(sHCCO) / W_sp(sHCCO)
!c(sC2H3) = Y(sC2H3) / W_sp(sC2H3)
!c(sCH2CHO) = Y(sCH2CHO) / W_sp(sCH2CHO)
!c(sC3H5XAXC3H5) = Y(sC3H5XAXC3H5) / W_sp(sC3H5XAXC3H5)
    c(sC3H6) = Y(sC3H6) / W_sp(sC3H6)
    c(sC3H5O) = Y(sC3H5O) / W_sp(sC3H5O)
!c(sIXC3H7) = Y(sIXC3H7) / W_sp(sIXC3H7)
!c(sNXC3H7) = Y(sNXC3H7) / W_sp(sNXC3H7)
    c(sC3H8) = Y(sC3H8) / W_sp(sC3H8)
!c(sIXC3H7O2) = Y(sIXC3H7O2) / W_sp(sIXC3H7O2)
!c(sNXC3H7O2) = Y(sNXC3H7O2) / W_sp(sNXC3H7O2)

    bla = 0.0_pr
    DO K = 1, nspec
      bla  = bla + c(K)
    END DO
    bla = P/(bla*T*8.31451_pr)
    DO K = 1, nspec
      c(K) = MAX(c(K),1.0e-60_pr) * bla
    END DO
  END SUBROUTINE YtoC

! --- Thirdbodies --- !
  SUBROUTINE get_thirdbodies(M,c)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nspec) :: c
    REAL(pr), DIMENSION(13) :: M

    M(mM7) = (1_pr)*c(sH2)&
                           + (4_pr)*c(sH2O)&
                           + (2_pr)*c(sCO2)&
                           + (1_pr)*c(sCO)&
                           + sum(c)
    M(mM5) = sum(c)
    M(mM2) = (1.5_pr)*c(sH2)&
                           + (2.8_pr)*c(sCO2)&
                           + (0.9_pr)*c(sCO)&
                           + (11_pr)*c(sH2O)&
                           + sum(c)
    M(mM9) = (5_pr)*c(sH2O)&
                           + (2.8_pr)*c(sCO2)&
                           + (0.9_pr)*c(sCO)&
                           + (1.5_pr)*c(sH2)&
                           + sum(c)
    M(mM4) = (1.5_pr)*c(sH2)&
                           + (11_pr)*c(sH2O)&
                           + (0.9_pr)*c(sCO)&
                           + (2.8_pr)*c(sCO2)&
                           + sum(c)
    M(mM3) = (1.5_pr)*c(sH2)&
                           + (2.8_pr)*c(sCO2)&
                           + (0.9_pr)*c(sCO)&
                           + (11_pr)*c(sH2O)&
                           + sum(c)
    M(mM1) = (1.5_pr)*c(sH2)&
                           + (11_pr)*c(sH2O)&
                           + (0.9_pr)*c(sCO)&
                           + (2.8_pr)*c(sCO2)&
                           + sum(c)
    M(mM8) = (1.5_pr)*c(sH2)&
                           + (2.8_pr)*c(sCO2)&
                           + (0.9_pr)*c(sCO)&
                           + (11_pr)*c(sH2O)&
                           + sum(c)
    M(mM12) = sum(c)
    M(mM6) = (1_pr)*c(sCO)&
                           + (2_pr)*c(sCO2)&
                           + (4_pr)*c(sH2O)&
                           + (1_pr)*c(sH2)&
                           + sum(c)
    M(mM11) = (1_pr)*c(sH2)&
                           + (4_pr)*c(sH2O)&
                           + (2_pr)*c(sCO2)&
                           + (1_pr)*c(sCO)&
                           + sum(c)
    M(mM0) = (11_pr)*c(sH2O)&
                           + (2.8_pr)*c(sCO2)&
                           + (0.9_pr)*c(sCO)&
                           + (1.5_pr)*c(sH2)&
                           + sum(c)

  END SUBROUTINE get_thirdbodies

! --- Rate coefficients --- !
  SUBROUTINE get_rate_coefficients(k,M,Tloc,Ploc)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nreac) :: k
    REAL(pr), DIMENSION(13) :: M
    REAL(pr) :: Tloc,Ploc

    REAL(pr) :: kH5f_0, kH5f_inf, FCH5f
    REAL(pr) :: kH8f_0, kH8f_inf, FCH8f
    REAL(pr) :: kH25f_0, kH25f_inf, FCH25f
    REAL(pr) :: kH29_0, kH29_inf, FCH29
    REAL(pr) :: kH42_0, kH42_inf, FCH42
    REAL(pr) :: kH59_0, kH59_inf, FCH59
    REAL(pr) :: kH73_0, kH73_inf, FCH73
    REAL(pr) :: kH74_0, kH74_inf, FCH74
    REAL(pr) :: kH80f_0, kH80f_inf, FCH80f
    REAL(pr) :: kH5b_0, kH5b_inf, FCH5b
    REAL(pr) :: kH8b_0, kH8b_inf, FCH8b
    REAL(pr) :: kH25b_0, kH25b_inf, FCH25b
    REAL(pr) :: kH80b_0, kH80b_inf, FCH80b

! Rate coefficients

    k(rH1) = (6.17000000e+03_pr)*Tloc**(-0.500_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH2) = (4.72000000e+06_pr)*Tloc**(-1.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH3f) = (5.08000000e-02_pr)*Tloc**(2.670_pr)*&
                    exp(-(2.633e+04_pr)/(8.314_pr*Tloc))

    k(rH3b) = (2.23100000e-02_pr)*Tloc**(2.670_pr)*&
                    exp(-(1.756e+04_pr)/(8.314_pr*Tloc))

    k(rH4) = (2.25000000e+10_pr)*Tloc**(-2.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    kH5f_0 = (3.04100000e+18_pr)*Tloc**(-4.630_pr)*&
                    exp(-(8.573e+03_pr)/(8.314_pr*Tloc))

    kH5f_inf = (1.23600000e+08_pr)*Tloc**(-0.370_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    FCH5f = (5.300e-01_pr)*&
                    exp(-Tloc/(1.000e+02_pr)) + (4.700e-01_pr)*&
                    exp(-Tloc/(2.000e+03_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+15_pr)/Tloc)

    k(rH5f) =&
                    getlindratecoeff&
                    (Tloc, kH5f_0, kH5f_inf, FCH5f, M(mM2), Ploc )

    k(rH6f) = (2.16000000e+02_pr)*Tloc**(1.510_pr)*&
                    exp(-(1.435e+04_pr)/(8.314_pr*Tloc))

    k(rH6b) = (9.35200000e+02_pr)*Tloc**(1.510_pr)*&
                    exp(-(7.774e+04_pr)/(8.314_pr*Tloc))

    k(rH7f) = (2.97000000e+00_pr)*Tloc**(2.020_pr)*&
                    exp(-(5.607e+04_pr)/(8.314_pr*Tloc))

    k(rH7b) = (3.01300000e-01_pr)*Tloc**(2.020_pr)*&
                    exp(-(-1.611e+04_pr)/(8.314_pr*Tloc))

    kH8f_0 = (3.50000000e+04_pr)*Tloc**(-0.410_pr)*&
                    exp(-(-4.669e+03_pr)/(8.314_pr*Tloc))

    kH8f_inf = (1.47500000e+06_pr)*Tloc**(0.600_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    FCH8f = (5.000e-01_pr)*&
                    exp(-Tloc/(1.000e-30_pr)) + (5.000e-01_pr)*&
                    exp(-Tloc/(1.000e+30_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+100_pr)/Tloc)

    k(rH8f) =&
                    getlindratecoeff&
                    (Tloc, kH8f_0, kH8f_inf, FCH8f, M(mM4), Ploc )

    k(rH9f) = (1.97000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(6.920e+04_pr)/(8.314_pr*Tloc))

    k(rH9b) = (1.55500000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.778e+03_pr)/(8.314_pr*Tloc))

    k(rH10f) = (2.89000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(-2.092e+03_pr)/(8.314_pr*Tloc))

    k(rH10b) = (6.88800000e+09_pr)*Tloc**(-0.330_pr)*&
                    exp(-(3.018e+05_pr)/(8.314_pr*Tloc))

    k(rH11f) = (3.25000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH11b) = (7.85700000e+08_pr)*Tloc**(-0.330_pr)*&
                    exp(-(2.318e+05_pr)/(8.314_pr*Tloc))

    k(rH12) = (1.66000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.431e+03_pr)/(8.314_pr*Tloc))

    k(rH13) = (1.30000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(-6.816e+03_pr)/(8.314_pr*Tloc))

    k(rH14) = (7.08000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.255e+03_pr)/(8.314_pr*Tloc))

    k(rH15) = (4.20000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(5.012e+04_pr)/(8.314_pr*Tloc))

    k(rH16f) = (4.82000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.326e+04_pr)/(8.314_pr*Tloc))

    k(rH16b) = (1.87500000e+06_pr)*Tloc**(0.330_pr)*&
                    exp(-(1.015e+05_pr)/(8.314_pr*Tloc))

    k(rH17) = (9.55000000e+00_pr)*Tloc**(2.000_pr)*&
                    exp(-(1.661e+04_pr)/(8.314_pr*Tloc))

    k(rH18) = (1.00000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH19f) = (5.80000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(4.000e+04_pr)/(8.314_pr*Tloc))

    k(rH19b) = (9.77100000e+07_pr)*Tloc**(0.330_pr)*&
                    exp(-(1.716e+05_pr)/(8.314_pr*Tloc))

    k(rH20) = (2.41000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.661e+04_pr)/(8.314_pr*Tloc))

    k(rH21) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH22) = (7.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH23) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH24f) = (7.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH24b) = (2.48200000e+11_pr)*Tloc**(-0.890_pr)*&
                    exp(-(6.749e+04_pr)/(8.314_pr*Tloc))

    kH25f_0 = (1.13500000e+24_pr)*Tloc**(-5.246_pr)*&
                    exp(-(7.134e+03_pr)/(8.314_pr*Tloc))

    kH25f_inf = (9.21400000e+10_pr)*Tloc**(-1.170_pr)*&
                    exp(-(2.660e+03_pr)/(8.314_pr*Tloc))

    FCH25f = (5.950e-01_pr)*&
                    exp(-Tloc/(1.120e+03_pr)) + (4.050e-01_pr)*&
                    exp(-Tloc/(6.960e+01_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+15_pr)/Tloc)

    k(rH25f) =&
                    getlindratecoeff&
                    (Tloc, kH25f_0, kH25f_inf, FCH25f, M(mM6), Ploc )

    k(rH26f) = (1.99500000e+12_pr)*Tloc**(-1.570_pr)*&
                    exp(-(1.222e+05_pr)/(8.314_pr*Tloc))

    k(rH26b) = (3.58500000e+12_pr)*Tloc**(-1.590_pr)*&
                    exp(-(-6.824e+03_pr)/(8.314_pr*Tloc))

    k(rH27) = (1.10000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH28f) = (2.65000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(9.146e+03_pr)/(8.314_pr*Tloc))

    k(rH28b) = (3.23600000e+04_pr)*Tloc**(0.890_pr)*&
                    exp(-(5.067e+03_pr)/(8.314_pr*Tloc))

    kH29_0 = (3.31000000e+18_pr)*Tloc**(-4.000_pr)*&
                    exp(-(8.820e+03_pr)/(8.314_pr*Tloc))

    kH29_inf = (2.13800000e+09_pr)*Tloc**(-0.400_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    FCH29 = (1.000e+00_pr)*&
                    exp(-Tloc/(1.000e-15_pr)) + (0.000e+00_pr)*&
                    exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)*&
                    exp(-(4.000e+01_pr)/Tloc)

    k(rH29) =&
                    getlindratecoeff&
                    (Tloc, kH29_0, kH29_inf, FCH29, M(mM7), Ploc )

    k(rH30f) = (6.84000000e+06_pr)*Tloc**(0.100_pr)*&
                    exp(-(4.435e+04_pr)/(8.314_pr*Tloc))

    k(rH31) = (7.47000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(5.962e+04_pr)/(8.314_pr*Tloc))

    k(rH32) = (8.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH33) = (3.60000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH34) = (2.25000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.799e+04_pr)/(8.314_pr*Tloc))

    k(rH35) = (2.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH36) = (3.36500000e+05_pr)*Tloc**(-0.330_pr)*&
                    exp(-(1.047e+04_pr)/(8.314_pr*Tloc))

    k(rH37f) = (1.93000000e-01_pr)*Tloc**(2.400_pr)*&
                    exp(-(8.812e+03_pr)/(8.314_pr*Tloc))

    k(rH37b) = (3.19900000e-02_pr)*Tloc**(2.400_pr)*&
                    exp(-(7.021e+04_pr)/(8.314_pr*Tloc))

    k(rH38f) = (1.72700000e-02_pr)*Tloc**(3.000_pr)*&
                    exp(-(3.441e+04_pr)/(8.314_pr*Tloc))

    k(rH38b) = (6.61000000e-04_pr)*Tloc**(3.000_pr)*&
                    exp(-(3.240e+04_pr)/(8.314_pr*Tloc))

    k(rH39f) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH39b) = (5.42900000e+09_pr)*Tloc**(-0.890_pr)*&
                    exp(-(6.548e+04_pr)/(8.314_pr*Tloc))

    k(rH40f) = (3.15000000e+06_pr)*Tloc**(0.500_pr)*&
                    exp(-(4.305e+04_pr)/(8.314_pr*Tloc))

    k(rH40b) = (5.29600000e+04_pr)*Tloc**(0.500_pr)*&
                    exp(-(3.228e+04_pr)/(8.314_pr*Tloc))

    k(rH41) = (3.01000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(9.623e+04_pr)/(8.314_pr*Tloc))

    kH42_0 = (1.35000000e+12_pr)*Tloc**(-2.788_pr)*&
                    exp(-(1.754e+04_pr)/(8.314_pr*Tloc))

    kH42_inf = (1.80000000e+04_pr)*Tloc**(0.000_pr)*&
                    exp(-(9.975e+03_pr)/(8.314_pr*Tloc))

    FCH42 =  + (1.000e+00_pr)*&
                    exp(-(0.000e+00_pr)/Tloc)

    k(rH42) =&
                    getlindratecoeff&
                    (Tloc, kH42_0, kH42_inf, FCH42, M(mM8), Ploc )

    k(rH43f) = (1.40000000e-01_pr)*Tloc**(1.950_pr)*&
                    exp(-(-5.636e+03_pr)/(8.314_pr*Tloc))

    k(rH43b) = (1.56800000e+01_pr)*Tloc**(1.950_pr)*&
                    exp(-(8.782e+04_pr)/(8.314_pr*Tloc))

    k(rH44f) = (1.62000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.996e+05_pr)/(8.314_pr*Tloc))

    k(rH44b) = (1.43300000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(2.256e+05_pr)/(8.314_pr*Tloc))

    k(rH45) = (3.02000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH46) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH47) = (7.34000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH48f) = (1.86000000e+11_pr)*Tloc**(-1.000_pr)*&
                    exp(-(7.113e+04_pr)/(8.314_pr*Tloc))

    k(rH48b) = (6.46700000e+01_pr)*Tloc**(0.000_pr)*&
                    exp(-(-1.849e+03_pr)/(8.314_pr*Tloc))

    k(rH49) = (1.02000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH50) = (1.21000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH51) = (7.58000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.715e+03_pr)/(8.314_pr*Tloc))

    k(rH52) = (5.82000000e-09_pr)*Tloc**(4.530_pr)*&
                    exp(-(2.743e+04_pr)/(8.314_pr*Tloc))

    k(rH53) = (3.43000000e+03_pr)*Tloc**(1.180_pr)*&
                    exp(-(-1.870e+03_pr)/(8.314_pr*Tloc))

    k(rH54) = (9.33400000e+02_pr)*Tloc**(1.500_pr)*&
                    exp(-(1.245e+04_pr)/(8.314_pr*Tloc))

    k(rH55) = (3.63600000e-12_pr)*Tloc**(5.420_pr)*&
                    exp(-(4.176e+03_pr)/(8.314_pr*Tloc))

    k(rH56) = (4.16000000e+05_pr)*Tloc**(0.570_pr)*&
                    exp(-(1.156e+04_pr)/(8.314_pr*Tloc))

    k(rH57) = (2.05000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.630e+05_pr)/(8.314_pr*Tloc))

    k(rH58) = (5.50000000e+04_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.014e+04_pr)/(8.314_pr*Tloc))

    kH59_0 = (2.34400000e+19_pr)*Tloc**(-2.700_pr)*&
                    exp(-(1.280e+05_pr)/(8.314_pr*Tloc))

    kH59_inf = (5.45000000e+13_pr)*Tloc**(0.000_pr)*&
                    exp(-(5.648e+04_pr)/(8.314_pr*Tloc))

    FCH59 =  + (1.000e+00_pr)*&
                    exp(-(0.000e+00_pr)/Tloc)

    k(rH59) =&
                    getlindratecoeff&
                    (Tloc, kH59_0, kH59_inf, FCH59, M(mM5), Ploc )

    k(rH60) = (3.00000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH61) = (1.75000000e+04_pr)*Tloc**(0.000_pr)*&
                    exp(-(-1.370e+04_pr)/(8.314_pr*Tloc))

    k(rH62) = (1.40000000e+10_pr)*Tloc**(-1.610_pr)*&
                    exp(-(7.782e+03_pr)/(8.314_pr*Tloc))

    k(rH63f) = (4.34300000e+21_pr)*Tloc**(-3.420_pr)*&
                    exp(-(1.275e+05_pr)/(8.314_pr*Tloc))

    k(rH63b) = (5.44000000e+13_pr)*Tloc**(-3.300_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH64) = (1.99000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(4.878e+04_pr)/(8.314_pr*Tloc))

    k(rH65) = (7.00000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(-4.184e+03_pr)/(8.314_pr*Tloc))

    k(rH66f) = (6.31000000e+14_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.770e+05_pr)/(8.314_pr*Tloc))

    k(rH66b) = (1.16600000e+05_pr)*Tloc**(0.600_pr)*&
                    exp(-(-7.410e+03_pr)/(8.314_pr*Tloc))

    k(rH67) = (1.43000000e+01_pr)*Tloc**(2.000_pr)*&
                    exp(-(7.950e+03_pr)/(8.314_pr*Tloc))

    k(rH68) = (2.00000000e+02_pr)*Tloc**(1.500_pr)*&
                    exp(-(1.259e+05_pr)/(8.314_pr*Tloc))

    k(rH69) = (2.12000000e-12_pr)*Tloc**(6.000_pr)*&
                    exp(-(3.968e+04_pr)/(8.314_pr*Tloc))

    k(rH70) = (3.50000000e+08_pr)*Tloc**(-0.610_pr)*&
                    exp(-(2.201e+04_pr)/(8.314_pr*Tloc))

    k(rH71) = (1.70000000e+23_pr)*Tloc**(-5.310_pr)*&
                    exp(-(2.720e+04_pr)/(8.314_pr*Tloc))

    k(rH72) = (2.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.046e+04_pr)/(8.314_pr*Tloc))

    kH73_0 = (1.16400000e+33_pr)*Tloc**(-6.821_pr)*&
                    exp(-(1.862e+05_pr)/(8.314_pr*Tloc))

    kH73_inf = (1.60600000e+10_pr)*Tloc**(1.028_pr)*&
                    exp(-(1.695e+05_pr)/(8.314_pr*Tloc))

    FCH73 = (0.000e+00_pr)*&
                    exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)*&
                    exp(-Tloc/(6.750e+02_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+15_pr)/Tloc)

    k(rH73) =&
                    getlindratecoeff&
                    (Tloc, kH73_0, kH73_inf, FCH73, M(mM11), Ploc )

    kH74_0 = (1.50000000e+09_pr)*Tloc**(0.000_pr)*&
                    exp(-(2.320e+05_pr)/(8.314_pr*Tloc))

    kH74_inf = (1.80000000e+13_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.180e+05_pr)/(8.314_pr*Tloc))

    FCH74 =  + (1.000e+00_pr)*&
                    exp(-(0.000e+00_pr)/Tloc)

    k(rH74) =&
                    getlindratecoeff&
                    (Tloc, kH74_0, kH74_inf, FCH74, M(mM5), Ploc )

    k(rH75f) = (8.42000000e-09_pr)*Tloc**(4.620_pr)*&
                    exp(-(1.081e+04_pr)/(8.314_pr*Tloc))

    k(rH75b) = (5.72300000e-07_pr)*Tloc**(3.790_pr)*&
                    exp(-(1.353e+04_pr)/(8.314_pr*Tloc))

    k(rH76) = (2.23000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(7.192e+04_pr)/(8.314_pr*Tloc))

    k(rH77f) = (2.05000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(2.490e+04_pr)/(8.314_pr*Tloc))

    k(rH77b) = (6.03300000e+09_pr)*Tloc**(-0.830_pr)*&
                    exp(-(9.104e+04_pr)/(8.314_pr*Tloc))

    k(rH78f) = (6.62000000e-06_pr)*Tloc**(3.700_pr)*&
                    exp(-(3.975e+04_pr)/(8.314_pr*Tloc))

    k(rH78b) = (1.44000000e-06_pr)*Tloc**(4.020_pr)*&
                    exp(-(2.290e+04_pr)/(8.314_pr*Tloc))

    k(rH79) = (1.02000000e+01_pr)*Tloc**(1.880_pr)*&
                    exp(-(7.490e+02_pr)/(8.314_pr*Tloc))

    kH80f_0 = (1.11200000e+22_pr)*Tloc**(-5.000_pr)*&
                    exp(-(1.861e+04_pr)/(8.314_pr*Tloc))

    kH80f_inf = (1.08100000e+06_pr)*Tloc**(0.450_pr)*&
                    exp(-(7.623e+03_pr)/(8.314_pr*Tloc))

    FCH80f = (0.000e+00_pr)*&
                    exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)*&
                    exp(-Tloc/(9.500e+01_pr)) + (1.000e+00_pr)*&
                    exp(-(2.000e+02_pr)/Tloc)

    k(rH80f) =&
                    getlindratecoeff&
                    (Tloc, kH80f_0, kH80f_inf, FCH80f, M(mM12), Ploc )

    k(rH81) = (3.39000000e+00_pr)*Tloc**(1.880_pr)*&
                    exp(-(7.490e+02_pr)/(8.314_pr*Tloc))

    k(rH82) = (2.67900000e+02_pr)*Tloc**(0.890_pr)*&
                    exp(-(-8.042e+03_pr)/(8.314_pr*Tloc))

    k(rH83) = (1.95000000e+07_pr)*Tloc**(-0.500_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH84) = (5.83100000e+05_pr)*Tloc**(0.599_pr)*&
                    exp(-(-1.219e+04_pr)/(8.314_pr*Tloc))

    k(rH85f) = (1.22000000e+24_pr)*Tloc**(-5.760_pr)*&
                    exp(-(4.226e+04_pr)/(8.314_pr*Tloc))

    k(rH85b) = (1.25900000e+24_pr)*Tloc**(-5.630_pr)*&
                    exp(-(9.330e+04_pr)/(8.314_pr*Tloc))

    k(rH86f) = (5.54000000e-04_pr)*Tloc**(3.500_pr)*&
                    exp(-(2.162e+04_pr)/(8.314_pr*Tloc))

    k(rH86b) = (1.35500000e-07_pr)*Tloc**(4.060_pr)*&
                    exp(-(3.706e+04_pr)/(8.314_pr*Tloc))

    k(rH87) = (5.80000000e+01_pr)*Tloc**(1.730_pr)*&
                    exp(-(4.853e+03_pr)/(8.314_pr*Tloc))

    k(rH88) = (1.20000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH89) = (1.51000000e-13_pr)*Tloc**(6.000_pr)*&
                    exp(-(2.530e+04_pr)/(8.314_pr*Tloc))

    k(rH90) = (1.30000000e+01_pr)*Tloc**(2.130_pr)*&
                    exp(-(2.172e+04_pr)/(8.314_pr*Tloc))

    k(rH91) = (1.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH92) = (8.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH93f) = (1.10000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH93b) = (2.04600000e+06_pr)*Tloc**(0.890_pr)*&
                    exp(-(1.164e+05_pr)/(8.314_pr*Tloc))

    k(rH94) = (2.40000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(-3.573e+03_pr)/(8.314_pr*Tloc))

    k(rH95) = (2.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.757e+04_pr)/(8.314_pr*Tloc))

    k(rH96) = (4.88700000e+50_pr)*Tloc**(-12.250_pr)*&
                    exp(-(1.175e+05_pr)/(8.314_pr*Tloc))

    k(rH97) = (7.00000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(-4.184e+03_pr)/(8.314_pr*Tloc))

    k(rH98) = (2.39700000e+48_pr)*Tloc**(-9.900_pr)*&
                    exp(-(3.434e+05_pr)/(8.314_pr*Tloc))

    k(rH99) = (7.00000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(-4.184e+03_pr)/(8.314_pr*Tloc))

    k(rH100) = (9.72000000e+23_pr)*Tloc**(-5.710_pr)*&
                    exp(-(8.975e+04_pr)/(8.314_pr*Tloc))

    k(rH101) = (7.14000000e+09_pr)*Tloc**(-1.210_pr)*&
                    exp(-(8.807e+04_pr)/(8.314_pr*Tloc))

    k(rH102) = (3.33200000e+04_pr)*Tloc**(0.340_pr)*&
                    exp(-(-2.326e+03_pr)/(8.314_pr*Tloc))

    k(rH103f) = (6.30000000e+02_pr)*Tloc**(1.900_pr)*&
                    exp(-(7.611e+04_pr)/(8.314_pr*Tloc))

    k(rH103b) = (1.09700000e+02_pr)*Tloc**(1.890_pr)*&
                    exp(-(6.628e+04_pr)/(8.314_pr*Tloc))

    k(rH104) = (3.12000000e+00_pr)*Tloc**(2.000_pr)*&
                    exp(-(-1.247e+03_pr)/(8.314_pr*Tloc))

    k(rH105) = (5.24000000e+05_pr)*Tloc**(0.700_pr)*&
                    exp(-(2.462e+04_pr)/(8.314_pr*Tloc))

    k(rH106f) = (1.73000000e-01_pr)*Tloc**(2.500_pr)*&
                    exp(-(1.043e+04_pr)/(8.314_pr*Tloc))

    k(rH106b) = (7.93300000e-02_pr)*Tloc**(2.510_pr)*&
                    exp(-(8.167e+04_pr)/(8.314_pr*Tloc))

    k(rH107f) = (4.83000000e+27_pr)*Tloc**(-5.810_pr)*&
                    exp(-(7.740e+04_pr)/(8.314_pr*Tloc))

    k(rH107b) = (2.31300000e+27_pr)*Tloc**(-5.900_pr)*&
                    exp(-(1.323e+05_pr)/(8.314_pr*Tloc))

    k(rH108) = (1.58000000e+01_pr)*Tloc**(1.760_pr)*&
                    exp(-(-5.088e+03_pr)/(8.314_pr*Tloc))

    k(rH109f) = (2.73000000e+62_pr)*Tloc**(-13.280_pr)*&
                    exp(-(5.155e+05_pr)/(8.314_pr*Tloc))

    k(rH109b) = (4.71200000e+53_pr)*Tloc**(-13.190_pr)*&
                    exp(-(1.236e+05_pr)/(8.314_pr*Tloc))

    k(rH110) = (2.21000000e-06_pr)*Tloc**(3.500_pr)*&
                    exp(-(2.374e+04_pr)/(8.314_pr*Tloc))

    k(rH111) = (4.50000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(2.100e+04_pr)/(8.314_pr*Tloc))

    k(rH112f) = (8.56900000e+18_pr)*Tloc**(-1.570_pr)*&
                    exp(-(1.688e+05_pr)/(8.314_pr*Tloc))

    k(rH112b) = (1.30000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(6.527e+03_pr)/(8.314_pr*Tloc))

    k(rH113) = (3.00000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.255e+04_pr)/(8.314_pr*Tloc))

    k(rH114f) = (2.66700000e+15_pr)*Tloc**(-0.640_pr)*&
                    exp(-(1.541e+05_pr)/(8.314_pr*Tloc))

    k(rH114b) = (1.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.046e+04_pr)/(8.314_pr*Tloc))

    k(rH115f) = (2.28400000e+14_pr)*Tloc**(-0.550_pr)*&
                    exp(-(1.188e+05_pr)/(8.314_pr*Tloc))

    k(rH115b) = (4.10000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.014e+04_pr)/(8.314_pr*Tloc))

    k(rH116) = (2.08000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH117) = (2.81000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(2.176e+04_pr)/(8.314_pr*Tloc))

    k(rH118) = (4.67000000e+01_pr)*Tloc**(1.610_pr)*&
                    exp(-(-1.460e+02_pr)/(8.314_pr*Tloc))

    k(rH119) = (1.29000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(4.853e+04_pr)/(8.314_pr*Tloc))

    k(rH120) = (1.68000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(8.548e+04_pr)/(8.314_pr*Tloc))

    k(rH121) = (5.60000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(7.406e+04_pr)/(8.314_pr*Tloc))

    k(rH122) = (1.05400000e+04_pr)*Tloc**(0.970_pr)*&
                    exp(-(6.636e+03_pr)/(8.314_pr*Tloc))

    k(rH123) = (3.98000000e+05_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.975e+04_pr)/(8.314_pr*Tloc))

    k(rH124f) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
                    exp(-(1.987e+05_pr)/(8.314_pr*Tloc))

    k(rH124b) = (2.08000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH125) = (1.13000000e+08_pr)*Tloc**(0.000_pr)*&
                    exp(-(3.284e+04_pr)/(8.314_pr*Tloc))

    k(rH126f) = (3.97200000e+00_pr)*Tloc**(2.750_pr)*&
                    exp(-(2.827e+04_pr)/(8.314_pr*Tloc))

    k(rH127f) = (1.30000000e+00_pr)*Tloc**(2.400_pr)*&
                    exp(-(1.871e+04_pr)/(8.314_pr*Tloc))

    k(rH127b) = (4.70900000e-01_pr)*Tloc**(2.150_pr)*&
                    exp(-(5.096e+04_pr)/(8.314_pr*Tloc))

    k(rH128) = (2.02800000e+12_pr)*Tloc**(0.090_pr)*&
                    exp(-(9.858e+04_pr)/(8.314_pr*Tloc))

    k(rH129f) = (2.80300000e+17_pr)*Tloc**(-0.620_pr)*&
                    exp(-(1.508e+05_pr)/(8.314_pr*Tloc))

    k(rH129b) = (7.54000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    k(rH130f) = (3.36400000e+19_pr)*Tloc**(-1.320_pr)*&
                    exp(-(1.496e+05_pr)/(8.314_pr*Tloc))

    k(rH130b) = (4.52000000e+06_pr)*Tloc**(0.000_pr)*&
                    exp(-(0.000e+00_pr)/(8.314_pr*Tloc))

    kH5b_0 = (4.07313004e+30_pr)*Tloc**(-5.83_pr)*&
                    exp(-(2.2816e+05_pr)/(8.314_pr*Tloc))

    kH5b_inf = (1.65550435e+20_pr)*Tloc**(-1.57_pr)*&
                    exp(-(2.1959e+05_pr)/(8.314_pr*Tloc))

    FCH5b = (5.300e-01_pr)*&
                    exp(-Tloc/(1.000e+02_pr)) + (4.700e-01_pr)*&
                    exp(-Tloc/(2.000e+03_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+15_pr)/Tloc)

    k(rH5b) =&
                    getlindratecoeff&
                    (Tloc, kH5b_0, kH5b_inf, FCH5b, M(mM2), Ploc )

    kH8b_0 = (5.93645620e+10_pr)*Tloc**(-0.45_pr)*&
                    exp(-(1.9684e+05_pr)/(8.314_pr*Tloc))

    kH8b_inf = (2.50179226e+12_pr)*Tloc**(0.56_pr)*&
                    exp(-(2.0151e+05_pr)/(8.314_pr*Tloc))

    FCH8b = (5.000e-01_pr)*&
                    exp(-Tloc/(1.000e-30_pr)) + (5.000e-01_pr)*&
                    exp(-Tloc/(1.000e+30_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+100_pr)/Tloc)

    k(rH8b) =&
                    getlindratecoeff&
                    (Tloc, kH8b_0, kH8b_inf, FCH8b, M(mM4), Ploc )

    kH25b_0 = (7.03162509e+38_pr)*Tloc**(-6.81_pr)*&
                    exp(-(3.9196e+05_pr)/(8.314_pr*Tloc))

    kH25b_inf = (5.70831662e+25_pr)*Tloc**(-2.73_pr)*&
                    exp(-(3.8749e+05_pr)/(8.314_pr*Tloc))

    FCH25b = (5.950e-01_pr)*&
                    exp(-Tloc/(1.120e+03_pr)) + (4.050e-01_pr)*&
                    exp(-Tloc/(6.960e+01_pr)) + (1.000e+00_pr)*&
                    exp(-(1.000e+15_pr)/Tloc)

    k(rH25b) =&
                    getlindratecoeff&
                    (Tloc, kH25b_0, kH25b_inf, FCH25b, M(mM6), Ploc )

    k(rH30b) = (1.35039738e+11_pr)*Tloc**(-0.81_pr)*&
                    exp(-(6.7179e+03_pr)/(8.314_pr*Tloc))

    kH80b_0 = (3.30072675e+27_pr)*Tloc**(-4.90_pr)*&
                    exp(-(1.7164e+05_pr)/(8.314_pr*Tloc))

    kH80b_inf = (3.20871009e+11_pr)*Tloc**(0.55_pr)*&
                    exp(-(1.6065e+05_pr)/(8.314_pr*Tloc))

    FCH80b = (0.000e+00_pr)*&
                    exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)*&
                    exp(-Tloc/(9.500e+01_pr)) + (1.000e+00_pr)*&
                    exp(-(2.000e+02_pr)/Tloc)

    k(rH80b) =&
                    getlindratecoeff&
                    (Tloc, kH80b_0, kH80b_inf, FCH80b, M(mM12), Ploc )

    k(rH126b) = (9.85582206e-03_pr)*Tloc**(3.24_pr)*&
                    exp(-(4.3192e+04_pr)/(8.314_pr*Tloc))

  END SUBROUTINE get_rate_coefficients

! --- Reaction rates --- !
  SUBROUTINE get_reaction_rates(w,k,m,c,cqss)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nspec) :: c
    REAL(pr), DIMENSION(nqss) :: cqss
    REAL(pr), DIMENSION(nreac) :: w,k
    REAL(pr), DIMENSION(13) :: m

    w(rH1) = k(rH1)  * c(sO)**2_pr  * m(mM0)
    w(rH2) = k(rH2)  * c(sO) * c(sH) * m(mM1)
    w(rH3f) = k(rH3f)  * c(sO) * c(sH2)
    w(rH3b) = k(rH3b)  * c(sH) * c(sOH)
    w(rH4) = k(rH4)  * c(sH) * c(sOH) * m(mM3)
    w(rH5f) = k(rH5f)  * c(sOH)**2_pr
    w(rH6f) = k(rH6f)  * c(sOH) * c(sH2)
    w(rH6b) = k(rH6b)  * c(sH) * c(sH2O)
    w(rH7f) = k(rH7f)  * c(sO) * c(sH2O)
    w(rH7b) = k(rH7b)  * c(sOH)**2_pr
    w(rH8f) = k(rH8f)  * c(sH) * c(sO2)
    w(rH9f) = k(rH9f)  * c(sH) * c(sO2)
    w(rH9b) = k(rH9b)  * c(sO) * c(sOH)
    w(rH10f) = k(rH10f)  * c(sHO2) * c(sOH)
    w(rH10b) = k(rH10b)  * c(sH2O) * c(sO2)
    w(rH11f) = k(rH11f)  * c(sHO2) * c(sO)
    w(rH11b) = k(rH11b)  * c(sOH) * c(sO2)
    w(rH12) = k(rH12)  * c(sHO2) * c(sH)
    w(rH13) = k(rH13)  * c(sHO2)**2_pr
    w(rH14) = k(rH14)  * c(sHO2) * c(sH)
    w(rH15) = k(rH15)  * c(sHO2)**2_pr
    w(rH16f) = k(rH16f)  * c(sH2O2) * c(sH)
    w(rH16b) = k(rH16b)  * c(sH2) * c(sHO2)
    w(rH17) = k(rH17)  * c(sH2O2) * c(sO)
    w(rH18) = k(rH18)  * c(sH2O2) * c(sOH)
    w(rH19f) = k(rH19f)  * c(sH2O2) * c(sOH)
    w(rH19b) = k(rH19b)  * c(sH2O) * c(sHO2)
    w(rH20) = k(rH20)  * c(sH2O2) * c(sH)
    w(rH21) = k(rH21)  * cqss(sCH2GSGXCH2-nspec) * c(sO)
    w(rH22) = k(rH22)  * cqss(sCH2GSGXCH2-nspec) * c(sO2)
    w(rH23) = k(rH23)  * cqss(sCH2GSGXCH2-nspec) * c(sOH)
    w(rH24f) = k(rH24f)  * cqss(sCH2GSGXCH2-nspec) * c(sH2)
    w(rH24b) = k(rH24b)  * c(sCH3) * c(sH)
    w(rH25f) = k(rH25f)  * c(sCH3)**2_pr
    w(rH26f) = k(rH26f)  * c(sCH3) * c(sO2)
    w(rH26b) = k(rH26b)  * cqss(sCH3O-nspec) * c(sO)
    w(rH27) = k(rH27)  * c(sCH3) * c(sHO2)
    w(rH28f) = k(rH28f)  * c(sCH3) * c(sOH)
    w(rH28b) = k(rH28b)  * cqss(sCH2GSGXCH2-nspec) * c(sH2O)
    w(rH29) = k(rH29)  * c(sCH3) * c(sH)
    w(rH30f) = k(rH30f)  * c(sCH3)**2_pr
    w(rH31) = k(rH31)  * c(sCH3) * c(sO2)
    w(rH32) = k(rH32)  * c(sCH3) * c(sO)
    w(rH33) = k(rH33)  * c(sCH3) * c(sHO2)
    w(rH34) = k(rH34)  * c(sCH3) * c(sOH)
    w(rH35) = k(rH35)  * cqss(sCH2GSGXCH2-nspec) * c(sCH3)
    w(rH36) = k(rH36)  * c(sCH3) * c(sH2O2)
    w(rH37f) = k(rH37f)  * c(sCH4) * c(sOH)
    w(rH37b) = k(rH37b)  * c(sCH3) * c(sH2O)
    w(rH38f) = k(rH38f)  * c(sCH4) * c(sH)
    w(rH38b) = k(rH38b)  * c(sCH3) * c(sH2)
    w(rH39f) = k(rH39f)  * cqss(sCH2GSGXCH2-nspec) * c(sCH4)
    w(rH39b) = k(rH39b)  * c(sCH3)**2_pr
    w(rH40f) = k(rH40f)  * c(sCH4) * c(sO)
    w(rH40b) = k(rH40b)  * c(sCH3) * c(sOH)
    w(rH41) = k(rH41)  * c(sCO) * c(sHO2)
    w(rH42) = k(rH42)  * c(sCO) * c(sO)
    w(rH43f) = k(rH43f)  * c(sCO) * c(sOH)
    w(rH43b) = k(rH43b)  * c(sCO2) * c(sH)
    w(rH44f) = k(rH44f)  * c(sCO) * c(sO2)
    w(rH44b) = k(rH44b)  * c(sCO2) * c(sO)
    w(rH45) = k(rH45)  * cqss(sHCO-nspec) * c(sO)
    w(rH46) = k(rH46)  * cqss(sHCO-nspec) * c(sO)
    w(rH47) = k(rH47)  * cqss(sHCO-nspec) * c(sH)
    w(rH48f) = k(rH48f)  * cqss(sHCO-nspec) * m(mM9)
    w(rH48b) = k(rH48b)  * c(sH) * c(sCO) * m(mM9)
    w(rH49) = k(rH49)  * cqss(sHCO-nspec) * c(sOH)
    w(rH50) = k(rH50)  * cqss(sHCO-nspec) * c(sCH3)
    w(rH51) = k(rH51)  * cqss(sHCO-nspec) * c(sO2)
    w(rH52) = k(rH52)  * c(sCH2O) * c(sHO2)
    w(rH53) = k(rH53)  * c(sCH2O) * c(sOH)
    w(rH54) = k(rH54)  * c(sCH2O) * c(sH)
    w(rH55) = k(rH55)  * c(sCH2O) * c(sCH3)
    w(rH56) = k(rH56)  * c(sCH2O) * c(sO)
    w(rH57) = k(rH57)  * c(sCH2O) * c(sO2)
    w(rH58) = k(rH58)  * cqss(sCH3O-nspec) * c(sO2)
    w(rH59) = k(rH59)  * cqss(sCH3O-nspec)
    w(rH60) = k(rH60)  * cqss(sCH2GSGXCH2-nspec) * c(sCO2)
    w(rH61) = k(rH61)  * c(sCH3O2) * c(sHO2)
    w(rH62) = k(rH62)  * c(sCH3O2)**2_pr
    w(rH63f) = k(rH63f)  * c(sCH3O2) * m(mM5)
    w(rH63b) = k(rH63b)  * c(sCH3) * c(sO2) * m(mM5)
    w(rH64) = k(rH64)  * c(sCH3O2) * c(sCH2O)
    w(rH65) = k(rH65)  * c(sCH3O2) * c(sCH3)
    w(rH66f) = k(rH66f)  * c(sCH3O2H)
    w(rH66b) = k(rH66b)  * cqss(sCH3O-nspec) * c(sOH)
    w(rH67) = k(rH67)  * c(sC2H2) * c(sO)
    w(rH68) = k(rH68)  * c(sC2H2) * c(sO2)
    w(rH69) = k(rH69)  * cqss(sC2H3-nspec) * c(sO2)
    w(rH70) = k(rH70)  * cqss(sC2H3-nspec) * c(sO2)
    w(rH71) = k(rH71)  * cqss(sC2H3-nspec) * c(sO2)
    w(rH72) = k(rH72)  * cqss(sC2H3-nspec) * c(sH)
    w(rH73) = k(rH73)  * cqss(sC2H3-nspec)
    w(rH74) = k(rH74)  * c(sC2H4)
    w(rH75f) = k(rH75f)  * c(sC2H4) * c(sH)
    w(rH75b) = k(rH75b)  * cqss(sC2H3-nspec) * c(sH2)
    w(rH76) = k(rH76)  * c(sC2H4) * c(sCH3O2)
    w(rH77f) = k(rH77f)  * c(sC2H4) * c(sOH)
    w(rH77b) = k(rH77b)  * cqss(sC2H3-nspec) * c(sH2O)
    w(rH78f) = k(rH78f)  * c(sC2H4) * c(sCH3)
    w(rH78b) = k(rH78b)  * cqss(sC2H3-nspec) * c(sCH4)
    w(rH79) = k(rH79)  * c(sC2H4) * c(sO)
    w(rH80f) = k(rH80f)  * c(sH) * c(sC2H4)
    w(rH81) = k(rH81)  * c(sC2H4) * c(sO)
    w(rH82) = k(rH82)  * cqss(sC2H5-nspec) * c(sHO2)
    w(rH83) = k(rH83)  * c(sCH3) * cqss(sC2H5-nspec)
    w(rH84) = k(rH84)  * c(sH) * cqss(sC2H5-nspec)
    w(rH85f) = k(rH85f)  * cqss(sC2H5-nspec) * c(sO2)
    w(rH85b) = k(rH85b)  * c(sC2H4) * c(sHO2)
    w(rH86f) = k(rH86f)  * c(sC2H6) * c(sH)
    w(rH86b) = k(rH86b)  * cqss(sC2H5-nspec) * c(sH2)
    w(rH87) = k(rH87)  * c(sC2H6) * c(sOH)
    w(rH88) = k(rH88)  * cqss(sCH2GSGXCH2-nspec) * c(sC2H6)
    w(rH89) = k(rH89)  * c(sC2H6) * c(sCH3)
    w(rH90) = k(rH90)  * c(sC2H6) * c(sO)
    w(rH91) = k(rH91)  * cqss(sHCCO-nspec) * c(sOH)
    w(rH92) = k(rH92)  * cqss(sHCCO-nspec) * c(sO)
    w(rH93f) = k(rH93f)  * cqss(sHCCO-nspec) * c(sH)
    w(rH93b) = k(rH93b)  * cqss(sCH2GSGXCH2-nspec) * c(sCO)
    w(rH94) = k(rH94)  * cqss(sHCCO-nspec) * c(sO2)
    w(rH95) = k(rH95)  * cqss(sCH2CHO-nspec) * c(sO2)
    w(rH96) = k(rH96)  * cqss(sC3H5XAXC3H5-nspec) * c(sH)
    w(rH97) = k(rH97)  * cqss(sC3H5XAXC3H5-nspec) * c(sHO2)
    w(rH98) = k(rH98)  * cqss(sC3H5XAXC3H5-nspec)
    w(rH99) = k(rH99)  * cqss(sC3H5XAXC3H5-nspec) * c(sCH3O2)
    w(rH100) = k(rH100)  * cqss(sC3H5XAXC3H5-nspec) * c(sO2)
    w(rH101) = k(rH101)  * cqss(sC3H5XAXC3H5-nspec) * c(sO2)
    w(rH102) = k(rH102)  * cqss(sC3H5XAXC3H5-nspec) * c(sHO2)
    w(rH103f) = k(rH103f)  * cqss(sC3H5XAXC3H5-nspec) * c(sCH2O)
    w(rH103b) = k(rH103b)  * c(sC3H6) * cqss(sHCO-nspec)
    w(rH104) = k(rH104)  * c(sC3H6) * c(sOH)
    w(rH105) = k(rH105)  * c(sC3H6) * c(sO)
    w(rH106f) = k(rH106f)  * c(sC3H6) * c(sH)
    w(rH106b) = k(rH106b)  * cqss(sC3H5XAXC3H5-nspec) * c(sH2)
    w(rH107f) = k(rH107f)  * c(sC3H6) * c(sH)
    w(rH107b) = k(rH107b)  * c(sC2H4) * c(sCH3)
    w(rH108) = k(rH108)  * c(sC3H6) * c(sO)
    w(rH109f) = k(rH109f)  * c(sC3H6)
    w(rH109b) = k(rH109b)  * cqss(sC2H3-nspec) * c(sCH3)
    w(rH110) = k(rH110)  * c(sC3H6) * c(sCH3)
    w(rH111) = k(rH111)  * cqss(sIXC3H7-nspec) * c(sO2)
    w(rH112f) = k(rH112f)  * cqss(sIXC3H7-nspec)
    w(rH112b) = k(rH112b)  * c(sH) * c(sC3H6)
    w(rH113) = k(rH113)  * cqss(sNXC3H7-nspec) * c(sO2)
    w(rH114f) = k(rH114f)  * cqss(sNXC3H7-nspec)
    w(rH114b) = k(rH114b)  * c(sH) * c(sC3H6)
    w(rH115f) = k(rH115f)  * cqss(sNXC3H7-nspec)
    w(rH115b) = k(rH115b)  * c(sCH3) * c(sC2H4)
    w(rH116) = k(rH116)  * cqss(sNXC3H7-nspec) * c(sHO2)
    w(rH117) = k(rH117)  * c(sC3H8) * c(sO)
    w(rH118) = k(rH118)  * c(sC3H8) * c(sOH)
    w(rH119) = k(rH119)  * c(sCH3) * c(sC3H8)
    w(rH120) = k(rH120)  * c(sC3H8) * c(sHO2)
    w(rH121) = k(rH121)  * c(sC3H8) * c(sHO2)
    w(rH122) = k(rH122)  * c(sC3H8) * c(sOH)
    w(rH123) = k(rH123)  * c(sCH3) * c(sC3H8)
    w(rH124f) = k(rH124f)  * c(sC3H8) * c(sO2)
    w(rH124b) = k(rH124b)  * cqss(sIXC3H7-nspec) * c(sHO2)
    w(rH125) = k(rH125)  * c(sC3H8) * c(sO)
    w(rH126f) = k(rH126f)  * c(sH) * c(sC3H8)
    w(rH127f) = k(rH127f)  * c(sH) * c(sC3H8)
    w(rH127b) = k(rH127b)  * c(sH2) * cqss(sIXC3H7-nspec)
    w(rH128) = k(rH128)  * c(sC3H5O)
    w(rH129f) = k(rH129f)  * cqss(sIXC3H7O2-nspec)
    w(rH129b) = k(rH129b)  * cqss(sIXC3H7-nspec) * c(sO2)
    w(rH130f) = k(rH130f)  * cqss(sNXC3H7O2-nspec)
    w(rH130b) = k(rH130b)  * cqss(sNXC3H7-nspec) * c(sO2)
    w(rH5b) = k(rH5b)  * c(sH2O2)
    w(rH8b) = k(rH8b)  * c(sHO2)
    w(rH25b) = k(rH25b)  * c(sC2H6)
    w(rH30b) = k(rH30b)  * c(sH) * cqss(sC2H5-nspec)
    w(rH80b) = k(rH80b)  * cqss(sC2H5-nspec)
    w(rH126b) = k(rH126b)  * c(sH2) * cqss(sNXC3H7-nspec)

  END SUBROUTINE get_reaction_rates

! --- Production rates --- !
  SUBROUTINE get_production_rates(cdot,w)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nspec) :: cdot
    REAL(pr), DIMENSION(nreac) :: w

    cdot(sN2) = 0.0_pr


    cdot(sO) = 0.0_pr - 2_pr *&
                   w(rH1) -&
                   w(rH2) -&
                   w(rH3f) +&
                   w(rH3b) -&
                   w(rH7f) +&
                   w(rH7b) +&
                   w(rH9f) -&
                   w(rH9b) -&
                   w(rH11f) +&
                   w(rH11b) -&
                   w(rH17) -&
                   w(rH21) +&
                   w(rH26f) -&
                   w(rH26b) -&
                   w(rH32) -&
                   w(rH40f) +&
                   w(rH40b) -&
                   w(rH42) +&
                   w(rH44f) -&
                   w(rH44b) -&
                   w(rH45) -&
                   w(rH46) -&
                   w(rH56) -&
                   w(rH67) +&
                   w(rH70) -&
                   w(rH79) -&
                   w(rH81) -&
                   w(rH90) -&
                   w(rH92) -&
                   w(rH105) -&
                   w(rH108) -&
                   w(rH117) -&
                   w(rH125)


    cdot(sO2) = 0.0_pr +&
                   w(rH1) -&
                   w(rH8f) -&
                   w(rH9f) +&
                   w(rH9b) +&
                   w(rH10f) -&
                   w(rH10b) +&
                   w(rH11f) -&
                   w(rH11b) +&
                   w(rH12) +&
                   w(rH13) +&
                   w(rH15) -&
                   w(rH22) -&
                   w(rH26f) +&
                   w(rH26b) -&
                   w(rH31) +&
                   w(rH33) -&
                   w(rH44f) +&
                   w(rH44b) -&
                   w(rH51) -&
                   w(rH57) -&
                   w(rH58) +&
                   w(rH61) +&
                   w(rH62) +&
                   w(rH63f) -&
                   w(rH63b) -&
                   w(rH68) -&
                   w(rH69) -&
                   w(rH70) -&
                   w(rH71) +&
                   w(rH82) -&
                   w(rH85f) +&
                   w(rH85b) -&
                   w(rH94) -&
                   w(rH95) -&
                   w(rH100) -&
                   w(rH101) +&
                   w(rH102) -&
                   w(rH111) -&
                   w(rH113) +&
                   w(rH116) -&
                   w(rH124f) +&
                   w(rH124b) +&
                   w(rH129f) -&
                   w(rH129b) +&
                   w(rH130f) -&
                   w(rH130b) +&
                   w(rH8b)


    cdot(sH) = 0.0_pr -&
                   w(rH2) +&
                   w(rH3f) -&
                   w(rH3b) -&
                   w(rH4) +&
                   w(rH6f) -&
                   w(rH6b) -&
                   w(rH8f) -&
                   w(rH9f) +&
                   w(rH9b) -&
                   w(rH12) -&
                   w(rH14) -&
                   w(rH16f) +&
                   w(rH16b) -&
                   w(rH20) + 2_pr *&
                   w(rH21) +&
                   w(rH22) +&
                   w(rH23) +&
                   w(rH24f) -&
                   w(rH24b) -&
                   w(rH29) +&
                   w(rH30f) +&
                   w(rH32) +&
                   w(rH35) -&
                   w(rH38f) +&
                   w(rH38b) +&
                   w(rH43f) -&
                   w(rH43b) +&
                   w(rH46) -&
                   w(rH47) +&
                   w(rH48f) -&
                   w(rH48b) -&
                   w(rH54) +&
                   w(rH59) +&
                   w(rH67) -&
                   w(rH72) +&
                   w(rH73) -&
                   w(rH75f) +&
                   w(rH75b) -&
                   w(rH80f) +&
                   w(rH81) -&
                   w(rH84) -&
                   w(rH86f) +&
                   w(rH86b) +&
                   w(rH92) -&
                   w(rH93f) +&
                   w(rH93b) -&
                   w(rH96) -&
                   w(rH106f) +&
                   w(rH106b) -&
                   w(rH107f)
    cdot(sH) = cdot(sH) +&
                   w(rH107b) +&
                   w(rH112f) -&
                   w(rH112b) +&
                   w(rH114f) -&
                   w(rH114b) -&
                   w(rH126f) -&
                   w(rH127f) +&
                   w(rH127b) +&
                   w(rH8b) -&
                   w(rH30b) +&
                   w(rH80b) +&
                   w(rH126b)


    cdot(sOH) = 0.0_pr +&
                   w(rH2) +&
                   w(rH3f) -&
                   w(rH3b) -&
                   w(rH4) - 2_pr *&
                   w(rH5f) -&
                   w(rH6f) +&
                   w(rH6b) + 2_pr *&
                   w(rH7f) - 2_pr *&
                   w(rH7b) +&
                   w(rH9f) -&
                   w(rH9b) -&
                   w(rH10f) +&
                   w(rH10b) +&
                   w(rH11f) -&
                   w(rH11b) + 2_pr *&
                   w(rH14) +&
                   w(rH17) -&
                   w(rH18) -&
                   w(rH19f) +&
                   w(rH19b) +&
                   w(rH20) +&
                   w(rH22) -&
                   w(rH23) +&
                   w(rH27) -&
                   w(rH28f) +&
                   w(rH28b) +&
                   w(rH31) -&
                   w(rH34) -&
                   w(rH37f) +&
                   w(rH37b) +&
                   w(rH40f) -&
                   w(rH40b) +&
                   w(rH41) -&
                   w(rH43f) +&
                   w(rH43b) +&
                   w(rH45) -&
                   w(rH49) -&
                   w(rH53) +&
                   w(rH56) +&
                   w(rH66f) -&
                   w(rH66b) +&
                   w(rH68) -&
                   w(rH77f) +&
                   w(rH77b) -&
                   w(rH87) +&
                   w(rH90) -&
                   w(rH91) +&
                   w(rH95) +&
                   w(rH97) +&
                   w(rH100)
    cdot(sOH) = cdot(sOH) -&
                   w(rH104) +&
                   w(rH105) +&
                   w(rH117) -&
                   w(rH118) -&
                   w(rH122) +&
                   w(rH125) + 2_pr *&
                   w(rH5b)


    cdot(sH2) = 0.0_pr -&
                   w(rH3f) +&
                   w(rH3b) -&
                   w(rH6f) +&
                   w(rH6b) +&
                   w(rH12) +&
                   w(rH16f) -&
                   w(rH16b) -&
                   w(rH24f) +&
                   w(rH24b) +&
                   w(rH34) +&
                   w(rH38f) -&
                   w(rH38b) +&
                   w(rH47) +&
                   w(rH54) +&
                   w(rH72) +&
                   w(rH74) +&
                   w(rH75f) -&
                   w(rH75b) +&
                   w(rH86f) -&
                   w(rH86b) +&
                   w(rH106f) -&
                   w(rH106b) +&
                   w(rH126f) +&
                   w(rH127f) -&
                   w(rH127b) -&
                   w(rH126b)


    cdot(sH2O) = 0.0_pr +&
                   w(rH4) +&
                   w(rH6f) -&
                   w(rH6b) -&
                   w(rH7f) +&
                   w(rH7b) +&
                   w(rH10f) -&
                   w(rH10b) +&
                   w(rH18) +&
                   w(rH19f) -&
                   w(rH19b) +&
                   w(rH20) +&
                   w(rH28f) -&
                   w(rH28b) +&
                   w(rH37f) -&
                   w(rH37b) +&
                   w(rH49) +&
                   w(rH53) +&
                   w(rH77f) -&
                   w(rH77b) +&
                   w(rH87) +&
                   w(rH104) +&
                   w(rH118) +&
                   w(rH122)


    cdot(sH2O2) = 0.0_pr +&
                   w(rH5f) +&
                   w(rH13) +&
                   w(rH15) -&
                   w(rH16f) +&
                   w(rH16b) -&
                   w(rH17) -&
                   w(rH18) -&
                   w(rH19f) +&
                   w(rH19b) -&
                   w(rH20) -&
                   w(rH36) +&
                   w(rH52) +&
                   w(rH120) +&
                   w(rH121) -&
                   w(rH5b)


    cdot(sHO2) = 0.0_pr +&
                   w(rH8f) -&
                   w(rH10f) +&
                   w(rH10b) -&
                   w(rH11f) +&
                   w(rH11b) -&
                   w(rH12) - 2_pr *&
                   w(rH13) -&
                   w(rH14) - 2_pr *&
                   w(rH15) +&
                   w(rH16f) -&
                   w(rH16b) +&
                   w(rH17) +&
                   w(rH18) +&
                   w(rH19f) -&
                   w(rH19b) -&
                   w(rH27) -&
                   w(rH33) +&
                   w(rH36) -&
                   w(rH41) +&
                   w(rH51) -&
                   w(rH52) +&
                   w(rH57) +&
                   w(rH58) -&
                   w(rH61) +&
                   w(rH69) -&
                   w(rH82) +&
                   w(rH85f) -&
                   w(rH85b) -&
                   w(rH97) -&
                   w(rH102) +&
                   w(rH111) +&
                   w(rH113) -&
                   w(rH116) -&
                   w(rH120) -&
                   w(rH121) +&
                   w(rH124f) -&
                   w(rH124b) -&
                   w(rH8b)

!cdot(sCH2GSGXCH2) = 0.0_pr


    cdot(sCO) = 0.0_pr +&
                   w(rH21) +&
                   w(rH22) -&
                   w(rH41) -&
                   w(rH42) -&
                   w(rH43f) +&
                   w(rH43b) -&
                   w(rH44f) +&
                   w(rH44b) +&
                   w(rH45) +&
                   w(rH47) +&
                   w(rH48f) -&
                   w(rH48b) +&
                   w(rH49) +&
                   w(rH50) +&
                   w(rH51) +&
                   w(rH60) + 2_pr *&
                   w(rH92) +&
                   w(rH93f) -&
                   w(rH93b) +&
                   w(rH95)


    cdot(sCH2O) = 0.0_pr +&
                   w(rH23) +&
                   w(rH31) +&
                   w(rH32) +&
                   w(rH34) -&
                   w(rH52) -&
                   w(rH53) -&
                   w(rH54) -&
                   w(rH55) -&
                   w(rH56) -&
                   w(rH57) +&
                   w(rH58) +&
                   w(rH59) +&
                   w(rH60) -&
                   w(rH64) +&
                   w(rH71) +&
                   w(rH95) +&
                   w(rH100) +&
                   w(rH101) -&
                   w(rH103f) +&
                   w(rH103b) +&
                   w(rH128)


    cdot(sCH3) = 0.0_pr +&
                   w(rH24f) -&
                   w(rH24b) - 2_pr *&
                   w(rH25f) -&
                   w(rH26f) +&
                   w(rH26b) -&
                   w(rH27) -&
                   w(rH28f) +&
                   w(rH28b) -&
                   w(rH29) - 2_pr *&
                   w(rH30f) -&
                   w(rH31) -&
                   w(rH32) -&
                   w(rH33) -&
                   w(rH34) -&
                   w(rH35) -&
                   w(rH36) +&
                   w(rH37f) -&
                   w(rH37b) +&
                   w(rH38f) -&
                   w(rH38b) + 2_pr *&
                   w(rH39f) - 2_pr *&
                   w(rH39b) +&
                   w(rH40f) -&
                   w(rH40b) -&
                   w(rH50) -&
                   w(rH55) +&
                   w(rH63f) -&
                   w(rH63b) -&
                   w(rH65) -&
                   w(rH78f) +&
                   w(rH78b) +&
                   w(rH79) -&
                   w(rH83) +&
                   w(rH88) -&
                   w(rH89) +&
                   w(rH98) +&
                   w(rH107f) -&
                   w(rH107b) +&
                   w(rH109f) -&
                   w(rH109b) -&
                   w(rH110) +&
                   w(rH115f) -&
                   w(rH115b) -&
                   w(rH119) -&
                   w(rH123) + 2_pr *&
                   w(rH25b) + 2_pr *&
                   w(rH30b)


    cdot(sC2H6) = 0.0_pr +&
                   w(rH25f) +&
                   w(rH82) +&
                   w(rH84) -&
                   w(rH86f) +&
                   w(rH86b) -&
                   w(rH87) -&
                   w(rH88) -&
                   w(rH89) -&
                   w(rH90) -&
                   w(rH25b)

!cdot(sCH3O) = 0.0_pr


    cdot(sCH4) = 0.0_pr +&
                   w(rH29) +&
                   w(rH33) +&
                   w(rH36) -&
                   w(rH37f) +&
                   w(rH37b) -&
                   w(rH38f) +&
                   w(rH38b) -&
                   w(rH39f) +&
                   w(rH39b) -&
                   w(rH40f) +&
                   w(rH40b) +&
                   w(rH50) +&
                   w(rH55) +&
                   w(rH78f) -&
                   w(rH78b) +&
                   w(rH83) +&
                   w(rH89) +&
                   w(rH110) +&
                   w(rH119) +&
                   w(rH123)

!cdot(sC2H5) = 0.0_pr


    cdot(sC2H4) = 0.0_pr +&
                   w(rH35) -&
                   w(rH74) -&
                   w(rH75f) +&
                   w(rH75b) -&
                   w(rH76) -&
                   w(rH77f) +&
                   w(rH77b) -&
                   w(rH78f) +&
                   w(rH78b) -&
                   w(rH79) -&
                   w(rH80f) -&
                   w(rH81) +&
                   w(rH83) +&
                   w(rH85f) -&
                   w(rH85b) +&
                   w(rH107f) -&
                   w(rH107b) +&
                   w(rH115f) -&
                   w(rH115b) +&
                   w(rH80b)


    cdot(sCO2) = 0.0_pr +&
                   w(rH41) +&
                   w(rH42) +&
                   w(rH43f) -&
                   w(rH43b) +&
                   w(rH44f) -&
                   w(rH44b) +&
                   w(rH46) -&
                   w(rH60) +&
                   w(rH94)

!cdot(sHCO) = 0.0_pr


    cdot(sCH3O2) = 0.0_pr -&
                   w(rH61) - 2_pr *&
                   w(rH62) -&
                   w(rH63f) +&
                   w(rH63b) -&
                   w(rH64) -&
                   w(rH65) -&
                   w(rH76) -&
                   w(rH99)


    cdot(sCH3O2H) = 0.0_pr +&
                   w(rH61) +&
                   w(rH64) -&
                   w(rH66f) +&
                   w(rH66b) +&
                   w(rH76)


    cdot(sC2H2) = 0.0_pr -&
                   w(rH67) -&
                   w(rH68) +&
                   w(rH69) +&
                   w(rH72) +&
                   w(rH73) +&
                   w(rH74) +&
                   w(rH98) +&
                   w(rH100)

!cdot(sHCCO) = 0.0_pr

!cdot(sC2H3) = 0.0_pr

!cdot(sCH2CHO) = 0.0_pr

!cdot(sC3H5XAXC3H5) = 0.0_pr


    cdot(sC3H6) = 0.0_pr +&
                   w(rH96) +&
                   w(rH102) +&
                   w(rH103f) -&
                   w(rH103b) -&
                   w(rH104) -&
                   w(rH105) -&
                   w(rH106f) +&
                   w(rH106b) -&
                   w(rH107f) +&
                   w(rH107b) -&
                   w(rH108) -&
                   w(rH109f) +&
                   w(rH109b) -&
                   w(rH110) +&
                   w(rH111) +&
                   w(rH112f) -&
                   w(rH112b) +&
                   w(rH113) +&
                   w(rH114f) -&
                   w(rH114b)


    cdot(sC3H5O) = 0.0_pr +&
                   w(rH97) +&
                   w(rH99) -&
                   w(rH128)

!cdot(sIXC3H7) = 0.0_pr

!cdot(sNXC3H7) = 0.0_pr


    cdot(sC3H8) = 0.0_pr +&
                   w(rH116) -&
                   w(rH117) -&
                   w(rH118) -&
                   w(rH119) -&
                   w(rH120) -&
                   w(rH121) -&
                   w(rH122) -&
                   w(rH123) -&
                   w(rH124f) +&
                   w(rH124b) -&
                   w(rH125) -&
                   w(rH126f) -&
                   w(rH127f) +&
                   w(rH127b) +&
                   w(rH126b)

!cdot(sIXC3H7O2) = 0.0_pr

!cdot(sNXC3H7O2) = 0.0_pr


  END SUBROUTINE get_production_rates

! --- Actual reactions --- !
  SUBROUTINE reaction_expressions
    IMPLICIT NONE

    reacexp(1) = '2 O + M0 -> O2 + M0'
    reacexp(2) = 'O + H + M1 -> OH + M1'
    reacexp(3) = 'O + H2 -> H + OH'
    reacexp(4) = 'H + OH -> O + H2'
    reacexp(5) = 'H + OH + M3 -> H2O + M3'
    reacexp(6) = '2 OH + M2 -> H2O2 + M2'
    reacexp(7) = 'OH + H2 -> H + H2O'
    reacexp(8) = 'H + H2O -> OH + H2'
    reacexp(9) = 'O + H2O -> 2 OH'
    reacexp(10) = '2 OH -> O + H2O'
    reacexp(11) = 'H + O2 + M4 -> HO2 + M4'
    reacexp(12) = 'H + O2 -> O + OH'
    reacexp(13) = 'O + OH -> H + O2'
    reacexp(14) = 'HO2 + OH -> H2O + O2'
    reacexp(15) = 'H2O + O2 -> HO2 + OH'
    reacexp(16) = 'HO2 + O -> OH + O2'
    reacexp(17) = 'OH + O2 -> HO2 + O'
    reacexp(18) = 'HO2 + H -> H2 + O2'
    reacexp(19) = '2 HO2 -> H2O2 + O2'
    reacexp(20) = 'HO2 + H -> 2 OH'
    reacexp(21) = '2 HO2 -> H2O2 + O2'
    reacexp(22) = 'H2O2 + H -> H2 + HO2'
    reacexp(23) = 'H2 + HO2 -> H2O2 + H'
    reacexp(24) = 'H2O2 + O -> OH + HO2'
    reacexp(25) = 'H2O2 + OH -> H2O + HO2'
    reacexp(26) = 'H2O2 + OH -> H2O + HO2'
    reacexp(27) = 'H2O + HO2 -> H2O2 + OH'
    reacexp(28) = 'H2O2 + H -> H2O + OH'
    reacexp(29) = 'CH2GSG-CH2 + O -> CO + 2 H'
    reacexp(30) = 'CH2GSG-CH2 + O2 -> CO + OH + H'
    reacexp(31) = 'CH2GSG-CH2 + OH -> CH2O + H'
    reacexp(32) = 'CH2GSG-CH2 + H2 -> CH3 + H'
    reacexp(33) = 'CH3 + H -> CH2GSG-CH2 + H2'
    reacexp(34) = '2 CH3 + M6 -> C2H6 + M6'
    reacexp(35) = 'CH3 + O2 -> CH3O + O'
    reacexp(36) = 'CH3O + O -> CH3 + O2'
    reacexp(37) = 'CH3 + HO2 -> CH3O + OH'
    reacexp(38) = 'CH3 + OH -> CH2GSG-CH2 + H2O'
    reacexp(39) = 'CH2GSG-CH2 + H2O -> CH3 + OH'
    reacexp(40) = 'CH3 + H + M7 -> CH4 + M7'
    reacexp(41) = '2 CH3 -> H + C2H5'
    reacexp(42) = 'CH3 + O2 -> CH2O + OH'
    reacexp(43) = 'CH3 + O -> CH2O + H'
    reacexp(44) = 'CH3 + HO2 -> CH4 + O2'
    reacexp(45) = 'CH3 + OH -> CH2O + H2'
    reacexp(46) = 'CH2GSG-CH2 + CH3 -> C2H4 + H'
    reacexp(47) = 'CH3 + H2O2 -> CH4 + HO2'
    reacexp(48) = 'CH4 + OH -> CH3 + H2O'
    reacexp(49) = 'CH3 + H2O -> CH4 + OH'
    reacexp(50) = 'CH4 + H -> CH3 + H2'
    reacexp(51) = 'CH3 + H2 -> CH4 + H'
    reacexp(52) = 'CH2GSG-CH2 + CH4 -> 2 CH3'
    reacexp(53) = '2 CH3 -> CH2GSG-CH2 + CH4'
    reacexp(54) = 'CH4 + O -> CH3 + OH'
    reacexp(55) = 'CH3 + OH -> CH4 + O'
    reacexp(56) = 'CO + HO2 -> CO2 + OH'
    reacexp(57) = 'CO + O + M8 -> CO2 + M8'
    reacexp(58) = 'CO + OH -> CO2 + H'
    reacexp(59) = 'CO2 + H -> CO + OH'
    reacexp(60) = 'CO + O2 -> CO2 + O'
    reacexp(61) = 'CO2 + O -> CO + O2'
    reacexp(62) = 'HCO + O -> CO + OH'
    reacexp(63) = 'HCO + O -> CO2 + H'
    reacexp(64) = 'HCO + H -> CO + H2'
    reacexp(65) = 'HCO + M9 -> H + CO + M9'
    reacexp(66) = 'H + CO + M9 -> HCO + M9'
    reacexp(67) = 'HCO + OH -> CO + H2O'
    reacexp(68) = 'HCO + CH3 -> CH4 + CO'
    reacexp(69) = 'HCO + O2 -> CO + HO2'
    reacexp(70) = 'CH2O + HO2 -> HCO + H2O2'
    reacexp(71) = 'CH2O + OH -> HCO + H2O'
    reacexp(72) = 'CH2O + H -> HCO + H2'
    reacexp(73) = 'CH2O + CH3 -> HCO + CH4'
    reacexp(74) = 'CH2O + O -> HCO + OH'
    reacexp(75) = 'CH2O + O2 -> HCO + HO2'
    reacexp(76) = 'CH3O + O2 -> CH2O + HO2'
    reacexp(77) = 'CH3O + M5 -> CH2O + H + M5'
    reacexp(78) = 'CH2GSG-CH2 + CO2 -> CH2O + CO'
    reacexp(79) = 'CH3O2 + HO2 -> CH3O2H + O2'
    reacexp(80) = '2 CH3O2 -> O2 + 2 CH3O'
    reacexp(81) = 'CH3O2 + M5 -> CH3 + O2 + M5'
    reacexp(82) = 'CH3 + O2 + M5 -> CH3O2 + M5'
    reacexp(83) = 'CH3O2 + CH2O -> CH3O2H + HCO'
    reacexp(84) = 'CH3O2 + CH3 -> 2 CH3O'
    reacexp(85) = 'CH3O2H -> CH3O + OH'
    reacexp(86) = 'CH3O + OH -> CH3O2H'
    reacexp(87) = 'C2H2 + O -> HCCO + H'
    reacexp(88) = 'C2H2 + O2 -> HCCO + OH'
    reacexp(89) = 'C2H3 + O2 -> C2H2 + HO2'
    reacexp(90) = 'C2H3 + O2 -> CH2CHO + O'
    reacexp(91) = 'C2H3 + O2 -> CH2O + HCO'
    reacexp(92) = 'C2H3 + H -> C2H2 + H2'
    reacexp(93) = 'C2H3 + M11 -> H + C2H2 + M11'
    reacexp(94) = 'C2H4 + M5 -> C2H2 + H2 + M5'
    reacexp(95) = 'C2H4 + H -> C2H3 + H2'
    reacexp(96) = 'C2H3 + H2 -> C2H4 + H'
    reacexp(97) = 'C2H4 + CH3O2 -> C2H3 + CH3O2H'
    reacexp(98) = 'C2H4 + OH -> C2H3 + H2O'
    reacexp(99) = 'C2H3 + H2O -> C2H4 + OH'
    reacexp(100) = 'C2H4 + CH3 -> C2H3 + CH4'
    reacexp(101) = 'C2H3 + CH4 -> C2H4 + CH3'
    reacexp(102) = 'C2H4 + O -> CH3 + HCO'
    reacexp(103) = 'H + C2H4 + M12 -> C2H5 + M12'
    reacexp(104) = 'C2H4 + O -> CH2CHO + H'
    reacexp(105) = 'C2H5 + HO2 -> C2H6 + O2'
    reacexp(106) = 'CH3 + C2H5 -> CH4 + C2H4'
    reacexp(107) = 'H + C2H5 -> C2H6'
    reacexp(108) = 'C2H5 + O2 -> C2H4 + HO2'
    reacexp(109) = 'C2H4 + HO2 -> C2H5 + O2'
    reacexp(110) = 'C2H6 + H -> C2H5 + H2'
    reacexp(111) = 'C2H5 + H2 -> C2H6 + H'
    reacexp(112) = 'C2H6 + OH -> C2H5 + H2O'
    reacexp(113) = 'CH2GSG-CH2 + C2H6 -> CH3 + C2H5'
    reacexp(114) = 'C2H6 + CH3 -> C2H5 + CH4'
    reacexp(115) = 'C2H6 + O -> C2H5 + OH'
    reacexp(116) = 'HCCO + OH -> 2 HCO'
    reacexp(117) = 'HCCO + O -> H + 2 CO'
    reacexp(118) = 'HCCO + H -> CH2GSG-CH2 + CO'
    reacexp(119) = 'CH2GSG-CH2 + CO -> HCCO + H'
    reacexp(120) = 'HCCO + O2 -> CO2 + HCO'
    reacexp(121) = 'CH2CHO + O2 -> CH2O + CO + OH'
    reacexp(122) = 'C3H5-A-C3H5 + H -> C3H6'
    reacexp(123) = 'C3H5-A-C3H5 + HO2 -> C3H5O + OH'
    reacexp(124) = 'C3H5-A-C3H5 -> C2H2 + CH3'
    reacexp(125) = 'C3H5-A-C3H5 + CH3O2 -> C3H5O + CH3O'
    reacexp(126) = 'C3H5-A-C3H5 + O2 -> C2H2 + CH2O + OH'
    reacexp(127) = 'C3H5-A-C3H5 + O2 -> CH2CHO + CH2O'
    reacexp(128) = 'C3H5-A-C3H5 + HO2 -> C3H6 + O2'
    reacexp(129) = 'C3H5-A-C3H5 + CH2O -> C3H6 + HCO'
    reacexp(130) = 'C3H6 + HCO -> C3H5-A-C3H5 + CH2O'
    reacexp(131) = 'C3H6 + OH -> C3H5-A-C3H5 + H2O'
    reacexp(132) = 'C3H6 + O -> C3H5-A-C3H5 + OH'
    reacexp(133) = 'C3H6 + H -> C3H5-A-C3H5 + H2'
    reacexp(134) = 'C3H5-A-C3H5 + H2 -> C3H6 + H'
    reacexp(135) = 'C3H6 + H -> C2H4 + CH3'
    reacexp(136) = 'C2H4 + CH3 -> C3H6 + H'
    reacexp(137) = 'C3H6 + O -> C2H5 + HCO'
    reacexp(138) = 'C3H6 -> C2H3 + CH3'
    reacexp(139) = 'C2H3 + CH3 -> C3H6'
    reacexp(140) = 'C3H6 + CH3 -> C3H5-A-C3H5 + CH4'
    reacexp(141) = 'I-C3H7 + O2 -> C3H6 + HO2'
    reacexp(142) = 'I-C3H7 -> H + C3H6'
    reacexp(143) = 'H + C3H6 -> I-C3H7'
    reacexp(144) = 'N-C3H7 + O2 -> C3H6 + HO2'
    reacexp(145) = 'N-C3H7 -> H + C3H6'
    reacexp(146) = 'H + C3H6 -> N-C3H7'
    reacexp(147) = 'N-C3H7 -> CH3 + C2H4'
    reacexp(148) = 'CH3 + C2H4 -> N-C3H7'
    reacexp(149) = 'N-C3H7 + HO2 -> C3H8 + O2'
    reacexp(150) = 'C3H8 + O -> I-C3H7 + OH'
    reacexp(151) = 'C3H8 + OH -> I-C3H7 + H2O'
    reacexp(152) = 'CH3 + C3H8 -> CH4 + N-C3H7'
    reacexp(153) = 'C3H8 + HO2 -> N-C3H7 + H2O2'
    reacexp(154) = 'C3H8 + HO2 -> I-C3H7 + H2O2'
    reacexp(155) = 'C3H8 + OH -> N-C3H7 + H2O'
    reacexp(156) = 'CH3 + C3H8 -> CH4 + I-C3H7'
    reacexp(157) = 'C3H8 + O2 -> I-C3H7 + HO2'
    reacexp(158) = 'I-C3H7 + HO2 -> C3H8 + O2'
    reacexp(159) = 'C3H8 + O -> N-C3H7 + OH'
    reacexp(160) = 'H + C3H8 -> H2 + N-C3H7'
    reacexp(161) = 'H + C3H8 -> H2 + I-C3H7'
    reacexp(162) = 'H2 + I-C3H7 -> H + C3H8'
    reacexp(163) = 'C3H5O -> C2H3 + CH2O'
    reacexp(164) = 'I-C3H7O2 -> I-C3H7 + O2'
    reacexp(165) = 'I-C3H7 + O2 -> I-C3H7O2'
    reacexp(166) = 'N-C3H7O2 -> N-C3H7 + O2'
    reacexp(167) = 'N-C3H7 + O2 -> N-C3H7O2'
    reacexp(168) = 'Reverse of 2 OH + M2 -> H2O2 + M2'
    reacexp(169) = 'Reverse of H + O2 + M4 -> HO2 + M4'
    reacexp(170) = 'Reverse of 2 CH3 + M6 -> C2H6 + M6'
    reacexp(171) = 'Reverse of 2 CH3 -> H + C2H5'
    reacexp(172) = 'Reverse of H + C2H4 + M12 -> C2H5 + M12'
    reacexp(173) = 'Reverse of H + C3H8 -> H2 + N-C3H7'

  END SUBROUTINE reaction_expressions

! --- Forward/Backward link --- !
  SUBROUTINE reverse_reactions
    IMPLICIT NONE

    fofb = 0

! Attach corresponding forward reaction to each backward reaction
    fofb(4) = 3
    fofb(8) = 7
    fofb(10) = 9
    fofb(13) = 12
    fofb(15) = 14
    fofb(17) = 16
    fofb(23) = 22
    fofb(27) = 26
    fofb(33) = 32
    fofb(36) = 35
    fofb(39) = 38
    fofb(49) = 48
    fofb(51) = 50
    fofb(53) = 52
    fofb(55) = 54
    fofb(59) = 58
    fofb(61) = 60
    fofb(66) = 65
    fofb(82) = 81
    fofb(86) = 85
    fofb(96) = 95
    fofb(99) = 98
    fofb(101) = 100
    fofb(109) = 108
    fofb(111) = 110
    fofb(119) = 118
    fofb(130) = 129
    fofb(134) = 133
    fofb(136) = 135
    fofb(139) = 138
    fofb(143) = 142
    fofb(146) = 145
    fofb(148) = 147
    fofb(158) = 157
    fofb(162) = 161
    fofb(165) = 164
    fofb(167) = 166
    fofb(168) = 6
    fofb(169) = 11
    fofb(170) = 34
    fofb(171) = 41
    fofb(172) = 103
    fofb(173) = 160

  END SUBROUTINE reverse_reactions

! --- Evaluation of QSS concentrations --- !
  SUBROUTINE get_QSS(cqss,c,k,M)
    IMPLICIT NONE

    REAL(pr), DIMENSION(nqss) :: cqss
    REAL(pr), DIMENSION(nspec) :: c
    REAL(pr), DIMENSION(nreac) :: k
    REAL(pr), DIMENSION(13) :: M

    REAL(pr) :: IXC3H7O2_denom1&
                   , IXC3H7O2_ct1&
                   , IXC3H7O2_denom2&
                   , IXC3H7O2_ct2
    REAL(pr) :: IXC3H7_denom1&
                   , IXC3H7_ct1&
                   , IXC3H7_denom2&
                   , IXC3H7_ct2&
                   , IXC3H7_IXC3H7O2&
                   , IXC3H7_IXC3H7O2_coeff
    REAL(pr) :: NXC3H7O2_denom1&
                   , NXC3H7O2_ct1&
                   , NXC3H7O2_denom2&
                   , NXC3H7O2_ct2&
                   , NXC3H7O2_NXC3H7&
                   , NXC3H7O2_NXC3H7_coeff
    REAL(pr) :: NXC3H7_denom1&
                   , NXC3H7_ct1&
                   , NXC3H7_denom2&
                   , NXC3H7_ct2
    REAL(pr) :: HCCO_denom1&
                   , HCCO_ct1&
                   , HCCO_denom2&
                   , HCCO_ct2&
                   , HCCO_CH2GSGXCH2&
                   , HCCO_CH2GSGXCH2_coeff
    REAL(pr) :: C3H5XAXC3H5_denom1&
                   , C3H5XAXC3H5_ct1&
                   , C3H5XAXC3H5_denom2&
                   , C3H5XAXC3H5_ct2
    REAL(pr) :: HCO_denom1&
                   , HCO_ct1&
                   , HCO_denom2&
                   , HCO_ct2&
                   , HCO_C3H5XAXC3H5&
                   , HCO_C3H5XAXC3H5_coeff
    REAL(pr) :: CH2GSGXCH2_denom1&
                   , CH2GSGXCH2_ct1&
                   , CH2GSGXCH2_denom2&
                   , CH2GSGXCH2_ct2


! c(sIXC3H7) c(sIXC3H7O2) (coupled)  --------------------

! Primary denominators-----------------------

    IXC3H7_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH111) *  c(sO2)&
                         + k(rH112f)&
                         + k(rH124b) *  c(sHO2)&
                         + k(rH127b) *  c(sH2)&
                         + k(rH129b) *  c(sO2)&
                         )

    IXC3H7O2_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH129f)&
                         )


! Primary constant parts -----------------------

    IXC3H7_ct1 = ( 0.0_pr&
                          + k(rH112b) * c(sH) * c(sC3H6)&
                          + k(rH117) * c(sC3H8) * c(sO)&
                          + k(rH118) * c(sC3H8) * c(sOH)&
                          + k(rH121) * c(sC3H8) * c(sHO2)&
                          + k(rH123) * c(sCH3) * c(sC3H8)&
                          + k(rH124f) * c(sC3H8) * c(sO2)&
                          + k(rH127f) * c(sH) * c(sC3H8)                )

    IXC3H7O2_ct1 = ( 0.0_pr                )


! IXC3H7 ---------------------------------------

    IXC3H7_denom2 = tiny(1.0_pr) + ( IXC3H7_denom1                 )

    IXC3H7_ct2 = ( IXC3H7_ct1&
                          ) / IXC3H7_denom2

    IXC3H7_IXC3H7O2 = ( 0.0_pr&
                          + k(rH129f)&
                          ) / IXC3H7_denom2



    IXC3H7_IXC3H7O2_coeff = ( 0.0_pr&
                          + k(rH129b)  * c(sO2)&
                          )


! IXC3H7O2 ---------------------------------------

    IXC3H7O2_denom2 = tiny(1.0_pr) + ( IXC3H7O2_denom1&
                          - IXC3H7_IXC3H7O2_coeff * IXC3H7_IXC3H7O2                )

    IXC3H7O2_ct2 = ( IXC3H7O2_ct1&
                          + IXC3H7_IXC3H7O2_coeff * IXC3H7_ct2&
                          ) / IXC3H7O2_denom2




! Reconstruction ------------------------------------

    cqss(sIXC3H7O2-nspec) = ( IXC3H7O2_ct2                 )

    cqss(sIXC3H7-nspec) = ( IXC3H7_ct2&
                          + IXC3H7_IXC3H7O2 * cqss(sIXC3H7O2-nspec)                )

! c(sNXC3H7O2) c(sNXC3H7) (coupled)  --------------------

! Primary denominators-----------------------

    NXC3H7O2_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH130f)&
                         )

    NXC3H7_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH113) *  c(sO2)&
                         + k(rH114f)&
                         + k(rH115f)&
                         + k(rH116) *  c(sHO2)&
                         + k(rH126b) *  c(sH2)&
                         + k(rH130b) *  c(sO2)&
                         )


! Primary constant parts -----------------------

    NXC3H7O2_ct1 = ( 0.0_pr                )

    NXC3H7_ct1 = ( 0.0_pr&
                          + k(rH114b) * c(sH) * c(sC3H6)&
                          + k(rH115b) * c(sCH3) * c(sC2H4)&
                          + k(rH119) * c(sCH3) * c(sC3H8)&
                          + k(rH120) * c(sC3H8) * c(sHO2)&
                          + k(rH122) * c(sC3H8) * c(sOH)&
                          + k(rH125) * c(sC3H8) * c(sO)&
                          + k(rH126f) * c(sH) * c(sC3H8)                )


! NXC3H7O2 ---------------------------------------

    NXC3H7O2_denom2 = tiny(1.0_pr) + ( NXC3H7O2_denom1                 )

    NXC3H7O2_ct2 = ( NXC3H7O2_ct1&
                          ) / NXC3H7O2_denom2

    NXC3H7O2_NXC3H7 = ( 0.0_pr&
                          + k(rH130b)  * c(sO2)&
                          ) / NXC3H7O2_denom2



    NXC3H7O2_NXC3H7_coeff = ( 0.0_pr&
                          + k(rH130f)&
                          )


! NXC3H7 ---------------------------------------

    NXC3H7_denom2 = tiny(1.0_pr) + ( NXC3H7_denom1&
                          - NXC3H7O2_NXC3H7_coeff * NXC3H7O2_NXC3H7                )

    NXC3H7_ct2 = ( NXC3H7_ct1&
                          + NXC3H7O2_NXC3H7_coeff * NXC3H7O2_ct2&
                          ) / NXC3H7_denom2




! Reconstruction ------------------------------------

    cqss(sNXC3H7-nspec) = ( NXC3H7_ct2                 )

    cqss(sNXC3H7O2-nspec) = ( NXC3H7O2_ct2&
                          + NXC3H7O2_NXC3H7 * cqss(sNXC3H7-nspec)                )

! c(sHCCO) c(sCH2GSGXCH2) (coupled)  --------------------

! Primary denominators-----------------------

    HCCO_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH91) *  c(sOH)&
                         + k(rH92) *  c(sO)&
                         + k(rH93f) *  c(sH)&
                         + k(rH94) *  c(sO2)&
                         )

    CH2GSGXCH2_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH21) *  c(sO)&
                         + k(rH22) *  c(sO2)&
                         + k(rH23) *  c(sOH)&
                         + k(rH24f) *  c(sH2)&
                         + k(rH28b) *  c(sH2O)&
                         + k(rH35) *  c(sCH3)&
                         + k(rH39f) *  c(sCH4)&
                         + k(rH60) *  c(sCO2)&
                         + k(rH88) *  c(sC2H6)&
                         + k(rH93b) *  c(sCO)&
                         )


! Primary constant parts -----------------------

    HCCO_ct1 = ( 0.0_pr&
                          + k(rH67) * c(sC2H2) * c(sO)&
                          + k(rH68) * c(sC2H2) * c(sO2)                )

    CH2GSGXCH2_ct1 = ( 0.0_pr&
                          + k(rH24b) * c(sCH3) * c(sH)&
                          + k(rH28f) * c(sCH3) * c(sOH)&
                          + k(rH39b) * c(sCH3) * c(sCH3)                )


! HCCO ---------------------------------------

    HCCO_denom2 = tiny(1.0_pr) + ( HCCO_denom1                 )

    HCCO_ct2 = ( HCCO_ct1&
                          ) / HCCO_denom2

    HCCO_CH2GSGXCH2 = ( 0.0_pr&
                          + k(rH93b)  * c(sCO)&
                          ) / HCCO_denom2



    HCCO_CH2GSGXCH2_coeff = ( 0.0_pr&
                          + k(rH93f)  * c(sH)&
                          )


! CH2GSGXCH2 ---------------------------------------

    CH2GSGXCH2_denom2 = tiny(1.0_pr) + ( CH2GSGXCH2_denom1&
                          - HCCO_CH2GSGXCH2_coeff * HCCO_CH2GSGXCH2                )

    CH2GSGXCH2_ct2 = ( CH2GSGXCH2_ct1&
                          + HCCO_CH2GSGXCH2_coeff * HCCO_ct2&
                          ) / CH2GSGXCH2_denom2




! Reconstruction ------------------------------------

    cqss(sCH2GSGXCH2-nspec) = ( CH2GSGXCH2_ct2                 )

    cqss(sHCCO-nspec) = ( HCCO_ct2&
                          + HCCO_CH2GSGXCH2 * cqss(sCH2GSGXCH2-nspec)                )

! cqss(sC2H5) (uncoupled) --------------------

    cqss(sC2H5-nspec) = ( 0.0_pr&
                          + k(rH30f) * c(sCH3) * c(sCH3)&
                          + k(rH80f) * c(sH) * c(sC2H4)&
                          + k(rH85b) * c(sC2H4) * c(sHO2)&
                          + k(rH86f) * c(sC2H6) * c(sH)&
                          + k(rH87) * c(sC2H6) * c(sOH)&
                          + k(rH88) * cqss(sCH2GSGXCH2-nspec) *c(sC2H6)&
                          + k(rH89) * c(sC2H6) * c(sCH3)&
                          + k(rH90) * c(sC2H6) * c(sO)&
                          + k(rH108) * c(sC3H6) * c(sO)&
                          ) / ( tiny(1.0_pr) + (&
                          + k(rH30b) *  c(sH)&
                          + k(rH80b)&
                          + k(rH82) *  c(sHO2)&
                          + k(rH83) *  c(sCH3)&
                          + k(rH84) *  c(sH)&
                          + k(rH85f) *  c(sO2)&
                          + k(rH86b) *  c(sH2)                 ) )


! cqss(sC2H3) (uncoupled) --------------------

    cqss(sC2H3-nspec) = ( 0.0_pr&
                          + k(rH75f) * c(sC2H4) * c(sH)&
                          + k(rH76) * c(sC2H4) * c(sCH3O2)&
                          + k(rH77f) * c(sC2H4) * c(sOH)&
                          + k(rH78f) * c(sC2H4) * c(sCH3)&
                          + k(rH109f) * c(sC3H6)&
                          + k(rH128) * c(sC3H5O)&
                          ) / ( tiny(1.0_pr) + (&
                          + k(rH69) *  c(sO2)&
                          + k(rH70) *  c(sO2)&
                          + k(rH71) *  c(sO2)&
                          + k(rH72) *  c(sH)&
                          + k(rH73)&
                          + k(rH75b) *  c(sH2)&
                          + k(rH77b) *  c(sH2O)&
                          + k(rH78b) *  c(sCH4)&
                          + k(rH109b) *  c(sCH3)                 ) )


! c(sHCO) c(sC3H5XAXC3H5) (coupled)  --------------------

! Primary denominators-----------------------

    HCO_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH45) *  c(sO)&
                         + k(rH46) *  c(sO)&
                         + k(rH47) *  c(sH)&
                         + k(rH48f) * M(mM9)&
                         + k(rH49) *  c(sOH)&
                         + k(rH50) *  c(sCH3)&
                         + k(rH51) *  c(sO2)&
                         + k(rH103b) *  c(sC3H6)&
                         )

    C3H5XAXC3H5_denom1 = tiny(1.0_pr) + ( 0.0_pr&
                         + k(rH96) *  c(sH)&
                         + k(rH97) *  c(sHO2)&
                         + k(rH98)&
                         + k(rH99) *  c(sCH3O2)&
                         + k(rH100) *  c(sO2)&
                         + k(rH101) *  c(sO2)&
                         + k(rH102) *  c(sHO2)&
                         + k(rH103f) *  c(sCH2O)&
                         + k(rH106b) *  c(sH2)&
                         )


! Primary constant parts -----------------------

    HCO_ct1 = ( 0.0_pr&
                          + k(rH48b) * c(sH) * c(sCO) * M(mM9)&
                          + k(rH52) * c(sCH2O) * c(sHO2)&
                          + k(rH53) * c(sCH2O) * c(sOH)&
                          + k(rH54) * c(sCH2O) * c(sH)&
                          + k(rH55) * c(sCH2O) * c(sCH3)&
                          + k(rH56) * c(sCH2O) * c(sO)&
                          + k(rH57) * c(sCH2O) * c(sO2)&
                          + k(rH64) * c(sCH3O2) * c(sCH2O)&
                          + k(rH71) * cqss(sC2H3-nspec) *c(sO2)&
                          + k(rH79) * c(sC2H4) * c(sO)&
                          + k(rH91) * 2.0 * cqss(sHCCO-nspec) *c(sOH)&
                          + k(rH94) * cqss(sHCCO-nspec) *c(sO2)&
                          + k(rH108) * c(sC3H6) * c(sO)                )

    C3H5XAXC3H5_ct1 = ( 0.0_pr&
                          + k(rH104) * c(sC3H6) * c(sOH)&
                          + k(rH105) * c(sC3H6) * c(sO)&
                          + k(rH106f) * c(sC3H6) * c(sH)&
                          + k(rH110) * c(sC3H6) * c(sCH3)                )


! HCO ---------------------------------------

    HCO_denom2 = tiny(1.0_pr) + ( HCO_denom1                 )

    HCO_ct2 = ( HCO_ct1&
                          ) / HCO_denom2

    HCO_C3H5XAXC3H5 = ( 0.0_pr&
                          + k(rH103f)  * c(sCH2O)&
                          ) / HCO_denom2



    HCO_C3H5XAXC3H5_coeff = ( 0.0_pr&
                          + k(rH103b) * c(sC3H6)&
                          )


! C3H5XAXC3H5 ---------------------------------------

    C3H5XAXC3H5_denom2 = tiny(1.0_pr) + ( C3H5XAXC3H5_denom1&
                          - HCO_C3H5XAXC3H5_coeff * HCO_C3H5XAXC3H5                )

    C3H5XAXC3H5_ct2 = ( C3H5XAXC3H5_ct1&
                          + HCO_C3H5XAXC3H5_coeff * HCO_ct2&
                          ) / C3H5XAXC3H5_denom2




! Reconstruction ------------------------------------

    cqss(sC3H5XAXC3H5-nspec) = ( C3H5XAXC3H5_ct2                 )

    cqss(sHCO-nspec) = ( HCO_ct2&
                          + HCO_C3H5XAXC3H5 * cqss(sC3H5XAXC3H5-nspec)                )

! cqss(sCH2CHO) (uncoupled) --------------------

    cqss(sCH2CHO-nspec) = ( 0.0_pr&
                          + k(rH70) * cqss(sC2H3-nspec) *c(sO2)&
                          + k(rH81) * c(sC2H4) * c(sO)&
                          + k(rH101) * cqss(sC3H5XAXC3H5-nspec) *c(sO2)&
                          ) / ( tiny(1.0_pr) + (&
                          + k(rH95) *  c(sO2)                 ) )


! cqss(sCH3O) (uncoupled) --------------------

    cqss(sCH3O-nspec) = ( 0.0_pr&
                          + k(rH26f) * c(sCH3) * c(sO2)&
                          + k(rH27) * c(sCH3) * c(sHO2)&
                          + k(rH62) * 2.0 * c(sCH3O2) * c(sCH3O2)&
                          + k(rH65) * 2.0 * c(sCH3O2) * c(sCH3)&
                          + k(rH66f) * c(sCH3O2H)&
                          + k(rH99) * cqss(sC3H5XAXC3H5-nspec) *c(sCH3O2)&
                          ) / ( tiny(1.0_pr) + (&
                          + k(rH26b) *  c(sO)&
                          + k(rH58) *  c(sO2)&
                          + k(rH59)&
                          + k(rH66b) *  c(sOH)                 ) )


  END SUBROUTINE get_QSS
END MODULE mC3H82217312FC

!subroutine externe au module pour CANTERA
SUBROUTINE customkinetics( P, T, Y, WDOT )
    use mC3H82217312FC

    IMPLICIT NONE

    REAL(pr), DIMENSION(nspec) :: c, Y,  WDOT
    real(pr), dimension(nspec + nqss) :: H,dH,Cp
    REAL(pr), DIMENSION(nqss) :: cqss
    REAL(pr), DIMENSION(nreac) :: w,k
    REAL(pr), DIMENSION(13) :: m

    REAL(pr) :: T,P

    CALL YtoC(c,P,T,Y)
    CALL get_thirdbodies ( M,c )
    CALL get_rate_coefficients ( k,M,T,P )
    CALL get_QSS ( cqss,c,k,M )
    CALL get_reaction_rates ( w,k,M,c,cqss )
    CALL get_production_rates ( WDOT,w )

    RETURN
END SUBROUTINE customkinetics