!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE 
!>    Module for calculating the analytical source terms for 
!!    @file    
!!    @authors 
!!    @date    
!!    @since   
!!    @note    
!--------------------------------------------------------------------------------------------------

!--------------------------------------------------------------------------------------------------
!     MODULE mod_customkinetics
!>    Module for calculating the analytical source terms for 
!!    Generated by the YARC2 fortran generator 
!!    @authors 
!!    @date    
!--------------------------------------------------------------------------------------------------
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 non-qss and qss species and reactions
  integer, parameter :: nspec = 28
  integer, parameter :: nqss = 14
  integer, parameter :: nreac = 205
  integer, parameter :: nreac_reverse = 66
  integer, parameter :: isc_T = 1
  integer, parameter :: neq = nspec + 1
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 3
  integer, parameter :: nFO = 13
  integer, parameter :: nTB_reverse = 1
  integer, parameter :: nFO_reverse = 7
  
  ! Index of species
  integer, parameter :: sN2 = 1
  integer, parameter :: sO = 2
  integer, parameter :: sH2 = 3
  integer, parameter :: sH = 4
  integer, parameter :: sOH = 5
  integer, parameter :: sH2O = 6
  integer, parameter :: sO2 = 7
  integer, parameter :: sHO2 = 8
  integer, parameter :: sCO = 9
  integer, parameter :: sCH2O = 10
  integer, parameter :: sCH3 = 11
  integer, parameter :: sCO2 = 12
  integer, parameter :: sCH4 = 13
  integer, parameter :: sC2H4 = 14
  integer, parameter :: sC2H5 = 15
  integer, parameter :: sC2H2 = 16
  integer, parameter :: sC3H3 = 17
  integer, parameter :: sAC3H5 = 18
  integer, parameter :: sC2H6 = 19
  integer, parameter :: sPC3H4 = 20
  integer, parameter :: sAC3H4 = 21
  integer, parameter :: sC6H6 = 22
  integer, parameter :: sC5H5 = 23
  integer, parameter :: sC3H6 = 24
  integer, parameter :: sC5H6 = 25
  integer, parameter :: sC10H8 = 26
  integer, parameter :: sC7H8 = 27
  integer, parameter :: sC8H6 = 28

  integer, parameter :: sqssSCH2 = 1
  integer, parameter :: sqssTCH2 = 2
  integer, parameter :: sqssCH = 3
  integer, parameter :: sqssHCO = 4
  integer, parameter :: sqssC2H3 = 5
  integer, parameter :: sqssC2H = 6
  integer, parameter :: sqssHCCO = 7
  integer, parameter :: sqssNC3H7 = 8
  integer, parameter :: sqssC6H5 = 9
  integer, parameter :: sqssC8H7 = 10
  integer, parameter :: sqssC7H7 = 11
  integer, parameter :: sqssC7H6O = 12
  integer, parameter :: sqssC8H5 = 13
  integer, parameter :: sqssC10H7 = 14

  ! 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 :: r1b = 206
  integer, parameter :: r2b = 207
  integer, parameter :: r5b = 208
  integer, parameter :: r6b = 209
  integer, parameter :: r8b = 210
  integer, parameter :: r17b = 211
  integer, parameter :: r22b = 212
  integer, parameter :: r23b = 213
  integer, parameter :: r27b = 214
  integer, parameter :: r31b = 215
  integer, parameter :: r33b = 216
  integer, parameter :: r37b = 217
  integer, parameter :: r38b = 218
  integer, parameter :: r45b = 219
  integer, parameter :: r47b = 220
  integer, parameter :: r48b = 221
  integer, parameter :: r49b = 222
  integer, parameter :: r51b = 223
  integer, parameter :: r53b = 224
  integer, parameter :: r54b = 225
  integer, parameter :: r56b = 226
  integer, parameter :: r61b = 227
  integer, parameter :: r62b = 228
  integer, parameter :: r65b = 229
  integer, parameter :: r70b = 230
  integer, parameter :: r75b = 231
  integer, parameter :: r76b = 232
  integer, parameter :: r83b = 233
  integer, parameter :: r92b = 234
  integer, parameter :: r93b = 235
  integer, parameter :: r101b = 236
  integer, parameter :: r102b = 237
  integer, parameter :: r108b = 238
  integer, parameter :: r120b = 239
  integer, parameter :: r121b = 240
  integer, parameter :: r126b = 241
  integer, parameter :: r127b = 242
  integer, parameter :: r128b = 243
  integer, parameter :: r129b = 244
  integer, parameter :: r130b = 245
  integer, parameter :: r136b = 246
  integer, parameter :: r137b = 247
  integer, parameter :: r138b = 248
  integer, parameter :: r142b = 249
  integer, parameter :: r148b = 250
  integer, parameter :: r149b = 251
  integer, parameter :: r150b = 252
  integer, parameter :: r151b = 253
  integer, parameter :: r164b = 254
  integer, parameter :: r175b = 255
  integer, parameter :: r176b = 256
  integer, parameter :: r182b = 257
  integer, parameter :: r183b = 258
  integer, parameter :: r187b = 259
  integer, parameter :: r188b = 260
  integer, parameter :: r189b = 261
  integer, parameter :: r196b = 262
  integer, parameter :: r197b = 263
  integer, parameter :: r198b = 264
  integer, parameter :: r199b = 265
  integer, parameter :: r200b = 266
  integer, parameter :: r201b = 267
  integer, parameter :: r202b = 268
  integer, parameter :: r203b = 269
  integer, parameter :: r204b = 270
  integer, parameter :: r205b = 271
  
  ! Index of third body species
  integer, parameter :: mM3 = 1
  integer, parameter :: mM4 = 2
  integer, parameter :: mM61 = 3
  
  integer, parameter :: mM7 = 4
  integer, parameter :: mM30 = 5
  integer, parameter :: mM33 = 6
  integer, parameter :: mM36 = 7
  integer, parameter :: mM52 = 8
  integer, parameter :: mM76 = 9
  integer, parameter :: mM84 = 10
  integer, parameter :: mM93 = 11
  integer, parameter :: mM101 = 12
  integer, parameter :: mM102 = 13
  integer, parameter :: mM108 = 14
  integer, parameter :: mM113 = 15
  integer, parameter :: mM148 = 16
  
  ! Index of third body reactions
  integer, parameter :: TBr3f = 1
  integer, parameter :: TBr4f = 2
  integer, parameter :: TBr61f = 3
  integer, parameter :: TBr61b = 4
  
  integer, parameter :: FOr7f = 1
  integer, parameter :: FOr30f = 2
  integer, parameter :: FOr33f = 3
  integer, parameter :: FOr36f = 4
  integer, parameter :: FOr52f = 5
  integer, parameter :: FOr76f = 6
  integer, parameter :: FOr84f = 7
  integer, parameter :: FOr93f = 8
  integer, parameter :: FOr101f = 9
  integer, parameter :: FOr102f = 10
  integer, parameter :: FOr108f = 11
  integer, parameter :: FOr113f = 12
  integer, parameter :: FOr148f = 13
  integer, parameter :: FOr33b = 14
  integer, parameter :: FOr76b = 15
  integer, parameter :: FOr93b = 16
  integer, parameter :: FOr101b = 17
  integer, parameter :: FOr102b = 18
  integer, parameter :: FOr108b = 19
  integer, parameter :: FOr148b = 20
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.02801348_pr, & ! N2
       0.0159994_pr, & ! O
       0.0020158800000000003_pr, & ! H2
       0.0010079400000000001_pr, & ! H
       0.01700734_pr, & ! OH
       0.01801528_pr, & ! H2O
       0.0319988_pr, & ! O2
       0.03300674_pr, & ! HO2
       0.028010399999999998_pr, & ! CO
       0.03002628_pr, & ! CH2O
       0.015034819999999999_pr, & ! CH3
       0.0440098_pr, & ! CO2
       0.01604276_pr, & ! CH4
       0.028053759999999997_pr, & ! C2H4
       0.0290617_pr, & ! C2H5
       0.02603788_pr, & ! C2H2
       0.03905682_pr, & ! C3H3
       0.0410727_pr, & ! AC3H5
       0.030069639999999998_pr, & ! C2H6
       0.04006476_pr, & ! PC3H4
       0.04006476_pr, & ! AC3H4
       0.07811364_pr, & ! C6H6
       0.06509469999999999_pr, & ! C5H5
       0.04208064_pr, & ! C3H6
       0.06610263999999999_pr, & ! C5H6
       0.12817351999999999_pr, & ! C10H8
       0.09214051999999999_pr, & ! C7H8
       0.10213564_pr & ! C8H6
      !0.014026879999999999_pr, & ! SCH2
      !0.014026879999999999_pr, & ! TCH2
      !0.01301894_pr, & ! CH
      !0.02901834_pr, & ! HCO
      !0.027045819999999998_pr, & ! C2H3
      !0.02502994_pr, & ! C2H
      !0.04102934_pr, & ! HCCO
      !0.04308858_pr, & ! NC3H7
      !0.0771057_pr, & ! C6H5
      !0.10314358_pr, & ! C8H7
      !0.09113258_pr, & ! C7H7
      !0.10612403999999999_pr, & ! C7H6O
      !0.10112769999999999_pr, & ! C8H5
      !0.12716558_pr & ! C10H7
       /)
  
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 = k0 * conc / max(kinf, 1.0e-60_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(mM3) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM4) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (5.3_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM61) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (-1.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM7) = (-0.25_pr)*c(sH2) &
         + (10.89_pr)*c(sH2O) &
         + (1.18_pr)*c(sCO2) &
         + (0.09_pr)*c(sCO) &
         + (-0.15_pr)*c(sO2) &
         + sum(c)

    M(mM30) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM33) = (1.0_pr)*c(sCO2) &
         + (0.5_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (5.0_pr)*c(sH2O) &
         + (2.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM36) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM52) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM76) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM84) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM93) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM101) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM102) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM108) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM113) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + sum(c)

    M(mM148) = (2.6_pr)*c(sCO2) &
         + (0.75_pr)*c(sCO) &
         + (1.0_pr)*c(sH2) &
         + (11.0_pr)*c(sH2O) &
         + (1.0_pr)*c(sCH4) &
         + (2.0_pr)*c(sC2H6) &
         + 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) = (15000000.0_pr)*exp((-2509.9816_pr)*R_T_inv)
    k(r2f) = (0.045899999999999996_pr)*exp((-26189.999040000002_pr)*R_T_inv + T_log *(2.7_pr))
    k(r3f) = (9430000.0_pr)*exp(T_log *(-1.0_pr))
    k(r4f) = (44000000000.0_pr)*exp(T_log *(-2.0_pr))
    k(r5f) = (173.0_pr)*exp((-14349.99032_pr)*R_T_inv + T_log *(1.51_pr))
    k(r6f) = (0.039700000000000006_pr)*exp((8829.997280000001_pr)*R_T_inv + T_log *(2.4_pr))
    k_0(FOr7f) = (63300000.0_pr)*exp(T_log *(-1.4_pr))
    k_inf(FOr7f) = (5120000.0_pr)*exp(T_log *(0.44_pr))
    FC(FOr7f) = ((1.0_pr - 0.5_pr)*exp(-Tloc/(0.0_pr))) + (0.5_pr*exp(-Tloc/(10000000000.0_pr))) + exp(-(10000000000.0)/Tloc)
    k(r7f) = getlindratecoeff(Tloc,k_0(FOr7f),k_inf(FOr7f),FC(FOr7f), M(mM7),Ploc)
    k(r8f) = (26400000000.0_pr)*exp((-71300.00424000001_pr)*R_T_inv + T_log *(-0.67_pr))
    k(r9f) = (3.649_pr)*exp((4570.01584_pr)*R_T_inv + T_log *(2.07_pr))
    k(r10f) = (3970000.0_pr)*exp((-2810.0162400000004_pr)*R_T_inv)
    k(r11f) = (74900000.0_pr)*exp((-2660.01984_pr)*R_T_inv)
    k(r12f) = (40000000.0_pr)
    k(r13f) = (23800000.0_pr)*exp((2089.99168_pr)*R_T_inv)
    k(r14f) = (10000000000.0_pr)*exp((-72510.01704_pr)*R_T_inv)
    k(r15f) = (57000000.0_pr)
    k(r16f) = (30000000.0_pr)
    k(r17f) = (108000000.0_pr)*exp((-13009.98064_pr)*R_T_inv)
    k(r18f) = (5710000.0_pr)*exp((3160.0078399999998_pr)*R_T_inv)
    k(r19f) = (67100000.0_pr)
    k(r20f) = (80000000.0_pr)
    k(r21f) = (20000000.0_pr)
    k(r22f) = (11.3_pr)*exp((-12549.99168_pr)*R_T_inv + T_log *(2.0_pr))
    k(r23f) = (0.5_pr)*exp((-30249.98528_pr)*R_T_inv + T_log *(2.0_pr))
    k(r24f) = (5800000.0_pr)*exp((-6280.016640000001_pr)*R_T_inv)
    k(r25f) = (2400000.0_pr)*exp((-6280.016640000001_pr)*R_T_inv)
    k(r26f) = (5000000.0_pr)*exp((-6280.016640000001_pr)*R_T_inv)
    k(r27f) = (70000000.0_pr)
    k(r28f) = (28000000.0_pr)
    k(r29f) = (12000000.0_pr)
    k_0(FOr30f) = (1.88e+26_pr)*exp((-21089.99592_pr)*R_T_inv + T_log *(-6.36_pr))
    k_inf(FOr30f) = (4.82e+11_pr)*exp((-4790.01056_pr)*R_T_inv + T_log *(-1.16_pr))
    FC(FOr30f) = ((1.0_pr - 0.6027_pr)*exp(-Tloc/(208.0_pr))) &
               + (0.6027_pr*exp(-Tloc/(3921.9999999999995_pr))) + exp(-(10180.0)/Tloc)
    k(r30f) = getlindratecoeff(Tloc,k_0(FOr30f),k_inf(FOr30f),FC(FOr30f), M(mM30),Ploc)
    k(r31f) = (30000000.0_pr)
    k(r32f) = (68200.0_pr)*exp((3909.9898399999997_pr)*R_T_inv + T_log *(0.25_pr))
    k_0(FOr33f) = (3.47e+26_pr)*exp((-21229.99256_pr)*R_T_inv + T_log *(-6.3_pr))
    k_inf(FOr33f) = (69200000.0_pr)*exp(T_log *(0.18_pr))
    k(r33f) = getlindratecoeff(Tloc,k_0(FOr33f),k_inf(FOr33f),1.0_pr, M(mM33),Ploc)
    k(r34f) = (50600000.0_pr)
    k(r35f) = (33700000.0_pr)
    k_0(FOr36f) = (4e+24_pr)*exp((-13140.01936_pr)*R_T_inv + T_log *(-5.92_pr))
    k_inf(FOr36f) = (2.79e+12_pr)*exp((-5569.99184_pr)*R_T_inv + T_log *(-1.43_pr))
    FC(FOr36f) = ((1.0_pr - 0.412_pr)*exp(-Tloc/(195.0_pr))) + (0.412_pr*exp(-Tloc/(5900.0_pr))) + exp(-(6394.0)/Tloc)
    k(r36f) = getlindratecoeff(Tloc,k_0(FOr36f),k_inf(FOr36f),FC(FOr36f), M(mM36),Ploc)
    k(r37f) = (56.0_pr)*exp((-22679.9996_pr)*R_T_inv + T_log *(1.6_pr))
    k(r38f) = (644000000000.0_pr)*exp((-5929.983200000001_pr)*R_T_inv + T_log *(-1.34_pr))
    k(r39f) = (13800000.0_pr)*exp((-127529.99360000002_pr)*R_T_inv)
    k(r40f) = (587000.0_pr)*exp((-57909.99088_pr)*R_T_inv)
    k(r41f) = (10000000.0_pr)
    k(r42f) = (3610000.0_pr)
    k(r43f) = (30000000.0_pr)
    k(r44f) = (100000000.0_pr)
    k(r45f) = (6840000.0_pr)*exp((-44349.9816_pr)*R_T_inv + T_log *(0.1_pr))
    k(r46f) = (2632000.0_pr)*exp((-57160.008879999994_pr)*R_T_inv + T_log *(-0.06_pr))
    k(r47f) = (660.0_pr)*exp((-45359.9992_pr)*R_T_inv + T_log *(1.62_pr))
    k(r48f) = (1020.0_pr)*exp((-35980.015120000004_pr)*R_T_inv + T_log *(1.5_pr))
    k(r49f) = (100.0_pr)*exp((-13049.97968_pr)*R_T_inv + T_log *(1.6_pr))
    k(r50f) = (60000000.0_pr)
    k(r51f) = (2.46_pr)*exp((-34600.0064_pr)*R_T_inv + T_log *(2.0_pr))
    k_0(FOr52f) = (1.17e+12_pr)*exp((-17539.99744_pr)*R_T_inv + T_log *(-2.79_pr))
    k_inf(FOr52f) = (13600.0_pr)*exp((-9980.01152_pr)*R_T_inv)
    FC(FOr52f) = ((1.0_pr - 1.0_pr)*exp(-Tloc/(1.0_pr))) + (1.0_pr*exp(-Tloc/(10000000.0_pr))) + exp(-(10000000.0)/Tloc)
    k(r52f) = getlindratecoeff(Tloc,k_0(FOr52f),k_inf(FOr52f),FC(FOr52f), M(mM52),Ploc)
    k(r53f) = (800000.0_pr)*exp((-30760.01488_pr)*R_T_inv + T_log *(0.14_pr))
    k(r54f) = (87800.0_pr)*exp((69.99832_pr)*R_T_inv + T_log *(0.03_pr))
    k(r55f) = (30100000.0_pr)*exp((-96229.99168_pr)*R_T_inv)
    k(r56f) = (9000000.0_pr)
    k(r57f) = (120000000.0_pr)
    k(r58f) = (30000000.0_pr)
    k(r59f) = (30000000.0_pr)
    k(r60f) = (30200000.0_pr)
    k(r61f) = (187000000000.0_pr)*exp((-71130.00832_pr)*R_T_inv + T_log *(-1.0_pr))
    k(r62f) = (2240000000000.0_pr)*exp((-71130.00832_pr)*R_T_inv + T_log *(-1.0_pr))
    k(r63f) = (12000.0_pr)*exp((3040.01072_pr)*R_T_inv + T_log *(0.81_pr))
    k(r64f) = (26500000.0_pr)
    k(r65f) = (57.4_pr)*exp((-11470.0176_pr)*R_T_inv + T_log *(1.9_pr))
    k(r66f) = (39000000.0_pr)*exp((-14809.979280000001_pr)*R_T_inv)
    k(r67f) = (3430.0_pr)*exp((1869.99696_pr)*R_T_inv + T_log *(1.18_pr))
    k(r68f) = (0.00332_pr)*exp((-24519.99728_pr)*R_T_inv + T_log *(2.81_pr))
    k(r69f) = (190000000.0_pr)*exp((-66070.00424000001_pr)*R_T_inv)
    k(r70f) = (7000000.0_pr)
    k(r71f) = (14000000.0_pr)
    k(r72f) = (50000000.0_pr)
    k(r73f) = (20000000.0_pr)
    k(r74f) = (10000000.0_pr)*exp((3160.0078399999998_pr)*R_T_inv)
    k(r75f) = (3.31_pr)*exp((-3769.9932_pr)*R_T_inv + T_log *(2.26_pr))
    k_0(FOr76f) = (6.34e+19_pr)*exp((-15819.99688_pr)*R_T_inv + T_log *(-4.66_pr))
    k_inf(FOr76f) = (17100.0_pr)*exp((-11329.97912_pr)*R_T_inv + T_log *(1.27_pr))
    FC(FOr76f) = ((1.0_pr - 0.2122_pr)*exp(-Tloc/(1.0_pr))) + (0.2122_pr*exp(-Tloc/(-10210.0_pr)))
    k(r76f) = getlindratecoeff(Tloc,k_0(FOr76f),k_inf(FOr76f),FC(FOr76f), M(mM76),Ploc)
    k(r77f) = (8.1_pr)*exp((-7950.018399999999_pr)*R_T_inv + T_log *(2.0_pr))
    k(r78f) = (12.5_pr)*exp((-7950.018399999999_pr)*R_T_inv + T_log *(2.0_pr))
    k(r79f) = (33240000000.0_pr)*exp((-128440.0136_pr)*R_T_inv + T_log *(-0.44_pr))
    k(r80f) = (2.63_pr)*exp((-71380.00232_pr)*R_T_inv + T_log *(2.14_pr))
    k(r81f) = (7.53_pr)*exp((-8809.99776_pr)*R_T_inv + T_log *(1.55_pr))
    k(r82f) = (1280.0_pr)*exp((-10789.99208_pr)*R_T_inv + T_log *(0.73_pr))
    k(r83f) = (190000000.0_pr)
    k_0(FOr84f) = (1.4e+18_pr)*exp((-13890.00136_pr)*R_T_inv + T_log *(-3.86_pr))
    k_inf(FOr84f) = (6080000.0_pr)*exp((-1170.01376_pr)*R_T_inv + T_log *(0.27_pr))
    FC(FOr84f) = ((1.0_pr - 0.782_pr)*exp(-Tloc/(207.49999999999997_pr))) + (0.782_pr*exp(-Tloc/(2663.0_pr))) + exp(-(6095.0)/Tloc)
    k(r84f) = getlindratecoeff(Tloc,k_0(FOr84f),k_inf(FOr84f),FC(FOr84f), M(mM84),Ploc)
    k(r85f) = (30000000.0_pr)
    k(r86f) = (10300000.0_pr)*exp((1789.9988799999999_pr)*R_T_inv + T_log *(0.21_pr))
    k(r87f) = (5000000.0_pr)
    k(r88f) = (1.34_pr)*exp((1610.0031999999999_pr)*R_T_inv + T_log *(1.61_pr))
    k(r89f) = (303000.0_pr)*exp((-49.998799999999996_pr)*R_T_inv + T_log *(0.29_pr))
    k(r90f) = (45800000000.0_pr)*exp((-4249.98168_pr)*R_T_inv + T_log *(-1.39_pr))
    k(r91f) = (9030000.0_pr)*exp((3200.0068800000004_pr)*R_T_inv)
    k(r92f) = (1930000000000.0_pr)*exp((-32089.982959999998_pr)*R_T_inv + T_log *(-1.25_pr))
    k_0(FOr93f) = (2.03e+27_pr)*exp((-24140.0064_pr)*R_T_inv + T_log *(-6.64_pr))
    k_inf(FOr93f) = (1370.0_pr)*exp((-5669.98944_pr)*R_T_inv + T_log *(1.46_pr))
    FC(FOr93f) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299.0_pr))) + (-0.569_pr*exp(-Tloc/(-9147.0_pr))) + exp(-(152.4)/Tloc)
    k(r93f) = getlindratecoeff(Tloc,k_0(FOr93f),k_inf(FOr93f),FC(FOr93f), M(mM93),Ploc)
    k(r94f) = (1.33_pr)*exp((-51209.98432_pr)*R_T_inv + T_log *(2.53_pr))
    k(r95f) = (7660.0_pr)*exp((-4770.01104_pr)*R_T_inv + T_log *(0.88_pr))
    k(r96f) = (0.0715_pr)*exp((-3889.9903200000003_pr)*R_T_inv + T_log *(2.47_pr))
    k(r97f) = (389.0_pr)*exp((-3709.9946400000003_pr)*R_T_inv + T_log *(1.36_pr))
    k(r98f) = (1.3100000000000002e-07_pr)*exp((3599.9972799999996_pr)*R_T_inv + T_log *(4.2_pr))
    k(r99f) = (3.7500000000000004e+30_pr)*exp((-29540.002319999996_pr)*R_T_inv + T_log *(-7.8_pr))
    k(r100f) = (0.227_pr)*exp((-38489.996719999996_pr)*R_T_inv + T_log *(2.0_pr))
    k_0(FOr101f) = (3e+51_pr)*exp((-76020.01648_pr)*R_T_inv + T_log *(-14.6_pr))
    k_inf(FOr101f) = (2.55_pr)*exp((-23850.01336_pr)*R_T_inv + T_log *(1.6_pr))
    FC(FOr101f) = ((1.0_pr - 0.1894_pr)*exp(-Tloc/(277.0_pr))) + (0.1894_pr*exp(-Tloc/(8748.0_pr))) + exp(-(7891.0)/Tloc)
    k(r101f) = getlindratecoeff(Tloc,k_0(FOr101f),k_inf(FOr101f),FC(FOr101f), M(mM101),Ploc)
    k_0(FOr102f) = (1.99e+29_pr)*exp((-27969.99816_pr)*R_T_inv + T_log *(-7.08_pr))
    k_inf(FOr102f) = (5.21e+11_pr)*exp((-6610.00872_pr)*R_T_inv + T_log *(-0.99_pr))
    FC(FOr102f) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125.0_pr))) + (0.8422_pr*exp(-Tloc/(2219.0_pr))) + exp(-(6882.0)/Tloc)
    k(r102f) = getlindratecoeff(Tloc,k_0(FOr102f),k_inf(FOr102f),FC(FOr102f), M(mM102),Ploc)
    k(r103f) = (31700000.0_pr)*exp((1650.00224_pr)*R_T_inv + T_log *(0.03_pr))
    k(r104f) = (19.2_pr)*exp((8510.00496_pr)*R_T_inv + T_log *(1.02_pr))
    k(r105f) = (120000000.0_pr)
    k(r106f) = (300000.0_pr)
    k(r107f) = (31000000.0_pr)
    k_0(FOr108f) = (3.72e+59_pr)*exp((-425010.00872_pr)*R_T_inv + T_log *(-13.14_pr))
    k_inf(FOr108f) = (1.88e+50_pr)*exp((-449120.01584_pr)*R_T_inv + T_log *(-9.72_pr))
    FC(FOr108f) = ((1.0_pr - 0.39_pr)*exp(-Tloc/(100.0_pr))) + (0.39_pr*exp(-Tloc/(1900.0_pr))) + exp(-(6000.0)/Tloc)
    k(r108f) = getlindratecoeff(Tloc,k_0(FOr108f),k_inf(FOr108f),FC(FOr108f), M(mM108),Ploc)
    k(r109f) = (0.17_pr)*exp((-24020.009280000002_pr)*R_T_inv + T_log *(2.7_pr))
    k(r110f) = (89.8_pr)*exp((-23810.014319999995_pr)*R_T_inv + T_log *(1.92_pr))
    k(r111f) = (1.61_pr)*exp((-3100.0092799999998_pr)*R_T_inv + T_log *(2.22_pr))
    k(r112f) = (843000000.0_pr)*exp((-93119.98264_pr)*R_T_inv)
    k_0(FOr113f) = (1.213e+29_pr)*exp((-320779.99984_pr)*R_T_inv + T_log *(-5.18_pr))
    k_inf(FOr113f) = (2.255e+20_pr)*exp((-312679.98504_pr)*R_T_inv + T_log *(-1.44_pr))
    FC(FOr113f) = ((1.0_pr - 0.5757_pr)*exp(-Tloc/(237.00000000000003_pr))) &
                + (0.5757_pr*exp(-Tloc/(1652.0_pr))) + exp(-(5069.0)/Tloc)
    k(r113f) = getlindratecoeff(Tloc,k_0(FOr113f),k_inf(FOr113f),FC(FOr113f), M(mM113),Ploc)
    k(r114f) = (100000000.0_pr)
    k(r115f) = (100000000.0_pr)
    k(r116f) = (42000.0_pr)*exp((-3569.998_pr)*R_T_inv)
    k(r117f) = (100000.0_pr)*exp((-12549.99168_pr)*R_T_inv)
    k(r118f) = (50000000.0_pr)
    k(r119f) = (10000000.0_pr)
    k(r120f) = (7.939999999999999e+23_pr)*exp((-20340.01392_pr)*R_T_inv + T_log *(-5.06_pr))
    k(r121f) = (3.16e+23_pr)*exp((-19709.9872_pr)*R_T_inv + T_log *(-5.0_pr))
    k(r122f) = (1280.0_pr)*exp((-10789.99208_pr)*R_T_inv + T_log *(0.73_pr))
    k(r123f) = (69500000.0_pr)
    k(r124f) = (0.17_pr)*exp((-6280.016640000001_pr)*R_T_inv + T_log *(1.7_pr))
    k(r125f) = (800000.0_pr)
    k(r126f) = (1.8699999999999998e+40_pr)*exp((-70309.98616000001_pr)*R_T_inv + T_log *(-9.84_pr))
    k(r127f) = (5.770000000000001e+31_pr)*exp((-131820.01616_pr)*R_T_inv + T_log *(-7.0_pr))
    k(r128f) = (235000.0_pr)*exp((-41820.000479999995_pr)*R_T_inv)
    k(r129f) = (3460000.0_pr)*exp((-22859.995280000003_pr)*R_T_inv + T_log *(0.44_pr))
    k(r130f) = (0.085_pr)*exp((-24020.009280000002_pr)*R_T_inv + T_log *(2.7_pr))
    k(r131f) = (0.805_pr)*exp((-3100.0092799999998_pr)*R_T_inv + T_log *(2.22_pr))
    k(r132f) = (422000000.0_pr)*exp((-93119.98264_pr)*R_T_inv)
    k(r133f) = (4.05_pr)*exp((-7950.018399999999_pr)*R_T_inv + T_log *(2.0_pr))
    k(r134f) = (6.25_pr)*exp((-7950.018399999999_pr)*R_T_inv + T_log *(2.0_pr))
    k(r135f) = (1280.0_pr)*exp((-10789.99208_pr)*R_T_inv + T_log *(0.73_pr))
    k(r136f) = (89500000.0_pr)*exp((-47070.0_pr)*R_T_inv + T_log *(-0.02_pr))
    k(r137f) = (7.76e+39_pr)*exp((-328219.98864000005_pr)*R_T_inv + T_log *(-7.8_pr))
    k(r138f) = (2470000000.0_pr)*exp((-26929.98128_pr)*R_T_inv + T_log *(-0.33_pr))
    k(r139f) = (1.33_pr)*exp((-51209.98432_pr)*R_T_inv + T_log *(2.53_pr))
    k(r140f) = (1.3100000000000002e-07_pr)*exp((3599.9972799999996_pr)*R_T_inv + T_log *(4.2_pr))
    k(r141f) = (0.227_pr)*exp((-38489.996719999996_pr)*R_T_inv + T_log *(2.0_pr))
    k(r142f) = (2.01e+43_pr)*exp((-82099.99608_pr)*R_T_inv + T_log *(-10.77_pr))
    k(r143f) = (0.00956_pr)*exp((-13770.00424_pr)*R_T_inv + T_log *(2.8_pr))
    k(r144f) = (6030000.0_pr)
    k(r145f) = (486000.0_pr)*exp((549.9867999999999_pr)*R_T_inv + T_log *(-0.32_pr))
    k(r146f) = (971000000000000.0_pr)*exp((-104520.00192_pr)*R_T_inv + T_log *(-2.7_pr))
    k(r147f) = (2.16e+33_pr)*exp((-99799.98968000001_pr)*R_T_inv + T_log *(-7.74_pr))
    k_0(FOr148f) = (6.26e+26_pr)*exp((-29290.00832_pr)*R_T_inv + T_log *(-6.66_pr))
    k_inf(FOr148f) = (306000000.0_pr)*exp((-16870.01352_pr)*R_T_inv + T_log *(-0.37_pr))
    FC(FOr148f) = ((1.0_pr - 1.0_pr)*exp(-Tloc/(1000.0_pr))) + (1.0_pr*exp(-Tloc/(1310.0_pr))) + exp(-(48100.0)/Tloc)
    k(r148f) = getlindratecoeff(Tloc,k_0(FOr148f),k_inf(FOr148f),FC(FOr148f), M(mM148),Ploc)
    k(r149f) = (4.04e+42_pr)*exp((-467899.98352_pr)*R_T_inv + T_log *(-7.67_pr))
    k(r150f) = (5.93e+48_pr)*exp((-98530.02016000001_pr)*R_T_inv + T_log *(-11.76_pr))
    k(r151f) = (8000000000000000.0_pr)*exp((-46780.00696_pr)*R_T_inv + T_log *(-2.39_pr))
    k(r152f) = (120.0_pr)*exp((-1370.00896_pr)*R_T_inv + T_log *(1.6_pr))
    k(r153f) = (35.0_pr)*exp((4069.986_pr)*R_T_inv + T_log *(1.6_pr))
    k(r154f) = (0.66_pr)*exp((-28269.990959999996_pr)*R_T_inv + T_log *(2.54_pr))
    k(r155f) = (0.0965_pr)*exp((-15550.003359999999_pr)*R_T_inv + T_log *(2.68_pr))
    k(r156f) = (200.0_pr)*exp((-2249.98784_pr)*R_T_inv + T_log *(1.46_pr))
    k(r157f) = (4.5199999999999997e-07_pr)*exp((-29929.992959999996_pr)*R_T_inv + T_log *(3.65_pr))
    k(r158f) = (2561000.0_pr)*exp((-13039.97992_pr)*R_T_inv + T_log *(0.06_pr))
    k(r159f) = (6.39e+23_pr)*exp((-147299.97936000003_pr)*R_T_inv + T_log *(-4.03_pr))
    k(r160f) = (70000000.0_pr)
    k(r161f) = (6.776e+23_pr)*exp((-48780.0008_pr)*R_T_inv + T_log *(-4.7_pr))
    k(r162f) = (1.19e+27_pr)*exp((-56069.9932_pr)*R_T_inv + T_log *(-6.52_pr))
    k(r163f) = (30200000.0_pr)
    k(r164f) = (1.73e+68_pr)*exp((-486899.98776_pr)*R_T_inv + T_log *(-15.16_pr))
    k(r165f) = (28000000.0_pr)*exp((-9449.9824_pr)*R_T_inv)
    k(r166f) = (660000000.0_pr)*exp((-51650.0156_pr)*R_T_inv)
    k(r167f) = (0.047700000000000006_pr)*exp((-4630.014399999999_pr)*R_T_inv + T_log *(2.71_pr))
    k(r168f) = (3.08_pr)*exp(T_log *(2.0_pr))
    k(r169f) = (1.8000000000000002e-07_pr)*exp(T_log *(4.0_pr))
    k(r170f) = (26000000.0_pr)*exp((-25610.012959999996_pr)*R_T_inv)
    k(r171f) = (30000000.0_pr)*exp((-37570.018800000005_pr)*R_T_inv)
    k(r172f) = (100000000.0_pr)
    k(r173f) = (30000000.0_pr)
    k(r174f) = (30000000.0_pr)
    k(r175f) = (3.29_pr)*exp((-13230.0172_pr)*R_T_inv + T_log *(2.05_pr))
    k(r176f) = (1.29e+61_pr)*exp((-619589.98304_pr)*R_T_inv + T_log *(-12.48_pr))
    k(r177f) = (602.0_pr)*exp((-68419.98968000001_pr)*R_T_inv + T_log *(1.8_pr))
    k(r178f) = (0.00040300000000000004_pr)*exp((-6089.979359999999_pr)*R_T_inv + T_log *(3.33_pr))
    k(r179f) = (2.752e-08_pr)*exp((-57060.01128_pr)*R_T_inv + T_log *(4.46_pr))
    k(r180f) = (22200000.0_pr)*exp((-18960.0052_pr)*R_T_inv)
    k(r181f) = (0.000132_pr)*exp((-23389.982560000004_pr)*R_T_inv + T_log *(3.25_pr))
    k(r182f) = (820000000000000.0_pr)*exp((-337550.01576_pr)*R_T_inv)
    k(r183f) = (5.83e+61_pr)*exp((-285890.0004_pr)*R_T_inv + T_log *(-14.15_pr))
    k(r184f) = (331000000.0_pr)
    k(r185f) = (10600000000.0_pr)*exp((-10559.997599999999_pr)*R_T_inv + T_log *(-0.94_pr))
    k(r186f) = (4.32e+33_pr)*exp((-99799.98968000001_pr)*R_T_inv + T_log *(-7.74_pr))
    k(r187f) = (2.31_pr)*exp((-17420.000319999996_pr)*R_T_inv + T_log *(2.17_pr))
    k(r188f) = (1.25e+18_pr)*exp((-396589.97952_pr)*R_T_inv + T_log *(-0.6_pr))
    k(r189f) = (2.16e+29_pr)*exp((-460929.98344_pr)*R_T_inv + T_log *(-3.58_pr))
    k(r190f) = (6.47e-06_pr)*exp((-14159.99488_pr)*R_T_inv + T_log *(3.98_pr))
    k(r191f) = (16200000.0_pr)*exp((-11590.01472_pr)*R_T_inv)
    k(r192f) = (422000000.0_pr)*exp((-93119.98264_pr)*R_T_inv)
    k(r193f) = (2610000000000000.0_pr)*exp((-337019.98664_pr)*R_T_inv + T_log *(0.15_pr))
    k(r194f) = (2050.0_pr)*exp((-10060.0096_pr)*R_T_inv + T_log *(1.16_pr))
    k(r195f) = (2.72_pr)*exp((-24769.991280000002_pr)*R_T_inv + T_log *(1.77_pr))
    k(r196f) = (0.0134_pr)*exp((-5369.99664_pr)*R_T_inv + T_log *(2.5_pr))
    k(r197f) = (3.6200000000000002e+22_pr)*exp((-99849.98848_pr)*R_T_inv + T_log *(-4.24_pr))
    k(r198f) = (2.1e+60_pr)*exp((-619549.984_pr)*R_T_inv + T_log *(-12.4_pr))
    k(r199f) = (132.0_pr)*exp((-70379.98448_pr)*R_T_inv + T_log *(1.88_pr))
    k(r200f) = (0.000134_pr)*exp((-6089.979359999999_pr)*R_T_inv + T_log *(3.33_pr))
    k(r201f) = (1.34e+17_pr)*exp((-172540.00120000003_pr)*R_T_inv + T_log *(-0.86_pr))
    k(r202f) = (70910000.0_pr)*exp((-29300.00808_pr)*R_T_inv + T_log *(-0.26_pr))
    k(r203f) = (8.6e+60_pr)*exp((-619549.984_pr)*R_T_inv + T_log *(-12.48_pr))
    k(r204f) = (265.0_pr)*exp((-71529.99872000002_pr)*R_T_inv + T_log *(1.87_pr))
    k(r205f) = (0.000963_pr)*exp((-18299.979199999998_pr)*R_T_inv + T_log *(3.02_pr))
    k(r1b) = (17729290.57642423_pr)*exp((-40271.85521197939_pr)*R_T_inv + T_log *(-0.12352841395486043_pr))
    k(r2b) = (0.022250519438231278_pr)*exp((-19900.25454070816_pr)*R_T_inv + T_log *(2.6908001930042906_pr))
    k(r5b) = (440.1782345473235_pr)*exp((-76020.25839665145_pr)*R_T_inv + T_log *(1.5920758925606673_pr))
    k(r6b) = (0.20837492457607695_pr)*exp((-59130.015228558455_pr)*R_T_inv + T_log *(2.491275712915376_pr))
    k(r8b) = (134265599.92056495_pr)*exp((-1342.102561897014_pr)*R_T_inv + T_log *(-0.3530558573624177_pr))
    k(r17b) = (731155417.4829124_pr)*exp((-879.8905518579306_pr)*R_T_inv + T_log *(-0.18088610001311003_pr))
    k(r22b) = (4.246927330556553_pr)*exp((-86350.34964081323_pr)*R_T_inv + T_log *(2.2629620314496512_pr))
    k(r23b) = (21.224414220257525_pr)*exp((-59260.37912534472_pr)*R_T_inv + T_log *(1.7782952171844644_pr))
    k(r27b) = (3512073963.061712_pr)*exp((-66772.26715985876_pr)*R_T_inv + T_log *(-0.34523313741984507_pr))
    k(r31b) = (35458579.34225838_pr)*exp((-37761.873578128936_pr)*R_T_inv + T_log *(-0.12352840716760596_pr))
    k_0(FOr33b) = (9.53636000109e+32_pr)*exp((-458465.466789_pr)*R_T_inv + T_log *(-5.8277005324_pr))
    k_inf(FOr33b) = (1.90177374962e+14_pr)*exp((-437235.473604_pr)*R_T_inv + T_log *(0.652299592632_pr))
    k(r33b) = getlindratecoeff(Tloc,k_0(FOr33b),k_inf(FOr33b),1.0_pr, M(mM33),Ploc)
    k(r37b) = (3.356640619771351_pr)*exp((-55339.87393280361_pr)*R_T_inv + T_log *(1.903780655129222_pr))
    k(r38b) = (32658975742.977776_pr)*exp((-827.983877978271_pr)*R_T_inv + T_log *(-0.9126909227024629_pr))
    k(r45b) = (61128097067.41484_pr)*exp((-4464.386899160203_pr)*R_T_inv + T_log *(-0.7159500745832347_pr))
    k(r47b) = (33.42867788060587_pr)*exp((-40724.971805766494_pr)*R_T_inv + T_log *(1.5635215188671316_pr))
    k(r48b) = (25.043952352686045_pr)*exp((-25055.24316953238_pr)*R_T_inv + T_log *(1.4343217229101575_pr))
    k(r49b) = (12.887175681807777_pr)*exp((-70085.2203970672_pr)*R_T_inv + T_log *(1.6255974043520434_pr))
    k(r51b) = (5.289029676794421_pr)*exp((-58975.37273243503_pr)*R_T_inv + T_log *(1.7218167594321432_pr))
    k(r53b) = (752754777491.3809_pr)*exp((-136838.3465694232_pr)*R_T_inv + T_log *(-0.9134395457245326_pr))
    k(r54b) = (82614856363.484_pr)*exp((-106008.33352628726_pr)*R_T_inv + T_log *(-1.0234395771598224_pr))
    k(r56b) = (10637573.509423345_pr)*exp((-37761.87355985963_pr)*R_T_inv + T_log *(-0.12352840350099344_pr))
    k(r61b) = (333279.14118726214_pr)*exp((-7462.737306276717_pr)*R_T_inv + T_log *(-1.0522126253409803_pr))
    k(r62b) = (3992220.2419204456_pr)*exp((-7462.737224414295_pr)*R_T_inv + T_log *(-1.0522126091681607_pr))
    k(r65b) = (1.6754590000232583_pr)*exp((-77640.90615077084_pr)*R_T_inv + T_log *(1.9858963092417024_pr))
    k(r70b) = (8273668.813509312_pr)*exp((-37761.873602245134_pr)*R_T_inv + T_log *(-0.12352841198991399_pr))
    k(r75b) = (1831.4030540628216_pr)*exp((-125628.39629633297_pr)*R_T_inv + T_log *(1.8751721823583902_pr))
    k_0(FOr76b) = (5.05500301487e+23_pr)*exp((-159505.418959_pr)*R_T_inv + T_log *(-4.06608812934_pr))
    k_inf(FOr76b) = (136341548.576_pr)*exp((-155015.401124_pr)*R_T_inv + T_log *(1.86391188567_pr))
    FC(FOr76b) = ((1.0_pr - 0.2122_pr)*exp(-Tloc/(1.0_pr))) + (0.2122_pr*exp(-Tloc/(-10210.0_pr)))
    k(r76b) = getlindratecoeff(Tloc,k_0(FOr76b),k_inf(FOr76b),FC(FOr76b), M(mM76),Ploc)
    k(r83b) = (17238180146751.81_pr)*exp((-104729.92548368075_pr)*R_T_inv + T_log *(-1.3080136227775843_pr))
    k(r92b) = (2.0696032153577708e+20_pr)*exp((-99688.90744280803_pr)*R_T_inv + T_log *(-2.978026526719999_pr))
    k_0(FOr93b) = (6.23917690211e+32_pr)*exp((-173773.804876_pr)*R_T_inv + T_log *(-6.53244020377_pr))
    k_inf(FOr93b) = (421067646.691_pr)*exp((-155303.787984_pr)*R_T_inv + T_log *(1.56755978268_pr))
    FC(FOr93b) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299.0_pr))) + (-0.569_pr*exp(-Tloc/(-9147.0_pr))) + exp(-(152.4)/Tloc)
    k(r93b) = getlindratecoeff(Tloc,k_0(FOr93b),k_inf(FOr93b),FC(FOr93b), M(mM93),Ploc)
    k_0(FOr101b) = (1.59592131765e+63_pr)*exp((-176791.070174_pr)*R_T_inv + T_log *(-15.9354800842_pr))
    k_inf(FOr101b) = (1.35653343924e+12_pr)*exp((-124621.06721_pr)*R_T_inv + T_log *(0.264519884531_pr))
    FC(FOr101b) = ((1.0_pr - 0.1894_pr)*exp(-Tloc/(277.0_pr))) + (0.1894_pr*exp(-Tloc/(8748.0_pr))) + exp(-(7891.0)/Tloc)
    k(r101b) = getlindratecoeff(Tloc,k_0(FOr101b),k_inf(FOr101b),FC(FOr101b), M(mM101),Ploc)
    k_0(FOr102b) = (9.64784500939e+36_pr)*exp((-446599.312153_pr)*R_T_inv + T_log *(-6.8773781844_pr))
    k_inf(FOr102b) = (2.52589236131e+19_pr)*exp((-425239.322522_pr)*R_T_inv + T_log *(-0.787378146066_pr))
    FC(FOr102b) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125.0_pr))) + (0.8422_pr*exp(-Tloc/(2219.0_pr))) + exp(-(6882.0)/Tloc)
    k(r102b) = getlindratecoeff(Tloc,k_0(FOr102b),k_inf(FOr102b),FC(FOr102b), M(mM102),Ploc)
    k_0(FOr108b) = (8.58580929884e+47_pr)*exp((-46266.2899098_pr)*R_T_inv + T_log *(-12.5266718371_pr))
    k_inf(FOr108b) = (4.33906480362e+38_pr)*exp((-70376.2970124_pr)*R_T_inv + T_log *(-9.10667183375_pr))
    FC(FOr108b) = ((1.0_pr - 0.39_pr)*exp(-Tloc/(100.0_pr))) + (0.39_pr*exp(-Tloc/(1900.0_pr))) + exp(-(6000.0)/Tloc)
    k(r108b) = getlindratecoeff(Tloc,k_0(FOr108b),k_inf(FOr108b),FC(FOr108b), M(mM108),Ploc)
    k(r120b) = (1.46534506175671e+30_pr)*exp((-389561.49330285477_pr)*R_T_inv + T_log *(-4.558931852346789_pr))
    k(r121b) = (1.780726319419713e+30_pr)*exp((-383809.1854230486_pr)*R_T_inv + T_log *(-4.572476349264602_pr))
    k(r126b) = (1.766672011437014e+55_pr)*exp((-667591.9831686344_pr)*R_T_inv + T_log *(-10.0741497423499_pr))
    k(r127b) = (8.216557984874969e+38_pr)*exp((-256049.9808476945_pr)*R_T_inv + T_log *(-7.428217499753913_pr))
    k(r128b) = (4.88702418541112e+17_pr)*exp((-347771.4212910579_pr)*R_T_inv + T_log *(-0.3161205409612194_pr))
    k(r129b) = (144.3154826853027_pr)*exp((-48281.30441142663_pr)*R_T_inv + T_log *(1.3175333231249051_pr))
    k(r130b) = (0.006411030595949381_pr)*exp((-87398.9765759684_pr)*R_T_inv + T_log *(2.6147528696899225_pr))
    k(r136b) = (1222.5573062172061_pr)*exp((-77613.59037579132_pr)*R_T_inv + T_log *(0.9310778031902854_pr))
    k(r137b) = (2.5413884356224326e+39_pr)*exp((-333342.26983105286_pr)*R_T_inv + T_log *(-7.726455509416434_pr))
    k(r138b) = (808921352.865255_pr)*exp((-32052.262498609212_pr)*R_T_inv + T_log *(-0.25645551495324115_pr))
    k(r142b) = (2.347492362835063e+50_pr)*exp((-323927.9330143372_pr)*R_T_inv + T_log *(-10.95303685226093_pr))
    k_0(FOr148b) = (4.74236001288e+32_pr)*exp((-164497.074934_pr)*R_T_inv + T_log *(-6.58837622385_pr))
    k_inf(FOr148b) = (2.31815074823e+14_pr)*exp((-152077.080228_pr)*R_T_inv + T_log *(-0.298376242578_pr))
    FC(FOr148b) = ((1.0_pr - 1.0_pr)*exp(-Tloc/(1000.0_pr))) + (1.0_pr*exp(-Tloc/(1310.0_pr))) + exp(-(48100.0)/Tloc)
    k(r148b) = getlindratecoeff(Tloc,k_0(FOr148b),k_inf(FOr148b),FC(FOr148b), M(mM148),Ploc)
    k(r149b) = (3.503876166348925e+29_pr)*exp((-39119.12191174787_pr)*R_T_inv + T_log *(-6.574220020123834_pr))
    k(r150b) = (6.376134983840494e+53_pr)*exp((-459711.95724727033_pr)*R_T_inv + T_log *(-11.12775344532129_pr))
    k(r151b) = (11392523160.47694_pr)*exp((-81216.01983118497_pr)*R_T_inv + T_log *(-0.9828961297086806_pr))
    k(r164b) = (6.809139754447457e+60_pr)*exp((-141565.05622061025_pr)*R_T_inv + T_log *(-15.03692365185586_pr))
    k(r175b) = (50061850534939.56_pr)*exp((-189796.01202681824_pr)*R_T_inv + T_log *(0.48997256714460014_pr))
    k(r176b) = (1.9444173759417796e+53_pr)*exp((-146537.95070946927_pr)*R_T_inv + T_log *(-12.674067755908245_pr))
    k(r182b) = (6.331619898516303_pr)*exp((-50661.87187748445_pr)*R_T_inv + T_log *(1.6119535021135627_pr))
    k(r183b) = (1.8029840064293446e+53_pr)*exp((-219059.30144182494_pr)*R_T_inv + T_log *(-12.199999259932445_pr))
    k(r187b) = (1.7486198551326773e-05_pr)*exp((-52202.686753638045_pr)*R_T_inv + T_log *(3.4058671578737547_pr))
    k(r188b) = (46117939350178.12_pr)*exp((-25151.33257022035_pr)*R_T_inv + T_log *(-1.5082013355170203_pr))
    k(r189b) = (2.4645459288773108e+16_pr)*exp((-22660.637477583638_pr)*R_T_inv + T_log *(-2.5382005845933087_pr))
    k(r196b) = (393194741040.2435_pr)*exp((-380303.2716716242_pr)*R_T_inv + T_log *(1.8349846921461765_pr))
    k(r197b) = (4.49494055065156e+32_pr)*exp((-378525.39204913995_pr)*R_T_inv + T_log *(-5.64192621967152_pr))
    k(r198b) = (7.290696451093596e+52_pr)*exp((-135820.90026390165_pr)*R_T_inv + T_log *(-12.568445754131329_pr))
    k(r199b) = (0.6378985628862892_pr)*exp((-19251.346874145205_pr)*R_T_inv + T_log *(2.1273753719381148_pr))
    k(r200b) = (1.647651406574653e-06_pr)*exp((-16631.610417983342_pr)*R_T_inv + T_log *(3.6594511479340737_pr))
    k(r201b) = (60348597353.743645_pr)*exp((-43026.75087625521_pr)*R_T_inv + T_log *(-0.8205146468469415_pr))
    k(r202b) = (1802310702339366.5_pr)*exp((-301635.3741169714_pr)*R_T_inv + T_log *(-0.8856508534641376_pr))
    k(r203b) = (1.5523562570439065e+53_pr)*exp((-108905.55911481557_pr)*R_T_inv + T_log *(-12.648324890460305_pr))
    k(r204b) = (0.665834920538346_pr)*exp((6513.980374038968_pr)*R_T_inv + T_log *(2.1174963032604324_pr))
    k(r205b) = (6.156446062933493e-06_pr)*exp((-1926.2691613862103_pr)*R_T_inv + T_log *(3.3495720014168837_pr))
  
    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(sN2) * cqss(sqssSCH2) 
    w(r2f) = k(r2f) * c(sH2) * c(sO) 
    w(r3f) = k(r3f) * c(sH) * c(sO) * m(mM3)
    w(r4f) = k(r4f) * c(sH) * c(sOH) * m(mM4)
    w(r5f) = k(r5f) * c(sH2) * c(sOH) 
    w(r6f) = k(r6f) * c(sOH)**2.0_pr 
    w(r7f) = k(r7f) * c(sH) * c(sO2) 
    w(r8f) = k(r8f) * c(sH) * c(sO2) 
    w(r9f) = k(r9f) * c(sH) * c(sHO2) 
    w(r10f) = k(r10f) * c(sH) * c(sHO2) 
    w(r11f) = k(r11f) * c(sH) * c(sHO2) 
    w(r12f) = k(r12f) * c(sHO2) * c(sO) 
    w(r13f) = k(r13f) * c(sHO2) * c(sOH) 
    w(r14f) = k(r14f) * c(sHO2) * c(sOH) 
    w(r15f) = k(r15f) * cqss(sqssCH) * c(sO) 
    w(r16f) = k(r16f) * cqss(sqssCH) * c(sOH) 
    w(r17f) = k(r17f) * c(sH2) * cqss(sqssCH) 
    w(r18f) = k(r18f) * c(sH2O) * cqss(sqssCH) 
    w(r19f) = k(r19f) * cqss(sqssCH) * c(sO2) 
    w(r20f) = k(r20f) * cqss(sqssTCH2) * c(sO) 
    w(r21f) = k(r21f) * cqss(sqssTCH2) * c(sOH) 
    w(r22f) = k(r22f) * cqss(sqssTCH2) * c(sOH) 
    w(r23f) = k(r23f) * c(sH2) * cqss(sqssTCH2) 
    w(r24f) = k(r24f) * cqss(sqssTCH2) * c(sO2) 
    w(r25f) = k(r25f) * cqss(sqssTCH2) * c(sO2) 
    w(r26f) = k(r26f) * cqss(sqssTCH2) * c(sO2) 
    w(r27f) = k(r27f) * c(sH2) * cqss(sqssSCH2) 
    w(r28f) = k(r28f) * cqss(sqssSCH2) * c(sO2) 
    w(r29f) = k(r29f) * cqss(sqssSCH2) * c(sO2) 
    w(r30f) = k(r30f) * c(sH2O) * cqss(sqssSCH2) 
    w(r31f) = k(r31f) * c(sH2O) * cqss(sqssSCH2) 
    w(r32f) = k(r32f) * c(sH2O) * cqss(sqssSCH2) 
    w(r33f) = k(r33f) * c(sH) * c(sCH3) 
    w(r34f) = k(r34f) * c(sCH3) * c(sO) 
    w(r35f) = k(r35f) * c(sCH3) * c(sO) 
    w(r36f) = k(r36f) * c(sCH3) * c(sOH) 
    w(r37f) = k(r37f) * c(sCH3) * c(sOH) 
    w(r38f) = k(r38f) * c(sCH3) * c(sOH) 
    w(r39f) = k(r39f) * c(sCH3) * c(sO2) 
    w(r40f) = k(r40f) * c(sCH3) * c(sO2) 
    w(r41f) = k(r41f) * c(sCH3) * c(sHO2) 
    w(r42f) = k(r42f) * c(sCH3) * c(sHO2) 
    w(r43f) = k(r43f) * cqss(sqssCH) * c(sCH3) 
    w(r44f) = k(r44f) * cqss(sqssTCH2) * c(sCH3) 
    w(r45f) = k(r45f) * c(sCH3)**2.0_pr 
    w(r46f) = k(r46f) * c(sCH3)**2.0_pr 
    w(r47f) = k(r47f) * c(sH) * c(sCH4) 
    w(r48f) = k(r48f) * c(sCH4) * c(sO) 
    w(r49f) = k(r49f) * c(sCH4) * c(sOH) 
    w(r50f) = k(r50f) * cqss(sqssCH) * c(sCH4) 
    w(r51f) = k(r51f) * cqss(sqssTCH2) * c(sCH4) 
    w(r52f) = k(r52f) * c(sCO) * c(sO) 
    w(r53f) = k(r53f) * c(sCO) * c(sOH) 
    w(r54f) = k(r54f) * c(sCO) * c(sOH) 
    w(r55f) = k(r55f) * c(sCO) * c(sHO2) 
    w(r56f) = k(r56f) * c(sCO) * cqss(sqssSCH2) 
    w(r57f) = k(r57f) * c(sH) * cqss(sqssHCO) 
    w(r58f) = k(r58f) * cqss(sqssHCO) * c(sO) 
    w(r59f) = k(r59f) * cqss(sqssHCO) * c(sO) 
    w(r60f) = k(r60f) * cqss(sqssHCO) * c(sOH) 
    w(r61f) = k(r61f) * cqss(sqssHCO) * m(mM61)
    w(r62f) = k(r62f) * c(sH2O) * cqss(sqssHCO) 
    w(r63f) = k(r63f) * cqss(sqssHCO) * c(sO2) 
    w(r64f) = k(r64f) * c(sCH3) * cqss(sqssHCO) 
    w(r65f) = k(r65f) * c(sCH2O) * c(sH) 
    w(r66f) = k(r66f) * c(sCH2O) * c(sO) 
    w(r67f) = k(r67f) * c(sCH2O) * c(sOH) 
    w(r68f) = k(r68f) * c(sCH2O) * c(sCH3) 
    w(r69f) = k(r69f) * cqss(sqssCH) * c(sCO2) 
    w(r70f) = k(r70f) * c(sCO2) * cqss(sqssSCH2) 
    w(r71f) = k(r71f) * c(sCO2) * cqss(sqssSCH2) 
    w(r72f) = k(r72f) * cqss(sqssC2H) * c(sO) 
    w(r73f) = k(r73f) * cqss(sqssC2H) * c(sOH) 
    w(r74f) = k(r74f) * cqss(sqssC2H) * c(sO2) 
    w(r75f) = k(r75f) * c(sH2) * cqss(sqssC2H) 
    w(r76f) = k(r76f) * c(sH) * c(sC2H2) 
    w(r77f) = k(r77f) * c(sC2H2) * c(sO) 
    w(r78f) = k(r78f) * c(sC2H2) * c(sO) 
    w(r79f) = k(r79f) * c(sC2H2) * c(sO) 
    w(r80f) = k(r80f) * c(sC2H2) * c(sOH) 
    w(r81f) = k(r81f) * c(sC2H2) * c(sOH) 
    w(r82f) = k(r82f) * c(sC2H2) * c(sOH) 
    w(r83f) = k(r83f) * cqss(sqssSCH2) * c(sC2H2) 
    w(r84f) = k(r84f) * c(sH) * cqss(sqssC2H3) 
    w(r85f) = k(r85f) * c(sH) * cqss(sqssC2H3) 
    w(r86f) = k(r86f) * cqss(sqssC2H3) * c(sO) 
    w(r87f) = k(r87f) * cqss(sqssC2H3) * c(sOH) 
    w(r88f) = k(r88f) * cqss(sqssC2H3) * c(sO2) 
    w(r89f) = k(r89f) * cqss(sqssC2H3) * c(sO2) 
    w(r90f) = k(r90f) * cqss(sqssC2H3) * c(sO2) 
    w(r91f) = k(r91f) * c(sCH3) * cqss(sqssC2H3) 
    w(r92f) = k(r92f) * c(sCH3) * cqss(sqssC2H3) 
    w(r93f) = k(r93f) * c(sH) * c(sC2H4) 
    w(r94f) = k(r94f) * c(sH) * c(sC2H4) 
    w(r95f) = k(r95f) * c(sC2H4) * c(sO) 
    w(r96f) = k(r96f) * c(sC2H4) * c(sO) 
    w(r97f) = k(r97f) * c(sC2H4) * c(sO) 
    w(r98f) = k(r98f) * c(sC2H4) * c(sOH) 
    w(r99f) = k(r99f) * c(sC2H4) * c(sOH) 
    w(r100f) = k(r100f) * c(sCH3) * c(sC2H4) 
    w(r101f) = k(r101f) * c(sCH3) * c(sC2H4) 
    w(r102f) = k(r102f) * c(sH) * c(sC2H5) 
    w(r103f) = k(r103f) * c(sC2H5) * c(sO) 
    w(r104f) = k(r104f) * c(sC2H5) * c(sO2) 
    w(r105f) = k(r105f) * c(sC2H5) * cqss(sqssHCO) 
    w(r106f) = k(r106f) * c(sC2H5) * c(sHO2) 
    w(r107f) = k(r107f) * c(sC2H5) * c(sHO2) 
    w(r108f) = k(r108f) * c(sC2H6) 
    w(r109f) = k(r109f) * c(sC2H6) * c(sH) 
    w(r110f) = k(r110f) * c(sC2H6) * c(sO) 
    w(r111f) = k(r111f) * c(sC2H6) * c(sOH) 
    w(r112f) = k(r112f) * c(sC2H6) * c(sCH3) 
    w(r113f) = k(r113f) * cqss(sqssHCCO) 
    w(r114f) = k(r114f) * c(sH) * cqss(sqssHCCO) 
    w(r115f) = k(r115f) * cqss(sqssHCCO) * c(sO) 
    w(r116f) = k(r116f) * cqss(sqssHCCO) * c(sO2) 
    w(r117f) = k(r117f) * c(sC2H2) * cqss(sqssHCCO) 
    w(r118f) = k(r118f) * c(sCH3) * cqss(sqssHCCO) 
    w(r119f) = k(r119f) * cqss(sqssHCCO) * c(sOH) 
    w(r120f) = k(r120f) * c(sH) * c(sC3H3) 
    w(r121f) = k(r121f) * c(sH) * c(sC3H3) 
    w(r122f) = k(r122f) * c(sC3H3) * c(sOH) 
    w(r123f) = k(r123f) * c(sC3H3) * c(sO) 
    w(r124f) = k(r124f) * c(sC3H3) * c(sO2) 
    w(r125f) = k(r125f) * c(sHO2) * c(sC3H3) 
    w(r126f) = k(r126f) * c(sC3H3)**2.0_pr 
    w(r127f) = k(r127f) * c(sC3H3)**2.0_pr 
    w(r128f) = k(r128f) * c(sC2H2) * c(sC3H3) 
    w(r129f) = k(r129f) * c(sH) * c(sPC3H4) 
    w(r130f) = k(r130f) * c(sH) * c(sPC3H4) 
    w(r131f) = k(r131f) * c(sPC3H4) * c(sOH) 
    w(r132f) = k(r132f) * c(sCH3) * c(sPC3H4) 
    w(r133f) = k(r133f) * c(sPC3H4) * c(sO) 
    w(r134f) = k(r134f) * c(sPC3H4) * c(sO) 
    w(r135f) = k(r135f) * c(sPC3H4) * c(sOH) 
    w(r136f) = k(r136f) * c(sH) * c(sAC3H4) 
    w(r137f) = k(r137f) * c(sAC3H4) 
    w(r138f) = k(r138f) * c(sH) * c(sAC3H4) 
    w(r139f) = k(r139f) * c(sH) * c(sAC3H4) 
    w(r140f) = k(r140f) * c(sAC3H4) * c(sOH) 
    w(r141f) = k(r141f) * c(sAC3H4) * c(sCH3) 
    w(r142f) = k(r142f) * c(sH) * c(sAC3H4) 
    w(r143f) = k(r143f) * c(sH) * c(sAC3H5) 
    w(r144f) = k(r144f) * c(sAC3H5) * c(sOH) 
    w(r145f) = k(r145f) * c(sAC3H5) * c(sCH3) 
    w(r146f) = k(r146f) * c(sAC3H5) * c(sO2) 
    w(r147f) = k(r147f) * c(sAC3H5) * c(sC3H3) 
    w(r148f) = k(r148f) * c(sC3H6) * c(sH) 
    w(r149f) = k(r149f) * c(sC3H6) 
    w(r150f) = k(r150f) * c(sH) * c(sAC3H5) 
    w(r151f) = k(r151f) * c(sC3H6) * c(sH) 
    w(r152f) = k(r152f) * c(sC3H6) * c(sO) 
    w(r153f) = k(r153f) * c(sC3H6) * c(sO) 
    w(r154f) = k(r154f) * c(sC3H6) * c(sH) 
    w(r155f) = k(r155f) * c(sC3H6) * c(sO) 
    w(r156f) = k(r156f) * c(sC3H6) * c(sOH) 
    w(r157f) = k(r157f) * c(sC3H6) * c(sCH3) 
    w(r158f) = k(r158f) * c(sHO2) * c(sC5H5) 
    w(r159f) = k(r159f) * c(sC5H5)**2.0_pr 
    w(r160f) = k(r160f) * c(sC5H5) * c(sO) 
    w(r161f) = k(r161f) * c(sHO2) * c(sC5H5) 
    w(r162f) = k(r162f) * c(sHO2) * c(sC5H5) 
    w(r163f) = k(r163f) * c(sC5H5) * c(sOH) 
    w(r164f) = k(r164f) * c(sC5H6) 
    w(r165f) = k(r165f) * c(sH) * c(sC5H6) 
    w(r166f) = k(r166f) * c(sH) * c(sC5H6) 
    w(r167f) = k(r167f) * c(sC5H6) * c(sO) 
    w(r168f) = k(r168f) * c(sC5H6) * c(sOH) 
    w(r169f) = k(r169f) * c(sCH3) * c(sC5H6) 
    w(r170f) = k(r170f) * c(sO2) * cqss(sqssC6H5) 
    w(r171f) = k(r171f) * c(sO2) * cqss(sqssC6H5) 
    w(r172f) = k(r172f) * cqss(sqssC6H5) * c(sO) 
    w(r173f) = k(r173f) * cqss(sqssC6H5) * c(sOH) 
    w(r174f) = k(r174f) * c(sHO2) * cqss(sqssC6H5) 
    w(r175f) = k(r175f) * c(sC2H2) * cqss(sqssC6H5) 
    w(r176f) = k(r176f) * c(sC6H6) 
    w(r177f) = k(r177f) * c(sH) * c(sC6H6) 
    w(r178f) = k(r178f) * c(sC6H6) * c(sOH) 
    w(r179f) = k(r179f) * c(sCH3) * c(sC6H6) 
    w(r180f) = k(r180f) * c(sC6H6) * c(sO) 
    w(r181f) = k(r181f) * c(sC6H6) * c(sOH) 
    w(r182f) = k(r182f) * cqss(sqssC7H7) 
    w(r183f) = k(r183f) * c(sH) * cqss(sqssC7H7) 
    w(r184f) = k(r184f) * cqss(sqssC7H7) * c(sO) 
    w(r185f) = k(r185f) * c(sHO2) * cqss(sqssC7H7) 
    w(r186f) = k(r186f) * c(sC3H3) * cqss(sqssC7H7) 
    w(r187f) = k(r187f) * c(sH) * c(sC7H8) 
    w(r188f) = k(r188f) * c(sC7H8) 
    w(r189f) = k(r189f) * c(sC7H8) 
    w(r190f) = k(r190f) * c(sH) * c(sC7H8) 
    w(r191f) = k(r191f) * c(sC7H8) * c(sOH) 
    w(r192f) = k(r192f) * c(sC7H8) * c(sCH3) 
    w(r193f) = k(r193f) * cqss(sqssC7H6O) 
    w(r194f) = k(r194f) * c(sH) * cqss(sqssC7H6O) 
    w(r195f) = k(r195f) * cqss(sqssC7H6O) * c(sCH3) 
    w(r196f) = k(r196f) * cqss(sqssC8H5) * c(sC2H2) 
    w(r197f) = k(r197f) * cqss(sqssC8H5) * c(sC2H4) 
    w(r198f) = k(r198f) * c(sC8H6) 
    w(r199f) = k(r199f) * c(sH) * c(sC8H6) 
    w(r200f) = k(r200f) * c(sC8H6) * c(sOH) 
    w(r201f) = k(r201f) * cqss(sqssC8H7) 
    w(r202f) = k(r202f) * cqss(sqssC8H7) * c(sC2H2) 
    w(r203f) = k(r203f) * c(sC10H8) 
    w(r204f) = k(r204f) * c(sH) * c(sC10H8) 
    w(r205f) = k(r205f) * c(sC10H8) * c(sOH) 
    w(r1b) = k(r1b) * cqss(sqssTCH2) * c(sN2) 
    w(r2b) = k(r2b) * c(sH) * c(sOH) 
    w(r5b) = k(r5b) * c(sH) * c(sH2O) 
    w(r6b) = k(r6b) * c(sH2O) * c(sO) 
    w(r8b) = k(r8b) * c(sO) * c(sOH) 
    w(r17b) = k(r17b) * c(sH) * cqss(sqssTCH2) 
    w(r22b) = k(r22b) * c(sH2O) * cqss(sqssCH) 
    w(r23b) = k(r23b) * c(sH) * c(sCH3) 
    w(r27b) = k(r27b) * c(sH) * c(sCH3) 
    w(r31b) = k(r31b) * cqss(sqssTCH2) * c(sH2O) 
    w(r33b) = k(r33b) * c(sCH4) 
    w(r37b) = k(r37b) * cqss(sqssTCH2) * c(sH2O) 
    w(r38b) = k(r38b) * c(sH2O) * cqss(sqssSCH2) 
    w(r45b) = k(r45b) * c(sH) * c(sC2H5) 
    w(r47b) = k(r47b) * c(sH2) * c(sCH3) 
    w(r48b) = k(r48b) * c(sCH3) * c(sOH) 
    w(r49b) = k(r49b) * c(sH2O) * c(sCH3) 
    w(r51b) = k(r51b) * c(sCH3)**2.0_pr 
    w(r53b) = k(r53b) * c(sH) * c(sCO2) 
    w(r54b) = k(r54b) * c(sH) * c(sCO2) 
    w(r56b) = k(r56b) * cqss(sqssTCH2) * c(sCO) 
    w(r61b) = k(r61b) * c(sH) * c(sCO) * m(mM61)
    w(r62b) = k(r62b) * c(sH) * c(sH2O) * c(sCO) 
    w(r65b) = k(r65b) * c(sH2) * cqss(sqssHCO) 
    w(r70b) = k(r70b) * cqss(sqssTCH2) * c(sCO2) 
    w(r75b) = k(r75b) * c(sH) * c(sC2H2) 
    w(r76b) = k(r76b) * cqss(sqssC2H3) 
    w(r83b) = k(r83b) * c(sH) * c(sC3H3) 
    w(r92b) = k(r92b) * c(sH) * c(sAC3H5) 
    w(r93b) = k(r93b) * c(sC2H5) 
    w(r101b) = k(r101b) * cqss(sqssNC3H7) 
    w(r102b) = k(r102b) * c(sC2H6) 
    w(r108b) = k(r108b) * c(sCH3)**2.0_pr 
    w(r120b) = k(r120b) * c(sPC3H4) 
    w(r121b) = k(r121b) * c(sAC3H4) 
    w(r126b) = k(r126b) * c(sC6H6) 
    w(r127b) = k(r127b) * c(sH) * cqss(sqssC6H5) 
    w(r128b) = k(r128b) * c(sC5H5) 
    w(r129b) = k(r129b) * c(sCH3) * c(sC2H2) 
    w(r130b) = k(r130b) * c(sH2) * c(sC3H3) 
    w(r136b) = k(r136b) * c(sCH3) * c(sC2H2) 
    w(r137b) = k(r137b) * c(sPC3H4) 
    w(r138b) = k(r138b) * c(sH) * c(sPC3H4) 
    w(r142b) = k(r142b) * c(sAC3H5) 
    w(r148b) = k(r148b) * cqss(sqssNC3H7) 
    w(r149b) = k(r149b) * c(sCH3) * cqss(sqssC2H3) 
    w(r150b) = k(r150b) * c(sC3H6) 
    w(r151b) = k(r151b) * c(sCH3) * c(sC2H4) 
    w(r164b) = k(r164b) * c(sH) * c(sC5H5) 
    w(r175b) = k(r175b) * cqss(sqssC8H7) 
    w(r176b) = k(r176b) * c(sH) * cqss(sqssC6H5) 
    w(r182b) = k(r182b) * c(sC2H2) * c(sC5H5) 
    w(r183b) = k(r183b) * c(sCH3) * cqss(sqssC6H5) 
    w(r187b) = k(r187b) * c(sCH3) * c(sC6H6) 
    w(r188b) = k(r188b) * c(sH) * cqss(sqssC7H7) 
    w(r189b) = k(r189b) * c(sCH3) * cqss(sqssC6H5) 
    w(r196b) = k(r196b) * cqss(sqssC10H7) 
    w(r197b) = k(r197b) * c(sH) * c(sC10H8) 
    w(r198b) = k(r198b) * cqss(sqssC8H5) * c(sH) 
    w(r199b) = k(r199b) * cqss(sqssC8H5) * c(sH2) 
    w(r200b) = k(r200b) * cqss(sqssC8H5) * c(sH2O) 
    w(r201b) = k(r201b) * c(sH) * c(sC8H6) 
    w(r202b) = k(r202b) * c(sH) * c(sC10H8) 
    w(r203b) = k(r203b) * c(sH) * cqss(sqssC10H7) 
    w(r204b) = k(r204b) * c(sH2) * cqss(sqssC10H7) 
    w(r205b) = k(r205b) * cqss(sqssC10H7) * 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(sN2) = 0.0_pr 

    cdot(sO) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r6f) &
         - w(r6b) &
         + w(r8f) &
         - w(r8b) &
         + w(r10f) &
         - w(r12f) &
         - w(r15f) &
         + w(r19f) &
         - w(r20f) &
         + w(r25f) &
         - w(r34f) &
         - w(r35f) &
         + w(r39f) &
         - w(r48f) &
         + w(r48b) &
         - w(r52f) &
         - w(r58f) &
         - w(r59f) &
         - w(r66f) &
         - w(r72f) &
         - w(r77f) &
         - w(r78f) &
         - w(r79f) &
         - w(r86f) &
         + w(r89f) &
         - w(r95f) &
         - w(r96f) &
         - w(r97f) &
         - w(r103f) &
         - w(r110f) &
         - w(r115f) &
         - w(r123f) &
         - w(r133f) &
         - w(r134f) &
         - w(r152f) &
         - w(r153f) &
         - w(r155f) &
         - w(r160f) &
         - w(r167f) &
         + w(r170f) &
         - w(r172f) &
         - w(r180f) &
         - w(r184f) 

    cdot(sH2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         - w(r5f) &
         + w(r5b) &
         + w(r9f) &
         - w(r17f) &
         + w(r17b) &
         - w(r23f) &
         + w(r23b) &
         - w(r27f) &
         + w(r27b) &
         + w(r30f) &
         + w(r32f) &
         + w(r35f) &
         + w(r36f) &
         + w(r47f) &
         - w(r47b) &
         + w(r57f) &
         + w(r65f) &
         - w(r65b) &
         - w(r75f) &
         + w(r75b) &
         + w(r85f) &
         + w(r94f) &
         + w(r109f) &
         + w(r130f) &
         - w(r130b) &
         + w(r139f) &
         + w(r143f) &
         + w(r154f) &
         + w(r165f) &
         + w(r177f) &
         + w(r190f) &
         + w(r194f) &
         + w(r199f) &
         - w(r199b) &
         + w(r204f) &
         - w(r204b) 

    cdot(sH) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         - w(r3f) &
         - w(r4f) &
         + w(r5f) &
         - w(r5b) &
         - w(r7f) &
         - w(r8f) &
         + w(r8b) &
         - w(r9f) &
         - w(r10f) &
         - w(r11f) &
         + w(r15f) &
         + w(r16f) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         + w(r20f) &
         + w(r21f) &
         + w(r23f) &
         - w(r23b) &
         + 2.0_pr * w(r24f) &
         + w(r26f) &
         + w(r27f) &
         - w(r27b) &
         + w(r28f) &
         - w(r33f) &
         + w(r33b) &
         + w(r34f) &
         + w(r35f) &
         + w(r39f) &
         + w(r41f) &
         + w(r43f) &
         + w(r44f) &
         + w(r45f) &
         - w(r45b) &
         - w(r47f) &
         + w(r47b) &
         + w(r50f) &
         + w(r53f) &
         - w(r53b) &
         + w(r54f) &
         - w(r54b) &
         - w(r57f) &
         + w(r59f) &
         + w(r61f) &
         - w(r61b) &
         + w(r62f) &
         - w(r62b) &
         - w(r65f) &
         + w(r65b) &
         + w(r73f) &
         + w(r75f) &
         - w(r75b) &
         - w(r76f) &
         + w(r76b) &
         + w(r77f) &
         + w(r83f) &
         - w(r83b) &
         - w(r84f) &
         - w(r85f) &
         + w(r92f) &
         - w(r92b) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         + w(r95f) &
         - w(r102f) &
         + w(r102b) &
         - w(r109f) &
         - w(r114f) &
         + w(r115f) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r121b) &
         + w(r123f) &
         + w(r127f) &
         - w(r127b) &
         - w(r129f) &
         + w(r129b) &
         - w(r130f) &
         + w(r130b) &
         - w(r136f) &
         + w(r136b) &
         - w(r139f) &
         - w(r142f) &
         + w(r142b) &
         - w(r143f) &
         + 2.0_pr * w(r147f) &
         - w(r148f) &
         + w(r148b) &
         - w(r150f) &
         + w(r150b) &
         - w(r151f) &
         + w(r151b) &
         - w(r154f) &
         + 2.0_pr * w(r159f) &
         + w(r163f) &
         + w(r164f) &
         - w(r164b) &
         - w(r165f) &
         - w(r166f) &
         + w(r171f) &
         + w(r176f) &
         - w(r176b) &
         - w(r177f) &
         + w(r180f) &
         + w(r181f) &
         - w(r183f) &
         + w(r183b) &
         + w(r184f) &
         + w(r185f) &
         + 2.0_pr * w(r186f) &
         - w(r187f) &
         + w(r187b) &
         + w(r188f) &
         - w(r188b) &
         - w(r190f) &
         - w(r194f) &
         + w(r197f) &
         - w(r197b) &
         + w(r198f) &
         - w(r198b) &
         - w(r199f) &
         + w(r199b) &
         + w(r201f) &
         - w(r201b) &
         + w(r202f) &
         - w(r202b) &
         + w(r203f) &
         - w(r203b) &
         - w(r204f) &
         + w(r204b) 

    cdot(sOH) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         + w(r3f) &
         - w(r4f) &
         - w(r5f) &
         + w(r5b) &
         -2.0_pr * w(r6f) &
         + 2.0_pr * w(r6b) &
         + w(r8f) &
         - w(r8b) &
         + 2.0_pr * w(r11f) &
         + w(r12f) &
         - w(r13f) &
         - w(r14f) &
         - w(r16f) &
         - w(r21f) &
         - w(r22f) &
         + w(r22b) &
         + w(r26f) &
         + w(r28f) &
         - w(r36f) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         + w(r40f) &
         + w(r41f) &
         + w(r48f) &
         - w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r53f) &
         + w(r53b) &
         - w(r54f) &
         + w(r54b) &
         + w(r55f) &
         + w(r58f) &
         - w(r60f) &
         + w(r66f) &
         - w(r67f) &
         - w(r73f) &
         + w(r79f) &
         - w(r80f) &
         - w(r81f) &
         - w(r82f) &
         - w(r87f) &
         - w(r98f) &
         - w(r99f) &
         + w(r107f) &
         + w(r110f) &
         - w(r111f) &
         + w(r116f) &
         - w(r119f) &
         - w(r122f) &
         + w(r125f) &
         - w(r131f) &
         - w(r135f) &
         - w(r140f) &
         - w(r144f) &
         + w(r146f) &
         + w(r155f) &
         - w(r156f) &
         + w(r161f) &
         - w(r163f) &
         + w(r167f) &
         - w(r168f) &
         - w(r173f) &
         + w(r174f) &
         - w(r178f) &
         - w(r181f) &
         + w(r185f) &
         - w(r191f) &
         - w(r200f) &
         + w(r200b) &
         - w(r205f) &
         + w(r205b) 

    cdot(sH2O) = 0.0_pr &
         + w(r4f) &
         + w(r5f) &
         - w(r5b) &
         + w(r6f) &
         - w(r6b) &
         + w(r10f) &
         + w(r13f) &
         + w(r14f) &
         - w(r18f) &
         + w(r22f) &
         - w(r22b) &
         + w(r29f) &
         - w(r30f) &
         - w(r32f) &
         + w(r37f) &
         - w(r37b) &
         + w(r38f) &
         - w(r38b) &
         + w(r49f) &
         - w(r49b) &
         + w(r60f) &
         + w(r67f) &
         + w(r80f) &
         + w(r87f) &
         + w(r98f) &
         + w(r111f) &
         + w(r131f) &
         + w(r140f) &
         + w(r144f) &
         + w(r156f) &
         + w(r162f) &
         + w(r168f) &
         + w(r178f) &
         + w(r191f) &
         + w(r200f) &
         - w(r200b) &
         + w(r205f) &
         - w(r205b) 

    cdot(sO2) = 0.0_pr &
         - w(r7f) &
         - w(r8f) &
         + w(r8b) &
         + w(r9f) &
         + w(r12f) &
         + w(r13f) &
         + w(r14f) &
         - w(r19f) &
         - w(r24f) &
         - w(r25f) &
         - w(r26f) &
         - w(r28f) &
         - w(r29f) &
         - w(r39f) &
         - w(r40f) &
         + w(r42f) &
         - w(r63f) &
         - w(r74f) &
         - w(r88f) &
         - w(r89f) &
         - w(r90f) &
         - w(r104f) &
         + w(r106f) &
         - w(r116f) &
         - w(r124f) &
         - w(r146f) &
         + w(r158f) &
         - w(r170f) &
         - w(r171f) 

    cdot(sHO2) = 0.0_pr &
         + w(r7f) &
         - w(r9f) &
         - w(r10f) &
         - w(r11f) &
         - w(r12f) &
         - w(r13f) &
         - w(r14f) &
         - w(r41f) &
         - w(r42f) &
         - w(r55f) &
         + w(r63f) &
         + w(r88f) &
         + w(r104f) &
         - w(r106f) &
         - w(r107f) &
         - w(r125f) &
         - w(r158f) &
         - w(r161f) &
         - w(r162f) &
         - w(r174f) &
         - w(r185f) 

    cdot(sCO) = 0.0_pr &
         + w(r15f) &
         + w(r26f) &
         + w(r28f) &
         + w(r29f) &
         + w(r35f) &
         - w(r52f) &
         - w(r53f) &
         + w(r53b) &
         - w(r54f) &
         + w(r54b) &
         - w(r55f) &
         + w(r57f) &
         + w(r58f) &
         + w(r60f) &
         + w(r61f) &
         - w(r61b) &
         + w(r62f) &
         - w(r62b) &
         + w(r63f) &
         + w(r64f) &
         + w(r69f) &
         + w(r71f) &
         + w(r72f) &
         + w(r74f) &
         + w(r78f) &
         + w(r81f) &
         + w(r82f) &
         + w(r86f) &
         + w(r89f) &
         + w(r95f) &
         + w(r105f) &
         + w(r113f) &
         + w(r114f) &
         + 2.0_pr * w(r115f) &
         + 2.0_pr * w(r116f) &
         + w(r117f) &
         + w(r118f) &
         + w(r122f) &
         + w(r123f) &
         + 2.0_pr * w(r124f) &
         + w(r125f) &
         + w(r134f) &
         + w(r135f) &
         + w(r152f) &
         + w(r160f) &
         + w(r161f) &
         + w(r162f) &
         + w(r163f) &
         + w(r170f) &
         + 2.0_pr * w(r171f) &
         + w(r172f) &
         + w(r173f) &
         + w(r174f) &
         + w(r180f) &
         + w(r181f) &
         + w(r194f) &
         + w(r195f) 

    cdot(sCH2O) = 0.0_pr &
         + w(r18f) &
         + w(r21f) &
         + w(r25f) &
         + w(r30f) &
         + w(r32f) &
         + w(r34f) &
         + w(r36f) &
         + w(r39f) &
         + w(r40f) &
         + w(r41f) &
         - w(r65f) &
         + w(r65b) &
         - w(r66f) &
         - w(r67f) &
         - w(r68f) &
         + w(r71f) &
         + w(r90f) &
         + w(r96f) &
         + w(r99f) &
         + w(r103f) &
         + w(r107f) &
         + w(r146f) 

    cdot(sCH3) = 0.0_pr &
         + w(r23f) &
         - w(r23b) &
         + w(r27f) &
         - w(r27b) &
         - w(r33f) &
         + w(r33b) &
         - w(r34f) &
         - w(r35f) &
         - w(r36f) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r39f) &
         - w(r40f) &
         - w(r41f) &
         - w(r42f) &
         - w(r43f) &
         - w(r44f) &
         -2.0_pr * w(r45f) &
         + 2.0_pr * w(r45b) &
         -2.0_pr * w(r46f) &
         + w(r47f) &
         - w(r47b) &
         + w(r48f) &
         - w(r48b) &
         + w(r49f) &
         - w(r49b) &
         + 2.0_pr * w(r51f) &
         - 2.0_pr * w(r51b) &
         - w(r64f) &
         - w(r68f) &
         + w(r81f) &
         + w(r82f) &
         + w(r86f) &
         + w(r89f) &
         - w(r91f) &
         - w(r92f) &
         + w(r92b) &
         + w(r95f) &
         + w(r97f) &
         + w(r99f) &
         - w(r100f) &
         - w(r101f) &
         + w(r101b) &
         + w(r103f) &
         + w(r107f) &
         + 2.0_pr * w(r108f) &
         - 2.0_pr * w(r108b) &
         - w(r112f) &
         - w(r118f) &
         + w(r124f) &
         + w(r129f) &
         - w(r129b) &
         - w(r132f) &
         + w(r133f) &
         + w(r136f) &
         - w(r136b) &
         - w(r141f) &
         - w(r145f) &
         + w(r149f) &
         - w(r149b) &
         + w(r151f) &
         - w(r151b) &
         + 2.0_pr * w(r152f) &
         - w(r157f) &
         - w(r169f) &
         - w(r179f) &
         + w(r183f) &
         - w(r183b) &
         + w(r187f) &
         - w(r187b) &
         + w(r189f) &
         - w(r189b) &
         - w(r192f) &
         - w(r195f) 

    cdot(sCO2) = 0.0_pr &
         + w(r24f) &
         + w(r52f) &
         + w(r53f) &
         - w(r53b) &
         + w(r54f) &
         - w(r54b) &
         + w(r55f) &
         + w(r59f) &
         - w(r69f) &
         - w(r71f) 

    cdot(sCH4) = 0.0_pr &
         + w(r33f) &
         - w(r33b) &
         + w(r42f) &
         + w(r46f) &
         - w(r47f) &
         + w(r47b) &
         - w(r48f) &
         + w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         - w(r51f) &
         + w(r51b) &
         + w(r64f) &
         + w(r68f) &
         + w(r91f) &
         + w(r100f) &
         + w(r112f) &
         + w(r132f) &
         + w(r141f) &
         + w(r145f) &
         + w(r157f) &
         + w(r169f) &
         + w(r179f) &
         + w(r192f) &
         + w(r195f) 

    cdot(sC2H4) = 0.0_pr &
         + w(r44f) &
         + w(r50f) &
         + w(r84f) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         - w(r95f) &
         - w(r96f) &
         - w(r97f) &
         - w(r98f) &
         - w(r99f) &
         - w(r100f) &
         - w(r101f) &
         + w(r101b) &
         + w(r104f) &
         + w(r118f) &
         + w(r122f) &
         + w(r134f) &
         + w(r151f) &
         - w(r151b) &
         - w(r197f) &
         + w(r197b) 

    cdot(sC2H5) = 0.0_pr &
         + w(r45f) &
         - w(r45b) &
         + w(r93f) &
         - w(r93b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         - w(r104f) &
         - w(r105f) &
         - w(r106f) &
         - w(r107f) &
         + w(r109f) &
         + w(r110f) &
         + w(r111f) &
         + w(r112f) &
         + w(r135f) &
         + w(r153f) 

    cdot(sC2H2) = 0.0_pr &
         + w(r75f) &
         - w(r75b) &
         - w(r76f) &
         + w(r76b) &
         - w(r77f) &
         - w(r78f) &
         - w(r79f) &
         - w(r80f) &
         - w(r81f) &
         - w(r82f) &
         - w(r83f) &
         + w(r83b) &
         + w(r85f) &
         + w(r87f) &
         + w(r88f) &
         + w(r91f) &
         - w(r117f) &
         + w(r123f) &
         - w(r128f) &
         + w(r128b) &
         + w(r129f) &
         - w(r129b) &
         + w(r136f) &
         - w(r136b) &
         + w(r146f) &
         + w(r160f) &
         + w(r161f) &
         + 2.0_pr * w(r162f) &
         + w(r163f) &
         + w(r166f) &
         + 2.0_pr * w(r171f) &
         - w(r175f) &
         + w(r175b) &
         + w(r182f) &
         - w(r182b) &
         - w(r196f) &
         + w(r196b) &
         - w(r202f) &
         + w(r202b) 

    cdot(sC3H3) = 0.0_pr &
         + w(r83f) &
         - w(r83b) &
         + w(r117f) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r121b) &
         - w(r122f) &
         - w(r123f) &
         - w(r124f) &
         - w(r125f) &
         -2.0_pr * w(r126f) &
         + 2.0_pr * w(r126b) &
         -2.0_pr * w(r127f) &
         + 2.0_pr * w(r127b) &
         - w(r128f) &
         + w(r128b) &
         + w(r130f) &
         - w(r130b) &
         + w(r131f) &
         + w(r132f) &
         + w(r139f) &
         + w(r140f) &
         + w(r141f) &
         - w(r147f) &
         - w(r186f) 

    cdot(sAC3H5) = 0.0_pr &
         + w(r92f) &
         - w(r92b) &
         + w(r142f) &
         - w(r142b) &
         - w(r143f) &
         - w(r144f) &
         - w(r145f) &
         - w(r146f) &
         - w(r147f) &
         - w(r150f) &
         + w(r150b) &
         + w(r154f) &
         + w(r155f) &
         + w(r156f) &
         + w(r157f) &
         + w(r166f) 

    cdot(sC2H6) = 0.0_pr &
         + w(r102f) &
         - w(r102b) &
         + w(r105f) &
         + w(r106f) &
         - w(r108f) &
         + w(r108b) &
         - w(r109f) &
         - w(r110f) &
         - w(r111f) &
         - w(r112f) 

    cdot(sPC3H4) = 0.0_pr &
         + w(r120f) &
         - w(r120b) &
         - w(r129f) &
         + w(r129b) &
         - w(r130f) &
         + w(r130b) &
         - w(r131f) &
         - w(r132f) &
         - w(r133f) &
         - w(r134f) &
         - w(r135f) &
         + w(r137f) &
         - w(r137b) &
         + w(r138f) &
         - w(r138b) 

    cdot(sAC3H4) = 0.0_pr &
         + w(r121f) &
         - w(r121b) &
         - w(r136f) &
         + w(r136b) &
         - w(r137f) &
         + w(r137b) &
         - w(r138f) &
         + w(r138b) &
         - w(r139f) &
         - w(r140f) &
         - w(r141f) &
         - w(r142f) &
         + w(r142b) &
         + w(r143f) &
         + w(r144f) &
         + w(r145f) 

    cdot(sC6H6) = 0.0_pr &
         + w(r126f) &
         - w(r126b) &
         + w(r147f) &
         - w(r176f) &
         + w(r176b) &
         - w(r177f) &
         - w(r178f) &
         - w(r179f) &
         - w(r180f) &
         - w(r181f) &
         + w(r187f) &
         - w(r187b) 

    cdot(sC5H5) = 0.0_pr &
         + w(r128f) &
         - w(r128b) &
         - w(r158f) &
         -2.0_pr * w(r159f) &
         - w(r160f) &
         - w(r161f) &
         - w(r162f) &
         - w(r163f) &
         + w(r164f) &
         - w(r164b) &
         + w(r165f) &
         + w(r167f) &
         + w(r168f) &
         + w(r169f) &
         + w(r170f) &
         + w(r172f) &
         + w(r174f) &
         + w(r180f) &
         + w(r182f) &
         - w(r182b) 

    cdot(sC3H6) = 0.0_pr &
         - w(r148f) &
         + w(r148b) &
         - w(r149f) &
         + w(r149b) &
         + w(r150f) &
         - w(r150b) &
         - w(r151f) &
         + w(r151b) &
         - w(r152f) &
         - w(r153f) &
         - w(r154f) &
         - w(r155f) &
         - w(r156f) &
         - w(r157f) 

    cdot(sC5H6) = 0.0_pr &
         + w(r158f) &
         - w(r164f) &
         + w(r164b) &
         - w(r165f) &
         - w(r166f) &
         - w(r167f) &
         - w(r168f) &
         - w(r169f) &
         + w(r173f) &
         + w(r181f) 

    cdot(sC10H8) = 0.0_pr &
         + w(r159f) &
         + w(r186f) &
         + w(r197f) &
         - w(r197b) &
         + w(r202f) &
         - w(r202b) &
         - w(r203f) &
         + w(r203b) &
         - w(r204f) &
         + w(r204b) &
         - w(r205f) &
         + w(r205b) 

    cdot(sC7H8) = 0.0_pr &
         - w(r187f) &
         + w(r187b) &
         - w(r188f) &
         + w(r188b) &
         - w(r189f) &
         + w(r189b) &
         - w(r190f) &
         - w(r191f) &
         - w(r192f) 

    cdot(sC8H6) = 0.0_pr &
         - w(r198f) &
         + w(r198b) &
         - w(r199f) &
         + w(r199b) &
         - w(r200f) &
         + w(r200b) &
         + w(r201f) &
         - w(r201b) 

  
    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) :: SCH2_ct
    real(pr) :: SCH2_num
    real(pr) :: SCH2_denom
    real(pr) :: SCH2_TCH2
    real(pr) :: SCH2_HCCO
    real(pr) :: TCH2_ct
    real(pr) :: TCH2_num
    real(pr) :: TCH2_denom
    real(pr) :: TCH2_SCH2
    real(pr) :: TCH2_CH
    real(pr) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_TCH2
    real(pr) :: CH_C2H
    real(pr) :: CH_HCCO
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_TCH2
    real(pr) :: HCO_CH
    real(pr) :: HCO_C2H3
    real(pr) :: HCO_C2H
    real(pr) :: HCO_HCCO
    real(pr) :: HCO_C7H6O
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_CH
    real(pr) :: C2H_ct
    real(pr) :: C2H_num
    real(pr) :: C2H_denom
    real(pr) :: HCCO_ct
    real(pr) :: HCCO_num
    real(pr) :: HCCO_denom
    real(pr) :: HCCO_C2H
    real(pr) :: NC3H7_ct
    real(pr) :: NC3H7_num
    real(pr) :: NC3H7_denom
    real(pr) :: C6H5_ct
    real(pr) :: C6H5_num
    real(pr) :: C6H5_denom
    real(pr) :: C6H5_C8H7
    real(pr) :: C6H5_C7H7
    real(pr) :: C6H5_C7H6O
    real(pr) :: C8H7_ct
    real(pr) :: C8H7_num
    real(pr) :: C8H7_denom
    real(pr) :: C8H7_C6H5
    real(pr) :: C7H7_ct
    real(pr) :: C7H7_num
    real(pr) :: C7H7_denom
    real(pr) :: C7H7_C6H5
    real(pr) :: C7H6O_ct
    real(pr) :: C7H6O_num
    real(pr) :: C7H6O_denom
    real(pr) :: C7H6O_C7H7
    real(pr) :: C8H5_ct
    real(pr) :: C8H5_num
    real(pr) :: C8H5_denom
    real(pr) :: C8H5_C10H7
    real(pr) :: C10H7_ct
    real(pr) :: C10H7_num
    real(pr) :: C10H7_denom
    real(pr) :: C10H7_C8H5
  
    C7H7_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r182f)&
             + k(r183f)* c(sH) &
             + k(r184f)* c(sO) &
             + k(r185f)* c(sHO2) &
             + k(r186f)* c(sC3H3) &
             + k(r188b)* c(sH)   )

    C7H7_num = ( 0.0_pr &
             + k(r188f)* c(sC7H8) &
             + k(r190f)* c(sH) * c(sC7H8) &
             + k(r191f)* c(sC7H8) * c(sOH) &
             + k(r192f)* c(sC7H8) * c(sCH3) &
             + k(r182b)* c(sC2H2) * c(sC5H5)   )

    C7H7_ct = C7H7_num / C7H7_denom

    C7H7_C6H5 = - ( 0.0_pr &
          + k(r183f) * c(sCH3)  ) / C7H7_denom

    C6H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r170f)* c(sO2) &
             + k(r171f)* c(sO2) &
             + k(r172f)* c(sO) &
             + k(r173f)* c(sOH) &
             + k(r174f)* c(sHO2) &
             + k(r175f)* c(sC2H2) &
             + k(r127b)* c(sH) &
             + k(r176b)* c(sH) &
             + k(r183b)* c(sCH3) &
             + k(r189b)* c(sCH3)   )

    C6H5_num = ( 0.0_pr &
             + k(r127f) *c(sC3H3)** 2.0_pr&
             + k(r176f)* c(sC6H6) &
             + k(r177f)* c(sH) * c(sC6H6) &
             + k(r178f)* c(sC6H6) * c(sOH) &
             + k(r179f)* c(sCH3) * c(sC6H6) &
             + k(r189f)* c(sC7H8)   )

    C6H5_ct = C6H5_num / C6H5_denom

    C6H5_C8H7 = - ( 0.0_pr &
          + k(r175f)  ) / C6H5_denom

    C6H5_C7H7 = - ( 0.0_pr &
          + k(r183f) * c(sH)  ) / C6H5_denom

    C6H5_C7H6O = - ( 0.0_pr &
          + k(r193f) &
          + k(r194f) * c(sH) &
          + k(r195f) * c(sCH3)  ) / C6H5_denom

    C7H6O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r193f)&
             + k(r194f)* c(sH) &
             + k(r195f)* c(sCH3)   )

    C7H6O_num = ( 0.0_pr   )

    C7H6O_ct = C7H6O_num / C7H6O_denom

    C7H6O_C7H7 = - ( 0.0_pr &
          + k(r184f) * c(sO) &
          + k(r185f) * c(sHO2)  ) / C7H6O_denom

    C8H7_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r201f)&
             + k(r202f)* c(sC2H2) &
             + k(r175b)  )

    C8H7_num = ( 0.0_pr &
             + k(r201b)* c(sH) * c(sC8H6) &
             + k(r202b)* c(sH) * c(sC10H8)   )

    C8H7_ct = C8H7_num / C8H7_denom

    C8H7_C6H5 = - ( 0.0_pr &
          + k(r175f) * c(sC2H2)  ) / C8H7_denom

    cqss(sqssC7H7) = ( & 
             (C7H7_ct) - (& 
             (& 
             (C6H5_ct) - (C8H7_ct) * (C6H5_C8H7)) - (C7H6O_ct) * (C6H5_C7H6O)) * (C7H7_C6H5) / (& 
                (1.0_pr) - (C8H7_C6H5) * (C6H5_C8H7)) ) / ( & 
                (1.0_pr) - (& 
                (C6H5_C7H7) - (C7H6O_C7H7) * (C6H5_C7H6O)) * (C7H7_C6H5) / (& 
                (1.0_pr) - (C8H7_C6H5) * (C6H5_C8H7)) )  

    cqss(sqssC6H5) = (& 
             (& 
             (C6H5_ct) - (C8H7_ct) * (C6H5_C8H7)) - (C7H6O_ct) * (C6H5_C7H6O)& 
              - (& 
                (C6H5_C7H7) - (C7H6O_C7H7) * (C6H5_C7H6O)) * cqss(sqssC7H7)) / (& 
                (1.0_pr) - (C8H7_C6H5) * (C6H5_C8H7))  

    cqss(sqssC7H6O) = C7H6O_ct& 
              - (C7H6O_C7H7) * cqss(sqssC7H7)  

    cqss(sqssC8H7) = C8H7_ct& 
              - (C8H7_C6H5) * cqss(sqssC6H5)  

    NC3H7_ct = 0.0_pr

    cqss(sqssNC3H7) = NC3H7_ct  

    C10H7_ct = 0.0_pr

    C8H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r196f)* c(sC2H2) &
             + k(r197f)* c(sC2H4) &
             + k(r198b)* c(sH) &
             + k(r199b)* c(sH2) &
             + k(r200b)* c(sH2O)   )

    C8H5_num = ( 0.0_pr &
             + k(r198f)* c(sC8H6) &
             + k(r199f)* c(sH) * c(sC8H6) &
             + k(r200f)* c(sC8H6) * c(sOH) &
             + k(r197b)* c(sH) * c(sC10H8)   )

    C8H5_ct = C8H5_num / C8H5_denom

    C8H5_C10H7 = - ( 0.0_pr &
          + k(r196f)  ) / C8H5_denom

    cqss(sqssC10H7) = ( & 
             (C10H7_ct) - (C8H5_ct) * (C10H7_C8H5) ) / ( & 
                (1.0_pr) - (C8H5_C10H7) * (C10H7_C8H5) )  

    cqss(sqssC8H5) = C8H5_ct& 
              - (C8H5_C10H7) * cqss(sqssC10H7)  

    C2H_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r72f)* c(sO) &
             + k(r73f)* c(sOH) &
             + k(r74f)* c(sO2) &
             + k(r75f)* c(sH2)   )

    C2H_num = ( 0.0_pr &
             + k(r79f)* c(sC2H2) * c(sO) &
             + k(r80f)* c(sC2H2) * c(sOH) &
             + k(r75b)* c(sH) * c(sC2H2)   )

    C2H_ct = C2H_num / C2H_denom

    cqss(sqssC2H) = C2H_ct  

    HCCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r113f)&
             + k(r114f)* c(sH) &
             + k(r115f)* c(sO) &
             + k(r116f)* c(sO2) &
             + k(r117f)* c(sC2H2) &
             + k(r118f)* c(sCH3) &
             + k(r119f)* c(sOH)   )

    HCCO_num = ( 0.0_pr &
             + k(r73f)* cqss(sqssC2H) * c(sOH) &
             + k(r77f)* c(sC2H2) * c(sO) &
             + k(r133f)* c(sPC3H4) * c(sO)   )

    HCCO_ct = HCCO_num / HCCO_denom

    cqss(sqssHCCO) = HCCO_ct  

    SCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r1f)* c(sN2) &
             + k(r27f)* c(sH2) &
             + k(r28f)* c(sO2) &
             + k(r29f)* c(sO2) &
             + k(r30f)* c(sH2O) &
             + k(r31f)* c(sH2O) &
             + k(r32f)* c(sH2O) &
             + k(r56f)* c(sCO) &
             + k(r70f)* c(sCO2) &
             + k(r71f)* c(sCO2) &
             + k(r83f)* c(sC2H2) &
             + k(r38b)* c(sH2O)   )

    SCH2_num = ( 0.0_pr &
             + k(r38f)* c(sCH3) * c(sOH) &
             + k(r46f) *c(sCH3)** 2.0_pr&
             + k(r114f)* c(sH) * cqss(sqssHCCO) &
             + k(r27b)* c(sH) * c(sCH3) &
             + k(r83b)* c(sH) * c(sC3H3)   )

    SCH2_ct = SCH2_num / SCH2_denom

    SCH2_TCH2 = - ( 0.0_pr &
          + k(r1b) * c(sN2) &
          + k(r31b) * c(sH2O) &
          + k(r56b) * c(sCO) &
          + k(r70b) * c(sCO2)  ) / SCH2_denom

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r15f)* c(sO) &
             + k(r16f)* c(sOH) &
             + k(r17f)* c(sH2) &
             + k(r18f)* c(sH2O) &
             + k(r19f)* c(sO2) &
             + k(r43f)* c(sCH3) &
             + k(r50f)* c(sCH4) &
             + k(r69f)* c(sCO2) &
             + k(r22b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r72f)* cqss(sqssC2H) * c(sO) &
             + k(r113f)* cqss(sqssHCCO)   )

    CH_ct = CH_num / CH_denom

    CH_TCH2 = - ( 0.0_pr &
          + k(r22f) * c(sOH) &
          + k(r17b) * c(sH)  ) / CH_denom

    TCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r20f)* c(sO) &
             + k(r21f)* c(sOH) &
             + k(r22f)* c(sOH) &
             + k(r23f)* c(sH2) &
             + k(r24f)* c(sO2) &
             + k(r25f)* c(sO2) &
             + k(r26f)* c(sO2) &
             + k(r44f)* c(sCH3) &
             + k(r51f)* c(sCH4) &
             + k(r1b)* c(sN2) &
             + k(r17b)* c(sH) &
             + k(r31b)* c(sH2O) &
             + k(r37b)* c(sH2O) &
             + k(r56b)* c(sCO) &
             + k(r70b)* c(sCO2)   )

    TCH2_num = ( 0.0_pr &
             + k(r37f)* c(sCH3) * c(sOH) &
             + k(r78f)* c(sC2H2) * c(sO) &
             + k(r96f)* c(sC2H4) * c(sO) &
             + k(r23b)* c(sH) * c(sCH3) &
             + k(r51b) *c(sCH3)** 2.0_pr  )

    TCH2_ct = TCH2_num / TCH2_denom

    TCH2_SCH2 = - ( 0.0_pr &
          + k(r1f) * c(sN2) &
          + k(r31f) * c(sH2O) &
          + k(r56f) * c(sCO) &
          + k(r70f) * c(sCO2)  ) / TCH2_denom

    TCH2_CH = - ( 0.0_pr &
          + k(r17f) * c(sH2) &
          + k(r22f) * c(sH2O)  ) / TCH2_denom

    cqss(sqssSCH2) = ( & 
             (& 
             (SCH2_ct) - (TCH2_ct) * (SCH2_TCH2)) - (& 
             (CH_ct) - (TCH2_ct) * (CH_TCH2)) * (& 
                 - (TCH2_CH) * (SCH2_TCH2)) / (& 
                (1.0_pr) - (TCH2_CH) * (CH_TCH2)) ) / ( & 
                (1.0_pr) - (TCH2_SCH2) * (SCH2_TCH2) )  

    cqss(sqssCH) = (& 
             (CH_ct) - (TCH2_ct) * (CH_TCH2)& 
              - (& 
                 - (TCH2_SCH2) * (CH_TCH2)) * cqss(sqssSCH2)) / (& 
                (1.0_pr) - (TCH2_CH) * (CH_TCH2))  

    cqss(sqssTCH2) = TCH2_ct& 
              - (TCH2_SCH2) * cqss(sqssSCH2)& 
              - (TCH2_CH) * cqss(sqssCH)  

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r84f)* c(sH) &
             + k(r85f)* c(sH) &
             + k(r86f)* c(sO) &
             + k(r87f)* c(sOH) &
             + k(r88f)* c(sO2) &
             + k(r89f)* c(sO2) &
             + k(r90f)* c(sO2) &
             + k(r91f)* c(sCH3) &
             + k(r92f)* c(sCH3) &
             + k(r76b)&
             + k(r149b)* c(sCH3)   )

    C2H3_num = ( 0.0_pr &
             + k(r43f)* cqss(sqssCH) * c(sCH3) &
             + k(r76f)* c(sH) * c(sC2H2) &
             + k(r94f)* c(sH) * c(sC2H4) &
             + k(r98f)* c(sC2H4) * c(sOH) &
             + k(r100f)* c(sCH3) * c(sC2H4) &
             + k(r125f)* c(sHO2) * c(sC3H3) &
             + k(r149f)* c(sC3H6) &
             + k(r160f)* c(sC5H5) * c(sO) &
             + k(r161f)* c(sHO2) * c(sC5H5) &
             + k(r163f)* c(sC5H5) * c(sOH) &
             + k(r92b)* c(sH) * c(sAC3H5)   )

    C2H3_ct = C2H3_num / C2H3_denom

    cqss(sqssC2H3) = C2H3_ct  

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r57f)* c(sH) &
             + k(r58f)* c(sO) &
             + k(r59f)* c(sO) &
             + k(r60f)* c(sOH) &
             + k(r61f)* m(mM61) &
             + k(r62f)* c(sH2O) &
             + k(r63f)* c(sO2) &
             + k(r64f)* c(sCH3) &
             + k(r105f)* c(sC2H5) &
             + k(r65b)* c(sH2)   )

    HCO_num = ( 0.0_pr &
             + k(r16f)* cqss(sqssCH) * c(sOH) &
             + k(r19f)* cqss(sqssCH) * c(sO2) &
             + k(r20f)* cqss(sqssTCH2) * c(sO) &
             + k(r65f)* c(sCH2O) * c(sH) &
             + k(r66f)* c(sCH2O) * c(sO) &
             + k(r67f)* c(sCH2O) * c(sOH) &
             + k(r68f)* c(sCH2O) * c(sCH3) &
             + k(r69f)* cqss(sqssCH) * c(sCO2) &
             + k(r74f)* cqss(sqssC2H) * c(sO2) &
             + k(r90f)* cqss(sqssC2H3) * c(sO2) &
             + k(r97f)* c(sC2H4) * c(sO) &
     +2.0_pr * k(r119f)* cqss(sqssHCCO) * c(sOH) &
             + k(r153f)* c(sC3H6) * c(sO) &
             + k(r193f)* cqss(sqssC7H6O) &
             + k(r61b)* c(sH) * c(sCO) * m(mM61) &
             + k(r62b)* c(sH) * c(sH2O) * c(sCO)   )

    HCO_ct = HCO_num / HCO_denom

    cqss(sqssHCO) = HCO_ct  

    do index = 1, nqss
      cqss(index) = max(cqss(index), 1e-60_pr)
    end do
    
    return
  end subroutine get_qss  
  
    
  ! ----------------------------------------------- !
  ! Mass fractions to concentrations                !
  ! ----------------------------------------------- !
  subroutine y2c(y, W_sp, rho, c)
    implicit none

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

    integer :: k
    
    ! Conversion
    c = y * rho / W_sp

    ! Concentrations clipping
    do k =1, nspec
      c(k) = max(c(k), 1e-60_pr)
    end do

        return
    end subroutine y2c

end module mod_customkinetics

! ----------------------------------------------- !
! Cantera routine                                 !
! ----------------------------------------------- !
subroutine customkinetics(P, T, rho, 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, rho, 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)
      
! Convert to kmol
  wdot = wdot / 1000.0_pr

  return
end subroutine customkinetics