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

!--------------------------------------------------------------------------------------------------
!     MODULE mod_C3H8_20_114_10_QM
!>    Module for calculating the analytical source terms for the
!C3H8_21_170_12_QM scheme.
!!    @authors Q. Male
!!    @date    10/12/2018
!--------------------------------------------------------------------------------------------------
module mod_C3H8_20_114_10_QM
      USE mod_param_defs
      implicit none

      integer, parameter :: nspec = 20
      integer, parameter :: nreac = 114
      integer, parameter :: isc_T = 1
      integer, parameter :: neq = nspec + 1

      ! QSS variables
      ! Number of QSS species
      integer, parameter :: nqss = 10
      ! 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 = 31
      ! 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 :: sO = 3
      integer, parameter :: sOH = 4
      integer, parameter :: sH = 5
      integer, parameter :: sH2O = 6
      integer, parameter :: sH2O2 = 7
      integer, parameter :: sO2 = 8
      integer, parameter :: sHO2 = 9
      integer, parameter :: sCH2O = 10
      integer, parameter :: sCO2 = 11
      integer, parameter :: sCO = 12
      integer, parameter :: sCH3 = 13
      integer, parameter :: sC2H6 = 14
      integer, parameter :: sC3H6 = 15
      integer, parameter :: sC2H4 = 16
      integer, parameter :: sCH2CO = 17
      integer, parameter :: sC3H8 = 18
      integer, parameter :: sOC3H5OOH = 19
      integer, parameter :: sC4H8 = 20
      integer, parameter :: sTXCH2 = 21
      integer, parameter :: sSXCH2 = 22
      integer, parameter :: sHCO = 23
      integer, parameter :: sC2H3 = 24
      integer, parameter :: sCH2CHO = 25
      integer, parameter :: sC2H5 = 26
      integer, parameter :: sC3H5 = 27
      integer, parameter :: sIXC3H7 = 28
      integer, parameter :: sNXC3H7 = 29
      integer, parameter :: sC3H6OOH = 30

      ! Index of reactions
      integer, parameter :: r1f = 1
      integer, parameter :: r2f = 2
      integer, parameter :: r3f = 3
      integer, parameter :: r4f = 4
      integer, parameter :: r5f = 5
      integer, parameter :: r6f = 6
      integer, parameter :: r7f = 7
      integer, parameter :: r8f = 8
      integer, parameter :: r9f = 9
      integer, parameter :: r10 = 10
      integer, parameter :: r11 = 11
      integer, parameter :: r12 = 12
      integer, parameter :: r13 = 13
      integer, parameter :: r14 = 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 :: r22f = 22
      integer, parameter :: r23 = 23
      integer, parameter :: r24f = 24
      integer, parameter :: r25 = 25
      integer, parameter :: r26 = 26
      integer, parameter :: r27 = 27
      integer, parameter :: r28f = 28
      integer, parameter :: r29 = 29
      integer, parameter :: r30 = 30
      integer, parameter :: r31 = 31
      integer, parameter :: r32f = 32
      integer, parameter :: r33 = 33
      integer, parameter :: r34 = 34
      integer, parameter :: r35 = 35
      integer, parameter :: r36 = 36
      integer, parameter :: r37 = 37
      integer, parameter :: r38 = 38
      integer, parameter :: r39 = 39
      integer, parameter :: r40 = 40
      integer, parameter :: r41 = 41
      integer, parameter :: r42 = 42
      integer, parameter :: r43f = 43
      integer, parameter :: r44 = 44
      integer, parameter :: r45f = 45
      integer, parameter :: r46f = 46
      integer, parameter :: r47 = 47
      integer, parameter :: r48 = 48
      integer, parameter :: r49 = 49
      integer, parameter :: r50 = 50
      integer, parameter :: r51 = 51
      integer, parameter :: r52 = 52
      integer, parameter :: r53 = 53
      integer, parameter :: r54 = 54
      integer, parameter :: r55 = 55
      integer, parameter :: r56 = 56
      integer, parameter :: r57 = 57
      integer, parameter :: r58f = 58
      integer, parameter :: r59 = 59
      integer, parameter :: r60 = 60
      integer, parameter :: r61f = 61
      integer, parameter :: r62 = 62
      integer, parameter :: r63f = 63
      integer, parameter :: r64 = 64
      integer, parameter :: r65 = 65
      integer, parameter :: r66 = 66
      integer, parameter :: r67f = 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 :: r77 = 77
      integer, parameter :: r78 = 78
      integer, parameter :: r79 = 79
      integer, parameter :: r80f = 80
      integer, parameter :: r81 = 81
      integer, parameter :: r82 = 82
      integer, parameter :: r83 = 83
      integer, parameter :: r84 = 84
      integer, parameter :: r85 = 85
      integer, parameter :: r86 = 86
      integer, parameter :: r87 = 87
      integer, parameter :: r88 = 88
      integer, parameter :: r89f = 89
      integer, parameter :: r90 = 90
      integer, parameter :: r1b = 91
      integer, parameter :: r2b = 92
      integer, parameter :: r3b = 93
      integer, parameter :: r4b = 94
      integer, parameter :: r5b = 95
      integer, parameter :: r6b = 96
      integer, parameter :: r7b = 97
      integer, parameter :: r8b = 98
      integer, parameter :: r9b = 99
      integer, parameter :: r22b = 100
      integer, parameter :: r24b = 101
      integer, parameter :: r28b = 102
      integer, parameter :: r32b = 103
      integer, parameter :: r43b = 104
      integer, parameter :: r45b = 105
      integer, parameter :: r46b = 106
      integer, parameter :: r58b = 107
      integer, parameter :: r61b = 108
      integer, parameter :: r63b = 109
      integer, parameter :: r67b = 110
      integer, parameter :: r72b = 111
      integer, parameter :: r73b = 112
      integer, parameter :: r80b = 113
      integer, parameter :: r89b = 114

      ! Index of third bodies
      integer, parameter :: mM31 = 1
      integer, parameter :: mM7 = 2
      integer, parameter :: mM11 = 3
      integer, parameter :: mM0 = 4
      integer, parameter :: mM30 = 5
      integer, parameter :: mM8 = 6
      integer, parameter :: mM29 = 7
      integer, parameter :: mM6 = 8
      integer, parameter :: mM16 = 9
      integer, parameter :: mM15 = 10
      integer, parameter :: mM12 = 11
      integer, parameter :: mM5 = 12
      integer, parameter :: mM2 = 13
      integer, parameter :: mM28 = 14

      contains

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

      ! if (irank.eq.iroot) print*, 'Using mechanism ./Spec_30_Reac_114_QSS_10.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

      ! 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) = 1

      ! Name of group of species
      ppname(1) = trim(gname(sN2))
      ppname(2) = trim(gname(sH2))
      ppname(3) = trim(gname(sO))
      ppname(4) = trim(gname(sOH))
      ppname(5) = trim(gname(sH))
      ppname(6) = trim(gname(sH2O))
      ppname(7) = trim(gname(sH2O2))
      ppname(8) = trim(gname(sO2))
      ppname(9) = trim(gname(sHO2))
      ppname(10) = trim(gname(sTXCH2))
      ppname(11) = trim(gname(sCH2O))
      ppname(12) = trim(gname(sCO2))
      ppname(13) = trim(gname(sCO))
      ppname(14) = trim(gname(sSXCH2))
      ppname(15) = trim(gname(sCH3))
      ppname(16) = trim(gname(sC2H6))
      ppname(17) = trim(gname(sHCO))
      ppname(18) = trim(gname(sC2H3))
      ppname(19) = trim(gname(sCH2CHO))
      ppname(20) = trim(gname(sC3H6))
      ppname(21) = trim(gname(sC2H4))
      ppname(22) = trim(gname(sC2H5))
      ppname(23) = trim(gname(sCH2CO))
      ppname(24) = trim(gname(sC3H5))
      ppname(25) = trim(gname(sIXC3H7))
      ppname(26) = trim(gname(sNXC3H7))
      ppname(27) = trim(gname(sC3H6OOH))
      ppname(28) = trim(gname(sC3H8))
      ppname(29) = trim(gname(sOC3H5OOH))
      ppname(30) = trim(gname(sC4H8))
      ppname(31) = '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(sO) = 0.016_pr
      W_sp(sOH) = 0.017008_pr
      W_sp(sH) = 0.001008_pr
      W_sp(sH2O) = 0.018016_pr
      W_sp(sH2O2) = 0.034016_pr
      W_sp(sO2) = 0.032_pr
      W_sp(sHO2) = 0.033008_pr
      W_sp(sTXCH2) = 0.014026_pr
      W_sp(sCH2O) = 0.030026_pr
      W_sp(sCO2) = 0.04401_pr
      W_sp(sCO) = 0.02801_pr
      W_sp(sSXCH2) = 0.014026_pr
      W_sp(sCH3) = 0.015034_pr
      W_sp(sC2H6) = 0.030068_pr
      W_sp(sHCO) = 0.029018_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(sC2H5) = 0.02906_pr
      W_sp(sCH2CO) = 0.042036_pr
      W_sp(sC3H5) = 0.04107_pr
      W_sp(sIXC3H7) = 0.043086_pr
      W_sp(sNXC3H7) = 0.043086_pr
      W_sp(sC3H6OOH) = 0.075086_pr
      W_sp(sC3H8) = 0.044094_pr
      W_sp(sOC3H5OOH) = 0.090078_pr
      W_sp(sC4H8) = 0.056104_pr

      end subroutine molar_mass

      ! Species names
      subroutine species_name
      implicit none

      gname(sN2) = 'N2'
      gname(sH2) = 'H2'
      gname(sO) = 'O'
      gname(sOH) = 'OH'
      gname(sH) = 'H'
      gname(sH2O) = 'H2O'
      gname(sH2O2) = 'H2O2'
      gname(sO2) = 'O2'
      gname(sHO2) = 'HO2'
      gname(sTXCH2) = 'T-CH2'
      gname(sCH2O) = 'CH2O'
      gname(sCO2) = 'CO2'
      gname(sCO) = 'CO'
      gname(sSXCH2) = 'S-CH2'
      gname(sCH3) = 'CH3'
      gname(sC2H6) = 'C2H6'
      gname(sHCO) = 'HCO'
      gname(sC2H3) = 'C2H3'
      gname(sCH2CHO) = 'CH2CHO'
      gname(sC3H6) = 'C3H6'
      gname(sC2H4) = 'C2H4'
      gname(sC2H5) = 'C2H5'
      gname(sCH2CO) = 'CH2CO'
      gname(sC3H5) = 'C3H5'
      gname(sIXC3H7) = 'I-C3H7'
      gname(sNXC3H7) = 'N-C3H7'
      gname(sC3H6OOH) = 'C3H6OOH'
      gname(sC3H8) = 'C3H8'
      gname(sOC3H5OOH) = 'OC3H5OOH'
      gname(sC4H8) = 'C4H8'

      end subroutine species_name

      ! Reaction names
      subroutine reaction_name
      implicit none

      rname(r1f) = '1f'
      rname(r2f) = '2f'
      rname(r3f) = '3f'
      rname(r4f) = '4f'
      rname(r5f) = '5f'
      rname(r6f) = '6f'
      rname(r7f) = '7f'
      rname(r8f) = '8f'
      rname(r9f) = '9f'
      rname(r10) = '10'
      rname(r11) = '11'
      rname(r12) = '12'
      rname(r13) = '13'
      rname(r14) = '14'
      rname(r15) = '15'
      rname(r16) = '16'
      rname(r17) = '17'
      rname(r18) = '18'
      rname(r19) = '19'
      rname(r20) = '20'
      rname(r21) = '21'
      rname(r22f) = '22f'
      rname(r23) = '23'
      rname(r24f) = '24f'
      rname(r25) = '25'
      rname(r26) = '26'
      rname(r27) = '27'
      rname(r28f) = '28f'
      rname(r29) = '29'
      rname(r30) = '30'
      rname(r31) = '31'
      rname(r32f) = '32f'
      rname(r33) = '33'
      rname(r34) = '34'
      rname(r35) = '35'
      rname(r36) = '36'
      rname(r37) = '37'
      rname(r38) = '38'
      rname(r39) = '39'
      rname(r40) = '40'
      rname(r41) = '41'
      rname(r42) = '42'
      rname(r43f) = '43f'
      rname(r44) = '44'
      rname(r45f) = '45f'
      rname(r46f) = '46f'
      rname(r47) = '47'
      rname(r48) = '48'
      rname(r49) = '49'
      rname(r50) = '50'
      rname(r51) = '51'
      rname(r52) = '52'
      rname(r53) = '53'
      rname(r54) = '54'
      rname(r55) = '55'
      rname(r56) = '56'
      rname(r57) = '57'
      rname(r58f) = '58f'
      rname(r59) = '59'
      rname(r60) = '60'
      rname(r61f) = '61f'
      rname(r62) = '62'
      rname(r63f) = '63f'
      rname(r64) = '64'
      rname(r65) = '65'
      rname(r66) = '66'
      rname(r67f) = '67f'
      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(r77) = '77'
      rname(r78) = '78'
      rname(r79) = '79'
      rname(r80f) = '80f'
      rname(r81) = '81'
      rname(r82) = '82'
      rname(r83) = '83'
      rname(r84) = '84'
      rname(r85) = '85'
      rname(r86) = '86'
      rname(r87) = '87'
      rname(r88) = '88'
      rname(r89f) = '89f'
      rname(r90) = '90'
      rname(r1b) = '1b'
      rname(r2b) = '2b'
      rname(r3b) = '3b'
      rname(r4b) = '4b'
      rname(r5b) = '5b'
      rname(r6b) = '6b'
      rname(r7b) = '7b'
      rname(r8b) = '8b'
      rname(r9b) = '9b'
      rname(r22b) = '22b'
      rname(r24b) = '24b'
      rname(r28b) = '28b'
      rname(r32b) = '32b'
      rname(r43b) = '43b'
      rname(r45b) = '45b'
      rname(r46b) = '46b'
      rname(r58b) = '58b'
      rname(r61b) = '61b'
      rname(r63b) = '63b'
      rname(r67b) = '67b'
      rname(r72b) = '72b'
      rname(r73b) = '73b'
      rname(r80b) = '80b'
      rname(r89b) = '89b'

      end subroutine reaction_name

      ! List of QSS species
      subroutine QSS_list
      implicit none

      iqss(sN2) = 0
      iqss(sH2) = 0
      iqss(sO) = 0
      iqss(sOH) = 0
      iqss(sH) = 0
      iqss(sH2O) = 0
      iqss(sH2O2) = 0
      iqss(sO2) = 0
      iqss(sHO2) = 0
      iqss(sTXCH2) = 1
      iqss(sCH2O) = 0
      iqss(sCO2) = 0
      iqss(sCO) = 0
      iqss(sSXCH2) = 1
      iqss(sCH3) = 0
      iqss(sC2H6) = 0
      iqss(sHCO) = 1
      iqss(sC2H3) = 1
      iqss(sCH2CHO) = 1
      iqss(sC3H6) = 0
      iqss(sC2H4) = 0
      iqss(sC2H5) = 1
      iqss(sCH2CO) = 0
      iqss(sC3H5) = 1
      iqss(sIXC3H7) = 1
      iqss(sNXC3H7) = 1
      iqss(sC3H6OOH) = 1
      iqss(sC3H8) = 0
      iqss(sOC3H5OOH) = 0
      iqss(sC4H8) = 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(sO) = Y(sO) / W_sp(sO)
      c(sOH) = Y(sOH) / W_sp(sOH)
      c(sH) = Y(sH) / W_sp(sH)
      c(sH2O) = Y(sH2O) / W_sp(sH2O)
      c(sH2O2) = Y(sH2O2) / W_sp(sH2O2)
      c(sO2) = Y(sO2) / W_sp(sO2)
      c(sHO2) = Y(sHO2) / W_sp(sHO2)
      !c(sTXCH2) = Y(sTXCH2) / W_sp(sTXCH2)
      c(sCH2O) = Y(sCH2O) / W_sp(sCH2O)
      c(sCO2) = Y(sCO2) / W_sp(sCO2)
      c(sCO) = Y(sCO) / W_sp(sCO)
      !c(sSXCH2) = Y(sSXCH2) / W_sp(sSXCH2)
      c(sCH3) = Y(sCH3) / W_sp(sCH3)
      c(sC2H6) = Y(sC2H6) / W_sp(sC2H6)
      !c(sHCO) = Y(sHCO) / W_sp(sHCO)
      !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(sC2H5) = Y(sC2H5) / W_sp(sC2H5)
      c(sCH2CO) = Y(sCH2CO) / W_sp(sCH2CO)
      !c(sC3H5) = Y(sC3H5) / W_sp(sC3H5)
      !c(sIXC3H7) = Y(sIXC3H7) / W_sp(sIXC3H7)
      !c(sNXC3H7) = Y(sNXC3H7) / W_sp(sNXC3H7)
      !c(sC3H6OOH) = Y(sC3H6OOH) / W_sp(sC3H6OOH)
      c(sC3H8) = Y(sC3H8) / W_sp(sC3H8)
      c(sOC3H5OOH) = Y(sOC3H5OOH) / W_sp(sOC3H5OOH)
      c(sC4H8) = Y(sC4H8) / W_sp(sC4H8)

      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(15) :: M

      M(mM31) = (5_pr)*c(sH2O) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM7) = (3_pr)*c(sCO2) &
                 + (11_pr)*c(sH2O) &
                 + (1_pr)*c(sCO) &
                 + (1.5_pr)*c(sH2) &
                 + sum(c) 
      M(mM11) = (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + (5_pr)*c(sH2O) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM0) = sum(c) 
      M(mM30) = (5_pr)*c(sH2O) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM8) = (11_pr)*c(sH2O) &
                 + (1.5_pr)*c(sCO) &
                 + (0.9_pr)*c(sH2) &
                 + (1.5_pr)*c(sCO2) &
                 + sum(c) 
      M(mM29) = (1_pr)*c(sCO2) &
                 + (5_pr)*c(sH2O) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + sum(c) 
      M(mM6) = (0.5_pr)*c(sCO) &
                 + (1.5_pr)*c(sH2) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sCO2) &
                 + (5_pr)*c(sH2O2) &
                 + sum(c) 
      M(mM16) = (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + (5_pr)*c(sH2O) &
                 + (1_pr)*c(sCO2) &
                 + sum(c) 
      M(mM15) = (1_pr)*c(sCO2) &
                 + (1_pr)*c(sH2) &
                 + (0.5_pr)*c(sCO) &
                 + (5_pr)*c(sH2O) &
                 + sum(c) 
      M(mM12) = (2.6_pr)*c(sCO2) &
                 + (14.4_pr)*c(sH2O) &
                 + (1.4_pr)*c(sH2) &
                 + (0.8_pr)*c(sCO) &
                 + sum(c) 
      M(mM5) = (1.5_pr)*c(sH2) &
                 + (0.2_pr)*c(sCO) &
                 + (15_pr)*c(sH2O) &
                 + (0.5_pr)*c(sC2H6) &
                 + (1.4_pr)*c(sCO2) &
                 + sum(c) 
      M(mM2) = (2.8_pr)*c(sCO2) &
                 + (11_pr)*c(sH2O) &
                 + (1.5_pr)*c(sH2) &
                 + (0.9_pr)*c(sCO) &
                 + sum(c) 
      M(mM28) = (1_pr)*c(sCO2) &
                 + (5_pr)*c(sH2O) &
                 + (2_pr)*c(sC2H6) &
                 + (1_pr)*c(sH2) &
                 + (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(15) :: M
      real(pr) :: Tloc,Ploc

      real(pr) :: k4f_0, k4f_inf, FC4f
      real(pr) :: k6f_0, k6f_inf, FC6f
      real(pr) :: k26_0, k26_inf, FC26
      real(pr) :: k29_0, k29_inf, FC29
      real(pr) :: k39_0, k39_inf, FC39
      real(pr) :: k46f_0, k46f_inf, FC46f
      real(pr) :: k63f_0, k63f_inf, FC63f
      real(pr) :: k67f_0, k67f_inf, FC67f
      real(pr) :: k73f_0, k73f_inf, FC73f
      real(pr) :: k74_0, k74_inf, FC74
      real(pr) :: k80f_0, k80f_inf, FC80f
      real(pr) :: k4b_0, k4b_inf, FC4b
      real(pr) :: k6b_0, k6b_inf, FC6b
      real(pr) :: k46b_0, k46b_inf, FC46b
      real(pr) :: k63b_0, k63b_inf, FC63b
      real(pr) :: k67b_0, k67b_inf, FC67b
      real(pr) :: k73b_0, k73b_inf, FC73b
      real(pr) :: k80b_0, k80b_inf, FC80b

      ! Rate coefficients

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

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

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

        k4f_0 = (2.76000000e+13_pr)*Tloc**(-3.200_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        k4f_inf = (9.55000000e+07_pr)*Tloc**(-0.270_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC4f = (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(r4f) =  &
          getlindratecoeff & 
        (Tloc, k4f_0, k4f_inf, FC4f, M(mM6), Ploc ) 

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

        k6f_0 = (5.75000000e+07_pr)*Tloc**(-1.400_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        k6f_inf = (4.65000000e+06_pr)*Tloc**(0.440_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC6f = (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(r6f) =  &
          getlindratecoeff & 
        (Tloc, k6f_0, k6f_inf, FC6f, M(mM5), Ploc ) 

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

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

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

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

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

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

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

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

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

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

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

        k(r18) = (2.63000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(6.240e+03_pr)/(8.314_pr*Tloc)) 

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

        k(r20) = (6.58000000e+06_pr)*Tloc**(0.000_pr)*&
          exp(-(6.240e+03_pr)/(8.314_pr*Tloc)) 

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

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

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

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

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

        k26_0 = (1.27000000e+29_pr)*Tloc**(-7.000_pr)*&
          exp(-(1.156e+04_pr)/(8.314_pr*Tloc)) 
 
        k26_inf = (1.81000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC26 = (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(r26) =  &
          getlindratecoeff & 
        (Tloc, k26_0, k26_inf, FC26, M(mM11), Ploc ) 

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

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

        k29_0 = (1.55000000e+12_pr)*Tloc**(-2.790_pr)*&
          exp(-(1.754e+04_pr)/(8.314_pr*Tloc)) 
 
        k29_inf = (1.80000000e+05_pr)*Tloc**(0.000_pr)* &
          exp(-(9.975e+03_pr)/(8.314_pr*Tloc)) 
 
        FC29 = (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(r29) =  &
          getlindratecoeff & 
        (Tloc, k29_0, k29_inf, FC29, M(mM7), Ploc ) 

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

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

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

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

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

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

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

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

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

        k39_0 = (4.27000000e+46_pr)*Tloc**(-11.940_pr)*&
          exp(-(4.088e+04_pr)/(8.314_pr*Tloc)) 
 
        k39_inf = (2.50000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC39 = (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(r39) =  &
          getlindratecoeff & 
        (Tloc, k39_0, k39_inf, FC39, M(mM29), Ploc ) 

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

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

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

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

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

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

        k46f_0 = (3.99000000e+27_pr)*Tloc**(-4.990_pr)*&
          exp(-(1.674e+05_pr)/(8.314_pr*Tloc)) 
 
        k46f_inf = (1.11000000e+10_pr)*Tloc**(1.037_pr)* &
          exp(-(1.538e+05_pr)/(8.314_pr*Tloc)) 
 
        FC46f = (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(r46f) =  &
          getlindratecoeff & 
        (Tloc, k46f_0, k46f_inf, FC46f, M(mM15), Ploc ) 

        k(r47) = (4.83200000e+12_pr)*Tloc**(-1.176_pr)*&
          exp(-(2.466e+04_pr)/(8.314_pr*Tloc)) 

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

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

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

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

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

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

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

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

        k(r56) = (1.00000000e+08_pr)*Tloc**(0.000_pr)*&
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 

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

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

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

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

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

        k(r62) = (5.67600000e+26_pr)*Tloc**(-4.718_pr)*&
          exp(-(1.420e+05_pr)/(8.314_pr*Tloc)) 

        k63f_0 = (1.33000000e+48_pr)*Tloc**(-12.000_pr)*&
          exp(-(2.497e+04_pr)/(8.314_pr*Tloc)) 
 
        k63f_inf = (2.00000000e+08_pr)*Tloc**(0.000_pr)* &
          exp(-(0.000e+00_pr)/(8.314_pr*Tloc)) 
 
        FC63f = (9.800e-01_pr)* &
          exp(-Tloc/(1.097e+03_pr)) + (2.000e-02_pr)* &
          exp(-Tloc/(1.097e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.860e+03_pr)/Tloc) 
 
        k(r63f) =  &
          getlindratecoeff & 
        (Tloc, k63f_0, k63f_inf, FC63f, M(mM28), Ploc ) 

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

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

        k(r66) = (1.70000000e-01_pr)*Tloc**(2.500_pr)*&
          exp(-(1.043e+04_pr)/(8.314_pr*Tloc)) 

        k67f_0 = (8.70000000e+30_pr)*Tloc**(-7.500_pr)*&
          exp(-(1.980e+04_pr)/(8.314_pr*Tloc)) 
 
        k67f_inf = (1.33000000e+07_pr)*Tloc**(0.000_pr)* &
          exp(-(6.530e+03_pr)/(8.314_pr*Tloc)) 
 
        FC67f = (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(r67f) =  &
          getlindratecoeff & 
        (Tloc, k67f_0, k67f_inf, FC67f, M(mM30), Ploc ) 

        k(r68) = (3.10000000e+00_pr)*Tloc**(2.000_pr)*&
          exp(-(-1.248e+03_pr)/(8.314_pr*Tloc)) 

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

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

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

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

        k73f_0 = (5.49000000e+43_pr)*Tloc**(-10.000_pr)*&
          exp(-(1.497e+05_pr)/(8.314_pr*Tloc)) 
 
        k73f_inf = (1.23000000e+13_pr)*Tloc**(-0.100_pr)* &
          exp(-(1.264e+05_pr)/(8.314_pr*Tloc)) 
 
        FC73f = (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(r73f) =  &
          getlindratecoeff & 
        (Tloc, k73f_0, k73f_inf, FC73f, M(mM0), Ploc ) 

        k74_0 = (2.15800000e+33_pr)*Tloc**(-6.949_pr)*&
          exp(-(1.698e+05_pr)/(8.314_pr*Tloc)) 
 
        k74_inf = (4.58600000e+13_pr)*Tloc**(-0.289_pr)* &
          exp(-(1.541e+05_pr)/(8.314_pr*Tloc)) 
 
        FC74 = (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(r74) =  &
          getlindratecoeff & 
        (Tloc, k74_0, k74_inf, FC74, M(mM31), Ploc ) 

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

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

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

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

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

        k80f_0 = (7.83000000e+12_pr)*Tloc**(0.000_pr)*&
          exp(-(2.719e+05_pr)/(8.314_pr*Tloc)) 
 
        k80f_inf = (1.10000000e+17_pr)*Tloc**(0.000_pr)* &
          exp(-(3.531e+05_pr)/(8.314_pr*Tloc)) 
 
        FC80f = (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(r80f) =  &
          getlindratecoeff & 
        (Tloc, k80f_0, k80f_inf, FC80f, M(mM0), Ploc ) 

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

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

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

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

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

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

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

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

        k(r89f) = (1.00000000e+16_pr)*Tloc**(0.000_pr)*&
          exp(-(3.050e+05_pr)/(8.314_pr*Tloc)) 

        k(r90) = (6.60000000e-01_pr)*Tloc**(2.540_pr)*&
          exp(-(2.830e+04_pr)/(8.314_pr*Tloc)) 

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

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

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

        k4b_0 = (1.19344641e+25_pr)*Tloc**(-4.25_pr)*&
          exp(-(2.1486e+05_pr)/(8.314_pr*Tloc)) 
 
        k4b_inf = (4.12949754e+19_pr)*Tloc**(-1.32_pr)* &
          exp(-(2.1486e+05_pr)/(8.314_pr*Tloc)) 
 
        FC4b = (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(r4b) =  &
          getlindratecoeff & 
        (Tloc, k4b_0, k4b_inf, FC4b, M(mM6), Ploc ) 

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

        k6b_0 = (1.02533717e+14_pr)*Tloc**(-1.44_pr)*&
          exp(-(2.0482e+05_pr)/(8.314_pr*Tloc)) 
 
        k6b_inf = (8.29185711e+12_pr)*Tloc**(0.40_pr)* &
          exp(-(2.0482e+05_pr)/(8.314_pr*Tloc)) 
 
        FC6b = (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(r6b) =  &
          getlindratecoeff & 
        (Tloc, k6b_0, k6b_inf, FC6b, M(mM5), Ploc ) 

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

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

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

        k(r22b) = (4.21142619e+06_pr)*Tloc**(-0.06_pr)*&
          exp(-(3.7280e+04_pr)/(8.314_pr*Tloc)) 

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

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

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

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

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

        k46b_0 = (1.72614099e+21_pr)*Tloc**(-4.84_pr)*&
          exp(-(1.5915e+04_pr)/(8.314_pr*Tloc)) 
 
        k46b_inf = (4.80204635e+03_pr)*Tloc**(1.19_pr)* &
          exp(-(2.3148e+03_pr)/(8.314_pr*Tloc)) 
 
        FC46b = (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(r46b) =  &
          getlindratecoeff & 
        (Tloc, k46b_0, k46b_inf, FC46b, M(mM15), Ploc ) 

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

        k(r61b) = (6.03518455e+06_pr)*Tloc**(0.38_pr)*&
          exp(-(1.6341e+05_pr)/(8.314_pr*Tloc)) 

        k63b_0 = (5.38095569e+54_pr)*Tloc**(-11.66_pr)*&
          exp(-(3.9320e+05_pr)/(8.314_pr*Tloc)) 
 
        k63b_inf = (8.09166270e+14_pr)*Tloc**(0.34_pr)* &
          exp(-(3.6823e+05_pr)/(8.314_pr*Tloc)) 
 
        FC63b = (9.800e-01_pr)* &
          exp(-Tloc/(1.097e+03_pr)) + (2.000e-02_pr)* &
          exp(-Tloc/(1.097e+03_pr)) + (1.000e+00_pr)* &
          exp(-(6.860e+03_pr)/Tloc) 
 
        k(r63b) =  &
          getlindratecoeff & 
        (Tloc, k63b_0, k63b_inf, FC63b, M(mM28), Ploc ) 

        k67b_0 = (2.85140703e+36_pr)*Tloc**(-7.39_pr)*&
          exp(-(1.6819e+05_pr)/(8.314_pr*Tloc)) 
 
        k67b_inf = (4.35904753e+12_pr)*Tloc**(0.11_pr)* &
          exp(-(1.5492e+05_pr)/(8.314_pr*Tloc)) 
 
        FC67b = (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(r67b) =  &
          getlindratecoeff & 
        (Tloc, k67b_0, k67b_inf, FC67b, M(mM30), Ploc ) 

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

        k73b_0 = (7.39650348e+30_pr)*Tloc**(-8.33_pr)*&
          exp(-(4.3656e+04_pr)/(8.314_pr*Tloc)) 
 
        k73b_inf = (1.65714012e+00_pr)*Tloc**(1.57_pr)* &
          exp(-(2.0356e+04_pr)/(8.314_pr*Tloc)) 
 
        FC73b = (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(r73b) =  &
          getlindratecoeff & 
        (Tloc, k73b_0, k73b_inf, FC73b, M(mM0), Ploc ) 

        k80b_0 = (1.08695906e-03_pr)*Tloc**(1.79_pr)*&
          exp(-(-1.0465e+05_pr)/(8.314_pr*Tloc)) 
 
        k80b_inf = (1.52701784e+01_pr)*Tloc**(1.79_pr)* &
          exp(-(-2.3452e+04_pr)/(8.314_pr*Tloc)) 
 
        FC80b = (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(r80b) =  &
          getlindratecoeff & 
        (Tloc, k80b_0, k80b_inf, FC80b, M(mM0), Ploc ) 

        k(r89b) = (6.66353595e+02_pr)*Tloc**(1.32_pr)*&
          exp(-(-1.9116e+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(15) :: m

      w(r1f) = k(r1f)  * c(sH2) * c(sO)
      w(r2f) = k(r2f)  * c(sH2) * c(sOH)
      w(r3f) = k(r3f)  * c(sH) * c(sOH) * m(mM2)
      w(r4f) = k(r4f)  * c(sOH)**2_pr 
      w(r5f) = k(r5f)  * c(sH2O) * c(sO)
      w(r6f) = k(r6f)  * c(sH) * c(sO2)
      w(r7f) = k(r7f)  * c(sH) * c(sO2)
      w(r8f) = k(r8f)  * c(sHO2) * c(sOH)
      w(r9f) = k(r9f)  * c(sHO2) * c(sOH)
      w(r10) = k(r10)  * c(sHO2) * c(sH)
      w(r11) = k(r11)  * c(sHO2) * c(sH)
      w(r12) = k(r12)  * c(sHO2)**2_pr 
      w(r13) = k(r13)  * c(sHO2) * c(sO)
      w(r14) = k(r14)  * c(sHO2)**2_pr 
      w(r15) = k(r15)  * c(sH2O2) * c(sOH)
      w(r16) = k(r16)  * c(sH2O2) * c(sOH)
      w(r17) = k(r17)  * cqss(sTXCH2-nspec) * c(sOH)
      w(r18) = k(r18)  * cqss(sTXCH2-nspec) * c(sO2)
      w(r19) = k(r19)  * cqss(sTXCH2-nspec) * c(sO)
      w(r20) = k(r20)  * cqss(sTXCH2-nspec) * c(sO2)
      w(r21) = k(r21)  * cqss(sTXCH2-nspec) * c(sO)
      w(r22f) = k(r22f)  * cqss(sSXCH2-nspec) * m(mM12)
      w(r23) = k(r23)  * cqss(sSXCH2-nspec) * c(sO2)
      w(r24f) = k(r24f)  * c(sCH3) * c(sOH)
      w(r25) = k(r25)  * c(sCH3) * c(sO2)
      w(r26) = k(r26)  * c(sCH3)**2_pr 
      w(r27) = k(r27)  * c(sCH3) * c(sO)
      w(r28f) = k(r28f)  * c(sCO) * c(sOH)
      w(r29) = k(r29)  * c(sCO) * c(sO)
      w(r30) = k(r30)  * c(sCO) * c(sHO2)
      w(r31) = k(r31)  * cqss(sHCO-nspec) * c(sOH)
      w(r32f) = k(r32f)  * cqss(sHCO-nspec) * m(mM8)
      w(r33) = k(r33)  * cqss(sHCO-nspec) * c(sO2)
      w(r34) = k(r34)  * c(sCH2O) * c(sHO2)
      w(r35) = k(r35)  * c(sCH2O) * c(sOH)
      w(r36) = k(r36)  * c(sCH2O) * c(sH)
      w(r37) = k(r37)  * c(sCH2O) * c(sO)
      w(r38) = k(r38)  * cqss(sC2H3-nspec) * c(sO2)
      w(r39) = k(r39)  * cqss(sC2H3-nspec) * c(sCH3)
      w(r40) = k(r40)  * cqss(sC2H3-nspec) * c(sO2)
      w(r41) = k(r41)  * c(sH) * cqss(sC2H3-nspec) * m(mM16)
      w(r42) = k(r42)  * c(sC2H4) * c(sO)
      w(r43f) = k(r43f)  * c(sC2H4) * c(sH)
      w(r44) = k(r44)  * c(sC2H4) * c(sO)
      w(r45f) = k(r45f)  * c(sC2H4) * c(sOH)
      w(r46f) = k(r46f)  * cqss(sC2H5-nspec)
      w(r47) = k(r47)  * c(sH) * cqss(sC2H5-nspec)
      w(r48) = k(r48)  * cqss(sC2H5-nspec) * c(sO2)
      w(r49) = k(r49)  * c(sC2H6) * c(sO)
      w(r50) = k(r50)  * c(sC2H6) * c(sH)
      w(r51) = k(r51)  * c(sC2H6) * c(sOH)
      w(r52) = k(r52)  * c(sCH2CO) * c(sO)
      w(r53) = k(r53)  * c(sCH2CO) * c(sCH3)
      w(r54) = k(r54)  * c(sCH2CO) * c(sH)
      w(r55) = k(r55)  * cqss(sCH2CHO-nspec) * c(sO2)
      w(r56) = k(r56)  * cqss(sCH2CHO-nspec) * c(sO)
      w(r57) = k(r57)  * cqss(sCH2CHO-nspec) * c(sCH3)
      w(r58f) = k(r58f)  * cqss(sCH2CHO-nspec)
      w(r59) = k(r59)  * cqss(sCH2CHO-nspec) * c(sHO2)
      w(r60) = k(r60)  * cqss(sCH2CHO-nspec) * c(sOH)
      w(r61f) = k(r61f)  * cqss(sC3H5-nspec) * c(sHO2)
      w(r62) = k(r62)  * c(sH) * cqss(sC3H5-nspec)
      w(r63f) = k(r63f)  * cqss(sC3H5-nspec) * c(sH)
      w(r64) = k(r64)  * cqss(sC3H5-nspec) * c(sHO2)
      w(r65) = k(r65)  * c(sC3H6) * c(sH)
      w(r66) = k(r66)  * c(sC3H6) * c(sH)
      w(r67f) = k(r67f)  * c(sC3H6) * c(sH)
      w(r68) = k(r68)  * c(sC3H6) * c(sOH)
      w(r69) = k(r69)  * c(sC3H6) * c(sO)
      w(r70) = k(r70)  * c(sC3H6) * c(sO)
      w(r71) = k(r71)  * cqss(sIXC3H7-nspec) * c(sO2)
      w(r72f) = k(r72f)  * cqss(sNXC3H7-nspec) * c(sO2)
      w(r73f) = k(r73f)  * cqss(sNXC3H7-nspec)
      w(r74) = k(r74)  * cqss(sNXC3H7-nspec)
      w(r75) = k(r75)  * cqss(sNXC3H7-nspec) * c(sO2)
      w(r76) = k(r76)  * c(sC3H8) * c(sH)
      w(r77) = k(r77)  * c(sC3H8) * c(sO)
      w(r78) = k(r78)  * c(sC3H8) * c(sOH)
      w(r79) = k(r79)  * c(sC3H8) * c(sO)
      w(r80f) = k(r80f)  * c(sC3H8)
      w(r81) = k(r81)  * c(sC3H8) * c(sO2)
      w(r82) = k(r82)  * c(sC3H8) * c(sHO2)
      w(r83) = k(r83)  * c(sC3H8) * c(sH)
      w(r84) = k(r84)  * c(sC3H8) * c(sO2)
      w(r85) = k(r85)  * c(sC3H8) * c(sHO2)
      w(r86) = k(r86)  * cqss(sC3H6OOH-nspec) * c(sO2)
      w(r87) = k(r87)  * cqss(sC3H6OOH-nspec)
      w(r88) = k(r88)  * c(sOC3H5OOH)
      w(r89f) = k(r89f)  * c(sC4H8)
      w(r90) = k(r90)  * c(sC4H8) * c(sH)
      w(r1b) = k(r1b)  * c(sOH) * c(sH)
      w(r2b) = k(r2b)  * c(sH2O) * c(sH)
      w(r3b) = k(r3b)  * c(sH2O) * m(mM2)
      w(r4b) = k(r4b)  * c(sH2O2)
      w(r5b) = k(r5b)  * c(sOH)**2_pr 
      w(r6b) = k(r6b)  * c(sHO2)
      w(r7b) = k(r7b)  * c(sOH) * c(sO)
      w(r8b) = k(r8b)  * c(sH2O) * c(sO2)
      w(r9b) = k(r9b)  * c(sH2O) * c(sO2)
      w(r22b) = k(r22b)  * cqss(sTXCH2-nspec) * m(mM12)
      w(r24b) = k(r24b)  * cqss(sSXCH2-nspec) * c(sH2O)
      w(r28b) = k(r28b)  * c(sCO2) * c(sH)
      w(r32b) = k(r32b)  * c(sCO) * c(sH) * m(mM8)
      w(r43b) = k(r43b)  * cqss(sC2H3-nspec) * c(sH2)
      w(r45b) = k(r45b)  * cqss(sC2H3-nspec) * c(sH2O)
      w(r46b) = k(r46b)  * c(sC2H4) * c(sH)
      w(r58b) = k(r58b)  * c(sCH2CO) * c(sH)
      w(r61b) = k(r61b)  * c(sC3H6) * c(sO2)
      w(r63b) = k(r63b)  * c(sC3H6)
      w(r67b) = k(r67b)  * cqss(sIXC3H7-nspec)
      w(r72b) = k(r72b)  * cqss(sC3H6OOH-nspec)
      w(r73b) = k(r73b)  * c(sCH3) * c(sC2H4)
      w(r80b) = k(r80b)  * c(sCH3) * cqss(sC2H5-nspec)
      w(r89b) = k(r89b)  * cqss(sC3H5-nspec) * c(sCH3)

      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(r1f) - &
         w(r2f) + &
         w(r10) + &
         w(r18) + &
         w(r19) + &
         w(r36) + &
         w(r43f) + &
         w(r50) + &
         w(r66) + &
         w(r76) + &
         w(r83) + &
         w(r90) + &
         w(r1b) + &
         w(r2b) - &
         w(r43b) 


        cdot(sO) = 0.0_pr - &
         w(r1f) - &
         w(r5f) + &
         w(r7f) - &
         w(r13) - &
         w(r19) - &
         w(r21) - &
         w(r27) - &
         w(r29) - &
         w(r37) + &
         w(r38) - &
         w(r42) - &
         w(r44) - &
         w(r49) - &
         w(r52) - &
         w(r56) - &
         w(r69) - &
         w(r70) - &
         w(r77) - &
         w(r79) + &
         w(r1b) + &
         w(r5b) - &
         w(r7b) 


        cdot(sOH) = 0.0_pr + &
         w(r1f) - &
         w(r2f) - &
         w(r3f) - 2_pr * &
         w(r4f) + 2_pr * &
         w(r5f) + &
         w(r7f) - &
         w(r8f) - &
         w(r9f) + 2_pr * &
         w(r11) + &
         w(r13) - &
         w(r15) - &
         w(r16) - &
         w(r17) + &
         w(r20) + &
         w(r23) - &
         w(r24f) + &
         w(r25) - &
         w(r28f) + &
         w(r30) - &
         w(r31) - &
         w(r35) + &
         w(r37) - &
         w(r45f) + &
         w(r49) - &
         w(r51) + &
         w(r55) + &
         w(r59) - &
         w(r60) + &
         w(r64) - &
         w(r68) + &
         w(r77) - &
         w(r78) + &
         w(r79) + &
         w(r86) + &
         w(r88) - &
         w(r1b) + &
         w(r2b) + &
         w(r3b) + 2_pr * &
         w(r4b) - 2_pr * &
         w(r5b) - &
         w(r7b) + &
         w(r8b) + &
         w(r9b) + &
         w(r24b) + &
         w(r28b) + &
         w(r45b) 


        cdot(sH) = 0.0_pr + &
         w(r1f) + &
         w(r2f) - &
         w(r3f) - &
         w(r6f) - &
         w(r7f) - &
         w(r10) - &
         w(r11) + &
         w(r17) + &
         w(r20) + 2_pr * &
         w(r21) + &
         w(r23) + &
         w(r27) + &
         w(r28f) + &
         w(r32f) - &
         w(r36) - &
         w(r41) + &
         w(r42) - &
         w(r43f) + &
         w(r46f) - &
         w(r47) - &
         w(r50) - &
         w(r54) + &
         w(r57) + &
         w(r58f) - &
         w(r62) - &
         w(r63f) - &
         w(r65) - &
         w(r66) - &
         w(r67f) + &
         w(r69) + &
         w(r74) - &
         w(r76) - &
         w(r83) - &
         w(r90) - &
         w(r1b) - &
         w(r2b) + &
         w(r3b) + &
         w(r6b) + &
         w(r7b) - &
         w(r28b) - &
         w(r32b) + &
         w(r43b) - &
         w(r46b) - &
         w(r58b) + &
         w(r63b) + &
         w(r67b) 


        cdot(sH2O) = 0.0_pr + &
         w(r2f) + &
         w(r3f) - &
         w(r5f) + &
         w(r8f) + &
         w(r9f) + &
         w(r15) + &
         w(r16) + &
         w(r24f) + &
         w(r31) + &
         w(r35) + &
         w(r45f) + &
         w(r51) + &
         w(r60) + &
         w(r68) + &
         w(r78) - &
         w(r2b) - &
         w(r3b) + &
         w(r5b) - &
         w(r8b) - &
         w(r9b) - &
         w(r24b) - &
         w(r45b) 


        cdot(sH2O2) = 0.0_pr + &
         w(r4f) + &
         w(r12) + &
         w(r14) - &
         w(r15) - &
         w(r16) + &
         w(r34) + &
         w(r82) + &
         w(r85) - &
         w(r4b) 


        cdot(sO2) = 0.0_pr - &
         w(r6f) - &
         w(r7f) + &
         w(r8f) + &
         w(r9f) + &
         w(r10) + &
         w(r12) + &
         w(r13) + &
         w(r14) - &
         w(r18) - &
         w(r20) - &
         w(r23) - &
         w(r25) - &
         w(r33) - &
         w(r38) - &
         w(r40) - &
         w(r48) - &
         w(r55) + &
         w(r61f) - &
         w(r71) - &
         w(r72f) - &
         w(r75) - &
         w(r81) - &
         w(r84) - &
         w(r86) + &
         w(r6b) + &
         w(r7b) - &
         w(r8b) - &
         w(r9b) - &
         w(r61b) + &
         w(r72b) 


        cdot(sHO2) = 0.0_pr + &
         w(r6f) - &
         w(r8f) - &
         w(r9f) - &
         w(r10) - &
         w(r11) - 2_pr * &
         w(r12) - &
         w(r13) - 2_pr * &
         w(r14) + &
         w(r15) + &
         w(r16) - &
         w(r30) + &
         w(r33) - &
         w(r34) + &
         w(r48) - &
         w(r59) - &
         w(r61f) - &
         w(r64) + &
         w(r71) + &
         w(r75) + &
         w(r81) - &
         w(r82) + &
         w(r84) - &
         w(r85) + &
         w(r87) - &
         w(r6b) + &
         w(r8b) + &
         w(r9b) + &
         w(r61b) 

      !cdot(sTXCH2) = 0.0_pr


        cdot(sCH2O) = 0.0_pr + &
         w(r17) + &
         w(r25) + &
         w(r27) - &
         w(r34) - &
         w(r35) - &
         w(r36) - &
         w(r37) + &
         w(r40) + &
         w(r55) + &
         w(r56) + &
         w(r59) + &
         w(r64) + &
         w(r88) 


        cdot(sCO2) = 0.0_pr + &
         w(r18) + &
         w(r28f) + &
         w(r29) + &
         w(r30) + &
         w(r52) - &
         w(r28b) 


        cdot(sCO) = 0.0_pr + &
         w(r19) + &
         w(r20) + &
         w(r21) + &
         w(r23) - &
         w(r28f) - &
         w(r29) - &
         w(r30) + &
         w(r31) + &
         w(r32f) + &
         w(r33) + &
         w(r53) + &
         w(r54) + &
         w(r55) + &
         w(r57) + &
         w(r28b) - &
         w(r32b) 

      !cdot(sSXCH2) = 0.0_pr


        cdot(sCH3) = 0.0_pr - &
         w(r24f) - &
         w(r25) - 2_pr * &
         w(r26) - &
         w(r27) - &
         w(r39) + &
         w(r44) + 2_pr * &
         w(r47) - &
         w(r53) + &
         w(r54) - &
         w(r57) + &
         w(r62) + &
         w(r65) + &
         w(r69) + &
         w(r73f) + &
         w(r80f) + &
         w(r89f) + &
         w(r24b) - &
         w(r73b) - &
         w(r80b) - &
         w(r89b) 


        cdot(sC2H6) = 0.0_pr + &
         w(r26) - &
         w(r49) - &
         w(r50) - &
         w(r51) 

      !cdot(sHCO) = 0.0_pr

      !cdot(sC2H3) = 0.0_pr

      !cdot(sCH2CHO) = 0.0_pr


        cdot(sC3H6) = 0.0_pr + &
         w(r39) + &
         w(r61f) + &
         w(r63f) - &
         w(r65) - &
         w(r66) - &
         w(r67f) - &
         w(r68) - &
         w(r69) - &
         w(r70) + &
         w(r71) + &
         w(r74) + &
         w(r75) + &
         w(r87) - &
         w(r61b) - &
         w(r63b) + &
         w(r67b) 


        cdot(sC2H4) = 0.0_pr + &
         w(r41) - &
         w(r42) - &
         w(r43f) - &
         w(r44) - &
         w(r45f) + &
         w(r46f) + &
         w(r48) + &
         w(r65) + &
         w(r73f) + &
         w(r90) + &
         w(r43b) + &
         w(r45b) - &
         w(r46b) - &
         w(r73b) 

      !cdot(sC2H5) = 0.0_pr


        cdot(sCH2CO) = 0.0_pr - &
         w(r52) - &
         w(r53) - &
         w(r54) + &
         w(r58f) + &
         w(r60) + &
         w(r69) - &
         w(r58b) 

      !cdot(sC3H5) = 0.0_pr

      !cdot(sIXC3H7) = 0.0_pr

      !cdot(sNXC3H7) = 0.0_pr

      !cdot(sC3H6OOH) = 0.0_pr


        cdot(sC3H8) = 0.0_pr - &
         w(r76) - &
         w(r77) - &
         w(r78) - &
         w(r79) - &
         w(r80f) - &
         w(r81) - &
         w(r82) - &
         w(r83) - &
         w(r84) - &
         w(r85) + &
         w(r80b) 


        cdot(sOC3H5OOH) = 0.0_pr + &
         w(r86) - &
         w(r88) 


        cdot(sC4H8) = 0.0_pr - &
         w(r89f) - &
         w(r90) + &
         w(r89b) 


      end subroutine get_production_rates

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

      reacexp(1) = 'H2 + O -> OH + H'
      reacexp(2) = 'H2 + OH -> H2O + H'
      reacexp(3) = 'H + OH + M2 -> H2O + M2'
      reacexp(4) = '2 OH + M6 -> H2O2 + M6'
      reacexp(5) = 'H2O + O -> 2 OH'
      reacexp(6) = 'H + O2 + M5 -> HO2 + M5'
      reacexp(7) = 'H + O2 -> OH + O'
      reacexp(8) = 'HO2 + OH -> H2O + O2'
      reacexp(9) = 'HO2 + OH -> H2O + O2'
      reacexp(10) = 'HO2 + H -> H2 + O2'
      reacexp(11) = 'HO2 + H -> 2 OH'
      reacexp(12) = '2 HO2 -> H2O2 + O2'
      reacexp(13) = 'HO2 + O -> OH + O2'
      reacexp(14) = '2 HO2 -> H2O2 + O2'
      reacexp(15) = 'H2O2 + OH -> H2O + HO2'
      reacexp(16) = 'H2O2 + OH -> H2O + HO2'
      reacexp(17) = 'T-CH2 + OH -> CH2O + H'
      reacexp(18) = 'T-CH2 + O2 -> CO2 + H2'
      reacexp(19) = 'T-CH2 + O -> CO + H2'
      reacexp(20) = 'T-CH2 + O2 -> CO + OH + H'
      reacexp(21) = 'T-CH2 + O -> CO + 2 H'
      reacexp(22) = 'S-CH2 + M12 -> T-CH2 + M12'
      reacexp(23) = 'S-CH2 + O2 -> CO + OH + H'
      reacexp(24) = 'CH3 + OH -> S-CH2 + H2O'
      reacexp(25) = 'CH3 + O2 -> CH2O + OH'
      reacexp(26) = '2 CH3 + M11 -> C2H6 + M11'
      reacexp(27) = 'CH3 + O -> CH2O + H'
      reacexp(28) = 'CO + OH -> CO2 + H'
      reacexp(29) = 'CO + O + M7 -> CO2 + M7'
      reacexp(30) = 'CO + HO2 -> CO2 + OH'
      reacexp(31) = 'HCO + OH -> CO + H2O'
      reacexp(32) = 'HCO + M8 -> CO + H + M8'
      reacexp(33) = 'HCO + O2 -> CO + HO2'
      reacexp(34) = 'CH2O + HO2 -> HCO + H2O2'
      reacexp(35) = 'CH2O + OH -> HCO + H2O'
      reacexp(36) = 'CH2O + H -> HCO + H2'
      reacexp(37) = 'CH2O + O -> HCO + OH'
      reacexp(38) = 'C2H3 + O2 -> CH2CHO + O'
      reacexp(39) = 'C2H3 + CH3 + M29 -> C3H6 + M29'
      reacexp(40) = 'C2H3 + O2 -> CH2O + HCO'
      reacexp(41) = 'H + C2H3 + M16 -> C2H4 + M16'
      reacexp(42) = 'C2H4 + O -> CH2CHO + H'
      reacexp(43) = 'C2H4 + H -> C2H3 + H2'
      reacexp(44) = 'C2H4 + O -> CH3 + HCO'
      reacexp(45) = 'C2H4 + OH -> C2H3 + H2O'
      reacexp(46) = 'C2H5 + M15 -> C2H4 + H + M15'
      reacexp(47) = 'H + C2H5 -> 2 CH3'
      reacexp(48) = 'C2H5 + O2 -> C2H4 + HO2'
      reacexp(49) = 'C2H6 + O -> C2H5 + OH'
      reacexp(50) = 'C2H6 + H -> C2H5 + H2'
      reacexp(51) = 'C2H6 + OH -> C2H5 + H2O'
      reacexp(52) = 'CH2CO + O -> T-CH2 + CO2'
      reacexp(53) = 'CH2CO + CH3 -> C2H5 + CO'
      reacexp(54) = 'CH2CO + H -> CH3 + CO'
      reacexp(55) = 'CH2CHO + O2 -> CH2O + CO + OH'
      reacexp(56) = 'CH2CHO + O -> CH2O + HCO'
      reacexp(57) = 'CH2CHO + CH3 -> C2H5 + CO + H'
      reacexp(58) = 'CH2CHO -> CH2CO + H'
      reacexp(59) = 'CH2CHO + HO2 -> CH2O + HCO + OH'
      reacexp(60) = 'CH2CHO + OH -> CH2CO + H2O'
      reacexp(61) = 'C3H5 + HO2 -> C3H6 + O2'
      reacexp(62) = 'H + C3H5 -> C2H3 + CH3'
      reacexp(63) = 'C3H5 + H + M28 -> C3H6 + M28'
      reacexp(64) = 'C3H5 + HO2 -> OH + C2H3 + CH2O'
      reacexp(65) = 'C3H6 + H -> C2H4 + CH3'
      reacexp(66) = 'C3H6 + H -> C3H5 + H2'
      reacexp(67) = 'C3H6 + H + M30 -> I-C3H7 + M30'
      reacexp(68) = 'C3H6 + OH -> C3H5 + H2O'
      reacexp(69) = 'C3H6 + O -> CH2CO + CH3 + H'
      reacexp(70) = 'C3H6 + O -> C2H5 + HCO'
      reacexp(71) = 'I-C3H7 + O2 -> C3H6 + HO2'
      reacexp(72) = 'N-C3H7 + O2 -> C3H6OOH'
      reacexp(73) = 'N-C3H7 + M0 -> CH3 + C2H4 + M0'
      reacexp(74) = 'N-C3H7 + M31 -> C3H6 + H + M31'
      reacexp(75) = 'N-C3H7 + O2 -> C3H6 + HO2'
      reacexp(76) = 'C3H8 + H -> N-C3H7 + H2'
      reacexp(77) = 'C3H8 + O -> N-C3H7 + OH'
      reacexp(78) = 'C3H8 + OH -> N-C3H7 + H2O'
      reacexp(79) = 'C3H8 + O -> I-C3H7 + OH'
      reacexp(80) = 'C3H8 + M0 -> CH3 + C2H5 + M0'
      reacexp(81) = 'C3H8 + O2 -> N-C3H7 + HO2'
      reacexp(82) = 'C3H8 + HO2 -> I-C3H7 + H2O2'
      reacexp(83) = 'C3H8 + H -> I-C3H7 + H2'
      reacexp(84) = 'C3H8 + O2 -> I-C3H7 + HO2'
      reacexp(85) = 'C3H8 + HO2 -> N-C3H7 + H2O2'
      reacexp(86) = 'C3H6OOH + O2 -> OC3H5OOH + OH'
      reacexp(87) = 'C3H6OOH -> C3H6 + HO2'
      reacexp(88) = 'OC3H5OOH -> CH2CHO + CH2O + OH'
      reacexp(89) = 'C4H8 -> C3H5 + CH3'
      reacexp(90) = 'C4H8 + H -> H2 + C2H3 + C2H4'
      reacexp(91) = 'Reverse of H2 + O -> OH + H'
      reacexp(92) = 'Reverse of H2 + OH -> H2O + H'
      reacexp(93) = 'Reverse of H + OH + M2 -> H2O + M2'
      reacexp(94) = 'Reverse of 2 OH + M6 -> H2O2 + M6'
      reacexp(95) = 'Reverse of H2O + O -> 2 OH'
      reacexp(96) = 'Reverse of H + O2 + M5 -> HO2 + M5'
      reacexp(97) = 'Reverse of H + O2 -> OH + O'
      reacexp(98) = 'Reverse of HO2 + OH -> H2O + O2'
      reacexp(99) = 'Reverse of HO2 + OH -> H2O + O2'
      reacexp(100) = 'Reverse of S-CH2 + M12 -> T-CH2 + M12'
      reacexp(101) = 'Reverse of CH3 + OH -> S-CH2 + H2O'
      reacexp(102) = 'Reverse of CO + OH -> CO2 + H'
      reacexp(103) = 'Reverse of HCO + M8 -> CO + H + M8'
      reacexp(104) = 'Reverse of C2H4 + H -> C2H3 + H2'
      reacexp(105) = 'Reverse of C2H4 + OH -> C2H3 + H2O'
      reacexp(106) = 'Reverse of C2H5 + M15 -> C2H4 + H + M15'
      reacexp(107) = 'Reverse of CH2CHO -> CH2CO + H'
      reacexp(108) = 'Reverse of C3H5 + HO2 -> C3H6 + O2'
      reacexp(109) = 'Reverse of C3H5 + H + M28 -> C3H6 + M28'
      reacexp(110) = 'Reverse of C3H6 + H + M30 -> I-C3H7 + M30'
      reacexp(111) = 'Reverse of N-C3H7 + O2 -> C3H6OOH'
      reacexp(112) = 'Reverse of N-C3H7 + M0 -> CH3 + C2H4 + M0'
      reacexp(113) = 'Reverse of C3H8 + M0 -> CH3 + C2H5 + M0'
      reacexp(114) = 'Reverse of C4H8 -> C3H5 + CH3'

      end subroutine reaction_expressions

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

      fofb = 0

! Attach corresponding forward reaction to each backward reaction
      fofb(91) = 1
      fofb(92) = 2
      fofb(93) = 3
      fofb(94) = 4
      fofb(95) = 5
      fofb(96) = 6
      fofb(97) = 7
      fofb(98) = 8
      fofb(99) = 9
      fofb(100) = 22
      fofb(101) = 24
      fofb(102) = 28
      fofb(103) = 32
      fofb(104) = 43
      fofb(105) = 45
      fofb(106) = 46
      fofb(107) = 58
      fofb(108) = 61
      fofb(109) = 63
      fofb(110) = 67
      fofb(111) = 72
      fofb(112) = 73
      fofb(113) = 80
      fofb(114) = 89

      end subroutine reverse_reactions

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

      real(pr), dimension(nqss) :: cqss
      real(pr), dimension(nspec) :: c
      real(pr), dimension(nreac) :: k
      real(pr), dimension(15) :: M
      real(pr) :: rho

      real(pr), parameter :: W_sTXCH2 = 0.014026_pr
      real(pr), parameter :: W_sSXCH2 = 0.014026_pr
      real(pr), parameter :: W_sHCO = 0.029018_pr
      real(pr), parameter :: W_sC2H3 = 0.027044_pr
      real(pr), parameter :: W_sCH2CHO = 0.043044_pr
      real(pr), parameter :: W_sC2H5 = 0.02906_pr
      real(pr), parameter :: W_sC3H5 = 0.04107_pr
      real(pr), parameter :: W_sIXC3H7 = 0.043086_pr
      real(pr), parameter :: W_sNXC3H7 = 0.043086_pr
      real(pr), parameter :: W_sC3H6OOH = 0.075086_pr

      real(pr) :: NXC3H7_denom1&
         , NXC3H7_ct1&
         , NXC3H7_denom2&
         , NXC3H7_ct2
      real(pr) :: TXCH2_denom1&
         , TXCH2_ct1&
         , TXCH2_denom2&
         , TXCH2_ct2&
         , TXCH2_SXCH2&
         , TXCH2_SXCH2_coeff
      real(pr) :: SXCH2_denom1&
         , SXCH2_ct1&
         , SXCH2_denom2&
         , SXCH2_ct2
      real(pr) :: C3H6OOH_denom1&
         , C3H6OOH_ct1&
         , C3H6OOH_denom2&
         , C3H6OOH_ct2&
         , C3H6OOH_NXC3H7&
         , C3H6OOH_NXC3H7_coeff


  ! c(sTXCH2) c(sSXCH2) (coupled)  --------------------

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

	 TXCH2_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r17) *  c(sOH) & 
               + k(r18) *  c(sO2) & 
               + k(r19) *  c(sO) & 
               + k(r20) *  c(sO2) & 
               + k(r21) *  c(sO) & 
               + k(r22b) * M(mM12)  & 
                )

	 SXCH2_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r22f) * M(mM12)  & 
               + k(r23) *  c(sO2) & 
               + k(r24b) *  c(sH2O) & 
                )


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

	 TXCH2_ct1 = ( 0.0_pr &
                + k(r52) * c(sCH2CO) * c(sO)                )

	 SXCH2_ct1 = ( 0.0_pr &
                + k(r24f) * c(sCH3) * c(sOH)                )


  ! TXCH2 ---------------------------------------

	 TXCH2_denom2 = tiny(1.0_pr) + ( TXCH2_denom1                 )

	 TXCH2_ct2 = ( TXCH2_ct1 &
                ) / TXCH2_denom2  

	 TXCH2_SXCH2 = ( 0.0_pr &
                + k(r22f)  * M(mM12) &
                ) / TXCH2_denom2 



   	 TXCH2_SXCH2_coeff = ( 0.0_pr &
                + k(r22b)  * M(mM12) &
                )


  ! SXCH2 ---------------------------------------

	 SXCH2_denom2 = tiny(1.0_pr) + ( SXCH2_denom1  &
                - TXCH2_SXCH2_coeff * TXCH2_SXCH2                )

	 SXCH2_ct2 = ( SXCH2_ct1 &
                + TXCH2_SXCH2_coeff * TXCH2_ct2 &
                ) / SXCH2_denom2  




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

	 cqss(sSXCH2-nspec) = ( SXCH2_ct2                 ) 
   !cqss(sSXCH2-nspec) = min( cqss(sSXCH2-nspec), rho/W_sSXCH2 )

	 cqss(sTXCH2-nspec) = ( TXCH2_ct2 &
                + TXCH2_SXCH2 * cqss(sSXCH2-nspec)                ) 

   cqss(sSXCH2-nspec) = min( cqss(sSXCH2-nspec), rho/W_sSXCH2 )
   cqss(sTXCH2-nspec) = min( cqss(sTXCH2-nspec), rho/W_sTXCH2 )

  ! cqss(sC3H5) (uncoupled) --------------------

	 cqss(sC3H5-nspec) = ( 0.0_pr &
                + k(r61b) * c(sO2) * c(sC3H6) &
                + k(r63b) * c(sC3H6) &
                + k(r66) * c(sC3H6) * c(sH) &
                + k(r68) * c(sC3H6) * c(sOH) &
                + k(r89f) * c(sC4H8) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r61f) *  c(sHO2) &
               + k(r62) *  c(sH) &
               + k(r63f) *  c(sH) &
               + k(r64) *  c(sHO2) &
               + k(r89b) *  c(sCH3)                 ) )

   cqss(sC3H5-nspec) = min( cqss(sC3H5-nspec), rho/W_sC3H5 )


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

	 cqss(sC2H3-nspec) = ( 0.0_pr &
                + k(r43f) * c(sC2H4) * c(sH) &
                + k(r45f) * c(sC2H4) * c(sOH) &
                + k(r62) * c(sH) * cqss(sC3H5-nspec) &
                + k(r64) * cqss(sC3H5-nspec) *c(sHO2) &
                + k(r90) * c(sC4H8) * c(sH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r38) *  c(sO2) &
               + k(r39) *  c(sCH3) &
               + k(r40) *  c(sO2) &
               + k(r41) *  c(sH) * M(mM16) &
               + k(r43b) *  c(sH2) &
               + k(r45b) *  c(sH2O)                 ) )
   cqss(sC2H3-nspec) = min( cqss(sC2H3-nspec), rho/W_sC2H3 )


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

	 cqss(sCH2CHO-nspec) = ( 0.0_pr &
                + k(r38) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r42) * c(sC2H4) * c(sO) &
                + k(r58b) * c(sH) * c(sCH2CO) &
                + k(r88) * c(sOC3H5OOH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r55) *  c(sO2) &
               + k(r56) *  c(sO) &
               + k(r57) *  c(sCH3) &
               + k(r58f)  &
               + k(r59) *  c(sHO2) &
               + k(r60) *  c(sOH)                 ) )

   cqss(sCH2CHO-nspec) = min( cqss(sCH2CHO-nspec), rho/W_sCH2CHO )


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

	 cqss(sC2H5-nspec) = ( 0.0_pr &
                + k(r46b) * c(sH) * c(sC2H4) &
                + k(r49) * c(sC2H6) * c(sO) &
                + k(r50) * c(sC2H6) * c(sH) &
                + k(r51) * c(sC2H6) * c(sOH) &
                + k(r53) * c(sCH2CO) * c(sCH3) &
                + k(r57) * cqss(sCH2CHO-nspec) *c(sCH3) &
                + k(r70) * c(sC3H6) * c(sO) &
                + k(r80f) * c(sC3H8) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r46f)  &
               + k(r47) *  c(sH) &
               + k(r48) *  c(sO2) &
               + k(r80b) *  c(sCH3)                 ) )

   cqss(sC2H5-nspec) = min( cqss(sC2H5-nspec), rho/W_sC2H5 )


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

	 cqss(sHCO-nspec) = ( 0.0_pr &
                + k(r32b) * c(sH) * c(sCO) * M(mM8) &
                + k(r34) * c(sCH2O) * c(sHO2) &
                + k(r35) * c(sCH2O) * c(sOH) &
                + k(r36) * c(sCH2O) * c(sH) &
                + k(r37) * c(sCH2O) * c(sO) &
                + k(r40) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r44) * c(sC2H4) * c(sO) &
                + k(r56) * cqss(sCH2CHO-nspec) *c(sO) &
                + k(r59) * cqss(sCH2CHO-nspec) *c(sHO2) &
                + k(r70) * c(sC3H6) * c(sO) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r31) *  c(sOH) &
               + k(r32f) * M(mM8)  &
               + k(r33) *  c(sO2)                 ) )

   cqss(sHCO-nspec) = min( cqss(sHCO-nspec), rho/W_sHCO )


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

	 cqss(sIXC3H7-nspec) = ( 0.0_pr &
                + k(r67f) * c(sC3H6) * c(sH) &
                + k(r79) * c(sC3H8) * c(sO) &
                + k(r82) * c(sC3H8) * c(sHO2) &
                + k(r83) * c(sC3H8) * c(sH) &
                + k(r84) * c(sC3H8) * c(sO2) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r67b)  &
               + k(r71) *  c(sO2)                 ) )

   cqss(sIXC3H7-nspec) = min( cqss(sIXC3H7-nspec), rho/W_sIXC3H7 )


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

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

	 C3H6OOH_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r72b)  & 
               + k(r86) *  c(sO2) & 
               + k(r87)  & 
                )

	 NXC3H7_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r72f) *  c(sO2) & 
               + k(r73f)  & 
               + k(r74)  & 
               + k(r75) *  c(sO2) & 
                )


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

	 C3H6OOH_ct1 = ( 0.0_pr                )

	 NXC3H7_ct1 = ( 0.0_pr &
                + k(r73b) * c(sC2H4) * c(sCH3) &
                + k(r76) * c(sC3H8) * c(sH) &
                + k(r77) * c(sC3H8) * c(sO) &
                + k(r78) * c(sC3H8) * c(sOH) &
                + k(r81) * c(sC3H8) * c(sO2) &
                + k(r85) * c(sC3H8) * c(sHO2)                )


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

	 C3H6OOH_denom2 = tiny(1.0_pr) + ( C3H6OOH_denom1                 )

	 C3H6OOH_ct2 = ( C3H6OOH_ct1 &
                ) / C3H6OOH_denom2  

	 C3H6OOH_NXC3H7 = ( 0.0_pr &
                + k(r72f)  * c(sO2) &
                ) / C3H6OOH_denom2 



   	 C3H6OOH_NXC3H7_coeff = ( 0.0_pr &
                + k(r72b)  &
                )


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

	 NXC3H7_denom2 = tiny(1.0_pr) + ( NXC3H7_denom1  &
                - C3H6OOH_NXC3H7_coeff * C3H6OOH_NXC3H7                )

	 NXC3H7_ct2 = ( NXC3H7_ct1 &
                + C3H6OOH_NXC3H7_coeff * C3H6OOH_ct2 &
                ) / NXC3H7_denom2  




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

	 cqss(sNXC3H7-nspec) = ( NXC3H7_ct2                 ) 
   !cqss(sNXC3H7-nspec) = min( cqss(sNXC3H7-nspec), rho/W_sNXC3H7 )

	 cqss(sC3H6OOH-nspec) = ( C3H6OOH_ct2 &
                + C3H6OOH_NXC3H7 * cqss(sNXC3H7-nspec)                ) 
               
   cqss(sNXC3H7-nspec) = min( cqss(sNXC3H7-nspec), rho/W_sNXC3H7 )
   cqss(sC3H6OOH-nspec) = min( cqss(sC3H6OOH-nspec), rho/W_sC3H6OOH )

      end subroutine get_QSS_clipped

! --- 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(15) :: M

      real(pr) :: NXC3H7_denom1&
         , NXC3H7_ct1&
         , NXC3H7_denom2&
         , NXC3H7_ct2
      real(pr) :: TXCH2_denom1&
         , TXCH2_ct1&
         , TXCH2_denom2&
         , TXCH2_ct2&
         , TXCH2_SXCH2&
         , TXCH2_SXCH2_coeff
      real(pr) :: SXCH2_denom1&
         , SXCH2_ct1&
         , SXCH2_denom2&
         , SXCH2_ct2
      real(pr) :: C3H6OOH_denom1&
         , C3H6OOH_ct1&
         , C3H6OOH_denom2&
         , C3H6OOH_ct2&
         , C3H6OOH_NXC3H7&
         , C3H6OOH_NXC3H7_coeff


  ! c(sTXCH2) c(sSXCH2) (coupled)  --------------------

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

	 TXCH2_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r17) *  c(sOH) & 
               + k(r18) *  c(sO2) & 
               + k(r19) *  c(sO) & 
               + k(r20) *  c(sO2) & 
               + k(r21) *  c(sO) & 
               + k(r22b) * M(mM12)  & 
                )

	 SXCH2_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r22f) * M(mM12)  & 
               + k(r23) *  c(sO2) & 
               + k(r24b) *  c(sH2O) & 
                )


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

	 TXCH2_ct1 = ( 0.0_pr &
                + k(r52) * c(sCH2CO) * c(sO)                )

	 SXCH2_ct1 = ( 0.0_pr &
                + k(r24f) * c(sCH3) * c(sOH)                )


  ! TXCH2 ---------------------------------------

	 TXCH2_denom2 = tiny(1.0_pr) + ( TXCH2_denom1                 )

	 TXCH2_ct2 = ( TXCH2_ct1 &
                ) / TXCH2_denom2  

	 TXCH2_SXCH2 = ( 0.0_pr &
                + k(r22f)  * M(mM12) &
                ) / TXCH2_denom2 



   	 TXCH2_SXCH2_coeff = ( 0.0_pr &
                + k(r22b)  * M(mM12) &
                )


  ! SXCH2 ---------------------------------------

	 SXCH2_denom2 = tiny(1.0_pr) + ( SXCH2_denom1  &
                - TXCH2_SXCH2_coeff * TXCH2_SXCH2                )

	 SXCH2_ct2 = ( SXCH2_ct1 &
                + TXCH2_SXCH2_coeff * TXCH2_ct2 &
                ) / SXCH2_denom2  




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

	 cqss(sSXCH2-nspec) = ( SXCH2_ct2                 ) 

	 cqss(sTXCH2-nspec) = ( TXCH2_ct2 &
                + TXCH2_SXCH2 * cqss(sSXCH2-nspec)                ) 

  ! cqss(sC3H5) (uncoupled) --------------------

	 cqss(sC3H5-nspec) = ( 0.0_pr &
                + k(r61b) * c(sO2) * c(sC3H6) &
                + k(r63b) * c(sC3H6) &
                + k(r66) * c(sC3H6) * c(sH) &
                + k(r68) * c(sC3H6) * c(sOH) &
                + k(r89f) * c(sC4H8) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r61f) *  c(sHO2) &
               + k(r62) *  c(sH) &
               + k(r63f) *  c(sH) &
               + k(r64) *  c(sHO2) &
               + k(r89b) *  c(sCH3)                 ) )


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

	 cqss(sC2H3-nspec) = ( 0.0_pr &
                + k(r43f) * c(sC2H4) * c(sH) &
                + k(r45f) * c(sC2H4) * c(sOH) &
                + k(r62) * c(sH) * cqss(sC3H5-nspec) &
                + k(r64) * cqss(sC3H5-nspec) *c(sHO2) &
                + k(r90) * c(sC4H8) * c(sH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r38) *  c(sO2) &
               + k(r39) *  c(sCH3) &
               + k(r40) *  c(sO2) &
               + k(r41) *  c(sH) * M(mM16) &
               + k(r43b) *  c(sH2) &
               + k(r45b) *  c(sH2O)                 ) )


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

	 cqss(sCH2CHO-nspec) = ( 0.0_pr &
                + k(r38) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r42) * c(sC2H4) * c(sO) &
                + k(r58b) * c(sH) * c(sCH2CO) &
                + k(r88) * c(sOC3H5OOH) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r55) *  c(sO2) &
               + k(r56) *  c(sO) &
               + k(r57) *  c(sCH3) &
               + k(r58f)  &
               + k(r59) *  c(sHO2) &
               + k(r60) *  c(sOH)                 ) )


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

	 cqss(sC2H5-nspec) = ( 0.0_pr &
                + k(r46b) * c(sH) * c(sC2H4) &
                + k(r49) * c(sC2H6) * c(sO) &
                + k(r50) * c(sC2H6) * c(sH) &
                + k(r51) * c(sC2H6) * c(sOH) &
                + k(r53) * c(sCH2CO) * c(sCH3) &
                + k(r57) * cqss(sCH2CHO-nspec) *c(sCH3) &
                + k(r70) * c(sC3H6) * c(sO) &
                + k(r80f) * c(sC3H8) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r46f)  &
               + k(r47) *  c(sH) &
               + k(r48) *  c(sO2) &
               + k(r80b) *  c(sCH3)                 ) )


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

	 cqss(sHCO-nspec) = ( 0.0_pr &
                + k(r32b) * c(sH) * c(sCO) * M(mM8) &
                + k(r34) * c(sCH2O) * c(sHO2) &
                + k(r35) * c(sCH2O) * c(sOH) &
                + k(r36) * c(sCH2O) * c(sH) &
                + k(r37) * c(sCH2O) * c(sO) &
                + k(r40) * cqss(sC2H3-nspec) *c(sO2) &
                + k(r44) * c(sC2H4) * c(sO) &
                + k(r56) * cqss(sCH2CHO-nspec) *c(sO) &
                + k(r59) * cqss(sCH2CHO-nspec) *c(sHO2) &
                + k(r70) * c(sC3H6) * c(sO) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r31) *  c(sOH) &
               + k(r32f) * M(mM8)  &
               + k(r33) *  c(sO2)                 ) )


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

	 cqss(sIXC3H7-nspec) = ( 0.0_pr &
                + k(r67f) * c(sC3H6) * c(sH) &
                + k(r79) * c(sC3H8) * c(sO) &
                + k(r82) * c(sC3H8) * c(sHO2) &
                + k(r83) * c(sC3H8) * c(sH) &
                + k(r84) * c(sC3H8) * c(sO2) &
                ) / ( tiny(1.0_pr) + ( &
               + k(r67b)  &
               + k(r71) *  c(sO2)                 ) )


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

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

	 C3H6OOH_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r72b)  & 
               + k(r86) *  c(sO2) & 
               + k(r87)  & 
                )

	 NXC3H7_denom1 = tiny(1.0_pr) + ( 0.0_pr &
               + k(r72f) *  c(sO2) & 
               + k(r73f)  & 
               + k(r74)  & 
               + k(r75) *  c(sO2) & 
                )


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

	 C3H6OOH_ct1 = ( 0.0_pr                )

	 NXC3H7_ct1 = ( 0.0_pr &
                + k(r73b) * c(sC2H4) * c(sCH3) &
                + k(r76) * c(sC3H8) * c(sH) &
                + k(r77) * c(sC3H8) * c(sO) &
                + k(r78) * c(sC3H8) * c(sOH) &
                + k(r81) * c(sC3H8) * c(sO2) &
                + k(r85) * c(sC3H8) * c(sHO2)                )


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

	 C3H6OOH_denom2 = tiny(1.0_pr) + ( C3H6OOH_denom1                 )

	 C3H6OOH_ct2 = ( C3H6OOH_ct1 &
                ) / C3H6OOH_denom2  

	 C3H6OOH_NXC3H7 = ( 0.0_pr &
                + k(r72f)  * c(sO2) &
                ) / C3H6OOH_denom2 



   	 C3H6OOH_NXC3H7_coeff = ( 0.0_pr &
                + k(r72b)  &
                )


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

	 NXC3H7_denom2 = tiny(1.0_pr) + ( NXC3H7_denom1  &
                - C3H6OOH_NXC3H7_coeff * C3H6OOH_NXC3H7                )

	 NXC3H7_ct2 = ( NXC3H7_ct1 &
                + C3H6OOH_NXC3H7_coeff * C3H6OOH_ct2 &
                ) / NXC3H7_denom2  




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

	 cqss(sNXC3H7-nspec) = ( NXC3H7_ct2                 ) 

	 cqss(sC3H6OOH-nspec) = ( C3H6OOH_ct2 &
                + C3H6OOH_NXC3H7 * cqss(sNXC3H7-nspec)                ) 

      end subroutine get_QSS

      SUBROUTINE GET_CQSS ( press, tempe, rho, C, CQSS )

        implicit none

        real(pr), intent(in) :: press, tempe, rho
        real(pr), dimension(nspec), intent(in) :: C
        real(pr), dimension(nqss), intent(out) :: CQSS

        real(pr), dimension(nreac) :: k
        real(pr), dimension(15) :: M

        CALL get_thirdbodies ( M, C )
        CALL get_rate_coefficients ( k, M, tempe, press )
        CALL get_QSS_clipped ( CQSS, C, k, M, rho )
      
      END SUBROUTINE

      !subroutine interne au module pour AVBP
      subroutine INTERNALAVBP(P, T, rho, c, WDOT, rrate )

      USE mod_input_param_defs, ONLY: nrrate_ARC, irrate_ARC
      USE mod_solut, ONLY: Lrrate_ARC
      USE mod_error, ONLY: print_error_and_quit

      implicit none

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

      real(pr) :: T,P,rho

      REAL(pr), DIMENSION(nrrate_ARC) :: rrate
      INTEGER :: i
      CHARACTER(LEN=strl) :: message

      CALL get_thirdbodies(M,c)
      CALL get_rate_coefficients(k,M,T,P)
      !CALL get_QSS(cqss,c,k,M)
      CALL get_QSS_clipped(cqss,c,k,M,rho)
      CALL get_reaction_rates(w,k,M,c,cqss)
      CALL get_production_rates(WDOT,w)

!     ARC reaction rate wanted
      IF ( irrate_ARC==2 ) THEN
        IF ( nrrate_ARC>nreac ) THEN
          WRITE ( message, '(A,I0)' ) "Maximum number of reaction rate demanded exceeded, check if the number of reactions in the name of your f90 file is the same as in your  limit is :", nreac
          CALL print_error_and_quit ( err_message=message )
        END IF
      ELSE IF ( irrate_ARC==1 ) THEN
        DO i=1,nrrate_ARC
          IF ( Lrrate_ARC(i)>nreac ) THEN
            WRITE ( message, '(A,I0)' ) "Maximum number of reaction rate demanded exceeded, limit is :", nreac
            CALL print_error_and_quit ( err_message=message )
          END IF
        END DO
      END IF
  
      rrate = zero
      DO i=1,nrrate_ARC
        SELECT CASE ( irrate_ARC )
        CASE (2)
          rrate(i) = w(i)
        CASE (1)
          rrate(i) = w(Lrrate_ARC(i))
        CASE DEFAULT
        END SELECT
      END DO


      end subroutine INTERNALAVBP

end module mod_C3H8_20_114_10_QM

!subroutine externe au module pour AVBP
SUBROUTINE C3H8_20_114_10_QM ( nnode,neqs,rhoinv,press,tempe,w_spec,wmol,source_spec,rrate_ARC )
      USE mod_param_defs
      USE mod_C3H8_20_114_10_QM, ONLY: INTERNALAVBP
      USE mod_input_param_defs, ONLY: chemical_ramp_niter_start, chemical_ramp_niter_end, Chemical_relax_start, Chemical_relax_end, nrrate_ARC
      USE mod_misc_defs, ONLY: nit_tot

      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
      REAL(pr), DIMENSION(1:nrrate_ARC,1:nnode), INTENT(OUT) :: rrate_ARC
      ! LOCAL
      REAL(pr), DIMENSION(1:neqs) :: C
      REAL(pr), DIMENSION(1:neqs) :: source_spec_buf
      REAL(pr), DIMENSION(1:nrrate_ARC) :: rrate_buf
      REAL(pr) :: P,T,rho
      REAL(pr) :: coef, Chemical_relax
      INTEGER :: n,k

!     Compute chemical relaxation ( CHARLELIE TIPS )
      IF ( chemical_ramp_niter_end > chemical_ramp_niter_start .AND. &
           nit_tot<chemical_ramp_niter_end ) THEN
         coef = MIN( 1.0d0, MAX( 0.0d0, DBLE(chemical_ramp_niter_end-nit_tot) / DBLE(chemical_ramp_niter_end-chemical_ramp_niter_start) ) )
         Chemical_relax = ( Chemical_relax_start - Chemical_relax_end ) * coef + Chemical_relax_end
         Chemical_relax = 1.0_pr / Chemical_relax
      ELSE
         Chemical_relax = 1.0_pr
      END IF

      DO n=1,nnode
        P = press(n)
        T = tempe(n)
        rho = one/rhoinv(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), verysmall)
        END DO

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

        !   Get the reaction rate
        DO k = 1,nrrate_ARC
          rrate_ARC(k,n) = rrate_buf(k)
        END DO

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

SUBROUTINE get_YQSS_C3H8_20_114_10_QM ( neqs, nqss, rho, press, tempe, w_spec, wmol, YQSS )

      USE mod_param_defs
      USE mod_C3H8_20_114_10_QM, ONLY : GET_CQSS

      IMPLICIT NONE

      ! IN/OUT
      INTEGER, INTENT(IN) :: neqs, nqss
      REAL(pr), INTENT(IN) :: rho, press, tempe
      REAL(pr), DIMENSION(1:neqs), INTENT(IN) :: wmol, w_spec
      REAL(pr), DIMENSION(1:nqss), INTENT(OUT) :: YQSS
      ! LOCAL
      INTEGER :: k
      REAL(pr), DIMENSION(1:neqs) :: C
      REAL(pr), DIMENSION(1:nqss) :: CQSS

        !Activity concentrations
        DO k=1,neqs
          C(k) = min( w_spec(k), rho )
          C(k) = C(k) / wmol(k)
          C(k) = max(C(k), verysmall)
        END DO

        CALL GET_CQSS ( press, tempe, rho, C, CQSS )

        DO k=1,neqs
          YQSS(k) = CQSS(k) * wmol(k) / rho
        END DO

END SUBROUTINE

SUBROUTINE get_reacname_C3H8_20_114_10_QM ( reacexp_ARC )

    USE mod_input_param_defs, ONLY: nrrate_ARC, irrate_ARC
    USE mod_solut, ONLY: Lrrate_ARC
    USE mod_C3H8_20_114_10_QM, ONLY: reacexp, reaction_expressions

    IMPLICIT NONE

    CHARACTER(LEN=65), DIMENSION(nrrate_ARC) :: reacexp_ARC
    INTEGER :: i


    CALL reaction_expressions
    DO i=1,nrrate_ARC
      IF ( irrate_ARC==1 ) THEN
        reacexp_ARC(i) = reacexp(Lrrate_ARC(i))
      ELSE
        reacexp_ARC(i) = reacexp(i)
      END IF
    END DO


END SUBROUTINE get_reacname_C3H8_20_114_10_QM