!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE ./database/S22R143QSS14_0/S22R143QSS14_0.f90
!>    @file ./database/S22R143QSS14_0/S22R143QSS14_0.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 = 4

  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 22
  integer, parameter :: nqss = 14
  integer, parameter :: nreac = 143
  integer, parameter :: nreac_reverse = 143

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

  ! Number of thirdbodies
  integer, parameter :: nTB = 6
  integer, parameter :: nFO = 10
  integer, parameter :: nTB_reverse = 6
  integer, parameter :: nFO_reverse = 10

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

  ! 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 :: sCH3OH = 14
  integer, parameter :: sC2H2 = 15
  integer, parameter :: sC2H4 = 16
  integer, parameter :: sC2H6 = 17
  integer, parameter :: sNO = 18
  integer, parameter :: sNO2 = 19
  integer, parameter :: sN2O = 20
  integer, parameter :: sHCN = 21
  integer, parameter :: sN2 = 22

  integer, parameter :: sqssCH = 1
  integer, parameter :: sqssCH2 = 2
  integer, parameter :: sqssCH2GSG = 3
  integer, parameter :: sqssHCO = 4
  integer, parameter :: sqssCH3O = 5
  integer, parameter :: sqssC2H3 = 6
  integer, parameter :: sqssC2H5 = 7
  integer, parameter :: sqssN = 8
  integer, parameter :: sqssNH = 9
  integer, parameter :: sqssNH2 = 10
  integer, parameter :: sqssNNH = 11
  integer, parameter :: sqssHNO = 12
  integer, parameter :: sqssHOCN = 13
  integer, parameter :: sqssHNCO = 14

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r39f = 39
  integer, parameter :: r40f = 40
  integer, parameter :: r41f = 41
  integer, parameter :: r42f = 42
  integer, parameter :: r43f = 43
  integer, parameter :: r44f = 44
  integer, parameter :: r45f = 45
  integer, parameter :: r46f = 46
  integer, parameter :: r47f = 47
  integer, parameter :: r48f = 48
  integer, parameter :: r49f = 49
  integer, parameter :: r50f = 50
  integer, parameter :: r51f = 51
  integer, parameter :: r52f = 52
  integer, parameter :: r53f = 53
  integer, parameter :: r54f = 54
  integer, parameter :: r55f = 55
  integer, parameter :: r56f = 56
  integer, parameter :: r57f = 57
  integer, parameter :: r58f = 58
  integer, parameter :: r59f = 59
  integer, parameter :: r60f = 60
  integer, parameter :: r61f = 61
  integer, parameter :: r62f = 62
  integer, parameter :: r63f = 63
  integer, parameter :: r64f = 64
  integer, parameter :: r65f = 65
  integer, parameter :: r66f = 66
  integer, parameter :: r67f = 67
  integer, parameter :: r68f = 68
  integer, parameter :: r69f = 69
  integer, parameter :: r70f = 70
  integer, parameter :: r71f = 71
  integer, parameter :: r72f = 72
  integer, parameter :: r73f = 73
  integer, parameter :: r74f = 74
  integer, parameter :: r75f = 75
  integer, parameter :: r76f = 76
  integer, parameter :: r77f = 77
  integer, parameter :: r78f = 78
  integer, parameter :: r79f = 79
  integer, parameter :: r80f = 80
  integer, parameter :: r81f = 81
  integer, parameter :: r82f = 82
  integer, parameter :: r83f = 83
  integer, parameter :: r84f = 84
  integer, parameter :: r85f = 85
  integer, parameter :: r86f = 86
  integer, parameter :: r87f = 87
  integer, parameter :: r88f = 88
  integer, parameter :: r89f = 89
  integer, parameter :: r90f = 90
  integer, parameter :: r91f = 91
  integer, parameter :: r92f = 92
  integer, parameter :: r93f = 93
  integer, parameter :: r94f = 94
  integer, parameter :: r95f = 95
  integer, parameter :: r96f = 96
  integer, parameter :: r97f = 97
  integer, parameter :: r98f = 98
  integer, parameter :: r99f = 99
  integer, parameter :: r100f = 100
  integer, parameter :: r101f = 101
  integer, parameter :: r102f = 102
  integer, parameter :: r103f = 103
  integer, parameter :: r104f = 104
  integer, parameter :: r105f = 105
  integer, parameter :: r106f = 106
  integer, parameter :: r107f = 107
  integer, parameter :: r108f = 108
  integer, parameter :: r109f = 109
  integer, parameter :: r110f = 110
  integer, parameter :: r111f = 111
  integer, parameter :: r112f = 112
  integer, parameter :: r113f = 113
  integer, parameter :: r114f = 114
  integer, parameter :: r115f = 115
  integer, parameter :: r116f = 116
  integer, parameter :: r117f = 117
  integer, parameter :: r118f = 118
  integer, parameter :: r119f = 119
  integer, parameter :: r120f = 120
  integer, parameter :: r121f = 121
  integer, parameter :: r122f = 122
  integer, parameter :: r123f = 123
  integer, parameter :: r124f = 124
  integer, parameter :: r125f = 125
  integer, parameter :: r126f = 126
  integer, parameter :: r127f = 127
  integer, parameter :: r128f = 128
  integer, parameter :: r129f = 129
  integer, parameter :: r130f = 130
  integer, parameter :: r131f = 131
  integer, parameter :: r132f = 132
  integer, parameter :: r133f = 133
  integer, parameter :: r134f = 134
  integer, parameter :: r135f = 135
  integer, parameter :: r136f = 136
  integer, parameter :: r137f = 137
  integer, parameter :: r138f = 138
  integer, parameter :: r139f = 139
  integer, parameter :: r140f = 140
  integer, parameter :: r141f = 141
  integer, parameter :: r142f = 142
  integer, parameter :: r143f = 143
  integer, parameter :: r1b = 144
  integer, parameter :: r2b = 145
  integer, parameter :: r3b = 146
  integer, parameter :: r4b = 147
  integer, parameter :: r5b = 148
  integer, parameter :: r6b = 149
  integer, parameter :: r7b = 150
  integer, parameter :: r8b = 151
  integer, parameter :: r9b = 152
  integer, parameter :: r10b = 153
  integer, parameter :: r11b = 154
  integer, parameter :: r12b = 155
  integer, parameter :: r13b = 156
  integer, parameter :: r14b = 157
  integer, parameter :: r15b = 158
  integer, parameter :: r16b = 159
  integer, parameter :: r17b = 160
  integer, parameter :: r18b = 161
  integer, parameter :: r19b = 162
  integer, parameter :: r20b = 163
  integer, parameter :: r21b = 164
  integer, parameter :: r22b = 165
  integer, parameter :: r23b = 166
  integer, parameter :: r24b = 167
  integer, parameter :: r25b = 168
  integer, parameter :: r26b = 169
  integer, parameter :: r27b = 170
  integer, parameter :: r28b = 171
  integer, parameter :: r29b = 172
  integer, parameter :: r30b = 173
  integer, parameter :: r31b = 174
  integer, parameter :: r32b = 175
  integer, parameter :: r33b = 176
  integer, parameter :: r34b = 177
  integer, parameter :: r35b = 178
  integer, parameter :: r36b = 179
  integer, parameter :: r37b = 180
  integer, parameter :: r38b = 181
  integer, parameter :: r39b = 182
  integer, parameter :: r40b = 183
  integer, parameter :: r41b = 184
  integer, parameter :: r42b = 185
  integer, parameter :: r43b = 186
  integer, parameter :: r44b = 187
  integer, parameter :: r45b = 188
  integer, parameter :: r46b = 189
  integer, parameter :: r47b = 190
  integer, parameter :: r48b = 191
  integer, parameter :: r49b = 192
  integer, parameter :: r50b = 193
  integer, parameter :: r51b = 194
  integer, parameter :: r52b = 195
  integer, parameter :: r53b = 196
  integer, parameter :: r54b = 197
  integer, parameter :: r55b = 198
  integer, parameter :: r56b = 199
  integer, parameter :: r57b = 200
  integer, parameter :: r58b = 201
  integer, parameter :: r59b = 202
  integer, parameter :: r60b = 203
  integer, parameter :: r61b = 204
  integer, parameter :: r62b = 205
  integer, parameter :: r63b = 206
  integer, parameter :: r64b = 207
  integer, parameter :: r65b = 208
  integer, parameter :: r66b = 209
  integer, parameter :: r67b = 210
  integer, parameter :: r68b = 211
  integer, parameter :: r69b = 212
  integer, parameter :: r70b = 213
  integer, parameter :: r71b = 214
  integer, parameter :: r72b = 215
  integer, parameter :: r73b = 216
  integer, parameter :: r74b = 217
  integer, parameter :: r75b = 218
  integer, parameter :: r76b = 219
  integer, parameter :: r77b = 220
  integer, parameter :: r78b = 221
  integer, parameter :: r79b = 222
  integer, parameter :: r80b = 223
  integer, parameter :: r81b = 224
  integer, parameter :: r82b = 225
  integer, parameter :: r83b = 226
  integer, parameter :: r84b = 227
  integer, parameter :: r85b = 228
  integer, parameter :: r86b = 229
  integer, parameter :: r87b = 230
  integer, parameter :: r88b = 231
  integer, parameter :: r89b = 232
  integer, parameter :: r90b = 233
  integer, parameter :: r91b = 234
  integer, parameter :: r92b = 235
  integer, parameter :: r93b = 236
  integer, parameter :: r94b = 237
  integer, parameter :: r95b = 238
  integer, parameter :: r96b = 239
  integer, parameter :: r97b = 240
  integer, parameter :: r98b = 241
  integer, parameter :: r99b = 242
  integer, parameter :: r100b = 243
  integer, parameter :: r101b = 244
  integer, parameter :: r102b = 245
  integer, parameter :: r103b = 246
  integer, parameter :: r104b = 247
  integer, parameter :: r105b = 248
  integer, parameter :: r106b = 249
  integer, parameter :: r107b = 250
  integer, parameter :: r108b = 251
  integer, parameter :: r109b = 252
  integer, parameter :: r110b = 253
  integer, parameter :: r111b = 254
  integer, parameter :: r112b = 255
  integer, parameter :: r113b = 256
  integer, parameter :: r114b = 257
  integer, parameter :: r115b = 258
  integer, parameter :: r116b = 259
  integer, parameter :: r117b = 260
  integer, parameter :: r118b = 261
  integer, parameter :: r119b = 262
  integer, parameter :: r120b = 263
  integer, parameter :: r121b = 264
  integer, parameter :: r122b = 265
  integer, parameter :: r123b = 266
  integer, parameter :: r124b = 267
  integer, parameter :: r125b = 268
  integer, parameter :: r126b = 269
  integer, parameter :: r127b = 270
  integer, parameter :: r128b = 271
  integer, parameter :: r129b = 272
  integer, parameter :: r130b = 273
  integer, parameter :: r131b = 274
  integer, parameter :: r132b = 275
  integer, parameter :: r133b = 276
  integer, parameter :: r134b = 277
  integer, parameter :: r135b = 278
  integer, parameter :: r136b = 279
  integer, parameter :: r137b = 280
  integer, parameter :: r138b = 281
  integer, parameter :: r139b = 282
  integer, parameter :: r140b = 283
  integer, parameter :: r141b = 284
  integer, parameter :: r142b = 285
  integer, parameter :: r143b = 286

  ! Index of third body species
  integer, parameter :: mM14 = 1
  integer, parameter :: mM20 = 2
  integer, parameter :: mM93 = 3
  integer, parameter :: mM107 = 4
  integer, parameter :: mM122 = 5
  integer, parameter :: mM125 = 6

  integer, parameter :: mM24 = 7
  integer, parameter :: mM26 = 8
  integer, parameter :: mM29 = 9
  integer, parameter :: mM35 = 10
  integer, parameter :: mM37 = 11
  integer, parameter :: mM39 = 12
  integer, parameter :: mM42 = 13
  integer, parameter :: mM50 = 14
  integer, parameter :: mM87 = 15
  integer, parameter :: mM105 = 16

  ! Index of third body reactions
  integer, parameter :: TBr14f = 1
  integer, parameter :: TBr20f = 2
  integer, parameter :: TBr93f = 3
  integer, parameter :: TBr107f = 4
  integer, parameter :: TBr122f = 5
  integer, parameter :: TBr125f = 6
  integer, parameter :: TBr14b = 7
  integer, parameter :: TBr20b = 8
  integer, parameter :: TBr93b = 9
  integer, parameter :: TBr107b = 10
  integer, parameter :: TBr122b = 11
  integer, parameter :: TBr125b = 12

  ! Index of fall off reactions
  integer, parameter :: FOr24f = 1
  integer, parameter :: FOr26f = 2
  integer, parameter :: FOr29f = 3
  integer, parameter :: FOr35f = 4
  integer, parameter :: FOr37f = 5
  integer, parameter :: FOr39f = 6
  integer, parameter :: FOr42f = 7
  integer, parameter :: FOr50f = 8
  integer, parameter :: FOr87f = 9
  integer, parameter :: FOr105f = 10
  integer, parameter :: FOr24b = 11
  integer, parameter :: FOr26b = 12
  integer, parameter :: FOr29b = 13
  integer, parameter :: FOr35b = 14
  integer, parameter :: FOr37b = 15
  integer, parameter :: FOr39b = 16
  integer, parameter :: FOr42b = 17
  integer, parameter :: FOr50b = 18
  integer, parameter :: FOr87b = 19
  integer, parameter :: FOr105b = 20

  ! 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.03204216_pr, & ! CH3OH
       0.02603788_pr, & ! C2H2
       0.02805376_pr, & ! C2H4
       0.03006964_pr, & ! C2H6
       0.03000614_pr, & ! NO
       0.04600554_pr, & ! NO2
       0.04401288_pr, & ! N2O
       0.02702568_pr, & ! HCN
       0.02801348_pr & ! N2
      !0.01301894_pr, & ! CH
      !0.01402688_pr, & ! CH2
      !0.01402688_pr, & ! CH2GSG
      !0.02901834_pr, & ! HCO
      !0.03103422_pr, & ! CH3O
      !0.02704582_pr, & ! C2H3
      !0.0290617_pr, & ! C2H5
      !0.01400674_pr, & ! N
      !0.01501468_pr, & ! NH
      !0.01602262_pr, & ! NH2
      !0.02902142_pr, & ! NNH
      !0.03101408_pr, & ! HNO
      !0.04302508_pr, & ! HOCN
      !0.04302508_pr & ! HNCO
       /)

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(mM14) = (-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(mM20) = (1_pr)*c(sCH4) &
         + (2.65_pr)*c(sH2O) &
         + sum(c)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  end subroutine get_thirdbodies

  ! ----------------------------------------------- !
  ! Evaluate rate coefficients                      !
  ! ----------------------------------------------- !
  subroutine get_rate_coefficients(k,M,Tloc,Ploc)
    implicit none

    real(pr), dimension(nreac + nreac_reverse) :: k
    real(pr), dimension(nFO + nFO_reverse) :: k_0
    real(pr), dimension(nFO + nFO_reverse) :: k_inf
    real(pr), dimension(nFO + nFO_reverse) :: FC
    real(pr), dimension(nTB + nFO) :: M
    real(pr) :: Tloc,Ploc,R_T_inv,T_log

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

    k(r1f) = (0.05_pr)*exp((-26317.36_pr)*R_T_inv + T_log *(2.67_pr))
    k(r2f) = (20000000_pr)
    k(r3f) = (57000000_pr)
    k(r4f) = (80000000_pr)
    k(r5f) = (84300000_pr)
    k(r6f) = (1020_pr)*exp((-35982.4_pr)*R_T_inv + T_log *(1.5_pr))
    k(r7f) = (30000000_pr)
    k(r8f) = (39000000_pr)*exp((-14811.36_pr)*R_T_inv)
    k(r9f) = (10.2_pr)*exp((-7949.6_pr)*R_T_inv + T_log *(2_pr))
    k(r10f) = (19.2_pr)*exp((-920.48_pr)*R_T_inv + T_log *(1.83_pr))
    k(r11f) = (132000000_pr)
    k(r12f) = (89.8_pr)*exp((-23806.96_pr)*R_T_inv + T_log *(1.92_pr))
    k(r13f) = (100000000_pr)*exp((-167360_pr)*R_T_inv)
    k(r14f) = (2800000_pr)*exp(T_log *(-0.86_pr))
    k(r15f) = (300000000_pr)*exp(T_log *(-1.72_pr))
    k(r16f) = (9380000_pr)*exp(T_log *(-0.76_pr))
    k(r17f) = (375000000_pr)*exp(T_log *(-1.72_pr))
    k(r18f) = (83000000_pr)*exp((-60303.992_pr)*R_T_inv)
    k(r19f) = (60000000_pr)*exp(T_log *(-1.25_pr))
    k(r20f) = (22000000000_pr)*exp(T_log *(-2_pr))
    k(r21f) = (3970000_pr)*exp((-2807.464_pr)*R_T_inv)
    k(r22f) = (28000000_pr)*exp((-4468.512_pr)*R_T_inv)
    k(r23f) = (134000000_pr)*exp((-2656.84_pr)*R_T_inv)
    k_0(FOr24f) = (3.2e+15_pr)*exp((-5146.32_pr)*R_T_inv + T_log *(-3.14_pr))
    k_inf(FOr24f) = (25000000000_pr)*exp(T_log *(-0.8_pr))
    FC(FOr24f) = ((1.0_pr - 0.68_pr)*exp(-Tloc/(78_pr))) + (0.68_pr*exp(-Tloc/(1995_pr)))&
              + exp(-(5590_pr)/Tloc)
    k(r24f) = getlindratecoeff(Tloc,k_0(FOr24f),k_inf(FOr24f),FC(FOr24f), M(mM24),Ploc)
    k(r25f) = (30000000_pr)
    k_0(FOr26f) = (2.477e+21_pr)*exp((-10208.96_pr)*R_T_inv + T_log *(-4.76_pr))
    k_inf(FOr26f) = (12700000000_pr)*exp((-1602.472_pr)*R_T_inv + T_log *(-0.63_pr))
    FC(FOr26f) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r26f) = getlindratecoeff(Tloc,k_0(FOr26f),k_inf(FOr26f),FC(FOr26f), M(mM26),Ploc)
    k(r27f) = (660_pr)*exp((-45354.56_pr)*R_T_inv + T_log *(1.62_pr))
    k(r28f) = (73400000_pr)
    k_0(FOr29f) = (2.2e+18_pr)*exp((-23263.04_pr)*R_T_inv + T_log *(-4.8_pr))
    k_inf(FOr29f) = (540000_pr)*exp((-10878.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr29f) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r29f) = getlindratecoeff(Tloc,k_0(FOr29f),k_inf(FOr29f),FC(FOr29f), M(mM29),Ploc)
    k(r30f) = (23000_pr)*exp((-13702.6_pr)*R_T_inv + T_log *(1.05_pr))
    k(r31f) = (20000000_pr)
    k(r32f) = (32000000_pr)
    k(r33f) = (16000000_pr)
    k(r34f) = (4.2_pr)*exp((-20376.08_pr)*R_T_inv + T_log *(2.1_pr))
    k_0(FOr35f) = (3.8e+28_pr)*exp((-30208.48_pr)*R_T_inv + T_log *(-7.27_pr))
    k_inf(FOr35f) = (5600000_pr)*exp((-10041.6_pr)*R_T_inv)
    FC(FOr35f) = ((1.0_pr - 0.7507_pr)*exp(-Tloc/(98.5_pr))) + (0.7507_pr*exp(-Tloc/(1302_pr)))&
              + exp(-(4167_pr)/Tloc)
    k(r35f) = getlindratecoeff(Tloc,k_0(FOr35f),k_inf(FOr35f),FC(FOr35f), M(mM35),Ploc)
    k(r36f) = (30000000_pr)
    k_0(FOr37f) = (1.2e+30_pr)*exp((-29162.48_pr)*R_T_inv + T_log *(-7.62_pr))
    k_inf(FOr37f) = (1080000_pr)*exp((-7614.88_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr37f) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r37f) = getlindratecoeff(Tloc,k_0(FOr37f),k_inf(FOr37f),FC(FOr37f), M(mM37),Ploc)
    k(r38f) = (1.325_pr)*exp((-51212.16_pr)*R_T_inv + T_log *(2.53_pr))
    k_0(FOr39f) = (1.99e+29_pr)*exp((-27970.04_pr)*R_T_inv + T_log *(-7.08_pr))
    k_inf(FOr39f) = (521000000000_pr)*exp((-6610.72_pr)*R_T_inv + T_log *(-0.99_pr))
    FC(FOr39f) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r39f) = getlindratecoeff(Tloc,k_0(FOr39f),k_inf(FOr39f),FC(FOr39f), M(mM39),Ploc)
    k(r40f) = (115_pr)*exp((-31505.52_pr)*R_T_inv + T_log *(1.9_pr))
    k(r41f) = (216_pr)*exp((-14351.12_pr)*R_T_inv + T_log *(1.51_pr))
    k_0(FOr42f) = (2300000_pr)*exp((7112.8_pr)*R_T_inv + T_log *(-0.9_pr))
    k_inf(FOr42f) = (74000000_pr)*exp(T_log *(-0.37_pr))
    FC(FOr42f) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r42f) = getlindratecoeff(Tloc,k_0(FOr42f),k_inf(FOr42f),FC(FOr42f), M(mM42),Ploc)
    k(r43f) = (0.0357_pr)*exp((8828.24_pr)*R_T_inv + T_log *(2.4_pr))
    k(r44f) = (29000000_pr)*exp((2092_pr)*R_T_inv)
    k(r45f) = (1750000_pr)*exp((-1338.88_pr)*R_T_inv)
    k(r46f) = (580000000_pr)*exp((-39999.04_pr)*R_T_inv)
    k(r47f) = (30000000_pr)
    k(r48f) = (20000000_pr)
    k(r49f) = (11.3_pr)*exp((-12552_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr50f) = (2.7e+26_pr)*exp((-12970.4_pr)*R_T_inv + T_log *(-6.3_pr))
    k_inf(FOr50f) = (63000000_pr)
    FC(FOr50f) = ((1.0_pr - 0.2105_pr)*exp(-Tloc/(83.5_pr))) + (0.2105_pr*exp(-Tloc/(5398_pr)))&
              + exp(-(8370_pr)/Tloc)
    k(r50f) = getlindratecoeff(Tloc,k_0(FOr50f),k_inf(FOr50f),FC(FOr50f), M(mM50),Ploc)
    k(r51f) = (56_pr)*exp((-22677.28_pr)*R_T_inv + T_log *(1.6_pr))
    k(r52f) = (25010000_pr)
    k(r53f) = (100_pr)*exp((-13054.08_pr)*R_T_inv + T_log *(1.6_pr))
    k(r54f) = (47.6_pr)*exp((-292.88_pr)*R_T_inv + T_log *(1.228_pr))
    k(r55f) = (50000000_pr)
    k(r56f) = (3430_pr)*exp((1870.248_pr)*R_T_inv + T_log *(1.18_pr))
    k(r57f) = (5000000_pr)
    k(r58f) = (6.3_pr)*exp((-6276_pr)*R_T_inv + T_log *(2_pr))
    k(r59f) = (3.6_pr)*exp((-10460_pr)*R_T_inv + T_log *(2_pr))
    k(r60f) = (3.54_pr)*exp((-3640.08_pr)*R_T_inv + T_log *(2.12_pr))
    k(r61f) = (130000_pr)*exp((6819.92_pr)*R_T_inv)
    k(r62f) = (420000000_pr)*exp((-50208_pr)*R_T_inv)
    k(r63f) = (1000000_pr)
    k(r64f) = (20000000_pr)
    k(r65f) = (150000000_pr)*exp((-98742.4_pr)*R_T_inv)
    k(r66f) = (1000000_pr)*exp((-33472_pr)*R_T_inv)
    k(r67f) = (33000000_pr)
    k(r68f) = (110.7_pr)*exp((-6987.28_pr)*R_T_inv + T_log *(1.79_pr))
    k(r69f) = (17130000_pr)*exp((3158.92_pr)*R_T_inv)
    k(r70f) = (30000000_pr)
    k(r71f) = (60000000_pr)
    k(r72f) = (3400000_pr)*exp((-2886.96_pr)*R_T_inv)
    k(r73f) = (13200000_pr)*exp((-6276_pr)*R_T_inv)
    k(r74f) = (0.5_pr)*exp((-30250.32_pr)*R_T_inv + T_log *(2_pr))
    k(r75f) = (40000000_pr)
    k(r76f) = (15000000_pr)*exp((-2510.4_pr)*R_T_inv)
    k(r77f) = (28000000_pr)
    k(r78f) = (12000000_pr)
    k(r79f) = (70000000_pr)
    k(r80f) = (30000000_pr)
    k(r81f) = (16000000_pr)*exp((2384.88_pr)*R_T_inv)
    k(r82f) = (9000000_pr)
    k(r83f) = (14000000_pr)
    k(r84f) = (26750000_pr)*exp((-120499.2_pr)*R_T_inv)
    k(r85f) = (36000_pr)*exp((-37404.96_pr)*R_T_inv)
    k(r86f) = (0.0245_pr)*exp((-21673.12_pr)*R_T_inv + T_log *(2.47_pr))
    k_0(FOr87f) = (1.77e+38_pr)*exp((-26024.48_pr)*R_T_inv + T_log *(-9.67_pr))
    k_inf(FOr87f) = (21200000000_pr)*exp((-2594.08_pr)*R_T_inv + T_log *(-0.97_pr))
    FC(FOr87f) = ((1.0_pr - 0.5325_pr)*exp(-Tloc/(151_pr))) + (0.5325_pr*exp(-Tloc/(1038_pr)))&
              + exp(-(4970_pr)/Tloc)
    k(r87f) = getlindratecoeff(Tloc,k_0(FOr87f),k_inf(FOr87f),FC(FOr87f), M(mM87),Ploc)
    k(r88f) = (4990000_pr)*exp((-44350.4_pr)*R_T_inv + T_log *(0.1_pr))
    k(r89f) = (26480000_pr)
    k(r90f) = (0.00332_pr)*exp((-24518.24_pr)*R_T_inv + T_log *(2.81_pr))
    k(r91f) = (10_pr)*exp((-41588.96_pr)*R_T_inv + T_log *(1.5_pr))
    k(r92f) = (2244000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r93f) = (187000000000_pr)*exp((-71128_pr)*R_T_inv + T_log *(-1_pr))
    k(r94f) = (7600000_pr)*exp((-1673.6_pr)*R_T_inv)
    k(r95f) = (4.28e-19_pr)*exp((14769.52_pr)*R_T_inv + T_log *(7.6_pr))
    k(r96f) = (3980000_pr)*exp((1004.16_pr)*R_T_inv)
    k(r97f) = (840000_pr)*exp((-16213_pr)*R_T_inv)
    k(r98f) = (35000000_pr)*exp((-1380.72_pr)*R_T_inv)
    k(r99f) = (2650000_pr)*exp((-26777.6_pr)*R_T_inv)
    k(r100f) = (73330000_pr)*exp((-4686.08_pr)*R_T_inv)
    k(r101f) = (1400000_pr)*exp((-45229.04_pr)*R_T_inv)
    k(r102f) = (29000000_pr)*exp((-96859.6_pr)*R_T_inv)
    k(r103f) = (440000000_pr)*exp((-78993.92_pr)*R_T_inv)
    k(r104f) = (2000000_pr)*exp((-88115.04_pr)*R_T_inv)
    k_0(FOr105f) = (620000000_pr)*exp((-234722.4_pr)*R_T_inv)
    k_inf(FOr105f) = (130000000000_pr)*exp((-249450.08_pr)*R_T_inv)
    k(r105f) = getlindratecoeff(Tloc,k_0(FOr105f),k_inf(FOr105f),1.0_pr, M(mM105),Ploc)
    k(r106f) = (2110000_pr)*exp((2008.32_pr)*R_T_inv)
    k(r107f) = (106000000_pr)*exp(T_log *(-1.41_pr))
    k(r108f) = (3900000_pr)*exp((1004.16_pr)*R_T_inv)
    k(r109f) = (132000000_pr)*exp((-1506.24_pr)*R_T_inv)
    k(r110f) = (50000000_pr)
    k(r111f) = (32000000_pr)*exp((-1380.72_pr)*R_T_inv)
    k(r112f) = (20000000_pr)
    k(r113f) = (2000_pr)*exp(T_log *(1.2_pr))
    k(r114f) = (0.461_pr)*exp((-27196_pr)*R_T_inv + T_log *(2_pr))
    k(r115f) = (1.28_pr)*exp((-418.4_pr)*R_T_inv + T_log *(1.5_pr))
    k(r116f) = (20000000_pr)*exp((-57948.4_pr)*R_T_inv)
    k(r117f) = (416000000_pr)*exp(T_log *(-0.45_pr))
    k(r118f) = (7000000_pr)
    k(r119f) = (46000000_pr)
    k(r120f) = (40000000_pr)*exp((-15271.6_pr)*R_T_inv)
    k(r121f) = (90_pr)*exp((1924.64_pr)*R_T_inv + T_log *(1.5_pr))
    k(r122f) = (130000000_pr)*exp((-20836.32_pr)*R_T_inv + T_log *(-0.11_pr))
    k(r123f) = (70000000_pr)
    k(r124f) = (50000000_pr)
    k(r125f) = (89500000_pr)*exp((-3096.16_pr)*R_T_inv + T_log *(-1.32_pr))
    k(r126f) = (25000000_pr)
    k(r127f) = (450000_pr)*exp((-2761.44_pr)*R_T_inv + T_log *(0.72_pr))
    k(r128f) = (13_pr)*exp((3974.8_pr)*R_T_inv + T_log *(1.9_pr))
    k(r129f) = (10000000_pr)*exp((-54392_pr)*R_T_inv)
    k(r130f) = (0.002767_pr)*exp((-20836.32_pr)*R_T_inv + T_log *(2.64_pr))
    k(r131f) = (1.1_pr)*exp((-55940.08_pr)*R_T_inv + T_log *(2.03_pr))
    k(r132f) = (0.0044_pr)*exp((-26777.6_pr)*R_T_inv + T_log *(2.26_pr))
    k(r133f) = (285.7_pr)*exp((-85353.6_pr)*R_T_inv + T_log *(1.1_pr))
    k(r134f) = (50000000_pr)
    k(r135f) = (310000000000_pr)*exp((-5313.68_pr)*R_T_inv + T_log *(-1.38_pr))
    k(r136f) = (290000000_pr)*exp((-3179.84_pr)*R_T_inv + T_log *(-0.69_pr))
    k(r137f) = (310000000000_pr)*exp((-5313.68_pr)*R_T_inv + T_log *(-1.38_pr))
    k(r138f) = (96000000_pr)*exp((-120499.2_pr)*R_T_inv)
    k(r139f) = (98_pr)*exp((-35564_pr)*R_T_inv + T_log *(1.41_pr))
    k(r140f) = (22.5_pr)*exp((-15899.2_pr)*R_T_inv + T_log *(1.7_pr))
    k(r141f) = (1550000_pr)*exp((-28660.4_pr)*R_T_inv)
    k(r142f) = (20_pr)*exp((-8368_pr)*R_T_inv + T_log *(2_pr))
    k(r143f) = (3700000_pr)*exp((376.56_pr)*R_T_inv + T_log *(0.15_pr))
    k(r1b) = (0.0254247930841329_pr)*exp((-20056.9432765112_pr)*R_T_inv + T_log *(2.65432763263115_pr))
    k(r2b) = (4096322.2461569_pr)*exp((-223095.984733352_pr)*R_T_inv + T_log *(0.233057403308661_pr))
    k(r3b) = (3286918582.30116_pr)*exp((-739228.309911701_pr)*R_T_inv + T_log *(-0.0294196747253529_pr))
    k(r4b) = (3423044661.20991_pr)*exp((-382560.288192676_pr)*R_T_inv + T_log *(-0.241942753349928_pr))
    k(r5b) = (3210332194.82639_pr)*exp((-287639.146017133_pr)*R_T_inv + T_log *(-0.119343188392549_pr))
    k(r6b) = (8.70418331850672_pr)*exp((-24207.3965707474_pr)*R_T_inv + T_log *(1.57729390856095_pr))
    k(r7b) = (3452574730326.69_pr)*exp((-468733.033799254_pr)*R_T_inv + T_log *(-0.701289799834911_pr))
    k(r8b) = (300164.092491249_pr)*exp((-74318.0916139794_pr)*R_T_inv + T_log *(0.159120574981163_pr))
    k(r9b) = (6.53872675348412e-05_pr)*exp((-199778.208530644_pr)*R_T_inv + T_log *(3.2042419033545_pr))
    k(r10b) = (0.00217233826275361_pr)*exp((-110702.751088735_pr)*R_T_inv + T_log *(2.42416702058187_pr))
    k(r11b) = (338656.375244391_pr)*exp((-326811.314533663_pr)*R_T_inv + T_log *(0.765172758243467_pr))
    k(r12b) = (0.0705029160430455_pr)*exp((-31136.281193286_pr)*R_T_inv + T_log *(2.20149566413092_pr))
    k(r13b) = (3757768.886961_pr)*exp((-3770.74696628824_pr)*R_T_inv + T_log *(-0.0739368446182744_pr))
    k(r14b) = (1203475020331.91_pr)*exp((-203377.238357849_pr)*R_T_inv + T_log *(-0.72246045035514_pr))
    k(r15b) = (128943750136136_pr)*exp((-203377.238346907_pr)*R_T_inv + T_log *(-1.58246044826867_pr))
    k(r16b) = (4031641090047.29_pr)*exp((-203377.238318793_pr)*R_T_inv + T_log *(-0.622460442903258_pr))
    k(r17b) = (161179682331002_pr)*exp((-203377.238324021_pr)*R_T_inv + T_log *(-1.58246044390852_pr))
    k(r18b) = (364195.722066508_pr)*exp((9745.45480291323_pr)*R_T_inv + T_log *(0.336899610132186_pr))
    k(r19b) = (10387403541152.2_pr)*exp((-432733.639828718_pr)*R_T_inv + T_log *(-0.863730682383367_pr))
    k(r20b) = (1.39203180686837e+16_pr)*exp((-494626.313492865_pr)*R_T_inv + T_log *(-1.58068609519003_pr))
    k(r21b) = (25644.5277057966_pr)*exp((-224007.092366566_pr)*R_T_inv + T_log *(0.618673958709982_pr))
    k(r22b) = (11278068.2093299_pr)*exp((-233824.913440838_pr)*R_T_inv + T_log *(0.248729773730298_pr))
    k(r23b) = (120427.497043328_pr)*exp((-155703.377923252_pr)*R_T_inv + T_log *(0.569957014788095_pr))
    k_0(FOr24b) = (4.11243768549599e+22_pr)*exp((-467033.953869435_pr)*R_T_inv + T_log *(-3.05112322838557_pr))
    k_inf(FOr24b) = (3.21284191304571e+17_pr)*exp((-461887.633863014_pr)*R_T_inv + T_log *(-0.711123227217508_pr))
    FC(FOr24b) = ((1.0_pr - 0.68_pr)*exp(-Tloc/(78_pr))) + (0.68_pr*exp(-Tloc/(1995_pr)))&
              + exp(-(5590_pr)/Tloc)
    k(r24b) = getlindratecoeff(Tloc,k_0(FOr24b),k_inf(FOr24b),FC(FOr24b), M(mM24),Ploc)
    k(r25b) = (4145450.42700321_pr)*exp((-49747.5995677441_pr)*R_T_inv + T_log *(0.0889959698879727_pr))
    k_0(FOr26b) = (2.55529633652742e+28_pr)*exp((-448457.18644934_pr)*R_T_inv + T_log *(-4.46669694218075_pr))
    k_inf(FOr26b) = (1.31014374591391e+17_pr)*exp((-439850.698387665_pr)*R_T_inv + T_log *(-0.336696930470406_pr))
    FC(FOr26b) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r26b) = getlindratecoeff(Tloc,k_0(FOr26b),k_inf(FOr26b),FC(FOr26b), M(mM26),Ploc)
    k(r27b) = (11.0760371472792_pr)*exp((-39839.973348218_pr)*R_T_inv + T_log *(1.71296626577427_pr))
    k(r28b) = (12980909.0067341_pr)*exp((-368724.51077046_pr)*R_T_inv + T_log *(0.409403289556344_pr))
    k_0(FOr29b) = (6.57774969059808e+26_pr)*exp((-111922.241734812_pr)*R_T_inv + T_log *(-5.4652751972734_pr))
    k_inf(FOr29b) = (161453859298397_pr)*exp((-99537.6017487496_pr)*R_T_inv + T_log *(-0.211275199929747_pr))
    FC(FOr29b) = ((1.0_pr - 0.758_pr)*exp(-Tloc/(94_pr))) + (0.758_pr*exp(-Tloc/(1555_pr)))&
              + exp(-(4200_pr)/Tloc)
    k(r29b) = getlindratecoeff(Tloc,k_0(FOr29b),k_inf(FOr29b),FC(FOr29b), M(mM29),Ploc)
    k(r30b) = (348.124453233529_pr)*exp((-79469.748315777_pr)*R_T_inv + T_log *(1.22479294652206_pr))
    k(r31b) = (11580.6009523209_pr)*exp((-344074.43810964_pr)*R_T_inv + T_log *(1.05154451197381_pr))
    k(r32b) = (247.409362386158_pr)*exp((-50174.8754119911_pr)*R_T_inv + T_log *(1.155215324728_pr))
    k(r33b) = (5.78402079735284_pr)*exp((-45222.4437617245_pr)*R_T_inv + T_log *(1.59353684343726_pr))
    k(r34b) = (19.9485854062454_pr)*exp((-16832.0676782678_pr)*R_T_inv + T_log *(1.64769314674002_pr))
    k_0(FOr35b) = (6.71682647796678e+32_pr)*exp((-174383.064570079_pr)*R_T_inv + T_log *(-6.78385967515134_pr))
    k_inf(FOr35b) = (98984818799.3674_pr)*exp((-154216.184622781_pr)*R_T_inv + T_log *(0.486140314812089_pr))
    FC(FOr35b) = ((1.0_pr - 0.7507_pr)*exp(-Tloc/(98.5_pr))) + (0.7507_pr*exp(-Tloc/(1302_pr)))&
              + exp(-(4167_pr)/Tloc)
    k(r35b) = getlindratecoeff(Tloc,k_0(FOr35b),k_inf(FOr35b),FC(FOr35b), M(mM35),Ploc)
    k(r36b) = (293830170.879477_pr)*exp((-288559.055120315_pr)*R_T_inv + T_log *(-0.099870980399811_pr))
    k_0(FOr37b) = (6.05300727443141e+35_pr)*exp((-179099.928055276_pr)*R_T_inv + T_log *(-7.57952936416865_pr))
    k_inf(FOr37b) = (544770651174.815_pr)*exp((-157552.32805083_pr)*R_T_inv + T_log *(0.494470636683265_pr))
    FC(FOr37b) = ((1.0_pr - 0.9753_pr)*exp(-Tloc/(210_pr))) + (0.9753_pr*exp(-Tloc/(984_pr)))&
              + exp(-(4374_pr)/Tloc)
    k(r37b) = getlindratecoeff(Tloc,k_0(FOr37b),k_inf(FOr37b),FC(FOr37b), M(mM37),Ploc)
    k(r38b) = (0.00568838952942243_pr)*exp((-20177.2840583028_pr)*R_T_inv + T_log *(2.72659195448168_pr))
    k_0(FOr39b) = (2.23133869975817e+37_pr)*exp((-447113.941893211_pr)*R_T_inv + T_log *(-6.99089871038336_pr))
    k_inf(FOr39b) = (5.84184589554799e+19_pr)*exp((-425754.621816626_pr)*R_T_inv + T_log *(-0.900898695715539_pr))
    FC(FOr39b) = ((1.0_pr - 0.8422_pr)*exp(-Tloc/(125_pr))) + (0.8422_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r39b) = getlindratecoeff(Tloc,k_0(FOr39b),k_inf(FOr39b),FC(FOr39b), M(mM39),Ploc)
    k(r40b) = (0.177558371375167_pr)*exp((-45095.2579039898_pr)*R_T_inv + T_log *(2.19716803393134_pr))
    k(r41b) = (789.449766778048_pr)*exp((-76243.7937041867_pr)*R_T_inv + T_log *(1.54304457952421_pr))
    k_0(FOr42b) = (5.02169009475452e+16_pr)*exp((-204775.243594019_pr)*R_T_inv + T_log *(-1.5709627798132_pr))
    k_inf(FOr42b) = (1.61567422842255e+18_pr)*exp((-211888.043604252_pr)*R_T_inv + T_log *(-1.04096278177281_pr))
    FC(FOr42b) = ((1.0_pr - 0.7346_pr)*exp(-Tloc/(94_pr))) + (0.7346_pr*exp(-Tloc/(1756_pr)))&
              + exp(-(5182_pr)/Tloc)
    k(r42b) = getlindratecoeff(Tloc,k_0(FOr42b),k_inf(FOr42b),FC(FOr42b), M(mM42),Ploc)
    k(r43b) = (0.256596952930681_pr)*exp((-59324.8503507529_pr)*R_T_inv + T_log *(2.44871696144176_pr))
    k(r44b) = (42691897.977904_pr)*exp((-289157.075191882_pr)*R_T_inv + T_log *(0.281774344172452_pr))
    k(r45b) = (56431.5014460642_pr)*exp((-131030.611994645_pr)*R_T_inv + T_log *(0.520319666156543_pr))
    k(r46b) = (18703012.0941078_pr)*exp((-169690.772001502_pr)*R_T_inv + T_log *(0.520319664844188_pr))
    k(r47b) = (19237092158.1421_pr)*exp((-376764.216735932_pr)*R_T_inv + T_log *(-0.423150763255335_pr))
    k(r48b) = (111188150084.124_pr)*exp((-323053.556676437_pr)*R_T_inv + T_log *(-0.401063346960796_pr))
    k(r49b) = (5.41958765873606_pr)*exp((-86501.1619910283_pr)*R_T_inv + T_log *(2.22992493642473_pr))
    k_0(FOr50b) = (1.27288910334069e+36_pr)*exp((-399073.17676565_pr)*R_T_inv + T_log *(-6.6166391589299_pr))
    k_inf(FOr50b) = (2.97007530608937e+17_pr)*exp((-386102.77693578_pr)*R_T_inv + T_log *(-0.316639191371164_pr))
    FC(FOr50b) = ((1.0_pr - 0.2105_pr)*exp(-Tloc/(83.5_pr))) + (0.2105_pr*exp(-Tloc/(5398_pr)))&
              + exp(-(8370_pr)/Tloc)
    k(r50b) = getlindratecoeff(Tloc,k_0(FOr50b),k_inf(FOr50b),FC(FOr50b), M(mM50),Ploc)
    k(r51b) = (2.75718047847825_pr)*exp((-55415.9596902191_pr)*R_T_inv + T_log *(1.93043712045487_pr))
    k(r52b) = (1169384.80179942_pr)*exp((4952.43158327489_pr)*R_T_inv + T_log *(0.438321506017214_pr))
    k(r53b) = (6.13354027395456_pr)*exp((-69432.1670581406_pr)*R_T_inv + T_log *(1.72601084421888_pr))
    k(r54b) = (60916128.074486_pr)*exp((-106561.819711313_pr)*R_T_inv + T_log *(0.132979285857758_pr))
    k(r55b) = (32318395.6279773_pr)*exp((-430617.184464_pr)*R_T_inv + T_log *(0.442447871054845_pr))
    k(r56b) = (189.745536043524_pr)*exp((-125789.574025541_pr)*R_T_inv + T_log *(1.3878375250205_pr))
    k(r57b) = (10581.3686743477_pr)*exp((-405967.11179929_pr)*R_T_inv + T_log *(1.08458909416224_pr))
    k(r58b) = (109.363927967479_pr)*exp((-64624.6613616757_pr)*R_T_inv + T_log *(1.5807377301802_pr))
    k(r59b) = (0.0564867639159602_pr)*exp((-41317.7977537572_pr)*R_T_inv + T_log *(2.22963653564526_pr))
    k(r60b) = (0.0199764071984045_pr)*exp((-79122.4917010152_pr)*R_T_inv + T_log *(2.45021259592765_pr))
    k(r61b) = (5934815.78501792_pr)*exp((-154737.42317011_pr)*R_T_inv + T_log *(-0.238545316692981_pr))
    k(r62b) = (19174020443.0614_pr)*exp((-211765.343177818_pr)*R_T_inv + T_log *(-0.238545318168657_pr))
    k(r63b) = (24001375.9008422_pr)*exp((-234870.987967411_pr)*R_T_inv + T_log *(0.155763531811789_pr))
    k(r64b) = (2324795342.77303_pr)*exp((-102871.662545334_pr)*R_T_inv + T_log *(-0.585258316376593_pr))
    k(r65b) = (172519225118.498_pr)*exp((-358057.877697701_pr)*R_T_inv + T_log *(-0.525063711376048_pr))
    k(r66b) = (1715512.48155282_pr)*exp((-31440.0900173551_pr)*R_T_inv + T_log *(-0.3124821385899_pr))
    k(r67b) = (92851502.6547889_pr)*exp((-306714.77005697_pr)*R_T_inv + T_log *(-0.0862511767858241_pr))
    k(r68b) = (843.588339643383_pr)*exp((5069.20828815179_pr)*R_T_inv + T_log *(1.59311964330268_pr))
    k(r69b) = (198562893457.911_pr)*exp((-245945.474720976_pr)*R_T_inv + T_log *(-0.630988290198015_pr))
    k(r70b) = (371171029181.557_pr)*exp((-229686.652975132_pr)*R_T_inv + T_log *(-0.836398195438195_pr))
    k(r71b) = (2901824264758.6_pr)*exp((-255206.94230168_pr)*R_T_inv + T_log *(-0.940023891166179_pr))
    k(r72b) = (1703.61603264645_pr)*exp((-273382.236975515_pr)*R_T_inv + T_log *(0.671869960268547_pr))
    k(r73b) = (2478296.80813578_pr)*exp((-318786.84148497_pr)*R_T_inv + T_log *(0.0949568385640545_pr))
    k(r74b) = (37.1162104712948_pr)*exp((-59404.3140823287_pr)*R_T_inv + T_log *(1.70260744613117_pr))
    k(r75b) = (15127127558929.5_pr)*exp((-272778.017298047_pr)*R_T_inv + T_log *(-0.836109811037985_pr))
    k(r76b) = (15795183.5151293_pr)*exp((-40201.5112773376_pr)*R_T_inv + T_log *(-0.10788438632958_pr))
    k(r77b) = (5.65490756311237_pr)*exp((-286192.825901994_pr)*R_T_inv + T_log *(0.0102060045718551_pr))
    k(r78b) = (1533465.02741625_pr)*exp((-780819.13728662_pr)*R_T_inv + T_log *(0.429520311885355_pr))
    k(r79b) = (5471734975.89459_pr)*exp((-66845.1053157836_pr)*R_T_inv + T_log *(-0.405276931918243_pr))
    k(r80b) = (31590367.5651516_pr)*exp((-37691.1112890249_pr)*R_T_inv + T_log *(-0.107884388560024_pr))
    k(r81b) = (20988792.1360098_pr)*exp((-58945.6386350438_pr)*R_T_inv + T_log *(-0.31231066057136_pr))
    k(r82b) = (9477109.89087751_pr)*exp((-37691.1112614549_pr)*R_T_inv + T_log *(-0.107884383296248_pr))
    k(r83b) = (64041.9613293096_pr)*exp((-254475.728243319_pr)*R_T_inv + T_log *(0.586072980634045_pr))
    k(r84b) = (15181490333.5883_pr)*exp((-274.877825532746_pr)*R_T_inv + T_log *(-0.818315722140235_pr))
    k(r85b) = (6015.63700460991_pr)*exp((-254994.6591732_pr)*R_T_inv + T_log *(0.217556429584783_pr))
    k(r86b) = (0.0128806689579625_pr)*exp((-94986.7649345585_pr)*R_T_inv + T_log *(2.86430882229322_pr))
    k_0(FOr87b) = (2.94593209180122e+50_pr)*exp((-405996.213326916_pr)*R_T_inv + T_log *(-10.4654146475174_pr))
    k_inf(FOr87b) = (3.5284609821569e+22_pr)*exp((-382565.813315892_pr)*R_T_inv + T_log *(-1.76541464551057_pr))
    FC(FOr87b) = ((1.0_pr - 0.5325_pr)*exp(-Tloc/(151_pr))) + (0.5325_pr*exp(-Tloc/(1038_pr)))&
              + exp(-(4970_pr)/Tloc)
    k(r87b) = getlindratecoeff(Tloc,k_0(FOr87b),k_inf(FOr87b),FC(FOr87b), M(mM87),Ploc)
    k(r88b) = (74069188324.061_pr)*exp((-5178.23142742055_pr)*R_T_inv + T_log *(-0.784515935987335_pr))
    k(r89b) = (279052880.797679_pr)*exp((-374239.097344983_pr)*R_T_inv + T_log *(0.316437038349758_pr))
    k(r90b) = (0.00299436201457904_pr)*exp((-95799.9749593436_pr)*R_T_inv + T_log *(2.89182668231866_pr))
    k(r91b) = (2830.23393996079_pr)*exp((-43559.5343119905_pr)*R_T_inv + T_log *(0.954726884373294_pr))
    k(r92b) = (2292324.81950883_pr)*exp((-7118.87096036709_pr)*R_T_inv + T_log *(-0.976866031663883_pr))
    k(r93b) = (191027.067697868_pr)*exp((-7118.8709583992_pr)*R_T_inv + T_log *(-0.976866031247588_pr))
    k(r94b) = (3336921.14896229_pr)*exp((-141041.709316715_pr)*R_T_inv + T_log *(0.160673518251964_pr))
    k(r95b) = (6.15273414847798e-22_pr)*exp((-99948.51661372_pr)*R_T_inv + T_log *(8.40281474868019_pr))
    k(r96b) = (34468.9092810795_pr)*exp((-363663.102923362_pr)*R_T_inv + T_log *(0.630803863874551_pr))
    k(r97b) = (715761.554213764_pr)*exp((-69652.7902740043_pr)*R_T_inv + T_log *(0.0970689192465346_pr))
    k(r98b) = (66558901.675827_pr)*exp((-315723.671647508_pr)*R_T_inv + T_log *(0.101015676191457_pr))
    k(r99b) = (316700.81678875_pr)*exp((-158619.425373855_pr)*R_T_inv + T_log *(0.0652464198256787_pr))
    k(r100b) = (1997230607.05143_pr)*exp((-206577.352122738_pr)*R_T_inv + T_log *(-0.271653180067665_pr))
    k(r101b) = (370.06436527558_pr)*exp((-373203.934434293_pr)*R_T_inv + T_log *(0.922547636543731_pr))
    k(r102b) = (481.739669057493_pr)*exp((-242333.368116912_pr)*R_T_inv + T_log *(0.886778388317402_pr))
    k(r103b) = (510.338877976389_pr)*exp((-336919.367629062_pr)*R_T_inv + T_log *(1.2594472471691_pr))
    k(r104b) = (2581.16097757486_pr)*exp((-192993.949671993_pr)*R_T_inv + T_log *(0.689490238788894_pr))
    k_0(FOr105b) = (0.00816872889424917_pr)*exp((-66174.6245106584_pr)*R_T_inv + T_log *(0.888850299437564_pr))
    k_inf(FOr105b) = (1.71279808335031_pr)*exp((-80902.3045468246_pr)*R_T_inv + T_log *(0.888850292573385_pr))
    k(r105b) = getlindratecoeff(Tloc,k_0(FOr105b),k_inf(FOr105b),1.0_pr, M(mM105),Ploc)
    k(r106b) = (94448574.4668766_pr)*exp((-30690.6089777002_pr)*R_T_inv + T_log *(-0.285482326766079_pr))
    k(r107b) = (4.64773221864783e+17_pr)*exp((-306125.614105506_pr)*R_T_inv + T_log *(-1.89484238101318_pr))
    k(r108b) = (17844.9685199676_pr)*exp((-189392.895775797_pr)*R_T_inv + T_log *(0.518539726172049_pr))
    k(r109b) = (2650.21964121852_pr)*exp((-121853.848989255_pr)*R_T_inv + T_log *(0.855439333231929_pr))
    k(r110b) = (1696248728.07463_pr)*exp((-297721.791640726_pr)*R_T_inv + T_log *(-0.210116002804543_pr))
    k(r111b) = (78385458.6156598_pr)*exp((-103471.6561548_pr)*R_T_inv + T_log *(0.0772095611956895_pr))
    k(r112b) = (28285575955.7012_pr)*exp((-76920.9827783366_pr)*R_T_inv + T_log *(-0.4964683277201_pr))
    k(r113b) = (17905.4922741144_pr)*exp((-163983.609850814_pr)*R_T_inv + T_log *(1.31025414225761_pr))
    k(r114b) = (2.86083442193115_pr)*exp((-34067.5360055196_pr)*R_T_inv + T_log *(1.84043127672516_pr))
    k(r115b) = (0.190540024408296_pr)*exp((-228090.744787866_pr)*R_T_inv + T_log *(1.62678361692442_pr))
    k(r116b) = (7739168111.2882_pr)*exp((-72976.7090803851_pr)*R_T_inv + T_log *(-0.52951290840586_pr))
    k(r117b) = (849568515748704_pr)*exp((-152248.023460562_pr)*R_T_inv + T_log *(-1.54689437903518_pr))
    k(r118b) = (1009216.31499357_pr)*exp((-42906.2742410266_pr)*R_T_inv + T_log *(0.0989875769874346_pr))
    k(r119b) = (9379487441.90869_pr)*exp((-119827.257053019_pr)*R_T_inv + T_log *(-0.397480757168506_pr))
    k(r120b) = (11341193.4118403_pr)*exp((-64438.2909101801_pr)*R_T_inv + T_log *(0.114659954753329_pr))
    k(r121b) = (93.2635772048696_pr)*exp((-109134.72468505_pr)*R_T_inv + T_log *(1.64770452077274_pr))
    k(r122b) = (185.156281255755_pr)*exp((-53650.6230804689_pr)*R_T_inv + T_log *(-0.023528967901363_pr))
    k(r123b) = (3705313.86265956_pr)*exp((-49114.0546899451_pr)*R_T_inv + T_log *(0.294515191764427_pr))
    k(r124b) = (12328793.8277553_pr)*exp((-465547.9424885_pr)*R_T_inv + T_log *(0.472740430152157_pr))
    k(r125b) = (328459447921874_pr)*exp((-208768.574317172_pr)*R_T_inv + T_log *(-1.23575538891207_pr))
    k(r126b) = (599686.764099662_pr)*exp((-220800.808790735_pr)*R_T_inv + T_log *(0.286352338632028_pr))
    k(r127b) = (21228.0239640656_pr)*exp((-229822.665551505_pr)*R_T_inv + T_log *(1.02202469887733_pr))
    k(r128b) = (2.24135753264699_pr)*exp((-284979.099226512_pr)*R_T_inv + T_log *(2.23506928396711_pr))
    k(r129b) = (1171171.05250748_pr)*exp((-52096.824046566_pr)*R_T_inv + T_log *(0.0532949374131235_pr))
    k(r130b) = (3.5675387850885e-06_pr)*exp((-156228.709481558_pr)*R_T_inv + T_log *(3.31343099013666_pr))
    k(r131b) = (2213.44442431723_pr)*exp((-24182.4018487794_pr)*R_T_inv + T_log *(1.43180309412836_pr))
    k(r132b) = (37.0618409156963_pr)*exp((-100447.858407247_pr)*R_T_inv + T_log *(1.52670938755944_pr))
    k(r133b) = (198064.994944848_pr)*exp((-77124.7780512641_pr)*R_T_inv + T_log *(0.506249487881391_pr))
    k(r134b) = (65918238673.7806_pr)*exp((-306114.129701307_pr)*R_T_inv + T_log *(-0.492734836216117_pr))
    k(r135b) = (2.29707976912738e+17_pr)*exp((-390894.139682752_pr)*R_T_inv + T_log *(-2.42481746122604_pr))
    k(r136b) = (25511621437.5286_pr)*exp((-315090.041262718_pr)*R_T_inv + T_log *(-1.00152684634946_pr))
    k(r137b) = (2.41885301211438e+17_pr)*exp((-428585.250934983_pr)*R_T_inv + T_log *(-2.5327018427609_pr))
    k(r138b) = (415803618.962138_pr)*exp((-465148.080636379_pr)*R_T_inv + T_log *(0.0189103348074121_pr))
    k(r139b) = (19.197146860908_pr)*exp((-203555.070849181_pr)*R_T_inv + T_log *(1.72170087625737_pr))
    k(r140b) = (2.38880954946344e-05_pr)*exp((-34715.0568300791_pr)*R_T_inv + T_log *(3.00773402619827_pr))
    k(r141b) = (2105988.98775727_pr)*exp((-153745.196622529_pr)*R_T_inv + T_log *(0.212713296597358_pr))
    k(r142b) = (83.719835460184_pr)*exp((-113795.936552612_pr)*R_T_inv + T_log *(1.8649062945505_pr))
    k(r143b) = (119424859.319102_pr)*exp((-484270.919378232_pr)*R_T_inv + T_log *(-0.135787486725005_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) * c(sCH3) * c(sO)
    w(r6f) = k(r6f) * c(sCH4) * c(sO)
    w(r7f) = k(r7f) * cqss(sqssHCO) * c(sO)
    w(r8f) = k(r8f) * c(sCH2O) * c(sO)
    w(r9f) = k(r9f) * c(sC2H2) * c(sO)
    w(r10f) = k(r10f) * c(sC2H4) * c(sO)
    w(r11f) = k(r11f) * cqss(sqssC2H5) * c(sO)
    w(r12f) = k(r12f) * c(sC2H6) * c(sO)
    w(r13f) = k(r13f) * c(sCH2O) * c(sO2)
    w(r14f) = k(r14f) * c(sH) * c(sO2) * m(mM14)
    w(r15f) = k(r15f) * c(sH) * c(sO2)
    w(r16f) = k(r16f) * c(sH) * c(sH2O) * c(sO2)
    w(r17f) = k(r17f) * c(sH) * c(sN2) * c(sO2)
    w(r18f) = k(r18f) * c(sH) * c(sO2)
    w(r19f) = k(r19f) * c(sH)**2.0_pr * c(sH2O)
    w(r20f) = k(r20f) * c(sH) * c(sOH) * m(mM20)
    w(r21f) = k(r21f) * c(sH) * c(sHO2)
    w(r22f) = k(r22f) * c(sH) * c(sHO2)
    w(r23f) = k(r23f) * c(sH) * c(sHO2)
    w(r24f) = k(r24f) * cqss(sqssCH2) * c(sH)
    w(r25f) = k(r25f) * cqss(sqssCH2GSG) * c(sH)
    w(r26f) = k(r26f) * c(sCH3) * c(sH)
    w(r27f) = k(r27f) * c(sCH4) * c(sH)
    w(r28f) = k(r28f) * c(sH) * cqss(sqssHCO)
    w(r29f) = k(r29f) * c(sCH2O) * c(sH)
    w(r30f) = k(r30f) * c(sCH2O) * c(sH)
    w(r31f) = k(r31f) * cqss(sqssCH3O) * c(sH)
    w(r32f) = k(r32f) * cqss(sqssCH3O) * c(sH)
    w(r33f) = k(r33f) * cqss(sqssCH3O) * c(sH)
    w(r34f) = k(r34f) * c(sCH3OH) * c(sH)
    w(r35f) = k(r35f) * c(sC2H2) * c(sH)
    w(r36f) = k(r36f) * cqss(sqssC2H3) * c(sH)
    w(r37f) = k(r37f) * c(sC2H4) * c(sH)
    w(r38f) = k(r38f) * c(sC2H4) * c(sH)
    w(r39f) = k(r39f) * cqss(sqssC2H5) * c(sH)
    w(r40f) = k(r40f) * c(sC2H6) * c(sH)
    w(r41f) = k(r41f) * c(sH2) * c(sOH)
    w(r42f) = k(r42f) * c(sOH)**2.0_pr
    w(r43f) = k(r43f) * c(sOH)**2.0_pr
    w(r44f) = k(r44f) * c(sHO2) * c(sOH)
    w(r45f) = k(r45f) * c(sH2O2) * c(sOH)
    w(r46f) = k(r46f) * c(sH2O2) * c(sOH)
    w(r47f) = k(r47f) * cqss(sqssCH) * c(sOH)
    w(r48f) = k(r48f) * cqss(sqssCH2) * c(sOH)
    w(r49f) = k(r49f) * cqss(sqssCH2) * c(sOH)
    w(r50f) = k(r50f) * c(sCH3) * c(sOH)
    w(r51f) = k(r51f) * c(sCH3) * c(sOH)
    w(r52f) = k(r52f) * c(sCH3) * c(sOH)
    w(r53f) = k(r53f) * c(sCH4) * c(sOH)
    w(r54f) = k(r54f) * c(sCO) * c(sOH)
    w(r55f) = k(r55f) * cqss(sqssHCO) * c(sOH)
    w(r56f) = k(r56f) * c(sCH2O) * c(sOH)
    w(r57f) = k(r57f) * cqss(sqssCH3O) * c(sOH)
    w(r58f) = k(r58f) * c(sCH3OH) * c(sOH)
    w(r59f) = k(r59f) * c(sC2H4) * c(sOH)
    w(r60f) = k(r60f) * c(sC2H6) * c(sOH)
    w(r61f) = k(r61f) * c(sHO2)**2.0_pr
    w(r62f) = k(r62f) * c(sHO2)**2.0_pr
    w(r63f) = k(r63f) * c(sCH3) * c(sHO2)
    w(r64f) = k(r64f) * c(sCH3) * c(sHO2)
    w(r65f) = k(r65f) * c(sCO) * c(sHO2)
    w(r66f) = k(r66f) * c(sCH2O) * c(sHO2)
    w(r67f) = k(r67f) * cqss(sqssCH) * c(sO2)
    w(r68f) = k(r68f) * cqss(sqssCH) * c(sH2)
    w(r69f) = k(r69f) * cqss(sqssCH) * c(sH2O)
    w(r70f) = k(r70f) * cqss(sqssCH) * c(sCH3)
    w(r71f) = k(r71f) * cqss(sqssCH) * c(sCH4)
    w(r72f) = k(r72f) * cqss(sqssCH) * c(sCO2)
    w(r73f) = k(r73f) * cqss(sqssCH2) * c(sO2)
    w(r74f) = k(r74f) * cqss(sqssCH2) * c(sH2)
    w(r75f) = k(r75f) * cqss(sqssCH2) * c(sCH3)
    w(r76f) = k(r76f) * cqss(sqssCH2GSG) * c(sN2)
    w(r77f) = k(r77f) * cqss(sqssCH2GSG) * c(sO2)
    w(r78f) = k(r78f) * cqss(sqssCH2GSG) * c(sO2)
    w(r79f) = k(r79f) * cqss(sqssCH2GSG) * c(sH2)
    w(r80f) = k(r80f) * cqss(sqssCH2GSG) * c(sH2O)
    w(r81f) = k(r81f) * cqss(sqssCH2GSG) * c(sCH4)
    w(r82f) = k(r82f) * cqss(sqssCH2GSG) * c(sCO)
    w(r83f) = k(r83f) * cqss(sqssCH2GSG) * c(sCO2)
    w(r84f) = k(r84f) * c(sCH3) * c(sO2)
    w(r85f) = k(r85f) * c(sCH3) * c(sO2)
    w(r86f) = k(r86f) * c(sCH3) * c(sH2O2)
    w(r87f) = k(r87f) * c(sCH3)**2.0_pr
    w(r88f) = k(r88f) * c(sCH3)**2.0_pr
    w(r89f) = k(r89f) * c(sCH3) * cqss(sqssHCO)
    w(r90f) = k(r90f) * c(sCH2O) * c(sCH3)
    w(r91f) = k(r91f) * c(sCH3) * c(sCH3OH)
    w(r92f) = k(r92f) * c(sH2O) * cqss(sqssHCO)
    w(r93f) = k(r93f) * cqss(sqssHCO) * m(mM93)
    w(r94f) = k(r94f) * cqss(sqssHCO) * c(sO2)
    w(r95f) = k(r95f) * cqss(sqssCH3O) * c(sO2)
    w(r96f) = k(r96f) * cqss(sqssC2H3) * c(sO2)
    w(r97f) = k(r97f) * cqss(sqssC2H5) * c(sO2)
    w(r98f) = k(r98f) * cqss(sqssN) * c(sNO)
    w(r99f) = k(r99f) * cqss(sqssN) * c(sO2)
    w(r100f) = k(r100f) * cqss(sqssN) * c(sOH)
    w(r101f) = k(r101f) * c(sN2O) * c(sO)
    w(r102f) = k(r102f) * c(sN2O) * c(sO)
    w(r103f) = k(r103f) * c(sH) * c(sN2O)
    w(r104f) = k(r104f) * c(sN2O) * c(sOH)
    w(r105f) = k(r105f) * c(sN2O)
    w(r106f) = k(r106f) * c(sHO2) * c(sNO)
    w(r107f) = k(r107f) * c(sNO) * c(sO) * m(mM107)
    w(r108f) = k(r108f) * c(sNO2) * c(sO)
    w(r109f) = k(r109f) * c(sH) * c(sNO2)
    w(r110f) = k(r110f) * cqss(sqssNH) * c(sO)
    w(r111f) = k(r111f) * c(sH) * cqss(sqssNH)
    w(r112f) = k(r112f) * cqss(sqssNH) * c(sOH)
    w(r113f) = k(r113f) * cqss(sqssNH) * c(sOH)
    w(r114f) = k(r114f) * cqss(sqssNH) * c(sO2)
    w(r115f) = k(r115f) * cqss(sqssNH) * c(sO2)
    w(r116f) = k(r116f) * c(sH2O) * cqss(sqssNH)
    w(r117f) = k(r117f) * cqss(sqssNH) * c(sNO)
    w(r118f) = k(r118f) * cqss(sqssNH2) * c(sO)
    w(r119f) = k(r119f) * cqss(sqssNH2) * c(sO)
    w(r120f) = k(r120f) * c(sH) * cqss(sqssNH2)
    w(r121f) = k(r121f) * cqss(sqssNH2) * c(sOH)
    w(r122f) = k(r122f) * cqss(sqssNNH) * m(mM122)
    w(r123f) = k(r123f) * cqss(sqssNNH) * c(sO)
    w(r124f) = k(r124f) * c(sH) * cqss(sqssNNH)
    w(r125f) = k(r125f) * c(sH) * c(sNO) * m(mM125)
    w(r126f) = k(r126f) * cqss(sqssHNO) * c(sO)
    w(r127f) = k(r127f) * c(sH) * cqss(sqssHNO)
    w(r128f) = k(r128f) * cqss(sqssHNO) * c(sOH)
    w(r129f) = k(r129f) * cqss(sqssHNO) * c(sO2)
    w(r130f) = k(r130f) * c(sHCN) * c(sO)
    w(r131f) = k(r131f) * c(sHCN) * c(sOH)
    w(r132f) = k(r132f) * c(sHCN) * c(sOH)
    w(r133f) = k(r133f) * cqss(sqssCH) * c(sN2)
    w(r134f) = k(r134f) * cqss(sqssCH) * c(sNO)
    w(r135f) = k(r135f) * cqss(sqssCH2) * c(sNO)
    w(r136f) = k(r136f) * cqss(sqssCH2) * c(sNO)
    w(r137f) = k(r137f) * cqss(sqssCH2GSG) * c(sNO)
    w(r138f) = k(r138f) * c(sCH3) * c(sNO)
    w(r139f) = k(r139f) * cqss(sqssHNCO) * c(sO)
    w(r140f) = k(r140f) * c(sH) * cqss(sqssHNCO)
    w(r141f) = k(r141f) * cqss(sqssHNCO) * c(sOH)
    w(r142f) = k(r142f) * c(sH) * cqss(sqssHOCN)
    w(r143f) = k(r143f) * c(sCH3) * cqss(sqssN)
    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(sCH2O) * c(sH)
    w(r6b) = k(r6b) * c(sCH3) * c(sOH)
    w(r7b) = k(r7b) * c(sCO2) * c(sH)
    w(r8b) = k(r8b) * cqss(sqssHCO) * c(sOH)
    w(r9b) = k(r9b) * cqss(sqssCH2) * c(sCO)
    w(r10b) = k(r10b) * c(sCH3) * cqss(sqssHCO)
    w(r11b) = k(r11b) * c(sCH2O) * c(sCH3)
    w(r12b) = k(r12b) * cqss(sqssC2H5) * c(sOH)
    w(r13b) = k(r13b) * cqss(sqssHCO) * c(sHO2)
    w(r14b) = k(r14b) * c(sHO2) * m(mM14)
    w(r15b) = k(r15b) * c(sHO2) * c(sO2)
    w(r16b) = k(r16b) * c(sH2O) * c(sHO2)
    w(r17b) = k(r17b) * c(sHO2) * c(sN2)
    w(r18b) = k(r18b) * c(sO) * c(sOH)
    w(r19b) = k(r19b) * c(sH2) * c(sH2O)
    w(r20b) = k(r20b) * c(sH2O) * m(mM20)
    w(r21b) = k(r21b) * c(sH2O) * c(sO)
    w(r22b) = k(r22b) * c(sH2) * c(sO2)
    w(r23b) = k(r23b) * c(sOH)**2.0_pr
    w(r24b) = k(r24b) * c(sCH3)
    w(r25b) = k(r25b) * cqss(sqssCH) * c(sH2)
    w(r26b) = k(r26b) * c(sCH4)
    w(r27b) = k(r27b) * c(sCH3) * c(sH2)
    w(r28b) = k(r28b) * c(sCO) * c(sH2)
    w(r29b) = k(r29b) * cqss(sqssCH3O)
    w(r30b) = k(r30b) * c(sH2) * cqss(sqssHCO)
    w(r31b) = k(r31b) * c(sCH2O) * c(sH2)
    w(r32b) = k(r32b) * c(sCH3) * c(sOH)
    w(r33b) = k(r33b) * cqss(sqssCH2GSG) * c(sH2O)
    w(r34b) = k(r34b) * cqss(sqssCH3O) * c(sH2)
    w(r35b) = k(r35b) * cqss(sqssC2H3)
    w(r36b) = k(r36b) * c(sC2H2) * c(sH2)
    w(r37b) = k(r37b) * cqss(sqssC2H5)
    w(r38b) = k(r38b) * cqss(sqssC2H3) * c(sH2)
    w(r39b) = k(r39b) * c(sC2H6)
    w(r40b) = k(r40b) * cqss(sqssC2H5) * c(sH2)
    w(r41b) = k(r41b) * c(sH) * c(sH2O)
    w(r42b) = k(r42b) * c(sH2O2)
    w(r43b) = k(r43b) * c(sH2O) * c(sO)
    w(r44b) = k(r44b) * c(sH2O) * c(sO2)
    w(r45b) = k(r45b) * c(sH2O) * c(sHO2)
    w(r46b) = k(r46b) * c(sH2O) * c(sHO2)
    w(r47b) = k(r47b) * c(sH) * cqss(sqssHCO)
    w(r48b) = k(r48b) * c(sCH2O) * c(sH)
    w(r49b) = k(r49b) * cqss(sqssCH) * c(sH2O)
    w(r50b) = k(r50b) * c(sCH3OH)
    w(r51b) = k(r51b) * cqss(sqssCH2) * c(sH2O)
    w(r52b) = k(r52b) * cqss(sqssCH2GSG) * c(sH2O)
    w(r53b) = k(r53b) * c(sCH3) * c(sH2O)
    w(r54b) = k(r54b) * c(sCO2) * c(sH)
    w(r55b) = k(r55b) * c(sCO) * c(sH2O)
    w(r56b) = k(r56b) * c(sH2O) * cqss(sqssHCO)
    w(r57b) = k(r57b) * c(sCH2O) * c(sH2O)
    w(r58b) = k(r58b) * cqss(sqssCH3O) * c(sH2O)
    w(r59b) = k(r59b) * cqss(sqssC2H3) * c(sH2O)
    w(r60b) = k(r60b) * cqss(sqssC2H5) * c(sH2O)
    w(r61b) = k(r61b) * c(sH2O2) * c(sO2)
    w(r62b) = k(r62b) * c(sH2O2) * c(sO2)
    w(r63b) = k(r63b) * c(sCH4) * c(sO2)
    w(r64b) = k(r64b) * cqss(sqssCH3O) * c(sOH)
    w(r65b) = k(r65b) * c(sCO2) * c(sOH)
    w(r66b) = k(r66b) * c(sH2O2) * cqss(sqssHCO)
    w(r67b) = k(r67b) * cqss(sqssHCO) * c(sO)
    w(r68b) = k(r68b) * cqss(sqssCH2) * c(sH)
    w(r69b) = k(r69b) * c(sCH2O) * c(sH)
    w(r70b) = k(r70b) * cqss(sqssC2H3) * c(sH)
    w(r71b) = k(r71b) * c(sC2H4) * c(sH)
    w(r72b) = k(r72b) * c(sCO) * cqss(sqssHCO)
    w(r73b) = k(r73b) * cqss(sqssHCO) * c(sOH)
    w(r74b) = k(r74b) * c(sCH3) * c(sH)
    w(r75b) = k(r75b) * c(sC2H4) * c(sH)
    w(r76b) = k(r76b) * cqss(sqssCH2) * c(sN2)
    w(r77b) = k(r77b) * c(sCO) * c(sH) * c(sOH)
    w(r78b) = k(r78b) * c(sCO) * c(sH2O)
    w(r79b) = k(r79b) * c(sCH3) * c(sH)
    w(r80b) = k(r80b) * cqss(sqssCH2) * c(sH2O)
    w(r81b) = k(r81b) * c(sCH3)**2.0_pr
    w(r82b) = k(r82b) * cqss(sqssCH2) * c(sCO)
    w(r83b) = k(r83b) * c(sCH2O) * c(sCO)
    w(r84b) = k(r84b) * cqss(sqssCH3O) * c(sO)
    w(r85b) = k(r85b) * c(sCH2O) * c(sOH)
    w(r86b) = k(r86b) * c(sCH4) * c(sHO2)
    w(r87b) = k(r87b) * c(sC2H6)
    w(r88b) = k(r88b) * cqss(sqssC2H5) * c(sH)
    w(r89b) = k(r89b) * c(sCH4) * c(sCO)
    w(r90b) = k(r90b) * c(sCH4) * cqss(sqssHCO)
    w(r91b) = k(r91b) * cqss(sqssCH3O) * c(sCH4)
    w(r92b) = k(r92b) * c(sCO) * c(sH) * c(sH2O)
    w(r93b) = k(r93b) * c(sCO) * c(sH) * m(mM93)
    w(r94b) = k(r94b) * c(sCO) * c(sHO2)
    w(r95b) = k(r95b) * c(sCH2O) * c(sHO2)
    w(r96b) = k(r96b) * c(sCH2O) * cqss(sqssHCO)
    w(r97b) = k(r97b) * c(sC2H4) * c(sHO2)
    w(r98b) = k(r98b) * c(sN2) * c(sO)
    w(r99b) = k(r99b) * c(sNO) * c(sO)
    w(r100b) = k(r100b) * c(sH) * c(sNO)
    w(r101b) = k(r101b) * c(sN2) * c(sO2)
    w(r102b) = k(r102b) * c(sNO)**2.0_pr
    w(r103b) = k(r103b) * c(sN2) * c(sOH)
    w(r104b) = k(r104b) * c(sHO2) * c(sN2)
    w(r105b) = k(r105b) * c(sN2) * c(sO)
    w(r106b) = k(r106b) * c(sNO2) * c(sOH)
    w(r107b) = k(r107b) * c(sNO2) * m(mM107)
    w(r108b) = k(r108b) * c(sNO) * c(sO2)
    w(r109b) = k(r109b) * c(sNO) * c(sOH)
    w(r110b) = k(r110b) * c(sH) * c(sNO)
    w(r111b) = k(r111b) * c(sH2) * cqss(sqssN)
    w(r112b) = k(r112b) * c(sH) * cqss(sqssHNO)
    w(r113b) = k(r113b) * c(sH2O) * cqss(sqssN)
    w(r114b) = k(r114b) * cqss(sqssHNO) * c(sO)
    w(r115b) = k(r115b) * c(sNO) * c(sOH)
    w(r116b) = k(r116b) * c(sH2) * cqss(sqssHNO)
    w(r117b) = k(r117b) * c(sH) * c(sN2O)
    w(r118b) = k(r118b) * cqss(sqssNH) * c(sOH)
    w(r119b) = k(r119b) * c(sH) * cqss(sqssHNO)
    w(r120b) = k(r120b) * c(sH2) * cqss(sqssNH)
    w(r121b) = k(r121b) * c(sH2O) * cqss(sqssNH)
    w(r122b) = k(r122b) * c(sH) * c(sN2) * m(mM122)
    w(r123b) = k(r123b) * cqss(sqssNH) * c(sNO)
    w(r124b) = k(r124b) * c(sH2) * c(sN2)
    w(r125b) = k(r125b) * cqss(sqssHNO) * m(mM125)
    w(r126b) = k(r126b) * c(sNO) * c(sOH)
    w(r127b) = k(r127b) * c(sH2) * c(sNO)
    w(r128b) = k(r128b) * c(sH2O) * c(sNO)
    w(r129b) = k(r129b) * c(sHO2) * c(sNO)
    w(r130b) = k(r130b) * c(sCO) * cqss(sqssNH)
    w(r131b) = k(r131b) * c(sH) * cqss(sqssHOCN)
    w(r132b) = k(r132b) * c(sH) * cqss(sqssHNCO)
    w(r133b) = k(r133b) * c(sHCN) * cqss(sqssN)
    w(r134b) = k(r134b) * c(sHCN) * c(sO)
    w(r135b) = k(r135b) * c(sH) * cqss(sqssHNCO)
    w(r136b) = k(r136b) * c(sHCN) * c(sOH)
    w(r137b) = k(r137b) * c(sH) * cqss(sqssHNCO)
    w(r138b) = k(r138b) * c(sH2O) * c(sHCN)
    w(r139b) = k(r139b) * c(sCO2) * cqss(sqssNH)
    w(r140b) = k(r140b) * c(sCO) * cqss(sqssNH2)
    w(r141b) = k(r141b) * c(sCO2) * cqss(sqssNH2)
    w(r142b) = k(r142b) * c(sH) * cqss(sqssHNCO)
    w(r143b) = k(r143b) * c(sH2) * c(sHCN)

    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(r19f) &
         - w(r19b) &
         + w(r22f) &
         - w(r22b) &
         + w(r25f) &
         - w(r25b) &
         + w(r27f) &
         - w(r27b) &
         + w(r28f) &
         - w(r28b) &
         + w(r30f) &
         - w(r30b) &
         + w(r31f) &
         - w(r31b) &
         + w(r34f) &
         - w(r34b) &
         + w(r36f) &
         - w(r36b) &
         + w(r38f) &
         - w(r38b) &
         + w(r40f) &
         - w(r40b) &
         - w(r41f) &
         + w(r41b) &
         - w(r68f) &
         + w(r68b) &
         - w(r74f) &
         + w(r74b) &
         - w(r79f) &
         + w(r79b) &
         + w(r111f) &
         - w(r111b) &
         + w(r116f) &
         - w(r116b) &
         + w(r120f) &
         - w(r120b) &
         + w(r124f) &
         - w(r124b) &
         + w(r127f) &
         - w(r127b) &
         + w(r143f) &
         - w(r143b)

    cdot(sH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r3f) &
         - w(r3b) &
         + w(r4f) &
         - w(r4b) &
         + w(r5f) &
         - w(r5b) &
         + w(r7f) &
         - w(r7b) &
         - w(r14f) &
         + w(r14b) &
         - w(r15f) &
         + w(r15b) &
         - w(r16f) &
         + w(r16b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         -2.0_pr * w(r19f) &
         + 2.0_pr * w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) &
         - w(r25f) &
         + w(r25b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         - w(r28f) &
         + w(r28b) &
         - w(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(r47f) &
         - w(r47b) &
         + w(r48f) &
         - w(r48b) &
         + w(r54f) &
         - w(r54b) &
         + w(r68f) &
         - w(r68b) &
         + w(r69f) &
         - w(r69b) &
         + w(r70f) &
         - w(r70b) &
         + w(r71f) &
         - w(r71b) &
         + w(r74f) &
         - w(r74b) &
         + w(r75f) &
         - w(r75b) &
         + w(r77f) &
         - w(r77b) &
         + w(r79f) &
         - w(r79b) &
         + w(r88f) &
         - w(r88b) &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         + w(r100f) &
         - w(r100b) &
         - w(r103f) &
         + w(r103b) &
         - w(r109f) &
         + w(r109b) &
         + w(r110f) &
         - w(r110b) &
         - w(r111f) &
         + w(r111b) &
         + w(r112f) &
         - w(r112b) &
         + w(r117f) &
         - w(r117b) &
         + w(r119f) &
         - w(r119b) &
         - w(r120f) &
         + w(r120b) &
         + w(r122f) &
         - w(r122b) &
         - w(r124f) &
         + w(r124b) &
         - w(r125f) &
         + w(r125b) &
         - w(r127f) &
         + w(r127b) &
         + w(r131f) &
         - w(r131b) &
         + w(r132f) &
         - w(r132b) &
         + w(r135f) &
         - w(r135b) &
         + w(r137f) &
         - w(r137b) &
         - w(r140f) &
         + w(r140b)

    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(r18f) &
         - w(r18b) &
         + w(r21f) &
         - w(r21b) &
         + w(r43f) &
         - w(r43b) &
         + w(r67f) &
         - w(r67b) &
         + w(r84f) &
         - w(r84b) &
         + w(r98f) &
         - w(r98b) &
         + w(r99f) &
         - w(r99b) &
         - w(r101f) &
         + w(r101b) &
         - w(r102f) &
         + w(r102b) &
         + w(r105f) &
         - w(r105b) &
         - w(r107f) &
         + w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r110f) &
         + w(r110b) &
         + w(r114f) &
         - w(r114b) &
         - w(r118f) &
         + w(r118b) &
         - w(r119f) &
         + w(r119b) &
         - w(r123f) &
         + w(r123b) &
         - w(r126f) &
         + w(r126b) &
         - w(r130f) &
         + w(r130b) &
         + w(r134f) &
         - w(r134b) &
         - w(r139f) &
         + w(r139b)

    cdot(sO2) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         - 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(r22f) &
         - w(r22b) &
         + w(r44f) &
         - w(r44b) &
         + w(r61f) &
         - w(r61b) &
         + w(r62f) &
         - w(r62b) &
         + w(r63f) &
         - w(r63b) &
         - w(r67f) &
         + w(r67b) &
         - w(r73f) &
         + w(r73b) &
         - w(r77f) &
         + w(r77b) &
         - w(r78f) &
         + w(r78b) &
         - w(r84f) &
         + w(r84b) &
         - w(r85f) &
         + w(r85b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         - w(r96f) &
         + w(r96b) &
         - w(r97f) &
         + w(r97b) &
         - w(r99f) &
         + w(r99b) &
         + w(r101f) &
         - w(r101b) &
         + w(r108f) &
         - w(r108b) &
         - w(r114f) &
         + w(r114b) &
         - w(r115f) &
         + w(r115b) &
         - w(r129f) &
         + w(r129b)

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r6f) &
         - w(r6b) &
         + w(r8f) &
         - w(r8b) &
         + w(r12f) &
         - w(r12b) &
         + w(r18f) &
         - w(r18b) &
         - w(r20f) &
         + w(r20b) &
         + 2.0_pr * w(r23f) &
         - 2.0_pr * w(r23b) &
         + w(r32f) &
         - w(r32b) &
         - w(r41f) &
         + w(r41b) &
         -2.0_pr * w(r42f) &
         + 2.0_pr * w(r42b) &
         -2.0_pr * w(r43f) &
         + 2.0_pr * w(r43b) &
         - w(r44f) &
         + w(r44b) &
         - w(r45f) &
         + w(r45b) &
         - w(r46f) &
         + w(r46b) &
         - w(r47f) &
         + w(r47b) &
         - w(r48f) &
         + w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         - w(r54f) &
         + w(r54b) &
         - w(r55f) &
         + w(r55b) &
         - w(r56f) &
         + w(r56b) &
         - w(r57f) &
         + w(r57b) &
         - w(r58f) &
         + w(r58b) &
         - w(r59f) &
         + w(r59b) &
         - w(r60f) &
         + w(r60b) &
         + w(r64f) &
         - w(r64b) &
         + w(r65f) &
         - w(r65b) &
         + w(r73f) &
         - w(r73b) &
         + w(r77f) &
         - w(r77b) &
         + w(r85f) &
         - w(r85b) &
         - w(r100f) &
         + w(r100b) &
         + w(r103f) &
         - w(r103b) &
         - w(r104f) &
         + w(r104b) &
         + w(r106f) &
         - w(r106b) &
         + w(r109f) &
         - w(r109b) &
         - w(r112f) &
         + w(r112b) &
         - w(r113f) &
         + w(r113b) &
         + w(r115f) &
         - w(r115b) &
         + w(r118f) &
         - w(r118b) &
         - w(r121f) &
         + w(r121b) &
         + w(r126f) &
         - w(r126b) &
         - w(r128f) &
         + w(r128b) &
         - w(r131f) &
         + w(r131b) &
         - w(r132f) &
         + w(r132b) &
         + w(r136f) &
         - w(r136b) &
         - w(r141f) &
         + w(r141b)

    cdot(sH2O) = 0.0_pr &
         + w(r20f) &
         - w(r20b) &
         + w(r21f) &
         - w(r21b) &
         + w(r33f) &
         - w(r33b) &
         + w(r41f) &
         - w(r41b) &
         + w(r43f) &
         - w(r43b) &
         + w(r44f) &
         - w(r44b) &
         + w(r45f) &
         - w(r45b) &
         + w(r46f) &
         - w(r46b) &
         + w(r49f) &
         - w(r49b) &
         + w(r51f) &
         - w(r51b) &
         + w(r52f) &
         - w(r52b) &
         + w(r53f) &
         - w(r53b) &
         + w(r55f) &
         - w(r55b) &
         + w(r56f) &
         - w(r56b) &
         + w(r57f) &
         - w(r57b) &
         + w(r58f) &
         - w(r58b) &
         + w(r59f) &
         - w(r59b) &
         + w(r60f) &
         - w(r60b) &
         - w(r69f) &
         + w(r69b) &
         + w(r78f) &
         - w(r78b) &
         + w(r113f) &
         - w(r113b) &
         - w(r116f) &
         + w(r116b) &
         + w(r121f) &
         - w(r121b) &
         + w(r128f) &
         - w(r128b) &
         + w(r138f) &
         - w(r138b)

    cdot(sHO2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         + w(r15f) &
         - w(r15b) &
         + w(r16f) &
         - w(r16b) &
         + w(r17f) &
         - w(r17b) &
         - w(r21f) &
         + w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         - w(r44f) &
         + w(r44b) &
         + w(r45f) &
         - w(r45b) &
         + w(r46f) &
         - w(r46b) &
         -2.0_pr * w(r61f) &
         + 2.0_pr * w(r61b) &
         -2.0_pr * w(r62f) &
         + 2.0_pr * w(r62b) &
         - w(r63f) &
         + w(r63b) &
         - w(r64f) &
         + w(r64b) &
         - w(r65f) &
         + w(r65b) &
         - w(r66f) &
         + w(r66b) &
         + w(r86f) &
         - w(r86b) &
         + w(r94f) &
         - w(r94b) &
         + w(r95f) &
         - w(r95b) &
         + w(r97f) &
         - w(r97b) &
         + w(r104f) &
         - w(r104b) &
         - w(r106f) &
         + w(r106b) &
         + w(r129f) &
         - w(r129b)

    cdot(sH2O2) = 0.0_pr &
         + w(r42f) &
         - w(r42b) &
         - w(r45f) &
         + w(r45b) &
         - w(r46f) &
         + w(r46b) &
         + w(r61f) &
         - w(r61b) &
         + w(r62f) &
         - w(r62b) &
         + w(r66f) &
         - w(r66b) &
         - w(r86f) &
         + w(r86b)

    cdot(sCH3) = 0.0_pr &
         - w(r5f) &
         + w(r5b) &
         + w(r6f) &
         - w(r6b) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         + w(r24f) &
         - w(r24b) &
         - w(r26f) &
         + w(r26b) &
         + w(r27f) &
         - w(r27b) &
         + w(r32f) &
         - w(r32b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         + w(r53f) &
         - w(r53b) &
         - w(r63f) &
         + w(r63b) &
         - w(r64f) &
         + w(r64b) &
         - w(r70f) &
         + w(r70b) &
         + w(r74f) &
         - w(r74b) &
         - w(r75f) &
         + w(r75b) &
         + w(r79f) &
         - w(r79b) &
         + 2.0_pr * w(r81f) &
         - 2.0_pr * w(r81b) &
         - w(r84f) &
         + w(r84b) &
         - w(r85f) &
         + w(r85b) &
         - w(r86f) &
         + w(r86b) &
         -2.0_pr * w(r87f) &
         + 2.0_pr * w(r87b) &
         -2.0_pr * w(r88f) &
         + 2.0_pr * w(r88b) &
         - w(r89f) &
         + w(r89b) &
         - w(r90f) &
         + w(r90b) &
         - w(r91f) &
         + w(r91b) &
         - w(r138f) &
         + w(r138b) &
         - w(r143f) &
         + w(r143b)

    cdot(sCH4) = 0.0_pr &
         - w(r6f) &
         + w(r6b) &
         + w(r26f) &
         - w(r26b) &
         - w(r27f) &
         + w(r27b) &
         - w(r53f) &
         + w(r53b) &
         + w(r63f) &
         - w(r63b) &
         - w(r71f) &
         + w(r71b) &
         - w(r81f) &
         + w(r81b) &
         + w(r86f) &
         - w(r86b) &
         + w(r89f) &
         - w(r89b) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         - w(r91b)

    cdot(sCO) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         + w(r9f) &
         - w(r9b) &
         + w(r28f) &
         - w(r28b) &
         - w(r54f) &
         + w(r54b) &
         + w(r55f) &
         - w(r55b) &
         - w(r65f) &
         + w(r65b) &
         + w(r72f) &
         - w(r72b) &
         + w(r77f) &
         - w(r77b) &
         + w(r78f) &
         - w(r78b) &
         + w(r83f) &
         - w(r83b) &
         + w(r89f) &
         - w(r89b) &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         + w(r94f) &
         - w(r94b) &
         + w(r130f) &
         - w(r130b) &
         + w(r140f) &
         - w(r140b)

    cdot(sCO2) = 0.0_pr &
         + w(r7f) &
         - w(r7b) &
         + w(r54f) &
         - w(r54b) &
         + w(r65f) &
         - w(r65b) &
         - w(r72f) &
         + w(r72b) &
         - w(r83f) &
         + w(r83b) &
         + w(r139f) &
         - w(r139b) &
         + w(r141f) &
         - w(r141b)

    cdot(sCH2O) = 0.0_pr &
         + w(r5f) &
         - w(r5b) &
         - w(r8f) &
         + w(r8b) &
         + w(r11f) &
         - w(r11b) &
         - w(r13f) &
         + w(r13b) &
         - w(r29f) &
         + w(r29b) &
         - w(r30f) &
         + w(r30b) &
         + w(r31f) &
         - w(r31b) &
         + w(r48f) &
         - w(r48b) &
         - w(r56f) &
         + w(r56b) &
         + w(r57f) &
         - w(r57b) &
         - w(r66f) &
         + w(r66b) &
         + w(r69f) &
         - w(r69b) &
         + w(r83f) &
         - w(r83b) &
         + w(r85f) &
         - w(r85b) &
         - w(r90f) &
         + w(r90b) &
         + w(r95f) &
         - w(r95b) &
         + w(r96f) &
         - w(r96b)

    cdot(sCH3OH) = 0.0_pr &
         - w(r34f) &
         + w(r34b) &
         + w(r50f) &
         - w(r50b) &
         - w(r58f) &
         + w(r58b) &
         - w(r91f) &
         + w(r91b)

    cdot(sC2H2) = 0.0_pr &
         - w(r9f) &
         + w(r9b) &
         - w(r35f) &
         + w(r35b) &
         + w(r36f) &
         - w(r36b)

    cdot(sC2H4) = 0.0_pr &
         - w(r10f) &
         + w(r10b) &
         - w(r37f) &
         + w(r37b) &
         - w(r38f) &
         + w(r38b) &
         - w(r59f) &
         + w(r59b) &
         + w(r71f) &
         - w(r71b) &
         + w(r75f) &
         - w(r75b) &
         + w(r97f) &
         - w(r97b)

    cdot(sC2H6) = 0.0_pr &
         - w(r12f) &
         + w(r12b) &
         + w(r39f) &
         - w(r39b) &
         - w(r40f) &
         + w(r40b) &
         - w(r60f) &
         + w(r60b) &
         + w(r87f) &
         - w(r87b)

    cdot(sNO) = 0.0_pr &
         - w(r98f) &
         + w(r98b) &
         + w(r99f) &
         - w(r99b) &
         + w(r100f) &
         - w(r100b) &
         + 2.0_pr * w(r102f) &
         - 2.0_pr * w(r102b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         + w(r107b) &
         + w(r108f) &
         - w(r108b) &
         + w(r109f) &
         - w(r109b) &
         + w(r110f) &
         - w(r110b) &
         + w(r115f) &
         - w(r115b) &
         - w(r117f) &
         + w(r117b) &
         + w(r123f) &
         - w(r123b) &
         - w(r125f) &
         + w(r125b) &
         + w(r126f) &
         - w(r126b) &
         + w(r127f) &
         - w(r127b) &
         + w(r128f) &
         - w(r128b) &
         + w(r129f) &
         - w(r129b) &
         - w(r134f) &
         + w(r134b) &
         - w(r135f) &
         + w(r135b) &
         - w(r136f) &
         + w(r136b) &
         - w(r137f) &
         + w(r137b) &
         - w(r138f) &
         + w(r138b)

    cdot(sNO2) = 0.0_pr &
         + w(r106f) &
         - w(r106b) &
         + w(r107f) &
         - w(r107b) &
         - w(r108f) &
         + w(r108b) &
         - w(r109f) &
         + w(r109b)

    cdot(sN2O) = 0.0_pr &
         - w(r101f) &
         + w(r101b) &
         - w(r102f) &
         + w(r102b) &
         - w(r103f) &
         + w(r103b) &
         - w(r104f) &
         + w(r104b) &
         - w(r105f) &
         + w(r105b) &
         + w(r117f) &
         - w(r117b)

    cdot(sHCN) = 0.0_pr &
         - w(r130f) &
         + w(r130b) &
         - w(r131f) &
         + w(r131b) &
         - w(r132f) &
         + w(r132b) &
         + w(r133f) &
         - w(r133b) &
         + w(r134f) &
         - w(r134b) &
         + w(r136f) &
         - w(r136b) &
         + w(r138f) &
         - w(r138b) &
         + w(r143f) &
         - w(r143b)

    cdot(sN2) = 0.0_pr &
         + w(r98f) &
         - w(r98b) &
         + w(r101f) &
         - w(r101b) &
         + w(r103f) &
         - w(r103b) &
         + w(r104f) &
         - w(r104b) &
         + w(r105f) &
         - w(r105b) &
         + w(r122f) &
         - w(r122b) &
         + w(r124f) &
         - w(r124b) &
         - w(r133f) &
         + w(r133b)


    return
  end subroutine get_production_rates


  ! ----------------------------------------------- !
  ! Evaluation of QSS concentrations                !
  ! ----------------------------------------------- !
  subroutine get_qss(cqss,c,k,M)
    implicit none

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

    integer :: index

    real(pr) :: CH_ct
    real(pr) :: CH_num
    real(pr) :: CH_denom
    real(pr) :: CH_CH2
    real(pr) :: CH_CH2GSG
    real(pr) :: CH_HCO
    real(pr) :: CH_C2H3
    real(pr) :: CH_N
    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_HNCO
    real(pr) :: CH2GSG_ct
    real(pr) :: CH2GSG_num
    real(pr) :: CH2GSG_denom
    real(pr) :: CH2GSG_CH
    real(pr) :: CH2GSG_CH2
    real(pr) :: CH2GSG_CH3O
    real(pr) :: CH2GSG_HNCO
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_CH
    real(pr) :: HCO_CH2
    real(pr) :: HCO_C2H3
    real(pr) :: CH3O_ct
    real(pr) :: CH3O_num
    real(pr) :: CH3O_denom
    real(pr) :: CH3O_CH2GSG
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_CH
    real(pr) :: C2H3_HCO
    real(pr) :: C2H5_ct
    real(pr) :: C2H5_num
    real(pr) :: C2H5_denom
    real(pr) :: N_ct
    real(pr) :: N_num
    real(pr) :: N_denom
    real(pr) :: N_CH
    real(pr) :: N_NH
    real(pr) :: NH_ct
    real(pr) :: NH_num
    real(pr) :: NH_denom
    real(pr) :: NH_N
    real(pr) :: NH_NH2
    real(pr) :: NH_NNH
    real(pr) :: NH_HNO
    real(pr) :: NH_HNCO
    real(pr) :: NH2_ct
    real(pr) :: NH2_num
    real(pr) :: NH2_denom
    real(pr) :: NH2_NH
    real(pr) :: NH2_HNO
    real(pr) :: NH2_HNCO
    real(pr) :: NNH_ct
    real(pr) :: NNH_num
    real(pr) :: NNH_denom
    real(pr) :: NNH_NH
    real(pr) :: HNO_ct
    real(pr) :: HNO_num
    real(pr) :: HNO_denom
    real(pr) :: HNO_NH
    real(pr) :: HNO_NH2
    real(pr) :: HOCN_ct
    real(pr) :: HOCN_num
    real(pr) :: HOCN_denom
    real(pr) :: HOCN_HNCO
    real(pr) :: HNCO_ct
    real(pr) :: HNCO_num
    real(pr) :: HNCO_denom
    real(pr) :: HNCO_CH2
    real(pr) :: HNCO_CH2GSG
    real(pr) :: HNCO_NH
    real(pr) :: HNCO_NH2
    real(pr) :: HNCO_HOCN

    real(pr) :: a_al_a
    real(pr) :: a_by_a
    real(pr) :: a_ct_a
    real(pr) :: a_cy_a
    real(pr) :: a_dg_a
    real(pr) :: a_dl_a
    real(pr) :: a_fg_a
    real(pr) :: a_fl_a
    real(pr) :: a_ad_b
    real(pr) :: a_aq_b
    real(pr) :: a_bd_b
    real(pr) :: a_bj_b
    real(pr) :: a_bk_b
    real(pr) :: a_bx_b
    real(pr) :: a_dd_b
    real(pr) :: a_dx_b
    real(pr) :: a_ed_b
    real(pr) :: a_ej_b
    real(pr) :: a_ek_b
    real(pr) :: a_eq_b
    real(pr) :: a_ew_b
    real(pr) :: a_ex_b
    real(pr) :: a_aq_c
    real(pr) :: a_bd_c
    real(pr) :: a_bf_c
    real(pr) :: a_bh_c
    real(pr) :: a_bl_c
    real(pr) :: a_by_c
    real(pr) :: a_cd_c
    real(pr) :: a_cf_c
    real(pr) :: a_ch_c
    real(pr) :: a_cl_c
    real(pr) :: a_cu_c
    real(pr) :: a_cy_c
    real(pr) :: a_dd_c
    real(pr) :: a_df_c
    real(pr) :: a_dh_c
    real(pr) :: a_dl_c
    real(pr) :: a_ed_c
    real(pr) :: a_eq_c
    real(pr) :: a_fd_c
    real(pr) :: a_ff_c
    real(pr) :: a_fh_c
    real(pr) :: a_fl_c
    real(pr) :: a_bd_d
    real(pr) :: a_cd_d
    real(pr) :: a_dd_d
    real(pr) :: a_ed_d
    real(pr) :: a_eq_d
    real(pr) :: a_fd_d
    real(pr) :: a_bx_e
    real(pr) :: a_ch_e
    real(pr) :: a_ck_e
    real(pr) :: a_cu_e
    real(pr) :: a_dh_e
    real(pr) :: a_dk_e
    real(pr) :: a_dx_e
    real(pr) :: a_eh_e
    real(pr) :: a_ek_e
    real(pr) :: a_eu_e
    real(pr) :: a_ex_e
    real(pr) :: a_fh_e
    real(pr) :: a_fk_e
    real(pr) :: a_ck_f
    real(pr) :: a_cl_f
    real(pr) :: a_cy_f
    real(pr) :: a_dk_f
    real(pr) :: a_dl_f
    real(pr) :: a_dx_f
    real(pr) :: a_ek_f
    real(pr) :: a_ex_f
    real(pr) :: a_ey_f
    real(pr) :: a_fk_f
    real(pr) :: a_fl_f
    real(pr) :: a_dh_h
    real(pr) :: a_dl_h
    real(pr) :: a_eh_h
    real(pr) :: a_eu_h
    real(pr) :: a_ey_h
    real(pr) :: a_fh_h
    real(pr) :: a_fl_h
    real(pr) :: a_el_i
    real(pr) :: a_ey_i
    real(pr) :: a_fl_i
    real(pr) :: a_ek_j
    real(pr) :: a_ex_j
    real(pr) :: a_fk_j
    real(pr) :: a_ex_k
    real(pr) :: a_fk_k
    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_I_A
    real(pr) :: A_J_A
    real(pr) :: A_K_A
    real(pr) :: A_L_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_I_B
    real(pr) :: A_J_B
    real(pr) :: A_K_B
    real(pr) :: A_L_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_I_C
    real(pr) :: A_J_C
    real(pr) :: A_K_C
    real(pr) :: A_L_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_I_D
    real(pr) :: A_J_D
    real(pr) :: A_K_D
    real(pr) :: A_L_D
    real(pr) :: A_E_E
    real(pr) :: A_F_E
    real(pr) :: A_G_E
    real(pr) :: A_H_E
    real(pr) :: A_I_E
    real(pr) :: A_J_E
    real(pr) :: A_K_E
    real(pr) :: A_L_E
    real(pr) :: A_F_F
    real(pr) :: A_G_F
    real(pr) :: A_H_F
    real(pr) :: A_I_F
    real(pr) :: A_J_F
    real(pr) :: A_K_F
    real(pr) :: A_L_F
    real(pr) :: A_G_G
    real(pr) :: A_H_G
    real(pr) :: A_I_G
    real(pr) :: A_J_G
    real(pr) :: A_K_G
    real(pr) :: A_L_G
    real(pr) :: A_H_H
    real(pr) :: A_I_H
    real(pr) :: A_J_H
    real(pr) :: A_K_H
    real(pr) :: A_L_H
    real(pr) :: A_I_I
    real(pr) :: A_J_I
    real(pr) :: A_K_I
    real(pr) :: A_L_I
    real(pr) :: A_J_J
    real(pr) :: A_K_J
    real(pr) :: A_L_J
    real(pr) :: A_K_K
    real(pr) :: A_L_K
    real(pr) :: A_L_L

    CH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r3f)* c(sO) &
             + k(r47f)* c(sOH) &
             + k(r67f)* c(sO2) &
             + k(r68f)* c(sH2) &
             + k(r69f)* c(sH2O) &
             + k(r70f)* c(sCH3) &
             + k(r71f)* c(sCH4) &
             + k(r72f)* c(sCO2) &
             + k(r133f)* c(sN2) &
             + k(r134f)* c(sNO) &
             + k(r25b)* c(sH2) &
             + k(r49b)* c(sH2O)   )

    CH_num = ( 0.0_pr &
             + k(r3b)* c(sCO) * c(sH) &
             + k(r69b)* c(sCH2O) * c(sH) &
             + k(r71b)* c(sC2H4) * c(sH) &
             + k(r134b)* c(sHCN) * c(sO)   )

    CH_ct = CH_num / CH_denom

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

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

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

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

    CH_N = - ( 0.0_pr &
          + k(r133b) * c(sHCN)  ) / CH_denom

    NH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r110f)* c(sO) &
             + k(r111f)* c(sH) &
             + k(r112f)* c(sOH) &
             + k(r113f)* c(sOH) &
             + k(r114f)* c(sO2) &
             + k(r115f)* c(sO2) &
             + k(r116f)* c(sH2O) &
             + k(r117f)* c(sNO) &
             + k(r118b)* c(sOH) &
             + k(r120b)* c(sH2) &
             + k(r121b)* c(sH2O) &
             + k(r123b)* c(sNO) &
             + k(r130b)* c(sCO) &
             + k(r139b)* c(sCO2)   )

    NH_num = ( 0.0_pr &
             + k(r130f)* c(sHCN) * c(sO) &
             + k(r110b)* c(sH) * c(sNO) &
             + k(r115b)* c(sNO) * c(sOH) &
             + k(r117b)* c(sH) * c(sN2O)   )

    NH_ct = NH_num / NH_denom

    NH_N = - ( 0.0_pr &
          + k(r111b) * c(sH2) &
          + k(r113b) * c(sH2O)  ) / NH_denom

    NH_NH2 = - ( 0.0_pr &
          + k(r118f) * c(sO) &
          + k(r120f) * c(sH) &
          + k(r121f) * c(sOH)  ) / NH_denom

    NH_NNH = - ( 0.0_pr &
          + k(r123f) * c(sO)  ) / NH_denom

    NH_HNO = - ( 0.0_pr &
          + k(r112b) * c(sH) &
          + k(r114b) * c(sO) &
          + k(r116b) * c(sH2)  ) / NH_denom

    NH_HNCO = - ( 0.0_pr &
          + k(r139f) * c(sO)  ) / NH_denom

    HNO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r126f)* c(sO) &
             + k(r127f)* c(sH) &
             + k(r128f)* c(sOH) &
             + k(r129f)* c(sO2) &
             + k(r112b)* c(sH) &
             + k(r114b)* c(sO) &
             + k(r116b)* c(sH2) &
             + k(r119b)* c(sH) &
             + k(r125b)* m(mM125)   )

    HNO_num = ( 0.0_pr &
             + k(r125f)* c(sH) * c(sNO) * m(mM125) &
             + k(r126b)* c(sNO) * c(sOH) &
             + k(r127b)* c(sH2) * c(sNO) &
             + k(r128b)* c(sH2O) * c(sNO) &
             + k(r129b)* c(sHO2) * c(sNO)   )

    HNO_ct = HNO_num / HNO_denom

    HNO_NH = - ( 0.0_pr &
          + k(r112f) * c(sOH) &
          + k(r114f) * c(sO2) &
          + k(r116f) * c(sH2O)  ) / HNO_denom

    HNO_NH2 = - ( 0.0_pr &
          + k(r119f) * c(sO)  ) / HNO_denom

    NNH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r122f)* m(mM122) &
             + k(r123f)* c(sO) &
             + k(r124f)* c(sH)   )

    NNH_num = ( 0.0_pr &
             + k(r122b)* c(sH) * c(sN2) * m(mM122) &
             + k(r124b)* c(sH2) * c(sN2)   )

    NNH_ct = NNH_num / NNH_denom

    NNH_NH = - ( 0.0_pr &
          + k(r123b) * c(sNO)  ) / NNH_denom

    CH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r4f)* c(sO) &
             + k(r24f)* c(sH) &
             + k(r48f)* c(sOH) &
             + k(r49f)* c(sOH) &
             + k(r73f)* c(sO2) &
             + k(r74f)* c(sH2) &
             + k(r75f)* c(sCH3) &
             + k(r135f)* c(sNO) &
             + k(r136f)* c(sNO) &
             + k(r9b)* c(sCO) &
             + k(r51b)* c(sH2O) &
             + k(r68b)* c(sH) &
             + k(r76b)* c(sN2) &
             + k(r80b)* c(sH2O) &
             + k(r82b)* c(sCO)   )

    CH2_num = ( 0.0_pr &
             + k(r9f)* c(sC2H2) * c(sO) &
             + k(r51f)* c(sCH3) * c(sOH) &
             + k(r24b)* c(sCH3) &
             + k(r48b)* c(sCH2O) * c(sH) &
             + k(r74b)* c(sCH3) * c(sH) &
             + k(r75b)* c(sC2H4) * c(sH) &
             + k(r136b)* c(sHCN) * c(sOH)   )

    CH2_ct = CH2_num / CH2_denom

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

    CH2_CH2GSG = - ( 0.0_pr &
          + k(r76f) * c(sN2) &
          + k(r80f) * c(sH2O) &
          + k(r82f) * c(sCO)  ) / CH2_denom

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

    CH2_HNCO = - ( 0.0_pr &
          + k(r135b) * c(sH)  ) / CH2_denom

    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r7f)* c(sO) &
             + k(r28f)* c(sH) &
             + k(r55f)* c(sOH) &
             + k(r89f)* c(sCH3) &
             + k(r92f)* c(sH2O) &
             + k(r93f)* m(mM93) &
             + k(r94f)* c(sO2) &
             + k(r4b)* c(sH) &
             + k(r8b)* c(sOH) &
             + k(r10b)* c(sCH3) &
             + k(r13b)* c(sHO2) &
             + k(r30b)* c(sH2) &
             + k(r47b)* c(sH) &
             + k(r56b)* c(sH2O) &
             + k(r66b)* c(sH2O2) &
             + k(r67b)* c(sO) &
             + k(r72b)* c(sCO) &
             + k(r73b)* c(sOH) &
             + k(r90b)* c(sCH4) &
             + k(r96b)* c(sCH2O)   )

    HCO_num = ( 0.0_pr &
             + k(r8f)* c(sCH2O) * c(sO) &
             + k(r10f)* c(sC2H4) * c(sO) &
             + k(r13f)* c(sCH2O) * c(sO2) &
             + k(r30f)* c(sCH2O) * c(sH) &
             + k(r56f)* c(sCH2O) * c(sOH) &
             + k(r66f)* c(sCH2O) * c(sHO2) &
             + k(r90f)* c(sCH2O) * c(sCH3) &
             + k(r7b)* c(sCO2) * c(sH) &
             + k(r28b)* c(sCO) * c(sH2) &
             + k(r55b)* c(sCO) * c(sH2O) &
             + k(r89b)* c(sCH4) * c(sCO) &
             + k(r92b)* c(sCO) * c(sH) * c(sH2O) &
             + k(r93b)* c(sCO) * c(sH) * m(mM93) &
             + k(r94b)* c(sCO) * c(sHO2)   )

    HCO_ct = HCO_num / HCO_denom

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

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

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

    CH3O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r31f)* c(sH) &
             + k(r32f)* c(sH) &
             + k(r33f)* c(sH) &
             + k(r57f)* c(sOH) &
             + k(r95f)* c(sO2) &
             + k(r29b)&
             + k(r34b)* c(sH2) &
             + k(r58b)* c(sH2O) &
             + k(r64b)* c(sOH) &
             + k(r84b)* c(sO) &
             + k(r91b)* c(sCH4)   )

    CH3O_num = ( 0.0_pr &
             + k(r29f)* c(sCH2O) * c(sH) &
             + k(r34f)* c(sCH3OH) * c(sH) &
             + k(r58f)* c(sCH3OH) * c(sOH) &
             + k(r64f)* c(sCH3) * c(sHO2) &
             + k(r84f)* c(sCH3) * c(sO2) &
             + k(r91f)* c(sCH3) * c(sCH3OH) &
             + k(r31b)* c(sCH2O) * c(sH2) &
             + k(r32b)* c(sCH3) * c(sOH) &
             + k(r57b)* c(sCH2O) * c(sH2O) &
             + k(r95b)* c(sCH2O) * c(sHO2)   )

    CH3O_ct = CH3O_num / CH3O_denom

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

    N_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r98f)* c(sNO) &
             + k(r99f)* c(sO2) &
             + k(r100f)* c(sOH) &
             + k(r143f)* c(sCH3) &
             + k(r111b)* c(sH2) &
             + k(r113b)* c(sH2O) &
             + k(r133b)* c(sHCN)   )

    N_num = ( 0.0_pr &
             + k(r98b)* c(sN2) * c(sO) &
             + k(r99b)* c(sNO) * c(sO) &
             + k(r100b)* c(sH) * c(sNO) &
             + k(r143b)* c(sH2) * c(sHCN)   )

    N_ct = N_num / N_denom

    N_CH = - ( 0.0_pr &
          + k(r133f) * c(sN2)  ) / N_denom

    N_NH = - ( 0.0_pr &
          + k(r111f) * c(sH) &
          + k(r113f) * c(sOH)  ) / N_denom

    HNCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r139f)* c(sO) &
             + k(r140f)* c(sH) &
             + k(r141f)* c(sOH) &
             + k(r132b)* c(sH) &
             + k(r135b)* c(sH) &
             + k(r137b)* c(sH) &
             + k(r142b)* c(sH)   )

    HNCO_num = ( 0.0_pr &
             + k(r132f)* c(sHCN) * c(sOH)   )

    HNCO_ct = HNCO_num / HNCO_denom

    HNCO_CH2 = - ( 0.0_pr &
          + k(r135f) * c(sNO)  ) / HNCO_denom

    HNCO_CH2GSG = - ( 0.0_pr &
          + k(r137f) * c(sNO)  ) / HNCO_denom

    HNCO_NH = - ( 0.0_pr &
          + k(r139b) * c(sCO2)  ) / HNCO_denom

    HNCO_NH2 = - ( 0.0_pr &
          + k(r140b) * c(sCO) &
          + k(r141b) * c(sCO2)  ) / HNCO_denom

    HNCO_HOCN = - ( 0.0_pr &
          + k(r142f) * c(sH)  ) / HNCO_denom

    HOCN_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r142f)* c(sH) &
             + k(r131b)* c(sH)   )

    HOCN_num = ( 0.0_pr &
             + k(r131f)* c(sHCN) * c(sOH)   )

    HOCN_ct = HOCN_num / HOCN_denom

    HOCN_HNCO = - ( 0.0_pr &
          + k(r142b) * c(sH)  ) / HOCN_denom

    CH2GSG_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r25f)* c(sH) &
             + k(r76f)* c(sN2) &
             + k(r77f)* c(sO2) &
             + k(r78f)* c(sO2) &
             + k(r79f)* c(sH2) &
             + k(r80f)* c(sH2O) &
             + k(r81f)* c(sCH4) &
             + k(r82f)* c(sCO) &
             + k(r83f)* c(sCO2) &
             + k(r137f)* c(sNO) &
             + k(r33b)* c(sH2O) &
             + k(r52b)* c(sH2O)   )

    CH2GSG_num = ( 0.0_pr &
             + k(r52f)* c(sCH3) * c(sOH) &
             + k(r77b)* c(sCO) * c(sH) * c(sOH) &
             + k(r78b)* c(sCO) * c(sH2O) &
             + k(r79b)* c(sCH3) * c(sH) &
             + k(r81b) *c(sCH3)** 2.0_pr&
             + k(r83b)* c(sCH2O) * c(sCO)   )

    CH2GSG_ct = CH2GSG_num / CH2GSG_denom

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

    CH2GSG_CH2 = - ( 0.0_pr &
          + k(r76b) * c(sN2) &
          + k(r80b) * c(sH2O) &
          + k(r82b) * c(sCO)  ) / CH2GSG_denom

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

    CH2GSG_HNCO = - ( 0.0_pr &
          + k(r137b) * c(sH)  ) / CH2GSG_denom

    NH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r118f)* c(sO) &
             + k(r119f)* c(sO) &
             + k(r120f)* c(sH) &
             + k(r121f)* c(sOH) &
             + k(r140b)* c(sCO) &
             + k(r141b)* c(sCO2)   )

    NH2_num = ( 0.0_pr   )

    NH2_ct = NH2_num / NH2_denom

    NH2_NH = - ( 0.0_pr &
          + k(r118b) * c(sOH) &
          + k(r120b) * c(sH2) &
          + k(r121b) * c(sH2O)  ) / NH2_denom

    NH2_HNO = - ( 0.0_pr &
          + k(r119b) * c(sH)  ) / NH2_denom

    NH2_HNCO = - ( 0.0_pr &
          + k(r140f) * c(sH) &
          + k(r141f) * c(sOH)  ) / NH2_denom

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r36f)* c(sH) &
             + k(r96f)* c(sO2) &
             + k(r35b)&
             + k(r38b)* c(sH2) &
             + k(r59b)* c(sH2O) &
             + k(r70b)* c(sH)   )

    C2H3_num = ( 0.0_pr &
             + k(r35f)* c(sC2H2) * c(sH) &
             + k(r38f)* c(sC2H4) * c(sH) &
             + k(r59f)* c(sC2H4) * c(sOH) &
             + k(r36b)* c(sC2H2) * c(sH2)   )

    C2H3_ct = C2H3_num / C2H3_denom

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

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

    a_al_a = CH2GSG_CH
    a_by_a = N_CH
    a_ct_a =  (1.0_pr) - (C2H3_HCO) * (HCO_C2H3)
    a_cy_a =  (HCO_CH) - (C2H3_CH) * (HCO_C2H3)
    a_dg_a = CH2_HCO
    a_dl_a = CH2_CH
    a_fg_a =  (CH_HCO) - (C2H3_HCO) * (CH_C2H3)
    a_fl_a =  (1.0_pr) - (C2H3_CH) * (CH_C2H3)
    a_ad_b = CH2GSG_HNCO
    a_aq_b = HOCN_HNCO
    a_bd_b =  (1.0_pr) - (NH2_HNCO) * (HNCO_NH2)
    a_bj_b =  - (NH2_HNO) * (HNCO_NH2)
    a_bk_b =  (HNCO_NH) - (NH2_NH) * (HNCO_NH2)
    a_bx_b = N_NH
    a_dd_b = CH2_HNCO
    a_dx_b = NNH_NH
    a_ed_b =  - (NH2_HNCO) * (HNO_NH2)
    a_ej_b =  (1.0_pr) - (NH2_HNO) * (HNO_NH2)
    a_ek_b =  (HNO_NH) - (NH2_NH) * (HNO_NH2)
    a_eq_b =  (NH_HNCO) - (NH2_HNCO) * (NH_NH2)
    a_ew_b =  (NH_HNO) - (NH2_HNO) * (NH_NH2)
    a_ex_b =  (1.0_pr) - (NH2_NH) * (NH_NH2)
    a_aq_c = a_aq_b
    a_bd_c =  (a_bd_b) - (a_ad_b) * (HNCO_CH2GSG)
    a_bf_c =  - (CH2GSG_CH3O) * (HNCO_CH2GSG)
    a_bh_c =  (HNCO_CH2) - (CH2GSG_CH2) * (HNCO_CH2GSG)
    a_bl_c =  - (a_al_a) * (HNCO_CH2GSG)
    a_by_c = a_by_a
    a_cd_c =  - (a_ad_b) * (CH3O_CH2GSG)
    a_cf_c =  (1.0_pr) - (CH2GSG_CH3O) * (CH3O_CH2GSG)
    a_ch_c =  - (CH2GSG_CH2) * (CH3O_CH2GSG)
    a_cl_c =  - (a_al_a) * (CH3O_CH2GSG)
    a_cu_c = HCO_CH2
    a_cy_c = a_cy_a
    a_dd_c =  (a_dd_b) - (a_ad_b) * (CH2_CH2GSG)
    a_df_c =  - (CH2GSG_CH3O) * (CH2_CH2GSG)
    a_dh_c =  (1.0_pr) - (CH2GSG_CH2) * (CH2_CH2GSG)
    a_dl_c =  (a_dl_a) - (a_al_a) * (CH2_CH2GSG)
    a_ed_c = a_ed_b
    a_eq_c = a_eq_b
    a_fd_c =  - (a_ad_b) * (CH_CH2GSG)
    a_ff_c =  - (CH2GSG_CH3O) * (CH_CH2GSG)
    a_fh_c =  (CH_CH2) - (CH2GSG_CH2) * (CH_CH2GSG)
    a_fl_c =  (a_fl_a) - (a_al_a) * (CH_CH2GSG)
    a_bd_d =  (a_bd_c) - (a_aq_c) * (HNCO_HOCN)
    a_cd_d = a_cd_c
    a_dd_d = a_dd_c
    a_ed_d = a_ed_c
    a_eq_d = a_eq_c
    a_fd_d = a_fd_c
    a_bx_e = a_bx_b
    a_ch_e =  (a_ch_c) - (a_bh_c) * (a_cd_d) / (a_bd_d)
    a_ck_e =  - (a_bk_b) * (a_cd_d) / (a_bd_d)
    a_cu_e = a_cu_c
    a_dh_e =  (a_dh_c) - (a_bh_c) * (a_dd_d) / (a_bd_d)
    a_dk_e =  - (a_bk_b) * (a_dd_d) / (a_bd_d)
    a_dx_e = a_dx_b
    a_eh_e =  - (a_bh_c) * (a_ed_d) / (a_bd_d)
    a_ek_e =  (a_ek_b) - (a_bk_b) * (a_ed_d) / (a_bd_d)
    a_eu_e =  - (a_bh_c) * (a_eq_d) / (a_bd_d)
    a_ex_e =  (a_ex_b) - (a_bk_b) * (a_eq_d) / (a_bd_d)
    a_fh_e =  (a_fh_c) - (a_bh_c) * (a_fd_d) / (a_bd_d)
    a_fk_e =  - (a_bk_b) * (a_fd_d) / (a_bd_d)
    a_ck_f = a_ck_e
    a_cl_f = a_cl_c
    a_cy_f = a_cy_c
    a_dk_f = a_dk_e
    a_dl_f = a_dl_c
    a_dx_f = a_dx_e
    a_ek_f = a_ek_e
    a_ex_f =  (a_ex_e) - (a_bx_e) * (NH_N)
    a_ey_f =  - (a_by_c) * (NH_N)
    a_fk_f =  (a_fk_e) - (a_bx_e) * (CH_N)
    a_fl_f =  (a_fl_c) - (a_by_c) * (CH_N)
    a_dh_h =  (a_dh_e) - (a_cu_e) * (a_dg_a) / (a_ct_a)
    a_dl_h =  (a_dl_f) - (a_cy_f) * (a_dg_a) / (a_ct_a)
    a_eh_h = a_eh_e
    a_eu_h = a_eu_e
    a_ey_h = a_ey_f
    a_fh_h =  (a_fh_e) - (a_cu_e) * (a_fg_a) / (a_ct_a)
    a_fl_h =  (a_fl_f) - (a_cy_f) * (a_fg_a) / (a_ct_a)
    a_el_i =  - (a_dl_h) * (a_eh_h) / (a_dh_h)
    a_ey_i =  (a_ey_h) - (a_dl_h) * (a_eu_h) / (a_dh_h)
    a_fl_i =  (a_fl_h) - (a_dl_h) * (a_fh_h) / (a_dh_h)
    a_ek_j = a_ek_f
    a_ex_j =  (a_ex_f) - (a_dx_f) * (NH_NNH)
    a_fk_j = a_fk_f
    a_ex_k =  (a_ex_j) - (a_ek_j) * (a_ew_b) / (a_ej_b)
    a_fk_k = a_fk_j
    A_A_A = NH2_ct
    A_B_A = CH2GSG_ct
    A_C_A = HOCN_ct
    A_D_A = HNCO_ct
    A_E_A = N_ct
    A_F_A = CH3O_ct
    A_G_A =  (HCO_ct) - (C2H3_ct) * (HCO_C2H3)
    A_H_A = CH2_ct
    A_I_A = NNH_ct
    A_J_A = HNO_ct
    A_K_A = NH_ct
    A_L_A =  (CH_ct) - (C2H3_ct) * (CH_C2H3)
    A_B_B = A_B_A
    A_C_B = A_C_A
    A_D_B =  (A_D_A) - (A_A_A) * (HNCO_NH2)
    A_E_B = A_E_A
    A_F_B = A_F_A
    A_G_B = A_G_A
    A_H_B = A_H_A
    A_I_B = A_I_A
    A_J_B =  (A_J_A) - (A_A_A) * (HNO_NH2)
    A_K_B =  (A_K_A) - (A_A_A) * (NH_NH2)
    A_L_B = A_L_A
    A_C_C = A_C_B
    A_D_C =  (A_D_B) - (A_B_B) * (HNCO_CH2GSG)
    A_E_C = A_E_B
    A_F_C =  (A_F_B) - (A_B_B) * (CH3O_CH2GSG)
    A_G_C = A_G_B
    A_H_C =  (A_H_B) - (A_B_B) * (CH2_CH2GSG)
    A_I_C = A_I_B
    A_J_C = A_J_B
    A_K_C = A_K_B
    A_L_C =  (A_L_B) - (A_B_B) * (CH_CH2GSG)
    A_D_D =  (A_D_C) - (A_C_C) * (HNCO_HOCN)
    A_E_D = A_E_C
    A_F_D = A_F_C
    A_G_D = A_G_C
    A_H_D = A_H_C
    A_I_D = A_I_C
    A_J_D = A_J_C
    A_K_D = A_K_C
    A_L_D = A_L_C
    A_E_E = A_E_D
    A_F_E =  (A_F_D) - (A_D_D) * (a_cd_d) / (a_bd_d)
    A_G_E = A_G_D
    A_H_E =  (A_H_D) - (A_D_D) * (a_dd_d) / (a_bd_d)
    A_I_E = A_I_D
    A_J_E =  (A_J_D) - (A_D_D) * (a_ed_d) / (a_bd_d)
    A_K_E =  (A_K_D) - (A_D_D) * (a_eq_d) / (a_bd_d)
    A_L_E =  (A_L_D) - (A_D_D) * (a_fd_d) / (a_bd_d)
    A_F_F = A_F_E
    A_G_F = A_G_E
    A_H_F = A_H_E
    A_I_F = A_I_E
    A_J_F = A_J_E
    A_K_F =  (A_K_E) - (A_E_E) * (NH_N)
    A_L_F =  (A_L_E) - (A_E_E) * (CH_N)
    A_G_G = A_G_F
    A_H_G =  (A_H_F) - (A_F_F) * (a_df_c) / (a_cf_c)
    A_I_G = A_I_F
    A_J_G = A_J_F
    A_K_G = A_K_F
    A_L_G =  (A_L_F) - (A_F_F) * (a_ff_c) / (a_cf_c)
    A_H_H =  (A_H_G) - (A_G_G) * (a_dg_a) / (a_ct_a)
    A_I_H = A_I_G
    A_J_H = A_J_G
    A_K_H = A_K_G
    A_L_H =  (A_L_G) - (A_G_G) * (a_fg_a) / (a_ct_a)
    A_I_I = A_I_H
    A_J_I =  (A_J_H) - (A_H_H) * (a_eh_h) / (a_dh_h)
    A_K_I =  (A_K_H) - (A_H_H) * (a_eu_h) / (a_dh_h)
    A_L_I =  (A_L_H) - (A_H_H) * (a_fh_h) / (a_dh_h)
    A_J_J = A_J_I
    A_K_J =  (A_K_I) - (A_I_I) * (NH_NNH)
    A_L_J = A_L_I
    A_K_K =  (A_K_J) - (A_J_J) * (a_ew_b) / (a_ej_b)
    A_L_K = A_L_J
    A_L_L =  (A_L_K) - (A_K_K) * (a_fk_k) / (a_ex_k)

    cqss(sqssCH) = ( A_L_L ) / ( a_fl_i )

    cqss(sqssNH) = (A_K_K - (a_ey_i) * cqss(sqssCH)) / (a_ex_k)

    cqss(sqssHNO) = (A_J_J - (a_el_i) * cqss(sqssCH) - (a_ek_j) &
                 * cqss(sqssNH)) / (a_ej_b)

    cqss(sqssNNH) = A_I_I - (a_dx_f) * cqss(sqssNH)

    cqss(sqssCH2) = (A_H_H - (a_dl_h) * cqss(sqssCH) - (a_dk_f) &
                 * cqss(sqssNH)) / (a_dh_h)

    cqss(sqssHCO) = (A_G_G - (a_cy_f) * cqss(sqssCH) - (a_cu_e) &
                 * cqss(sqssCH2)) / (a_ct_a)

    cqss(sqssCH3O) = (A_F_F - (a_cl_f) * cqss(sqssCH) - (a_ck_f) &
                 * cqss(sqssNH) - (a_ch_e) * cqss(sqssCH2)) &
                 / (a_cf_c)

    cqss(sqssN) = A_E_E - (a_by_c) * cqss(sqssCH) - (a_bx_e) &
                 * cqss(sqssNH)

    cqss(sqssHNCO) = (A_D_D - (a_bl_c) * cqss(sqssCH) - (a_bk_b) &
                 * cqss(sqssNH) - (a_bj_b) * cqss(sqssHNO) - (a_bh_c) &
                 * cqss(sqssCH2) - (a_bf_c) * cqss(sqssCH3O)) &
                 / (a_bd_d)

    cqss(sqssHOCN) = A_C_C - (a_aq_c) * cqss(sqssHNCO)

    cqss(sqssCH2GSG) = A_B_B - (a_al_a) * cqss(sqssCH) - (CH2GSG_CH2) &
                 * cqss(sqssCH2) - (CH2GSG_CH3O) * cqss(sqssCH3O) &
                 - (a_ad_b) * cqss(sqssHNCO)

    cqss(sqssNH2) = A_A_A - (NH2_NH) * cqss(sqssNH) - (NH2_HNO) &
                 * cqss(sqssHNO) - (NH2_HNCO) * cqss(sqssHNCO)

    cqss(sqssC2H3) = C2H3_ct - (C2H3_CH) * cqss(sqssCH) - (C2H3_HCO) &
                 * cqss(sqssHCO)

    C2H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r11f)* c(sO) &
             + k(r39f)* c(sH) &
             + k(r97f)* c(sO2) &
             + k(r12b)* c(sOH) &
             + k(r37b)&
             + k(r40b)* c(sH2) &
             + k(r60b)* c(sH2O) &
             + k(r88b)* c(sH)   )

    C2H5_num = ( 0.0_pr &
             + k(r12f)* c(sC2H6) * c(sO) &
             + k(r37f)* c(sC2H4) * c(sH) &
             + k(r40f)* c(sC2H6) * c(sH) &
             + k(r60f)* c(sC2H6) * c(sOH) &
             + k(r88f) *c(sCH3)** 2.0_pr&
             + k(r11b)* c(sCH2O) * c(sCH3) &
             + k(r39b)* c(sC2H6) &
             + k(r97b)* 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