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

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

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

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

  integer, parameter :: sqssHO2 = 1
  integer, parameter :: sqssCH3OH = 2
  integer, parameter :: sqssCH3O = 3
  integer, parameter :: sqssCH2OH = 4
  integer, parameter :: sqssHCO = 5
  integer, parameter :: sqssC2H5 = 6
  integer, parameter :: sqssC2H4 = 7
  integer, parameter :: sqssC2H3 = 8
  integer, parameter :: sqssCH2CO = 9

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r39f = 39
  integer, parameter :: r40f = 40
  integer, parameter :: r41f = 41
  integer, parameter :: r42f = 42
  integer, parameter :: r43f = 43
  integer, parameter :: r44f = 44
  integer, parameter :: r45f = 45
  integer, parameter :: r46f = 46
  integer, parameter :: r47f = 47
  integer, parameter :: r48f = 48
  integer, parameter :: r49f = 49
  integer, parameter :: r50f = 50
  integer, parameter :: r51f = 51
  integer, parameter :: r52f = 52
  integer, parameter :: r53f = 53
  integer, parameter :: r54f = 54
  integer, parameter :: r55f = 55
  integer, parameter :: r56f = 56
  integer, parameter :: r57f = 57
  integer, parameter :: r58f = 58
  integer, parameter :: r59f = 59
  integer, parameter :: r60f = 60
  integer, parameter :: r61f = 61
  integer, parameter :: r62f = 62
  integer, parameter :: r63f = 63
  integer, parameter :: r64f = 64
  integer, parameter :: r65f = 65
  integer, parameter :: r66f = 66
  integer, parameter :: r67f = 67
  integer, parameter :: r68f = 68
  integer, parameter :: r69f = 69
  integer, parameter :: r70f = 70
  integer, parameter :: r71f = 71
  integer, parameter :: r72f = 72
  integer, parameter :: r73f = 73
  integer, parameter :: r74f = 74
  integer, parameter :: r75f = 75
  integer, parameter :: r76f = 76
  integer, parameter :: r77f = 77
  integer, parameter :: r78f = 78
  integer, parameter :: r79f = 79
  integer, parameter :: r80f = 80
  integer, parameter :: r81f = 81
  integer, parameter :: r82f = 82
  integer, parameter :: r83f = 83
  integer, parameter :: r84f = 84
  integer, parameter :: r85f = 85
  integer, parameter :: r86f = 86
  integer, parameter :: r87f = 87
  integer, parameter :: r88f = 88
  integer, parameter :: r89f = 89
  integer, parameter :: r90f = 90
  integer, parameter :: r91f = 91
  integer, parameter :: r92f = 92
  integer, parameter :: r93f = 93
  integer, parameter :: r94f = 94
  integer, parameter :: r95f = 95
  integer, parameter :: r96f = 96
  integer, parameter :: r97f = 97
  integer, parameter :: r98f = 98
  integer, parameter :: r99f = 99
  integer, parameter :: r100f = 100
  integer, parameter :: r101f = 101
  integer, parameter :: r102f = 102
  integer, parameter :: r103f = 103
  integer, parameter :: r104f = 104
  integer, parameter :: r105f = 105
  integer, parameter :: r106f = 106
  integer, parameter :: r107f = 107
  integer, parameter :: r108f = 108
  integer, parameter :: r109f = 109
  integer, parameter :: r110f = 110
  integer, parameter :: r111f = 111
  integer, parameter :: r112f = 112
  integer, parameter :: r113f = 113
  integer, parameter :: r114f = 114
  integer, parameter :: r115f = 115
  integer, parameter :: r116f = 116
  integer, parameter :: r117f = 117
  integer, parameter :: r118f = 118
  integer, parameter :: r119f = 119
  integer, parameter :: r120f = 120
  integer, parameter :: r121f = 121
  integer, parameter :: r122f = 122
  integer, parameter :: r123f = 123
  integer, parameter :: r124f = 124
  integer, parameter :: r125f = 125
  integer, parameter :: r126f = 126
  integer, parameter :: r127f = 127
  integer, parameter :: r128f = 128
  integer, parameter :: r129f = 129
  integer, parameter :: r130f = 130
  integer, parameter :: r131f = 131
  integer, parameter :: r132f = 132
  integer, parameter :: r133f = 133
  integer, parameter :: r134f = 134
  integer, parameter :: r135f = 135
  integer, parameter :: r136f = 136
  integer, parameter :: r137f = 137
  integer, parameter :: r138f = 138
  integer, parameter :: r1b = 139
  integer, parameter :: r2b = 140
  integer, parameter :: r3b = 141
  integer, parameter :: r4b = 142
  integer, parameter :: r5b = 143
  integer, parameter :: r6b = 144
  integer, parameter :: r7b = 145
  integer, parameter :: r8b = 146
  integer, parameter :: r9b = 147
  integer, parameter :: r10b = 148
  integer, parameter :: r11b = 149
  integer, parameter :: r12b = 150
  integer, parameter :: r13b = 151
  integer, parameter :: r14b = 152
  integer, parameter :: r15b = 153
  integer, parameter :: r16b = 154
  integer, parameter :: r17b = 155
  integer, parameter :: r18b = 156
  integer, parameter :: r19b = 157
  integer, parameter :: r20b = 158
  integer, parameter :: r22b = 159
  integer, parameter :: r23b = 160
  integer, parameter :: r24b = 161
  integer, parameter :: r25b = 162
  integer, parameter :: r26b = 163
  integer, parameter :: r27b = 164
  integer, parameter :: r28b = 165
  integer, parameter :: r29b = 166
  integer, parameter :: r32b = 167
  integer, parameter :: r33b = 168
  integer, parameter :: r34b = 169
  integer, parameter :: r35b = 170
  integer, parameter :: r36b = 171
  integer, parameter :: r37b = 172
  integer, parameter :: r38b = 173
  integer, parameter :: r39b = 174
  integer, parameter :: r40b = 175
  integer, parameter :: r41b = 176
  integer, parameter :: r42b = 177
  integer, parameter :: r43b = 178
  integer, parameter :: r44b = 179
  integer, parameter :: r45b = 180
  integer, parameter :: r46b = 181
  integer, parameter :: r47b = 182
  integer, parameter :: r48b = 183
  integer, parameter :: r49b = 184
  integer, parameter :: r50b = 185
  integer, parameter :: r51b = 186
  integer, parameter :: r52b = 187
  integer, parameter :: r53b = 188
  integer, parameter :: r54b = 189
  integer, parameter :: r55b = 190
  integer, parameter :: r56b = 191
  integer, parameter :: r57b = 192
  integer, parameter :: r58b = 193
  integer, parameter :: r59b = 194
  integer, parameter :: r60b = 195
  integer, parameter :: r61b = 196
  integer, parameter :: r62b = 197
  integer, parameter :: r63b = 198
  integer, parameter :: r64b = 199
  integer, parameter :: r65b = 200
  integer, parameter :: r66b = 201
  integer, parameter :: r67b = 202
  integer, parameter :: r68b = 203
  integer, parameter :: r69b = 204
  integer, parameter :: r70b = 205
  integer, parameter :: r71b = 206
  integer, parameter :: r72b = 207
  integer, parameter :: r73b = 208
  integer, parameter :: r74b = 209
  integer, parameter :: r75b = 210
  integer, parameter :: r76b = 211
  integer, parameter :: r77b = 212
  integer, parameter :: r78b = 213
  integer, parameter :: r79b = 214
  integer, parameter :: r80b = 215
  integer, parameter :: r81b = 216
  integer, parameter :: r82b = 217
  integer, parameter :: r83b = 218
  integer, parameter :: r84b = 219
  integer, parameter :: r85b = 220
  integer, parameter :: r86b = 221
  integer, parameter :: r87b = 222
  integer, parameter :: r88b = 223
  integer, parameter :: r89b = 224
  integer, parameter :: r90b = 225
  integer, parameter :: r91b = 226
  integer, parameter :: r92b = 227
  integer, parameter :: r93b = 228
  integer, parameter :: r94b = 229
  integer, parameter :: r95b = 230
  integer, parameter :: r98b = 231
  integer, parameter :: r99b = 232
  integer, parameter :: r100b = 233
  integer, parameter :: r101b = 234
  integer, parameter :: r102b = 235
  integer, parameter :: r103b = 236
  integer, parameter :: r104b = 237
  integer, parameter :: r105b = 238
  integer, parameter :: r106b = 239
  integer, parameter :: r107b = 240
  integer, parameter :: r108b = 241
  integer, parameter :: r109b = 242
  integer, parameter :: r110b = 243
  integer, parameter :: r111b = 244
  integer, parameter :: r112b = 245
  integer, parameter :: r113b = 246
  integer, parameter :: r117b = 247
  integer, parameter :: r119b = 248
  integer, parameter :: r120b = 249
  integer, parameter :: r121b = 250
  integer, parameter :: r122b = 251
  integer, parameter :: r123b = 252
  integer, parameter :: r124b = 253
  integer, parameter :: r127b = 254
  integer, parameter :: r131b = 255
  integer, parameter :: r138b = 256
  
  ! Index of third body species
  integer, parameter :: mM4 = 1
  integer, parameter :: mM63 = 2
  integer, parameter :: mM82 = 3
  integer, parameter :: mM117 = 4
  integer, parameter :: mM136 = 5
  
  integer, parameter :: mM11 = 6
  integer, parameter :: mM12 = 7
  integer, parameter :: mM15 = 8
  integer, parameter :: mM28 = 9
  integer, parameter :: mM40 = 10
  integer, parameter :: mM41 = 11
  integer, parameter :: mM57 = 12
  integer, parameter :: mM58 = 13
  integer, parameter :: mM70 = 14
  integer, parameter :: mM71 = 15
  integer, parameter :: mM72 = 16
  integer, parameter :: mM73 = 17
  integer, parameter :: mM78 = 18
  integer, parameter :: mM81 = 19
  integer, parameter :: mM89 = 20
  integer, parameter :: mM108 = 21
  
  ! Index of third body reactions
  integer, parameter :: TBr4f = 1
  integer, parameter :: TBr63f = 2
  integer, parameter :: TBr82f = 3
  integer, parameter :: TBr117f = 4
  integer, parameter :: TBr136f = 5
  integer, parameter :: TBr4b = 6
  integer, parameter :: TBr63b = 7
  integer, parameter :: TBr82b = 8
  integer, parameter :: TBr117b = 9
  
  ! Index of fall off reactions
  integer, parameter :: FOr11f = 1
  integer, parameter :: FOr12f = 2
  integer, parameter :: FOr15f = 3
  integer, parameter :: FOr28f = 4
  integer, parameter :: FOr40f = 5
  integer, parameter :: FOr41f = 6
  integer, parameter :: FOr57f = 7
  integer, parameter :: FOr58f = 8
  integer, parameter :: FOr70f = 9
  integer, parameter :: FOr71f = 10
  integer, parameter :: FOr72f = 11
  integer, parameter :: FOr73f = 12
  integer, parameter :: FOr78f = 13
  integer, parameter :: FOr81f = 14
  integer, parameter :: FOr89f = 15
  integer, parameter :: FOr108f = 16
  integer, parameter :: FOr11b = 17
  integer, parameter :: FOr12b = 18
  integer, parameter :: FOr15b = 19
  integer, parameter :: FOr28b = 20
  integer, parameter :: FOr40b = 21
  integer, parameter :: FOr41b = 22
  integer, parameter :: FOr57b = 23
  integer, parameter :: FOr58b = 24
  integer, parameter :: FOr70b = 25
  integer, parameter :: FOr71b = 26
  integer, parameter :: FOr72b = 27
  integer, parameter :: FOr73b = 28
  integer, parameter :: FOr78b = 29
  integer, parameter :: FOr81b = 30
  integer, parameter :: FOr89b = 31
  integer, parameter :: FOr108b = 32
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.02801348_pr, & ! N2
       0.00201588_pr, & ! H2
       0.00100794_pr, & ! H
       0.0319988_pr, & ! O2
       0.0159994_pr, & ! O
       0.01801528_pr, & ! H2O
       0.01700734_pr, & ! OH
       0.0280104_pr, & ! CO
       0.0440098_pr, & ! CO2
       0.01604276_pr, & ! CH4
       0.01503482_pr, & ! CH3
       0.03002628_pr, & ! CH2O
       0.03006964_pr, & ! C2H6
       0.02603788_pr, & ! C2H2
       0.01402688_pr & ! CH2GL1G
      !0.03300674_pr, & ! HO2
      !0.03204216_pr, & ! CH3OH
      !0.03103422_pr, & ! CH3O
      !0.03103422_pr, & ! CH2OH
      !0.02901834_pr, & ! HCO
      !0.0290617_pr, & ! C2H5
      !0.02805376_pr, & ! C2H4
      !0.02704582_pr, & ! C2H3
      !0.04203728_pr & ! CH2CO
       /)
  
contains

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

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

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

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

  end function getlindratecoeff

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

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

    M(mM4) = (1_pr)*c(sCH4) &
         + (2.65_pr)*c(sH2O) &
         + sum(c)

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

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

    M(mM117) = sum(c)

    M(mM136) = sum(c)

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

    M(mM12) = (0.75_pr)*c(sCO) &
         + (2.6_pr)*c(sCO2) &
         + (11_pr)*c(sH2O) &
         + sum(c)

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

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

    M(mM40) = sum(c)

    M(mM41) = sum(c)

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

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

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

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

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

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

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

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

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

    M(mM108) = (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.0508_pr)*exp((-26325.728_pr)*R_T_inv + T_log *(2.67_pr))
    k(r2f) = (43800000_pr)*exp((-29246.16_pr)*R_T_inv)
    k(r3f) = (114000000_pr)*exp((-63956.624_pr)*R_T_inv)
    k(r4f) = (35000000000_pr)*exp(T_log *(-2_pr))
    k(r5f) = (67_pr)*exp((-62704.7712_pr)*R_T_inv + T_log *(1.704_pr))
    k(r6f) = (70790000_pr)*exp((-1234.28_pr)*R_T_inv)
    k(r7f) = (11402_pr)*exp((-2317.01552_pr)*R_T_inv + T_log *(1.0827_pr))
    k(r8f) = (32500000_pr)
    k(r9f) = (7000000_pr)*exp((4572.94464_pr)*R_T_inv)
    k(r10f) = (450000000_pr)*exp((-45729.4464_pr)*R_T_inv)
    k_0(FOr11f) = (17370000_pr)*exp(T_log *(-1.23_pr))
    k_inf(FOr11f) = (4650000_pr)*exp(T_log *(0.44_pr))
    FC(FOr11f) = ((1.0_pr - 0.67_pr)*exp(-Tloc/(1e-30_pr))) + (0.67_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r11f) = getlindratecoeff(Tloc,k_0(FOr11f),k_inf(FOr11f),FC(FOr11f), M(mM11),Ploc)
    k_0(FOr12f) = (1173000000000_pr)*exp((-17535.144_pr)*R_T_inv + T_log *(-2.79_pr))
    k_inf(FOr12f) = (13620_pr)*exp((-9974.656_pr)*R_T_inv)
    k(r12f) = getlindratecoeff(Tloc,k_0(FOr12f),k_inf(FOr12f),1.0_pr, M(mM12),Ploc)
    k(r13f) = (0.07015_pr)*exp((1488.2488_pr)*R_T_inv + T_log *(2.053_pr))
    k(r14f) = (5757000_pr)*exp((-1388.2512_pr)*R_T_inv + T_log *(-0.664_pr))
    k_0(FOr15f) = (2.477e+21_pr)*exp((-10208.96_pr)*R_T_inv + T_log *(-4.76_pr))
    k_inf(FOr15f) = (12700000000_pr)*exp((-1602.472_pr)*R_T_inv + T_log *(-0.63_pr))
    FC(FOr15f) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r15f) = getlindratecoeff(Tloc,k_0(FOr15f),k_inf(FOr15f),FC(FOr15f), M(mM15),Ploc)
    k(r16f) = (0.614_pr)*exp((-40112.008_pr)*R_T_inv + T_log *(2.5_pr))
    k(r17f) = (1020_pr)*exp((-35982.4_pr)*R_T_inv + T_log *(1.5_pr))
    k(r18f) = (0.0583_pr)*exp((-9162.96_pr)*R_T_inv + T_log *(2.6_pr))
    k(r19f) = (0.116_pr)*exp((12644.048_pr)*R_T_inv + T_log *(2.23_pr))
    k(r20f) = (13.26516_pr)*exp((-36632.94924_pr)*R_T_inv + T_log *(1.785451_pr))
    k(r21f) = (0.1746412_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r22f) = (0.07484624_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r23f) = (0.0935578_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r24f) = (0.0935578_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r25f) = (0.4366031_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r26f) = (0.1871156_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r27f) = (0.08732062_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k_0(FOr28f) = (1.72555e+16_pr)*exp((-7177.58924_pr)*R_T_inv + T_log *(-3.354549_pr))
    k_inf(FOr28f) = (134808600000_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-1.014549_pr))
    FC(FOr28f) = ((1.0_pr - 0.68_pr)*exp(-Tloc/(78_pr))) + (0.68_pr*exp(-Tloc/(1995_pr)))&
              + exp(-(5590_pr)/Tloc)
    k(r28f) = getlindratecoeff(Tloc,k_0(FOr28f),k_inf(FOr28f),FC(FOr28f), M(mM28),Ploc)
    k(r29f) = (57158830_pr)*exp((-8307.26924_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r30f) = (14235780_pr)*exp((-8307.26924_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r31f) = (269617100_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r32f) = (7546000_pr)*exp((-118490.88_pr)*R_T_inv)
    k(r33f) = (2.641e-06_pr)*exp((-33911.32_pr)*R_T_inv + T_log *(3.283_pr))
    k(r34f) = (55400000_pr)*exp((569.024_pr)*R_T_inv + T_log *(0.05_pr))
    k(r35f) = (528200000000_pr)*exp((-7414.048_pr)*R_T_inv + T_log *(-1.518_pr))
    k(r36f) = (46860_pr)*exp((-14920.144_pr)*R_T_inv + T_log *(0.833_pr))
    k(r37f) = (1230_pr)*exp((-49998.8_pr)*R_T_inv + T_log *(1.011_pr))
    k(r38f) = (0.04293_pr)*exp((-16726.7952_pr)*R_T_inv + T_log *(2.568_pr))
    k(r39f) = (1000000_pr)*exp((2876.5_pr)*R_T_inv + T_log *(0.269_pr))
    k_0(FOr40f) = (1.5e+37_pr)*exp((-409999.3648_pr)*R_T_inv + T_log *(-6.995_pr))
    k_inf(FOr40f) = (2.084e+18_pr)*exp((-387189.8704_pr)*R_T_inv + T_log *(-0.615_pr))
    FC(FOr40f) = ((1.0_pr - -0.4748_pr)*exp(-Tloc/(35580_pr))) + (-0.4748_pr*exp(-Tloc/(1116_pr)))&
              + exp(-(9023_pr)/Tloc)
    k(r40f) = getlindratecoeff(Tloc,k_0(FOr40f),k_inf(FOr40f),FC(FOr40f), M(mM40),Ploc)
    k_0(FOr41f) = (1.43e+41_pr)*exp((-415961.1464_pr)*R_T_inv + T_log *(-8.227_pr))
    k_inf(FOr41f) = (3.121e+18_pr)*exp((-383723.008_pr)*R_T_inv + T_log *(-1.017_pr))
    FC(FOr41f) = ((1.0_pr - 2.545_pr)*exp(-Tloc/(3290_pr))) + (2.545_pr*exp(-Tloc/(47320_pr)))&
              + exp(-(47110_pr)/Tloc)
    k(r41f) = getlindratecoeff(Tloc,k_0(FOr41f),k_inf(FOr41f),FC(FOr41f), M(mM41),Ploc)
    k(r42f) = (0.199_pr)*exp((-43095.2_pr)*R_T_inv + T_log *(2.56_pr))
    k(r43f) = (0.307_pr)*exp((-22760.96_pr)*R_T_inv + T_log *(2.55_pr))
    k(r44f) = (0.0388_pr)*exp((-12886.72_pr)*R_T_inv + T_log *(2.5_pr))
    k(r45f) = (0.388_pr)*exp((-12886.72_pr)*R_T_inv + T_log *(2.5_pr))
    k(r46f) = (0.00015_pr)*exp((3192.392_pr)*R_T_inv + T_log *(3.03_pr))
    k(r47f) = (0.0308_pr)*exp((3375.2328_pr)*R_T_inv + T_log *(2.65_pr))
    k(r48f) = (0.00322_pr)*exp((-35896.628_pr)*R_T_inv + T_log *(2.425_pr))
    k(r49f) = (1510000000_pr)*exp(T_log *(-1_pr))
    k(r50f) = (241000000_pr)*exp((-20991.128_pr)*R_T_inv)
    k(r51f) = (6000000_pr)
    k(r52f) = (24000000_pr)
    k(r53f) = (42000000_pr)
    k(r54f) = (4.38e-25_pr)*exp((23016.184_pr)*R_T_inv + T_log *(9.5_pr))
    k(r55f) = (20000000_pr)
    k(r56f) = (12000000_pr)
    k_0(FOr57f) = (1350000000000_pr)*exp((-5962.2_pr)*R_T_inv + T_log *(-2.57_pr))
    k_inf(FOr57f) = (1090000_pr)*exp((1087.84_pr)*R_T_inv + T_log *(0.48_pr))
    FC(FOr57f) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r57f) = getlindratecoeff(Tloc,k_0(FOr57f),k_inf(FOr57f),FC(FOr57f), M(mM57),Ploc)
    k_0(FOr58f) = (5.07e+15_pr)*exp((-352912.032_pr)*R_T_inv + T_log *(-3.42_pr))
    k_inf(FOr58f) = (43_pr)*exp((-333046.4_pr)*R_T_inv + T_log *(1.5_pr))
    FC(FOr58f) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r58f) = getlindratecoeff(Tloc,k_0(FOr58f),k_inf(FOr58f),FC(FOr58f), M(mM58),Ploc)
    k(r59f) = (6260_pr)*exp((-9455.84_pr)*R_T_inv + T_log *(1.15_pr))
    k(r60f) = (57.4_pr)*exp((-11464.16_pr)*R_T_inv + T_log *(1.9_pr))
    k(r61f) = (78.2_pr)*exp((4414.12_pr)*R_T_inv + T_log *(1.63_pr))
    k(r62f) = (3.83e-05_pr)*exp((-18041.408_pr)*R_T_inv + T_log *(3.36_pr))
    k(r63f) = (570000_pr)*exp((-62216.08_pr)*R_T_inv + T_log *(0.66_pr))
    k(r64f) = (7580000_pr)*exp((-1715.44_pr)*R_T_inv)
    k(r65f) = (30200000_pr)
    k(r66f) = (73400000_pr)
    k(r67f) = (30110000_pr)
    k(r68f) = (26500000_pr)
    k(r69f) = (30000000_pr)
    k_0(FOr70f) = (1.27e+20_pr)*exp((-27321.52_pr)*R_T_inv + T_log *(-4.82_pr))
    k_inf(FOr70f) = (540000_pr)*exp((-15062.4_pr)*R_T_inv + T_log *(0.454_pr))
    FC(FOr70f) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r70f) = getlindratecoeff(Tloc,k_0(FOr70f),k_inf(FOr70f),FC(FOr70f), M(mM70),Ploc)
    k_0(FOr71f) = (1.867e+19_pr)*exp((-101700.488_pr)*R_T_inv + T_log *(-3_pr))
    k_inf(FOr71f) = (68000000000000_pr)*exp((-109495.28_pr)*R_T_inv)
    FC(FOr71f) = ((1.0_pr - 0.9_pr)*exp(-Tloc/(2500_pr))) + (0.9_pr*exp(-Tloc/(1300_pr)))&
              + exp(-(1e+99_pr)/Tloc)
    k(r71f) = getlindratecoeff(Tloc,k_0(FOr71f),k_inf(FOr71f),FC(FOr71f), M(mM71),Ploc)
    k_0(FOr72f) = (8.054e+19_pr)*exp((-4107.0144_pr)*R_T_inv + T_log *(-3.75_pr))
    k_inf(FOr72f) = (2277000000_pr)*exp((-731.7816_pr)*R_T_inv + T_log *(-0.69_pr))
    FC(FOr72f) = ((1.0_pr - 0_pr)*exp(-Tloc/(570_pr))) + (0_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r72f) = getlindratecoeff(Tloc,k_0(FOr72f),k_inf(FOr72f),FC(FOr72f), M(mM72),Ploc)
    k_0(FOr73f) = (1.99e+29_pr)*exp((-27970.04_pr)*R_T_inv + T_log *(-7.08_pr))
    k_inf(FOr73f) = (521000000000_pr)*exp((-6610.72_pr)*R_T_inv + T_log *(-0.99_pr))
    FC(FOr73f) = ((1.0_pr - 0.842_pr)*exp(-Tloc/(125_pr))) + (0.842_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r73f) = getlindratecoeff(Tloc,k_0(FOr73f),k_inf(FOr73f),FC(FOr73f), M(mM73),Ploc)
    k(r74f) = (3.55_pr)*exp((-24392.72_pr)*R_T_inv + T_log *(2.4_pr))
    k(r75f) = (115_pr)*exp((-31505.52_pr)*R_T_inv + T_log *(1.9_pr))
    k(r76f) = (14.8_pr)*exp((-3974.8_pr)*R_T_inv + T_log *(1.9_pr))
    k(r77f) = (5.55e-10_pr)*exp((-13518.504_pr)*R_T_inv + T_log *(4.72_pr))
    k_0(FOr78f) = (1.419e+27_pr)*exp((-24137.496_pr)*R_T_inv + T_log *(-6.642_pr))
    k_inf(FOr78f) = (956.9_pr)*exp((-5669.32_pr)*R_T_inv + T_log *(1.463_pr))
    FC(FOr78f) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299_pr))) + (-0.569_pr*exp(-Tloc/(-9147_pr)))&
              + exp(-(152.4_pr)/Tloc)
    k(r78f) = getlindratecoeff(Tloc,k_0(FOr78f),k_inf(FOr78f),FC(FOr78f), M(mM78),Ploc)
    k(r79f) = (0.0118_pr)*exp((12221.464_pr)*R_T_inv + T_log *(2.45_pr))
    k(r80f) = (310000000_pr)*exp((-55950.54_pr)*R_T_inv + T_log *(-0.362_pr))
    k_0(FOr81f) = (1.4e+18_pr)*exp((-13890.88_pr)*R_T_inv + T_log *(-3.86_pr))
    k_inf(FOr81f) = (6080000_pr)*exp((-1171.52_pr)*R_T_inv + T_log *(0.27_pr))
    FC(FOr81f) = ((1.0_pr - 0.782_pr)*exp(-Tloc/(207.5_pr))) + (0.782_pr*exp(-Tloc/(2663_pr)))&
              + exp(-(6095_pr)/Tloc)
    k(r81f) = getlindratecoeff(Tloc,k_0(FOr81f),k_inf(FOr81f),FC(FOr81f), M(mM81),Ploc)
    k(r82f) = (26100000000_pr)*exp((-283771.432_pr)*R_T_inv)
    k(r83f) = (50.7_pr)*exp((-54182.8_pr)*R_T_inv + T_log *(1.93_pr))
    k(r84f) = (0.0223_pr)*exp((-9269.652_pr)*R_T_inv + T_log *(2.745_pr))
    k(r85f) = (0.000976_pr)*exp((-63379.232_pr)*R_T_inv + T_log *(2.947_pr))
    k(r86f) = (8.13e-11_pr)*exp((-36968.9872_pr)*R_T_inv + T_log *(4.417_pr))
    k(r87f) = (7.453_pr)*exp((-765.672_pr)*R_T_inv + T_log *(1.88_pr))
    k(r88f) = (0.1247437_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k_0(FOr89f) = (6.346e+19_pr)*exp((-15815.52_pr)*R_T_inv + T_log *(-4.664_pr))
    k_inf(FOr89f) = (17100_pr)*exp((-11334.456_pr)*R_T_inv + T_log *(1.266_pr))
    FC(FOr89f) = ((1.0_pr - 0.788_pr)*exp(-Tloc/(-10200_pr))) + (0.788_pr*exp(-Tloc/(1e-30_pr)))
    k(r89f) = getlindratecoeff(Tloc,k_0(FOr89f),k_inf(FOr89f),FC(FOr89f), M(mM89),Ploc)
    k(r90f) = (21.5_pr)*exp((-14087.528_pr)*R_T_inv + T_log *(1.19_pr))
    k(r91f) = (4.6e-05_pr)*exp((2061.8752_pr)*R_T_inv + T_log *(2.76_pr))
    k(r92f) = (0.00106_pr)*exp((-25857.12_pr)*R_T_inv + T_log *(2.39_pr))
    k(r93f) = (5.26e-07_pr)*exp((-7434.968_pr)*R_T_inv + T_log *(3.01_pr))
    k(r94f) = (2.73e+29_pr)*exp((-49454.88_pr)*R_T_inv + T_log *(-7.32_pr))
    k(r95f) = (6080000000_pr)*exp((-2701.6088_pr)*R_T_inv + T_log *(-1.31_pr))
    k(r96f) = (6.36e+29_pr)*exp((-49454.88_pr)*R_T_inv + T_log *(-7.32_pr))
    k(r97f) = (14200000000_pr)*exp((-2701.6088_pr)*R_T_inv + T_log *(-1.31_pr))
    k(r98f) = (103000_pr)*exp((3128.7952_pr)*R_T_inv + T_log *(-0.33_pr))
    k(r99f) = (5.77e+15_pr)*exp((-19966.048_pr)*R_T_inv + T_log *(-3.54_pr))
    k(r100f) = (7.25e+25_pr)*exp((-43680.96_pr)*R_T_inv + T_log *(-6.7_pr))
    k(r101f) = (53200000_pr)*exp((-1868.9928_pr)*R_T_inv + T_log *(-1.14_pr))
    k(r102f) = (96352500_pr)
    k(r103f) = (30110000_pr)
    k(r104f) = (739.5_pr)*exp((-10342.848_pr)*R_T_inv + T_log *(1.28_pr))
    k(r105f) = (10_pr)*exp((-25104_pr)*R_T_inv + T_log *(2_pr))
    k(r106f) = (7.528_pr)*exp((-8811.504_pr)*R_T_inv + T_log *(1.55_pr))
    k(r107f) = (1277_pr)*exp((-10790.536_pr)*R_T_inv + T_log *(0.73_pr))
    k_0(FOr108f) = (3000000000_pr)*exp((-238488_pr)*R_T_inv)
    k_inf(FOr108f) = (30000000000000_pr)*exp((-297064_pr)*R_T_inv)
    k(r108f) = getlindratecoeff(Tloc,k_0(FOr108f),k_inf(FOr108f),1.0_pr, M(mM108),Ploc)
    k(r109f) = (30000000_pr)
    k(r110f) = (777_pr)*exp((-11631.52_pr)*R_T_inv + T_log *(1.45_pr))
    k(r111f) = (2000000_pr)*exp((4225.84_pr)*R_T_inv)
    k(r112f) = (4.16e+18_pr)*exp((-40835.84_pr)*R_T_inv + T_log *(-3.19_pr))
    k(r113f) = (4920_pr)*exp((-2509.43768_pr)*R_T_inv + T_log *(1.02_pr))
    k(r114f) = (10000000_pr)*exp((-62760_pr)*R_T_inv)
    k(r115f) = (20000000_pr)*exp((-9623.2_pr)*R_T_inv)
    k(r116f) = (10000000_pr)
    k(r117f) = (10000_pr)
    k(r118f) = (100000000_pr)*exp((-100416_pr)*R_T_inv)
    k(r119f) = (107846900_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r120f) = (269617100_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r121f) = (161770300_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r122f) = (593157.7_pr)*exp((-6215.26924_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r123f) = (0.1871156_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r124f) = (0.268199_pr)*exp((17630.744216_pr)*R_T_inv + T_log *(1.988656_pr))
    k(r125f) = (500000000_pr)*exp((-133888_pr)*R_T_inv)
    k(r126f) = (647081100_pr)*exp((-2031.2675664_pr)*R_T_inv + T_log *(-0.2145491_pr))
    k(r127f) = (30000000000000_pr)*exp((-328444_pr)*R_T_inv)
    k(r128f) = (20000000_pr)*exp((-20920_pr)*R_T_inv)
    k(r129f) = (100000_pr)
    k(r130f) = (5000000_pr)*exp((-62760_pr)*R_T_inv)
    k(r131f) = (500000_pr)*exp((-27196_pr)*R_T_inv)
    k(r132f) = (300000_pr)*exp((-108784_pr)*R_T_inv)
    k(r133f) = (100000000_pr)*exp((-200832_pr)*R_T_inv)
    k(r134f) = (300000_pr)*exp((-112968_pr)*R_T_inv)
    k(r135f) = (100000000_pr)*exp((-179912_pr)*R_T_inv)
    k(r136f) = (50000_pr)
    k(r137f) = (15000000_pr)
    k(r138f) = (10.83_pr)*exp((-36743.72064_pr)*R_T_inv + T_log *(2_pr))
    k(r1b) = (0.0270125153693532_pr)*exp((-20345.1738760158_pr)*R_T_inv + T_log *(2.64832403899644_pr))
    k(r2b) = (284406184.196714_pr)*exp((-91248.0139363202_pr)*R_T_inv + T_log *(-0.0442320199425153_pr))
    k(r3b) = (440038.685960121_pr)*exp((5920.90113184252_pr)*R_T_inv + T_log *(0.354141305726111_pr))
    k(r4b) = (6.68351087968594e+16_pr)*exp((-495069.107515644_pr)*R_T_inv + T_log *(-1.72920006958837_pr))
    k(r5b) = (5.4866995655307_pr)*exp((5277.63688551364_pr)*R_T_inv + T_log *(1.72655606236567_pr))
    k(r6b) = (39742.8681842056_pr)*exp((-154234.588683103_pr)*R_T_inv + T_log *(0.633270526293657_pr))
    k(r7b) = (3118.75340569907_pr)*exp((-231175.403510548_pr)*R_T_inv + T_log *(1.38350517297583_pr))
    k(r8b) = (4726990.61993216_pr)*exp((-222877.833974141_pr)*R_T_inv + T_log *(0.279129192740927_pr))
    k(r9b) = (12432627.4077715_pr)*exp((-286287.297212063_pr)*R_T_inv + T_log *(0.256573167246556_pr))
    k(r10b) = (799240369.647056_pr)*exp((-336589.688283625_pr)*R_T_inv + T_log *(0.256573161359728_pr))
    k_0(FOr11b) = (18675472030069.4_pr)*exp((-204208.865719999_pr)*R_T_inv + T_log *(-1.21577324658642_pr))
    k_inf(FOr11b) = (4999477905722.04_pr)*exp((-204208.865604176_pr)*R_T_inv + T_log *(0.454226773906472_pr))
    FC(FOr11b) = ((1.0_pr - 0.67_pr)*exp(-Tloc/(1e-30_pr))) + (0.67_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r11b) = getlindratecoeff(Tloc,k_0(FOr11b),k_inf(FOr11b),FC(FOr11b), M(mM11),Ploc)
    k_0(FOr12b) = (3.18037647702451e+23_pr)*exp((-550827.478381066_pr)*R_T_inv + T_log *(-3.63255340053932_pr))
    k_inf(FOr12b) = (3.69281562074678e+15_pr)*exp((-543266.990389345_pr)*R_T_inv + T_log *(-0.842553399135274_pr))
    k(r12b) = getlindratecoeff(Tloc,k_0(FOr12b),k_inf(FOr12b),1.0_pr, M(mM12),Ploc)
    k(r13b) = (121628.458041892_pr)*exp((-104717.38614579_pr)*R_T_inv + T_log *(0.917090607268091_pr))
    k(r14b) = (9981682452645.99_pr)*exp((-107593.886136258_pr)*R_T_inv + T_log *(-1.79990939109523_pr))
    k_0(FOr15b) = (3.55966249336337e+29_pr)*exp((-449910.022521636_pr)*R_T_inv + T_log *(-4.82105412018431_pr))
    k_inf(FOr15b) = (1.82510033376441e+18_pr)*exp((-441303.534866241_pr)*R_T_inv + T_log *(-0.691054182340573_pr))
    FC(FOr15b) = ((1.0_pr - 0.783_pr)*exp(-Tloc/(74_pr))) + (0.783_pr*exp(-Tloc/(2941_pr)))&
              + exp(-(6964_pr)/Tloc)
    k(r15b) = getlindratecoeff(Tloc,k_0(FOr15b),k_inf(FOr15b),FC(FOr15b), M(mM15),Ploc)
    k(r16b) = (0.00125648475023843_pr)*exp((-33478.1990166204_pr)*R_T_inv + T_log *(2.87608607731348_pr))
    k(r17b) = (1.10991655856455_pr)*exp((-23368.0368797288_pr)*R_T_inv + T_log *(1.85441011832776_pr))
    k(r18b) = (0.000774679937512472_pr)*exp((-64531.0049185985_pr)*R_T_inv + T_log *(2.93185406400405_pr))
    k(r19b) = (15.5049085026644_pr)*exp((-222848.149091627_pr)*R_T_inv + T_log *(2.15471907516999_pr))
    k(r20b) = (8.55376078647189_pr)*exp((-76325.7852488525_pr)*R_T_inv + T_log *(1.7150021341284_pr))
    k(r22b) = (0.0218169204160266_pr)*exp((-742114.415019099_pr)*R_T_inv + T_log *(2.45349676597155_pr))
    k(r23b) = (1.08806343525867_pr)*exp((-749990.089069758_pr)*R_T_inv + T_log *(2.14358698707589_pr))
    k(r24b) = (21.3051750330677_pr)*exp((-381687.402378174_pr)*R_T_inv + T_log *(1.56511894301452_pr))
    k(r25b) = (137.576028590874_pr)*exp((-28695.9011120244_pr)*R_T_inv + T_log *(1.5421209983148_pr))
    k(r26b) = (5650.65924786369_pr)*exp((-322882.664065485_pr)*R_T_inv + T_log *(1.40528995624174_pr))
    k(r27b) = (0.00152089261784454_pr)*exp((-216677.029366013_pr)*R_T_inv + T_log *(2.54119930585758_pr))
    k_0(FOr28b) = (1.59902508054954e+24_pr)*exp((-486571.488284621_pr)*R_T_inv + T_log *(-3.48605207668191_pr))
    k_inf(FOr28b) = (1.24923914046432e+19_pr)*exp((-481425.167043674_pr)*R_T_inv + T_log *(-1.14605215337338_pr))
    FC(FOr28b) = ((1.0_pr - 0.68_pr)*exp(-Tloc/(78_pr))) + (0.68_pr*exp(-Tloc/(1995_pr)))&
              + exp(-(5590_pr)/Tloc)
    k(r28b) = getlindratecoeff(Tloc,k_0(FOr28b),k_inf(FOr28b),FC(FOr28b), M(mM28),Ploc)
    k(r29b) = (50242869.4842222_pr)*exp((-337747.89068388_pr)*R_T_inv + T_log *(-0.28394484866537_pr))
    k(r32b) = (6897078696.29131_pr)*exp((1639.0656183633_pr)*R_T_inv + T_log *(-0.882623673913189_pr))
    k(r33b) = (5.19504096539813e-07_pr)*exp((-252240.003839691_pr)*R_T_inv + T_log *(3.47863427753638_pr))
    k(r34b) = (2823217993.07441_pr)*exp((-287637.185103268_pr)*R_T_inv + T_log *(-0.108507052499265_pr))
    k(r35b) = (10884472418.0727_pr)*exp((-23089.2566375588_pr)*R_T_inv + T_log *(-1.11569702249245_pr))
    k(r36b) = (125596760.184155_pr)*exp((-682.723825590269_pr)*R_T_inv + T_log *(0.0235446986626907_pr))
    k(r37b) = (291251297.340444_pr)*exp((253.619773138916_pr)*R_T_inv + T_log *(-0.225765092305187_pr))
    k(r38b) = (0.000884646729248711_pr)*exp((-32402.0038399276_pr)*R_T_inv + T_log *(2.97030297709788_pr))
    k(r39b) = (132938171.240877_pr)*exp((-99871.3882517169_pr)*R_T_inv + T_log *(-0.334494463385433_pr))
    k_0(FOr40b) = (7.08968109392954e+26_pr)*exp((-23137.8995271779_pr)*R_T_inv + T_log *(-6.47138034090048_pr))
    k_inf(FOr40b) = (98498953.6478154_pr)*exp((-328.40248741695_pr)*R_T_inv + T_log *(-0.0913798826133818_pr))
    FC(FOr40b) = ((1.0_pr - -0.4748_pr)*exp(-Tloc/(35580_pr))) + (-0.4748_pr*exp(-Tloc/(1116_pr)))&
              + exp(-(9023_pr)/Tloc)
    k(r40b) = getlindratecoeff(Tloc,k_0(FOr40b),k_inf(FOr40b),FC(FOr40b), M(mM40),Ploc)
    k_0(FOr41b) = (1.39277349970087e+29_pr)*exp((-44774.8898177605_pr)*R_T_inv + T_log *(-7.30107737222089_pr))
    k_inf(FOr41b) = (3039751.26313933_pr)*exp((-12536.7511131587_pr)*R_T_inv + T_log *(-0.0910773190259202_pr))
    FC(FOr41b) = ((1.0_pr - 2.545_pr)*exp(-Tloc/(3290_pr))) + (2.545_pr*exp(-Tloc/(47320_pr)))&
              + exp(-(47110_pr)/Tloc)
    k(r41b) = getlindratecoeff(Tloc,k_0(FOr41b),k_inf(FOr41b),FC(FOr41b), M(mM41),Ploc)
    k(r42b) = (0.65497184004023_pr)*exp((-39048.5678966186_pr)*R_T_inv + T_log *(2.16188661530813_pr))
    k(r43b) = (0.0114372638711732_pr)*exp((-54729.3289222633_pr)*R_T_inv + T_log *(2.57919616339815_pr))
    k(r44b) = (0.0679051319171975_pr)*exp((-2859.5337808263_pr)*R_T_inv + T_log *(2.08021065277604_pr))
    k(r45b) = (0.0076862925500728_pr)*exp((-38874.5349756581_pr)*R_T_inv + T_log *(2.50752017122843_pr))
    k(r46b) = (0.00320572123360217_pr)*exp((-54762.8298673442_pr)*R_T_inv + T_log *(2.58765459020606_pr))
    k(r47b) = (0.00745074149636707_pr)*exp((-90594.9902135496_pr)*R_T_inv + T_log *(2.6349641171484_pr))
    k(r48b) = (5.17888892525164_pr)*exp((-38483.8049220901_pr)*R_T_inv + T_log *(1.6508005309345_pr))
    k(r49b) = (197393852.201591_pr)*exp((-79565.794585665_pr)*R_T_inv + T_log *(-0.628180811219007_pr))
    k(r50b) = (31504582.3604037_pr)*exp((-100556.922600545_pr)*R_T_inv + T_log *(0.371819186085481_pr))
    k(r51b) = (214539.82949832_pr)*exp((-308424.182583617_pr)*R_T_inv + T_log *(0.672624360223304_pr))
    k(r52b) = (5572279.07458404_pr)*exp((-370426.03654224_pr)*R_T_inv + T_log *(0.628392336959043_pr))
    k(r53b) = (798559.528582209_pr)*exp((-302443.62847706_pr)*R_T_inv + T_log *(0.650948396160461_pr))
    k(r54b) = (6.4810448701502e-28_pr)*exp((-92564.6116593055_pr)*R_T_inv + T_log *(10.2991287282051_pr))
    k(r55b) = (8094.70252911178_pr)*exp((-344439.183615062_pr)*R_T_inv + T_log *(1.09993390729457_pr))
    k(r56b) = (2373357.88972685_pr)*exp((-351072.992499742_pr)*R_T_inv + T_log *(0.723847848187727_pr))
    k_0(FOr57b) = (2.79956727712364e+19_pr)*exp((-374244.161340104_pr)*R_T_inv + T_log *(-2.43647303305687_pr))
    k_inf(FOr57b) = (22603900222366.7_pr)*exp((-367194.120908147_pr)*R_T_inv + T_log *(0.613527043593833_pr))
    FC(FOr57b) = ((1.0_pr - 0.7824_pr)*exp(-Tloc/(271_pr))) + (0.7824_pr*exp(-Tloc/(2755_pr)))&
              + exp(-(6570_pr)/Tloc)
    k(r57b) = getlindratecoeff(Tloc,k_0(FOr57b),k_inf(FOr57b),FC(FOr57b), M(mM57),Ploc)
    k_0(FOr58b) = (2.05872237832027e+24_pr)*exp((-352891.309687275_pr)*R_T_inv + T_log *(-3.86494160203652_pr))
    k_inf(FOr58b) = (17460562461.8546_pr)*exp((-333025.677599799_pr)*R_T_inv + T_log *(1.05505841349847_pr))
    FC(FOr58b) = ((1.0_pr - 0.932_pr)*exp(-Tloc/(197_pr))) + (0.932_pr*exp(-Tloc/(1540_pr)))&
              + exp(-(10300_pr)/Tloc)
    k(r58b) = getlindratecoeff(Tloc,k_0(FOr58b),k_inf(FOr58b),FC(FOr58b), M(mM58),Ploc)
    k(r59b) = (47.2052451606734_pr)*exp((-68260.5782120931_pr)*R_T_inv + T_log *(1.3098290044007_pr))
    k(r60b) = (0.814000523602517_pr)*exp((-76249.4492400561_pr)*R_T_inv + T_log *(2.08150550714807_pr))
    k(r61b) = (7.20089154986084_pr)*exp((-122373.026245856_pr)*R_T_inv + T_log *(1.76727295083412_pr))
    k(r62b) = (0.000265414442501381_pr)*exp((-89460.5093517373_pr)*R_T_inv + T_log *(3.16541888256575_pr))
    k(r63b) = (0.0989852573209113_pr)*exp((2548.48985998606_pr)*R_T_inv + T_log *(0.923436607890063_pr))
    k(r64b) = (1415261.02707123_pr)*exp((-141159.735711162_pr)*R_T_inv + T_log *(0.27766338762487_pr))
    k(r65b) = (820116.756876822_pr)*exp((-362322.12972001_pr)*R_T_inv + T_log *(0.55679257392235_pr))
    k(r66b) = (3748550.52369545_pr)*exp((-368302.68363775_pr)*R_T_inv + T_log *(0.578468570470829_pr))
    k(r67b) = (9984882.65241001_pr)*exp((-430304.53754352_pr)*R_T_inv + T_log *(0.534236558017052_pr))
    k(r68b) = (661339396.332628_pr)*exp((-374936.492730935_pr)*R_T_inv + T_log *(0.20238247447273_pr))
    k(r69b) = (1412529373431.91_pr)*exp((-468527.764606979_pr)*R_T_inv + T_log *(-0.579116808239996_pr))
    k_0(FOr70b) = (1.04452486182596e+27_pr)*exp((-151964.591078037_pr)*R_T_inv + T_log *(-5.17759242540463_pr))
    k_inf(FOr70b) = (4441286692333.91_pr)*exp((-139705.471057965_pr)*R_T_inv + T_log *(0.0964075781535183_pr))
    FC(FOr70b) = ((1.0_pr - 0.7187_pr)*exp(-Tloc/(103_pr))) + (0.7187_pr*exp(-Tloc/(1291_pr)))&
              + exp(-(4160_pr)/Tloc)
    k(r70b) = getlindratecoeff(Tloc,k_0(FOr70b),k_inf(FOr70b),FC(FOr70b), M(mM70),Ploc)
    k_0(FOr71b) = (25694696089.0061_pr)*exp((-13072.4179712031_pr)*R_T_inv + T_log *(-2.21509803082354_pr))
    k_inf(FOr71b) = (93585.385930726_pr)*exp((-20867.2098974968_pr)*R_T_inv + T_log *(0.784901982487235_pr))
    FC(FOr71b) = ((1.0_pr - 0.9_pr)*exp(-Tloc/(2500_pr))) + (0.9_pr*exp(-Tloc/(1300_pr)))&
              + exp(-(1e+99_pr)/Tloc)
    k(r71b) = getlindratecoeff(Tloc,k_0(FOr71b),k_inf(FOr71b),FC(FOr71b), M(mM71),Ploc)
    k_0(FOr72b) = (1.04980509386092e+33_pr)*exp((-384971.434484873_pr)*R_T_inv + T_log *(-4.82216580254776_pr))
    k_inf(FOr72b) = (2.96797438205629e+22_pr)*exp((-381596.20179663_pr)*R_T_inv + T_log *(-1.76216582259599_pr))
    FC(FOr72b) = ((1.0_pr - 0_pr)*exp(-Tloc/(570_pr))) + (0_pr*exp(-Tloc/(1e+30_pr)))&
              + exp(-(1e+30_pr)/Tloc)
    k(r72b) = getlindratecoeff(Tloc,k_0(FOr72b),k_inf(FOr72b),FC(FOr72b), M(mM72),Ploc)
    k_0(FOr73b) = (5.14229612938501e+37_pr)*exp((-449030.290716163_pr)*R_T_inv + T_log *(-7.16781993143772_pr))
    k_inf(FOr73b) = (1.34629944718175e+20_pr)*exp((-427670.970617583_pr)*R_T_inv + T_log *(-1.07781991266864_pr))
    FC(FOr73b) = ((1.0_pr - 0.842_pr)*exp(-Tloc/(125_pr))) + (0.842_pr*exp(-Tloc/(2219_pr)))&
              + exp(-(6882_pr)/Tloc)
    k(r73b) = getlindratecoeff(Tloc,k_0(FOr73b),k_inf(FOr73b),FC(FOr73b), M(mM73),Ploc)
    k(r74b) = (0.00214831185579349_pr)*exp((-30419.168690166_pr)*R_T_inv + T_log *(2.78117592841327_pr))
    k(r75b) = (0.130877667139832_pr)*exp((-43512.522893964_pr)*R_T_inv + T_log *(2.30285187523546_pr))
    k(r76b) = (0.109369018638292_pr)*exp((-77983.6567838101_pr)*R_T_inv + T_log *(2.25861986443174_pr))
    k(r77b) = (3.08653942532647e-10_pr)*exp((-32159.3158644512_pr)*R_T_inv + T_log *(4.74676580047897_pr))
    k_0(FOr78b) = (2.98251299762534e+33_pr)*exp((-173580.77080033_pr)*R_T_inv + T_log *(-6.72871759581829_pr))
    k_inf(FOr78b) = (2011251904.28879_pr)*exp((-155112.594740908_pr)*R_T_inv + T_log *(1.37628241491903_pr))
    FC(FOr78b) = ((1.0_pr - -0.569_pr)*exp(-Tloc/(299_pr))) + (-0.569_pr*exp(-Tloc/(-9147_pr)))&
              + exp(-(152.4_pr)/Tloc)
    k(r78b) = getlindratecoeff(Tloc,k_0(FOr78b),k_inf(FOr78b),FC(FOr78b), M(mM78),Ploc)
    k(r79b) = (0.806813550839828_pr)*exp((-278036.337863386_pr)*R_T_inv + T_log *(2.47566099866483_pr))
    k(r80b) = (15637043525798.5_pr)*exp((-15754.7092090043_pr)*R_T_inv + T_log *(-1.34634584345038_pr))
    k_0(FOr81b) = (5.83595097084938e+26_pr)*exp((-475779.76525933_pr)*R_T_inv + T_log *(-3.99965554899756_pr))
    k_inf(FOr81b) = (2.53448096719882e+15_pr)*exp((-463060.408366002_pr)*R_T_inv + T_log *(0.130343895940619_pr))
    FC(FOr81b) = ((1.0_pr - 0.782_pr)*exp(-Tloc/(207.5_pr))) + (0.782_pr*exp(-Tloc/(2663_pr)))&
              + exp(-(6095_pr)/Tloc)
    k(r81b) = getlindratecoeff(Tloc,k_0(FOr81b),k_inf(FOr81b),FC(FOr81b), M(mM81),Ploc)
    k(r82b) = (709.258441650225_pr)*exp((-107187.457968666_pr)*R_T_inv + T_log *(0.00496995391731932_pr))
    k(r83b) = (0.0357681188181754_pr)*exp((-25361.1683975375_pr)*R_T_inv + T_log *(2.38468748439084_pr))
    k(r84b) = (0.000102154614520924_pr)*exp((-42449.8744422983_pr)*R_T_inv + T_log *(3.15545544645163_pr))
    k(r85b) = (0.000336472187115447_pr)*exp((-41191.4094919261_pr)*R_T_inv + T_log *(3.02560138802983_pr))
    k(r86b) = (2.80278567460698e-11_pr)*exp((-14781.1646751653_pr)*R_T_inv + T_log *(4.49560139101445_pr))
    k(r87b) = (0.000763159931003861_pr)*exp((-110329.025244877_pr)*R_T_inv + T_log *(2.48559424859911_pr))
    k(r88b) = (277420.898822613_pr)*exp((-272124.048943009_pr)*R_T_inv + T_log *(0.95952472846128_pr))
    k_0(FOr89b) = (1.6475713443409e+24_pr)*exp((-163577.89612925_pr)*R_T_inv + T_log *(-4.21428845608793_pr))
    k_inf(FOr89b) = (443956332.356998_pr)*exp((-159096.832099461_pr)*R_T_inv + T_log *(1.71571154914153_pr))
    FC(FOr89b) = ((1.0_pr - 0.788_pr)*exp(-Tloc/(-10200_pr))) + (0.788_pr*exp(-Tloc/(1e-30_pr)))
    k(r89b) = getlindratecoeff(Tloc,k_0(FOr89b),k_inf(FOr89b),FC(FOr89b), M(mM89),Ploc)
    k(r90b) = (890.360905277887_pr)*exp((-70534.0176164249_pr)*R_T_inv + T_log *(0.75451520490111_pr))
    k(r91b) = (0.00190495814659086_pr)*exp((-54384.614389779_pr)*R_T_inv + T_log *(2.32451520962494_pr))
    k(r92b) = (0.0165583544413694_pr)*exp((-333632.933069291_pr)*R_T_inv + T_log *(1.98854151208298_pr))
    k(r93b) = (8.2166911980476e-06_pr)*exp((-315210.780922358_pr)*R_T_inv + T_log *(2.60854153837122_pr))
    k(r94b) = (1.46581166430223e+28_pr)*exp((-412149.102752964_pr)*R_T_inv + T_log *(-6.95178298750225_pr))
    k(r95b) = (326451822.939102_pr)*exp((-365395.831540798_pr)*R_T_inv + T_log *(-0.94178298526351_pr))
    k(r98b) = (697829.348586601_pr)*exp((-383428.927695747_pr)*R_T_inv + T_log *(-0.483248344698474_pr))
    k(r99b) = (3.90919933047341e+16_pr)*exp((-406523.77089201_pr)*R_T_inv + T_log *(-3.69324834367174_pr))
    k(r100b) = (3.59662516745663e+27_pr)*exp((-586696.738355013_pr)*R_T_inv + T_log *(-6.75239276185443_pr))
    k(r101b) = (2639179154.3965_pr)*exp((-544884.77126834_pr)*R_T_inv + T_log *(-1.19239278201826_pr))
    k(r102b) = (1091416664.49225_pr)*exp((-285304.877556_pr)*R_T_inv + T_log *(-0.134679613206917_pr))
    k(r103b) = (2214640326.89386_pr)*exp((-347306.731377978_pr)*R_T_inv + T_log *(-0.178911612135578_pr))
    k(r104b) = (0.00153574202163739_pr)*exp((-185398.993654445_pr)*R_T_inv + T_log *(2.59058994106078_pr))
    k(r105b) = (0.0450857977514186_pr)*exp((-108101.806160688_pr)*R_T_inv + T_log *(2.71314817099975_pr))
    k(r106b) = (5057.96615697832_pr)*exp((-107140.518630986_pr)*R_T_inv + T_log *(0.950755806119762_pr))
    k(r107b) = (1.5715416935799_pr)*exp((-238153.881586271_pr)*R_T_inv + T_log *(1.61573081632511_pr))
    k_0(FOr108b) = (5.92968920988823e-05_pr)*exp((111871.56859445_pr)*R_T_inv + T_log *(1.61647817544697_pr))
    k_inf(FOr108b) = (0.592968989022686_pr)*exp((53295.5685112857_pr)*R_T_inv + T_log *(1.61647816053834_pr))
    k(r108b) = getlindratecoeff(Tloc,k_0(FOr108b),k_inf(FOr108b),1.0_pr, M(mM108),Ploc)
    k(r109b) = (121407763123.891_pr)*exp((-377653.338039168_pr)*R_T_inv + T_log *(-0.755599764020485_pr))
    k(r110b) = (0.00142317865083233_pr)*exp((-140665.850444999_pr)*R_T_inv + T_log *(2.93497509966357_pr))
    k(r111b) = (9818.49631506828_pr)*exp((-110571.070978953_pr)*R_T_inv + T_log *(0.675519668622081_pr))
    k(r112b) = (9.24777926438917e+17_pr)*exp((-58340.8537784271_pr)*R_T_inv + T_log *(-3.18184752639103_pr))
    k(r113b) = (1093.72773558332_pr)*exp((-20014.4514555522_pr)*R_T_inv + T_log *(1.02815247412094_pr))
    k(r117b) = (2785387303049.38_pr)*exp((-274086.390569059_pr)*R_T_inv + T_log *(-0.339914501472034_pr))
    k(r119b) = (1735524986.82584_pr)*exp((-416660.596440834_pr)*R_T_inv + T_log *(-0.0826155138966575_pr))
    k(r120b) = (3135599634.73594_pr)*exp((-769652.099137342_pr)*R_T_inv + T_log *(-0.0596178136936434_pr))
    k(r121b) = (4885262789347.18_pr)*exp((-342544.675877501_pr)*R_T_inv + T_log *(-0.797915148892618_pr))
    k(r122b) = (10331.2221530176_pr)*exp((-240523.042462496_pr)*R_T_inv + T_log *(0.337994269658532_pr))
    k(r123b) = (7.39964231226105e-06_pr)*exp((-316922.832459067_pr)*R_T_inv + T_log *(1.82855556166483_pr))
    k(r124b) = (0.172942573343235_pr)*exp((-22062.0920581944_pr)*R_T_inv + T_log *(1.91820708809775_pr))
    k(r127b) = (16243.8433295383_pr)*exp((-191007.727726374_pr)*R_T_inv + T_log *(0.489569536169723_pr))
    k(r131b) = (3661130.85351629_pr)*exp((-183654.055593167_pr)*R_T_inv + T_log *(0.100855558813585_pr))
    k(r138b) = (0.00406273036585156_pr)*exp((-1941.53506758441_pr)*R_T_inv + T_log *(2.43301149708684_pr))
  
    return
  end subroutine get_rate_coefficients

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

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

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

    real(pr), dimension(nspec) :: cdot
    real(pr), dimension(nreac + nreac_reverse) :: w

    cdot(sN2) = 0.0_pr 

    cdot(sH2) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         - w(r2f) &
         + w(r2b) &
         + w(r7f) &
         - w(r7b) &
         + w(r16f) &
         - w(r16b) &
         + w(r23f) &
         - w(r23b) &
         - w(r25f) &
         + w(r25b) &
         + w(r42f) &
         - w(r42b) &
         + w(r43f) &
         - w(r43b) &
         + w(r51f) &
         - w(r51b) &
         + w(r55f) &
         - w(r55b) &
         - w(r58f) &
         + w(r58b) &
         + w(r60f) &
         - w(r60b) &
         + w(r66f) &
         - w(r66b) &
         + w(r75f) &
         - w(r75b) &
         + w(r82f) &
         - w(r82b) &
         + w(r83f) &
         - w(r83b) &
         + w(r102f) &
         - w(r102b) &
         + w(r120f) &
         - w(r120b) &
         + w(r125f) &
         + w(r127f) &
         - w(r127b) &
         + w(r129f) 

    cdot(sH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r2f) &
         - w(r2b) &
         - w(r3f) &
         + w(r3b) &
         - w(r4f) &
         + w(r4b) &
         - w(r6f) &
         + w(r6b) &
         - w(r7f) &
         + w(r7b) &
         - w(r11f) &
         + w(r11b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         - w(r15f) &
         + w(r15b) &
         - w(r16f) &
         + w(r16b) &
         + w(r21f) &
         + w(r24f) &
         - w(r24b) &
         + w(r25f) &
         - w(r25b) &
         + w(r26f) &
         - w(r26b) &
         - w(r28f) &
         + w(r28b) &
         + 2.0_pr * w(r30f) &
         + 2.0_pr * w(r31f) &
         + w(r34f) &
         - w(r34b) &
         + w(r36f) &
         - w(r36b) &
         + w(r37f) &
         - w(r37b) &
         - w(r42f) &
         + w(r42b) &
         - w(r43f) &
         + w(r43b) &
         - w(r51f) &
         + w(r51b) &
         - w(r55f) &
         + w(r55b) &
         - w(r57f) &
         + w(r57b) &
         - w(r60f) &
         + w(r60b) &
         + w(r63f) &
         - w(r63b) &
         - w(r66f) &
         + w(r66b) &
         + w(r69f) &
         - w(r69b) &
         - w(r70f) &
         + w(r70b) &
         + w(r71f) &
         - w(r71b) &
         - w(r73f) &
         + w(r73b) &
         - w(r75f) &
         + w(r75b) &
         - w(r78f) &
         + w(r78b) &
         + w(r80f) &
         - w(r80b) &
         - w(r81f) &
         + w(r81b) &
         - w(r83f) &
         + w(r83b) &
         + w(r88f) &
         - w(r88b) &
         - w(r89f) &
         + w(r89b) &
         + w(r96f) &
         + w(r97f) &
         - w(r102f) &
         + w(r102b) &
         + w(r106f) &
         - w(r106b) &
         + w(r109f) &
         - w(r109b) &
         - w(r110f) &
         + w(r110b) &
         + w(r121f) &
         - w(r121b) &
         + 2.0_pr * w(r123f) &
         - 2.0_pr * w(r123b) &
         + 2.0_pr * w(r126f) &
         + 2.0_pr * w(r128f) &
         + w(r129f) 

    cdot(sO2) = 0.0_pr &
         - w(r3f) &
         + w(r3b) &
         + w(r7f) &
         - w(r7b) &
         + w(r8f) &
         - w(r8b) &
         + w(r9f) &
         - w(r9b) &
         + w(r10f) &
         - w(r10b) &
         - w(r11f) &
         + w(r11b) &
         + w(r19f) &
         - w(r19b) &
         - w(r21f) &
         - w(r22f) &
         + w(r22b) &
         - w(r29f) &
         + w(r29b) &
         - w(r30f) &
         - w(r32f) &
         + w(r32b) &
         - w(r33f) &
         + w(r33b) &
         - w(r49f) &
         + w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r54f) &
         + w(r54b) &
         - w(r64f) &
         + w(r64b) &
         - w(r90f) &
         + w(r90b) &
         - w(r91f) &
         + w(r91b) &
         - w(r92f) &
         + w(r92b) &
         - w(r93f) &
         + w(r93b) &
         - w(r94f) &
         + w(r94b) &
         - w(r95f) &
         + w(r95b) &
         - w(r96f) &
         - w(r97f) &
         - w(r98f) &
         + w(r98b) &
         - w(r99f) &
         + w(r99b) &
         - w(r100f) &
         + w(r100b) &
         - w(r101f) &
         + w(r101b) &
         - w(r118f) &
         - w(r132f) &
         - w(r133f) &
         - w(r134f) &
         - w(r135f) 

    cdot(sO) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         + w(r3f) &
         - w(r3b) &
         - w(r5f) &
         + w(r5b) &
         - w(r8f) &
         + w(r8b) &
         - w(r12f) &
         + w(r12b) &
         - w(r17f) &
         + w(r17b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) &
         - w(r31f) &
         + w(r32f) &
         - w(r32b) &
         - w(r34f) &
         + w(r34b) &
         - w(r44f) &
         + w(r44b) &
         - w(r45f) &
         + w(r45b) &
         - w(r53f) &
         + w(r53b) &
         - w(r59f) &
         + w(r59b) &
         - w(r65f) &
         + w(r65b) &
         - w(r69f) &
         + w(r69b) &
         - w(r74f) &
         + w(r74b) &
         - w(r87f) &
         + w(r87b) &
         - w(r104f) &
         + w(r104b) &
         - w(r109f) &
         + w(r109b) &
         - w(r114f) &
         - w(r115f) &
         - w(r117f) &
         + w(r117b) &
         - w(r120f) &
         + w(r120b) &
         - w(r123f) &
         + w(r123b) &
         - w(r128f) &
         - w(r136f) &
         - w(r138f) &
         + w(r138b) 

    cdot(sH2O) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         + w(r4f) &
         - w(r4b) &
         - w(r5f) &
         + w(r5b) &
         + w(r9f) &
         - w(r9b) &
         + w(r10f) &
         - w(r10b) &
         + w(r18f) &
         - w(r18b) &
         + w(r22f) &
         - w(r22b) &
         + w(r35f) &
         - w(r35b) &
         + w(r38f) &
         - w(r38b) &
         + w(r41f) &
         - w(r41b) &
         + w(r46f) &
         - w(r46b) &
         + w(r47f) &
         - w(r47b) &
         + w(r52f) &
         - w(r52b) &
         + w(r61f) &
         - w(r61b) &
         + w(r67f) &
         - w(r67b) &
         + w(r76f) &
         - w(r76b) &
         + w(r84f) &
         - w(r84b) &
         + w(r103f) &
         - w(r103b) 

    cdot(sOH) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         - w(r2f) &
         + w(r2b) &
         + w(r3f) &
         - w(r3b) &
         - w(r4f) &
         + w(r4b) &
         + 2.0_pr * w(r5f) &
         - 2.0_pr * w(r5b) &
         + 2.0_pr * w(r6f) &
         - 2.0_pr * w(r6b) &
         + w(r8f) &
         - w(r8b) &
         - w(r9f) &
         + w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r13f) &
         + w(r13b) &
         - w(r14f) &
         + w(r14b) &
         + w(r17f) &
         - w(r17b) &
         - w(r18f) &
         + w(r18b) &
         + w(r21f) &
         - w(r26f) &
         + w(r26b) &
         + w(r29f) &
         - w(r29b) &
         + w(r33f) &
         - w(r33b) &
         - 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(r44f) &
         - w(r44b) &
         + w(r45f) &
         - w(r45b) &
         - w(r46f) &
         + w(r46b) &
         - w(r47f) &
         + w(r47b) &
         - w(r52f) &
         + w(r52b) &
         + w(r53f) &
         - w(r53b) &
         + w(r59f) &
         - w(r59b) &
         - w(r61f) &
         + w(r61b) &
         + w(r65f) &
         - w(r65b) &
         - w(r67f) &
         + w(r67b) &
         + w(r74f) &
         - w(r74b) &
         - w(r76f) &
         + w(r76b) &
         - w(r84f) &
         + w(r84b) &
         + w(r92f) &
         - w(r92b) &
         + w(r93f) &
         - w(r93b) &
         - w(r103f) &
         + w(r103b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         + w(r107b) &
         - w(r111f) &
         + w(r111b) &
         - w(r116f) &
         - w(r117f) &
         + w(r117b) &
         - w(r121f) &
         + w(r121b) &
         - w(r129f) &
         + w(r138f) &
         - w(r138b) 

    cdot(sCO) = 0.0_pr &
         - w(r12f) &
         + w(r12b) &
         - w(r13f) &
         + w(r13b) &
         - w(r14f) &
         + w(r14b) &
         + w(r21f) &
         + w(r22f) &
         - w(r22b) &
         + w(r23f) &
         - w(r23b) &
         + w(r27f) &
         - w(r27b) &
         + w(r31f) &
         - w(r58f) &
         + w(r58b) &
         + w(r63f) &
         - w(r63b) &
         + w(r64f) &
         - w(r64b) &
         + w(r65f) &
         - w(r65b) &
         + w(r66f) &
         - w(r66b) &
         + w(r67f) &
         - w(r67b) &
         + w(r68f) &
         - w(r68b) &
         + w(r96f) &
         + w(r97f) &
         + w(r98f) &
         - w(r98b) &
         + w(r99f) &
         - w(r99b) &
         + w(r104f) &
         - w(r104b) &
         + w(r105f) &
         - w(r105b) &
         + w(r107f) &
         - w(r107b) &
         + w(r108f) &
         - w(r108b) &
         + w(r110f) &
         - w(r110b) &
         + w(r111f) &
         - w(r111b) &
         + w(r119f) &
         - w(r119b) &
         + w(r120f) &
         - w(r120b) &
         + w(r122f) &
         - w(r122b) &
         + w(r123f) &
         - w(r123b) &
         - w(r131f) &
         + w(r131b) &
         + w(r132f) 

    cdot(sCO2) = 0.0_pr &
         + w(r12f) &
         - w(r12b) &
         + w(r13f) &
         - w(r13b) &
         + w(r14f) &
         - w(r14b) &
         - w(r27f) &
         + w(r27b) &
         + w(r30f) &
         + w(r69f) &
         - w(r69b) &
         + w(r100f) &
         - w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r116f) &
         - w(r122f) &
         + w(r122b) &
         + w(r128f) &
         + w(r129f) &
         + w(r131f) &
         - w(r131b) 

    cdot(sCH4) = 0.0_pr &
         + w(r15f) &
         - w(r15b) &
         - w(r16f) &
         + w(r16b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         + w(r19f) &
         - w(r19b) &
         - w(r20f) &
         + w(r20b) &
         + w(r48f) &
         - w(r48b) &
         + w(r56f) &
         - w(r56b) &
         + w(r62f) &
         - w(r62b) &
         + w(r68f) &
         - w(r68b) &
         + w(r77f) &
         - w(r77b) &
         + w(r79f) &
         - w(r79b) &
         + w(r85f) &
         - w(r85b) &
         + w(r86f) &
         - w(r86b) &
         - w(r124f) &
         + w(r124b) &
         + w(r137f) 

    cdot(sCH3) = 0.0_pr &
         - w(r15f) &
         + w(r15b) &
         + w(r16f) &
         - w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r18f) &
         - w(r18b) &
         - w(r19f) &
         + w(r19b) &
         + 2.0_pr * w(r20f) &
         - 2.0_pr * w(r20b) &
         + w(r25f) &
         - w(r25b) &
         + w(r28f) &
         - w(r28b) &
         - 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(r48f) &
         + w(r48b) &
         - w(r56f) &
         + w(r56b) &
         - w(r62f) &
         + w(r62b) &
         - w(r68f) &
         + w(r68b) &
         -2.0_pr * w(r72f) &
         + 2.0_pr * w(r72b) &
         - w(r77f) &
         + w(r77b) &
         - w(r79f) &
         + w(r79b) &
         -2.0_pr * w(r80f) &
         + 2.0_pr * w(r80b) &
         - w(r85f) &
         + w(r85b) &
         - w(r86f) &
         + w(r86b) &
         + w(r87f) &
         - w(r87b) &
         - w(r88f) &
         + w(r88b) &
         + w(r100f) &
         - w(r100b) &
         + w(r101f) &
         - w(r101b) &
         + w(r107f) &
         - w(r107b) &
         + w(r110f) &
         - w(r110b) &
         + w(r112f) &
         - w(r112b) &
         + w(r113f) &
         - w(r113b) &
         + w(r116f) &
         + w(r119f) &
         - w(r119b) &
         + 2.0_pr * w(r124f) &
         - 2.0_pr * w(r124b) &
         -2.0_pr * w(r125f) &
         + w(r131f) &
         - w(r131b) &
         - w(r136f) &
         - w(r137f) 

    cdot(sCH2O) = 0.0_pr &
         + w(r26f) &
         - w(r26b) &
         + w(r27f) &
         - w(r27b) &
         + w(r33f) &
         - w(r33b) &
         + w(r34f) &
         - w(r34b) &
         + 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(r61f) &
         + w(r61b) &
         - w(r62f) &
         + w(r62b) &
         - w(r70f) &
         + w(r70b) &
         + w(r71f) &
         - w(r71b) &
         + w(r94f) &
         - w(r94b) &
         + w(r95f) &
         - w(r95b) &
         + w(r96f) &
         + w(r97f) &
         + w(r118f) &
         + w(r121f) &
         - w(r121b) &
         + w(r122f) &
         - w(r122b) &
         - w(r128f) &
         - w(r129f) &
         + w(r130f) &
         + w(r132f) &
         + 2.0_pr * w(r133f) &
         + w(r137f) 

    cdot(sC2H6) = 0.0_pr &
         + w(r72f) &
         - w(r72b) &
         + w(r73f) &
         - w(r73b) &
         - w(r74f) &
         + w(r74b) &
         - w(r75f) &
         + w(r75b) &
         - w(r76f) &
         + w(r76b) &
         - w(r77f) &
         + w(r77b) &
         - w(r127f) &
         + w(r127b) 

    cdot(sC2H2) = 0.0_pr &
         + w(r82f) &
         - w(r82b) &
         - w(r89f) &
         + w(r89b) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         - w(r91b) &
         + w(r102f) &
         - w(r102b) &
         + w(r103f) &
         - w(r103b) &
         - w(r104f) &
         + w(r104b) &
         - w(r105f) &
         + w(r105b) &
         - w(r106f) &
         + w(r106b) &
         - w(r107f) &
         + w(r107b) &
         - w(r114f) &
         + w(r126f) &
         - w(r130f) &
         - w(r132f) &
         - w(r134f) 

    cdot(sCH2GL1G) = 0.0_pr &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         - 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(r31f) &
         + w(r35f) &
         - w(r35b) &
         + w(r38f) &
         - w(r38b) &
         + w(r41f) &
         - w(r41b) &
         - w(r88f) &
         + w(r88b) &
         + w(r104f) &
         - w(r104b) &
         + w(r108f) &
         - w(r108b) &
         - w(r112f) &
         + w(r112b) &
         - w(r113f) &
         + w(r113b) &
         - w(r119f) &
         + w(r119b) &
         - w(r120f) &
         + w(r120b) &
         - w(r121f) &
         + w(r121b) &
         - w(r122f) &
         + w(r122b) &
         - w(r123f) &
         + w(r123b) &
         - w(r124f) &
         + w(r124b) &
         -2.0_pr * w(r126f) 

  
    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) :: HO2_ct
    real(pr) :: HO2_num
    real(pr) :: HO2_denom
    real(pr) :: HO2_CH3O
    real(pr) :: HO2_CH2OH
    real(pr) :: HO2_HCO
    real(pr) :: HO2_C2H3
    real(pr) :: CH3OH_ct
    real(pr) :: CH3OH_num
    real(pr) :: CH3OH_denom
    real(pr) :: CH3OH_CH3O
    real(pr) :: CH3OH_CH2OH
    real(pr) :: CH3O_ct
    real(pr) :: CH3O_num
    real(pr) :: CH3O_denom
    real(pr) :: CH3O_HO2
    real(pr) :: CH3O_CH3OH
    real(pr) :: CH3O_C2H5
    real(pr) :: CH3O_C2H4
    real(pr) :: CH3O_C2H3
    real(pr) :: CH2OH_ct
    real(pr) :: CH2OH_num
    real(pr) :: CH2OH_denom
    real(pr) :: CH2OH_HO2
    real(pr) :: CH2OH_CH3OH
    real(pr) :: CH2OH_CH2CO
    real(pr) :: HCO_ct
    real(pr) :: HCO_num
    real(pr) :: HCO_denom
    real(pr) :: HCO_HO2
    real(pr) :: HCO_C2H4
    real(pr) :: HCO_C2H3
    real(pr) :: HCO_CH2CO
    real(pr) :: C2H5_ct
    real(pr) :: C2H5_num
    real(pr) :: C2H5_denom
    real(pr) :: C2H5_C2H4
    real(pr) :: C2H4_ct
    real(pr) :: C2H4_num
    real(pr) :: C2H4_denom
    real(pr) :: C2H4_HCO
    real(pr) :: C2H4_C2H5
    real(pr) :: C2H4_C2H3
    real(pr) :: C2H3_ct
    real(pr) :: C2H3_num
    real(pr) :: C2H3_denom
    real(pr) :: C2H3_HO2
    real(pr) :: C2H3_CH3O
    real(pr) :: C2H3_HCO
    real(pr) :: C2H3_C2H4
    real(pr) :: C2H3_CH2CO
    real(pr) :: CH2CO_ct
    real(pr) :: CH2CO_num
    real(pr) :: CH2CO_denom
    real(pr) :: CH2CO_CH2OH
    real(pr) :: CH2CO_C2H3
  
    HCO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r57f)* c(sH) &
             + k(r63f)* m(mM63) &
             + k(r64f)* c(sO2) &
             + k(r65f)* c(sO) &
             + k(r66f)* c(sH) &
             + k(r67f)* c(sOH) &
             + k(r68f)* c(sCH3) &
             + k(r69f)* c(sO) &
             + k(r105f)* c(sC2H2) &
             + k(r119f)* c(sCH2GL1G) &
             + k(r24b)* c(sH) &
             + k(r29b)* c(sOH) &
             + k(r59b)* c(sOH) &
             + k(r60b)* c(sH2) &
             + k(r61b)* c(sH2O) &
             + k(r62b)* c(sCH4) &
             + k(r87b)* c(sCH3) &
             + k(r94b)* c(sCH2O) &
             + k(r95b)* c(sCH2O)   )

    HCO_num = ( 0.0_pr &
             + k(r24f)* c(sCH2GL1G) * c(sO) &
             + k(r29f)* c(sCH2GL1G) * c(sO2) &
             + k(r59f)* c(sCH2O) * c(sO) &
             + k(r60f)* c(sCH2O) * c(sH) &
             + k(r61f)* c(sCH2O) * c(sOH) &
             + k(r62f)* c(sCH2O) * c(sCH3) &
     +2.0_pr * k(r134f)* c(sC2H2) * c(sO2) &
             + k(r57b)* c(sCH2O) &
             + k(r63b)* c(sCO) * c(sH) * m(mM63) &
             + k(r65b)* c(sCO) * c(sOH) &
             + k(r66b)* c(sCO) * c(sH2) &
             + k(r67b)* c(sCO) * c(sH2O) &
             + k(r68b)* c(sCH4) * c(sCO) &
             + k(r69b)* c(sCO2) * c(sH) &
             + k(r119b)* c(sCH3) * c(sCO)   )

    HCO_ct = HCO_num / HCO_denom

    HCO_HO2 = - ( 0.0_pr &
          + k(r130f) * c(sC2H2) &
          + k(r64b) * c(sCO)  ) / HCO_denom

    HCO_C2H4 = - ( 0.0_pr &
          + k(r87f) * c(sO) &
          + k(r135f) * c(sO2)  ) / HCO_denom

    HCO_C2H3 = - ( 0.0_pr &
          + k(r94f) * c(sO2) &
          + k(r95f) * c(sO2) &
          + k(r105b) * c(sCO)  ) / HCO_denom

    HCO_CH2CO = - ( 0.0_pr &
          + k(r115f)  * 2.0_pr * c(sO)  ) / HCO_denom

    CH3O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r54f)* c(sO2) &
             + k(r55f)* c(sH) &
             + k(r56f)* c(sCH3) &
             + k(r71f)&
             + k(r131f)* c(sCO) &
             + k(r32b)* c(sO) &
             + k(r37b)* c(sH) &
             + k(r39b)* c(sOH) &
             + k(r42b)* c(sH2) &
             + k(r44b)* c(sOH) &
             + k(r46b)* c(sH2O) &
             + k(r48b)* c(sCH4) &
             + k(r98b)* c(sCO) &
             + k(r99b)* c(sCO)   )

    CH3O_num = ( 0.0_pr &
             + k(r32f)* c(sCH3) * c(sO2) &
             + k(r37f)* c(sCH3) * c(sOH) &
             + k(r136f)* c(sCH3) * c(sO) * m(mM136) &
             + k(r55b)* c(sCH2O) * c(sH2) &
             + k(r56b)* c(sCH2O) * c(sCH4) &
             + k(r71b)* c(sCH2O) * c(sH) &
             + k(r131b)* c(sCH3) * c(sCO2)   )

    CH3O_ct = CH3O_num / CH3O_denom

    CH3O_HO2 = - ( 0.0_pr &
          + k(r39f) * c(sCH3) &
          + k(r54b) * c(sCH2O)  ) / CH3O_denom

    CH3O_CH3OH = - ( 0.0_pr &
          + k(r42f) * c(sH) &
          + k(r44f) * c(sO) &
          + k(r46f) * c(sOH) &
          + k(r48f) * c(sCH3)  ) / CH3O_denom

    CH3O_C2H5 = - ( 0.0_pr &
          + k(r118f) * c(sO2)  ) / CH3O_denom

    CH3O_C2H4 = - ( 0.0_pr &
          + k(r135f) * c(sO2)  ) / CH3O_denom

    CH3O_C2H3 = - ( 0.0_pr &
          + k(r98f) * c(sO2) &
          + k(r99f) * c(sO2)  ) / CH3O_denom

    C2H4_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r78f)* c(sH) &
             + k(r82f)* m(mM82) &
             + k(r83f)* c(sH) &
             + k(r84f)* c(sOH) &
             + k(r85f)* c(sCH3) &
             + k(r86f)* c(sCH3) &
             + k(r87f)* c(sO) &
             + k(r112f)* c(sCH2GL1G) &
             + k(r113f)* c(sCH2GL1G) &
             + k(r133f)* c(sO2) &
             + k(r135f)* c(sO2) &
             + k(r138f)* c(sO) &
             + k(r79b)* c(sCH4) &
             + k(r81b)&
             + k(r88b)* c(sH) &
             + k(r127b)* c(sH2)   )

    C2H4_num = ( 0.0_pr &
             + k(r88f)* c(sCH2GL1G) * c(sCH3) &
             + k(r125f) *c(sCH3)** 2.0_pr&
             + k(r127f)* c(sC2H6) &
             + k(r82b)* c(sC2H2) * c(sH2) * m(mM82)   )

    C2H4_ct = C2H4_num / C2H4_denom

    C2H4_HCO = - ( 0.0_pr &
          + k(r87b) * c(sCH3)  ) / C2H4_denom

    C2H4_C2H5 = - ( 0.0_pr &
          + k(r79f) * c(sCH3) &
          + k(r78b)  ) / C2H4_denom

    C2H4_C2H3 = - ( 0.0_pr &
          + k(r81f) * c(sH) &
          + k(r83b) * c(sH2) &
          + k(r84b) * c(sH2O) &
          + k(r85b) * c(sCH4) &
          + k(r86b) * c(sCH4) &
          + k(r112b) * c(sCH3) &
          + k(r113b) * c(sCH3) &
          + k(r138b) * c(sOH)  ) / C2H4_denom

    CH3OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r40f)&
             + k(r41f)&
             + k(r42f)* c(sH) &
             + k(r43f)* c(sH) &
             + k(r44f)* c(sO) &
             + k(r45f)* c(sO) &
             + k(r46f)* c(sOH) &
             + k(r47f)* c(sOH) &
             + k(r48f)* c(sCH3)   )

    CH3OH_num = ( 0.0_pr &
             + k(r40b)* c(sCH3) * c(sOH) &
             + k(r41b)* c(sCH2GL1G) * c(sH2O)   )

    CH3OH_ct = CH3OH_num / CH3OH_denom

    CH3OH_CH3O = - ( 0.0_pr &
          + k(r42b) * c(sH2) &
          + k(r44b) * c(sOH) &
          + k(r46b) * c(sH2O) &
          + k(r48b) * c(sCH4)  ) / CH3OH_denom

    CH3OH_CH2OH = - ( 0.0_pr &
          + k(r43b) * c(sH2) &
          + k(r45b) * c(sOH) &
          + k(r47b) * c(sH2O)  ) / CH3OH_denom

    C2H3_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r81f)* c(sH) &
             + k(r90f)* c(sO2) &
             + k(r91f)* c(sO2) &
             + k(r92f)* c(sO2) &
             + k(r93f)* c(sO2) &
             + k(r94f)* c(sO2) &
             + k(r95f)* c(sO2) &
             + k(r96f)* c(sO2) &
             + k(r97f)* c(sO2) &
             + k(r98f)* c(sO2) &
             + k(r99f)* c(sO2) &
             + k(r100f)* c(sO2) &
             + k(r101f)* c(sO2) &
             + k(r102f)* c(sH) &
             + k(r103f)* c(sOH) &
             + k(r109f)* c(sO) &
             + k(r83b)* c(sH2) &
             + k(r84b)* c(sH2O) &
             + k(r85b)* c(sCH4) &
             + k(r86b)* c(sCH4) &
             + k(r89b)&
             + k(r105b)* c(sCO) &
             + k(r112b)* c(sCH3) &
             + k(r113b)* c(sCH3) &
             + k(r138b)* c(sOH)   )

    C2H3_num = ( 0.0_pr &
             + k(r89f)* c(sC2H2) * c(sH) &
             + k(r100b)* c(sCH3) * c(sCO2) &
             + k(r101b)* c(sCH3) * c(sCO2) &
             + k(r102b)* c(sC2H2) * c(sH2) &
             + k(r103b)* c(sC2H2) * c(sH2O)   )

    C2H3_ct = C2H3_num / C2H3_denom

    C2H3_HO2 = - ( 0.0_pr &
          + k(r90b) * c(sC2H2) &
          + k(r91b) * c(sC2H2)  ) / C2H3_denom

    C2H3_CH3O = - ( 0.0_pr &
          + k(r98b) * c(sCO) &
          + k(r99b) * c(sCO)  ) / C2H3_denom

    C2H3_HCO = - ( 0.0_pr &
          + k(r105f) * c(sC2H2) &
          + k(r94b) * c(sCH2O) &
          + k(r95b) * c(sCH2O)  ) / C2H3_denom

    C2H3_C2H4 = - ( 0.0_pr &
          + k(r83f) * c(sH) &
          + k(r84f) * c(sOH) &
          + k(r85f) * c(sCH3) &
          + k(r86f) * c(sCH3) &
          + k(r112f) * c(sCH2GL1G) &
          + k(r113f) * c(sCH2GL1G) &
          + k(r138f) * c(sO) &
          + k(r81b)  ) / C2H3_denom

    C2H3_CH2CO = - ( 0.0_pr &
          + k(r92b) * c(sOH) &
          + k(r93b) * c(sOH) &
          + k(r109b) * c(sH)  ) / C2H3_denom

    C2H5_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r73f)* c(sH) &
             + k(r79f)* c(sCH3) &
             + k(r118f)* c(sO2) &
             + k(r74b)* c(sOH) &
             + k(r75b)* c(sH2) &
             + k(r76b)* c(sH2O) &
             + k(r77b)* c(sCH4) &
             + k(r78b)&
             + k(r80b)* c(sH)   )

    C2H5_num = ( 0.0_pr &
             + k(r74f)* c(sC2H6) * c(sO) &
             + k(r75f)* c(sC2H6) * c(sH) &
             + k(r76f)* c(sC2H6) * c(sOH) &
             + k(r77f)* c(sC2H6) * c(sCH3) &
             + k(r80f) *c(sCH3)** 2.0_pr&
             + k(r73b)* c(sC2H6)   )

    C2H5_ct = C2H5_num / C2H5_denom

    C2H5_C2H4 = - ( 0.0_pr &
          + k(r78f) * c(sH) &
          + k(r79b) * c(sCH4)  ) / C2H5_denom

    CH2CO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r108f)&
             + k(r110f)* c(sH) &
             + k(r111f)* c(sOH) &
             + k(r115f)* c(sO) &
             + k(r116f)* c(sOH) &
             + k(r92b)* c(sOH) &
             + k(r93b)* c(sOH) &
             + k(r106b)* c(sH) &
             + k(r109b)* c(sH)   )

    CH2CO_num = ( 0.0_pr &
             + k(r106f)* c(sC2H2) * c(sOH) &
             + k(r114f)* c(sC2H2) * c(sO) &
             + k(r108b)* c(sCH2GL1G) * c(sCO) &
             + k(r110b)* c(sCH3) * c(sCO)   )

    CH2CO_ct = CH2CO_num / CH2CO_denom

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

    CH2CO_C2H3 = - ( 0.0_pr &
          + k(r92f) * c(sO2) &
          + k(r93f) * c(sO2) &
          + k(r109f) * c(sO)  ) / CH2CO_denom

    CH2OH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r49f)* c(sO2) &
             + k(r50f)* c(sO2) &
             + k(r51f)* c(sH) &
             + k(r52f)* c(sOH) &
             + k(r53f)* c(sO) &
             + k(r137f)* c(sCH3) &
             + k(r36b)* c(sH) &
             + k(r43b)* c(sH2) &
             + k(r45b)* c(sOH) &
             + k(r47b)* c(sH2O) &
             + k(r70b)&
             + k(r111b)* c(sCO)   )

    CH2OH_num = ( 0.0_pr &
             + k(r36f)* c(sCH3) * c(sOH) &
             + k(r70f)* c(sCH2O) * c(sH) &
             + k(r51b)* c(sCH2O) * c(sH2) &
             + k(r52b)* c(sCH2O) * c(sH2O) &
             + k(r53b)* c(sCH2O) * c(sOH)   )

    CH2OH_ct = CH2OH_num / CH2OH_denom

    CH2OH_HO2 = - ( 0.0_pr &
          + k(r49b) * c(sCH2O) &
          + k(r50b) * c(sCH2O)  ) / CH2OH_denom

    CH2OH_CH3OH = - ( 0.0_pr &
          + k(r43f) * c(sH) &
          + k(r45f) * c(sO) &
          + k(r47f) * c(sOH)  ) / CH2OH_denom

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

    HO2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r6f)* c(sH) &
             + k(r7f)* c(sH) &
             + k(r8f)* c(sO) &
             + k(r9f)* c(sOH) &
             + k(r10f)* c(sOH) &
             + k(r19f)* c(sCH3) &
             + k(r39f)* c(sCH3) &
             + k(r130f)* c(sC2H2) &
             + k(r11b)&
             + k(r49b)* c(sCH2O) &
             + k(r50b)* c(sCH2O) &
             + k(r54b)* c(sCH2O) &
             + k(r64b)* c(sCO) &
             + k(r90b)* c(sC2H2) &
             + k(r91b)* c(sC2H2) &
             + k(r117b)* m(mM117)   )

    HO2_num = ( 0.0_pr &
             + k(r11f)* c(sH) * c(sO2) &
             + k(r117f)* c(sO) * c(sOH) * m(mM117) &
             + k(r6b) *c(sOH)** 2.0_pr&
             + k(r7b)* c(sH2) * c(sO2) &
             + k(r8b)* c(sO2) * c(sOH) &
             + k(r9b)* c(sH2O) * c(sO2) &
             + k(r10b)* c(sH2O) * c(sO2) &
             + k(r19b)* c(sCH4) * c(sO2)   )

    HO2_ct = HO2_num / HO2_denom

    HO2_CH3O = - ( 0.0_pr &
          + k(r54f) * c(sO2) &
          + k(r39b) * c(sOH)  ) / HO2_denom

    HO2_CH2OH = - ( 0.0_pr &
          + k(r49f) * c(sO2) &
          + k(r50f) * c(sO2)  ) / HO2_denom

    HO2_HCO = - ( 0.0_pr &
          + k(r64f) * c(sO2)  ) / HO2_denom

    HO2_C2H3 = - ( 0.0_pr &
          + k(r90f) * c(sO2) &
          + k(r91f) * c(sO2)  ) / HO2_denom

    cqss(sqssHCO) = (  ( ( ( ( ( ( (HCO_ct) - (HO2_ct) * (HCO_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (HCO_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (HCO_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (HCO_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( ( (HCO_C2H3) &
                 - (HO2_C2H3) * (HCO_HO2)) - (CH2CO_C2H3) * ( (HCO_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( ( ( (CH3OH_ct) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) - ( (CH2CO_ct) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( - (CH2OH_CH3OH) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - ( ( (C2H4_ct) - (C2H5_ct) * (C2H4_C2H5)) &
                 - ( ( ( (C2H3_ct) - (HO2_ct) * (C2H3_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * (C2H4_C2H3) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( (HCO_C2H4) - (C2H3_C2H4) * ( ( (HCO_C2H3) &
                 - (HO2_C2H3) * (HCO_HO2)) - (CH2CO_C2H3) * ( (HCO_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))))) &
                 - ( ( ( ( ( ( ( (CH3O_ct) - (HO2_ct) * (CH3O_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - (C2H5_ct) * (CH3O_C2H5)) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( ( (CH3O_C2H3) &
                 - (HO2_C2H3) * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( ( ( (CH3OH_ct) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) - ( (CH2CO_ct) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( (CH3O_CH3OH) - (CH2OH_CH3OH) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - ( ( (C2H4_ct) - (C2H5_ct) * (C2H4_C2H5)) &
                 - ( ( ( (C2H3_ct) - (HO2_ct) * (C2H3_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * (C2H4_C2H3) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( ( (CH3O_C2H4) - (C2H5_C2H4) * (CH3O_C2H5)) &
                 - (C2H3_C2H4) * ( ( (CH3O_C2H3) - (HO2_C2H3) &
                 * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))))) &
                 * ( ( ( - (HO2_CH3O) * (HCO_HO2)) - ( (C2H3_CH3O) &
                 - (HO2_CH3O) * (C2H3_HO2)) * ( ( (HCO_C2H3) &
                 - (HO2_C2H3) * (HCO_HO2)) - (CH2CO_C2H3) * ( (HCO_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( (CH3OH_CH3O) - ( (C2H3_CH3O) - (HO2_CH3O) &
                 * (C2H3_HO2)) * ( - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * (CH3OH_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) / ( (1.0_pr) - (CH2OH_CH2CO) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2)))) / ( ( (1.0_pr) - (HO2_C2H3) &
                 * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( - (CH2OH_CH3OH) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( ( (1.0_pr) - (HO2_CH3O) * (CH3O_HO2)) &
                 - ( (C2H3_CH3O) - (HO2_CH3O) * (C2H3_HO2)) &
                 * ( ( (CH3O_C2H3) - (HO2_C2H3) * (CH3O_HO2)) &
                 - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( (CH3OH_CH3O) - ( (C2H3_CH3O) - (HO2_CH3O) &
                 * (C2H3_HO2)) * ( - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * (CH3OH_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) / ( (1.0_pr) - (CH2OH_CH2CO) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2)))) / ( ( (1.0_pr) - (HO2_C2H3) &
                 * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( (CH3O_CH3OH) - (CH2OH_CH3OH) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 ) / (  ( ( (1.0_pr) - (HO2_HCO) * (HCO_HO2)) &
                 - ( (C2H3_HCO) - (HO2_HCO) * (C2H3_HO2)) * ( ( (HCO_C2H3) &
                 - (HO2_C2H3) * (HCO_HO2)) - (CH2CO_C2H3) * ( (HCO_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( (C2H4_HCO) - ( (C2H3_HCO) - (HO2_HCO) * (C2H3_HO2)) &
                 * (C2H4_C2H3) / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) &
                 - (CH2CO_C2H3) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))))) * ( (HCO_C2H4) &
                 - (C2H3_C2H4) * ( ( (HCO_C2H3) - (HO2_C2H3) &
                 * (HCO_HO2)) - (CH2CO_C2H3) * ( (HCO_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (HCO_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 )  

    cqss(sqssCH3O) = ( ( ( ( ( ( ( (CH3O_ct) - (HO2_ct) * (CH3O_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - (C2H5_ct) * (CH3O_C2H5)) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( ( (CH3O_C2H3) &
                 - (HO2_C2H3) * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( ( ( (CH3OH_ct) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) - ( (CH2CO_ct) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( (CH3O_CH3OH) - (CH2OH_CH3OH) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - ( ( (C2H4_ct) - (C2H5_ct) * (C2H4_C2H5)) &
                 - ( ( ( (C2H3_ct) - (HO2_ct) * (C2H3_HO2)) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * (C2H4_C2H3) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( ( (CH3O_C2H4) - (C2H5_C2H4) * (CH3O_C2H5)) &
                 - (C2H3_C2H4) * ( ( (CH3O_C2H3) - (HO2_C2H3) &
                 * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( ( ( - (HO2_HCO) * (CH3O_HO2)) - ( (C2H3_HCO) &
                 - (HO2_HCO) * (C2H3_HO2)) * ( ( (CH3O_C2H3) &
                 - (HO2_C2H3) * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( (C2H4_HCO) - ( (C2H3_HCO) - (HO2_HCO) * (C2H3_HO2)) &
                 * (C2H4_C2H3) / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) &
                 - (CH2CO_C2H3) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))))) * ( ( (CH3O_C2H4) &
                 - (C2H5_C2H4) * (CH3O_C2H5)) - (C2H3_C2H4) &
                 * ( ( (CH3O_C2H3) - (HO2_C2H3) * (CH3O_HO2)) &
                 - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))))) &
                 * cqss(sqssHCO)) / ( ( ( (1.0_pr) - (HO2_CH3O) &
                 * (CH3O_HO2)) - ( (C2H3_CH3O) - (HO2_CH3O) &
                 * (C2H3_HO2)) * ( ( (CH3O_C2H3) - (HO2_C2H3) &
                 * (CH3O_HO2)) - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (CH3O_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 - ( (CH3OH_CH3O) - ( (C2H3_CH3O) - (HO2_CH3O) &
                 * (C2H3_HO2)) * ( - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) &
                 * (CH3OH_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) / ( (1.0_pr) - (CH2OH_CH2CO) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2)))) / ( ( (1.0_pr) - (HO2_C2H3) &
                 * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * ( (CH3O_CH3OH) - (CH2OH_CH3OH) * ( - (HO2_CH2OH) &
                 * (CH3O_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH3OH) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))  

    cqss(sqssC2H4) = ( ( (C2H4_ct) - (C2H5_ct) * (C2H4_C2H5)) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * (C2H4_C2H3) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - ( (C2H4_HCO) - ( (C2H3_HCO) - (HO2_HCO) * (C2H3_HO2)) &
                 * (C2H4_C2H3) / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) &
                 - (CH2CO_C2H3) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))))) * cqss(sqssHCO) &
                 - ( - ( (C2H3_CH3O) - (HO2_CH3O) * (C2H3_HO2)) &
                 * (C2H4_C2H3) / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) &
                 - (CH2CO_C2H3) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))))) * cqss(sqssCH3O)) &
                 / ( ( (1.0_pr) - (C2H5_C2H4) * (C2H4_C2H5)) &
                 - (C2H3_C2H4) * (C2H4_C2H3) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))))  

    cqss(sqssCH3OH) = ( ( ( (CH3OH_ct) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) - ( (CH2CO_ct) &
                 - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) - ( ( ( (C2H3_ct) &
                 - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) - (HO2_ct) &
                 * (CH2OH_HO2)) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 - ( - ( (C2H3_HCO) - (HO2_HCO) * (C2H3_HO2)) &
                 * ( - (CH2CO_C2H3) * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * cqss(sqssHCO) - ( (CH3OH_CH3O) - ( (C2H3_CH3O) &
                 - (HO2_CH3O) * (C2H3_HO2)) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * cqss(sqssCH3O) - ( - (C2H3_C2H4) * ( - (CH2CO_C2H3) &
                 * ( - (CH2OH_CH2CO) * (CH3OH_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2)))) / ( ( (1.0_pr) &
                 - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) * ( (C2H3_CH2CO) &
                 - (CH2OH_CH2CO) * ( - (HO2_CH2OH) * (C2H3_HO2)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))) &
                 * cqss(sqssC2H4)) / ( (1.0_pr) - (CH2OH_CH3OH) &
                 * (CH3OH_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2)))  

    cqss(sqssC2H3) = ( ( ( (C2H3_ct) - (HO2_ct) * (C2H3_HO2)) - ( (CH2OH_ct) &
                 - (HO2_ct) * (CH2OH_HO2)) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 - ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2))) * ( (C2H3_CH2CO) - (CH2OH_CH2CO) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) / ( (1.0_pr) &
                 - (CH2OH_CH2CO) * (CH2CO_CH2OH) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) - ( (C2H3_HCO) &
                 - (HO2_HCO) * (C2H3_HO2)) * cqss(sqssHCO) - ( (C2H3_CH3O) &
                 - (HO2_CH3O) * (C2H3_HO2)) * cqss(sqssCH3O) &
                 - (C2H3_C2H4) * cqss(sqssC2H4) - ( - (CH2OH_CH3OH) &
                 * ( - (HO2_CH2OH) * (C2H3_HO2)) / ( (1.0_pr) &
                 - (HO2_CH2OH) * (CH2OH_HO2))) * cqss(sqssCH3OH)) &
                 / ( ( (1.0_pr) - (HO2_C2H3) * (C2H3_HO2)) - (CH2CO_C2H3) &
                 * ( (C2H3_CH2CO) - (CH2OH_CH2CO) * ( - (HO2_CH2OH) &
                 * (C2H3_HO2)) / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))))  

    cqss(sqssC2H5) = C2H5_ct - (C2H5_C2H4) * cqss(sqssC2H4)  

    cqss(sqssCH2CO) = ( (CH2CO_ct) - ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2)) &
                 * (CH2CO_CH2OH) / ( (1.0_pr) - (HO2_CH2OH) &
                 * (CH2OH_HO2)) - ( - (CH2OH_CH3OH) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))) &
                 * cqss(sqssCH3OH) - (CH2CO_C2H3) * cqss(sqssC2H3)) &
                 / ( (1.0_pr) - (CH2OH_CH2CO) * (CH2CO_CH2OH) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2)))  

    cqss(sqssCH2OH) = ( (CH2OH_ct) - (HO2_ct) * (CH2OH_HO2) - ( - (HO2_HCO) &
                 * (CH2OH_HO2)) * cqss(sqssHCO) - ( - (HO2_CH3O) &
                 * (CH2OH_HO2)) * cqss(sqssCH3O) - (CH2OH_CH3OH) &
                 * cqss(sqssCH3OH) - ( - (HO2_C2H3) * (CH2OH_HO2)) &
                 * cqss(sqssC2H3) - (CH2OH_CH2CO) * cqss(sqssCH2CO)) &
                 / ( (1.0_pr) - (HO2_CH2OH) * (CH2OH_HO2))  

    cqss(sqssHO2) = HO2_ct - (HO2_HCO) * cqss(sqssHCO) - (HO2_CH3O) &
                 * cqss(sqssCH3O) - (HO2_C2H3) * cqss(sqssC2H3) &
                 - (HO2_CH2OH) * cqss(sqssCH2OH)  

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

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

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

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

    ! Conversion
    c = y * rho / W_sp

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

        return
    end subroutine y2c

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

  use mod_customkinetics
  implicit none

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

  real(pr) :: P, T, rho

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

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

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

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

  ! Evaluate production rates
  call get_production_rates(wdot,w)

  return
end subroutine customkinetics