!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE ./database/S16R129QSS10_2/S16R129QSS10_2.f90
!>    @file ./database/S16R129QSS10_2/S16R129QSS10_2.f90
!!    Module for calculating the analytical source terms in Cantera
!!    @details
!!    @authors Quentin Cazeres
!!    @date    2020/09/02
!!    @since   2.3
!!    @note
!--------------------------------------------------------------------------------------------------

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

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

  ! Ideal gas constant
  real(pr), parameter :: Rcst = 8.3144621_pr

  ! Use semi-implicit formulation for select species
  logical, parameter :: semiImplicit = .False.

  ! Number of elements in the chemical system
  integer, parameter :: ne = 3

  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 16
  integer, parameter :: nqss = 10
  integer, parameter :: nreac = 129
  integer, parameter :: nreac_reverse = 121

  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp

  ! Number of thirdbodies
  integer, parameter :: nTB = 4
  integer, parameter :: nFO = 15
  integer, parameter :: nTB_reverse = 4
  integer, parameter :: nFO_reverse = 15

  ! Index of elements
  integer, parameter :: eH = 1
  integer, parameter :: eO = 2
  integer, parameter :: eC = 3

  ! 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 :: sC2H4 = 14
  integer, parameter :: sC2H6 = 15
  integer, parameter :: sN2 = 16

  integer, parameter :: sqssCH = 1
  integer, parameter :: sqssCH2 = 2
  integer, parameter :: sqssCH2GSG = 3
  integer, parameter :: sqssHCO = 4
  integer, parameter :: sqssCH2OH = 5
  integer, parameter :: sqssCH3O = 6
  integer, parameter :: sqssCH3OH = 7
  integer, parameter :: sqssC2H2 = 8
  integer, parameter :: sqssC2H5 = 9
  integer, parameter :: sqssCH2CHO = 10

  ! 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 :: r1b = 130
  integer, parameter :: r2b = 131
  integer, parameter :: r3b = 132
  integer, parameter :: r4b = 133
  integer, parameter :: r5b = 134
  integer, parameter :: r6b = 135
  integer, parameter :: r7b = 136
  integer, parameter :: r8b = 137
  integer, parameter :: r9b = 138
  integer, parameter :: r10b = 139
  integer, parameter :: r11b = 140
  integer, parameter :: r12b = 141
  integer, parameter :: r13b = 142
  integer, parameter :: r14b = 143
  integer, parameter :: r15b = 144
  integer, parameter :: r16b = 145
  integer, parameter :: r17b = 146
  integer, parameter :: r18b = 147
  integer, parameter :: r19b = 148
  integer, parameter :: r20b = 149
  integer, parameter :: r21b = 150
  integer, parameter :: r22b = 151
  integer, parameter :: r23b = 152
  integer, parameter :: r24b = 153
  integer, parameter :: r25b = 154
  integer, parameter :: r26b = 155
  integer, parameter :: r27b = 156
  integer, parameter :: r28b = 157
  integer, parameter :: r29b = 158
  integer, parameter :: r30b = 159
  integer, parameter :: r31b = 160
  integer, parameter :: r32b = 161
  integer, parameter :: r33b = 162
  integer, parameter :: r34b = 163
  integer, parameter :: r35b = 164
  integer, parameter :: r36b = 165
  integer, parameter :: r37b = 166
  integer, parameter :: r38b = 167
  integer, parameter :: r39b = 168
  integer, parameter :: r40b = 169
  integer, parameter :: r41b = 170
  integer, parameter :: r42b = 171
  integer, parameter :: r43b = 172
  integer, parameter :: r44b = 173
  integer, parameter :: r45b = 174
  integer, parameter :: r46b = 175
  integer, parameter :: r47b = 176
  integer, parameter :: r48b = 177
  integer, parameter :: r49b = 178
  integer, parameter :: r50b = 179
  integer, parameter :: r51b = 180
  integer, parameter :: r52b = 181
  integer, parameter :: r53b = 182
  integer, parameter :: r54b = 183
  integer, parameter :: r55b = 184
  integer, parameter :: r56b = 185
  integer, parameter :: r57b = 186
  integer, parameter :: r58b = 187
  integer, parameter :: r59b = 188
  integer, parameter :: r60b = 189
  integer, parameter :: r61b = 190
  integer, parameter :: r62b = 191
  integer, parameter :: r63b = 192
  integer, parameter :: r64b = 193
  integer, parameter :: r65b = 194
  integer, parameter :: r66b = 195
  integer, parameter :: r67b = 196
  integer, parameter :: r68b = 197
  integer, parameter :: r69b = 198
  integer, parameter :: r70b = 199
  integer, parameter :: r71b = 200
  integer, parameter :: r72b = 201
  integer, parameter :: r73b = 202
  integer, parameter :: r74b = 203
  integer, parameter :: r75b = 204
  integer, parameter :: r76b = 205
  integer, parameter :: r77b = 206
  integer, parameter :: r78b = 207
  integer, parameter :: r79b = 208
  integer, parameter :: r80b = 209
  integer, parameter :: r81b = 210
  integer, parameter :: r82b = 211
  integer, parameter :: r83b = 212
  integer, parameter :: r84b = 213
  integer, parameter :: r85b = 214
  integer, parameter :: r86b = 215
  integer, parameter :: r88b = 216
  integer, parameter :: r89b = 217
  integer, parameter :: r90b = 218
  integer, parameter :: r91b = 219
  integer, parameter :: r92b = 220
  integer, parameter :: r93b = 221
  integer, parameter :: r94b = 222
  integer, parameter :: r95b = 223
  integer, parameter :: r96b = 224
  integer, parameter :: r97b = 225
  integer, parameter :: r98b = 226
  integer, parameter :: r99b = 227
  integer, parameter :: r100b = 228
  integer, parameter :: r101b = 229
  integer, parameter :: r102b = 230
  integer, parameter :: r103b = 231
  integer, parameter :: r104b = 232
  integer, parameter :: r105b = 233
  integer, parameter :: r106b = 234
  integer, parameter :: r107b = 235
  integer, parameter :: r108b = 236
  integer, parameter :: r109b = 237
  integer, parameter :: r110b = 238
  integer, parameter :: r111b = 239
  integer, parameter :: r112b = 240
  integer, parameter :: r113b = 241
  integer, parameter :: r114b = 242
  integer, parameter :: r115b = 243
  integer, parameter :: r116b = 244
  integer, parameter :: r117b = 245
  integer, parameter :: r119b = 246
  integer, parameter :: r120b = 247
  integer, parameter :: r122b = 248
  integer, parameter :: r124b = 249
  integer, parameter :: r129b = 250

  ! Index of third body species
  integer, parameter :: mM21 = 1
  integer, parameter :: mM26 = 2
  integer, parameter :: mM28 = 3
  integer, parameter :: mM112 = 4

  integer, parameter :: mM9 = 5
  integer, parameter :: mM33 = 6
  integer, parameter :: mM35 = 7
  integer, parameter :: mM37 = 8
  integer, parameter :: mM39 = 9
  integer, parameter :: mM40 = 10
  integer, parameter :: mM50 = 11
  integer, parameter :: mM51 = 12
  integer, parameter :: mM54 = 13
  integer, parameter :: mM56 = 14
  integer, parameter :: mM64 = 15
  integer, parameter :: mM95 = 16
  integer, parameter :: mM105 = 17
  integer, parameter :: mM116 = 18
  integer, parameter :: mM122 = 19

  ! Index of third body reactions
  integer, parameter :: TBr21f = 1
  integer, parameter :: TBr26f = 2
  integer, parameter :: TBr28f = 3
  integer, parameter :: TBr112f = 4
  integer, parameter :: TBr21b = 5
  integer, parameter :: TBr26b = 6
  integer, parameter :: TBr28b = 7
  integer, parameter :: TBr112b = 8

  ! Index of fall off reactions
  integer, parameter :: FOr9f = 1
  integer, parameter :: FOr33f = 2
  integer, parameter :: FOr35f = 3
  integer, parameter :: FOr37f = 4
  integer, parameter :: FOr39f = 5
  integer, parameter :: FOr40f = 6
  integer, parameter :: FOr50f = 7
  integer, parameter :: FOr51f = 8
  integer, parameter :: FOr54f = 9
  integer, parameter :: FOr56f = 10
  integer, parameter :: FOr64f = 11
  integer, parameter :: FOr95f = 12
  integer, parameter :: FOr105f = 13
  integer, parameter :: FOr116f = 14
  integer, parameter :: FOr122f = 15
  integer, parameter :: FOr9b = 16
  integer, parameter :: FOr33b = 17
  integer, parameter :: FOr35b = 18
  integer, parameter :: FOr37b = 19
  integer, parameter :: FOr39b = 20
  integer, parameter :: FOr40b = 21
  integer, parameter :: FOr50b = 22
  integer, parameter :: FOr51b = 23
  integer, parameter :: FOr54b = 24
  integer, parameter :: FOr56b = 25
  integer, parameter :: FOr64b = 26
  integer, parameter :: FOr95b = 27
  integer, parameter :: FOr105b = 28
  integer, parameter :: FOr116b = 29
  integer, parameter :: FOr122b = 30

  ! 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.02805376_pr, & ! C2H4
       0.03006964_pr, & ! C2H6
       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.03103422_pr, & ! CH3O
      !0.03204216_pr, & ! CH3OH
      !0.02603788_pr, & ! C2H2
      !0.0290617_pr, & ! C2H5
      !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 = abs(k0) * conc / max(abs(kinf), tiny(1.0_pr)) + tiny(1.0_pr)
    lgknull = log10(k0)
    f = (lgknull+ccoeff)/(ntmp-dcoeff*(lgknull+ccoeff))
    f = fc**(1.0_pr / ( f * f + 1.0_pr ))
    getlindratecoeff = kinf * f * k0 / ( 1.0_pr + k0 )

  end function getlindratecoeff

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

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

    M(mM21) = (-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(mM26) = (1_pr)*c(sCH4) &
         + (-1_pr)*c(sCO2) &
         + (-1_pr)*c(sH2) &
         + (-1_pr)*c(sH2O) &
         + sum(c)

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

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

    M(mM9) = (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(mM33) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

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

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

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

    M(mM40) = (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) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

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

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

    M(mM56) = (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) = (1_pr)*c(sCH4) &
         + (0.5_pr)*c(sCO) &
         + (1_pr)*c(sCO2) &
         + (1_pr)*c(sH2) &
         + (5_pr)*c(sH2O) &
         + sum(c)

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

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

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

    M(mM122) = (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) = (0.0387_pr)*exp((-26191.84_pr)*R_T_inv + T_log *(2.7_pr))
    k(r2f) = (20000000_pr)
    k(r3f) = (57000000_pr)
    k(r4f) = (80000000_pr)
    k(r5f) = (15000000_pr)
    k(r6f) = (15000000_pr)
    k(r7f) = (50600000_pr)
    k(r8f) = (1020_pr)*exp((-35982.4_pr)*R_T_inv + T_log *(1.5_pr))
    k_0(FOr9f) = (602_pr)*exp((-12552_pr)*R_T_inv)
    k_inf(FOr9f) = (18000_pr)*exp((-9978.84_pr)*R_T_inv)
    k(r9f) = getlindratecoeff(Tloc,k_0(FOr9f),k_inf(FOr9f),1.0_pr, M(mM9),Ploc)
    k(r10f) = (30000000_pr)
    k(r11f) = (30000000_pr)
    k(r12f) = (39000000_pr)*exp((-14811.36_pr)*R_T_inv)
    k(r13f) = (10000000_pr)
    k(r14f) = (0.388_pr)*exp((-12970.4_pr)*R_T_inv + T_log *(2.5_pr))
    k(r15f) = (6.94_pr)*exp((-7949.6_pr)*R_T_inv + T_log *(2_pr))
    k(r16f) = (12.5_pr)*exp((-920.48_pr)*R_T_inv + T_log *(1.83_pr))
    k(r17f) = (22400000_pr)
    k(r18f) = (89.8_pr)*exp((-23806.96_pr)*R_T_inv + T_log *(1.92_pr))
    k(r19f) = (2500000_pr)*exp((-199995.2_pr)*R_T_inv)
    k(r20f) = (100000000_pr)*exp((-167360_pr)*R_T_inv)
    k(r21f) = (2800000_pr)*exp(T_log *(-0.86_pr))
    k(r22f) = (20800000_pr)*exp(T_log *(-1.24_pr))
    k(r23f) = (11260000_pr)*exp(T_log *(-0.76_pr))
    k(r24f) = (26000000_pr)*exp(T_log *(-1.24_pr))
    k(r25f) = (26500000000_pr)*exp((-71299.544_pr)*R_T_inv + T_log *(-0.6707_pr))
    k(r26f) = (1000000_pr)*exp(T_log *(-1_pr))
    k(r27f) = (60000000_pr)*exp(T_log *(-1.25_pr))
    k(r28f) = (22000000000_pr)*exp(T_log *(-2_pr))
    k(r29f) = (3970000_pr)*exp((-2807.464_pr)*R_T_inv)
    k(r30f) = (44800000_pr)*exp((-4468.512_pr)*R_T_inv)
    k(r31f) = (84000000_pr)*exp((-2656.84_pr)*R_T_inv)
    k(r32f) = (12.1_pr)*exp((-21756.8_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr33f) = (104000000000000_pr)*exp((-6694.4_pr)*R_T_inv + T_log *(-2.76_pr))
    k_inf(FOr33f) = (600000000_pr)
    FC(FOr33f) = ((1.0_pr - 0.562_pr)*exp(-Tloc/(91_pr))) + (0.562_pr*exp(-Tloc/(5836_pr)))&
              + exp(-(8552_pr)/Tloc)
    k(r33f) = getlindratecoeff(Tloc,k_0(FOr33f),k_inf(FOr33f),FC(FOr33f), M(mM33),Ploc)
    k(r34f) = (30000000_pr)
    k_0(FOr35f) = (2.62e+21_pr)*exp((-10208.96_pr)*R_T_inv + T_log *(-4.76_pr))
    k_inf(FOr35f) = (13900000000_pr)*exp((-2242.624_pr)*R_T_inv + T_log *(-0.534_pr))
    FC(FOr35f) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r35f) = getlindratecoeff(Tloc,k_0(FOr35f),k_inf(FOr35f),FC(FOr35f), M(mM35),Ploc)
    k(r36f) = (660_pr)*exp((-45354.56_pr)*R_T_inv + T_log *(1.62_pr))
    k_0(FOr37f) = (2470000000000_pr)*exp((-1778.2_pr)*R_T_inv + T_log *(-2.57_pr))
    k_inf(FOr37f) = (1090000_pr)*exp((1087.84_pr)*R_T_inv + T_log *(0.48_pr))
    FC(FOr37f) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r37f) = getlindratecoeff(Tloc,k_0(FOr37f),k_inf(FOr37f),FC(FOr37f), M(mM37),Ploc)
    k(r38f) = (73400000_pr)
    k_0(FOr39f) = (1.27e+20_pr)*exp((-27321.52_pr)*R_T_inv + T_log *(-4.82_pr))
    k_inf(FOr39f) = (540000_pr)*exp((-15062.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr39f) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r39f) = getlindratecoeff(Tloc,k_0(FOr39f),k_inf(FOr39f),FC(FOr39f), M(mM39),Ploc)
    k_0(FOr40f) = (2.2e+18_pr)*exp((-23263.04_pr)*R_T_inv + T_log *(-4.8_pr))
    k_inf(FOr40f) = (540000_pr)*exp((-10878.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr40f) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r40f) = getlindratecoeff(Tloc,k_0(FOr40f),k_inf(FOr40f),FC(FOr40f), M(mM40),Ploc)
    k(r41f) = (57.4_pr)*exp((-11472.528_pr)*R_T_inv + T_log *(1.9_pr))
    k(r42f) = (20000000_pr)
    k(r43f) = (165000_pr)*exp((1188.256_pr)*R_T_inv + T_log *(0.65_pr))
    k(r44f) = (32800000_pr)*exp((-2552.24_pr)*R_T_inv + T_log *(-0.09_pr))
    k(r45f) = (20000000_pr)
    k(r46f) = (1500000_pr)*exp((460.24_pr)*R_T_inv + T_log *(0.5_pr))
    k(r47f) = (262000000_pr)*exp((-4476.88_pr)*R_T_inv + T_log *(-0.23_pr))
    k(r48f) = (17_pr)*exp((-20376.08_pr)*R_T_inv + T_log *(2.1_pr))
    k(r49f) = (4.2_pr)*exp((-20376.08_pr)*R_T_inv + T_log *(2.1_pr))
    k_0(FOr50f) = (6e+29_pr)*exp((-29162.48_pr)*R_T_inv + T_log *(-7.62_pr))
    k_inf(FOr50f) = (540000_pr)*exp((-7614.88_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr50f) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r50f) = getlindratecoeff(Tloc,k_0(FOr50f),k_inf(FOr50f),FC(FOr50f), M(mM50),Ploc)
    k_0(FOr51f) = (1.99e+29_pr)*exp((-27970.04_pr)*R_T_inv + T_log *(-7.08_pr))
    k_inf(FOr51f) = (521000000000_pr)*exp((-6610.72_pr)*R_T_inv + T_log *(-0.99_pr))
    FC(FOr51f) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r51f) = getlindratecoeff(Tloc,k_0(FOr51f),k_inf(FOr51f),FC(FOr51f), M(mM51),Ploc)
    k(r52f) = (2000000_pr)
    k(r53f) = (115_pr)*exp((-31505.52_pr)*R_T_inv + T_log *(1.9_pr))
    k_0(FOr54f) = (5.07e+15_pr)*exp((-352920.4_pr)*R_T_inv + T_log *(-3.42_pr))
    k_inf(FOr54f) = (43_pr)*exp((-333046.4_pr)*R_T_inv + T_log *(1.5_pr))
    FC(FOr54f) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r54f) = getlindratecoeff(Tloc,k_0(FOr54f),k_inf(FOr54f),FC(FOr54f), M(mM54),Ploc)
    k(r55f) = (216_pr)*exp((-14351.12_pr)*R_T_inv + T_log *(1.51_pr))
    k_0(FOr56f) = (2300000_pr)*exp((7112.8_pr)*R_T_inv + T_log *(-0.9_pr))
    k_inf(FOr56f) = (74000000_pr)*exp(T_log *(-0.37_pr))
    FC(FOr56f) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r56f) = getlindratecoeff(Tloc,k_0(FOr56f),k_inf(FOr56f),FC(FOr56f), M(mM56),Ploc)
    k(r57f) = (0.0357_pr)*exp((8828.24_pr)*R_T_inv + T_log *(2.4_pr))
    k(r58f) = (14500000_pr)*exp((2092_pr)*R_T_inv)
    k(r59f) = (1700000000000_pr)*exp((-123051.44_pr)*R_T_inv)
    k(r60f) = (30000000_pr)
    k(r61f) = (20000000_pr)
    k(r62f) = (11.3_pr)*exp((-12552_pr)*R_T_inv + T_log *(2_pr))
    k(r63f) = (30000000_pr)
    k_0(FOr64f) = (4e+24_pr)*exp((-13137.76_pr)*R_T_inv + T_log *(-5.92_pr))
    k_inf(FOr64f) = (2790000000000_pr)*exp((-5564.72_pr)*R_T_inv + T_log *(-1.43_pr))
    FC(FOr64f) = ((1.0_pr - 0.412_pr)*exp(-Tloc/(195_pr))) + (0.412_pr*exp(-Tloc/(5900_pr)))&
              + exp(-(6394_pr)/Tloc)
    k(r64f) = getlindratecoeff(Tloc,k_0(FOr64f),k_inf(FOr64f),FC(FOr64f), M(mM64),Ploc)
    k(r65f) = (56_pr)*exp((-22677.28_pr)*R_T_inv + T_log *(1.6_pr))
    k(r66f) = (644000000000_pr)*exp((-5928.728_pr)*R_T_inv + T_log *(-1.34_pr))
    k(r67f) = (100_pr)*exp((-13054.08_pr)*R_T_inv + T_log *(1.6_pr))
    k(r68f) = (47.6_pr)*exp((-292.88_pr)*R_T_inv + T_log *(1.228_pr))
    k(r69f) = (50000000_pr)
    k(r70f) = (3430_pr)*exp((1870.248_pr)*R_T_inv + T_log *(1.18_pr))
    k(r71f) = (5000000_pr)
    k(r72f) = (1.44_pr)*exp((3514.56_pr)*R_T_inv + T_log *(2_pr))
    k(r73f) = (6.3_pr)*exp((-6276_pr)*R_T_inv + T_log *(2_pr))
    k(r74f) = (4.83e-10_pr)*exp((8368_pr)*R_T_inv + T_log *(4_pr))
    k(r75f) = (3.54_pr)*exp((-3640.08_pr)*R_T_inv + T_log *(2.12_pr))
    k(r76f) = (420000000_pr)*exp((-50208_pr)*R_T_inv)
    k(r77f) = (20000000_pr)
    k(r78f) = (1000000_pr)
    k(r79f) = (37800000_pr)
    k(r80f) = (150000000_pr)*exp((-98742.4_pr)*R_T_inv)
    k(r81f) = (5.6_pr)*exp((-50208_pr)*R_T_inv + T_log *(2_pr))
    k(r82f) = (67100000_pr)
    k(r83f) = (108000000_pr)*exp((-13012.24_pr)*R_T_inv)
    k(r84f) = (5710000_pr)*exp((3158.92_pr)*R_T_inv)
    k(r85f) = (60000000_pr)
    k(r86f) = (190000000_pr)*exp((-66073.728_pr)*R_T_inv)
    k(r87f) = (5000000_pr)*exp((-6276_pr)*R_T_inv)
    k(r88f) = (0.5_pr)*exp((-30250.32_pr)*R_T_inv + T_log *(2_pr))
    k(r89f) = (40000000_pr)
    k(r90f) = (2.46_pr)*exp((-34601.68_pr)*R_T_inv + T_log *(2_pr))
    k(r91f) = (15000000_pr)*exp((-2510.4_pr)*R_T_inv)
    k(r92f) = (28000000_pr)
    k(r93f) = (12000000_pr)
    k(r94f) = (70000000_pr)
    k_0(FOr95f) = (1.88e+26_pr)*exp((-21087.36_pr)*R_T_inv + T_log *(-6.36_pr))
    k_inf(FOr95f) = (482000000000_pr)*exp((-4790.68_pr)*R_T_inv + T_log *(-1.16_pr))
    FC(FOr95f) = ((1.0_pr - 0.6027_pr)*exp(-Tloc/(208_pr))) + (0.6027_pr*exp(-Tloc/(3922_pr)))&
              + exp(-(10180_pr)/Tloc)
    k(r95f) = getlindratecoeff(Tloc,k_0(FOr95f),k_inf(FOr95f),FC(FOr95f), M(mM95),Ploc)
    k(r96f) = (30000000_pr)
    k(r97f) = (12000000_pr)*exp((2384.88_pr)*R_T_inv)
    k(r98f) = (16000000_pr)*exp((2384.88_pr)*R_T_inv)
    k(r99f) = (9000000_pr)
    k(r100f) = (7000000_pr)
    k(r101f) = (14000000_pr)
    k(r102f) = (35600000_pr)*exp((-127528.32_pr)*R_T_inv)
    k(r103f) = (2310000_pr)*exp((-84997.96_pr)*R_T_inv)
    k(r104f) = (0.0245_pr)*exp((-21673.12_pr)*R_T_inv + T_log *(2.47_pr))
    k_0(FOr105f) = (3.4e+29_pr)*exp((-11556.208_pr)*R_T_inv + T_log *(-7.03_pr))
    k_inf(FOr105f) = (67700000000_pr)*exp((-2736.336_pr)*R_T_inv + T_log *(-1.18_pr))
    FC(FOr105f) = ((1.0_pr - 0.619_pr)*exp(-Tloc/(73.2_pr))) + (0.619_pr*exp(-Tloc/(1180_pr)))&
              + exp(-(9999_pr)/Tloc)
    k(r105f) = getlindratecoeff(Tloc,k_0(FOr105f),k_inf(FOr105f),FC(FOr105f), M(mM105),Ploc)
    k(r106f) = (6840000_pr)*exp((-44350.4_pr)*R_T_inv + T_log *(0.1_pr))
    k(r107f) = (26480000_pr)
    k(r108f) = (0.00332_pr)*exp((-24518.24_pr)*R_T_inv + T_log *(2.81_pr))
    k(r109f) = (10_pr)*exp((-41588.96_pr)*R_T_inv + T_log *(1.5_pr))
    k(r110f) = (6.14_pr)*exp((-43722.8_pr)*R_T_inv + T_log *(1.74_pr))
    k(r111f) = (1500000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r112f) = (187000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r113f) = (13450000_pr)*exp((-1673.6_pr)*R_T_inv)
    k(r114f) = (18000000_pr)*exp((-3765.6_pr)*R_T_inv)
    k(r115f) = (4.28e-19_pr)*exp((14769.52_pr)*R_T_inv + T_log *(7.6_pr))
    k_0(FOr116f) = (1.58e+45_pr)*exp((-409195.2_pr)*R_T_inv + T_log *(-9.3_pr))
    k_inf(FOr116f) = (8000000000000_pr)*exp((-363045.68_pr)*R_T_inv + T_log *(0.44_pr))
    FC(FOr116f) = ((1.0_pr - 0.7345_pr)*exp(-Tloc/(180_pr))) + (0.7345_pr*exp(-Tloc/(1035_pr)))&
              + exp(-(5417_pr)/Tloc)
    k(r116f) = getlindratecoeff(Tloc,k_0(FOr116f),k_inf(FOr116f),FC(FOr116f), M(mM116),Ploc)
    k(r117f) = (840000_pr)*exp((-16213_pr)*R_T_inv)
    k(r118f) = (33700000_pr)
    k(r119f) = (6.7_pr)*exp((-920.48_pr)*R_T_inv + T_log *(1.83_pr))
    k(r120f) = (5000000000_pr)*exp((-72508.72_pr)*R_T_inv)
    k(r121f) = (8000_pr)*exp((7342.92_pr)*R_T_inv + T_log *(0.5_pr))
    k_0(FOr122f) = (48200000000000_pr)*exp((-2468.56_pr)*R_T_inv + T_log *(-2.8_pr))
    k_inf(FOr122f) = (1970000_pr)*exp((1548.08_pr)*R_T_inv + T_log *(0.43_pr))
    FC(FOr122f) = ((1.0_pr - 0.578_pr)*exp(-Tloc/(122_pr))) + (0.578_pr*exp(-Tloc/(2535_pr)))&
              + exp(-(9365_pr)/Tloc)
    k(r122f) = getlindratecoeff(Tloc,k_0(FOr122f),k_inf(FOr122f),FC(FOr122f), M(mM122),Ploc)
    k(r123f) = (5800000_pr)*exp((-6276_pr)*R_T_inv)
    k(r124f) = (2400000_pr)*exp((-6276_pr)*R_T_inv)
    k(r125f) = (68200_pr)*exp((3912.04_pr)*R_T_inv + T_log *(0.25_pr))
    k(r126f) = (150000000_pr)
    k(r127f) = (18100_pr)
    k(r128f) = (23500_pr)
    k(r129f) = (22000000_pr)
    k(r1b) = (0.0206329471029809_pr)*exp((-18166.1522806783_pr)*R_T_inv + T_log *(2.67792789463785_pr))
    k(r2b) = (3011784.28657983_pr)*exp((-221106.574347875_pr)*R_T_inv + T_log *(0.274363527908898_pr))
    k(r3b) = (2650359855.1806_pr)*exp((-739091.916935856_pr)*R_T_inv + T_log *(-0.000506956884162345_pr))
    k(r4b) = (3952908618.60814_pr)*exp((-382650.966331323_pr)*R_T_inv + T_log *(-0.261298831611412_pr))
    k(r5b) = (80410533.6817515_pr)*exp((-788725.617543665_pr)*R_T_inv + T_log *(0.11285818119583_pr))
    k(r6b) = (663934446.502167_pr)*exp((-420239.943003471_pr)*R_T_inv + T_log *(-0.347447999313118_pr))
    k(r7b) = (2598937783.74877_pr)*exp((-287827.834672731_pr)*R_T_inv + T_log *(-0.159567518165995_pr))
    k(r8b) = (1.6196842604001_pr)*exp((-21349.5453478993_pr)*R_T_inv + T_log *(1.80326072383876_pr))
    k_0(FOr9b) = (152263054883873_pr)*exp((-545805.078078398_pr)*R_T_inv + T_log *(-0.833125395867326_pr))
    k_inf(FOr9b) = (4.55271590087635e+15_pr)*exp((-543231.918074133_pr)*R_T_inv + T_log *(-0.833125395124001_pr))
    k(r9b) = getlindratecoeff(Tloc,k_0(FOr9b),k_inf(FOr9b),1.0_pr, M(mM9),Ploc)
    k(r10b) = (1937132.30041074_pr)*exp((-360459.986796635_pr)*R_T_inv + T_log *(0.43823407956681_pr))
    k(r11b) = (3388094297002.73_pr)*exp((-468721.011341084_pr)*R_T_inv + T_log *(-0.69876121480493_pr))
    k(r12b) = (130049.624236063_pr)*exp((-71994.7127952712_pr)*R_T_inv + T_log *(0.271516409392677_pr))
    k(r13b) = (126190.532556803_pr)*exp((-302192.897291954_pr)*R_T_inv + T_log *(0.705086465434834_pr))
    k(r14b) = (0.00887054759646716_pr)*exp((-34627.174054719_pr)*R_T_inv + T_log *(2.48400824495903_pr))
    k(r15b) = (2.97353381645469e-05_pr)*exp((-199523.552595021_pr)*R_T_inv + T_log *(3.2583902331948_pr))
    k(r16b) = (0.000578286461866481_pr)*exp((-110137.848571894_pr)*R_T_inv + T_log *(2.54437389488512_pr))
    k(r17b) = (40757.8438955762_pr)*exp((-326594.351820143_pr)*R_T_inv + T_log *(0.811360287107584_pr))
    k(r18b) = (0.0322505812737488_pr)*exp((-28847.2023510494_pr)*R_T_inv + T_log *(2.30659367941472_pr))
    k(r19b) = (16638595957.5674_pr)*exp((-236321.847275149_pr)*R_T_inv + T_log *(-0.780968881256363_pr))
    k(r20b) = (2214374.19227209_pr)*exp((-3436.77937385951_pr)*R_T_inv + T_log *(-0.00284722293883299_pr))
    k(r21b) = (2688846218230.72_pr)*exp((-203885.479118053_pr)*R_T_inv + T_log *(-0.830493617224921_pr))
    k(r22b) = (19974286260286.7_pr)*exp((-203885.479120324_pr)*R_T_inv + T_log *(-1.21049361767258_pr))
    k(r23b) = (10813003290149.5_pr)*exp((-203885.47913709_pr)*R_T_inv + T_log *(-0.730493620641542_pr))
    k(r24b) = (24967857855156.7_pr)*exp((-203885.479121343_pr)*R_T_inv + T_log *(-1.21049361782343_pr))
    k(r25b) = (100838485.591073_pr)*exp((634.83330177955_pr)*R_T_inv + T_log *(-0.314673580885157_pr))
    k(r26b) = (271238628429.747_pr)*exp((-433017.741249994_pr)*R_T_inv + T_log *(-0.674057995545302_pr))
    k(r27b) = (16274316836688.4_pr)*exp((-433017.741211138_pr)*R_T_inv + T_log *(-0.924057988602065_pr))
    k(r28b) = (3.71189323627178e+16_pr)*exp((-497041.49105575_pr)*R_T_inv + T_log *(-1.71245710875496_pr))
    k(r29b) = (26542.130037223_pr)*exp((-224029.098573096_pr)*R_T_inv + T_log *(0.614062938719847_pr))
    k(r30b) = (12653817.441446_pr)*exp((-233600.774038532_pr)*R_T_inv + T_log *(0.296435638410145_pr))
    k(r31b) = (48134.1808569423_pr)*exp((-151829.037030551_pr)*R_T_inv + T_log *(0.63038994966306_pr))
    k(r32b) = (0.0874533271090287_pr)*exp((-89430.2464641844_pr)*R_T_inv + T_log *(2.51410249200596_pr))
    k_0(FOr33b) = (4.33882456584479e+21_pr)*exp((-469326.232368474_pr)*R_T_inv + T_log *(-2.82937781922589_pr))
    k_inf(FOr33b) = (2.50316969383812e+16_pr)*exp((-462631.832855263_pr)*R_T_inv + T_log *(-0.0693779062300217_pr))
    FC(FOr33b) = ((1.0_pr - 0.562_pr)*exp(-Tloc/(91_pr))) + (0.562_pr*exp(-Tloc/(5836_pr)))&
              + exp(-(8552_pr)/Tloc)
    k(r33b) = getlindratecoeff(Tloc,k_0(FOr33b),k_inf(FOr33b),FC(FOr33b), M(mM33),Ploc)
    k(r34b) = (3458685.06005046_pr)*exp((-49633.6975949257_pr)*R_T_inv + T_log *(0.11336572085656_pr))
    k_0(FOr35b) = (2.38601356409496e+29_pr)*exp((-449833.868128927_pr)*R_T_inv + T_log *(-4.75939081505185_pr))
    k_inf(FOr35b) = (1.26586262085018e+18_pr)*exp((-441867.532394821_pr)*R_T_inv + T_log *(-0.5333908625885_pr))
    FC(FOr35b) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r35b) = getlindratecoeff(Tloc,k_0(FOr35b),k_inf(FOr35b),FC(FOr35b), M(mM35),Ploc)
    k(r36b) = (1.96572584270787_pr)*exp((-38747.3915781591_pr)*R_T_inv + T_log *(1.94533309150892_pr))
    k_0(FOr37b) = (1.07115991608643e+20_pr)*exp((-369586.900759822_pr)*R_T_inv + T_log *(-2.53764651466939_pr))
    k_inf(FOr37b) = (47269813076636.7_pr)*exp((-366720.860806242_pr)*R_T_inv + T_log *(0.512353477032747_pr))
    FC(FOr37b) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r37b) = getlindratecoeff(Tloc,k_0(FOr37b),k_inf(FOr37b),FC(FOr37b), M(mM37),Ploc)
    k(r38b) = (8889632.04242039_pr)*exp((-368485.674511254_pr)*R_T_inv + T_log *(0.460306185862355_pr))
    k_0(FOr39b) = (1.45538796893112e+27_pr)*exp((-150120.676227656_pr)*R_T_inv + T_log *(-5.22121656434435_pr))
    k_inf(FOr39b) = (6188264438792.04_pr)*exp((-137861.55630219_pr)*R_T_inv + T_log *(0.0527834223377546_pr))
    FC(FOr39b) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r39b) = getlindratecoeff(Tloc,k_0(FOr39b),k_inf(FOr39b),FC(FOr39b), M(mM39),Ploc)
    k_0(FOr40b) = (7.29297525428035e+25_pr)*exp((-116015.724655564_pr)*R_T_inv + T_log *(-5.01206208381703_pr))
    k_inf(FOr40b) = (17900937359160.5_pr)*exp((-103631.084578208_pr)*R_T_inv + T_log *(0.241937929989484_pr))
    FC(FOr40b) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r40b) = getlindratecoeff(Tloc,k_0(FOr40b),k_inf(FOr40b),FC(FOr40b), M(mM40),Ploc)
    k(r41b) = (0.359009541928841_pr)*exp((-76681.5683530109_pr)*R_T_inv + T_log *(2.19358854334518_pr))
    k(r42b) = (473376.247239249_pr)*exp((-310218.584999948_pr)*R_T_inv + T_log *(0.727158572800848_pr))
    k(r43b) = (40.5382809898171_pr)*exp((-13176.8065809668_pr)*R_T_inv + T_log *(1.51465399035595_pr))
    k(r44b) = (363.816248265675_pr)*exp((-13737.9845523932_pr)*R_T_inv + T_log *(1.21772387655136_pr))
    k(r45b) = (163643.755960311_pr)*exp((-340265.056531946_pr)*R_T_inv + T_log *(0.538004099418597_pr))
    k(r46b) = (127.398891234893_pr)*exp((-43951.2941607302_pr)*R_T_inv + T_log *(1.17549950851338_pr))
    k(r47b) = (1004.62165708942_pr)*exp((-45709.0961199229_pr)*R_T_inv + T_log *(0.888569396867642_pr))
    k(r48b) = (0.72898561734933_pr)*exp((-50058.544541109_pr)*R_T_inv + T_log *(2.10607986268624_pr))
    k(r49b) = (0.520984449469503_pr)*exp((-20012.0703126441_pr)*R_T_inv + T_log *(2.29523481134728_pr))
    k_0(FOr50b) = (6.61573998624193e+35_pr)*exp((-179594.197493771_pr)*R_T_inv + T_log *(-7.68463290374203_pr))
    k_inf(FOr50b) = (595416608253.056_pr)*exp((-158046.597505444_pr)*R_T_inv + T_log *(0.389367094187727_pr))
    FC(FOr50b) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r50b) = getlindratecoeff(Tloc,k_0(FOr50b),k_inf(FOr50b),FC(FOr50b), M(mM50),Ploc)
    k_0(FOr51b) = (8.01297302882961e+37_pr)*exp((-447921.851137878_pr)*R_T_inv + T_log *(-7.16272377287274_pr))
    k_inf(FOr51b) = (2.09786885309244e+20_pr)*exp((-426562.531149954_pr)*R_T_inv + T_log *(-1.07272377504382_pr))
    FC(FOr51b) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r51b) = getlindratecoeff(Tloc,k_0(FOr51b),k_inf(FOr51b),FC(FOr51b), M(mM51),Ploc)
    k(r52b) = (491987.816427886_pr)*exp((-282586.023740342_pr)*R_T_inv + T_log *(0.390574911039666_pr))
    k(r53b) = (0.0774655754696939_pr)*exp((-44571.4500721255_pr)*R_T_inv + T_log *(2.30866578450944_pr))
    k_0(FOr54b) = (1.81542210644361e+24_pr)*exp((-352243.426292636_pr)*R_T_inv + T_log *(-3.84795270834916_pr))
    k_inf(FOr54b) = (15397072034.6088_pr)*exp((-332369.426335381_pr)*R_T_inv + T_log *(1.07204728390438_pr))
    FC(FOr54b) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r54b) = getlindratecoeff(Tloc,k_0(FOr54b),k_inf(FOr54b),FC(FOr54b), M(mM54),Ploc)
    k(r55b) = (1343.61549379748_pr)*exp((-78374.8697995786_pr)*R_T_inv + T_log *(1.47160088792134_pr))
    k_0(FOr56b) = (1.50630915377225e+17_pr)*exp((-209059.297764783_pr)*R_T_inv + T_log *(-1.71855043877653_pr))
    k_inf(FOr56b) = (4.84638578819488e+18_pr)*exp((-216172.097737142_pr)*R_T_inv + T_log *(-1.18855043382193_pr))
    FC(FOr56b) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r56b) = getlindratecoeff(Tloc,k_0(FOr56b),k_inf(FOr56b),FC(FOr56b), M(mM56),Ploc)
    k(r57b) = (0.416523260153799_pr)*exp((-63221.1976723598_pr)*R_T_inv + T_log *(2.38367296625554_pr))
    k(r58b) = (25476094.0731401_pr)*exp((-291064.011903378_pr)*R_T_inv + T_log *(0.258036514637478_pr))
    k(r59b) = (76429514883.233_pr)*exp((-254748.636201363_pr)*R_T_inv + T_log *(0.475703390885456_pr))
    k(r60b) = (21603050324.5977_pr)*exp((-378631.933031913_pr)*R_T_inv + T_log *(-0.438741597813832_pr))
    k(r61b) = (296355038458.805_pr)*exp((-325467.613580337_pr)*R_T_inv + T_log *(-0.532815248978477_pr))
    k(r62b) = (9.04653626254542_pr)*exp((-88620.470764269_pr)*R_T_inv + T_log *(2.16111576922577_pr))
    k(r63b) = (398208644105.312_pr)*exp((-363056.590175417_pr)*R_T_inv + T_log *(-0.618964402961445_pr))
    k_0(FOr64b) = (1.02982092602499e+35_pr)*exp((-402107.976846437_pr)*R_T_inv + T_log *(-6.46479232758987_pr))
    k_inf(FOr64b) = (7.1830010870007e+22_pr)*exp((-394534.93685917_pr)*R_T_inv + T_log *(-1.9747923299281_pr))
    FC(FOr64b) = ((1.0_pr - 0.412_pr)*exp(-Tloc/(195_pr))) + (0.412_pr*exp(-Tloc/(5900_pr)))&
              + exp(-(6394_pr)/Tloc)
    k(r64b) = getlindratecoeff(Tloc,k_0(FOr64b),k_inf(FOr64b),FC(FOr64b), M(mM64),Ploc)
    k(r65b) = (2.26475978321703_pr)*exp((-57086.9388117409_pr)*R_T_inv + T_log *(1.95692068869599_pr))
    k(r66b) = (29074531941.6556_pr)*exp((-2749.40997430892_pr)*R_T_inv + T_log *(-0.896930114313579_pr))
    k(r67b) = (1.85268277893599_pr)*exp((-70470.6627348669_pr)*R_T_inv + T_log *(1.88693374032154_pr))
    k(r68b) = (83253629.4091096_pr)*exp((-108553.904571856_pr)*R_T_inv + T_log *(0.091004700687006_pr))
    k(r69b) = (37668554.3825608_pr)*exp((-432509.424263376_pr)*R_T_inv + T_log *(0.421907082201754_pr))
    k(r70b) = (133.447220570628_pr)*exp((-127362.542256477_pr)*R_T_inv + T_log *(1.43518941292305_pr))
    k(r71b) = (736152.386811317_pr)*exp((-374242.334801735_pr)*R_T_inv + T_log *(0.688759460302075_pr))
    k(r72b) = (0.384106918717434_pr)*exp((-90191.6515410707_pr)*R_T_inv + T_log *(1.9676812440437_pr))
    k(r73b) = (4.86112950549529_pr)*exp((-69935.7399761492_pr)*R_T_inv + T_log *(2.15683572326573_pr))
    k(r74b) = (5.97031810434878e-13_pr)*exp((-220845.731461799_pr)*R_T_inv + T_log *(4.8851425096356_pr))
    k(r75b) = (0.0148332201689907_pr)*exp((-80729.7599017697_pr)*R_T_inv + T_log *(2.4902666671264_pr))
    k(r76b) = (16413527988.2691_pr)*exp((-211666.815684592_pr)*R_T_inv + T_log *(-0.217666873144441_pr))
    k(r77b) = (169819130.10676_pr)*exp((-474639.810606179_pr)*R_T_inv + T_log *(0.0975747015016859_pr))
    k(r78b) = (94833894.4594763_pr)*exp((-235739.428846986_pr)*R_T_inv + T_log *(-0.0288971684811179_pr))
    k(r79b) = (255030246.231868_pr)*exp((-104760.662836931_pr)*R_T_inv + T_log *(-0.0451095529287638_pr))
    k(r80b) = (150335589016.601_pr)*exp((-356175.621639624_pr)*R_T_inv + T_log *(-0.506605356325564_pr))
    k(r81b) = (4.84608755868487_pr)*exp((-47743.5940526735_pr)*R_T_inv + T_log *(1.77948602246129_pr))
    k(r82b) = (183864048.116401_pr)*exp((-306697.555786768_pr)*R_T_inv + T_log *(-0.0827151889224419_pr))
    k(r83b) = (839152989.801403_pr)*exp((-967.519074679188_pr)*R_T_inv + T_log *(-0.199514888362545_pr))
    k(r84b) = (105685296949.529_pr)*exp((-246240.22286528_pr)*R_T_inv + T_log *(-0.693931026872802_pr))
    k(r85b) = (1483003917134.29_pr)*exp((-254781.709886815_pr)*R_T_inv + T_log *(-0.84985478302744_pr))
    k(r86b) = (78226.022243568_pr)*exp((-336444.636537656_pr)*R_T_inv + T_log *(0.698253687627095_pr))
    k(r88b) = (76.9054706556152_pr)*exp((-59864.4111493516_pr)*R_T_inv + T_log *(1.60468017076929_pr))
    k(r89b) = (42722211321261.3_pr)*exp((-273433.597687662_pr)*R_T_inv + T_log *(-0.975672713663772_pr))
    k(r90b) = (1.12694372463172_pr)*exp((-57608.6042196744_pr)*R_T_inv + T_log *(1.93001299935225_pr))
    k(r91b) = (13436879.2635952_pr)*exp((-40099.3766731866_pr)*R_T_inv + T_log *(-0.0861491680249256_pr))
    k(r92b) = (2.10575841485329_pr)*exp((-283773.498865759_pr)*R_T_inv + T_log *(0.142942617602521_pr))
    k(r93b) = (1522665.76273333_pr)*exp((-780814.989970731_pr)*R_T_inv + T_log *(0.430485500686761_pr))
    k(r94b) = (9644781971.51236_pr)*exp((-67203.0678035724_pr)*R_T_inv + T_log *(-0.481468993862963_pr))
    k_0(FOr95b) = (1.07209364519654e+38_pr)*exp((-413236.894899233_pr)*R_T_inv + T_log *(-7.34786221802827_pr))
    k_inf(FOr95b) = (2.74866554880215e+23_pr)*exp((-396940.214879663_pr)*R_T_inv + T_log *(-2.14786221456293_pr))
    FC(FOr95b) = ((1.0_pr - 0.6027_pr)*exp(-Tloc/(208_pr))) + (0.6027_pr*exp(-Tloc/(3922_pr)))&
              + exp(-(10180_pr)/Tloc)
    k(r95b) = getlindratecoeff(Tloc,k_0(FOr95b),k_inf(FOr95b),FC(FOr95b), M(mM95),Ploc)
    k(r96b) = (26873759.5767095_pr)*exp((-37588.976701611_pr)*R_T_inv + T_log *(-0.0861491731019821_pr))
    k(r97b) = (11481064735815_pr)*exp((-308637.694413429_pr)*R_T_inv + T_log *(-1.06182189105891_pr))
    k(r98b) = (6565900.08745883_pr)*exp((-58211.0209078635_pr)*R_T_inv + T_log *(-0.156136171362597_pr))
    k(r99b) = (8062127.74046857_pr)*exp((-37588.9766896408_pr)*R_T_inv + T_log *(-0.0861491709648683_pr))
    k(r100b) = (6270543.46080477_pr)*exp((-37588.9766505078_pr)*R_T_inv + T_log *(-0.0861491639703141_pr))
    k(r101b) = (106248.129477802_pr)*exp((-254795.565689549_pr)*R_T_inv + T_log *(0.518030880928361_pr))
    k(r102b) = (1594982998.5902_pr)*exp((-11182.4085518218_pr)*R_T_inv + T_log *(-0.319473091975016_pr))
    k(r103b) = (451479.264790068_pr)*exp((-300891.417407832_pr)*R_T_inv + T_log *(0.196458894259829_pr))
    k(r104b) = (0.0594534570494178_pr)*exp((-95953.7333171516_pr)*R_T_inv + T_log *(2.65876967705197_pr))
    k_0(FOr105b) = (3.86457479015067e+42_pr)*exp((-392741.50201074_pr)*R_T_inv + T_log *(-8.08365158176097_pr))
    k_inf(FOr105b) = (7.69505079828506e+23_pr)*exp((-383921.630049687_pr)*R_T_inv + T_log *(-2.23365158863998_pr))
    FC(FOr105b) = ((1.0_pr - 0.619_pr)*exp(-Tloc/(73.2_pr))) + (0.619_pr*exp(-Tloc/(1180_pr)))&
              + exp(-(9999_pr)/Tloc)
    k(r105b) = getlindratecoeff(Tloc,k_0(FOr105b),k_inf(FOr105b),FC(FOr105b), M(mM105),Ploc)
    k(r106b) = (193080444393.979_pr)*exp((-5583.88286148333_pr)*R_T_inv + T_log *(-0.870927806802739_pr))
    k(r107b) = (1076777092.20914_pr)*exp((-375092.841430875_pr)*R_T_inv + T_log *(0.134973359047523_pr))
    k(r108b) = (0.00697191919383231_pr)*exp((-96334.4473747459_pr)*R_T_inv + T_log *(2.77825569853129_pr))
    k(r109b) = (416.481339051414_pr)*exp((-47832.1171904263_pr)*R_T_inv + T_log *(1.36990199191127_pr))
    k(r110b) = (1.38867128153506_pr)*exp((-63395.8971107513_pr)*R_T_inv + T_log *(1.82333293673591_pr))
    k(r111b) = (669772.703777375_pr)*exp((-6595.93325996503_pr)*R_T_inv + T_log *(-0.865635818444764_pr))
    k(r112b) = (83498.3299169219_pr)*exp((-6595.93325565694_pr)*R_T_inv + T_log *(-0.86563581768805_pr))
    k(r113b) = (5767218.49050619_pr)*exp((-141027.01238183_pr)*R_T_inv + T_log *(0.163870563639617_pr))
    k(r114b) = (1508361.36160335_pr)*exp((-84851.9229240185_pr)*R_T_inv + T_log *(0.430722941090067_pr))
    k(r115b) = (1.23985141340484e-20_pr)*exp((-96363.2744596223_pr)*R_T_inv + T_log *(7.8415684670767_pr))
    k_0(FOr116b) = (4.95251324369475e+37_pr)*exp((-232692.458092495_pr)*R_T_inv + T_log *(-9.31433232832737_pr))
    k_inf(FOr116b) = (250760.205416341_pr)*exp((-186542.938212177_pr)*R_T_inv + T_log *(0.425667650335722_pr))
    FC(FOr116b) = ((1.0_pr - 0.7345_pr)*exp(-Tloc/(180_pr))) + (0.7345_pr*exp(-Tloc/(1035_pr)))&
              + exp(-(5417_pr)/Tloc)
    k(r116b) = getlindratecoeff(Tloc,k_0(FOr116b),k_inf(FOr116b),FC(FOr116b), M(mM116),Ploc)
    k(r117b) = (731577.001623313_pr)*exp((-69666.761628575_pr)*R_T_inv + T_log *(0.0941392857083723_pr))
    k(r119b) = (98.0159588954126_pr)*exp((-59644.661503801_pr)*R_T_inv + T_log *(1.33474456422871_pr))
    k(r120b) = (8784859901.51405_pr)*exp((-365664.731893228_pr)*R_T_inv + T_log *(0.258036516473117_pr))
    k_0(FOr122b) = (1.56243937309683e+22_pr)*exp((-453055.671434165_pr)*R_T_inv + T_log *(-3.06889270590257_pr))
    k_inf(FOr122b) = (638590524332406_pr)*exp((-449039.031620948_pr)*R_T_inv + T_log *(0.161107262122784_pr))
    FC(FOr122b) = ((1.0_pr - 0.578_pr)*exp(-Tloc/(122_pr))) + (0.578_pr*exp(-Tloc/(2535_pr)))&
              + exp(-(9365_pr)/Tloc)
    k(r122b) = getlindratecoeff(Tloc,k_0(FOr122b),k_inf(FOr122b),FC(FOr122b), M(mM122),Ploc)
    k(r124b) = (135323695.282269_pr)*exp((-259809.236020762_pr)*R_T_inv + T_log *(-0.176788783905835_pr))
    k(r129b) = (69.5718754103485_pr)*exp((-50493.18708114_pr)*R_T_inv + T_log *(1.20962932832788_pr))

    return
  end subroutine get_rate_coefficients

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

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

    w(r1f) = k(r1f) * c(sH2) * c(sO)
    w(r2f) = k(r2f) * c(sHO2) * c(sO)
    w(r3f) = k(r3f) * cqss(sqssCH) * c(sO)
    w(r4f) = k(r4f) * cqss(sqssCH2) * c(sO)
    w(r5f) = k(r5f) * cqss(sqssCH2GSG) * c(sO)
    w(r6f) = k(r6f) * cqss(sqssCH2GSG) * c(sO)
    w(r7f) = k(r7f) * c(sCH3) * c(sO)
    w(r8f) = k(r8f) * c(sCH4) * c(sO)
    w(r9f) = k(r9f) * c(sCO) * c(sO)
    w(r10f) = k(r10f) * cqss(sqssHCO) * c(sO)
    w(r11f) = k(r11f) * cqss(sqssHCO) * c(sO)
    w(r12f) = k(r12f) * c(sCH2O) * c(sO)
    w(r13f) = k(r13f) * cqss(sqssCH2OH) * c(sO)
    w(r14f) = k(r14f) * cqss(sqssCH3OH) * c(sO)
    w(r15f) = k(r15f) * cqss(sqssC2H2) * c(sO)
    w(r16f) = k(r16f) * c(sC2H4) * c(sO)
    w(r17f) = k(r17f) * cqss(sqssC2H5) * c(sO)
    w(r18f) = k(r18f) * c(sC2H6) * c(sO)
    w(r19f) = k(r19f) * c(sCO) * c(sO2)
    w(r20f) = k(r20f) * c(sCH2O) * c(sO2)
    w(r21f) = k(r21f) * c(sH) * c(sO2) * m(mM21)
    w(r22f) = k(r22f) * c(sH) * c(sO2)
    w(r23f) = k(r23f) * c(sH) * c(sH2O) * c(sO2)
    w(r24f) = k(r24f) * c(sH) * c(sN2) * c(sO2)
    w(r25f) = k(r25f) * c(sH) * c(sO2)
    w(r26f) = k(r26f) * c(sH)**2.0_pr * m(mM26)
    w(r27f) = k(r27f) * c(sH)**2.0_pr * c(sH2O)
    w(r28f) = k(r28f) * c(sH) * c(sOH) * m(mM28)
    w(r29f) = k(r29f) * c(sH) * c(sHO2)
    w(r30f) = k(r30f) * c(sH) * c(sHO2)
    w(r31f) = k(r31f) * c(sH) * c(sHO2)
    w(r32f) = k(r32f) * c(sH) * c(sH2O2)
    w(r33f) = k(r33f) * cqss(sqssCH2) * c(sH)
    w(r34f) = k(r34f) * cqss(sqssCH2GSG) * c(sH)
    w(r35f) = k(r35f) * c(sCH3) * c(sH)
    w(r36f) = k(r36f) * c(sCH4) * c(sH)
    w(r37f) = k(r37f) * c(sH) * cqss(sqssHCO)
    w(r38f) = k(r38f) * c(sH) * cqss(sqssHCO)
    w(r39f) = k(r39f) * c(sCH2O) * c(sH)
    w(r40f) = k(r40f) * c(sCH2O) * c(sH)
    w(r41f) = k(r41f) * c(sCH2O) * c(sH)
    w(r42f) = k(r42f) * cqss(sqssCH2OH) * c(sH)
    w(r43f) = k(r43f) * cqss(sqssCH2OH) * c(sH)
    w(r44f) = k(r44f) * cqss(sqssCH2OH) * c(sH)
    w(r45f) = k(r45f) * cqss(sqssCH3O) * c(sH)
    w(r46f) = k(r46f) * cqss(sqssCH3O) * c(sH)
    w(r47f) = k(r47f) * cqss(sqssCH3O) * c(sH)
    w(r48f) = k(r48f) * cqss(sqssCH3OH) * c(sH)
    w(r49f) = k(r49f) * cqss(sqssCH3OH) * c(sH)
    w(r50f) = k(r50f) * c(sC2H4) * c(sH)
    w(r51f) = k(r51f) * cqss(sqssC2H5) * c(sH)
    w(r52f) = k(r52f) * cqss(sqssC2H5) * c(sH)
    w(r53f) = k(r53f) * c(sC2H6) * c(sH)
    w(r54f) = k(r54f) * c(sCO) * c(sH2)
    w(r55f) = k(r55f) * c(sH2) * c(sOH)
    w(r56f) = k(r56f) * c(sOH)**2.0_pr
    w(r57f) = k(r57f) * c(sOH)**2.0_pr
    w(r58f) = k(r58f) * c(sHO2) * c(sOH)
    w(r59f) = k(r59f) * c(sH2O2) * c(sOH)
    w(r60f) = k(r60f) * cqss(sqssCH) * c(sOH)
    w(r61f) = k(r61f) * cqss(sqssCH2) * c(sOH)
    w(r62f) = k(r62f) * cqss(sqssCH2) * c(sOH)
    w(r63f) = k(r63f) * cqss(sqssCH2GSG) * c(sOH)
    w(r64f) = k(r64f) * c(sCH3) * c(sOH)
    w(r65f) = k(r65f) * c(sCH3) * c(sOH)
    w(r66f) = k(r66f) * c(sCH3) * c(sOH)
    w(r67f) = k(r67f) * c(sCH4) * c(sOH)
    w(r68f) = k(r68f) * c(sCO) * c(sOH)
    w(r69f) = k(r69f) * cqss(sqssHCO) * c(sOH)
    w(r70f) = k(r70f) * c(sCH2O) * c(sOH)
    w(r71f) = k(r71f) * cqss(sqssCH2OH) * c(sOH)
    w(r72f) = k(r72f) * cqss(sqssCH3OH) * c(sOH)
    w(r73f) = k(r73f) * cqss(sqssCH3OH) * c(sOH)
    w(r74f) = k(r74f) * cqss(sqssC2H2) * c(sOH)
    w(r75f) = k(r75f) * c(sC2H6) * c(sOH)
    w(r76f) = k(r76f) * c(sHO2)**2.0_pr
    w(r77f) = k(r77f) * cqss(sqssCH2) * c(sHO2)
    w(r78f) = k(r78f) * c(sCH3) * c(sHO2)
    w(r79f) = k(r79f) * c(sCH3) * c(sHO2)
    w(r80f) = k(r80f) * c(sCO) * c(sHO2)
    w(r81f) = k(r81f) * c(sCH2O) * c(sHO2)
    w(r82f) = k(r82f) * cqss(sqssCH) * c(sO2)
    w(r83f) = k(r83f) * cqss(sqssCH) * c(sH2)
    w(r84f) = k(r84f) * cqss(sqssCH) * c(sH2O)
    w(r85f) = k(r85f) * cqss(sqssCH) * c(sCH4)
    w(r86f) = k(r86f) * cqss(sqssCH) * c(sCO2)
    w(r87f) = k(r87f) * cqss(sqssCH2) * c(sO2)
    w(r88f) = k(r88f) * cqss(sqssCH2) * c(sH2)
    w(r89f) = k(r89f) * cqss(sqssCH2) * c(sCH3)
    w(r90f) = k(r90f) * cqss(sqssCH2) * c(sCH4)
    w(r91f) = k(r91f) * cqss(sqssCH2GSG) * c(sN2)
    w(r92f) = k(r92f) * cqss(sqssCH2GSG) * c(sO2)
    w(r93f) = k(r93f) * cqss(sqssCH2GSG) * c(sO2)
    w(r94f) = k(r94f) * cqss(sqssCH2GSG) * c(sH2)
    w(r95f) = k(r95f) * cqss(sqssCH2GSG) * c(sH2O)
    w(r96f) = k(r96f) * cqss(sqssCH2GSG) * c(sH2O)
    w(r97f) = k(r97f) * cqss(sqssCH2GSG) * c(sCH3)
    w(r98f) = k(r98f) * cqss(sqssCH2GSG) * c(sCH4)
    w(r99f) = k(r99f) * cqss(sqssCH2GSG) * c(sCO)
    w(r100f) = k(r100f) * cqss(sqssCH2GSG) * c(sCO2)
    w(r101f) = k(r101f) * cqss(sqssCH2GSG) * c(sCO2)
    w(r102f) = k(r102f) * c(sCH3) * c(sO2)
    w(r103f) = k(r103f) * c(sCH3) * c(sO2)
    w(r104f) = k(r104f) * c(sCH3) * c(sH2O2)
    w(r105f) = k(r105f) * c(sCH3)**2.0_pr
    w(r106f) = k(r106f) * c(sCH3)**2.0_pr
    w(r107f) = k(r107f) * c(sCH3) * cqss(sqssHCO)
    w(r108f) = k(r108f) * c(sCH2O) * c(sCH3)
    w(r109f) = k(r109f) * c(sCH3) * cqss(sqssCH3OH)
    w(r110f) = k(r110f) * c(sC2H6) * c(sCH3)
    w(r111f) = k(r111f) * c(sH2O) * cqss(sqssHCO)
    w(r112f) = k(r112f) * cqss(sqssHCO) * m(mM112)
    w(r113f) = k(r113f) * cqss(sqssHCO) * c(sO2)
    w(r114f) = k(r114f) * cqss(sqssCH2OH) * c(sO2)
    w(r115f) = k(r115f) * cqss(sqssCH3O) * c(sO2)
    w(r116f) = k(r116f) * c(sC2H4)
    w(r117f) = k(r117f) * cqss(sqssC2H5) * c(sO2)
    w(r118f) = k(r118f) * c(sCH3) * c(sO)
    w(r119f) = k(r119f) * c(sC2H4) * c(sO)
    w(r120f) = k(r120f) * c(sHO2) * c(sOH)
    w(r121f) = k(r121f) * c(sCH3) * c(sOH)
    w(r122f) = k(r122f) * cqss(sqssCH) * c(sH2)
    w(r123f) = k(r123f) * cqss(sqssCH2) * c(sO2)
    w(r124f) = k(r124f) * cqss(sqssCH2) * c(sO2)
    w(r125f) = k(r125f) * cqss(sqssCH2GSG) * c(sH2O)
    w(r126f) = k(r126f) * cqss(sqssCH2CHO) * c(sO)
    w(r127f) = k(r127f) * cqss(sqssCH2CHO) * c(sO2)
    w(r128f) = k(r128f) * cqss(sqssCH2CHO) * c(sO2)
    w(r129f) = k(r129f) * cqss(sqssCH2CHO) * c(sH)
    w(r1b) = k(r1b) * c(sH) * c(sOH)
    w(r2b) = k(r2b) * c(sO2) * c(sOH)
    w(r3b) = k(r3b) * c(sCO) * c(sH)
    w(r4b) = k(r4b) * c(sH) * cqss(sqssHCO)
    w(r5b) = k(r5b) * c(sCO) * c(sH2)
    w(r6b) = k(r6b) * c(sH) * cqss(sqssHCO)
    w(r7b) = k(r7b) * c(sCH2O) * c(sH)
    w(r8b) = k(r8b) * c(sCH3) * c(sOH)
    w(r9b) = k(r9b) * c(sCO2)
    w(r10b) = k(r10b) * c(sCO) * c(sOH)
    w(r11b) = k(r11b) * c(sCO2) * c(sH)
    w(r12b) = k(r12b) * cqss(sqssHCO) * c(sOH)
    w(r13b) = k(r13b) * c(sCH2O) * c(sOH)
    w(r14b) = k(r14b) * cqss(sqssCH2OH) * c(sOH)
    w(r15b) = k(r15b) * cqss(sqssCH2) * c(sCO)
    w(r16b) = k(r16b) * c(sCH3) * cqss(sqssHCO)
    w(r17b) = k(r17b) * c(sCH2O) * c(sCH3)
    w(r18b) = k(r18b) * cqss(sqssC2H5) * c(sOH)
    w(r19b) = k(r19b) * c(sCO2) * c(sO)
    w(r20b) = k(r20b) * cqss(sqssHCO) * c(sHO2)
    w(r21b) = k(r21b) * c(sHO2) * m(mM21)
    w(r22b) = k(r22b) * c(sHO2) * c(sO2)
    w(r23b) = k(r23b) * c(sH2O) * c(sHO2)
    w(r24b) = k(r24b) * c(sHO2) * c(sN2)
    w(r25b) = k(r25b) * c(sO) * c(sOH)
    w(r26b) = k(r26b) * c(sH2) * m(mM26)
    w(r27b) = k(r27b) * c(sH2) * c(sH2O)
    w(r28b) = k(r28b) * c(sH2O) * m(mM28)
    w(r29b) = k(r29b) * c(sH2O) * c(sO)
    w(r30b) = k(r30b) * c(sH2) * c(sO2)
    w(r31b) = k(r31b) * c(sOH)**2.0_pr
    w(r32b) = k(r32b) * c(sH2) * c(sHO2)
    w(r33b) = k(r33b) * c(sCH3)
    w(r34b) = k(r34b) * cqss(sqssCH) * c(sH2)
    w(r35b) = k(r35b) * c(sCH4)
    w(r36b) = k(r36b) * c(sCH3) * c(sH2)
    w(r37b) = k(r37b) * c(sCH2O)
    w(r38b) = k(r38b) * c(sCO) * c(sH2)
    w(r39b) = k(r39b) * cqss(sqssCH2OH)
    w(r40b) = k(r40b) * cqss(sqssCH3O)
    w(r41b) = k(r41b) * c(sH2) * cqss(sqssHCO)
    w(r42b) = k(r42b) * c(sCH2O) * c(sH2)
    w(r43b) = k(r43b) * c(sCH3) * c(sOH)
    w(r44b) = k(r44b) * cqss(sqssCH2GSG) * c(sH2O)
    w(r45b) = k(r45b) * c(sCH2O) * c(sH2)
    w(r46b) = k(r46b) * c(sCH3) * c(sOH)
    w(r47b) = k(r47b) * cqss(sqssCH2GSG) * c(sH2O)
    w(r48b) = k(r48b) * cqss(sqssCH2OH) * c(sH2)
    w(r49b) = k(r49b) * cqss(sqssCH3O) * c(sH2)
    w(r50b) = k(r50b) * cqss(sqssC2H5)
    w(r51b) = k(r51b) * c(sC2H6)
    w(r52b) = k(r52b) * c(sC2H4) * c(sH2)
    w(r53b) = k(r53b) * cqss(sqssC2H5) * c(sH2)
    w(r54b) = k(r54b) * c(sCH2O)
    w(r55b) = k(r55b) * c(sH) * c(sH2O)
    w(r56b) = k(r56b) * c(sH2O2)
    w(r57b) = k(r57b) * c(sH2O) * c(sO)
    w(r58b) = k(r58b) * c(sH2O) * c(sO2)
    w(r59b) = k(r59b) * c(sH2O) * c(sHO2)
    w(r60b) = k(r60b) * c(sH) * cqss(sqssHCO)
    w(r61b) = k(r61b) * c(sCH2O) * c(sH)
    w(r62b) = k(r62b) * cqss(sqssCH) * c(sH2O)
    w(r63b) = k(r63b) * c(sCH2O) * c(sH)
    w(r64b) = k(r64b) * cqss(sqssCH3OH)
    w(r65b) = k(r65b) * cqss(sqssCH2) * c(sH2O)
    w(r66b) = k(r66b) * cqss(sqssCH2GSG) * c(sH2O)
    w(r67b) = k(r67b) * c(sCH3) * c(sH2O)
    w(r68b) = k(r68b) * c(sCO2) * c(sH)
    w(r69b) = k(r69b) * c(sCO) * c(sH2O)
    w(r70b) = k(r70b) * c(sH2O) * cqss(sqssHCO)
    w(r71b) = k(r71b) * c(sCH2O) * c(sH2O)
    w(r72b) = k(r72b) * cqss(sqssCH2OH) * c(sH2O)
    w(r73b) = k(r73b) * cqss(sqssCH3O) * c(sH2O)
    w(r74b) = k(r74b) * c(sCH3) * c(sCO)
    w(r75b) = k(r75b) * cqss(sqssC2H5) * c(sH2O)
    w(r76b) = k(r76b) * c(sH2O2) * c(sO2)
    w(r77b) = k(r77b) * c(sCH2O) * c(sOH)
    w(r78b) = k(r78b) * c(sCH4) * c(sO2)
    w(r79b) = k(r79b) * cqss(sqssCH3O) * c(sOH)
    w(r80b) = k(r80b) * c(sCO2) * c(sOH)
    w(r81b) = k(r81b) * c(sH2O2) * cqss(sqssHCO)
    w(r82b) = k(r82b) * cqss(sqssHCO) * c(sO)
    w(r83b) = k(r83b) * cqss(sqssCH2) * c(sH)
    w(r84b) = k(r84b) * c(sCH2O) * c(sH)
    w(r85b) = k(r85b) * c(sC2H4) * c(sH)
    w(r86b) = k(r86b) * c(sCO) * cqss(sqssHCO)
    w(r88b) = k(r88b) * c(sCH3) * c(sH)
    w(r89b) = k(r89b) * c(sC2H4) * c(sH)
    w(r90b) = k(r90b) * c(sCH3)**2.0_pr
    w(r91b) = k(r91b) * cqss(sqssCH2) * c(sN2)
    w(r92b) = k(r92b) * c(sCO) * c(sH) * c(sOH)
    w(r93b) = k(r93b) * c(sCO) * c(sH2O)
    w(r94b) = k(r94b) * c(sCH3) * c(sH)
    w(r95b) = k(r95b) * cqss(sqssCH3OH)
    w(r96b) = k(r96b) * cqss(sqssCH2) * c(sH2O)
    w(r97b) = k(r97b) * c(sC2H4) * c(sH)
    w(r98b) = k(r98b) * c(sCH3)**2.0_pr
    w(r99b) = k(r99b) * cqss(sqssCH2) * c(sCO)
    w(r100b) = k(r100b) * cqss(sqssCH2) * c(sCO2)
    w(r101b) = k(r101b) * c(sCH2O) * c(sCO)
    w(r102b) = k(r102b) * cqss(sqssCH3O) * c(sO)
    w(r103b) = k(r103b) * c(sCH2O) * c(sOH)
    w(r104b) = k(r104b) * c(sCH4) * c(sHO2)
    w(r105b) = k(r105b) * c(sC2H6)
    w(r106b) = k(r106b) * cqss(sqssC2H5) * c(sH)
    w(r107b) = k(r107b) * c(sCH4) * c(sCO)
    w(r108b) = k(r108b) * c(sCH4) * cqss(sqssHCO)
    w(r109b) = k(r109b) * cqss(sqssCH3O) * c(sCH4)
    w(r110b) = k(r110b) * cqss(sqssC2H5) * c(sCH4)
    w(r111b) = k(r111b) * c(sCO) * c(sH) * c(sH2O)
    w(r112b) = k(r112b) * c(sCO) * c(sH) * m(mM112)
    w(r113b) = k(r113b) * c(sCO) * c(sHO2)
    w(r114b) = k(r114b) * c(sCH2O) * c(sHO2)
    w(r115b) = k(r115b) * c(sCH2O) * c(sHO2)
    w(r116b) = k(r116b) * cqss(sqssC2H2) * c(sH2)
    w(r117b) = k(r117b) * c(sC2H4) * c(sHO2)
    w(r119b) = k(r119b) * cqss(sqssCH2CHO) * c(sH)
    w(r120b) = k(r120b) * c(sH2O) * c(sO2)
    w(r122b) = k(r122b) * c(sCH3)
    w(r124b) = k(r124b) * c(sCH2O) * c(sO)
    w(r129b) = k(r129b) * c(sCH3) * cqss(sqssHCO)

    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(r1f) &
         + w(r1b) &
         + w(r5f) &
         - w(r5b) &
         + w(r26f) &
         - w(r26b) &
         + w(r27f) &
         - w(r27b) &
         + w(r30f) &
         - w(r30b) &
         + w(r32f) &
         - w(r32b) &
         + w(r34f) &
         - w(r34b) &
         + w(r36f) &
         - w(r36b) &
         + w(r38f) &
         - w(r38b) &
         + w(r41f) &
         - w(r41b) &
         + w(r42f) &
         - w(r42b) &
         + w(r45f) &
         - w(r45b) &
         + w(r48f) &
         - w(r48b) &
         + w(r49f) &
         - w(r49b) &
         + w(r52f) &
         - w(r52b) &
         + w(r53f) &
         - w(r53b) &
         - w(r54f) &
         + w(r54b) &
         - w(r55f) &
         + w(r55b) &
         - w(r83f) &
         + w(r83b) &
         - w(r88f) &
         + w(r88b) &
         - w(r94f) &
         + w(r94b) &
         + w(r116f) &
         - w(r116b) &
         + w(r118f) &
         + w(r121f) &
         - w(r122f) &
         + w(r122b) &
         + w(r125f)

    cdot(sH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r3f) &
         - w(r3b) &
         + w(r4f) &
         - w(r4b) &
         + w(r6f) &
         - w(r6b) &
         + w(r7f) &
         - w(r7b) &
         + w(r11f) &
         - w(r11b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) &
         - w(r25f) &
         + w(r25b) &
         -2.0_pr * w(r26f) &
         + 2.0_pr * w(r26b) &
         -2.0_pr * w(r27f) &
         + 2.0_pr * w(r27b) &
         - 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(r35f) &
         + w(r35b) &
         - w(r36f) &
         + w(r36b) &
         - w(r37f) &
         + 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(r55f) &
         - w(r55b) &
         + w(r60f) &
         - w(r60b) &
         + w(r61f) &
         - w(r61b) &
         + w(r63f) &
         - w(r63b) &
         + w(r68f) &
         - w(r68b) &
         + w(r83f) &
         - w(r83b) &
         + w(r84f) &
         - w(r84b) &
         + w(r85f) &
         - w(r85b) &
         + w(r87f) &
         + w(r88f) &
         - w(r88b) &
         + w(r89f) &
         - w(r89b) &
         + w(r92f) &
         - w(r92b) &
         + w(r94f) &
         - w(r94b) &
         + w(r97f) &
         - w(r97b) &
         + w(r106f) &
         - w(r106b) &
         + w(r111f) &
         - w(r111b) &
         + w(r112f) &
         - w(r112b) &
         + w(r118f) &
         + w(r119f) &
         - w(r119b) &
         + 2.0_pr * w(r123f) &
         + w(r126f) &
         - w(r129f) &
         + w(r129b)

    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(r25f) &
         - w(r25b) &
         + w(r29f) &
         - w(r29b) &
         + w(r57f) &
         - w(r57b) &
         + w(r82f) &
         - w(r82b) &
         + w(r102f) &
         - w(r102b) &
         - w(r118f) &
         - w(r119f) &
         + w(r119b) &
         + w(r124f) &
         - w(r124b) &
         - w(r126f)

    cdot(sO2) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         - 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(r30f) &
         - w(r30b) &
         + w(r58f) &
         - w(r58b) &
         + w(r76f) &
         - w(r76b) &
         + w(r78f) &
         - w(r78b) &
         - w(r82f) &
         + w(r82b) &
         - w(r87f) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r113f) &
         + w(r113b) &
         - w(r114f) &
         + w(r114b) &
         - w(r115f) &
         + w(r115b) &
         - w(r117f) &
         + w(r117b) &
         + w(r120f) &
         - w(r120b) &
         - w(r123f) &
         - w(r124f) &
         + w(r124b) &
         - w(r127f) &
         - w(r128f)

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r8f) &
         - w(r8b) &
         + w(r10f) &
         - w(r10b) &
         + w(r12f) &
         - w(r12b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         + w(r18f) &
         - w(r18b) &
         + w(r25f) &
         - w(r25b) &
         - w(r28f) &
         + w(r28b) &
         + 2.0_pr * w(r31f) &
         - 2.0_pr * w(r31b) &
         + w(r43f) &
         - w(r43b) &
         + w(r46f) &
         - w(r46b) &
         - w(r55f) &
         + w(r55b) &
         -2.0_pr * w(r56f) &
         + 2.0_pr * w(r56b) &
         -2.0_pr * w(r57f) &
         + 2.0_pr * 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(r72f) &
         + w(r72b) &
         - w(r73f) &
         + w(r73b) &
         - w(r74f) &
         + w(r74b) &
         - w(r75f) &
         + w(r75b) &
         + w(r77f) &
         - w(r77b) &
         + w(r79f) &
         - w(r79b) &
         + w(r80f) &
         - w(r80b) &
         + w(r87f) &
         + w(r92f) &
         - w(r92b) &
         + w(r103f) &
         - w(r103b) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r127f) &
         + w(r128f)

    cdot(sH2O) = 0.0_pr &
         + w(r28f) &
         - w(r28b) &
         + w(r29f) &
         - w(r29b) &
         + w(r44f) &
         - w(r44b) &
         + w(r47f) &
         - w(r47b) &
         + w(r55f) &
         - w(r55b) &
         + w(r57f) &
         - w(r57b) &
         + w(r58f) &
         - w(r58b) &
         + w(r59f) &
         - w(r59b) &
         + w(r62f) &
         - w(r62b) &
         + w(r65f) &
         - w(r65b) &
         + w(r66f) &
         - w(r66b) &
         + w(r67f) &
         - w(r67b) &
         + w(r69f) &
         - w(r69b) &
         + w(r70f) &
         - w(r70b) &
         + w(r71f) &
         - w(r71b) &
         + w(r72f) &
         - w(r72b) &
         + w(r73f) &
         - w(r73b) &
         + w(r75f) &
         - w(r75b) &
         - w(r84f) &
         + w(r84b) &
         + w(r93f) &
         - w(r93b) &
         - w(r95f) &
         + w(r95b) &
         + w(r120f) &
         - w(r120b) &
         - w(r125f)

    cdot(sHO2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         + w(r20f) &
         - w(r20b) &
         + w(r21f) &
         - w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r23f) &
         - w(r23b) &
         + w(r24f) &
         - w(r24b) &
         - w(r29f) &
         + w(r29b) &
         - w(r30f) &
         + w(r30b) &
         - w(r31f) &
         + w(r31b) &
         + w(r32f) &
         - w(r32b) &
         - w(r58f) &
         + w(r58b) &
         + w(r59f) &
         - w(r59b) &
         -2.0_pr * w(r76f) &
         + 2.0_pr * w(r76b) &
         - w(r77f) &
         + w(r77b) &
         - w(r78f) &
         + w(r78b) &
         - w(r79f) &
         + w(r79b) &
         - w(r80f) &
         + w(r80b) &
         - w(r81f) &
         + w(r81b) &
         + w(r104f) &
         - w(r104b) &
         + w(r113f) &
         - w(r113b) &
         + w(r114f) &
         - w(r114b) &
         + w(r115f) &
         - w(r115b) &
         + w(r117f) &
         - w(r117b) &
         - w(r120f) &
         + w(r120b)

    cdot(sH2O2) = 0.0_pr &
         - w(r32f) &
         + w(r32b) &
         + w(r56f) &
         - w(r56b) &
         - w(r59f) &
         + w(r59b) &
         + w(r76f) &
         - w(r76b) &
         + w(r81f) &
         - w(r81b) &
         - w(r104f) &
         + w(r104b)

    cdot(sCH3) = 0.0_pr &
         - w(r7f) &
         + w(r7b) &
         + w(r8f) &
         - w(r8b) &
         + w(r16f) &
         - w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r33f) &
         - w(r33b) &
         - w(r35f) &
         + w(r35b) &
         + w(r36f) &
         - w(r36b) &
         + w(r43f) &
         - w(r43b) &
         + w(r46f) &
         - w(r46b) &
         - w(r64f) &
         + w(r64b) &
         - w(r65f) &
         + w(r65b) &
         - w(r66f) &
         + w(r66b) &
         + w(r67f) &
         - w(r67b) &
         + w(r74f) &
         - w(r74b) &
         - w(r78f) &
         + w(r78b) &
         - w(r79f) &
         + w(r79b) &
         + w(r88f) &
         - w(r88b) &
         - w(r89f) &
         + w(r89b) &
         + 2.0_pr * w(r90f) &
         - 2.0_pr * w(r90b) &
         + w(r94f) &
         - w(r94b) &
         - w(r97f) &
         + w(r97b) &
         + 2.0_pr * w(r98f) &
         - 2.0_pr * w(r98b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         -2.0_pr * w(r105f) &
         + 2.0_pr * w(r105b) &
         -2.0_pr * w(r106f) &
         + 2.0_pr * w(r106b) &
         - w(r107f) &
         + w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r109f) &
         + w(r109b) &
         - w(r110f) &
         + w(r110b) &
         - w(r118f) &
         - w(r121f) &
         + w(r122f) &
         - w(r122b) &
         + w(r129f) &
         - w(r129b)

    cdot(sCH4) = 0.0_pr &
         - w(r8f) &
         + w(r8b) &
         + w(r35f) &
         - w(r35b) &
         - w(r36f) &
         + w(r36b) &
         - w(r67f) &
         + w(r67b) &
         + w(r78f) &
         - w(r78b) &
         - w(r85f) &
         + w(r85b) &
         - w(r90f) &
         + w(r90b) &
         - w(r98f) &
         + w(r98b) &
         + w(r104f) &
         - w(r104b) &
         + w(r107f) &
         - w(r107b) &
         + w(r108f) &
         - w(r108b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b)

    cdot(sCO) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         + w(r5f) &
         - w(r5b) &
         - w(r9f) &
         + w(r9b) &
         + w(r10f) &
         - w(r10b) &
         + w(r15f) &
         - w(r15b) &
         - w(r19f) &
         + w(r19b) &
         + w(r38f) &
         - w(r38b) &
         - w(r54f) &
         + w(r54b) &
         - w(r68f) &
         + w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r74f) &
         - w(r74b) &
         - w(r80f) &
         + w(r80b) &
         + w(r86f) &
         - w(r86b) &
         + w(r87f) &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         + w(r101f) &
         - w(r101b) &
         + w(r107f) &
         - w(r107b) &
         + w(r111f) &
         - w(r111b) &
         + w(r112f) &
         - w(r112b) &
         + w(r113f) &
         - w(r113b) &
         + w(r118f) &
         + w(r127f)

    cdot(sCO2) = 0.0_pr &
         + w(r9f) &
         - w(r9b) &
         + w(r11f) &
         - w(r11b) &
         + w(r19f) &
         - w(r19b) &
         + w(r68f) &
         - w(r68b) &
         + w(r80f) &
         - w(r80b) &
         - w(r86f) &
         + w(r86b) &
         - w(r101f) &
         + w(r101b) &
         + w(r123f) &
         + w(r126f)

    cdot(sCH2O) = 0.0_pr &
         + w(r7f) &
         - w(r7b) &
         - w(r12f) &
         + w(r12b) &
         + w(r13f) &
         - w(r13b) &
         + w(r17f) &
         - w(r17b) &
         - w(r20f) &
         + w(r20b) &
         + w(r37f) &
         - w(r37b) &
         - w(r39f) &
         + w(r39b) &
         - w(r40f) &
         + w(r40b) &
         - w(r41f) &
         + w(r41b) &
         + w(r42f) &
         - w(r42b) &
         + w(r45f) &
         - w(r45b) &
         + w(r54f) &
         - w(r54b) &
         + w(r61f) &
         - w(r61b) &
         + w(r63f) &
         - w(r63b) &
         - w(r70f) &
         + w(r70b) &
         + w(r71f) &
         - w(r71b) &
         + w(r77f) &
         - w(r77b) &
         - w(r81f) &
         + w(r81b) &
         + w(r84f) &
         - w(r84b) &
         + w(r101f) &
         - w(r101b) &
         + w(r103f) &
         - w(r103b) &
         - w(r108f) &
         + w(r108b) &
         + w(r114f) &
         - w(r114b) &
         + w(r115f) &
         - w(r115b) &
         + w(r121f) &
         + w(r124f) &
         - w(r124b) &
         + w(r125f) &
         + w(r127f)

    cdot(sC2H4) = 0.0_pr &
         - w(r16f) &
         + w(r16b) &
         - w(r50f) &
         + w(r50b) &
         + w(r52f) &
         - w(r52b) &
         + w(r85f) &
         - w(r85b) &
         + w(r89f) &
         - w(r89b) &
         + w(r97f) &
         - w(r97b) &
         - w(r116f) &
         + w(r116b) &
         + w(r117f) &
         - w(r117b) &
         - w(r119f) &
         + w(r119b)

    cdot(sC2H6) = 0.0_pr &
         - w(r18f) &
         + w(r18b) &
         + w(r51f) &
         - w(r51b) &
         - w(r53f) &
         + w(r53b) &
         - w(r75f) &
         + w(r75b) &
         + w(r105f) &
         - w(r105b) &
         - w(r110f) &
         + w(r110b)

    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) :: 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_C2H2
    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_CH3O
    real(pr) :: CH2GSG_CH3OH
    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_CH2CHO
    real(pr) :: CH2OH_ct
    real(pr) :: CH2OH_num
    real(pr) :: CH2OH_denom
    real(pr) :: CH2OH_CH2GSG
    real(pr) :: CH2OH_CH3OH
    real(pr) :: CH3O_ct
    real(pr) :: CH3O_num
    real(pr) :: CH3O_denom
    real(pr) :: CH3O_CH2GSG
    real(pr) :: CH3O_CH3OH
    real(pr) :: CH3OH_ct
    real(pr) :: CH3OH_num
    real(pr) :: CH3OH_denom
    real(pr) :: CH3OH_CH2GSG
    real(pr) :: CH3OH_CH2OH
    real(pr) :: CH3OH_CH3O
    real(pr) :: C2H2_ct
    real(pr) :: C2H2_num
    real(pr) :: C2H2_denom
    real(pr) :: C2H2_CH2
    real(pr) :: C2H5_ct
    real(pr) :: C2H5_num
    real(pr) :: C2H5_denom
    real(pr) :: CH2CHO_ct
    real(pr) :: CH2CHO_num
    real(pr) :: CH2CHO_denom
    real(pr) :: CH2CHO_HCO

    real(pr) :: a_j_a
    real(pr) :: a_k_a
    real(pr) :: a_l_a
    real(pr) :: a_m_a
    real(pr) :: a_n_a
    real(pr) :: a_q_a
    real(pr) :: a_s_a
    real(pr) :: a_t_a
    real(pr) :: a_u_a
    real(pr) :: a_v_a
    real(pr) :: a_w_a
    real(pr) :: a_z_a
    real(pr) :: a_ab_a
    real(pr) :: a_ac_a
    real(pr) :: a_ad_a
    real(pr) :: a_ae_a
    real(pr) :: a_af_a
    real(pr) :: a_ai_a
    real(pr) :: a_ak_a
    real(pr) :: a_al_a
    real(pr) :: a_am_a
    real(pr) :: a_an_a
    real(pr) :: a_ao_a
    real(pr) :: a_ar_a
    real(pr) :: a_at_a
    real(pr) :: a_au_a
    real(pr) :: a_av_a
    real(pr) :: a_aw_a
    real(pr) :: a_ax_a
    real(pr) :: a_ba_a
    real(pr) :: a_bd_a
    real(pr) :: a_bs_a
    real(pr) :: a_bu_a
    real(pr) :: a_bv_a
    real(pr) :: a_bw_a
    real(pr) :: a_bx_a
    real(pr) :: a_by_a
    real(pr) :: a_cb_a
    real(pr) :: a_u_b
    real(pr) :: a_v_b
    real(pr) :: a_ad_b
    real(pr) :: a_ae_b
    real(pr) :: a_am_b
    real(pr) :: a_an_b
    real(pr) :: a_av_b
    real(pr) :: a_aw_b
    real(pr) :: a_bw_b
    real(pr) :: a_bx_b
    real(pr) :: a_af_c
    real(pr) :: a_ag_c
    real(pr) :: a_ah_c
    real(pr) :: a_ai_c
    real(pr) :: a_ao_c
    real(pr) :: a_ap_c
    real(pr) :: a_aq_c
    real(pr) :: a_ar_c
    real(pr) :: a_ax_c
    real(pr) :: a_ay_c
    real(pr) :: a_bz_c
    real(pr) :: a_ba_c
    real(pr) :: a_bg_c
    real(pr) :: a_bh_c
    real(pr) :: a_bi_c
    real(pr) :: a_bj_c
    real(pr) :: a_br_c
    real(pr) :: a_bs_c
    real(pr) :: a_by_c
    real(pr) :: a_cz_c
    real(pr) :: a_ca_c
    real(pr) :: a_cb_c
    real(pr) :: a_bj_f
    real(pr) :: a_bs_f
    real(pr) :: a_cb_f
    real(pr) :: a_cb_h
    real(pr) :: A_A_A
    real(pr) :: A_B_A
    real(pr) :: A_C_A
    real(pr) :: A_D_A
    real(pr) :: A_E_A
    real(pr) :: A_F_A
    real(pr) :: A_G_A
    real(pr) :: A_H_A
    real(pr) :: A_B_B
    real(pr) :: A_C_B
    real(pr) :: A_D_B
    real(pr) :: A_E_B
    real(pr) :: A_F_B
    real(pr) :: A_G_B
    real(pr) :: A_H_B
    real(pr) :: A_C_C
    real(pr) :: A_D_C
    real(pr) :: A_E_C
    real(pr) :: A_F_C
    real(pr) :: A_G_C
    real(pr) :: A_H_C
    real(pr) :: A_D_D
    real(pr) :: A_E_D
    real(pr) :: A_F_D
    real(pr) :: A_G_D
    real(pr) :: A_H_D
    real(pr) :: A_E_E
    real(pr) :: A_F_E
    real(pr) :: A_G_E
    real(pr) :: A_H_E
    real(pr) :: A_F_F
    real(pr) :: A_G_F
    real(pr) :: A_H_F
    real(pr) :: A_G_G
    real(pr) :: A_H_G
    real(pr) :: A_H_H

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r10f)* c(sO) &
             + k(r11f)* c(sO) &
             + k(r37f)* c(sH) &
             + k(r38f)* c(sH) &
             + k(r69f)* c(sOH) &
             + k(r107f)* c(sCH3) &
             + k(r111f)* c(sH2O) &
             + k(r112f)* m(mM112) &
             + k(r113f)* c(sO2) &
             + k(r4b)* c(sH) &
             + k(r6b)* c(sH) &
             + k(r12b)* c(sOH) &
             + k(r16b)* c(sCH3) &
             + k(r20b)* c(sHO2) &
             + k(r41b)* c(sH2) &
             + k(r60b)* c(sH) &
             + k(r70b)* c(sH2O) &
             + k(r81b)* c(sH2O2) &
             + k(r82b)* c(sO) &
             + k(r86b)* c(sCO) &
             + k(r108b)* c(sCH4) &
             + k(r129b)* c(sCH3)   )

    HCO_num = ( 0.0_pr &
             + k(r12f)* c(sCH2O) * c(sO) &
             + k(r16f)* c(sC2H4) * c(sO) &
             + k(r20f)* c(sCH2O) * c(sO2) &
             + k(r41f)* c(sCH2O) * c(sH) &
             + k(r70f)* c(sCH2O) * c(sOH) &
             + k(r81f)* c(sCH2O) * c(sHO2) &
             + k(r108f)* c(sCH2O) * c(sCH3) &
             + k(r10b)* c(sCO) * c(sOH) &
             + k(r11b)* c(sCO2) * c(sH) &
             + k(r37b)* c(sCH2O) &
             + k(r38b)* c(sCO) * c(sH2) &
             + k(r69b)* c(sCO) * c(sH2O) &
             + k(r107b)* c(sCH4) * c(sCO) &
             + k(r111b)* c(sCO) * c(sH) * c(sH2O) &
             + k(r112b)* c(sCO) * c(sH) * m(mM112) &
             + k(r113b)* c(sCO) * c(sHO2)   )

    HCO_ct = HCO_num / HCO_denom

    HCO_CH = - ( 0.0_pr &
          + k(r60f) * c(sOH) &
          + k(r82f) * c(sO2) &
          + k(r86f) * c(sCO2)  ) / HCO_denom

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

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

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

    CH2CHO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r126f)* c(sO) &
             + k(r127f)* c(sO2) &
             + k(r128f)* c(sO2) &
             + k(r129f)* c(sH) &
             + k(r119b)* c(sH)   )

    CH2CHO_num = ( 0.0_pr &
             + k(r119f)* c(sC2H4) * c(sO)   )

    CH2CHO_ct = CH2CHO_num / CH2CHO_denom

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

    C2H2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r15f)* c(sO) &
             + k(r74f)* c(sOH) &
             + k(r116b)* c(sH2)   )

    C2H2_num = ( 0.0_pr &
             + k(r116f)* c(sC2H4) &
             + k(r74b)* c(sCH3) * c(sCO)   )

    C2H2_ct = C2H2_num / C2H2_denom

    C2H2_CH2 = - ( 0.0_pr &
          + k(r15b) * c(sCO)  ) / C2H2_denom

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r3f)* c(sO) &
             + k(r60f)* c(sOH) &
             + k(r82f)* c(sO2) &
             + k(r83f)* c(sH2) &
             + k(r84f)* c(sH2O) &
             + k(r85f)* c(sCH4) &
             + k(r86f)* c(sCO2) &
             + k(r122f)* c(sH2) &
             + k(r34b)* c(sH2) &
             + k(r62b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r3b)* c(sCO) * c(sH) &
             + k(r84b)* c(sCH2O) * c(sH) &
             + k(r85b)* c(sC2H4) * c(sH) &
             + k(r122b)* c(sCH3)   )

    CH_ct = CH_num / CH_denom

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

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

    CH_HCO = - ( 0.0_pr &
          + k(r60b) * c(sH) &
          + k(r82b) * c(sO) &
          + k(r86b) * c(sCO)  ) / CH_denom

    CH2OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r13f)* c(sO) &
             + k(r42f)* c(sH) &
             + k(r43f)* c(sH) &
             + k(r44f)* c(sH) &
             + k(r71f)* c(sOH) &
             + k(r114f)* c(sO2) &
             + k(r14b)* c(sOH) &
             + k(r39b)&
             + k(r48b)* c(sH2) &
             + k(r72b)* c(sH2O)   )

    CH2OH_num = ( 0.0_pr &
             + k(r39f)* c(sCH2O) * c(sH) &
             + k(r13b)* c(sCH2O) * c(sOH) &
             + k(r42b)* c(sCH2O) * c(sH2) &
             + k(r43b)* c(sCH3) * c(sOH) &
             + k(r71b)* c(sCH2O) * c(sH2O) &
             + k(r114b)* c(sCH2O) * c(sHO2)   )

    CH2OH_ct = CH2OH_num / CH2OH_denom

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

    CH2OH_CH3OH = - ( 0.0_pr &
          + k(r14f) * c(sO) &
          + k(r48f) * c(sH) &
          + k(r72f) * c(sOH)  ) / CH2OH_denom

    CH3O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r45f)* c(sH) &
             + k(r46f)* c(sH) &
             + k(r47f)* c(sH) &
             + k(r115f)* c(sO2) &
             + k(r40b)&
             + k(r49b)* c(sH2) &
             + k(r73b)* c(sH2O) &
             + k(r79b)* c(sOH) &
             + k(r102b)* c(sO) &
             + k(r109b)* c(sCH4)   )

    CH3O_num = ( 0.0_pr &
             + k(r40f)* c(sCH2O) * c(sH) &
             + k(r79f)* c(sCH3) * c(sHO2) &
             + k(r102f)* c(sCH3) * c(sO2) &
             + k(r45b)* c(sCH2O) * c(sH2) &
             + k(r46b)* c(sCH3) * c(sOH) &
             + k(r115b)* c(sCH2O) * c(sHO2)   )

    CH3O_ct = CH3O_num / CH3O_denom

    CH3O_CH2GSG = - ( 0.0_pr &
          + k(r47b) * c(sH2O)  ) / CH3O_denom

    CH3O_CH3OH = - ( 0.0_pr &
          + k(r49f) * c(sH) &
          + k(r73f) * c(sOH) &
          + k(r109f) * c(sCH3)  ) / CH3O_denom

    CH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r4f)* c(sO) &
             + k(r33f)* c(sH) &
             + k(r61f)* c(sOH) &
             + k(r62f)* c(sOH) &
             + k(r77f)* c(sHO2) &
             + k(r87f)* c(sO2) &
             + k(r88f)* c(sH2) &
             + k(r89f)* c(sCH3) &
             + k(r90f)* c(sCH4) &
             + k(r123f)* c(sO2) &
             + k(r124f)* c(sO2) &
             + k(r15b)* c(sCO) &
             + k(r65b)* c(sH2O) &
             + k(r83b)* c(sH) &
             + k(r91b)* c(sN2) &
             + k(r96b)* c(sH2O) &
             + k(r99b)* c(sCO) &
             + k(r100b)* c(sCO2)   )

    CH2_num = ( 0.0_pr &
             + k(r65f)* c(sCH3) * c(sOH) &
             + k(r33b)* c(sCH3) &
             + k(r61b)* c(sCH2O) * c(sH) &
             + k(r77b)* c(sCH2O) * c(sOH) &
             + k(r88b)* c(sCH3) * c(sH) &
             + k(r89b)* c(sC2H4) * c(sH) &
             + k(r90b) *c(sCH3)** 2.0_pr&
             + k(r124b)* c(sCH2O) * c(sO)   )

    CH2_ct = CH2_num / CH2_denom

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

    CH2_CH2GSG = - ( 0.0_pr &
          + k(r91f) * c(sN2) &
          + k(r96f) * c(sH2O) &
          + k(r99f) * c(sCO) &
          + k(r100f) * c(sCO2)  ) / CH2_denom

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

    CH2_C2H2 = - ( 0.0_pr &
          + k(r15f) * c(sO)  ) / CH2_denom

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

    CH3OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r14f)* c(sO) &
             + k(r48f)* c(sH) &
             + k(r49f)* c(sH) &
             + k(r72f)* c(sOH) &
             + k(r73f)* c(sOH) &
             + k(r109f)* c(sCH3) &
             + k(r64b)&
             + k(r95b)  )

    CH3OH_num = ( 0.0_pr &
             + k(r64f)* c(sCH3) * c(sOH)   )

    CH3OH_ct = CH3OH_num / CH3OH_denom

    CH3OH_CH2GSG = - ( 0.0_pr &
          + k(r95f) * c(sH2O)  ) / CH3OH_denom

    CH3OH_CH2OH = - ( 0.0_pr &
          + k(r14b) * c(sOH) &
          + k(r48b) * c(sH2) &
          + k(r72b) * c(sH2O)  ) / CH3OH_denom

    CH3OH_CH3O = - ( 0.0_pr &
          + k(r49b) * c(sH2) &
          + k(r73b) * c(sH2O) &
          + k(r109b) * c(sCH4)  ) / CH3OH_denom

    CH2GSG_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r5f)* c(sO) &
             + k(r6f)* c(sO) &
             + k(r34f)* c(sH) &
             + k(r63f)* c(sOH) &
             + k(r91f)* c(sN2) &
             + k(r92f)* c(sO2) &
             + k(r93f)* c(sO2) &
             + k(r94f)* c(sH2) &
             + k(r95f)* c(sH2O) &
             + k(r96f)* c(sH2O) &
             + k(r97f)* c(sCH3) &
             + k(r98f)* c(sCH4) &
             + k(r99f)* c(sCO) &
             + k(r100f)* c(sCO2) &
             + k(r101f)* c(sCO2) &
             + k(r125f)* c(sH2O) &
             + k(r44b)* c(sH2O) &
             + k(r47b)* c(sH2O) &
             + k(r66b)* c(sH2O)   )

    CH2GSG_num = ( 0.0_pr &
             + k(r66f)* c(sCH3) * c(sOH) &
             + k(r5b)* c(sCO) * c(sH2) &
             + k(r63b)* c(sCH2O) * c(sH) &
             + k(r92b)* c(sCO) * c(sH) * c(sOH) &
             + k(r93b)* c(sCO) * c(sH2O) &
             + k(r94b)* c(sCH3) * c(sH) &
             + k(r97b)* c(sC2H4) * c(sH) &
             + k(r98b) *c(sCH3)** 2.0_pr&
             + k(r101b)* c(sCH2O) * c(sCO)   )

    CH2GSG_ct = CH2GSG_num / CH2GSG_denom

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

    CH2GSG_CH2 = - ( 0.0_pr &
          + k(r91b) * c(sN2) &
          + k(r96b) * c(sH2O) &
          + k(r99b) * c(sCO) &
          + k(r100b) * c(sCO2)  ) / CH2GSG_denom

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

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

    CH2GSG_CH3O = - ( 0.0_pr &
          + k(r47f) * c(sH)  ) / CH2GSG_denom

    CH2GSG_CH3OH = - ( 0.0_pr &
          + k(r95b)  ) / CH2GSG_denom

    a_j_a =  (1.0_pr) - (CH2GSG_CH3OH) * (CH3OH_CH2GSG)
    a_k_a =  - (CH2GSG_CH2) * (CH3OH_CH2GSG)
    a_l_a =  (CH3OH_CH3O) - (CH2GSG_CH3O) * (CH3OH_CH2GSG)
    a_m_a =  (CH3OH_CH2OH) - (CH2GSG_CH2OH) * (CH3OH_CH2GSG)
    a_n_a =  - (CH2GSG_CH) * (CH3OH_CH2GSG)
    a_q_a =  - (CH2GSG_HCO) * (CH3OH_CH2GSG)
    a_s_a =  - (CH2GSG_CH3OH) * (CH2_CH2GSG)
    a_t_a =  (1.0_pr) - (CH2GSG_CH2) * (CH2_CH2GSG)
    a_u_a =  - (CH2GSG_CH3O) * (CH2_CH2GSG)
    a_v_a =  - (CH2GSG_CH2OH) * (CH2_CH2GSG)
    a_w_a =  (CH2_CH) - (CH2GSG_CH) * (CH2_CH2GSG)
    a_z_a =  (CH2_HCO) - (CH2GSG_HCO) * (CH2_CH2GSG)
    a_ab_a =  (CH3O_CH3OH) - (CH2GSG_CH3OH) * (CH3O_CH2GSG)
    a_ac_a =  - (CH2GSG_CH2) * (CH3O_CH2GSG)
    a_ad_a =  (1.0_pr) - (CH2GSG_CH3O) * (CH3O_CH2GSG)
    a_ae_a =  - (CH2GSG_CH2OH) * (CH3O_CH2GSG)
    a_af_a =  - (CH2GSG_CH) * (CH3O_CH2GSG)
    a_ai_a =  - (CH2GSG_HCO) * (CH3O_CH2GSG)
    a_ak_a =  (CH2OH_CH3OH) - (CH2GSG_CH3OH) * (CH2OH_CH2GSG)
    a_al_a =  - (CH2GSG_CH2) * (CH2OH_CH2GSG)
    a_am_a =  - (CH2GSG_CH3O) * (CH2OH_CH2GSG)
    a_an_a =  (1.0_pr) - (CH2GSG_CH2OH) * (CH2OH_CH2GSG)
    a_ao_a =  - (CH2GSG_CH) * (CH2OH_CH2GSG)
    a_ar_a =  - (CH2GSG_HCO) * (CH2OH_CH2GSG)
    a_at_a =  - (CH2GSG_CH3OH) * (CH_CH2GSG)
    a_au_a =  (CH_CH2) - (CH2GSG_CH2) * (CH_CH2GSG)
    a_av_a =  - (CH2GSG_CH3O) * (CH_CH2GSG)
    a_aw_a =  - (CH2GSG_CH2OH) * (CH_CH2GSG)
    a_ax_a =  (1.0_pr) - (CH2GSG_CH) * (CH_CH2GSG)
    a_ba_a =  (CH_HCO) - (CH2GSG_HCO) * (CH_CH2GSG)
    a_bd_a = C2H2_CH2
    a_bs_a = CH2CHO_HCO
    a_bu_a =  - (CH2GSG_CH3OH) * (HCO_CH2GSG)
    a_bv_a =  (HCO_CH2) - (CH2GSG_CH2) * (HCO_CH2GSG)
    a_bw_a =  - (CH2GSG_CH3O) * (HCO_CH2GSG)
    a_bx_a =  - (CH2GSG_CH2OH) * (HCO_CH2GSG)
    a_by_a =  (HCO_CH) - (CH2GSG_CH) * (HCO_CH2GSG)
    a_cb_a =  (1.0_pr) - (CH2GSG_HCO) * (HCO_CH2GSG)
    a_u_b =  (a_u_a) - (a_l_a) * (a_s_a) / (a_j_a)
    a_v_b =  (a_v_a) - (a_m_a) * (a_s_a) / (a_j_a)
    a_ad_b =  (a_ad_a) - (a_l_a) * (a_ab_a) / (a_j_a)
    a_ae_b =  (a_ae_a) - (a_m_a) * (a_ab_a) / (a_j_a)
    a_am_b =  (a_am_a) - (a_l_a) * (a_ak_a) / (a_j_a)
    a_an_b =  (a_an_a) - (a_m_a) * (a_ak_a) / (a_j_a)
    a_av_b =  (a_av_a) - (a_l_a) * (a_at_a) / (a_j_a)
    a_aw_b =  (a_aw_a) - (a_m_a) * (a_at_a) / (a_j_a)
    a_bw_b =  (a_bw_a) - (a_l_a) * (a_bu_a) / (a_j_a)
    a_bx_b =  (a_bx_a) - (a_m_a) * (a_bu_a) / (a_j_a)
    a_af_c =  (a_af_a) - (a_w_a) * (a_ac_a) / (a_t_a)
    a_ag_c =  - (CH2_C2H2) * (a_ac_a) / (a_t_a)
    a_ah_c =  - (CH2_CH2CHO) * (a_ac_a) / (a_t_a)
    a_ai_c =  (a_ai_a) - (a_z_a) * (a_ac_a) / (a_t_a)
    a_ao_c =  (a_ao_a) - (a_w_a) * (a_al_a) / (a_t_a)
    a_ap_c =  - (CH2_C2H2) * (a_al_a) / (a_t_a)
    a_aq_c =  - (CH2_CH2CHO) * (a_al_a) / (a_t_a)
    a_ar_c =  (a_ar_a) - (a_z_a) * (a_al_a) / (a_t_a)
    a_ax_c =  (a_ax_a) - (a_w_a) * (a_au_a) / (a_t_a)
    a_ay_c =  - (CH2_C2H2) * (a_au_a) / (a_t_a)
    a_bz_c =  - (CH2_CH2CHO) * (a_au_a) / (a_t_a)
    a_ba_c =  (a_ba_a) - (a_z_a) * (a_au_a) / (a_t_a)
    a_bg_c =  - (a_w_a) * (a_bd_a) / (a_t_a)
    a_bh_c =  (1.0_pr) - (CH2_C2H2) * (a_bd_a) / (a_t_a)
    a_bi_c =  - (CH2_CH2CHO) * (a_bd_a) / (a_t_a)
    a_bj_c =  - (a_z_a) * (a_bd_a) / (a_t_a)
    a_br_c = 1.0_pr
    a_bs_c = a_bs_a
    a_by_c =  (a_by_a) - (a_w_a) * (a_bv_a) / (a_t_a)
    a_cz_c =  - (CH2_C2H2) * (a_bv_a) / (a_t_a)
    a_ca_c =  (HCO_CH2CHO) - (CH2_CH2CHO) * (a_bv_a) / (a_t_a)
    a_cb_c =  (a_cb_a) - (a_z_a) * (a_bv_a) / (a_t_a)
    a_bj_f =  (a_bj_c) - (a_ba_c) * (a_bg_c) / (a_ax_c)
    a_bs_f = a_bs_c
    a_cb_f =  (a_cb_c) - (a_ba_c) * (a_by_c) / (a_ax_c)
    a_cb_h =  (a_cb_f) - (a_bs_f) * (a_ca_c) / (a_br_c)
    A_A_A =  (CH3OH_ct) - (CH2GSG_ct) * (CH3OH_CH2GSG)
    A_B_A =  (CH2_ct) - (CH2GSG_ct) * (CH2_CH2GSG)
    A_C_A =  (CH3O_ct) - (CH2GSG_ct) * (CH3O_CH2GSG)
    A_D_A =  (CH2OH_ct) - (CH2GSG_ct) * (CH2OH_CH2GSG)
    A_E_A =  (CH_ct) - (CH2GSG_ct) * (CH_CH2GSG)
    A_F_A = C2H2_ct
    A_G_A = CH2CHO_ct
    A_H_A =  (HCO_ct) - (CH2GSG_ct) * (HCO_CH2GSG)
    A_B_B =  (A_B_A) - (A_A_A) * (a_s_a) / (a_j_a)
    A_C_B =  (A_C_A) - (A_A_A) * (a_ab_a) / (a_j_a)
    A_D_B =  (A_D_A) - (A_A_A) * (a_ak_a) / (a_j_a)
    A_E_B =  (A_E_A) - (A_A_A) * (a_at_a) / (a_j_a)
    A_F_B = A_F_A
    A_G_B = A_G_A
    A_H_B =  (A_H_A) - (A_A_A) * (a_bu_a) / (a_j_a)
    A_C_C =  (A_C_B) - (A_B_B) * (a_ac_a) / (a_t_a)
    A_D_C =  (A_D_B) - (A_B_B) * (a_al_a) / (a_t_a)
    A_E_C =  (A_E_B) - (A_B_B) * (a_au_a) / (a_t_a)
    A_F_C =  (A_F_B) - (A_B_B) * (a_bd_a) / (a_t_a)
    A_G_C = A_G_B
    A_H_C =  (A_H_B) - (A_B_B) * (a_bv_a) / (a_t_a)
    A_D_D =  (A_D_C) - (A_C_C) * (a_am_b) / (a_ad_b)
    A_E_D =  (A_E_C) - (A_C_C) * (a_av_b) / (a_ad_b)
    A_F_D = A_F_C
    A_G_D = A_G_C
    A_H_D =  (A_H_C) - (A_C_C) * (a_bw_b) / (a_ad_b)
    A_E_E =  (A_E_D) - (A_D_D) * (a_aw_b) / (a_an_b)
    A_F_E = A_F_D
    A_G_E = A_G_D
    A_H_E =  (A_H_D) - (A_D_D) * (a_bx_b) / (a_an_b)
    A_F_F =  (A_F_E) - (A_E_E) * (a_bg_c) / (a_ax_c)
    A_G_F = A_G_E
    A_H_F =  (A_H_E) - (A_E_E) * (a_by_c) / (a_ax_c)
    A_G_G = A_G_F
    A_H_G =  (A_H_F) - (A_F_F) * (a_cz_c) / (a_bh_c)
    A_H_H =  (A_H_G) - (A_G_G) * (a_ca_c) / (a_br_c)

    cqss(sqssHCO) = ( A_H_H ) / ( a_cb_h )

    cqss(sqssCH2CHO) = (A_G_G - (a_bs_f) * cqss(sqssHCO)) / (a_br_c)

    cqss(sqssC2H2) = (A_F_F - (a_bj_f) * cqss(sqssHCO) - (a_bi_c) &
                 * cqss(sqssCH2CHO)) / (a_bh_c)

    cqss(sqssCH) = (A_E_E - (a_ba_c) * cqss(sqssHCO) - (a_bz_c) &
                 * cqss(sqssCH2CHO) - (a_ay_c) * cqss(sqssC2H2)) &
                 / (a_ax_c)

    cqss(sqssCH2OH) = (A_D_D - (a_ar_c) * cqss(sqssHCO) - (a_aq_c) &
                 * cqss(sqssCH2CHO) - (a_ap_c) * cqss(sqssC2H2) &
                 - (a_ao_c) * cqss(sqssCH)) / (a_an_b)

    cqss(sqssCH3O) = (A_C_C - (a_ai_c) * cqss(sqssHCO) - (a_ah_c) &
                 * cqss(sqssCH2CHO) - (a_ag_c) * cqss(sqssC2H2) &
                 - (a_af_c) * cqss(sqssCH) - (a_ae_b) * cqss(sqssCH2OH)) &
                 / (a_ad_b)

    cqss(sqssCH2) = (A_B_B - (a_z_a) * cqss(sqssHCO) - (CH2_CH2CHO) &
                 * cqss(sqssCH2CHO) - (CH2_C2H2) * cqss(sqssC2H2) &
                 - (a_w_a) * cqss(sqssCH) - (a_v_b) * cqss(sqssCH2OH) &
                 - (a_u_b) * cqss(sqssCH3O)) / (a_t_a)

    cqss(sqssCH3OH) = (A_A_A - (a_q_a) * cqss(sqssHCO) - (a_n_a) &
                 * cqss(sqssCH) - (a_m_a) * cqss(sqssCH2OH) &
                 - (a_l_a) * cqss(sqssCH3O) - (a_k_a) * cqss(sqssCH2)) &
                 / (a_j_a)

    cqss(sqssCH2GSG) = CH2GSG_ct - (CH2GSG_HCO) * cqss(sqssHCO) - (CH2GSG_CH) &
                 * cqss(sqssCH) - (CH2GSG_CH2OH) * cqss(sqssCH2OH) &
                 - (CH2GSG_CH3O) * cqss(sqssCH3O) - (CH2GSG_CH2) &
                 * cqss(sqssCH2) - (CH2GSG_CH3OH) * cqss(sqssCH3OH)

    C2H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r17f)* c(sO) &
             + k(r51f)* c(sH) &
             + k(r52f)* c(sH) &
             + k(r117f)* c(sO2) &
             + k(r18b)* c(sOH) &
             + k(r50b)&
             + k(r53b)* c(sH2) &
             + k(r75b)* c(sH2O) &
             + k(r106b)* c(sH) &
             + k(r110b)* c(sCH4)   )

    C2H5_num = ( 0.0_pr &
             + k(r18f)* c(sC2H6) * c(sO) &
             + k(r50f)* c(sC2H4) * c(sH) &
             + k(r53f)* c(sC2H6) * c(sH) &
             + k(r75f)* c(sC2H6) * c(sOH) &
             + k(r106f) *c(sCH3)** 2.0_pr&
             + k(r110f)* c(sC2H6) * c(sCH3) &
             + k(r17b)* c(sCH2O) * c(sCH3) &
             + k(r51b)* c(sC2H6) &
             + k(r52b)* c(sC2H4) * c(sH2) &
             + k(r117b)* c(sC2H4) * c(sHO2)   )

    C2H5_ct = C2H5_num / C2H5_denom

    cqss(sqssC2H5) = C2H5_ct

    cqss = max(cqss, tiny(1.0_pr))
    cqss = min(cqss, 1e03_pr)

    return
  end subroutine get_qss


  ! ----------------------------------------------- !
  ! Mass fractions to concentrations                !
  ! ----------------------------------------------- !
  subroutine y2c(y, W_sp, P, T, c)
    implicit none

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

    integer :: k

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

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

    ! Conversion
    c = y * rho / W_sp

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

        return
    end subroutine y2c

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

  use mod_customkinetics
  implicit none

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

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

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

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

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

  ! Evaluate production rates
  call get_production_rates(wdot,w)


  return
end subroutine customkinetics