!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE C3H8_SanDiego_ARC21_153_11_QM.f90
!>    Module for calculating the analytical source terms for the C3H8_SanDiego_ARC21_153_11_QM scheme.
!!    @file    C3H8_SanDiego_ARC21_153_11_QM.f90
!!    @authors Q. Male
!!    @date    22/01/2018
!!    @since   
!!    @note    
!--------------------------------------------------------------------------------------------------

!--------------------------------------------------------------------------------------------------
!     MODULE mod_C3H8_SanDiego_ARC21_153_11_QM
!>    Module for calculating the analytical source terms for the C3H8_SanDiego_ARC21_153_11_QM scheme.
!!    @authors Q. Male
!!    @date    22/01/2018
!--------------------------------------------------------------------------------------------------
module mod_C3H8_SanDiego_ARC21_153_11_QM
      USE mod_param_defs
      implicit none

      integer, parameter :: nspec = 21
      integer, parameter :: nreac = 153
      integer, parameter :: isc_T = 1
      integer, parameter :: neq = nspec + 1

      ! QSS variables
      ! Number of QSS species
      integer, parameter :: nqss = 11
      ! QSS species
      integer, dimension(nspec + nqss) :: iqss

      real(pr), dimension(nspec + nqss) :: W_sp,Cp_sp,h_sp,dh_sp
      character(len=15), dimension(nspec + nqss) :: gname
      character(len=5), dimension(nreac) :: rname

      ! Post processing variables
      ! Number of groups
      integer, parameter :: ng = 33
      ! Max number of species in groups
      integer, parameter :: maxppn = 1
      ! Number of species in each group
      integer, dimension(ng) :: ppn 
      ! Species in each group
      integer, dimension(ng,maxppn) :: pp 
      ! Name of species in each group
      character(len=30), dimension(ng) :: ppname 
      ! Actual expression of each reaction
      character(len=65), dimension(nreac) :: reacexp
      ! Link between backward and forward rates
      integer, dimension(nreac) :: fofb

      ! Index of species
      integer, parameter :: sN2 = 1
      integer, parameter :: sH2 = 2
      integer, parameter :: sH = 3
      integer, parameter :: sOH = 4
      integer, parameter :: sH2O = 5
      integer, parameter :: sHO2 = 6
      integer, parameter :: sH2O2 = 7
      integer, parameter :: sO2 = 8
      integer, parameter :: sCH2O = 9
      integer, parameter :: sCO = 10
      integer, parameter :: sCH3 = 11
      integer, parameter :: sCH4 = 12
      integer, parameter :: sC2H6 = 13
      integer, parameter :: sCO2 = 14
      integer, parameter :: sC2H2 = 15
      integer, parameter :: sCH2CO = 16
      integer, parameter :: sC3H6 = 17
      integer, parameter :: sC2H4 = 18
      integer, parameter :: sC2H4O = 19
      integer, parameter :: sC3H8 = 20
      integer, parameter :: sOC3H5OOH = 21
      integer, parameter :: sO = 22
      integer, parameter :: sSXCH2 = 23
      integer, parameter :: sC2H5 = 24
      integer, parameter :: sCH3O = 25
      integer, parameter :: sHCO = 26
      integer, parameter :: sCH2OH = 27
      integer, parameter :: sC2H3 = 28
      integer, parameter :: sCH2CHO = 29
      integer, parameter :: sNXC3H7 = 30
      integer, parameter :: sIXC3H7 = 31
      integer, parameter :: sC3H6OOH = 32

      ! Index of reactions
      integer, parameter :: r1 = 1
      integer, parameter :: r2f = 2
      integer, parameter :: r3 = 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 :: r11 = 11
      integer, parameter :: r12 = 12
      integer, parameter :: r13f = 13
      integer, parameter :: r14f = 14
      integer, parameter :: r15 = 15
      integer, parameter :: r16 = 16
      integer, parameter :: r17 = 17
      integer, parameter :: r18 = 18
      integer, parameter :: r19 = 19
      integer, parameter :: r20 = 20
      integer, parameter :: r21 = 21
      integer, parameter :: r22 = 22
      integer, parameter :: r23f = 23
      integer, parameter :: r24 = 24
      integer, parameter :: r25f = 25
      integer, parameter :: r26f = 26
      integer, parameter :: r27 = 27
      integer, parameter :: r28f = 28
      integer, parameter :: r29 = 29
      integer, parameter :: r30 = 30
      integer, parameter :: r31 = 31
      integer, parameter :: r32 = 32
      integer, parameter :: r33f = 33
      integer, parameter :: r34f = 34
      integer, parameter :: r35f = 35
      integer, parameter :: r36 = 36
      integer, parameter :: r37 = 37
      integer, parameter :: r38f = 38
      integer, parameter :: r39 = 39
      integer, parameter :: r40f = 40
      integer, parameter :: r41 = 41
      integer, parameter :: r42 = 42
      integer, parameter :: r43 = 43
      integer, parameter :: r44 = 44
      integer, parameter :: r45 = 45
      integer, parameter :: r46 = 46
      integer, parameter :: r47 = 47
      integer, parameter :: r48 = 48
      integer, parameter :: r49 = 49
      integer, parameter :: r50f = 50
      integer, parameter :: r51 = 51
      integer, parameter :: r52 = 52
      integer, parameter :: r53 = 53
      integer, parameter :: r54 = 54
      integer, parameter :: r55 = 55
      integer, parameter :: r56f = 56
      integer, parameter :: r57f = 57
      integer, parameter :: r58 = 58
      integer, parameter :: r59 = 59
      integer, parameter :: r60f = 60
      integer, parameter :: r61 = 61
      integer, parameter :: r62 = 62
      integer, parameter :: r63f = 63
      integer, parameter :: r64 = 64
      integer, parameter :: r65f = 65
      integer, parameter :: r66 = 66
      integer, parameter :: r67 = 67
      integer, parameter :: r68 = 68
      integer, parameter :: r69 = 69
      integer, parameter :: r70 = 70
      integer, parameter :: r71 = 71
      integer, parameter :: r72f = 72
      integer, parameter :: r73f = 73
      integer, parameter :: r74 = 74
      integer, parameter :: r75 = 75
      integer, parameter :: r76 = 76
      integer, parameter :: r77f = 77
      integer, parameter :: r80 = 78
      integer, parameter :: r81f = 79
      integer, parameter :: r82 = 80
      integer, parameter :: r83 = 81
      integer, parameter :: r84 = 82
      integer, parameter :: r85f = 83
      integer, parameter :: r86 = 84
      integer, parameter :: r87 = 85
      integer, parameter :: r88 = 86
      integer, parameter :: r89 = 87
      integer, parameter :: r90f = 88
      integer, parameter :: r91 = 89
      integer, parameter :: r92 = 90
      integer, parameter :: r94 = 91
      integer, parameter :: r95f = 92
      integer, parameter :: r96 = 93
      integer, parameter :: r97f = 94
      integer, parameter :: r98 = 95
      integer, parameter :: r99 = 96
      integer, parameter :: r100f = 97
      integer, parameter :: r101 = 98
      integer, parameter :: r102f = 99
      integer, parameter :: r103f = 100
      integer, parameter :: r104 = 101
      integer, parameter :: r105 = 102
      integer, parameter :: r106f = 103
      integer, parameter :: r107f = 104
      integer, parameter :: r108 = 105
      integer, parameter :: r109 = 106
      integer, parameter :: r110f = 107
      integer, parameter :: r111 = 108
      integer, parameter :: r112 = 109
      integer, parameter :: r113f = 110
      integer, parameter :: r114 = 111
      integer, parameter :: r115 = 112
      integer, parameter :: r116 = 113
      integer, parameter :: r2b = 114
      integer, parameter :: r4b = 115
      integer, parameter :: r5b = 116
      integer, parameter :: r6b = 117
      integer, parameter :: r7b = 118
      integer, parameter :: r8b = 119
      integer, parameter :: r9b = 120
      integer, parameter :: r10b = 121
      integer, parameter :: r13b = 122
      integer, parameter :: r14b = 123
      integer, parameter :: r23b = 124
      integer, parameter :: r25b = 125
      integer, parameter :: r26b = 126
      integer, parameter :: r28b = 127
      integer, parameter :: r33b = 128
      integer, parameter :: r34b = 129
      integer, parameter :: r35b = 130
      integer, parameter :: r38b = 131
      integer, parameter :: r40b = 132
      integer, parameter :: r50b = 133
      integer, parameter :: r56b = 134
      integer, parameter :: r57b = 135
      integer, parameter :: r60b = 136
      integer, parameter :: r63b = 137
      integer, parameter :: r65b = 138
      integer, parameter :: r72b = 139
      integer, parameter :: r73b = 140
      integer, parameter :: r77b = 141
      integer, parameter :: r81b = 142
      integer, parameter :: r85b = 143
      integer, parameter :: r90b = 144
      integer, parameter :: r95b = 145
      integer, parameter :: r97b = 146
      integer, parameter :: r100b = 147
      integer, parameter :: r102b = 148
      integer, parameter :: r103b = 149
      integer, parameter :: r106b = 150
      integer, parameter :: r107b = 151
      integer, parameter :: r110b = 152
      integer, parameter :: r113b = 153

      ! Index of third bodies
      integer, parameter :: mM32 = 1
      integer, parameter :: mM5 = 2
      integer, parameter :: mM12 = 3
      integer, parameter :: mM16 = 4
      integer, parameter :: mM10 = 5
      integer, parameter :: mM4 = 6
      integer, parameter :: mM20 = 7
      integer, parameter :: mM9 = 8
      integer, parameter :: mM30 = 9
      integer, parameter :: mM7 = 10
      integer, parameter :: mM14 = 11
      integer, parameter :: mM19 = 12
      integer, parameter :: mM21 = 13
      integer, parameter :: mM1 = 14
      integer, parameter :: mM17 = 15
      integer, parameter :: mM2 = 16
      integer, parameter :: mM8 = 17
      integer, parameter :: mM18 = 18
      integer, parameter :: mM0 = 19
      integer, parameter :: mM6 = 20
      integer, parameter :: mM11 = 21
      integer, parameter :: mM15 = 22
      integer, parameter :: mM31 = 23

      contains

      ! Name of mechanism
      ! subroutine which_mechanism
      ! use parallel
      !implicit none

      ! if (irank.eq.iroot) print*, 'Using mechanism ./Spec_32_Reac_156_QSS_11.mech'

      ! return
      ! end subroutine which_mechanism

      ! Subroutine to define groups for post-processing
      subroutine pp_data
      implicit none

      ! Number of species in each group
      ppn(1) = 1
      ppn(2) = 1
      ppn(3) = 1
      ppn(4) = 1
      ppn(5) = 1
      ppn(6) = 1
      ppn(7) = 1
      ppn(8) = 1
      ppn(9) = 1
      ppn(10) = 1
      ppn(11) = 1
      ppn(12) = 1
      ppn(13) = 1
      ppn(14) = 1
      ppn(15) = 1
      ppn(16) = 1
      ppn(17) = 1
      ppn(18) = 1
      ppn(19) = 1
      ppn(20) = 1
      ppn(21) = 1
      ppn(22) = 1
      ppn(23) = 1
      ppn(24) = 1
      ppn(25) = 1
      ppn(26) = 1
      ppn(27) = 1
      ppn(28) = 1
      ppn(29) = 1
      ppn(30) = 1
      ppn(31) = 1
      ppn(32) = 1
      ppn(33) = 1

      ! Indices of species in each group
      pp(1,1) = 1
      pp(2,1) = 2
      pp(3,1) = 3
      pp(4,1) = 4
      pp(5,1) = 5
      pp(6,1) = 6
      pp(7,1) = 7
      pp(8,1) = 8
      pp(9,1) = 9
      pp(10,1) = 10
      pp(11,1) = 11
      pp(12,1) = 12
      pp(13,1) = 13
      pp(14,1) = 14
      pp(15,1) = 15
      pp(16,1) = 16
      pp(17,1) = 17
      pp(18,1) = 18
      pp(19,1) = 19
      pp(20,1) = 20
      pp(21,1) = 21
      pp(22,1) = 22
      pp(23,1) = 23
      pp(24,1) = 24
      pp(25,1) = 25
      pp(26,1) = 26
      pp(27,1) = 27
      pp(28,1) = 28
      pp(29,1) = 29
      pp(30,1) = 30
      pp(31,1) = 31
      pp(32,1) = 32
      pp(33,1) = 1

      ! Name of group of species
      ppname(1) = trim(gname(sN2))
      ppname(2) = trim(gname(sH2))
      ppname(3) = trim(gname(sH))
      ppname(4) = trim(gname(sO))
      ppname(5) = trim(gname(sOH))
      ppname(6) = trim(gname(sH2O))
      ppname(7) = trim(gname(sHO2))
      ppname(8) = trim(gname(sH2O2))
      ppname(9) = trim(gname(sO2))
      ppname(10) = trim(gname(sSXCH2))
      ppname(11) = trim(gname(sCH2O))
      ppname(12) = trim(gname(sCO))
      ppname(13) = trim(gname(sCH3))
      ppname(14) = trim(gname(sCH4))
      ppname(15) = trim(gname(sC2H5))
      ppname(16) = trim(gname(sC2H6))
      ppname(17) = trim(gname(sCH3O))
      ppname(18) = trim(gname(sCO2))
      ppname(19) = trim(gname(sHCO))
      ppname(20) = trim(gname(sCH2OH))
      ppname(21) = trim(gname(sC2H2))
      ppname(22) = trim(gname(sCH2CO))
      ppname(23) = trim(gname(sC2H3))
      ppname(24) = trim(gname(sCH2CHO))
      ppname(25) = trim(gname(sC3H6))
      ppname(26) = trim(gname(sC2H4))
      ppname(27) = trim(gname(sC2H4O))
      ppname(28) = trim(gname(sNXC3H7))
      ppname(29) = trim(gname(sIXC3H7))
      ppname(30) = trim(gname(sC3H6OOH))
      ppname(31) = trim(gname(sC3H8))
      ppname(32) = trim(gname(sOC3H5OOH))
      ppname(33) = 'N2X'

      end subroutine pp_data

      ! Molar mass
      subroutine molar_mass
      implicit none

      W_sp(sN2) = 0.02802_pr
      W_sp(sH2) = 0.002016_pr
      W_sp(sH) = 0.001008_pr
      W_sp(sO) = 0.016_pr
      W_sp(sOH) = 0.017008_pr
      W_sp(sH2O) = 0.018016_pr
      W_sp(sHO2) = 0.033008_pr
      W_sp(sH2O2) = 0.034016_pr
      W_sp(sO2) = 0.032_pr
      W_sp(sSXCH2) = 0.014026_pr
      W_sp(sCH2O) = 0.030026_pr
      W_sp(sCO) = 0.02801_pr
      W_sp(sCH3) = 0.015034_pr
      W_sp(sCH4) = 0.016042_pr
      W_sp(sC2H5) = 0.02906_pr
      W_sp(sC2H6) = 0.030068_pr
      W_sp(sCH3O) = 0.031034_pr
      W_sp(sCO2) = 0.04401_pr
      W_sp(sHCO) = 0.029018_pr
      W_sp(sCH2OH) = 0.031034_pr
      W_sp(sC2H2) = 0.026036_pr
      W_sp(sCH2CO) = 0.042036_pr
      W_sp(sC2H3) = 0.027044_pr
      W_sp(sCH2CHO) = 0.043044_pr
      W_sp(sC3H6) = 0.042078_pr
      W_sp(sC2H4) = 0.028052_pr
      W_sp(sC2H4O) = 0.044052_pr
      W_sp(sNXC3H7) = 0.043086_pr
      W_sp(sIXC3H7) = 0.043086_pr
      W_sp(sC3H6OOH) = 0.075086_pr
      W_sp(sC3H8) = 0.044094_pr
      W_sp(sOC3H5OOH) = 0.090078_pr

      end subroutine molar_mass

      ! Species names
      subroutine species_name
      implicit none

      gname(sN2) = 'N2'
      gname(sH2) = 'H2'
      gname(sH) = 'H'
      gname(sO) = 'O'
      gname(sOH) = 'OH'
      gname(sH2O) = 'H2O'
      gname(sHO2) = 'HO2'
      gname(sH2O2) = 'H2O2'
      gname(sO2) = 'O2'
      gname(sSXCH2) = 'S-CH2'
      gname(sCH2O) = 'CH2O'
      gname(sCO) = 'CO'
      gname(sCH3) = 'CH3'
      gname(sCH4) = 'CH4'
      gname(sC2H5) = 'C2H5'
      gname(sC2H6) = 'C2H6'
      gname(sCH3O) = 'CH3O'
      gname(sCO2) = 'CO2'
      gname(sHCO) = 'HCO'
      gname(sCH2OH) = 'CH2OH'
      gname(sC2H2) = 'C2H2'
      gname(sCH2CO) = 'CH2CO'
      gname(sC2H3) = 'C2H3'
      gname(sCH2CHO) = 'CH2CHO'
      gname(sC3H6) = 'C3H6'
      gname(sC2H4) = 'C2H4'
      gname(sC2H4O) = 'C2H4O'
      gname(sNXC3H7) = 'N-C3H7'
      gname(sIXC3H7) = 'I-C3H7'
      gname(sC3H6OOH) = 'C3H6OOH'
      gname(sC3H8) = 'C3H8'
      gname(sOC3H5OOH) = 'OC3H5OOH'

      end subroutine species_name

      ! Reaction names
      subroutine reaction_name
      implicit none

      rname(r1) = '1'
      rname(r2f) = '2f'
      rname(r3) = '3'
      rname(r4f) = '4f'
      rname(r5f) = '5f'
      rname(r6f) = '6f'
      rname(r7f) = '7f'
      rname(r8f) = '8f'
      rname(r9f) = '9f'
      rname(r10f) = '10f'
      rname(r11) = '11'
      rname(r12) = '12'
      rname(r13f) = '13f'
      rname(r14f) = '14f'
      rname(r15) = '15'
      rname(r16) = '16'
      rname(r17) = '17'
      rname(r18) = '18'
      rname(r19) = '19'
      rname(r20) = '20'
      rname(r21) = '21'
      rname(r22) = '22'
      rname(r23f) = '23f'
      rname(r24) = '24'
      rname(r25f) = '25f'
      rname(r26f) = '26f'
      rname(r27) = '27'
      rname(r28f) = '28f'
      rname(r29) = '29'
      rname(r30) = '30'
      rname(r31) = '31'
      rname(r32) = '32'
      rname(r33f) = '33f'
      rname(r34f) = '34f'
      rname(r35f) = '35f'
      rname(r36) = '36'
      rname(r37) = '37'
      rname(r38f) = '38f'
      rname(r39) = '39'
      rname(r40f) = '40f'
      rname(r41) = '41'
      rname(r42) = '42'
      rname(r43) = '43'
      rname(r44) = '44'
      rname(r45) = '45'
      rname(r46) = '46'
      rname(r47) = '47'
      rname(r48) = '48'
      rname(r49) = '49'
      rname(r50f) = '50f'
      rname(r51) = '51'
      rname(r52) = '52'
      rname(r53) = '53'
      rname(r54) = '54'
      rname(r55) = '55'
      rname(r56f) = '56f'
      rname(r57f) = '57f'
      rname(r58) = '58'
      rname(r59) = '59'
      rname(r60f) = '60f'
      rname(r61) = '61'
      rname(r62) = '62'
      rname(r63f) = '63f'
      rname(r64) = '64'
      rname(r65f) = '65f'
      rname(r66) = '66'
      rname(r67) = '67'
      rname(r68) = '68'
      rname(r69) = '69'
      rname(r70) = '70'
      rname(r71) = '71'
      rname(r72f) = '72f'
      rname(r73f) = '73f'
      rname(r74) = '74'
      rname(r75) = '75'
      rname(r76) = '76'
      rname(r77f) = '77f'
      rname(r80) = '80'
      rname(r81f) = '81f'
      rname(r82) = '82'
      rname(r83) = '83'
      rname(r84) = '84'
      rname(r85f) = '85f'
      rname(r86) = '86'
      rname(r87) = '87'
      rname(r88) = '88'
      rname(r89) = '89'
      rname(r90f) = '90f'
      rname(r91) = '91'
      rname(r92) = '92'
      rname(r94) = '94'
      rname(r95f) = '95f'
      rname(r96) = '96'
      rname(r97f) = '97f'
      rname(r98) = '98'
      rname(r99) = '99'
      rname(r100f) = '100f'
      rname(r101) = '101'
      rname(r102f) = '102f'
      rname(r103f) = '103f'
      rname(r104) = '104'
      rname(r105) = '105'
      rname(r106f) = '106f'
      rname(r107f) = '107f'
      rname(r108) = '108'
      rname(r109) = '109'
      rname(r110f) = '110f'
      rname(r111) = '111'
      rname(r112) = '112'
      rname(r113f) = '113f'
      rname(r114) = '114'
      rname(r115) = '115'
      rname(r116) = '116'
      rname(r2b) = '2b'
      rname(r4b) = '4b'
      rname(r5b) = '5b'
      rname(r6b) = '6b'
      rname(r7b) = '7b'
      rname(r8b) = '8b'
      rname(r9b) = '9b'
      rname(r10b) = '10b'
      rname(r13b) = '13b'
      rname(r14b) = '14b'
      rname(r23b) = '23b'
      rname(r25b) = '25b'
      rname(r26b) = '26b'
      rname(r28b) = '28b'
      rname(r33b) = '33b'
      rname(r34b) = '34b'
      rname(r35b) = '35b'
      rname(r38b) = '38b'
      rname(r40b) = '40b'
      rname(r50b) = '50b'
      rname(r56b) = '56b'
      rname(r57b) = '57b'
      rname(r60b) = '60b'
      rname(r63b) = '63b'
      rname(r65b) = '65b'
      rname(r72b) = '72b'
      rname(r73b) = '73b'
      rname(r77b) = '77b'
      rname(r81b) = '81b'
      rname(r85b) = '85b'
      rname(r90b) = '90b'
      rname(r95b) = '95b'
      rname(r97b) = '97b'
      rname(r100b) = '100b'
      rname(r102b) = '102b'
      rname(r103b) = '103b'
      rname(r106b) = '106b'
      rname(r107b) = '107b'
      rname(r110b) = '110b'
      rname(r113b) = '113b'

      end subroutine reaction_name

      ! List of QSS species
      subroutine QSS_list
      implicit none

      iqss(sN2) = 0
      iqss(sH2) = 0
      iqss(sH) = 0
      iqss(sO) = 1
      iqss(sOH) = 0
      iqss(sH2O) = 0
      iqss(sHO2) = 0
      iqss(sH2O2) = 0
      iqss(sO2) = 0
      iqss(sSXCH2) = 1
      iqss(sCH2O) = 0
      iqss(sCO) = 0
      iqss(sCH3) = 0
      iqss(sCH4) = 0
      iqss(sC2H5) = 1
      iqss(sC2H6) = 0
      iqss(sCH3O) = 1
      iqss(sCO2) = 0
      iqss(sHCO) = 1
      iqss(sCH2OH) = 1
      iqss(sC2H2) = 0
      iqss(sCH2CO) = 0
      iqss(sC2H3) = 1
      iqss(sCH2CHO) = 1
      iqss(sC3H6) = 0
      iqss(sC2H4) = 0
      iqss(sC2H4O) = 0
      iqss(sNXC3H7) = 1
      iqss(sIXC3H7) = 1
      iqss(sC3H6OOH) = 1
      iqss(sC3H8) = 0
      iqss(sOC3H5OOH) = 0

      end subroutine QSS_list

      ! 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), parameter :: R = 8.31434_pr
      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 / ( R * Tloc )
      end if
      ntmp = 0.75_pr - 1.27_pr * dlog10( fc )

      ccoeff = - 0.4_pr - 0.67_pr * dlog10( fc )
      dcoeff = 0.14_pr
      k0 = k0 * conc / max(kinf, 1.0e-60_pr)
      lgknull = dlog10(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

      subroutine YtoC(c,P,T,Y)
      implicit none

      real(pr), dimension(nspec) :: c,Y
      real(pr) :: bla, P, T
      integer :: K

      call molar_mass

      c(sN2) = Y(sN2) / W_sp(sN2)
      c(sH2) = Y(sH2) / W_sp(sH2)
      c(sH) = Y(sH) / W_sp(sH)
      !c(sO) = Y(sO) / W_sp(sO)
      c(sOH) = Y(sOH) / W_sp(sOH)
      c(sH2O) = Y(sH2O) / W_sp(sH2O)
      c(sHO2) = Y(sHO2) / W_sp(sHO2)
      c(sH2O2) = Y(sH2O2) / W_sp(sH2O2)
      c(sO2) = Y(sO2) / W_sp(sO2)
      !c(sSXCH2) = Y(sSXCH2) / W_sp(sSXCH2)
      c(sCH2O) = Y(sCH2O) / W_sp(sCH2O)
      c(sCO) = Y(sCO) / W_sp(sCO)
      c(sCH3) = Y(sCH3) / W_sp(sCH3)
      c(sCH4) = Y(sCH4) / W_sp(sCH4)
      !c(sC2H5) = Y(sC2H5) / W_sp(sC2H5)
      c(sC2H6) = Y(sC2H6) / W_sp(sC2H6)
      !c(sCH3O) = Y(sCH3O) / W_sp(sCH3O)
      c(sCO2) = Y(sCO2) / W_sp(sCO2)
      !c(sHCO) = Y(sHCO) / W_sp(sHCO)
      !c(sCH2OH) = Y(sCH2OH) / W_sp(sCH2OH)
      c(sC2H2) = Y(sC2H2) / W_sp(sC2H2)
      c(sCH2CO) = Y(sCH2CO) / W_sp(sCH2CO)
      !c(sC2H3) = Y(sC2H3) / W_sp(sC2H3)
      !c(sCH2CHO) = Y(sCH2CHO) / W_sp(sCH2CHO)
      c(sC3H6) = Y(sC3H6) / W_sp(sC3H6)
      c(sC2H4) = Y(sC2H4) / W_sp(sC2H4)
      c(sC2H4O) = Y(sC2H4O) / W_sp(sC2H4O)
      !c(sNXC3H7) = Y(sNXC3H7) / W_sp(sNXC3H7)
      !c(sIXC3H7) = Y(sIXC3H7) / W_sp(sIXC3H7)
      !c(sC3H6OOH) = Y(sC3H6OOH) / W_sp(sC3H6OOH)
      c(sC3H8) = Y(sC3H8) / W_sp(sC3H8)
      c(sOC3H5OOH) = Y(sOC3H5OOH) / W_sp(sOC3H5OOH)

      bla = 0.0_pr
      DO K = 1, nspec
        bla  = bla + c(K)
      ENDDO
      bla = P/(bla*T*8.31451_pr)
      DO K = 1, nspec
        c(K) = max(c(K),1.0e-60_pr) * bla
      ENDDO
      end subroutine YtoC

! --- Thirdbodies --- !
      subroutine get_thirdbodies(M,c)
      implicit none

      real(pr), dimension(nspec) :: c
      real(pr), dimension(24) :: M

      M(mM32) = (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 
      M(mM5) = (1.5_pr)*c(sH2) &
                 + (2.8_pr)*c(sCO2) &
                 + (0.9_pr)*c(sCO) &
                 + (11_pr)*c(sH2O) &
                 + sum(c) 
      M(mM12) = (1_pr)*c(sCO2) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sH2) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sCH4) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 
      M(mM16) = (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + sum(c) 
      M(mM10) = (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (2_pr)*c(sC2H6) &
                 + sum(c) 
      M(mM4) = (1.5_pr)*c(sH2) &
                 + (2.8_pr)*c(sCO2) &
                 + (0.9_pr)*c(sCO) &
                 + (11_pr)*c(sH2O) &
                 + sum(c) 
      M(mM20) = (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM9) = (1.5_pr)*c(sCO) &
                 + (11_pr)*c(sH2O) &
                 + (0.9_pr)*c(sH2) &
                 + (1.5_pr)*c(sCO2) &
                 + sum(c) 
      M(mM30) = (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 
      M(mM7) = (1.5_pr)*c(sH2) &
                 + (5_pr)*c(sH2O2) &
                 + (1_pr)*c(sCO2) &
                 + (0.5_pr)*c(sCO) &
                 + (5_pr)*c(sH2O) &
                 + sum(c) 
      M(mM14) = (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + sum(c) 
      M(mM19) = (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM21) = (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM1) = (11_pr)*c(sH2O) &
                 + (0.9_pr)*c(sCO) &
                 + (2.8_pr)*c(sCO2) &
                 + (1.5_pr)*c(sH2) &
                 + sum(c) 
      M(mM17) = (1_pr)*c(sH2) &
                 + (1_pr)*c(sCO2) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + sum(c) 
      M(mM2) = (1.5_pr)*c(sH2) &
                 + (2.8_pr)*c(sCO2) &
                 + (0.9_pr)*c(sCO) &
                 + (11_pr)*c(sH2O) &
                 + sum(c) 
      M(mM8) = (3_pr)*c(sCO2) &
                 + (1.5_pr)*c(sH2) &
                 + (11_pr)*c(sH2O) &
                 + (1_pr)*c(sCO) &
                 + sum(c) 
      M(mM18) = (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + sum(c) 
      M(mM0) = sum(c) 
      M(mM6) = (0.5_pr)*c(sC2H6) &
                 + (1.5_pr)*c(sH2) &
                 + (1.4_pr)*c(sCO2) &
                 + (0.2_pr)*c(sCO) &
                 + (15_pr)*c(sH2O) &
                 + sum(c) 
      M(mM11) = (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sCH4) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 
      M(mM15) = (5_pr)*c(sH2O) &
                 + (1_pr)*c(sCH4) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (2_pr)*c(sC2H6) &
                 + sum(c) 
      M(mM31) = (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sCH4) &
                 + (5_pr)*c(sH2O) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 

      end subroutine get_thirdbodies

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

      real(pr), dimension(nreac) :: k
      real(pr), dimension(24) :: M
      real(pr) :: Tloc,Ploc

      real(pr) :: k6f_0, k6f_inf, FC6f
      real(pr) :: k10f_0, k10f_inf, FC10f
      real(pr) :: k26f_0, k26f_inf, FC26f
      real(pr) :: k33f_0, k33f_inf, FC33f
      real(pr) :: k37_0, k37_inf, FC37
      real(pr) :: k49_0, k49_inf, FC49
      real(pr) :: k63f_0, k63f_inf, FC63f
      real(pr) :: k65f_0, k65f_inf, FC65f
      real(pr) :: k74_0, k74_inf, FC74
      real(pr) :: k77f_0, k77f_inf, FC77f
      real(pr) :: k95f_0, k95f_inf, FC95f
      real(pr) :: k97f_0, k97f_inf, FC97f
      real(pr) :: k103f_0, k103f_inf, FC103f
      real(pr) :: k110f_0, k110f_inf, FC110f
      real(pr) :: k6b_0, k6b_inf, FC6b
      real(pr) :: k10b_0, k10b_inf, FC10b
      real(pr) :: k26b_0, k26b_inf, FC26b
      real(pr) :: k33b_0, k33b_inf, FC33b
      real(pr) :: k63b_0, k63b_inf, FC63b
      real(pr) :: k65b_0, k65b_inf, FC65b
      real(pr) :: k77b_0, k77b_inf, FC77b
      real(pr) :: k95b_0, k95b_inf, FC95b
      real(pr) :: k97b_0, k97b_inf, FC97b
      real(pr) :: k103b_0, k103b_inf, FC103b
      real(pr) :: k110b_0, k110b_inf, FC110b

      ! Rate coefficients

        k(r1) = (3.03900000e+11_pr)*Tloc**(-0.650_pr)*&
          exp(-(4.331e+05_pr)/(8.314_pr*Tloc)) 

        k(r2f) = (5.06000000e-02_pr)*Tloc**(2.670_pr)*&
          exp(-(2.632e+04_pr)/(8.314_pr*Tloc)) 

        k(r3) = (4.71000000e+06_pr)*Tloc**(-1.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r4f) = (4.00000000e+10_pr)*Tloc**(-2.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r5f) = (8.00000000e+03_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k6f_0 = (2.76000000e+13_pr)*Tloc**(-3.200_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        k6f_inf = (9.55000000e+07_pr)*Tloc**(-0.270_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC6f = (4.300e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (5.700e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r6f) =  &
          getlindratecoeff & 
        (Tloc, k6f_0, k6f_inf, FC6f, M(mM7), Ploc ) 

        k(r7f) = (1.17000000e+03_pr)*Tloc**(1.300_pr)*&
          exp(-(1.521e+04_pr)/(8.314_pr*Tloc)) 

        k(r8f) = (7.00000000e-01_pr)*Tloc**(2.330_pr)*&
          exp(-(6.087e+04_pr)/(8.314_pr*Tloc)) 

        k(r9f) = (3.52000000e+10_pr)*Tloc**(-0.700_pr)*&
          exp(-(7.142e+04_pr)/(8.314_pr*Tloc)) 

        k10f_0 = (5.75000000e+07_pr)*Tloc**(-1.400_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        k10f_inf = (4.65000000e+06_pr)*Tloc**(0.440_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC10f = (5.000e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (5.000e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r10f) =  &
          getlindratecoeff & 
        (Tloc, k10f_0, k10f_inf, FC10f, M(mM6), Ploc ) 

        k(r11) = (2.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r12) = (7.08000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.234e+03_pr)/(8.314_pr*Tloc)) 

        k(r13f) = (4.50000000e+08_pr)*Tloc**(0.000_pr)*&
          exp(-(4.573e+04_pr)/(8.314_pr*Tloc)) 

        k(r14f) = (7.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(-4.580e+03_pr)/(8.314_pr*Tloc)) 

        k(r15) = (3.10000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(7.200e+03_pr)/(8.314_pr*Tloc)) 

        k(r16) = (1.94000000e+05_pr)*Tloc**(0.000_pr)*&
          exp(-(-5.895e+03_pr)/(8.314_pr*Tloc)) 

        k(r17) = (1.03000000e+08_pr)*Tloc**(0.000_pr)*&
          exp(-(4.620e+04_pr)/(8.314_pr*Tloc)) 

        k(r18) = (1.66000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(3.443e+03_pr)/(8.314_pr*Tloc)) 

        k(r19) = (1.74000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(6.000e+03_pr)/(8.314_pr*Tloc)) 

        k(r20) = (7.59000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(3.043e+04_pr)/(8.314_pr*Tloc)) 

        k(r21) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r22) = (3.13000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r23f) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.047e+04_pr)/(8.314_pr*Tloc)) 

        k(r24) = (1.40800000e+05_pr)*Tloc**(0.149_pr)*&
          exp(-(5.430e+02_pr)/(8.314_pr*Tloc)) 

        k(r25f) = (3.16000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(6.150e+04_pr)/(8.314_pr*Tloc)) 

        k26f_0 = (1.27000000e+29_pr)*Tloc**(-7.000_pr)*&
          exp(-(1.156e+04_pr)/(8.314_pr*Tloc)) 
 
        k26f_inf = (1.81000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC26f = (3.800e-01_pr)* &
          exp(-Tloc/(7.300e+01_pr)) + (6.200e-01_pr)* &
          exp(-Tloc/(1.200e+03_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r26f) =  &
          getlindratecoeff & 
        (Tloc, k26f_0, k26f_inf, FC26f, M(mM12), Ploc ) 

        k(r27) = (8.43000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r28f) = (1.55000000e+08_pr)*Tloc**(0.000_pr)*&
          exp(-(5.640e+04_pr)/(8.314_pr*Tloc)) 

        k(r29) = (3.30000000e+05_pr)*Tloc**(0.000_pr)*&
          exp(-(3.741e+04_pr)/(8.314_pr*Tloc)) 

        k(r30) = (5.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r31) = (1.10000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.164e+05_pr)/(8.314_pr*Tloc)) 

        k(r32) = (1.52800000e+06_pr)*Tloc**(-0.095_pr)*&
          exp(-(2.720e+04_pr)/(8.314_pr*Tloc)) 

        k33f_0 = (2.47000000e+21_pr)*Tloc**(-4.760_pr)*&
          exp(-(1.021e+04_pr)/(8.314_pr*Tloc)) 
 
        k33f_inf = (1.27000000e+10_pr)*Tloc**(-0.630_pr)* &
          exp(-(1.602e+03_pr)/(8.314_pr*Tloc)) 
 
        FC33f = (2.170e-01_pr)* &
          exp(-Tloc/(7.400e+01_pr)) + (7.830e-01_pr)* &
          exp(-Tloc/(2.941e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.964e+03_pr)/Tloc) 
 
        k(r33f) =  &
          getlindratecoeff & 
        (Tloc, k33f_0, k33f_inf, FC33f, M(mM11), Ploc ) 

        k(r34f) = (1.30000000e-02_pr)*Tloc**(3.000_pr)*&
          exp(-(3.363e+04_pr)/(8.314_pr*Tloc)) 

        k(r35f) = (1.60000000e+01_pr)*Tloc**(1.830_pr)*&
          exp(-(1.164e+04_pr)/(8.314_pr*Tloc)) 

        k(r36) = (1.90000000e+03_pr)*Tloc**(1.440_pr)*&
          exp(-(3.630e+04_pr)/(8.314_pr*Tloc)) 

        k37_0 = (1.55000000e+12_pr)*Tloc**(-2.790_pr)*&
          exp(-(1.754e+04_pr)/(8.314_pr*Tloc)) 
 
        k37_inf = (1.80000000e+05_pr)*Tloc**(0.000_pr)* &
          exp(-(9.975e+03_pr)/(8.314_pr*Tloc)) 
 
        FC37 = (0.000e+00_pr)* &
          exp(-Tloc/(1.000e+00_pr)) + (1.000e+00_pr)* &
          exp(-Tloc/(1.000e+07_pr)) + (1.000e+00_pr)* &
          exp(-(1.000e+07_pr)/Tloc) 
 
        k(r37) =  &
          getlindratecoeff & 
        (Tloc, k37_0, k37_inf, FC37, M(mM8), Ploc ) 

        k(r38f) = (4.40000000e+00_pr)*Tloc**(1.500_pr)*&
          exp(-(-3.100e+03_pr)/(8.314_pr*Tloc)) 

        k(r39) = (2.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(9.600e+04_pr)/(8.314_pr*Tloc)) 

        k(r40f) = (1.86000000e+11_pr)*Tloc**(-1.000_pr)*&
          exp(-(7.113e+04_pr)/(8.314_pr*Tloc)) 

        k(r41) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r42) = (7.58000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(1.715e+03_pr)/(8.314_pr*Tloc)) 

        k(r43) = (5.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r44) = (5.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r45) = (3.50000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.470e+04_pr)/(8.314_pr*Tloc)) 

        k(r46) = (3.90000000e+04_pr)*Tloc**(0.890_pr)*&
          exp(-(1.700e+03_pr)/(8.314_pr*Tloc)) 

        k(r47) = (4.11000000e-02_pr)*Tloc**(2.500_pr)*&
          exp(-(4.272e+04_pr)/(8.314_pr*Tloc)) 

        k(r48) = (5.74000000e+01_pr)*Tloc**(1.900_pr)*&
          exp(-(1.150e+04_pr)/(8.314_pr*Tloc)) 

        k49_0 = (1.89500000e+20_pr)*Tloc**(-2.669_pr)*&
          exp(-(3.714e+05_pr)/(8.314_pr*Tloc)) 
 
        k49_inf = (1.53000000e+14_pr)*Tloc**(0.381_pr)* &
          exp(-(3.685e+05_pr)/(8.314_pr*Tloc)) 
 
        FC49 = (2.176e-01_pr)* &
          exp(-Tloc/(2.710e+02_pr)) + (7.824e-01_pr)* &
          exp(-Tloc/(2.755e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.570e+03_pr)/Tloc) 
 
        k(r49) =  &
          getlindratecoeff & 
        (Tloc, k49_0, k49_inf, FC49, M(mM10), Ploc ) 

        k(r50f) = (1.00000000e+08_pr)*Tloc**(0.000_pr)*&
          exp(-(8.000e+04_pr)/(8.314_pr*Tloc)) 

        k(r51) = (7.78000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(5.654e+04_pr)/(8.314_pr*Tloc)) 

        k(r52) = (4.28000000e-19_pr)*Tloc**(7.600_pr)*&
          exp(-(-1.480e+04_pr)/(8.314_pr*Tloc)) 

        k(r53) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r54) = (2.40000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r55) = (5.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.051e+05_pr)/(8.314_pr*Tloc)) 

        k(r56f) = (2.50000000e+11_pr)*Tloc**(-0.930_pr)*&
          exp(-(2.145e+04_pr)/(8.314_pr*Tloc)) 

        k(r57f) = (5.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r58) = (3.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r59) = (4.60000000e+09_pr)*Tloc**(-0.540_pr)*&
          exp(-(1.880e+05_pr)/(8.314_pr*Tloc)) 

        k(r60f) = (1.90000000e+01_pr)*Tloc**(1.700_pr)*&
          exp(-(4.180e+03_pr)/(8.314_pr*Tloc)) 

        k(r61) = (5.19000000e+09_pr)*Tloc**(-1.260_pr)*&
          exp(-(1.386e+04_pr)/(8.314_pr*Tloc)) 

        k(r62) = (7.00000000e+08_pr)*Tloc**(-0.611_pr)*&
          exp(-(2.202e+04_pr)/(8.314_pr*Tloc)) 

        k63f_0 = (4.27000000e+46_pr)*Tloc**(-11.940_pr)*&
          exp(-(4.088e+04_pr)/(8.314_pr*Tloc)) 
 
        k63f_inf = (2.50000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC63f = (8.250e-01_pr)* &
          exp(-Tloc/(1.341e+03_pr)) + (1.750e-01_pr)* &
          exp(-Tloc/(6.000e+04_pr)) + (1.000e+00_pr)* &
          exp(-(1.014e+04_pr)/Tloc) 
 
        k(r63f) =  &
          getlindratecoeff & 
        (Tloc, k63f_0, k63f_inf, FC63f, M(mM30), Ploc ) 

        k(r64) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k65f_0 = (1.51000000e+08_pr)*Tloc**(0.100_pr)*&
          exp(-(1.368e+05_pr)/(8.314_pr*Tloc)) 
 
        k65f_inf = (6.38000000e+09_pr)*Tloc**(1.000_pr)* &
          exp(-(1.574e+05_pr)/(8.314_pr*Tloc)) 
 
        FC65f = (7.000e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (3.000e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r65f) =  &
          getlindratecoeff & 
        (Tloc, k65f_0, k65f_inf, FC65f, M(mM19), Ploc ) 

        k(r66) = (1.70000000e+23_pr)*Tloc**(-5.312_pr)*&
          exp(-(2.721e+04_pr)/(8.314_pr*Tloc)) 

        k(r67) = (4.67800000e+02_pr)*Tloc**(0.136_pr)*&
          exp(-(-6.248e+04_pr)/(8.314_pr*Tloc)) 

        k(r68) = (1.21000000e+00_pr)*Tloc**(2.080_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r69) = (2.23000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(7.192e+04_pr)/(8.314_pr*Tloc)) 

        k(r70) = (2.25000000e+00_pr)*Tloc**(2.080_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r71) = (3.50000000e+10_pr)*Tloc**(0.000_pr)*&
          exp(-(2.993e+05_pr)/(8.314_pr*Tloc)) 

        k(r72f) = (4.49000000e+01_pr)*Tloc**(2.120_pr)*&
          exp(-(5.590e+04_pr)/(8.314_pr*Tloc)) 

        k(r73f) = (5.53000000e-01_pr)*Tloc**(2.310_pr)*&
          exp(-(1.240e+04_pr)/(8.314_pr*Tloc)) 

        k74_0 = (3.94400000e+27_pr)*Tloc**(-6.221_pr)*&
          exp(-(2.685e+04_pr)/(8.314_pr*Tloc)) 
 
        k74_inf = (7.12300000e+11_pr)*Tloc**(-1.021_pr)* &
          exp(-(6.025e+03_pr)/(8.314_pr*Tloc)) 
 
        FC74 = (1.600e-01_pr)* &
          exp(-Tloc/(1.250e+02_pr)) + (8.400e-01_pr)* &
          exp(-Tloc/(2.219e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.882e+03_pr)/Tloc) 
 
        k(r74) =  &
          getlindratecoeff & 
        (Tloc, k74_0, k74_inf, FC74, M(mM15), Ploc ) 

        k(r75) = (7.50000000e+08_pr)*Tloc**(-1.000_pr)*&
          exp(-(2.008e+04_pr)/(8.314_pr*Tloc)) 

        k(r76) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k77f_0 = (3.99000000e+27_pr)*Tloc**(-4.990_pr)*&
          exp(-(1.674e+05_pr)/(8.314_pr*Tloc)) 
 
        k77f_inf = (1.11000000e+10_pr)*Tloc**(1.037_pr)* &
          exp(-(1.538e+05_pr)/(8.314_pr*Tloc)) 
 
        FC77f = (8.320e-01_pr)* &
          exp(-Tloc/(1.200e+03_pr)) + (1.680e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r77f) =  &
          getlindratecoeff & 
        (Tloc, k77f_0, k77f_inf, FC77f, M(mM16), Ploc ) 

        k(r80) = (2.20000000e+01_pr)*Tloc**(1.900_pr)*&
          exp(-(4.700e+03_pr)/(8.314_pr*Tloc)) 

        k(r81f) = (5.40000000e-04_pr)*Tloc**(3.500_pr)*&
          exp(-(2.180e+04_pr)/(8.314_pr*Tloc)) 

        k(r82) = (5.50000000e-07_pr)*Tloc**(4.000_pr)*&
          exp(-(3.470e+04_pr)/(8.314_pr*Tloc)) 

        k(r83) = (1.40000000e-06_pr)*Tloc**(4.300_pr)*&
          exp(-(1.160e+04_pr)/(8.314_pr*Tloc)) 

        k(r84) = (9.00000000e+04_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r85f) = (1.50000000e+03_pr)*Tloc**(1.430_pr)*&
          exp(-(1.125e+04_pr)/(8.314_pr*Tloc)) 

        k(r86) = (1.02000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r87) = (5.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r88) = (3.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r89) = (4.90000000e+08_pr)*Tloc**(-0.500_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r90f) = (1.04700000e+37_pr)*Tloc**(-7.189_pr)*&
          exp(-(1.855e+05_pr)/(8.314_pr*Tloc)) 

        k(r91) = (3.00000000e+04_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r92) = (7.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r94) = (4.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(7.116e+04_pr)/(8.314_pr*Tloc)) 

        k95f_0 = (6.26000000e+26_pr)*Tloc**(-6.660_pr)*&
          exp(-(2.929e+04_pr)/(8.314_pr*Tloc)) 
 
        k95f_inf = (1.33000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(1.364e+04_pr)/(8.314_pr*Tloc)) 
 
        FC95f = (0.000e+00_pr)* &
          exp(-Tloc/(1.000e+03_pr)) + (1.000e+00_pr)* &
          exp(-Tloc/(1.310e+03_pr)) + (1.000e+00_pr)* &
          exp(-(4.810e+04_pr)/Tloc) 
 
        k(r95f) =  &
          getlindratecoeff & 
        (Tloc, k95f_0, k95f_inf, FC95f, M(mM32), Ploc ) 

        k(r96) = (1.20000000e+02_pr)*Tloc**(1.650_pr)*&
          exp(-(1.370e+03_pr)/(8.314_pr*Tloc)) 

        k97f_0 = (8.70000000e+30_pr)*Tloc**(-7.500_pr)*&
          exp(-(1.980e+04_pr)/(8.314_pr*Tloc)) 
 
        k97f_inf = (1.33000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(6.530e+03_pr)/(8.314_pr*Tloc)) 
 
        FC97f = (0.000e+00_pr)* &
          exp(-Tloc/(1.000e+03_pr)) + (1.000e+00_pr)* &
          exp(-Tloc/(6.454e+02_pr)) + (1.000e+00_pr)* &
          exp(-(6.844e+03_pr)/Tloc) 
 
        k(r97f) =  &
          getlindratecoeff & 
        (Tloc, k97f_0, k97f_inf, FC97f, M(mM31), Ploc ) 

        k(r98) = (1.60000000e+16_pr)*Tloc**(-2.390_pr)*&
          exp(-(4.680e+04_pr)/(8.314_pr*Tloc)) 

        k(r99) = (3.50000000e+01_pr)*Tloc**(1.650_pr)*&
          exp(-(-4.070e+03_pr)/(8.314_pr*Tloc)) 

        k(r100f) = (1.30000000e+05_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k(r101) = (3.50000000e+10_pr)*Tloc**(-1.600_pr)*&
          exp(-(1.464e+04_pr)/(8.314_pr*Tloc)) 

        k(r102f) = (2.00000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

        k103f_0 = (5.49000000e+43_pr)*Tloc**(-10.000_pr)*&
          exp(-(1.497e+05_pr)/(8.314_pr*Tloc)) 
 
        k103f_inf = (1.23000000e+13_pr)*Tloc**(-0.100_pr)* &
          exp(-(1.264e+05_pr)/(8.314_pr*Tloc)) 
 
        FC103f = (2.170e+00_pr)* &
          exp(-Tloc/(2.510e+02_pr)) + (-1.170e+00_pr)* &
          exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)* &
          exp(-(1.185e+03_pr)/Tloc) 
 
        k(r103f) =  &
          getlindratecoeff & 
        (Tloc, k103f_0, k103f_inf, FC103f, M(mM0), Ploc ) 

        k(r104) = (4.76000000e-02_pr)*Tloc**(2.550_pr)*&
          exp(-(6.900e+04_pr)/(8.314_pr*Tloc)) 

        k(r105) = (9.64000000e-03_pr)*Tloc**(2.600_pr)*&
          exp(-(5.823e+04_pr)/(8.314_pr*Tloc)) 

        k(r106f) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(1.987e+05_pr)/(8.314_pr*Tloc)) 

        k(r107f) = (1.33000000e+00_pr)*Tloc**(2.540_pr)*&
          exp(-(2.829e+04_pr)/(8.314_pr*Tloc)) 

        k(r108) = (4.00000000e+07_pr)*Tloc**(0.000_pr)*&
          exp(-(2.131e+05_pr)/(8.314_pr*Tloc)) 

        k(r109) = (1.90000000e-01_pr)*Tloc**(2.680_pr)*&
          exp(-(1.556e+04_pr)/(8.314_pr*Tloc)) 

        k110f_0 = (7.83000000e+12_pr)*Tloc**(0.000_pr)*&
          exp(-(2.719e+05_pr)/(8.314_pr*Tloc)) 
 
        k110f_inf = (1.10000000e+17_pr)*Tloc**(0.000_pr)* &
          exp(-(3.531e+05_pr)/(8.314_pr*Tloc)) 
 
        FC110f = (2.400e-01_pr)* &
          exp(-Tloc/(1.900e+03_pr)) + (7.600e-01_pr)* &
          exp(-Tloc/(3.800e+01_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r110f) =  &
          getlindratecoeff & 
        (Tloc, k110f_0, k110f_inf, FC110f, M(mM0), Ploc ) 

        k(r111) = (1.00000000e+04_pr)*Tloc**(1.000_pr)*&
          exp(-(6.694e+03_pr)/(8.314_pr*Tloc)) 

        k(r112) = (4.76000000e-02_pr)*Tloc**(2.710_pr)*&
          exp(-(8.817e+03_pr)/(8.314_pr*Tloc)) 

        k(r113f) = (1.30000000e+00_pr)*Tloc**(2.400_pr)*&
          exp(-(1.871e+04_pr)/(8.314_pr*Tloc)) 

        k(r114) = (1.50000000e+02_pr)*Tloc**(0.000_pr)*&
          exp(-(-2.929e+04_pr)/(8.314_pr*Tloc)) 

        k(r115) = (2.50000000e+35_pr)*Tloc**(-8.300_pr)*&
          exp(-(9.205e+04_pr)/(8.314_pr*Tloc)) 

        k(r116) = (1.00000000e+15_pr)*Tloc**(0.000_pr)*&
          exp(-(1.799e+05_pr)/(8.314_pr*Tloc)) 

        k(r2b) = (3.08572107e-02_pr)*Tloc**(2.63_pr)*&
          exp(-(2.0241e+04_pr)/(8.314_pr*Tloc)) 

        k(r4b) = (1.34762266e+17_pr)*Tloc**(-1.79_pr)*&
          exp(-(4.9633e+05_pr)/(8.314_pr*Tloc)) 

        k(r5b) = (6.98293820e+12_pr)*Tloc**(-0.47_pr)*&
          exp(-(2.7562e+05_pr)/(8.314_pr*Tloc)) 

        k6b_0 = (1.19344641e+25_pr)*Tloc**(-4.25_pr)*&
          exp(-(2.1486e+05_pr)/(8.314_pr*Tloc)) 
 
        k6b_inf = (4.12949754e+19_pr)*Tloc**(-1.32_pr)* &
          exp(-(2.1486e+05_pr)/(8.314_pr*Tloc)) 
 
        FC6b = (4.300e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (5.700e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r6b) =  &
          getlindratecoeff & 
        (Tloc, k6b_0, k6b_inf, FC6b, M(mM7), Ploc ) 

        k(r7b) = (1.45933806e+04_pr)*Tloc**(1.18_pr)*&
          exp(-(7.8338e+04_pr)/(8.314_pr*Tloc)) 

        k(r8b) = (3.42242660e-02_pr)*Tloc**(2.41_pr)*&
          exp(-(-8.3362e+03_pr)/(8.314_pr*Tloc)) 

        k(r9b) = (7.19106655e+07_pr)*Tloc**(-0.27_pr)*&
          exp(-(6.2578e+02_pr)/(8.314_pr*Tloc)) 

        k10b_0 = (1.02533717e+14_pr)*Tloc**(-1.44_pr)*&
          exp(-(2.0482e+05_pr)/(8.314_pr*Tloc)) 
 
        k10b_inf = (8.29185711e+12_pr)*Tloc**(0.40_pr)* &
          exp(-(2.0482e+05_pr)/(8.314_pr*Tloc)) 
 
        FC10b = (5.000e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (5.000e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r10b) =  &
          getlindratecoeff & 
        (Tloc, k10b_0, k10b_inf, FC10b, M(mM6), Ploc ) 

        k(r13b) = (8.50201705e+08_pr)*Tloc**(0.25_pr)*&
          exp(-(3.3724e+05_pr)/(8.314_pr*Tloc)) 

        k(r14b) = (1.32253599e+07_pr)*Tloc**(0.25_pr)*&
          exp(-(2.8693e+05_pr)/(8.314_pr*Tloc)) 

        k(r23b) = (1.37132417e+06_pr)*Tloc**(0.48_pr)*&
          exp(-(5.2179e+03_pr)/(8.314_pr*Tloc)) 

        k(r25b) = (5.11170724e+12_pr)*Tloc**(-1.18_pr)*&
          exp(-(2.4677e+04_pr)/(8.314_pr*Tloc)) 

        k26b_0 = (3.51669517e+43_pr)*Tloc**(-8.44_pr)*&
          exp(-(3.9663e+05_pr)/(8.314_pr*Tloc)) 
 
        k26b_inf = (5.01198288e+21_pr)*Tloc**(-1.44_pr)* &
          exp(-(3.8507e+05_pr)/(8.314_pr*Tloc)) 
 
        FC26b = (3.800e-01_pr)* &
          exp(-Tloc/(7.300e+01_pr)) + (6.200e-01_pr)* &
          exp(-Tloc/(1.200e+03_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r26b) =  &
          getlindratecoeff & 
        (Tloc, k26b_0, k26b_inf, FC26b, M(mM12), Ploc ) 

        k(r28b) = (4.26031577e+05_pr)*Tloc**(0.60_pr)*&
          exp(-(-1.1980e+04_pr)/(8.314_pr*Tloc)) 

        k33b_0 = (1.75987289e+30_pr)*Tloc**(-5.00_pr)*&
          exp(-(4.5275e+05_pr)/(8.314_pr*Tloc)) 
 
        k33b_inf = (9.04873914e+18_pr)*Tloc**(-0.87_pr)* &
          exp(-(4.4415e+05_pr)/(8.314_pr*Tloc)) 
 
        FC33b = (2.170e-01_pr)* &
          exp(-Tloc/(7.400e+01_pr)) + (7.830e-01_pr)* &
          exp(-Tloc/(2.941e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.964e+03_pr)/Tloc) 
 
        k(r33b) =  &
          getlindratecoeff & 
        (Tloc, k33b_0, k33b_inf, FC33b, M(mM11), Ploc ) 

        k(r34b) = (4.92830154e-06_pr)*Tloc**(3.57_pr)*&
          exp(-(2.4293e+04_pr)/(8.314_pr*Tloc)) 

        k(r35b) = (7.56560999e-02_pr)*Tloc**(2.28_pr)*&
          exp(-(6.5430e+04_pr)/(8.314_pr*Tloc)) 

        k(r38b) = (2.40355366e+07_pr)*Tloc**(0.22_pr)*&
          exp(-(1.0458e+05_pr)/(8.314_pr*Tloc)) 

        k(r40b) = (2.89407746e+04_pr)*Tloc**(-0.74_pr)*&
          exp(-(5.2287e+03_pr)/(8.314_pr*Tloc)) 

        k(r50b) = (3.78931373e+05_pr)*Tloc**(0.56_pr)*&
          exp(-(1.1461e+05_pr)/(8.314_pr*Tloc)) 

        k(r56b) = (4.35956032e+07_pr)*Tloc**(-0.02_pr)*&
          exp(-(3.4605e+04_pr)/(8.314_pr*Tloc)) 

        k(r57b) = (1.27426086e+06_pr)*Tloc**(0.30_pr)*&
          exp(-(7.9735e+04_pr)/(8.314_pr*Tloc)) 

        k(r60b) = (1.22889550e+05_pr)*Tloc**(0.96_pr)*&
          exp(-(1.0301e+05_pr)/(8.314_pr*Tloc)) 

        k63b_0 = (6.15094823e+61_pr)*Tloc**(-13.48_pr)*&
          exp(-(4.7310e+05_pr)/(8.314_pr*Tloc)) 
 
        k63b_inf = (3.60125775e+22_pr)*Tloc**(-1.54_pr)* &
          exp(-(4.3222e+05_pr)/(8.314_pr*Tloc)) 
 
        FC63b = (8.250e-01_pr)* &
          exp(-Tloc/(1.341e+03_pr)) + (1.750e-01_pr)* &
          exp(-Tloc/(6.000e+04_pr)) + (1.000e+00_pr)* &
          exp(-(1.014e+04_pr)/Tloc) 
 
        k(r63b) =  &
          getlindratecoeff & 
        (Tloc, k63b_0, k63b_inf, FC63b, M(mM30), Ploc ) 

        k65b_0 = (5.26609134e+02_pr)*Tloc**(-0.04_pr)*&
          exp(-(-1.0169e+04_pr)/(8.314_pr*Tloc)) 
 
        k65b_inf = (2.22501078e+04_pr)*Tloc**(0.86_pr)* &
          exp(-(1.0431e+04_pr)/(8.314_pr*Tloc)) 
 
        FC65b = (7.000e-01_pr)* &
          exp(-Tloc/(1.000e+30_pr)) + (3.000e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r65b) =  &
          getlindratecoeff & 
        (Tloc, k65b_0, k65b_inf, FC65b, M(mM19), Ploc ) 

        k(r72b) = (1.56275683e-02_pr)*Tloc**(2.63_pr)*&
          exp(-(2.2357e+04_pr)/(8.314_pr*Tloc)) 

        k(r73b) = (2.40071300e-03_pr)*Tloc**(2.70_pr)*&
          exp(-(4.1984e+04_pr)/(8.314_pr*Tloc)) 

        k77b_0 = (1.72614099e+21_pr)*Tloc**(-4.84_pr)*&
          exp(-(1.5915e+04_pr)/(8.314_pr*Tloc)) 
 
        k77b_inf = (4.80204635e+03_pr)*Tloc**(1.19_pr)* &
          exp(-(2.3148e+03_pr)/(8.314_pr*Tloc)) 
 
        FC77b = (8.320e-01_pr)* &
          exp(-Tloc/(1.200e+03_pr)) + (1.680e-01_pr)* &
          exp(-Tloc/(1.000e-30_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r77b) =  &
          getlindratecoeff & 
        (Tloc, k77b_0, k77b_inf, FC77b, M(mM16), Ploc ) 

        k(r81b) = (8.52079101e-08_pr)*Tloc**(4.08_pr)*&
          exp(-(3.3110e+04_pr)/(8.314_pr*Tloc)) 

        k(r85b) = (4.79124352e-04_pr)*Tloc**(2.99_pr)*&
          exp(-(1.4047e+05_pr)/(8.314_pr*Tloc)) 

        k(r90b) = (5.16504907e+28_pr)*Tloc**(-6.62_pr)*&
          exp(-(2.5864e+04_pr)/(8.314_pr*Tloc)) 

        k95b_0 = (2.50322114e+33_pr)*Tloc**(-6.97_pr)*&
          exp(-(1.6986e+05_pr)/(8.314_pr*Tloc)) 
 
        k95b_inf = (5.31834523e+13_pr)*Tloc**(-0.31_pr)* &
          exp(-(1.5421e+05_pr)/(8.314_pr*Tloc)) 
 
        FC95b = (0.000e+00_pr)* &
          exp(-Tloc/(1.000e+03_pr)) + (1.000e+00_pr)* &
          exp(-Tloc/(1.310e+03_pr)) + (1.000e+00_pr)* &
          exp(-(4.810e+04_pr)/Tloc) 
 
        k(r95b) =  &
          getlindratecoeff & 
        (Tloc, k95b_0, k95b_inf, FC95b, M(mM32), Ploc ) 

        k97b_0 = (2.85140703e+36_pr)*Tloc**(-7.39_pr)*&
          exp(-(1.6819e+05_pr)/(8.314_pr*Tloc)) 
 
        k97b_inf = (4.35904753e+12_pr)*Tloc**(0.11_pr)* &
          exp(-(1.5492e+05_pr)/(8.314_pr*Tloc)) 
 
        FC97b = (0.000e+00_pr)* &
          exp(-Tloc/(1.000e+03_pr)) + (1.000e+00_pr)* &
          exp(-Tloc/(6.454e+02_pr)) + (1.000e+00_pr)* &
          exp(-(6.844e+03_pr)/Tloc) 
 
        k(r97b) =  &
          getlindratecoeff & 
        (Tloc, k97b_0, k97b_inf, FC97b, M(mM31), Ploc ) 

        k(r100b) = (7.07297699e+05_pr)*Tloc**(-0.15_pr)*&
          exp(-(5.6430e+04_pr)/(8.314_pr*Tloc)) 

        k(r102b) = (1.56621153e+19_pr)*Tloc**(-1.89_pr)*&
          exp(-(9.7992e+04_pr)/(8.314_pr*Tloc)) 

        k103b_0 = (7.39650348e+30_pr)*Tloc**(-8.33_pr)*&
          exp(-(4.3656e+04_pr)/(8.314_pr*Tloc)) 
 
        k103b_inf = (1.65714012e+00_pr)*Tloc**(1.57_pr)* &
          exp(-(2.0356e+04_pr)/(8.314_pr*Tloc)) 
 
        FC103b = (2.170e+00_pr)* &
          exp(-Tloc/(2.510e+02_pr)) + (-1.170e+00_pr)* &
          exp(-Tloc/(1.000e-15_pr)) + (1.000e+00_pr)* &
          exp(-(1.185e+03_pr)/Tloc) 
 
        k(r103b) =  &
          getlindratecoeff & 
        (Tloc, k103b_0, k103b_inf, FC103b, M(mM0), Ploc ) 

        k(r106b) = (2.60600114e+03_pr)*Tloc**(0.64_pr)*&
          exp(-(-1.0652e+04_pr)/(8.314_pr*Tloc)) 

        k(r107b) = (1.60136592e-04_pr)*Tloc**(3.14_pr)*&
          exp(-(3.9503e+04_pr)/(8.314_pr*Tloc)) 

        k110b_0 = (1.08695906e-03_pr)*Tloc**(1.79_pr)*&
          exp(-(-1.0465e+05_pr)/(8.314_pr*Tloc)) 
 
        k110b_inf = (1.52701784e+01_pr)*Tloc**(1.79_pr)* &
          exp(-(-2.3452e+04_pr)/(8.314_pr*Tloc)) 
 
        FC110b = (2.400e-01_pr)* &
          exp(-Tloc/(1.900e+03_pr)) + (7.600e-01_pr)* &
          exp(-Tloc/(3.800e+01_pr)) + (0.000e+00_pr)* &
          exp(-(0.000e+00_pr)/Tloc) 
 
        k(r110b) =  &
          getlindratecoeff & 
        (Tloc, k110b_0, k110b_inf, FC110b, M(mM0), Ploc ) 

        k(r113b) = (1.28291349e-05_pr)*Tloc**(3.41_pr)*&
          exp(-(3.7742e+04_pr)/(8.314_pr*Tloc)) 

      end subroutine get_rate_coefficients

! --- 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) :: w,k
      real(pr), dimension(24) :: m

      w(r1) = k(r1)  * c(sH2) * m(mM1)
      w(r2f) = k(r2f)  * c(sH2) * cqss(sO-nspec)
      w(r3) = k(r3)  * c(sH) * cqss(sO-nspec) * m(mM4)
      w(r4f) = k(r4f)  * c(sH) * c(sOH) * m(mM2)
      w(r5f) = k(r5f)  * cqss(sO-nspec) * c(sOH) * m(mM5)
      w(r6f) = k(r6f)  * c(sOH)**2_pr 
      w(r7f) = k(r7f)  * c(sH2) * c(sOH)
      w(r8f) = k(r8f)  * c(sH2O) * cqss(sO-nspec)
      w(r9f) = k(r9f)  * c(sH) * c(sO2)
      w(r10f) = k(r10f)  * c(sH) * c(sO2)
      w(r11) = k(r11)  * c(sHO2) * cqss(sO-nspec)
      w(r12) = k(r12)  * c(sHO2) * c(sH)
      w(r13f) = k(r13f)  * c(sHO2) * c(sOH)
      w(r14f) = k(r14f)  * c(sHO2) * c(sOH)
      w(r15) = k(r15)  * c(sHO2) * c(sH)
      w(r16) = k(r16)  * c(sHO2)**2_pr 
      w(r17) = k(r17)  * c(sHO2)**2_pr 
      w(r18) = k(r18)  * c(sHO2) * c(sH)
      w(r19) = k(r19)  * c(sH2O2) * c(sOH)
      w(r20) = k(r20)  * c(sH2O2) * c(sOH)
      w(r21) = k(r21)  * cqss(sSXCH2-nspec) * c(sOH)
      w(r22) = k(r22)  * cqss(sSXCH2-nspec) * c(sO2)
      w(r23f) = k(r23f)  * c(sCH3) * c(sOH)
      w(r24) = k(r24)  * c(sHO2) * c(sCH3)
      w(r25f) = k(r25f)  * c(sCH3)**2_pr 
      w(r26f) = k(r26f)  * c(sCH3)**2_pr 
      w(r27) = k(r27)  * c(sCH3) * cqss(sO-nspec)
      w(r28f) = k(r28f)  * c(sCH3) * c(sH)
      w(r29) = k(r29)  * c(sCH3) * c(sO2)
      w(r30) = k(r30)  * c(sCH3) * c(sHO2)
      w(r31) = k(r31)  * c(sCH3) * c(sO2)
      w(r32) = k(r32)  * c(sH2O2) * c(sCH3)
      w(r33f) = k(r33f)  * c(sH) * c(sCH3)
      w(r34f) = k(r34f)  * c(sCH4) * c(sH)
      w(r35f) = k(r35f)  * c(sCH4) * c(sOH)
      w(r36) = k(r36)  * c(sCH4) * cqss(sO-nspec)
      w(r37) = k(r37)  * c(sCO) * cqss(sO-nspec)
      w(r38f) = k(r38f)  * c(sCO) * c(sOH)
      w(r39) = k(r39)  * c(sCO) * c(sHO2)
      w(r40f) = k(r40f)  * cqss(sHCO-nspec) * m(mM9)
      w(r41) = k(r41)  * cqss(sHCO-nspec) * c(sOH)
      w(r42) = k(r42)  * cqss(sHCO-nspec) * c(sO2)
      w(r43) = k(r43)  * cqss(sHCO-nspec) * c(sH)
      w(r44) = k(r44)  * cqss(sHCO-nspec) * c(sCH3)
      w(r45) = k(r45)  * c(sCH2O) * cqss(sO-nspec)
      w(r46) = k(r46)  * c(sCH2O) * c(sOH)
      w(r47) = k(r47)  * c(sCH2O) * c(sHO2)
      w(r48) = k(r48)  * c(sCH2O) * c(sH)
      w(r49) = k(r49)  * c(sCH2O)
      w(r50f) = k(r50f)  * cqss(sCH3O-nspec) * m(mM21)
      w(r51) = k(r51)  * cqss(sCH3O-nspec) * m(mM14)
      w(r52) = k(r52)  * cqss(sCH3O-nspec) * c(sO2)
      w(r53) = k(r53)  * cqss(sCH2OH-nspec) * c(sH)
      w(r54) = k(r54)  * cqss(sCH2OH-nspec) * c(sOH)
      w(r55) = k(r55)  * cqss(sCH2OH-nspec) * m(mM20)
      w(r56f) = k(r56f)  * cqss(sCH2OH-nspec) * c(sH)
      w(r57f) = k(r57f)  * cqss(sCH2OH-nspec) * c(sO2)
      w(r58) = k(r58)  * cqss(sSXCH2-nspec) * c(sCO2)
      w(r59) = k(r59)  * c(sC2H2) * c(sO2)
      w(r60f) = k(r60f)  * c(sC2H2) * c(sOH)
      w(r61) = k(r61)  * cqss(sC2H3-nspec) * c(sO2)
      w(r62) = k(r62)  * cqss(sC2H3-nspec) * c(sO2)
      w(r63f) = k(r63f)  * cqss(sC2H3-nspec) * c(sCH3)
      w(r64) = k(r64)  * cqss(sC2H3-nspec) * c(sH)
      w(r65f) = k(r65f)  * cqss(sC2H3-nspec)
      w(r66) = k(r66)  * cqss(sC2H3-nspec) * c(sO2)
      w(r67) = k(r67)  * c(sH) * cqss(sC2H3-nspec) * m(mM17)
      w(r68) = k(r68)  * c(sC2H4) * cqss(sO-nspec)
      w(r69) = k(r69)  * c(sC2H4) * c(sHO2)
      w(r70) = k(r70)  * c(sC2H4) * cqss(sO-nspec)
      w(r71) = k(r71)  * c(sC2H4) * m(mM18)
      w(r72f) = k(r72f)  * c(sC2H4) * c(sH)
      w(r73f) = k(r73f)  * c(sC2H4) * c(sOH)
      w(r74) = k(r74)  * c(sH) * cqss(sC2H5-nspec)
      w(r75) = k(r75)  * cqss(sC2H5-nspec) * c(sO2)
      w(r76) = k(r76)  * cqss(sC2H5-nspec) * c(sH)
      w(r77f) = k(r77f)  * cqss(sC2H5-nspec)
      w(r80) = k(r80)  * c(sC2H6) * c(sOH)
      w(r81f) = k(r81f)  * c(sC2H6) * c(sH)
      w(r82) = k(r82)  * c(sC2H6) * c(sCH3)
      w(r83) = k(r83)  * c(sC2H6) * cqss(sO-nspec)
      w(r84) = k(r84)  * c(sCH2CO) * c(sCH3)
      w(r85f) = k(r85f)  * c(sCH2CO) * c(sH)
      w(r86) = k(r86)  * c(sCH2CO) * c(sOH)
      w(r87) = k(r87)  * cqss(sCH2CHO-nspec) * c(sH)
      w(r88) = k(r88)  * cqss(sCH2CHO-nspec) * c(sOH)
      w(r89) = k(r89)  * cqss(sCH2CHO-nspec) * c(sCH3)
      w(r90f) = k(r90f)  * cqss(sCH2CHO-nspec)
      w(r91) = k(r91)  * cqss(sCH2CHO-nspec) * c(sO2)
      w(r92) = k(r92)  * cqss(sCH2CHO-nspec) * c(sHO2)
      w(r94) = k(r94)  * c(sC2H4O) * c(sHO2)
      w(r95f) = k(r95f)  * c(sH) * c(sC3H6)
      w(r96) = k(r96)  * c(sC3H6) * cqss(sO-nspec)
      w(r97f) = k(r97f)  * c(sC3H6) * c(sH)
      w(r98) = k(r98)  * c(sC3H6) * c(sH)
      w(r99) = k(r99)  * c(sC3H6) * cqss(sO-nspec)
      w(r100f) = k(r100f)  * cqss(sIXC3H7-nspec) * c(sO2)
      w(r101) = k(r101)  * cqss(sNXC3H7-nspec) * c(sO2)
      w(r102f) = k(r102f)  * cqss(sNXC3H7-nspec) * c(sO2)
      w(r103f) = k(r103f)  * cqss(sNXC3H7-nspec)
      w(r104) = k(r104)  * c(sC3H8) * c(sHO2)
      w(r105) = k(r105)  * c(sC3H8) * c(sHO2)
      w(r106f) = k(r106f)  * c(sC3H8) * c(sO2)
      w(r107f) = k(r107f)  * c(sC3H8) * c(sH)
      w(r108) = k(r108)  * c(sC3H8) * c(sO2)
      w(r109) = k(r109)  * c(sC3H8) * cqss(sO-nspec)
      w(r110f) = k(r110f)  * c(sC3H8)
      w(r111) = k(r111)  * c(sC3H8) * c(sOH)
      w(r112) = k(r112)  * c(sC3H8) * cqss(sO-nspec)
      w(r113f) = k(r113f)  * c(sC3H8) * c(sH)
      w(r114) = k(r114)  * cqss(sC3H6OOH-nspec) * c(sO2)
      w(r115) = k(r115)  * cqss(sC3H6OOH-nspec)
      w(r116) = k(r116)  * c(sOC3H5OOH)
      w(r2b) = k(r2b)  * c(sOH) * c(sH)
      w(r4b) = k(r4b)  * c(sH2O) * m(mM2)
      w(r5b) = k(r5b)  * c(sHO2) * m(mM5)
      w(r6b) = k(r6b)  * c(sH2O2)
      w(r7b) = k(r7b)  * c(sH2O) * c(sH)
      w(r8b) = k(r8b)  * c(sOH)**2_pr 
      w(r9b) = k(r9b)  * c(sOH) * cqss(sO-nspec)
      w(r10b) = k(r10b)  * c(sHO2)
      w(r13b) = k(r13b)  * c(sH2O) * c(sO2)
      w(r14b) = k(r14b)  * c(sH2O) * c(sO2)
      w(r23b) = k(r23b)  * cqss(sSXCH2-nspec) * c(sH2O)
      w(r25b) = k(r25b)  * cqss(sC2H5-nspec) * c(sH)
      w(r26b) = k(r26b)  * c(sC2H6)
      w(r28b) = k(r28b)  * cqss(sSXCH2-nspec) * c(sH2)
      w(r33b) = k(r33b)  * c(sCH4)
      w(r34b) = k(r34b)  * c(sH2) * c(sCH3)
      w(r35b) = k(r35b)  * c(sH2O) * c(sCH3)
      w(r38b) = k(r38b)  * c(sCO2) * c(sH)
      w(r40b) = k(r40b)  * c(sCO) * c(sH) * m(mM9)
      w(r50b) = k(r50b)  * cqss(sCH2OH-nspec) * m(mM21)
      w(r56b) = k(r56b)  * c(sCH3) * c(sOH)
      w(r57b) = k(r57b)  * c(sCH2O) * c(sHO2)
      w(r60b) = k(r60b)  * c(sCH2CO) * c(sH)
      w(r63b) = k(r63b)  * c(sC3H6)
      w(r65b) = k(r65b)  * c(sC2H2) * c(sH)
      w(r72b) = k(r72b)  * cqss(sC2H3-nspec) * c(sH2)
      w(r73b) = k(r73b)  * cqss(sC2H3-nspec) * c(sH2O)
      w(r77b) = k(r77b)  * c(sC2H4) * c(sH)
      w(r81b) = k(r81b)  * cqss(sC2H5-nspec) * c(sH2)
      w(r85b) = k(r85b)  * c(sCH3) * c(sCO)
      w(r90b) = k(r90b)  * c(sCH2CO) * c(sH)
      w(r95b) = k(r95b)  * cqss(sNXC3H7-nspec)
      w(r97b) = k(r97b)  * cqss(sIXC3H7-nspec)
      w(r100b) = k(r100b)  * c(sC3H6) * c(sHO2)
      w(r102b) = k(r102b)  * cqss(sC3H6OOH-nspec)
      w(r103b) = k(r103b)  * c(sCH3) * c(sC2H4)
      w(r106b) = k(r106b)  * cqss(sIXC3H7-nspec) * c(sHO2)
      w(r107b) = k(r107b)  * cqss(sNXC3H7-nspec) * c(sH2)
      w(r110b) = k(r110b)  * c(sCH3) * cqss(sC2H5-nspec)
      w(r113b) = k(r113b)  * cqss(sIXC3H7-nspec) * c(sH2)

      end subroutine get_reaction_rates

! --- Production rates --- !
      subroutine get_production_rates(cdot,w)
      implicit none

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

        cdot(sN2) = 0.0_pr 


        cdot(sH2) = 0.0_pr - &
         w(r1) - &
         w(r2f) - &
         w(r7f) + &
         w(r18) + &
         w(r28f) + &
         w(r34f) + &
         w(r43) + &
         w(r48) + &
         w(r53) + &
         w(r64) + &
         w(r71) + &
         w(r72f) + &
         w(r76) + &
         w(r81f) + &
         w(r107f) + &
         w(r113f) + &
         w(r2b) + &
         w(r7b) - &
         w(r28b) - &
         w(r34b) - &
         w(r72b) - &
         w(r81b) - &
         w(r107b) - &
         w(r113b) 


        cdot(sH) = 0.0_pr + 2_pr * &
         w(r1) + &
         w(r2f) - &
         w(r3) - &
         w(r4f) + &
         w(r7f) - &
         w(r9f) - &
         w(r10f) - &
         w(r12) - &
         w(r15) - &
         w(r18) + &
         w(r21) + &
         w(r22) + &
         w(r25f) + &
         w(r27) - &
         w(r28f) - &
         w(r33f) - &
         w(r34f) + &
         w(r38f) + &
         w(r40f) - &
         w(r43) - &
         w(r48) + &
         w(r49) + &
         w(r51) - &
         w(r53) + &
         w(r55) - &
         w(r56f) + &
         w(r60f) - &
         w(r64) + &
         w(r65f) - &
         w(r67) + &
         w(r68) - &
         w(r72f) - &
         w(r74) - &
         w(r76) + &
         w(r77f) - &
         w(r81f) - &
         w(r85f) - &
         w(r87) + &
         w(r89) + &
         w(r90f) - &
         w(r95f) + &
         w(r96) - &
         w(r97f) - &
         w(r98) - &
         w(r107f) - &
         w(r113f) - &
         w(r2b) + &
         w(r4b) - &
         w(r7b) + &
         w(r9b) 
    cdot(sH) = cdot(sH) + &
         w(r10b) - &
         w(r25b) + &
         w(r28b) + &
         w(r33b) + &
         w(r34b) - &
         w(r38b) - &
         w(r40b) + &
         w(r56b) - &
         w(r60b) - &
         w(r65b) + &
         w(r72b) - &
         w(r77b) + &
         w(r81b) + &
         w(r85b) - &
         w(r90b) + &
         w(r95b) + &
         w(r97b) + &
         w(r107b) + &
         w(r113b) 

      !cdot(sO) = 0.0_pr


        cdot(sOH) = 0.0_pr + &
         w(r2f) + &
         w(r3) - &
         w(r4f) - &
         w(r5f) - 2_pr * &
         w(r6f) - &
         w(r7f) + 2_pr * &
         w(r8f) + &
         w(r9f) + &
         w(r11) + 2_pr * &
         w(r12) - &
         w(r13f) - &
         w(r14f) - &
         w(r19) - &
         w(r20) - &
         w(r21) + &
         w(r22) - &
         w(r23f) + &
         w(r29) + &
         w(r30) - &
         w(r35f) + &
         w(r36) - &
         w(r38f) + &
         w(r39) - &
         w(r41) + &
         w(r45) - &
         w(r46) - &
         w(r54) + &
         w(r56f) - &
         w(r60f) + &
         w(r69) - &
         w(r73f) - &
         w(r80) + &
         w(r83) - &
         w(r86) - &
         w(r88) + &
         w(r91) + &
         w(r92) + &
         w(r109) - &
         w(r111) + &
         w(r112) + &
         w(r114) + &
         w(r116) - &
         w(r2b) + &
         w(r4b) + &
         w(r5b) + 2_pr * &
         w(r6b) + &
         w(r7b) - 2_pr * &
         w(r8b) - &
         w(r9b) + &
         w(r13b) 
    cdot(sOH) = cdot(sOH) + &
         w(r14b) + &
         w(r23b) + &
         w(r35b) + &
         w(r38b) - &
         w(r56b) + &
         w(r60b) + &
         w(r73b) 


        cdot(sH2O) = 0.0_pr + &
         w(r4f) + &
         w(r7f) - &
         w(r8f) + &
         w(r13f) + &
         w(r14f) + &
         w(r15) + &
         w(r19) + &
         w(r20) + &
         w(r23f) + &
         w(r35f) + &
         w(r41) + &
         w(r46) + &
         w(r54) + &
         w(r73f) + &
         w(r80) + &
         w(r88) + &
         w(r111) - &
         w(r4b) - &
         w(r7b) + &
         w(r8b) - &
         w(r13b) - &
         w(r14b) - &
         w(r23b) - &
         w(r35b) - &
         w(r73b) 


        cdot(sHO2) = 0.0_pr + &
         w(r5f) + &
         w(r10f) - &
         w(r11) - &
         w(r12) - &
         w(r13f) - &
         w(r14f) - &
         w(r15) - 2_pr * &
         w(r16) - 2_pr * &
         w(r17) - &
         w(r18) + &
         w(r19) + &
         w(r20) - &
         w(r24) - &
         w(r30) + &
         w(r32) - &
         w(r39) + &
         w(r42) - &
         w(r47) + &
         w(r52) + &
         w(r57f) + &
         w(r61) - &
         w(r69) + &
         w(r75) - &
         w(r92) - &
         w(r94) + &
         w(r100f) + &
         w(r101) - &
         w(r104) - &
         w(r105) + &
         w(r106f) + &
         w(r108) + &
         w(r115) - &
         w(r5b) - &
         w(r10b) + &
         w(r13b) + &
         w(r14b) - &
         w(r57b) - &
         w(r100b) - &
         w(r106b) 


        cdot(sH2O2) = 0.0_pr + &
         w(r6f) + &
         w(r16) + &
         w(r17) - &
         w(r19) - &
         w(r20) - &
         w(r32) + &
         w(r47) + &
         w(r94) + &
         w(r104) + &
         w(r105) - &
         w(r6b) 


        cdot(sO2) = 0.0_pr - &
         w(r9f) - &
         w(r10f) + &
         w(r11) + &
         w(r13f) + &
         w(r14f) + &
         w(r16) + &
         w(r17) + &
         w(r18) - &
         w(r22) + &
         w(r24) - &
         w(r29) - &
         w(r31) - &
         w(r42) - &
         w(r52) - &
         w(r57f) - &
         w(r59) - &
         w(r61) - &
         w(r62) - &
         w(r66) - &
         w(r75) - &
         w(r91) - &
         w(r100f) - &
         w(r101) - &
         w(r102f) - &
         w(r106f) - &
         w(r108) - &
         w(r114) + &
         w(r9b) + &
         w(r10b) - &
         w(r13b) - &
         w(r14b) + &
         w(r57b) + &
         w(r100b) + &
         w(r102b) + &
         w(r106b) 

      !cdot(sSXCH2) = 0.0_pr


        cdot(sCH2O) = 0.0_pr + &
         w(r21) + &
         w(r27) + &
         w(r29) - &
         w(r45) - &
         w(r46) - &
         w(r47) - &
         w(r48) - &
         w(r49) + &
         w(r51) + &
         w(r52) + &
         w(r53) + &
         w(r54) + &
         w(r55) + &
         w(r57f) + &
         w(r58) + &
         w(r59) + &
         w(r66) + &
         w(r91) + &
         w(r92) + &
         w(r116) - &
         w(r57b) 


        cdot(sCO) = 0.0_pr + &
         w(r22) - &
         w(r37) - &
         w(r38f) - &
         w(r39) + &
         w(r40f) + &
         w(r41) + &
         w(r42) + &
         w(r43) + &
         w(r44) + &
         w(r58) + &
         w(r59) + &
         w(r84) + &
         w(r85f) + &
         w(r86) + &
         w(r89) + &
         w(r91) + &
         w(r94) + &
         w(r38b) - &
         w(r40b) - &
         w(r85b) 


        cdot(sCH3) = 0.0_pr - &
         w(r23f) - &
         w(r24) - 2_pr * &
         w(r25f) - 2_pr * &
         w(r26f) - &
         w(r27) - &
         w(r28f) - &
         w(r29) - &
         w(r30) - &
         w(r31) - &
         w(r32) - &
         w(r33f) + &
         w(r34f) + &
         w(r35f) + &
         w(r36) - &
         w(r44) + &
         w(r56f) - &
         w(r63f) + &
         w(r70) - &
         w(r82) - &
         w(r84) + &
         w(r85f) + &
         w(r87) - &
         w(r89) + &
         w(r94) + &
         w(r96) + &
         w(r98) + &
         w(r103f) + &
         w(r110f) + &
         w(r23b) + 2_pr * &
         w(r25b) + 2_pr * &
         w(r26b) + &
         w(r28b) + &
         w(r33b) - &
         w(r34b) - &
         w(r35b) - &
         w(r56b) + &
         w(r63b) - &
         w(r85b) - &
         w(r103b) - &
         w(r110b) 


        cdot(sCH4) = 0.0_pr + &
         w(r24) + &
         w(r32) + &
         w(r33f) - &
         w(r34f) - &
         w(r35f) - &
         w(r36) + &
         w(r44) + &
         w(r82) - &
         w(r33b) + &
         w(r34b) + &
         w(r35b) 

      !cdot(sC2H5) = 0.0_pr


        cdot(sC2H6) = 0.0_pr + &
         w(r26f) + &
         w(r74) - &
         w(r80) - &
         w(r81f) - &
         w(r82) - &
         w(r83) - &
         w(r26b) + &
         w(r81b) 

      !cdot(sCH3O) = 0.0_pr


        cdot(sCO2) = 0.0_pr + &
         w(r37) + &
         w(r38f) + &
         w(r39) - &
         w(r58) - &
         w(r38b) 

      !cdot(sHCO) = 0.0_pr

      !cdot(sCH2OH) = 0.0_pr


        cdot(sC2H2) = 0.0_pr - &
         w(r59) - &
         w(r60f) + &
         w(r61) + &
         w(r64) + &
         w(r65f) + &
         w(r71) + &
         w(r60b) - &
         w(r65b) 


        cdot(sCH2CO) = 0.0_pr + &
         w(r60f) - &
         w(r84) - &
         w(r85f) - &
         w(r86) + &
         w(r88) + &
         w(r90f) + &
         w(r96) - &
         w(r60b) + &
         w(r85b) - &
         w(r90b) 

      !cdot(sC2H3) = 0.0_pr

      !cdot(sCH2CHO) = 0.0_pr


        cdot(sC3H6) = 0.0_pr + &
         w(r63f) - &
         w(r95f) - &
         w(r96) - &
         w(r97f) - &
         w(r98) - &
         w(r99) + &
         w(r100f) + &
         w(r101) + &
         w(r115) - &
         w(r63b) + &
         w(r95b) + &
         w(r97b) - &
         w(r100b) 


        cdot(sC2H4) = 0.0_pr + &
         w(r67) - &
         w(r68) - &
         w(r69) - &
         w(r70) - &
         w(r71) - &
         w(r72f) - &
         w(r73f) + &
         w(r75) + &
         w(r76) + &
         w(r77f) + &
         w(r98) + &
         w(r103f) + &
         w(r72b) + &
         w(r73b) - &
         w(r77b) - &
         w(r103b) 


        cdot(sC2H4O) = 0.0_pr + &
         w(r69) - &
         w(r94) 

      !cdot(sNXC3H7) = 0.0_pr

      !cdot(sIXC3H7) = 0.0_pr

      !cdot(sC3H6OOH) = 0.0_pr


        cdot(sC3H8) = 0.0_pr - &
         w(r104) - &
         w(r105) - &
         w(r106f) - &
         w(r107f) - &
         w(r108) - &
         w(r109) - &
         w(r110f) - &
         w(r111) - &
         w(r112) - &
         w(r113f) + &
         w(r106b) + &
         w(r107b) + &
         w(r110b) + &
         w(r113b) 


        cdot(sOC3H5OOH) = 0.0_pr + &
         w(r114) - &
         w(r116) 


      end subroutine get_production_rates

! --- Actual reactions --- !
      subroutine reaction_expressions
      implicit none

      reacexp(1) = 'H2 + M1 -> 2 H + M1'
      reacexp(2) = 'H2 + O -> OH + H'
      reacexp(3) = 'H + O + M4 -> OH + M4'
      reacexp(4) = 'H + OH + M2 -> H2O + M2'
      reacexp(5) = 'O + OH + M5 -> HO2 + M5'
      reacexp(6) = '2 OH + M7 -> H2O2 + M7'
      reacexp(7) = 'H2 + OH -> H2O + H'
      reacexp(8) = 'H2O + O -> 2 OH'
      reacexp(9) = 'H + O2 -> OH + O'
      reacexp(10) = 'H + O2 + M6 -> HO2 + M6'
      reacexp(11) = 'HO2 + O -> OH + O2'
      reacexp(12) = 'HO2 + H -> 2 OH'
      reacexp(13) = 'HO2 + OH -> H2O + O2'
      reacexp(14) = 'HO2 + OH -> H2O + O2'
      reacexp(15) = 'HO2 + H -> H2O + O'
      reacexp(16) = '2 HO2 -> H2O2 + O2'
      reacexp(17) = '2 HO2 -> H2O2 + O2'
      reacexp(18) = 'HO2 + H -> H2 + O2'
      reacexp(19) = 'H2O2 + OH -> H2O + HO2'
      reacexp(20) = 'H2O2 + OH -> H2O + HO2'
      reacexp(21) = 'S-CH2 + OH -> CH2O + H'
      reacexp(22) = 'S-CH2 + O2 -> CO + OH + H'
      reacexp(23) = 'CH3 + OH -> S-CH2 + H2O'
      reacexp(24) = 'HO2 + CH3 -> O2 + CH4'
      reacexp(25) = '2 CH3 -> C2H5 + H'
      reacexp(26) = '2 CH3 + M12 -> C2H6 + M12'
      reacexp(27) = 'CH3 + O -> CH2O + H'
      reacexp(28) = 'CH3 + H -> S-CH2 + H2'
      reacexp(29) = 'CH3 + O2 -> CH2O + OH'
      reacexp(30) = 'CH3 + HO2 -> CH3O + OH'
      reacexp(31) = 'CH3 + O2 -> CH3O + O'
      reacexp(32) = 'H2O2 + CH3 -> HO2 + CH4'
      reacexp(33) = 'H + CH3 + M11 -> CH4 + M11'
      reacexp(34) = 'CH4 + H -> H2 + CH3'
      reacexp(35) = 'CH4 + OH -> H2O + CH3'
      reacexp(36) = 'CH4 + O -> CH3 + OH'
      reacexp(37) = 'CO + O + M8 -> CO2 + M8'
      reacexp(38) = 'CO + OH -> CO2 + H'
      reacexp(39) = 'CO + HO2 -> CO2 + OH'
      reacexp(40) = 'HCO + M9 -> CO + H + M9'
      reacexp(41) = 'HCO + OH -> CO + H2O'
      reacexp(42) = 'HCO + O2 -> CO + HO2'
      reacexp(43) = 'HCO + H -> CO + H2'
      reacexp(44) = 'HCO + CH3 -> CO + CH4'
      reacexp(45) = 'CH2O + O -> HCO + OH'
      reacexp(46) = 'CH2O + OH -> HCO + H2O'
      reacexp(47) = 'CH2O + HO2 -> HCO + H2O2'
      reacexp(48) = 'CH2O + H -> HCO + H2'
      reacexp(49) = 'CH2O + M10 -> HCO + H + M10'
      reacexp(50) = 'CH3O + M21 -> CH2OH + M21'
      reacexp(51) = 'CH3O + M14 -> CH2O + H + M14'
      reacexp(52) = 'CH3O + O2 -> CH2O + HO2'
      reacexp(53) = 'CH2OH + H -> CH2O + H2'
      reacexp(54) = 'CH2OH + OH -> CH2O + H2O'
      reacexp(55) = 'CH2OH + M20 -> CH2O + H + M20'
      reacexp(56) = 'CH2OH + H -> CH3 + OH'
      reacexp(57) = 'CH2OH + O2 -> CH2O + HO2'
      reacexp(58) = 'S-CH2 + CO2 -> CO + CH2O'
      reacexp(59) = 'C2H2 + O2 -> CH2O + CO'
      reacexp(60) = 'C2H2 + OH -> CH2CO + H'
      reacexp(61) = 'C2H3 + O2 -> C2H2 + HO2'
      reacexp(62) = 'C2H3 + O2 -> CH2CHO + O'
      reacexp(63) = 'C2H3 + CH3 + M30 -> C3H6 + M30'
      reacexp(64) = 'C2H3 + H -> C2H2 + H2'
      reacexp(65) = 'C2H3 + M19 -> C2H2 + H + M19'
      reacexp(66) = 'C2H3 + O2 -> CH2O + HCO'
      reacexp(67) = 'H + C2H3 + M17 -> C2H4 + M17'
      reacexp(68) = 'C2H4 + O -> CH2CHO + H'
      reacexp(69) = 'C2H4 + HO2 -> C2H4O + OH'
      reacexp(70) = 'C2H4 + O -> CH3 + HCO'
      reacexp(71) = 'C2H4 + M18 -> C2H2 + H2 + M18'
      reacexp(72) = 'C2H4 + H -> C2H3 + H2'
      reacexp(73) = 'C2H4 + OH -> C2H3 + H2O'
      reacexp(74) = 'H + C2H5 + M15 -> C2H6 + M15'
      reacexp(75) = 'C2H5 + O2 -> C2H4 + HO2'
      reacexp(76) = 'C2H5 + H -> C2H4 + H2'
      reacexp(77) = 'C2H5 + M16 -> C2H4 + H + M16'
      reacexp(78) = 'C2H6 + OH -> C2H5 + H2O'
      reacexp(79) = 'C2H6 + H -> C2H5 + H2'
      reacexp(80) = 'C2H6 + CH3 -> C2H5 + CH4'
      reacexp(81) = 'C2H6 + O -> C2H5 + OH'
      reacexp(82) = 'CH2CO + CH3 -> C2H5 + CO'
      reacexp(83) = 'CH2CO + H -> CH3 + CO'
      reacexp(84) = 'CH2CO + OH -> CH2OH + CO'
      reacexp(85) = 'CH2CHO + H -> CH3 + HCO'
      reacexp(86) = 'CH2CHO + OH -> CH2CO + H2O'
      reacexp(87) = 'CH2CHO + CH3 -> C2H5 + CO + H'
      reacexp(88) = 'CH2CHO -> CH2CO + H'
      reacexp(89) = 'CH2CHO + O2 -> CH2O + CO + OH'
      reacexp(90) = 'CH2CHO + HO2 -> CH2O + HCO + OH'
      reacexp(91) = 'C2H4O + HO2 -> CH3 + CO + H2O2'
      reacexp(92) = 'H + C3H6 + M32 -> N-C3H7 + M32'
      reacexp(93) = 'C3H6 + O -> CH2CO + CH3 + H'
      reacexp(94) = 'C3H6 + H + M31 -> I-C3H7 + M31'
      reacexp(95) = 'C3H6 + H -> C2H4 + CH3'
      reacexp(96) = 'C3H6 + O -> C2H5 + HCO'
      reacexp(97) = 'I-C3H7 + O2 -> C3H6 + HO2'
      reacexp(98) = 'N-C3H7 + O2 -> C3H6 + HO2'
      reacexp(99) = 'N-C3H7 + O2 -> C3H6OOH'
      reacexp(100) = 'N-C3H7 + M0 -> CH3 + C2H4 + M0'
      reacexp(101) = 'C3H8 + HO2 -> N-C3H7 + H2O2'
      reacexp(102) = 'C3H8 + HO2 -> I-C3H7 + H2O2'
      reacexp(103) = 'C3H8 + O2 -> I-C3H7 + HO2'
      reacexp(104) = 'C3H8 + H -> N-C3H7 + H2'
      reacexp(105) = 'C3H8 + O2 -> N-C3H7 + HO2'
      reacexp(106) = 'C3H8 + O -> N-C3H7 + OH'
      reacexp(107) = 'C3H8 + M0 -> CH3 + C2H5 + M0'
      reacexp(108) = 'C3H8 + OH -> N-C3H7 + H2O'
      reacexp(109) = 'C3H8 + O -> I-C3H7 + OH'
      reacexp(110) = 'C3H8 + H -> I-C3H7 + H2'
      reacexp(111) = 'C3H6OOH + O2 -> OC3H5OOH + OH'
      reacexp(112) = 'C3H6OOH -> C3H6 + HO2'
      reacexp(113) = 'OC3H5OOH -> CH2CHO + CH2O + OH'
      reacexp(114) = 'Reverse of H2 + O -> OH + H'
      reacexp(115) = 'Reverse of H + OH + M2 -> H2O + M2'
      reacexp(116) = 'Reverse of O + OH + M5 -> HO2 + M5'
      reacexp(117) = 'Reverse of 2 OH + M7 -> H2O2 + M7'
      reacexp(118) = 'Reverse of H2 + OH -> H2O + H'
      reacexp(119) = 'Reverse of H2O + O -> 2 OH'
      reacexp(120) = 'Reverse of H + O2 -> OH + O'
      reacexp(121) = 'Reverse of H + O2 + M6 -> HO2 + M6'
      reacexp(122) = 'Reverse of HO2 + OH -> H2O + O2'
      reacexp(123) = 'Reverse of HO2 + OH -> H2O + O2'
      reacexp(124) = 'Reverse of CH3 + OH -> S-CH2 + H2O'
      reacexp(125) = 'Reverse of 2 CH3 -> C2H5 + H'
      reacexp(126) = 'Reverse of 2 CH3 + M12 -> C2H6 + M12'
      reacexp(127) = 'Reverse of CH3 + H -> S-CH2 + H2'
      reacexp(128) = 'Reverse of H + CH3 + M11 -> CH4 + M11'
      reacexp(129) = 'Reverse of CH4 + H -> H2 + CH3'
      reacexp(130) = 'Reverse of CH4 + OH -> H2O + CH3'
      reacexp(131) = 'Reverse of CO + OH -> CO2 + H'
      reacexp(132) = 'Reverse of HCO + M9 -> CO + H + M9'
      reacexp(133) = 'Reverse of CH3O + M21 -> CH2OH + M21'
      reacexp(134) = 'Reverse of CH2OH + H -> CH3 + OH'
      reacexp(135) = 'Reverse of CH2OH + O2 -> CH2O + HO2'
      reacexp(136) = 'Reverse of C2H2 + OH -> CH2CO + H'
      reacexp(137) = 'Reverse of C2H3 + CH3 + M30 -> C3H6 + M30'
      reacexp(138) = 'Reverse of C2H3 + M19 -> C2H2 + H + M19'
      reacexp(139) = 'Reverse of C2H4 + H -> C2H3 + H2'
      reacexp(140) = 'Reverse of C2H4 + OH -> C2H3 + H2O'
      reacexp(141) = 'Reverse of C2H5 + M16 -> C2H4 + H + M16'
      reacexp(142) = 'Reverse of C2H6 + H -> C2H5 + H2'
      reacexp(143) = 'Reverse of CH2CO + H -> CH3 + CO'
      reacexp(144) = 'Reverse of CH2CHO -> CH2CO + H'
      reacexp(145) = 'Reverse of H + C3H6 + M32 -> N-C3H7 + M32'
      reacexp(146) = 'Reverse of C3H6 + H + M31 -> I-C3H7 + M31'
      reacexp(147) = 'Reverse of I-C3H7 + O2 -> C3H6 + HO2'
      reacexp(148) = 'Reverse of N-C3H7 + O2 -> C3H6OOH'
      reacexp(149) = 'Reverse of N-C3H7 + M0 -> CH3 + C2H4 + M0'
      reacexp(150) = 'Reverse of C3H8 + O2 -> I-C3H7 + HO2'
      reacexp(151) = 'Reverse of C3H8 + H -> N-C3H7 + H2'
      reacexp(152) = 'Reverse of C3H8 + M0 -> CH3 + C2H5 + M0'
      reacexp(153) = 'Reverse of C3H8 + H -> I-C3H7 + H2'

      end subroutine reaction_expressions

! --- Forward/Backward link --- !
      subroutine reverse_reactions
      implicit none

      fofb = 0

! Attach corresponding forward reaction to each backward reaction
      fofb(114) = 2
      fofb(115) = 4
      fofb(116) = 5
      fofb(117) = 6
      fofb(118) = 7
      fofb(119) = 8
      fofb(120) = 9
      fofb(121) = 10
      fofb(122) = 13
      fofb(123) = 14
      fofb(124) = 23
      fofb(125) = 25
      fofb(126) = 26
      fofb(127) = 28
      fofb(128) = 33
      fofb(129) = 34
      fofb(130) = 35
      fofb(131) = 38
      fofb(132) = 40
      fofb(133) = 50
      fofb(134) = 56
      fofb(135) = 57
      fofb(136) = 60
      fofb(137) = 63
      fofb(138) = 65
      fofb(139) = 72
      fofb(140) = 73
      fofb(141) = 77
      fofb(142) = 79
      fofb(143) = 83
      fofb(144) = 88
      fofb(145) = 92
      fofb(146) = 94
      fofb(147) = 97
      fofb(148) = 99
      fofb(149) = 100
      fofb(150) = 103
      fofb(151) = 104
      fofb(152) = 107
      fofb(153) = 110

      end subroutine reverse_reactions

! --- 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) :: k
      real(pr), dimension(24) :: M

      real(pr) :: C3H6OOH_denom1&
         , C3H6OOH_ct1&
         , C3H6OOH_denom2&
         , C3H6OOH_ct2
      real(pr) :: NXC3H7_denom1&
         , NXC3H7_ct1&
         , NXC3H7_denom2&
         , NXC3H7_ct2&
         , NXC3H7_C3H6OOH&
         , NXC3H7_C3H6OOH_coeff
      real(pr) :: CH2OH_denom1&
         , CH2OH_ct1&
         , CH2OH_denom2&
         , CH2OH_ct2
      real(pr) :: CH3O_denom1&
         , CH3O_ct1&
         , CH3O_denom2&
         , CH3O_ct2&
         , CH3O_CH2OH&
         , CH3O_CH2OH_coeff


  ! c(sCH3O) c(sCH2OH) (coupled)  --------------------

  ! Primary denominators-----------------------

	 CH3O_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r50f) * M(mM21)  & 
               + k(r51) * M(mM14)  & 
               + k(r52) *  c(sO2) & 
                )

	 CH2OH_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r50b) * M(mM21)  & 
               + k(r53) *  c(sH) & 
               + k(r54) *  c(sOH) & 
               + k(r55) * M(mM20)  & 
               + k(r56f) *  c(sH) & 
               + k(r57f) *  c(sO2) & 
                )


  ! Primary constant parts -----------------------

	 CH3O_ct1 = ( 0.0_pr &
                + k(r30) * c(sCH3) * c(sHO2) &
                + k(r31) * c(sCH3) * c(sO2)                )

	 CH2OH_ct1 = ( 0.0_pr &
                + k(r56b) * c(sOH) * c(sCH3) &
                + k(r57b) * c(sHO2) * c(sCH2O) &
                + k(r86) * c(sCH2CO) * c(sOH)                )


  ! CH3O ---------------------------------------

	 CH3O_denom2 = tiny(1.0_pr) + ( CH3O_denom1                 )

	 CH3O_ct2 = ( CH3O_ct1 &
                ) / CH3O_denom2  

	 CH3O_CH2OH = ( 0.0_pr &
                + k(r50b)  * M(mM21) &
                ) / CH3O_denom2 



   	 CH3O_CH2OH_coeff = ( 0.0_pr &
                + k(r50f)  * M(mM21) &
                )


  ! CH2OH ---------------------------------------

	 CH2OH_denom2 = tiny(1.0_pr) + ( CH2OH_denom1  &
                - CH3O_CH2OH_coeff * CH3O_CH2OH                )

	 CH2OH_ct2 = ( CH2OH_ct1 &
                + CH3O_CH2OH_coeff * CH3O_ct2 &
                ) / CH2OH_denom2  




  ! Reconstruction ------------------------------------

	 cqss(sCH2OH-nspec) = ( CH2OH_ct2                 ) 

	 cqss(sCH3O-nspec) = ( CH3O_ct2 &
                + CH3O_CH2OH * cqss(sCH2OH-nspec)                ) 

  ! cqss(sC2H3) (uncoupled) --------------------

	 cqss(sC2H3-nspec) = ( 0.0_pr &
                + k(r63b) * c(sC3H6) &
                + k(r65b) * c(sH) * c(sC2H2) &
                + k(r72f) * c(sC2H4) * c(sH) &
                + k(r73f) * c(sC2H4) * c(sOH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r61) *  c(sO2) &
               + k(r62) *  c(sO2) &
               + k(r63f) *  c(sCH3) &
               + k(r64) *  c(sH) &
               + k(r65f)  &
               + k(r66) *  c(sO2) &
               + k(r67) *  c(sH) * M(mM17) &
               + k(r72b) *  c(sH2) &
               + k(r73b) *  c(sH2O)                 ) )


  ! cqss(sO) (uncoupled) --------------------

	 cqss(sO-nspec) = ( 0.0_pr &
                + k(r2b) * c(sH) * c(sOH) &
                + k(r5b) * c(sHO2) * M(mM5) &
                + k(r8b) * c(sOH) * c(sOH) &
                + k(r9f) * c(sH) * c(sO2) &
                + k(r15) * c(sHO2) * c(sH) &
                + k(r31) * c(sCH3) * c(sO2) &
                + k(r62) * cqss(sC2H3-nspec) *c(sO2) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r2f) *  c(sH2) &
               + k(r3) *  c(sH) * M(mM4) &
               + k(r5f) *  c(sOH) * M(mM5) &
               + k(r8f) *  c(sH2O) &
               + k(r9b) *  c(sOH) &
               + k(r11) *  c(sHO2) &
               + k(r27) *  c(sCH3) &
               + k(r36) *  c(sCH4) &
               + k(r37) *  c(sCO) &
               + k(r45) *  c(sCH2O) &
               + k(r68) *  c(sC2H4) &
               + k(r70) *  c(sC2H4) &
               + k(r83) *  c(sC2H6) &
               + k(r96) *  c(sC3H6) &
               + k(r99) *  c(sC3H6) &
               + k(r109) *  c(sC3H8) &
               + k(r112) *  c(sC3H8)                 ) )


  ! cqss(sCH2CHO) (uncoupled) --------------------

	 cqss(sCH2CHO-nspec) = ( 0.0_pr &
                + k(r62) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r68) * c(sC2H4) * cqss(sO-nspec) &
                + k(r90b) * c(sH) * c(sCH2CO) &
                + k(r116) * c(sOC3H5OOH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r87) *  c(sH) &
               + k(r88) *  c(sOH) &
               + k(r89) *  c(sCH3) &
               + k(r90f)  &
               + k(r91) *  c(sO2) &
               + k(r92) *  c(sHO2)                 ) )


  ! cqss(sHCO) (uncoupled) --------------------

	 cqss(sHCO-nspec) = ( 0.0_pr &
                + k(r40b) * c(sH) * c(sCO) * M(mM9) &
                + k(r45) * c(sCH2O) * cqss(sO-nspec) &
                + k(r46) * c(sCH2O) * c(sOH) &
                + k(r47) * c(sCH2O) * c(sHO2) &
                + k(r48) * c(sCH2O) * c(sH) &
                + k(r49) * c(sCH2O) &
                + k(r66) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r70) * c(sC2H4) * cqss(sO-nspec) &
                + k(r87) * cqss(sCH2CHO-nspec) *c(sH) &
                + k(r92) * cqss(sCH2CHO-nspec) *c(sHO2) &
                + k(r99) * c(sC3H6) * cqss(sO-nspec) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r40f) * M(mM9)  &
               + k(r41) *  c(sOH) &
               + k(r42) *  c(sO2) &
               + k(r43) *  c(sH) &
               + k(r44) *  c(sCH3)                 ) )


  ! cqss(sIXC3H7) (uncoupled) --------------------

	 cqss(sIXC3H7-nspec) = ( 0.0_pr &
                + k(r97f) * c(sC3H6) * c(sH) &
                + k(r100b) * c(sHO2) * c(sC3H6) &
                + k(r105) * c(sC3H8) * c(sHO2) &
                + k(r106f) * c(sC3H8) * c(sO2) &
                + k(r112) * c(sC3H8) * cqss(sO-nspec) &
                + k(r113f) * c(sC3H8) * c(sH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r97b)  &
               + k(r100f) *  c(sO2) &
               + k(r106b) *  c(sHO2) &
               + k(r113b) *  c(sH2)                 ) )


  ! c(sNXC3H7) c(sC3H6OOH) (coupled)  --------------------

  ! Primary denominators-----------------------

	 NXC3H7_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r95b)  & 
               + k(r101) *  c(sO2) & 
               + k(r102f) *  c(sO2) & 
               + k(r103f)  & 
               + k(r107b) *  c(sH2) & 
                )

	 C3H6OOH_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r102b)  & 
               + k(r114) *  c(sO2) & 
               + k(r115)  & 
                )


  ! Primary constant parts -----------------------

	 NXC3H7_ct1 = ( 0.0_pr &
                + k(r95f) * c(sH) * c(sC3H6) &
                + k(r103b) * c(sC2H4) * c(sCH3) &
                + k(r104) * c(sC3H8) * c(sHO2) &
                + k(r107f) * c(sC3H8) * c(sH) &
                + k(r108) * c(sC3H8) * c(sO2) &
                + k(r109) * c(sC3H8) * cqss(sO-nspec) &
                + k(r111) * c(sC3H8) * c(sOH)                )

	 C3H6OOH_ct1 = ( 0.0_pr                )


  ! NXC3H7 ---------------------------------------

	 NXC3H7_denom2 = tiny(1.0_pr) + ( NXC3H7_denom1                 )

	 NXC3H7_ct2 = ( NXC3H7_ct1 &
                ) / NXC3H7_denom2  

	 NXC3H7_C3H6OOH = ( 0.0_pr &
                + k(r102b)  &
                ) / NXC3H7_denom2 



   	 NXC3H7_C3H6OOH_coeff = ( 0.0_pr &
                + k(r102f)  * c(sO2) &
                )


  ! C3H6OOH ---------------------------------------

	 C3H6OOH_denom2 = tiny(1.0_pr) + ( C3H6OOH_denom1  &
                - NXC3H7_C3H6OOH_coeff * NXC3H7_C3H6OOH                )

	 C3H6OOH_ct2 = ( C3H6OOH_ct1 &
                + NXC3H7_C3H6OOH_coeff * NXC3H7_ct2 &
                ) / C3H6OOH_denom2  




  ! Reconstruction ------------------------------------

	 cqss(sC3H6OOH-nspec) = ( C3H6OOH_ct2                 ) 

	 cqss(sNXC3H7-nspec) = ( NXC3H7_ct2 &
                + NXC3H7_C3H6OOH * cqss(sC3H6OOH-nspec)                ) 

  ! cqss(sC2H5) (uncoupled) --------------------

	 cqss(sC2H5-nspec) = ( 0.0_pr &
                + k(r25f) * c(sCH3) * c(sCH3) &
                + k(r77b) * c(sH) * c(sC2H4) &
                + k(r80) * c(sC2H6) * c(sOH) &
                + k(r81f) * c(sC2H6) * c(sH) &
                + k(r82) * c(sC2H6) * c(sCH3) &
                + k(r83) * c(sC2H6) * cqss(sO-nspec) &
                + k(r84) * c(sCH2CO) * c(sCH3) &
                + k(r89) * cqss(sCH2CHO-nspec) *c(sCH3) &
                + k(r99) * c(sC3H6) * cqss(sO-nspec) &
                + k(r110f) * c(sC3H8) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r25b) *  c(sH) &
               + k(r74) *  c(sH) &
               + k(r75) *  c(sO2) &
               + k(r76) *  c(sH) &
               + k(r77f)  &
               + k(r81b) *  c(sH2) &
               + k(r110b) *  c(sCH3)                 ) )


  ! cqss(sSXCH2) (uncoupled) --------------------

	 cqss(sSXCH2-nspec) = ( 0.0_pr &
                + k(r23f) * c(sCH3) * c(sOH) &
                + k(r28f) * c(sCH3) * c(sH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r21) *  c(sOH) &
               + k(r22) *  c(sO2) &
               + k(r23b) *  c(sH2O) &
               + k(r28b) *  c(sH2) &
               + k(r58) *  c(sCO2)                 ) )


      end subroutine get_QSS

      !subroutine interne au module pour AVBP
      subroutine INTERNALAVBP(P, T, c, WDOT)
      implicit none

      real(pr), dimension(nspec) :: c, WDOT
      real(pr), dimension(nqss) :: cqss
      real(pr), dimension(nreac) :: w,k
      real(pr), dimension(24) :: m

      real(pr) :: T,P

      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)
      CALL get_production_rates(WDOT,w)
      end subroutine INTERNALAVBP

end module mod_C3H8_SanDiego_ARC21_153_11_QM

!subroutine externe au module pour AVBP
SUBROUTINE C3H8_SanDiego_ARC21_153_11_QM ( nnode,neqs,rhoinv,press,tempe,w_spec,wmol,source_spec )
      USE mod_param_defs
      USE mod_C3H8_SanDiego_ARC21_153_11_QM, ONLY: INTERNALAVBP
      IMPLICIT NONE
      ! IN/OUT
      INTEGER, INTENT(IN) :: nnode, neqs
      REAL(pr), DIMENSION(1:neqs), INTENT(IN) :: wmol
      REAL(pr), DIMENSION(1:nnode), INTENT(IN) :: rhoinv,press,tempe
      REAL(pr), DIMENSION(1:neqs,1:nnode), INTENT(IN) :: w_spec
      REAL(pr), DIMENSION(1:neqs,1:nnode), INTENT(OUT) ::source_spec
      ! LOCAL
      REAL(pr), DIMENSION(1:neqs) :: C
      REAL(pr), DIMENSION(1:neqs) :: source_spec_buf
      REAL(pr) :: P,T
      INTEGER :: n,k

      DO n=1,nnode
        P = press(n)
        T = tempe(n)

        !Activity concentrations
        DO k=1,neqs
          C(k) = min( w_spec(k,n), one/rhoinv(n) )
          C(k) = C(k) / wmol(k)
          C(k) = max(C(k),1.0e-10_pr)
        END DO

        !   The internal YARC routine is called
        CALL INTERNALAVBP ( P, T, C, source_spec_buf )

        !   Conversion from mol/m3/s to kg/m3/s
        DO k=1,neqs
          source_spec(k,n) = source_spec_buf(k) * wmol(k)
        END DO
      END DO
END SUBROUTINE C3H8_SanDiego_ARC21_153_11_QM