!
! ********* C1-C3 MECHANISM *********
! Ranzi, E., Frassoldati, A., Stagni, A., Pelucchi, M., Cuoci, A., Faravelli, T., Reduced kinetic schemes of complex reaction systems: Fossil and biomass-derived transportation fuels (2014) International Journal of Chemical Kinetics, 46 (9), pp. 512-542, DOI: 10.1002/kin.20867
! Ranzi, E., Cavallotti, C., Cuoci, A., Frassoldati, A., Pelucchi, M., Faravelli, T., New reaction classes in the kinetic modeling of low temperature oxidation of n-alkanes (2015) Combustion and Flame, 162 (5), pp. 1679-1691, DOI: 10.1016/j.combustflame.2014.11.030
! Bagheri, G., Ranzi, E., Pelucchi, M., Parente, A., Frassoldati, A., Faravelli, T., Comprehensive kinetic study of combustion technologies for low environmental impact: MILD and OXY-fuel combustion of methane (2020) Combustion and Flame 212, pp. 142-155, DOI: 10.1016/j.combustflame.2019.10.014
!
! ********* NOx SUB-MECHANISM *********
! Song, Y., Marrodán, L., Vin, N., Herbinet, O., Assaf, E., Fittschen, C., Stagni, A., Faravelli, T. & Battin-Leclerc, F., The sensitizing effects of NO2 and NO on methane low temperature oxidation in a jet stirred reactor. (2019) Proceedings of the Combustion Institute, 37(1), 667-675. DOI: 10.1016/j.proci.2018.06.115
!
! ********* REDUCED USING THE ARCANE SOFTWARE *********
! Cazères, Q., Pepiot, P., Riber, E., & Cuenot, B., A fully automatic procedure for the analytical reduction of chemical kinetics mechanisms for Computational Fluid Dynamics applications. (2021) Fuel, 303, 121247. DOI: 10.1016/j.fuel.2021.121247
!
! ********* REDUCED AT ETH ZURICH BY QUENTIN MALE *********
! Malé, Q., Pandey, K., & Noiray, N., The LEAF concept operated with hydrogen: Flame topology and NOx formation. (2024) Proceedings of the Combustion Institute.
!
!--------------------------------------------------------------------------------------------------
!     Copyright (c) CERFACS (all rights reserved)
!--------------------------------------------------------------------------------------------------
!     FILE ARC_S12R38QSS4.f90
!>    @file ARC_S12R38QSS4.f90
!!    Module for calculating the ARC source terms.
!!    @details 
!!    @authors Q. Malé
!!    @date    2023/10/18
!!    @since   
!!    @note    
!--------------------------------------------------------------------------------------------------
!
!--------------------------------------------------------------------------------------------------
!     MODULE mod_customkinetics
!>    @details Generated by ARCANE custom kinetics routine to compute the chemical source terms.
!!    @authors Q. Cazères, P. Pepiot
!!    @date    2019/01/24
!--------------------------------------------------------------------------------------------------
module mod_customkinetics
  implicit none

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

  ! Ideal gas constant
  real(pr), parameter :: Rcst = 8.3144621_pr
  
  ! Number of elements in the chemical system 
  integer, parameter :: ne = 3
  
  ! Number of non-qss and qss species and reactions
  integer, parameter :: nspec = 12
  integer, parameter :: nqss = 4
  integer, parameter :: nreac = 38
  integer, parameter :: nreac_reverse = 38
  
  ! Actual expression of each reaction
  character(len=65), dimension(nreac + nreac_reverse) :: reacexp
  
  ! Number of thirdbodies
  integer, parameter :: nTB = 3
  integer, parameter :: nFO = 3
  integer, parameter :: nTB_reverse = 3
  integer, parameter :: nFO_reverse = 3
  
  ! Index of elements
  integer, parameter :: eN = 1
  integer, parameter :: eH = 2
  integer, parameter :: eO = 3

  ! Index of species
  integer, parameter :: sN2 = 1
  integer, parameter :: sH2 = 2
  integer, parameter :: sH = 3
  integer, parameter :: sO2 = 4
  integer, parameter :: sO = 5
  integer, parameter :: sH2O = 6
  integer, parameter :: sOH = 7
  integer, parameter :: sHO2 = 8
  integer, parameter :: sNO = 9
  integer, parameter :: sNO2 = 10
  integer, parameter :: sN = 11
  integer, parameter :: sNH = 12

  integer, parameter :: sqssN2O = 1
  integer, parameter :: sqssHONO = 2
  integer, parameter :: sqssNNH = 3
  integer, parameter :: sqssNH2 = 4

  ! Index of reactions
  integer, parameter :: r1f = 1
  integer, parameter :: r2f = 2
  integer, parameter :: r3f = 3
  integer, parameter :: r4f = 4
  integer, parameter :: r5f = 5
  integer, parameter :: r6f = 6
  integer, parameter :: r7f = 7
  integer, parameter :: r8f = 8
  integer, parameter :: r9f = 9
  integer, parameter :: r10f = 10
  integer, parameter :: r11f = 11
  integer, parameter :: r12f = 12
  integer, parameter :: r13f = 13
  integer, parameter :: r14f = 14
  integer, parameter :: r15f = 15
  integer, parameter :: r16f = 16
  integer, parameter :: r17f = 17
  integer, parameter :: r18f = 18
  integer, parameter :: r19f = 19
  integer, parameter :: r20f = 20
  integer, parameter :: r21f = 21
  integer, parameter :: r22f = 22
  integer, parameter :: r23f = 23
  integer, parameter :: r24f = 24
  integer, parameter :: r25f = 25
  integer, parameter :: r26f = 26
  integer, parameter :: r27f = 27
  integer, parameter :: r28f = 28
  integer, parameter :: r29f = 29
  integer, parameter :: r30f = 30
  integer, parameter :: r31f = 31
  integer, parameter :: r32f = 32
  integer, parameter :: r33f = 33
  integer, parameter :: r34f = 34
  integer, parameter :: r35f = 35
  integer, parameter :: r36f = 36
  integer, parameter :: r37f = 37
  integer, parameter :: r38f = 38
  integer, parameter :: r1b = 39
  integer, parameter :: r2b = 40
  integer, parameter :: r3b = 41
  integer, parameter :: r4b = 42
  integer, parameter :: r5b = 43
  integer, parameter :: r6b = 44
  integer, parameter :: r7b = 45
  integer, parameter :: r8b = 46
  integer, parameter :: r9b = 47
  integer, parameter :: r10b = 48
  integer, parameter :: r11b = 49
  integer, parameter :: r12b = 50
  integer, parameter :: r13b = 51
  integer, parameter :: r14b = 52
  integer, parameter :: r15b = 53
  integer, parameter :: r16b = 54
  integer, parameter :: r17b = 55
  integer, parameter :: r18b = 56
  integer, parameter :: r19b = 57
  integer, parameter :: r20b = 58
  integer, parameter :: r21b = 59
  integer, parameter :: r22b = 60
  integer, parameter :: r23b = 61
  integer, parameter :: r24b = 62
  integer, parameter :: r25b = 63
  integer, parameter :: r26b = 64
  integer, parameter :: r27b = 65
  integer, parameter :: r28b = 66
  integer, parameter :: r29b = 67
  integer, parameter :: r30b = 68
  integer, parameter :: r31b = 69
  integer, parameter :: r32b = 70
  integer, parameter :: r33b = 71
  integer, parameter :: r34b = 72
  integer, parameter :: r35b = 73
  integer, parameter :: r36b = 74
  integer, parameter :: r37b = 75
  integer, parameter :: r38b = 76
  
  ! Index of third body species
  integer, parameter :: mM1 = 1
  integer, parameter :: mM5 = 2
  integer, parameter :: mM13 = 3
  
  integer, parameter :: mM12 = 4
  integer, parameter :: mM29 = 5
  integer, parameter :: mM36 = 6
  
  ! Index of third body reactions
  integer, parameter :: TBr1f = 1
  integer, parameter :: TBr5f = 2
  integer, parameter :: TBr13f = 3
  integer, parameter :: TBr1b = 4
  integer, parameter :: TBr5b = 5
  integer, parameter :: TBr13b = 6
  
  ! Index of fall off reactions
  integer, parameter :: FOr12f = 1
  integer, parameter :: FOr29f = 2
  integer, parameter :: FOr36f = 3
  integer, parameter :: FOr12b = 4
  integer, parameter :: FOr29b = 5
  integer, parameter :: FOr36b = 6
  
  ! Molar mass
  real(pr), parameter, dimension(nspec) :: W_sp =(/ &
       2.801348e-02_pr, & ! N2
       2.015880e-03_pr, & ! H2
       1.007940e-03_pr, & ! H
       3.199880e-02_pr, & ! O2
       1.599940e-02_pr, & ! O
       1.801528e-02_pr, & ! H2O
       1.700734e-02_pr, & ! OH
       3.300674e-02_pr, & ! HO2
       3.000614e-02_pr, & ! NO
       4.600554e-02_pr, & ! NO2
       1.400674e-02_pr, & ! N
       1.501468e-02_pr & ! NH
      !4.401288e-02_pr, & ! N2O
      !4.701348e-02_pr, & ! HONO
      !2.902142e-02_pr, & ! NNH
      !1.602262e-02_pr & ! NH2
       /)
  
contains

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

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

    if (concin.gt.0.0_pr) then
       conc = concin
    else
       conc = Ploc / ( Rcst * Tloc )
    end if
    
    redP = abs(k0) * conc / max(abs(kinf), tiny(1.0_pr)) + tiny(1.0_pr)

    ntmp = 0.75_pr - 1.27_pr * log10( fc )
    ccoeff = - 0.4_pr - 0.67_pr * log10( fc )
    dcoeff = 0.14_pr
    lgknull = log10(redP)
    f = (lgknull+ccoeff)/(ntmp-dcoeff*(lgknull+ccoeff))
    f = fc**(1.0_pr / ( f * f + 1.0_pr ))
    
    getlindratecoeff = kinf * f * redP / ( 1.0_pr + redP )

  end function getlindratecoeff

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

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

    M(mM1) = (1.500000e+00_pr)*c(sH2) &
         + (1.100000e+01_pr)*c(sH2O) &
         + sum(c)

    M(mM5) = (-2.700000e-01_pr)*c(sH2) &
         + (2.650000e+00_pr)*c(sH2O) &
         + sum(c)

    M(mM13) = sum(c)

    M(mM12) = (3.000000e-01_pr)*c(sH2) &
         + (9.000000e+00_pr)*c(sH2O) &
         + sum(c)

    M(mM29) = (9.000000e+00_pr)*c(sH2O) &
         + (7.000000e-01_pr)*c(sN2) &
         + (5.000000e-01_pr)*c(sO2) &
         + sum(c)

    M(mM36) = (1.100000e+01_pr)*c(sH2O) &
         + (7.000000e-01_pr)*c(sN2) &
         + (4.000000e-01_pr)*c(sO2) &
         + sum(c)

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

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

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

    k(r1f) = (4.577000e+13_pr)*exp((-4.368096e+05_pr)*R_T_inv + (-1.400000e+00_pr)*T_log )
    k(r2f) = (5.080000e-02_pr)*exp((-2.632573e+04_pr)*R_T_inv + (2.670000e+00_pr)*T_log )
    k(r3f) = (4.380000e+07_pr)*exp((-2.924616e+04_pr)*R_T_inv )
    k(r4f) = (1.140000e+08_pr)*exp((-6.395662e+04_pr)*R_T_inv )
    k(r5f) = (3.500000e+10_pr)*exp((-2.000000e+00_pr)*T_log )
    k(r6f) = (6.700000e+01_pr)*exp((-6.270477e+04_pr)*R_T_inv + (1.704000e+00_pr)*T_log )
    k(r7f) = (7.079000e+07_pr)*exp((-1.234280e+03_pr)*R_T_inv )
    k(r8f) = (1.140200e+04_pr)*exp((-2.317016e+03_pr)*R_T_inv + (1.083000e+00_pr)*T_log )
    k(r9f) = (3.250000e+07_pr)
    k(r10f) = (7.000000e+06_pr)*exp((4.572945e+03_pr)*R_T_inv )
    k(r11f) = (4.500000e+08_pr)*exp((-4.572945e+04_pr)*R_T_inv )
    k_0(FOr12f) =(1.740000e+07_pr)*exp((-1.230000e+00_pr)*T_log )
    k_inf(FOr12f) =(4.650000e+06_pr)*exp((4.400000e-01_pr)*T_log )
    FC(FOr12f) = ((1.0_pr - 6.700000e-01_pr)*exp(-Tloc/(1.000000e-30_pr)) + (6.700000e-01_pr)*exp(-Tloc/(1.000000e+30_pr)))&
              + exp(-(1.000000e+30_pr)/Tloc)
    k(r12f) = getlindratecoeff(Tloc, k_0(FOr12f), k_inf(FOr12f),FC(FOr12f), M(mM12),Ploc)
    k(r13f) = (1.000000e+04_pr)
    k(r14f) = (4.000000e+07_pr)*exp((-1.527160e+04_pr)*R_T_inv )
    k(r15f) = (6.200000e+09_pr)*exp((-1.053322e+04_pr)*R_T_inv + (-1.150000e+00_pr)*T_log )
    k(r16f) = (3.500000e+07_pr)*exp((-7.232881e+03_pr)*R_T_inv )
    k(r17f) = (7.000000e+07_pr)
    k(r18f) = (1.570000e+01_pr)*exp((2.426720e+03_pr)*R_T_inv + (1.740000e+00_pr)*T_log )
    k(r19f) = (1.000000e+07_pr)
    k(r20f) = (2.010000e+09_pr)*exp((-2.372328e+04_pr)*R_T_inv + (-1.380000e+00_pr)*T_log )
    k(r21f) = (1.800000e+08_pr)*exp((1.020896e+03_pr)*R_T_inv + (-3.510000e-01_pr)*T_log )
    k(r22f) = (2.830000e+07_pr)
    k(r23f) = (9.027000e+03_pr)*exp((-2.719600e+04_pr)*R_T_inv + (1.000000e+00_pr)*T_log )
    k(r24f) = (4.280000e+07_pr)*exp((-6.568880e+03_pr)*R_T_inv )
    k(r25f) = (1.000000e+09_pr)
    k(r26f) = (1.900000e+08_pr)*exp((9.204800e+01_pr)*R_T_inv + (-2.740000e-01_pr)*T_log )
    k(r27f) = (5.200000e+05_pr)*exp((1.711256e+03_pr)*R_T_inv + (3.880000e-01_pr)*T_log )
    k(r28f) = (2.110000e+06_pr)*exp((2.008320e+03_pr)*R_T_inv )
    k_0(FOr29f) =(4.720000e+12_pr)*exp((-6.485200e+03_pr)*R_T_inv + (-2.870000e+00_pr)*T_log )
    k_inf(FOr29f) =(1.300000e+09_pr)*exp((-7.500000e-01_pr)*T_log )
    FC(FOr29f) = ((1.0_pr - 7.500000e-01_pr)*exp(-Tloc/(1.000000e+03_pr)) + (7.500000e-01_pr)*exp(-Tloc/(1.000000e+05_pr)))&
              + exp(-(1.000000e+30_pr)/Tloc)
    k(r29f) = getlindratecoeff(Tloc, k_0(FOr29f), k_inf(FOr29f),FC(FOr29f), M(mM29),Ploc)
    k(r30f) = (3.090000e+17_pr)*exp((-6.782264e+03_pr)*R_T_inv + (-4.170000e+00_pr)*T_log )
    k(r31f) = (1.890000e-03_pr)*exp((-5.952577e+03_pr)*R_T_inv + (2.830000e+00_pr)*T_log )
    k(r32f) = (4.300000e+03_pr)*exp((-1.702888e+04_pr)*R_T_inv + (9.800000e-01_pr)*T_log )
    k(r33f) = (1.700000e+06_pr)*exp((2.175680e+03_pr)*R_T_inv )
    k(r34f) = (2.000000e+05_pr)*exp((4.426672e+03_pr)*R_T_inv + (8.400000e-01_pr)*T_log )
    k(r35f) = (3.920000e+06_pr)*exp((9.957920e+02_pr)*R_T_inv )
    k_0(FOr36f) =(6.020000e+08_pr)*exp((-2.403457e+05_pr)*R_T_inv )
    k_inf(FOr36f) =(9.900000e+10_pr)*exp((-2.422578e+05_pr)*R_T_inv )
    redP = abs(k_0(FOr36f)) * M(mM36) / max(abs(k_inf(FOr36f)), tiny(1.0_pr)) + tiny(1.0_pr)
    k(r36f) = k_inf(FOr36f) * redP / ( 1.0_pr + redP )
    k(r37f) = (5.000000e+08_pr)*exp((-7.573040e+04_pr)*R_T_inv )
    k(r38f) = (6.620000e+07_pr)*exp((-1.114199e+05_pr)*R_T_inv )
    k(r1b) = (2.598477e+08_pr)*exp((-4.054129e+03_pr)*R_T_inv + (-1.784129e+00_pr)*T_log )
    k(r2b) = (2.502565e-02_pr)*exp((-2.013635e+04_pr)*R_T_inv + (2.658110e+00_pr)*T_log )
    k(r3b) = (1.527374e+08_pr)*exp((-9.105995e+04_pr)*R_T_inv + (3.945820e-02_pr)*T_log )
    k(r4b) = (4.965349e+05_pr)*exp((6.110298e+03_pr)*R_T_inv + (3.373898e-01_pr)*T_log )
    k(r5b) = (2.149816e+16_pr)*exp((-4.945693e+05_pr)*R_T_inv + (-1.576413e+00_pr)*T_log )
    k(r6b) = (9.465100e+00_pr)*exp((5.298392e+03_pr)*R_T_inv + (1.652652e+00_pr)*T_log )
    k(r7b) = (4.315121e+04_pr)*exp((-1.538968e+05_pr)*R_T_inv + (6.251890e-01_pr)*T_log )
    k(r8b) = (3.239183e+03_pr)*exp((-2.312358e+05_pr)*R_T_inv + (1.382689e+00_pr)*T_log )
    k(r9b) = (4.548409e+06_pr)*exp((-2.227294e+05_pr)*R_T_inv + (2.877992e-01_pr)*T_log )
    k(r10b) = (6.934637e+06_pr)*exp((-2.861597e+05_pr)*R_T_inv + (3.391470e-01_pr)*T_log )
    k(r11b) = (4.457981e+08_pr)*exp((-3.364621e+05_pr)*R_T_inv + (3.391470e-01_pr)*T_log )
    k_0(FOr12b) =(1.078839e+13_pr)*exp((-2.038366e+05_pr)*R_T_inv + (-1.145560e+00_pr)*T_log )
    k_inf(FOr12b) =(2.883105e+12_pr)*exp((-2.038366e+05_pr)*R_T_inv + (5.244405e-01_pr)*T_log )
    FC(FOr12b) = ((1.0_pr - 6.700000e-01_pr)*exp(-Tloc/(1.000000e-30_pr)) + (6.700000e-01_pr)*exp(-Tloc/(1.000000e+30_pr)))&
              + exp(-(1.000000e+30_pr)/Tloc)
    k(r12b) = getlindratecoeff(Tloc, k_0(FOr12b), k_inf(FOr12b),FC(FOr12b), M(mM12),Ploc)
    k(r13b) = (1.423517e+12_pr)*exp((-2.739036e+05_pr)*R_T_inv + (-2.529494e-01_pr)*T_log )
    k(r14b) = (1.853949e+07_pr)*exp((-5.997735e+04_pr)*R_T_inv + (4.196312e-02_pr)*T_log )
    k(r15b) = (1.214733e+10_pr)*exp((-1.373899e+03_pr)*R_T_inv + (-1.298055e+00_pr)*T_log )
    k(r16b) = (8.669342e+07_pr)*exp((-1.110774e+05_pr)*R_T_inv + (7.427599e-02_pr)*T_log )
    k(r17b) = (2.400132e+09_pr)*exp((-2.995616e+05_pr)*R_T_inv + (-2.112206e-01_pr)*T_log )
    k(r18b) = (1.356091e+02_pr)*exp((-1.632316e+05_pr)*R_T_inv + (1.853734e+00_pr)*T_log )
    k(r19b) = (6.960099e+08_pr)*exp((-3.057510e+05_pr)*R_T_inv + (-1.993309e-01_pr)*T_log )
    k(r20b) = (3.001774e+08_pr)*exp((-2.532180e+05_pr)*R_T_inv + (-1.253831e+00_pr)*T_log )
    k(r21b) = (3.704061e+14_pr)*exp((-1.530077e+05_pr)*R_T_inv + (-1.448680e+00_pr)*T_log )
    k(r22b) = (7.952135e+08_pr)*exp((-2.019065e+05_pr)*R_T_inv + (-2.736069e-01_pr)*T_log )
    k(r23b) = (1.104805e+03_pr)*exp((-1.590356e+05_pr)*R_T_inv + (1.063783e+00_pr)*T_log )
    k(r24b) = (8.017840e+07_pr)*exp((-3.208938e+05_pr)*R_T_inv + (1.034297e-01_pr)*T_log )
    k(r25b) = (3.070011e+03_pr)*exp((-3.308966e+04_pr)*R_T_inv + (-2.829555e-02_pr)*T_log )
    k(r26b) = (4.556481e+13_pr)*exp((-2.016122e+05_pr)*R_T_inv + (-1.193552e+00_pr)*T_log )
    k(r27b) = (6.060014e+04_pr)*exp((-4.596444e+04_pr)*R_T_inv + (5.661281e-01_pr)*T_log )
    k(r28b) = (6.566553e+07_pr)*exp((-3.035686e+04_pr)*R_T_inv + (-2.325011e-01_pr)*T_log )
    k_0(FOr29b) =(2.091027e+22_pr)*exp((-3.127540e+05_pr)*R_T_inv + (-3.355450e+00_pr)*T_log )
    k_inf(FOr29b) =(5.759183e+18_pr)*exp((-3.062688e+05_pr)*R_T_inv + (-1.235450e+00_pr)*T_log )
    FC(FOr29b) = ((1.0_pr - 7.500000e-01_pr)*exp(-Tloc/(1.000000e+03_pr)) + (7.500000e-01_pr)*exp(-Tloc/(1.000000e+05_pr)))&
              + exp(-(1.000000e+30_pr)/Tloc)
    k(r29b) = getlindratecoeff(Tloc, k_0(FOr29b), k_inf(FOr29b),FC(FOr29b), M(mM29),Ploc)
    k(r30b) = (4.602619e+29_pr)*exp((-2.162417e+05_pr)*R_T_inv + (-5.342724e+00_pr)*T_log )
    k(r31b) = (1.141067e-05_pr)*exp((-1.089513e+05_pr)*R_T_inv + (3.529163e+00_pr)*T_log )
    k(r32b) = (1.773189e-03_pr)*exp((-3.021387e+05_pr)*R_T_inv + (2.576311e+00_pr)*T_log )
    k(r33b) = (3.579064e+04_pr)*exp((-1.626368e+05_pr)*R_T_inv + (7.386209e-01_pr)*T_log )
    k(r34b) = (3.917383e+00_pr)*exp((-1.158707e+05_pr)*R_T_inv + (1.697690e+00_pr)*T_log )
    k(r35b) = (1.762817e+04_pr)*exp((-1.893685e+05_pr)*R_T_inv + (5.203003e-01_pr)*T_log )
    k_0(FOr36b) =(7.706559e-03_pr)*exp((-7.173106e+04_pr)*R_T_inv + (8.912561e-01_pr)*T_log )
    k_inf(FOr36b) =(1.267358e+00_pr)*exp((-7.364315e+04_pr)*R_T_inv + (8.912561e-01_pr)*T_log )
    redP = abs(k_0(FOr36b)) * M(mM36) / max(abs(k_inf(FOr36b)), tiny(1.0_pr)) + tiny(1.0_pr)
    k(r36b) = k_inf(FOr36b) * redP / ( 1.0_pr + redP )
    k(r37b) = (5.554150e+02_pr)*exp((-3.336819e+05_pr)*R_T_inv + (1.263496e+00_pr)*T_log )
    k(r38b) = (1.103035e+03_pr)*exp((-2.569530e+05_pr)*R_T_inv + (8.864592e-01_pr)*T_log )
  
    return
  end subroutine get_rate_coefficients

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

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

    w(r1f) = k(r1f) * c(sH2) * m(mM1)
    w(r2f) = k(r2f) * c(sH2) * c(sO) 
    w(r3f) = k(r3f) * c(sH2) * c(sOH) 
    w(r4f) = k(r4f) * c(sH) * c(sO2) 
    w(r5f) = k(r5f) * c(sH) * c(sOH) * m(mM5)
    w(r6f) = k(r6f) * c(sH2O) * c(sO) 
    w(r7f) = k(r7f) * c(sH) * c(sHO2) 
    w(r8f) = k(r8f) * c(sH) * c(sHO2) 
    w(r9f) = k(r9f) * c(sHO2) * c(sO) 
    w(r10f) = k(r10f) * c(sHO2) * c(sOH) 
    w(r11f) = k(r11f) * c(sHO2) * c(sOH) 
    w(r12f) = k(r12f) * c(sH) * c(sO2) 
    w(r13f) = k(r13f) * c(sO) * c(sOH) * m(mM13)
    w(r14f) = k(r14f) * c(sH) * cqss(sqssNH2) 
    w(r15f) = k(r15f) * cqss(sqssNH2) * c(sNO) 
    w(r16f) = k(r16f) * c(sH) * c(sNH) 
    w(r17f) = k(r17f) * c(sNH) * c(sO) 
    w(r18f) = k(r18f) * c(sNH) * c(sOH) 
    w(r19f) = k(r19f) * c(sNH) * c(sOH) 
    w(r20f) = k(r20f) * c(sNH) * c(sO2) 
    w(r21f) = k(r21f) * c(sNH) * c(sNO) 
    w(r22f) = k(r22f) * c(sN) * c(sOH) 
    w(r23f) = k(r23f) * c(sN) * c(sO2) 
    w(r24f) = k(r24f) * c(sN) * c(sNO) 
    w(r25f) = k(r25f) * cqss(sqssNNH) 
    w(r26f) = k(r26f) * cqss(sqssNNH) * c(sO) 
    w(r27f) = k(r27f) * cqss(sqssNNH) * c(sO) 
    w(r28f) = k(r28f) * c(sHO2) * c(sNO) 
    w(r29f) = k(r29f) * c(sNO) * c(sO) 
    w(r30f) = k(r30f) * c(sNO) * c(sOH) 
    w(r31f) = k(r31f) * c(sH) * cqss(sqssHONO) 
    w(r32f) = k(r32f) * c(sH) * cqss(sqssHONO) 
    w(r33f) = k(r33f) * cqss(sqssHONO) * c(sOH) 
    w(r34f) = k(r34f) * c(sH) * c(sNO2) 
    w(r35f) = k(r35f) * c(sNO2) * c(sO) 
    w(r36f) = k(r36f) * cqss(sqssN2O) 
    w(r37f) = k(r37f) * c(sH) * cqss(sqssN2O) 
    w(r38f) = k(r38f) * cqss(sqssN2O) * c(sO) 
    w(r1b) = k(r1b) * c(sH)**2.0_pr * m(mM1)
    w(r2b) = k(r2b) * c(sH) * c(sOH) 
    w(r3b) = k(r3b) * c(sH) * c(sH2O) 
    w(r4b) = k(r4b) * c(sO) * c(sOH) 
    w(r5b) = k(r5b) * c(sH2O) * m(mM5)
    w(r6b) = k(r6b) * c(sOH)**2.0_pr 
    w(r7b) = k(r7b) * c(sOH)**2.0_pr 
    w(r8b) = k(r8b) * c(sH2) * c(sO2) 
    w(r9b) = k(r9b) * c(sO2) * c(sOH) 
    w(r10b) = k(r10b) * c(sH2O) * c(sO2) 
    w(r11b) = k(r11b) * c(sH2O) * c(sO2) 
    w(r12b) = k(r12b) * c(sHO2) 
    w(r13b) = k(r13b) * c(sHO2) * m(mM13)
    w(r14b) = k(r14b) * c(sH2) * c(sNH) 
    w(r15b) = k(r15b) * cqss(sqssNNH) * c(sOH) 
    w(r16b) = k(r16b) * c(sH2) * c(sN) 
    w(r17b) = k(r17b) * c(sH) * c(sNO) 
    w(r18b) = k(r18b) * c(sH2O) * c(sN) 
    w(r19b) = k(r19b) * c(sH2) * c(sNO) 
    w(r20b) = k(r20b) * c(sNO) * c(sOH) 
    w(r21b) = k(r21b) * c(sH) * cqss(sqssN2O) 
    w(r22b) = k(r22b) * c(sH) * c(sNO) 
    w(r23b) = k(r23b) * c(sNO) * c(sO) 
    w(r24b) = k(r24b) * c(sN2) * c(sO) 
    w(r25b) = k(r25b) * c(sH) * c(sN2) 
    w(r26b) = k(r26b) * c(sH) * cqss(sqssN2O) 
    w(r27b) = k(r27b) * c(sNH) * c(sNO) 
    w(r28b) = k(r28b) * c(sNO2) * c(sOH) 
    w(r29b) = k(r29b) * c(sNO2) 
    w(r30b) = k(r30b) * cqss(sqssHONO) 
    w(r31b) = k(r31b) * c(sH2) * c(sNO2) 
    w(r32b) = k(r32b) * c(sH2O) * c(sNO) 
    w(r33b) = k(r33b) * c(sH2O) * c(sNO2) 
    w(r34b) = k(r34b) * c(sNO) * c(sOH) 
    w(r35b) = k(r35b) * c(sNO) * c(sO2) 
    w(r36b) = k(r36b) * c(sN2) * c(sO) 
    w(r37b) = k(r37b) * c(sN2) * c(sOH) 
    w(r38b) = k(r38b) * c(sNO)**2.0_pr 
  
    return
  end subroutine get_reaction_rates
  
  ! ----------------------------------------------- !
  ! Evaluate production rates                       !
  ! ----------------------------------------------- !
  subroutine get_production_rates(cdot,w)
    implicit none

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

    cdot(sN2) = 0.0_pr &
         + w(r24f) &
         - w(r24b) &
         + w(r25f) &
         - w(r25b) &
         + w(r36f) &
         - w(r36b) &
         + w(r37f) &
         - w(r37b) 

    cdot(sH2) = 0.0_pr &
         - w(r1f) &
         + w(r1b) &
         - w(r2f) &
         + w(r2b) &
         - w(r3f) &
         + w(r3b) &
         + w(r8f) &
         - w(r8b) &
         + w(r14f) &
         - w(r14b) &
         + w(r16f) &
         - w(r16b) &
         + w(r19f) &
         - w(r19b) &
         + w(r31f) &
         - w(r31b) 

    cdot(sH) = 0.0_pr &
         + 2.0_pr * w(r1f) &
         - 2.0_pr * w(r1b) &
         + w(r2f) &
         - w(r2b) &
         + w(r3f) &
         - w(r3b) &
         - w(r4f) &
         + w(r4b) &
         - w(r5f) &
         + w(r5b) &
         - w(r7f) &
         + w(r7b) &
         - w(r8f) &
         + w(r8b) &
         - w(r12f) &
         + w(r12b) &
         - w(r14f) &
         + w(r14b) &
         - w(r16f) &
         + w(r16b) &
         + w(r17f) &
         - w(r17b) &
         + w(r21f) &
         - w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r25f) &
         - w(r25b) &
         + w(r26f) &
         - w(r26b) &
         - w(r31f) &
         + w(r31b) &
         - w(r32f) &
         + w(r32b) &
         - w(r34f) &
         + w(r34b) &
         - w(r37f) &
         + w(r37b) 

    cdot(sO2) = 0.0_pr &
         - w(r4f) &
         + w(r4b) &
         + w(r8f) &
         - w(r8b) &
         + w(r9f) &
         - w(r9b) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         - w(r12f) &
         + w(r12b) &
         - w(r20f) &
         + w(r20b) &
         - w(r23f) &
         + w(r23b) &
         + w(r35f) &
         - w(r35b) 

    cdot(sO) = 0.0_pr &
         - w(r2f) &
         + w(r2b) &
         + w(r4f) &
         - w(r4b) &
         - w(r6f) &
         + w(r6b) &
         - w(r9f) &
         + w(r9b) &
         - w(r13f) &
         + w(r13b) &
         - w(r17f) &
         + w(r17b) &
         + w(r23f) &
         - w(r23b) &
         + w(r24f) &
         - w(r24b) &
         - w(r26f) &
         + w(r26b) &
         - w(r27f) &
         + w(r27b) &
         - w(r29f) &
         + w(r29b) &
         - w(r35f) &
         + w(r35b) &
         + w(r36f) &
         - w(r36b) &
         - w(r38f) &
         + w(r38b) 

    cdot(sH2O) = 0.0_pr &
         + w(r3f) &
         - w(r3b) &
         + w(r5f) &
         - w(r5b) &
         - w(r6f) &
         + w(r6b) &
         + w(r10f) &
         - w(r10b) &
         + w(r11f) &
         - w(r11b) &
         + w(r18f) &
         - w(r18b) &
         + w(r32f) &
         - w(r32b) &
         + w(r33f) &
         - w(r33b) 

    cdot(sOH) = 0.0_pr &
         + w(r2f) &
         - w(r2b) &
         - w(r3f) &
         + w(r3b) &
         + w(r4f) &
         - w(r4b) &
         - w(r5f) &
         + w(r5b) &
         + 2.0_pr * w(r6f) &
         - 2.0_pr * w(r6b) &
         + 2.0_pr * w(r7f) &
         - 2.0_pr * w(r7b) &
         + w(r9f) &
         - w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         - w(r13f) &
         + w(r13b) &
         + w(r15f) &
         - w(r15b) &
         - w(r18f) &
         + w(r18b) &
         - w(r19f) &
         + w(r19b) &
         + w(r20f) &
         - w(r20b) &
         - w(r22f) &
         + w(r22b) &
         + w(r28f) &
         - w(r28b) &
         - w(r30f) &
         + w(r30b) &
         - w(r33f) &
         + w(r33b) &
         + w(r34f) &
         - w(r34b) &
         + w(r37f) &
         - w(r37b) 

    cdot(sHO2) = 0.0_pr &
         - w(r7f) &
         + w(r7b) &
         - w(r8f) &
         + w(r8b) &
         - w(r9f) &
         + w(r9b) &
         - w(r10f) &
         + w(r10b) &
         - w(r11f) &
         + w(r11b) &
         + w(r12f) &
         - w(r12b) &
         + w(r13f) &
         - w(r13b) &
         - w(r28f) &
         + w(r28b) 

    cdot(sNO) = 0.0_pr &
         - w(r15f) &
         + w(r15b) &
         + w(r17f) &
         - w(r17b) &
         + w(r19f) &
         - w(r19b) &
         + w(r20f) &
         - w(r20b) &
         - w(r21f) &
         + w(r21b) &
         + w(r22f) &
         - w(r22b) &
         + w(r23f) &
         - w(r23b) &
         - w(r24f) &
         + w(r24b) &
         + w(r27f) &
         - w(r27b) &
         - w(r28f) &
         + w(r28b) &
         - w(r29f) &
         + w(r29b) &
         - w(r30f) &
         + w(r30b) &
         + w(r32f) &
         - w(r32b) &
         + w(r34f) &
         - w(r34b) &
         + w(r35f) &
         - w(r35b) &
         + 2.0_pr * w(r38f) &
         - 2.0_pr * w(r38b) 

    cdot(sNO2) = 0.0_pr &
         + w(r28f) &
         - w(r28b) &
         + w(r29f) &
         - w(r29b) &
         + w(r31f) &
         - w(r31b) &
         + w(r33f) &
         - w(r33b) &
         - w(r34f) &
         + w(r34b) &
         - w(r35f) &
         + w(r35b) 

    cdot(sN) = 0.0_pr &
         + w(r16f) &
         - w(r16b) &
         + w(r18f) &
         - w(r18b) &
         - w(r22f) &
         + w(r22b) &
         - w(r23f) &
         + w(r23b) &
         - w(r24f) &
         + w(r24b) 

    cdot(sNH) = 0.0_pr &
         + w(r14f) &
         - w(r14b) &
         - w(r16f) &
         + w(r16b) &
         - w(r17f) &
         + w(r17b) &
         - w(r18f) &
         + w(r18b) &
         - w(r19f) &
         + w(r19b) &
         - w(r20f) &
         + w(r20b) &
         - w(r21f) &
         + w(r21b) &
         + w(r27f) &
         - w(r27b) 

  
    return
  end subroutine get_production_rates


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

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

    integer :: index

    real(pr) :: N2O_ct
    real(pr) :: N2O_num
    real(pr) :: N2O_denom
    real(pr) :: N2O_NNH
    real(pr) :: HONO_ct
    real(pr) :: HONO_num
    real(pr) :: HONO_denom
    real(pr) :: NNH_ct
    real(pr) :: NNH_num
    real(pr) :: NNH_denom
    real(pr) :: NNH_N2O
    real(pr) :: NNH_NH2
    real(pr) :: NH2_ct
    real(pr) :: NH2_num
    real(pr) :: NH2_denom
    real(pr) :: NH2_NNH
  
    real(pr) :: a_d_a
    real(pr) :: a_g_a
    real(pr) :: a_h_b
    real(pr) :: A_A_A
    real(pr) :: A_B_A
    real(pr) :: A_B_B

    NH2_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r14f)* c(sH) &
             + k(r15f)* c(sNO)   )

    NH2_num = ( 0.0_pr &
             + k(r14b)* c(sH2) * c(sNH)   )

    NH2_ct = NH2_num / NH2_denom

    NH2_NNH = - ( 0.0_pr &
          + k(r15b) * c(sOH)  ) / NH2_denom

    NNH_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r25f)&
             + k(r26f)* c(sO) &
             + k(r27f)* c(sO) &
             + k(r15b)* c(sOH)   )

    NNH_num = ( 0.0_pr &
             + k(r25b)* c(sH) * c(sN2) &
             + k(r27b)* c(sNH) * c(sNO)   )

    NNH_ct = NNH_num / NNH_denom

    NNH_N2O = - ( 0.0_pr &
          + k(r26b) * c(sH)  ) / NNH_denom

    NNH_NH2 = - ( 0.0_pr &
          + k(r15f) * c(sNO)  ) / NNH_denom

    N2O_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r36f)&
             + k(r37f)* c(sH) &
             + k(r38f)* c(sO) &
             + k(r21b)* c(sH) &
             + k(r26b)* c(sH)   )

    N2O_num = ( 0.0_pr &
             + k(r21f)* c(sNH) * c(sNO) &
             + k(r36b)* c(sN2) * c(sO) &
             + k(r37b)* c(sN2) * c(sOH) &
             + k(r38b) *c(sNO)** 2.0_pr  )

    N2O_ct = N2O_num / N2O_denom

    N2O_NNH = - ( 0.0_pr &
          + k(r26f) * c(sO)  ) / N2O_denom

    a_d_a =  (1.0_pr) - (N2O_NNH) * (NNH_N2O)
    a_g_a = NH2_NNH
    a_h_b =  (1.0_pr) - (NNH_NH2) * (a_g_a) / (a_d_a)
    A_A_A =  (NNH_ct) - (N2O_ct) * (NNH_N2O)
    A_B_A = NH2_ct
    A_B_B =  (A_B_A) - (A_A_A) * (a_g_a) / (a_d_a)

    cqss(sqssNH2) = ( A_B_B ) / ( a_h_b )  

    cqss(sqssNNH) = (A_A_A - (NNH_NH2) * cqss(sqssNH2)) / (a_d_a)  

    cqss(sqssN2O) = N2O_ct - (N2O_NNH) * cqss(sqssNNH)  

    HONO_denom = tiny(1.0_pr) + ( 0.0_pr &
             + k(r31f)* c(sH) &
             + k(r32f)* c(sH) &
             + k(r33f)* c(sOH) &
             + k(r30b)  )

    HONO_num = ( 0.0_pr &
             + k(r30f)* c(sNO) * c(sOH) &
             + k(r31b)* c(sH2) * c(sNO2) &
             + k(r32b)* c(sH2O) * c(sNO) &
             + k(r33b)* c(sH2O) * c(sNO2)   )

    HONO_ct = HONO_num / HONO_denom

    cqss(sqssHONO) = HONO_ct  

    cqss = max(cqss, tiny(1.0_pr))
    cqss = min(cqss, 1e03_pr)

    ! Necessary for very low concentration variations
    ! Needs to be integrated only for Cantera
    do index=1,nqss
      if (cqss(index) /= cqss(index)) then
        cqss(index) = tiny(1.0_pr)
      end if
    end do

    return
  end subroutine get_qss


  ! ----------------------------------------------- !
  ! Mass fractions to concentrations                !
  ! ----------------------------------------------- !
  subroutine y2c(y, W_sp, P, T, c)
    implicit none

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

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

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

    ! Conversion
    c = y * rho / W_sp

    return
  end subroutine y2c

end module mod_customkinetics

! ---------------------------------------------------- !
! routine to get the species source terms from the ARC !
!                                                      !
! INPUTS                                               !
!  P: Pressure                                         !
!  T: Temperature                                      !
!  y: species mass fractions                           !
! OUTPUT                                               !
!  wdot: species source terms                          !
!                                                      !
! ---------------------------------------------------- !
subroutine customkinetics(P, T, y, wdot)

  use mod_customkinetics
  implicit none

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

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

  ! Evaluate QSS concentrations and reaction rates
  call get_thirdbodies(M,c)
  call get_rate_coefficients(k, M, T, P)
  call get_qss(cqss, c, k, M)
  call get_reaction_rates(w, k, M, c, cqss)

  ! Evaluate production rates
  call get_production_rates(wdot, w)

  return
end subroutine customkinetics

