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

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

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

  ! Ideal gas constant
  real(pr), parameter :: Rcst = 8.3144621_pr
  
  ! Use semi-implicit formulation for select species
  logical, parameter :: semiImplicit = .False.
  
  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 20
  integer, parameter :: nqss = 9
  integer, parameter :: nreac = 166
  integer, parameter :: nreac_reverse = 158
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 5
  integer, parameter :: nFO = 18
  integer, parameter :: nTB_reverse = 5
  integer, parameter :: nFO_reverse = 18
  
  ! Index of species
  integer, parameter :: sH2 = 1
  integer, parameter :: sH = 2
  integer, parameter :: sO = 3
  integer, parameter :: sO2 = 4
  integer, parameter :: sOH = 5
  integer, parameter :: sH2O = 6
  integer, parameter :: sHO2 = 7
  integer, parameter :: sH2O2 = 8
  integer, parameter :: sCH3 = 9
  integer, parameter :: sCH4 = 10
  integer, parameter :: sCO = 11
  integer, parameter :: sCO2 = 12
  integer, parameter :: sCH2O = 13
  integer, parameter :: sCH3O = 14
  integer, parameter :: sCH3OH = 15
  integer, parameter :: sC2H2 = 16
  integer, parameter :: sC2H4 = 17
  integer, parameter :: sC2H6 = 18
  integer, parameter :: sCH2CO = 19
  integer, parameter :: sN2 = 20

  integer, parameter :: sqssCH = 1
  integer, parameter :: sqssCH2 = 2
  integer, parameter :: sqssCH2GSG = 3
  integer, parameter :: sqssHCO = 4
  integer, parameter :: sqssCH2OH = 5
  integer, parameter :: sqssC2H3 = 6
  integer, parameter :: sqssC2H5 = 7
  integer, parameter :: sqssHCCO = 8
  integer, parameter :: sqssCH2CHO = 9

  ! 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 :: r1b = 167
  integer, parameter :: r2b = 168
  integer, parameter :: r3b = 169
  integer, parameter :: r4b = 170
  integer, parameter :: r5b = 171
  integer, parameter :: r6b = 172
  integer, parameter :: r7b = 173
  integer, parameter :: r8b = 174
  integer, parameter :: r9b = 175
  integer, parameter :: r10b = 176
  integer, parameter :: r11b = 177
  integer, parameter :: r12b = 178
  integer, parameter :: r13b = 179
  integer, parameter :: r14b = 180
  integer, parameter :: r15b = 181
  integer, parameter :: r16b = 182
  integer, parameter :: r17b = 183
  integer, parameter :: r18b = 184
  integer, parameter :: r19b = 185
  integer, parameter :: r20b = 186
  integer, parameter :: r21b = 187
  integer, parameter :: r22b = 188
  integer, parameter :: r23b = 189
  integer, parameter :: r24b = 190
  integer, parameter :: r25b = 191
  integer, parameter :: r26b = 192
  integer, parameter :: r27b = 193
  integer, parameter :: r28b = 194
  integer, parameter :: r29b = 195
  integer, parameter :: r30b = 196
  integer, parameter :: r31b = 197
  integer, parameter :: r32b = 198
  integer, parameter :: r33b = 199
  integer, parameter :: r34b = 200
  integer, parameter :: r35b = 201
  integer, parameter :: r36b = 202
  integer, parameter :: r37b = 203
  integer, parameter :: r38b = 204
  integer, parameter :: r39b = 205
  integer, parameter :: r40b = 206
  integer, parameter :: r41b = 207
  integer, parameter :: r42b = 208
  integer, parameter :: r43b = 209
  integer, parameter :: r44b = 210
  integer, parameter :: r45b = 211
  integer, parameter :: r46b = 212
  integer, parameter :: r47b = 213
  integer, parameter :: r48b = 214
  integer, parameter :: r49b = 215
  integer, parameter :: r50b = 216
  integer, parameter :: r51b = 217
  integer, parameter :: r52b = 218
  integer, parameter :: r53b = 219
  integer, parameter :: r54b = 220
  integer, parameter :: r55b = 221
  integer, parameter :: r56b = 222
  integer, parameter :: r57b = 223
  integer, parameter :: r58b = 224
  integer, parameter :: r59b = 225
  integer, parameter :: r60b = 226
  integer, parameter :: r61b = 227
  integer, parameter :: r62b = 228
  integer, parameter :: r63b = 229
  integer, parameter :: r64b = 230
  integer, parameter :: r65b = 231
  integer, parameter :: r66b = 232
  integer, parameter :: r67b = 233
  integer, parameter :: r68b = 234
  integer, parameter :: r69b = 235
  integer, parameter :: r70b = 236
  integer, parameter :: r71b = 237
  integer, parameter :: r72b = 238
  integer, parameter :: r73b = 239
  integer, parameter :: r74b = 240
  integer, parameter :: r75b = 241
  integer, parameter :: r76b = 242
  integer, parameter :: r77b = 243
  integer, parameter :: r78b = 244
  integer, parameter :: r79b = 245
  integer, parameter :: r80b = 246
  integer, parameter :: r81b = 247
  integer, parameter :: r82b = 248
  integer, parameter :: r83b = 249
  integer, parameter :: r84b = 250
  integer, parameter :: r85b = 251
  integer, parameter :: r86b = 252
  integer, parameter :: r87b = 253
  integer, parameter :: r88b = 254
  integer, parameter :: r89b = 255
  integer, parameter :: r90b = 256
  integer, parameter :: r91b = 257
  integer, parameter :: r92b = 258
  integer, parameter :: r93b = 259
  integer, parameter :: r94b = 260
  integer, parameter :: r95b = 261
  integer, parameter :: r96b = 262
  integer, parameter :: r97b = 263
  integer, parameter :: r98b = 264
  integer, parameter :: r99b = 265
  integer, parameter :: r100b = 266
  integer, parameter :: r101b = 267
  integer, parameter :: r102b = 268
  integer, parameter :: r103b = 269
  integer, parameter :: r104b = 270
  integer, parameter :: r105b = 271
  integer, parameter :: r106b = 272
  integer, parameter :: r107b = 273
  integer, parameter :: r108b = 274
  integer, parameter :: r109b = 275
  integer, parameter :: r110b = 276
  integer, parameter :: r111b = 277
  integer, parameter :: r112b = 278
  integer, parameter :: r113b = 279
  integer, parameter :: r114b = 280
  integer, parameter :: r116b = 281
  integer, parameter :: r117b = 282
  integer, parameter :: r118b = 283
  integer, parameter :: r119b = 284
  integer, parameter :: r120b = 285
  integer, parameter :: r121b = 286
  integer, parameter :: r122b = 287
  integer, parameter :: r123b = 288
  integer, parameter :: r124b = 289
  integer, parameter :: r125b = 290
  integer, parameter :: r126b = 291
  integer, parameter :: r127b = 292
  integer, parameter :: r128b = 293
  integer, parameter :: r129b = 294
  integer, parameter :: r130b = 295
  integer, parameter :: r131b = 296
  integer, parameter :: r132b = 297
  integer, parameter :: r133b = 298
  integer, parameter :: r134b = 299
  integer, parameter :: r135b = 300
  integer, parameter :: r136b = 301
  integer, parameter :: r137b = 302
  integer, parameter :: r138b = 303
  integer, parameter :: r139b = 304
  integer, parameter :: r140b = 305
  integer, parameter :: r141b = 306
  integer, parameter :: r142b = 307
  integer, parameter :: r143b = 308
  integer, parameter :: r144b = 309
  integer, parameter :: r145b = 310
  integer, parameter :: r146b = 311
  integer, parameter :: r147b = 312
  integer, parameter :: r148b = 313
  integer, parameter :: r149b = 314
  integer, parameter :: r150b = 315
  integer, parameter :: r151b = 316
  integer, parameter :: r153b = 317
  integer, parameter :: r154b = 318
  integer, parameter :: r157b = 319
  integer, parameter :: r159b = 320
  integer, parameter :: r160b = 321
  integer, parameter :: r161b = 322
  integer, parameter :: r165b = 323
  integer, parameter :: r166b = 324
  
  ! Index of third body species
  integer, parameter :: mM1 = 1
  integer, parameter :: mM30 = 2
  integer, parameter :: mM35 = 3
  integer, parameter :: mM38 = 4
  integer, parameter :: mM144 = 5
  
  integer, parameter :: mM11 = 6
  integer, parameter :: mM44 = 7
  integer, parameter :: mM46 = 8
  integer, parameter :: mM49 = 9
  integer, parameter :: mM50 = 10
  integer, parameter :: mM55 = 11
  integer, parameter :: mM62 = 12
  integer, parameter :: mM64 = 13
  integer, parameter :: mM66 = 14
  integer, parameter :: mM72 = 15
  integer, parameter :: mM74 = 16
  integer, parameter :: mM83 = 17
  integer, parameter :: mM112 = 18
  integer, parameter :: mM119 = 19
  integer, parameter :: mM124 = 20
  integer, parameter :: mM135 = 21
  integer, parameter :: mM149 = 22
  integer, parameter :: mM161 = 23
  
  ! Index of third body reactions
  integer, parameter :: TBr1f = 1
  integer, parameter :: TBr30f = 2
  integer, parameter :: TBr35f = 3
  integer, parameter :: TBr38f = 4
  integer, parameter :: TBr144f = 5
  integer, parameter :: TBr1b = 6
  integer, parameter :: TBr30b = 7
  integer, parameter :: TBr35b = 8
  integer, parameter :: TBr38b = 9
  integer, parameter :: TBr144b = 10
  
  ! Index of fall off reactions
  integer, parameter :: FOr11f = 1
  integer, parameter :: FOr44f = 2
  integer, parameter :: FOr46f = 3
  integer, parameter :: FOr49f = 4
  integer, parameter :: FOr50f = 5
  integer, parameter :: FOr55f = 6
  integer, parameter :: FOr62f = 7
  integer, parameter :: FOr64f = 8
  integer, parameter :: FOr66f = 9
  integer, parameter :: FOr72f = 10
  integer, parameter :: FOr74f = 11
  integer, parameter :: FOr83f = 12
  integer, parameter :: FOr112f = 13
  integer, parameter :: FOr119f = 14
  integer, parameter :: FOr124f = 15
  integer, parameter :: FOr135f = 16
  integer, parameter :: FOr149f = 17
  integer, parameter :: FOr161f = 18
  integer, parameter :: FOr11b = 19
  integer, parameter :: FOr44b = 20
  integer, parameter :: FOr46b = 21
  integer, parameter :: FOr49b = 22
  integer, parameter :: FOr50b = 23
  integer, parameter :: FOr55b = 24
  integer, parameter :: FOr62b = 25
  integer, parameter :: FOr64b = 26
  integer, parameter :: FOr66b = 27
  integer, parameter :: FOr72b = 28
  integer, parameter :: FOr74b = 29
  integer, parameter :: FOr83b = 30
  integer, parameter :: FOr112b = 31
  integer, parameter :: FOr119b = 32
  integer, parameter :: FOr124b = 33
  integer, parameter :: FOr135b = 34
  integer, parameter :: FOr149b = 35
  integer, parameter :: FOr161b = 36
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.00201588_pr, & ! H2
       0.00100794_pr, & ! H
       0.0159994_pr, & ! O
       0.0319988_pr, & ! O2
       0.01700734_pr, & ! OH
       0.01801528_pr, & ! H2O
       0.03300674_pr, & ! HO2
       0.03401468_pr, & ! H2O2
       0.01503482_pr, & ! CH3
       0.01604276_pr, & ! CH4
       0.0280104_pr, & ! CO
       0.0440098_pr, & ! CO2
       0.03002628_pr, & ! CH2O
       0.03103422_pr, & ! CH3O
       0.03204216_pr, & ! CH3OH
       0.02603788_pr, & ! C2H2
       0.02805376_pr, & ! C2H4
       0.03006964_pr, & ! C2H6
       0.04203728_pr, & ! CH2CO
       0.02801348_pr & ! N2
      !0.01301894_pr, & ! CH
      !0.01402688_pr, & ! CH2
      !0.01402688_pr, & ! CH2GSG
      !0.02901834_pr, & ! HCO
      !0.03103422_pr, & ! CH2OH
      !0.02704582_pr, & ! C2H3
      !0.0290617_pr, & ! C2H5
      !0.04102934_pr, & ! HCCO
      !0.04304522_pr & ! CH2CHO
       /)
  
contains

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

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

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

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

  end function getlindratecoeff

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

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

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

    M(mM30) = (0.5_pr)*c(sC2H6) &
         + (-0.25_pr)*c(sCO) &
         + (0.5_pr)*c(sCO2) &
         + (-1_pr)*c(sH2O) &
         + (-1_pr)*c(sN2) &
         + (-1_pr)*c(sO2) &
         + sum(c)

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

    M(mM38) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (-0.27_pr)*c(sH2) &
         + (2.65_pr)*c(sH2O) &
         + sum(c)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    k(r1f) = (500000_pr)*exp((-0_pr)*R_T_inv + T_log *(-1_pr))
    k(r2f) = (0.0387_pr)*exp((-26191.84_pr)*R_T_inv + T_log *(2.7_pr))
    k(r3f) = (20000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r4f) = (9.63_pr)*exp((-16736_pr)*R_T_inv + T_log *(2_pr))
    k(r5f) = (57000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r6f) = (80000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r7f) = (15000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r8f) = (15000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r9f) = (50600000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r10f) = (1020_pr)*exp((-35982.4_pr)*R_T_inv + T_log *(1.5_pr))
    k_0(FOr11f) = (602_pr)*exp((-12552_pr)*R_T_inv + T_log *(0_pr))
    k_inf(FOr11f) = (18000_pr)*exp((-9978.84_pr)*R_T_inv + T_log *(0_pr))
    k(r11f) = getlindratecoeff(Tloc,k_0(FOr11f),k_inf(FOr11f),1.0_pr, M(mM11),Ploc)
    k(r12f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r13f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r14f) = (39000000_pr)*exp((-14811.36_pr)*R_T_inv + T_log *(0_pr))
    k(r15f) = (10000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r16f) = (10000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r17f) = (0.388_pr)*exp((-12970.4_pr)*R_T_inv + T_log *(2.5_pr))
    k(r18f) = (0.13_pr)*exp((-20920_pr)*R_T_inv + T_log *(2.5_pr))
    k(r19f) = (13.5_pr)*exp((-7949.6_pr)*R_T_inv + T_log *(2_pr))
    k(r20f) = (6.94_pr)*exp((-7949.6_pr)*R_T_inv + T_log *(2_pr))
    k(r21f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r22f) = (12.5_pr)*exp((-920.48_pr)*R_T_inv + T_log *(1.83_pr))
    k(r23f) = (22400000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r24f) = (89.8_pr)*exp((-23806.96_pr)*R_T_inv + T_log *(1.92_pr))
    k(r25f) = (100000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r26f) = (10000000_pr)*exp((-33472_pr)*R_T_inv + T_log *(0_pr))
    k(r27f) = (1750000_pr)*exp((-5648.4_pr)*R_T_inv + T_log *(0_pr))
    k(r28f) = (2500000_pr)*exp((-199995.2_pr)*R_T_inv + T_log *(0_pr))
    k(r29f) = (100000000_pr)*exp((-167360_pr)*R_T_inv + T_log *(0_pr))
    k(r30f) = (2800000_pr)*exp((-0_pr)*R_T_inv + T_log *(-0.86_pr))
    k(r31f) = (20800000_pr)*exp((-0_pr)*R_T_inv + T_log *(-1.24_pr))
    k(r32f) = (11260000_pr)*exp((-0_pr)*R_T_inv + T_log *(-0.76_pr))
    k(r33f) = (26000000_pr)*exp((-0_pr)*R_T_inv + T_log *(-1.24_pr))
    k(r34f) = (26500000000_pr)*exp((-71299.544_pr)*R_T_inv + T_log *(-0.6707_pr))
    k(r35f) = (1000000_pr)*exp((-0_pr)*R_T_inv + T_log *(-1_pr))
    k(r36f) = (90000_pr)*exp((-0_pr)*R_T_inv + T_log *(-0.6_pr))
    k(r37f) = (60000000_pr)*exp((-0_pr)*R_T_inv + T_log *(-1.25_pr))
    k(r38f) = (22000000000_pr)*exp((-0_pr)*R_T_inv + T_log *(-2_pr))
    k(r39f) = (3970000_pr)*exp((-2807.464_pr)*R_T_inv + T_log *(0_pr))
    k(r40f) = (44800000_pr)*exp((-4468.512_pr)*R_T_inv + T_log *(0_pr))
    k(r41f) = (84000000_pr)*exp((-2656.84_pr)*R_T_inv + T_log *(0_pr))
    k(r42f) = (12.1_pr)*exp((-21756.8_pr)*R_T_inv + T_log *(2_pr))
    k(r43f) = (10000000_pr)*exp((-15062.4_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr44f) = (104000000000000_pr)*exp((-6694.4_pr)*R_T_inv + T_log *(-2.76_pr))
    k_inf(FOr44f) = (600000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    FC(FOr44f) = ((1.0_pr - 0.562_pr)*exp(-Tloc/(91_pr))) + (0.562_pr*exp(-Tloc/(5836_pr)))&
              + exp(-(8552_pr)/Tloc)
    k(r44f) = getlindratecoeff(Tloc,k_0(FOr44f),k_inf(FOr44f),FC(FOr44f), M(mM44),Ploc)
    k(r45f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr46f) = (2.62e+21_pr)*exp((-10208.96_pr)*R_T_inv + T_log *(-4.76_pr))
    k_inf(FOr46f) = (13900000000_pr)*exp((-2242.624_pr)*R_T_inv + T_log *(-0.534_pr))
    FC(FOr46f) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r46f) = getlindratecoeff(Tloc,k_0(FOr46f),k_inf(FOr46f),FC(FOr46f), M(mM46),Ploc)
    k(r47f) = (660_pr)*exp((-45354.56_pr)*R_T_inv + T_log *(1.62_pr))
    k(r48f) = (73400000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr49f) = (1.27e+20_pr)*exp((-27321.52_pr)*R_T_inv + T_log *(-4.82_pr))
    k_inf(FOr49f) = (540000_pr)*exp((-15062.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr49f) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r49f) = getlindratecoeff(Tloc,k_0(FOr49f),k_inf(FOr49f),FC(FOr49f), M(mM49),Ploc)
    k_0(FOr50f) = (2.2e+18_pr)*exp((-23263.04_pr)*R_T_inv + T_log *(-4.8_pr))
    k_inf(FOr50f) = (540000_pr)*exp((-10878.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr50f) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r50f) = getlindratecoeff(Tloc,k_0(FOr50f),k_inf(FOr50f),FC(FOr50f), M(mM50),Ploc)
    k(r51f) = (57.4_pr)*exp((-11472.528_pr)*R_T_inv + T_log *(1.9_pr))
    k(r52f) = (20000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r53f) = (165000_pr)*exp((1188.256_pr)*R_T_inv + T_log *(0.65_pr))
    k(r54f) = (32800000_pr)*exp((-2552.24_pr)*R_T_inv + T_log *(-0.09_pr))
    k_0(FOr55f) = (4.66e+29_pr)*exp((-58910.72_pr)*R_T_inv + T_log *(-7.44_pr))
    k_inf(FOr55f) = (2430000_pr)*exp((-209.2_pr)*R_T_inv + T_log *(0.515_pr))
    FC(FOr55f) = ((1.0_pr - 0.7_pr)*exp(-Tloc/(100_pr))) + (0.7_pr*exp(-Tloc/(90000_pr)))&
              + exp(-(10000_pr)/Tloc)
    k(r55f) = getlindratecoeff(Tloc,k_0(FOr55f),k_inf(FOr55f),FC(FOr55f), M(mM55),Ploc)
    k(r56f) = (41.5_pr)*exp((-8050.016_pr)*R_T_inv + T_log *(1.63_pr))
    k(r57f) = (20000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r58f) = (1500000_pr)*exp((460.24_pr)*R_T_inv + T_log *(0.5_pr))
    k(r59f) = (262000000_pr)*exp((-4476.88_pr)*R_T_inv + T_log *(-0.23_pr))
    k(r60f) = (17_pr)*exp((-20376.08_pr)*R_T_inv + T_log *(2.1_pr))
    k(r61f) = (4.2_pr)*exp((-20376.08_pr)*R_T_inv + T_log *(2.1_pr))
    k_0(FOr62f) = (3.8e+28_pr)*exp((-30208.48_pr)*R_T_inv + T_log *(-7.27_pr))
    k_inf(FOr62f) = (5600000_pr)*exp((-10041.6_pr)*R_T_inv + T_log *(0_pr))
    FC(FOr62f) = ((1.0_pr - 0.7507_pr)*exp(-Tloc/(98.5_pr))) + (0.7507_pr*exp(-Tloc/(1302_pr)))&
              + exp(-(4167_pr)/Tloc)
    k(r62f) = getlindratecoeff(Tloc,k_0(FOr62f),k_inf(FOr62f),FC(FOr62f), M(mM62),Ploc)
    k(r63f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr64f) = (6e+29_pr)*exp((-29162.48_pr)*R_T_inv + T_log *(-7.62_pr))
    k_inf(FOr64f) = (540000_pr)*exp((-7614.88_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr64f) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r64f) = getlindratecoeff(Tloc,k_0(FOr64f),k_inf(FOr64f),FC(FOr64f), M(mM64),Ploc)
    k(r65f) = (1.325_pr)*exp((-51212.16_pr)*R_T_inv + T_log *(2.53_pr))
    k_0(FOr66f) = (1.99e+29_pr)*exp((-27970.04_pr)*R_T_inv + T_log *(-7.08_pr))
    k_inf(FOr66f) = (521000000000_pr)*exp((-6610.72_pr)*R_T_inv + T_log *(-0.99_pr))
    FC(FOr66f) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r66f) = getlindratecoeff(Tloc,k_0(FOr66f),k_inf(FOr66f),FC(FOr66f), M(mM66),Ploc)
    k(r67f) = (2000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r68f) = (115_pr)*exp((-31505.52_pr)*R_T_inv + T_log *(1.9_pr))
    k(r69f) = (100000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r70f) = (50000000_pr)*exp((-33472_pr)*R_T_inv + T_log *(0_pr))
    k(r71f) = (11300000_pr)*exp((-14342.752_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr72f) = (5.07e+15_pr)*exp((-352920.4_pr)*R_T_inv + T_log *(-3.42_pr))
    k_inf(FOr72f) = (43_pr)*exp((-333046.4_pr)*R_T_inv + T_log *(1.5_pr))
    FC(FOr72f) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r72f) = getlindratecoeff(Tloc,k_0(FOr72f),k_inf(FOr72f),FC(FOr72f), M(mM72),Ploc)
    k(r73f) = (216_pr)*exp((-14351.12_pr)*R_T_inv + T_log *(1.51_pr))
    k_0(FOr74f) = (2300000_pr)*exp((7112.8_pr)*R_T_inv + T_log *(-0.9_pr))
    k_inf(FOr74f) = (74000000_pr)*exp((-0_pr)*R_T_inv + T_log *(-0.37_pr))
    FC(FOr74f) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r74f) = getlindratecoeff(Tloc,k_0(FOr74f),k_inf(FOr74f),FC(FOr74f), M(mM74),Ploc)
    k(r75f) = (0.0357_pr)*exp((8828.24_pr)*R_T_inv + T_log *(2.4_pr))
    k(r76f) = (14500000_pr)*exp((2092_pr)*R_T_inv + T_log *(0_pr))
    k(r77f) = (2000000_pr)*exp((-1786.568_pr)*R_T_inv + T_log *(0_pr))
    k(r78f) = (1700000000000_pr)*exp((-123051.44_pr)*R_T_inv + T_log *(0_pr))
    k(r79f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r80f) = (20000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r81f) = (11.3_pr)*exp((-12552_pr)*R_T_inv + T_log *(2_pr))
    k(r82f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr83f) = (4e+24_pr)*exp((-13137.76_pr)*R_T_inv + T_log *(-5.92_pr))
    k_inf(FOr83f) = (2790000000000_pr)*exp((-5564.72_pr)*R_T_inv + T_log *(-1.43_pr))
    FC(FOr83f) = ((1.0_pr - 0.412_pr)*exp(-Tloc/(195_pr))) + (0.412_pr*exp(-Tloc/(5900_pr)))&
              + exp(-(6394_pr)/Tloc)
    k(r83f) = getlindratecoeff(Tloc,k_0(FOr83f),k_inf(FOr83f),FC(FOr83f), M(mM83),Ploc)
    k(r84f) = (56_pr)*exp((-22677.28_pr)*R_T_inv + T_log *(1.6_pr))
    k(r85f) = (644000000000_pr)*exp((-5928.728_pr)*R_T_inv + T_log *(-1.34_pr))
    k(r86f) = (100_pr)*exp((-13054.08_pr)*R_T_inv + T_log *(1.6_pr))
    k(r87f) = (47.6_pr)*exp((-292.88_pr)*R_T_inv + T_log *(1.228_pr))
    k(r88f) = (50000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r89f) = (3430_pr)*exp((1870.248_pr)*R_T_inv + T_log *(1.18_pr))
    k(r90f) = (5000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r91f) = (5000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r92f) = (1.44_pr)*exp((3514.56_pr)*R_T_inv + T_log *(2_pr))
    k(r93f) = (6.3_pr)*exp((-6276_pr)*R_T_inv + T_log *(2_pr))
    k(r94f) = (2.18e-10_pr)*exp((4184_pr)*R_T_inv + T_log *(4.5_pr))
    k(r95f) = (4.83e-10_pr)*exp((8368_pr)*R_T_inv + T_log *(4_pr))
    k(r96f) = (5000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r97f) = (3.6_pr)*exp((-10460_pr)*R_T_inv + T_log *(2_pr))
    k(r98f) = (3.54_pr)*exp((-3640.08_pr)*R_T_inv + T_log *(2.12_pr))
    k(r99f) = (7500000_pr)*exp((-8368_pr)*R_T_inv + T_log *(0_pr))
    k(r100f) = (130000_pr)*exp((6819.92_pr)*R_T_inv + T_log *(0_pr))
    k(r101f) = (420000000_pr)*exp((-50208_pr)*R_T_inv + T_log *(0_pr))
    k(r102f) = (20000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r103f) = (1000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r104f) = (37800000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r105f) = (150000000_pr)*exp((-98742.4_pr)*R_T_inv + T_log *(0_pr))
    k(r106f) = (5.6_pr)*exp((-50208_pr)*R_T_inv + T_log *(2_pr))
    k(r107f) = (67100000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r108f) = (108000000_pr)*exp((-13012.24_pr)*R_T_inv + T_log *(0_pr))
    k(r109f) = (5710000_pr)*exp((3158.92_pr)*R_T_inv + T_log *(0_pr))
    k(r110f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r111f) = (60000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr112f) = (2.69e+16_pr)*exp((-8100.224_pr)*R_T_inv + T_log *(-3.74_pr))
    k_inf(FOr112f) = (50000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    FC(FOr112f) = ((1.0_pr - 0.5757_pr)*exp(-Tloc/(237_pr))) + (0.5757_pr*exp(-Tloc/(1652_pr)))&
              + exp(-(5069_pr)/Tloc)
    k(r112f) = getlindratecoeff(Tloc,k_0(FOr112f),k_inf(FOr112f),FC(FOr112f), M(mM112),Ploc)
    k(r113f) = (190000000_pr)*exp((-66073.728_pr)*R_T_inv + T_log *(0_pr))
    k(r114f) = (94600000_pr)*exp((2154.76_pr)*R_T_inv + T_log *(0_pr))
    k(r115f) = (5000000_pr)*exp((-6276_pr)*R_T_inv + T_log *(0_pr))
    k(r116f) = (0.5_pr)*exp((-30250.32_pr)*R_T_inv + T_log *(2_pr))
    k(r117f) = (40000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r118f) = (2.46_pr)*exp((-34601.68_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr119f) = (2.69e+21_pr)*exp((-29685.48_pr)*R_T_inv + T_log *(-5.11_pr))
    k_inf(FOr119f) = (810000_pr)*exp((-18869.84_pr)*R_T_inv + T_log *(0.5_pr))
    FC(FOr119f) = ((1.0_pr - 0.5907_pr)*exp(-Tloc/(275_pr))) + (0.5907_pr*exp(-Tloc/(1226_pr)))&
              + exp(-(5185_pr)/Tloc)
    k(r119f) = getlindratecoeff(Tloc,k_0(FOr119f),k_inf(FOr119f),FC(FOr119f), M(mM119),Ploc)
    k(r120f) = (15000000_pr)*exp((-2510.4_pr)*R_T_inv + T_log *(0_pr))
    k(r121f) = (28000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r122f) = (12000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r123f) = (70000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k_0(FOr124f) = (1.88e+26_pr)*exp((-21087.36_pr)*R_T_inv + T_log *(-6.36_pr))
    k_inf(FOr124f) = (482000000000_pr)*exp((-4790.68_pr)*R_T_inv + T_log *(-1.16_pr))
    FC(FOr124f) = ((1.0_pr - 0.6027_pr)*exp(-Tloc/(208_pr))) + (0.6027_pr*exp(-Tloc/(3922_pr)))&
              + exp(-(10180_pr)/Tloc)
    k(r124f) = getlindratecoeff(Tloc,k_0(FOr124f),k_inf(FOr124f),FC(FOr124f), M(mM124),Ploc)
    k(r125f) = (30000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r126f) = (12000000_pr)*exp((2384.88_pr)*R_T_inv + T_log *(0_pr))
    k(r127f) = (16000000_pr)*exp((2384.88_pr)*R_T_inv + T_log *(0_pr))
    k(r128f) = (9000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r129f) = (7000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r130f) = (14000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r131f) = (40000000_pr)*exp((2301.2_pr)*R_T_inv + T_log *(0_pr))
    k(r132f) = (35600000_pr)*exp((-127528.32_pr)*R_T_inv + T_log *(0_pr))
    k(r133f) = (2310000_pr)*exp((-84997.96_pr)*R_T_inv + T_log *(0_pr))
    k(r134f) = (0.0245_pr)*exp((-21673.12_pr)*R_T_inv + T_log *(2.47_pr))
    k_0(FOr135f) = (3.4e+29_pr)*exp((-11556.208_pr)*R_T_inv + T_log *(-7.03_pr))
    k_inf(FOr135f) = (67700000000_pr)*exp((-2736.336_pr)*R_T_inv + T_log *(-1.18_pr))
    FC(FOr135f) = ((1.0_pr - 0.619_pr)*exp(-Tloc/(73.2_pr))) + (0.619_pr*exp(-Tloc/(1180_pr)))&
              + exp(-(9999_pr)/Tloc)
    k(r135f) = getlindratecoeff(Tloc,k_0(FOr135f),k_inf(FOr135f),FC(FOr135f), M(mM135),Ploc)
    k(r136f) = (6840000_pr)*exp((-44350.4_pr)*R_T_inv + T_log *(0.1_pr))
    k(r137f) = (26480000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r138f) = (0.00332_pr)*exp((-24518.24_pr)*R_T_inv + T_log *(2.81_pr))
    k(r139f) = (30_pr)*exp((-41588.96_pr)*R_T_inv + T_log *(1.5_pr))
    k(r140f) = (10_pr)*exp((-41588.96_pr)*R_T_inv + T_log *(1.5_pr))
    k(r141f) = (0.227_pr)*exp((-38492.8_pr)*R_T_inv + T_log *(2_pr))
    k(r142f) = (6.14_pr)*exp((-43722.8_pr)*R_T_inv + T_log *(1.74_pr))
    k(r143f) = (1500000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r144f) = (187000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r145f) = (13450000_pr)*exp((-1673.6_pr)*R_T_inv + T_log *(0_pr))
    k(r146f) = (18000000_pr)*exp((-3765.6_pr)*R_T_inv + T_log *(0_pr))
    k(r147f) = (4.28e-19_pr)*exp((14769.52_pr)*R_T_inv + T_log *(7.6_pr))
    k(r148f) = (45800000000_pr)*exp((-4246.76_pr)*R_T_inv + T_log *(-1.39_pr))
    k_0(FOr149f) = (1.58e+45_pr)*exp((-409195.2_pr)*R_T_inv + T_log *(-9.3_pr))
    k_inf(FOr149f) = (8000000000000_pr)*exp((-363045.68_pr)*R_T_inv + T_log *(0.44_pr))
    FC(FOr149f) = ((1.0_pr - 0.7345_pr)*exp(-Tloc/(180_pr))) + (0.7345_pr*exp(-Tloc/(1035_pr)))&
              + exp(-(5417_pr)/Tloc)
    k(r149f) = getlindratecoeff(Tloc,k_0(FOr149f),k_inf(FOr149f),FC(FOr149f), M(mM149),Ploc)
    k(r150f) = (840000_pr)*exp((-16213_pr)*R_T_inv + T_log *(0_pr))
    k(r151f) = (3200000_pr)*exp((-3573.136_pr)*R_T_inv + T_log *(0_pr))
    k(r152f) = (33700000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r153f) = (6.7_pr)*exp((-920.48_pr)*R_T_inv + T_log *(1.83_pr))
    k(r154f) = (5000000000_pr)*exp((-72508.72_pr)*R_T_inv + T_log *(0_pr))
    k(r155f) = (8000_pr)*exp((7342.92_pr)*R_T_inv + T_log *(0.5_pr))
    k(r156f) = (5800000_pr)*exp((-6276_pr)*R_T_inv + T_log *(0_pr))
    k(r157f) = (2400000_pr)*exp((-6276_pr)*R_T_inv + T_log *(0_pr))
    k(r158f) = (68200_pr)*exp((3912.04_pr)*R_T_inv + T_log *(0.25_pr))
    k(r159f) = (303000_pr)*exp((-46.024_pr)*R_T_inv + T_log *(0.29_pr))
    k(r160f) = (1.337_pr)*exp((1606.656_pr)*R_T_inv + T_log *(1.61_pr))
    k_0(FOr161f) = (1.012e+30_pr)*exp((-16125.136_pr)*R_T_inv + T_log *(-7.63_pr))
    k_inf(FOr161f) = (486500_pr)*exp((7342.92_pr)*R_T_inv + T_log *(0.422_pr))
    FC(FOr161f) = ((1.0_pr - 0.465_pr)*exp(-Tloc/(201_pr))) + (0.465_pr*exp(-Tloc/(1773_pr)))&
              + exp(-(5333_pr)/Tloc)
    k(r161f) = getlindratecoeff(Tloc,k_0(FOr161f),k_inf(FOr161f),FC(FOr161f), M(mM161),Ploc)
    k(r162f) = (150000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r163f) = (18100_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r164f) = (23500_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r165f) = (22000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r166f) = (11000000_pr)*exp((-0_pr)*R_T_inv + T_log *(0_pr))
    k(r1b) = (77654354800.6578_pr)*exp((-425035.941446839_pr)*R_T_inv + T_log *(-0.705734671876919_pr))
    k(r2b) = (0.020711486094766_pr)*exp((-18168.1250206447_pr)*R_T_inv + T_log *(2.67740604415194_pr))
    k(r3b) = (2921026.85047749_pr)*exp((-221089.458553749_pr)*R_T_inv + T_log *(0.27853070856989_pr))
    k(r4b) = (0.0370460100167413_pr)*exp((-76384.9465235186_pr)*R_T_inv + T_log *(2.49231966612624_pr))
    k(r5b) = (2545377683.25348_pr)*exp((-739065.845262563_pr)*R_T_inv + T_log *(0.00489705077585355_pr))
    k(r6b) = (3972879879.95113_pr)*exp((-382652.015363043_pr)*R_T_inv + T_log *(-0.262036531165194_pr))
    k(r7b) = (77331278.5661296_pr)*exp((-788703.34955258_pr)*R_T_inv + T_log *(0.118163721712543_pr))
    k(r8b) = (655739839.118918_pr)*exp((-420230.742822539_pr)*R_T_inv + T_log *(-0.345821927004021_pr))
    k(r9b) = (2649502503.69646_pr)*exp((-287836.718042318_pr)*R_T_inv + T_log *(-0.162246817165187_pr))
    k(r10b) = (1.32757125025021_pr)*exp((-21232.1347562493_pr)*R_T_inv + T_log *(1.83016824693274_pr))
    k_0(FOr11b) = (165832885694719_pr)*exp((-545853.407002819_pr)*R_T_inv + T_log *(-0.844735733613503_pr))
    k_inf(FOr11b) = (4.95845665144248e+15_pr)*exp((-543280.246747439_pr)*R_T_inv + T_log *(-0.844735688537068_pr))
    k(r11b) = getlindratecoeff(Tloc,k_0(FOr11b),k_inf(FOr11b),1.0_pr, M(mM11),Ploc)
    k(r12b) = (1893407.97207187_pr)*exp((-360448.890124488_pr)*R_T_inv + T_log *(0.441391983213283_pr))
    k(r13b) = (3358323969035.05_pr)*exp((-468714.355731817_pr)*R_T_inv + T_log *(-0.697609088210368_pr))
    k(r14b) = (119613.617390403_pr)*exp((-71947.3233354498_pr)*R_T_inv + T_log *(0.28289158761451_pr))
    k(r15b) = (131156.678235481_pr)*exp((-302216.616998152_pr)*R_T_inv + T_log *(0.699891154241341_pr))
    k(r16b) = (41088.024310921_pr)*exp((-332206.055471112_pr)*R_T_inv + T_log *(0.524092187918471_pr))
    k(r17b) = (0.00784561153720865_pr)*exp((-34557.0446160389_pr)*R_T_inv + T_log *(2.50068841357956_pr))
    k(r18b) = (0.00839099629785807_pr)*exp((-12517.2060603269_pr)*R_T_inv + T_log *(2.67648739419993_pr))
    k(r19b) = (34.425971560248_pr)*exp((-89876.6795560546_pr)*R_T_inv + T_log *(1.84921923312609_pr))
    k(r20b) = (2.84005468976049e-05_pr)*exp((-199497.539519457_pr)*R_T_inv + T_log *(3.26463599023883_pr))
    k(r21b) = (163304335013.597_pr)*exp((-379586.586383181_pr)*R_T_inv + T_log *(-0.625650501027362_pr))
    k(r22b) = (0.000528132947309131_pr)*exp((-110087.274141435_pr)*R_T_inv + T_log *(2.55673460427883_pr))
    k(r23b) = (39533.8024700913_pr)*exp((-326578.49327088_pr)*R_T_inv + T_log *(0.815547905891994_pr))
    k(r24b) = (0.029834748586851_pr)*exp((-28803.1796501584_pr)*R_T_inv + T_log *(2.3171844499336_pr))
    k(r25b) = (0.00323860187147881_pr)*exp((-427685.824304008_pr)*R_T_inv + T_log *(1.30050682572713_pr))
    k(r26b) = (11716.0658310934_pr)*exp((-15886.7999454295_pr)*R_T_inv + T_log *(0.452044089793627_pr))
    k(r27b) = (5835.96939427452_pr)*exp((-205949.525469324_pr)*R_T_inv + T_log *(0.728459784304265_pr))
    k(r28b) = (16787969429.9696_pr)*exp((-236325.25707961_pr)*R_T_inv + T_log *(-0.78223221538397_pr))
    k(r29b) = (2099958.33529657_pr)*exp((-3406.50514375938_pr)*R_T_inv + T_log *(0.00436081517845624_pr))
    k(r30b) = (2977474913454_pr)*exp((-203946.482593931_pr)*R_T_inv + T_log *(-0.844265327199912_pr))
    k(r31b) = (22118385838470.1_pr)*exp((-203946.482619027_pr)*R_T_inv + T_log *(-1.22426533170906_pr))
    k(r32b) = (11973704134340.6_pr)*exp((-203946.482682537_pr)*R_T_inv + T_log *(-0.744265342872238_pr))
    k(r33b) = (27647988031676.5_pr)*exp((-203946.482770268_pr)*R_T_inv + T_log *(-1.22426535864353_pr))
    k(r34b) = (100328807.717269_pr)*exp((635.864533327733_pr)*R_T_inv + T_log *(-0.313931145463416_pr))
    k(r35b) = (290198724209.349_pr)*exp((-433059.656399876_pr)*R_T_inv + T_log *(-0.683140711293222_pr))
    k(r36b) = (26117885254.1253_pr)*exp((-433059.656406948_pr)*R_T_inv + T_log *(-0.283140711343168_pr))
    k(r37b) = (17411922435263.5_pr)*exp((-433059.656356855_pr)*R_T_inv + T_log *(-0.933140703723797_pr))
    k(r38b) = (4.2040214411069e+16_pr)*exp((-497116.466748875_pr)*R_T_inv + T_log *(-1.72925846822621_pr))
    k(r39b) = (27009.8734999967_pr)*exp((-224042.039802433_pr)*R_T_inv + T_log *(0.611775682478179_pr))
    k(r40b) = (12225970.3666321_pr)*exp((-233581.685647792_pr)*R_T_inv + T_log *(0.301124644121301_pr))
    k(r41b) = (46447.7600307324_pr)*exp((-151810.890291503_pr)*R_T_inv + T_log *(0.635299516067345_pr))
    k(r42b) = (0.0869761645463356_pr)*exp((-89429.4615470536_pr)*R_T_inv + T_log *(2.51491361425604_pr))
    k(r43b) = (261.726248328552_pr)*exp((-295945.922469627_pr)*R_T_inv + T_log *(1.10409532366667_pr))
    k_0(FOr44b) = (4.99476042263987e+21_pr)*exp((-469409.675261043_pr)*R_T_inv + T_log *(-2.84841594333284_pr))
    k_inf(FOr44b) = (2.8815737210813e+16_pr)*exp((-462715.270485629_pr)*R_T_inv + T_log *(-0.0884150950219076_pr))
    FC(FOr44b) = ((1.0_pr - 0.562_pr)*exp(-Tloc/(91_pr))) + (0.562_pr*exp(-Tloc/(5836_pr)))&
              + exp(-(8552_pr)/Tloc)
    k(r44b) = getlindratecoeff(Tloc,k_0(FOr44b),k_inf(FOr44b),FC(FOr44b), M(mM44),Ploc)
    k(r45b) = (3463429.66981424_pr)*exp((-49637.5019156926_pr)*R_T_inv + T_log *(0.113267102167906_pr))
    k_0(FOr46b) = (3.12636203702489e+29_pr)*exp((-449995.16640668_pr)*R_T_inv + T_log *(-4.79590286897603_pr))
    k_inf(FOr46b) = (1.65863867570345e+18_pr)*exp((-442028.828745086_pr)*R_T_inv + T_log *(-0.569902573139977_pr))
    FC(FOr46b) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r46b) = getlindratecoeff(Tloc,k_0(FOr46b),k_inf(FOr46b),FC(FOr46b), M(mM46),Ploc)
    k(r47b) = (1.60509899443808_pr)*exp((-38628.0106234292_pr)*R_T_inv + T_log *(1.97276204915205_pr))
    k(r48b) = (8656028.67291496_pr)*exp((-368472.605058181_pr)*R_T_inv + T_log *(0.463985946036272_pr))
    k_0(FOr49b) = (1.50386582713963e+27_pr)*exp((-150140.84439361_pr)*R_T_inv + T_log *(-5.22562581644432_pr))
    k_inf(FOr49b) = (6394390882939.78_pr)*exp((-137881.724479564_pr)*R_T_inv + T_log *(0.048374168294516_pr))
    FC(FOr49b) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r49b) = getlindratecoeff(Tloc,k_0(FOr49b),k_inf(FOr49b),FC(FOr49b), M(mM49),Ploc)
    k_0(FOr50b) = (8.31578323990648e+25_pr)*exp((-116092.925852392_pr)*R_T_inv + T_log *(-5.02982683814769_pr))
    k_inf(FOr50b) = (20411471245921_pr)*exp((-103708.285970061_pr)*R_T_inv + T_log *(0.224173140894518_pr))
    FC(FOr50b) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r50b) = getlindratecoeff(Tloc,k_0(FOr50b),k_inf(FOr50b),FC(FOr50b), M(mM50),Ploc)
    k(r51b) = (0.328948431745645_pr)*exp((-76632.206665515_pr)*R_T_inv + T_log *(2.2054854817067_pr))
    k(r52b) = (490139.943386981_pr)*exp((-310240.331959117_pr)*R_T_inv + T_log *(0.7224851136321_pr))
    k(r53b) = (41.3295324676424_pr)*exp((-13191.6428841233_pr)*R_T_inv + T_log *(1.51213798438914_pr))
    k(r54b) = (371.352666220173_pr)*exp((-13754.603077139_pr)*R_T_inv + T_log *(1.2150808256898_pr))
    k_0(FOr55b) = (1.12127343712802e+36_pr)*exp((-492349.45535876_pr)*R_T_inv + T_log *(-7.32222206101444_pr))
    k_inf(FOr55b) = (5846984571650.62_pr)*exp((-433647.935455226_pr)*R_T_inv + T_log *(0.632777921864153_pr))
    FC(FOr55b) = ((1.0_pr - 0.7_pr)*exp(-Tloc/(100_pr))) + (0.7_pr*exp(-Tloc/(90000_pr)))&
              + exp(-(10000_pr)/Tloc)
    k(r55b) = getlindratecoeff(Tloc,k_0(FOr55b),k_inf(FOr55b),FC(FOr55b), M(mM55),Ploc)
    k(r56b) = (13.0008880387493_pr)*exp((-38039.4545949244_pr)*R_T_inv + T_log *(1.45420101240349_pr))
    k(r57b) = (153548.267319829_pr)*exp((-340229.770384813_pr)*R_T_inv + T_log *(0.546686155125088_pr))
    k(r58b) = (117.704397538656_pr)*exp((-43909.0974090943_pr)*R_T_inv + T_log *(1.18633900865301_pr))
    k(r59b) = (929.263533042616_pr)*exp((-45668.6816354582_pr)*R_T_inv + T_log *(0.899281844320309_pr))
    k(r60b) = (0.642308449553108_pr)*exp((-49986.4394844425_pr)*R_T_inv + T_log *(2.1232823887815_pr))
    k(r61b) = (0.506546290695425_pr)*exp((-19997.0009862452_pr)*R_T_inv + T_log *(2.29908135945212_pr))
    k_0(FOr62b) = (2.35979139542732e+33_pr)*exp((-175170.114727901_pr)*R_T_inv + T_log *(-6.95290903647189_pr))
    k_inf(FOr62b) = (347758723791.028_pr)*exp((-155003.234711258_pr)*R_T_inv + T_log *(0.317090966586302_pr))
    FC(FOr62b) = ((1.0_pr - 0.7507_pr)*exp(-Tloc/(98.5_pr))) + (0.7507_pr*exp(-Tloc/(1302_pr)))&
              + exp(-(4167_pr)/Tloc)
    k(r62b) = getlindratecoeff(Tloc,k_0(FOr62b),k_inf(FOr62b),FC(FOr62b), M(mM62),Ploc)
    k(r63b) = (140209257.829813_pr)*exp((-288098.105796075_pr)*R_T_inv + T_log *(-0.000246601911888883_pr))
    k_0(FOr64b) = (7.2734990847763e+35_pr)*exp((-179650.758418867_pr)*R_T_inv + T_log *(-7.69743946117778_pr))
    k_inf(FOr64b) = (654614792187.994_pr)*exp((-158103.158278364_pr)*R_T_inv + T_log *(0.376560563683083_pr))
    FC(FOr64b) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r64b) = getlindratecoeff(Tloc,k_0(FOr64b),k_inf(FOr64b),FC(FOr64b), M(mM64),Ploc)
    k(r65b) = (0.00208599420281629_pr)*exp((-19549.9789483594_pr)*R_T_inv + T_log *(2.86159150202553_pr))
    k_0(FOr66b) = (9.30256861130987e+37_pr)*exp((-448009.761866563_pr)*R_T_inv + T_log *(-7.18291913369542_pr))
    k_inf(FOr66b) = (2.43549596567901e+20_pr)*exp((-426650.441674565_pr)*R_T_inv + T_log *(-1.09291909951663_pr))
    FC(FOr66b) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r66b) = getlindratecoeff(Tloc,k_0(FOr66b),k_inf(FOr66b),FC(FOr66b), M(mM66),Ploc)
    k(r67b) = (478777.187797747_pr)*exp((-282571.378137728_pr)*R_T_inv + T_log *(0.394298722048155_pr))
    k(r68b) = (0.0713910027736444_pr)*exp((-44525.4544468034_pr)*R_T_inv + T_log *(2.31977843779398_pr))
    k(r69b) = (182.301398562242_pr)*exp((-72042.1326168749_pr)*R_T_inv + T_log *(1.49920213246028_pr))
    k(r70b) = (109459.013092112_pr)*exp((-23910.5149520728_pr)*R_T_inv + T_log *(0.474638040475586_pr))
    k(r71b) = (6.56992817698278_pr)*exp((-144057.745834803_pr)*R_T_inv + T_log *(1.48477955493198_pr))
    k_0(FOr72b) = (2.17703277422895e+24_pr)*exp((-352347.772785561_pr)*R_T_inv + T_log *(-3.87261215789313_pr))
    k_inf(FOr72b) = (18463982465.4624_pr)*exp((-332473.772641431_pr)*R_T_inv + T_log *(1.04738786724042_pr))
    FC(FOr72b) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r72b) = getlindratecoeff(Tloc,k_0(FOr72b),k_inf(FOr72b),FC(FOr72b), M(mM72),Ploc)
    k(r73b) = (1422.33089320161_pr)*exp((-78407.930600131_pr)*R_T_inv + T_log *(1.46388219842263_pr))
    k_0(FOr74b) = (1.67927937788963e+17_pr)*exp((-209120.144545601_pr)*R_T_inv + T_log *(-1.73335383879549_pr))
    k_inf(FOr74b) = (5.4028997401273e+18_pr)*exp((-216232.944663931_pr)*R_T_inv + T_log *(-1.2033538597521_pr))
    FC(FOr74b) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r74b) = getlindratecoeff(Tloc,k_0(FOr74b),k_inf(FOr74b),FC(FOr74b), M(mM74),Ploc)
    k(r75b) = (0.439253023654174_pr)*exp((-63252.2854909635_pr)*R_T_inv + T_log *(2.37647616969559_pr))
    k(r76b) = (26056747.432318_pr)*exp((-291077.984221164_pr)*R_T_inv + T_log *(0.255006847177944_pr))
    k(r77b) = (94665.5022652033_pr)*exp((-133516.040103392_pr)*R_T_inv + T_log *(0.468795820501492_pr))
    k(r78b) = (80465685035.8442_pr)*exp((-254780.91217683_pr)*R_T_inv + T_log *(0.468795807406884_pr))
    k(r79b) = (21226392940.7599_pr)*exp((-378616.955968607_pr)*R_T_inv + T_log *(-0.436495087983854_pr))
    k(r80b) = (323839262704.935_pr)*exp((-325516.052154305_pr)*R_T_inv + T_log *(-0.544928141776289_pr))
    k(r81b) = (9.7585591424852_pr)*exp((-88667.5852375846_pr)*R_T_inv + T_log *(2.15093466406397_pr))
    k(r82b) = (427607683391.282_pr)*exp((-363094.779409074_pr)*R_T_inv + T_log *(-0.628713500260053_pr))
    k_0(FOr83b) = (1.22654710600691e+35_pr)*exp((-402207.158039837_pr)*R_T_inv + T_log *(-6.48856108558369_pr))
    k_inf(FOr83b) = (8.5551675388514e+22_pr)*exp((-394634.118164237_pr)*R_T_inv + T_log *(-1.99856110802375_pr))
    FC(FOr83b) = ((1.0_pr - 0.412_pr)*exp(-Tloc/(195_pr))) + (0.412_pr*exp(-Tloc/(5900_pr)))&
              + exp(-(6394_pr)/Tloc)
    k(r83b) = getlindratecoeff(Tloc,k_0(FOr83b),k_inf(FOr83b),FC(FOr83b), M(mM83),Ploc)
    k(r84b) = (2.22817302912156_pr)*exp((-57078.4714305577_pr)*R_T_inv + T_log *(1.95915748485357_pr))
    k(r85b) = (29108649727.9622_pr)*exp((-2751.1922583827_pr)*R_T_inv + T_log *(-0.89705717012382_pr))
    k(r86b) = (1.60141577618367_pr)*exp((-70384.3400666882_pr)*R_T_inv + T_log *(1.90664444894449_pr))
    k(r87b) = (84427796.929894_pr)*exp((-108558.345808898_pr)*R_T_inv + T_log *(0.0889988938610438_pr))
    k(r88b) = (38827514.6367383_pr)*exp((-432529.415772079_pr)*R_T_inv + T_log *(0.417868123691779_pr))
    k(r89b) = (129.436515761268_pr)*exp((-127346.241051639_pr)*R_T_inv + T_log *(1.4393677177393_pr))
    k(r90b) = (806876.175718482_pr)*exp((-374297.142386748_pr)*R_T_inv + T_log *(0.676367343410574_pr))
    k(r91b) = (252773.672006485_pr)*exp((-404286.581019818_pr)*R_T_inv + T_log *(0.500568348716207_pr))
    k(r92b) = (0.358264776705963_pr)*exp((-90152.6101083563_pr)*R_T_inv + T_log *(1.97716458296435_pr))
    k(r93b) = (5.0033085347215_pr)*exp((-69953.7315763695_pr)*R_T_inv + T_log *(2.15296355948785_pr))
    k(r94b) = (4.74490121757084e-07_pr)*exp((-95328.279694001_pr)*R_T_inv + T_log *(3.89717512911257_pr))
    k(r95b) = (6.11223450183379e-13_pr)*exp((-220859.273557687_pr)*R_T_inv + T_log *(4.88195467894354_pr))
    k(r96b) = (153858806.852544_pr)*exp((-352154.832241754_pr)*R_T_inv + T_log *(-0.0463494709264068_pr))
    k(r97b) = (0.0373204256811633_pr)*exp((-42854.6295100157_pr)*R_T_inv + T_log *(2.28547370727389_pr))
    k(r98b) = (0.0144709086748178_pr)*exp((-80716.825058759_pr)*R_T_inv + T_log *(2.49366063414898_pr))
    k(r99b) = (108115.917088564_pr)*exp((-62863.3254811453_pr)*R_T_inv + T_log *(0.428520251638122_pr))
    k(r100b) = (4935530.05506712_pr)*exp((-154620.592110153_pr)*R_T_inv + T_log *(-0.213788971889509_pr))
    k(r101b) = (15945558054.8193_pr)*exp((-211648.512083266_pr)*R_T_inv + T_log *(-0.213788967132282_pr))
    k(r102b) = (179066727.716594_pr)*exp((-474670.10233691_pr)*R_T_inv + T_log *(0.0903714034147012_pr))
    k(r103b) = (112215069.373567_pr)*exp((-235839.729341859_pr)*R_T_inv + T_log *(-0.051638523588781_pr))
    k(r104b) = (266364516.623079_pr)*exp((-104784.713742335_pr)*R_T_inv + T_log *(-0.0510396455618949_pr))
    k(r105b) = (147114358947.814_pr)*exp((-356161.915781535_pr)*R_T_inv + T_log *(-0.503701533850294_pr))
    k(r106b) = (4.46466558222771_pr)*exp((-47695.016882915_pr)*R_T_inv + T_log *(1.79057190877595_pr))
    k(r107b) = (179745115.217238_pr)*exp((-306681.547195364_pr)*R_T_inv + T_log *(-0.0797261897499784_pr))
    k(r108b) = (823499556.846774_pr)*exp((-953.465324075955_pr)*R_T_inv + T_log *(-0.197052458918754_pr))
    k(r109b) = (107060236726.102_pr)*exp((-246241.546649827_pr)*R_T_inv + T_log *(-0.695862757828971_pr))
    k(r110b) = (423292392205.071_pr)*exp((-229764.265791704_pr)*R_T_inv + T_log *(-0.854232144995402_pr))
    k(r111b) = (1307767687909.99_pr)*exp((-254699.896567169_pr)*R_T_inv + T_log *(-0.833061441953099_pr))
    k_0(FOr112b) = (3.70913447723497e+28_pr)*exp((-319480.24573168_pr)*R_T_inv + T_log *(-5.03560992176423_pr))
    k_inf(FOr112b) = (6.89430159622198e+19_pr)*exp((-311380.021687711_pr)*R_T_inv + T_log *(-1.29560991412606_pr))
    FC(FOr112b) = ((1.0_pr - 0.5757_pr)*exp(-Tloc/(237_pr))) + (0.5757_pr*exp(-Tloc/(1652_pr)))&
              + exp(-(5069_pr)/Tloc)
    k(r112b) = getlindratecoeff(Tloc,k_0(FOr112b),k_inf(FOr112b),FC(FOr112b), M(mM112),Ploc)
    k(r113b) = (75793.1792985743_pr)*exp((-336425.218299929_pr)*R_T_inv + T_log *(0.702505993305569_pr))
    k(r114b) = (138762773622580_pr)*exp((-319359.37410918_pr)*R_T_inv + T_log *(-1.31763582467194_pr))
    k(r116b) = (82.7476982454678_pr)*exp((-59905.9390200582_pr)*R_T_inv + T_log *(1.59472473988448_pr))
    k(r117b) = (47015644659177.8_pr)*exp((-273485.221568028_pr)*R_T_inv + T_log *(-0.988771196752537_pr))
    k(r118b) = (0.990098514064229_pr)*exp((-57530.7487796671_pr)*R_T_inv + T_log *(1.94748693869814_pr))
    k_0(FOr119b) = (2.22203794760432e+35_pr)*exp((-362685.761506126_pr)*R_T_inv + T_log *(-6.6831955124041_pr))
    k_inf(FOr119b) = (6.69089370531227e+19_pr)*exp((-351870.121372284_pr)*R_T_inv + T_log *(-1.07319548852198_pr))
    FC(FOr119b) = ((1.0_pr - 0.5907_pr)*exp(-Tloc/(275_pr))) + (0.5907_pr*exp(-Tloc/(1226_pr)))&
              + exp(-(5185_pr)/Tloc)
    k(r119b) = getlindratecoeff(Tloc,k_0(FOr119b),k_inf(FOr119b),FC(FOr119b), M(mM119),Ploc)
    k(r120b) = (13204318.9554577_pr)*exp((-40089.1272645726_pr)*R_T_inv + T_log *(-0.0837853611206428_pr))
    k(r121b) = (1.88323929171178_pr)*exp((-283708.283379444_pr)*R_T_inv + T_log *(0.158073506963395_pr))
    k(r122b) = (1542308.05230658_pr)*exp((-780824.750414087_pr)*R_T_inv + T_log *(0.428814990427002_pr))
    k(r123b) = (10197852531.7205_pr)*exp((-67234.3463222137_pr)*R_T_inv + T_log *(-0.489060627960015_pr))
    k_0(FOr124b) = (1.27539841126381e+38_pr)*exp((-413334.293736671_pr)*R_T_inv + T_log *(-7.37150390785924_pr))
    k_inf(FOr124b) = (3.26990388142892e+23_pr)*exp((-397037.61361243_pr)*R_T_inv + T_log *(-2.17150388580472_pr))
    FC(FOr124b) = ((1.0_pr - 0.6027_pr)*exp(-Tloc/(208_pr))) + (0.6027_pr*exp(-Tloc/(3922_pr)))&
              + exp(-(10180_pr)/Tloc)
    k(r124b) = getlindratecoeff(Tloc,k_0(FOr124b),k_inf(FOr124b),FC(FOr124b), M(mM124),Ploc)
    k(r125b) = (26408645.9643953_pr)*exp((-37578.7274873188_pr)*R_T_inv + T_log *(-0.0837854007120996_pr))
    k(r126b) = (12416189421966.1_pr)*exp((-308679.068718033_pr)*R_T_inv + T_log *(-1.07255653763348_pr))
    k(r127b) = (5668758.1883892_pr)*exp((-58122.9158762821_pr)*R_T_inv + T_log *(-0.136298392590972_pr))
    k(r128b) = (7922591.44736081_pr)*exp((-37578.7272715189_pr)*R_T_inv + T_log *(-0.0837853623328468_pr))
    k(r129b) = (6162015.39075275_pr)*exp((-37578.7272500854_pr)*R_T_inv + T_log *(-0.0837853585567989_pr))
    k(r130b) = (112505.53160049_pr)*exp((-254829.313764051_pr)*R_T_inv + T_log *(0.51028757610441_pr))
    k(r131b) = (3617565.12794363_pr)*exp((-77953.080833721_pr)*R_T_inv + T_log *(-0.0692822014646095_pr))
    k(r132b) = (1717625312.26455_pr)*exp((-11223.5742068129_pr)*R_T_inv + T_log *(-0.329570179682377_pr))
    k(r133b) = (457937.113370923_pr)*exp((-300899.269964827_pr)*R_T_inv + T_log *(0.194521957558832_pr))
    k(r134b) = (0.0724141718784036_pr)*exp((-96072.3317630534_pr)*R_T_inv + T_log *(2.63215141959656_pr))
    k_0(FOr135b) = (4.71543604905195e+42_pr)*exp((-392854.154450694_pr)*R_T_inv + T_log *(-8.11071382359718_pr))
    k_inf(FOr135b) = (9.38926808677002e+23_pr)*exp((-384034.282667138_pr)*R_T_inv + T_log *(-2.26071386198154_pr))
    FC(FOr135b) = ((1.0_pr - 0.619_pr)*exp(-Tloc/(73.2_pr))) + (0.619_pr*exp(-Tloc/(1180_pr)))&
              + exp(-(9999_pr)/Tloc)
    k(r135b) = getlindratecoeff(Tloc,k_0(FOr135b),k_inf(FOr135b),FC(FOr135b), M(mM135),Ploc)
    k(r136b) = (202931386115.706_pr)*exp((-5608.62455897509_pr)*R_T_inv + T_log *(-0.877794684821185_pr))
    k(r137b) = (1284051174.98956_pr)*exp((-375199.153798522_pr)*R_T_inv + T_log *(0.111224013417628_pr))
    k(r138b) = (0.00782341689641042_pr)*exp((-96404.4686440941_pr)*R_T_inv + T_log *(2.76272332885858_pr))
    k(r139b) = (466.078011851944_pr)*exp((-77925.8697756884_pr)*R_T_inv + T_log *(1.17052018098384_pr))
    k(r140b) = (495.920975040689_pr)*exp((-47936.431247477_pr)*R_T_inv + T_log *(1.34631915706947_pr))
    k(r141b) = (0.14694868623737_pr)*exp((-13557.1692150551_pr)*R_T_inv + T_log *(1.97882929857391_pr))
    k(r142b) = (1.56731624327277_pr)*exp((-63469.2846966943_pr)*R_T_inv + T_log *(1.80701623736264_pr))
    k(r143b) = (609562.815452753_pr)*exp((-6540.94879155969_pr)*R_T_inv + T_log *(-0.852873366505309_pr))
    k(r144b) = (75992.1632677637_pr)*exp((-6540.94878161098_pr)*R_T_inv + T_log *(-0.85287336468662_pr))
    k(r145b) = (5812187.58020879_pr)*exp((-141033.031471993_pr)*R_T_inv + T_log *(0.162861291149088_pr))
    k(r146b) = (1616431.37239286_pr)*exp((-84892.7582943105_pr)*R_T_inv + T_log *(0.421360472540468_pr))
    k(r147b) = (1.20407461091903e-20_pr)*exp((-96347.0767916354_pr)*R_T_inv + T_log *(7.84556150201477_pr))
    k(r148b) = (455297057.816222_pr)*exp((-369000.759751352_pr)*R_T_inv + T_log *(-0.777740915683757_pr))
    k_0(FOr149b) = (4.00557021696959e+37_pr)*exp((-232571.384321892_pr)*R_T_inv + T_log *(-9.2854994794244_pr))
    k_inf(FOr149b) = (202813.630183864_pr)*exp((-186421.864134675_pr)*R_T_inv + T_log *(0.454500553565961_pr))
    FC(FOr149b) = ((1.0_pr - 0.7345_pr)*exp(-Tloc/(180_pr))) + (0.7345_pr*exp(-Tloc/(1035_pr)))&
              + exp(-(5417_pr)/Tloc)
    k(r149b) = getlindratecoeff(Tloc,k_0(FOr149b),k_inf(FOr149b),FC(FOr149b), M(mM149),Ploc)
    k(r150b) = (736846.917780792_pr)*exp((-69671.2042485753_pr)*R_T_inv + T_log *(0.0931741209560246_pr))
    k(r151b) = (3.92362278486906e-07_pr)*exp((-359323.551660622_pr)*R_T_inv + T_log *(1.65727570159031_pr))
    k(r153b) = (89.6008362249231_pr)*exp((-59591.9320788026_pr)*R_T_inv + T_log *(1.34689709490345_pr))
    k(r154b) = (8985084870.45986_pr)*exp((-365678.704184587_pr)*R_T_inv + T_log *(0.25500685369819_pr))
    k(r157b) = (147126347.466701_pr)*exp((-259856.643519386_pr)*R_T_inv + T_log *(-0.18815926495398_pr))
    k(r159b) = (18207997.6008078_pr)*exp((-26467.963643811_pr)*R_T_inv + T_log *(-0.145331607623535_pr))
    k(r160b) = (22.8945202979817_pr)*exp((-57378.191997019_pr)*R_T_inv + T_log *(1.30864368596726_pr))
    k_0(FOr161b) = (4.58289120939449e+35_pr)*exp((-159931.839068909_pr)*R_T_inv + T_log *(-7.50218460369188_pr))
    k_inf(FOr161b) = (220313885889.93_pr)*exp((-136463.783053088_pr)*R_T_inv + T_log *(0.549815399099047_pr))
    FC(FOr161b) = ((1.0_pr - 0.465_pr)*exp(-Tloc/(201_pr))) + (0.465_pr*exp(-Tloc/(1773_pr)))&
              + exp(-(5333_pr)/Tloc)
    k(r161b) = getlindratecoeff(Tloc,k_0(FOr161b),k_inf(FOr161b),FC(FOr161b), M(mM161),Ploc)
    k(r165b) = (69.5054141490795_pr)*exp((-50495.3420056011_pr)*R_T_inv + T_log *(1.2098375194734_pr))
    k(r166b) = (7049026.25220497_pr)*exp((-289252.953292496_pr)*R_T_inv + T_log *(0.189043898481464_pr))
  
    return
  end subroutine get_rate_coefficients

  ! ----------------------------------------------- !
  ! Evaluate reaction rates                         !
  ! ----------------------------------------------- !
  subroutine get_reaction_rates(w,k,m,c,cqss)
    implicit none

    real(pr), dimension(nspec) :: c
    real(pr), dimension(nqss) :: cqss
    real(pr), dimension(nreac + nreac_reverse) :: w,k
    real(pr), dimension(nTB + nFO) :: m

    w(r1f) = k(r1f) * c(sH) * c(sO) * m(mM1)
    w(r2f) = k(r2f) * c(sH2) * c(sO) 
    w(r3f) = k(r3f) * c(sHO2) * c(sO) 
    w(r4f) = k(r4f) * c(sH2O2) * c(sO) 
    w(r5f) = k(r5f) * cqss(sqssCH) * c(sO) 
    w(r6f) = k(r6f) * cqss(sqssCH2) * c(sO) 
    w(r7f) = k(r7f) * cqss(sqssCH2GSG) * c(sO) 
    w(r8f) = k(r8f) * cqss(sqssCH2GSG) * c(sO) 
    w(r9f) = k(r9f) * c(sCH3) * c(sO) 
    w(r10f) = k(r10f) * c(sCH4) * c(sO) 
    w(r11f) = k(r11f) * c(sCO) * c(sO) 
    w(r12f) = k(r12f) * cqss(sqssHCO) * c(sO) 
    w(r13f) = k(r13f) * cqss(sqssHCO) * c(sO) 
    w(r14f) = k(r14f) * c(sCH2O) * c(sO) 
    w(r15f) = k(r15f) * cqss(sqssCH2OH) * c(sO) 
    w(r16f) = k(r16f) * c(sCH3O) * c(sO) 
    w(r17f) = k(r17f) * c(sCH3OH) * c(sO) 
    w(r18f) = k(r18f) * c(sCH3OH) * c(sO) 
    w(r19f) = k(r19f) * c(sC2H2) * c(sO) 
    w(r20f) = k(r20f) * c(sC2H2) * c(sO) 
    w(r21f) = k(r21f) * cqss(sqssC2H3) * c(sO) 
    w(r22f) = k(r22f) * c(sC2H4) * c(sO) 
    w(r23f) = k(r23f) * cqss(sqssC2H5) * c(sO) 
    w(r24f) = k(r24f) * c(sC2H6) * c(sO) 
    w(r25f) = k(r25f) * cqss(sqssHCCO) * c(sO) 
    w(r26f) = k(r26f) * c(sCH2CO) * c(sO) 
    w(r27f) = k(r27f) * c(sCH2CO) * c(sO) 
    w(r28f) = k(r28f) * c(sCO) * c(sO2) 
    w(r29f) = k(r29f) * c(sCH2O) * c(sO2) 
    w(r30f) = k(r30f) * c(sH) * c(sO2) * m(mM30)
    w(r31f) = k(r31f) * c(sH) * c(sO2) 
    w(r32f) = k(r32f) * c(sH) * c(sH2O) * c(sO2) 
    w(r33f) = k(r33f) * c(sH) * c(sN2) * c(sO2) 
    w(r34f) = k(r34f) * c(sH) * c(sO2) 
    w(r35f) = k(r35f) * c(sH)**2.0_pr * m(mM35)
    w(r36f) = k(r36f) * c(sH)**2.0_pr * c(sH2) 
    w(r37f) = k(r37f) * c(sH)**2.0_pr * c(sH2O) 
    w(r38f) = k(r38f) * c(sH) * c(sOH) * m(mM38)
    w(r39f) = k(r39f) * c(sH) * c(sHO2) 
    w(r40f) = k(r40f) * c(sH) * c(sHO2) 
    w(r41f) = k(r41f) * c(sH) * c(sHO2) 
    w(r42f) = k(r42f) * c(sH) * c(sH2O2) 
    w(r43f) = k(r43f) * c(sH) * c(sH2O2) 
    w(r44f) = k(r44f) * cqss(sqssCH2) * c(sH) 
    w(r45f) = k(r45f) * cqss(sqssCH2GSG) * c(sH) 
    w(r46f) = k(r46f) * c(sCH3) * c(sH) 
    w(r47f) = k(r47f) * c(sCH4) * c(sH) 
    w(r48f) = k(r48f) * c(sH) * cqss(sqssHCO) 
    w(r49f) = k(r49f) * c(sCH2O) * c(sH) 
    w(r50f) = k(r50f) * c(sCH2O) * c(sH) 
    w(r51f) = k(r51f) * c(sCH2O) * c(sH) 
    w(r52f) = k(r52f) * cqss(sqssCH2OH) * c(sH) 
    w(r53f) = k(r53f) * cqss(sqssCH2OH) * c(sH) 
    w(r54f) = k(r54f) * cqss(sqssCH2OH) * c(sH) 
    w(r55f) = k(r55f) * c(sCH3O) * c(sH) 
    w(r56f) = k(r56f) * c(sCH3O) * c(sH) 
    w(r57f) = k(r57f) * c(sCH3O) * c(sH) 
    w(r58f) = k(r58f) * c(sCH3O) * c(sH) 
    w(r59f) = k(r59f) * c(sCH3O) * c(sH) 
    w(r60f) = k(r60f) * c(sCH3OH) * c(sH) 
    w(r61f) = k(r61f) * c(sCH3OH) * c(sH) 
    w(r62f) = k(r62f) * c(sC2H2) * c(sH) 
    w(r63f) = k(r63f) * cqss(sqssC2H3) * c(sH) 
    w(r64f) = k(r64f) * c(sC2H4) * c(sH) 
    w(r65f) = k(r65f) * c(sC2H4) * c(sH) 
    w(r66f) = k(r66f) * cqss(sqssC2H5) * c(sH) 
    w(r67f) = k(r67f) * cqss(sqssC2H5) * c(sH) 
    w(r68f) = k(r68f) * c(sC2H6) * c(sH) 
    w(r69f) = k(r69f) * c(sH) * cqss(sqssHCCO) 
    w(r70f) = k(r70f) * c(sCH2CO) * c(sH) 
    w(r71f) = k(r71f) * c(sCH2CO) * c(sH) 
    w(r72f) = k(r72f) * c(sCO) * c(sH2) 
    w(r73f) = k(r73f) * c(sH2) * c(sOH) 
    w(r74f) = k(r74f) * c(sOH)**2.0_pr 
    w(r75f) = k(r75f) * c(sOH)**2.0_pr 
    w(r76f) = k(r76f) * c(sHO2) * c(sOH) 
    w(r77f) = k(r77f) * c(sH2O2) * c(sOH) 
    w(r78f) = k(r78f) * c(sH2O2) * c(sOH) 
    w(r79f) = k(r79f) * cqss(sqssCH) * c(sOH) 
    w(r80f) = k(r80f) * cqss(sqssCH2) * c(sOH) 
    w(r81f) = k(r81f) * cqss(sqssCH2) * c(sOH) 
    w(r82f) = k(r82f) * cqss(sqssCH2GSG) * c(sOH) 
    w(r83f) = k(r83f) * c(sCH3) * c(sOH) 
    w(r84f) = k(r84f) * c(sCH3) * c(sOH) 
    w(r85f) = k(r85f) * c(sCH3) * c(sOH) 
    w(r86f) = k(r86f) * c(sCH4) * c(sOH) 
    w(r87f) = k(r87f) * c(sCO) * c(sOH) 
    w(r88f) = k(r88f) * cqss(sqssHCO) * c(sOH) 
    w(r89f) = k(r89f) * c(sCH2O) * c(sOH) 
    w(r90f) = k(r90f) * cqss(sqssCH2OH) * c(sOH) 
    w(r91f) = k(r91f) * c(sCH3O) * c(sOH) 
    w(r92f) = k(r92f) * c(sCH3OH) * c(sOH) 
    w(r93f) = k(r93f) * c(sCH3OH) * c(sOH) 
    w(r94f) = k(r94f) * c(sC2H2) * c(sOH) 
    w(r95f) = k(r95f) * c(sC2H2) * c(sOH) 
    w(r96f) = k(r96f) * cqss(sqssC2H3) * c(sOH) 
    w(r97f) = k(r97f) * c(sC2H4) * c(sOH) 
    w(r98f) = k(r98f) * c(sC2H6) * c(sOH) 
    w(r99f) = k(r99f) * c(sCH2CO) * c(sOH) 
    w(r100f) = k(r100f) * c(sHO2)**2.0_pr 
    w(r101f) = k(r101f) * c(sHO2)**2.0_pr 
    w(r102f) = k(r102f) * cqss(sqssCH2) * c(sHO2) 
    w(r103f) = k(r103f) * c(sCH3) * c(sHO2) 
    w(r104f) = k(r104f) * c(sCH3) * c(sHO2) 
    w(r105f) = k(r105f) * c(sCO) * c(sHO2) 
    w(r106f) = k(r106f) * c(sCH2O) * c(sHO2) 
    w(r107f) = k(r107f) * cqss(sqssCH) * c(sO2) 
    w(r108f) = k(r108f) * cqss(sqssCH) * c(sH2) 
    w(r109f) = k(r109f) * cqss(sqssCH) * c(sH2O) 
    w(r110f) = k(r110f) * cqss(sqssCH) * c(sCH3) 
    w(r111f) = k(r111f) * cqss(sqssCH) * c(sCH4) 
    w(r112f) = k(r112f) * cqss(sqssCH) * c(sCO) 
    w(r113f) = k(r113f) * cqss(sqssCH) * c(sCO2) 
    w(r114f) = k(r114f) * cqss(sqssCH) * c(sCH2O) 
    w(r115f) = k(r115f) * cqss(sqssCH2) * c(sO2) 
    w(r116f) = k(r116f) * cqss(sqssCH2) * c(sH2) 
    w(r117f) = k(r117f) * cqss(sqssCH2) * c(sCH3) 
    w(r118f) = k(r118f) * cqss(sqssCH2) * c(sCH4) 
    w(r119f) = k(r119f) * cqss(sqssCH2) * c(sCO) 
    w(r120f) = k(r120f) * cqss(sqssCH2GSG) * c(sN2) 
    w(r121f) = k(r121f) * cqss(sqssCH2GSG) * c(sO2) 
    w(r122f) = k(r122f) * cqss(sqssCH2GSG) * c(sO2) 
    w(r123f) = k(r123f) * cqss(sqssCH2GSG) * c(sH2) 
    w(r124f) = k(r124f) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r125f) = k(r125f) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r126f) = k(r126f) * cqss(sqssCH2GSG) * c(sCH3) 
    w(r127f) = k(r127f) * cqss(sqssCH2GSG) * c(sCH4) 
    w(r128f) = k(r128f) * cqss(sqssCH2GSG) * c(sCO) 
    w(r129f) = k(r129f) * cqss(sqssCH2GSG) * c(sCO2) 
    w(r130f) = k(r130f) * cqss(sqssCH2GSG) * c(sCO2) 
    w(r131f) = k(r131f) * c(sC2H6) * cqss(sqssCH2GSG) 
    w(r132f) = k(r132f) * c(sCH3) * c(sO2) 
    w(r133f) = k(r133f) * c(sCH3) * c(sO2) 
    w(r134f) = k(r134f) * c(sCH3) * c(sH2O2) 
    w(r135f) = k(r135f) * c(sCH3)**2.0_pr 
    w(r136f) = k(r136f) * c(sCH3)**2.0_pr 
    w(r137f) = k(r137f) * c(sCH3) * cqss(sqssHCO) 
    w(r138f) = k(r138f) * c(sCH2O) * c(sCH3) 
    w(r139f) = k(r139f) * c(sCH3) * c(sCH3OH) 
    w(r140f) = k(r140f) * c(sCH3) * c(sCH3OH) 
    w(r141f) = k(r141f) * c(sC2H4) * c(sCH3) 
    w(r142f) = k(r142f) * c(sC2H6) * c(sCH3) 
    w(r143f) = k(r143f) * c(sH2O) * cqss(sqssHCO) 
    w(r144f) = k(r144f) * cqss(sqssHCO) * m(mM144)
    w(r145f) = k(r145f) * cqss(sqssHCO) * c(sO2) 
    w(r146f) = k(r146f) * cqss(sqssCH2OH) * c(sO2) 
    w(r147f) = k(r147f) * c(sCH3O) * c(sO2) 
    w(r148f) = k(r148f) * cqss(sqssC2H3) * c(sO2) 
    w(r149f) = k(r149f) * c(sC2H4) 
    w(r150f) = k(r150f) * cqss(sqssC2H5) * c(sO2) 
    w(r151f) = k(r151f) * cqss(sqssHCCO) * c(sO2) 
    w(r152f) = k(r152f) * c(sCH3) * c(sO) 
    w(r153f) = k(r153f) * c(sC2H4) * c(sO) 
    w(r154f) = k(r154f) * c(sHO2) * c(sOH) 
    w(r155f) = k(r155f) * c(sCH3) * c(sOH) 
    w(r156f) = k(r156f) * cqss(sqssCH2) * c(sO2) 
    w(r157f) = k(r157f) * cqss(sqssCH2) * c(sO2) 
    w(r158f) = k(r158f) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r159f) = k(r159f) * cqss(sqssC2H3) * c(sO2) 
    w(r160f) = k(r160f) * cqss(sqssC2H3) * c(sO2) 
    w(r161f) = k(r161f) * c(sCH2CO) * c(sH) 
    w(r162f) = k(r162f) * cqss(sqssCH2CHO) * c(sO) 
    w(r163f) = k(r163f) * cqss(sqssCH2CHO) * c(sO2) 
    w(r164f) = k(r164f) * cqss(sqssCH2CHO) * c(sO2) 
    w(r165f) = k(r165f) * cqss(sqssCH2CHO) * c(sH) 
    w(r166f) = k(r166f) * cqss(sqssCH2CHO) * c(sH) 
    w(r1b) = k(r1b) * c(sOH) * m(mM1)
    w(r2b) = k(r2b) * c(sH) * c(sOH) 
    w(r3b) = k(r3b) * c(sO2) * c(sOH) 
    w(r4b) = k(r4b) * c(sHO2) * c(sOH) 
    w(r5b) = k(r5b) * c(sCO) * c(sH) 
    w(r6b) = k(r6b) * c(sH) * cqss(sqssHCO) 
    w(r7b) = k(r7b) * c(sCO) * c(sH2) 
    w(r8b) = k(r8b) * c(sH) * cqss(sqssHCO) 
    w(r9b) = k(r9b) * c(sCH2O) * c(sH) 
    w(r10b) = k(r10b) * c(sCH3) * c(sOH) 
    w(r11b) = k(r11b) * c(sCO2) 
    w(r12b) = k(r12b) * c(sCO) * c(sOH) 
    w(r13b) = k(r13b) * c(sCO2) * c(sH) 
    w(r14b) = k(r14b) * cqss(sqssHCO) * c(sOH) 
    w(r15b) = k(r15b) * c(sCH2O) * c(sOH) 
    w(r16b) = k(r16b) * c(sCH2O) * c(sOH) 
    w(r17b) = k(r17b) * cqss(sqssCH2OH) * c(sOH) 
    w(r18b) = k(r18b) * c(sCH3O) * c(sOH) 
    w(r19b) = k(r19b) * c(sH) * cqss(sqssHCCO) 
    w(r20b) = k(r20b) * cqss(sqssCH2) * c(sCO) 
    w(r21b) = k(r21b) * c(sCH2CO) * c(sH) 
    w(r22b) = k(r22b) * c(sCH3) * cqss(sqssHCO) 
    w(r23b) = k(r23b) * c(sCH2O) * c(sCH3) 
    w(r24b) = k(r24b) * cqss(sqssC2H5) * c(sOH) 
    w(r25b) = k(r25b) * c(sCO) * c(sH) 
    w(r26b) = k(r26b) * cqss(sqssHCCO) * c(sOH) 
    w(r27b) = k(r27b) * cqss(sqssCH2) * c(sCO2) 
    w(r28b) = k(r28b) * c(sCO2) * c(sO) 
    w(r29b) = k(r29b) * cqss(sqssHCO) * c(sHO2) 
    w(r30b) = k(r30b) * c(sHO2) * m(mM30)
    w(r31b) = k(r31b) * c(sHO2) * c(sO2) 
    w(r32b) = k(r32b) * c(sH2O) * c(sHO2) 
    w(r33b) = k(r33b) * c(sHO2) * c(sN2) 
    w(r34b) = k(r34b) * c(sO) * c(sOH) 
    w(r35b) = k(r35b) * c(sH2) * m(mM35)
    w(r36b) = k(r36b) * c(sH2)**2.0_pr 
    w(r37b) = k(r37b) * c(sH2) * c(sH2O) 
    w(r38b) = k(r38b) * c(sH2O) * m(mM38)
    w(r39b) = k(r39b) * c(sH2O) * c(sO) 
    w(r40b) = k(r40b) * c(sH2) * c(sO2) 
    w(r41b) = k(r41b) * c(sOH)**2.0_pr 
    w(r42b) = k(r42b) * c(sH2) * c(sHO2) 
    w(r43b) = k(r43b) * c(sH2O) * c(sOH) 
    w(r44b) = k(r44b) * c(sCH3) 
    w(r45b) = k(r45b) * cqss(sqssCH) * c(sH2) 
    w(r46b) = k(r46b) * c(sCH4) 
    w(r47b) = k(r47b) * c(sCH3) * c(sH2) 
    w(r48b) = k(r48b) * c(sCO) * c(sH2) 
    w(r49b) = k(r49b) * cqss(sqssCH2OH) 
    w(r50b) = k(r50b) * c(sCH3O) 
    w(r51b) = k(r51b) * c(sH2) * cqss(sqssHCO) 
    w(r52b) = k(r52b) * c(sCH2O) * c(sH2) 
    w(r53b) = k(r53b) * c(sCH3) * c(sOH) 
    w(r54b) = k(r54b) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r55b) = k(r55b) * c(sCH3OH) 
    w(r56b) = k(r56b) * cqss(sqssCH2OH) * c(sH) 
    w(r57b) = k(r57b) * c(sCH2O) * c(sH2) 
    w(r58b) = k(r58b) * c(sCH3) * c(sOH) 
    w(r59b) = k(r59b) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r60b) = k(r60b) * cqss(sqssCH2OH) * c(sH2) 
    w(r61b) = k(r61b) * c(sCH3O) * c(sH2) 
    w(r62b) = k(r62b) * cqss(sqssC2H3) 
    w(r63b) = k(r63b) * c(sC2H2) * c(sH2) 
    w(r64b) = k(r64b) * cqss(sqssC2H5) 
    w(r65b) = k(r65b) * cqss(sqssC2H3) * c(sH2) 
    w(r66b) = k(r66b) * c(sC2H6) 
    w(r67b) = k(r67b) * c(sC2H4) * c(sH2) 
    w(r68b) = k(r68b) * cqss(sqssC2H5) * c(sH2) 
    w(r69b) = k(r69b) * cqss(sqssCH2GSG) * c(sCO) 
    w(r70b) = k(r70b) * c(sH2) * cqss(sqssHCCO) 
    w(r71b) = k(r71b) * c(sCH3) * c(sCO) 
    w(r72b) = k(r72b) * c(sCH2O) 
    w(r73b) = k(r73b) * c(sH) * c(sH2O) 
    w(r74b) = k(r74b) * c(sH2O2) 
    w(r75b) = k(r75b) * c(sH2O) * c(sO) 
    w(r76b) = k(r76b) * c(sH2O) * c(sO2) 
    w(r77b) = k(r77b) * c(sH2O) * c(sHO2) 
    w(r78b) = k(r78b) * c(sH2O) * c(sHO2) 
    w(r79b) = k(r79b) * c(sH) * cqss(sqssHCO) 
    w(r80b) = k(r80b) * c(sCH2O) * c(sH) 
    w(r81b) = k(r81b) * cqss(sqssCH) * c(sH2O) 
    w(r82b) = k(r82b) * c(sCH2O) * c(sH) 
    w(r83b) = k(r83b) * c(sCH3OH) 
    w(r84b) = k(r84b) * cqss(sqssCH2) * c(sH2O) 
    w(r85b) = k(r85b) * cqss(sqssCH2GSG) * c(sH2O) 
    w(r86b) = k(r86b) * c(sCH3) * c(sH2O) 
    w(r87b) = k(r87b) * c(sCO2) * c(sH) 
    w(r88b) = k(r88b) * c(sCO) * c(sH2O) 
    w(r89b) = k(r89b) * c(sH2O) * cqss(sqssHCO) 
    w(r90b) = k(r90b) * c(sCH2O) * c(sH2O) 
    w(r91b) = k(r91b) * c(sCH2O) * c(sH2O) 
    w(r92b) = k(r92b) * cqss(sqssCH2OH) * c(sH2O) 
    w(r93b) = k(r93b) * c(sCH3O) * c(sH2O) 
    w(r94b) = k(r94b) * c(sCH2CO) * c(sH) 
    w(r95b) = k(r95b) * c(sCH3) * c(sCO) 
    w(r96b) = k(r96b) * c(sC2H2) * c(sH2O) 
    w(r97b) = k(r97b) * cqss(sqssC2H3) * c(sH2O) 
    w(r98b) = k(r98b) * cqss(sqssC2H5) * c(sH2O) 
    w(r99b) = k(r99b) * c(sH2O) * cqss(sqssHCCO) 
    w(r100b) = k(r100b) * c(sH2O2) * c(sO2) 
    w(r101b) = k(r101b) * c(sH2O2) * c(sO2) 
    w(r102b) = k(r102b) * c(sCH2O) * c(sOH) 
    w(r103b) = k(r103b) * c(sCH4) * c(sO2) 
    w(r104b) = k(r104b) * c(sCH3O) * c(sOH) 
    w(r105b) = k(r105b) * c(sCO2) * c(sOH) 
    w(r106b) = k(r106b) * c(sH2O2) * cqss(sqssHCO) 
    w(r107b) = k(r107b) * cqss(sqssHCO) * c(sO) 
    w(r108b) = k(r108b) * cqss(sqssCH2) * c(sH) 
    w(r109b) = k(r109b) * c(sCH2O) * c(sH) 
    w(r110b) = k(r110b) * cqss(sqssC2H3) * c(sH) 
    w(r111b) = k(r111b) * c(sC2H4) * c(sH) 
    w(r112b) = k(r112b) * cqss(sqssHCCO) 
    w(r113b) = k(r113b) * c(sCO) * cqss(sqssHCO) 
    w(r114b) = k(r114b) * c(sCH2CO) * c(sH) 
    w(r116b) = k(r116b) * c(sCH3) * c(sH) 
    w(r117b) = k(r117b) * c(sC2H4) * c(sH) 
    w(r118b) = k(r118b) * c(sCH3)**2.0_pr 
    w(r119b) = k(r119b) * c(sCH2CO) 
    w(r120b) = k(r120b) * cqss(sqssCH2) * c(sN2) 
    w(r121b) = k(r121b) * c(sCO) * c(sH) * c(sOH) 
    w(r122b) = k(r122b) * c(sCO) * c(sH2O) 
    w(r123b) = k(r123b) * c(sCH3) * c(sH) 
    w(r124b) = k(r124b) * c(sCH3OH) 
    w(r125b) = k(r125b) * cqss(sqssCH2) * c(sH2O) 
    w(r126b) = k(r126b) * c(sC2H4) * c(sH) 
    w(r127b) = k(r127b) * c(sCH3)**2.0_pr 
    w(r128b) = k(r128b) * cqss(sqssCH2) * c(sCO) 
    w(r129b) = k(r129b) * cqss(sqssCH2) * c(sCO2) 
    w(r130b) = k(r130b) * c(sCH2O) * c(sCO) 
    w(r131b) = k(r131b) * cqss(sqssC2H5) * c(sCH3) 
    w(r132b) = k(r132b) * c(sCH3O) * c(sO) 
    w(r133b) = k(r133b) * c(sCH2O) * c(sOH) 
    w(r134b) = k(r134b) * c(sCH4) * c(sHO2) 
    w(r135b) = k(r135b) * c(sC2H6) 
    w(r136b) = k(r136b) * cqss(sqssC2H5) * c(sH) 
    w(r137b) = k(r137b) * c(sCH4) * c(sCO) 
    w(r138b) = k(r138b) * c(sCH4) * cqss(sqssHCO) 
    w(r139b) = k(r139b) * cqss(sqssCH2OH) * c(sCH4) 
    w(r140b) = k(r140b) * c(sCH3O) * c(sCH4) 
    w(r141b) = k(r141b) * cqss(sqssC2H3) * c(sCH4) 
    w(r142b) = k(r142b) * cqss(sqssC2H5) * c(sCH4) 
    w(r143b) = k(r143b) * c(sCO) * c(sH) * c(sH2O) 
    w(r144b) = k(r144b) * c(sCO) * c(sH) * m(mM144)
    w(r145b) = k(r145b) * c(sCO) * c(sHO2) 
    w(r146b) = k(r146b) * c(sCH2O) * c(sHO2) 
    w(r147b) = k(r147b) * c(sCH2O) * c(sHO2) 
    w(r148b) = k(r148b) * c(sCH2O) * cqss(sqssHCO) 
    w(r149b) = k(r149b) * c(sC2H2) * c(sH2) 
    w(r150b) = k(r150b) * c(sC2H4) * c(sHO2) 
    w(r151b) = k(r151b) * c(sCO) * c(sOH) 
    w(r153b) = k(r153b) * cqss(sqssCH2CHO) * c(sH) 
    w(r154b) = k(r154b) * c(sH2O) * c(sO2) 
    w(r157b) = k(r157b) * c(sCH2O) * c(sO) 
    w(r159b) = k(r159b) * cqss(sqssCH2CHO) * c(sO) 
    w(r160b) = k(r160b) * c(sC2H2) * c(sHO2) 
    w(r161b) = k(r161b) * cqss(sqssCH2CHO) 
    w(r165b) = k(r165b) * c(sCH3) * cqss(sqssHCO) 
    w(r166b) = k(r166b) * c(sCH2CO) * c(sH2) 
  
    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(sH2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         + w(r7f) &
         - w(r7b) &
         + w(r35f) &
         - w(r35b) &
         + w(r36f) &
         - w(r36b) &
         + w(r37f) &
         - w(r37b) &
         + w(r40f) &
         - w(r40b) &
         + w(r42f) &
         - w(r42b) &
         + w(r45f) &
         - w(r45b) &
         + w(r47f) &
         - w(r47b) &
         + w(r48f) &
         - w(r48b) &
         + w(r51f) &
         - w(r51b) &
         + w(r52f) &
         - w(r52b) &
         + w(r57f) &
         - w(r57b) &
         + w(r60f) &
         - w(r60b) &
         + w(r61f) &
         - w(r61b) &
         + w(r63f) &
         - w(r63b) &
         + w(r65f) &
         - w(r65b) &
         + w(r67f) &
         - w(r67b) &
         + w(r68f) &
         - w(r68b) &
         + w(r70f) &
         - w(r70b) &
         - w(r72f) &
         + w(r72b) &
         - w(r73f) &
         + w(r73b) &
         - w(r108f) &
         + w(r108b) &
         - w(r116f) &
         + w(r116b) &
         - w(r123f) &
         + w(r123b) &
         + w(r149f) &
         - w(r149b) &
         + w(r152f) &
         + w(r155f) &
         + w(r158f) &
         + w(r166f) &
         - w(r166b) 

    cdot(sH) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r5f) &
         - w(r5b) &
         + w(r6f) &
         - w(r6b) &
         + w(r8f) &
         - w(r8b) &
         + w(r9f) &
         - w(r9b) &
         + w(r13f) &
         - w(r13b) &
         + w(r19f) &
         - w(r19b) &
         + w(r21f) &
         - w(r21b) &
         + w(r25f) &
         - w(r25b) &
         - w(r30f) &
         + w(r30b) &
         - w(r31f) &
         + w(r31b) &
         - w(r32f) &
         + w(r32b) &
         - w(r33f) &
         + w(r33b) &
         - w(r34f) &
         + w(r34b) &
         -2.0_pr * w(r35f) &
         + 2.0_pr * w(r35b) &
         -2.0_pr * w(r36f) &
         + 2.0_pr * w(r36b) &
         -2.0_pr * w(r37f) &
         + 2.0_pr * w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r39f) &
         + w(r39b) &
         - 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(r47f) &
         + w(r47b) &
         - w(r48f) &
         + w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         - w(r54f) &
         + w(r54b) &
         - w(r55f) &
         + w(r55b) &
         - w(r57f) &
         + w(r57b) &
         - w(r58f) &
         + w(r58b) &
         - w(r59f) &
         + w(r59b) &
         - w(r60f) &
         + w(r60b) &
         - w(r61f) &
         + w(r61b) &
         - w(r62f) &
         + w(r62b) &
         - w(r63f) &
         + w(r63b) &
         - w(r64f) &
         + w(r64b) &
         - 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(r73f) &
         - w(r73b) &
         + w(r79f) &
         - w(r79b) &
         + w(r80f) &
         - w(r80b) &
         + w(r82f) &
         - w(r82b) &
         + w(r87f) &
         - w(r87b) &
         + w(r94f) &
         - w(r94b) &
         + w(r108f) &
         - w(r108b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b) &
         + w(r111f) &
         - w(r111b) &
         + w(r114f) &
         - w(r114b) &
         + w(r115f) &
         + w(r116f) &
         - w(r116b) &
         + w(r117f) &
         - w(r117b) &
         + w(r121f) &
         - w(r121b) &
         + w(r123f) &
         - w(r123b) &
         + w(r126f) &
         - w(r126b) &
         + w(r136f) &
         - w(r136b) &
         + w(r143f) &
         - w(r143b) &
         + w(r144f) &
         - w(r144b) &
         + w(r152f) &
         + w(r153f) &
         - w(r153b) &
         + 2.0_pr * w(r156f) &
         - w(r161f) &
         + w(r161b) &
         + w(r162f) &
         - w(r165f) &
         + w(r165b) &
         - w(r166f) &
         + w(r166b) 

    cdot(sO) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r3b) &
         - w(r4f) &
         + w(r4b) &
         - w(r5f) &
         + w(r5b) &
         - w(r6f) &
         + w(r6b) &
         - 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) &
         - w(r14f) &
         + w(r14b) &
         - w(r15f) &
         + w(r15b) &
         - w(r16f) &
         + w(r16b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) &
         - w(r25f) &
         + w(r25b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         + w(r28f) &
         - w(r28b) &
         + w(r34f) &
         - w(r34b) &
         + w(r39f) &
         - w(r39b) &
         + w(r75f) &
         - w(r75b) &
         + w(r107f) &
         - w(r107b) &
         + w(r132f) &
         - w(r132b) &
         - w(r152f) &
         - w(r153f) &
         + w(r153b) &
         + w(r157f) &
         - w(r157b) &
         + w(r159f) &
         - w(r159b) &
         - w(r162f) 

    cdot(sO2) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         - w(r28f) &
         + w(r28b) &
         - w(r29f) &
         + w(r29b) &
         - w(r30f) &
         + w(r30b) &
         - w(r31f) &
         + w(r31b) &
         - w(r32f) &
         + w(r32b) &
         - w(r33f) &
         + w(r33b) &
         - w(r34f) &
         + w(r34b) &
         + w(r40f) &
         - w(r40b) &
         + w(r76f) &
         - w(r76b) &
         + w(r100f) &
         - w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r103f) &
         - w(r103b) &
         - w(r107f) &
         + w(r107b) &
         - w(r115f) &
         - w(r121f) &
         + w(r121b) &
         - w(r122f) &
         + w(r122b) &
         - w(r132f) &
         + w(r132b) &
         - w(r133f) &
         + w(r133b) &
         - w(r145f) &
         + w(r145b) &
         - w(r146f) &
         + w(r146b) &
         - w(r147f) &
         + w(r147b) &
         - w(r148f) &
         + w(r148b) &
         - w(r150f) &
         + w(r150b) &
         - w(r151f) &
         + w(r151b) &
         + w(r154f) &
         - w(r154b) &
         - w(r156f) &
         - w(r157f) &
         + w(r157b) &
         - w(r159f) &
         + w(r159b) &
         - w(r160f) &
         + w(r160b) &
         - w(r163f) &
         - w(r164f) 

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r3f) &
         - w(r3b) &
         + w(r4f) &
         - w(r4b) &
         + w(r10f) &
         - w(r10b) &
         + w(r12f) &
         - w(r12b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r16f) &
         - w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         + w(r24f) &
         - w(r24b) &
         + w(r26f) &
         - w(r26b) &
         + w(r34f) &
         - w(r34b) &
         - w(r38f) &
         + w(r38b) &
         + 2.0_pr * w(r41f) &
         - 2.0_pr * w(r41b) &
         + w(r43f) &
         - w(r43b) &
         + w(r53f) &
         - w(r53b) &
         + w(r58f) &
         - w(r58b) &
         - w(r73f) &
         + w(r73b) &
         -2.0_pr * w(r74f) &
         + 2.0_pr * w(r74b) &
         -2.0_pr * w(r75f) &
         + 2.0_pr * w(r75b) &
         - 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(r82f) &
         + w(r82b) &
         - w(r83f) &
         + w(r83b) &
         - w(r84f) &
         + w(r84b) &
         - w(r85f) &
         + w(r85b) &
         - w(r86f) &
         + w(r86b) &
         - w(r87f) &
         + w(r87b) &
         - w(r88f) &
         + w(r88b) &
         - w(r89f) &
         + w(r89b) &
         - w(r90f) &
         + w(r90b) &
         - w(r91f) &
         + w(r91b) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         - w(r96f) &
         + w(r96b) &
         - w(r97f) &
         + w(r97b) &
         - w(r98f) &
         + w(r98b) &
         - w(r99f) &
         + w(r99b) &
         + w(r102f) &
         - w(r102b) &
         + w(r104f) &
         - w(r104b) &
         + w(r105f) &
         - w(r105b) &
         + w(r115f) &
         + w(r121f) &
         - w(r121b) &
         + w(r133f) &
         - w(r133b) &
         + w(r151f) &
         - w(r151b) &
         - w(r154f) &
         + w(r154b) &
         - w(r155f) &
         + w(r163f) &
         + w(r164f) 

    cdot(sH2O) = 0.0_pr &
         + w(r38f) &
         - w(r38b) &
         + w(r39f) &
         - w(r39b) &
         + w(r43f) &
         - w(r43b) &
         + w(r54f) &
         - w(r54b) &
         + w(r59f) &
         - w(r59b) &
         + w(r73f) &
         - w(r73b) &
         + w(r75f) &
         - w(r75b) &
         + w(r76f) &
         - w(r76b) &
         + w(r77f) &
         - w(r77b) &
         + w(r78f) &
         - w(r78b) &
         + w(r81f) &
         - w(r81b) &
         + w(r84f) &
         - w(r84b) &
         + w(r85f) &
         - w(r85b) &
         + w(r86f) &
         - w(r86b) &
         + w(r88f) &
         - w(r88b) &
         + w(r89f) &
         - w(r89b) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         - w(r91b) &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         + w(r96f) &
         - w(r96b) &
         + w(r97f) &
         - w(r97b) &
         + w(r98f) &
         - w(r98b) &
         + w(r99f) &
         - w(r99b) &
         - w(r109f) &
         + w(r109b) &
         + w(r122f) &
         - w(r122b) &
         - w(r124f) &
         + w(r124b) &
         + w(r154f) &
         - w(r154b) &
         - w(r158f) 

    cdot(sHO2) = 0.0_pr &
         - w(r3f) &
         + w(r3b) &
         + w(r4f) &
         - w(r4b) &
         + w(r29f) &
         - w(r29b) &
         + w(r30f) &
         - w(r30b) &
         + w(r31f) &
         - w(r31b) &
         + w(r32f) &
         - w(r32b) &
         + w(r33f) &
         - w(r33b) &
         - w(r39f) &
         + w(r39b) &
         - w(r40f) &
         + w(r40b) &
         - w(r41f) &
         + w(r41b) &
         + w(r42f) &
         - w(r42b) &
         - w(r76f) &
         + w(r76b) &
         + w(r77f) &
         - w(r77b) &
         + w(r78f) &
         - w(r78b) &
         -2.0_pr * w(r100f) &
         + 2.0_pr * w(r100b) &
         -2.0_pr * w(r101f) &
         + 2.0_pr * w(r101b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         - w(r105f) &
         + w(r105b) &
         - w(r106f) &
         + w(r106b) &
         + w(r134f) &
         - w(r134b) &
         + w(r145f) &
         - w(r145b) &
         + w(r146f) &
         - w(r146b) &
         + w(r147f) &
         - w(r147b) &
         + w(r150f) &
         - w(r150b) &
         - w(r154f) &
         + w(r154b) &
         + w(r160f) &
         - w(r160b) 

    cdot(sH2O2) = 0.0_pr &
         - w(r4f) &
         + w(r4b) &
         - w(r42f) &
         + w(r42b) &
         - w(r43f) &
         + w(r43b) &
         + w(r74f) &
         - w(r74b) &
         - w(r77f) &
         + w(r77b) &
         - w(r78f) &
         + w(r78b) &
         + w(r100f) &
         - w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r106f) &
         - w(r106b) &
         - w(r134f) &
         + w(r134b) 

    cdot(sCH3) = 0.0_pr &
         - w(r9f) &
         + w(r9b) &
         + w(r10f) &
         - w(r10b) &
         + w(r22f) &
         - w(r22b) &
         + w(r23f) &
         - w(r23b) &
         + w(r44f) &
         - w(r44b) &
         - w(r46f) &
         + w(r46b) &
         + w(r47f) &
         - w(r47b) &
         + w(r53f) &
         - w(r53b) &
         + w(r58f) &
         - w(r58b) &
         + w(r71f) &
         - w(r71b) &
         - w(r83f) &
         + w(r83b) &
         - w(r84f) &
         + w(r84b) &
         - w(r85f) &
         + w(r85b) &
         + w(r86f) &
         - w(r86b) &
         + w(r95f) &
         - w(r95b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         - w(r110f) &
         + w(r110b) &
         + w(r116f) &
         - w(r116b) &
         - w(r117f) &
         + w(r117b) &
         + 2.0_pr * w(r118f) &
         - 2.0_pr * w(r118b) &
         + w(r123f) &
         - w(r123b) &
         - w(r126f) &
         + w(r126b) &
         + 2.0_pr * w(r127f) &
         - 2.0_pr * w(r127b) &
         + w(r131f) &
         - w(r131b) &
         - w(r132f) &
         + w(r132b) &
         - w(r133f) &
         + w(r133b) &
         - w(r134f) &
         + w(r134b) &
         -2.0_pr * w(r135f) &
         + 2.0_pr * w(r135b) &
         -2.0_pr * w(r136f) &
         + 2.0_pr * 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(r152f) &
         - w(r155f) &
         + w(r165f) &
         - w(r165b) 

    cdot(sCH4) = 0.0_pr &
         - w(r10f) &
         + w(r10b) &
         + w(r46f) &
         - w(r46b) &
         - w(r47f) &
         + w(r47b) &
         - w(r86f) &
         + w(r86b) &
         + w(r103f) &
         - w(r103b) &
         - w(r111f) &
         + w(r111b) &
         - w(r118f) &
         + w(r118b) &
         - w(r127f) &
         + w(r127b) &
         + w(r134f) &
         - w(r134b) &
         + w(r137f) &
         - w(r137b) &
         + w(r138f) &
         - w(r138b) &
         + w(r139f) &
         - w(r139b) &
         + w(r140f) &
         - w(r140b) &
         + w(r141f) &
         - w(r141b) &
         + w(r142f) &
         - w(r142b) 

    cdot(sCO) = 0.0_pr &
         + w(r5f) &
         - w(r5b) &
         + w(r7f) &
         - w(r7b) &
         - w(r11f) &
         + w(r11b) &
         + w(r12f) &
         - w(r12b) &
         + w(r20f) &
         - w(r20b) &
         + 2.0_pr * w(r25f) &
         - 2.0_pr * w(r25b) &
         - w(r28f) &
         + w(r28b) &
         + w(r48f) &
         - w(r48b) &
         + w(r69f) &
         - w(r69b) &
         + w(r71f) &
         - w(r71b) &
         - w(r72f) &
         + w(r72b) &
         - w(r87f) &
         + w(r87b) &
         + w(r88f) &
         - w(r88b) &
         + w(r95f) &
         - w(r95b) &
         - w(r105f) &
         + w(r105b) &
         - w(r112f) &
         + w(r112b) &
         + w(r113f) &
         - w(r113b) &
         + w(r115f) &
         - w(r119f) &
         + w(r119b) &
         + w(r121f) &
         - w(r121b) &
         + w(r122f) &
         - w(r122b) &
         + w(r130f) &
         - w(r130b) &
         + w(r137f) &
         - w(r137b) &
         + w(r143f) &
         - w(r143b) &
         + w(r144f) &
         - w(r144b) &
         + w(r145f) &
         - w(r145b) &
         + 2.0_pr * w(r151f) &
         - 2.0_pr * w(r151b) &
         + w(r152f) &
         + w(r163f) 

    cdot(sCO2) = 0.0_pr &
         + w(r11f) &
         - w(r11b) &
         + w(r13f) &
         - w(r13b) &
         + w(r27f) &
         - w(r27b) &
         + w(r28f) &
         - w(r28b) &
         + w(r87f) &
         - w(r87b) &
         + w(r105f) &
         - w(r105b) &
         - w(r113f) &
         + w(r113b) &
         - w(r130f) &
         + w(r130b) &
         + w(r156f) &
         + w(r162f) 

    cdot(sCH2O) = 0.0_pr &
         + w(r9f) &
         - w(r9b) &
         - w(r14f) &
         + w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r16f) &
         - w(r16b) &
         + w(r23f) &
         - w(r23b) &
         - w(r29f) &
         + w(r29b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         + w(r52f) &
         - w(r52b) &
         + w(r57f) &
         - w(r57b) &
         + w(r72f) &
         - w(r72b) &
         + w(r80f) &
         - w(r80b) &
         + w(r82f) &
         - w(r82b) &
         - w(r89f) &
         + w(r89b) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         - w(r91b) &
         + w(r102f) &
         - w(r102b) &
         - w(r106f) &
         + w(r106b) &
         + w(r109f) &
         - w(r109b) &
         - w(r114f) &
         + w(r114b) &
         + w(r130f) &
         - w(r130b) &
         + w(r133f) &
         - w(r133b) &
         - w(r138f) &
         + w(r138b) &
         + w(r146f) &
         - w(r146b) &
         + w(r147f) &
         - w(r147b) &
         + w(r148f) &
         - w(r148b) &
         + w(r155f) &
         + w(r157f) &
         - w(r157b) &
         + w(r158f) &
         + w(r163f) 

    cdot(sCH3O) = 0.0_pr &
         - w(r16f) &
         + w(r16b) &
         + w(r18f) &
         - w(r18b) &
         + w(r50f) &
         - w(r50b) &
         - w(r55f) &
         + w(r55b) &
         - w(r56f) &
         + w(r56b) &
         - w(r57f) &
         + w(r57b) &
         - w(r58f) &
         + w(r58b) &
         - w(r59f) &
         + w(r59b) &
         + w(r61f) &
         - w(r61b) &
         - w(r91f) &
         + w(r91b) &
         + w(r93f) &
         - w(r93b) &
         + w(r104f) &
         - w(r104b) &
         + w(r132f) &
         - w(r132b) &
         + w(r140f) &
         - w(r140b) &
         - w(r147f) &
         + w(r147b) 

    cdot(sCH3OH) = 0.0_pr &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         + w(r55f) &
         - w(r55b) &
         - w(r60f) &
         + w(r60b) &
         - w(r61f) &
         + w(r61b) &
         + w(r83f) &
         - w(r83b) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         + w(r124f) &
         - w(r124b) &
         - w(r139f) &
         + w(r139b) &
         - w(r140f) &
         + w(r140b) 

    cdot(sC2H2) = 0.0_pr &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r62f) &
         + w(r62b) &
         + w(r63f) &
         - w(r63b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         + w(r96f) &
         - w(r96b) &
         + w(r149f) &
         - w(r149b) &
         + w(r160f) &
         - w(r160b) 

    cdot(sC2H4) = 0.0_pr &
         - w(r22f) &
         + w(r22b) &
         - w(r64f) &
         + w(r64b) &
         - w(r65f) &
         + w(r65b) &
         + w(r67f) &
         - w(r67b) &
         - w(r97f) &
         + w(r97b) &
         + w(r111f) &
         - w(r111b) &
         + w(r117f) &
         - w(r117b) &
         + w(r126f) &
         - w(r126b) &
         - w(r141f) &
         + w(r141b) &
         - w(r149f) &
         + w(r149b) &
         + w(r150f) &
         - w(r150b) &
         - w(r153f) &
         + w(r153b) 

    cdot(sC2H6) = 0.0_pr &
         - w(r24f) &
         + w(r24b) &
         + w(r66f) &
         - w(r66b) &
         - w(r68f) &
         + w(r68b) &
         - w(r98f) &
         + w(r98b) &
         - w(r131f) &
         + w(r131b) &
         + w(r135f) &
         - w(r135b) &
         - w(r142f) &
         + w(r142b) 

    cdot(sCH2CO) = 0.0_pr &
         + w(r21f) &
         - w(r21b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         - w(r70f) &
         + w(r70b) &
         - w(r71f) &
         + w(r71b) &
         + w(r94f) &
         - w(r94b) &
         - w(r99f) &
         + w(r99b) &
         + w(r114f) &
         - w(r114b) &
         + w(r119f) &
         - w(r119b) &
         - w(r161f) &
         + w(r161b) &
         + w(r166f) &
         - w(r166b) 

    cdot(sN2) = 0.0_pr 

  
    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) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_CH2
    real(pr) :: CH_CH2GSG
    real(pr) :: CH_HCO
    real(pr) :: CH_C2H3
    real(pr) :: CH_HCCO
    real(pr) :: CH2_ct
    real(pr) :: CH2_num
    real(pr) :: CH2_denom
    real(pr) :: CH2_CH
    real(pr) :: CH2_CH2GSG
    real(pr) :: CH2_HCO
    real(pr) :: CH2_CH2CHO
    real(pr) :: CH2GSG_ct
    real(pr) :: CH2GSG_num
    real(pr) :: CH2GSG_denom
    real(pr) :: CH2GSG_CH
    real(pr) :: CH2GSG_CH2
    real(pr) :: CH2GSG_HCO
    real(pr) :: CH2GSG_CH2OH
    real(pr) :: CH2GSG_C2H5
    real(pr) :: CH2GSG_HCCO
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_CH
    real(pr) :: HCO_CH2
    real(pr) :: HCO_CH2GSG
    real(pr) :: HCO_C2H3
    real(pr) :: HCO_CH2CHO
    real(pr) :: CH2OH_ct
    real(pr) :: CH2OH_num
    real(pr) :: CH2OH_denom
    real(pr) :: CH2OH_CH2GSG
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_CH
    real(pr) :: C2H3_HCO
    real(pr) :: C2H3_CH2CHO
    real(pr) :: C2H5_ct
    real(pr) :: C2H5_num
    real(pr) :: C2H5_denom
    real(pr) :: C2H5_CH2GSG
    real(pr) :: HCCO_ct
    real(pr) :: HCCO_num
    real(pr) :: HCCO_denom
    real(pr) :: HCCO_CH
    real(pr) :: HCCO_CH2GSG
    real(pr) :: CH2CHO_ct
    real(pr) :: CH2CHO_num
    real(pr) :: CH2CHO_denom
    real(pr) :: CH2CHO_HCO
    real(pr) :: CH2CHO_C2H3
  
    C2H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r23f)* c(sO) &
             + k(r66f)* c(sH) &
             + k(r67f)* c(sH) &
             + k(r150f)* c(sO2) &
             + k(r24b)* c(sOH) &
             + k(r64b)&
             + k(r68b)* c(sH2) &
             + k(r98b)* c(sH2O) &
             + k(r131b)* c(sCH3) &
             + k(r136b)* c(sH) &
             + k(r142b)* c(sCH4)   )

    C2H5_num = ( 0.0_pr &
             + k(r24f)* c(sC2H6) * c(sO) &
             + k(r64f)* c(sC2H4) * c(sH) &
             + k(r68f)* c(sC2H6) * c(sH) &
             + k(r98f)* c(sC2H6) * c(sOH) &
             + k(r136f) *c(sCH3)** 2.0_pr&
             + k(r142f)* c(sC2H6) * c(sCH3) &
             + k(r23b)* c(sCH2O) * c(sCH3) &
             + k(r66b)* c(sC2H6) &
             + k(r67b)* c(sC2H4) * c(sH2) &
             + k(r150b)* c(sC2H4) * c(sHO2)   )

    C2H5_ct = C2H5_num / C2H5_denom

    C2H5_CH2GSG = - ( 0.0_pr &
          + k(r131f) * c(sC2H6)  ) / C2H5_denom

    CH2GSG_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r7f)* c(sO) &
             + k(r8f)* c(sO) &
             + k(r45f)* c(sH) &
             + k(r82f)* c(sOH) &
             + k(r120f)* c(sN2) &
             + k(r121f)* c(sO2) &
             + k(r122f)* c(sO2) &
             + k(r123f)* c(sH2) &
             + k(r124f)* c(sH2O) &
             + k(r125f)* c(sH2O) &
             + k(r126f)* c(sCH3) &
             + k(r127f)* c(sCH4) &
             + k(r128f)* c(sCO) &
             + k(r129f)* c(sCO2) &
             + k(r130f)* c(sCO2) &
             + k(r131f)* c(sC2H6) &
             + k(r158f)* c(sH2O) &
             + k(r54b)* c(sH2O) &
             + k(r59b)* c(sH2O) &
             + k(r69b)* c(sCO) &
             + k(r85b)* c(sH2O)   )

    CH2GSG_num = ( 0.0_pr &
             + k(r59f)* c(sCH3O) * c(sH) &
             + k(r85f)* c(sCH3) * c(sOH) &
             + k(r7b)* c(sCO) * c(sH2) &
             + k(r82b)* c(sCH2O) * c(sH) &
             + k(r121b)* c(sCO) * c(sH) * c(sOH) &
             + k(r122b)* c(sCO) * c(sH2O) &
             + k(r123b)* c(sCH3) * c(sH) &
             + k(r124b)* c(sCH3OH) &
             + k(r126b)* c(sC2H4) * c(sH) &
             + k(r127b) *c(sCH3)** 2.0_pr&
             + k(r130b)* c(sCH2O) * c(sCO)   )

    CH2GSG_ct = CH2GSG_num / CH2GSG_denom

    CH2GSG_CH = - ( 0.0_pr &
          + k(r45b) * c(sH2)  ) / CH2GSG_denom

    CH2GSG_CH2 = - ( 0.0_pr &
          + k(r120b) * c(sN2) &
          + k(r125b) * c(sH2O) &
          + k(r128b) * c(sCO) &
          + k(r129b) * c(sCO2)  ) / CH2GSG_denom

    CH2GSG_HCO = - ( 0.0_pr &
          + k(r8b) * c(sH)  ) / CH2GSG_denom

    CH2GSG_CH2OH = - ( 0.0_pr &
          + k(r54f) * c(sH)  ) / CH2GSG_denom

    CH2GSG_C2H5 = - ( 0.0_pr &
          + k(r131b) * c(sCH3)  ) / CH2GSG_denom

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

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r21f)* c(sO) &
             + k(r63f)* c(sH) &
             + k(r96f)* c(sOH) &
             + k(r148f)* c(sO2) &
             + k(r159f)* c(sO2) &
             + k(r160f)* c(sO2) &
             + k(r62b)&
             + k(r65b)* c(sH2) &
             + k(r97b)* c(sH2O) &
             + k(r110b)* c(sH) &
             + k(r141b)* c(sCH4)   )

    C2H3_num = ( 0.0_pr &
             + k(r62f)* c(sC2H2) * c(sH) &
             + k(r65f)* c(sC2H4) * c(sH) &
             + k(r97f)* c(sC2H4) * c(sOH) &
             + k(r141f)* c(sC2H4) * c(sCH3) &
             + k(r21b)* c(sCH2CO) * c(sH) &
             + k(r63b)* c(sC2H2) * c(sH2) &
             + k(r96b)* c(sC2H2) * c(sH2O) &
             + k(r160b)* c(sC2H2) * c(sHO2)   )

    C2H3_ct = C2H3_num / C2H3_denom

    C2H3_CH = - ( 0.0_pr &
          + k(r110f) * c(sCH3)  ) / C2H3_denom

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

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

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r5f)* c(sO) &
             + k(r79f)* c(sOH) &
             + k(r107f)* c(sO2) &
             + k(r108f)* c(sH2) &
             + k(r109f)* c(sH2O) &
             + k(r110f)* c(sCH3) &
             + k(r111f)* c(sCH4) &
             + k(r112f)* c(sCO) &
             + k(r113f)* c(sCO2) &
             + k(r114f)* c(sCH2O) &
             + k(r45b)* c(sH2) &
             + k(r81b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r5b)* c(sCO) * c(sH) &
             + k(r109b)* c(sCH2O) * c(sH) &
             + k(r111b)* c(sC2H4) * c(sH) &
             + k(r114b)* c(sCH2CO) * c(sH)   )

    CH_ct = CH_num / CH_denom

    CH_CH2 = - ( 0.0_pr &
          + k(r81f) * c(sOH) &
          + k(r108b) * c(sH)  ) / CH_denom

    CH_CH2GSG = - ( 0.0_pr &
          + k(r45f) * c(sH)  ) / CH_denom

    CH_HCO = - ( 0.0_pr &
          + k(r79b) * c(sH) &
          + k(r107b) * c(sO) &
          + k(r113b) * c(sCO)  ) / CH_denom

    CH_C2H3 = - ( 0.0_pr &
          + k(r110b) * c(sH)  ) / CH_denom

    CH_HCCO = - ( 0.0_pr &
          + k(r112b)  ) / CH_denom

    CH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r6f)* c(sO) &
             + k(r44f)* c(sH) &
             + k(r80f)* c(sOH) &
             + k(r81f)* c(sOH) &
             + k(r102f)* c(sHO2) &
             + k(r115f)* c(sO2) &
             + k(r116f)* c(sH2) &
             + k(r117f)* c(sCH3) &
             + k(r118f)* c(sCH4) &
             + k(r119f)* c(sCO) &
             + k(r156f)* c(sO2) &
             + k(r157f)* c(sO2) &
             + k(r20b)* c(sCO) &
             + k(r27b)* c(sCO2) &
             + k(r84b)* c(sH2O) &
             + k(r108b)* c(sH) &
             + k(r120b)* c(sN2) &
             + k(r125b)* c(sH2O) &
             + k(r128b)* c(sCO) &
             + k(r129b)* c(sCO2)   )

    CH2_num = ( 0.0_pr &
             + k(r20f)* c(sC2H2) * c(sO) &
             + k(r27f)* c(sCH2CO) * c(sO) &
             + k(r84f)* c(sCH3) * c(sOH) &
             + k(r44b)* c(sCH3) &
             + k(r80b)* c(sCH2O) * c(sH) &
             + k(r102b)* c(sCH2O) * c(sOH) &
             + k(r116b)* c(sCH3) * c(sH) &
             + k(r117b)* c(sC2H4) * c(sH) &
             + k(r118b) *c(sCH3)** 2.0_pr&
             + k(r119b)* c(sCH2CO) &
             + k(r157b)* c(sCH2O) * c(sO)   )

    CH2_ct = CH2_num / CH2_denom

    CH2_CH = - ( 0.0_pr &
          + k(r108f) * c(sH2) &
          + k(r81b) * c(sH2O)  ) / CH2_denom

    CH2_CH2GSG = - ( 0.0_pr &
          + k(r120f) * c(sN2) &
          + k(r125f) * c(sH2O) &
          + k(r128f) * c(sCO) &
          + k(r129f) * c(sCO2)  ) / CH2_denom

    CH2_HCO = - ( 0.0_pr &
          + k(r6b) * c(sH)  ) / CH2_denom

    CH2_CH2CHO = - ( 0.0_pr &
          + k(r162f) * c(sO)  ) / CH2_denom

    CH2OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r15f)* c(sO) &
             + k(r52f)* c(sH) &
             + k(r53f)* c(sH) &
             + k(r54f)* c(sH) &
             + k(r90f)* c(sOH) &
             + k(r146f)* c(sO2) &
             + k(r17b)* c(sOH) &
             + k(r49b)&
             + k(r56b)* c(sH) &
             + k(r60b)* c(sH2) &
             + k(r92b)* c(sH2O) &
             + k(r139b)* c(sCH4)   )

    CH2OH_num = ( 0.0_pr &
             + k(r17f)* c(sCH3OH) * c(sO) &
             + k(r49f)* c(sCH2O) * c(sH) &
             + k(r56f)* c(sCH3O) * c(sH) &
             + k(r60f)* c(sCH3OH) * c(sH) &
             + k(r92f)* c(sCH3OH) * c(sOH) &
             + k(r139f)* c(sCH3) * c(sCH3OH) &
             + k(r15b)* c(sCH2O) * c(sOH) &
             + k(r52b)* c(sCH2O) * c(sH2) &
             + k(r53b)* c(sCH3) * c(sOH) &
             + k(r90b)* c(sCH2O) * c(sH2O) &
             + k(r146b)* c(sCH2O) * c(sHO2)   )

    CH2OH_ct = CH2OH_num / CH2OH_denom

    CH2OH_CH2GSG = - ( 0.0_pr &
          + k(r54b) * c(sH2O)  ) / CH2OH_denom

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r12f)* c(sO) &
             + k(r13f)* c(sO) &
             + k(r48f)* c(sH) &
             + k(r88f)* c(sOH) &
             + k(r137f)* c(sCH3) &
             + k(r143f)* c(sH2O) &
             + k(r144f)* m(mM144) &
             + k(r145f)* c(sO2) &
             + k(r6b)* c(sH) &
             + k(r8b)* c(sH) &
             + k(r14b)* c(sOH) &
             + k(r22b)* c(sCH3) &
             + k(r29b)* c(sHO2) &
             + k(r51b)* c(sH2) &
             + k(r79b)* c(sH) &
             + k(r89b)* c(sH2O) &
             + k(r106b)* c(sH2O2) &
             + k(r107b)* c(sO) &
             + k(r113b)* c(sCO) &
             + k(r138b)* c(sCH4) &
             + k(r148b)* c(sCH2O) &
             + k(r165b)* c(sCH3)   )

    HCO_num = ( 0.0_pr &
             + k(r14f)* c(sCH2O) * c(sO) &
             + k(r22f)* c(sC2H4) * c(sO) &
             + k(r29f)* c(sCH2O) * c(sO2) &
             + k(r51f)* c(sCH2O) * c(sH) &
             + k(r89f)* c(sCH2O) * c(sOH) &
             + k(r106f)* c(sCH2O) * c(sHO2) &
             + k(r138f)* c(sCH2O) * c(sCH3) &
             + k(r12b)* c(sCO) * c(sOH) &
             + k(r13b)* c(sCO2) * c(sH) &
             + k(r48b)* c(sCO) * c(sH2) &
             + k(r88b)* c(sCO) * c(sH2O) &
             + k(r137b)* c(sCH4) * c(sCO) &
             + k(r143b)* c(sCO) * c(sH) * c(sH2O) &
             + k(r144b)* c(sCO) * c(sH) * m(mM144) &
             + k(r145b)* c(sCO) * c(sHO2)   )

    HCO_ct = HCO_num / HCO_denom

    HCO_CH = - ( 0.0_pr &
          + k(r79f) * c(sOH) &
          + k(r107f) * c(sO2) &
          + k(r113f) * c(sCO2)  ) / HCO_denom

    HCO_CH2 = - ( 0.0_pr &
          + k(r6f) * c(sO)  ) / HCO_denom

    HCO_CH2GSG = - ( 0.0_pr &
          + k(r8f) * c(sO)  ) / HCO_denom

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

    HCO_CH2CHO = - ( 0.0_pr &
          + k(r164f)  * 2.0_pr * c(sO2) &
          + k(r165f) * c(sH)  ) / HCO_denom

    HCCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r25f)* c(sO) &
             + k(r69f)* c(sH) &
             + k(r151f)* c(sO2) &
             + k(r19b)* c(sH) &
             + k(r26b)* c(sOH) &
             + k(r70b)* c(sH2) &
             + k(r99b)* c(sH2O) &
             + k(r112b)  )

    HCCO_num = ( 0.0_pr &
             + k(r19f)* c(sC2H2) * c(sO) &
             + k(r26f)* c(sCH2CO) * c(sO) &
             + k(r70f)* c(sCH2CO) * c(sH) &
             + k(r99f)* c(sCH2CO) * c(sOH) &
             + k(r25b) *c(sCO)** 2.0_pr* c(sH) &
             + k(r151b) *c(sCO)** 2.0_pr* c(sOH)   )

    HCCO_ct = HCCO_num / HCCO_denom

    HCCO_CH = - ( 0.0_pr &
          + k(r112f) * c(sCO)  ) / HCCO_denom

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

    CH2CHO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r162f)* c(sO) &
             + k(r163f)* c(sO2) &
             + k(r164f)* c(sO2) &
             + k(r165f)* c(sH) &
             + k(r166f)* c(sH) &
             + k(r153b)* c(sH) &
             + k(r159b)* c(sO) &
             + k(r161b)  )

    CH2CHO_num = ( 0.0_pr &
             + k(r153f)* c(sC2H4) * c(sO) &
             + k(r161f)* c(sCH2CO) * c(sH) &
             + k(r166b)* c(sCH2CO) * c(sH2)   )

    CH2CHO_ct = CH2CHO_num / CH2CHO_denom

    CH2CHO_HCO = - ( 0.0_pr &
          + k(r165b) * c(sCH3)  ) / CH2CHO_denom

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

    cqss(sqssC2H5) = (  (C2H5_ct) - ( ( ( ( ( ( (CH2GSG_ct) - (HCCO_ct) &
                 * (CH2GSG_HCCO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * (CH2GSG_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - (CH2OH_ct) &
                 * (CH2GSG_CH2OH)) - ( ( (CH2_ct) - (CH2CHO_ct) &
                 * (CH2_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) - ( ( ( (CH_ct) &
                 - (HCCO_ct) * (CH_HCCO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( ( (CH2_ct) - (CH2CHO_ct) &
                 * (CH2_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) * ( ( ( (CH2GSG_CH) - (HCCO_CH) &
                 * (CH2GSG_HCCO)) - (HCO_CH) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) - ( ( ( ( (C2H3_ct) - (CH2CHO_ct) &
                 * (C2H3_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( ( (CH2_ct) - (CH2CHO_ct) &
                 * (CH2_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 - ( ( ( (CH_ct) - (HCCO_ct) * (CH_HCCO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( ( (CH2_ct) &
                 - (CH2CHO_ct) * (CH2_CH2CHO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * ( ( (C2H3_CH) &
                 - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 / ( ( ( (1.0_pr) - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) * ( ( - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * (CH2GSG_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH_C2H3) &
                 - ( (HCO_C2H3) - (CH2CHO_C2H3) * (HCO_CH2CHO)) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( ( ( (CH2GSG_CH) - (HCCO_CH) * (CH2GSG_HCCO)) &
                 - (HCO_CH) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH2GSG_CH2) &
                 - (HCO_CH2) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) / ( ( ( (1.0_pr) - (CH2CHO_C2H3) &
                 * (C2H3_CH2CHO)) - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH_C2H3) - ( (HCO_C2H3) &
                 - (CH2CHO_C2H3) * (HCO_CH2CHO)) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( ( (C2H3_CH) - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( - (HCO_CH2) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) - (HCCO_CH) &
                 * (CH_HCCO)) - (HCO_CH) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH2_CH) &
                 - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))))) &
                 * (C2H5_CH2GSG) / ( ( ( ( ( (1.0_pr) - (HCCO_CH2GSG) &
                 * (CH2GSG_HCCO)) - (HCO_CH2GSG) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - (CH2OH_CH2GSG) * (CH2GSG_CH2OH)) - ( (CH2_CH2GSG) &
                 - (HCO_CH2GSG) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) - ( ( ( (CH_CH2GSG) &
                 - (HCCO_CH2GSG) * (CH_HCCO)) - (HCO_CH2GSG) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * ( ( ( (CH2GSG_CH) &
                 - (HCCO_CH) * (CH2GSG_HCCO)) - (HCO_CH) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) ) / (  (1.0_pr) - (CH2GSG_C2H5) &
                 * (C2H5_CH2GSG) / ( ( ( ( ( (1.0_pr) - (HCCO_CH2GSG) &
                 * (CH2GSG_HCCO)) - (HCO_CH2GSG) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - (CH2OH_CH2GSG) * (CH2GSG_CH2OH)) - ( (CH2_CH2GSG) &
                 - (HCO_CH2GSG) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) - ( ( ( (CH_CH2GSG) &
                 - (HCCO_CH2GSG) * (CH_HCCO)) - (HCO_CH2GSG) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * ( ( ( (CH2GSG_CH) &
                 - (HCCO_CH) * (CH2GSG_HCCO)) - (HCO_CH) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) )  

    cqss(sqssCH2GSG) = ( ( ( ( ( ( (CH2GSG_ct) - (HCCO_ct) * (CH2GSG_HCCO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - (CH2OH_ct) * (CH2GSG_CH2OH)) &
                 - ( ( (CH2_ct) - (CH2CHO_ct) * (CH2_CH2CHO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH2GSG_CH2) - (HCO_CH2) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) - ( ( ( (CH_ct) - (HCCO_ct) &
                 * (CH_HCCO)) - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( ( (CH2_ct) - (CH2CHO_ct) * (CH2_CH2CHO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 * ( ( ( (CH2GSG_CH) - (HCCO_CH) * (CH2GSG_HCCO)) &
                 - (HCO_CH) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH2GSG_CH2) &
                 - (HCO_CH2) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) - ( ( ( ( (C2H3_ct) - (CH2CHO_ct) &
                 * (C2H3_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( ( (CH2_ct) - (CH2CHO_ct) &
                 * (CH2_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 - ( ( ( (CH_ct) - (HCCO_ct) * (CH_HCCO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( ( (CH2_ct) &
                 - (CH2CHO_ct) * (CH2_CH2CHO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * ( ( (C2H3_CH) &
                 - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 / ( ( ( (1.0_pr) - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) * ( ( - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * (CH2GSG_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH_C2H3) &
                 - ( (HCO_C2H3) - (CH2CHO_C2H3) * (HCO_CH2CHO)) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( ( ( (CH2GSG_CH) - (HCCO_CH) * (CH2GSG_HCCO)) &
                 - (HCO_CH) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH2GSG_CH2) &
                 - (HCO_CH2) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))) / ( ( ( (1.0_pr) - (CH2CHO_C2H3) &
                 * (C2H3_CH2CHO)) - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH_C2H3) - ( (HCO_C2H3) &
                 - (CH2CHO_C2H3) * (HCO_CH2CHO)) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( ( (C2H3_CH) - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( - (HCO_CH2) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) - (HCCO_CH) &
                 * (CH_HCCO)) - (HCO_CH) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH2_CH) &
                 - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))))) &
                 - (CH2GSG_C2H5) * cqss(sqssC2H5)) / ( ( ( ( ( (1.0_pr) &
                 - (HCCO_CH2GSG) * (CH2GSG_HCCO)) - (HCO_CH2GSG) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - (CH2OH_CH2GSG) * (CH2GSG_CH2OH)) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH2GSG_CH2) &
                 - (HCO_CH2) * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) - ( ( ( (CH_CH2GSG) &
                 - (HCCO_CH2GSG) * (CH_HCCO)) - (HCO_CH2GSG) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * ( ( ( (CH2GSG_CH) &
                 - (HCCO_CH) * (CH2GSG_HCCO)) - (HCO_CH) * (CH2GSG_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH2GSG_CH2) - (HCO_CH2) &
                 * (CH2GSG_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))))  

    cqss(sqssC2H3) = ( ( ( ( (C2H3_ct) - (CH2CHO_ct) * (C2H3_CH2CHO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( ( (CH2_ct) - (CH2CHO_ct) * (CH2_CH2CHO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( - (HCO_CH2) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) - ( ( ( (CH_ct) &
                 - (HCCO_ct) * (CH_HCCO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( ( (CH2_ct) - (CH2CHO_ct) &
                 * (CH2_CH2CHO)) - ( (HCO_ct) - (CH2CHO_ct) &
                 * (HCO_CH2CHO)) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) * ( ( (C2H3_CH) - (HCO_CH) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 / ( ( ( (1.0_pr) - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) - ( ( ( - (HCO_CH2GSG) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH2_CH2GSG) &
                 - (HCO_CH2GSG) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( - (HCO_CH2) * ( (C2H3_HCO) &
                 - (CH2CHO_HCO) * (C2H3_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 - ( ( ( (CH_CH2GSG) - (HCCO_CH2GSG) * (CH_HCCO)) &
                 - (HCO_CH2GSG) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH2GSG) - (HCO_CH2GSG) &
                 * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))) &
                 * ( ( (C2H3_CH) - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( - (HCO_CH2) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) - (HCCO_CH) &
                 * (CH_HCCO)) - (HCO_CH) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH2_CH) &
                 - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))))) &
                 * cqss(sqssCH2GSG)) / ( ( ( (1.0_pr) - (CH2CHO_C2H3) &
                 * (C2H3_CH2CHO)) - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH_C2H3) - ( (HCO_C2H3) &
                 - (CH2CHO_C2H3) * (HCO_CH2CHO)) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( ( (C2H3_CH) - (HCO_CH) * ( (C2H3_HCO) - (CH2CHO_HCO) &
                 * (C2H3_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( - (HCO_CH2) &
                 * ( (C2H3_HCO) - (CH2CHO_HCO) * (C2H3_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO)))) / ( ( ( (1.0_pr) - (HCCO_CH) &
                 * (CH_HCCO)) - (HCO_CH) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( (CH2_CH) &
                 - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))))  

    cqss(sqssCH) = ( ( ( (CH_ct) - (HCCO_ct) * (CH_HCCO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( ( (CH2_ct) &
                 - (CH2CHO_ct) * (CH2_CH2CHO)) - ( (HCO_ct) &
                 - (CH2CHO_ct) * (HCO_CH2CHO)) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) - ( ( ( (CH_CH2GSG) &
                 - (HCCO_CH2GSG) * (CH_HCCO)) - (HCO_CH2GSG) &
                 * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * ( (CH_CH2) &
                 - (HCO_CH2) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO)))) * cqss(sqssCH2GSG) &
                 - ( (CH_C2H3) - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * (CH_HCO) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * cqss(sqssC2H3)) / ( ( ( (1.0_pr) &
                 - (HCCO_CH) * (CH_HCCO)) - (HCO_CH) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * ( (CH_CH2) - (HCO_CH2) * (CH_HCO) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))) &
                 / ( (1.0_pr) - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))))  

    cqss(sqssCH2) = ( ( (CH2_ct) - (CH2CHO_ct) * (CH2_CH2CHO)) &
                 - ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO)) &
                 * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)) &
                 - ( (CH2_CH2GSG) - (HCO_CH2GSG) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * cqss(sqssCH2GSG) &
                 - ( ( - (CH2CHO_C2H3) * (CH2_CH2CHO)) - ( (HCO_C2H3) &
                 - (CH2CHO_C2H3) * (HCO_CH2CHO)) * ( (CH2_HCO) &
                 - (CH2CHO_HCO) * (CH2_CH2CHO)) / ( (1.0_pr) &
                 - (CH2CHO_HCO) * (HCO_CH2CHO))) * cqss(sqssC2H3) &
                 - ( (CH2_CH) - (HCO_CH) * ( (CH2_HCO) - (CH2CHO_HCO) &
                 * (CH2_CH2CHO)) / ( (1.0_pr) - (CH2CHO_HCO) &
                 * (HCO_CH2CHO))) * cqss(sqssCH)) / ( (1.0_pr) &
                 - (HCO_CH2) * ( (CH2_HCO) - (CH2CHO_HCO) * (CH2_CH2CHO)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO)))  

    cqss(sqssCH2OH) = CH2OH_ct - (CH2OH_CH2GSG) * cqss(sqssCH2GSG)  

    cqss(sqssHCO) = ( (HCO_ct) - (CH2CHO_ct) * (HCO_CH2CHO) - (HCO_CH2GSG) &
                 * cqss(sqssCH2GSG) - ( (HCO_C2H3) - (CH2CHO_C2H3) &
                 * (HCO_CH2CHO)) * cqss(sqssC2H3) - (HCO_CH) &
                 * cqss(sqssCH) - (HCO_CH2) * cqss(sqssCH2)) &
                 / ( (1.0_pr) - (CH2CHO_HCO) * (HCO_CH2CHO))  

    cqss(sqssHCCO) = HCCO_ct - (HCCO_CH2GSG) * cqss(sqssCH2GSG) &
                 - (HCCO_CH) * cqss(sqssCH)  

    cqss(sqssCH2CHO) = CH2CHO_ct - (CH2CHO_C2H3) * cqss(sqssC2H3) &
                 - (CH2CHO_HCO) * cqss(sqssHCO)  

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

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

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

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

    ! Conversion
    c = y * rho / W_sp

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

        return
    end subroutine y2c

end module mod_customkinetics

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

  use mod_customkinetics
  implicit none

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

  real(pr) :: P, T, rho

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

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

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

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

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