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

!--------------------------------------------------------------------------------------------------
!     MODULE mod_customkinetics
!>    @details Generated by ARCANE custom kinetics routine to compute the chemical source terms.
!!    @authors Q. Cazères, P. Pepiot, J. Wirtz
!!    @date    24/01/2019
!--------------------------------------------------------------------------------------------------
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 = 21
  integer, parameter :: nqss = 14
  integer, parameter :: nreac = 161
  integer, parameter :: nreac_reverse = 161
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 11
  integer, parameter :: nFO = 16
  integer, parameter :: nTB_reverse = 11
  integer, parameter :: nFO_reverse = 16
  
  ! Index of species
  integer, parameter :: sN2 = 1
  integer, parameter :: sH = 2
  integer, parameter :: sO2 = 3
  integer, parameter :: sOH = 4
  integer, parameter :: sO = 5
  integer, parameter :: sH2 = 6
  integer, parameter :: sH2O = 7
  integer, parameter :: sHO2 = 8
  integer, parameter :: sH2O2 = 9
  integer, parameter :: sCO = 10
  integer, parameter :: sCO2 = 11
  integer, parameter :: sCH3 = 12
  integer, parameter :: sCH4 = 13
  integer, parameter :: sCH2O = 14
  integer, parameter :: sC2H4 = 15
  integer, parameter :: sC2H6 = 16
  integer, parameter :: sC2H2 = 17
  integer, parameter :: sCH2CO = 18
  integer, parameter :: sC3H6 = 19
  integer, parameter :: sC3H8 = 20
  integer, parameter :: sOC3H5OOH = 21

  integer, parameter :: sqssHCO = 1
  integer, parameter :: sqssTXCH2 = 2
  integer, parameter :: sqssSXCH2 = 3
  integer, parameter :: sqssCH3O = 4
  integer, parameter :: sqssC2H5 = 5
  integer, parameter :: sqssCH = 6
  integer, parameter :: sqssC2H3 = 7
  integer, parameter :: sqssCH2CHO = 8
  integer, parameter :: sqssHCCO = 9
  integer, parameter :: sqssCH2OH = 10
  integer, parameter :: sqssC3H5 = 11
  integer, parameter :: sqssIXC3H7 = 12
  integer, parameter :: sqssNXC3H7 = 13
  integer, parameter :: sqssC3H6OOH = 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 :: r1b = 162
  integer, parameter :: r2b = 163
  integer, parameter :: r3b = 164
  integer, parameter :: r4b = 165
  integer, parameter :: r5b = 166
  integer, parameter :: r6b = 167
  integer, parameter :: r7b = 168
  integer, parameter :: r8b = 169
  integer, parameter :: r9b = 170
  integer, parameter :: r10b = 171
  integer, parameter :: r11b = 172
  integer, parameter :: r12b = 173
  integer, parameter :: r13b = 174
  integer, parameter :: r14b = 175
  integer, parameter :: r15b = 176
  integer, parameter :: r16b = 177
  integer, parameter :: r17b = 178
  integer, parameter :: r18b = 179
  integer, parameter :: r19b = 180
  integer, parameter :: r20b = 181
  integer, parameter :: r21b = 182
  integer, parameter :: r22b = 183
  integer, parameter :: r23b = 184
  integer, parameter :: r24b = 185
  integer, parameter :: r25b = 186
  integer, parameter :: r26b = 187
  integer, parameter :: r27b = 188
  integer, parameter :: r28b = 189
  integer, parameter :: r29b = 190
  integer, parameter :: r30b = 191
  integer, parameter :: r31b = 192
  integer, parameter :: r32b = 193
  integer, parameter :: r33b = 194
  integer, parameter :: r34b = 195
  integer, parameter :: r35b = 196
  integer, parameter :: r36b = 197
  integer, parameter :: r37b = 198
  integer, parameter :: r38b = 199
  integer, parameter :: r39b = 200
  integer, parameter :: r40b = 201
  integer, parameter :: r41b = 202
  integer, parameter :: r42b = 203
  integer, parameter :: r43b = 204
  integer, parameter :: r44b = 205
  integer, parameter :: r45b = 206
  integer, parameter :: r46b = 207
  integer, parameter :: r47b = 208
  integer, parameter :: r48b = 209
  integer, parameter :: r49b = 210
  integer, parameter :: r50b = 211
  integer, parameter :: r51b = 212
  integer, parameter :: r52b = 213
  integer, parameter :: r53b = 214
  integer, parameter :: r54b = 215
  integer, parameter :: r55b = 216
  integer, parameter :: r56b = 217
  integer, parameter :: r57b = 218
  integer, parameter :: r58b = 219
  integer, parameter :: r59b = 220
  integer, parameter :: r60b = 221
  integer, parameter :: r61b = 222
  integer, parameter :: r62b = 223
  integer, parameter :: r63b = 224
  integer, parameter :: r64b = 225
  integer, parameter :: r65b = 226
  integer, parameter :: r66b = 227
  integer, parameter :: r67b = 228
  integer, parameter :: r68b = 229
  integer, parameter :: r69b = 230
  integer, parameter :: r70b = 231
  integer, parameter :: r71b = 232
  integer, parameter :: r72b = 233
  integer, parameter :: r73b = 234
  integer, parameter :: r74b = 235
  integer, parameter :: r75b = 236
  integer, parameter :: r76b = 237
  integer, parameter :: r77b = 238
  integer, parameter :: r78b = 239
  integer, parameter :: r79b = 240
  integer, parameter :: r80b = 241
  integer, parameter :: r81b = 242
  integer, parameter :: r82b = 243
  integer, parameter :: r83b = 244
  integer, parameter :: r84b = 245
  integer, parameter :: r85b = 246
  integer, parameter :: r86b = 247
  integer, parameter :: r87b = 248
  integer, parameter :: r88b = 249
  integer, parameter :: r89b = 250
  integer, parameter :: r90b = 251
  integer, parameter :: r91b = 252
  integer, parameter :: r92b = 253
  integer, parameter :: r93b = 254
  integer, parameter :: r94b = 255
  integer, parameter :: r95b = 256
  integer, parameter :: r96b = 257
  integer, parameter :: r97b = 258
  integer, parameter :: r98b = 259
  integer, parameter :: r99b = 260
  integer, parameter :: r100b = 261
  integer, parameter :: r101b = 262
  integer, parameter :: r102b = 263
  integer, parameter :: r103b = 264
  integer, parameter :: r104b = 265
  integer, parameter :: r105b = 266
  integer, parameter :: r106b = 267
  integer, parameter :: r107b = 268
  integer, parameter :: r108b = 269
  integer, parameter :: r109b = 270
  integer, parameter :: r110b = 271
  integer, parameter :: r111b = 272
  integer, parameter :: r112b = 273
  integer, parameter :: r113b = 274
  integer, parameter :: r114b = 275
  integer, parameter :: r115b = 276
  integer, parameter :: r116b = 277
  integer, parameter :: r117b = 278
  integer, parameter :: r118b = 279
  integer, parameter :: r119b = 280
  integer, parameter :: r120b = 281
  integer, parameter :: r121b = 282
  integer, parameter :: r122b = 283
  integer, parameter :: r123b = 284
  integer, parameter :: r124b = 285
  integer, parameter :: r125b = 286
  integer, parameter :: r126b = 287
  integer, parameter :: r127b = 288
  integer, parameter :: r128b = 289
  integer, parameter :: r129b = 290
  integer, parameter :: r130b = 291
  integer, parameter :: r131b = 292
  integer, parameter :: r132b = 293
  integer, parameter :: r133b = 294
  integer, parameter :: r134b = 295
  integer, parameter :: r135b = 296
  integer, parameter :: r136b = 297
  integer, parameter :: r137b = 298
  integer, parameter :: r138b = 299
  integer, parameter :: r139b = 300
  integer, parameter :: r140b = 301
  integer, parameter :: r141b = 302
  integer, parameter :: r142b = 303
  integer, parameter :: r143b = 304
  integer, parameter :: r144b = 305
  integer, parameter :: r145b = 306
  integer, parameter :: r146b = 307
  integer, parameter :: r147b = 308
  integer, parameter :: r148b = 309
  integer, parameter :: r149b = 310
  integer, parameter :: r150b = 311
  integer, parameter :: r151b = 312
  integer, parameter :: r152b = 313
  integer, parameter :: r153b = 314
  integer, parameter :: r154b = 315
  integer, parameter :: r155b = 316
  integer, parameter :: r156b = 317
  integer, parameter :: r157b = 318
  integer, parameter :: r158b = 319
  integer, parameter :: r159b = 320
  integer, parameter :: r160b = 321
  integer, parameter :: r161b = 322
  
  ! Index of third body species
  integer, parameter :: mM5 = 1
  integer, parameter :: mM6 = 2
  integer, parameter :: mM7 = 3
  integer, parameter :: mM8 = 4
  integer, parameter :: mM28 = 5
  integer, parameter :: mM61 = 6
  integer, parameter :: mM78 = 7
  integer, parameter :: mM95 = 8
  integer, parameter :: mM96 = 9
  integer, parameter :: mM119 = 10
  integer, parameter :: mM120 = 11
  
  integer, parameter :: mM9 = 12
  integer, parameter :: mM16 = 13
  integer, parameter :: mM24 = 14
  integer, parameter :: mM35 = 15
  integer, parameter :: mM56 = 16
  integer, parameter :: mM57 = 17
  integer, parameter :: mM83 = 18
  integer, parameter :: mM89 = 19
  integer, parameter :: mM98 = 20
  integer, parameter :: mM131 = 21
  integer, parameter :: mM135 = 22
  integer, parameter :: mM138 = 23
  integer, parameter :: mM141 = 24
  integer, parameter :: mM153 = 25
  integer, parameter :: mM155 = 26
  integer, parameter :: mM156 = 27
  
  ! Index of third body reactions
  integer, parameter :: TBr5f = 1
  integer, parameter :: TBr6f = 2
  integer, parameter :: TBr7f = 3
  integer, parameter :: TBr8f = 4
  integer, parameter :: TBr28f = 5
  integer, parameter :: TBr61f = 6
  integer, parameter :: TBr78f = 7
  integer, parameter :: TBr95f = 8
  integer, parameter :: TBr96f = 9
  integer, parameter :: TBr119f = 10
  integer, parameter :: TBr120f = 11
  integer, parameter :: TBr5b = 12
  integer, parameter :: TBr6b = 13
  integer, parameter :: TBr7b = 14
  integer, parameter :: TBr8b = 15
  integer, parameter :: TBr28b = 16
  integer, parameter :: TBr61b = 17
  integer, parameter :: TBr78b = 18
  integer, parameter :: TBr95b = 19
  integer, parameter :: TBr96b = 20
  integer, parameter :: TBr119b = 21
  integer, parameter :: TBr120b = 22
  
  ! Index of fall off reactions
  integer, parameter :: FOr9f = 1
  integer, parameter :: FOr16f = 2
  integer, parameter :: FOr24f = 3
  integer, parameter :: FOr35f = 4
  integer, parameter :: FOr56f = 5
  integer, parameter :: FOr57f = 6
  integer, parameter :: FOr83f = 7
  integer, parameter :: FOr89f = 8
  integer, parameter :: FOr98f = 9
  integer, parameter :: FOr131f = 10
  integer, parameter :: FOr135f = 11
  integer, parameter :: FOr138f = 12
  integer, parameter :: FOr141f = 13
  integer, parameter :: FOr153f = 14
  integer, parameter :: FOr155f = 15
  integer, parameter :: FOr156f = 16
  integer, parameter :: FOr9b = 17
  integer, parameter :: FOr16b = 18
  integer, parameter :: FOr24b = 19
  integer, parameter :: FOr35b = 20
  integer, parameter :: FOr56b = 21
  integer, parameter :: FOr57b = 22
  integer, parameter :: FOr83b = 23
  integer, parameter :: FOr89b = 24
  integer, parameter :: FOr98b = 25
  integer, parameter :: FOr131b = 26
  integer, parameter :: FOr135b = 27
  integer, parameter :: FOr138b = 28
  integer, parameter :: FOr141b = 29
  integer, parameter :: FOr153b = 30
  integer, parameter :: FOr155b = 31
  integer, parameter :: FOr156b = 32
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.02801348_pr, & ! N2
       0.00100794_pr, & ! H
       0.0319988_pr, & ! O2
       0.01700734_pr, & ! OH
       0.0159994_pr, & ! O
       0.00201588_pr, & ! H2
       0.01801528_pr, & ! H2O
       0.03300674_pr, & ! HO2
       0.03401468_pr, & ! H2O2
       0.0280104_pr, & ! CO
       0.0440098_pr, & ! CO2
       0.01503482_pr, & ! CH3
       0.01604276_pr, & ! CH4
       0.03002628_pr, & ! CH2O
       0.02805376_pr, & ! C2H4
       0.03006964_pr, & ! C2H6
       0.02603788_pr, & ! C2H2
       0.04203728_pr, & ! CH2CO
       0.04208064_pr, & ! C3H6
       0.04409652_pr, & ! C3H8
       0.09007884_pr & ! OC3H5OOH
      !0.02901834_pr, & ! HCO
      !0.01402688_pr, & ! TXCH2
      !0.01402688_pr, & ! SXCH2
      !0.03103422_pr, & ! CH3O
      !0.0290617_pr, & ! C2H5
      !0.01301894_pr, & ! CH
      !0.02704582_pr, & ! C2H3
      !0.04304522_pr, & ! CH2CHO
      !0.04102934_pr, & ! HCCO
      !0.03103422_pr, & ! CH2OH
      !0.0410727_pr, & ! C3H5
      !0.04308858_pr, & ! IXC3H7
      !0.04308858_pr, & ! NXC3H7
      !0.07508738_pr & ! C3H6OOH
       /)
  
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(mM5) = (0.9_pr)*c(sCO) &
         + (2.8_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM6) = (0.9_pr)*c(sCO) &
         + (2.8_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM7) = (0.9_pr)*c(sCO) &
         + (2.8_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM8) = (0.9_pr)*c(sCO) &
         + (2.8_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM28) = (1.5_pr)*c(sCO) &
         + (1.5_pr)*c(sCO2) &
         + (0.9_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM61) = (0.8_pr)*c(sCO) &
         + (2.6_pr)*c(sCO2) &
         + (1.4_pr)*c(sH2) &
         + (14.4_pr)*c(sH2O) &
         + sum(c)

    M(mM78) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM95) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM96) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM119) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM120) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM9) = (0.5_pr)*c(sC2H6) &
         + (0.2_pr)*c(sCO) &
         + (1.4_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (15_pr)*c(sH2O) &
         + sum(c)

    M(mM16) = (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + (5_pr)*c(sH2O2) &
         + sum(c)

    M(mM24) = (1_pr)*c(sCO) &
         + (3_pr)*c(sCO2) &
         + (1.5_pr)*c(sH2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

    M(mM35) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM56) = (3_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (15_pr)*c(sH2O) &
         + sum(c)

    M(mM57) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM83) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM89) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM98) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM131) = sum(c)

    M(mM135) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM138) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM141) = sum(c)

    M(mM153) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

    M(mM155) = sum(c)

    M(mM156) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

  end subroutine get_thirdbodies
  
  ! ----------------------------------------------- !
  ! Evaluate rate coefficients                      !
  ! ----------------------------------------------- !  
  subroutine get_rate_coefficients(k,M,Tloc,Ploc)
    implicit none
      
    real(pr), dimension(nreac + nreac_reverse) :: k
    real(pr), dimension(nFO + nFO_reverse) :: k_0
    real(pr), dimension(nFO + nFO_reverse) :: k_inf
    real(pr), dimension(nFO + nFO_reverse) :: FC
    real(pr), dimension(nTB + nFO) :: M
    real(pr) :: Tloc,Ploc,R_T_inv,T_log

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

    k(r1f) = (35200000000_pr)*exp((-71420.00136_pr)*R_T_inv + T_log *(-0.7_pr))
    k(r2f) = (0.0506_pr)*exp((-26319.99592_pr)*R_T_inv + T_log *(2.67_pr))
    k(r3f) = (1170_pr)*exp((-15210.01152_pr)*R_T_inv + T_log *(1.3_pr))
    k(r4f) = (0.7_pr)*exp((-60870.00352_pr)*R_T_inv + T_log *(2.33_pr))
    k(r5f) = (1300000_pr)*exp(T_log *(-1_pr))
    k(r6f) = (40000000000_pr)*exp(T_log *(-2_pr))
    k(r7f) = (6170_pr)*exp(T_log *(-0.5_pr))
    k(r8f) = (4710000_pr)*exp(T_log *(-1_pr))
    k_0(FOr9f) = (57500000_pr)*exp(T_log *(-1.4_pr))
    k_inf(FOr9f) = (4650000_pr)*exp(T_log *(0.44_pr))
    FC(FOr9f) = ((1.0_pr - 0.5_pr)*exp(-Tloc/(0_pr))) + (0.5_pr*exp(-Tloc/(1e+30_pr)))
    k(r9f) = getlindratecoeff(Tloc,k_0(FOr9f),k_inf(FOr9f),FC(FOr9f), M(mM9),Ploc)
    k(r10f) = (70800000_pr)*exp((-1233.98712_pr)*R_T_inv)
    k(r11f) = (16600000_pr)*exp((-3443.0136_pr)*R_T_inv)
    k(r12f) = (31000000_pr)*exp((-7199.99456_pr)*R_T_inv)
    k(r13f) = (20000000_pr)
    k(r14f) = (7000000_pr)*exp((4580.0156_pr)*R_T_inv)
    k(r15f) = (450000000_pr)*exp((-45729.99032_pr)*R_T_inv)
    k_0(FOr16f) = (27600000000000_pr)*exp(T_log *(-3.2_pr))
    k_inf(FOr16f) = (95500000_pr)*exp(T_log *(-0.27_pr))
    FC(FOr16f) = ((1.0_pr - 0.57_pr)*exp(-Tloc/(1e+30_pr))) + (0.57_pr*exp(-Tloc/(0_pr)))
    k(r16f) = getlindratecoeff(Tloc,k_0(FOr16f),k_inf(FOr16f),FC(FOr16f), M(mM16),Ploc)
    k(r17f) = (103000000_pr)*exp((-46200.02088_pr)*R_T_inv)
    k(r18f) = (194000_pr)*exp((5895.00496_pr)*R_T_inv)
    k(r19f) = (23000000_pr)*exp((-33263.0092_pr)*R_T_inv)
    k(r20f) = (10000000_pr)*exp((-15000.01656_pr)*R_T_inv)
    k(r21f) = (1740000_pr)*exp((-5999.98152_pr)*R_T_inv)
    k(r22f) = (75900000_pr)*exp((-30429.98096_pr)*R_T_inv)
    k(r23f) = (9.63_pr)*exp((-16700.0176_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr24f) = (1550000000000_pr)*exp((-17535.01848_pr)*R_T_inv + T_log *(-2.79_pr))
    k_inf(FOr24f) = (180000_pr)*exp((-9974.99072_pr)*R_T_inv)
    FC(FOr24f) = ((1.0_pr - 1_pr)*exp(-Tloc/(1_pr))) + (1_pr*exp(-Tloc/(10000000_pr)))&
              + exp(-(10000000_pr)/Tloc)
    k(r24f) = getlindratecoeff(Tloc,k_0(FOr24f),k_inf(FOr24f),FC(FOr24f), M(mM24),Ploc)
    k(r25f) = (4.4_pr)*exp((3100.00928_pr)*R_T_inv + T_log *(1.5_pr))
    k(r26f) = (20000000_pr)*exp((-95999.9972_pr)*R_T_inv)
    k(r27f) = (1000000_pr)*exp((-199577.0092_pr)*R_T_inv)
    k(r28f) = (186000000000_pr)*exp((-71130.00832_pr)*R_T_inv + T_log *(-1_pr))
    k(r29f) = (50000000_pr)
    k(r30f) = (30000000_pr)
    k(r31f) = (30000000_pr)
    k(r32f) = (30000000_pr)
    k(r33f) = (7580000_pr)*exp((-1714.97976_pr)*R_T_inv)
    k(r34f) = (50000000_pr)
    k_0(FOr35f) = (1350000000000_pr)*exp((-1777.9908_pr)*R_T_inv + T_log *(-2.57_pr))
    k_inf(FOr35f) = (1090000_pr)*exp((1088.00736_pr)*R_T_inv + T_log *(0.48_pr))
    FC(FOr35f) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r35f) = getlindratecoeff(Tloc,k_0(FOr35f),k_inf(FOr35f),FC(FOr35f), M(mM35),Ploc)
    k(r36f) = (57.4_pr)*exp((-11500.01688_pr)*R_T_inv + T_log *(1.9_pr))
    k(r37f) = (35000000_pr)*exp((-14699.98192_pr)*R_T_inv)
    k(r38f) = (39000_pr)*exp((-1700.00104_pr)*R_T_inv + T_log *(0.89_pr))
    k(r39f) = (60000000_pr)*exp((-170180.016_pr)*R_T_inv)
    k(r40f) = (0.0411_pr)*exp((-42720.02072_pr)*R_T_inv + T_log *(2.5_pr))
    k(r41f) = (0.013_pr)*exp((-33629.98784_pr)*R_T_inv + T_log *(3_pr))
    k(r42f) = (16_pr)*exp((-11640.01352_pr)*R_T_inv + T_log *(1.83_pr))
    k(r43f) = (1900_pr)*exp((-36300.00744_pr)*R_T_inv + T_log *(1.44_pr))
    k(r44f) = (39800000_pr)*exp((-238030.01936_pr)*R_T_inv)
    k(r45f) = (9030000_pr)*exp((-103099.99416_pr)*R_T_inv)
    k(r46f) = (180000000_pr)*exp((-63199.98944_pr)*R_T_inv)
    k(r47f) = (155000000_pr)*exp((-56399.98528_pr)*R_T_inv)
    k(r48f) = (40000000_pr)*exp((-10469.99976_pr)*R_T_inv)
    k(r49f) = (84300000_pr)
    k(r50f) = (42200000_pr)
    k(r51f) = (5000000_pr)
    k(r52f) = (330000_pr)*exp((-37409.9808_pr)*R_T_inv)
    k(r53f) = (11000000_pr)*exp((-116399.00552_pr)*R_T_inv)
    k(r54f) = (100000000_pr)*exp((-133900.00808_pr)*R_T_inv)
    k(r55f) = (31600000_pr)*exp((-61499.9884_pr)*R_T_inv)
    k_0(FOr56f) = (1.59e+21_pr)*exp((-10174.98592_pr)*R_T_inv + T_log *(-4.761_pr))
    k_inf(FOr56f) = (135100000_pr)*exp((-367.02048_pr)*R_T_inv + T_log *(0.091_pr))
    FC(FOr56f) = ((1.0_pr - 0.834_pr)*exp(-Tloc/(36.8_pr))) + (0.834_pr*exp(-Tloc/(778_pr)))&
              + exp(-(2464_pr)/Tloc)
    k(r56f) = getlindratecoeff(Tloc,k_0(FOr56f),k_inf(FOr56f),FC(FOr56f), M(mM56),Ploc)
    k_0(FOr57f) = (1.27e+29_pr)*exp((-11560.01544_pr)*R_T_inv + T_log *(-7_pr))
    k_inf(FOr57f) = (18100000_pr)
    FC(FOr57f) = ((1.0_pr - 0.62_pr)*exp(-Tloc/(73_pr))) + (0.62_pr*exp(-Tloc/(1200_pr)))
    k(r57f) = getlindratecoeff(Tloc,k_0(FOr57f),k_inf(FOr57f),FC(FOr57f), M(mM57),Ploc)
    k(r58f) = (30000000_pr)
    k(r59f) = (31300000_pr)
    k(r60f) = (3000000_pr)
    k(r61f) = (6000000_pr)
    k(r62f) = (6020000_pr)*exp((7479.98784_pr)*R_T_inv)
    k(r63f) = (25000000_pr)
    k(r64f) = (11.3_pr)*exp((-12549.99168_pr)*R_T_inv + T_log *(2_pr))
    k(r65f) = (80000000_pr)
    k(r66f) = (40000000_pr)
    k(r67f) = (2630000_pr)*exp((-6240.0176_pr)*R_T_inv)
    k(r68f) = (6580000_pr)*exp((-6240.0176_pr)*R_T_inv)
    k(r69f) = (40000000_pr)
    k(r70f) = (177000_pr)*exp((1999.99384_pr)*R_T_inv + T_log *(0.76_pr))
    k(r71f) = (1170000000_pr)*exp(T_log *(-0.75_pr))
    k(r72f) = (4.8e-05_pr)*exp((13500.01072_pr)*R_T_inv + T_log *(3.22_pr))
    k(r73f) = (20000000_pr)
    k(r74f) = (16000000_pr)
    k(r75f) = (5000000_pr)
    k(r76f) = (10000000_pr)
    k(r77f) = (4.28e-19_pr)*exp((14799.97952_pr)*R_T_inv + T_log *(7.6_pr))
    k(r78f) = (77800000_pr)*exp((-56539.98192_pr)*R_T_inv)
    k(r79f) = (0.00054_pr)*exp((-21800.02072_pr)*R_T_inv + T_log *(3.5_pr))
    k(r80f) = (1.4e-06_pr)*exp((-11600.01448_pr)*R_T_inv + T_log *(4.3_pr))
    k(r81f) = (22_pr)*exp((-4700.01272_pr)*R_T_inv + T_log *(1.9_pr))
    k(r82f) = (5.5e-07_pr)*exp((-34700.004_pr)*R_T_inv + T_log *(4_pr))
    k_0(FOr83f) = (4.9e+36_pr)*exp((-448525.0092_pr)*R_T_inv + T_log *(-6.43_pr))
    k_inf(FOr83f) = (8.85e+20_pr)*exp((-427699.986_pr)*R_T_inv + T_log *(-1.23_pr))
    FC(FOr83f) = ((1.0_pr - 0.84_pr)*exp(-Tloc/(125_pr))) + (0.84_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r83f) = getlindratecoeff(Tloc,k_0(FOr83f),k_inf(FOr83f),FC(FOr83f), M(mM83),Ploc)
    k(r84f) = (13200000_pr)*exp((-85646.01976_pr)*R_T_inv)
    k(r85f) = (30000000_pr)
    k(r86f) = (30600000_pr)
    k(r87f) = (42400000_pr)
    k(r88f) = (750000000_pr)*exp((-20082.9908_pr)*R_T_inv + T_log *(-1_pr))
    k_0(FOr89f) = (3.99e+27_pr)*exp((-167360_pr)*R_T_inv + T_log *(-4.99_pr))
    k_inf(FOr89f) = (11100000000_pr)*exp((-153839.98976_pr)*R_T_inv + T_log *(1.037_pr))
    FC(FOr89f) = ((1.0_pr - 0.168_pr)*exp(-Tloc/(1200_pr))) + (0.168_pr*exp(-Tloc/(0_pr)))
    k(r89f) = getlindratecoeff(Tloc,k_0(FOr89f),k_inf(FOr89f),FC(FOr89f), M(mM89),Ploc)
    k(r90f) = (44.9_pr)*exp((-55899.99728_pr)*R_T_inv + T_log *(2.12_pr))
    k(r91f) = (0.553_pr)*exp((-12399.99528_pr)*R_T_inv + T_log *(2.31_pr))
    k(r92f) = (2.25_pr)*exp(T_log *(2.08_pr))
    k(r93f) = (1.21_pr)*exp(T_log *(2.08_pr))
    k(r94f) = (42200000_pr)*exp((-241095.00856_pr)*R_T_inv)
    k(r95f) = (260000000000_pr)*exp((-404041.01408_pr)*R_T_inv)
    k(r96f) = (35000000000_pr)*exp((-299290.01352_pr)*R_T_inv)
    k(r97f) = (40000000_pr)
    k_0(FOr98f) = (151000000_pr)*exp((-136775.00184_pr)*R_T_inv + T_log *(0.1_pr))
    k_inf(FOr98f) = (6380000000_pr)*exp((-157429.98728_pr)*R_T_inv + T_log *(1_pr))
    FC(FOr98f) = ((1.0_pr - 0.3_pr)*exp(-Tloc/(1e+30_pr))) + (0.3_pr*exp(-Tloc/(0_pr)))
    k(r98f) = getlindratecoeff(Tloc,k_0(FOr98f),k_inf(FOr98f),FC(FOr98f), M(mM98),Ploc)
    k(r99f) = (1.7e+23_pr)*exp((-27209.01224_pr)*R_T_inv + T_log *(-5.312_pr))
    k(r100f) = (700000000_pr)*exp((-22018.00712_pr)*R_T_inv + T_log *(-0.611_pr))
    k(r101f) = (5190000000_pr)*exp((-13860.00208_pr)*R_T_inv + T_log *(-1.26_pr))
    k(r102f) = (400000000_pr)*exp((-44600.01744_pr)*R_T_inv)
    k(r103f) = (160000000_pr)*exp((-41400.01056_pr)*R_T_inv)
    k(r104f) = (4600000000_pr)*exp((-188000.00672_pr)*R_T_inv + T_log *(-0.54_pr))
    k(r105f) = (19_pr)*exp((-4179.98336_pr)*R_T_inv + T_log *(1.7_pr))
    k(r106f) = (1500_pr)*exp((-11249.98104_pr)*R_T_inv + T_log *(1.43_pr))
    k(r107f) = (20000000_pr)*exp((-9600.02064_pr)*R_T_inv)
    k(r108f) = (10000000_pr)*exp((-8370.00832_pr)*R_T_inv)
    k(r109f) = (90000_pr)
    k(r110f) = (150000000_pr)
    k(r111f) = (2000000_pr)
    k(r112f) = (96400000_pr)
    k(r113f) = (28.8_pr)*exp((-4189.98312_pr)*R_T_inv + T_log *(1.7_pr))
    k(r114f) = (14_pr)*exp((-4189.98312_pr)*R_T_inv + T_log *(1.7_pr))
    k(r115f) = (30000000_pr)
    k(r116f) = (250000000000_pr)*exp((-21450.99144_pr)*R_T_inv + T_log *(-0.93_pr))
    k(r117f) = (24000000_pr)
    k(r118f) = (5000000_pr)
    k(r119f) = (50000000_pr)*exp((-105099.988_pr)*R_T_inv)
    k(r120f) = (100000000_pr)*exp((-80000.00464_pr)*R_T_inv)
    k(r121f) = (10200000_pr)
    k(r122f) = (1.047e+37_pr)*exp((-185519.98256_pr)*R_T_inv + T_log *(-7.189_pr))
    k(r123f) = (50000000_pr)
    k(r124f) = (20000000_pr)
    k(r125f) = (100000000_pr)
    k(r126f) = (30000000_pr)
    k(r127f) = (30000_pr)
    k(r128f) = (490000000_pr)*exp(T_log *(-0.5_pr))
    k(r129f) = (7000000_pr)
    k(r130f) = (1.17e+43_pr)*exp((-183258.9908_pr)*R_T_inv + T_log *(-9.8_pr))
    k_0(FOr131f) = (0.002_pr)*exp(T_log *(1_pr))
    k_inf(FOr131f) = (600_pr)
    FC(FOr131f) = ((1.0_pr - 0.5_pr)*exp(-Tloc/(1e+30_pr))) + (0.5_pr*exp(-Tloc/(0_pr)))
    k(r131f) = getlindratecoeff(Tloc,k_0(FOr131f),k_inf(FOr131f),FC(FOr131f), M(mM131),Ploc)
    k(r132f) = (3.1_pr)*exp((1248.00352_pr)*R_T_inv + T_log *(2_pr))
    k(r133f) = (120_pr)*exp((-1370.00896_pr)*R_T_inv + T_log *(1.65_pr))
    k(r134f) = (0.17_pr)*exp((-10430.00072_pr)*R_T_inv + T_log *(2.5_pr))
    k_0(FOr135f) = (1.33e+48_pr)*exp((-24969.98648_pr)*R_T_inv + T_log *(-12_pr))
    k_inf(FOr135f) = (200000000_pr)
    FC(FOr135f) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(1097_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r135f) = getlindratecoeff(Tloc,k_0(FOr135f),k_inf(FOr135f),FC(FOr135f), M(mM135),Ploc)
    k(r136f) = (2660000_pr)
    k(r137f) = (3000000_pr)
    k_0(FOr138f) = (4.27e+46_pr)*exp((-40881.98952_pr)*R_T_inv + T_log *(-11.94_pr))
    k_inf(FOr138f) = (25000000_pr)
    FC(FOr138f) = ((1.0_pr - 0.175_pr)*exp(-Tloc/(1341_pr))) + (0.175_pr*exp(-Tloc/(60000_pr)))&
              + exp(-(10140_pr)/Tloc)
    k(r138f) = getlindratecoeff(Tloc,k_0(FOr138f),k_inf(FOr138f),FC(FOr138f), M(mM138),Ploc)
    k(r139f) = (1.6e+16_pr)*exp((-46800.00648_pr)*R_T_inv + T_log *(-2.39_pr))
    k(r140f) = (1.5e+18_pr)*exp((-77900.0132_pr)*R_T_inv + T_log *(-2.83_pr))
    k_0(FOr141f) = (7830000000000_pr)*exp((-271875.98528_pr)*R_T_inv)
    k_inf(FOr141f) = (1.1e+17_pr)*exp((-353100.01912_pr)*R_T_inv)
    FC(FOr141f) = ((1.0_pr - 0.76_pr)*exp(-Tloc/(1900_pr))) + (0.76_pr*exp(-Tloc/(38_pr)))
    k(r141f) = getlindratecoeff(Tloc,k_0(FOr141f),k_inf(FOr141f),FC(FOr141f), M(mM141),Ploc)
    k(r142f) = (40000000_pr)*exp((-198740_pr)*R_T_inv)
    k(r143f) = (40000000_pr)*exp((-213099.99008_pr)*R_T_inv)
    k(r144f) = (1.3_pr)*exp((-18706.99872_pr)*R_T_inv + T_log *(2.4_pr))
    k(r145f) = (1.33_pr)*exp((-28289.99048_pr)*R_T_inv + T_log *(2.54_pr))
    k(r146f) = (0.0476_pr)*exp((-8816.98504_pr)*R_T_inv + T_log *(2.71_pr))
    k(r147f) = (0.19_pr)*exp((-15557.9948_pr)*R_T_inv + T_log *(2.68_pr))
    k(r148f) = (10000_pr)*exp((-6693.9816_pr)*R_T_inv + T_log *(1_pr))
    k(r149f) = (20_pr)*exp((417.9816_pr)*R_T_inv + T_log *(-1.6_pr))
    k(r150f) = (0.00964_pr)*exp((-58229.9832_pr)*R_T_inv + T_log *(2.6_pr))
    k(r151f) = (0.0476_pr)*exp((-69000.0176_pr)*R_T_inv + T_log *(2.55_pr))
    k(r152f) = (8.4e-09_pr)*exp((-36300.00744_pr)*R_T_inv + T_log *(4.2_pr))
    k_0(FOr153f) = (8.7e+30_pr)*exp((-19798.98088_pr)*R_T_inv + T_log *(-7.5_pr))
    k_inf(FOr153f) = (13300000_pr)*exp((-6530.01064_pr)*R_T_inv)
    FC(FOr153f) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(645.4_pr)))&
              + exp(-(6844_pr)/Tloc)
    k(r153f) = getlindratecoeff(Tloc,k_0(FOr153f),k_inf(FOr153f),FC(FOr153f), M(mM153),Ploc)
    k(r154f) = (130000_pr)
    k_0(FOr155f) = (5.49e+43_pr)*exp((-149703.98024_pr)*R_T_inv + T_log *(-10_pr))
    k_inf(FOr155f) = (12300000000000_pr)*exp((-126400.02072_pr)*R_T_inv + T_log *(-0.1_pr))
    FC(FOr155f) = ((1.0_pr - -1.17_pr)*exp(-Tloc/(251_pr))) + (-1.17_pr*exp(-Tloc/(0_pr)))&
              + exp(-(1185_pr)/Tloc)
    k(r155f) = getlindratecoeff(Tloc,k_0(FOr155f),k_inf(FOr155f),FC(FOr155f), M(mM155),Ploc)
    k_0(FOr156f) = (6.26e+26_pr)*exp((-29288_pr)*R_T_inv + T_log *(-6.66_pr))
    k_inf(FOr156f) = (13300000_pr)*exp((-13640.00736_pr)*R_T_inv)
    FC(FOr156f) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(1310_pr)))&
              + exp(-(48100_pr)/Tloc)
    k(r156f) = getlindratecoeff(Tloc,k_0(FOr156f),k_inf(FOr156f),FC(FOr156f), M(mM156),Ploc)
    k(r157f) = (35000000000_pr)*exp((-14644_pr)*R_T_inv + T_log *(-1.6_pr))
    k(r158f) = (2000000_pr)
    k(r159f) = (2.5e+35_pr)*exp((-92048_pr)*R_T_inv + T_log *(-8.3_pr))
    k(r160f) = (150_pr)*exp((29288_pr)*R_T_inv)
    k(r161f) = (1e+15_pr)*exp((-179912_pr)*R_T_inv)
    k(r1b) = (132796281.261502_pr)*exp((-1277.45895479505_pr)*R_T_inv + T_log *(-0.342724717448034_pr))
    k(r2b) = (0.0269844924429187_pr)*exp((-20089.1470108947_pr)*R_T_inv + T_log *(2.64791248174178_pr))
    k(r3b) = (7731.59879497729_pr)*exp((-77473.9560444158_pr)*R_T_inv + T_log *(1.25337576194212_pr))
    k(r4b) = (0.0564908783135058_pr)*exp((7624.78991727421_pr)*R_T_inv + T_log *(2.35453672032186_pr))
    k(r5b) = (377258279297.854_pr)*exp((-433059.656278805_pr)*R_T_inv + T_log *(-0.683140689925532_pr))
    k(r6b) = (7.67076920981103e+16_pr)*exp((-495323.600888425_pr)*R_T_inv + T_log *(-1.72976494280662_pr))
    k(r7b) = (253105316174.772_pr)*exp((-496971.350080239_pr)*R_T_inv + T_log *(-0.562503547413838_pr))
    k(r8b) = (728920210106.337_pr)*exp((-426828.80722891_pr)*R_T_inv + T_log *(-0.705228184512192_pr))
    k_0(FOr9b) = (61144583709585.7_pr)*exp((-203946.482708775_pr)*R_T_inv + T_log *(-1.38426534756023_pr))
    k_inf(FOr9b) = (4944733982618.58_pr)*exp((-203946.482425346_pr)*R_T_inv + T_log *(0.455734702773399_pr))
    FC(FOr9b) = ((1.0_pr - 0.5_pr)*exp(-Tloc/(0_pr))) + (0.5_pr*exp(-Tloc/(1e+30_pr)))
    k(r9b) = getlindratecoeff(Tloc,k_0(FOr9b),k_inf(FOr9b),FC(FOr9b), M(mM9),Ploc)
    k(r10b) = (38872.7788501702_pr)*exp((-153973.769485512_pr)*R_T_inv + T_log *(0.63631240301974_pr))
    k(r11b) = (4530157.24624011_pr)*exp((-232556.18701885_pr)*R_T_inv + T_log *(0.301124684688439_pr))
    k(r12b) = (210908.338555442_pr)*exp((-228434.570385431_pr)*R_T_inv + T_log *(0.611775678530261_pr))
    k(r13b) = (2910710.49137187_pr)*exp((-222882.324660551_pr)*R_T_inv + T_log *(0.279037139624366_pr))
    k(r14b) = (12623703.6624577_pr)*exp((-286797.102530792_pr)*R_T_inv + T_log *(0.254500412512219_pr))
    k(r15b) = (811523802.126965_pr)*exp((-337107.108449863_pr)*R_T_inv + T_log *(0.254500413247054_pr))
    k_0(FOr16b) = (2.0294452119005e+24_pr)*exp((-212647.212423917_pr)*R_T_inv + T_log *(-4.034366717336_pr))
    k_inf(FOr16b) = (7.02217601061399e+18_pr)*exp((-212647.212575054_pr)*R_T_inv + T_log *(-1.10436674424188_pr))
    FC(FOr16b) = ((1.0_pr - 0.57_pr)*exp(-Tloc/(1e+30_pr))) + (0.57_pr*exp(-Tloc/(0_pr)))
    k(r16b) = getlindratecoeff(Tloc,k_0(FOr16b),k_inf(FOr16b),FC(FOr16b), M(mM16),Ploc)
    k(r17b) = (3910458470.58983_pr)*exp((-207640.533000139_pr)*R_T_inv + T_log *(-0.21378897325161_pr))
    k(r18b) = (7365329.31417772_pr)*exp((-155545.507135003_pr)*R_T_inv + T_log *(-0.213788969198506_pr))
    k(r19b) = (165326.596432793_pr)*exp((-100935.670758475_pr)*R_T_inv + T_log *(0.514913612148235_pr))
    k(r20b) = (260.801838723788_pr)*exp((-297676.404976674_pr)*R_T_inv + T_log *(1.10460178241294_pr))
    k(r21b) = (82650.8922154039_pr)*exp((-135936.587536768_pr)*R_T_inv + T_log *(0.46828938577721_pr))
    k(r22b) = (3605288.89246454_pr)*exp((-160366.586971498_pr)*R_T_inv + T_log *(0.468289386740391_pr))
    k(r23b) = (0.0369151761648382_pr)*exp((-78141.8303045318_pr)*R_T_inv + T_log *(2.49282608429797_pr))
    k_0(FOr24b) = (4.26978282376581e+23_pr)*exp((-550836.425353072_pr)*R_T_inv + T_log *(-3.6347357099319_pr))
    k_inf(FOr24b) = (4.9584565347813e+16_pr)*exp((-543276.397454135_pr)*R_T_inv + T_log *(-0.844735685404702_pr))
    FC(FOr24b) = ((1.0_pr - 1_pr)*exp(-Tloc/(1_pr))) + (1_pr*exp(-Tloc/(10000000_pr)))&
              + exp(-(10000000_pr)/Tloc)
    k(r24b) = getlindratecoeff(Tloc,k_0(FOr24b),k_inf(FOr24b),FC(FOr24b), M(mM24),Ploc)
    k(r25b) = (7831908.28998695_pr)*exp((-103372.59021153_pr)*R_T_inv + T_log *(0.360492500217439_pr))
    k(r26b) = (19545972395.3478_pr)*exp((-355212.37912014_pr)*R_T_inv + T_log *(-0.503195108388735_pr))
    k(r27b) = (6715187840.07503_pr)*exp((-235907.066285719_pr)*R_T_inv + T_log *(-0.782232216746173_pr))
    k(r28b) = (75585.7910444025_pr)*exp((-6542.95713027912_pr)*R_T_inv + T_log *(-0.852873369811746_pr))
    k(r29b) = (5896479.95672707_pr)*exp((-368472.605387643_pr)*R_T_inv + T_log *(0.463985887488034_pr))
    k(r30b) = (1886720.64430552_pr)*exp((-362241.756121677_pr)*R_T_inv + T_log *(0.441898432752029_pr))
    k(r31b) = (3358323969035.05_pr)*exp((-468714.355731817_pr)*R_T_inv + T_log *(-0.697609088210368_pr))
    k(r32b) = (23379076.234311_pr)*exp((-430736.549612428_pr)*R_T_inv + T_log *(0.417361702583794_pr))
    k(r33b) = (3275567.37112066_pr)*exp((-141074.411220096_pr)*R_T_inv + T_log *(0.162861293259345_pr))
    k(r34b) = (2424583919.56353_pr)*exp((-375199.158618956_pr)*R_T_inv + T_log *(0.111223157275263_pr))
    k_0(FOr35b) = (6.83617932315801e+19_pr)*exp((-369677.968717578_pr)*R_T_inv + T_log *(-2.55862622487201_pr))
    k_inf(FOr35b) = (55195818114615.9_pr)*exp((-366811.970556178_pr)*R_T_inv + T_log *(0.491373775430381_pr))
    FC(FOr35b) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r35b) = getlindratecoeff(Tloc,k_0(FOr35b),k_inf(FOr35b),FC(FOr35b), M(mM35),Ploc)
    k(r36b) = (0.328948144125278_pr)*exp((-76659.6948989427_pr)*R_T_inv + T_log *(2.20548559494605_pr))
    k(r37b) = (106966.433844106_pr)*exp((-73628.8113509731_pr)*R_T_inv + T_log *(0.283398020795668_pr))
    k(r38b) = (1476.94327869162_pr)*exp((-129123.623947019_pr)*R_T_inv + T_log *(1.14886129320022_pr))
    k(r39b) = (1259971.79001692_pr)*exp((-6226.51928210437_pr)*R_T_inv + T_log *(0.00436114605122311_pr))
    k(r40b) = (0.0327674568621369_pr)*exp((-40207.0376152626_pr)*R_T_inv + T_log *(2.29057190666917_pr))
    k(r41b) = (3.16155490236012e-05_pr)*exp((-26903.4375857454_pr)*R_T_inv + T_log *(3.35276220125743_pr))
    k(r42b) = (0.257134754120746_pr)*exp((-67177.4077448366_pr)*R_T_inv + T_log *(2.13613797121124_pr))
    k(r43b) = (2.46419278943352_pr)*exp((-23342.6082238415_pr)*R_T_inv + T_log *(1.77067469229196_pr))
    k(r44b) = (354678.876513258_pr)*exp((-2190.29579342555_pr)*R_T_inv + T_log *(0.0516374969192687_pr))
    k(r45b) = (3055135.53162912_pr)*exp((-28700.78275748_pr)*R_T_inv + T_log *(-0.162151483887863_pr))
    k(r46b) = (1087643.514692_pr)*exp((-33544.3703984629_pr)*R_T_inv + T_log *(0.405275263938278_pr))
    k(r47b) = (1063949.50852216_pr)*exp((10834.3610289363_pr)*R_T_inv + T_log *(0.489060625715628_pr))
    k(r48b) = (1814398.76019352_pr)*exp((-5499.59794269533_pr)*R_T_inv + T_log *(0.442436393289239_pr))
    k(r49b) = (4414091659.68467_pr)*exp((-287836.717963824_pr)*R_T_inv + T_log *(-0.162246803781_pr))
    k(r50b) = (49601489718957.2_pr)*exp((-273485.221341383_pr)*R_T_inv + T_log *(-0.98877115645228_pr))
    k(r51b) = (637960904.132582_pr)*exp((-102927.361251273_pr)*R_T_inv + T_log *(-0.597885473735203_pr))
    k(r52b) = (65188.5142937241_pr)*exp((-255104.156567395_pr)*R_T_inv + T_log *(0.195028442799429_pr))
    k(r53b) = (9643789134.52329_pr)*exp((3555.95792746122_pr)*R_T_inv + T_log *(-0.876922606000626_pr))
    k(r54b) = (710225867164.596_pr)*exp((-377729.610627335_pr)*R_T_inv + T_log *(-0.583495936471742_pr))
    k(r55b) = (937519251085.236_pr)*exp((-22758.2129423332_pr)*R_T_inv + T_log *(-0.977794682399029_pr))
    k_0(FOr56b) = (1.89729627341367e+29_pr)*exp((-449961.192416399_pr)*R_T_inv + T_log *(-4.79690288452907_pr))
    k_inf(FOr56b) = (1.61209881432389e+16_pr)*exp((-440153.224081238_pr)*R_T_inv + T_log *(0.055097629497815_pr))
    FC(FOr56b) = ((1.0_pr - 0.834_pr)*exp(-Tloc/(36.8_pr))) + (0.834_pr*exp(-Tloc/(778_pr)))&
              + exp(-(2464_pr)/Tloc)
    k(r56b) = getlindratecoeff(Tloc,k_0(FOr56b),k_inf(FOr56b),FC(FOr56b), M(mM56),Ploc)
    k_0(FOr57b) = (1.76135397213646e+42_pr)*exp((-392857.96185687_pr)*R_T_inv + T_log *(-8.08071381760312_pr))
    k_inf(FOr57b) = (2.51027661550021e+20_pr)*exp((-381297.94655639_pr)*R_T_inv + T_log *(-1.08071384255436_pr))
    FC(FOr57b) = ((1.0_pr - 0.62_pr)*exp(-Tloc/(73_pr))) + (0.62_pr*exp(-Tloc/(1200_pr)))
    k(r57b) = getlindratecoeff(Tloc,k_0(FOr57b),k_inf(FOr57b),FC(FOr57b), M(mM57),Ploc)
    k(r58b) = (429123425011.447_pr)*exp((-361301.913610816_pr)*R_T_inv + T_log *(-0.629219986648906_pr))
    k(r59b) = (2.09775316629084_pr)*exp((-285501.148015687_pr)*R_T_inv + T_log *(0.158580205141468_pr))
    k(r60b) = (24108.3338598667_pr)*exp((-254829.31393606_pr)*R_T_inv + T_log *(0.510287545656794_pr))
    k(r61b) = (5281727.41819362_pr)*exp((-37578.7272419531_pr)*R_T_inv + T_log *(-0.0837853570858638_pr))
    k(r62b) = (789508.712214611_pr)*exp((-4578.78691948731_pr)*R_T_inv + T_log *(0.197052444005639_pr))
    k(r63b) = (406233810472.778_pr)*exp((-323723.18607299_pr)*R_T_inv + T_log *(-0.545434576801976_pr))
    k(r64b) = (9.79314821322167_pr)*exp((-86872.7109642109_pr)*R_T_inv + T_log *(2.15042820583261_pr))
    k(r65b) = (1614.48101139884_pr)*exp((-318064.964657318_pr)*R_T_inv + T_log *(-0.114909987385297_pr))
    k(r66b) = (234260017.34198_pr)*exp((-751124.620599222_pr)*R_T_inv + T_log *(0.201949383183794_pr))
    k(r67b) = (103431350606.543_pr)*exp((-793694.695456978_pr)*R_T_inv + T_log *(-0.580282864049124_pr))
    k(r68b) = (0.500970013327124_pr)*exp((-254162.439232151_pr)*R_T_inv + T_log *(0.242365404687159_pr))
    k(r69b) = (1786233196.29809_pr)*exp((-739065.84654749_pr)*R_T_inv + T_log *(0.00489681341698843_pr))
    k(r70b) = (474141.502694601_pr)*exp((-304681.553573017_pr)*R_T_inv + T_log *(0.680273771757537_pr))
    k(r71b) = (21937040843472.4_pr)*exp((-249400.466801886_pr)*R_T_inv + T_log *(-1.44586278507641_pr))
    k(r72b) = (1.91477497536833e-08_pr)*exp((-256851.479549289_pr)*R_T_inv + T_log *(3.92250599877951_pr))
    k(r73b) = (8450.22813620844_pr)*exp((-343879.988217087_pr)*R_T_inv + T_log *(1.09403855151573_pr))
    k(r74b) = (3.12306549145816_pr)*exp((-44842.0194498925_pr)*R_T_inv + T_log *(1.67663424306421_pr))
    k(r75b) = (13960.2072603002_pr)*exp((-406143.932697785_pr)*R_T_inv + T_log *(1.04741432069848_pr))
    k(r76b) = (2253.21260683798_pr)*exp((-337649.139302167_pr)*R_T_inv + T_log *(1.07195103404382_pr))
    k(r77b) = (6.62638769187194e-22_pr)*exp((-99966.8349790883_pr)*R_T_inv + T_log *(8.39291391955515_pr))
    k(r78b) = (0.113272001011377_pr)*exp((32639.6862035955_pr)*R_T_inv + T_log *(0.77717925194639_pr))
    k(r79b) = (3.3522730674248e-07_pr)*exp((-34819.9551439832_pr)*R_T_inv + T_log *(3.91977844192167_pr))
    k(r80b) = (4.63486837946379e-10_pr)*exp((-18389.1000030037_pr)*R_T_inv + T_log *(4.69769092217837_pr))
    k(r81b) = (0.0902509479646934_pr)*exp((-79983.8916934692_pr)*R_T_inv + T_log *(2.27315419919512_pr))
    k(r82b) = (1.40394781305603e-07_pr)*exp((-54446.4887159319_pr)*R_T_inv + T_log *(4.06701623380999_pr))
    k_0(FOr83b) = (1.04820509793461e+28_pr)*exp((-28485.2873216619_pr)*R_T_inv + T_log *(-6.32708086422556_pr))
    k_inf(FOr83b) = (1893187216307_pr)*exp((-7660.26430369486_pr)*R_T_inv + T_log *(-1.12708089543282_pr))
    FC(FOr83b) = ((1.0_pr - 0.84_pr)*exp(-Tloc/(125_pr))) + (0.84_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r83b) = getlindratecoeff(Tloc,k_0(FOr83b),k_inf(FOr83b),FC(FOr83b), M(mM83),Ploc)
    k(r84b) = (1139999.2406062_pr)*exp((-30993.2924427981_pr)*R_T_inv + T_log *(-0.0951351378992994_pr))
    k(r85b) = (7181656.91689163_pr)*exp((-282571.378045101_pr)*R_T_inv + T_log *(0.394298738281577_pr))
    k(r86b) = (3906506.23026556_pr)*exp((-276340.529068546_pr)*R_T_inv + T_log *(0.372211232483691_pr))
    k(r87b) = (74831.8637260831_pr)*exp((-326578.493498927_pr)*R_T_inv + T_log *(0.815547865415727_pr))
    k(r88b) = (657899025.617812_pr)*exp((-73541.1950395366_pr)*R_T_inv + T_log *(-0.906825877441833_pr))
    k_0(FOr89b) = (3.29140063609936e+21_pr)*exp((-16871.7215662511_pr)*R_T_inv + T_log *(-4.91256053619032_pr))
    k_inf(FOr89b) = (9156.52963860351_pr)*exp((-3351.71145961802_pr)*R_T_inv + T_log *(1.11443944210628_pr))
    FC(FOr89b) = ((1.0_pr - 0.168_pr)*exp(-Tloc/(1200_pr))) + (0.168_pr*exp(-Tloc/(0_pr)))
    k(r89b) = getlindratecoeff(Tloc,k_0(FOr89b),k_inf(FOr89b),FC(FOr89b), M(mM89),Ploc)
    k(r90b) = (0.0706876778789452_pr)*exp((-24237.8164928914_pr)*R_T_inv + T_log *(2.45159145576011_pr))
    k(r91b) = (0.00575315089204133_pr)*exp((-43001.7586957207_pr)*R_T_inv + T_log *(2.59496727391991_pr))
    k(r92b) = (9.50639220899262e-05_pr)*exp((-109166.794076378_pr)*R_T_inv + T_log *(2.80673461577417_pr))
    k(r93b) = (1682.40260231383_pr)*exp((-72559.5501013055_pr)*R_T_inv + T_log *(1.07473406140143_pr))
    k(r94b) = (243446.867739692_pr)*exp((19680.3463974452_pr)*R_T_inv + T_log *(0.0304669039364054_pr))
    k(r95b) = (1410.50409798784_pr)*exp((60680.8245756321_pr)*R_T_inv + T_log *(0.0147324272816156_pr))
    k(r96b) = (887.317456061081_pr)*exp((-122666.204097675_pr)*R_T_inv + T_log *(0.0144994089223173_pr))
    k(r97b) = (186914864.175687_pr)*exp((-288097.985331569_pr)*R_T_inv + T_log *(-0.000225204160703346_pr))
    k_0(FOr98b) = (2431.57166700907_pr)*exp((8186.63264593708_pr)*R_T_inv + T_log *(-0.217091006438007_pr))
    k_inf(FOr98b) = (102737.892098096_pr)*exp((-12468.3525353712_pr)*R_T_inv + T_log *(0.682909039892877_pr))
    FC(FOr98b) = ((1.0_pr - 0.3_pr)*exp(-Tloc/(1e+30_pr))) + (0.3_pr*exp(-Tloc/(0_pr)))
    k(r98b) = getlindratecoeff(Tloc,k_0(FOr98b),k_inf(FOr98b),FC(FOr98b), M(mM98),Ploc)
    k(r99b) = (1.68996722741507e+21_pr)*exp((-391963.011982765_pr)*R_T_inv + T_log *(-4.6997409143664_pr))
    k(r100b) = (4373457216322.75_pr)*exp((-62328.0447565959_pr)*R_T_inv + T_log *(-1.56849463587637_pr))
    k(r101b) = (88872520088.8422_pr)*exp((-72844.8500714125_pr)*R_T_inv + T_log *(-1.5613563129369_pr))
    k(r102b) = (1020028847.10885_pr)*exp((-126527.097039789_pr)*R_T_inv + T_log *(-0.150780774488602_pr))
    k(r103b) = (654.767669340815_pr)*exp((-232947.950099406_pr)*R_T_inv + T_log *(1.26463598678372_pr))
    k(r104b) = (1153995.48848784_pr)*exp((-633128.589468261_pr)*R_T_inv + T_log *(0.536476774105956_pr))
    k(r105b) = (41501.2163252567_pr)*exp((-101899.396954885_pr)*R_T_inv + T_log *(1.09666869665431_pr))
    k(r106b) = (0.000872114396363547_pr)*exp((-140964.974905931_pr)*R_T_inv + T_log *(2.9147795495044_pr))
    k(r107b) = (66696.7871412729_pr)*exp((-209901.146044741_pr)*R_T_inv + T_log *(0.728459795878215_pr))
    k(r108b) = (11674.6869204283_pr)*exp((7422.32566639724_pr)*R_T_inv + T_log *(0.452550527542136_pr))
    k(r109b) = (1552.45080132631_pr)*exp((-90973.2184560916_pr)*R_T_inv + T_log *(0.506984858711749_pr))
    k(r110b) = (273.452098067941_pr)*exp((-72042.1326175416_pr)*R_T_inv + T_log *(1.49920213235818_pr))
    k(r111b) = (0.00102991441166647_pr)*exp((-65444.0680564133_pr)*R_T_inv + T_log *(0.858608415816998_pr))
    k(r112b) = (0.00312201211986639_pr)*exp((-427685.824286051_pr)*R_T_inv + T_log *(1.30050682930822_pr))
    k(r113b) = (3.51878974588193e-12_pr)*exp((-361733.265050855_pr)*R_T_inv + T_log *(3.35778210316527_pr))
    k(r114b) = (3.04469511228076e-06_pr)*exp((-468205.864589144_pr)*R_T_inv + T_log *(2.21827459516449_pr))
    k(r115b) = (1119812.78963549_pr)*exp((-307864.987095841_pr)*R_T_inv + T_log *(0.666729020151121_pr))
    k(r116b) = (95041686.7382007_pr)*exp((-35248.4116492695_pr)*R_T_inv + T_log *(-0.12311169243674_pr))
    k(r117b) = (5919962.08464534_pr)*exp((-370128.931733372_pr)*R_T_inv + T_log *(0.620104762024867_pr))
    k(r118b) = (683894.173070168_pr)*exp((-78751.8135416192_pr)*R_T_inv + T_log *(0.365604359527175_pr))
    k(r119b) = (6.43129945518984_pr)*exp((20094.6812071124_pr)*R_T_inv + T_log *(0.349869714229434_pr))
    k(r120b) = (1131916.12747205_pr)*exp((-116015.005723146_pr)*R_T_inv + T_log *(0.42730953787722_pr))
    k(r121b) = (15599.4143789583_pr)*exp((-115917.573746136_pr)*R_T_inv + T_log *(0.677891226115172_pr))
    k(r122b) = (2.22372130358313e+29_pr)*exp((-27825.181487297_pr)*R_T_inv + T_log *(-6.79465236620169_pr))
    k(r123b) = (1.51935291667142_pr)*exp((-36607.2439932304_pr)*R_T_inv + T_log *(1.7320005512336_pr))
    k(r124b) = (123270.498370355_pr)*exp((-275364.855296911_pr)*R_T_inv + T_log *(0.71120692777944_pr))
    k(r125b) = (159.111847342632_pr)*exp((-324443.962104549_pr)*R_T_inv + T_log *(1.5697537217619_pr))
    k(r126b) = (1221894.59272035_pr)*exp((-337628.799626305_pr)*R_T_inv + T_log *(0.664582724362025_pr))
    k(r127b) = (7.31803058668378e-11_pr)*exp((-189714.368477031_pr)*R_T_inv + T_log *(2.07415564079724_pr))
    k(r128b) = (0.179516817912891_pr)*exp((66721.5825322113_pr)*R_T_inv + T_log *(0.401332477588984_pr))
    k(r129b) = (3.9514244253821e-08_pr)*exp((-50354.93702576_pr)*R_T_inv + T_log *(1.91129434601947_pr))
    k(r130b) = (1.44478000272462e+29_pr)*exp((-155279.183574078_pr)*R_T_inv + T_log *(-7.92087281327732_pr))
    k_0(FOr131b) = (16136039377.6447_pr)*exp((-207907.720268506_pr)*R_T_inv + T_log *(-0.438266960552671_pr))
    k_inf(FOr131b) = (4.84081209567616e+15_pr)*exp((-207907.720310142_pr)*R_T_inv + T_log *(-1.43826696815772_pr))
    FC(FOr131b) = ((1.0_pr - 0.5_pr)*exp(-Tloc/(1e+30_pr))) + (0.5_pr*exp(-Tloc/(0_pr)))
    k(r131b) = getlindratecoeff(Tloc,k_0(FOr131b),k_inf(FOr131b),FC(FOr131b), M(mM131),Ploc)
    k(r132b) = (4.42528681312071_pr)*exp((-127528.275534118_pr)*R_T_inv + T_log *(1.80919830591043_pr))
    k(r133b) = (3.74266985094159e-09_pr)*exp((48536.8121271288_pr)*R_T_inv + T_log *(2.3186706892314_pr))
    k(r134b) = (0.0367235980894567_pr)*exp((-76942.3352883342_pr)*R_T_inv + T_log *(2.35582253724501_pr))
    k_0(FOr135b) = (1.78669648684422e+54_pr)*exp((-391517.308223941_pr)*R_T_inv + T_log *(-11.5389632330121_pr))
    k_inf(FOr135b) = (268676301280138_pr)*exp((-366547.322118902_pr)*R_T_inv + T_log *(0.461036700387404_pr))
    FC(FOr135b) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(1097_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r135b) = getlindratecoeff(Tloc,k_0(FOr135b),k_inf(FOr135b),FC(FOr135b), M(mM135),Ploc)
    k(r136b) = (3360397.25050102_pr)*exp((-162600.839017673_pr)*R_T_inv + T_log *(0.445302117722778_pr))
    k(r137b) = (4.28954497330382e-09_pr)*exp((49198.3926359801_pr)*R_T_inv + T_log *(1.93465174193552_pr))
    k_0(FOr138b) = (7.45253516140443e+60_pr)*exp((-470375.396943466_pr)*R_T_inv + T_log *(-13.2343211820906_pr))
    k_inf(FOr138b) = (4.36331139555938e+21_pr)*exp((-429493.407492946_pr)*R_T_inv + T_log *(-1.29432119425044_pr))
    FC(FOr138b) = ((1.0_pr - 0.175_pr)*exp(-Tloc/(1341_pr))) + (0.175_pr*exp(-Tloc/(60000_pr)))&
              + exp(-(10140_pr)/Tloc)
    k(r138b) = getlindratecoeff(Tloc,k_0(FOr138b),k_inf(FOr138b),FC(FOr138b), M(mM138),Ploc)
    k(r139b) = (16898262943.3053_pr)*exp((-82028.4364398054_pr)*R_T_inv + T_log *(-1.11041101920375_pr))
    k(r140b) = (1.94880452024033e+26_pr)*exp((-140846.098885359_pr)*R_T_inv + T_log *(-4.58535794995739_pr))
    k_0(FOr141b) = (0.0077314857003756_pr)*exp((102199.758973207_pr)*R_T_inv + T_log *(1.55501313283501_pr))
    k_inf(FOr141b) = (108.615976771575_pr)*exp((20975.7254244849_pr)*R_T_inv + T_log *(1.55501318271976_pr))
    FC(FOr141b) = ((1.0_pr - 0.76_pr)*exp(-Tloc/(1900_pr))) + (0.76_pr*exp(-Tloc/(38_pr)))
    k(r141b) = getlindratecoeff(Tloc,k_0(FOr141b),k_inf(FOr141b),FC(FOr141b), M(mM141),Ploc)
    k(r142b) = (3799.50885917369_pr)*exp((10233.5485911267_pr)*R_T_inv + T_log *(0.592728223741924_pr))
    k(r143b) = (71451.1165254303_pr)*exp((3045.65431901644_pr)*R_T_inv + T_log *(0.12941443375889_pr))
    k(r144b) = (3.36989346444262e-05_pr)*exp((-38846.6238165371_pr)*R_T_inv + T_log *(3.29385286064002_pr))
    k(r145b) = (0.000648344731095291_pr)*exp((-41257.5197121933_pr)*R_T_inv + T_log *(2.97053908026138_pr))
    k(r146b) = (6.58026677341199e-07_pr)*exp((-22725.7612146148_pr)*R_T_inv + T_log *(3.5817653447125_pr))
    k(r147b) = (4.9393712654647e-05_pr)*exp((-22294.6751031144_pr)*R_T_inv + T_log *(3.08845156551798_pr))
    k(r148b) = (32.2134911863125_pr)*exp((-81925.4553597226_pr)*R_T_inv + T_log *(1.38391484163582_pr))
    k(r149b) = (0.00342599128506903_pr)*exp((-81985.5880059171_pr)*R_T_inv + T_log *(-0.75277137480976_pr))
    k(r150b) = (3.47644182257949e-05_pr)*exp((-10696.9467416049_pr)*R_T_inv + T_log *(2.97893924787738_pr))
    k(r151b) = (0.00322809628272877_pr)*exp((-14294.885394503_pr)*R_T_inv + T_log *(2.46562544688082_pr))
    k(r152b) = (1.57964978452216e-07_pr)*exp((-29127.911576728_pr)*R_T_inv + T_log *(3.73668621940584_pr))
    k_0(FOr153b) = (1.00764996273815e+36_pr)*exp((-166671.402371657_pr)*R_T_inv + T_log *(-7.27587006998354_pr))
    k_inf(FOr153b) = (1540430152200.01_pr)*exp((-153402.43201168_pr)*R_T_inv + T_log *(0.2241299510989_pr))
    FC(FOr153b) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(645.4_pr)))&
              + exp(-(6844_pr)/Tloc)
    k(r153b) = getlindratecoeff(Tloc,k_0(FOr153b),k_inf(FOr153b),FC(FOr153b), M(mM153),Ploc)
    k(r154b) = (1193556.74025728_pr)*exp((-57074.0612347226_pr)*R_T_inv + T_log *(-0.208395280564794_pr))
    k_0(FOr155b) = (2.66208826392825e+31_pr)*exp((-45232.0845593488_pr)*R_T_inv + T_log *(-8.48122716635297_pr))
    k_inf(FOr155b) = (5.96424149675058_pr)*exp((-21928.1250426598_pr)*R_T_inv + T_log *(1.41877283295164_pr))
    FC(FOr155b) = ((1.0_pr - -1.17_pr)*exp(-Tloc/(251_pr))) + (-1.17_pr*exp(-Tloc/(0_pr)))&
              + exp(-(1185_pr)/Tloc)
    k(r155b) = getlindratecoeff(Tloc,k_0(FOr155b),k_inf(FOr155b),FC(FOr155b), M(mM155),Ploc)
    k_0(FOr156b) = (1.36347225388333e+33_pr)*exp((-168988.325634904_pr)*R_T_inv + T_log *(-6.89918385190694_pr))
    k_inf(FOr156b) = (28968336439253.9_pr)*exp((-153340.332899562_pr)*R_T_inv + T_log *(-0.239183834977203_pr))
    FC(FOr156b) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(1310_pr)))&
              + exp(-(48100_pr)/Tloc)
    k(r156b) = getlindratecoeff(Tloc,k_0(FOr156b),k_inf(FOr156b),FC(FOr156b), M(mM156),Ploc)
    k(r157b) = (17087801774.2972_pr)*exp((-78890.1570342813_pr)*R_T_inv + T_log *(-1.34508148866925_pr))
    k(r158b) = (5.58615470025163e+18_pr)*exp((-97014.4457455474_pr)*R_T_inv + T_log *(-1.76086716283669_pr))
    k(r159b) = (4.36993694896037e+22_pr)*exp((-59279.7112422598_pr)*R_T_inv + T_log *(-6.28421431750243_pr))
    k(r160b) = (1140366.78820435_pr)*exp((-221066.586171088_pr)*R_T_inv + T_log *(-1.00545065805924_pr))
    k(r161b) = (2.36624007483949e-14_pr)*exp((51817.744117877_pr)*R_T_inv + T_log *(3.83212843860929_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(sH) * c(sO2) 
    w(r2f) = k(r2f) * c(sH2) * c(sO) 
    w(r3f) = k(r3f) * c(sH2) * c(sOH) 
    w(r4f) = k(r4f) * c(sH2O) * c(sO) 
    w(r5f) = k(r5f) * c(sH)**2.0_pr * m(mM5)
    w(r6f) = k(r6f) * c(sH) * c(sOH) * m(mM6)
    w(r7f) = k(r7f) * c(sO)**2.0_pr * m(mM7)
    w(r8f) = k(r8f) * c(sH) * c(sO) * m(mM8)
    w(r9f) = k(r9f) * c(sH) * c(sO2) 
    w(r10f) = k(r10f) * c(sH) * c(sHO2) 
    w(r11f) = k(r11f) * c(sH) * c(sHO2) 
    w(r12f) = k(r12f) * c(sH) * c(sHO2) 
    w(r13f) = k(r13f) * c(sHO2) * c(sO) 
    w(r14f) = k(r14f) * c(sHO2) * c(sOH) 
    w(r15f) = k(r15f) * c(sHO2) * c(sOH) 
    w(r16f) = k(r16f) * c(sOH)**2.0_pr 
    w(r17f) = k(r17f) * c(sHO2)**2.0_pr 
    w(r18f) = k(r18f) * c(sHO2)**2.0_pr 
    w(r19f) = k(r19f) * c(sH) * c(sH2O2) 
    w(r20f) = k(r20f) * c(sH) * c(sH2O2) 
    w(r21f) = k(r21f) * c(sH2O2) * c(sOH) 
    w(r22f) = k(r22f) * c(sH2O2) * c(sOH) 
    w(r23f) = k(r23f) * c(sH2O2) * c(sO) 
    w(r24f) = k(r24f) * c(sCO) * c(sO) 
    w(r25f) = k(r25f) * c(sCO) * c(sOH) 
    w(r26f) = k(r26f) * c(sCO) * c(sHO2) 
    w(r27f) = k(r27f) * c(sCO) * c(sO2) 
    w(r28f) = k(r28f) * cqss(sqssHCO) * m(mM28)
    w(r29f) = k(r29f) * c(sH) * cqss(sqssHCO) 
    w(r30f) = k(r30f) * cqss(sqssHCO) * c(sO) 
    w(r31f) = k(r31f) * cqss(sqssHCO) * c(sO) 
    w(r32f) = k(r32f) * cqss(sqssHCO) * c(sOH) 
    w(r33f) = k(r33f) * cqss(sqssHCO) * c(sO2) 
    w(r34f) = k(r34f) * c(sCH3) * cqss(sqssHCO) 
    w(r35f) = k(r35f) * c(sH) * cqss(sqssHCO) 
    w(r36f) = k(r36f) * c(sCH2O) * c(sH) 
    w(r37f) = k(r37f) * c(sCH2O) * c(sO) 
    w(r38f) = k(r38f) * c(sCH2O) * c(sOH) 
    w(r39f) = k(r39f) * c(sCH2O) * c(sO2) 
    w(r40f) = k(r40f) * c(sCH2O) * c(sHO2) 
    w(r41f) = k(r41f) * c(sCH4) * c(sH) 
    w(r42f) = k(r42f) * c(sCH4) * c(sOH) 
    w(r43f) = k(r43f) * c(sCH4) * c(sO) 
    w(r44f) = k(r44f) * c(sCH4) * c(sO2) 
    w(r45f) = k(r45f) * c(sCH4) * c(sHO2) 
    w(r46f) = k(r46f) * c(sCH3) * c(sH) 
    w(r47f) = k(r47f) * c(sCH3) * c(sH) 
    w(r48f) = k(r48f) * c(sCH3) * c(sOH) 
    w(r49f) = k(r49f) * c(sCH3) * c(sO) 
    w(r50f) = k(r50f) * c(sCH3) * cqss(sqssTXCH2) 
    w(r51f) = k(r51f) * c(sCH3) * c(sHO2) 
    w(r52f) = k(r52f) * c(sCH3) * c(sO2) 
    w(r53f) = k(r53f) * c(sCH3) * c(sO2) 
    w(r54f) = k(r54f) * c(sCH3)**2.0_pr 
    w(r55f) = k(r55f) * c(sCH3)**2.0_pr 
    w(r56f) = k(r56f) * c(sCH3) * c(sH) 
    w(r57f) = k(r57f) * c(sCH3)**2.0_pr 
    w(r58f) = k(r58f) * c(sOH) * cqss(sqssSXCH2) 
    w(r59f) = k(r59f) * c(sO2) * cqss(sqssSXCH2) 
    w(r60f) = k(r60f) * c(sCO2) * cqss(sqssSXCH2) 
    w(r61f) = k(r61f) * cqss(sqssSXCH2) * m(mM61)
    w(r62f) = k(r62f) * c(sH) * cqss(sqssTXCH2) 
    w(r63f) = k(r63f) * c(sOH) * cqss(sqssTXCH2) 
    w(r64f) = k(r64f) * c(sOH) * cqss(sqssTXCH2) 
    w(r65f) = k(r65f) * c(sO) * cqss(sqssTXCH2) 
    w(r66f) = k(r66f) * c(sO) * cqss(sqssTXCH2) 
    w(r67f) = k(r67f) * c(sO2) * cqss(sqssTXCH2) 
    w(r68f) = k(r68f) * c(sO2) * cqss(sqssTXCH2) 
    w(r69f) = k(r69f) * cqss(sqssCH) * c(sO) 
    w(r70f) = k(r70f) * cqss(sqssCH) * c(sO2) 
    w(r71f) = k(r71f) * cqss(sqssCH) * c(sH2O) 
    w(r72f) = k(r72f) * cqss(sqssCH) * c(sCO2) 
    w(r73f) = k(r73f) * cqss(sqssCH3O) * c(sH) 
    w(r74f) = k(r74f) * cqss(sqssCH3O) * c(sH) 
    w(r75f) = k(r75f) * cqss(sqssCH3O) * c(sOH) 
    w(r76f) = k(r76f) * cqss(sqssCH3O) * c(sO) 
    w(r77f) = k(r77f) * cqss(sqssCH3O) * c(sO2) 
    w(r78f) = k(r78f) * cqss(sqssCH3O) * m(mM78)
    w(r79f) = k(r79f) * c(sC2H6) * c(sH) 
    w(r80f) = k(r80f) * c(sC2H6) * c(sO) 
    w(r81f) = k(r81f) * c(sC2H6) * c(sOH) 
    w(r82f) = k(r82f) * c(sC2H6) * c(sCH3) 
    w(r83f) = k(r83f) * c(sC2H6) 
    w(r84f) = k(r84f) * c(sC2H6) * c(sHO2) 
    w(r85f) = k(r85f) * cqss(sqssC2H5) * c(sH) 
    w(r86f) = k(r86f) * cqss(sqssC2H5) * c(sO) 
    w(r87f) = k(r87f) * cqss(sqssC2H5) * c(sO) 
    w(r88f) = k(r88f) * cqss(sqssC2H5) * c(sO2) 
    w(r89f) = k(r89f) * cqss(sqssC2H5) 
    w(r90f) = k(r90f) * c(sC2H4) * c(sH) 
    w(r91f) = k(r91f) * c(sC2H4) * c(sOH) 
    w(r92f) = k(r92f) * c(sC2H4) * c(sO) 
    w(r93f) = k(r93f) * c(sC2H4) * c(sO) 
    w(r94f) = k(r94f) * c(sC2H4) * c(sO2) 
    w(r95f) = k(r95f) * c(sC2H4) * m(mM95)
    w(r96f) = k(r96f) * c(sC2H4) * m(mM96)
    w(r97f) = k(r97f) * cqss(sqssC2H3) * c(sH) 
    w(r98f) = k(r98f) * cqss(sqssC2H3) 
    w(r99f) = k(r99f) * cqss(sqssC2H3) * c(sO2) 
    w(r100f) = k(r100f) * cqss(sqssC2H3) * c(sO2) 
    w(r101f) = k(r101f) * cqss(sqssC2H3) * c(sO2) 
    w(r102f) = k(r102f) * c(sC2H2) * c(sO) 
    w(r103f) = k(r103f) * c(sC2H2) * c(sO) 
    w(r104f) = k(r104f) * c(sC2H2) * c(sO2) 
    w(r105f) = k(r105f) * c(sC2H2) * c(sOH) 
    w(r106f) = k(r106f) * c(sCH2CO) * c(sH) 
    w(r107f) = k(r107f) * c(sCH2CO) * c(sO) 
    w(r108f) = k(r108f) * c(sCH2CO) * c(sO) 
    w(r109f) = k(r109f) * c(sCH2CO) * c(sCH3) 
    w(r110f) = k(r110f) * c(sH) * cqss(sqssHCCO) 
    w(r111f) = k(r111f) * cqss(sqssHCCO) * c(sOH) 
    w(r112f) = k(r112f) * cqss(sqssHCCO) * c(sO) 
    w(r113f) = k(r113f) * cqss(sqssHCCO) * c(sO2) 
    w(r114f) = k(r114f) * cqss(sqssHCCO) * c(sO2) 
    w(r115f) = k(r115f) * cqss(sqssCH2OH) * c(sH) 
    w(r116f) = k(r116f) * cqss(sqssCH2OH) * c(sH) 
    w(r117f) = k(r117f) * cqss(sqssCH2OH) * c(sOH) 
    w(r118f) = k(r118f) * cqss(sqssCH2OH) * c(sO2) 
    w(r119f) = k(r119f) * cqss(sqssCH2OH) * m(mM119)
    w(r120f) = k(r120f) * cqss(sqssCH3O) * m(mM120)
    w(r121f) = k(r121f) * c(sCH2CO) * c(sOH) 
    w(r122f) = k(r122f) * cqss(sqssCH2CHO) 
    w(r123f) = k(r123f) * cqss(sqssCH2CHO) * c(sH) 
    w(r124f) = k(r124f) * cqss(sqssCH2CHO) * c(sH) 
    w(r125f) = k(r125f) * cqss(sqssCH2CHO) * c(sO) 
    w(r126f) = k(r126f) * cqss(sqssCH2CHO) * c(sOH) 
    w(r127f) = k(r127f) * cqss(sqssCH2CHO) * c(sO2) 
    w(r128f) = k(r128f) * cqss(sqssCH2CHO) * c(sCH3) 
    w(r129f) = k(r129f) * cqss(sqssCH2CHO) * c(sHO2) 
    w(r130f) = k(r130f) * cqss(sqssCH2CHO) 
    w(r131f) = k(r131f) * c(sC2H2) * c(sCH3) 
    w(r132f) = k(r132f) * c(sC3H6) * c(sOH) 
    w(r133f) = k(r133f) * c(sC3H6) * c(sO) 
    w(r134f) = k(r134f) * c(sC3H6) * c(sH) 
    w(r135f) = k(r135f) * cqss(sqssC3H5) * c(sH) 
    w(r136f) = k(r136f) * cqss(sqssC3H5) * c(sHO2) 
    w(r137f) = k(r137f) * cqss(sqssC3H5) * c(sHO2) 
    w(r138f) = k(r138f) * cqss(sqssC2H3) * c(sCH3) 
    w(r139f) = k(r139f) * c(sC3H6) * c(sH) 
    w(r140f) = k(r140f) * cqss(sqssC2H3) * c(sCH3) 
    w(r141f) = k(r141f) * c(sC3H8) 
    w(r142f) = k(r142f) * c(sC3H8) * c(sO2) 
    w(r143f) = k(r143f) * c(sC3H8) * c(sO2) 
    w(r144f) = k(r144f) * c(sC3H8) * c(sH) 
    w(r145f) = k(r145f) * c(sC3H8) * c(sH) 
    w(r146f) = k(r146f) * c(sC3H8) * c(sO) 
    w(r147f) = k(r147f) * c(sC3H8) * c(sO) 
    w(r148f) = k(r148f) * c(sC3H8) * c(sOH) 
    w(r149f) = k(r149f) * c(sC3H8) * c(sOH) 
    w(r150f) = k(r150f) * c(sC3H8) * c(sHO2) 
    w(r151f) = k(r151f) * c(sC3H8) * c(sHO2) 
    w(r152f) = k(r152f) * c(sC3H8) * cqss(sqssIXC3H7) 
    w(r153f) = k(r153f) * c(sC3H6) * c(sH) 
    w(r154f) = k(r154f) * cqss(sqssIXC3H7) * c(sO2) 
    w(r155f) = k(r155f) * cqss(sqssNXC3H7) 
    w(r156f) = k(r156f) * c(sC3H6) * c(sH) 
    w(r157f) = k(r157f) * cqss(sqssNXC3H7) * c(sO2) 
    w(r158f) = k(r158f) * cqss(sqssNXC3H7) * c(sO2) 
    w(r159f) = k(r159f) * cqss(sqssC3H6OOH) 
    w(r160f) = k(r160f) * cqss(sqssC3H6OOH) * c(sO2) 
    w(r161f) = k(r161f) * c(sOC3H5OOH) 
    w(r1b) = k(r1b) * c(sO) * c(sOH) 
    w(r2b) = k(r2b) * c(sH) * c(sOH) 
    w(r3b) = k(r3b) * c(sH) * c(sH2O) 
    w(r4b) = k(r4b) * c(sOH)**2.0_pr 
    w(r5b) = k(r5b) * c(sH2) * m(mM5)
    w(r6b) = k(r6b) * c(sH2O) * m(mM6)
    w(r7b) = k(r7b) * c(sO2) * m(mM7)
    w(r8b) = k(r8b) * c(sOH) * m(mM8)
    w(r9b) = k(r9b) * c(sHO2) 
    w(r10b) = k(r10b) * c(sOH)**2.0_pr 
    w(r11b) = k(r11b) * c(sH2) * c(sO2) 
    w(r12b) = k(r12b) * c(sH2O) * c(sO) 
    w(r13b) = k(r13b) * c(sO2) * c(sOH) 
    w(r14b) = k(r14b) * c(sH2O) * c(sO2) 
    w(r15b) = k(r15b) * c(sH2O) * c(sO2) 
    w(r16b) = k(r16b) * c(sH2O2) 
    w(r17b) = k(r17b) * c(sH2O2) * c(sO2) 
    w(r18b) = k(r18b) * c(sH2O2) * c(sO2) 
    w(r19b) = k(r19b) * c(sH2) * c(sHO2) 
    w(r20b) = k(r20b) * c(sH2O) * c(sOH) 
    w(r21b) = k(r21b) * c(sH2O) * c(sHO2) 
    w(r22b) = k(r22b) * c(sH2O) * c(sHO2) 
    w(r23b) = k(r23b) * c(sHO2) * c(sOH) 
    w(r24b) = k(r24b) * c(sCO2) 
    w(r25b) = k(r25b) * c(sCO2) * c(sH) 
    w(r26b) = k(r26b) * c(sCO2) * c(sOH) 
    w(r27b) = k(r27b) * c(sCO2) * c(sO) 
    w(r28b) = k(r28b) * c(sCO) * c(sH) * m(mM28)
    w(r29b) = k(r29b) * c(sCO) * c(sH2) 
    w(r30b) = k(r30b) * c(sCO) * c(sOH) 
    w(r31b) = k(r31b) * c(sCO2) * c(sH) 
    w(r32b) = k(r32b) * c(sCO) * c(sH2O) 
    w(r33b) = k(r33b) * c(sCO) * c(sHO2) 
    w(r34b) = k(r34b) * c(sCH4) * c(sCO) 
    w(r35b) = k(r35b) * c(sCH2O) 
    w(r36b) = k(r36b) * c(sH2) * cqss(sqssHCO) 
    w(r37b) = k(r37b) * cqss(sqssHCO) * c(sOH) 
    w(r38b) = k(r38b) * c(sH2O) * cqss(sqssHCO) 
    w(r39b) = k(r39b) * cqss(sqssHCO) * c(sHO2) 
    w(r40b) = k(r40b) * c(sH2O2) * cqss(sqssHCO) 
    w(r41b) = k(r41b) * c(sCH3) * c(sH2) 
    w(r42b) = k(r42b) * c(sCH3) * c(sH2O) 
    w(r43b) = k(r43b) * c(sCH3) * c(sOH) 
    w(r44b) = k(r44b) * c(sCH3) * c(sHO2) 
    w(r45b) = k(r45b) * c(sCH3) * c(sH2O2) 
    w(r46b) = k(r46b) * c(sH2) * cqss(sqssTXCH2) 
    w(r47b) = k(r47b) * c(sH2) * cqss(sqssSXCH2) 
    w(r48b) = k(r48b) * c(sH2O) * cqss(sqssSXCH2) 
    w(r49b) = k(r49b) * c(sCH2O) * c(sH) 
    w(r50b) = k(r50b) * c(sC2H4) * c(sH) 
    w(r51b) = k(r51b) * cqss(sqssCH3O) * c(sOH) 
    w(r52b) = k(r52b) * c(sCH2O) * c(sOH) 
    w(r53b) = k(r53b) * cqss(sqssCH3O) * c(sO) 
    w(r54b) = k(r54b) * c(sC2H4) * c(sH2) 
    w(r55b) = k(r55b) * cqss(sqssC2H5) * c(sH) 
    w(r56b) = k(r56b) * c(sCH4) 
    w(r57b) = k(r57b) * c(sC2H6) 
    w(r58b) = k(r58b) * c(sCH2O) * c(sH) 
    w(r59b) = k(r59b) * c(sCO) * c(sH) * c(sOH) 
    w(r60b) = k(r60b) * c(sCH2O) * c(sCO) 
    w(r61b) = k(r61b) * cqss(sqssTXCH2) * m(mM61)
    w(r62b) = k(r62b) * cqss(sqssCH) * c(sH2) 
    w(r63b) = k(r63b) * c(sCH2O) * c(sH) 
    w(r64b) = k(r64b) * cqss(sqssCH) * c(sH2O) 
    w(r65b) = k(r65b) * c(sCO) * c(sH) 
    w(r66b) = k(r66b) * c(sCO) * c(sH2) 
    w(r67b) = k(r67b) * c(sCO2) * c(sH2) 
    w(r68b) = k(r68b) * c(sCO) * c(sH) * c(sOH) 
    w(r69b) = k(r69b) * c(sCO) * c(sH) 
    w(r70b) = k(r70b) * cqss(sqssHCO) * c(sO) 
    w(r71b) = k(r71b) * c(sCH2O) * c(sH) 
    w(r72b) = k(r72b) * c(sCO) * cqss(sqssHCO) 
    w(r73b) = k(r73b) * c(sCH2O) * c(sH2) 
    w(r74b) = k(r74b) * c(sH2O) * cqss(sqssSXCH2) 
    w(r75b) = k(r75b) * c(sCH2O) * c(sH2O) 
    w(r76b) = k(r76b) * c(sCH2O) * c(sOH) 
    w(r77b) = k(r77b) * c(sCH2O) * c(sHO2) 
    w(r78b) = k(r78b) * c(sCH2O) * c(sH) * m(mM78)
    w(r79b) = k(r79b) * cqss(sqssC2H5) * c(sH2) 
    w(r80b) = k(r80b) * cqss(sqssC2H5) * c(sOH) 
    w(r81b) = k(r81b) * cqss(sqssC2H5) * c(sH2O) 
    w(r82b) = k(r82b) * cqss(sqssC2H5) * c(sCH4) 
    w(r83b) = k(r83b) * cqss(sqssC2H5) * c(sH) 
    w(r84b) = k(r84b) * cqss(sqssC2H5) * c(sH2O2) 
    w(r85b) = k(r85b) * c(sC2H4) * c(sH2) 
    w(r86b) = k(r86b) * c(sC2H4) * c(sOH) 
    w(r87b) = k(r87b) * c(sCH2O) * c(sCH3) 
    w(r88b) = k(r88b) * c(sC2H4) * c(sHO2) 
    w(r89b) = k(r89b) * c(sC2H4) * c(sH) 
    w(r90b) = k(r90b) * cqss(sqssC2H3) * c(sH2) 
    w(r91b) = k(r91b) * cqss(sqssC2H3) * c(sH2O) 
    w(r92b) = k(r92b) * c(sCH3) * cqss(sqssHCO) 
    w(r93b) = k(r93b) * cqss(sqssCH2CHO) * c(sH) 
    w(r94b) = k(r94b) * cqss(sqssC2H3) * c(sHO2) 
    w(r95b) = k(r95b) * cqss(sqssC2H3) * c(sH) * m(mM95)
    w(r96b) = k(r96b) * c(sC2H2) * c(sH2) * m(mM96)
    w(r97b) = k(r97b) * c(sC2H2) * c(sH2) 
    w(r98b) = k(r98b) * c(sC2H2) * c(sH) 
    w(r99b) = k(r99b) * c(sCH2O) * cqss(sqssHCO) 
    w(r100b) = k(r100b) * cqss(sqssCH2CHO) * c(sO) 
    w(r101b) = k(r101b) * c(sC2H2) * c(sHO2) 
    w(r102b) = k(r102b) * c(sH) * cqss(sqssHCCO) 
    w(r103b) = k(r103b) * c(sCO) * cqss(sqssTXCH2) 
    w(r104b) = k(r104b) * c(sCH2O) * c(sCO) 
    w(r105b) = k(r105b) * c(sCH2CO) * c(sH) 
    w(r106b) = k(r106b) * c(sCH3) * c(sCO) 
    w(r107b) = k(r107b) * c(sCO2) * cqss(sqssTXCH2) 
    w(r108b) = k(r108b) * cqss(sqssHCCO) * c(sOH) 
    w(r109b) = k(r109b) * cqss(sqssC2H5) * c(sCO) 
    w(r110b) = k(r110b) * c(sCO) * cqss(sqssSXCH2) 
    w(r111b) = k(r111b) * c(sCO) * c(sH) * cqss(sqssHCO) 
    w(r112b) = k(r112b) * c(sCO) * c(sH) 
    w(r113b) = k(r113b) * c(sCO) * c(sOH) 
    w(r114b) = k(r114b) * c(sCO) * c(sCO2) * c(sH) 
    w(r115b) = k(r115b) * c(sCH2O) * c(sH2) 
    w(r116b) = k(r116b) * c(sCH3) * c(sOH) 
    w(r117b) = k(r117b) * c(sCH2O) * c(sH2O) 
    w(r118b) = k(r118b) * c(sCH2O) * c(sHO2) 
    w(r119b) = k(r119b) * c(sCH2O) * c(sH) * m(mM119)
    w(r120b) = k(r120b) * cqss(sqssCH2OH) * m(mM120)
    w(r121b) = k(r121b) * cqss(sqssCH2OH) * c(sCO) 
    w(r122b) = k(r122b) * c(sCH2CO) * c(sH) 
    w(r123b) = k(r123b) * c(sCH3) * cqss(sqssHCO) 
    w(r124b) = k(r124b) * c(sCH2CO) * c(sH2) 
    w(r125b) = k(r125b) * c(sCH2O) * cqss(sqssHCO) 
    w(r126b) = k(r126b) * c(sCH2CO) * c(sH2O) 
    w(r127b) = k(r127b) * c(sCH2O) * c(sCO) * c(sOH) 
    w(r128b) = k(r128b) * cqss(sqssC2H5) * c(sCO) * c(sH) 
    w(r129b) = k(r129b) * c(sCH2O) * cqss(sqssHCO) * c(sOH) 
    w(r130b) = k(r130b) * c(sCH3) * c(sCO) 
    w(r131b) = k(r131b) * cqss(sqssC3H5) 
    w(r132b) = k(r132b) * cqss(sqssC3H5) * c(sH2O) 
    w(r133b) = k(r133b) * c(sCH2CO) * c(sCH3) * c(sH) 
    w(r134b) = k(r134b) * cqss(sqssC3H5) * c(sH2) 
    w(r135b) = k(r135b) * c(sC3H6) 
    w(r136b) = k(r136b) * c(sC3H6) * c(sO2) 
    w(r137b) = k(r137b) * cqss(sqssC2H3) * c(sCH2O) * c(sOH) 
    w(r138b) = k(r138b) * c(sC3H6) 
    w(r139b) = k(r139b) * c(sC2H4) * c(sCH3) 
    w(r140b) = k(r140b) * cqss(sqssC3H5) * c(sH) 
    w(r141b) = k(r141b) * cqss(sqssC2H5) * c(sCH3) 
    w(r142b) = k(r142b) * c(sHO2) * cqss(sqssIXC3H7) 
    w(r143b) = k(r143b) * c(sHO2) * cqss(sqssNXC3H7) 
    w(r144b) = k(r144b) * c(sH2) * cqss(sqssIXC3H7) 
    w(r145b) = k(r145b) * c(sH2) * cqss(sqssNXC3H7) 
    w(r146b) = k(r146b) * cqss(sqssIXC3H7) * c(sOH) 
    w(r147b) = k(r147b) * cqss(sqssNXC3H7) * c(sOH) 
    w(r148b) = k(r148b) * c(sH2O) * cqss(sqssNXC3H7) 
    w(r149b) = k(r149b) * c(sH2O) * cqss(sqssIXC3H7) 
    w(r150b) = k(r150b) * c(sH2O2) * cqss(sqssIXC3H7) 
    w(r151b) = k(r151b) * c(sH2O2) * cqss(sqssNXC3H7) 
    w(r152b) = k(r152b) * c(sC3H8) * cqss(sqssNXC3H7) 
    w(r153b) = k(r153b) * cqss(sqssIXC3H7) 
    w(r154b) = k(r154b) * c(sC3H6) * c(sHO2) 
    w(r155b) = k(r155b) * c(sC2H4) * c(sCH3) 
    w(r156b) = k(r156b) * cqss(sqssNXC3H7) 
    w(r157b) = k(r157b) * c(sC3H6) * c(sHO2) 
    w(r158b) = k(r158b) * cqss(sqssC3H6OOH) 
    w(r159b) = k(r159b) * c(sC3H6) * c(sHO2) 
    w(r160b) = k(r160b) * c(sOC3H5OOH) * c(sOH) 
    w(r161b) = k(r161b) * cqss(sqssCH2CHO) * c(sCH2O) * c(sOH) 
  
    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(sH) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r3f) &
         - w(r3b) &
         -2.0_pr * w(r5f) &
         + 2.0_pr * w(r5b) &
         - w(r6f) &
         + w(r6b) &
         - w(r8f) &
         + w(r8b) &
         - w(r9f) &
         + w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         - w(r12f) &
         + w(r12b) &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         + w(r25f) &
         - w(r25b) &
         + w(r28f) &
         - w(r28b) &
         - w(r29f) &
         + w(r29b) &
         + w(r31f) &
         - w(r31b) &
         - w(r35f) &
         + w(r35b) &
         - w(r36f) &
         + w(r36b) &
         - w(r41f) &
         + w(r41b) &
         - w(r46f) &
         + w(r46b) &
         - w(r47f) &
         + w(r47b) &
         + w(r49f) &
         - w(r49b) &
         + w(r50f) &
         - w(r50b) &
         + w(r55f) &
         - w(r55b) &
         - w(r56f) &
         + w(r56b) &
         + w(r58f) &
         - w(r58b) &
         + w(r59f) &
         - w(r59b) &
         - w(r62f) &
         + w(r62b) &
         + w(r63f) &
         - w(r63b) &
         + 2.0_pr * w(r65f) &
         - 2.0_pr * w(r65b) &
         + w(r68f) &
         - w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r71f) &
         - w(r71b) &
         - w(r73f) &
         + w(r73b) &
         - w(r74f) &
         + w(r74b) &
         + w(r78f) &
         - w(r78b) &
         - w(r79f) &
         + w(r79b) &
         + w(r83f) &
         - w(r83b) &
         - w(r85f) &
         + w(r85b) &
         + w(r89f) &
         - w(r89b) &
         - w(r90f) &
         + w(r90b) &
         + w(r93f) &
         - w(r93b) &
         + w(r95f) &
         - w(r95b) &
         - w(r97f) &
         + w(r97b) &
         + w(r98f) &
         - w(r98b) &
         + w(r102f) &
         - w(r102b) &
         + w(r105f) &
         - w(r105b) &
         - w(r106f) &
         + w(r106b) &
         - w(r110f) &
         + w(r110b) &
         + w(r111f) &
         - w(r111b) &
         + w(r112f) &
         - w(r112b) &
         + w(r114f) &
         - w(r114b) &
         - w(r115f) &
         + w(r115b) &
         - w(r116f) &
         + w(r116b) &
         + w(r119f) &
         - w(r119b) &
         + w(r122f) &
         - w(r122b) &
         - w(r123f) &
         + w(r123b) &
         - w(r124f) &
         + w(r124b) &
         + w(r128f) &
         - w(r128b) &
         + w(r133f) &
         - w(r133b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r139f) &
         + w(r139b) &
         + w(r140f) &
         - w(r140b) &
         - w(r144f) &
         + w(r144b) &
         - w(r145f) &
         + w(r145b) &
         - w(r153f) &
         + w(r153b) &
         - w(r156f) &
         + w(r156b) 

    cdot(sO2) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         + w(r7f) &
         - w(r7b) &
         - w(r9f) &
         + w(r9b) &
         + w(r11f) &
         - w(r11b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         - w(r27f) &
         + w(r27b) &
         - w(r33f) &
         + w(r33b) &
         - w(r39f) &
         + w(r39b) &
         - w(r44f) &
         + w(r44b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         - w(r59f) &
         + w(r59b) &
         - w(r67f) &
         + w(r67b) &
         - w(r68f) &
         + w(r68b) &
         - w(r70f) &
         + w(r70b) &
         - w(r77f) &
         + w(r77b) &
         - w(r88f) &
         + w(r88b) &
         - w(r94f) &
         + w(r94b) &
         - w(r99f) &
         + w(r99b) &
         - w(r100f) &
         + w(r100b) &
         - w(r101f) &
         + w(r101b) &
         - w(r104f) &
         + w(r104b) &
         - w(r113f) &
         + w(r113b) &
         - w(r114f) &
         + w(r114b) &
         - w(r118f) &
         + w(r118b) &
         - w(r127f) &
         + w(r127b) &
         + w(r136f) &
         - w(r136b) &
         - w(r142f) &
         + w(r142b) &
         - w(r143f) &
         + w(r143b) &
         - w(r154f) &
         + w(r154b) &
         - w(r157f) &
         + w(r157b) &
         - w(r158f) &
         + w(r158b) &
         - w(r160f) &
         + w(r160b) 

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         - w(r3f) &
         + w(r3b) &
         + 2.0_pr * w(r4f) &
         - 2.0_pr * w(r4b) &
         - w(r6f) &
         + w(r6b) &
         + w(r8f) &
         - w(r8b) &
         + 2.0_pr * w(r10f) &
         - 2.0_pr * w(r10b) &
         + w(r13f) &
         - w(r13b) &
         - w(r14f) &
         + w(r14b) &
         - w(r15f) &
         + w(r15b) &
         -2.0_pr * w(r16f) &
         + 2.0_pr * w(r16b) &
         + w(r20f) &
         - w(r20b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         + w(r23f) &
         - w(r23b) &
         - w(r25f) &
         + w(r25b) &
         + w(r26f) &
         - w(r26b) &
         + w(r30f) &
         - w(r30b) &
         - w(r32f) &
         + w(r32b) &
         + w(r37f) &
         - w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r42f) &
         + w(r42b) &
         + w(r43f) &
         - w(r43b) &
         - w(r48f) &
         + w(r48b) &
         + w(r51f) &
         - w(r51b) &
         + w(r52f) &
         - w(r52b) &
         - w(r58f) &
         + w(r58b) &
         + w(r59f) &
         - w(r59b) &
         - w(r63f) &
         + w(r63b) &
         - w(r64f) &
         + w(r64b) &
         + w(r68f) &
         - w(r68b) &
         - w(r75f) &
         + w(r75b) &
         + w(r76f) &
         - w(r76b) &
         + w(r80f) &
         - w(r80b) &
         - w(r81f) &
         + w(r81b) &
         + w(r86f) &
         - w(r86b) &
         - w(r91f) &
         + w(r91b) &
         - w(r105f) &
         + w(r105b) &
         + w(r108f) &
         - w(r108b) &
         - w(r111f) &
         + w(r111b) &
         + w(r113f) &
         - w(r113b) &
         + w(r116f) &
         - w(r116b) &
         - w(r117f) &
         + w(r117b) &
         - w(r121f) &
         + w(r121b) &
         - w(r126f) &
         + w(r126b) &
         + w(r127f) &
         - w(r127b) &
         + w(r129f) &
         - w(r129b) &
         - w(r132f) &
         + w(r132b) &
         + w(r137f) &
         - w(r137b) &
         + w(r146f) &
         - w(r146b) &
         + w(r147f) &
         - w(r147b) &
         - w(r148f) &
         + w(r148b) &
         - w(r149f) &
         + w(r149b) &
         + w(r160f) &
         - w(r160b) &
         + w(r161f) &
         - w(r161b) 

    cdot(sO) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         - w(r2f) &
         + w(r2b) &
         - w(r4f) &
         + w(r4b) &
         -2.0_pr * w(r7f) &
         + 2.0_pr * w(r7b) &
         - w(r8f) &
         + w(r8b) &
         + w(r12f) &
         - w(r12b) &
         - w(r13f) &
         + w(r13b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) &
         + w(r27f) &
         - w(r27b) &
         - w(r30f) &
         + w(r30b) &
         - w(r31f) &
         + w(r31b) &
         - w(r37f) &
         + w(r37b) &
         - w(r43f) &
         + w(r43b) &
         - w(r49f) &
         + w(r49b) &
         + w(r53f) &
         - w(r53b) &
         - w(r65f) &
         + w(r65b) &
         - w(r66f) &
         + w(r66b) &
         - w(r69f) &
         + w(r69b) &
         + w(r70f) &
         - w(r70b) &
         - w(r76f) &
         + w(r76b) &
         - w(r80f) &
         + w(r80b) &
         - w(r86f) &
         + w(r86b) &
         - w(r87f) &
         + w(r87b) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         + w(r100f) &
         - w(r100b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r107f) &
         + w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r112f) &
         + w(r112b) &
         - w(r125f) &
         + w(r125b) &
         - w(r133f) &
         + w(r133b) &
         - w(r146f) &
         + w(r146b) &
         - w(r147f) &
         + w(r147b) 

    cdot(sH2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r3b) &
         + w(r5f) &
         - w(r5b) &
         + w(r11f) &
         - w(r11b) &
         + w(r19f) &
         - w(r19b) &
         + w(r29f) &
         - w(r29b) &
         + w(r36f) &
         - w(r36b) &
         + w(r41f) &
         - w(r41b) &
         + w(r46f) &
         - w(r46b) &
         + w(r47f) &
         - w(r47b) &
         + w(r54f) &
         - w(r54b) &
         + w(r62f) &
         - w(r62b) &
         + w(r66f) &
         - w(r66b) &
         + w(r67f) &
         - w(r67b) &
         + w(r73f) &
         - w(r73b) &
         + w(r79f) &
         - w(r79b) &
         + w(r85f) &
         - w(r85b) &
         + w(r90f) &
         - w(r90b) &
         + w(r96f) &
         - w(r96b) &
         + w(r97f) &
         - w(r97b) &
         + w(r115f) &
         - w(r115b) &
         + w(r124f) &
         - w(r124b) &
         + w(r134f) &
         - w(r134b) &
         + w(r144f) &
         - w(r144b) &
         + w(r145f) &
         - w(r145b) 

    cdot(sH2O) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         - w(r4f) &
         + w(r4b) &
         + w(r6f) &
         - w(r6b) &
         + w(r12f) &
         - w(r12b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r20f) &
         - w(r20b) &
         + w(r21f) &
         - w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r32f) &
         - w(r32b) &
         + w(r38f) &
         - w(r38b) &
         + w(r42f) &
         - w(r42b) &
         + w(r48f) &
         - w(r48b) &
         + w(r64f) &
         - w(r64b) &
         - w(r71f) &
         + w(r71b) &
         + w(r74f) &
         - w(r74b) &
         + w(r75f) &
         - w(r75b) &
         + w(r81f) &
         - w(r81b) &
         + w(r91f) &
         - w(r91b) &
         + w(r117f) &
         - w(r117b) &
         + w(r126f) &
         - w(r126b) &
         + w(r132f) &
         - w(r132b) &
         + w(r148f) &
         - w(r148b) &
         + w(r149f) &
         - w(r149b) 

    cdot(sHO2) = 0.0_pr &
         + w(r9f) &
         - w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         - w(r12f) &
         + w(r12b) &
         - w(r13f) &
         + w(r13b) &
         - w(r14f) &
         + w(r14b) &
         - w(r15f) &
         + w(r15b) &
         -2.0_pr * w(r17f) &
         + 2.0_pr * w(r17b) &
         -2.0_pr * w(r18f) &
         + 2.0_pr * w(r18b) &
         + w(r19f) &
         - w(r19b) &
         + w(r21f) &
         - w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r23f) &
         - w(r23b) &
         - w(r26f) &
         + w(r26b) &
         + w(r33f) &
         - w(r33b) &
         + w(r39f) &
         - w(r39b) &
         - w(r40f) &
         + w(r40b) &
         + w(r44f) &
         - w(r44b) &
         - w(r45f) &
         + w(r45b) &
         - w(r51f) &
         + w(r51b) &
         + w(r77f) &
         - w(r77b) &
         - w(r84f) &
         + w(r84b) &
         + w(r88f) &
         - w(r88b) &
         + w(r94f) &
         - w(r94b) &
         + w(r101f) &
         - w(r101b) &
         + w(r118f) &
         - w(r118b) &
         - w(r129f) &
         + w(r129b) &
         - w(r136f) &
         + w(r136b) &
         - w(r137f) &
         + w(r137b) &
         + w(r142f) &
         - w(r142b) &
         + w(r143f) &
         - w(r143b) &
         - w(r150f) &
         + w(r150b) &
         - w(r151f) &
         + w(r151b) &
         + w(r154f) &
         - w(r154b) &
         + w(r157f) &
         - w(r157b) &
         + w(r159f) &
         - w(r159b) 

    cdot(sH2O2) = 0.0_pr &
         + w(r16f) &
         - w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         + w(r40f) &
         - w(r40b) &
         + w(r45f) &
         - w(r45b) &
         + w(r84f) &
         - w(r84b) &
         + w(r150f) &
         - w(r150b) &
         + w(r151f) &
         - w(r151b) 

    cdot(sCO) = 0.0_pr &
         - w(r24f) &
         + w(r24b) &
         - w(r25f) &
         + w(r25b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         + w(r28f) &
         - w(r28b) &
         + w(r29f) &
         - w(r29b) &
         + w(r30f) &
         - w(r30b) &
         + w(r32f) &
         - w(r32b) &
         + w(r33f) &
         - w(r33b) &
         + w(r34f) &
         - w(r34b) &
         + w(r59f) &
         - w(r59b) &
         + w(r60f) &
         - w(r60b) &
         + w(r65f) &
         - w(r65b) &
         + w(r66f) &
         - w(r66b) &
         + w(r68f) &
         - w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r72f) &
         - w(r72b) &
         + w(r103f) &
         - w(r103b) &
         + w(r104f) &
         - w(r104b) &
         + w(r106f) &
         - w(r106b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b) &
         + w(r111f) &
         - w(r111b) &
         + 2.0_pr * w(r112f) &
         - 2.0_pr * w(r112b) &
         + 2.0_pr * w(r113f) &
         - 2.0_pr * w(r113b) &
         + w(r114f) &
         - w(r114b) &
         + w(r121f) &
         - w(r121b) &
         + w(r127f) &
         - w(r127b) &
         + w(r128f) &
         - w(r128b) &
         + w(r130f) &
         - w(r130b) 

    cdot(sCO2) = 0.0_pr &
         + w(r24f) &
         - w(r24b) &
         + w(r25f) &
         - w(r25b) &
         + w(r26f) &
         - w(r26b) &
         + w(r27f) &
         - w(r27b) &
         + w(r31f) &
         - w(r31b) &
         - w(r60f) &
         + w(r60b) &
         + w(r67f) &
         - w(r67b) &
         - w(r72f) &
         + w(r72b) &
         + w(r107f) &
         - w(r107b) &
         + w(r114f) &
         - w(r114b) 

    cdot(sCH3) = 0.0_pr &
         - w(r34f) &
         + w(r34b) &
         + w(r41f) &
         - w(r41b) &
         + w(r42f) &
         - w(r42b) &
         + w(r43f) &
         - w(r43b) &
         + w(r44f) &
         - w(r44b) &
         + w(r45f) &
         - w(r45b) &
         - w(r46f) &
         + w(r46b) &
         - w(r47f) &
         + w(r47b) &
         - w(r48f) &
         + w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         -2.0_pr * w(r54f) &
         + 2.0_pr * w(r54b) &
         -2.0_pr * w(r55f) &
         + 2.0_pr * w(r55b) &
         - w(r56f) &
         + w(r56b) &
         -2.0_pr * w(r57f) &
         + 2.0_pr * w(r57b) &
         - w(r82f) &
         + w(r82b) &
         + w(r87f) &
         - w(r87b) &
         + w(r92f) &
         - w(r92b) &
         + w(r106f) &
         - w(r106b) &
         - w(r109f) &
         + w(r109b) &
         + w(r116f) &
         - w(r116b) &
         + w(r123f) &
         - w(r123b) &
         - w(r128f) &
         + w(r128b) &
         + w(r130f) &
         - w(r130b) &
         - w(r131f) &
         + w(r131b) &
         + w(r133f) &
         - w(r133b) &
         - w(r138f) &
         + w(r138b) &
         + w(r139f) &
         - w(r139b) &
         - w(r140f) &
         + w(r140b) &
         + w(r141f) &
         - w(r141b) &
         + w(r155f) &
         - w(r155b) 

    cdot(sCH4) = 0.0_pr &
         + w(r34f) &
         - w(r34b) &
         - w(r41f) &
         + w(r41b) &
         - w(r42f) &
         + w(r42b) &
         - w(r43f) &
         + w(r43b) &
         - w(r44f) &
         + w(r44b) &
         - w(r45f) &
         + w(r45b) &
         + w(r56f) &
         - w(r56b) &
         + w(r82f) &
         - w(r82b) 

    cdot(sCH2O) = 0.0_pr &
         + w(r35f) &
         - w(r35b) &
         - w(r36f) &
         + w(r36b) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r39f) &
         + w(r39b) &
         - w(r40f) &
         + w(r40b) &
         + w(r49f) &
         - w(r49b) &
         + w(r52f) &
         - w(r52b) &
         + w(r58f) &
         - w(r58b) &
         + w(r60f) &
         - w(r60b) &
         + w(r63f) &
         - w(r63b) &
         + w(r71f) &
         - w(r71b) &
         + w(r73f) &
         - w(r73b) &
         + w(r75f) &
         - w(r75b) &
         + w(r76f) &
         - w(r76b) &
         + w(r77f) &
         - w(r77b) &
         + w(r78f) &
         - w(r78b) &
         + w(r87f) &
         - w(r87b) &
         + w(r99f) &
         - w(r99b) &
         + w(r104f) &
         - w(r104b) &
         + w(r115f) &
         - w(r115b) &
         + w(r117f) &
         - w(r117b) &
         + w(r118f) &
         - w(r118b) &
         + w(r119f) &
         - w(r119b) &
         + w(r125f) &
         - w(r125b) &
         + w(r127f) &
         - w(r127b) &
         + w(r129f) &
         - w(r129b) &
         + w(r137f) &
         - w(r137b) &
         + w(r161f) &
         - w(r161b) 

    cdot(sC2H4) = 0.0_pr &
         + w(r50f) &
         - w(r50b) &
         + w(r54f) &
         - w(r54b) &
         + w(r85f) &
         - w(r85b) &
         + w(r86f) &
         - w(r86b) &
         + w(r88f) &
         - w(r88b) &
         + w(r89f) &
         - w(r89b) &
         - w(r90f) &
         + w(r90b) &
         - w(r91f) &
         + w(r91b) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         - w(r96f) &
         + w(r96b) &
         + w(r139f) &
         - w(r139b) &
         + w(r155f) &
         - w(r155b) 

    cdot(sC2H6) = 0.0_pr &
         + w(r57f) &
         - w(r57b) &
         - w(r79f) &
         + w(r79b) &
         - w(r80f) &
         + w(r80b) &
         - w(r81f) &
         + w(r81b) &
         - w(r82f) &
         + w(r82b) &
         - w(r83f) &
         + w(r83b) &
         - w(r84f) &
         + w(r84b) 

    cdot(sC2H2) = 0.0_pr &
         + w(r96f) &
         - w(r96b) &
         + w(r97f) &
         - w(r97b) &
         + w(r98f) &
         - w(r98b) &
         + w(r101f) &
         - w(r101b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         - w(r105f) &
         + w(r105b) &
         - w(r131f) &
         + w(r131b) 

    cdot(sCH2CO) = 0.0_pr &
         + w(r105f) &
         - w(r105b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         + w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r109f) &
         + w(r109b) &
         - w(r121f) &
         + w(r121b) &
         + w(r122f) &
         - w(r122b) &
         + w(r124f) &
         - w(r124b) &
         + w(r126f) &
         - w(r126b) &
         + w(r133f) &
         - w(r133b) 

    cdot(sC3H6) = 0.0_pr &
         - w(r132f) &
         + w(r132b) &
         - w(r133f) &
         + w(r133b) &
         - w(r134f) &
         + w(r134b) &
         + w(r135f) &
         - w(r135b) &
         + w(r136f) &
         - w(r136b) &
         + w(r138f) &
         - w(r138b) &
         - w(r139f) &
         + w(r139b) &
         - w(r153f) &
         + w(r153b) &
         + w(r154f) &
         - w(r154b) &
         - w(r156f) &
         + w(r156b) &
         + w(r157f) &
         - w(r157b) &
         + w(r159f) &
         - w(r159b) 

    cdot(sC3H8) = 0.0_pr &
         - w(r141f) &
         + w(r141b) &
         - w(r142f) &
         + w(r142b) &
         - w(r143f) &
         + w(r143b) &
         - w(r144f) &
         + w(r144b) &
         - w(r145f) &
         + w(r145b) &
         - w(r146f) &
         + w(r146b) &
         - w(r147f) &
         + w(r147b) &
         - w(r148f) &
         + w(r148b) &
         - w(r149f) &
         + w(r149b) &
         - w(r150f) &
         + w(r150b) &
         - w(r151f) &
         + w(r151b) 

    cdot(sOC3H5OOH) = 0.0_pr &
         + w(r160f) &
         - w(r160b) &
         - w(r161f) &
         + w(r161b) 

  
    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) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_CH
    real(pr) :: HCO_C2H3
    real(pr) :: HCO_CH2CHO
    real(pr) :: HCO_HCCO
    real(pr) :: TXCH2_ct
    real(pr) :: TXCH2_num
    real(pr) :: TXCH2_denom
    real(pr) :: TXCH2_SXCH2
    real(pr) :: TXCH2_CH
    real(pr) :: SXCH2_ct
    real(pr) :: SXCH2_num
    real(pr) :: SXCH2_denom
    real(pr) :: SXCH2_TXCH2
    real(pr) :: SXCH2_CH3O
    real(pr) :: SXCH2_HCCO
    real(pr) :: CH3O_ct
    real(pr) :: CH3O_num
    real(pr) :: CH3O_denom
    real(pr) :: CH3O_SXCH2
    real(pr) :: CH3O_CH2OH
    real(pr) :: C2H5_ct
    real(pr) :: C2H5_num
    real(pr) :: C2H5_denom
    real(pr) :: C2H5_CH2CHO
    real(pr) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_HCO
    real(pr) :: CH_TXCH2
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_HCO
    real(pr) :: C2H3_CH2CHO
    real(pr) :: C2H3_C3H5
    real(pr) :: CH2CHO_ct
    real(pr) :: CH2CHO_num
    real(pr) :: CH2CHO_denom
    real(pr) :: CH2CHO_HCO
    real(pr) :: CH2CHO_C2H5
    real(pr) :: CH2CHO_C2H3
    real(pr) :: HCCO_ct
    real(pr) :: HCCO_num
    real(pr) :: HCCO_denom
    real(pr) :: HCCO_HCO
    real(pr) :: HCCO_SXCH2
    real(pr) :: CH2OH_ct
    real(pr) :: CH2OH_num
    real(pr) :: CH2OH_denom
    real(pr) :: CH2OH_CH3O
    real(pr) :: C3H5_ct
    real(pr) :: C3H5_num
    real(pr) :: C3H5_denom
    real(pr) :: C3H5_C2H3
    real(pr) :: IXC3H7_ct
    real(pr) :: IXC3H7_num
    real(pr) :: IXC3H7_denom
    real(pr) :: IXC3H7_NXC3H7
    real(pr) :: NXC3H7_ct
    real(pr) :: NXC3H7_num
    real(pr) :: NXC3H7_denom
    real(pr) :: NXC3H7_IXC3H7
    real(pr) :: NXC3H7_C3H6OOH
    real(pr) :: C3H6OOH_ct
    real(pr) :: C3H6OOH_num
    real(pr) :: C3H6OOH_denom
    real(pr) :: C3H6OOH_NXC3H7
  
    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r28f)* m(mM28) &
             + k(r29f)* c(sH) &
             + k(r30f)* c(sO) &
             + k(r31f)* c(sO) &
             + k(r32f)* c(sOH) &
             + k(r33f)* c(sO2) &
             + k(r34f)* c(sCH3) &
             + k(r35f)* c(sH) &
             + k(r36b)* c(sH2) &
             + k(r37b)* c(sOH) &
             + k(r38b)* c(sH2O) &
             + k(r39b)* c(sHO2) &
             + k(r40b)* c(sH2O2) &
             + k(r70b)* c(sO) &
             + k(r72b)* c(sCO) &
             + k(r92b)* c(sCH3) &
             + k(r99b)* c(sCH2O) &
             + k(r111b)* c(sCO) * c(sH) &
             + k(r123b)* c(sCH3) &
             + k(r125b)* c(sCH2O) &
             + k(r129b)* c(sCH2O) * c(sOH)   )

    HCO_num = ( 0.0_pr &
             + k(r36f)* c(sCH2O) * c(sH) &
             + k(r37f)* c(sCH2O) * c(sO) &
             + k(r38f)* c(sCH2O) * c(sOH) &
             + k(r39f)* c(sCH2O) * c(sO2) &
             + k(r40f)* c(sCH2O) * c(sHO2) &
             + k(r92f)* c(sC2H4) * c(sO) &
             + k(r28b)* c(sCO) * c(sH) * m(mM28) &
             + k(r29b)* c(sCO) * c(sH2) &
             + k(r30b)* c(sCO) * c(sOH) &
             + k(r31b)* c(sCO2) * c(sH) &
             + k(r32b)* c(sCO) * c(sH2O) &
             + k(r33b)* c(sCO) * c(sHO2) &
             + k(r34b)* c(sCH4) * c(sCO) &
             + k(r35b)* c(sCH2O)   )

    HCO_ct = HCO_num / HCO_denom

    HCO_CH = - ( 0.0_pr &
          + k(r70f) * c(sO2) &
          + k(r72f) * c(sCO2)  ) / HCO_denom

    HCO_C2H3 = - ( 0.0_pr &
          + k(r99f) * c(sO2)  ) / HCO_denom

    HCO_CH2CHO = - ( 0.0_pr &
          + k(r123f) * c(sH) &
          + k(r125f) * c(sO) &
          + k(r129f) * c(sHO2)  ) / HCO_denom

    HCO_HCCO = - ( 0.0_pr &
          + k(r111f) * c(sOH)  ) / HCO_denom

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r97f)* c(sH) &
             + k(r98f)&
             + k(r99f)* c(sO2) &
             + k(r100f)* c(sO2) &
             + k(r101f)* c(sO2) &
             + k(r138f)* c(sCH3) &
             + k(r140f)* c(sCH3) &
             + k(r90b)* c(sH2) &
             + k(r91b)* c(sH2O) &
             + k(r94b)* c(sHO2) &
             + k(r95b)* c(sH) * m(mM95) &
             + k(r137b)* c(sCH2O) * c(sOH)   )

    C2H3_num = ( 0.0_pr &
             + k(r90f)* c(sC2H4) * c(sH) &
             + k(r91f)* c(sC2H4) * c(sOH) &
             + k(r94f)* c(sC2H4) * c(sO2) &
             + k(r95f)* c(sC2H4) * m(mM95) &
             + k(r97b)* c(sC2H2) * c(sH2) &
             + k(r98b)* c(sC2H2) * c(sH) &
             + k(r101b)* c(sC2H2) * c(sHO2) &
             + k(r138b)* c(sC3H6)   )

    C2H3_ct = C2H3_num / C2H3_denom

    C2H3_HCO = - ( 0.0_pr &
          + k(r99b) * c(sCH2O)  ) / C2H3_denom

    C2H3_CH2CHO = - ( 0.0_pr &
          + k(r100b) * c(sO)  ) / C2H3_denom

    C2H3_C3H5 = - ( 0.0_pr &
          + k(r137f) * c(sHO2) &
          + k(r140b) * c(sH)  ) / C2H3_denom

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r69f)* c(sO) &
             + k(r70f)* c(sO2) &
             + k(r71f)* c(sH2O) &
             + k(r72f)* c(sCO2) &
             + k(r62b)* c(sH2) &
             + k(r64b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r69b)* c(sCO) * c(sH) &
             + k(r71b)* c(sCH2O) * c(sH)   )

    CH_ct = CH_num / CH_denom

    CH_HCO = - ( 0.0_pr &
          + k(r70b) * c(sO) &
          + k(r72b) * c(sCO)  ) / CH_denom

    CH_TXCH2 = - ( 0.0_pr &
          + k(r62f) * c(sH) &
          + k(r64f) * c(sOH)  ) / CH_denom

    CH2CHO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r122f)&
             + k(r123f)* c(sH) &
             + k(r124f)* c(sH) &
             + k(r125f)* c(sO) &
             + k(r126f)* c(sOH) &
             + k(r127f)* c(sO2) &
             + k(r128f)* c(sCH3) &
             + k(r129f)* c(sHO2) &
             + k(r130f)&
             + k(r93b)* c(sH) &
             + k(r100b)* c(sO) &
             + k(r161b)* c(sCH2O) * c(sOH)   )

    CH2CHO_num = ( 0.0_pr &
             + k(r93f)* c(sC2H4) * c(sO) &
             + k(r161f)* c(sOC3H5OOH) &
             + k(r122b)* c(sCH2CO) * c(sH) &
             + k(r124b)* c(sCH2CO) * c(sH2) &
             + k(r126b)* c(sCH2CO) * c(sH2O) &
             + k(r127b)* c(sCH2O) * c(sCO) * c(sOH) &
             + k(r130b)* c(sCH3) * c(sCO)   )

    CH2CHO_ct = CH2CHO_num / CH2CHO_denom

    CH2CHO_HCO = - ( 0.0_pr &
          + k(r123b) * c(sCH3) &
          + k(r125b) * c(sCH2O) &
          + k(r129b) * c(sCH2O) * c(sOH)  ) / CH2CHO_denom

    CH2CHO_C2H5 = - ( 0.0_pr &
          + k(r128b) * c(sCO) * c(sH)  ) / CH2CHO_denom

    CH2CHO_C2H3 = - ( 0.0_pr &
          + k(r100f) * c(sO2)  ) / CH2CHO_denom

    TXCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r50f)* c(sCH3) &
             + k(r62f)* c(sH) &
             + k(r63f)* c(sOH) &
             + k(r64f)* c(sOH) &
             + k(r65f)* c(sO) &
             + k(r66f)* c(sO) &
             + k(r67f)* c(sO2) &
             + k(r68f)* c(sO2) &
             + k(r46b)* c(sH2) &
             + k(r61b)* m(mM61) &
             + k(r103b)* c(sCO) &
             + k(r107b)* c(sCO2)   )

    TXCH2_num = ( 0.0_pr &
             + k(r46f)* c(sCH3) * c(sH) &
             + k(r103f)* c(sC2H2) * c(sO) &
             + k(r107f)* c(sCH2CO) * c(sO) &
             + k(r50b)* c(sC2H4) * c(sH) &
             + k(r63b)* c(sCH2O) * c(sH) &
             + k(r65b)* c(sCO)  *c(sH)** 2.0_pr&
             + k(r66b)* c(sCO) * c(sH2) &
             + k(r67b)* c(sCO2) * c(sH2) &
             + k(r68b)* c(sCO) * c(sH) * c(sOH)   )

    TXCH2_ct = TXCH2_num / TXCH2_denom

    TXCH2_SXCH2 = - ( 0.0_pr &
          + k(r61f) * m(mM61)  ) / TXCH2_denom

    TXCH2_CH = - ( 0.0_pr &
          + k(r62b) * c(sH2) &
          + k(r64b) * c(sH2O)  ) / TXCH2_denom

    SXCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r58f)* c(sOH) &
             + k(r59f)* c(sO2) &
             + k(r60f)* c(sCO2) &
             + k(r61f)* m(mM61) &
             + k(r47b)* c(sH2) &
             + k(r48b)* c(sH2O) &
             + k(r74b)* c(sH2O) &
             + k(r110b)* c(sCO)   )

    SXCH2_num = ( 0.0_pr &
             + k(r47f)* c(sCH3) * c(sH) &
             + k(r48f)* c(sCH3) * c(sOH) &
             + k(r58b)* c(sCH2O) * c(sH) &
             + k(r59b)* c(sCO) * c(sH) * c(sOH) &
             + k(r60b)* c(sCH2O) * c(sCO)   )

    SXCH2_ct = SXCH2_num / SXCH2_denom

    SXCH2_TXCH2 = - ( 0.0_pr &
          + k(r61b) * m(mM61)  ) / SXCH2_denom

    SXCH2_CH3O = - ( 0.0_pr &
          + k(r74f) * c(sH)  ) / SXCH2_denom

    SXCH2_HCCO = - ( 0.0_pr &
          + k(r110f) * c(sH)  ) / SXCH2_denom

    C2H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r85f)* c(sH) &
             + k(r86f)* c(sO) &
             + k(r87f)* c(sO) &
             + k(r88f)* c(sO2) &
             + k(r89f)&
             + k(r55b)* c(sH) &
             + k(r79b)* c(sH2) &
             + k(r80b)* c(sOH) &
             + k(r81b)* c(sH2O) &
             + k(r82b)* c(sCH4) &
             + k(r83b)* c(sH) &
             + k(r84b)* c(sH2O2) &
             + k(r109b)* c(sCO) &
             + k(r128b)* c(sCO) * c(sH) &
             + k(r141b)* c(sCH3)   )

    C2H5_num = ( 0.0_pr &
             + k(r55f) *c(sCH3)** 2.0_pr&
             + k(r79f)* c(sC2H6) * c(sH) &
             + k(r80f)* c(sC2H6) * c(sO) &
             + k(r81f)* c(sC2H6) * c(sOH) &
             + k(r82f)* c(sC2H6) * c(sCH3) &
             + k(r83f)* c(sC2H6) &
             + k(r84f)* c(sC2H6) * c(sHO2) &
             + k(r109f)* c(sCH2CO) * c(sCH3) &
             + k(r141f)* c(sC3H8) &
             + k(r85b)* c(sC2H4) * c(sH2) &
             + k(r86b)* c(sC2H4) * c(sOH) &
             + k(r87b)* c(sCH2O) * c(sCH3) &
             + k(r88b)* c(sC2H4) * c(sHO2) &
             + k(r89b)* c(sC2H4) * c(sH)   )

    C2H5_ct = C2H5_num / C2H5_denom

    C2H5_CH2CHO = - ( 0.0_pr &
          + k(r128f) * c(sCH3)  ) / C2H5_denom

    C3H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r135f)* c(sH) &
             + k(r136f)* c(sHO2) &
             + k(r137f)* c(sHO2) &
             + k(r131b)&
             + k(r132b)* c(sH2O) &
             + k(r134b)* c(sH2) &
             + k(r140b)* c(sH)   )

    C3H5_num = ( 0.0_pr &
             + k(r131f)* c(sC2H2) * c(sCH3) &
             + k(r132f)* c(sC3H6) * c(sOH) &
             + k(r134f)* c(sC3H6) * c(sH) &
             + k(r135b)* c(sC3H6) &
             + k(r136b)* c(sC3H6) * c(sO2)   )

    C3H5_ct = C3H5_num / C3H5_denom

    C3H5_C2H3 = - ( 0.0_pr &
          + k(r140f) * c(sCH3) &
          + k(r137b) * c(sCH2O) * c(sOH)  ) / C3H5_denom

    HCCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r110f)* c(sH) &
             + k(r111f)* c(sOH) &
             + k(r112f)* c(sO) &
             + k(r113f)* c(sO2) &
             + k(r114f)* c(sO2) &
             + k(r102b)* c(sH) &
             + k(r108b)* c(sOH)   )

    HCCO_num = ( 0.0_pr &
             + k(r102f)* c(sC2H2) * c(sO) &
             + k(r108f)* c(sCH2CO) * c(sO) &
             + k(r112b) *c(sCO)** 2.0_pr* c(sH) &
             + k(r113b) *c(sCO)** 2.0_pr* c(sOH) &
             + k(r114b)* c(sCO) * c(sCO2) * c(sH)   )

    HCCO_ct = HCCO_num / HCCO_denom

    HCCO_HCO = - ( 0.0_pr &
          + k(r111b) * c(sCO) * c(sH)  ) / HCCO_denom

    HCCO_SXCH2 = - ( 0.0_pr &
          + k(r110b) * c(sCO)  ) / HCCO_denom

    CH3O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r73f)* c(sH) &
             + k(r74f)* c(sH) &
             + k(r75f)* c(sOH) &
             + k(r76f)* c(sO) &
             + k(r77f)* c(sO2) &
             + k(r78f)* m(mM78) &
             + k(r120f)* m(mM120) &
             + k(r51b)* c(sOH) &
             + k(r53b)* c(sO)   )

    CH3O_num = ( 0.0_pr &
             + k(r51f)* c(sCH3) * c(sHO2) &
             + k(r53f)* c(sCH3) * c(sO2) &
             + k(r73b)* c(sCH2O) * c(sH2) &
             + k(r75b)* c(sCH2O) * c(sH2O) &
             + k(r76b)* c(sCH2O) * c(sOH) &
             + k(r77b)* c(sCH2O) * c(sHO2) &
             + k(r78b)* c(sCH2O) * c(sH) * m(mM78)   )

    CH3O_ct = CH3O_num / CH3O_denom

    CH3O_SXCH2 = - ( 0.0_pr &
          + k(r74b) * c(sH2O)  ) / CH3O_denom

    CH3O_CH2OH = - ( 0.0_pr &
          + k(r120b) * m(mM120)  ) / CH3O_denom

    CH2OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r115f)* c(sH) &
             + k(r116f)* c(sH) &
             + k(r117f)* c(sOH) &
             + k(r118f)* c(sO2) &
             + k(r119f)* m(mM119) &
             + k(r120b)* m(mM120) &
             + k(r121b)* c(sCO)   )

    CH2OH_num = ( 0.0_pr &
             + k(r121f)* c(sCH2CO) * c(sOH) &
             + k(r115b)* c(sCH2O) * c(sH2) &
             + k(r116b)* c(sCH3) * c(sOH) &
             + k(r117b)* c(sCH2O) * c(sH2O) &
             + k(r118b)* c(sCH2O) * c(sHO2) &
             + k(r119b)* c(sCH2O) * c(sH) * m(mM119)   )

    CH2OH_ct = CH2OH_num / CH2OH_denom

    CH2OH_CH3O = - ( 0.0_pr &
          + k(r120f) * m(mM120)  ) / CH2OH_denom

    cqss(sqssHCO) = (  ( ( ( ( ( (HCO_ct) - (HCCO_ct) * (HCO_HCCO)) &
                 - ( ( (SXCH2_ct) - ( (CH3O_ct) - (CH2OH_ct) &
                 * (CH3O_CH2OH)) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_ct) &
                 * (SXCH2_HCCO)) * ( - (HCCO_SXCH2) * (HCO_HCCO)) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO))) - ( (TXCH2_ct) &
                 - ( ( (SXCH2_ct) - ( (CH3O_ct) - (CH2OH_ct) &
                 * (CH3O_CH2OH)) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_ct) &
                 * (SXCH2_HCCO)) * (TXCH2_SXCH2) / ( ( (1.0_pr) &
                 - (CH3O_SXCH2) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_SXCH2) &
                 * (SXCH2_HCCO))) * ( - (SXCH2_TXCH2) * ( - (HCCO_SXCH2) &
                 * (HCO_HCCO)) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO))) &
                 / ( (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO)))) - ( (CH2CHO_ct) &
                 - (C2H5_ct) * (CH2CHO_C2H5)) * (HCO_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5))) &
                 - ( (CH_ct) - ( (TXCH2_ct) - ( ( (SXCH2_ct) &
                 - ( (CH3O_ct) - (CH2OH_ct) * (CH3O_CH2OH)) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_ct) * (SXCH2_HCCO)) &
                 * (TXCH2_SXCH2) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO))) &
                 * (CH_TXCH2) / ( (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO)))) * ( (HCO_CH) &
                 - (TXCH2_CH) * ( - (SXCH2_TXCH2) * ( - (HCCO_SXCH2) &
                 * (HCO_HCCO)) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO))) &
                 / ( (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO)))) / ( (1.0_pr) &
                 - (TXCH2_CH) * (CH_TXCH2) / ( (1.0_pr) - (SXCH2_TXCH2) &
                 * (TXCH2_SXCH2) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO))))) &
                 - ( ( (C2H3_ct) - (C3H5_ct) * (C2H3_C3H5)) &
                 - ( (CH2CHO_ct) - (C2H5_ct) * (CH2CHO_C2H5)) &
                 * (C2H3_CH2CHO) / ( (1.0_pr) - (C2H5_CH2CHO) &
                 * (CH2CHO_C2H5))) * ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO) / ( (1.0_pr) - (C2H5_CH2CHO) &
                 * (CH2CHO_C2H5))) / ( ( (1.0_pr) - (C3H5_C2H3) &
                 * (C2H3_C3H5)) - (CH2CHO_C2H3) * (C2H3_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5))) &
                 ) / (  ( ( ( (1.0_pr) - (HCCO_HCO) * (HCO_HCCO)) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO) / ( (1.0_pr) &
                 - (C2H5_CH2CHO) * (CH2CHO_C2H5))) - (CH_HCO) &
                 * ( (HCO_CH) - (TXCH2_CH) * ( - (SXCH2_TXCH2) &
                 * ( - (HCCO_SXCH2) * (HCO_HCCO)) / ( ( (1.0_pr) &
                 - (CH3O_SXCH2) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_SXCH2) &
                 * (SXCH2_HCCO))) / ( (1.0_pr) - (SXCH2_TXCH2) &
                 * (TXCH2_SXCH2) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO)))) &
                 / ( (1.0_pr) - (TXCH2_CH) * (CH_TXCH2) / ( (1.0_pr) &
                 - (SXCH2_TXCH2) * (TXCH2_SXCH2) / ( ( (1.0_pr) &
                 - (CH3O_SXCH2) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_SXCH2) &
                 * (SXCH2_HCCO))))) - ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO) / ( (1.0_pr) - (C2H5_CH2CHO) &
                 * (CH2CHO_C2H5))) * ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO) / ( (1.0_pr) - (C2H5_CH2CHO) &
                 * (CH2CHO_C2H5))) / ( ( (1.0_pr) - (C3H5_C2H3) &
                 * (C2H3_C3H5)) - (CH2CHO_C2H3) * (C2H3_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5))) &
                 )  

    cqss(sqssC2H3) = ( ( (C2H3_ct) - (C3H5_ct) * (C2H3_C3H5)) - ( (CH2CHO_ct) &
                 - (C2H5_ct) * (CH2CHO_C2H5)) * (C2H3_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5)) &
                 - ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5))) &
                 * cqss(sqssHCO)) / ( ( (1.0_pr) - (C3H5_C2H3) &
                 * (C2H3_C3H5)) - (CH2CHO_C2H3) * (C2H3_CH2CHO) &
                 / ( (1.0_pr) - (C2H5_CH2CHO) * (CH2CHO_C2H5)))  

    cqss(sqssCH) = ( (CH_ct) - ( (TXCH2_ct) - ( ( (SXCH2_ct) &
                 - ( (CH3O_ct) - (CH2OH_ct) * (CH3O_CH2OH)) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_ct) * (SXCH2_HCCO)) &
                 * (TXCH2_SXCH2) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO))) &
                 * (CH_TXCH2) / ( (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO))) - (CH_HCO) &
                 * cqss(sqssHCO)) / ( (1.0_pr) - (TXCH2_CH) &
                 * (CH_TXCH2) / ( (1.0_pr) - (SXCH2_TXCH2) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO))))  

    cqss(sqssCH2CHO) = ( (CH2CHO_ct) - (C2H5_ct) * (CH2CHO_C2H5) &
                 - (CH2CHO_HCO) * cqss(sqssHCO) - (CH2CHO_C2H3) &
                 * cqss(sqssC2H3)) / ( (1.0_pr) - (C2H5_CH2CHO) &
                 * (CH2CHO_C2H5))  

    cqss(sqssTXCH2) = ( (TXCH2_ct) - ( ( (SXCH2_ct) - ( (CH3O_ct) &
                 - (CH2OH_ct) * (CH3O_CH2OH)) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_ct) * (SXCH2_HCCO)) * (TXCH2_SXCH2) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO)) - (TXCH2_CH) &
                 * cqss(sqssCH)) / ( (1.0_pr) - (SXCH2_TXCH2) &
                 * (TXCH2_SXCH2) / ( ( (1.0_pr) - (CH3O_SXCH2) &
                 * (SXCH2_CH3O) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))) - (HCCO_SXCH2) * (SXCH2_HCCO)))  

    cqss(sqssSXCH2) = ( ( (SXCH2_ct) - ( (CH3O_ct) - (CH2OH_ct) &
                 * (CH3O_CH2OH)) * (SXCH2_CH3O) / ( (1.0_pr) &
                 - (CH2OH_CH3O) * (CH3O_CH2OH))) - (HCCO_ct) &
                 * (SXCH2_HCCO) - ( - (HCCO_HCO) * (SXCH2_HCCO)) &
                 * cqss(sqssHCO) - (SXCH2_TXCH2) * cqss(sqssTXCH2)) &
                 / ( ( (1.0_pr) - (CH3O_SXCH2) * (SXCH2_CH3O) &
                 / ( (1.0_pr) - (CH2OH_CH3O) * (CH3O_CH2OH))) &
                 - (HCCO_SXCH2) * (SXCH2_HCCO))  

    cqss(sqssC2H5) = C2H5_ct - (C2H5_CH2CHO) * cqss(sqssCH2CHO)  

    cqss(sqssC3H5) = C3H5_ct - (C3H5_C2H3) * cqss(sqssC2H3)  

    cqss(sqssHCCO) = HCCO_ct - (HCCO_HCO) * cqss(sqssHCO) - (HCCO_SXCH2) &
                 * cqss(sqssSXCH2)  

    cqss(sqssCH3O) = ( (CH3O_ct) - (CH2OH_ct) * (CH3O_CH2OH) - (CH3O_SXCH2) &
                 * cqss(sqssSXCH2)) / ( (1.0_pr) - (CH2OH_CH3O) &
                 * (CH3O_CH2OH))  

    cqss(sqssCH2OH) = CH2OH_ct - (CH2OH_CH3O) * cqss(sqssCH3O)  

    IXC3H7_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r152f)* c(sC3H8) &
             + k(r154f)* c(sO2) &
             + k(r142b)* c(sHO2) &
             + k(r144b)* c(sH2) &
             + k(r146b)* c(sOH) &
             + k(r149b)* c(sH2O) &
             + k(r150b)* c(sH2O2) &
             + k(r153b)  )

    IXC3H7_num = ( 0.0_pr &
             + k(r142f)* c(sC3H8) * c(sO2) &
             + k(r144f)* c(sC3H8) * c(sH) &
             + k(r146f)* c(sC3H8) * c(sO) &
             + k(r149f)* c(sC3H8) * c(sOH) &
             + k(r150f)* c(sC3H8) * c(sHO2) &
             + k(r153f)* c(sC3H6) * c(sH) &
             + k(r154b)* c(sC3H6) * c(sHO2)   )

    IXC3H7_ct = IXC3H7_num / IXC3H7_denom

    IXC3H7_NXC3H7 = - ( 0.0_pr &
          + k(r152b) * c(sC3H8)  ) / IXC3H7_denom

    NXC3H7_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r155f)&
             + k(r157f)* c(sO2) &
             + k(r158f)* c(sO2) &
             + k(r143b)* c(sHO2) &
             + k(r145b)* c(sH2) &
             + k(r147b)* c(sOH) &
             + k(r148b)* c(sH2O) &
             + k(r151b)* c(sH2O2) &
             + k(r152b)* c(sC3H8) &
             + k(r156b)  )

    NXC3H7_num = ( 0.0_pr &
             + k(r143f)* c(sC3H8) * c(sO2) &
             + k(r145f)* c(sC3H8) * c(sH) &
             + k(r147f)* c(sC3H8) * c(sO) &
             + k(r148f)* c(sC3H8) * c(sOH) &
             + k(r151f)* c(sC3H8) * c(sHO2) &
             + k(r156f)* c(sC3H6) * c(sH) &
             + k(r155b)* c(sC2H4) * c(sCH3) &
             + k(r157b)* c(sC3H6) * c(sHO2)   )

    NXC3H7_ct = NXC3H7_num / NXC3H7_denom

    NXC3H7_IXC3H7 = - ( 0.0_pr &
          + k(r152f) * c(sC3H8)  ) / NXC3H7_denom

    NXC3H7_C3H6OOH = - ( 0.0_pr &
          + k(r158b)  ) / NXC3H7_denom

    C3H6OOH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r159f)&
             + k(r160f)* c(sO2) &
             + k(r158b)  )

    C3H6OOH_num = ( 0.0_pr &
             + k(r159b)* c(sC3H6) * c(sHO2) &
             + k(r160b)* c(sOC3H5OOH) * c(sOH)   )

    C3H6OOH_ct = C3H6OOH_num / C3H6OOH_denom

    C3H6OOH_NXC3H7 = - ( 0.0_pr &
          + k(r158f) * c(sO2)  ) / C3H6OOH_denom

    cqss(sqssIXC3H7) = (  (IXC3H7_ct) - ( (NXC3H7_ct) - (C3H6OOH_ct) &
                 * (NXC3H7_C3H6OOH)) * (IXC3H7_NXC3H7) / ( (1.0_pr) &
                 - (C3H6OOH_NXC3H7) * (NXC3H7_C3H6OOH)) ) / (  (1.0_pr) &
                 - (NXC3H7_IXC3H7) * (IXC3H7_NXC3H7) / ( (1.0_pr) &
                 - (C3H6OOH_NXC3H7) * (NXC3H7_C3H6OOH)) )  

    cqss(sqssNXC3H7) = ( (NXC3H7_ct) - (C3H6OOH_ct) * (NXC3H7_C3H6OOH) &
                 - (NXC3H7_IXC3H7) * cqss(sqssIXC3H7)) / ( (1.0_pr) &
                 - (C3H6OOH_NXC3H7) * (NXC3H7_C3H6OOH))  

    cqss(sqssC3H6OOH) = C3H6OOH_ct - (C3H6OOH_NXC3H7) * cqss(sqssNXC3H7)  

    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, P, T, c)
    implicit none

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

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

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

    ! Conversion
    c = y * rho / W_sp

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

        return
    end subroutine y2c

end module mod_customkinetics

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

  use mod_customkinetics
  implicit none

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

  real(pr) :: P, T, rho

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

  ! Evaluate QSS concentrations and reaction rates
  call get_thirdbodies(M,c)

  call get_rate_coefficients(k, M, T, P)
  call get_qss(cqss, c, k, M)

  call get_reaction_rates(w, k, M, c, cqss)

  ! Evaluate production rates
  call get_production_rates(wdot,w)
  
  return
end subroutine customkinetics