!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE S33R182QSS14.f90
!>    @file S33R182QSS14.f90 
!!    Module for calculating the analytical source terms in Cantera  
!!    @details 
!!    @authors 
!!    @date    2020/10/06   
!!    @since     
!!    @note      
!--------------------------------------------------------------------------------------------------

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

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

  ! Ideal gas constant
  real(pr), parameter :: Rcst = 8.3144621_pr
  
  ! Use semi-implicit formulation for select species
  logical, parameter :: semiImplicit = .False.
  
  ! Number of elements in the chemical system 
  integer, parameter :: ne = 3
  
  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 33
  integer, parameter :: nqss = 0
  integer, parameter :: nreac = 182
  integer, parameter :: nreac_reverse = 47
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 0
  integer, parameter :: nFO = 6
  integer, parameter :: nTB_reverse = 0
  integer, parameter :: nFO_reverse = 6
  
  ! Index of elements
  integer, parameter :: eH = 1
  integer, parameter :: eC = 2
  integer, parameter :: eO = 3

  ! Index of species
  integer, parameter :: sH2O = 1
  integer, parameter :: sH2 = 2
  integer, parameter :: sCH4 = 3
  integer, parameter :: sC2H2 = 4
  integer, parameter :: sC2H4 = 5
  integer, parameter :: sC2H6 = 6
  integer, parameter :: sAC3H4 = 7
  integer, parameter :: sC3H6 = 8
  integer, parameter :: sC3H8 = 9
  integer, parameter :: sC4H6 = 10
  integer, parameter :: sNC4H8 = 11
  integer, parameter :: sNC4H10 = 12
  integer, parameter :: sCYC5H6 = 13
  integer, parameter :: sC5H8 = 14
  integer, parameter :: sNC5H10 = 15
  integer, parameter :: sC6H6 = 16
  integer, parameter :: sMCPTD = 17
  integer, parameter :: sH = 18
  integer, parameter :: sCH3 = 19
  integer, parameter :: sC2H3 = 20
  integer, parameter :: sC2H5 = 21
  integer, parameter :: sC3H3 = 22
  integer, parameter :: sCHCHCH3 = 23
  integer, parameter :: sCH2CHCH2 = 24
  integer, parameter :: sNC3H7 = 25
  integer, parameter :: sIC3H7 = 26
  integer, parameter :: sC4H5 = 27
  integer, parameter :: sCH2C3H5 = 28
  integer, parameter :: sSC4H7 = 29
  integer, parameter :: sNC4H9S = 30
  integer, parameter :: sNC4H9P = 31
  integer, parameter :: sCYC5H5 = 32
  integer, parameter :: sNC5H9X3 = 33

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r39f = 39
  integer, parameter :: r40f = 40
  integer, parameter :: r41f = 41
  integer, parameter :: r42f = 42
  integer, parameter :: r43f = 43
  integer, parameter :: r44f = 44
  integer, parameter :: r45f = 45
  integer, parameter :: r46f = 46
  integer, parameter :: r47f = 47
  integer, parameter :: r48f = 48
  integer, parameter :: r49f = 49
  integer, parameter :: r50f = 50
  integer, parameter :: r51f = 51
  integer, parameter :: r52f = 52
  integer, parameter :: r53f = 53
  integer, parameter :: r54f = 54
  integer, parameter :: r55f = 55
  integer, parameter :: r56f = 56
  integer, parameter :: r57f = 57
  integer, parameter :: r58f = 58
  integer, parameter :: r59f = 59
  integer, parameter :: r60f = 60
  integer, parameter :: r61f = 61
  integer, parameter :: r62f = 62
  integer, parameter :: r63f = 63
  integer, parameter :: r64f = 64
  integer, parameter :: r65f = 65
  integer, parameter :: r66f = 66
  integer, parameter :: r67f = 67
  integer, parameter :: r68f = 68
  integer, parameter :: r69f = 69
  integer, parameter :: r70f = 70
  integer, parameter :: r71f = 71
  integer, parameter :: r72f = 72
  integer, parameter :: r73f = 73
  integer, parameter :: r74f = 74
  integer, parameter :: r75f = 75
  integer, parameter :: r76f = 76
  integer, parameter :: r77f = 77
  integer, parameter :: r78f = 78
  integer, parameter :: r79f = 79
  integer, parameter :: r80f = 80
  integer, parameter :: r81f = 81
  integer, parameter :: r82f = 82
  integer, parameter :: r83f = 83
  integer, parameter :: r84f = 84
  integer, parameter :: r85f = 85
  integer, parameter :: r86f = 86
  integer, parameter :: r87f = 87
  integer, parameter :: r88f = 88
  integer, parameter :: r89f = 89
  integer, parameter :: r90f = 90
  integer, parameter :: r91f = 91
  integer, parameter :: r92f = 92
  integer, parameter :: r93f = 93
  integer, parameter :: r94f = 94
  integer, parameter :: r95f = 95
  integer, parameter :: r96f = 96
  integer, parameter :: r97f = 97
  integer, parameter :: r98f = 98
  integer, parameter :: r99f = 99
  integer, parameter :: r100f = 100
  integer, parameter :: r101f = 101
  integer, parameter :: r102f = 102
  integer, parameter :: r103f = 103
  integer, parameter :: r104f = 104
  integer, parameter :: r105f = 105
  integer, parameter :: r106f = 106
  integer, parameter :: r107f = 107
  integer, parameter :: r108f = 108
  integer, parameter :: r109f = 109
  integer, parameter :: r110f = 110
  integer, parameter :: r111f = 111
  integer, parameter :: r112f = 112
  integer, parameter :: r113f = 113
  integer, parameter :: r114f = 114
  integer, parameter :: r115f = 115
  integer, parameter :: r116f = 116
  integer, parameter :: r117f = 117
  integer, parameter :: r118f = 118
  integer, parameter :: r119f = 119
  integer, parameter :: r120f = 120
  integer, parameter :: r121f = 121
  integer, parameter :: r122f = 122
  integer, parameter :: r123f = 123
  integer, parameter :: r124f = 124
  integer, parameter :: r125f = 125
  integer, parameter :: r126f = 126
  integer, parameter :: r127f = 127
  integer, parameter :: r128f = 128
  integer, parameter :: r129f = 129
  integer, parameter :: r130f = 130
  integer, parameter :: r131f = 131
  integer, parameter :: r132f = 132
  integer, parameter :: r133f = 133
  integer, parameter :: r134f = 134
  integer, parameter :: r135f = 135
  integer, parameter :: r136f = 136
  integer, parameter :: r137f = 137
  integer, parameter :: r138f = 138
  integer, parameter :: r139f = 139
  integer, parameter :: r140f = 140
  integer, parameter :: r141f = 141
  integer, parameter :: r142f = 142
  integer, parameter :: r143f = 143
  integer, parameter :: r144f = 144
  integer, parameter :: r145f = 145
  integer, parameter :: r146f = 146
  integer, parameter :: r147f = 147
  integer, parameter :: r148f = 148
  integer, parameter :: r149f = 149
  integer, parameter :: r150f = 150
  integer, parameter :: r151f = 151
  integer, parameter :: r152f = 152
  integer, parameter :: r153f = 153
  integer, parameter :: r154f = 154
  integer, parameter :: r155f = 155
  integer, parameter :: r156f = 156
  integer, parameter :: r157f = 157
  integer, parameter :: r158f = 158
  integer, parameter :: r159f = 159
  integer, parameter :: r160f = 160
  integer, parameter :: r161f = 161
  integer, parameter :: r162f = 162
  integer, parameter :: r163f = 163
  integer, parameter :: r164f = 164
  integer, parameter :: r165f = 165
  integer, parameter :: r166f = 166
  integer, parameter :: r167f = 167
  integer, parameter :: r168f = 168
  integer, parameter :: r169f = 169
  integer, parameter :: r170f = 170
  integer, parameter :: r171f = 171
  integer, parameter :: r172f = 172
  integer, parameter :: r173f = 173
  integer, parameter :: r174f = 174
  integer, parameter :: r175f = 175
  integer, parameter :: r176f = 176
  integer, parameter :: r177f = 177
  integer, parameter :: r178f = 178
  integer, parameter :: r179f = 179
  integer, parameter :: r180f = 180
  integer, parameter :: r181f = 181
  integer, parameter :: r182f = 182
  integer, parameter :: r1b = 183
  integer, parameter :: r2b = 184
  integer, parameter :: r3b = 185
  integer, parameter :: r4b = 186
  integer, parameter :: r5b = 187
  integer, parameter :: r6b = 188
  integer, parameter :: r7b = 189
  integer, parameter :: r8b = 190
  integer, parameter :: r10b = 191
  integer, parameter :: r11b = 192
  integer, parameter :: r16b = 193
  integer, parameter :: r17b = 194
  integer, parameter :: r18b = 195
  integer, parameter :: r19b = 196
  integer, parameter :: r20b = 197
  integer, parameter :: r21b = 198
  integer, parameter :: r22b = 199
  integer, parameter :: r23b = 200
  integer, parameter :: r24b = 201
  integer, parameter :: r25b = 202
  integer, parameter :: r26b = 203
  integer, parameter :: r27b = 204
  integer, parameter :: r28b = 205
  integer, parameter :: r42b = 206
  integer, parameter :: r43b = 207
  integer, parameter :: r44b = 208
  integer, parameter :: r45b = 209
  integer, parameter :: r47b = 210
  integer, parameter :: r48b = 211
  integer, parameter :: r49b = 212
  integer, parameter :: r50b = 213
  integer, parameter :: r51b = 214
  integer, parameter :: r52b = 215
  integer, parameter :: r53b = 216
  integer, parameter :: r54b = 217
  integer, parameter :: r55b = 218
  integer, parameter :: r79b = 219
  integer, parameter :: r80b = 220
  integer, parameter :: r81b = 221
  integer, parameter :: r83b = 222
  integer, parameter :: r84b = 223
  integer, parameter :: r85b = 224
  integer, parameter :: r86b = 225
  integer, parameter :: r87b = 226
  integer, parameter :: r90b = 227
  integer, parameter :: r91b = 228
  integer, parameter :: r92b = 229
  
  ! Index of third body species
  integer, parameter :: mM1 = 1
  integer, parameter :: mM3 = 2
  integer, parameter :: mM4 = 3
  integer, parameter :: mM16 = 4
  integer, parameter :: mM17 = 5
  integer, parameter :: mM18 = 6
  
  ! Index of fall off reactions
  integer, parameter :: FOr1f = 1
  integer, parameter :: FOr3f = 2
  integer, parameter :: FOr4f = 3
  integer, parameter :: FOr16f = 4
  integer, parameter :: FOr17f = 5
  integer, parameter :: FOr18f = 6
  integer, parameter :: FOr1b = 7
  integer, parameter :: FOr3b = 8
  integer, parameter :: FOr4b = 9
  integer, parameter :: FOr16b = 10
  integer, parameter :: FOr17b = 11
  integer, parameter :: FOr18b = 12
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       0.01801528_pr, & ! H2O
       0.00201588_pr, & ! H2
       0.01604276_pr, & ! CH4
       0.02603788_pr, & ! C2H2
       0.02805376_pr, & ! C2H4
       0.03006964_pr, & ! C2H6
       0.04006476_pr, & ! AC3H4
       0.04208064_pr, & ! C3H6
       0.04409652_pr, & ! C3H8
       0.05409164_pr, & ! C4H6
       0.05610752_pr, & ! NC4H8
       0.0581234_pr, & ! NC4H10
       0.06610264_pr, & ! CYC5H6
       0.06811852_pr, & ! C5H8
       0.0701344_pr, & ! NC5H10
       0.07811364_pr, & ! C6H6
       0.08012952_pr, & ! MCPTD
       0.00100794_pr, & ! H
       0.01503482_pr, & ! CH3
       0.02704582_pr, & ! C2H3
       0.0290617_pr, & ! C2H5
       0.03905682_pr, & ! C3H3
       0.0410727_pr, & ! CHCHCH3
       0.0410727_pr, & ! CH2CHCH2
       0.04308858_pr, & ! NC3H7
       0.04308858_pr, & ! IC3H7
       0.0530837_pr, & ! C4H5
       0.05509958_pr, & ! CH2C3H5
       0.05509958_pr, & ! SC4H7
       0.05711546_pr, & ! NC4H9S
       0.05711546_pr, & ! NC4H9P
       0.0650947_pr, & ! CYC5H5
       0.06912646_pr & ! NC5H9X3
       /)
  
contains

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

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

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

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

  end function getlindratecoeff

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

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

    M(mM1) = (1_pr)*c(sH2) &
         + sum(c)

    M(mM3) = (2_pr)*c(sC2H6) &
         + (1_pr)*c(sCH4) &
         + (1_pr)*c(sH2) &
         + sum(c)

    M(mM4) = sum(c)

    M(mM16) = (1_pr)*c(sH2) &
         + sum(c)

    M(mM17) = (1_pr)*c(sH2) &
         + sum(c)

    M(mM18) = 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_0(FOr1f) = (2.33e+22_pr)*exp((5020.8_pr)*R_T_inv + T_log *(-5.03_pr))
    k_inf(FOr1f) = (25000000_pr)
    FC(FOr1f) = ((1.0_pr - 0.38_pr)*exp(-Tloc/(73_pr))) + (0.38_pr*exp(-Tloc/(1180_pr)))
    k(r1f) = getlindratecoeff(Tloc,k_0(FOr1f),k_inf(FOr1f),FC(FOr1f), M(mM1),Ploc)
    k(r2f) = (18100000_pr)
    k_0(FOr3f) = (1.33e+48_pr)*exp((-24969.2752_pr)*R_T_inv + T_log *(-12_pr))
    k_inf(FOr3f) = (100000000_pr)
    FC(FOr3f) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(10970_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r3f) = getlindratecoeff(Tloc,k_0(FOr3f),k_inf(FOr3f),FC(FOr3f), M(mM3),Ploc)
    k_0(FOr4f) = (6.8e+49_pr)*exp((-25104_pr)*R_T_inv + T_log *(-13.42_pr))
    k_inf(FOr4f) = (960000000_pr)*exp(T_log *(-0.5_pr))
    FC(FOr4f) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(1434_pr)))&
              + exp(-(5329_pr)/Tloc)
    k(r4f) = getlindratecoeff(Tloc,k_0(FOr4f),k_inf(FOr4f),FC(FOr4f), M(mM4),Ploc)
    k(r5f) = (2e+16_pr)*exp((-326352_pr)*R_T_inv)
    k(r6f) = (7500000000000_pr)*exp((-167360_pr)*R_T_inv)
    k(r7f) = (50000_pr)*exp((-117152_pr)*R_T_inv)
    k(r8f) = (13330000_pr)
    k(r9f) = (1000000_pr)
    k(r10f) = (40000000_pr)*exp((-195811.2_pr)*R_T_inv)
    k(r11f) = (57800000_pr)*exp((-215894.4_pr)*R_T_inv)
    k(r12f) = (25000_pr)*exp((-6276_pr)*R_T_inv)
    k(r13f) = (500000_pr)
    k(r14f) = (700000_pr)
    k(r15f) = (700000_pr)
    k_0(FOr16f) = (39000_pr)*exp((2343.04_pr)*R_T_inv)
    k_inf(FOr16f) = (10000000_pr)*exp((-11589.68_pr)*R_T_inv)
    k(r16f) = getlindratecoeff(Tloc,k_0(FOr16f),k_inf(FOr16f),1.0_pr, M(mM16),Ploc)
    k_0(FOr17f) = (4600000_pr)*exp((-4476.88_pr)*R_T_inv)
    k_inf(FOr17f) = (17700000_pr)*exp((-8828.24_pr)*R_T_inv)
    FC(FOr17f) = ((1.0_pr - 1_pr)*exp(-Tloc/(1e-15_pr))) + (1_pr*exp(-Tloc/(95_pr)))&
              + exp(-(200_pr)/Tloc)
    k(r17f) = getlindratecoeff(Tloc,k_0(FOr17f),k_inf(FOr17f),FC(FOr17f), M(mM17),Ploc)
    k_0(FOr18f) = (5.66e+21_pr)*exp((-18610.432_pr)*R_T_inv + T_log *(-5_pr))
    k_inf(FOr18f) = (120000_pr)*exp((-12581.288_pr)*R_T_inv + T_log *(0.69_pr))
    k(r18f) = getlindratecoeff(Tloc,k_0(FOr18f),k_inf(FOr18f),1.0_pr, M(mM18),Ploc)
    k(r19f) = (20000000000000_pr)*exp((-158992_pr)*R_T_inv)
    k(r20f) = (50000000000000_pr)*exp((-171544_pr)*R_T_inv)
    k(r21f) = (2000000000000_pr)*exp((-175728_pr)*R_T_inv)
    k(r22f) = (10000000000000_pr)*exp((-133888_pr)*R_T_inv)
    k(r23f) = (10000000000000_pr)*exp((-165268_pr)*R_T_inv)
    k(r24f) = (200000000000000_pr)*exp((-213384_pr)*R_T_inv)
    k(r25f) = (2500000000000_pr)*exp((-150624_pr)*R_T_inv)
    k(r26f) = (25000000000000_pr)*exp((-158992_pr)*R_T_inv)
    k(r27f) = (10000000000000_pr)*exp((-158992_pr)*R_T_inv)
    k(r28f) = (140000000_pr)*exp((-58576_pr)*R_T_inv)
    k(r29f) = (25000000_pr)*exp((-129704_pr)*R_T_inv)
    k(r30f) = (100000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r31f) = (300000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r32f) = (180000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r33f) = (200000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r34f) = (100000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r35f) = (50000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r36f) = (300000_pr)*exp((-75312_pr)*R_T_inv)
    k(r37f) = (200000_pr)*exp((-75312_pr)*R_T_inv)
    k(r38f) = (200000_pr)*exp((-25104_pr)*R_T_inv)
    k(r39f) = (200000_pr)*exp((-54392_pr)*R_T_inv)
    k(r40f) = (300000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r41f) = (300000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r42f) = (30000000000000_pr)*exp((-297064_pr)*R_T_inv)
    k(r43f) = (18000000000000_pr)*exp((-326352_pr)*R_T_inv)
    k(r44f) = (50000000000000_pr)*exp((-292880_pr)*R_T_inv)
    k(r45f) = (50000000000000_pr)*exp((-292880_pr)*R_T_inv)
    k(r46f) = (150000_pr)*exp((-94140_pr)*R_T_inv)
    k(r47f) = (10_pr)*exp((-17572.8_pr)*R_T_inv + T_log *(2_pr))
    k(r48f) = (50_pr)*exp((-20920_pr)*R_T_inv + T_log *(2_pr))
    k(r49f) = (40000000_pr)*exp((-66944_pr)*R_T_inv)
    k(r50f) = (30_pr)*exp((-41840_pr)*R_T_inv + T_log *(2_pr))
    k(r51f) = (143000000_pr)*exp((-43932_pr)*R_T_inv)
    k(r52f) = (127000000_pr)*exp((-43932_pr)*R_T_inv)
    k(r53f) = (32000000_pr)*exp((-33053.6_pr)*R_T_inv)
    k(r54f) = (1.5e+15_pr)*exp((-340996_pr)*R_T_inv)
    k(r55f) = (400000_pr)*exp((-92048_pr)*R_T_inv)
    k(r56f) = (1600000000_pr)*exp((-33158.2_pr)*R_T_inv)
    k(r57f) = (180000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r58f) = (12000000_pr)*exp((-106692_pr)*R_T_inv)
    k(r59f) = (1000000_pr)*exp((-12552_pr)*R_T_inv)
    k(r60f) = (1500000_pr)*exp((-12552_pr)*R_T_inv)
    k(r61f) = (500000_pr)*exp((-12552_pr)*R_T_inv)
    k(r62f) = (150000_pr)*exp((-54392_pr)*R_T_inv)
    k(r63f) = (10000000_pr)*exp((-12552_pr)*R_T_inv)
    k(r64f) = (15000000_pr)*exp((-8368_pr)*R_T_inv)
    k(r65f) = (6000000_pr)
    k(r66f) = (300000_pr)*exp((-125520_pr)*R_T_inv)
    k(r67f) = (300000_pr)*exp((-125520_pr)*R_T_inv)
    k(r68f) = (22.3_pr)*exp((-5732.08_pr)*R_T_inv + T_log *(1.31_pr))
    k(r69f) = (2000000_pr)*exp((-69036_pr)*R_T_inv)
    k(r70f) = (400000_pr)*exp((-62760_pr)*R_T_inv)
    k(r71f) = (300000_pr)*exp((-29288_pr)*R_T_inv)
    k(r72f) = (400000_pr)*exp((-54392_pr)*R_T_inv)
    k(r73f) = (5000000_pr)
    k(r74f) = (30000000_pr)*exp((-8368_pr)*R_T_inv)
    k(r75f) = (250000_pr)*exp((-25104_pr)*R_T_inv)
    k(r76f) = (1000000_pr)*exp((-12552_pr)*R_T_inv)
    k(r77f) = (100000_pr)*exp((-92048_pr)*R_T_inv)
    k(r78f) = (200000_pr)*exp((-25104_pr)*R_T_inv)
    k(r79f) = (3e+16_pr)*exp((-347272_pr)*R_T_inv)
    k(r80f) = (1e+16_pr)*exp((-301248_pr)*R_T_inv)
    k(r81f) = (5e+16_pr)*exp((-334720_pr)*R_T_inv)
    k(r82f) = (30000000000000_pr)*exp((-154808_pr)*R_T_inv)
    k(r83f) = (20000000000000_pr)*exp((-125520_pr)*R_T_inv)
    k(r84f) = (15000000000000_pr)*exp((-161084_pr)*R_T_inv)
    k(r85f) = (1000000000000_pr)*exp((-144348_pr)*R_T_inv)
    k(r86f) = (25000000000000_pr)*exp((-135980_pr)*R_T_inv)
    k(r87f) = (30000000000000_pr)*exp((-170707.2_pr)*R_T_inv)
    k(r88f) = (150000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r89f) = (150000_pr)*exp((-31798.4_pr)*R_T_inv)
    k(r90f) = (50000000000000_pr)*exp((-297064_pr)*R_T_inv)
    k(r91f) = (95000000_pr)*exp((-43932_pr)*R_T_inv)
    k(r92f) = (63500000_pr)*exp((-33053.6_pr)*R_T_inv)
    k(r93f) = (0.9496_pr)*exp((-35395.67768_pr)*R_T_inv + T_log *(2_pr))
    k(r94f) = (0.378_pr)*exp((-34093.49136_pr)*R_T_inv + T_log *(2_pr))
    k(r95f) = (0.4759_pr)*exp((-85363.47424_pr)*R_T_inv + T_log *(2_pr))
    k(r96f) = (0.378_pr)*exp((-34194.3676_pr)*R_T_inv + T_log *(2_pr))
    k(r97f) = (0.1609_pr)*exp((-61178.40616_pr)*R_T_inv + T_log *(2_pr))
    k(r98f) = (0.1685_pr)*exp((-98187.30872_pr)*R_T_inv + T_log *(2_pr))
    k(r99f) = (0.09473_pr)*exp((-66357.4032_pr)*R_T_inv + T_log *(2_pr))
    k(r100f) = (0.4748_pr)*exp((-46415.6224_pr)*R_T_inv + T_log *(2_pr))
    k(r101f) = (0.189_pr)*exp((-45021.26256_pr)*R_T_inv + T_log *(2_pr))
    k(r102f) = (19.25_pr)*exp((-43554.47768_pr)*R_T_inv + T_log *(2_pr))
    k(r103f) = (0.3122_pr)*exp((-47670.8224_pr)*R_T_inv + T_log *(2_pr))
    k(r104f) = (0.1839_pr)*exp((-62463.85648_pr)*R_T_inv + T_log *(2_pr))
    k(r105f) = (0.1925_pr)*exp((-99676.60352_pr)*R_T_inv + T_log *(2_pr))
    k(r106f) = (0.3513_pr)*exp((-31890.40616_pr)*R_T_inv + T_log *(2_pr))
    k(r107f) = (0.2166_pr)*exp((-77437.76488_pr)*R_T_inv + T_log *(2_pr))
    k(r108f) = (0.1218_pr)*exp((-47940.56488_pr)*R_T_inv + T_log *(2_pr))
    k(r109f) = (0.6105_pr)*exp((-31920.65648_pr)*R_T_inv + T_log *(2_pr))
    k(r110f) = (0.243_pr)*exp((-30665.45648_pr)*R_T_inv + T_log *(2_pr))
    k(r111f) = (0.1171_pr)*exp((-23502.90872_pr)*R_T_inv + T_log *(2_pr))
    k(r112f) = (0.06895_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r113f) = (0.0406_pr)*exp((-33770.69576_pr)*R_T_inv + T_log *(2_pr))
    k(r114f) = (0.2035_pr)*exp((-23737.00352_pr)*R_T_inv + T_log *(2_pr))
    k(r115f) = (0.08101_pr)*exp((-64988.23104_pr)*R_T_inv + T_log *(2_pr))
    k(r116f) = (0.08101_pr)*exp((-22579.20704_pr)*R_T_inv + T_log *(2_pr))
    k(r117f) = (0.102_pr)*exp((-61635.00608_pr)*R_T_inv + T_log *(2_pr))
    k(r118f) = (0.08101_pr)*exp((-22345.11224_pr)*R_T_inv + T_log *(2_pr))
    k(r119f) = (0.04452_pr)*exp((-30638.30232_pr)*R_T_inv + T_log *(2_pr))
    k(r120f) = (10.53_pr)*exp((-42353.1676_pr)*R_T_inv + T_log *(2_pr))
    k(r121f) = (0.1708_pr)*exp((-46341.7748_pr)*R_T_inv + T_log *(2_pr))
    k(r122f) = (0.1006_pr)*exp((-61178.40616_pr)*R_T_inv + T_log *(2_pr))
    k(r123f) = (0.1053_pr)*exp((-98284.71224_pr)*R_T_inv + T_log *(2_pr))
    k(r124f) = (0.2967_pr)*exp((-46266.00256_pr)*R_T_inv + T_log *(2_pr))
    k(r125f) = (0.3122_pr)*exp((-31445.56328_pr)*R_T_inv + T_log *(2_pr))
    k(r126f) = (0.1839_pr)*exp((-44038.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r127f) = (0.1925_pr)*exp((-77749.18_pr)*R_T_inv + T_log *(2_pr))
    k(r128f) = (0.5426_pr)*exp((-31442.46712_pr)*R_T_inv + T_log *(2_pr))
    k(r129f) = (0.07806_pr)*exp((-21170.2032_pr)*R_T_inv + T_log *(2_pr))
    k(r130f) = (0.04597_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r131f) = (0.04813_pr)*exp((-63267.89576_pr)*R_T_inv + T_log *(2_pr))
    k(r132f) = (57.76_pr)*exp((-43507.49136_pr)*R_T_inv + T_log *(2_pr))
    k(r133f) = (0.9367_pr)*exp((-47531.66256_pr)*R_T_inv + T_log *(2_pr))
    k(r134f) = (0.5516_pr)*exp((-62463.85648_pr)*R_T_inv + T_log *(2_pr))
    k(r135f) = (1.628_pr)*exp((-47455.89032_pr)*R_T_inv + T_log *(2_pr))
    k(r136f) = (12.03_pr)*exp((-17582.67424_pr)*R_T_inv + T_log *(2_pr))
    k(r137f) = (0.1952_pr)*exp((-23502.90872_pr)*R_T_inv + T_log *(2_pr))
    k(r138f) = (0.1149_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r139f) = (0.1203_pr)*exp((-61635.00608_pr)*R_T_inv + T_log *(2_pr))
    k(r140f) = (0.3391_pr)*exp((-23737.00352_pr)*R_T_inv + T_log *(2_pr))
    k(r141f) = (7.22_pr)*exp((-27302.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r142f) = (0.1171_pr)*exp((-31486.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r143f) = (0.06895_pr)*exp((-44038.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r144f) = (0.0722_pr)*exp((-77720.18488_pr)*R_T_inv + T_log *(2_pr))
    k(r145f) = (0.2035_pr)*exp((-31486.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r146f) = (0.3903_pr)*exp((-20381.47736_pr)*R_T_inv + T_log *(2_pr))
    k(r147f) = (0.2298_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r148f) = (0.2407_pr)*exp((-63819.97456_pr)*R_T_inv + T_log *(2_pr))
    k(r149f) = (0.6783_pr)*exp((-20381.47736_pr)*R_T_inv + T_log *(2_pr))
    k(r150f) = (24.07_pr)*exp((-11142.28488_pr)*R_T_inv + T_log *(2_pr))
    k(r151f) = (0.3903_pr)*exp((-14828.76544_pr)*R_T_inv + T_log *(2_pr))
    k(r152f) = (0.2298_pr)*exp((-26042.59672_pr)*R_T_inv + T_log *(2_pr))
    k(r153f) = (0.2407_pr)*exp((-56869.89032_pr)*R_T_inv + T_log *(2_pr))
    k(r154f) = (0.6783_pr)*exp((-14828.76544_pr)*R_T_inv + T_log *(2_pr))
    k(r155f) = (19.25_pr)*exp((-16529.18488_pr)*R_T_inv + T_log *(2_pr))
    k(r156f) = (0.3122_pr)*exp((-20381.47736_pr)*R_T_inv + T_log *(2_pr))
    k(r157f) = (0.1839_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r158f) = (0.1925_pr)*exp((-63819.97456_pr)*R_T_inv + T_log *(2_pr))
    k(r159f) = (0.2342_pr)*exp((-31569.24232_pr)*R_T_inv + T_log *(2_pr))
    k(r160f) = (0.1379_pr)*exp((-44038.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r161f) = (0.1444_pr)*exp((-77662.61304_pr)*R_T_inv + T_log *(2_pr))
    k(r162f) = (0.0812_pr)*exp((-39436.58488_pr)*R_T_inv + T_log *(2_pr))
    k(r163f) = (0.0812_pr)*exp((-48165.41304_pr)*R_T_inv + T_log *(2_pr))
    k(r164f) = (0.407_pr)*exp((-31575.3928_pr)*R_T_inv + T_log *(2_pr))
    k(r165f) = (0.162_pr)*exp((-75507.3928_pr)*R_T_inv + T_log *(2_pr))
    k(r166f) = (0.162_pr)*exp((-30320.1928_pr)*R_T_inv + T_log *(2_pr))
    k(r167f) = (0.204_pr)*exp((-77662.61304_pr)*R_T_inv + T_log *(2_pr))
    k(r168f) = (0.162_pr)*exp((-30314.04232_pr)*R_T_inv + T_log *(2_pr))
    k(r169f) = (0.0812_pr)*exp((-39854.98488_pr)*R_T_inv + T_log *(2_pr))
    k(r170f) = (0.08904_pr)*exp((-44818.21304_pr)*R_T_inv + T_log *(2_pr))
    k(r171f) = (0.1561_pr)*exp((-20863.34864_pr)*R_T_inv + T_log *(2_pr))
    k(r172f) = (0.09193_pr)*exp((-32041.36488_pr)*R_T_inv + T_log *(2_pr))
    k(r173f) = (0.09627_pr)*exp((-63482.70232_pr)*R_T_inv + T_log *(2_pr))
    k(r174f) = (0.05413_pr)*exp((-27750.38_pr)*R_T_inv + T_log *(2_pr))
    k(r175f) = (0.05413_pr)*exp((-35618.35016_pr)*R_T_inv + T_log *(2_pr))
    k(r176f) = (0.2713_pr)*exp((-20899.4984_pr)*R_T_inv + T_log *(2_pr))
    k(r177f) = (0.108_pr)*exp((-62150.76776_pr)*R_T_inv + T_log *(2_pr))
    k(r178f) = (0.108_pr)*exp((-19741.74376_pr)*R_T_inv + T_log *(2_pr))
    k(r179f) = (0.136_pr)*exp((-63482.70232_pr)*R_T_inv + T_log *(2_pr))
    k(r180f) = (0.108_pr)*exp((-19705.594_pr)*R_T_inv + T_log *(2_pr))
    k(r181f) = (0.05413_pr)*exp((-28139.78488_pr)*R_T_inv + T_log *(2_pr))
    k(r182f) = (0.05413_pr)*exp((-28919.01304_pr)*R_T_inv + T_log *(2_pr))
    k_0(FOr1b) = (2.33396874365701e+35_pr)*exp((-373940.431000222_pr)*R_T_inv + T_log *(-6.05404325307077_pr))
    k_inf(FOr1b) = (2.50425829536754e+20_pr)*exp((-378961.230995734_pr)*R_T_inv + T_log *(-1.02404325241078_pr))
    FC(FOr1b) = ((1.0_pr - 0.38_pr)*exp(-Tloc/(73_pr))) + (0.38_pr*exp(-Tloc/(1180_pr)))
    k(r1b) = getlindratecoeff(Tloc,k_0(FOr1b),k_inf(FOr1b),FC(FOr1b), M(mM1),Ploc)
    k(r2b) = (54587.3303484287_pr)*exp((-179994.252427512_pr)*R_T_inv + T_log *(1.03389144562243_pr))
    k_0(FOr3b) = (7.83504117360934e+53_pr)*exp((-384437.328426216_pr)*R_T_inv + T_log *(-11.2054855708216_pr))
    k_inf(FOr3b) = (58910081680886.8_pr)*exp((-359468.053197589_pr)*R_T_inv + T_log *(0.794514434330109_pr))
    FC(FOr3b) = ((1.0_pr - 0.02_pr)*exp(-Tloc/(1097_pr))) + (0.02_pr*exp(-Tloc/(10970_pr)))&
              + exp(-(6860_pr)/Tloc)
    k(r3b) = getlindratecoeff(Tloc,k_0(FOr3b),k_inf(FOr3b),FC(FOr3b), M(mM3),Ploc)
    k_0(FOr4b) = (6.62693111881404e+65_pr)*exp((-396188.076251446_pr)*R_T_inv + T_log *(-15.2453750940007_pr))
    k_inf(FOr4b) = (9.35566728949499e+24_pr)*exp((-371084.076238164_pr)*R_T_inv + T_log *(-2.3253750915958_pr))
    FC(FOr4b) = ((1.0_pr - 1_pr)*exp(-Tloc/(1000_pr))) + (1_pr*exp(-Tloc/(1434_pr)))&
              + exp(-(5329_pr)/Tloc)
    k(r4b) = getlindratecoeff(Tloc,k_0(FOr4b),k_inf(FOr4b),FC(FOr4b), M(mM4),Ploc)
    k(r5b) = (35882.1420616314_pr)*exp((-13572.7368598865_pr)*R_T_inv + T_log *(0.685910145174778_pr))
    k(r6b) = (0.0033721928926209_pr)*exp((-8687.6650417772_pr)*R_T_inv + T_log *(2.49227942043424_pr))
    k(r7b) = (9.61090660120252e+17_pr)*exp((-289417.754534428_pr)*R_T_inv + T_log *(-1.50491928469738_pr))
    k(r8b) = (20093167845.6705_pr)*exp((-281239.169188498_pr)*R_T_inv + T_log *(-0.387410421998399_pr))
    k(r10b) = (722772.077173024_pr)*exp((7216.60158853646_pr)*R_T_inv + T_log *(-0.163505731289003_pr))
    k(r11b) = (10111483.086977_pr)*exp((-7926.79186180594_pr)*R_T_inv + T_log *(-0.569700375176713_pr))
    k_0(FOr16b) = (415140982.982556_pr)*exp((-154097.211613388_pr)*R_T_inv + T_log *(0.631008702262352_pr))
    k_inf(FOr16b) = (106446402463.238_pr)*exp((-168029.931589869_pr)*R_T_inv + T_log *(0.631008706449389_pr))
    k(r16b) = getlindratecoeff(Tloc,k_0(FOr16b),k_inf(FOr16b),1.0_pr, M(mM16),Ploc)
    k_0(FOr17b) = (474061274975.222_pr)*exp((-155977.325063102_pr)*R_T_inv + T_log *(0.224814058468624_pr))
    k_inf(FOr17b) = (1824105401631.21_pr)*exp((-160328.685087344_pr)*R_T_inv + T_log *(0.22481405411987_pr))
    FC(FOr17b) = ((1.0_pr - 1_pr)*exp(-Tloc/(1e-15_pr))) + (1_pr*exp(-Tloc/(95_pr)))&
              + exp(-(200_pr)/Tloc)
    k(r17b) = getlindratecoeff(Tloc,k_0(FOr17b),k_inf(FOr17b),FC(FOr17b), M(mM17),Ploc)
    k_0(FOr18b) = (4.99774826300176e+29_pr)*exp((-271598.801809528_pr)*R_T_inv + T_log *(-5.7073525572355_pr))
    k_inf(FOr18b) = (10595927785180.8_pr)*exp((-265569.657471017_pr)*R_T_inv + T_log *(-0.0173524967753623_pr))
    k(r18b) = getlindratecoeff(Tloc,k_0(FOr18b),k_inf(FOr18b),1.0_pr, M(mM18),Ploc)
    k(r19b) = (797.113915553473_pr)*exp((-54137.4595947394_pr)*R_T_inv + T_log *(1.09451894375325_pr))
    k(r20b) = (4526434.47233418_pr)*exp((-14711.604934024_pr)*R_T_inv + T_log *(0.321737865141025_pr))
    k(r21b) = (408375714568.843_pr)*exp((-161267.743995236_pr)*R_T_inv + T_log *(0.0491701185767253_pr))
    k(r22b) = (4.45599648789442_pr)*exp((-27752.1098878131_pr)*R_T_inv + T_log *(1.55137247872864_pr))
    k(r23b) = (4433598.03570279_pr)*exp((-22895.86092952_pr)*R_T_inv + T_log *(0.272567748214871_pr))
    k(r24b) = (34535143.2025942_pr)*exp((-19810.036841356_pr)*R_T_inv + T_log *(0.28120595899214_pr))
    k(r25b) = (18128826636765.1_pr)*exp((-217915.380740639_pr)*R_T_inv + T_log *(-0.128249528124765_pr))
    k(r26b) = (0.0173355507581085_pr)*exp((-16883.9147151146_pr)*R_T_inv + T_log *(2.28888440802617_pr))
    k(r27b) = (12521634.1354508_pr)*exp((-32709.4176782456_pr)*R_T_inv + T_log *(0.152956413678047_pr))
    k(r28b) = (345529638696.139_pr)*exp((-18628.1905734005_pr)*R_T_inv + T_log *(-0.660202460456301_pr))
    k(r42b) = (19100.028725938_pr)*exp((-159637.136801569_pr)*R_T_inv + T_log *(0.465565612795451_pr))
    k(r43b) = (92150.1711432488_pr)*exp((-146878.199205491_pr)*R_T_inv + T_log *(0.23937701558639_pr))
    k(r44b) = (13190.8416553463_pr)*exp((-167041.85195304_pr)*R_T_inv + T_log *(0.648295195584218_pr))
    k(r45b) = (1052615.14277896_pr)*exp((-181673.770104686_pr)*R_T_inv + T_log *(0.312154934895535_pr))
    k(r47b) = (4.52044921834736_pr)*exp((-91087.3690290153_pr)*R_T_inv + T_log *(1.53202445594049_pr))
    k(r48b) = (0.37475132043283_pr)*exp((-88331.2390190665_pr)*R_T_inv + T_log *(2.05913491311201_pr))
    k(r49b) = (18063884.7140573_pr)*exp((-139052.037508936_pr)*R_T_inv + T_log *(-0.0238056815971398_pr))
    k(r50b) = (0.497901650896935_pr)*exp((-37143.2020287512_pr)*R_T_inv + T_log *(2.08294050348568_pr))
    k(r51b) = (9382.64882340023_pr)*exp((-58005.5818815662_pr)*R_T_inv + T_log *(0.690379668261717_pr))
    k(r52b) = (75570.0892234335_pr)*exp((-60465.9910695623_pr)*R_T_inv + T_log *(0.375727439193656_pr))
    k(r53b) = (93253.7460821326_pr)*exp((-64047.8470841944_pr)*R_T_inv + T_log *(0.326557318737782_pr))
    k(r54b) = (7808145.87753833_pr)*exp((5588.30964519984_pr)*R_T_inv + T_log *(0.397573133254674_pr))
    k(r55b) = (4.76575667825343e+19_pr)*exp((-387637.3063435_pr)*R_T_inv + T_log *(-1.23347998707587_pr))
    k(r79b) = (0.00124203954410609_pr)*exp((17958.6660477153_pr)*R_T_inv + T_log *(2.75040695556418_pr))
    k(r80b) = (37.9429370665294_pr)*exp((15848.5919726153_pr)*R_T_inv + T_log *(1.29991525025141_pr))
    k(r81b) = (19615.1859324796_pr)*exp((-19735.5916186733_pr)*R_T_inv + T_log *(0.97861765379411_pr))
    k(r83b) = (0.00683798547899473_pr)*exp((-24289.8559296036_pr)*R_T_inv + T_log *(2.33791467078461_pr))
    k(r84b) = (1956165.7745632_pr)*exp((-16876.6263472893_pr)*R_T_inv + T_log *(0.238887820298105_pr))
    k(r85b) = (562697750601.446_pr)*exp((-157274.353693947_pr)*R_T_inv + T_log *(0.139345018532055_pr))
    k(r86b) = (6.12374682063169_pr)*exp((-25535.0627265987_pr)*R_T_inv + T_log *(1.5799673783301_pr))
    k(r87b) = (6952811.07043107_pr)*exp((-13573.4725790357_pr)*R_T_inv + T_log *(0.0995428150392337_pr))
    k(r90b) = (2040.30447605918_pr)*exp((-170338.281566626_pr)*R_T_inv + T_log *(0.753104930950074_pr))
    k(r91b) = (29725.8412766681_pr)*exp((-61413.6551660327_pr)*R_T_inv + T_log *(0.514217120129066_pr))
    k(r92b) = (11180.4547179744_pr)*exp((-63461.6088825114_pr)*R_T_inv + T_log *(0.653562134638527_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(sCH3)**2.0_pr 
    w(r2f) = k(r2f) * c(sCH2CHCH2) * c(sH) 
    w(r3f) = k(r3f) * c(sCH2CHCH2) * c(sH) 
    w(r4f) = k(r4f) * c(sC2H5) * c(sCH3) 
    w(r5f) = k(r5f) * c(sNC4H8) 
    w(r6f) = k(r6f) * c(sC4H5) 
    w(r7f) = k(r7f) * c(sC2H2) * c(sC2H4) 
    w(r8f) = k(r8f) * c(sC2H3) * c(sCH3) 
    w(r9f) = k(r9f) * c(sC2H5) * c(sSC4H7) 
    w(r10f) = k(r10f) * c(sC2H2) * c(sC3H6) 
    w(r11f) = k(r11f) * c(sC2H4) * c(sC3H6) 
    w(r12f) = k(r12f) * c(sCH2CHCH2)**2.0_pr 
    w(r13f) = k(r13f) * c(sC3H3) * c(sCH2CHCH2) 
    w(r14f) = k(r14f) * c(sCH2CHCH2) * c(sSC4H7) 
    w(r15f) = k(r15f) * c(sSC4H7)**2.0_pr 
    w(r16f) = k(r16f) * c(sC2H2) * c(sH) 
    w(r17f) = k(r17f) * c(sC2H4) * c(sH) 
    w(r18f) = k(r18f) * c(sAC3H4) * c(sH) 
    w(r19f) = k(r19f) * c(sCHCHCH3) 
    w(r20f) = k(r20f) * c(sIC3H7) 
    w(r21f) = k(r21f) * c(sIC3H7) 
    w(r22f) = k(r22f) * c(sNC3H7) 
    w(r23f) = k(r23f) * c(sNC3H7) 
    w(r24f) = k(r24f) * c(sSC4H7) 
    w(r25f) = k(r25f) * c(sCH2C3H5) 
    w(r26f) = k(r26f) * c(sCH2C3H5) 
    w(r27f) = k(r27f) * c(sCH2C3H5) 
    w(r28f) = k(r28f) * c(sCH3)**2.0_pr 
    w(r29f) = k(r29f) * c(sCH3) * c(sCH4) 
    w(r30f) = k(r30f) * c(sC4H6) * c(sCH3) 
    w(r31f) = k(r31f) * c(sC2H2) * c(sC2H5) 
    w(r32f) = k(r32f) * c(sC2H5) * c(sC3H6) 
    w(r33f) = k(r33f) * c(sC2H4) * c(sNC3H7) 
    w(r34f) = k(r34f) * c(sC2H4) * c(sIC3H7) 
    w(r35f) = k(r35f) * c(sC3H6) * c(sIC3H7) 
    w(r36f) = k(r36f) * c(sC3H6) * c(sCH2CHCH2) 
    w(r37f) = k(r37f) * c(sC3H6) * c(sCH2CHCH2) 
    w(r38f) = k(r38f) * c(sC2H4) * c(sCHCHCH3) 
    w(r39f) = k(r39f) * c(sC2H4) * c(sSC4H7) 
    w(r40f) = k(r40f) * c(sC2H4) * c(sCH2C3H5) 
    w(r41f) = k(r41f) * c(sC3H6) * c(sCH2C3H5) 
    w(r42f) = k(r42f) * c(sC2H6) 
    w(r43f) = k(r43f) * c(sC3H6) 
    w(r44f) = k(r44f) * c(sC3H8) 
    w(r45f) = k(r45f) * c(sNC4H8) 
    w(r46f) = k(r46f) * c(sC2H2) * c(sC4H6) 
    w(r47f) = k(r47f) * c(sC3H6) * c(sH) 
    w(r48f) = k(r48f) * c(sAC3H4) * c(sH) 
    w(r49f) = k(r49f) * c(sAC3H4) * c(sCH3) 
    w(r50f) = k(r50f) * c(sCH4) * c(sH) 
    w(r51f) = k(r51f) * c(sC2H6) * c(sH) 
    w(r52f) = k(r52f) * c(sC3H8) * c(sH) 
    w(r53f) = k(r53f) * c(sC3H8) * c(sH) 
    w(r54f) = k(r54f) * c(sCYC5H6) 
    w(r55f) = k(r55f) * c(sAC3H4) * c(sC2H2) 
    w(r56f) = k(r56f) * c(sCYC5H6) * c(sH) 
    w(r57f) = k(r57f) * c(sCYC5H5) * c(sH2) 
    w(r58f) = k(r58f) * c(sCYC5H5) * c(sCYC5H6) 
    w(r59f) = k(r59f) * c(sCH3) * c(sCYC5H5) 
    w(r60f) = k(r60f) * c(sCH3) * c(sCYC5H5) 
    w(r61f) = k(r61f) * c(sAC3H4) * c(sC2H3) 
    w(r62f) = k(r62f) * c(sC2H4) * c(sSC4H7) 
    w(r63f) = k(r63f) * c(sCH3) * c(sCYC5H5) 
    w(r64f) = k(r64f) * c(sH) * c(sMCPTD) 
    w(r65f) = k(r65f) * c(sCH2CHCH2) * c(sCHCHCH3) 
    w(r66f) = k(r66f) * c(sC2H4) * c(sCYC5H6) 
    w(r67f) = k(r67f) * c(sC3H6) * c(sCYC5H6) 
    w(r68f) = k(r68f) * c(sC2H3) * c(sC4H6) 
    w(r69f) = k(r69f) * c(sAC3H4) * c(sCH2CHCH2) 
    w(r70f) = k(r70f) * c(sC4H6) * c(sCH2CHCH2) 
    w(r71f) = k(r71f) * c(sC2H4) * c(sC4H5) 
    w(r72f) = k(r72f) * c(sC3H6) * c(sSC4H7) 
    w(r73f) = k(r73f) * c(sC2H3) * c(sCH2CHCH2) 
    w(r74f) = k(r74f) * c(sC5H8) * c(sH) 
    w(r75f) = k(r75f) * c(sC2H3) * c(sC3H6) 
    w(r76f) = k(r76f) * c(sC2H3) * c(sC5H8) 
    w(r77f) = k(r77f) * c(sC2H4) * c(sCH2CHCH2) 
    w(r78f) = k(r78f) * c(sC3H6) * c(sCHCHCH3) 
    w(r79f) = k(r79f) * c(sNC4H10) 
    w(r80f) = k(r80f) * c(sNC5H10) 
    w(r81f) = k(r81f) * c(sNC5H10) 
    w(r82f) = k(r82f) * c(sNC5H9X3) 
    w(r83f) = k(r83f) * c(sNC4H9P) 
    w(r84f) = k(r84f) * c(sNC4H9P) 
    w(r85f) = k(r85f) * c(sNC4H9P) 
    w(r86f) = k(r86f) * c(sNC4H9S) 
    w(r87f) = k(r87f) * c(sNC4H9S) 
    w(r88f) = k(r88f) * c(sC2H5) * c(sNC4H8) 
    w(r89f) = k(r89f) * c(sC2H4) * c(sNC4H9P) 
    w(r90f) = k(r90f) * c(sNC4H10) 
    w(r91f) = k(r91f) * c(sH) * c(sNC4H10) 
    w(r92f) = k(r92f) * c(sH) * c(sNC4H10) 
    w(r93f) = k(r93f) * c(sC2H3) * c(sH2) 
    w(r94f) = k(r94f) * c(sC4H5) * c(sH2) 
    w(r95f) = k(r95f) * c(sH2) * c(sSC4H7) 
    w(r96f) = k(r96f) * c(sCHCHCH3) * c(sH2) 
    w(r97f) = k(r97f) * c(sC2H5) * c(sCH4) 
    w(r98f) = k(r98f) * c(sCH2CHCH2) * c(sCH4) 
    w(r99f) = k(r99f) * c(sCH4) * c(sIC3H7) 
    w(r100f) = k(r100f) * c(sC2H3) * c(sCH4) 
    w(r101f) = k(r101f) * c(sC4H5) * c(sCH4) 
    w(r102f) = k(r102f) * c(sC2H4) * c(sH) 
    w(r103f) = k(r103f) * c(sC2H4) * c(sCH3) 
    w(r104f) = k(r104f) * c(sC2H4) * c(sC2H5) 
    w(r105f) = k(r105f) * c(sC2H4) * c(sCH2CHCH2) 
    w(r106f) = k(r106f) * c(sC2H6) * c(sCH3) 
    w(r107f) = k(r107f) * c(sC2H6) * c(sCH2CHCH2) 
    w(r108f) = k(r108f) * c(sC2H6) * c(sIC3H7) 
    w(r109f) = k(r109f) * c(sC2H3) * c(sC2H6) 
    w(r110f) = k(r110f) * c(sC2H6) * c(sC4H5) 
    w(r111f) = k(r111f) * c(sC3H6) * c(sCH3) 
    w(r112f) = k(r112f) * c(sC2H5) * c(sC3H6) 
    w(r113f) = k(r113f) * c(sC3H6) * c(sIC3H7) 
    w(r114f) = k(r114f) * c(sC2H3) * c(sC3H6) 
    w(r115f) = k(r115f) * c(sC3H3) * c(sC3H6) 
    w(r116f) = k(r116f) * c(sC3H6) * c(sC4H5) 
    w(r117f) = k(r117f) * c(sC3H6) * c(sSC4H7) 
    w(r118f) = k(r118f) * c(sC3H6) * c(sCHCHCH3) 
    w(r119f) = k(r119f) * c(sC3H6) * c(sNC4H9S) 
    w(r120f) = k(r120f) * c(sC3H6) * c(sH) 
    w(r121f) = k(r121f) * c(sC3H6) * c(sCH3) 
    w(r122f) = k(r122f) * c(sC2H5) * c(sC3H6) 
    w(r123f) = k(r123f) * c(sC3H6) * c(sCH2CHCH2) 
    w(r124f) = k(r124f) * c(sC2H3) * c(sC3H6) 
    w(r125f) = k(r125f) * c(sC3H8) * c(sCH3) 
    w(r126f) = k(r126f) * c(sC2H5) * c(sC3H8) 
    w(r127f) = k(r127f) * c(sC3H8) * c(sCH2CHCH2) 
    w(r128f) = k(r128f) * c(sC2H3) * c(sC3H8) 
    w(r129f) = k(r129f) * c(sC3H8) * c(sCH3) 
    w(r130f) = k(r130f) * c(sC2H5) * c(sC3H8) 
    w(r131f) = k(r131f) * c(sC3H8) * c(sCH2CHCH2) 
    w(r132f) = k(r132f) * c(sC4H6) * c(sH) 
    w(r133f) = k(r133f) * c(sC4H6) * c(sCH3) 
    w(r134f) = k(r134f) * c(sC2H5) * c(sC4H6) 
    w(r135f) = k(r135f) * c(sC2H3) * c(sC4H6) 
    w(r136f) = k(r136f) * c(sH) * c(sNC4H8) 
    w(r137f) = k(r137f) * c(sCH3) * c(sNC4H8) 
    w(r138f) = k(r138f) * c(sC2H5) * c(sNC4H8) 
    w(r139f) = k(r139f) * c(sCH2CHCH2) * c(sNC4H8) 
    w(r140f) = k(r140f) * c(sC2H3) * c(sNC4H8) 
    w(r141f) = k(r141f) * c(sH) * c(sNC4H8) 
    w(r142f) = k(r142f) * c(sCH3) * c(sNC4H8) 
    w(r143f) = k(r143f) * c(sC2H5) * c(sNC4H8) 
    w(r144f) = k(r144f) * c(sCH2CHCH2) * c(sNC4H8) 
    w(r145f) = k(r145f) * c(sC2H3) * c(sNC4H8) 
    w(r146f) = k(r146f) * c(sCH3) * c(sCYC5H6) 
    w(r147f) = k(r147f) * c(sC2H5) * c(sCYC5H6) 
    w(r148f) = k(r148f) * c(sCH2CHCH2) * c(sCYC5H6) 
    w(r149f) = k(r149f) * c(sC2H3) * c(sCYC5H6) 
    w(r150f) = k(r150f) * c(sH) * c(sMCPTD) 
    w(r151f) = k(r151f) * c(sCH3) * c(sMCPTD) 
    w(r152f) = k(r152f) * c(sC2H5) * c(sMCPTD) 
    w(r153f) = k(r153f) * c(sCH2CHCH2) * c(sMCPTD) 
    w(r154f) = k(r154f) * c(sC2H3) * c(sMCPTD) 
    w(r155f) = k(r155f) * c(sH) * c(sNC5H10) 
    w(r156f) = k(r156f) * c(sCH3) * c(sNC5H10) 
    w(r157f) = k(r157f) * c(sC2H5) * c(sNC5H10) 
    w(r158f) = k(r158f) * c(sCH2CHCH2) * c(sNC5H10) 
    w(r159f) = k(r159f) * c(sCH3) * c(sNC4H10) 
    w(r160f) = k(r160f) * c(sC2H5) * c(sNC4H10) 
    w(r161f) = k(r161f) * c(sCH2CHCH2) * c(sNC4H10) 
    w(r162f) = k(r162f) * c(sNC3H7) * c(sNC4H10) 
    w(r163f) = k(r163f) * c(sIC3H7) * c(sNC4H10) 
    w(r164f) = k(r164f) * c(sC2H3) * c(sNC4H10) 
    w(r165f) = k(r165f) * c(sC3H3) * c(sNC4H10) 
    w(r166f) = k(r166f) * c(sC4H5) * c(sNC4H10) 
    w(r167f) = k(r167f) * c(sNC4H10) * c(sSC4H7) 
    w(r168f) = k(r168f) * c(sCHCHCH3) * c(sNC4H10) 
    w(r169f) = k(r169f) * c(sCH2C3H5) * c(sNC4H10) 
    w(r170f) = k(r170f) * c(sNC4H10) * c(sNC4H9S) 
    w(r171f) = k(r171f) * c(sCH3) * c(sNC4H10) 
    w(r172f) = k(r172f) * c(sC2H5) * c(sNC4H10) 
    w(r173f) = k(r173f) * c(sCH2CHCH2) * c(sNC4H10) 
    w(r174f) = k(r174f) * c(sNC3H7) * c(sNC4H10) 
    w(r175f) = k(r175f) * c(sIC3H7) * c(sNC4H10) 
    w(r176f) = k(r176f) * c(sC2H3) * c(sNC4H10) 
    w(r177f) = k(r177f) * c(sC3H3) * c(sNC4H10) 
    w(r178f) = k(r178f) * c(sC4H5) * c(sNC4H10) 
    w(r179f) = k(r179f) * c(sNC4H10) * c(sSC4H7) 
    w(r180f) = k(r180f) * c(sCHCHCH3) * c(sNC4H10) 
    w(r181f) = k(r181f) * c(sCH2C3H5) * c(sNC4H10) 
    w(r182f) = k(r182f) * c(sNC4H10) * c(sNC4H9P) 
    w(r1b) = k(r1b) * c(sC2H6) 
    w(r2b) = k(r2b) * c(sAC3H4) * c(sH2) 
    w(r3b) = k(r3b) * c(sC3H6) 
    w(r4b) = k(r4b) * c(sC3H8) 
    w(r5b) = k(r5b) * c(sCH2CHCH2) * c(sCH3) 
    w(r6b) = k(r6b) * c(sC2H2) * c(sC2H3) 
    w(r7b) = k(r7b) * c(sC4H6) 
    w(r8b) = k(r8b) * c(sC2H2) * c(sCH4) 
    w(r10b) = k(r10b) * c(sC2H3) * c(sCH2CHCH2) 
    w(r11b) = k(r11b) * c(sC2H5) * c(sCH2CHCH2) 
    w(r16b) = k(r16b) * c(sC2H3) 
    w(r17b) = k(r17b) * c(sC2H5) 
    w(r18b) = k(r18b) * c(sCH2CHCH2) 
    w(r19b) = k(r19b) * c(sC2H2) * c(sCH3) 
    w(r20b) = k(r20b) * c(sC3H6) * c(sH) 
    w(r21b) = k(r21b) * c(sNC3H7) 
    w(r22b) = k(r22b) * c(sC2H4) * c(sCH3) 
    w(r23b) = k(r23b) * c(sC3H6) * c(sH) 
    w(r24b) = k(r24b) * c(sC4H6) * c(sH) 
    w(r25b) = k(r25b) * c(sSC4H7) 
    w(r26b) = k(r26b) * c(sC2H3) * c(sC2H4) 
    w(r27b) = k(r27b) * c(sC4H6) * c(sH) 
    w(r28b) = k(r28b) * c(sC2H5) * c(sH) 
    w(r42b) = k(r42b) * c(sC2H4) * c(sH2) 
    w(r43b) = k(r43b) * c(sAC3H4) * c(sH2) 
    w(r44b) = k(r44b) * c(sC3H6) * c(sH2) 
    w(r45b) = k(r45b) * c(sC4H6) * c(sH2) 
    w(r47b) = k(r47b) * c(sCH2CHCH2) * c(sH2) 
    w(r48b) = k(r48b) * c(sC3H3) * c(sH2) 
    w(r49b) = k(r49b) * c(sC3H3) * c(sCH4) 
    w(r50b) = k(r50b) * c(sCH3) * c(sH2) 
    w(r51b) = k(r51b) * c(sC2H5) * c(sH2) 
    w(r52b) = k(r52b) * c(sH2) * c(sNC3H7) 
    w(r53b) = k(r53b) * c(sH2) * c(sIC3H7) 
    w(r54b) = k(r54b) * c(sCYC5H5) * c(sH) 
    w(r55b) = k(r55b) * c(sCYC5H6) 
    w(r79b) = k(r79b) * c(sC2H5)**2.0_pr 
    w(r80b) = k(r80b) * c(sC2H5) * c(sCH2CHCH2) 
    w(r81b) = k(r81b) * c(sCH3) * c(sSC4H7) 
    w(r83b) = k(r83b) * c(sC2H4) * c(sC2H5) 
    w(r84b) = k(r84b) * c(sH) * c(sNC4H8) 
    w(r85b) = k(r85b) * c(sNC4H9S) 
    w(r86b) = k(r86b) * c(sC3H6) * c(sCH3) 
    w(r87b) = k(r87b) * c(sH) * c(sNC4H8) 
    w(r90b) = k(r90b) * c(sH2) * c(sNC4H8) 
    w(r91b) = k(r91b) * c(sH2) * c(sNC4H9P) 
    w(r92b) = k(r92b) * c(sH2) * c(sNC4H9S) 
  
    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(sH2O) = 0.0_pr 

    cdot(sH2) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         + 0.3_pr * w(r37f) &
         + w(r42f) &
         - w(r42b) &
         + w(r43f) &
         - w(r43b) &
         + w(r44f) &
         - w(r44b) &
         + w(r45f) &
         - w(r45b) &
         + w(r46f) &
         + w(r47f) &
         - w(r47b) &
         + w(r48f) &
         - w(r48b) &
         + w(r50f) &
         - w(r50b) &
         + w(r51f) &
         - w(r51b) &
         + w(r52f) &
         - w(r52b) &
         + w(r53f) &
         - w(r53b) &
         + w(r56f) &
         - w(r57f) &
         + w(r62f) &
         + w(r65f) &
         + 0.15_pr * w(r68f) &
         + w(r69f) &
         + 0.3_pr * w(r70f) &
         + w(r71f) &
         + 1.5_pr * w(r72f) &
         + w(r90f) &
         - w(r90b) &
         + w(r91f) &
         - w(r91b) &
         + w(r92f) &
         - w(r92b) &
         - w(r93f) &
         - w(r94f) &
         - w(r95f) &
         - w(r96f) &
         + w(r102f) &
         + w(r120f) &
         + w(r132f) &
         + w(r136f) &
         + w(r141f) &
         + w(r150f) &
         + w(r155f) 

    cdot(sCH4) = 0.0_pr &
         + w(r8f) &
         - w(r8b) &
         - w(r29f) &
         + w(r49f) &
         - w(r49b) &
         - w(r50f) &
         + w(r50b) &
         - w(r97f) &
         - w(r98f) &
         - w(r99f) &
         - w(r100f) &
         - w(r101f) &
         + w(r103f) &
         + w(r106f) &
         + w(r111f) &
         + w(r121f) &
         + w(r125f) &
         + w(r129f) &
         + w(r133f) &
         + w(r137f) &
         + w(r142f) &
         + w(r146f) &
         + w(r151f) &
         + w(r156f) &
         + w(r159f) &
         + w(r171f) 

    cdot(sC2H2) = 0.0_pr &
         + w(r6f) &
         - w(r6b) &
         - w(r7f) &
         + w(r7b) &
         + w(r8f) &
         - w(r8b) &
         - w(r10f) &
         + w(r10b) &
         + w(r13f) &
         - w(r16f) &
         + w(r16b) &
         + w(r19f) &
         - w(r19b) &
         - w(r31f) &
         - w(r46f) &
         - w(r55f) &
         + w(r55b) &
         + w(r60f) 

    cdot(sC2H4) = 0.0_pr &
         - w(r7f) &
         + w(r7b) &
         - w(r11f) &
         + w(r11b) &
         - w(r17f) &
         + w(r17b) &
         + w(r22f) &
         - w(r22b) &
         + w(r26f) &
         - w(r26b) &
         + w(r30f) &
         + 0.5_pr * w(r32f) &
         - w(r33f) &
         - w(r34f) &
         + 0.8_pr * w(r36f) &
         - w(r38f) &
         - w(r39f) &
         - w(r40f) &
         + w(r42f) &
         - w(r42b) &
         - w(r62f) &
         - w(r66f) &
         - w(r71f) &
         + 0.25_pr * w(r72f) &
         + 0.05_pr * w(r73f) &
         + 0.3_pr * w(r74f) &
         - w(r77f) &
         + w(r83f) &
         - w(r83b) &
         + 0.55_pr * w(r88f) &
         - w(r89f) &
         + w(r93f) &
         + w(r100f) &
         - w(r102f) &
         - w(r103f) &
         - w(r104f) &
         - w(r105f) &
         + w(r109f) &
         + w(r114f) &
         + w(r124f) &
         + w(r128f) &
         + w(r135f) &
         + w(r140f) &
         + w(r145f) &
         + w(r149f) &
         + w(r154f) &
         + w(r164f) &
         + w(r176f) 

    cdot(sC2H6) = 0.0_pr &
         + w(r1f) &
         - w(r1b) &
         + w(r9f) &
         + w(r29f) &
         - w(r42f) &
         + w(r42b) &
         - w(r51f) &
         + w(r51b) &
         + w(r97f) &
         + w(r104f) &
         - w(r106f) &
         - w(r107f) &
         - w(r108f) &
         - w(r109f) &
         - w(r110f) &
         + w(r112f) &
         + w(r122f) &
         + w(r126f) &
         + w(r130f) &
         + w(r134f) &
         + w(r138f) &
         + w(r143f) &
         + w(r147f) &
         + w(r152f) &
         + w(r157f) &
         + w(r160f) &
         + w(r172f) 

    cdot(sAC3H4) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         + w(r12f) &
         - w(r18f) &
         + w(r18b) &
         + w(r43f) &
         - w(r43b) &
         - w(r48f) &
         + w(r48b) &
         - w(r49f) &
         + w(r49b) &
         - w(r55f) &
         + w(r55b) &
         - w(r61f) &
         - w(r69f) &
         + 0.05_pr * w(r73f) &
         + w(r115f) &
         + w(r165f) &
         + w(r177f) 

    cdot(sC3H6) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         + w(r12f) &
         + w(r14f) &
         + w(r20f) &
         - w(r20b) &
         + w(r23f) &
         - w(r23b) &
         - w(r32f) &
         + w(r33f) &
         + 0.8_pr * w(r34f) &
         - w(r36f) &
         - w(r37f) &
         + 0.85_pr * w(r39f) &
         + 0.2_pr * w(r40f) &
         - w(r41f) &
         - w(r43f) &
         + w(r43b) &
         + w(r44f) &
         - w(r44b) &
         - w(r47f) &
         + w(r47b) &
         - w(r67f) &
         - w(r72f) &
         + 0.2_pr * w(r74f) &
         - w(r75f) &
         + 0.67_pr * w(r77f) &
         - w(r78f) &
         + w(r86f) &
         - w(r86b) &
         + 0.45_pr * w(r88f) &
         + w(r89f) &
         + w(r96f) &
         + w(r98f) &
         + w(r105f) &
         + w(r107f) &
         - w(r111f) &
         - w(r112f) &
         - w(r113f) &
         - w(r114f) &
         - w(r115f) &
         - w(r116f) &
         - w(r117f) &
         - w(r119f) &
         - w(r120f) &
         - w(r121f) &
         - w(r122f) &
         - w(r124f) &
         + w(r127f) &
         + w(r131f) &
         + w(r139f) &
         + w(r144f) &
         + w(r148f) &
         + w(r153f) &
         + w(r158f) &
         + w(r161f) &
         + w(r168f) &
         + w(r173f) &
         + w(r180f) 

    cdot(sC3H8) = 0.0_pr &
         + w(r4f) &
         - w(r4b) &
         - w(r44f) &
         + w(r44b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         + w(r99f) &
         + w(r108f) &
         + w(r113f) &
         - w(r125f) &
         - w(r126f) &
         - w(r127f) &
         - w(r128f) &
         - w(r129f) &
         - w(r130f) &
         - w(r131f) &
         + w(r162f) &
         + w(r163f) &
         + w(r174f) &
         + w(r175f) 

    cdot(sC4H6) = 0.0_pr &
         + w(r7f) &
         - w(r7b) &
         + w(r9f) &
         + w(r13f) &
         + w(r14f) &
         + w(r15f) &
         + w(r24f) &
         - w(r24b) &
         + w(r27f) &
         - w(r27b) &
         - w(r30f) &
         + 0.2_pr * w(r36f) &
         + w(r38f) &
         + 0.15_pr * w(r39f) &
         + 0.8_pr * w(r40f) &
         + w(r41f) &
         + w(r45f) &
         - w(r45b) &
         - w(r46f) &
         + w(r60f) &
         - w(r68f) &
         - w(r70f) &
         + 0.5_pr * w(r74f) &
         + 0.6_pr * w(r75f) &
         + 0.5_pr * w(r78f) &
         + 0.65_pr * w(r82f) &
         + w(r94f) &
         + w(r101f) &
         + w(r110f) &
         + w(r116f) &
         - w(r132f) &
         - w(r133f) &
         - w(r134f) &
         - w(r135f) &
         + w(r166f) &
         + w(r178f) 

    cdot(sNC4H8) = 0.0_pr &
         - w(r5f) &
         + w(r5b) &
         + w(r15f) &
         + 0.5_pr * w(r32f) &
         + 0.2_pr * w(r34f) &
         + 0.7_pr * w(r37f) &
         - w(r45f) &
         + w(r45b) &
         + 0.35_pr * w(r78f) &
         + w(r84f) &
         - w(r84b) &
         + w(r87f) &
         - w(r87b) &
         - w(r88f) &
         + w(r90f) &
         - w(r90b) &
         + w(r95f) &
         + w(r117f) &
         - w(r136f) &
         - w(r137f) &
         - w(r138f) &
         - w(r139f) &
         - w(r140f) &
         - w(r141f) &
         - w(r142f) &
         - w(r143f) &
         - w(r144f) &
         - w(r145f) &
         + w(r167f) &
         + w(r169f) &
         + w(r179f) &
         + w(r181f) 

    cdot(sNC4H10) = 0.0_pr &
         - w(r79f) &
         + w(r79b) &
         - w(r90f) &
         + w(r90b) &
         - w(r91f) &
         + w(r91b) &
         - w(r92f) &
         + w(r92b) &
         + w(r119f) &
         - w(r159f) &
         - w(r160f) &
         - w(r161f) &
         - w(r162f) &
         - w(r163f) &
         - w(r164f) &
         - w(r165f) &
         - w(r166f) &
         - w(r167f) &
         - w(r168f) &
         - w(r169f) &
         - w(r171f) &
         - w(r172f) &
         - w(r173f) &
         - w(r174f) &
         - w(r175f) &
         - w(r176f) &
         - w(r177f) &
         - w(r178f) &
         - w(r179f) &
         - w(r180f) &
         - w(r181f) 

    cdot(sCYC5H6) = 0.0_pr &
         + 0.3_pr * w(r37f) &
         - w(r54f) &
         + w(r54b) &
         + w(r55f) &
         - w(r55b) &
         - w(r56f) &
         + w(r57f) &
         - w(r58f) &
         + w(r61f) &
         + w(r62f) &
         + w(r64f) &
         - w(r66f) &
         - w(r67f) &
         + 0.85_pr * w(r68f) &
         + 0.7_pr * w(r70f) &
         + 0.5_pr * w(r72f) &
         + w(r76f) &
         - w(r146f) &
         - w(r147f) &
         - w(r148f) &
         - w(r149f) 

    cdot(sC5H8) = 0.0_pr &
         + 0.95_pr * w(r73f) &
         - w(r74f) &
         + 0.4_pr * w(r75f) &
         - w(r76f) &
         + 0.33_pr * w(r77f) &
         + 0.15_pr * w(r78f) &
         + 0.35_pr * w(r82f) 

    cdot(sNC5H10) = 0.0_pr &
         - w(r80f) &
         + w(r80b) &
         - w(r81f) &
         + w(r81b) &
         - w(r155f) &
         - w(r156f) &
         - w(r157f) &
         - w(r158f) 

    cdot(sC6H6) = 0.0_pr &
         + w(r46f) &
         + w(r58f) &
         + w(r66f) &
         + w(r67f) &
         + 0.15_pr * w(r68f) &
         + w(r69f) &
         + 0.3_pr * w(r70f) &
         + w(r71f) &
         + 0.5_pr * w(r72f) &
         + w(r150f) &
         + w(r151f) &
         + w(r152f) &
         + w(r153f) &
         + w(r154f) 

    cdot(sMCPTD) = 0.0_pr &
         + w(r63f) &
         - w(r64f) &
         + w(r65f) &
         - w(r150f) &
         - w(r151f) &
         - w(r152f) &
         - w(r153f) &
         - w(r154f) 

    cdot(sH) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r3b) &
         - w(r16f) &
         + w(r16b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         + w(r20f) &
         - w(r20b) &
         + w(r23f) &
         - w(r23b) &
         + w(r24f) &
         - w(r24b) &
         + w(r27f) &
         - w(r27b) &
         + w(r28f) &
         - w(r28b) &
         + w(r29f) &
         - w(r47f) &
         + w(r47b) &
         - w(r48f) &
         + w(r48b) &
         - w(r50f) &
         + w(r50b) &
         - w(r51f) &
         + w(r51b) &
         - w(r52f) &
         + w(r52b) &
         - w(r53f) &
         + w(r53b) &
         + w(r54f) &
         - w(r54b) &
         - w(r56f) &
         + w(r57f) &
         + w(r61f) &
         - w(r64f) &
         + w(r66f) &
         + 0.15_pr * w(r68f) &
         + w(r69f) &
         + w(r71f) &
         - w(r74f) &
         + 0.4_pr * w(r75f) &
         + 0.33_pr * w(r77f) &
         + 0.35_pr * w(r82f) &
         + w(r84f) &
         - w(r84b) &
         + w(r87f) &
         - w(r87b) &
         - w(r91f) &
         + w(r91b) &
         - w(r92f) &
         + w(r92b) &
         + w(r93f) &
         + w(r94f) &
         + w(r95f) &
         + w(r96f) &
         - w(r102f) &
         - w(r120f) &
         - w(r132f) &
         - w(r136f) &
         - w(r141f) &
         + w(r151f) &
         + w(r152f) &
         + w(r153f) &
         + w(r154f) &
         - w(r155f) 

    cdot(sCH3) = 0.0_pr &
         -2.0_pr * w(r1f) &
         + 2.0_pr * w(r1b) &
         - w(r4f) &
         + w(r4b) &
         + w(r5f) &
         - w(r5b) &
         - w(r8f) &
         + w(r8b) &
         + w(r19f) &
         - w(r19b) &
         + w(r22f) &
         - w(r22b) &
         -2.0_pr * w(r28f) &
         + 2.0_pr * w(r28b) &
         - w(r29f) &
         - w(r30f) &
         + 0.5_pr * w(r32f) &
         + 0.2_pr * w(r34f) &
         + 0.3_pr * w(r37f) &
         + w(r38f) &
         - w(r49f) &
         + w(r49b) &
         + w(r50f) &
         - w(r50b) &
         - w(r59f) &
         - w(r60f) &
         + w(r62f) &
         - w(r63f) &
         + w(r64f) &
         + w(r66f) &
         + 2.0_pr * w(r67f) &
         + 0.85_pr * w(r68f) &
         + 0.3_pr * w(r70f) &
         + w(r72f) &
         + 0.5_pr * w(r74f) &
         + 0.6_pr * w(r75f) &
         + 0.15_pr * w(r78f) &
         + w(r81f) &
         - w(r81b) &
         + 0.65_pr * w(r82f) &
         + w(r86f) &
         - w(r86b) &
         + w(r97f) &
         + w(r98f) &
         + w(r99f) &
         + w(r100f) &
         + w(r101f) &
         - w(r103f) &
         - w(r106f) &
         - w(r111f) &
         - w(r121f) &
         - w(r125f) &
         - w(r129f) &
         - w(r133f) &
         - w(r137f) &
         - w(r142f) &
         - w(r146f) &
         - w(r151f) &
         - w(r156f) &
         - w(r159f) &
         - w(r171f) 

    cdot(sC2H3) = 0.0_pr &
         + w(r6f) &
         - w(r6b) &
         - w(r8f) &
         + w(r8b) &
         + w(r10f) &
         - w(r10b) &
         + w(r16f) &
         - w(r16b) &
         + w(r26f) &
         - w(r26b) &
         + 0.7_pr * w(r37f) &
         + w(r59f) &
         - w(r61f) &
         - w(r68f) &
         - w(r73f) &
         + 0.2_pr * w(r74f) &
         - w(r75f) &
         - w(r76f) &
         + 0.67_pr * w(r77f) &
         + 0.35_pr * w(r78f) &
         - w(r93f) &
         - w(r100f) &
         + w(r102f) &
         + w(r103f) &
         + w(r104f) &
         + w(r105f) &
         - w(r109f) &
         - w(r114f) &
         - w(r124f) &
         - w(r128f) &
         - w(r135f) &
         - w(r140f) &
         - w(r145f) &
         - w(r149f) &
         - w(r154f) &
         - w(r164f) &
         - w(r176f) 

    cdot(sC2H5) = 0.0_pr &
         - w(r4f) &
         + w(r4b) &
         - w(r9f) &
         + w(r11f) &
         - w(r11b) &
         + w(r17f) &
         - w(r17b) &
         + w(r28f) &
         - w(r28b) &
         - w(r31f) &
         - w(r32f) &
         + w(r33f) &
         + 0.8_pr * w(r34f) &
         + 0.2_pr * w(r36f) &
         + 0.15_pr * w(r39f) &
         + 0.8_pr * w(r40f) &
         + w(r51f) &
         - w(r51b) &
         + 0.7_pr * w(r70f) &
         + w(r76f) &
         + 0.5_pr * w(r78f) &
         + 2.0_pr * w(r79f) &
         - 2.0_pr * w(r79b) &
         + w(r80f) &
         - w(r80b) &
         + w(r83f) &
         - w(r83b) &
         - w(r88f) &
         - w(r97f) &
         - w(r104f) &
         + w(r106f) &
         + w(r107f) &
         + w(r108f) &
         + w(r109f) &
         + w(r110f) &
         - w(r112f) &
         - w(r122f) &
         - w(r126f) &
         - w(r130f) &
         - w(r134f) &
         - w(r138f) &
         - w(r143f) &
         - w(r147f) &
         - w(r152f) &
         - w(r157f) &
         - w(r160f) &
         - w(r172f) 

    cdot(sC3H3) = 0.0_pr &
         - w(r13f) &
         + w(r48f) &
         - w(r48b) &
         + w(r49f) &
         - w(r49b) &
         - w(r115f) &
         - w(r165f) &
         - w(r177f) 

    cdot(sCHCHCH3) = 0.0_pr &
         - w(r19f) &
         + w(r19b) &
         - w(r38f) &
         + 0.1_pr * w(r40f) &
         - w(r65f) &
         + 0.3_pr * w(r74f) &
         - w(r78f) &
         - w(r96f) &
         - w(r118f) &
         + w(r120f) &
         + w(r121f) &
         + w(r122f) &
         + w(r123f) &
         + w(r124f) &
         - w(r168f) &
         - w(r180f) 

    cdot(sCH2CHCH2) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r3b) &
         + w(r5f) &
         - w(r5b) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         -2.0_pr * w(r12f) &
         - w(r13f) &
         - w(r14f) &
         + w(r18f) &
         - w(r18b) &
         + w(r30f) &
         - w(r36f) &
         - w(r37f) &
         + 0.85_pr * w(r39f) &
         + 0.1_pr * w(r40f) &
         + w(r47f) &
         - w(r47b) &
         - w(r65f) &
         - w(r69f) &
         - w(r70f) &
         - w(r73f) &
         - w(r77f) &
         + w(r80f) &
         - w(r80b) &
         - w(r98f) &
         - w(r105f) &
         - w(r107f) &
         + w(r111f) &
         + w(r112f) &
         + w(r113f) &
         + w(r114f) &
         + w(r115f) &
         + w(r116f) &
         + w(r117f) &
         + w(r118f) &
         + w(r119f) &
         - w(r123f) &
         - w(r127f) &
         - w(r131f) &
         - w(r139f) &
         - w(r144f) &
         - w(r148f) &
         - w(r153f) &
         - w(r158f) &
         - w(r161f) &
         - w(r173f) 

    cdot(sNC3H7) = 0.0_pr &
         + w(r21f) &
         - w(r21b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         + 0.5_pr * w(r32f) &
         - w(r33f) &
         + w(r35f) &
         + w(r41f) &
         + w(r52f) &
         - w(r52b) &
         + 0.45_pr * w(r88f) &
         + w(r89f) &
         + w(r125f) &
         + w(r126f) &
         + w(r127f) &
         + w(r128f) &
         - w(r162f) &
         - w(r174f) 

    cdot(sIC3H7) = 0.0_pr &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         - w(r34f) &
         - w(r35f) &
         + w(r53f) &
         - w(r53b) &
         - w(r99f) &
         - w(r108f) &
         - w(r113f) &
         + w(r129f) &
         + w(r130f) &
         + w(r131f) &
         - w(r163f) &
         - w(r175f) 

    cdot(sC4H5) = 0.0_pr &
         - w(r6f) &
         + w(r6b) &
         + w(r58f) &
         + w(r59f) &
         - w(r71f) &
         - w(r94f) &
         - w(r101f) &
         - w(r110f) &
         - w(r116f) &
         + w(r132f) &
         + w(r133f) &
         + w(r134f) &
         + w(r135f) &
         - w(r166f) &
         - w(r178f) 

    cdot(sCH2C3H5) = 0.0_pr &
         - w(r25f) &
         + w(r25b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         + w(r31f) &
         + 0.8_pr * w(r36f) &
         - w(r40f) &
         - w(r41f) &
         + w(r141f) &
         + w(r142f) &
         + w(r143f) &
         + w(r144f) &
         + w(r145f) &
         - w(r169f) &
         - w(r181f) 

    cdot(sSC4H7) = 0.0_pr &
         - w(r9f) &
         - w(r14f) &
         -2.0_pr * w(r15f) &
         - w(r24f) &
         + w(r24b) &
         + w(r25f) &
         - w(r25b) &
         - w(r39f) &
         - w(r62f) &
         - w(r72f) &
         + w(r81f) &
         - w(r81b) &
         - w(r95f) &
         - w(r117f) &
         + w(r136f) &
         + w(r137f) &
         + w(r138f) &
         + w(r139f) &
         + w(r140f) &
         - w(r167f) &
         - w(r179f) 

    cdot(sNC4H9S) = 0.0_pr &
         + w(r85f) &
         - w(r85b) &
         - w(r86f) &
         + w(r86b) &
         - w(r87f) &
         + w(r87b) &
         + 0.32_pr * w(r88f) &
         + w(r92f) &
         - w(r92b) &
         - w(r119f) &
         - w(r170f) &
         + w(r171f) &
         + w(r172f) &
         + w(r173f) &
         + w(r174f) &
         + w(r175f) &
         + w(r176f) &
         + w(r177f) &
         + w(r178f) &
         + w(r179f) &
         + w(r180f) &
         + w(r181f) &
         + w(r182f) 

    cdot(sNC4H9P) = 0.0_pr &
         - w(r83f) &
         + w(r83b) &
         - w(r84f) &
         + w(r84b) &
         - w(r85f) &
         + w(r85b) &
         + 0.23_pr * w(r88f) &
         - w(r89f) &
         + w(r91f) &
         - w(r91b) &
         + w(r159f) &
         + w(r160f) &
         + w(r161f) &
         + w(r162f) &
         + w(r163f) &
         + w(r164f) &
         + w(r165f) &
         + w(r166f) &
         + w(r167f) &
         + w(r168f) &
         + w(r169f) &
         + w(r170f) &
         - w(r182f) 

    cdot(sCYC5H5) = 0.0_pr &
         + w(r54f) &
         - w(r54b) &
         + w(r56f) &
         - w(r57f) &
         - w(r58f) &
         - w(r59f) &
         - w(r60f) &
         - w(r63f) &
         + w(r146f) &
         + w(r147f) &
         + w(r148f) &
         + w(r149f) 

    cdot(sNC5H9X3) = 0.0_pr &
         - w(r82f) &
         + w(r155f) &
         + w(r156f) &
         + w(r157f) &
         + w(r158f) 

  
    return
  end subroutine get_production_rates


  ! ----------------------------------------------- !
  ! 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
    !  c = max(c, tiny(1.0_pr))

        return
    end subroutine y2c

end module mod_S33R182QSS14
! ----------------------------------------------- !
! Cantera routine                                 !
! ----------------------------------------------- !
subroutine S33R182QSS14(P, T, y, wdot)

  use mod_S33R182QSS14
  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_rate_coefficients(k, M, T, P)
  call get_reaction_rates(w, k, M, c, cqss)

  ! Evaluate production rates
  call get_production_rates(wdot,w)

  return
end subroutine S33R182QSS14