!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE reduced/reducedS31R224QSS24_1.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, J. Wirtz, P. Pepiot
!!    @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 = 31
  integer, parameter :: nqss = 24
  integer, parameter :: nreac = 220
  integer, parameter :: nreac_reverse = 166
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 3
  integer, parameter :: nFO = 11
  integer, parameter :: nTB_reverse = 3
  integer, parameter :: nFO_reverse = 11
  
  ! Index of species
  integer, parameter :: sN2 = 1
  integer, parameter :: sH2 = 2
  integer, parameter :: sH = 3
  integer, parameter :: sO2 = 4
  integer, parameter :: sO = 5
  integer, parameter :: sH2O = 6
  integer, parameter :: sOH = 7
  integer, parameter :: sH2O2 = 8
  integer, parameter :: sHO2 = 9
  integer, parameter :: sCO = 10
  integer, parameter :: sCO2 = 11
  integer, parameter :: sCH4 = 12
  integer, parameter :: sCH3 = 13
  integer, parameter :: sCH2O = 14
  integer, parameter :: sC2H6 = 15
  integer, parameter :: sC2H5 = 16
  integer, parameter :: sC2H4 = 17
  integer, parameter :: sC2H3 = 18
  integer, parameter :: sC2H2 = 19
  integer, parameter :: sC2H = 20
  integer, parameter :: sC3H6 = 21
  integer, parameter :: sC3H5XT = 22
  integer, parameter :: sC3H4XP = 23
  integer, parameter :: sC3H4XA = 24
  integer, parameter :: sC2H3CHO = 25
  integer, parameter :: sIC4H8 = 26
  integer, parameter :: sIC4H7 = 27
  integer, parameter :: sIC8H18 = 28
  integer, parameter :: sIC8H16 = 29
  integer, parameter :: sIC16H34 = 30
  integer, parameter :: sIC12H26 = 31

  integer, parameter :: sqssCH2 = 1
  integer, parameter :: sqssCH2GSG = 2
  integer, parameter :: sqssC = 3
  integer, parameter :: sqssCH = 4
  integer, parameter :: sqssCH3O = 5
  integer, parameter :: sqssCH2OH = 6
  integer, parameter :: sqssHCO = 7
  integer, parameter :: sqssCH3CO = 8
  integer, parameter :: sqssCH2CHO = 9
  integer, parameter :: sqssCH2CO = 10
  integer, parameter :: sqssHCCO = 11
  integer, parameter :: sqssC3H5XA = 12
  integer, parameter :: sqssC3H5O = 13
  integer, parameter :: sqssC3H3 = 14
  integer, parameter :: sqssCH3COCH2 = 15
  integer, parameter :: sqssIC4H9 = 16
  integer, parameter :: sqssTC4H9 = 17
  integer, parameter :: sqssIC4H7O = 18
  integer, parameter :: sqssNEOC5H11 = 19
  integer, parameter :: sqssIC8H17 = 20
  integer, parameter :: sqssIC16H33 = 21
  integer, parameter :: sqssIC12H25 = 22
  integer, parameter :: sqssC4H9O2GL1G = 23
  integer, parameter :: sqssC3H7GL1G = 24

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r39f = 39
  integer, parameter :: r40f = 40
  integer, parameter :: r41f = 41
  integer, parameter :: r42f = 42
  integer, parameter :: r43f = 43
  integer, parameter :: r44f = 44
  integer, parameter :: r45f = 45
  integer, parameter :: r46f = 46
  integer, parameter :: r47f = 47
  integer, parameter :: r48f = 48
  integer, parameter :: r49f = 49
  integer, parameter :: r50f = 50
  integer, parameter :: r51f = 51
  integer, parameter :: r52f = 52
  integer, parameter :: r53f = 53
  integer, parameter :: r54f = 54
  integer, parameter :: r55f = 55
  integer, parameter :: r56f = 56
  integer, parameter :: r57f = 57
  integer, parameter :: r58f = 58
  integer, parameter :: r59f = 59
  integer, parameter :: r60f = 60
  integer, parameter :: r61f = 61
  integer, parameter :: r62f = 62
  integer, parameter :: r63f = 63
  integer, parameter :: r64f = 64
  integer, parameter :: r65f = 65
  integer, parameter :: r66f = 66
  integer, parameter :: r67f = 67
  integer, parameter :: r68f = 68
  integer, parameter :: r69f = 69
  integer, parameter :: r70f = 70
  integer, parameter :: r71f = 71
  integer, parameter :: r72f = 72
  integer, parameter :: r73f = 73
  integer, parameter :: r74f = 74
  integer, parameter :: r75f = 75
  integer, parameter :: r76f = 76
  integer, parameter :: r77f = 77
  integer, parameter :: r78f = 78
  integer, parameter :: r79f = 79
  integer, parameter :: r80f = 80
  integer, parameter :: r81f = 81
  integer, parameter :: r82f = 82
  integer, parameter :: r83f = 83
  integer, parameter :: r84f = 84
  integer, parameter :: r85f = 85
  integer, parameter :: r86f = 86
  integer, parameter :: r87f = 87
  integer, parameter :: r88f = 88
  integer, parameter :: r89f = 89
  integer, parameter :: r90f = 90
  integer, parameter :: r91f = 91
  integer, parameter :: r92f = 92
  integer, parameter :: r93f = 93
  integer, parameter :: r94f = 94
  integer, parameter :: r95f = 95
  integer, parameter :: r96f = 96
  integer, parameter :: r97f = 97
  integer, parameter :: r98f = 98
  integer, parameter :: r99f = 99
  integer, parameter :: r100f = 100
  integer, parameter :: r101f = 101
  integer, parameter :: r102f = 102
  integer, parameter :: r103f = 103
  integer, parameter :: r104f = 104
  integer, parameter :: r105f = 105
  integer, parameter :: r106f = 106
  integer, parameter :: r107f = 107
  integer, parameter :: r108f = 108
  integer, parameter :: r109f = 109
  integer, parameter :: r110f = 110
  integer, parameter :: r111f = 111
  integer, parameter :: r112f = 112
  integer, parameter :: r113f = 113
  integer, parameter :: r114f = 114
  integer, parameter :: r115f = 115
  integer, parameter :: r116f = 116
  integer, parameter :: r117f = 117
  integer, parameter :: r118f = 118
  integer, parameter :: r119f = 119
  integer, parameter :: r120f = 120
  integer, parameter :: r121f = 121
  integer, parameter :: r122f = 122
  integer, parameter :: r123f = 123
  integer, parameter :: r124f = 124
  integer, parameter :: r125f = 125
  integer, parameter :: r126f = 126
  integer, parameter :: r127f = 127
  integer, parameter :: r128f = 128
  integer, parameter :: r129f = 129
  integer, parameter :: r130f = 130
  integer, parameter :: r131f = 131
  integer, parameter :: r132f = 132
  integer, parameter :: r133f = 133
  integer, parameter :: r134f = 134
  integer, parameter :: r135f = 135
  integer, parameter :: r136f = 136
  integer, parameter :: r137f = 137
  integer, parameter :: r138f = 138
  integer, parameter :: r139f = 139
  integer, parameter :: r140f = 140
  integer, parameter :: r141f = 141
  integer, parameter :: r142f = 142
  integer, parameter :: r143f = 143
  integer, parameter :: r144f = 144
  integer, parameter :: r145f = 145
  integer, parameter :: r146f = 146
  integer, parameter :: r147f = 147
  integer, parameter :: r148f = 148
  integer, parameter :: r149f = 149
  integer, parameter :: r150f = 150
  integer, parameter :: r151f = 151
  integer, parameter :: r152f = 152
  integer, parameter :: r153f = 153
  integer, parameter :: r154f = 154
  integer, parameter :: r155f = 155
  integer, parameter :: r156f = 156
  integer, parameter :: r157f = 157
  integer, parameter :: r158f = 158
  integer, parameter :: r159f = 159
  integer, parameter :: r160f = 160
  integer, parameter :: r161f = 161
  integer, parameter :: r162f = 162
  integer, parameter :: r163f = 163
  integer, parameter :: r164f = 164
  integer, parameter :: r165f = 165
  integer, parameter :: r166f = 166
  integer, parameter :: r167f = 167
  integer, parameter :: r168f = 168
  integer, parameter :: r169f = 169
  integer, parameter :: r170f = 170
  integer, parameter :: r171f = 171
  integer, parameter :: r172f = 172
  integer, parameter :: r173f = 173
  integer, parameter :: r174f = 174
  integer, parameter :: r175f = 175
  integer, parameter :: r176f = 176
  integer, parameter :: r177f = 177
  integer, parameter :: r178f = 178
  integer, parameter :: r179f = 179
  integer, parameter :: r180f = 180
  integer, parameter :: r181f = 181
  integer, parameter :: r182f = 182
  integer, parameter :: r183f = 183
  integer, parameter :: r184f = 184
  integer, parameter :: r185f = 185
  integer, parameter :: r186f = 186
  integer, parameter :: r187f = 187
  integer, parameter :: r188f = 188
  integer, parameter :: r189f = 189
  integer, parameter :: r190f = 190
  integer, parameter :: r191f = 191
  integer, parameter :: r192f = 192
  integer, parameter :: r193f = 193
  integer, parameter :: r194f = 194
  integer, parameter :: r195f = 195
  integer, parameter :: r196f = 196
  integer, parameter :: r197f = 197
  integer, parameter :: r198f = 198
  integer, parameter :: r199f = 199
  integer, parameter :: r200f = 200
  integer, parameter :: r201f = 201
  integer, parameter :: r202f = 202
  integer, parameter :: r203f = 203
  integer, parameter :: r204f = 204
  integer, parameter :: r205f = 205
  integer, parameter :: r206f = 206
  integer, parameter :: r207f = 207
  integer, parameter :: r208f = 208
  integer, parameter :: r209f = 209
  integer, parameter :: r210f = 210
  integer, parameter :: r211f = 211
  integer, parameter :: r212f = 212
  integer, parameter :: r213f = 213
  integer, parameter :: r214f = 214
  integer, parameter :: r215f = 215
  integer, parameter :: r216f = 216
  integer, parameter :: r217f = 217
  integer, parameter :: r218f = 218
  integer, parameter :: r219f = 219
  integer, parameter :: r220f = 220
  integer, parameter :: r1b = 221
  integer, parameter :: r2b = 222
  integer, parameter :: r3b = 223
  integer, parameter :: r4b = 224
  integer, parameter :: r5b = 225
  integer, parameter :: r6b = 226
  integer, parameter :: r7b = 227
  integer, parameter :: r8b = 228
  integer, parameter :: r9b = 229
  integer, parameter :: r10b = 230
  integer, parameter :: r11b = 231
  integer, parameter :: r12b = 232
  integer, parameter :: r13b = 233
  integer, parameter :: r14b = 234
  integer, parameter :: r15b = 235
  integer, parameter :: r16b = 236
  integer, parameter :: r17b = 237
  integer, parameter :: r18b = 238
  integer, parameter :: r19b = 239
  integer, parameter :: r20b = 240
  integer, parameter :: r21b = 241
  integer, parameter :: r22b = 242
  integer, parameter :: r23b = 243
  integer, parameter :: r25b = 244
  integer, parameter :: r26b = 245
  integer, parameter :: r27b = 246
  integer, parameter :: r29b = 247
  integer, parameter :: r30b = 248
  integer, parameter :: r31b = 249
  integer, parameter :: r32b = 250
  integer, parameter :: r33b = 251
  integer, parameter :: r34b = 252
  integer, parameter :: r35b = 253
  integer, parameter :: r36b = 254
  integer, parameter :: r37b = 255
  integer, parameter :: r38b = 256
  integer, parameter :: r39b = 257
  integer, parameter :: r40b = 258
  integer, parameter :: r41b = 259
  integer, parameter :: r42b = 260
  integer, parameter :: r43b = 261
  integer, parameter :: r44b = 262
  integer, parameter :: r45b = 263
  integer, parameter :: r46b = 264
  integer, parameter :: r47b = 265
  integer, parameter :: r48b = 266
  integer, parameter :: r49b = 267
  integer, parameter :: r50b = 268
  integer, parameter :: r52b = 269
  integer, parameter :: r53b = 270
  integer, parameter :: r54b = 271
  integer, parameter :: r55b = 272
  integer, parameter :: r56b = 273
  integer, parameter :: r57b = 274
  integer, parameter :: r58b = 275
  integer, parameter :: r59b = 276
  integer, parameter :: r60b = 277
  integer, parameter :: r61b = 278
  integer, parameter :: r62b = 279
  integer, parameter :: r63b = 280
  integer, parameter :: r64b = 281
  integer, parameter :: r65b = 282
  integer, parameter :: r66b = 283
  integer, parameter :: r67b = 284
  integer, parameter :: r68b = 285
  integer, parameter :: r69b = 286
  integer, parameter :: r70b = 287
  integer, parameter :: r71b = 288
  integer, parameter :: r72b = 289
  integer, parameter :: r73b = 290
  integer, parameter :: r76b = 291
  integer, parameter :: r77b = 292
  integer, parameter :: r78b = 293
  integer, parameter :: r79b = 294
  integer, parameter :: r80b = 295
  integer, parameter :: r81b = 296
  integer, parameter :: r82b = 297
  integer, parameter :: r83b = 298
  integer, parameter :: r84b = 299
  integer, parameter :: r85b = 300
  integer, parameter :: r86b = 301
  integer, parameter :: r87b = 302
  integer, parameter :: r88b = 303
  integer, parameter :: r90b = 304
  integer, parameter :: r92b = 305
  integer, parameter :: r93b = 306
  integer, parameter :: r94b = 307
  integer, parameter :: r95b = 308
  integer, parameter :: r96b = 309
  integer, parameter :: r97b = 310
  integer, parameter :: r98b = 311
  integer, parameter :: r99b = 312
  integer, parameter :: r100b = 313
  integer, parameter :: r101b = 314
  integer, parameter :: r102b = 315
  integer, parameter :: r103b = 316
  integer, parameter :: r104b = 317
  integer, parameter :: r105b = 318
  integer, parameter :: r106b = 319
  integer, parameter :: r107b = 320
  integer, parameter :: r108b = 321
  integer, parameter :: r109b = 322
  integer, parameter :: r110b = 323
  integer, parameter :: r111b = 324
  integer, parameter :: r112b = 325
  integer, parameter :: r113b = 326
  integer, parameter :: r114b = 327
  integer, parameter :: r115b = 328
  integer, parameter :: r116b = 329
  integer, parameter :: r117b = 330
  integer, parameter :: r118b = 331
  integer, parameter :: r119b = 332
  integer, parameter :: r120b = 333
  integer, parameter :: r121b = 334
  integer, parameter :: r122b = 335
  integer, parameter :: r123b = 336
  integer, parameter :: r124b = 337
  integer, parameter :: r125b = 338
  integer, parameter :: r126b = 339
  integer, parameter :: r127b = 340
  integer, parameter :: r128b = 341
  integer, parameter :: r129b = 342
  integer, parameter :: r130b = 343
  integer, parameter :: r131b = 344
  integer, parameter :: r132b = 345
  integer, parameter :: r133b = 346
  integer, parameter :: r134b = 347
  integer, parameter :: r135b = 348
  integer, parameter :: r136b = 349
  integer, parameter :: r137b = 350
  integer, parameter :: r138b = 351
  integer, parameter :: r139b = 352
  integer, parameter :: r140b = 353
  integer, parameter :: r141b = 354
  integer, parameter :: r142b = 355
  integer, parameter :: r143b = 356
  integer, parameter :: r144b = 357
  integer, parameter :: r145b = 358
  integer, parameter :: r146b = 359
  integer, parameter :: r147b = 360
  integer, parameter :: r148b = 361
  integer, parameter :: r149b = 362
  integer, parameter :: r150b = 363
  integer, parameter :: r151b = 364
  integer, parameter :: r152b = 365
  integer, parameter :: r153b = 366
  integer, parameter :: r154b = 367
  integer, parameter :: r155b = 368
  integer, parameter :: r156b = 369
  integer, parameter :: r157b = 370
  integer, parameter :: r158b = 371
  integer, parameter :: r159b = 372
  integer, parameter :: r160b = 373
  integer, parameter :: r161b = 374
  integer, parameter :: r165b = 375
  integer, parameter :: r166b = 376
  integer, parameter :: r167b = 377
  integer, parameter :: r168b = 378
  integer, parameter :: r169b = 379
  integer, parameter :: r176b = 380
  integer, parameter :: r178b = 381
  integer, parameter :: r179b = 382
  integer, parameter :: r180b = 383
  integer, parameter :: r181b = 384
  integer, parameter :: r196b = 385
  integer, parameter :: r198b = 386
  
  ! Index of third body species
  integer, parameter :: mM4 = 1
  integer, parameter :: mM47 = 2
  integer, parameter :: mM165 = 3
  
  integer, parameter :: mM6 = 4
  integer, parameter :: mM16 = 5
  integer, parameter :: mM19 = 6
  integer, parameter :: mM52 = 7
  integer, parameter :: mM53 = 8
  integer, parameter :: mM54 = 9
  integer, parameter :: mM58 = 10
  integer, parameter :: mM65 = 11
  integer, parameter :: mM86 = 12
  integer, parameter :: mM92 = 13
  integer, parameter :: mM178 = 14
  
  ! Index of third body reactions
  integer, parameter :: TBr4f = 1
  integer, parameter :: TBr47f = 2
  integer, parameter :: TBr165f = 3
  integer, parameter :: TBr4b = 4
  integer, parameter :: TBr47b = 5
  integer, parameter :: TBr165b = 6
  
  ! Index of fall off reactions
  integer, parameter :: FOr6f = 1
  integer, parameter :: FOr16f = 2
  integer, parameter :: FOr19f = 3
  integer, parameter :: FOr52f = 4
  integer, parameter :: FOr53f = 5
  integer, parameter :: FOr54f = 6
  integer, parameter :: FOr58f = 7
  integer, parameter :: FOr65f = 8
  integer, parameter :: FOr86f = 9
  integer, parameter :: FOr92f = 10
  integer, parameter :: FOr178f = 11
  integer, parameter :: FOr6b = 12
  integer, parameter :: FOr16b = 13
  integer, parameter :: FOr19b = 14
  integer, parameter :: FOr52b = 15
  integer, parameter :: FOr53b = 16
  integer, parameter :: FOr54b = 17
  integer, parameter :: FOr58b = 18
  integer, parameter :: FOr65b = 19
  integer, parameter :: FOr86b = 20
  integer, parameter :: FOr92b = 21
  integer, parameter :: FOr178b = 22
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.02801348_pr, & ! N2
       0.00201588_pr, & ! H2
       0.00100794_pr, & ! H
       0.0319988_pr, & ! O2
       0.0159994_pr, & ! O
       0.01801528_pr, & ! H2O
       0.01700734_pr, & ! OH
       0.03401468_pr, & ! H2O2
       0.03300674_pr, & ! HO2
       0.0280104_pr, & ! CO
       0.0440098_pr, & ! CO2
       0.01604276_pr, & ! CH4
       0.01503482_pr, & ! CH3
       0.03002628_pr, & ! CH2O
       0.03006964_pr, & ! C2H6
       0.0290617_pr, & ! C2H5
       0.02805376_pr, & ! C2H4
       0.02704582_pr, & ! C2H3
       0.02603788_pr, & ! C2H2
       0.02502994_pr, & ! C2H
       0.04208064_pr, & ! C3H6
       0.0410727_pr, & ! C3H5XT
       0.04006476_pr, & ! C3H4XP
       0.04006476_pr, & ! C3H4XA
       0.05606416_pr, & ! C2H3CHO
       0.05610752_pr, & ! IC4H8
       0.05509958_pr, & ! IC4H7
       0.11423092_pr, & ! IC8H18
       0.11221504_pr, & ! IC8H16
       0.22644596_pr, & ! IC16H34
       0.17033844_pr & ! IC12H26
      !0.01402688_pr, & ! CH2
      !0.01402688_pr, & ! CH2GSG
      !0.012011_pr, & ! C
      !0.01301894_pr, & ! CH
      !0.03103422_pr, & ! CH3O
      !0.03103422_pr, & ! CH2OH
      !0.02901834_pr, & ! HCO
      !0.04304522_pr, & ! CH3CO
      !0.04304522_pr, & ! CH2CHO
      !0.04203728_pr, & ! CH2CO
      !0.04102934_pr, & ! HCCO
      !0.0410727_pr, & ! C3H5XA
      !0.0570721_pr, & ! C3H5O
      !0.03905682_pr, & ! C3H3
      !0.0570721_pr, & ! CH3COCH2
      !0.05711546_pr, & ! IC4H9
      !0.05711546_pr, & ! TC4H9
      !0.07109898_pr, & ! IC4H7O
      !0.07114234_pr, & ! NEOC5H11
      !0.11322298_pr, & ! IC8H17
      !0.22543802_pr, & ! IC16H33
      !0.1693305_pr, & ! IC12H25
      !0.08911426_pr, & ! C4H9O2GL1G
      !0.04308858_pr & ! C3H7GL1G
       /)
  
contains

  ! ----------------------------------------------- !
  ! Subroutine for pressure dependent coefficients  !
  ! ----------------------------------------------- !
  real(pr) function getlindratecoeff(Tloc,k0,kinf,fc,concin,Ploc)
    implicit none

    real(pr) ::  Tloc,k0,kinf,fc,Ploc
    real(pr) :: ntmp,ccoeff,dcoeff,lgknull
    real(pr) :: f
    real(pr) :: conc, concin

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

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

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

    M(mM165) = sum(c)

    M(mM6) = (1.8_pr)*c(sCO) &
         + (0.6_pr)*c(sCO2) &
         + (6.65_pr)*c(sH2O) &
         + (0.5_pr)*c(sN2) &
         + (0.2_pr)*c(sO2) &
         + sum(c)

    M(mM16) = (0.9_pr)*c(sCO) &
         + (2.8_pr)*c(sCO2) &
         + (9_pr)*c(sH2O) &
         + sum(c)

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

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

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

    M(mM54) = (1_pr)*c(sCO) &
         + (2_pr)*c(sCO2) &
         + (4_pr)*c(sH2O) &
         + sum(c)

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

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

    M(mM86) = sum(c)

    M(mM92) = sum(c)

    M(mM178) = (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (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) = (0.0508_pr)*exp((-26325.728_pr)*R_T_inv + T_log *(2.67_pr))
    k(r2f) = (43800000_pr)*exp((-29246.16_pr)*R_T_inv)
    k(r3f) = (114000000_pr)*exp((-63956.624_pr)*R_T_inv)
    k(r4f) = (35000000000_pr)*exp(T_log *(-2_pr))
    k(r5f) = (67_pr)*exp((-62704.7712_pr)*R_T_inv + T_log *(1.704_pr))
    k_0(FOr6f) = (2.49e+18_pr)*exp((-203965.816_pr)*R_T_inv + T_log *(-2.3_pr))
    k_inf(FOr6f) = (2000000000000_pr)*exp((-203965.816_pr)*R_T_inv + T_log *(0.9_pr))
    FC(FOr6f) = ((1.0_pr - 0.43_pr)*exp(-Tloc/(1e-30_pr))) + (0.43_pr*exp(-Tloc/(1e+30_pr)))
    k(r6f) = getlindratecoeff(Tloc,k_0(FOr6f),k_inf(FOr6f),FC(FOr6f), M(mM6),Ploc)
    k(r7f) = (1740000_pr)*exp((-1330.512_pr)*R_T_inv)
    k(r8f) = (75900000_pr)*exp((-30413.496_pr)*R_T_inv)
    k(r9f) = (70790000_pr)*exp((-1234.28_pr)*R_T_inv)
    k(r10f) = (11402_pr)*exp((-2317.01552_pr)*R_T_inv + T_log *(1.083_pr))
    k(r11f) = (32500000_pr)
    k(r12f) = (7000000_pr)*exp((4572.94464_pr)*R_T_inv)
    k(r13f) = (450000000_pr)*exp((-45729.4464_pr)*R_T_inv)
    k(r14f) = (100000000_pr)*exp((-46195.04192_pr)*R_T_inv)
    k(r15f) = (190000_pr)*exp((5894.92128_pr)*R_T_inv)
    k_0(FOr16f) = (17400000_pr)*exp(T_log *(-1.23_pr))
    k_inf(FOr16f) = (4650000_pr)*exp(T_log *(0.44_pr))
    FC(FOr16f) = ((1.0_pr - 0.67_pr)*exp(-Tloc/(1e-30_pr))) + (0.67_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r16f) = getlindratecoeff(Tloc,k_0(FOr16f),k_inf(FOr16f),FC(FOr16f), M(mM16),Ploc)
    k(r17f) = (0.07015_pr)*exp((1488.2488_pr)*R_T_inv + T_log *(2.053_pr))
    k(r18f) = (5757000_pr)*exp((-1388.2512_pr)*R_T_inv + T_log *(-0.664_pr))
    k_0(FOr19f) = (2.48e+21_pr)*exp((-10208.96_pr)*R_T_inv + T_log *(-4.76_pr))
    k_inf(FOr19f) = (12700000000_pr)*exp((-1602.472_pr)*R_T_inv + T_log *(-0.63_pr))
    FC(FOr19f) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r19f) = getlindratecoeff(Tloc,k_0(FOr19f),k_inf(FOr19f),FC(FOr19f), M(mM19),Ploc)
    k(r20f) = (0.614_pr)*exp((-40112.008_pr)*R_T_inv + T_log *(2.5_pr))
    k(r21f) = (0.0583_pr)*exp((-9162.96_pr)*R_T_inv + T_log *(2.6_pr))
    k(r22f) = (0.116_pr)*exp((12644.048_pr)*R_T_inv + T_log *(2.23_pr))
    k(r23f) = (15000000_pr)*exp((-2510.4_pr)*R_T_inv)
    k(r24f) = (28000000_pr)
    k(r25f) = (12000000_pr)
    k(r26f) = (70000000_pr)
    k(r27f) = (10600000_pr)*exp((-6276_pr)*R_T_inv)
    k(r28f) = (2640000_pr)*exp((-6276_pr)*R_T_inv)
    k(r29f) = (30000000_pr)
    k(r30f) = (33000000_pr)
    k(r31f) = (110000000_pr)
    k(r32f) = (50000000_pr)
    k(r33f) = (7546000_pr)*exp((-118490.88_pr)*R_T_inv)
    k(r34f) = (2.641e-06_pr)*exp((-33911.32_pr)*R_T_inv + T_log *(3.283_pr))
    k(r35f) = (55400000_pr)*exp((569.024_pr)*R_T_inv + T_log *(0.05_pr))
    k(r36f) = (528200000000_pr)*exp((-7414.048_pr)*R_T_inv + T_log *(-1.518_pr))
    k(r37f) = (46860_pr)*exp((-14920.144_pr)*R_T_inv + T_log *(0.833_pr))
    k(r38f) = (0.04293_pr)*exp((-16726.7952_pr)*R_T_inv + T_log *(2.568_pr))
    k(r39f) = (1000000_pr)*exp((2876.5_pr)*R_T_inv + T_log *(0.269_pr))
    k(r40f) = (1510000000_pr)*exp(T_log *(-1_pr))
    k(r41f) = (241000000_pr)*exp((-20991.128_pr)*R_T_inv)
    k(r42f) = (4.38e-25_pr)*exp((23016.184_pr)*R_T_inv + T_log *(9.5_pr))
    k(r43f) = (8070000000_pr)*exp((-223509.28_pr)*R_T_inv)
    k(r44f) = (57.4_pr)*exp((-11464.16_pr)*R_T_inv + T_log *(1.9_pr))
    k(r45f) = (78.2_pr)*exp((4414.12_pr)*R_T_inv + T_log *(1.63_pr))
    k(r46f) = (0.0188_pr)*exp((-48199.68_pr)*R_T_inv + T_log *(2.7_pr))
    k(r47f) = (570000_pr)*exp((-62216.08_pr)*R_T_inv + T_log *(0.66_pr))
    k(r48f) = (7580000_pr)*exp((-1715.44_pr)*R_T_inv)
    k(r49f) = (73400000_pr)
    k(r50f) = (30110000_pr)
    k(r51f) = (30000000_pr)
    k_0(FOr52f) = (1.27e+20_pr)*exp((-27321.52_pr)*R_T_inv + T_log *(-4.82_pr))
    k_inf(FOr52f) = (540000_pr)*exp((-15062.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr52f) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r52f) = getlindratecoeff(Tloc,k_0(FOr52f),k_inf(FOr52f),FC(FOr52f), M(mM52),Ploc)
    k_0(FOr53f) = (1.87e+19_pr)*exp((-101700.488_pr)*R_T_inv + T_log *(-3_pr))
    k_inf(FOr53f) = (68000000000000_pr)*exp((-109495.28_pr)*R_T_inv)
    FC(FOr53f) = ((1.0_pr - 0.9_pr)*exp(-Tloc/(2500_pr))) + (0.9_pr*exp(-Tloc/(1300_pr)))&
              + exp(-(1e+99_pr)/Tloc)
    k(r53f) = getlindratecoeff(Tloc,k_0(FOr53f),k_inf(FOr53f),FC(FOr53f), M(mM53),Ploc)
    k_0(FOr54f) = (8.05e+19_pr)*exp((-4107.0144_pr)*R_T_inv + T_log *(-3.75_pr))
    k_inf(FOr54f) = (2277000000_pr)*exp((-731.7816_pr)*R_T_inv + T_log *(-0.69_pr))
    FC(FOr54f) = ((1.0_pr - 0_pr)*exp(-Tloc/(570_pr))) + (0_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r54f) = getlindratecoeff(Tloc,k_0(FOr54f),k_inf(FOr54f),FC(FOr54f), M(mM54),Ploc)
    k(r55f) = (3.55_pr)*exp((-24392.72_pr)*R_T_inv + T_log *(2.4_pr))
    k(r56f) = (115_pr)*exp((-31505.52_pr)*R_T_inv + T_log *(1.9_pr))
    k(r57f) = (5.55e-10_pr)*exp((-13518.504_pr)*R_T_inv + T_log *(4.72_pr))
    k_0(FOr58f) = (1.42e+27_pr)*exp((-24137.496_pr)*R_T_inv + T_log *(-6.642_pr))
    k_inf(FOr58f) = (956.9_pr)*exp((-5669.32_pr)*R_T_inv + T_log *(1.463_pr))
    FC(FOr58f) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299_pr))) + (-0.569_pr*exp(-Tloc/(-9147_pr)))&
              + exp(-(152.4_pr)/Tloc)
    k(r58f) = getlindratecoeff(Tloc,k_0(FOr58f),k_inf(FOr58f),FC(FOr58f), M(mM58),Ploc)
    k(r59f) = (310000000_pr)*exp((-55950.54_pr)*R_T_inv + T_log *(-0.362_pr))
    k(r60f) = (18.43_pr)*exp((3014.9904_pr)*R_T_inv + T_log *(1.13_pr))
    k(r61f) = (0.0223_pr)*exp((-9269.652_pr)*R_T_inv + T_log *(2.745_pr))
    k(r62f) = (0.000976_pr)*exp((-63379.232_pr)*R_T_inv + T_log *(2.947_pr))
    k(r63f) = (8.13e-11_pr)*exp((-36968.9872_pr)*R_T_inv + T_log *(4.417_pr))
    k(r64f) = (6.098_pr)*exp((-765.672_pr)*R_T_inv + T_log *(1.88_pr))
    k_0(FOr65f) = (6.35e+19_pr)*exp((-15815.52_pr)*R_T_inv + T_log *(-4.664_pr))
    k_inf(FOr65f) = (17100_pr)*exp((-11334.456_pr)*R_T_inv + T_log *(1.266_pr))
    FC(FOr65f) = ((1.0_pr - 0.788_pr)*exp(-Tloc/(-10200_pr))) + (0.788_pr*exp(-Tloc/(1e-30_pr)))
    k(r65f) = getlindratecoeff(Tloc,k_0(FOr65f),k_inf(FOr65f),FC(FOr65f), M(mM65),Ploc)
    k(r66f) = (710000000000000_pr)*exp((-27149.976_pr)*R_T_inv + T_log *(-2.65_pr))
    k(r67f) = (20200_pr)*exp((-160.6656_pr)*R_T_inv + T_log *(0.58_pr))
    k(r68f) = (21.5_pr)*exp((-14087.528_pr)*R_T_inv + T_log *(1.19_pr))
    k(r69f) = (4.6e-05_pr)*exp((2061.8752_pr)*R_T_inv + T_log *(2.76_pr))
    k(r70f) = (0.00106_pr)*exp((-25857.12_pr)*R_T_inv + T_log *(2.39_pr))
    k(r71f) = (5.26e-07_pr)*exp((-7434.968_pr)*R_T_inv + T_log *(3.01_pr))
    k(r72f) = (2.73e+29_pr)*exp((-49454.88_pr)*R_T_inv + T_log *(-7.32_pr))
    k(r73f) = (6080000000_pr)*exp((-2701.6088_pr)*R_T_inv + T_log *(-1.31_pr))
    k(r74f) = (6.36e+29_pr)*exp((-49454.88_pr)*R_T_inv + T_log *(-7.32_pr))
    k(r75f) = (14200000000_pr)*exp((-2701.6088_pr)*R_T_inv + T_log *(-1.31_pr))
    k(r76f) = (103000_pr)*exp((3128.7952_pr)*R_T_inv + T_log *(-0.33_pr))
    k(r77f) = (5.77e+15_pr)*exp((-19966.048_pr)*R_T_inv + T_log *(-3.54_pr))
    k(r78f) = (7.25e+25_pr)*exp((-43680.96_pr)*R_T_inv + T_log *(-6.7_pr))
    k(r79f) = (53200000_pr)*exp((-1868.9928_pr)*R_T_inv + T_log *(-1.14_pr))
    k(r80f) = (96353000_pr)
    k(r81f) = (30110000_pr)
    k(r82f) = (2958_pr)*exp((-10342.848_pr)*R_T_inv + T_log *(1.28_pr))
    k(r83f) = (7.528_pr)*exp((-8811.504_pr)*R_T_inv + T_log *(1.55_pr))
    k(r84f) = (50000000_pr)*exp((-6276_pr)*R_T_inv)
    k(r85f) = (6.45e+18_pr)*exp((-68770.316_pr)*R_T_inv + T_log *(-2.52_pr))
    k_0(FOr86f) = (9.52e+27_pr)*exp((-172799.2_pr)*R_T_inv + T_log *(-5.07_pr))
    k_inf(FOr86f) = (2930000000000_pr)*exp((-168615.2_pr)*R_T_inv + T_log *(0.29_pr))
    FC(FOr86f) = ((1.0_pr - 7.13e-17_pr)*exp(-Tloc/(1150_pr))) + (7.13e-17_pr*exp(-Tloc/(4990000000_pr)))&
              + exp(-(1790000000_pr)/Tloc)
    k(r86f) = getlindratecoeff(Tloc,k_0(FOr86f),k_inf(FOr86f),FC(FOr86f), M(mM86),Ploc)
    k(r87f) = (777_pr)*exp((-11631.52_pr)*R_T_inv + T_log *(1.45_pr))
    k(r88f) = (2000000_pr)*exp((4225.84_pr)*R_T_inv)
    k(r89f) = (100000000_pr)
    k(r90f) = (100000000_pr)
    k(r91f) = (4780000_pr)*exp((-4811.6_pr)*R_T_inv + T_log *(-0.142_pr))
    k_0(FOr92f) = (4.27e+46_pr)*exp((-40876.8432_pr)*R_T_inv + T_log *(-11.94_pr))
    k_inf(FOr92f) = (25000000_pr)
    FC(FOr92f) = ((1.0_pr - 0.175_pr)*exp(-Tloc/(1341_pr))) + (0.175_pr*exp(-Tloc/(60000_pr)))&
              + exp(-(10140_pr)/Tloc)
    k(r92f) = getlindratecoeff(Tloc,k_0(FOr92f),k_inf(FOr92f),FC(FOr92f), M(mM92),Ploc)
    k(r93f) = (2.13e+52_pr)*exp((-85353.6_pr)*R_T_inv + T_log *(-13.5_pr))
    k(r94f) = (1.3e+34_pr)*exp((-24534.1392_pr)*R_T_inv + T_log *(-8.77_pr))
    k(r95f) = (4.44e+29_pr)*exp((-58157.6_pr)*R_T_inv + T_log *(-6.55_pr))
    k(r96f) = (389000000_pr)*exp((-5178.1184_pr)*R_T_inv + T_log *(-0.42_pr))
    k(r97f) = (4.16e+18_pr)*exp((-40835.84_pr)*R_T_inv + T_log *(-3.19_pr))
    k(r98f) = (4920_pr)*exp((-2509.43768_pr)*R_T_inv + T_log *(1.02_pr))
    k(r99f) = (5.3e+23_pr)*exp((-60249.6_pr)*R_T_inv + T_log *(-4.57_pr))
    k(r100f) = (44800_pr)*exp((-5947.9744_pr)*R_T_inv + T_log *(0.6_pr))
    k(r101f) = (5.8e+75_pr)*exp((-560656_pr)*R_T_inv + T_log *(-17.2_pr))
    k(r102f) = (6.7e+54_pr)*exp((-476306.56_pr)*R_T_inv + T_log *(-11.8_pr))
    k(r103f) = (1.08e+71_pr)*exp((-522414.24_pr)*R_T_inv + T_log *(-15.9_pr))
    k(r104f) = (6.28e+42_pr)*exp((-410048.736_pr)*R_T_inv + T_log *(-8.51_pr))
    k(r105f) = (2.31e+54_pr)*exp((-91629.6_pr)*R_T_inv + T_log *(-14_pr))
    k(r106f) = (1.09e+34_pr)*exp((-26977.5952_pr)*R_T_inv + T_log *(-8.66_pr))
    k(r107f) = (2.98e+24_pr)*exp((-50208_pr)*R_T_inv + T_log *(-4.79_pr))
    k(r108f) = (817000_pr)*exp((-4956.3664_pr)*R_T_inv + T_log *(0.49_pr))
    k(r109f) = (2.01e+18_pr)*exp((-45605.6_pr)*R_T_inv + T_log *(-2.86_pr))
    k(r110f) = (997_pr)*exp((-6679.3376_pr)*R_T_inv + T_log *(1.36_pr))
    k(r111f) = (0.3644_pr)*exp((-18247.2608_pr)*R_T_inv + T_log *(2.455_pr))
    k(r112f) = (120000000000000_pr)*exp((-193267.7464_pr)*R_T_inv + T_log *(-1.67_pr))
    k(r113f) = (524000_pr)*exp((-24618.656_pr)*R_T_inv + T_log *(0.7_pr))
    k(r114f) = (4.46_pr)*exp((-4396.5472_pr)*R_T_inv + T_log *(2.072_pr))
    k(r115f) = (0.0001498_pr)*exp((-37277.348_pr)*R_T_inv + T_log *(3.381_pr))
    k(r116f) = (1.8_pr)*exp((-9352.0768_pr)*R_T_inv + T_log *(1.979_pr))
    k(r117f) = (7.45_pr)*exp((-765.672_pr)*R_T_inv + T_log *(1.88_pr))
    k(r118f) = (1.04e+43_pr)*exp((-64262.056_pr)*R_T_inv + T_log *(-11.5_pr))
    k(r119f) = (3.78e+22_pr)*exp((-23535.4184_pr)*R_T_inv + T_log *(-5.57_pr))
    k(r120f) = (2670000_pr)*exp((-22723.7224_pr)*R_T_inv + T_log *(0.47_pr))
    k(r121f) = (1e-16_pr)
    k(r122f) = (3.26e+55_pr)*exp((-84353.624_pr)*R_T_inv + T_log *(-14.94_pr))
    k(r123f) = (1.06e+24_pr)*exp((-23486.4656_pr)*R_T_inv + T_log *(-5.63_pr))
    k(r124f) = (7.67e+41_pr)*exp((-93579.344_pr)*R_T_inv + T_log *(-11.17_pr))
    k(r125f) = (2.6e+27_pr)*exp((-51948.544_pr)*R_T_inv + T_log *(-7.46_pr))
    k(r126f) = (9.86e+19_pr)*exp((-47092.5936_pr)*R_T_inv + T_log *(-3.751_pr))
    k(r127f) = (255000000000000_pr)*exp((-6550.8888_pr)*R_T_inv + T_log *(-2.608_pr))
    k(r128f) = (777000000000_pr)*exp((-9954.5728_pr)*R_T_inv + T_log *(-1.52_pr))
    k(r129f) = (1.07e+26_pr)*exp((-103022.2136_pr)*R_T_inv + T_log *(-4.56_pr))
    k(r130f) = (1.15e+25_pr)*exp((-84982.0608_pr)*R_T_inv + T_log *(-4.24_pr))
    k(r131f) = (20_pr)*exp((-20920_pr)*R_T_inv + T_log *(2_pr))
    k(r132f) = (2480000000_pr)*exp((-26928.224_pr)*R_T_inv + T_log *(-0.33_pr))
    k(r133f) = (1e-16_pr)
    k(r134f) = (4.67e+45_pr)*exp((-89286.56_pr)*R_T_inv + T_log *(-11.45_pr))
    k(r135f) = (3.32e+24_pr)*exp((-28923.992_pr)*R_T_inv + T_log *(-5.78_pr))
    k(r136f) = (1.9e+47_pr)*exp((-69981.584_pr)*R_T_inv + T_log *(-12.59_pr))
    k(r137f) = (2.81e+34_pr)*exp((-32844.4_pr)*R_T_inv + T_log *(-9.42_pr))
    k(r138f) = (126000000000000_pr)*exp((-62772.552_pr)*R_T_inv + T_log *(-1.83_pr))
    k(r139f) = (0.0123_pr)*exp((-26505.64_pr)*R_T_inv + T_log *(2.68_pr))
    k(r140f) = (2.87e+47_pr)*exp((-70512.952_pr)*R_T_inv + T_log *(-12.51_pr))
    k(r141f) = (6.93e+33_pr)*exp((-31204.272_pr)*R_T_inv + T_log *(-9.11_pr))
    k(r142f) = (3460000_pr)*exp((-22857.192_pr)*R_T_inv + T_log *(0.442_pr))
    k(r143f) = (20_pr)*exp((-4184_pr)*R_T_inv + T_log *(1.8_pr))
    k(r144f) = (20000000_pr)
    k(r145f) = (10000000_pr)
    k(r146f) = (35600000000_pr)*exp((-144766.4_pr)*R_T_inv + T_log *(0.88_pr))
    k(r147f) = (6.75e+44_pr)*exp((-155682.456_pr)*R_T_inv + T_log *(-10.07_pr))
    k(r148f) = (3.59e+88_pr)*exp((-534769.592_pr)*R_T_inv + T_log *(-20.99_pr))
    k(r149f) = (0.04395_pr)*exp((3460.5864_pr)*R_T_inv + T_log *(2.678_pr))
    k(r150f) = (31200000_pr)*exp((-156690.8_pr)*R_T_inv)
    k(r151f) = (1050000_pr)*exp((-24618.656_pr)*R_T_inv + T_log *(0.7_pr))
    k(r152f) = (0.04021_pr)*exp((-37017.94_pr)*R_T_inv + T_log *(2.5_pr))
    k(r153f) = (777000000000_pr)*exp((-9954.5728_pr)*R_T_inv + T_log *(-1.52_pr))
    k(r154f) = (1.07e+26_pr)*exp((-103022.2136_pr)*R_T_inv + T_log *(-4.56_pr))
    k(r155f) = (1.07e+24_pr)*exp((-79433.24_pr)*R_T_inv + T_log *(-4.16_pr))
    k(r156f) = (1.04e+43_pr)*exp((-64262.056_pr)*R_T_inv + T_log *(-11.5_pr))
    k(r157f) = (3.78e+22_pr)*exp((-23535.4184_pr)*R_T_inv + T_log *(-5.57_pr))
    k(r158f) = (3.26e+55_pr)*exp((-84353.624_pr)*R_T_inv + T_log *(-14.94_pr))
    k(r159f) = (1.06e+24_pr)*exp((-23486.4656_pr)*R_T_inv + T_log *(-5.63_pr))
    k(r160f) = (890000_pr)*exp((-22723.7224_pr)*R_T_inv + T_log *(0.47_pr))
    k(r161f) = (1.34_pr)*exp((-22723.7224_pr)*R_T_inv + T_log *(0.47_pr))
    k(r162f) = (1000000000_pr)*exp((-171544_pr)*R_T_inv)
    k(r163f) = (300000000_pr)*exp((-167360_pr)*R_T_inv)
    k(r164f) = (500000_pr)*exp((-12552_pr)*R_T_inv)
    k(r165f) = (10000_pr)
    k(r166f) = (100000000000000_pr)*exp((-129704_pr)*R_T_inv)
    k(r167f) = (50000000_pr)
    k(r168f) = (30000000_pr)
    k(r169f) = (7.94e+23_pr)*exp((-20338.424_pr)*R_T_inv + T_log *(-5.06_pr))
    k(r170f) = (10000000_pr)*exp((-12552_pr)*R_T_inv)
    k(r171f) = (20000000_pr)*exp((-20920_pr)*R_T_inv)
    k(r172f) = (100000000_pr)*exp((-179912_pr)*R_T_inv)
    k(r173f) = (8000000_pr)*exp((-92048_pr)*R_T_inv)
    k(r174f) = (50000000_pr)
    k(r175f) = (300000_pr)
    k(r176f) = (200000000000_pr)*exp((-133888_pr)*R_T_inv)
    k(r177f) = (9.9996_pr)*exp((13990.58472_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr178f) = (1.33e+48_pr)*exp((-24969.2752_pr)*R_T_inv + T_log *(-12_pr))
    k_inf(FOr178f) = (100000000_pr)
    FC(FOr178f) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(10970_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r178f) = getlindratecoeff(Tloc,k_0(FOr178f),k_inf(FOr178f),FC(FOr178f), M(mM178),Ploc)
    k(r179f) = (200000_pr)
    k(r180f) = (5.2e+60_pr)*exp((-381162.4_pr)*R_T_inv + T_log *(-13.93_pr))
    k(r181f) = (1800000_pr)*exp((-12552_pr)*R_T_inv)
    k(r182f) = (4.7_pr)*exp((-33712.70552_pr)*R_T_inv + T_log *(2_pr))
    k(r183f) = (18.8_pr)*exp((-33712.70552_pr)*R_T_inv + T_log *(2_pr))
    k(r184f) = (141_pr)*exp((-33712.70552_pr)*R_T_inv + T_log *(2_pr))
    k(r185f) = (24_pr)*exp((-39563.40192_pr)*R_T_inv + T_log *(2_pr))
    k(r186f) = (18.8_pr)*exp((-33712.70552_pr)*R_T_inv + T_log *(2_pr))
    k(r187f) = (12000_pr)*exp((-9656.0444_pr)*R_T_inv + T_log *(1_pr))
    k(r188f) = (18_pr)*exp((-18367.25792_pr)*R_T_inv + T_log *(2_pr))
    k(r189f) = (1500000_pr)
    k(r190f) = (1500000_pr)
    k(r191f) = (43407.1_pr)*exp((-111831.50008_pr)*R_T_inv + T_log *(2.509943_pr))
    k(r192f) = (8.726576e+17_pr)*exp((-133122.70456_pr)*R_T_inv + T_log *(-1.402041_pr))
    k(r193f) = (56554380000000_pr)*exp((-106786.68392_pr)*R_T_inv + T_log *(-0.7629124_pr))
    k(r194f) = (141385900000000_pr)*exp((-110970.68392_pr)*R_T_inv + T_log *(-0.7629124_pr))
    k(r195f) = (8.070114e+22_pr)*exp((-138804.6184_pr)*R_T_inv + T_log *(-3.078473_pr))
    k(r196f) = (40000000000000_pr)*exp((-128867.2_pr)*R_T_inv)
    k(r197f) = (2e+17_pr)*exp((-326352_pr)*R_T_inv)
    k(r198f) = (1e+16_pr)*exp((-297064_pr)*R_T_inv)
    k(r199f) = (25000000_pr)*exp((-8368_pr)*R_T_inv)
    k(r200f) = (7800000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r201f) = (21600000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r202f) = (15000000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r203f) = (15600000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r204f) = (200000000000000_pr)*exp((-158992_pr)*R_T_inv)
    k(r205f) = (134000_pr)*exp((-10460_pr)*R_T_inv)
    k(r206f) = (0.1_pr)*exp((4707_pr)*R_T_inv + T_log *(2.56_pr))
    k(r207f) = (27_pr)*exp((-16736_pr)*R_T_inv + T_log *(2_pr))
    k(r208f) = (18_pr)*exp((-16736_pr)*R_T_inv + T_log *(2_pr))
    k(r209f) = (9000_pr)*exp((-209.2_pr)*R_T_inv + T_log *(1_pr))
    k(r210f) = (500000_pr)*exp((-73220_pr)*R_T_inv)
    k(r211f) = (4e+17_pr)*exp((-324260_pr)*R_T_inv)
    k(r212f) = (20000000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r213f) = (2.5e+17_pr)*exp((-322168_pr)*R_T_inv)
    k(r214f) = (20000000000000_pr)*exp((-121336_pr)*R_T_inv)
    k(r215f) = (50000000000000_pr)*exp((-121336_pr)*R_T_inv)
    k(r216f) = (250000_pr)*exp((-14644_pr)*R_T_inv)
    k(r217f) = (39_pr)*exp((-16736_pr)*R_T_inv + T_log *(2_pr))
    k(r218f) = (19500_pr)*exp((-209.2_pr)*R_T_inv + T_log *(1_pr))
    k(r219f) = (48_pr)*exp((-16736_pr)*R_T_inv + T_log *(2_pr))
    k(r220f) = (24000_pr)*exp((-209.2_pr)*R_T_inv + T_log *(1_pr))
    k(r1b) = (0.0268129129822277_pr)*exp((-20179.429910244_pr)*R_T_inv + T_log *(2.64882382000562_pr))
    k(r2b) = (263852571.883227_pr)*exp((-91402.3739876693_pr)*R_T_inv + T_log *(-0.0340885799754926_pr))
    k(r3b) = (429670.340257317_pr)*exp((6199.28059980642_pr)*R_T_inv + T_log *(0.356901390784223_pr))
    k(r4b) = (6.18795866616868e+16_pr)*exp((-495233.816255088_pr)*R_T_inv + T_log *(-1.71856505265788_pr))
    k(r5b) = (5.87040322418415_pr)*exp((5597.74076760676_pr)*R_T_inv + T_log *(1.71691238054548_pr))
    k_0(FOr6b) = (21958746.1034916_pr)*exp((8433.90375290541_pr)*R_T_inv + T_log *(-1.40640375179644_pr))
    k_inf(FOr6b) = (17.637549006168_pr)*exp((8433.90367227979_pr)*R_T_inv + T_log *(1.79359623397226_pr))
    FC(FOr6b) = ((1.0_pr - 0.43_pr)*exp(-Tloc/(1e-30_pr))) + (0.43_pr*exp(-Tloc/(1e+30_pr)))
    k(r6b) = getlindratecoeff(Tloc,k_0(FOr6b),k_inf(FOr6b),FC(FOr6b), M(mM6),Ploc)
    k(r7b) = (81985.4226270192_pr)*exp((-131883.09104881_pr)*R_T_inv + T_log *(0.467598634100207_pr))
    k(r8b) = (3576260.41025708_pr)*exp((-160966.074994454_pr)*R_T_inv + T_log *(0.467598643767672_pr))
    k(r9b) = (23424.6012097025_pr)*exp((-153515.797680653_pr)*R_T_inv + T_log *(0.707432521026287_pr))
    k(r10b) = (1896.57729009101_pr)*exp((-230900.735934586_pr)*R_T_inv + T_log *(1.45470730232354_pr))
    k(r11b) = (2853338.36263151_pr)*exp((-222437.42252403_pr)*R_T_inv + T_log *(0.350531086813303_pr))
    k(r12b) = (7014143.61287375_pr)*exp((-286166.98956522_pr)*R_T_inv + T_log *(0.337618757255384_pr))
    k(r13b) = (450909233.8516_pr)*exp((-336469.380608213_pr)*R_T_inv + T_log *(0.33761875676606_pr))
    k(r14b) = (2126616576.06856_pr)*exp((-206382.397027698_pr)*R_T_inv + T_log *(-0.129979868640911_pr))
    k(r15b) = (4040572.61805334_pr)*exp((-154292.434033297_pr)*R_T_inv + T_log *(-0.129979904708465_pr))
    k_0(FOr16b) = (30700965633784.9_pr)*exp((-204493.882119868_pr)*R_T_inv + T_log *(-1.28618382217184_pr))
    k_inf(FOr16b) = (8204570130930.22_pr)*exp((-204493.882273679_pr)*R_T_inv + T_log *(0.383816150468114_pr))
    FC(FOr16b) = ((1.0_pr - 0.67_pr)*exp(-Tloc/(1e-30_pr))) + (0.67_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r16b) = getlindratecoeff(Tloc,k_0(FOr16b),k_inf(FOr16b),FC(FOr16b), M(mM16),Ploc)
    k(r17b) = (123636.34053641_pr)*exp((-104907.326437622_pr)*R_T_inv + T_log *(0.915171820757095_pr))
    k(r18b) = (10146463705640.4_pr)*exp((-107783.826453895_pr)*R_T_inv + T_log *(-1.8018281822741_pr))
    k_0(FOr19b) = (3.59990899730388e+29_pr)*exp((-449920.623953222_pr)*R_T_inv + T_log *(-4.82231424864122_pr))
    k_inf(FOr19b) = (1.84350280771922e+18_pr)*exp((-441314.136359572_pr)*R_T_inv + T_log *(-0.692314320656943_pr))
    FC(FOr19b) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r19b) = getlindratecoeff(Tloc,k_0(FOr19b),k_inf(FOr19b),FC(FOr19b), M(mM19),Ploc)
    k(r20b) = (0.00124142581424412_pr)*exp((-33477.9464875461_pr)*R_T_inv + T_log *(2.87783774495054_pr))
    k(r21b) = (0.000710081476063557_pr)*exp((-64685.1124875859_pr)*R_T_inv + T_log *(2.94374916306033_pr))
    k(r22b) = (9.54322922506894_pr)*exp((-222573.734365095_pr)*R_T_inv + T_log *(2.2238694805159_pr))
    k(r23b) = (9567213.75472332_pr)*exp((-39916.6129144906_pr)*R_T_inv + T_log *(-0.0400605594622435_pr))
    k(r25b) = (1270217.48854075_pr)*exp((-779608.576260738_pr)*R_T_inv + T_log *(0.454685011329288_pr))
    k(r26b) = (8760112621.08348_pr)*exp((-66235.5570469003_pr)*R_T_inv + T_log *(-0.468831459957723_pr))
    k(r27b) = (5786247.14314582_pr)*exp((-318528.945041323_pr)*R_T_inv + T_log *(-0.051474977324336_pr))
    k(r29b) = (4942560.35474627_pr)*exp((-12213.5306892411_pr)*R_T_inv + T_log *(0.165903503798759_pr))
    k(r30b) = (207154242.293403_pr)*exp((-306185.712090001_pr)*R_T_inv + T_log *(-0.196202238550358_pr))
    k(r31b) = (67045866.4794991_pr)*exp((-96656.1894991135_pr)*R_T_inv + T_log *(0.216970988762198_pr))
    k(r32b) = (25989323.5013698_pr)*exp((-577322.727089117_pr)*R_T_inv + T_log *(0.167135866658773_pr))
    k(r33b) = (6958708830.64264_pr)*exp((1656.14419284427_pr)*R_T_inv + T_log *(-0.883521887631753_pr))
    k(r34b) = (5.23652860870369e-07_pr)*exp((-252043.118864436_pr)*R_T_inv + T_log *(3.47735725614659_pr))
    k(r35b) = (2914435769.2793_pr)*exp((-287718.679745721_pr)*R_T_inv + T_log *(-0.112544185190428_pr))
    k(r36b) = (25425759359.6992_pr)*exp((-3334.70488018401_pr)*R_T_inv + T_log *(-1.08325710915205_pr))
    k(r37b) = (126735760.337053_pr)*exp((-886.783490260695_pr)*R_T_inv + T_log *(0.0225392904300487_pr))
    k(r38b) = (0.00131804649228375_pr)*exp((-50053.6651422291_pr)*R_T_inv + T_log *(2.96268230523803_pr))
    k(r39b) = (80962051.105837_pr)*exp((-99413.897865177_pr)*R_T_inv + T_log *(-0.263990718230607_pr))
    k(r40b) = (334545343.951444_pr)*exp((-79883.641942173_pr)*R_T_inv + T_log *(-0.702614599362838_pr))
    k(r41b) = (53394307.667233_pr)*exp((-100874.769731261_pr)*R_T_inv + T_log *(0.297385438180316_pr))
    k(r42b) = (1.07267371978274e-27_pr)*exp((-92825.2169186665_pr)*R_T_inv + T_log *(10.2273479883772_pr))
    k(r43b) = (680741988.729528_pr)*exp((-60217.361480601_pr)*R_T_inv + T_log *(-0.188768608357115_pr))
    k(r44b) = (0.805398631753326_pr)*exp((-76755.9627160875_pr)*R_T_inv + T_log *(2.08293855052652_pr))
    k(r45b) = (6.60986228516293_pr)*exp((-123033.895754165_pr)*R_T_inv + T_log *(1.77885013696305_pr))
    k(r46b) = (0.0337253216668472_pr)*exp((-45095.1166789496_pr)*R_T_inv + T_log *(2.38125150740495_pr))
    k(r47b) = (0.0980182552035039_pr)*exp((3068.31801111421_pr)*R_T_inv + T_log *(0.924785612195135_pr))
    k(r48b) = (2299874.80848002_pr)*exp((-140924.924258992_pr)*R_T_inv + T_log *(0.20860176251745_pr))
    k(r49b) = (3704421.33388246_pr)*exp((-367793.204365344_pr)*R_T_inv + T_log *(0.580309119577134_pr))
    k(r50b) = (9154239.861871_pr)*exp((-429949.418342745_pr)*R_T_inv + T_log *(0.546220541312997_pr))
    k_0(FOr52b) = (1.01141311054869e+27_pr)*exp((-151931.760368305_pr)*R_T_inv + T_log *(-5.17356925670052_pr))
    k_inf(FOr52b) = (4300496639801.57_pr)*exp((-139672.640359683_pr)*R_T_inv + T_log *(0.100430744830484_pr))
    FC(FOr52b) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r52b) = getlindratecoeff(Tloc,k_0(FOr52b),k_inf(FOr52b),FC(FOr52b), M(mM52),Ploc)
    k_0(FOr53b) = (25955672254.5419_pr)*exp((-13048.0066850761_pr)*R_T_inv + T_log *(-2.21646816915867_pr))
    k_inf(FOr53b) = (94384.2620487426_pr)*exp((-20842.7986807314_pr)*R_T_inv + T_log *(0.783531831836205_pr))
    FC(FOr53b) = ((1.0_pr - 0.9_pr)*exp(-Tloc/(2500_pr))) + (0.9_pr*exp(-Tloc/(1300_pr)))&
              + exp(-(1e+99_pr)/Tloc)
    k(r53b) = getlindratecoeff(Tloc,k_0(FOr53b),k_inf(FOr53b),FC(FOr53b), M(mM53),Ploc)
    k_0(FOr54b) = (1.11123147930692e+33_pr)*exp((-385100.131301753_pr)*R_T_inv + T_log *(-4.82877955786906_pr))
    k_inf(FOr54b) = (3.14319682346887e+22_pr)*exp((-381724.89831826_pr)*R_T_inv + T_log *(-1.76877952525485_pr))
    FC(FOr54b) = ((1.0_pr - 0_pr)*exp(-Tloc/(570_pr))) + (0_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r54b) = getlindratecoeff(Tloc,k_0(FOr54b),k_inf(FOr54b),FC(FOr54b), M(mM54),Ploc)
    k(r55b) = (0.00208024595964166_pr)*exp((-30207.6230252537_pr)*R_T_inv + T_log *(2.78474632748836_pr))
    k(r56b) = (0.127674443424973_pr)*exp((-43466.7212600404_pr)*R_T_inv + T_log *(2.30592248207337_pr))
    k(r57b) = (3.04752007462381e-10_pr)*exp((-32113.7665338237_pr)*R_T_inv + T_log *(4.74808477892966_pr))
    k_0(FOr58b) = (3.03267979124316e+33_pr)*exp((-173619.466062885_pr)*R_T_inv + T_log *(-6.73052245269828_pr))
    k_inf(FOr58b) = (2043641519.16723_pr)*exp((-155151.289978989_pr)*R_T_inv + T_log *(1.37447756229616_pr))
    FC(FOr58b) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299_pr))) + (-0.569_pr*exp(-Tloc/(-9147_pr)))&
              + exp(-(152.4_pr)/Tloc)
    k(r58b) = getlindratecoeff(Tloc,k_0(FOr58b),k_inf(FOr58b),FC(FOr58b), M(mM58),Ploc)
    k(r59b) = (16187661201211.9_pr)*exp((-15827.2554928885_pr)*R_T_inv + T_log *(-1.35038053226632_pr))
    k(r60b) = (15.2261456870438_pr)*exp((-51996.9217771678_pr)*R_T_inv + T_log *(1.16233860921817_pr))
    k(r61b) = (9.30041263292953e-05_pr)*exp((-42574.7994147688_pr)*R_T_inv + T_log *(3.16798747708243_pr))
    k(r62b) = (0.000334200846906747_pr)*exp((-41162.2269317263_pr)*R_T_inv + T_log *(3.02623831322536_pr))
    k(r63b) = (2.78386579188197e-11_pr)*exp((-14751.9821662472_pr)*R_T_inv + T_log *(4.4962383071406_pr))
    k(r64b) = (743.069962252979_pr)*exp((-72524.0291862548_pr)*R_T_inv + T_log *(1.13590253772062_pr))
    k_0(FOr65b) = (1.63324758795428e+24_pr)*exp((-163561.838228171_pr)*R_T_inv + T_log *(-4.21321884186752_pr))
    k_inf(FOr65b) = (439819349.011031_pr)*exp((-159080.774096075_pr)*R_T_inv + T_log *(1.71678118184675_pr))
    FC(FOr65b) = ((1.0_pr - 0.788_pr)*exp(-Tloc/(-10200_pr))) + (0.788_pr*exp(-Tloc/(1e-30_pr)))
    k(r65b) = getlindratecoeff(Tloc,k_0(FOr65b),k_inf(FOr65b),FC(FOr65b), M(mM65),Ploc)
    k(r66b) = (8.92361229065844e+17_pr)*exp((-63749.7931719251_pr)*R_T_inv + T_log *(-3.47309593465077_pr))
    k(r67b) = (25388304.1144654_pr)*exp((-36760.4827352798_pr)*R_T_inv + T_log *(-0.243095928124836_pr))
    k(r68b) = (1474.90146450007_pr)*exp((-70835.091991644_pr)*R_T_inv + T_log *(0.683035001752591_pr))
    k(r69b) = (0.00315560336468553_pr)*exp((-54685.6888457823_pr)*R_T_inv + T_log *(2.25303499226057_pr))
    k(r70b) = (0.0162234026239117_pr)*exp((-333413.249061909_pr)*R_T_inv + T_log *(1.9909106953077_pr))
    k(r71b) = (8.05048126304841e-06_pr)*exp((-314991.09709524_pr)*R_T_inv + T_log *(2.61091069003308_pr))
    k(r72b) = (1.52374998881344e+28_pr)*exp((-412691.200641824_pr)*R_T_inv + T_log *(-6.95681341334233_pr))
    k(r73b) = (339355345.209079_pr)*exp((-365937.929521142_pr)*R_T_inv + T_log *(-0.946813427154175_pr))
    k(r76b) = (712245.462129313_pr)*exp((-383475.608681635_pr)*R_T_inv + T_log *(-0.485559620475093_pr))
    k(r77b) = (3.98995786892095e+16_pr)*exp((-406570.451929645_pr)*R_T_inv + T_log *(-3.69555962964043_pr))
    k(r78b) = (3.61133457399296e+27_pr)*exp((-586672.05892159_pr)*R_T_inv + T_log *(-6.75296455889297_pr))
    k(r79b) = (2649972321.50734_pr)*exp((-544860.091698195_pr)*R_T_inv + T_log *(-1.19296455482641_pr))
    k(r80b) = (1099458773.33563_pr)*exp((-285331.283819581_pr)*R_T_inv + T_log *(-0.135257591868737_pr))
    k(r81b) = (2069721455.86566_pr)*exp((-347487.498174095_pr)*R_T_inv + T_log *(-0.169346236674931_pr))
    k(r82b) = (21488.6207428614_pr)*exp((-91779.5920268907_pr)*R_T_inv + T_log *(0.976597733288695_pr))
    k(r83b) = (5075.63984781476_pr)*exp((-107338.551301309_pr)*R_T_inv + T_log *(0.950443188357544_pr))
    k(r84b) = (266791.079194657_pr)*exp((-641607.354397729_pr)*R_T_inv + T_log *(0.737433034335535_pr))
    k(r85b) = (16577358.8699098_pr)*exp((-19287.725674161_pr)*R_T_inv + T_log *(-1.29753100587307_pr))
    k_0(FOr86b) = (1.38195287406732e+15_pr)*exp((-145863.602000382_pr)*R_T_inv + T_log *(-3.45638773310014_pr))
    k_inf(FOr86b) = (0.42532795503191_pr)*exp((-141679.60203859_pr)*R_T_inv + T_log *(1.9036122601444_pr))
    FC(FOr86b) = ((1.0_pr - 7.13e-17_pr)*exp(-Tloc/(1150_pr))) + (7.13e-17_pr*exp(-Tloc/(4990000000_pr)))&
              + exp(-(1790000000_pr)/Tloc)
    k(r86b) = getlindratecoeff(Tloc,k_0(FOr86b),k_inf(FOr86b),FC(FOr86b), M(mM86),Ploc)
    k(r87b) = (0.0014348159353116_pr)*exp((-140670.914623578_pr)*R_T_inv + T_log *(2.93395292456293_pr))
    k(r88b) = (9988.5396984917_pr)*exp((-110780.194013003_pr)*R_T_inv + T_log *(0.673492232451349_pr))
    k(r90b) = (72.2844179606601_pr)*exp((-73747.8428673762_pr)*R_T_inv + T_log *(1.63545364099356_pr))
    k_0(FOr92b) = (3.77899225018767e+60_pr)*exp((-461754.515910601_pr)*R_T_inv + T_log *(-13.2925194710715_pr))
    k_inf(FOr92b) = (2.21252463663223e+21_pr)*exp((-420877.672677596_pr)*R_T_inv + T_log *(-1.35251946519998_pr))
    FC(FOr92b) = ((1.0_pr - 0.175_pr)*exp(-Tloc/(1341_pr))) + (0.175_pr*exp(-Tloc/(60000_pr)))&
              + exp(-(10140_pr)/Tloc)
    k(r92b) = getlindratecoeff(Tloc,k_0(FOr92b),k_inf(FOr92b),FC(FOr92b), M(mM92),Ploc)
    k(r93b) = (1.63323810608449e+65_pr)*exp((-543615.763169701_pr)*R_T_inv + T_log *(-14.8642748715874_pr))
    k(r94b) = (9.96812032271198e+46_pr)*exp((-482796.302403421_pr)*R_T_inv + T_log *(-10.1342748770509_pr))
    k(r95b) = (6.32731173472858e+36_pr)*exp((-150343.215546423_pr)*R_T_inv + T_log *(-8.33222750890367_pr))
    k(r96b) = (5.54352310883203e+15_pr)*exp((-97363.7339451751_pr)*R_T_inv + T_log *(-2.2022275086623_pr))
    k(r97b) = (3.6042521083341e+17_pr)*exp((-78220.3305547411_pr)*R_T_inv + T_log *(-3.20175541722593_pr))
    k(r98b) = (426.272129451274_pr)*exp((-39893.9282433165_pr)*R_T_inv + T_log *(1.00824458120861_pr))
    k(r99b) = (8.71746679736693e+31_pr)*exp((-115050.724993761_pr)*R_T_inv + T_log *(-6.34047209200654_pr))
    k(r100b) = (7368726541737.97_pr)*exp((-60749.0993829691_pr)*R_T_inv + T_log *(-1.17047209006873_pr))
    k(r101b) = (6.55359924488163e+61_pr)*exp((-139778.327311792_pr)*R_T_inv + T_log *(-15.8474805328961_pr))
    k(r102b) = (7.57053732680227e+40_pr)*exp((-55428.8873383033_pr)*R_T_inv + T_log *(-10.4474805374731_pr))
    k(r103b) = (2.00719744128114e+65_pr)*exp((-156337.692339244_pr)*R_T_inv + T_log *(-16.317952630823_pr))
    k(r104b) = (1.16714804413738e+37_pr)*exp((-43972.1882769979_pr)*R_T_inv + T_log *(-8.92795262000112_pr))
    k(r105b) = (1.01582837781262e+62_pr)*exp((-532955.277758806_pr)*R_T_inv + T_log *(-13.9863439471663_pr))
    k(r106b) = (4.79330255096323e+41_pr)*exp((-468303.272930833_pr)*R_T_inv + T_log *(-8.64634394218893_pr))
    k(r107b) = (2.43551596642155e+26_pr)*exp((-125457.130126051_pr)*R_T_inv + T_log *(-5.19429658295112_pr))
    k(r108b) = (66772365.3285378_pr)*exp((-80205.4965048085_pr)*R_T_inv + T_log *(0.0857034208139897_pr))
    k(r109b) = (998748907896.382_pr)*exp((-66053.6050909217_pr)*R_T_inv + T_log *(-1.49382448360071_pr))
    k(r110b) = (0.000495399335383525_pr)*exp((-27127.3426933531_pr)*R_T_inv + T_log *(2.72617551602296_pr))
    k(r111b) = (0.198763602367995_pr)*exp((-85248.3154265063_pr)*R_T_inv + T_log *(2.35257089245903_pr))
    k(r112b) = (393505148865639_pr)*exp((-31685.080937123_pr)*R_T_inv + T_log *(-2.14413646702753_pr))
    k(r113b) = (150858.642329758_pr)*exp((-85473.4127590162_pr)*R_T_inv + T_log *(0.576394673538386_pr))
    k(r114b) = (14.6548254074002_pr)*exp((-133553.815947636_pr)*R_T_inv + T_log *(1.93548228930381_pr))
    k(r115b) = (9.99759911647743e-07_pr)*exp((-29029.2727237954_pr)*R_T_inv + T_log *(3.68286743633537_pr))
    k(r116b) = (0.0723675016914409_pr)*exp((-63260.2154245226_pr)*R_T_inv + T_log *(2.24677887174509_pr))
    k(r117b) = (0.000191681051448594_pr)*exp((-111800.541094326_pr)*R_T_inv + T_log *(2.70731555396415_pr))
    k(r118b) = (2.15681178783546e+49_pr)*exp((-205538.050666252_pr)*R_T_inv + T_log *(-11.5963088112973_pr))
    k(r119b) = (7.83918138307051e+28_pr)*exp((-164811.413073553_pr)*R_T_inv + T_log *(-5.66630881257756_pr))
    k(r120b) = (12.7892561176751_pr)*exp((-63774.7187732069_pr)*R_T_inv + T_log *(1.68096689830787_pr))
    k(r121b) = (4.78998344374439e-22_pr)*exp((-41050.9963543205_pr)*R_T_inv + T_log *(1.21096690156334_pr))
    k(r122b) = (6.76077538539168e+61_pr)*exp((-225629.618663366_pr)*R_T_inv + T_log *(-15.0363088107915_pr))
    k(r123b) = (2.19828899452781e+30_pr)*exp((-164762.460285456_pr)*R_T_inv + T_log *(-5.72630881464566_pr))
    k(r124b) = (3.32078122713074e+53_pr)*exp((-193804.342319828_pr)*R_T_inv + T_log *(-12.4772757141749_pr))
    k(r125b) = (1.12568862223624e+39_pr)*exp((-152173.542365621_pr)*R_T_inv + T_log *(-8.767275722332_pr))
    k(r126b) = (1.63584602153115e+21_pr)*exp((-81906.1642954428_pr)*R_T_inv + T_log *(-4.03682338054791_pr))
    k(r127b) = (473181759848.257_pr)*exp((-374433.406843521_pr)*R_T_inv + T_log *(-1.83632127873384_pr))
    k(r128b) = (4171459520524.85_pr)*exp((-61039.0248901059_pr)*R_T_inv + T_log *(-1.86630063059446_pr))
    k(r129b) = (13616785032.2027_pr)*exp((-10774.5531899344_pr)*R_T_inv + T_log *(-2.19268618677673_pr))
    k(r130b) = (2.1823739995559e+18_pr)*exp((-37119.1835274763_pr)*R_T_inv + T_log *(-4.46488680397383_pr))
    k(r131b) = (0.80822069983239_pr)*exp((-138062.019955325_pr)*R_T_inv + T_log *(2.07468392393173_pr))
    k(r132b) = (674054865.058928_pr)*exp((-31979.7347805199_pr)*R_T_inv + T_log *(-0.235164858112746_pr))
    k(r133b) = (2.71796306154855e-17_pr)*exp((-5051.51075313179_pr)*R_T_inv + T_log *(0.094835146873703_pr))
    k(r134b) = (1.60872221503267e+53_pr)*exp((-321994.741118883_pr)*R_T_inv + T_log *(-11.7519414774484_pr))
    k(r135b) = (1.14367403768075e+32_pr)*exp((-261632.173119431_pr)*R_T_inv + T_log *(-6.0819414774895_pr))
    k(r136b) = (8.00835018581059e+52_pr)*exp((-227440.634981048_pr)*R_T_inv + T_log *(-12.4876448924009_pr))
    k(r137b) = (1.1843929017728e+40_pr)*exp((-190303.451017254_pr)*R_T_inv + T_log *(-9.31764489881628_pr))
    k(r138b) = (1025987350.05572_pr)*exp((-92933.2899224551_pr)*R_T_inv + T_log *(-0.812250548067382_pr))
    k(r139b) = (1.00155912853939e-07_pr)*exp((-56666.3779578022_pr)*R_T_inv + T_log *(3.69774944561208_pr))
    k(r140b) = (4.45069459581437e+53_pr)*exp((-222920.492216762_pr)*R_T_inv + T_log *(-12.502480037198_pr))
    k(r141b) = (1.07467987142443e+40_pr)*exp((-183611.812187638_pr)*R_T_inv + T_log *(-9.10248003197475_pr))
    k(r142b) = (103.658285404572_pr)*exp((-47966.4191889858_pr)*R_T_inv + T_log *(1.36491430170789_pr))
    k(r143b) = (0.000767970625641617_pr)*exp((-499237.229736376_pr)*R_T_inv + T_log *(3.09663478263827_pr))
    k(r144b) = (829.074186482758_pr)*exp((-139121.326228441_pr)*R_T_inv + T_log *(1.2706914008602_pr))
    k(r145b) = (1487.67804543101_pr)*exp((-49569.4535813765_pr)*R_T_inv + T_log *(1.08580493033398_pr))
    k(r146b) = (141389805741.382_pr)*exp((-166168.408145878_pr)*R_T_inv + T_log *(0.904271119578402_pr))
    k(r147b) = (2.40425298942181e+30_pr)*exp((-42396.7920147055_pr)*R_T_inv + T_log *(-8.14809889500124_pr))
    k(r148b) = (7.83116713653773e+84_pr)*exp((-182714.6645408_pr)*R_T_inv + T_log *(-22.4465244879731_pr))
    k(r149b) = (16.9500277913285_pr)*exp((-139718.30261799_pr)*R_T_inv + T_log *(1.50291042028881_pr))
    k(r150b) = (12008519441.8943_pr)*exp((-9129.75472876181_pr)*R_T_inv + T_log *(-1.51270832203533_pr))
    k(r151b) = (35480846.3065202_pr)*exp((-99495.0330593657_pr)*R_T_inv + T_log *(-0.462177200932148_pr))
    k(r152b) = (6.04175176852994e+15_pr)*exp((-258811.285985566_pr)*R_T_inv + T_log *(-0.472795529137085_pr))
    k(r153b) = (1107108048.08054_pr)*exp((-45754.397988922_pr)*R_T_inv + T_log *(-0.382742024326596_pr))
    k(r154b) = (23672720828.2308_pr)*exp((-16526.0108992787_pr)*R_T_inv + T_log *(-2.37156639402373_pr))
    k(r155b) = (497694336.62899_pr)*exp((-66476.9143252097_pr)*R_T_inv + T_log *(-2.09049511064508_pr))
    k(r156b) = (7.85702779772584e+49_pr)*exp((-202484.175761909_pr)*R_T_inv + T_log *(-11.7896189892085_pr))
    k(r157b) = (2.85572745796542e+29_pr)*exp((-161757.538174134_pr)*R_T_inv + T_log *(-5.85961899134663_pr))
    k(r158b) = (9.78161689934835e+62_pr)*exp((-243977.751906337_pr)*R_T_inv + T_log *(-15.2053478694185_pr))
    k(r159b) = (3.18052585245848e+31_pr)*exp((-183110.593532338_pr)*R_T_inv + T_log *(-5.89534787398958_pr))
    k(r160b) = (0.0239492195875729_pr)*exp((-47660.178205628_pr)*R_T_inv + T_log *(2.10228211068689_pr))
    k(r161b) = (3.60583752773728e-08_pr)*exp((-47660.1781997025_pr)*R_T_inv + T_log *(2.10228211168262_pr))
    k(r165b) = (4681361326133.66_pr)*exp((-274649.786551834_pr)*R_T_inv + T_log *(-0.413085183123137_pr))
    k(r166b) = (295.909052623258_pr)*exp((4036.74099051907_pr)*R_T_inv + T_log *(0.958562291902718_pr))
    k(r167b) = (365472442.60408_pr)*exp((-750202.053543877_pr)*R_T_inv + T_log *(0.171932833658989_pr))
    k(r168b) = (586682564585.553_pr)*exp((-323263.346066304_pr)*R_T_inv + T_log *(-0.57013892317797_pr))
    k(r169b) = (9.4415902768573e+30_pr)*exp((-403481.733962642_pr)*R_T_inv + T_log *(-4.75841433989317_pr))
    k(r176b) = (11296155519.8552_pr)*exp((-156434.992337873_pr)*R_T_inv + T_log *(0.391143270739186_pr))
    k_0(FOr178b) = (7.15624633799689e+53_pr)*exp((-391045.822828036_pr)*R_T_inv + T_log *(-11.5820473631706_pr))
    k_inf(FOr178b) = (53806377516368.5_pr)*exp((-366076.547818252_pr)*R_T_inv + T_log *(0.417952602843039_pr))
    FC(FOr178b) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(10970_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r178b) = getlindratecoeff(Tloc,k_0(FOr178b),k_inf(FOr178b),FC(FOr178b), M(mM178),Ploc)
    k(r179b) = (26.6154508443533_pr)*exp((-119346.746432544_pr)*R_T_inv + T_log *(1.75846594918124_pr))
    k(r180b) = (1.91319740672143e+61_pr)*exp((-376110.889230704_pr)*R_T_inv + T_log *(-14.0248351438786_pr))
    k(r181b) = (105847.917053358_pr)*exp((-57421.7543364525_pr)*R_T_inv + T_log *(0.209164025348686_pr))
    k(r196b) = (0.0636327909798835_pr)*exp((-29990.308156632_pr)*R_T_inv + T_log *(2.07868858614449_pr))
    k(r198b) = (18.3835502613327_pr)*exp((-20966.4566034921_pr)*R_T_inv + T_log *(1.07627966878759_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(sH2) * c(sO) 
    w(r2f) = k(r2f) * c(sH2) * c(sOH) 
    w(r3f) = k(r3f) * c(sH) * c(sO2) 
    w(r4f) = k(r4f) * c(sH) * c(sOH) * m(mM4)
    w(r5f) = k(r5f) * c(sH2O) * c(sO) 
    w(r6f) = k(r6f) * c(sH2O2) 
    w(r7f) = k(r7f) * c(sH2O2) * c(sOH) 
    w(r8f) = k(r8f) * c(sH2O2) * c(sOH) 
    w(r9f) = k(r9f) * c(sH) * c(sHO2) 
    w(r10f) = k(r10f) * c(sH) * c(sHO2) 
    w(r11f) = k(r11f) * c(sHO2) * c(sO) 
    w(r12f) = k(r12f) * c(sHO2) * c(sOH) 
    w(r13f) = k(r13f) * c(sHO2) * c(sOH) 
    w(r14f) = k(r14f) * c(sHO2)**2.0_pr 
    w(r15f) = k(r15f) * c(sHO2)**2.0_pr 
    w(r16f) = k(r16f) * c(sH) * c(sO2) 
    w(r17f) = k(r17f) * c(sCO) * c(sOH) 
    w(r18f) = k(r18f) * c(sCO) * c(sOH) 
    w(r19f) = k(r19f) * c(sCH3) * c(sH) 
    w(r20f) = k(r20f) * c(sCH4) * c(sH) 
    w(r21f) = k(r21f) * c(sCH4) * c(sOH) 
    w(r22f) = k(r22f) * c(sCH3) * c(sHO2) 
    w(r23f) = k(r23f) * cqss(sqssCH2GSG) * c(sN2) 
    w(r24f) = k(r24f) * cqss(sqssCH2GSG) * c(sO2) 
    w(r25f) = k(r25f) * cqss(sqssCH2GSG) * c(sO2) 
    w(r26f) = k(r26f) * cqss(sqssCH2GSG) * c(sH2) 
    w(r27f) = k(r27f) * cqss(sqssCH2) * c(sO2) 
    w(r28f) = k(r28f) * cqss(sqssCH2) * c(sO2) 
    w(r29f) = k(r29f) * cqss(sqssCH2) * c(sH) 
    w(r30f) = k(r30f) * cqss(sqssCH) * c(sO2) 
    w(r31f) = k(r31f) * cqss(sqssCH) * c(sH) 
    w(r32f) = k(r32f) * cqss(sqssC) * c(sO2) 
    w(r33f) = k(r33f) * c(sCH3) * c(sO2) 
    w(r34f) = k(r34f) * c(sCH3) * c(sO2) 
    w(r35f) = k(r35f) * c(sCH3) * c(sO) 
    w(r36f) = k(r36f) * c(sCH3) * c(sOH) 
    w(r37f) = k(r37f) * c(sCH3) * c(sOH) 
    w(r38f) = k(r38f) * c(sCH3) * c(sOH) 
    w(r39f) = k(r39f) * c(sCH3) * c(sHO2) 
    w(r40f) = k(r40f) * cqss(sqssCH2OH) * c(sO2) 
    w(r41f) = k(r41f) * cqss(sqssCH2OH) * c(sO2) 
    w(r42f) = k(r42f) * cqss(sqssCH3O) * c(sO2) 
    w(r43f) = k(r43f) * c(sCH2O) * c(sO2) 
    w(r44f) = k(r44f) * c(sCH2O) * c(sH) 
    w(r45f) = k(r45f) * c(sCH2O) * c(sOH) 
    w(r46f) = k(r46f) * c(sCH2O) * c(sHO2) 
    w(r47f) = k(r47f) * cqss(sqssHCO) * m(mM47)
    w(r48f) = k(r48f) * cqss(sqssHCO) * c(sO2) 
    w(r49f) = k(r49f) * c(sH) * cqss(sqssHCO) 
    w(r50f) = k(r50f) * cqss(sqssHCO) * c(sOH) 
    w(r51f) = k(r51f) * cqss(sqssHCO) * c(sHO2) 
    w(r52f) = k(r52f) * c(sCH2O) * c(sH) 
    w(r53f) = k(r53f) * cqss(sqssCH3O) 
    w(r54f) = k(r54f) * c(sCH3)**2.0_pr 
    w(r55f) = k(r55f) * c(sC2H6) * c(sO) 
    w(r56f) = k(r56f) * c(sC2H6) * c(sH) 
    w(r57f) = k(r57f) * c(sC2H6) * c(sCH3) 
    w(r58f) = k(r58f) * c(sC2H4) * c(sH) 
    w(r59f) = k(r59f) * c(sCH3)**2.0_pr 
    w(r60f) = k(r60f) * c(sC2H5) * c(sO2) 
    w(r61f) = k(r61f) * c(sC2H4) * c(sOH) 
    w(r62f) = k(r62f) * c(sC2H4) * c(sCH3) 
    w(r63f) = k(r63f) * c(sC2H4) * c(sCH3) 
    w(r64f) = k(r64f) * c(sC2H4) * c(sO) 
    w(r65f) = k(r65f) * c(sC2H2) * c(sH) 
    w(r66f) = k(r66f) * c(sC2H3) * c(sO2) 
    w(r67f) = k(r67f) * c(sC2H3) * c(sO2) 
    w(r68f) = k(r68f) * c(sC2H3) * c(sO2) 
    w(r69f) = k(r69f) * c(sC2H3) * c(sO2) 
    w(r70f) = k(r70f) * c(sC2H3) * c(sO2) 
    w(r71f) = k(r71f) * c(sC2H3) * c(sO2) 
    w(r72f) = k(r72f) * c(sC2H3) * c(sO2) 
    w(r73f) = k(r73f) * c(sC2H3) * c(sO2) 
    w(r74f) = k(r74f) * c(sC2H3) * c(sO2) 
    w(r75f) = k(r75f) * c(sC2H3) * c(sO2) 
    w(r76f) = k(r76f) * c(sC2H3) * c(sO2) 
    w(r77f) = k(r77f) * c(sC2H3) * c(sO2) 
    w(r78f) = k(r78f) * c(sC2H3) * c(sO2) 
    w(r79f) = k(r79f) * c(sC2H3) * c(sO2) 
    w(r80f) = k(r80f) * c(sC2H3) * c(sH) 
    w(r81f) = k(r81f) * c(sC2H3) * c(sOH) 
    w(r82f) = k(r82f) * c(sC2H2) * c(sO) 
    w(r83f) = k(r83f) * c(sC2H2) * c(sOH) 
    w(r84f) = k(r84f) * c(sC2H) * c(sO2) 
    w(r85f) = k(r85f) * cqss(sqssCH3CO) 
    w(r86f) = k(r86f) * cqss(sqssCH2CHO) 
    w(r87f) = k(r87f) * cqss(sqssCH2CO) * c(sH) 
    w(r88f) = k(r88f) * cqss(sqssCH2CO) * c(sOH) 
    w(r89f) = k(r89f) * cqss(sqssHCCO) * c(sOH) 
    w(r90f) = k(r90f) * c(sH) * cqss(sqssHCCO) 
    w(r91f) = k(r91f) * cqss(sqssHCCO) * c(sO2) 
    w(r92f) = k(r92f) * c(sC2H3) * c(sCH3) 
    w(r93f) = k(r93f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r94f) = k(r94f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r95f) = k(r95f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r96f) = k(r96f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r97f) = k(r97f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r98f) = k(r98f) * c(sC2H4) * cqss(sqssCH2GSG) 
    w(r99f) = k(r99f) * c(sC2H3) * c(sCH3) 
    w(r100f) = k(r100f) * c(sC2H3) * c(sCH3) 
    w(r101f) = k(r101f) * c(sC3H6) 
    w(r102f) = k(r102f) * c(sC3H6) 
    w(r103f) = k(r103f) * c(sC3H6) 
    w(r104f) = k(r104f) * c(sC3H6) 
    w(r105f) = k(r105f) * c(sC3H5XT) * c(sH) 
    w(r106f) = k(r106f) * c(sC3H5XT) * c(sH) 
    w(r107f) = k(r107f) * c(sC3H5XT) * c(sH) 
    w(r108f) = k(r108f) * c(sC3H5XT) * c(sH) 
    w(r109f) = k(r109f) * c(sC3H5XT) * c(sH) 
    w(r110f) = k(r110f) * c(sC3H5XT) * c(sH) 
    w(r111f) = k(r111f) * c(sC3H6) * c(sH) 
    w(r112f) = k(r112f) * c(sC3H6) * c(sO2) 
    w(r113f) = k(r113f) * c(sC3H6) * c(sO) 
    w(r114f) = k(r114f) * c(sC3H6) * c(sOH) 
    w(r115f) = k(r115f) * c(sC3H6) * c(sH) 
    w(r116f) = k(r116f) * c(sC3H6) * c(sOH) 
    w(r117f) = k(r117f) * c(sC3H6) * c(sO) 
    w(r118f) = k(r118f) * c(sC3H6) * c(sH) 
    w(r119f) = k(r119f) * c(sC3H6) * c(sH) 
    w(r120f) = k(r120f) * c(sC3H6) * c(sH) 
    w(r121f) = k(r121f) * c(sC3H6) * c(sH) 
    w(r122f) = k(r122f) * c(sC3H6) * c(sH) 
    w(r123f) = k(r123f) * c(sC3H6) * c(sH) 
    w(r124f) = k(r124f) * c(sC2H4) * c(sCH3) 
    w(r125f) = k(r125f) * c(sC2H4) * c(sCH3) 
    w(r126f) = k(r126f) * c(sC3H5XT) * c(sO2) 
    w(r127f) = k(r127f) * c(sC3H5XT) * c(sO2) 
    w(r128f) = k(r128f) * cqss(sqssC3H5XA) * c(sHO2) 
    w(r129f) = k(r129f) * cqss(sqssC3H5O) 
    w(r130f) = k(r130f) * cqss(sqssC3H5O) 
    w(r131f) = k(r131f) * c(sC3H4XA) * c(sOH) 
    w(r132f) = k(r132f) * c(sC3H4XA) * c(sH) 
    w(r133f) = k(r133f) * c(sC3H4XA) * c(sH) 
    w(r134f) = k(r134f) * c(sC3H4XA) * c(sH) 
    w(r135f) = k(r135f) * c(sC3H4XA) * c(sH) 
    w(r136f) = k(r136f) * c(sC3H4XA) * c(sH) 
    w(r137f) = k(r137f) * c(sC3H4XA) * c(sH) 
    w(r138f) = k(r138f) * c(sC3H4XA) * c(sH) 
    w(r139f) = k(r139f) * c(sC3H4XA) * c(sH) 
    w(r140f) = k(r140f) * c(sC3H4XP) * c(sH) 
    w(r141f) = k(r141f) * c(sC3H4XP) * c(sH) 
    w(r142f) = k(r142f) * c(sC3H4XP) * c(sH) 
    w(r143f) = k(r143f) * c(sC3H4XA) * c(sO) 
    w(r144f) = k(r144f) * cqss(sqssC3H3) * c(sO) 
    w(r145f) = k(r145f) * cqss(sqssC3H3) * c(sOH) 
    w(r146f) = k(r146f) * cqss(sqssIC4H9) 
    w(r147f) = k(r147f) * cqss(sqssIC4H9) 
    w(r148f) = k(r148f) * c(sIC4H8) 
    w(r149f) = k(r149f) * c(sIC4H8) * c(sOH) 
    w(r150f) = k(r150f) * c(sIC4H8) * c(sO2) 
    w(r151f) = k(r151f) * c(sIC4H8) * c(sO) 
    w(r152f) = k(r152f) * c(sC3H4XA) * c(sCH3) 
    w(r153f) = k(r153f) * c(sHO2) * c(sIC4H7) 
    w(r154f) = k(r154f) * cqss(sqssIC4H7O) 
    w(r155f) = k(r155f) * cqss(sqssIC4H7O) 
    w(r156f) = k(r156f) * c(sH) * c(sIC4H8) 
    w(r157f) = k(r157f) * c(sH) * c(sIC4H8) 
    w(r158f) = k(r158f) * c(sH) * c(sIC4H8) 
    w(r159f) = k(r159f) * c(sH) * c(sIC4H8) 
    w(r160f) = k(r160f) * c(sH) * c(sIC4H8) 
    w(r161f) = k(r161f) * c(sH) * c(sIC4H8) 
    w(r162f) = k(r162f) * c(sC3H4XA) * c(sO2) 
    w(r163f) = k(r163f) * cqss(sqssCH2CO) * c(sO2) 
    w(r164f) = k(r164f) * cqss(sqssCH2CHO) * c(sO2) 
    w(r165f) = k(r165f) * c(sO) * c(sOH) * m(mM165)
    w(r166f) = k(r166f) * cqss(sqssCH3COCH2) 
    w(r167f) = k(r167f) * cqss(sqssCH2) * c(sO) 
    w(r168f) = k(r168f) * cqss(sqssCH2) * c(sOH) 
    w(r169f) = k(r169f) * cqss(sqssC3H3) * c(sH) 
    w(r170f) = k(r170f) * c(sC3H6) * c(sO) 
    w(r171f) = k(r171f) * c(sCH2O) * c(sO) 
    w(r172f) = k(r172f) * c(sC2H4) * c(sO2) 
    w(r173f) = k(r173f) * c(sIC4H7) * c(sO2) 
    w(r174f) = k(r174f) * c(sIC4H7) * c(sO) 
    w(r175f) = k(r175f) * c(sHO2) * c(sIC4H7) 
    w(r176f) = k(r176f) * cqss(sqssCH2CHO) 
    w(r177f) = k(r177f) * c(sC2H3CHO) * c(sOH) 
    w(r178f) = k(r178f) * cqss(sqssC3H5XA) * c(sH) 
    w(r179f) = k(r179f) * cqss(sqssC3H5XA) * c(sIC4H7) 
    w(r180f) = k(r180f) * c(sC3H4XP) 
    w(r181f) = k(r181f) * c(sO2) * cqss(sqssTC4H9) 
    w(r182f) = k(r182f) * c(sC2H2) * c(sO) 
    w(r183f) = k(r183f) * c(sC2H4) * c(sO) 
    w(r184f) = k(r184f) * c(sC3H4XA) * c(sO) 
    w(r185f) = k(r185f) * cqss(sqssCH2CO) * c(sH) 
    w(r186f) = k(r186f) * cqss(sqssCH2CO) * c(sO) 
    w(r187f) = k(r187f) * cqss(sqssCH2CO) * c(sOH) 
    w(r188f) = k(r188f) * c(sH) * c(sIC4H8) 
    w(r189f) = k(r189f) * cqss(sqssIC4H9) * c(sO2) 
    w(r190f) = k(r190f) * c(sO2) * cqss(sqssTC4H9) 
    w(r191f) = k(r191f) * cqss(sqssC4H9O2GL1G) 
    w(r192f) = k(r192f) * cqss(sqssC4H9O2GL1G) 
    w(r193f) = k(r193f) * cqss(sqssC4H9O2GL1G) 
    w(r194f) = k(r194f) * cqss(sqssC4H9O2GL1G) 
    w(r195f) = k(r195f) * cqss(sqssC4H9O2GL1G) 
    w(r196f) = k(r196f) * cqss(sqssNEOC5H11) 
    w(r197f) = k(r197f) * c(sIC8H18) 
    w(r198f) = k(r198f) * c(sIC8H16) 
    w(r199f) = k(r199f) * c(sH) * c(sIC8H16) 
    w(r200f) = k(r200f) * cqss(sqssIC8H17) 
    w(r201f) = k(r201f) * cqss(sqssIC8H17) 
    w(r202f) = k(r202f) * cqss(sqssIC8H17) 
    w(r203f) = k(r203f) * cqss(sqssIC8H17) 
    w(r204f) = k(r204f) * cqss(sqssIC8H17) 
    w(r205f) = k(r205f) * cqss(sqssIC8H17) * c(sO2) 
    w(r206f) = k(r206f) * c(sIC8H18) * c(sOH) 
    w(r207f) = k(r207f) * c(sH) * c(sIC8H18) 
    w(r208f) = k(r208f) * c(sH) * c(sIC8H16) 
    w(r209f) = k(r209f) * c(sIC8H16) * c(sOH) 
    w(r210f) = k(r210f) * c(sIC4H7) * c(sO2) 
    w(r211f) = k(r211f) * c(sIC12H26) 
    w(r212f) = k(r212f) * cqss(sqssIC12H25) 
    w(r213f) = k(r213f) * c(sIC16H34) 
    w(r214f) = k(r214f) * cqss(sqssIC16H33) 
    w(r215f) = k(r215f) * cqss(sqssIC16H33) 
    w(r216f) = k(r216f) * cqss(sqssIC16H33) * c(sO2) 
    w(r217f) = k(r217f) * c(sH) * c(sIC12H26) 
    w(r218f) = k(r218f) * c(sIC12H26) * c(sOH) 
    w(r219f) = k(r219f) * c(sH) * c(sIC16H34) 
    w(r220f) = k(r220f) * c(sIC16H34) * c(sOH) 
    w(r1b) = k(r1b) * c(sH) * c(sOH) 
    w(r2b) = k(r2b) * c(sH) * c(sH2O) 
    w(r3b) = k(r3b) * c(sO) * c(sOH) 
    w(r4b) = k(r4b) * c(sH2O) * m(mM4)
    w(r5b) = k(r5b) * c(sOH)**2.0_pr 
    w(r6b) = k(r6b) * c(sOH)**2.0_pr 
    w(r7b) = k(r7b) * c(sH2O) * c(sHO2) 
    w(r8b) = k(r8b) * c(sH2O) * c(sHO2) 
    w(r9b) = k(r9b) * c(sOH)**2.0_pr 
    w(r10b) = k(r10b) * c(sH2) * c(sO2) 
    w(r11b) = k(r11b) * c(sO2) * c(sOH) 
    w(r12b) = k(r12b) * c(sH2O) * c(sO2) 
    w(r13b) = k(r13b) * c(sH2O) * c(sO2) 
    w(r14b) = k(r14b) * c(sH2O2) * c(sO2) 
    w(r15b) = k(r15b) * c(sH2O2) * c(sO2) 
    w(r16b) = k(r16b) * c(sHO2) 
    w(r17b) = k(r17b) * c(sCO2) * c(sH) 
    w(r18b) = k(r18b) * c(sCO2) * c(sH) 
    w(r19b) = k(r19b) * c(sCH4) 
    w(r20b) = k(r20b) * c(sCH3) * c(sH2) 
    w(r21b) = k(r21b) * c(sCH3) * c(sH2O) 
    w(r22b) = k(r22b) * c(sCH4) * c(sO2) 
    w(r23b) = k(r23b) * cqss(sqssCH2) * c(sN2) 
    w(r25b) = k(r25b) * c(sCO) * c(sH2O) 
    w(r26b) = k(r26b) * c(sCH3) * c(sH) 
    w(r27b) = k(r27b) * cqss(sqssHCO) * c(sOH) 
    w(r29b) = k(r29b) * cqss(sqssCH) * c(sH2) 
    w(r30b) = k(r30b) * cqss(sqssHCO) * c(sO) 
    w(r31b) = k(r31b) * cqss(sqssC) * c(sH2) 
    w(r32b) = k(r32b) * c(sCO) * c(sO) 
    w(r33b) = k(r33b) * cqss(sqssCH3O) * c(sO) 
    w(r34b) = k(r34b) * c(sCH2O) * c(sOH) 
    w(r35b) = k(r35b) * c(sCH2O) * c(sH) 
    w(r36b) = k(r36b) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r37b) = k(r37b) * cqss(sqssCH2OH) * c(sH) 
    w(r38b) = k(r38b) * cqss(sqssCH2) * c(sH2O) 
    w(r39b) = k(r39b) * cqss(sqssCH3O) * c(sOH) 
    w(r40b) = k(r40b) * c(sCH2O) * c(sHO2) 
    w(r41b) = k(r41b) * c(sCH2O) * c(sHO2) 
    w(r42b) = k(r42b) * c(sCH2O) * c(sHO2) 
    w(r43b) = k(r43b) * cqss(sqssHCO) * c(sHO2) 
    w(r44b) = k(r44b) * c(sH2) * cqss(sqssHCO) 
    w(r45b) = k(r45b) * c(sH2O) * cqss(sqssHCO) 
    w(r46b) = k(r46b) * c(sH2O2) * cqss(sqssHCO) 
    w(r47b) = k(r47b) * c(sCO) * c(sH) * m(mM47)
    w(r48b) = k(r48b) * c(sCO) * c(sHO2) 
    w(r49b) = k(r49b) * c(sCO) * c(sH2) 
    w(r50b) = k(r50b) * c(sCO) * c(sH2O) 
    w(r52b) = k(r52b) * cqss(sqssCH2OH) 
    w(r53b) = k(r53b) * c(sCH2O) * c(sH) 
    w(r54b) = k(r54b) * c(sC2H6) 
    w(r55b) = k(r55b) * c(sC2H5) * c(sOH) 
    w(r56b) = k(r56b) * c(sC2H5) * c(sH2) 
    w(r57b) = k(r57b) * c(sC2H5) * c(sCH4) 
    w(r58b) = k(r58b) * c(sC2H5) 
    w(r59b) = k(r59b) * c(sC2H5) * c(sH) 
    w(r60b) = k(r60b) * c(sC2H4) * c(sHO2) 
    w(r61b) = k(r61b) * c(sC2H3) * c(sH2O) 
    w(r62b) = k(r62b) * c(sC2H3) * c(sCH4) 
    w(r63b) = k(r63b) * c(sC2H3) * c(sCH4) 
    w(r64b) = k(r64b) * cqss(sqssCH2CHO) * c(sH) 
    w(r65b) = k(r65b) * c(sC2H3) 
    w(r66b) = k(r66b) * cqss(sqssCH2CHO) * c(sO) 
    w(r67b) = k(r67b) * cqss(sqssCH2CHO) * c(sO) 
    w(r68b) = k(r68b) * c(sC2H2) * c(sHO2) 
    w(r69b) = k(r69b) * c(sC2H2) * c(sHO2) 
    w(r70b) = k(r70b) * cqss(sqssCH2CO) * c(sOH) 
    w(r71b) = k(r71b) * cqss(sqssCH2CO) * c(sOH) 
    w(r72b) = k(r72b) * c(sCH2O) * cqss(sqssHCO) 
    w(r73b) = k(r73b) * c(sCH2O) * cqss(sqssHCO) 
    w(r76b) = k(r76b) * cqss(sqssCH3O) * c(sCO) 
    w(r77b) = k(r77b) * cqss(sqssCH3O) * c(sCO) 
    w(r78b) = k(r78b) * c(sCH3) * c(sCO2) 
    w(r79b) = k(r79b) * c(sCH3) * c(sCO2) 
    w(r80b) = k(r80b) * c(sC2H2) * c(sH2) 
    w(r81b) = k(r81b) * c(sC2H2) * c(sH2O) 
    w(r82b) = k(r82b) * c(sH) * cqss(sqssHCCO) 
    w(r83b) = k(r83b) * cqss(sqssCH2CO) * c(sH) 
    w(r84b) = k(r84b) * c(sCO) * cqss(sqssHCO) 
    w(r85b) = k(r85b) * c(sCH3) * c(sCO) 
    w(r86b) = k(r86b) * c(sCH3) * c(sCO) 
    w(r87b) = k(r87b) * c(sCH3) * c(sCO) 
    w(r88b) = k(r88b) * cqss(sqssCH2OH) * c(sCO) 
    w(r90b) = k(r90b) * cqss(sqssCH2GSG) * c(sCO) 
    w(r92b) = k(r92b) * c(sC3H6) 
    w(r93b) = k(r93b) * c(sC3H6) 
    w(r94b) = k(r94b) * c(sC3H6) 
    w(r95b) = k(r95b) * cqss(sqssC3H5XA) * c(sH) 
    w(r96b) = k(r96b) * cqss(sqssC3H5XA) * c(sH) 
    w(r97b) = k(r97b) * c(sC2H3) * c(sCH3) 
    w(r98b) = k(r98b) * c(sC2H3) * c(sCH3) 
    w(r99b) = k(r99b) * cqss(sqssC3H5XA) * c(sH) 
    w(r100b) = k(r100b) * cqss(sqssC3H5XA) * c(sH) 
    w(r101b) = k(r101b) * c(sC2H3) * c(sCH3) 
    w(r102b) = k(r102b) * c(sC2H3) * c(sCH3) 
    w(r103b) = k(r103b) * cqss(sqssC3H5XA) * c(sH) 
    w(r104b) = k(r104b) * cqss(sqssC3H5XA) * c(sH) 
    w(r105b) = k(r105b) * c(sC3H6) 
    w(r106b) = k(r106b) * c(sC3H6) 
    w(r107b) = k(r107b) * cqss(sqssC3H5XA) * c(sH) 
    w(r108b) = k(r108b) * cqss(sqssC3H5XA) * c(sH) 
    w(r109b) = k(r109b) * c(sC2H3) * c(sCH3) 
    w(r110b) = k(r110b) * c(sC2H3) * c(sCH3) 
    w(r111b) = k(r111b) * cqss(sqssC3H5XA) * c(sH2) 
    w(r112b) = k(r112b) * cqss(sqssC3H5XA) * c(sHO2) 
    w(r113b) = k(r113b) * cqss(sqssC3H5XA) * c(sOH) 
    w(r114b) = k(r114b) * cqss(sqssC3H5XA) * c(sH2O) 
    w(r115b) = k(r115b) * c(sC3H5XT) * c(sH2) 
    w(r116b) = k(r116b) * c(sC3H5XT) * c(sH2O) 
    w(r117b) = k(r117b) * c(sC2H5) * cqss(sqssHCO) 
    w(r118b) = k(r118b) * cqss(sqssC3H7GL1G) 
    w(r119b) = k(r119b) * cqss(sqssC3H7GL1G) 
    w(r120b) = k(r120b) * c(sC2H4) * c(sCH3) 
    w(r121b) = k(r121b) * c(sC2H4) * c(sCH3) 
    w(r122b) = k(r122b) * cqss(sqssC3H7GL1G) 
    w(r123b) = k(r123b) * cqss(sqssC3H7GL1G) 
    w(r124b) = k(r124b) * cqss(sqssC3H7GL1G) 
    w(r125b) = k(r125b) * cqss(sqssC3H7GL1G) 
    w(r126b) = k(r126b) * cqss(sqssCH3COCH2) * c(sO) 
    w(r127b) = k(r127b) * c(sCH2O) * cqss(sqssCH3CO) 
    w(r128b) = k(r128b) * cqss(sqssC3H5O) * c(sOH) 
    w(r129b) = k(r129b) * c(sC2H3) * c(sCH2O) 
    w(r130b) = k(r130b) * c(sC2H3CHO) * c(sH) 
    w(r131b) = k(r131b) * cqss(sqssC3H3) * c(sH2O) 
    w(r132b) = k(r132b) * c(sC3H4XP) * c(sH) 
    w(r133b) = k(r133b) * c(sC3H4XP) * c(sH) 
    w(r134b) = k(r134b) * cqss(sqssC3H5XA) 
    w(r135b) = k(r135b) * cqss(sqssC3H5XA) 
    w(r136b) = k(r136b) * c(sC3H5XT) 
    w(r137b) = k(r137b) * c(sC3H5XT) 
    w(r138b) = k(r138b) * c(sC2H2) * c(sCH3) 
    w(r139b) = k(r139b) * c(sC2H2) * c(sCH3) 
    w(r140b) = k(r140b) * c(sC3H5XT) 
    w(r141b) = k(r141b) * c(sC3H5XT) 
    w(r142b) = k(r142b) * c(sC2H2) * c(sCH3) 
    w(r143b) = k(r143b) * c(sC2H4) * c(sCO) 
    w(r144b) = k(r144b) * c(sC2H) * c(sCH2O) 
    w(r145b) = k(r145b) * c(sC2H3) * cqss(sqssHCO) 
    w(r146b) = k(r146b) * cqss(sqssTC4H9) 
    w(r147b) = k(r147b) * c(sC3H6) * c(sCH3) 
    w(r148b) = k(r148b) * c(sH) * c(sIC4H7) 
    w(r149b) = k(r149b) * c(sH2O) * c(sIC4H7) 
    w(r150b) = k(r150b) * c(sHO2) * c(sIC4H7) 
    w(r151b) = k(r151b) * c(sIC4H7) * c(sOH) 
    w(r152b) = k(r152b) * c(sIC4H7) 
    w(r153b) = k(r153b) * cqss(sqssIC4H7O) * c(sOH) 
    w(r154b) = k(r154b) * c(sC3H5XT) * c(sCH2O) 
    w(r155b) = k(r155b) * c(sC3H6) * cqss(sqssHCO) 
    w(r156b) = k(r156b) * cqss(sqssIC4H9) 
    w(r157b) = k(r157b) * cqss(sqssIC4H9) 
    w(r158b) = k(r158b) * cqss(sqssTC4H9) 
    w(r159b) = k(r159b) * cqss(sqssTC4H9) 
    w(r160b) = k(r160b) * c(sC3H6) * c(sCH3) 
    w(r161b) = k(r161b) * c(sC3H6) * c(sCH3) 
    w(r165b) = k(r165b) * c(sHO2) * m(mM165)
    w(r166b) = k(r166b) * cqss(sqssCH2CO) * c(sCH3) 
    w(r167b) = k(r167b) * c(sCO) * c(sH2) 
    w(r168b) = k(r168b) * c(sCH2O) * c(sH) 
    w(r169b) = k(r169b) * c(sC3H4XP) 
    w(r176b) = k(r176b) * cqss(sqssCH3CO) 
    w(r178b) = k(r178b) * c(sC3H6) 
    w(r179b) = k(r179b) * c(sC3H4XA) * c(sIC4H8) 
    w(r180b) = k(r180b) * c(sC3H4XA) 
    w(r181b) = k(r181b) * c(sHO2) * c(sIC4H8) 
    w(r196b) = k(r196b) * c(sCH3) * c(sIC4H8) 
    w(r198b) = k(r198b) * c(sIC4H7) * cqss(sqssTC4H9) 
  
    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(sH2) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         - w(r2f) &
         + w(r2b) &
         + w(r10f) &
         - w(r10b) &
         + w(r20f) &
         - w(r20b) &
         - w(r26f) &
         + w(r26b) &
         + w(r29f) &
         - w(r29b) &
         + w(r31f) &
         - w(r31b) &
         + w(r44f) &
         - w(r44b) &
         + w(r49f) &
         - w(r49b) &
         + w(r56f) &
         - w(r56b) &
         + w(r80f) &
         - w(r80b) &
         + w(r89f) &
         + w(r111f) &
         - w(r111b) &
         + w(r115f) &
         - w(r115b) &
         + w(r167f) &
         - w(r167b) &
         + w(r185f) &
         + w(r188f) &
         + w(r207f) &
         + w(r208f) &
         + w(r217f) &
         + w(r219f) 

    cdot(sH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         - w(r3f) &
         + w(r3b) &
         - w(r4f) &
         + w(r4b) &
         - w(r9f) &
         + w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r16f) &
         + w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         + w(r24f) &
         + w(r26f) &
         - w(r26b) &
         + 2.0_pr * w(r28f) &
         - w(r29f) &
         + w(r29b) &
         - w(r31f) &
         + w(r31b) &
         + w(r35f) &
         - w(r35b) &
         + w(r37f) &
         - w(r37b) &
         - w(r44f) &
         + w(r44b) &
         + w(r47f) &
         - w(r47b) &
         - w(r49f) &
         + w(r49b) &
         + w(r51f) &
         - w(r52f) &
         + w(r52b) &
         + w(r53f) &
         - w(r53b) &
         - w(r56f) &
         + w(r56b) &
         - w(r58f) &
         + w(r58b) &
         + w(r59f) &
         - w(r59b) &
         + w(r64f) &
         - w(r64b) &
         - w(r65f) &
         + w(r65b) &
         + w(r74f) &
         + w(r75f) &
         - w(r80f) &
         + w(r80b) &
         + w(r82f) &
         - w(r82b) &
         + w(r83f) &
         - w(r83b) &
         - w(r87f) &
         + w(r87b) &
         - w(r90f) &
         + w(r90b) &
         + w(r91f) &
         + w(r95f) &
         - w(r95b) &
         + w(r96f) &
         - w(r96b) &
         + w(r99f) &
         - w(r99b) &
         + w(r100f) &
         - w(r100b) &
         + w(r103f) &
         - w(r103b) &
         + w(r104f) &
         - w(r104b) &
         - w(r105f) &
         + w(r105b) &
         - w(r106f) &
         + w(r106b) &
         - w(r109f) &
         + w(r109b) &
         - w(r110f) &
         + w(r110b) &
         - w(r111f) &
         + w(r111b) &
         - w(r115f) &
         + w(r115b) &
         - w(r118f) &
         + w(r118b) &
         - w(r119f) &
         + w(r119b) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r121b) &
         - w(r122f) &
         + w(r122b) &
         - w(r123f) &
         + w(r123b) &
         + w(r130f) &
         - w(r130b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r136f) &
         + w(r136b) &
         - w(r137f) &
         + w(r137b) &
         - w(r138f) &
         + w(r138b) &
         - w(r139f) &
         + w(r139b) &
         - w(r140f) &
         + w(r140b) &
         - w(r141f) &
         + w(r141b) &
         - w(r142f) &
         + w(r142b) &
         + w(r148f) &
         - w(r148b) &
         - w(r156f) &
         + w(r156b) &
         - w(r157f) &
         + w(r157b) &
         - w(r158f) &
         + w(r158b) &
         - w(r159f) &
         + w(r159b) &
         - w(r160f) &
         + w(r160b) &
         - w(r161f) &
         + w(r161b) &
         + w(r168f) &
         - w(r168b) &
         - w(r169f) &
         + w(r169b) &
         + 2.0_pr * w(r170f) &
         + 2.0_pr * w(r171f) &
         - w(r178f) &
         + w(r178b) &
         - w(r185f) &
         - w(r188f) &
         - w(r199f) &
         + w(r204f) &
         - w(r207f) &
         - w(r208f) &
         - w(r217f) &
         - w(r219f) 

    cdot(sO2) = 0.0_pr &
         - w(r3f) &
         + w(r3b) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         + w(r12f) &
         - w(r12b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         - w(r16f) &
         + w(r16b) &
         + w(r22f) &
         - w(r22b) &
         - w(r24f) &
         - w(r25f) &
         + w(r25b) &
         - w(r27f) &
         + w(r27b) &
         - w(r28f) &
         - w(r30f) &
         + w(r30b) &
         - w(r32f) &
         + w(r32b) &
         - w(r33f) &
         + w(r33b) &
         - w(r34f) &
         + w(r34b) &
         - w(r40f) &
         + w(r40b) &
         - w(r41f) &
         + w(r41b) &
         - w(r42f) &
         + w(r42b) &
         - w(r43f) &
         + w(r43b) &
         - w(r48f) &
         + w(r48b) &
         - w(r60f) &
         + w(r60b) &
         - w(r66f) &
         + w(r66b) &
         - w(r67f) &
         + w(r67b) &
         - w(r68f) &
         + w(r68b) &
         - w(r69f) &
         + w(r69b) &
         - w(r70f) &
         + w(r70b) &
         - w(r71f) &
         + w(r71b) &
         - w(r72f) &
         + w(r72b) &
         - w(r73f) &
         + w(r73b) &
         - w(r74f) &
         - w(r75f) &
         - w(r76f) &
         + w(r76b) &
         - w(r77f) &
         + w(r77b) &
         - w(r78f) &
         + w(r78b) &
         - w(r79f) &
         + w(r79b) &
         - w(r84f) &
         + w(r84b) &
         - w(r91f) &
         - w(r112f) &
         + w(r112b) &
         - w(r126f) &
         + w(r126b) &
         - w(r127f) &
         + w(r127b) &
         - w(r150f) &
         + w(r150b) &
         - w(r162f) &
         - w(r163f) &
         - w(r164f) &
         - w(r172f) &
         - w(r181f) &
         + w(r181b) &
         - w(r189f) &
         - w(r190f) &
         + w(r191f) &
         + w(r192f) &
         - w(r205f) &
         - w(r210f) &
         - w(r216f) 

    cdot(sO) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         + w(r3f) &
         - w(r3b) &
         - w(r5f) &
         + w(r5b) &
         - w(r11f) &
         + w(r11b) &
         + w(r30f) &
         - w(r30b) &
         + w(r32f) &
         - w(r32b) &
         + w(r33f) &
         - w(r33b) &
         - w(r35f) &
         + w(r35b) &
         - w(r55f) &
         + w(r55b) &
         - w(r64f) &
         + w(r64b) &
         + w(r66f) &
         - w(r66b) &
         + w(r67f) &
         - w(r67b) &
         - w(r82f) &
         + w(r82b) &
         - w(r113f) &
         + w(r113b) &
         - w(r117f) &
         + w(r117b) &
         + w(r126f) &
         - w(r126b) &
         - w(r143f) &
         + w(r143b) &
         - w(r144f) &
         + w(r144b) &
         - w(r151f) &
         + w(r151b) &
         - w(r165f) &
         + w(r165b) &
         - w(r167f) &
         + w(r167b) &
         - w(r170f) &
         - w(r171f) &
         - w(r174f) &
         - w(r182f) &
         - w(r183f) &
         - w(r184f) &
         - w(r186f) 

    cdot(sH2O) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         + w(r4f) &
         - w(r4b) &
         - w(r5f) &
         + w(r5b) &
         + w(r7f) &
         - w(r7b) &
         + w(r8f) &
         - w(r8b) &
         + w(r12f) &
         - w(r12b) &
         + w(r13f) &
         - w(r13b) &
         + w(r21f) &
         - w(r21b) &
         + w(r25f) &
         - w(r25b) &
         + w(r36f) &
         - w(r36b) &
         + w(r38f) &
         - w(r38b) &
         + w(r45f) &
         - w(r45b) &
         + w(r50f) &
         - w(r50b) &
         + w(r61f) &
         - w(r61b) &
         + w(r81f) &
         - w(r81b) &
         + w(r114f) &
         - w(r114b) &
         + w(r116f) &
         - w(r116b) &
         + w(r131f) &
         - w(r131b) &
         + w(r149f) &
         - w(r149b) &
         + w(r177f) &
         + w(r187f) &
         + w(r206f) &
         + w(r209f) &
         + w(r218f) &
         + w(r220f) 

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         - w(r2f) &
         + w(r2b) &
         + w(r3f) &
         - w(r3b) &
         - w(r4f) &
         + w(r4b) &
         + 2.0_pr * w(r5f) &
         - 2.0_pr * w(r5b) &
         + 2.0_pr * w(r6f) &
         - 2.0_pr * w(r6b) &
         - w(r7f) &
         + w(r7b) &
         - w(r8f) &
         + w(r8b) &
         + 2.0_pr * w(r9f) &
         - 2.0_pr * w(r9b) &
         + w(r11f) &
         - w(r11b) &
         - w(r12f) &
         + w(r12b) &
         - w(r13f) &
         + w(r13b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         - w(r21f) &
         + w(r21b) &
         + w(r24f) &
         + w(r27f) &
         - w(r27b) &
         + w(r34f) &
         - w(r34b) &
         - w(r36f) &
         + w(r36b) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         + w(r39f) &
         - w(r39b) &
         - w(r45f) &
         + w(r45b) &
         - w(r50f) &
         + w(r50b) &
         + w(r51f) &
         + w(r55f) &
         - w(r55b) &
         - w(r61f) &
         + w(r61b) &
         + w(r70f) &
         - w(r70b) &
         + w(r71f) &
         - w(r71b) &
         - w(r81f) &
         + w(r81b) &
         - w(r83f) &
         + w(r83b) &
         - w(r88f) &
         + w(r88b) &
         - w(r89f) &
         + w(r113f) &
         - w(r113b) &
         - w(r114f) &
         + w(r114b) &
         - w(r116f) &
         + w(r116b) &
         + w(r128f) &
         - w(r128b) &
         - w(r131f) &
         + w(r131b) &
         - w(r145f) &
         + w(r145b) &
         - w(r149f) &
         + w(r149b) &
         + w(r151f) &
         - w(r151b) &
         + w(r153f) &
         - w(r153b) &
         + w(r163f) &
         - w(r165f) &
         + w(r165b) &
         - w(r168f) &
         + w(r168b) &
         + w(r175f) &
         - w(r177f) &
         + w(r182f) &
         + w(r183f) &
         + w(r184f) &
         + w(r186f) &
         - w(r187f) &
         + w(r194f) &
         - w(r206f) &
         - w(r209f) &
         - w(r218f) &
         - w(r220f) 

    cdot(sH2O2) = 0.0_pr &
         - w(r6f) &
         + w(r6b) &
         - w(r7f) &
         + w(r7b) &
         - w(r8f) &
         + w(r8b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r46f) &
         - w(r46b) 

    cdot(sHO2) = 0.0_pr &
         + w(r7f) &
         - w(r7b) &
         + w(r8f) &
         - w(r8b) &
         - w(r9f) &
         + w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         - w(r12f) &
         + w(r12b) &
         - w(r13f) &
         + w(r13b) &
         -2.0_pr * w(r14f) &
         + 2.0_pr * w(r14b) &
         -2.0_pr * w(r15f) &
         + 2.0_pr * w(r15b) &
         + w(r16f) &
         - w(r16b) &
         - w(r22f) &
         + w(r22b) &
         - w(r39f) &
         + w(r39b) &
         + w(r40f) &
         - w(r40b) &
         + w(r41f) &
         - w(r41b) &
         + w(r42f) &
         - w(r42b) &
         + w(r43f) &
         - w(r43b) &
         - w(r46f) &
         + w(r46b) &
         + w(r48f) &
         - w(r48b) &
         - w(r51f) &
         + w(r60f) &
         - w(r60b) &
         + w(r68f) &
         - w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r112f) &
         - w(r112b) &
         - w(r128f) &
         + w(r128b) &
         + w(r150f) &
         - w(r150b) &
         - w(r153f) &
         + w(r153b) &
         + w(r164f) &
         + w(r165f) &
         - w(r165b) &
         - w(r175f) &
         + w(r181f) &
         - w(r181b) &
         + w(r193f) &
         + w(r195f) &
         + w(r205f) &
         + w(r216f) 

    cdot(sCO) = 0.0_pr &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         + w(r24f) &
         + w(r25f) &
         - w(r25b) &
         + w(r32f) &
         - w(r32b) &
         + w(r47f) &
         - w(r47b) &
         + w(r48f) &
         - w(r48b) &
         + w(r49f) &
         - w(r49b) &
         + w(r50f) &
         - w(r50b) &
         + w(r74f) &
         + w(r75f) &
         + w(r76f) &
         - w(r76b) &
         + w(r77f) &
         - w(r77b) &
         + w(r84f) &
         - w(r84b) &
         + w(r85f) &
         - w(r85b) &
         + w(r86f) &
         - w(r86b) &
         + w(r87f) &
         - w(r87b) &
         + w(r88f) &
         - w(r88b) &
         + 2.0_pr * w(r89f) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         + w(r143f) &
         - w(r143b) &
         + w(r163f) &
         + w(r167f) &
         - w(r167b) &
         + w(r177f) 

    cdot(sCO2) = 0.0_pr &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         + w(r28f) &
         + w(r51f) &
         + w(r78f) &
         - w(r78b) &
         + w(r79f) &
         - w(r79b) &
         + w(r91f) &
         + w(r171f) 

    cdot(sCH4) = 0.0_pr &
         + w(r19f) &
         - w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r57f) &
         - w(r57b) &
         + w(r62f) &
         - w(r62b) &
         + w(r63f) &
         - w(r63b) 

    cdot(sCH3) = 0.0_pr &
         - w(r19f) &
         + w(r19b) &
         + w(r20f) &
         - w(r20b) &
         + w(r21f) &
         - w(r21b) &
         - w(r22f) &
         + w(r22b) &
         + w(r26f) &
         - w(r26b) &
         - w(r33f) &
         + w(r33b) &
         - w(r34f) &
         + w(r34b) &
         - w(r35f) &
         + w(r35b) &
         - w(r36f) &
         + w(r36b) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r39f) &
         + w(r39b) &
         -2.0_pr * w(r54f) &
         + 2.0_pr * w(r54b) &
         - w(r57f) &
         + w(r57b) &
         -2.0_pr * w(r59f) &
         + 2.0_pr * w(r59b) &
         - w(r62f) &
         + w(r62b) &
         - w(r63f) &
         + w(r63b) &
         + w(r78f) &
         - w(r78b) &
         + w(r79f) &
         - w(r79b) &
         + w(r85f) &
         - w(r85b) &
         + w(r86f) &
         - w(r86b) &
         + w(r87f) &
         - w(r87b) &
         - w(r92f) &
         + w(r92b) &
         + w(r97f) &
         - w(r97b) &
         + w(r98f) &
         - w(r98b) &
         - w(r99f) &
         + w(r99b) &
         - w(r100f) &
         + w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r102f) &
         - w(r102b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b) &
         + w(r120f) &
         - w(r120b) &
         + w(r121f) &
         - w(r121b) &
         - w(r124f) &
         + w(r124b) &
         - w(r125f) &
         + w(r125b) &
         + w(r138f) &
         - w(r138b) &
         + w(r139f) &
         - w(r139b) &
         + w(r142f) &
         - w(r142b) &
         + w(r147f) &
         - w(r147b) &
         - w(r152f) &
         + w(r152b) &
         + w(r160f) &
         - w(r160b) &
         + w(r161f) &
         - w(r161b) &
         + w(r166f) &
         - w(r166b) &
         + w(r173f) &
         + w(r196f) &
         - w(r196b) &
         + w(r203f) &
         + 0.67_pr * w(r214f) 

    cdot(sCH2O) = 0.0_pr &
         + w(r34f) &
         - w(r34b) &
         + w(r35f) &
         - w(r35b) &
         + w(r40f) &
         - w(r40b) &
         + 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(r52f) &
         + w(r52b) &
         + w(r53f) &
         - w(r53b) &
         + w(r72f) &
         - w(r72b) &
         + w(r73f) &
         - w(r73b) &
         + w(r74f) &
         + w(r75f) &
         + w(r127f) &
         - w(r127b) &
         + w(r129f) &
         - w(r129b) &
         + w(r144f) &
         - w(r144b) &
         + w(r154f) &
         - w(r154b) &
         + w(r162f) &
         + w(r168f) &
         - w(r168b) &
         - w(r171f) &
         + w(r174f) &
         + w(r175f) &
         + w(r194f) &
         + w(r210f) 

    cdot(sC2H6) = 0.0_pr &
         + w(r54f) &
         - w(r54b) &
         - w(r55f) &
         + w(r55b) &
         - w(r56f) &
         + w(r56b) &
         - w(r57f) &
         + w(r57b) 

    cdot(sC2H5) = 0.0_pr &
         + w(r55f) &
         - w(r55b) &
         + w(r56f) &
         - w(r56b) &
         + w(r57f) &
         - w(r57b) &
         + w(r58f) &
         - w(r58b) &
         + w(r59f) &
         - w(r59b) &
         - w(r60f) &
         + w(r60b) &
         + w(r117f) &
         - w(r117b) 

    cdot(sC2H4) = 0.0_pr &
         - w(r58f) &
         + w(r58b) &
         + w(r60f) &
         - w(r60b) &
         - w(r61f) &
         + w(r61b) &
         - w(r62f) &
         + w(r62b) &
         - w(r63f) &
         + w(r63b) &
         - w(r64f) &
         + w(r64b) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         - w(r96f) &
         + w(r96b) &
         - w(r97f) &
         + w(r97b) &
         - w(r98f) &
         + w(r98b) &
         + w(r120f) &
         - w(r120b) &
         + w(r121f) &
         - w(r121b) &
         - w(r124f) &
         + w(r124b) &
         - w(r125f) &
         + w(r125b) &
         + w(r143f) &
         - w(r143b) &
         - w(r172f) &
         - w(r183f) 

    cdot(sC2H3) = 0.0_pr &
         + w(r61f) &
         - w(r61b) &
         + w(r62f) &
         - w(r62b) &
         + w(r63f) &
         - w(r63b) &
         + w(r65f) &
         - w(r65b) &
         - w(r66f) &
         + w(r66b) &
         - w(r67f) &
         + w(r67b) &
         - w(r68f) &
         + w(r68b) &
         - w(r69f) &
         + w(r69b) &
         - w(r70f) &
         + w(r70b) &
         - w(r71f) &
         + w(r71b) &
         - w(r72f) &
         + w(r72b) &
         - w(r73f) &
         + w(r73b) &
         - w(r74f) &
         - w(r75f) &
         - w(r76f) &
         + w(r76b) &
         - w(r77f) &
         + w(r77b) &
         - w(r78f) &
         + w(r78b) &
         - w(r79f) &
         + w(r79b) &
         - w(r80f) &
         + w(r80b) &
         - w(r81f) &
         + w(r81b) &
         - w(r92f) &
         + w(r92b) &
         + w(r97f) &
         - w(r97b) &
         + w(r98f) &
         - w(r98b) &
         - w(r99f) &
         + w(r99b) &
         - w(r100f) &
         + w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r102f) &
         - w(r102b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b) &
         + w(r129f) &
         - w(r129b) &
         + w(r145f) &
         - w(r145b) &
         + w(r177f) &
         + w(r183f) 

    cdot(sC2H2) = 0.0_pr &
         - w(r65f) &
         + w(r65b) &
         + w(r68f) &
         - w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r80f) &
         - w(r80b) &
         + w(r81f) &
         - w(r81b) &
         - w(r82f) &
         + w(r82b) &
         - w(r83f) &
         + w(r83b) &
         + w(r138f) &
         - w(r138b) &
         + w(r139f) &
         - w(r139b) &
         + w(r142f) &
         - w(r142b) &
         - w(r182f) 

    cdot(sC2H) = 0.0_pr &
         - w(r84f) &
         + w(r84b) &
         + w(r144f) &
         - w(r144b) &
         + w(r182f) 

    cdot(sC3H6) = 0.0_pr &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         + w(r94f) &
         - w(r94b) &
         - w(r101f) &
         + w(r101b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         + w(r105f) &
         - w(r105b) &
         + w(r106f) &
         - w(r106b) &
         - w(r111f) &
         + w(r111b) &
         - w(r112f) &
         + w(r112b) &
         - w(r113f) &
         + w(r113b) &
         - w(r114f) &
         + w(r114b) &
         - w(r115f) &
         + w(r115b) &
         - w(r116f) &
         + w(r116b) &
         - w(r117f) &
         + w(r117b) &
         - w(r118f) &
         + w(r118b) &
         - w(r119f) &
         + w(r119b) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r121b) &
         - w(r122f) &
         + w(r122b) &
         - w(r123f) &
         + w(r123b) &
         + w(r147f) &
         - w(r147b) &
         + w(r155f) &
         - w(r155b) &
         + w(r160f) &
         - w(r160b) &
         + w(r161f) &
         - w(r161b) &
         - w(r170f) &
         + w(r178f) &
         - w(r178b) &
         + w(r194f) &
         + w(r200f) &
         + 0.67_pr * w(r213f) &
         + 0.56_pr * w(r214f) 

    cdot(sC3H5XT) = 0.0_pr &
         - w(r105f) &
         + w(r105b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         + w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r109f) &
         + w(r109b) &
         - w(r110f) &
         + w(r110b) &
         + w(r115f) &
         - w(r115b) &
         + w(r116f) &
         - w(r116b) &
         - w(r126f) &
         + w(r126b) &
         - w(r127f) &
         + w(r127b) &
         + w(r136f) &
         - w(r136b) &
         + w(r137f) &
         - w(r137b) &
         + w(r140f) &
         - w(r140b) &
         + w(r141f) &
         - w(r141b) &
         + w(r154f) &
         - w(r154b) &
         + w(r174f) &
         + w(r175f) 

    cdot(sC3H4XP) = 0.0_pr &
         + w(r132f) &
         - w(r132b) &
         + w(r133f) &
         - w(r133b) &
         - w(r140f) &
         + w(r140b) &
         - w(r141f) &
         + w(r141b) &
         - w(r142f) &
         + w(r142b) &
         + w(r169f) &
         - w(r169b) &
         - w(r180f) &
         + w(r180b) 

    cdot(sC3H4XA) = 0.0_pr &
         - w(r131f) &
         + w(r131b) &
         - w(r132f) &
         + w(r132b) &
         - w(r133f) &
         + w(r133b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r136f) &
         + w(r136b) &
         - w(r137f) &
         + w(r137b) &
         - w(r138f) &
         + w(r138b) &
         - w(r139f) &
         + w(r139b) &
         - w(r143f) &
         + w(r143b) &
         - w(r152f) &
         + w(r152b) &
         - w(r162f) &
         + w(r173f) &
         + w(r179f) &
         - w(r179b) &
         + w(r180f) &
         - w(r180b) &
         - w(r184f) &
         + 0.5_pr * w(r208f) &
         + 0.5_pr * w(r209f) 

    cdot(sC2H3CHO) = 0.0_pr &
         + w(r130f) &
         - w(r130b) &
         + w(r170f) &
         - w(r177f) 

    cdot(sIC4H8) = 0.0_pr &
         - w(r148f) &
         + w(r148b) &
         - w(r149f) &
         + w(r149b) &
         - w(r150f) &
         + w(r150b) &
         - w(r151f) &
         + w(r151b) &
         - w(r156f) &
         + w(r156b) &
         - w(r157f) &
         + w(r157b) &
         - w(r158f) &
         + w(r158b) &
         - w(r159f) &
         + w(r159b) &
         - w(r160f) &
         + w(r160b) &
         - w(r161f) &
         + w(r161b) &
         + w(r179f) &
         - w(r179b) &
         + w(r181f) &
         - w(r181b) &
         - w(r188f) &
         + w(r193f) &
         + w(r195f) &
         + w(r196f) &
         - w(r196b) &
         + w(r201f) &
         + w(r202f) &
         + 0.25_pr * w(r203f) &
         + 0.5_pr * w(r208f) &
         + 0.5_pr * w(r209f) &
         + w(r212f) &
         + w(r213f) &
         + w(r214f) 

    cdot(sIC4H7) = 0.0_pr &
         + w(r148f) &
         - w(r148b) &
         + w(r149f) &
         - w(r149b) &
         + w(r150f) &
         - w(r150b) &
         + w(r151f) &
         - w(r151b) &
         + w(r152f) &
         - w(r152b) &
         - w(r153f) &
         + w(r153b) &
         - w(r173f) &
         - w(r174f) &
         - w(r175f) &
         - w(r179f) &
         + w(r179b) &
         + w(r188f) &
         + w(r198f) &
         - w(r198b) &
         + 0.5_pr * w(r208f) &
         + 0.5_pr * w(r209f) &
         - w(r210f) 

    cdot(sIC8H18) = 0.0_pr &
         - w(r197f) &
         - w(r206f) &
         - w(r207f) 

    cdot(sIC8H16) = 0.0_pr &
         - w(r198f) &
         + w(r198b) &
         - w(r199f) &
         + 0.75_pr * w(r203f) &
         + w(r204f) &
         + w(r205f) &
         - w(r208f) &
         - w(r209f) &
         + w(r214f) &
         + w(r215f) &
         + 2.0_pr * w(r216f) 

    cdot(sIC16H34) = 0.0_pr &
         - w(r213f) &
         - w(r219f) &
         - w(r220f) 

    cdot(sIC12H26) = 0.0_pr &
         - w(r211f) &
         - w(r217f) &
         - w(r218f) 

  
    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) :: CH2_ct
    real(pr) :: CH2_num
    real(pr) :: CH2_denom
    real(pr) :: CH2_CH2GSG
    real(pr) :: CH2_CH
    real(pr) :: CH2_HCO
    real(pr) :: CH2GSG_ct
    real(pr) :: CH2GSG_num
    real(pr) :: CH2GSG_denom
    real(pr) :: CH2GSG_CH2
    real(pr) :: CH2GSG_HCCO
    real(pr) :: CH2GSG_C3H5XA
    real(pr) :: C_ct
    real(pr) :: C_num
    real(pr) :: C_denom
    real(pr) :: C_CH
    real(pr) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_CH2
    real(pr) :: CH_C
    real(pr) :: CH_HCO
    real(pr) :: CH3O_ct
    real(pr) :: CH3O_num
    real(pr) :: CH3O_denom
    real(pr) :: CH2OH_ct
    real(pr) :: CH2OH_num
    real(pr) :: CH2OH_denom
    real(pr) :: CH2OH_CH2CO
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_CH2
    real(pr) :: HCO_CH
    real(pr) :: HCO_CH2CO
    real(pr) :: HCO_C3H3
    real(pr) :: HCO_IC4H7O
    real(pr) :: CH3CO_ct
    real(pr) :: CH3CO_num
    real(pr) :: CH3CO_denom
    real(pr) :: CH3CO_CH2CHO
    real(pr) :: CH2CHO_ct
    real(pr) :: CH2CHO_num
    real(pr) :: CH2CHO_denom
    real(pr) :: CH2CHO_CH3CO
    real(pr) :: CH2CO_ct
    real(pr) :: CH2CO_num
    real(pr) :: CH2CO_denom
    real(pr) :: CH2CO_CH2OH
    real(pr) :: CH2CO_CH2CHO
    real(pr) :: CH2CO_CH3COCH2
    real(pr) :: HCCO_ct
    real(pr) :: HCCO_num
    real(pr) :: HCCO_denom
    real(pr) :: HCCO_CH2GSG
    real(pr) :: HCCO_CH2CO
    real(pr) :: C3H5XA_ct
    real(pr) :: C3H5XA_num
    real(pr) :: C3H5XA_denom
    real(pr) :: C3H5XA_CH2GSG
    real(pr) :: C3H5XA_C3H5O
    real(pr) :: C3H5O_ct
    real(pr) :: C3H5O_num
    real(pr) :: C3H5O_denom
    real(pr) :: C3H5O_C3H5XA
    real(pr) :: C3H3_ct
    real(pr) :: C3H3_num
    real(pr) :: C3H3_denom
    real(pr) :: C3H3_HCO
    real(pr) :: CH3COCH2_ct
    real(pr) :: CH3COCH2_num
    real(pr) :: CH3COCH2_denom
    real(pr) :: CH3COCH2_CH2CO
    real(pr) :: IC4H9_ct
    real(pr) :: IC4H9_num
    real(pr) :: IC4H9_denom
    real(pr) :: IC4H9_TC4H9
    real(pr) :: IC4H9_IC8H17
    real(pr) :: IC4H9_C4H9O2GL1G
    real(pr) :: TC4H9_ct
    real(pr) :: TC4H9_num
    real(pr) :: TC4H9_denom
    real(pr) :: TC4H9_IC4H9
    real(pr) :: TC4H9_IC8H17
    real(pr) :: TC4H9_C4H9O2GL1G
    real(pr) :: IC4H7O_ct
    real(pr) :: IC4H7O_num
    real(pr) :: IC4H7O_denom
    real(pr) :: IC4H7O_HCO
    real(pr) :: NEOC5H11_ct
    real(pr) :: NEOC5H11_num
    real(pr) :: NEOC5H11_denom
    real(pr) :: NEOC5H11_IC8H17
    real(pr) :: IC8H17_ct
    real(pr) :: IC8H17_num
    real(pr) :: IC8H17_denom
    real(pr) :: IC8H17_IC16H33
    real(pr) :: IC8H17_IC12H25
    real(pr) :: IC16H33_ct
    real(pr) :: IC16H33_num
    real(pr) :: IC16H33_denom
    real(pr) :: IC12H25_ct
    real(pr) :: IC12H25_num
    real(pr) :: IC12H25_denom
    real(pr) :: C4H9O2GL1G_ct
    real(pr) :: C4H9O2GL1G_num
    real(pr) :: C4H9O2GL1G_denom
    real(pr) :: C4H9O2GL1G_IC4H9
    real(pr) :: C4H9O2GL1G_TC4H9
    real(pr) :: C3H7GL1G_ct
    real(pr) :: C3H7GL1G_num
    real(pr) :: C3H7GL1G_denom
  
    CH3O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r42f)* c(sO2) &
             + k(r53f)&
             + k(r33b)* c(sO) &
             + k(r39b)* c(sOH) &
             + k(r76b)* c(sCO) &
             + k(r77b)* c(sCO)   )

    CH3O_num = ( 0.0_pr &
             + k(r33f)* c(sCH3) * c(sO2) &
             + k(r39f)* c(sCH3) * c(sHO2) &
             + k(r76f)* c(sC2H3) * c(sO2) &
             + k(r77f)* c(sC2H3) * c(sO2) &
             + k(r172f)* c(sC2H4) * c(sO2) &
             + k(r42b)* c(sCH2O) * c(sHO2) &
             + k(r53b)* c(sCH2O) * c(sH)   )

    CH3O_ct = CH3O_num / CH3O_denom

    cqss(sqssCH3O) = CH3O_ct  

    CH3CO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r85f)&
             + k(r127b)* c(sCH2O) &
             + k(r176b)  )

    CH3CO_num = ( 0.0_pr &
             + k(r127f)* c(sC3H5XT) * c(sO2) &
             + k(r85b)* c(sCH3) * c(sCO)   )

    CH3CO_ct = CH3CO_num / CH3CO_denom

    CH3CO_CH2CHO = - ( 0.0_pr &
          + k(r176f)  ) / CH3CO_denom

    CH2CHO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r86f)&
             + k(r164f)* c(sO2) &
             + k(r176f)&
             + k(r64b)* c(sH) &
             + k(r66b)* c(sO) &
             + k(r67b)* c(sO)   )

    CH2CHO_num = ( 0.0_pr &
             + k(r64f)* c(sC2H4) * c(sO) &
             + k(r66f)* c(sC2H3) * c(sO2) &
             + k(r67f)* c(sC2H3) * c(sO2) &
             + k(r86b)* c(sCH3) * c(sCO)   )

    CH2CHO_ct = CH2CHO_num / CH2CHO_denom

    CH2CHO_CH3CO = - ( 0.0_pr &
          + k(r176b)  ) / CH2CHO_denom

    cqss(sqssCH3CO) = (  (CH3CO_ct) - (CH2CHO_ct) * (CH3CO_CH2CHO) &
                 ) / (  (1.0_pr) - (CH2CHO_CH3CO) * (CH3CO_CH2CHO) &
                 )  

    cqss(sqssCH2CHO) = CH2CHO_ct - (CH2CHO_CH3CO) * cqss(sqssCH3CO)  

    IC16H33_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r214f)&
             + k(r215f)&
             + k(r216f)* c(sO2)   )

    IC16H33_num = ( 0.0_pr &
             + k(r219f)* c(sH) * c(sIC16H34) &
             + k(r220f)* c(sIC16H34) * c(sOH)   )

    IC16H33_ct = IC16H33_num / IC16H33_denom

    cqss(sqssIC16H33) = IC16H33_ct  

    IC12H25_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r212f)  )

    IC12H25_num = ( 0.0_pr &
             + k(r217f)* c(sH) * c(sIC12H26) &
             + k(r218f)* c(sIC12H26) * c(sOH)   )

    IC12H25_ct = IC12H25_num / IC12H25_denom

    cqss(sqssIC12H25) = IC12H25_ct  

    C3H7GL1G_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r118b)&
             + k(r119b)&
             + k(r122b)&
             + k(r123b)&
             + k(r124b)&
             + k(r125b)  )

    C3H7GL1G_num = ( 0.0_pr &
             + k(r118f)* c(sC3H6) * c(sH) &
             + k(r119f)* c(sC3H6) * c(sH) &
             + k(r122f)* c(sC3H6) * c(sH) &
             + k(r123f)* c(sC3H6) * c(sH) &
             + k(r124f)* c(sC2H4) * c(sCH3) &
             + k(r125f)* c(sC2H4) * c(sCH3) &
     +0.2_pr * k(r197f)* c(sIC8H18)   )

    C3H7GL1G_ct = C3H7GL1G_num / C3H7GL1G_denom

    cqss(sqssC3H7GL1G) = C3H7GL1G_ct  

    CH2CO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r87f)* c(sH) &
             + k(r88f)* c(sOH) &
             + k(r163f)* c(sO2) &
             + k(r185f)* c(sH) &
             + k(r186f)* c(sO) &
             + k(r187f)* c(sOH) &
             + k(r70b)* c(sOH) &
             + k(r71b)* c(sOH) &
             + k(r83b)* c(sH) &
             + k(r166b)* c(sCH3)   )

    CH2CO_num = ( 0.0_pr &
             + k(r70f)* c(sC2H3) * c(sO2) &
             + k(r71f)* c(sC2H3) * c(sO2) &
             + k(r83f)* c(sC2H2) * c(sOH) &
             + k(r162f)* c(sC3H4XA) * c(sO2) &
             + k(r164f)* cqss(sqssCH2CHO) * c(sO2) &
             + k(r87b)* c(sCH3) * c(sCO)   )

    CH2CO_ct = CH2CO_num / CH2CO_denom

    CH2CO_CH2OH = - ( 0.0_pr &
          + k(r88b) * c(sCO)  ) / CH2CO_denom

    CH2CO_CH2CHO = - ( 0.0_pr &
          + k(r164f) * cqss(sqssCH2CHO) * c(sO2)  ) / CH2CO_denom

    CH2CO_CH3COCH2 = - ( 0.0_pr &
          + k(r166f)  ) / CH2CO_denom

    CH2OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r40f)* c(sO2) &
             + k(r41f)* c(sO2) &
             + k(r37b)* c(sH) &
             + k(r52b)&
             + k(r88b)* c(sCO)   )

    CH2OH_num = ( 0.0_pr &
             + k(r37f)* c(sCH3) * c(sOH) &
             + k(r52f)* c(sCH2O) * c(sH) &
             + k(r40b)* c(sCH2O) * c(sHO2) &
             + k(r41b)* c(sCH2O) * c(sHO2)   )

    CH2OH_ct = CH2OH_num / CH2OH_denom

    CH2OH_CH2CO = - ( 0.0_pr &
          + k(r88f) * c(sOH)  ) / CH2OH_denom

    CH3COCH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r166f)&
             + k(r126b)* c(sO)   )

    CH3COCH2_num = ( 0.0_pr &
             + k(r126f)* c(sC3H5XT) * c(sO2) &
             + k(r210f)* c(sIC4H7) * c(sO2)   )

    CH3COCH2_ct = CH3COCH2_num / CH3COCH2_denom

    CH3COCH2_CH2CO = - ( 0.0_pr &
          + k(r166b) * c(sCH3)  ) / CH3COCH2_denom

    cqss(sqssCH2CO) = (  ( (CH2CO_ct) - (CH3COCH2_ct) * (CH2CO_CH3COCH2)) &
                 - (CH2OH_ct) * (CH2CO_CH2OH) ) / (  ( (1.0_pr) &
                 - (CH3COCH2_CH2CO) * (CH2CO_CH3COCH2)) - (CH2OH_CH2CO) &
                 * (CH2CO_CH2OH) )  

    cqss(sqssCH2OH) = CH2OH_ct - (CH2OH_CH2CO) * cqss(sqssCH2CO)  

    cqss(sqssCH3COCH2) = CH3COCH2_ct - (CH3COCH2_CH2CO) * cqss(sqssCH2CO)  

    IC8H17_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r200f)&
             + k(r201f)&
             + k(r202f)&
             + k(r203f)&
             + k(r204f)&
             + k(r205f)* c(sO2)   )

    IC8H17_num = ( 0.0_pr &
             + k(r199f)* c(sH) * c(sIC8H16) &
             + k(r206f)* c(sIC8H18) * c(sOH) &
             + k(r207f)* c(sH) * c(sIC8H18) &
             + k(r211f)* c(sIC12H26) &
             + k(r212f)* cqss(sqssIC12H25) &
     +0.33_pr * k(r213f)* c(sIC16H34) &
             + k(r215f)* cqss(sqssIC16H33)   )

    IC8H17_ct = IC8H17_num / IC8H17_denom

    cqss(sqssIC8H17) = IC8H17_ct  

    CH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r27f)* c(sO2) &
             + k(r28f)* c(sO2) &
             + k(r29f)* c(sH) &
             + k(r167f)* c(sO) &
             + k(r168f)* c(sOH) &
             + k(r23b)* c(sN2) &
             + k(r38b)* c(sH2O)   )

    CH2_num = ( 0.0_pr &
             + k(r38f)* c(sCH3) * c(sOH) &
             + k(r167b)* c(sCO) * c(sH2) &
             + k(r168b)* c(sCH2O) * c(sH)   )

    CH2_ct = CH2_num / CH2_denom

    CH2_CH2GSG = - ( 0.0_pr &
          + k(r23f) * c(sN2)  ) / CH2_denom

    CH2_CH = - ( 0.0_pr &
          + k(r29b) * c(sH2)  ) / CH2_denom

    CH2_HCO = - ( 0.0_pr &
          + k(r27b) * c(sOH)  ) / CH2_denom

    C_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r32f)* c(sO2) &
             + k(r31b)* c(sH2)   )

    C_num = ( 0.0_pr &
             + k(r32b)* c(sCO) * c(sO)   )

    C_ct = C_num / C_denom

    C_CH = - ( 0.0_pr &
          + k(r31f) * c(sH)  ) / C_denom

    C3H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r144f)* c(sO) &
             + k(r145f)* c(sOH) &
             + k(r169f)* c(sH) &
             + k(r131b)* c(sH2O)   )

    C3H3_num = ( 0.0_pr &
             + k(r131f)* c(sC3H4XA) * c(sOH) &
             + k(r184f)* c(sC3H4XA) * c(sO) &
             + k(r144b)* c(sC2H) * c(sCH2O) &
             + k(r169b)* c(sC3H4XP)   )

    C3H3_ct = C3H3_num / C3H3_denom

    C3H3_HCO = - ( 0.0_pr &
          + k(r145b) * c(sC2H3)  ) / C3H3_denom

    C3H5O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r129f)&
             + k(r130f)&
             + k(r128b)* c(sOH)   )

    C3H5O_num = ( 0.0_pr &
             + k(r129b)* c(sC2H3) * c(sCH2O) &
             + k(r130b)* c(sC2H3CHO) * c(sH)   )

    C3H5O_ct = C3H5O_num / C3H5O_denom

    C3H5O_C3H5XA = - ( 0.0_pr &
          + k(r128f) * c(sHO2)  ) / C3H5O_denom

    HCCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r89f)* c(sOH) &
             + k(r90f)* c(sH) &
             + k(r91f)* c(sO2) &
             + k(r82b)* c(sH)   )

    HCCO_num = ( 0.0_pr &
             + k(r82f)* c(sC2H2) * c(sO) &
             + k(r185f)* cqss(sqssCH2CO) * c(sH) &
             + k(r186f)* cqss(sqssCH2CO) * c(sO) &
             + k(r187f)* cqss(sqssCH2CO) * c(sOH)   )

    HCCO_ct = HCCO_num / HCCO_denom

    HCCO_CH2GSG = - ( 0.0_pr &
          + k(r90b) * c(sCO)  ) / HCCO_denom

    HCCO_CH2CO = - ( 0.0_pr &
          + k(r185f) * cqss(sqssCH2CO) * c(sH) &
          + k(r186f) * cqss(sqssCH2CO) * c(sO) &
          + k(r187f) * cqss(sqssCH2CO) * c(sOH)  ) / HCCO_denom

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r47f)* m(mM47) &
             + k(r48f)* c(sO2) &
             + k(r49f)* c(sH) &
             + k(r50f)* c(sOH) &
             + k(r51f)* c(sHO2) &
             + k(r27b)* c(sOH) &
             + k(r30b)* c(sO) &
             + k(r43b)* c(sHO2) &
             + k(r44b)* c(sH2) &
             + k(r45b)* c(sH2O) &
             + k(r46b)* c(sH2O2) &
             + k(r72b)* c(sCH2O) &
             + k(r73b)* c(sCH2O) &
             + k(r84b)* c(sCO) &
             + k(r117b)* c(sC2H5) &
             + k(r145b)* c(sC2H3) &
             + k(r155b)* c(sC3H6)   )

    HCO_num = ( 0.0_pr &
             + k(r43f)* c(sCH2O) * c(sO2) &
             + k(r44f)* c(sCH2O) * c(sH) &
             + k(r45f)* c(sCH2O) * c(sOH) &
             + k(r46f)* c(sCH2O) * c(sHO2) &
             + k(r72f)* c(sC2H3) * c(sO2) &
             + k(r73f)* c(sC2H3) * c(sO2) &
             + k(r84f)* c(sC2H) * c(sO2) &
             + k(r117f)* c(sC3H6) * c(sO) &
             + k(r163f)* cqss(sqssCH2CO) * c(sO2) &
             + k(r172f)* c(sC2H4) * c(sO2) &
             + k(r47b)* c(sCO) * c(sH) * m(mM47) &
             + k(r48b)* c(sCO) * c(sHO2) &
             + k(r49b)* c(sCO) * c(sH2) &
             + k(r50b)* c(sCO) * c(sH2O)   )

    HCO_ct = HCO_num / HCO_denom

    HCO_CH2 = - ( 0.0_pr &
          + k(r27f) * c(sO2)  ) / HCO_denom

    HCO_CH = - ( 0.0_pr &
          + k(r30f) * c(sO2)  ) / HCO_denom

    HCO_CH2CO = - ( 0.0_pr &
          + k(r163f) * cqss(sqssCH2CO) * c(sO2)  ) / HCO_denom

    HCO_C3H3 = - ( 0.0_pr &
          + k(r145f) * c(sOH)  ) / HCO_denom

    HCO_IC4H7O = - ( 0.0_pr &
          + k(r155f)  ) / HCO_denom

    C3H5XA_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r128f)* c(sHO2) &
             + k(r178f)* c(sH) &
             + k(r179f)* c(sIC4H7) &
             + k(r95b)* c(sH) &
             + k(r96b)* c(sH) &
             + k(r99b)* c(sH) &
             + k(r100b)* c(sH) &
             + k(r103b)* c(sH) &
             + k(r104b)* c(sH) &
             + k(r107b)* c(sH) &
             + k(r108b)* c(sH) &
             + k(r111b)* c(sH2) &
             + k(r112b)* c(sHO2) &
             + k(r113b)* c(sOH) &
             + k(r114b)* c(sH2O) &
             + k(r134b)&
             + k(r135b)  )

    C3H5XA_num = ( 0.0_pr &
             + k(r99f)* c(sC2H3) * c(sCH3) &
             + k(r100f)* c(sC2H3) * c(sCH3) &
             + k(r103f)* c(sC3H6) &
             + k(r104f)* c(sC3H6) &
             + k(r107f)* c(sC3H5XT) * c(sH) &
             + k(r108f)* c(sC3H5XT) * c(sH) &
             + k(r111f)* c(sC3H6) * c(sH) &
             + k(r112f)* c(sC3H6) * c(sO2) &
             + k(r113f)* c(sC3H6) * c(sO) &
             + k(r114f)* c(sC3H6) * c(sOH) &
             + k(r134f)* c(sC3H4XA) * c(sH) &
             + k(r135f)* c(sC3H4XA) * c(sH) &
             + k(r178b)* c(sC3H6) &
             + k(r179b)* c(sC3H4XA) * c(sIC4H8)   )

    C3H5XA_ct = C3H5XA_num / C3H5XA_denom

    C3H5XA_CH2GSG = - ( 0.0_pr &
          + k(r95f) * c(sC2H4) &
          + k(r96f) * c(sC2H4)  ) / C3H5XA_denom

    C3H5XA_C3H5O = - ( 0.0_pr &
          + k(r128b) * c(sOH)  ) / C3H5XA_denom

    CH2GSG_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r23f)* c(sN2) &
             + k(r24f)* c(sO2) &
             + k(r25f)* c(sO2) &
             + k(r26f)* c(sH2) &
             + k(r93f)* c(sC2H4) &
             + k(r94f)* c(sC2H4) &
             + k(r95f)* c(sC2H4) &
             + k(r96f)* c(sC2H4) &
             + k(r97f)* c(sC2H4) &
             + k(r98f)* c(sC2H4) &
             + k(r36b)* c(sH2O) &
             + k(r90b)* c(sCO)   )

    CH2GSG_num = ( 0.0_pr &
             + k(r36f)* c(sCH3) * c(sOH) &
             + k(r25b)* c(sCO) * c(sH2O) &
             + k(r26b)* c(sCH3) * c(sH) &
             + k(r93b)* c(sC3H6) &
             + k(r94b)* c(sC3H6) &
             + k(r97b)* c(sC2H3) * c(sCH3) &
             + k(r98b)* c(sC2H3) * c(sCH3)   )

    CH2GSG_ct = CH2GSG_num / CH2GSG_denom

    CH2GSG_CH2 = - ( 0.0_pr &
          + k(r23b) * c(sN2)  ) / CH2GSG_denom

    CH2GSG_HCCO = - ( 0.0_pr &
          + k(r90f) * c(sH)  ) / CH2GSG_denom

    CH2GSG_C3H5XA = - ( 0.0_pr &
          + k(r95b) * c(sH) &
          + k(r96b) * c(sH)  ) / CH2GSG_denom

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r30f)* c(sO2) &
             + k(r31f)* c(sH) &
             + k(r29b)* c(sH2)   )

    CH_num = ( 0.0_pr   )

    CH_ct = CH_num / CH_denom

    CH_CH2 = - ( 0.0_pr &
          + k(r29f) * c(sH)  ) / CH_denom

    CH_C = - ( 0.0_pr &
          + k(r31b) * c(sH2)  ) / CH_denom

    CH_HCO = - ( 0.0_pr &
          + k(r30b) * c(sO)  ) / CH_denom

    IC4H7O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r154f)&
             + k(r155f)&
             + k(r153b)* c(sOH)   )

    IC4H7O_num = ( 0.0_pr &
             + k(r153f)* c(sHO2) * c(sIC4H7) &
             + k(r154b)* c(sC3H5XT) * c(sCH2O)   )

    IC4H7O_ct = IC4H7O_num / IC4H7O_denom

    IC4H7O_HCO = - ( 0.0_pr &
          + k(r155b) * c(sC3H6)  ) / IC4H7O_denom

    cqss(sqssCH2) = (  ( ( ( ( ( ( ( (CH2_ct) - (CH_ct) * (CH2_CH)) &
                 - (CH2GSG_ct) * (CH2_CH2GSG)) - ( (C3H5XA_ct) &
                 - (CH2GSG_ct) * (C3H5XA_CH2GSG)) * ( - (CH2GSG_C3H5XA) &
                 * (CH2_CH2GSG)) / ( (1.0_pr) - (CH2GSG_C3H5XA) &
                 * (C3H5XA_CH2GSG))) - ( ( (HCO_ct) - (IC4H7O_ct) &
                 * (HCO_IC4H7O)) - (CH_ct) * (HCO_CH)) * ( (CH2_HCO) &
                 - (CH_HCO) * (CH2_CH)) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) - ( ( (HCCO_ct) &
                 - (CH2GSG_ct) * (HCCO_CH2GSG)) - ( (C3H5XA_ct) &
                 - (CH2GSG_ct) * (C3H5XA_CH2GSG)) * ( - (CH2GSG_C3H5XA) &
                 * (HCCO_CH2GSG)) / ( (1.0_pr) - (CH2GSG_C3H5XA) &
                 * (C3H5XA_CH2GSG))) * ( - (CH2GSG_HCCO) * (CH2_CH2GSG)) &
                 / ( (1.0_pr) - (CH2GSG_HCCO) * (HCCO_CH2GSG))) &
                 - ( (C3H5O_ct) - ( (C3H5XA_ct) - (CH2GSG_ct) &
                 * (C3H5XA_CH2GSG)) * (C3H5O_C3H5XA) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG))) * ( - (C3H5XA_C3H5O) &
                 * ( - (CH2GSG_C3H5XA) * (CH2_CH2GSG)) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG))) / ( (1.0_pr) &
                 - (C3H5XA_C3H5O) * (C3H5O_C3H5XA) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG)))) - ( (C3H3_ct) &
                 - ( ( (HCO_ct) - (IC4H7O_ct) * (HCO_IC4H7O)) &
                 - (CH_ct) * (HCO_CH)) * (C3H3_HCO) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) &
                 * ( - (HCO_C3H3) * ( (CH2_HCO) - (CH_HCO) * (CH2_CH)) &
                 / ( ( (1.0_pr) - (IC4H7O_HCO) * (HCO_IC4H7O)) &
                 - (CH_HCO) * (HCO_CH))) / ( (1.0_pr) - (HCO_C3H3) &
                 * (C3H3_HCO) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH)))) - ( ( ( (C_ct) &
                 - (CH_ct) * (C_CH)) - ( ( (HCO_ct) - (IC4H7O_ct) &
                 * (HCO_IC4H7O)) - (CH_ct) * (HCO_CH)) * ( - (CH_HCO) &
                 * (C_CH)) / ( ( (1.0_pr) - (IC4H7O_HCO) * (HCO_IC4H7O)) &
                 - (CH_HCO) * (HCO_CH))) - ( (C3H3_ct) - ( ( (HCO_ct) &
                 - (IC4H7O_ct) * (HCO_IC4H7O)) - (CH_ct) * (HCO_CH)) &
                 * (C3H3_HCO) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) * ( - (HCO_C3H3) &
                 * ( - (CH_HCO) * (C_CH)) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) / ( (1.0_pr) &
                 - (HCO_C3H3) * (C3H3_HCO) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH)))) * ( - (CH_C) &
                 * (CH2_CH)) / ( (1.0_pr) - (CH_C) * (C_CH)) &
                 ) / (  ( ( (1.0_pr) - (CH_CH2) * (CH2_CH)) &
                 - (CH2GSG_CH2) * (CH2_CH2GSG)) - ( (HCO_CH2) &
                 - (CH_CH2) * (HCO_CH)) * ( (CH2_HCO) - (CH_HCO) &
                 * (CH2_CH)) / ( ( (1.0_pr) - (IC4H7O_HCO) * (HCO_IC4H7O)) &
                 - (CH_HCO) * (HCO_CH)) )  

    cqss(sqssC) = ( ( ( (C_ct) - (CH_ct) * (C_CH)) - ( ( (HCO_ct) &
                 - (IC4H7O_ct) * (HCO_IC4H7O)) - (CH_ct) * (HCO_CH)) &
                 * ( - (CH_HCO) * (C_CH)) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) - ( (C3H3_ct) &
                 - ( ( (HCO_ct) - (IC4H7O_ct) * (HCO_IC4H7O)) &
                 - (CH_ct) * (HCO_CH)) * (C3H3_HCO) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) &
                 * ( - (HCO_C3H3) * ( - (CH_HCO) * (C_CH)) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) &
                 / ( (1.0_pr) - (HCO_C3H3) * (C3H3_HCO) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) &
                 - ( ( - (CH_CH2) * (C_CH)) - ( (HCO_CH2) - (CH_CH2) &
                 * (HCO_CH)) * ( - (CH_HCO) * (C_CH)) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH))) &
                 * cqss(sqssCH2)) / ( (1.0_pr) - (CH_C) * (C_CH))  

    cqss(sqssC3H3) = ( (C3H3_ct) - ( ( (HCO_ct) - (IC4H7O_ct) * (HCO_IC4H7O)) &
                 - (CH_ct) * (HCO_CH)) * (C3H3_HCO) / ( ( (1.0_pr) &
                 - (IC4H7O_HCO) * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH)) &
                 - ( - ( (HCO_CH2) - (CH_CH2) * (HCO_CH)) * (C3H3_HCO) &
                 / ( ( (1.0_pr) - (IC4H7O_HCO) * (HCO_IC4H7O)) &
                 - (CH_HCO) * (HCO_CH))) * cqss(sqssCH2)) / ( (1.0_pr) &
                 - (HCO_C3H3) * (C3H3_HCO) / ( ( (1.0_pr) - (IC4H7O_HCO) &
                 * (HCO_IC4H7O)) - (CH_HCO) * (HCO_CH)))  

    cqss(sqssC3H5O) =  (C3H5O_ct) - ( (C3H5XA_ct) - (CH2GSG_ct) &
                 * (C3H5XA_CH2GSG)) * (C3H5O_C3H5XA) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG)) /  (1.0_pr) &
                 - (C3H5XA_C3H5O) * (C3H5O_C3H5XA) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG))  

    cqss(sqssHCCO) = ( ( (HCCO_ct) - (CH2GSG_ct) * (HCCO_CH2GSG)) &
                 - ( (C3H5XA_ct) - (CH2GSG_ct) * (C3H5XA_CH2GSG)) &
                 * ( - (CH2GSG_C3H5XA) * (HCCO_CH2GSG)) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG)) - ( - (CH2GSG_CH2) &
                 * (HCCO_CH2GSG)) * cqss(sqssCH2) - ( - (C3H5XA_C3H5O) &
                 * ( - (CH2GSG_C3H5XA) * (HCCO_CH2GSG)) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG))) * cqss(sqssC3H5O)) &
                 / ( (1.0_pr) - (CH2GSG_HCCO) * (HCCO_CH2GSG))  

    cqss(sqssHCO) = ( ( (HCO_ct) - (IC4H7O_ct) * (HCO_IC4H7O)) &
                 - (CH_ct) * (HCO_CH) - ( (HCO_CH2) - (CH_CH2) &
                 * (HCO_CH)) * cqss(sqssCH2) - ( - (CH_C) * (HCO_CH)) &
                 * cqss(sqssC) - (HCO_C3H3) * cqss(sqssC3H3)) &
                 / ( ( (1.0_pr) - (IC4H7O_HCO) * (HCO_IC4H7O)) &
                 - (CH_HCO) * (HCO_CH))  

    cqss(sqssC3H5XA) = ( (C3H5XA_ct) - (CH2GSG_ct) * (C3H5XA_CH2GSG) &
                 - ( - (CH2GSG_CH2) * (C3H5XA_CH2GSG)) * cqss(sqssCH2) &
                 - (C3H5XA_C3H5O) * cqss(sqssC3H5O) - ( - (CH2GSG_HCCO) &
                 * (C3H5XA_CH2GSG)) * cqss(sqssHCCO)) / ( (1.0_pr) &
                 - (CH2GSG_C3H5XA) * (C3H5XA_CH2GSG))  

    cqss(sqssCH2GSG) = CH2GSG_ct - (CH2GSG_CH2) * cqss(sqssCH2) - (CH2GSG_HCCO) &
                 * cqss(sqssHCCO) - (CH2GSG_C3H5XA) * cqss(sqssC3H5XA)  

    cqss(sqssCH) = CH_ct - (CH_CH2) * cqss(sqssCH2) - (CH_C) &
                 * cqss(sqssC) - (CH_HCO) * cqss(sqssHCO)  

    cqss(sqssIC4H7O) = IC4H7O_ct - (IC4H7O_HCO) * cqss(sqssHCO)  

    IC4H9_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r146f)&
             + k(r147f)&
             + k(r189f)* c(sO2) &
             + k(r156b)&
             + k(r157b)  )

    IC4H9_num = ( 0.0_pr &
             + k(r156f)* c(sH) * c(sIC4H8) &
             + k(r157f)* c(sH) * c(sIC4H8) &
     +0.8_pr * k(r197f)* c(sIC8H18) &
             + k(r202f)* cqss(sqssIC8H17) &
             + k(r147b)* c(sC3H6) * c(sCH3)   )

    IC4H9_ct = IC4H9_num / IC4H9_denom

    IC4H9_TC4H9 = - ( 0.0_pr &
          + k(r146b)  ) / IC4H9_denom

    IC4H9_C4H9O2GL1G = - ( 0.0_pr &
          + k(r191f)  ) / IC4H9_denom

    TC4H9_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r181f)* c(sO2) &
             + k(r190f)* c(sO2) &
             + k(r146b)&
             + k(r158b)&
             + k(r159b)&
             + k(r198b)* c(sIC4H7)   )

    TC4H9_num = ( 0.0_pr &
             + k(r158f)* c(sH) * c(sIC4H8) &
             + k(r159f)* c(sH) * c(sIC4H8) &
     +0.8_pr * k(r197f)* c(sIC8H18) &
             + k(r198f)* c(sIC8H16) &
             + k(r201f)* cqss(sqssIC8H17) &
             + k(r211f)* c(sIC12H26) &
             + k(r213f)* c(sIC16H34) &
             + k(r181b)* c(sHO2) * c(sIC4H8)   )

    TC4H9_ct = TC4H9_num / TC4H9_denom

    TC4H9_IC4H9 = - ( 0.0_pr &
          + k(r146f)  ) / TC4H9_denom

    TC4H9_C4H9O2GL1G = - ( 0.0_pr &
          + k(r192f)  ) / TC4H9_denom

    C4H9O2GL1G_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r191f)&
             + k(r192f)&
             + k(r193f)&
             + k(r194f)&
             + k(r195f)  )

    C4H9O2GL1G_num = ( 0.0_pr   )

    C4H9O2GL1G_ct = C4H9O2GL1G_num / C4H9O2GL1G_denom

    C4H9O2GL1G_IC4H9 = - ( 0.0_pr &
          + k(r189f) * c(sO2)  ) / C4H9O2GL1G_denom

    C4H9O2GL1G_TC4H9 = - ( 0.0_pr &
          + k(r190f) * c(sO2)  ) / C4H9O2GL1G_denom

    cqss(sqssIC4H9) = (  ( (IC4H9_ct) - (C4H9O2GL1G_ct) * (IC4H9_C4H9O2GL1G)) &
                 - ( (TC4H9_ct) - (C4H9O2GL1G_ct) * (TC4H9_C4H9O2GL1G)) &
                 * ( (IC4H9_TC4H9) - (C4H9O2GL1G_TC4H9) * (IC4H9_C4H9O2GL1G)) &
                 / ( (1.0_pr) - (C4H9O2GL1G_TC4H9) * (TC4H9_C4H9O2GL1G)) &
                 ) / (  ( (1.0_pr) - (C4H9O2GL1G_IC4H9) * (IC4H9_C4H9O2GL1G)) &
                 - ( (TC4H9_IC4H9) - (C4H9O2GL1G_IC4H9) * (TC4H9_C4H9O2GL1G)) &
                 * ( (IC4H9_TC4H9) - (C4H9O2GL1G_TC4H9) * (IC4H9_C4H9O2GL1G)) &
                 / ( (1.0_pr) - (C4H9O2GL1G_TC4H9) * (TC4H9_C4H9O2GL1G)) &
                 )  

    cqss(sqssTC4H9) = ( (TC4H9_ct) - (C4H9O2GL1G_ct) * (TC4H9_C4H9O2GL1G) &
                 - ( (TC4H9_IC4H9) - (C4H9O2GL1G_IC4H9) * (TC4H9_C4H9O2GL1G)) &
                 * cqss(sqssIC4H9)) / ( (1.0_pr) - (C4H9O2GL1G_TC4H9) &
                 * (TC4H9_C4H9O2GL1G))  

    cqss(sqssC4H9O2GL1G) = C4H9O2GL1G_ct - (C4H9O2GL1G_IC4H9) * cqss(sqssIC4H9) &
                 - (C4H9O2GL1G_TC4H9) * cqss(sqssTC4H9)  

    NEOC5H11_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r196f)  )

    NEOC5H11_num = ( 0.0_pr &
     +0.2_pr * k(r197f)* c(sIC8H18) &
             + k(r200f)* cqss(sqssIC8H17) &
     +0.5_pr * k(r208f)* c(sH) * c(sIC8H16) &
     +0.5_pr * k(r209f)* c(sIC8H16) * c(sOH) &
     +0.67_pr * k(r213f)* c(sIC16H34) &
     +0.33_pr * k(r214f)* cqss(sqssIC16H33) &
             + k(r196b)* c(sCH3) * c(sIC4H8)   )

    NEOC5H11_ct = NEOC5H11_num / NEOC5H11_denom

    cqss(sqssNEOC5H11) = NEOC5H11_ct  

    do index = 1, nqss
      cqss(index) = max(cqss(index), 1e-60_pr)
      cqss(index) = min(cqss(index), 1e60_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