!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE JetAPAH_29_277_15_LG.f90
!>    MODULE for calculating the analytical source terms for
!!    @file    File that contains the module and subroutine used to calculate the source terms for ethylene (C2H4) analytical mechanism
!!    @authors L. Gallen
!!    @date    27/01/2020
!!    @since   2.3
!!    @note    Analytical reduction using DRGEP and manual reduction techniques (Bisetti et al . CnF 2012)
!!             from Blanquart mechanism (Blanquart et al. 2009),
!!             QSS step using ARCANE (Q. Cazeres, P.Pepiot-Desjardins, J. Wirtz)
!!             Combination of previous ARC (C2H4_28_205_14_LG) and HyChem pyrolysis model(Wang et al. CnF 2019, Felden et al. CnF 2018)
!!             to predict PAH burning aviation jet fuel (Jet A1 : POSF10325)
!!             Validity range : 300 K
!!                              1 bar
!!                              Phi : 0.6-2.4
!--------------------------------------------------------------------------------------------------

!--------------------------------------------------------------------------------------------------
!     MODULE mod_customkinetics
!>    @details Generated by ARCANE custom kinetics routine to compute the chemical source terms.
!!    @authors Q. Cazères, P. Pepiot
!!    @date    2019/01/24
!--------------------------------------------------------------------------------------------------
module mod_customkinetics
  implicit none

  integer, parameter :: pr = selected_real_kind(15,307)

  ! Ideal gas constant
  real(pr), parameter :: Rcst = 8.3144621_pr
  
  ! Use semi-implicit formulation for select species
  logical, parameter :: semiImplicit = .False.
  
  ! Number of elements in the chemical system 
  integer, parameter :: ne = 3
  
  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 29
  integer, parameter :: nqss = 15
  integer, parameter :: nreac = 212
  integer, parameter :: nreac_reverse = 65
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 3
  integer, parameter :: nFO = 13
  integer, parameter :: nTB_reverse = 1
  integer, parameter :: nFO_reverse = 7
  
  ! Index of elements
  integer, parameter :: eC = 1
  integer, parameter :: eH = 2
  integer, parameter :: eO = 3

  ! Index of species
  integer, parameter :: sPOSF10325 = 1
  integer, parameter :: sN2 = 2
  integer, parameter :: sO = 3
  integer, parameter :: sH2 = 4
  integer, parameter :: sH = 5
  integer, parameter :: sOH = 6
  integer, parameter :: sH2O = 7
  integer, parameter :: sO2 = 8
  integer, parameter :: sHO2 = 9
  integer, parameter :: sCO = 10
  integer, parameter :: sCH2O = 11
  integer, parameter :: sCH3 = 12
  integer, parameter :: sCO2 = 13
  integer, parameter :: sCH4 = 14
  integer, parameter :: sC2H4 = 15
  integer, parameter :: sC2H5 = 16
  integer, parameter :: sC2H2 = 17
  integer, parameter :: sC3H3 = 18
  integer, parameter :: sAXC3H5 = 19
  integer, parameter :: sC2H6 = 20
  integer, parameter :: sPXC3H4 = 21
  integer, parameter :: sAXC3H4 = 22
  integer, parameter :: sA1 = 23
  integer, parameter :: sC5H5 = 24
  integer, parameter :: sC3H6 = 25
  integer, parameter :: sC5H6 = 26
  integer, parameter :: sA2 = 27
  integer, parameter :: sA1CH3 = 28
  integer, parameter :: sA1C2H = 29

  integer, parameter :: sqssSXCH2 = 1
  integer, parameter :: sqssCH = 2
  integer, parameter :: sqssHCO = 3
  integer, parameter :: sqssC2H3 = 4
  integer, parameter :: sqssHCCO = 5
  integer, parameter :: sqssC2H = 6
  integer, parameter :: sqssNXC3H7 = 7
  integer, parameter :: sqssA1C2H2 = 8
  integer, parameter :: sqssTXCH2 = 9
  integer, parameter :: sqssA1C2 = 10
  integer, parameter :: sqssA1CHO = 11
  integer, parameter :: sqssA1CH2 = 12
  integer, parameter :: sqssA1X = 13
  integer, parameter :: sqssA2X = 14
  integer, parameter :: sqssC4H8 = 15

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r39f = 39
  integer, parameter :: r40f = 40
  integer, parameter :: r41f = 41
  integer, parameter :: r42f = 42
  integer, parameter :: r43f = 43
  integer, parameter :: r44f = 44
  integer, parameter :: r45f = 45
  integer, parameter :: r46f = 46
  integer, parameter :: r47f = 47
  integer, parameter :: r48f = 48
  integer, parameter :: r49f = 49
  integer, parameter :: r50f = 50
  integer, parameter :: r51f = 51
  integer, parameter :: r52f = 52
  integer, parameter :: r53f = 53
  integer, parameter :: r54f = 54
  integer, parameter :: r55f = 55
  integer, parameter :: r56f = 56
  integer, parameter :: r57f = 57
  integer, parameter :: r58f = 58
  integer, parameter :: r59f = 59
  integer, parameter :: r60f = 60
  integer, parameter :: r61f = 61
  integer, parameter :: r62f = 62
  integer, parameter :: r63f = 63
  integer, parameter :: r64f = 64
  integer, parameter :: r65f = 65
  integer, parameter :: r66f = 66
  integer, parameter :: r67f = 67
  integer, parameter :: r68f = 68
  integer, parameter :: r69f = 69
  integer, parameter :: r70f = 70
  integer, parameter :: r71f = 71
  integer, parameter :: r72f = 72
  integer, parameter :: r73f = 73
  integer, parameter :: r74f = 74
  integer, parameter :: r75f = 75
  integer, parameter :: r76f = 76
  integer, parameter :: r77f = 77
  integer, parameter :: r78f = 78
  integer, parameter :: r79f = 79
  integer, parameter :: r80f = 80
  integer, parameter :: r81f = 81
  integer, parameter :: r82f = 82
  integer, parameter :: r83f = 83
  integer, parameter :: r84f = 84
  integer, parameter :: r85f = 85
  integer, parameter :: r86f = 86
  integer, parameter :: r87f = 87
  integer, parameter :: r88f = 88
  integer, parameter :: r89f = 89
  integer, parameter :: r90f = 90
  integer, parameter :: r91f = 91
  integer, parameter :: r92f = 92
  integer, parameter :: r93f = 93
  integer, parameter :: r94f = 94
  integer, parameter :: r95f = 95
  integer, parameter :: r96f = 96
  integer, parameter :: r97f = 97
  integer, parameter :: r98f = 98
  integer, parameter :: r99f = 99
  integer, parameter :: r100f = 100
  integer, parameter :: r101f = 101
  integer, parameter :: r102f = 102
  integer, parameter :: r103f = 103
  integer, parameter :: r104f = 104
  integer, parameter :: r105f = 105
  integer, parameter :: r106f = 106
  integer, parameter :: r107f = 107
  integer, parameter :: r108f = 108
  integer, parameter :: r109f = 109
  integer, parameter :: r110f = 110
  integer, parameter :: r111f = 111
  integer, parameter :: r112f = 112
  integer, parameter :: r113f = 113
  integer, parameter :: r114f = 114
  integer, parameter :: r115f = 115
  integer, parameter :: r116f = 116
  integer, parameter :: r117f = 117
  integer, parameter :: r118f = 118
  integer, parameter :: r119f = 119
  integer, parameter :: r120f = 120
  integer, parameter :: r121f = 121
  integer, parameter :: r122f = 122
  integer, parameter :: r123f = 123
  integer, parameter :: r124f = 124
  integer, parameter :: r125f = 125
  integer, parameter :: r126f = 126
  integer, parameter :: r127f = 127
  integer, parameter :: r128f = 128
  integer, parameter :: r129f = 129
  integer, parameter :: r130f = 130
  integer, parameter :: r131f = 131
  integer, parameter :: r132f = 132
  integer, parameter :: r133f = 133
  integer, parameter :: r134f = 134
  integer, parameter :: r135f = 135
  integer, parameter :: r136f = 136
  integer, parameter :: r137f = 137
  integer, parameter :: r138f = 138
  integer, parameter :: r139f = 139
  integer, parameter :: r140f = 140
  integer, parameter :: r141f = 141
  integer, parameter :: r142f = 142
  integer, parameter :: r143f = 143
  integer, parameter :: r144f = 144
  integer, parameter :: r145f = 145
  integer, parameter :: r146f = 146
  integer, parameter :: r147f = 147
  integer, parameter :: r148f = 148
  integer, parameter :: r149f = 149
  integer, parameter :: r150f = 150
  integer, parameter :: r151f = 151
  integer, parameter :: r152f = 152
  integer, parameter :: r153f = 153
  integer, parameter :: r154f = 154
  integer, parameter :: r155f = 155
  integer, parameter :: r156f = 156
  integer, parameter :: r157f = 157
  integer, parameter :: r158f = 158
  integer, parameter :: r159f = 159
  integer, parameter :: r160f = 160
  integer, parameter :: r161f = 161
  integer, parameter :: r162f = 162
  integer, parameter :: r163f = 163
  integer, parameter :: r164f = 164
  integer, parameter :: r165f = 165
  integer, parameter :: r166f = 166
  integer, parameter :: r167f = 167
  integer, parameter :: r168f = 168
  integer, parameter :: r169f = 169
  integer, parameter :: r170f = 170
  integer, parameter :: r171f = 171
  integer, parameter :: r172f = 172
  integer, parameter :: r173f = 173
  integer, parameter :: r174f = 174
  integer, parameter :: r175f = 175
  integer, parameter :: r176f = 176
  integer, parameter :: r177f = 177
  integer, parameter :: r178f = 178
  integer, parameter :: r179f = 179
  integer, parameter :: r180f = 180
  integer, parameter :: r181f = 181
  integer, parameter :: r182f = 182
  integer, parameter :: r183f = 183
  integer, parameter :: r184f = 184
  integer, parameter :: r185f = 185
  integer, parameter :: r186f = 186
  integer, parameter :: r187f = 187
  integer, parameter :: r188f = 188
  integer, parameter :: r189f = 189
  integer, parameter :: r190f = 190
  integer, parameter :: r191f = 191
  integer, parameter :: r192f = 192
  integer, parameter :: r193f = 193
  integer, parameter :: r194f = 194
  integer, parameter :: r195f = 195
  integer, parameter :: r196f = 196
  integer, parameter :: r197f = 197
  integer, parameter :: r198f = 198
  integer, parameter :: r199f = 199
  integer, parameter :: r200f = 200
  integer, parameter :: r201f = 201
  integer, parameter :: r202f = 202
  integer, parameter :: r203f = 203
  integer, parameter :: r204f = 204
  integer, parameter :: r205f = 205
  integer, parameter :: r206f = 206
  integer, parameter :: r207f = 207
  integer, parameter :: r208f = 208
  integer, parameter :: r209f = 209
  integer, parameter :: r210f = 210
  integer, parameter :: r211f = 211
  integer, parameter :: r212f = 212
  integer, parameter :: r6b = 213
  integer, parameter :: r7b = 214
  integer, parameter :: r10b = 215
  integer, parameter :: r11b = 216
  integer, parameter :: r13b = 217
  integer, parameter :: r25b = 218
  integer, parameter :: r26b = 219
  integer, parameter :: r30b = 220
  integer, parameter :: r34b = 221
  integer, parameter :: r36b = 222
  integer, parameter :: r40b = 223
  integer, parameter :: r41b = 224
  integer, parameter :: r48b = 225
  integer, parameter :: r50b = 226
  integer, parameter :: r51b = 227
  integer, parameter :: r52b = 228
  integer, parameter :: r54b = 229
  integer, parameter :: r56b = 230
  integer, parameter :: r57b = 231
  integer, parameter :: r59b = 232
  integer, parameter :: r64b = 233
  integer, parameter :: r65b = 234
  integer, parameter :: r68b = 235
  integer, parameter :: r73b = 236
  integer, parameter :: r78b = 237
  integer, parameter :: r79b = 238
  integer, parameter :: r86b = 239
  integer, parameter :: r96b = 240
  integer, parameter :: r97b = 241
  integer, parameter :: r105b = 242
  integer, parameter :: r106b = 243
  integer, parameter :: r112b = 244
  integer, parameter :: r125b = 245
  integer, parameter :: r126b = 246
  integer, parameter :: r131b = 247
  integer, parameter :: r132b = 248
  integer, parameter :: r133b = 249
  integer, parameter :: r134b = 250
  integer, parameter :: r135b = 251
  integer, parameter :: r141b = 252
  integer, parameter :: r142b = 253
  integer, parameter :: r143b = 254
  integer, parameter :: r147b = 255
  integer, parameter :: r153b = 256
  integer, parameter :: r154b = 257
  integer, parameter :: r155b = 258
  integer, parameter :: r156b = 259
  integer, parameter :: r171b = 260
  integer, parameter :: r182b = 261
  integer, parameter :: r183b = 262
  integer, parameter :: r189b = 263
  integer, parameter :: r190b = 264
  integer, parameter :: r194b = 265
  integer, parameter :: r195b = 266
  integer, parameter :: r196b = 267
  integer, parameter :: r203b = 268
  integer, parameter :: r204b = 269
  integer, parameter :: r205b = 270
  integer, parameter :: r206b = 271
  integer, parameter :: r207b = 272
  integer, parameter :: r208b = 273
  integer, parameter :: r209b = 274
  integer, parameter :: r210b = 275
  integer, parameter :: r211b = 276
  integer, parameter :: r212b = 277
  
  ! Index of third body species
  integer, parameter :: mM8 = 1
  integer, parameter :: mM9 = 2
  integer, parameter :: mM64 = 3
  
  integer, parameter :: mM12 = 4
  integer, parameter :: mM33 = 5
  integer, parameter :: mM36 = 6
  integer, parameter :: mM39 = 7
  integer, parameter :: mM55 = 8
  integer, parameter :: mM79 = 9
  integer, parameter :: mM87 = 10
  integer, parameter :: mM97 = 11
  integer, parameter :: mM105 = 12
  integer, parameter :: mM106 = 13
  integer, parameter :: mM112 = 14
  integer, parameter :: mM117 = 15
  integer, parameter :: mM153 = 16
  
  ! Index of third body reactions
  integer, parameter :: TBr8f = 1
  integer, parameter :: TBr9f = 2
  integer, parameter :: TBr64f = 3
  integer, parameter :: TBr64b = 4
  
  ! Index of fall off reactions
  integer, parameter :: FOr12f = 1
  integer, parameter :: FOr33f = 2
  integer, parameter :: FOr36f = 3
  integer, parameter :: FOr39f = 4
  integer, parameter :: FOr55f = 5
  integer, parameter :: FOr79f = 6
  integer, parameter :: FOr87f = 7
  integer, parameter :: FOr97f = 8
  integer, parameter :: FOr105f = 9
  integer, parameter :: FOr106f = 10
  integer, parameter :: FOr112f = 11
  integer, parameter :: FOr117f = 12
  integer, parameter :: FOr153f = 13
  integer, parameter :: FOr36b = 14
  integer, parameter :: FOr79b = 15
  integer, parameter :: FOr97b = 16
  integer, parameter :: FOr105b = 17
  integer, parameter :: FOr106b = 18
  integer, parameter :: FOr112b = 19
  integer, parameter :: FOr153b = 20
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       1.542957e-01_pr, & ! POSF10325
       2.801348e-02_pr, & ! N2
       1.599940e-02_pr, & ! O
       2.015880e-03_pr, & ! H2
       1.007940e-03_pr, & ! H
       1.700734e-02_pr, & ! OH
       1.801528e-02_pr, & ! H2O
       3.199880e-02_pr, & ! O2
       3.300674e-02_pr, & ! HO2
       2.801040e-02_pr, & ! CO
       3.002628e-02_pr, & ! CH2O
       1.503482e-02_pr, & ! CH3
       4.400980e-02_pr, & ! CO2
       1.604276e-02_pr, & ! CH4
       2.805376e-02_pr, & ! C2H4
       2.906170e-02_pr, & ! C2H5
       2.603788e-02_pr, & ! C2H2
       3.905682e-02_pr, & ! C3H3
       4.107270e-02_pr, & ! AXC3H5
       3.006964e-02_pr, & ! C2H6
       4.006476e-02_pr, & ! PXC3H4
       4.006476e-02_pr, & ! AXC3H4
       7.811364e-02_pr, & ! A1
       6.509470e-02_pr, & ! C5H5
       4.208064e-02_pr, & ! C3H6
       6.610264e-02_pr, & ! C5H6
       1.281735e-01_pr, & ! A2
       9.214052e-02_pr, & ! A1CH3
       1.021356e-01_pr & ! A1C2H
      !1.402688e-02_pr, & ! SXCH2
      !1.301894e-02_pr, & ! CH
      !2.901834e-02_pr, & ! HCO
      !2.704582e-02_pr, & ! C2H3
      !4.102934e-02_pr, & ! HCCO
      !2.502994e-02_pr, & ! C2H
      !4.308858e-02_pr, & ! NXC3H7
      !1.031436e-01_pr, & ! A1C2H2
      !1.402688e-02_pr, & ! TXCH2
      !1.011277e-01_pr, & ! A1C2
      !1.061240e-01_pr, & ! A1CHO
      !9.113258e-02_pr, & ! A1CH2
      !7.710570e-02_pr, & ! A1X
      !1.271656e-01_pr, & ! A2X
      !5.610752e-02_pr & ! C4H8
       /)
  
contains

  ! ----------------------------------------------- !
  ! 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) :: ntmp,ccoeff,dcoeff,lgknull
    real(pr) :: f
    real(pr) :: conc, concin

    if (concin.gt.0.0_pr) then
       conc = concin
    else
       conc = Ploc / ( Rcst * Tloc )
    end if

    ntmp = 0.75_pr - 1.27_pr * log10( fc )
    ccoeff = - 0.4_pr - 0.67_pr * log10( fc )
    dcoeff = 0.14_pr
    k0 = abs(k0) * conc / max(abs(kinf), tiny(1.0_pr)) + tiny(1.0_pr)
    lgknull = log10(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

  ! ----------------------------------------------- !
  ! Evaluate thirdbodies                            !
  ! ----------------------------------------------- !
  subroutine get_thirdbodies(M,c)
    implicit none

    real(pr), dimension(nspec) :: c
    real(pr), dimension(nTB + nFO) :: M

    M(mM8) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM9) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (5.300000e+00_pr)*c(sH2O) &
         + sum(c)

    M(mM64) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (-1.000000e+00_pr)*c(sH2O) &
         + sum(c)

    M(mM12) = (9.000000e-02_pr)*c(sCO) &
         + (1.180000e+00_pr)*c(sCO2) &
         + (-2.500000e-01_pr)*c(sH2) &
         + (1.089000e+01_pr)*c(sH2O) &
         + (-1.500000e-01_pr)*c(sO2) &
         + sum(c)

    M(mM33) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM36) = (2.000000e+00_pr)*c(sC2H6) &
         + (2.000000e+00_pr)*c(sCH4) &
         + (5.000000e-01_pr)*c(sCO) &
         + (1.000000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (5.000000e+00_pr)*c(sH2O) &
         + sum(c)

    M(mM39) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM55) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM79) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM87) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM97) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM105) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM106) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM112) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM117) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM153) = (2.000000e+00_pr)*c(sC2H6) &
         + (1.000000e+00_pr)*c(sCH4) &
         + (7.500000e-01_pr)*c(sCO) &
         + (2.600000e+00_pr)*c(sCO2) &
         + (1.000000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

  end subroutine get_thirdbodies
  
  ! ----------------------------------------------- !
  ! Evaluate rate coefficients                      !
  ! ----------------------------------------------- !
  subroutine get_rate_coefficients(k,M,Tloc,Ploc)
    implicit none

    real(pr), dimension(nreac + nreac_reverse) :: k
    real(pr), dimension(nFO + nFO_reverse) :: k_0
    real(pr), dimension(nFO + nFO_reverse) :: k_inf
    real(pr), dimension(nFO + nFO_reverse) :: FC
    real(pr), dimension(nTB + nFO) :: M
    real(pr) :: Tloc,Ploc,R_T_inv,T_log

    ! Rate coefficients
    R_T_inv = 1.0_pr/(Rcst*Tloc)
    T_log = log(Tloc)

    k(r1f) = (1.530000e+27_pr)*exp((-3.669242e+05_pr)*R_T_inv + (-2.580000e+00_pr)*T_log )
    k(r2f) = (7.660000e-08_pr)*exp((-5.417862e+03_pr)*R_T_inv + (4.760000e+00_pr)*T_log )
    k(r3f) = (3.170000e-13_pr)*exp((-2.405131e+04_pr)*R_T_inv + (5.950000e+00_pr)*T_log )
    k(r4f) = (2.960000e+03_pr)*exp((-8.920288e+02_pr)*R_T_inv + (1.020000e+00_pr)*T_log )
    k(r5f) = (8.940000e-05_pr)*exp((-3.201178e+03_pr)*R_T_inv + (3.860000e+00_pr)*T_log )
    k(r6f) = (1.500000e+07_pr)*exp((-2.509982e+03_pr)*R_T_inv )
    k(r7f) = (4.590000e-02_pr)*exp((-2.619000e+04_pr)*R_T_inv + (2.700000e+00_pr)*T_log )
    k(r8f) = (9.430000e+06_pr)*exp((-1.000000e+00_pr)*T_log )
    k(r9f) = (4.400000e+10_pr)*exp((-2.000000e+00_pr)*T_log )
    k(r10f) = (1.730000e+02_pr)*exp((-1.434999e+04_pr)*R_T_inv + (1.510000e+00_pr)*T_log )
    k(r11f) = (3.970000e-02_pr)*exp((8.829997e+03_pr)*R_T_inv + (2.400000e+00_pr)*T_log )
    k_0(FOr12f) =(6.330000e+07_pr)*exp((-1.400000e+00_pr)*T_log )
    k_inf(FOr12f) =(5.120000e+06_pr)*exp((4.400000e-01_pr)*T_log )
    FC(FOr12f) = ((1.0_pr - 5.000000e-01_pr)*exp(-Tloc/(0.000000e+00_pr))) + (5.000000e-01_pr*exp(-Tloc/(1.000000e+10_pr)))&
              + exp(-(1.000000e+10_pr)/Tloc)
    k(r12f) = getlindratecoeff(Tloc,k_0(FOr12f),k_inf(FOr12f),FC(FOr12f), M(mM12),Ploc)
    k(r13f) = (2.640000e+10_pr)*exp((-7.130000e+04_pr)*R_T_inv + (-6.700000e-01_pr)*T_log )
    k(r14f) = (3.649000e+00_pr)*exp((4.570016e+03_pr)*R_T_inv + (2.070000e+00_pr)*T_log )
    k(r15f) = (3.970000e+06_pr)*exp((-2.810016e+03_pr)*R_T_inv )
    k(r16f) = (7.490000e+07_pr)*exp((-2.660020e+03_pr)*R_T_inv )
    k(r17f) = (4.000000e+07_pr)
    k(r18f) = (2.380000e+07_pr)*exp((2.089992e+03_pr)*R_T_inv )
    k(r19f) = (1.000000e+10_pr)*exp((-7.251002e+04_pr)*R_T_inv )
    k(r20f) = (5.700000e+07_pr)
    k(r21f) = (5.710000e+06_pr)*exp((3.160008e+03_pr)*R_T_inv )
    k(r22f) = (6.710000e+07_pr)
    k(r23f) = (8.000000e+07_pr)
    k(r24f) = (2.000000e+07_pr)
    k(r25f) = (1.130000e+01_pr)*exp((-1.254999e+04_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r26f) = (5.000000e-01_pr)*exp((-3.024999e+04_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r27f) = (5.800000e+06_pr)*exp((-6.280017e+03_pr)*R_T_inv )
    k(r28f) = (2.400000e+06_pr)*exp((-6.280017e+03_pr)*R_T_inv )
    k(r29f) = (5.000000e+06_pr)*exp((-6.280017e+03_pr)*R_T_inv )
    k(r30f) = (7.000000e+07_pr)
    k(r31f) = (2.800000e+07_pr)
    k(r32f) = (1.200000e+07_pr)
    k_0(FOr33f) =(1.880000e+23_pr)*exp((-2.109000e+04_pr)*R_T_inv + (-6.360000e+00_pr)*T_log )
    k_inf(FOr33f) =(4.820000e+08_pr)*exp((-4.790011e+03_pr)*R_T_inv + (-1.160000e+00_pr)*T_log )
    FC(FOr33f) = ((1.0_pr - 6.027000e-01_pr)*exp(-Tloc/(2.080000e+02_pr))) + (6.027000e-01_pr*exp(-Tloc/(3.922000e+03_pr)))&
              + exp(-(1.018000e+04_pr)/Tloc)
    k(r33f) = getlindratecoeff(Tloc,k_0(FOr33f),k_inf(FOr33f),FC(FOr33f), M(mM33),Ploc)
    k(r34f) = (3.000000e+07_pr)
    k(r35f) = (6.820000e+04_pr)*exp((3.909990e+03_pr)*R_T_inv + (2.500000e-01_pr)*T_log )
    k_0(FOr36f) =(3.470000e+26_pr)*exp((-2.122999e+04_pr)*R_T_inv + (-6.300000e+00_pr)*T_log )
    k_inf(FOr36f) =(6.920000e+07_pr)*exp((1.800000e-01_pr)*T_log )
    k(r36f) = getlindratecoeff(Tloc,k_0(FOr36f),k_inf(FOr36f),1.0_pr, M(mM36),Ploc)
    k(r37f) = (5.060000e+07_pr)
    k(r38f) = (3.370000e+07_pr)
    k_0(FOr39f) =(4.000000e+21_pr)*exp((-1.314002e+04_pr)*R_T_inv + (-5.920000e+00_pr)*T_log )
    k_inf(FOr39f) =(2.790000e+09_pr)*exp((-5.569992e+03_pr)*R_T_inv + (-1.430000e+00_pr)*T_log )
    FC(FOr39f) = ((1.0_pr - 4.120000e-01_pr)*exp(-Tloc/(1.950000e+02_pr))) + (4.120000e-01_pr*exp(-Tloc/(5.900000e+03_pr)))&
              + exp(-(6.394000e+03_pr)/Tloc)
    k(r39f) = getlindratecoeff(Tloc,k_0(FOr39f),k_inf(FOr39f),FC(FOr39f), M(mM39),Ploc)
    k(r40f) = (5.600000e+01_pr)*exp((-2.268000e+04_pr)*R_T_inv + (1.600000e+00_pr)*T_log )
    k(r41f) = (6.440000e+11_pr)*exp((-5.929983e+03_pr)*R_T_inv + (-1.340000e+00_pr)*T_log )
    k(r42f) = (1.380000e+07_pr)*exp((-1.275300e+05_pr)*R_T_inv )
    k(r43f) = (5.870000e+05_pr)*exp((-5.790999e+04_pr)*R_T_inv )
    k(r44f) = (1.000000e+07_pr)
    k(r45f) = (3.610000e+06_pr)
    k(r46f) = (3.000000e+07_pr)
    k(r47f) = (1.000000e+08_pr)
    k(r48f) = (6.840000e+06_pr)*exp((-4.434998e+04_pr)*R_T_inv + (1.000000e-01_pr)*T_log )
    k(r49f) = (2.632000e+06_pr)*exp((-5.716001e+04_pr)*R_T_inv + (-6.000000e-02_pr)*T_log )
    k(r50f) = (6.600000e+02_pr)*exp((-4.536000e+04_pr)*R_T_inv + (1.620000e+00_pr)*T_log )
    k(r51f) = (1.020000e+03_pr)*exp((-3.598002e+04_pr)*R_T_inv + (1.500000e+00_pr)*T_log )
    k(r52f) = (1.000000e+02_pr)*exp((-1.304998e+04_pr)*R_T_inv + (1.600000e+00_pr)*T_log )
    k(r53f) = (6.000000e+07_pr)
    k(r54f) = (2.460000e+00_pr)*exp((-3.460001e+04_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k_0(FOr55f) =(1.170000e+12_pr)*exp((-1.754000e+04_pr)*R_T_inv + (-2.790000e+00_pr)*T_log )
    k_inf(FOr55f) =(1.360000e+04_pr)*exp((-9.980012e+03_pr)*R_T_inv )
    FC(FOr55f) = ((1.0_pr - 1.000000e+00_pr)*exp(-Tloc/(1.000000e+00_pr))) + (1.000000e+00_pr*exp(-Tloc/(1.000000e+07_pr)))&
              + exp(-(1.000000e+07_pr)/Tloc)
    k(r55f) = getlindratecoeff(Tloc,k_0(FOr55f),k_inf(FOr55f),FC(FOr55f), M(mM55),Ploc)
    k(r56f) = (8.000000e+05_pr)*exp((-3.076001e+04_pr)*R_T_inv + (1.400000e-01_pr)*T_log )
    k(r57f) = (8.780000e+04_pr)*exp((6.999832e+01_pr)*R_T_inv + (3.000000e-02_pr)*T_log )
    k(r58f) = (3.010000e+07_pr)*exp((-9.622999e+04_pr)*R_T_inv )
    k(r59f) = (9.000000e+06_pr)
    k(r60f) = (1.200000e+08_pr)
    k(r61f) = (3.000000e+07_pr)
    k(r62f) = (3.000000e+07_pr)
    k(r63f) = (3.020000e+07_pr)
    k(r64f) = (1.870000e+11_pr)*exp((-7.113001e+04_pr)*R_T_inv + (-1.000000e+00_pr)*T_log )
    k(r65f) = (2.240000e+12_pr)*exp((-7.113001e+04_pr)*R_T_inv + (-1.000000e+00_pr)*T_log )
    k(r66f) = (1.200000e+04_pr)*exp((3.040011e+03_pr)*R_T_inv + (8.100000e-01_pr)*T_log )
    k(r67f) = (2.650000e+07_pr)
    k(r68f) = (5.740000e+01_pr)*exp((-1.147002e+04_pr)*R_T_inv + (1.900000e+00_pr)*T_log )
    k(r69f) = (3.900000e+07_pr)*exp((-1.480998e+04_pr)*R_T_inv )
    k(r70f) = (3.430000e+03_pr)*exp((1.869997e+03_pr)*R_T_inv + (1.180000e+00_pr)*T_log )
    k(r71f) = (3.320000e-03_pr)*exp((-2.452000e+04_pr)*R_T_inv + (2.810000e+00_pr)*T_log )
    k(r72f) = (1.900000e+08_pr)*exp((-6.607000e+04_pr)*R_T_inv )
    k(r73f) = (7.000000e+06_pr)
    k(r74f) = (1.400000e+07_pr)
    k(r75f) = (5.000000e+07_pr)
    k(r76f) = (2.000000e+07_pr)
    k(r77f) = (1.000000e+07_pr)*exp((3.160008e+03_pr)*R_T_inv )
    k(r78f) = (3.310000e+00_pr)*exp((-3.769993e+03_pr)*R_T_inv + (2.260000e+00_pr)*T_log )
    k_0(FOr79f) =(6.340000e+19_pr)*exp((-1.582000e+04_pr)*R_T_inv + (-4.660000e+00_pr)*T_log )
    k_inf(FOr79f) =(1.710000e+04_pr)*exp((-1.132998e+04_pr)*R_T_inv + (1.270000e+00_pr)*T_log )
    FC(FOr79f) = ((1.0_pr - 2.122000e-01_pr)*exp(-Tloc/(1.000000e+00_pr))) + (2.122000e-01_pr*exp(-Tloc/(-1.021000e+04_pr)))
    k(r79f) = getlindratecoeff(Tloc,k_0(FOr79f),k_inf(FOr79f),FC(FOr79f), M(mM79),Ploc)
    k(r80f) = (8.100000e+00_pr)*exp((-7.950018e+03_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r81f) = (1.250000e+01_pr)*exp((-7.950018e+03_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r82f) = (3.324000e+10_pr)*exp((-1.284400e+05_pr)*R_T_inv + (-4.400000e-01_pr)*T_log )
    k(r83f) = (2.630000e+00_pr)*exp((-7.138000e+04_pr)*R_T_inv + (2.140000e+00_pr)*T_log )
    k(r84f) = (7.530000e+00_pr)*exp((-8.809998e+03_pr)*R_T_inv + (1.550000e+00_pr)*T_log )
    k(r85f) = (1.280000e+03_pr)*exp((-1.078999e+04_pr)*R_T_inv + (7.300000e-01_pr)*T_log )
    k(r86f) = (1.900000e+08_pr)
    k_0(FOr87f) =(1.400000e+18_pr)*exp((-1.389000e+04_pr)*R_T_inv + (-3.860000e+00_pr)*T_log )
    k_inf(FOr87f) =(6.080000e+06_pr)*exp((-1.170014e+03_pr)*R_T_inv + (2.700000e-01_pr)*T_log )
    FC(FOr87f) = ((1.0_pr - 7.820000e-01_pr)*exp(-Tloc/(2.075000e+02_pr))) + (7.820000e-01_pr*exp(-Tloc/(2.663000e+03_pr)))&
              + exp(-(6.095000e+03_pr)/Tloc)
    k(r87f) = getlindratecoeff(Tloc,k_0(FOr87f),k_inf(FOr87f),FC(FOr87f), M(mM87),Ploc)
    k(r88f) = (3.000000e+07_pr)
    k(r89f) = (1.030000e+07_pr)*exp((1.789999e+03_pr)*R_T_inv + (2.100000e-01_pr)*T_log )
    k(r90f) = (5.000000e+06_pr)
    k(r91f) = (1.340000e+00_pr)*exp((1.610003e+03_pr)*R_T_inv + (1.610000e+00_pr)*T_log )
    k(r92f) = (3.030000e+05_pr)*exp((-4.999880e+01_pr)*R_T_inv + (2.900000e-01_pr)*T_log )
    k(r93f) = (4.580000e+10_pr)*exp((-4.249982e+03_pr)*R_T_inv + (-1.390000e+00_pr)*T_log )
    k(r94f) = (9.000000e+07_pr)
    k(r95f) = (9.030000e+06_pr)*exp((3.200007e+03_pr)*R_T_inv )
    k(r96f) = (1.930000e+12_pr)*exp((-3.208998e+04_pr)*R_T_inv + (-1.250000e+00_pr)*T_log )
    k_0(FOr97f) =(2.030000e+27_pr)*exp((-2.414001e+04_pr)*R_T_inv + (-6.640000e+00_pr)*T_log )
    k_inf(FOr97f) =(1.370000e+03_pr)*exp((-5.669989e+03_pr)*R_T_inv + (1.460000e+00_pr)*T_log )
    FC(FOr97f) = ((1.0_pr - -5.690000e-01_pr)*exp(-Tloc/(2.990000e+02_pr))) + (-5.690000e-01_pr*exp(-Tloc/(-9.147000e+03_pr)))&
              + exp(-(1.524000e+02_pr)/Tloc)
    k(r97f) = getlindratecoeff(Tloc,k_0(FOr97f),k_inf(FOr97f),FC(FOr97f), M(mM97),Ploc)
    k(r98f) = (1.330000e+00_pr)*exp((-5.120998e+04_pr)*R_T_inv + (2.530000e+00_pr)*T_log )
    k(r99f) = (7.660000e+03_pr)*exp((-4.770011e+03_pr)*R_T_inv + (8.800000e-01_pr)*T_log )
    k(r100f) = (7.150000e-02_pr)*exp((-3.889990e+03_pr)*R_T_inv + (2.470000e+00_pr)*T_log )
    k(r101f) = (3.890000e+02_pr)*exp((-3.709995e+03_pr)*R_T_inv + (1.360000e+00_pr)*T_log )
    k(r102f) = (1.310000e-07_pr)*exp((3.599997e+03_pr)*R_T_inv + (4.200000e+00_pr)*T_log )
    k(r103f) = (3.750000e+30_pr)*exp((-2.954000e+04_pr)*R_T_inv + (-7.800000e+00_pr)*T_log )
    k(r104f) = (2.270000e-01_pr)*exp((-3.849000e+04_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k_0(FOr105f) =(3.000000e+51_pr)*exp((-7.602002e+04_pr)*R_T_inv + (-1.460000e+01_pr)*T_log )
    k_inf(FOr105f) =(2.550000e+00_pr)*exp((-2.385001e+04_pr)*R_T_inv + (1.600000e+00_pr)*T_log )
    FC(FOr105f) = ((1.0_pr - 1.894000e-01_pr)*exp(-Tloc/(2.770000e+02_pr))) + (1.894000e-01_pr*exp(-Tloc/(8.748000e+03_pr)))&
              + exp(-(7.891000e+03_pr)/Tloc)
    k(r105f) = getlindratecoeff(Tloc,k_0(FOr105f),k_inf(FOr105f),FC(FOr105f), M(mM105),Ploc)
    k_0(FOr106f) =(1.990000e+29_pr)*exp((-2.797000e+04_pr)*R_T_inv + (-7.080000e+00_pr)*T_log )
    k_inf(FOr106f) =(5.210000e+11_pr)*exp((-6.610009e+03_pr)*R_T_inv + (-9.900000e-01_pr)*T_log )
    FC(FOr106f) = ((1.0_pr - 8.422000e-01_pr)*exp(-Tloc/(1.250000e+02_pr))) + (8.422000e-01_pr*exp(-Tloc/(2.219000e+03_pr)))&
              + exp(-(6.882000e+03_pr)/Tloc)
    k(r106f) = getlindratecoeff(Tloc,k_0(FOr106f),k_inf(FOr106f),FC(FOr106f), M(mM106),Ploc)
    k(r107f) = (3.170000e+07_pr)*exp((1.650002e+03_pr)*R_T_inv + (3.000000e-02_pr)*T_log )
    k(r108f) = (1.920000e+01_pr)*exp((8.510005e+03_pr)*R_T_inv + (1.020000e+00_pr)*T_log )
    k(r109f) = (1.200000e+08_pr)
    k(r110f) = (3.000000e+05_pr)
    k(r111f) = (3.100000e+07_pr)
    k_0(FOr112f) =(3.720000e+59_pr)*exp((-4.250100e+05_pr)*R_T_inv + (-1.314000e+01_pr)*T_log )
    k_inf(FOr112f) =(1.880000e+50_pr)*exp((-4.491200e+05_pr)*R_T_inv + (-9.720000e+00_pr)*T_log )
    FC(FOr112f) = ((1.0_pr - 3.900000e-01_pr)*exp(-Tloc/(1.000000e+02_pr))) + (3.900000e-01_pr*exp(-Tloc/(1.900000e+03_pr)))&
              + exp(-(6.000000e+03_pr)/Tloc)
    k(r112f) = getlindratecoeff(Tloc,k_0(FOr112f),k_inf(FOr112f),FC(FOr112f), M(mM112),Ploc)
    k(r113f) = (1.700000e-01_pr)*exp((-2.402001e+04_pr)*R_T_inv + (2.700000e+00_pr)*T_log )
    k(r114f) = (8.980000e+01_pr)*exp((-2.381001e+04_pr)*R_T_inv + (1.920000e+00_pr)*T_log )
    k(r115f) = (1.610000e+00_pr)*exp((-3.100009e+03_pr)*R_T_inv + (2.220000e+00_pr)*T_log )
    k(r116f) = (8.430000e+08_pr)*exp((-9.311998e+04_pr)*R_T_inv )
    k_0(FOr117f) =(1.213000e+29_pr)*exp((-3.207800e+05_pr)*R_T_inv + (-5.180000e+00_pr)*T_log )
    k_inf(FOr117f) =(2.255000e+20_pr)*exp((-3.126800e+05_pr)*R_T_inv + (-1.440000e+00_pr)*T_log )
    FC(FOr117f) = ((1.0_pr - 5.757000e-01_pr)*exp(-Tloc/(2.370000e+02_pr))) + (5.757000e-01_pr*exp(-Tloc/(1.652000e+03_pr)))&
              + exp(-(5.069000e+03_pr)/Tloc)
    k(r117f) = getlindratecoeff(Tloc,k_0(FOr117f),k_inf(FOr117f),FC(FOr117f), M(mM117),Ploc)
    k(r118f) = (1.000000e+08_pr)
    k(r119f) = (1.000000e+08_pr)
    k(r120f) = (4.200000e+04_pr)*exp((-3.569998e+03_pr)*R_T_inv )
    k(r121f) = (1.000000e+07_pr)
    k(r122f) = (1.000000e+05_pr)*exp((-1.254999e+04_pr)*R_T_inv )
    k(r123f) = (5.000000e+07_pr)
    k(r124f) = (1.000000e+07_pr)
    k(r125f) = (7.940000e+23_pr)*exp((-2.034001e+04_pr)*R_T_inv + (-5.060000e+00_pr)*T_log )
    k(r126f) = (3.160000e+23_pr)*exp((-1.970999e+04_pr)*R_T_inv + (-5.000000e+00_pr)*T_log )
    k(r127f) = (1.280000e+03_pr)*exp((-1.078999e+04_pr)*R_T_inv + (7.300000e-01_pr)*T_log )
    k(r128f) = (6.950000e+07_pr)
    k(r129f) = (1.700000e-01_pr)*exp((-6.280017e+03_pr)*R_T_inv + (1.700000e+00_pr)*T_log )
    k(r130f) = (8.000000e+05_pr)
    k(r131f) = (1.870000e+40_pr)*exp((-7.030999e+04_pr)*R_T_inv + (-9.840000e+00_pr)*T_log )
    k(r132f) = (5.770000e+31_pr)*exp((-1.318200e+05_pr)*R_T_inv + (-7.000000e+00_pr)*T_log )
    k(r133f) = (2.350000e+05_pr)*exp((-4.182000e+04_pr)*R_T_inv )
    k(r134f) = (3.460000e+06_pr)*exp((-2.286000e+04_pr)*R_T_inv + (4.400000e-01_pr)*T_log )
    k(r135f) = (8.500000e-02_pr)*exp((-2.402001e+04_pr)*R_T_inv + (2.700000e+00_pr)*T_log )
    k(r136f) = (8.050000e-01_pr)*exp((-3.100009e+03_pr)*R_T_inv + (2.220000e+00_pr)*T_log )
    k(r137f) = (4.220000e+08_pr)*exp((-9.311998e+04_pr)*R_T_inv )
    k(r138f) = (4.050000e+00_pr)*exp((-7.950018e+03_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r139f) = (6.250000e+00_pr)*exp((-7.950018e+03_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r140f) = (1.280000e+03_pr)*exp((-1.078999e+04_pr)*R_T_inv + (7.300000e-01_pr)*T_log )
    k(r141f) = (8.950000e+07_pr)*exp((-4.707000e+04_pr)*R_T_inv + (-2.000000e-02_pr)*T_log )
    k(r142f) = (7.760000e+39_pr)*exp((-3.282200e+05_pr)*R_T_inv + (-7.800000e+00_pr)*T_log )
    k(r143f) = (2.470000e+09_pr)*exp((-2.692998e+04_pr)*R_T_inv + (-3.300000e-01_pr)*T_log )
    k(r144f) = (1.330000e+00_pr)*exp((-5.120998e+04_pr)*R_T_inv + (2.530000e+00_pr)*T_log )
    k(r145f) = (1.310000e-07_pr)*exp((3.599997e+03_pr)*R_T_inv + (4.200000e+00_pr)*T_log )
    k(r146f) = (2.270000e-01_pr)*exp((-3.849000e+04_pr)*R_T_inv + (2.000000e+00_pr)*T_log )
    k(r147f) = (2.010000e+43_pr)*exp((-8.210000e+04_pr)*R_T_inv + (-1.077000e+01_pr)*T_log )
    k(r148f) = (9.560000e-03_pr)*exp((-1.377000e+04_pr)*R_T_inv + (2.800000e+00_pr)*T_log )
    k(r149f) = (6.030000e+06_pr)
    k(r150f) = (4.860000e+05_pr)*exp((5.499868e+02_pr)*R_T_inv + (-3.200000e-01_pr)*T_log )
    k(r151f) = (9.710000e+14_pr)*exp((-1.045200e+05_pr)*R_T_inv + (-2.700000e+00_pr)*T_log )
    k(r152f) = (2.160000e+33_pr)*exp((-9.979999e+04_pr)*R_T_inv + (-7.740000e+00_pr)*T_log )
    k_0(FOr153f) =(6.260000e+26_pr)*exp((-2.929001e+04_pr)*R_T_inv + (-6.660000e+00_pr)*T_log )
    k_inf(FOr153f) =(3.060000e+08_pr)*exp((-1.687001e+04_pr)*R_T_inv + (-3.700000e-01_pr)*T_log )
    FC(FOr153f) = ((1.0_pr - 1.000000e+00_pr)*exp(-Tloc/(1.000000e+03_pr))) + (1.000000e+00_pr*exp(-Tloc/(1.310000e+03_pr)))&
              + exp(-(4.810000e+04_pr)/Tloc)
    k(r153f) = getlindratecoeff(Tloc,k_0(FOr153f),k_inf(FOr153f),FC(FOr153f), M(mM153),Ploc)
    k(r154f) = (4.040000e+42_pr)*exp((-4.679000e+05_pr)*R_T_inv + (-7.670000e+00_pr)*T_log )
    k(r155f) = (5.930000e+48_pr)*exp((-9.853002e+04_pr)*R_T_inv + (-1.176000e+01_pr)*T_log )
    k(r156f) = (8.000000e+15_pr)*exp((-4.678001e+04_pr)*R_T_inv + (-2.390000e+00_pr)*T_log )
    k(r157f) = (1.200000e+02_pr)*exp((-1.370009e+03_pr)*R_T_inv + (1.600000e+00_pr)*T_log )
    k(r158f) = (3.500000e+01_pr)*exp((4.069986e+03_pr)*R_T_inv + (1.600000e+00_pr)*T_log )
    k(r159f) = (6.600000e-01_pr)*exp((-2.826999e+04_pr)*R_T_inv + (2.540000e+00_pr)*T_log )
    k(r160f) = (9.650000e-02_pr)*exp((-1.555000e+04_pr)*R_T_inv + (2.680000e+00_pr)*T_log )
    k(r161f) = (2.000000e+02_pr)*exp((-2.249988e+03_pr)*R_T_inv + (1.460000e+00_pr)*T_log )
    k(r162f) = (4.520000e-07_pr)*exp((-2.992999e+04_pr)*R_T_inv + (3.650000e+00_pr)*T_log )
    k(r163f) = (7.680000e+06_pr)*exp((-6.190019e+03_pr)*R_T_inv + (1.100000e-01_pr)*T_log )
    k(r164f) = (5.000000e+15_pr)*exp((-2.970600e+05_pr)*R_T_inv )
    k(r165f) = (2.561000e+06_pr)*exp((-1.303998e+04_pr)*R_T_inv + (6.000000e-02_pr)*T_log )
    k(r166f) = (6.390000e+23_pr)*exp((-1.473000e+05_pr)*R_T_inv + (-4.030000e+00_pr)*T_log )
    k(r167f) = (7.000000e+07_pr)
    k(r168f) = (6.776000e+23_pr)*exp((-4.878000e+04_pr)*R_T_inv + (-4.700000e+00_pr)*T_log )
    k(r169f) = (1.190000e+27_pr)*exp((-5.606999e+04_pr)*R_T_inv + (-6.520000e+00_pr)*T_log )
    k(r170f) = (3.020000e+07_pr)
    k(r171f) = (1.730000e+68_pr)*exp((-4.869000e+05_pr)*R_T_inv + (-1.516000e+01_pr)*T_log )
    k(r172f) = (2.800000e+07_pr)*exp((-9.449982e+03_pr)*R_T_inv )
    k(r173f) = (6.600000e+08_pr)*exp((-5.165002e+04_pr)*R_T_inv )
    k(r174f) = (4.770000e-02_pr)*exp((-4.630014e+03_pr)*R_T_inv + (2.710000e+00_pr)*T_log )
    k(r175f) = (3.080000e+00_pr)*exp((2.000000e+00_pr)*T_log )
    k(r176f) = (1.800000e-07_pr)*exp((4.000000e+00_pr)*T_log )
    k(r177f) = (2.600000e+07_pr)*exp((-2.561001e+04_pr)*R_T_inv )
    k(r178f) = (3.000000e+07_pr)*exp((-3.757002e+04_pr)*R_T_inv )
    k(r179f) = (1.000000e+08_pr)
    k(r180f) = (3.000000e+07_pr)
    k(r181f) = (3.000000e+07_pr)
    k(r182f) = (3.290000e+00_pr)*exp((-1.323002e+04_pr)*R_T_inv + (2.050000e+00_pr)*T_log )
    k(r183f) = (1.290000e+61_pr)*exp((-6.195900e+05_pr)*R_T_inv + (-1.248000e+01_pr)*T_log )
    k(r184f) = (6.020000e+02_pr)*exp((-6.841999e+04_pr)*R_T_inv + (1.800000e+00_pr)*T_log )
    k(r185f) = (4.030000e-04_pr)*exp((-6.089979e+03_pr)*R_T_inv + (3.330000e+00_pr)*T_log )
    k(r186f) = (2.752000e-08_pr)*exp((-5.706001e+04_pr)*R_T_inv + (4.460000e+00_pr)*T_log )
    k(r187f) = (2.220000e+07_pr)*exp((-1.896001e+04_pr)*R_T_inv )
    k(r188f) = (1.320000e-04_pr)*exp((-2.338998e+04_pr)*R_T_inv + (3.250000e+00_pr)*T_log )
    k(r189f) = (8.200000e+14_pr)*exp((-3.375500e+05_pr)*R_T_inv )
    k(r190f) = (5.830000e+61_pr)*exp((-2.858900e+05_pr)*R_T_inv + (-1.415000e+01_pr)*T_log )
    k(r191f) = (3.310000e+08_pr)
    k(r192f) = (1.060000e+10_pr)*exp((-1.056000e+04_pr)*R_T_inv + (-9.400000e-01_pr)*T_log )
    k(r193f) = (4.320000e+33_pr)*exp((-9.979999e+04_pr)*R_T_inv + (-7.740000e+00_pr)*T_log )
    k(r194f) = (2.310000e+00_pr)*exp((-1.742000e+04_pr)*R_T_inv + (2.170000e+00_pr)*T_log )
    k(r195f) = (1.250000e+18_pr)*exp((-3.965900e+05_pr)*R_T_inv + (-6.000000e-01_pr)*T_log )
    k(r196f) = (2.160000e+29_pr)*exp((-4.609300e+05_pr)*R_T_inv + (-3.580000e+00_pr)*T_log )
    k(r197f) = (6.470000e-06_pr)*exp((-1.415999e+04_pr)*R_T_inv + (3.980000e+00_pr)*T_log )
    k(r198f) = (1.620000e+07_pr)*exp((-1.159001e+04_pr)*R_T_inv )
    k(r199f) = (4.220000e+08_pr)*exp((-9.311998e+04_pr)*R_T_inv )
    k(r200f) = (2.610000e+15_pr)*exp((-3.370200e+05_pr)*R_T_inv + (1.500000e-01_pr)*T_log )
    k(r201f) = (2.050000e+03_pr)*exp((-1.006001e+04_pr)*R_T_inv + (1.160000e+00_pr)*T_log )
    k(r202f) = (2.720000e+00_pr)*exp((-2.476999e+04_pr)*R_T_inv + (1.770000e+00_pr)*T_log )
    k(r203f) = (1.340000e-02_pr)*exp((-5.369997e+03_pr)*R_T_inv + (2.500000e+00_pr)*T_log )
    k(r204f) = (3.620000e+22_pr)*exp((-9.984999e+04_pr)*R_T_inv + (-4.240000e+00_pr)*T_log )
    k(r205f) = (2.100000e+60_pr)*exp((-6.195500e+05_pr)*R_T_inv + (-1.240000e+01_pr)*T_log )
    k(r206f) = (1.320000e+02_pr)*exp((-7.037998e+04_pr)*R_T_inv + (1.880000e+00_pr)*T_log )
    k(r207f) = (1.340000e-04_pr)*exp((-6.089979e+03_pr)*R_T_inv + (3.330000e+00_pr)*T_log )
    k(r208f) = (1.340000e+17_pr)*exp((-1.725400e+05_pr)*R_T_inv + (-8.600000e-01_pr)*T_log )
    k(r209f) = (7.091000e+07_pr)*exp((-2.930001e+04_pr)*R_T_inv + (-2.600000e-01_pr)*T_log )
    k(r210f) = (8.600000e+60_pr)*exp((-6.195500e+05_pr)*R_T_inv + (-1.248000e+01_pr)*T_log )
    k(r211f) = (2.650000e+02_pr)*exp((-7.153000e+04_pr)*R_T_inv + (1.870000e+00_pr)*T_log )
    k(r212f) = (9.630000e-04_pr)*exp((-1.829998e+04_pr)*R_T_inv + (3.020000e+00_pr)*T_log )
    k(r6b) = (1.580383e+07_pr)*exp((-4.020140e+04_pr)*R_T_inv + (-1.079592e-01_pr)*T_log )
    k(r7b) = (2.333543e-02_pr)*exp((-1.992948e+04_pr)*R_T_inv + (2.684354e+00_pr)*T_log )
    k(r10b) = (6.312372e+02_pr)*exp((-7.624174e+04_pr)*R_T_inv + (1.543273e+00_pr)*T_log )
    k(r11b) = (2.849272e-01_pr)*exp((-5.932228e+04_pr)*R_T_inv + (2.448919e+00_pr)*T_log )
    k(r13b) = (1.158872e+08_pr)*exp((-1.250730e+03_pr)*R_T_inv + (-3.331572e-01_pr)*T_log )
    k(r25b) = (5.411455e+00_pr)*exp((-8.649826e+04_pr)*R_T_inv + (2.230129e+00_pr)*T_log )
    k(r26b) = (3.593488e+01_pr)*exp((-5.958422e+04_pr)*R_T_inv + (1.707023e+00_pr)*T_log )
    k(r30b) = (5.300481e+09_pr)*exp((-6.702565e+04_pr)*R_T_inv + (-4.009361e-01_pr)*T_log )
    k(r34b) = (3.160766e+07_pr)*exp((-3.769142e+04_pr)*R_T_inv + (-1.079592e-01_pr)*T_log )
    k_0(FOr36b) =(3.789535e+33_pr)*exp((-4.593118e+05_pr)*R_T_inv + (-6.014518e+00_pr)*T_log )
    k_inf(FOr36b) =(7.557229e+14_pr)*exp((-4.380818e+05_pr)*R_T_inv + (4.654818e-01_pr)*T_log )
    k(r36b) = getlindratecoeff(Tloc,k_0(FOr36b),k_inf(FOr36b),1.0_pr, M(mM36),Ploc)
    k(r40b) = (2.843075e+00_pr)*exp((-5.523752e+04_pr)*R_T_inv + (1.926250e+00_pr)*T_log )
    k(r41b) = (3.103238e+10_pr)*exp((-7.960817e+02_pr)*R_T_inv + (-9.057907e-01_pr)*T_log )
    k(r48b) = (1.076035e+11_pr)*exp((-4.813845e+03_pr)*R_T_inv + (-7.924392e-01_pr)*T_log )
    k(r50b) = (1.044894e+01_pr)*exp((-4.001105e+04_pr)*R_T_inv + (1.720966e+00_pr)*T_log )
    k(r51b) = (8.209779e+00_pr)*exp((-2.437054e+04_pr)*R_T_inv + (1.585321e+00_pr)*T_log )
    k(r52b) = (5.776631e+00_pr)*exp((-6.959278e+04_pr)*R_T_inv + (1.734239e+00_pr)*T_log )
    k(r54b) = (2.799043e+00_pr)*exp((-5.858529e+04_pr)*R_T_inv + (1.807989e+00_pr)*T_log )
    k(r56b) = (1.022875e+12_pr)*exp((-1.370286e+05_pr)*R_T_inv + (-9.548937e-01_pr)*T_log )
    k(r57b) = (1.122605e+11_pr)*exp((-1.061985e+05_pr)*R_T_inv + (-1.064894e+00_pr)*T_log )
    k(r59b) = (9.482297e+06_pr)*exp((-3.769142e+04_pr)*R_T_inv + (-1.079592e-01_pr)*T_log )
    k(r64b) = (1.915154e+05_pr)*exp((-7.122291e+03_pr)*R_T_inv + (-9.772162e-01_pr)*T_log )
    k(r65b) = (2.294089e+06_pr)*exp((-7.122291e+03_pr)*R_T_inv + (-9.772162e-01_pr)*T_log )
    k(r68b) = (8.712558e-01_pr)*exp((-7.723870e+04_pr)*R_T_inv + (2.074405e+00_pr)*T_log )
    k(r73b) = (7.375120e+06_pr)*exp((-3.769142e+04_pr)*R_T_inv + (-1.079592e-01_pr)*T_log )
    k(r78b) = (2.495642e+03_pr)*exp((-1.258202e+05_pr)*R_T_inv + (1.833330e+00_pr)*T_log )
    k_0(FOr79b) =(1.116668e+24_pr)*exp((-1.599926e+05_pr)*R_T_inv + (-4.173371e+00_pr)*T_log )
    k_inf(FOr79b) =(3.011833e+08_pr)*exp((-1.555026e+05_pr)*R_T_inv + (1.756629e+00_pr)*T_log )
    FC(FOr79b) = ((1.0_pr - 2.122000e-01_pr)*exp(-Tloc/(1.000000e+00_pr))) + (2.122000e-01_pr*exp(-Tloc/(-1.021000e+04_pr)))
    k(r79b) = getlindratecoeff(Tloc,k_0(FOr79b),k_inf(FOr79b),FC(FOr79b), M(mM79),Ploc)
    k(r86b) = (1.669998e+13_pr)*exp((-1.047121e+05_pr)*R_T_inv + (-1.303670e+00_pr)*T_log )
    k(r96b) = (2.654199e+20_pr)*exp((-9.984391e+04_pr)*R_T_inv + (-3.011639e+00_pr)*T_log )
    k_0(FOr97b) =(1.021455e+33_pr)*exp((-1.740761e+05_pr)*R_T_inv + (-6.599194e+00_pr)*T_log )
    k_inf(FOr97b) =(6.893563e+08_pr)*exp((-1.556060e+05_pr)*R_T_inv + (1.500806e+00_pr)*T_log )
    FC(FOr97b) = ((1.0_pr - -5.690000e-01_pr)*exp(-Tloc/(2.990000e+02_pr))) + (-5.690000e-01_pr*exp(-Tloc/(-9.147000e+03_pr)))&
              + exp(-(1.524000e+02_pr)/Tloc)
    k(r97b) = getlindratecoeff(Tloc,k_0(FOr97b),k_inf(FOr97b),FC(FOr97b), M(mM97),Ploc)
    k_0(FOr105b) =(2.733940e+63_pr)*exp((-1.771235e+05_pr)*R_T_inv + (-1.600830e+01_pr)*T_log )
    k_inf(FOr105b) =(2.323849e+12_pr)*exp((-1.249535e+05_pr)*R_T_inv + (1.917019e-01_pr)*T_log )
    FC(FOr105b) = ((1.0_pr - 1.894000e-01_pr)*exp(-Tloc/(2.770000e+02_pr))) + (1.894000e-01_pr*exp(-Tloc/(8.748000e+03_pr)))&
              + exp(-(7.891000e+03_pr)/Tloc)
    k(r105b) = getlindratecoeff(Tloc,k_0(FOr105b),k_inf(FOr105b),FC(FOr105b), M(mM105),Ploc)
    k_0(FOr106b) =(2.222499e+37_pr)*exp((-4.471117e+05_pr)*R_T_inv + (-6.990355e+00_pr)*T_log )
    k_inf(FOr106b) =(5.818704e+19_pr)*exp((-4.257517e+05_pr)*R_T_inv + (-9.003551e-01_pr)*T_log )
    FC(FOr106b) = ((1.0_pr - 8.422000e-01_pr)*exp(-Tloc/(1.250000e+02_pr))) + (8.422000e-01_pr*exp(-Tloc/(2.219000e+03_pr)))&
              + exp(-(6.882000e+03_pr)/Tloc)
    k(r106b) = getlindratecoeff(Tloc,k_0(FOr106b),k_inf(FOr106b),FC(FOr106b), M(mM106),Ploc)
    k_0(FOr112b) =(2.117309e+47_pr)*exp((-4.540446e+04_pr)*R_T_inv + (-1.233721e+01_pr)*T_log )
    k_inf(FOr112b) =(1.070038e+38_pr)*exp((-6.951447e+04_pr)*R_T_inv + (-8.917206e+00_pr)*T_log )
    FC(FOr112b) = ((1.0_pr - 3.900000e-01_pr)*exp(-Tloc/(1.000000e+02_pr))) + (3.900000e-01_pr*exp(-Tloc/(1.900000e+03_pr)))&
              + exp(-(6.000000e+03_pr)/Tloc)
    k(r112b) = getlindratecoeff(Tloc,k_0(FOr112b),k_inf(FOr112b),FC(FOr112b), M(mM112),Ploc)
    k(r125b) = (4.065840e+30_pr)*exp((-3.901881e+05_pr)*R_T_inv + (-4.697097e+00_pr)*T_log )
    k(r126b) = (5.083203e+30_pr)*exp((-3.844538e+05_pr)*R_T_inv + (-4.714468e+00_pr)*T_log )
    k(r131b) = (1.462677e+56_pr)*exp((-6.688943e+05_pr)*R_T_inv + (-1.036019e+01_pr)*T_log )
    k(r132b) = (3.006920e+39_pr)*exp((-2.568508e+05_pr)*R_T_inv + (-7.603729e+00_pr)*T_log )
    k(r133b) = (2.112284e+18_pr)*exp((-3.486701e+05_pr)*R_T_inv + (-5.142990e-01_pr)*T_log )
    k(r134b) = (1.006436e+02_pr)*exp((-4.805837e+04_pr)*R_T_inv + (1.366279e+00_pr)*T_log )
    k(r135b) = (2.869948e-03_pr)*exp((-8.690478e+04_pr)*R_T_inv + (2.723545e+00_pr)*T_log )
    k(r141b) = (8.287287e+02_pr)*exp((-7.737267e+04_pr)*R_T_inv + (9.836499e-01_pr)*T_log )
    k(r142b) = (2.470251e+39_pr)*exp((-3.333243e+05_pr)*R_T_inv + (-7.722629e+00_pr)*T_log )
    k(r143b) = (7.862783e+08_pr)*exp((-3.203428e+04_pr)*R_T_inv + (-2.526289e-01_pr)*T_log )
    k(r147b) = (4.508130e+50_pr)*exp((-3.243292e+05_pr)*R_T_inv + (-1.104136e+01_pr)*T_log )
    k_0(FOr153b) =(7.188448e+32_pr)*exp((-1.647529e+05_pr)*R_T_inv + (-6.644673e+00_pr)*T_log )
    k_inf(FOr153b) =(3.513842e+14_pr)*exp((-1.523329e+05_pr)*R_T_inv + (-3.546725e-01_pr)*T_log )
    FC(FOr153b) = ((1.0_pr - 1.000000e+00_pr)*exp(-Tloc/(1.000000e+03_pr))) + (1.000000e+00_pr*exp(-Tloc/(1.310000e+03_pr)))&
              + exp(-(4.810000e+04_pr)/Tloc)
    k(r153b) = getlindratecoeff(Tloc,k_0(FOr153b),k_inf(FOr153b),FC(FOr153b), M(mM153),Ploc)
    k(r154b) = (1.267709e+29_pr)*exp((-3.849329e+04_pr)*R_T_inv + (-6.436628e+00_pr)*T_log )
    k(r155b) = (1.374169e+54_pr)*exp((-4.601828e+05_pr)*R_T_inv + (-1.123173e+01_pr)*T_log )
    k(r156b) = (1.008052e+10_pr)*exp((-8.113948e+04_pr)*R_T_inv + (-9.663744e-01_pr)*T_log )
    k(r171b) = (2.838062e+60_pr)*exp((-1.410233e+05_pr)*R_T_inv + (-1.491858e+01_pr)*T_log )
    k(r182b) = (6.758064e+13_pr)*exp((-1.899816e+05_pr)*R_T_inv + (4.493911e-01_pr)*T_log )
    k(r183b) = (8.594663e+52_pr)*exp((-1.460365e+05_pr)*R_T_inv + (-1.256354e+01_pr)*T_log )
    k(r189b) = (2.710074e+00_pr)*exp((-5.013494e+04_pr)*R_T_inv + (1.726656e+00_pr)*T_log )
    k(r190b) = (1.696682e+53_pr)*exp((-2.190202e+05_pr)*R_T_inv + (-1.219183e+01_pr)*T_log )
    k(r194b) = (1.405955e-05_pr)*exp((-5.206747e+04_pr)*R_T_inv + (3.435357e+00_pr)*T_log )
    k(r195b) = (1.741713e+13_pr)*exp((-2.455379e+04_pr)*R_T_inv + (-1.376359e+00_pr)*T_log )
    k(r196b) = (8.758954e+15_pr)*exp((-2.202397e+04_pr)*R_T_inv + (-2.398187e+00_pr)*T_log )
    k(r203b) = (1.642035e+12_pr)*exp((-3.811843e+05_pr)*R_T_inv + (1.641569e+00_pr)*T_log )
    k(r204b) = (7.792154e+32_pr)*exp((-3.788659e+05_pr)*R_T_inv + (-5.716326e+00_pr)*T_log )
    k(r205b) = (3.266888e+52_pr)*exp((-1.353278e+05_pr)*R_T_inv + (-1.245977e+01_pr)*T_log )
    k(r206b) = (3.550369e-01_pr)*exp((-1.889066e+04_pr)*R_T_inv + (2.206678e+00_pr)*T_log )
    k(r207b) = (1.315076e-06_pr)*exp((-1.649241e+04_pr)*R_T_inv + (3.689951e+00_pr)*T_log )
    k(r208b) = (3.763237e+10_pr)*exp((-4.273620e+04_pr)*R_T_inv + (-7.565933e-01_pr)*T_log )
    k(r209b) = (4.716479e+15_pr)*exp((-3.022288e+05_pr)*R_T_inv + (-1.015807e+00_pr)*T_log )
    k(r210b) = (6.922092e+52_pr)*exp((-1.084095e+05_pr)*R_T_inv + (-1.253899e+01_pr)*T_log )
    k(r211b) = (3.687818e-01_pr)*exp((6.877639e+03_pr)*R_T_inv + (2.197460e+00_pr)*T_log )
    k(r212b) = (4.889860e-06_pr)*exp((-1.784092e+03_pr)*R_T_inv + (3.380733e+00_pr)*T_log )
  
    return
  end subroutine get_rate_coefficients

  ! ----------------------------------------------- !
  ! Evaluate 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 + nreac_reverse) :: w,k
    real(pr), dimension(nTB + nFO) :: m

    w(r1f) = k(r1f) * c(sPOSF10325) 
    w(r2f) = k(r2f) * c(sH) * c(sPOSF10325) 
    w(r3f) = k(r3f) * c(sCH3) * c(sPOSF10325) 
    w(r4f) = k(r4f) * c(sOH) * c(sPOSF10325) 
    w(r5f) = k(r5f) * c(sO) * c(sPOSF10325) 
    w(r6f) = k(r6f) * c(sN2) * cqss(sqssSXCH2) 
    w(r7f) = k(r7f) * c(sH2) * c(sO) 
    w(r8f) = k(r8f) * c(sH) * c(sO) * m(mM8)
    w(r9f) = k(r9f) * c(sH) * c(sOH) * m(mM9)
    w(r10f) = k(r10f) * c(sH2) * c(sOH) 
    w(r11f) = k(r11f) * c(sOH)**2.0_pr 
    w(r12f) = k(r12f) * c(sH) * c(sO2) 
    w(r13f) = k(r13f) * c(sH) * c(sO2) 
    w(r14f) = k(r14f) * c(sH) * c(sHO2) 
    w(r15f) = k(r15f) * c(sH) * c(sHO2) 
    w(r16f) = k(r16f) * c(sH) * c(sHO2) 
    w(r17f) = k(r17f) * c(sHO2) * c(sO) 
    w(r18f) = k(r18f) * c(sHO2) * c(sOH) 
    w(r19f) = k(r19f) * c(sHO2) * c(sOH) 
    w(r20f) = k(r20f) * cqss(sqssCH) * c(sO) 
    w(r21f) = k(r21f) * cqss(sqssCH) * c(sH2O) 
    w(r22f) = k(r22f) * cqss(sqssCH) * c(sO2) 
    w(r23f) = k(r23f) * c(sO) * cqss(sqssTXCH2) 
    w(r24f) = k(r24f) * c(sOH) * cqss(sqssTXCH2) 
    w(r25f) = k(r25f) * c(sOH) * cqss(sqssTXCH2) 
    w(r26f) = k(r26f) * c(sH2) * cqss(sqssTXCH2) 
    w(r27f) = k(r27f) * c(sO2) * cqss(sqssTXCH2) 
    w(r28f) = k(r28f) * c(sO2) * cqss(sqssTXCH2) 
    w(r29f) = k(r29f) * c(sO2) * cqss(sqssTXCH2) 
    w(r30f) = k(r30f) * c(sH2) * cqss(sqssSXCH2) 
    w(r31f) = k(r31f) * c(sO2) * cqss(sqssSXCH2) 
    w(r32f) = k(r32f) * c(sO2) * cqss(sqssSXCH2) 
    w(r33f) = k(r33f) * c(sH2O) * cqss(sqssSXCH2) 
    w(r34f) = k(r34f) * c(sH2O) * cqss(sqssSXCH2) 
    w(r35f) = k(r35f) * c(sH2O) * cqss(sqssSXCH2) 
    w(r36f) = k(r36f) * c(sCH3) * c(sH) 
    w(r37f) = k(r37f) * c(sCH3) * c(sO) 
    w(r38f) = k(r38f) * c(sCH3) * c(sO) 
    w(r39f) = k(r39f) * c(sCH3) * c(sOH) 
    w(r40f) = k(r40f) * c(sCH3) * c(sOH) 
    w(r41f) = k(r41f) * c(sCH3) * c(sOH) 
    w(r42f) = k(r42f) * c(sCH3) * c(sO2) 
    w(r43f) = k(r43f) * c(sCH3) * c(sO2) 
    w(r44f) = k(r44f) * c(sCH3) * c(sHO2) 
    w(r45f) = k(r45f) * c(sCH3) * c(sHO2) 
    w(r46f) = k(r46f) * cqss(sqssCH) * c(sCH3) 
    w(r47f) = k(r47f) * c(sCH3) * cqss(sqssTXCH2) 
    w(r48f) = k(r48f) * c(sCH3)**2.0_pr 
    w(r49f) = k(r49f) * c(sCH3)**2.0_pr 
    w(r50f) = k(r50f) * c(sCH4) * c(sH) 
    w(r51f) = k(r51f) * c(sCH4) * c(sO) 
    w(r52f) = k(r52f) * c(sCH4) * c(sOH) 
    w(r53f) = k(r53f) * cqss(sqssCH) * c(sCH4) 
    w(r54f) = k(r54f) * c(sCH4) * cqss(sqssTXCH2) 
    w(r55f) = k(r55f) * c(sCO) * c(sO) 
    w(r56f) = k(r56f) * c(sCO) * c(sOH) 
    w(r57f) = k(r57f) * c(sCO) * c(sOH) 
    w(r58f) = k(r58f) * c(sCO) * c(sHO2) 
    w(r59f) = k(r59f) * c(sCO) * cqss(sqssSXCH2) 
    w(r60f) = k(r60f) * c(sH) * cqss(sqssHCO) 
    w(r61f) = k(r61f) * cqss(sqssHCO) * c(sO) 
    w(r62f) = k(r62f) * cqss(sqssHCO) * c(sO) 
    w(r63f) = k(r63f) * cqss(sqssHCO) * c(sOH) 
    w(r64f) = k(r64f) * cqss(sqssHCO) * m(mM64)
    w(r65f) = k(r65f) * c(sH2O) * cqss(sqssHCO) 
    w(r66f) = k(r66f) * cqss(sqssHCO) * c(sO2) 
    w(r67f) = k(r67f) * c(sCH3) * cqss(sqssHCO) 
    w(r68f) = k(r68f) * c(sCH2O) * c(sH) 
    w(r69f) = k(r69f) * c(sCH2O) * c(sO) 
    w(r70f) = k(r70f) * c(sCH2O) * c(sOH) 
    w(r71f) = k(r71f) * c(sCH2O) * c(sCH3) 
    w(r72f) = k(r72f) * cqss(sqssCH) * c(sCO2) 
    w(r73f) = k(r73f) * c(sCO2) * cqss(sqssSXCH2) 
    w(r74f) = k(r74f) * c(sCO2) * cqss(sqssSXCH2) 
    w(r75f) = k(r75f) * cqss(sqssC2H) * c(sO) 
    w(r76f) = k(r76f) * cqss(sqssC2H) * c(sOH) 
    w(r77f) = k(r77f) * cqss(sqssC2H) * c(sO2) 
    w(r78f) = k(r78f) * cqss(sqssC2H) * c(sH2) 
    w(r79f) = k(r79f) * c(sC2H2) * c(sH) 
    w(r80f) = k(r80f) * c(sC2H2) * c(sO) 
    w(r81f) = k(r81f) * c(sC2H2) * c(sO) 
    w(r82f) = k(r82f) * c(sC2H2) * c(sO) 
    w(r83f) = k(r83f) * c(sC2H2) * c(sOH) 
    w(r84f) = k(r84f) * c(sC2H2) * c(sOH) 
    w(r85f) = k(r85f) * c(sC2H2) * c(sOH) 
    w(r86f) = k(r86f) * c(sC2H2) * cqss(sqssSXCH2) 
    w(r87f) = k(r87f) * cqss(sqssC2H3) * c(sH) 
    w(r88f) = k(r88f) * cqss(sqssC2H3) * c(sH) 
    w(r89f) = k(r89f) * cqss(sqssC2H3) * c(sO) 
    w(r90f) = k(r90f) * cqss(sqssC2H3) * c(sOH) 
    w(r91f) = k(r91f) * cqss(sqssC2H3) * c(sO2) 
    w(r92f) = k(r92f) * cqss(sqssC2H3) * c(sO2) 
    w(r93f) = k(r93f) * cqss(sqssC2H3) * c(sO2) 
    w(r94f) = k(r94f) * cqss(sqssC2H3) * cqss(sqssHCO) 
    w(r95f) = k(r95f) * cqss(sqssC2H3) * c(sCH3) 
    w(r96f) = k(r96f) * cqss(sqssC2H3) * c(sCH3) 
    w(r97f) = k(r97f) * c(sC2H4) * c(sH) 
    w(r98f) = k(r98f) * c(sC2H4) * c(sH) 
    w(r99f) = k(r99f) * c(sC2H4) * c(sO) 
    w(r100f) = k(r100f) * c(sC2H4) * c(sO) 
    w(r101f) = k(r101f) * c(sC2H4) * c(sO) 
    w(r102f) = k(r102f) * c(sC2H4) * c(sOH) 
    w(r103f) = k(r103f) * c(sC2H4) * c(sOH) 
    w(r104f) = k(r104f) * c(sC2H4) * c(sCH3) 
    w(r105f) = k(r105f) * c(sC2H4) * c(sCH3) 
    w(r106f) = k(r106f) * c(sC2H5) * c(sH) 
    w(r107f) = k(r107f) * c(sC2H5) * c(sO) 
    w(r108f) = k(r108f) * c(sC2H5) * c(sO2) 
    w(r109f) = k(r109f) * c(sC2H5) * cqss(sqssHCO) 
    w(r110f) = k(r110f) * c(sC2H5) * c(sHO2) 
    w(r111f) = k(r111f) * c(sC2H5) * c(sHO2) 
    w(r112f) = k(r112f) * c(sC2H6) 
    w(r113f) = k(r113f) * c(sC2H6) * c(sH) 
    w(r114f) = k(r114f) * c(sC2H6) * c(sO) 
    w(r115f) = k(r115f) * c(sC2H6) * c(sOH) 
    w(r116f) = k(r116f) * c(sC2H6) * c(sCH3) 
    w(r117f) = k(r117f) * cqss(sqssHCCO) 
    w(r118f) = k(r118f) * c(sH) * cqss(sqssHCCO) 
    w(r119f) = k(r119f) * cqss(sqssHCCO) * c(sO) 
    w(r120f) = k(r120f) * cqss(sqssHCCO) * c(sO2) 
    w(r121f) = k(r121f) * cqss(sqssHCCO)**2.0_pr 
    w(r122f) = k(r122f) * c(sC2H2) * cqss(sqssHCCO) 
    w(r123f) = k(r123f) * c(sCH3) * cqss(sqssHCCO) 
    w(r124f) = k(r124f) * cqss(sqssHCCO) * c(sOH) 
    w(r125f) = k(r125f) * c(sC3H3) * c(sH) 
    w(r126f) = k(r126f) * c(sC3H3) * c(sH) 
    w(r127f) = k(r127f) * c(sC3H3) * c(sOH) 
    w(r128f) = k(r128f) * c(sC3H3) * c(sO) 
    w(r129f) = k(r129f) * c(sC3H3) * c(sO2) 
    w(r130f) = k(r130f) * c(sC3H3) * c(sHO2) 
    w(r131f) = k(r131f) * c(sC3H3)**2.0_pr 
    w(r132f) = k(r132f) * c(sC3H3)**2.0_pr 
    w(r133f) = k(r133f) * c(sC2H2) * c(sC3H3) 
    w(r134f) = k(r134f) * c(sH) * c(sPXC3H4) 
    w(r135f) = k(r135f) * c(sH) * c(sPXC3H4) 
    w(r136f) = k(r136f) * c(sOH) * c(sPXC3H4) 
    w(r137f) = k(r137f) * c(sCH3) * c(sPXC3H4) 
    w(r138f) = k(r138f) * c(sO) * c(sPXC3H4) 
    w(r139f) = k(r139f) * c(sO) * c(sPXC3H4) 
    w(r140f) = k(r140f) * c(sOH) * c(sPXC3H4) 
    w(r141f) = k(r141f) * c(sAXC3H4) * c(sH) 
    w(r142f) = k(r142f) * c(sAXC3H4) 
    w(r143f) = k(r143f) * c(sAXC3H4) * c(sH) 
    w(r144f) = k(r144f) * c(sAXC3H4) * c(sH) 
    w(r145f) = k(r145f) * c(sAXC3H4) * c(sOH) 
    w(r146f) = k(r146f) * c(sAXC3H4) * c(sCH3) 
    w(r147f) = k(r147f) * c(sAXC3H4) * c(sH) 
    w(r148f) = k(r148f) * c(sAXC3H5) * c(sH) 
    w(r149f) = k(r149f) * c(sAXC3H5) * c(sOH) 
    w(r150f) = k(r150f) * c(sAXC3H5) * c(sCH3) 
    w(r151f) = k(r151f) * c(sAXC3H5) * c(sO2) 
    w(r152f) = k(r152f) * c(sAXC3H5) * c(sC3H3) 
    w(r153f) = k(r153f) * c(sC3H6) * c(sH) 
    w(r154f) = k(r154f) * c(sC3H6) 
    w(r155f) = k(r155f) * c(sAXC3H5) * c(sH) 
    w(r156f) = k(r156f) * c(sC3H6) * c(sH) 
    w(r157f) = k(r157f) * c(sC3H6) * c(sO) 
    w(r158f) = k(r158f) * c(sC3H6) * c(sO) 
    w(r159f) = k(r159f) * c(sC3H6) * c(sH) 
    w(r160f) = k(r160f) * c(sC3H6) * c(sO) 
    w(r161f) = k(r161f) * c(sC3H6) * c(sOH) 
    w(r162f) = k(r162f) * c(sC3H6) * c(sCH3) 
    w(r163f) = k(r163f) * cqss(sqssC4H8) * c(sH) 
    w(r164f) = k(r164f) * cqss(sqssC4H8) 
    w(r165f) = k(r165f) * c(sC5H5) * c(sHO2) 
    w(r166f) = k(r166f) * c(sC5H5)**2.0_pr 
    w(r167f) = k(r167f) * c(sC5H5) * c(sO) 
    w(r168f) = k(r168f) * c(sC5H5) * c(sHO2) 
    w(r169f) = k(r169f) * c(sC5H5) * c(sHO2) 
    w(r170f) = k(r170f) * c(sC5H5) * c(sOH) 
    w(r171f) = k(r171f) * c(sC5H6) 
    w(r172f) = k(r172f) * c(sC5H6) * c(sH) 
    w(r173f) = k(r173f) * c(sC5H6) * c(sH) 
    w(r174f) = k(r174f) * c(sC5H6) * c(sO) 
    w(r175f) = k(r175f) * c(sC5H6) * c(sOH) 
    w(r176f) = k(r176f) * c(sC5H6) * c(sCH3) 
    w(r177f) = k(r177f) * cqss(sqssA1X) * c(sO2) 
    w(r178f) = k(r178f) * cqss(sqssA1X) * c(sO2) 
    w(r179f) = k(r179f) * cqss(sqssA1X) * c(sO) 
    w(r180f) = k(r180f) * cqss(sqssA1X) * c(sOH) 
    w(r181f) = k(r181f) * cqss(sqssA1X) * c(sHO2) 
    w(r182f) = k(r182f) * cqss(sqssA1X) * c(sC2H2) 
    w(r183f) = k(r183f) * c(sA1) 
    w(r184f) = k(r184f) * c(sA1) * c(sH) 
    w(r185f) = k(r185f) * c(sA1) * c(sOH) 
    w(r186f) = k(r186f) * c(sA1) * c(sCH3) 
    w(r187f) = k(r187f) * c(sA1) * c(sO) 
    w(r188f) = k(r188f) * c(sA1) * c(sOH) 
    w(r189f) = k(r189f) * cqss(sqssA1CH2) 
    w(r190f) = k(r190f) * cqss(sqssA1CH2) * c(sH) 
    w(r191f) = k(r191f) * cqss(sqssA1CH2) * c(sO) 
    w(r192f) = k(r192f) * cqss(sqssA1CH2) * c(sHO2) 
    w(r193f) = k(r193f) * cqss(sqssA1CH2) * c(sC3H3) 
    w(r194f) = k(r194f) * c(sA1CH3) * c(sH) 
    w(r195f) = k(r195f) * c(sA1CH3) 
    w(r196f) = k(r196f) * c(sA1CH3) 
    w(r197f) = k(r197f) * c(sA1CH3) * c(sH) 
    w(r198f) = k(r198f) * c(sA1CH3) * c(sOH) 
    w(r199f) = k(r199f) * c(sA1CH3) * c(sCH3) 
    w(r200f) = k(r200f) * cqss(sqssA1CHO) 
    w(r201f) = k(r201f) * cqss(sqssA1CHO) * c(sH) 
    w(r202f) = k(r202f) * cqss(sqssA1CHO) * c(sCH3) 
    w(r203f) = k(r203f) * cqss(sqssA1C2) * c(sC2H2) 
    w(r204f) = k(r204f) * cqss(sqssA1C2) * c(sC2H4) 
    w(r205f) = k(r205f) * c(sA1C2H) 
    w(r206f) = k(r206f) * c(sA1C2H) * c(sH) 
    w(r207f) = k(r207f) * c(sA1C2H) * c(sOH) 
    w(r208f) = k(r208f) * cqss(sqssA1C2H2) 
    w(r209f) = k(r209f) * cqss(sqssA1C2H2) * c(sC2H2) 
    w(r210f) = k(r210f) * c(sA2) 
    w(r211f) = k(r211f) * c(sA2) * c(sH) 
    w(r212f) = k(r212f) * c(sA2) * c(sOH) 
    w(r6b) = k(r6b) * c(sN2) * cqss(sqssTXCH2) 
    w(r7b) = k(r7b) * c(sH) * c(sOH) 
    w(r10b) = k(r10b) * c(sH) * c(sH2O) 
    w(r11b) = k(r11b) * c(sH2O) * c(sO) 
    w(r13b) = k(r13b) * c(sO) * c(sOH) 
    w(r25b) = k(r25b) * cqss(sqssCH) * c(sH2O) 
    w(r26b) = k(r26b) * c(sCH3) * c(sH) 
    w(r30b) = k(r30b) * c(sCH3) * c(sH) 
    w(r34b) = k(r34b) * c(sH2O) * cqss(sqssTXCH2) 
    w(r36b) = k(r36b) * c(sCH4) 
    w(r40b) = k(r40b) * c(sH2O) * cqss(sqssTXCH2) 
    w(r41b) = k(r41b) * c(sH2O) * cqss(sqssSXCH2) 
    w(r48b) = k(r48b) * c(sC2H5) * c(sH) 
    w(r50b) = k(r50b) * c(sCH3) * c(sH2) 
    w(r51b) = k(r51b) * c(sCH3) * c(sOH) 
    w(r52b) = k(r52b) * c(sCH3) * c(sH2O) 
    w(r54b) = k(r54b) * c(sCH3)**2.0_pr 
    w(r56b) = k(r56b) * c(sCO2) * c(sH) 
    w(r57b) = k(r57b) * c(sCO2) * c(sH) 
    w(r59b) = k(r59b) * c(sCO) * cqss(sqssTXCH2) 
    w(r64b) = k(r64b) * c(sCO) * c(sH) * m(mM64)
    w(r65b) = k(r65b) * c(sCO) * c(sH) * c(sH2O) 
    w(r68b) = k(r68b) * c(sH2) * cqss(sqssHCO) 
    w(r73b) = k(r73b) * c(sCO2) * cqss(sqssTXCH2) 
    w(r78b) = k(r78b) * c(sC2H2) * c(sH) 
    w(r79b) = k(r79b) * cqss(sqssC2H3) 
    w(r86b) = k(r86b) * c(sC3H3) * c(sH) 
    w(r96b) = k(r96b) * c(sAXC3H5) * c(sH) 
    w(r97b) = k(r97b) * c(sC2H5) 
    w(r105b) = k(r105b) * cqss(sqssNXC3H7) 
    w(r106b) = k(r106b) * c(sC2H6) 
    w(r112b) = k(r112b) * c(sCH3)**2.0_pr 
    w(r125b) = k(r125b) * c(sPXC3H4) 
    w(r126b) = k(r126b) * c(sAXC3H4) 
    w(r131b) = k(r131b) * c(sA1) 
    w(r132b) = k(r132b) * cqss(sqssA1X) * c(sH) 
    w(r133b) = k(r133b) * c(sC5H5) 
    w(r134b) = k(r134b) * c(sC2H2) * c(sCH3) 
    w(r135b) = k(r135b) * c(sC3H3) * c(sH2) 
    w(r141b) = k(r141b) * c(sC2H2) * c(sCH3) 
    w(r142b) = k(r142b) * c(sPXC3H4) 
    w(r143b) = k(r143b) * c(sH) * c(sPXC3H4) 
    w(r147b) = k(r147b) * c(sAXC3H5) 
    w(r153b) = k(r153b) * cqss(sqssNXC3H7) 
    w(r154b) = k(r154b) * cqss(sqssC2H3) * c(sCH3) 
    w(r155b) = k(r155b) * c(sC3H6) 
    w(r156b) = k(r156b) * c(sC2H4) * c(sCH3) 
    w(r171b) = k(r171b) * c(sC5H5) * c(sH) 
    w(r182b) = k(r182b) * cqss(sqssA1C2H2) 
    w(r183b) = k(r183b) * cqss(sqssA1X) * c(sH) 
    w(r189b) = k(r189b) * c(sC2H2) * c(sC5H5) 
    w(r190b) = k(r190b) * cqss(sqssA1X) * c(sCH3) 
    w(r194b) = k(r194b) * c(sA1) * c(sCH3) 
    w(r195b) = k(r195b) * cqss(sqssA1CH2) * c(sH) 
    w(r196b) = k(r196b) * cqss(sqssA1X) * c(sCH3) 
    w(r203b) = k(r203b) * cqss(sqssA2X) 
    w(r204b) = k(r204b) * c(sA2) * c(sH) 
    w(r205b) = k(r205b) * cqss(sqssA1C2) * c(sH) 
    w(r206b) = k(r206b) * cqss(sqssA1C2) * c(sH2) 
    w(r207b) = k(r207b) * cqss(sqssA1C2) * c(sH2O) 
    w(r208b) = k(r208b) * c(sA1C2H) * c(sH) 
    w(r209b) = k(r209b) * c(sA2) * c(sH) 
    w(r210b) = k(r210b) * cqss(sqssA2X) * c(sH) 
    w(r211b) = k(r211b) * cqss(sqssA2X) * c(sH2) 
    w(r212b) = k(r212b) * cqss(sqssA2X) * c(sH2O) 
  
    return
  end subroutine get_reaction_rates
  
  ! ----------------------------------------------- !
  ! Evaluate production rates                       !
  ! ----------------------------------------------- !
  subroutine get_production_rates(cdot,w)
    implicit none

    real(pr), dimension(nspec) :: cdot
    real(pr), dimension(nreac + nreac_reverse) :: w

    cdot(sPOSF10325) = 0.0_pr &
         - w(r1f) &
         - w(r2f) &
         - w(r3f) &
         - w(r4f) &
         - w(r5f) 

    cdot(sN2) = 0.0_pr 

    cdot(sO) = 0.0_pr &
         - w(r5f) &
         - w(r7f) &
         + w(r7b) &
         - w(r8f) &
         + w(r11f) &
         - w(r11b) &
         + w(r13f) &
         - w(r13b) &
         + w(r15f) &
         - w(r17f) &
         - w(r20f) &
         + w(r22f) &
         - w(r23f) &
         + w(r28f) &
         - w(r37f) &
         - w(r38f) &
         + w(r42f) &
         - w(r51f) &
         + w(r51b) &
         - w(r55f) &
         - w(r61f) &
         - w(r62f) &
         - w(r69f) &
         - w(r75f) &
         - w(r80f) &
         - w(r81f) &
         - w(r82f) &
         - w(r89f) &
         + w(r92f) &
         - w(r99f) &
         - w(r100f) &
         - w(r101f) &
         - w(r107f) &
         - w(r114f) &
         - w(r119f) &
         - w(r128f) &
         - w(r138f) &
         - w(r139f) &
         - w(r157f) &
         - w(r158f) &
         - w(r160f) &
         - w(r167f) &
         - w(r174f) &
         + w(r177f) &
         - w(r179f) &
         - w(r187f) &
         - w(r191f) 

    cdot(sH2) = 0.0_pr &
         + w(r2f) &
         - w(r7f) &
         + w(r7b) &
         - w(r10f) &
         + w(r10b) &
         + w(r14f) &
         - w(r26f) &
         + w(r26b) &
         - w(r30f) &
         + w(r30b) &
         + w(r33f) &
         + w(r35f) &
         + w(r38f) &
         + w(r39f) &
         + w(r50f) &
         - w(r50b) &
         + w(r60f) &
         + w(r68f) &
         - w(r68b) &
         - w(r78f) &
         + w(r78b) &
         + w(r88f) &
         + w(r98f) &
         + w(r113f) &
         + w(r135f) &
         - w(r135b) &
         + w(r144f) &
         + w(r148f) &
         + w(r159f) &
         + w(r172f) &
         + w(r184f) &
         + w(r197f) &
         + w(r201f) &
         + w(r206f) &
         - w(r206b) &
         + w(r211f) &
         - w(r211b) 

    cdot(sH) = 0.0_pr &
         + 0.5_pr * w(r1f) &
         -0.7_pr * w(r2f) &
         + 0.3_pr * w(r3f) &
         + 0.3_pr * w(r4f) &
         + 0.3_pr * w(r5f) &
         + w(r7f) &
         - w(r7b) &
         - w(r8f) &
         - w(r9f) &
         + w(r10f) &
         - w(r10b) &
         - w(r12f) &
         - w(r13f) &
         + w(r13b) &
         - w(r14f) &
         - w(r15f) &
         - w(r16f) &
         + w(r20f) &
         + w(r21f) &
         + w(r23f) &
         + w(r24f) &
         + w(r26f) &
         - w(r26b) &
         + 2.0_pr * w(r27f) &
         + w(r29f) &
         + w(r30f) &
         - w(r30b) &
         + w(r31f) &
         - w(r36f) &
         + w(r36b) &
         + w(r37f) &
         + w(r38f) &
         + w(r42f) &
         + w(r44f) &
         + w(r46f) &
         + w(r47f) &
         + w(r48f) &
         - w(r48b) &
         - w(r50f) &
         + w(r50b) &
         + w(r53f) &
         + w(r56f) &
         - w(r56b) &
         + w(r57f) &
         - w(r57b) &
         - w(r60f) &
         + w(r62f) &
         + w(r64f) &
         - w(r64b) &
         + w(r65f) &
         - w(r65b) &
         - w(r68f) &
         + w(r68b) &
         + w(r76f) &
         + w(r78f) &
         - w(r78b) &
         - w(r79f) &
         + w(r79b) &
         + w(r80f) &
         + w(r86f) &
         - w(r86b) &
         - w(r87f) &
         - w(r88f) &
         + w(r96f) &
         - w(r96b) &
         - w(r97f) &
         + w(r97b) &
         - w(r98f) &
         + w(r99f) &
         - w(r106f) &
         + w(r106b) &
         - w(r113f) &
         - w(r118f) &
         + w(r119f) &
         - w(r125f) &
         + w(r125b) &
         - w(r126f) &
         + w(r126b) &
         + w(r128f) &
         + w(r132f) &
         - w(r132b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r141f) &
         + w(r141b) &
         - w(r144f) &
         - w(r147f) &
         + w(r147b) &
         - w(r148f) &
         + 2.0_pr * w(r152f) &
         - w(r153f) &
         + w(r153b) &
         - w(r155f) &
         + w(r155b) &
         - w(r156f) &
         + w(r156b) &
         - w(r159f) &
         - w(r163f) &
         + 2.0_pr * w(r166f) &
         + w(r170f) &
         + w(r171f) &
         - w(r171b) &
         - w(r172f) &
         - w(r173f) &
         + w(r178f) &
         + w(r183f) &
         - w(r183b) &
         - w(r184f) &
         + w(r187f) &
         + w(r188f) &
         - w(r190f) &
         + w(r190b) &
         + w(r191f) &
         + w(r192f) &
         + 2.0_pr * w(r193f) &
         - w(r194f) &
         + w(r194b) &
         + w(r195f) &
         - w(r195b) &
         - w(r197f) &
         - w(r201f) &
         + w(r204f) &
         - w(r204b) &
         + w(r205f) &
         - w(r205b) &
         - w(r206f) &
         + w(r206b) &
         + w(r208f) &
         - w(r208b) &
         + w(r209f) &
         - w(r209b) &
         + w(r210f) &
         - w(r210b) &
         - w(r211f) &
         + w(r211b) 

    cdot(sOH) = 0.0_pr &
         - w(r4f) &
         + w(r5f) &
         + w(r7f) &
         - w(r7b) &
         + w(r8f) &
         - w(r9f) &
         - w(r10f) &
         + w(r10b) &
         -2.0_pr * w(r11f) &
         + 2.0_pr * w(r11b) &
         + w(r13f) &
         - w(r13b) &
         + 2.0_pr * w(r16f) &
         + w(r17f) &
         - w(r18f) &
         - w(r19f) &
         - w(r24f) &
         - w(r25f) &
         + w(r25b) &
         + w(r29f) &
         + w(r31f) &
         - w(r39f) &
         - w(r40f) &
         + w(r40b) &
         - w(r41f) &
         + w(r41b) &
         + w(r43f) &
         + w(r44f) &
         + w(r51f) &
         - w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r56f) &
         + w(r56b) &
         - w(r57f) &
         + w(r57b) &
         + w(r58f) &
         + w(r61f) &
         - w(r63f) &
         + w(r69f) &
         - w(r70f) &
         - w(r76f) &
         + w(r82f) &
         - w(r83f) &
         - w(r84f) &
         - w(r85f) &
         - w(r90f) &
         - w(r102f) &
         - w(r103f) &
         + w(r111f) &
         + w(r114f) &
         - w(r115f) &
         + w(r120f) &
         - w(r124f) &
         - w(r127f) &
         + w(r130f) &
         - w(r136f) &
         - w(r140f) &
         - w(r145f) &
         - w(r149f) &
         + w(r151f) &
         + w(r160f) &
         - w(r161f) &
         + w(r168f) &
         - w(r170f) &
         + w(r174f) &
         - w(r175f) &
         - w(r180f) &
         + w(r181f) &
         - w(r185f) &
         - w(r188f) &
         + w(r192f) &
         - w(r198f) &
         - w(r207f) &
         + w(r207b) &
         - w(r212f) &
         + w(r212b) 

    cdot(sH2O) = 0.0_pr &
         + w(r4f) &
         + w(r9f) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         + w(r15f) &
         + w(r18f) &
         + w(r19f) &
         - w(r21f) &
         + w(r25f) &
         - w(r25b) &
         + w(r32f) &
         - w(r33f) &
         - w(r35f) &
         + w(r40f) &
         - w(r40b) &
         + w(r41f) &
         - w(r41b) &
         + w(r52f) &
         - w(r52b) &
         + w(r63f) &
         + w(r70f) &
         + w(r83f) &
         + w(r90f) &
         + w(r102f) &
         + w(r115f) &
         + w(r136f) &
         + w(r145f) &
         + w(r149f) &
         + w(r161f) &
         + w(r169f) &
         + w(r175f) &
         + w(r185f) &
         + w(r198f) &
         + w(r207f) &
         - w(r207b) &
         + w(r212f) &
         - w(r212b) 

    cdot(sO2) = 0.0_pr &
         - w(r12f) &
         - w(r13f) &
         + w(r13b) &
         + w(r14f) &
         + w(r17f) &
         + w(r18f) &
         + w(r19f) &
         - w(r22f) &
         - w(r27f) &
         - w(r28f) &
         - w(r29f) &
         - w(r31f) &
         - w(r32f) &
         - w(r42f) &
         - w(r43f) &
         + w(r45f) &
         - w(r66f) &
         - w(r77f) &
         - w(r91f) &
         - w(r92f) &
         - w(r93f) &
         - w(r108f) &
         + w(r110f) &
         - w(r120f) &
         - w(r129f) &
         - w(r151f) &
         + w(r165f) &
         - w(r177f) &
         - w(r178f) 

    cdot(sHO2) = 0.0_pr &
         + w(r12f) &
         - w(r14f) &
         - w(r15f) &
         - w(r16f) &
         - w(r17f) &
         - w(r18f) &
         - w(r19f) &
         - w(r44f) &
         - w(r45f) &
         - w(r58f) &
         + w(r66f) &
         + w(r91f) &
         + w(r108f) &
         - w(r110f) &
         - w(r111f) &
         - w(r130f) &
         - w(r165f) &
         - w(r168f) &
         - w(r169f) &
         - w(r181f) &
         - w(r192f) 

    cdot(sCO) = 0.0_pr &
         + w(r20f) &
         + w(r29f) &
         + w(r31f) &
         + w(r32f) &
         + w(r38f) &
         - w(r55f) &
         - w(r56f) &
         + w(r56b) &
         - w(r57f) &
         + w(r57b) &
         - w(r58f) &
         + w(r60f) &
         + w(r61f) &
         + w(r63f) &
         + w(r64f) &
         - w(r64b) &
         + w(r65f) &
         - w(r65b) &
         + w(r66f) &
         + w(r67f) &
         + w(r72f) &
         + w(r74f) &
         + w(r75f) &
         + w(r77f) &
         + w(r81f) &
         + w(r84f) &
         + w(r85f) &
         + w(r89f) &
         + w(r92f) &
         + w(r94f) &
         + w(r99f) &
         + w(r109f) &
         + w(r117f) &
         + w(r118f) &
         + 2.0_pr * w(r119f) &
         + 2.0_pr * w(r120f) &
         + 2.0_pr * w(r121f) &
         + w(r122f) &
         + w(r123f) &
         + w(r127f) &
         + w(r128f) &
         + 2.0_pr * w(r129f) &
         + w(r130f) &
         + w(r139f) &
         + w(r140f) &
         + w(r157f) &
         + w(r167f) &
         + w(r168f) &
         + w(r169f) &
         + w(r170f) &
         + w(r177f) &
         + 2.0_pr * w(r178f) &
         + w(r179f) &
         + w(r180f) &
         + w(r181f) &
         + w(r187f) &
         + w(r188f) &
         + w(r201f) &
         + w(r202f) 

    cdot(sCH2O) = 0.0_pr &
         + w(r21f) &
         + w(r24f) &
         + w(r28f) &
         + w(r33f) &
         + w(r35f) &
         + w(r37f) &
         + w(r39f) &
         + w(r42f) &
         + w(r43f) &
         + w(r44f) &
         - w(r68f) &
         + w(r68b) &
         - w(r69f) &
         - w(r70f) &
         - w(r71f) &
         + w(r74f) &
         + w(r93f) &
         + w(r100f) &
         + w(r103f) &
         + w(r107f) &
         + w(r111f) &
         + w(r151f) 

    cdot(sCH3) = 0.0_pr &
         + 1.5_pr * w(r1f) &
         + 0.7_pr * w(r2f) &
         -0.30000000000000004_pr * w(r3f) &
         + 0.7_pr * w(r4f) &
         + 0.7_pr * w(r5f) &
         + w(r26f) &
         - w(r26b) &
         + w(r30f) &
         - w(r30b) &
         - w(r36f) &
         + w(r36b) &
         - w(r37f) &
         - w(r38f) &
         - w(r39f) &
         - w(r40f) &
         + w(r40b) &
         - w(r41f) &
         + w(r41b) &
         - w(r42f) &
         - w(r43f) &
         - w(r44f) &
         - w(r45f) &
         - w(r46f) &
         - w(r47f) &
         -2.0_pr * w(r48f) &
         + 2.0_pr * w(r48b) &
         -2.0_pr * w(r49f) &
         + w(r50f) &
         - w(r50b) &
         + w(r51f) &
         - w(r51b) &
         + w(r52f) &
         - w(r52b) &
         + 2.0_pr * w(r54f) &
         - 2.0_pr * w(r54b) &
         - w(r67f) &
         - w(r71f) &
         + w(r84f) &
         + w(r85f) &
         + w(r89f) &
         + w(r92f) &
         - w(r95f) &
         - w(r96f) &
         + w(r96b) &
         + w(r99f) &
         + w(r101f) &
         + w(r103f) &
         - w(r104f) &
         - w(r105f) &
         + w(r105b) &
         + w(r107f) &
         + w(r111f) &
         + 2.0_pr * w(r112f) &
         - 2.0_pr * w(r112b) &
         - w(r116f) &
         - w(r123f) &
         + w(r129f) &
         + w(r134f) &
         - w(r134b) &
         - w(r137f) &
         + w(r138f) &
         + w(r141f) &
         - w(r141b) &
         - w(r146f) &
         - w(r150f) &
         + w(r154f) &
         - w(r154b) &
         + w(r156f) &
         - w(r156b) &
         + 2.0_pr * w(r157f) &
         - w(r162f) &
         + w(r164f) &
         - w(r176f) &
         - w(r186f) &
         + w(r190f) &
         - w(r190b) &
         + w(r194f) &
         - w(r194b) &
         + w(r196f) &
         - w(r196b) &
         - w(r199f) &
         - w(r202f) 

    cdot(sCO2) = 0.0_pr &
         + w(r27f) &
         + w(r55f) &
         + w(r56f) &
         - w(r56b) &
         + w(r57f) &
         - w(r57b) &
         + w(r58f) &
         + w(r62f) &
         - w(r72f) &
         - w(r74f) 

    cdot(sCH4) = 0.0_pr &
         + 0.45_pr * w(r2f) &
         + 1.45_pr * w(r3f) &
         + 0.45_pr * w(r4f) &
         + 0.45_pr * w(r5f) &
         + w(r36f) &
         - w(r36b) &
         + w(r45f) &
         + w(r49f) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         - w(r54f) &
         + w(r54b) &
         + w(r67f) &
         + w(r71f) &
         + w(r95f) &
         + w(r104f) &
         + w(r116f) &
         + w(r137f) &
         + w(r146f) &
         + w(r150f) &
         + w(r162f) &
         + w(r176f) &
         + w(r186f) &
         + w(r199f) &
         + w(r202f) 

    cdot(sC2H4) = 0.0_pr &
         + 1.74268_pr * w(r1f) &
         + 1.59458_pr * w(r2f) &
         + 1.59458_pr * w(r3f) &
         + 1.59458_pr * w(r4f) &
         + 1.59458_pr * w(r5f) &
         + w(r47f) &
         + w(r53f) &
         + w(r87f) &
         + w(r94f) &
         - w(r97f) &
         + w(r97b) &
         - w(r98f) &
         - w(r99f) &
         - w(r100f) &
         - w(r101f) &
         - w(r102f) &
         - w(r103f) &
         - w(r104f) &
         - w(r105f) &
         + w(r105b) &
         + w(r108f) &
         + w(r123f) &
         + w(r127f) &
         + w(r139f) &
         + w(r156f) &
         - w(r156b) &
         + w(r163f) &
         - w(r204f) &
         + w(r204b) 

    cdot(sC2H5) = 0.0_pr &
         + w(r48f) &
         - w(r48b) &
         + w(r97f) &
         - w(r97b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         - w(r108f) &
         - w(r109f) &
         - w(r110f) &
         - w(r111f) &
         + w(r113f) &
         + w(r114f) &
         + w(r115f) &
         + w(r116f) &
         + w(r140f) &
         + w(r158f) &
         + w(r163f) 

    cdot(sC2H2) = 0.0_pr &
         + w(r78f) &
         - w(r78b) &
         - w(r79f) &
         + w(r79b) &
         - w(r80f) &
         - w(r81f) &
         - w(r82f) &
         - w(r83f) &
         - w(r84f) &
         - w(r85f) &
         - w(r86f) &
         + w(r86b) &
         + w(r88f) &
         + w(r90f) &
         + w(r91f) &
         + w(r95f) &
         + w(r121f) &
         - w(r122f) &
         + w(r128f) &
         - w(r133f) &
         + w(r133b) &
         + w(r134f) &
         - w(r134b) &
         + w(r141f) &
         - w(r141b) &
         + w(r151f) &
         + w(r167f) &
         + w(r168f) &
         + 2.0_pr * w(r169f) &
         + w(r170f) &
         + w(r173f) &
         + 2.0_pr * w(r178f) &
         - w(r182f) &
         + w(r182b) &
         + w(r189f) &
         - w(r189b) &
         - w(r203f) &
         + w(r203b) &
         - w(r209f) &
         + w(r209b) 

    cdot(sC3H3) = 0.0_pr &
         + w(r86f) &
         - w(r86b) &
         + w(r122f) &
         - w(r125f) &
         + w(r125b) &
         - w(r126f) &
         + w(r126b) &
         - w(r127f) &
         - w(r128f) &
         - w(r129f) &
         - w(r130f) &
         -2.0_pr * w(r131f) &
         + 2.0_pr * w(r131b) &
         -2.0_pr * w(r132f) &
         + 2.0_pr * w(r132b) &
         - w(r133f) &
         + w(r133b) &
         + w(r135f) &
         - w(r135b) &
         + w(r136f) &
         + w(r137f) &
         + w(r144f) &
         + w(r145f) &
         + w(r146f) &
         - w(r152f) &
         - w(r193f) 

    cdot(sAXC3H5) = 0.0_pr &
         + w(r96f) &
         - w(r96b) &
         + w(r147f) &
         - w(r147b) &
         - w(r148f) &
         - w(r149f) &
         - w(r150f) &
         - w(r151f) &
         - w(r152f) &
         - w(r155f) &
         + w(r155b) &
         + w(r159f) &
         + w(r160f) &
         + w(r161f) &
         + w(r162f) &
         + w(r164f) &
         + w(r173f) 

    cdot(sC2H6) = 0.0_pr &
         + w(r106f) &
         - w(r106b) &
         + w(r109f) &
         + w(r110f) &
         - w(r112f) &
         + w(r112b) &
         - w(r113f) &
         - w(r114f) &
         - w(r115f) &
         - w(r116f) 

    cdot(sPXC3H4) = 0.0_pr &
         + w(r125f) &
         - w(r125b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r136f) &
         - w(r137f) &
         - w(r138f) &
         - w(r139f) &
         - w(r140f) &
         + w(r142f) &
         - w(r142b) &
         + w(r143f) &
         - w(r143b) 

    cdot(sAXC3H4) = 0.0_pr &
         + w(r126f) &
         - w(r126b) &
         - w(r141f) &
         + w(r141b) &
         - w(r142f) &
         + w(r142b) &
         - w(r143f) &
         + w(r143b) &
         - w(r144f) &
         - w(r145f) &
         - w(r146f) &
         - w(r147f) &
         + w(r147b) &
         + w(r148f) &
         + w(r149f) &
         + w(r150f) 

    cdot(sA1) = 0.0_pr &
         + 0.17_pr * w(r1f) &
         + 0.2465_pr * w(r2f) &
         + 0.2465_pr * w(r3f) &
         + 0.2465_pr * w(r4f) &
         + 0.2465_pr * w(r5f) &
         + w(r131f) &
         - w(r131b) &
         + w(r152f) &
         - w(r183f) &
         + w(r183b) &
         - w(r184f) &
         - w(r185f) &
         - w(r186f) &
         - w(r187f) &
         - w(r188f) &
         + w(r194f) &
         - w(r194b) 

    cdot(sC5H5) = 0.0_pr &
         + w(r133f) &
         - w(r133b) &
         - w(r165f) &
         -2.0_pr * w(r166f) &
         - w(r167f) &
         - w(r168f) &
         - w(r169f) &
         - w(r170f) &
         + w(r171f) &
         - w(r171b) &
         + w(r172f) &
         + w(r174f) &
         + w(r175f) &
         + w(r176f) &
         + w(r177f) &
         + w(r179f) &
         + w(r181f) &
         + w(r187f) &
         + w(r189f) &
         - w(r189b) 

    cdot(sC3H6) = 0.0_pr &
         + 0.819058_pr * w(r1f) &
         + 0.749451_pr * w(r2f) &
         + 0.749451_pr * w(r3f) &
         + 0.749451_pr * w(r4f) &
         + 0.749451_pr * w(r5f) &
         - w(r153f) &
         + w(r153b) &
         - w(r154f) &
         + w(r154b) &
         + w(r155f) &
         - w(r155b) &
         - w(r156f) &
         + w(r156b) &
         - w(r157f) &
         - w(r158f) &
         - w(r159f) &
         - w(r160f) &
         - w(r161f) &
         - w(r162f) 

    cdot(sC5H6) = 0.0_pr &
         + w(r165f) &
         - w(r171f) &
         + w(r171b) &
         - w(r172f) &
         - w(r173f) &
         - w(r174f) &
         - w(r175f) &
         - w(r176f) &
         + w(r180f) &
         + w(r188f) 

    cdot(sA2) = 0.0_pr &
         + w(r166f) &
         + w(r193f) &
         + w(r204f) &
         - w(r204b) &
         + w(r209f) &
         - w(r209b) &
         - w(r210f) &
         + w(r210b) &
         - w(r211f) &
         + w(r211b) &
         - w(r212f) &
         + w(r212b) 

    cdot(sA1CH3) = 0.0_pr &
         + 0.163333_pr * w(r1f) &
         + 0.236833_pr * w(r2f) &
         + 0.236833_pr * w(r3f) &
         + 0.236833_pr * w(r4f) &
         + 0.236833_pr * w(r5f) &
         - w(r194f) &
         + w(r194b) &
         - w(r195f) &
         + w(r195b) &
         - w(r196f) &
         + w(r196b) &
         - w(r197f) &
         - w(r198f) &
         - w(r199f) 

    cdot(sA1C2H) = 0.0_pr &
         - w(r205f) &
         + w(r205b) &
         - w(r206f) &
         + w(r206b) &
         - w(r207f) &
         + w(r207b) &
         + w(r208f) &
         - w(r208b) 

  
    return
  end subroutine get_production_rates


  ! ----------------------------------------------- !
  ! 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 + nreac_reverse) :: k
    real(pr), dimension(nTB + nFO) :: M
    
    integer :: index

    real(pr) :: SXCH2_ct
    real(pr) :: SXCH2_num
    real(pr) :: SXCH2_denom
    real(pr) :: SXCH2_HCCO
    real(pr) :: SXCH2_TXCH2
    real(pr) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_HCCO
    real(pr) :: CH_C2H
    real(pr) :: CH_TXCH2
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_CH
    real(pr) :: HCO_C2H3
    real(pr) :: HCO_HCCO
    real(pr) :: HCO_C2H
    real(pr) :: HCO_TXCH2
    real(pr) :: HCO_A1CHO
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_CH
    real(pr) :: HCCO_ct
    real(pr) :: HCCO_num
    real(pr) :: HCCO_denom
    real(pr) :: HCCO_C2H
    real(pr) :: C2H_ct
    real(pr) :: C2H_num
    real(pr) :: C2H_denom
    real(pr) :: NXC3H7_ct
    real(pr) :: NXC3H7_num
    real(pr) :: NXC3H7_denom
    real(pr) :: A1C2H2_ct
    real(pr) :: A1C2H2_num
    real(pr) :: A1C2H2_denom
    real(pr) :: A1C2H2_A1X
    real(pr) :: TXCH2_ct
    real(pr) :: TXCH2_num
    real(pr) :: TXCH2_denom
    real(pr) :: TXCH2_SXCH2
    real(pr) :: TXCH2_CH
    real(pr) :: A1C2_ct
    real(pr) :: A1C2_num
    real(pr) :: A1C2_denom
    real(pr) :: A1C2_A2X
    real(pr) :: A1CHO_ct
    real(pr) :: A1CHO_num
    real(pr) :: A1CHO_denom
    real(pr) :: A1CHO_A1CH2
    real(pr) :: A1CH2_ct
    real(pr) :: A1CH2_num
    real(pr) :: A1CH2_denom
    real(pr) :: A1CH2_A1X
    real(pr) :: A1X_ct
    real(pr) :: A1X_num
    real(pr) :: A1X_denom
    real(pr) :: A1X_A1C2H2
    real(pr) :: A1X_A1CHO
    real(pr) :: A1X_A1CH2
    real(pr) :: A2X_ct
    real(pr) :: A2X_num
    real(pr) :: A2X_denom
    real(pr) :: A2X_A1C2
    real(pr) :: C4H8_ct
    real(pr) :: C4H8_num
    real(pr) :: C4H8_denom
  
    real(pr) :: a_f_a
    real(pr) :: a_j_a
    real(pr) :: a_j_b
    real(pr) :: a_n_b
    real(pr) :: a_o_c
    real(pr) :: b_c_a
    real(pr) :: c_e_a
    real(pr) :: c_h_a
    real(pr) :: c_h_b
    real(pr) :: A_A_A
    real(pr) :: A_B_A
    real(pr) :: A_C_A
    real(pr) :: A_B_B
    real(pr) :: A_C_B
    real(pr) :: A_C_C
    real(pr) :: B_A_A
    real(pr) :: C_A_A
    real(pr) :: C_B_A
    real(pr) :: C_B_B

    C2H_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r75f)* c(sO) &
             + k(r76f)* c(sOH) &
             + k(r77f)* c(sO2) &
             + k(r78f)* c(sH2)   )

    C2H_num = ( 0.0_pr &
             + k(r82f)* c(sC2H2) * c(sO) &
             + k(r83f)* c(sC2H2) * c(sOH) &
             + k(r78b)* c(sC2H2) * c(sH)   )

    C2H_ct = C2H_num / C2H_denom

    cqss(sqssC2H) = C2H_ct  

    NXC3H7_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r105b)&
             + k(r153b)  )

    NXC3H7_num = ( 0.0_pr &
             + k(r105f)* c(sC2H4) * c(sCH3) &
             + k(r153f)* c(sC3H6) * c(sH)   )

    NXC3H7_ct = NXC3H7_num / NXC3H7_denom

    cqss(sqssNXC3H7) = NXC3H7_ct  

    A1CHO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r200f)&
             + k(r201f)* c(sH) &
             + k(r202f)* c(sCH3)   )

    A1CHO_num = ( 0.0_pr   )

    A1CHO_ct = A1CHO_num / A1CHO_denom

    A1CHO_A1CH2 = - ( 0.0_pr &
          + k(r191f) * c(sO) &
          + k(r192f) * c(sHO2)  ) / A1CHO_denom

    A1X_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r177f)* c(sO2) &
             + k(r178f)* c(sO2) &
             + k(r179f)* c(sO) &
             + k(r180f)* c(sOH) &
             + k(r181f)* c(sHO2) &
             + k(r182f)* c(sC2H2) &
             + k(r132b)* c(sH) &
             + k(r183b)* c(sH) &
             + k(r190b)* c(sCH3) &
             + k(r196b)* c(sCH3)   )

    A1X_num = ( 0.0_pr &
             + k(r132f) *c(sC3H3)** 2.0_pr&
             + k(r183f)* c(sA1) &
             + k(r184f)* c(sA1) * c(sH) &
             + k(r185f)* c(sA1) * c(sOH) &
             + k(r186f)* c(sA1) * c(sCH3) &
             + k(r196f)* c(sA1CH3)   )

    A1X_ct = A1X_num / A1X_denom

    A1X_A1C2H2 = - ( 0.0_pr &
          + k(r182b)  ) / A1X_denom

    A1X_A1CHO = - ( 0.0_pr &
          + k(r200f) &
          + k(r201f) * c(sH) &
          + k(r202f) * c(sCH3)  ) / A1X_denom

    A1X_A1CH2 = - ( 0.0_pr &
          + k(r190f) * c(sH)  ) / A1X_denom

    A1CH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r189f)&
             + k(r190f)* c(sH) &
             + k(r191f)* c(sO) &
             + k(r192f)* c(sHO2) &
             + k(r193f)* c(sC3H3) &
             + k(r195b)* c(sH)   )

    A1CH2_num = ( 0.0_pr &
             + k(r195f)* c(sA1CH3) &
             + k(r197f)* c(sA1CH3) * c(sH) &
             + k(r198f)* c(sA1CH3) * c(sOH) &
             + k(r199f)* c(sA1CH3) * c(sCH3) &
             + k(r189b)* c(sC2H2) * c(sC5H5)   )

    A1CH2_ct = A1CH2_num / A1CH2_denom

    A1CH2_A1X = - ( 0.0_pr &
          + k(r190b) * c(sCH3)  ) / A1CH2_denom

    A1C2H2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r208f)&
             + k(r209f)* c(sC2H2) &
             + k(r182b)  )

    A1C2H2_num = ( 0.0_pr &
             + k(r208b)* c(sA1C2H) * c(sH) &
             + k(r209b)* c(sA2) * c(sH)   )

    A1C2H2_ct = A1C2H2_num / A1C2H2_denom

    A1C2H2_A1X = - ( 0.0_pr &
          + k(r182f) * c(sC2H2)  ) / A1C2H2_denom

    a_f_a = A1CH2_A1X
    a_j_a =  (1.0_pr) - (A1C2H2_A1X) * (A1X_A1C2H2)
    a_j_b =  (a_j_a) - (a_f_a) * (A1X_A1CH2)
    a_n_b =  - (a_f_a) * (A1CHO_A1CH2)
    a_o_c =  (1.0_pr) - (A1X_A1CHO) * (a_n_b) / (a_j_b)
    A_A_A = A1CH2_ct
    A_B_A =  (A1X_ct) - (A1C2H2_ct) * (A1X_A1C2H2)
    A_C_A = A1CHO_ct
    A_B_B =  (A_B_A) - (A_A_A) * (A1X_A1CH2)
    A_C_B =  (A_C_A) - (A_A_A) * (A1CHO_A1CH2)
    A_C_C =  (A_C_B) - (A_B_B) * (a_n_b) / (a_j_b)

    cqss(sqssA1CHO) = ( A_C_C ) / ( a_o_c )  

    cqss(sqssA1X) = (A_B_B - (A1X_A1CHO) * cqss(sqssA1CHO)) / (a_j_b)  

    cqss(sqssA1CH2) = A_A_A - (a_f_a) * cqss(sqssA1X)  

    cqss(sqssA1C2H2) = A1C2H2_ct - (A1C2H2_A1X) * cqss(sqssA1X)  

    A2X_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r203b)&
             + k(r210b)* c(sH) &
             + k(r211b)* c(sH2) &
             + k(r212b)* c(sH2O)   )

    A2X_num = ( 0.0_pr &
             + k(r210f)* c(sA2) &
             + k(r211f)* c(sA2) * c(sH) &
             + k(r212f)* c(sA2) * c(sOH)   )

    A2X_ct = A2X_num / A2X_denom

    A2X_A1C2 = - ( 0.0_pr &
          + k(r203f) * c(sC2H2)  ) / A2X_denom

    A1C2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r203f)* c(sC2H2) &
             + k(r204f)* c(sC2H4) &
             + k(r205b)* c(sH) &
             + k(r206b)* c(sH2) &
             + k(r207b)* c(sH2O)   )

    A1C2_num = ( 0.0_pr &
             + k(r205f)* c(sA1C2H) &
             + k(r206f)* c(sA1C2H) * c(sH) &
             + k(r207f)* c(sA1C2H) * c(sOH) &
             + k(r204b)* c(sA2) * c(sH)   )

    A1C2_ct = A1C2_num / A1C2_denom

    A1C2_A2X = - ( 0.0_pr &
          + k(r203b)  ) / A1C2_denom

    b_c_a =  (1.0_pr) - (A1C2_A2X) * (A2X_A1C2)
    B_A_A =  (A2X_ct) - (A1C2_ct) * (A2X_A1C2)

    cqss(sqssA2X) = ( B_A_A ) / ( b_c_a )  

    cqss(sqssA1C2) = A1C2_ct - (A1C2_A2X) * cqss(sqssA2X)  

    C4H8_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r163f)* c(sH) &
             + k(r164f)  )

    C4H8_num = ( 0.0_pr &
     +0.348535_pr * k(r1f)* c(sPOSF10325) &
     +0.318915_pr * k(r2f)* c(sH) * c(sPOSF10325) &
     +0.318915_pr * k(r3f)* c(sCH3) * c(sPOSF10325) &
     +0.318915_pr * k(r4f)* c(sOH) * c(sPOSF10325) &
     +0.318915_pr * k(r5f)* c(sO) * c(sPOSF10325)   )

    C4H8_ct = C4H8_num / C4H8_denom

    cqss(sqssC4H8) = C4H8_ct  

    HCCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r117f)&
             + k(r118f)* c(sH) &
             + k(r119f)* c(sO) &
             + k(r120f)* c(sO2) &
             + k(r121f)&
             + k(r122f)* c(sC2H2) &
             + k(r123f)* c(sCH3) &
             + k(r124f)* c(sOH)   )

    HCCO_num = ( 0.0_pr &
             + k(r76f)* cqss(sqssC2H) * c(sOH) &
             + k(r80f)* c(sC2H2) * c(sO) &
             + k(r138f)* c(sO) * c(sPXC3H4)   )

    HCCO_ct = HCCO_num / HCCO_denom

    cqss(sqssHCCO) = HCCO_ct  

    TXCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r23f)* c(sO) &
             + k(r24f)* c(sOH) &
             + k(r25f)* c(sOH) &
             + k(r26f)* c(sH2) &
             + k(r27f)* c(sO2) &
             + k(r28f)* c(sO2) &
             + k(r29f)* c(sO2) &
             + k(r47f)* c(sCH3) &
             + k(r54f)* c(sCH4) &
             + k(r6b)* c(sN2) &
             + k(r34b)* c(sH2O) &
             + k(r40b)* c(sH2O) &
             + k(r59b)* c(sCO) &
             + k(r73b)* c(sCO2)   )

    TXCH2_num = ( 0.0_pr &
             + k(r40f)* c(sCH3) * c(sOH) &
             + k(r81f)* c(sC2H2) * c(sO) &
             + k(r100f)* c(sC2H4) * c(sO) &
             + k(r26b)* c(sCH3) * c(sH) &
             + k(r54b) *c(sCH3)** 2.0_pr  )

    TXCH2_ct = TXCH2_num / TXCH2_denom

    TXCH2_SXCH2 = - ( 0.0_pr &
          + k(r6f) * c(sN2) &
          + k(r34f) * c(sH2O) &
          + k(r59f) * c(sCO) &
          + k(r73f) * c(sCO2)  ) / TXCH2_denom

    TXCH2_CH = - ( 0.0_pr &
          + k(r25b) * c(sH2O)  ) / TXCH2_denom

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r20f)* c(sO) &
             + k(r21f)* c(sH2O) &
             + k(r22f)* c(sO2) &
             + k(r46f)* c(sCH3) &
             + k(r53f)* c(sCH4) &
             + k(r72f)* c(sCO2) &
             + k(r25b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r75f)* cqss(sqssC2H) * c(sO) &
             + k(r117f)* cqss(sqssHCCO)   )

    CH_ct = CH_num / CH_denom

    CH_TXCH2 = - ( 0.0_pr &
          + k(r25f) * c(sOH)  ) / CH_denom

    SXCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r6f)* c(sN2) &
             + k(r30f)* c(sH2) &
             + k(r31f)* c(sO2) &
             + k(r32f)* c(sO2) &
             + k(r33f)* c(sH2O) &
             + k(r34f)* c(sH2O) &
             + k(r35f)* c(sH2O) &
             + k(r59f)* c(sCO) &
             + k(r73f)* c(sCO2) &
             + k(r74f)* c(sCO2) &
             + k(r86f)* c(sC2H2) &
             + k(r41b)* c(sH2O)   )

    SXCH2_num = ( 0.0_pr &
             + k(r41f)* c(sCH3) * c(sOH) &
             + k(r49f) *c(sCH3)** 2.0_pr&
             + k(r118f)* c(sH) * cqss(sqssHCCO) &
             + k(r30b)* c(sCH3) * c(sH) &
             + k(r86b)* c(sC3H3) * c(sH)   )

    SXCH2_ct = SXCH2_num / SXCH2_denom

    SXCH2_TXCH2 = - ( 0.0_pr &
          + k(r6b) * c(sN2) &
          + k(r34b) * c(sH2O) &
          + k(r59b) * c(sCO) &
          + k(r73b) * c(sCO2)  ) / SXCH2_denom

    c_e_a = CH_TXCH2
    c_h_a =  (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2)
    c_h_b =  (c_h_a) - (c_e_a) * (TXCH2_CH)
    C_A_A = CH_ct
    C_B_A =  (TXCH2_ct) - (SXCH2_ct) * (TXCH2_SXCH2)
    C_B_B =  (C_B_A) - (C_A_A) * (TXCH2_CH)

    cqss(sqssTXCH2) = ( C_B_B ) / ( c_h_b )  

    cqss(sqssCH) = C_A_A - (c_e_a) * cqss(sqssTXCH2)  

    cqss(sqssSXCH2) = SXCH2_ct - (SXCH2_TXCH2) * cqss(sqssTXCH2)  

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r87f)* c(sH) &
             + k(r88f)* c(sH) &
             + k(r89f)* c(sO) &
             + k(r90f)* c(sOH) &
             + k(r91f)* c(sO2) &
             + k(r92f)* c(sO2) &
             + k(r93f)* c(sO2) &
             + k(r94f)* cqss(sqssHCO) &
             + k(r95f)* c(sCH3) &
             + k(r96f)* c(sCH3) &
             + k(r79b)&
             + k(r154b)* c(sCH3)   )

    C2H3_num = ( 0.0_pr &
             + k(r46f)* cqss(sqssCH) * c(sCH3) &
             + k(r79f)* c(sC2H2) * c(sH) &
             + k(r98f)* c(sC2H4) * c(sH) &
             + k(r102f)* c(sC2H4) * c(sOH) &
             + k(r104f)* c(sC2H4) * c(sCH3) &
             + k(r130f)* c(sC3H3) * c(sHO2) &
             + k(r154f)* c(sC3H6) &
             + k(r167f)* c(sC5H5) * c(sO) &
             + k(r168f)* c(sC5H5) * c(sHO2) &
             + k(r170f)* c(sC5H5) * c(sOH) &
             + k(r96b)* c(sAXC3H5) * c(sH)   )

    C2H3_ct = C2H3_num / C2H3_denom

    cqss(sqssC2H3) = C2H3_ct  

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r60f)* c(sH) &
             + k(r61f)* c(sO) &
             + k(r62f)* c(sO) &
             + k(r63f)* c(sOH) &
             + k(r64f)* m(mM64) &
             + k(r65f)* c(sH2O) &
             + k(r66f)* c(sO2) &
             + k(r67f)* c(sCH3) &
             + k(r94f)* cqss(sqssC2H3) &
             + k(r109f)* c(sC2H5) &
             + k(r68b)* c(sH2)   )

    HCO_num = ( 0.0_pr &
             + k(r22f)* cqss(sqssCH) * c(sO2) &
             + k(r23f)* c(sO) * cqss(sqssTXCH2) &
             + k(r68f)* c(sCH2O) * c(sH) &
             + k(r69f)* c(sCH2O) * c(sO) &
             + k(r70f)* c(sCH2O) * c(sOH) &
             + k(r71f)* c(sCH2O) * c(sCH3) &
             + k(r72f)* cqss(sqssCH) * c(sCO2) &
             + k(r77f)* cqss(sqssC2H) * c(sO2) &
             + k(r93f)* cqss(sqssC2H3) * c(sO2) &
             + k(r101f)* c(sC2H4) * c(sO) &
     +2.0_pr * k(r124f)* cqss(sqssHCCO) * c(sOH) &
             + k(r158f)* c(sC3H6) * c(sO) &
             + k(r200f)* cqss(sqssA1CHO) &
             + k(r64b)* c(sCO) * c(sH) * m(mM64) &
             + k(r65b)* c(sCO) * c(sH) * c(sH2O)   )

    HCO_ct = HCO_num / HCO_denom

    cqss(sqssHCO) = HCO_ct  

    cqss = max(cqss, tiny(1.0_pr))
    cqss = min(cqss, 1e03_pr)
    
    ! Necessary for very low concentration variations
    ! Needs to be integrated only for Cantera
    do index=1,nqss
      if (cqss(index) /= cqss(index)) then
        cqss(index) = tiny(1.0_pr)
      end if
    end do
    
    return
  end subroutine get_qss  
  
    
  ! ----------------------------------------------- !
  ! Mass fractions to concentrations                !
  ! ----------------------------------------------- !
  subroutine y2c(y, W_sp, P, T, c)
    implicit none

    real(pr),dimension(nspec) :: W_sp
    real(pr),dimension(nspec) :: c, y
    real(pr) :: rho, P, T, inv_W_g

    integer :: k
    
    ! Gas molecular weight inverse
    inv_W_g = 0.0_pr
    do k =1, nspec
      inv_W_g = inv_W_g + y(k) / W_sp(k)
    end do

    ! Gas density
    rho = P / (Rcst * inv_W_g * T)

    ! Conversion
    c = y * rho / W_sp

    ! Concentrations clipping
    !do  k = 1, nspec
    !   if (c(k)<1.0e-20_pr) c(k) = 0.0_pr
    !end do

        return
    end subroutine y2c

end module mod_customkinetics
! ----------------------------------------------- !
! Cantera routine                                 !
! ----------------------------------------------- !
subroutine customkinetics(P, T, y, wdot)

  use mod_customkinetics
  implicit none

  real(pr), dimension(nspec) :: y, c, wdot, cdot
  real(pr), dimension(nqss) :: cqss
  real(pr), dimension(nreac + nreac_reverse) :: w,k
  real(pr), dimension(nTB + nFO) :: M
  real(pr) :: P, T, rho

  ! Convert to concentrations
  call y2c(y, W_sp, P, T, c)

  ! Evaluate QSS concentrations and reaction rates
  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)

  ! Evaluate production rates
  call get_production_rates(wdot,w)


  return
end subroutine customkinetics