#include atoms_red
#include CBMZ.spc
#include CBMZ.eqn

#include ../boxmox/wrapper
#include ../boxmox/setup

#LOOKATALL

#INLINE F90_GLOBAL

  ! global variables used by special rate constants of this mechanism

  ! air/water concentration (molecules/cm3)
  REAL(KIND=dp)      :: c_m, c_h2o
  ! relative humidity (not used)
  REAL(KIND=dp)      :: rh

#ENDINLINE

#INLINE F90_INIT

  ! 1E-03 * (1.0/28.97) * 6.022e23 * air density
  c_m    = 2.078702e+19 * 1.25

  ! 1E-03 * (1.0/18.02) * 6.022e23 * air density, then times kg(water)/kg(air)
  c_h2o  = 3.341842e+19 * 1.25 * 0.010

#ENDINLINE

#INLINE F90_RATES
!__________________________________________________

    REAL(KIND=dp) FUNCTION ARR3( A0,B0, TEMP )
    REAL(KIND=dp),INTENT(IN) :: TEMP
    REAL(KIND=dp),INTENT(IN):: A0,B0
    ARR3 = A0 * EXP(- B0 /TEMP )
    END FUNCTION ARR3
!__________________________________________________

    REAL(KIND=dp) FUNCTION ARR3MS( A0,B0,TEMP,C_M )
    REAL(KIND=dp), INTENT(IN) :: A0,B0
    REAL(KIND=dp), INTENT(IN) :: TEMP,C_M

    ARR3MS = C_M*A0 *(TEMP/300._dp)**(-B0)
    END FUNCTION ARR3MS
!__________________________________________________

    REAL(KIND=dp) FUNCTION TROEMS(k0_300K,n,kinf_300K,m,TEMP,C_M)

    INTRINSIC LOG10

    REAL(KIND=dp), INTENT(IN) :: TEMP      ! TEMPerature [K]
    REAL(KIND=dp), INTENT(IN) :: C_M      ! air concentration [molecules/cm3]
    REAL(KIND=dp), INTENT(IN) :: k0_300K   ! low pressure limit at 300 K
    REAL(KIND=dp), INTENT(IN) :: n         ! exponent for low pressure limit
    REAL(KIND=dp), INTENT(IN) :: kinf_300K ! high pressure limit at 300 K
    REAL(KIND=dp), INTENT(IN) :: m         ! exponent for high pressure limit
    REAL(KIND=dp)             :: zt_help, k0_T, kinf_T, k_ratio

    zt_help = TEMP/300._dp
    k0_T    = k0_300K   * zt_help**(n) * C_M ! k_0   at current T
    kinf_T  = kinf_300K * zt_help**(m)        ! k_inf at current T
    k_ratio = k0_T/kinf_T
    TROEMS   = k0_T/(1._dp+k_ratio)*0.6_dp**(1._dp/(1._dp+LOG10(k_ratio)**2))

    END FUNCTION TROEMS
!__________________________________________________

    REAL(KIND=dp) FUNCTION TROEEMS(A, B, k0_300K,n,kinf_300K,m,TEMP,C_M)

    INTRINSIC LOG10

    REAL(KIND=dp), INTENT(IN) :: TEMP      ! TEMPerature [K]
    REAL(KIND=dp), INTENT(IN) :: C_M      ! air concentration [molecules/cm3]
    REAL(KIND=dp), INTENT(IN) :: k0_300K   ! low pressure limit at 300 K
    REAL(KIND=dp), INTENT(IN) :: n         ! exponent for low pressure limit
    REAL(KIND=dp), INTENT(IN) :: kinf_300K ! high pressure limit at 300 K
    REAL(KIND=dp), INTENT(IN) :: m         ! exponent for high pressure limit
    REAL(KIND=dp), INTENT(IN) :: A, B
    REAL(KIND=dp)             :: zt_help, k0_T, kinf_T, k_ratio, troe


    zt_help = TEMP/300._dp
    k0_T    = k0_300K   * zt_help**(n) * C_M ! k_0   at current T
    kinf_T  = kinf_300K * zt_help**(m)        ! k_inf at current T
    k_ratio = k0_T/kinf_T
    troe   = k0_T/(1._dp+k_ratio)*0.6_dp**(1._dp/(1._dp+LOG10(k_ratio)**2))

    TROEEMS = A * EXP( - B / TEMP) * troe
    END FUNCTION TROEEMS

!__________________________________________________


   REAL(KIND=dp) FUNCTION k46( TEMP, C_M )
   REAL(KIND=dp), INTENT(IN) :: TEMP, C_M
   REAL(KIND=dp) :: k0, k2, k3
   k0=7.2E-15_dp * EXP(785._dp/TEMP)
   k2=4.1E-16_dp * EXP(1440._dp/TEMP)
   k3=1.9E-33_dp * EXP(725._dp/TEMP)
   k46=k0+k3/(1+k3/k2)
   END FUNCTION k46
!__________________________________________________
   REAL(KIND=dp) FUNCTION RK_HO_HNO3( TEMP, C_M )
   REAL(KIND=dp), INTENT(IN) :: TEMP, C_M
   REAL(KIND=dp) :: k1, k2, k3
   k1=7.2E-15_dp * EXP(785._dp/TEMP)
   k2=1.9E-33_dp * EXP(725._dp/TEMP)
   k3=4.1E-16_dp * EXP(1440._dp/TEMP)
   RK_HO_HNO3=k1+(C_M*k2)/(1+(C_M*k2)/k3)
   END FUNCTION RK_HO_HNO3
!__________________________________________________
   REAL(KIND=dp) FUNCTION RK_2HO2( TEMP, C_M )
   REAL(KIND=dp), INTENT(IN) :: TEMP, C_M
   REAL(KIND=dp) :: k1, k2, k3
   k1=2.3E-13_dp * EXP(600._dp/TEMP)
   k2=1.7E-33_dp * EXP(1000._dp/TEMP)
   RK_2HO2=k1+(C_M*k2)
   END FUNCTION RK_2HO2
!__________________________________________________
   REAL(KIND=dp) FUNCTION RK_2HO2_H2O( TEMP, C_M )
   REAL(KIND=dp), INTENT(IN) :: TEMP, C_M
   REAL(KIND=dp) :: k1, k2, k3
   k1=2.3E-13_dp * EXP(600._dp/TEMP)
   k2=1.7E-33_dp * EXP(1000._dp/TEMP) * C_M
   k3=1.4E-21_dp * EXP(2200._dp/TEMP)
   RK_2HO2_H2O=(k1+k2)*k3
   END FUNCTION RK_2HO2_H2O
!__________________________________________________
   REAL(KIND=dp) FUNCTION RK_CO_HO( TEMP, C_M )
   REAL(KIND=dp), INTENT(IN) :: TEMP, C_M
   RK_CO_HO =1.5e-13 * (1.0 + 8.18e-23 * TEMP * C_M)
   END FUNCTION RK_CO_HO
!__________________________________________________
   REAL(KIND=dp) FUNCTION peroxy(K,TEMP,C_M)
   REAL(KIND=dp), INTENT(IN) :: TEMP,C_M
   INTEGER :: nperox, I, J, K
   PARAMETER(nperox=10)
   REAL(KIND=dp) :: Aperox(nperox,nperox),Bperox(nperox,nperox)
   REAL(KIND=dp) :: RK_PEROX(nperox,nperox)
   REAL(KIND=dp) :: RK_PARAM(nperox),SPEROX(nperox)
!
   SPEROX(1)=C(ind_CH3O2)
   SPEROX(2)=C(ind_ETHP)
   SPEROX(3)=C(ind_RO2)
   SPEROX(4)=C(ind_C2O3)
   SPEROX(5)=C(ind_ANO2)
   SPEROX(6)=C(ind_NAP)
   SPEROX(7)=C(ind_ISOPP)
   SPEROX(8)=C(ind_ISOPN)
   SPEROX(9)=C(ind_ISOPO2)
   SPEROX(10)=C(ind_XO2)
!
   Aperox(1,1)=2.5e-13
   Aperox(2,2)=6.8e-14
   Aperox(3,3)=2.9e-12
   Aperox(4,4)=8.0e-12
   Aperox(5,5)=1.0e-12
   Aperox(6,6)=5.3e-16
   Aperox(7,7)=3.1e-14
   Aperox(8,8)=3.1e-14
   Aperox(9,9)=3.1e-14
   Aperox(10,10)=3.1e-14
   Bperox(1,1)=190.
   Bperox(2,2)=0.0
   Bperox(3,3)=500.
   Bperox(4,4)=0.0
   Bperox(5,5)=0.0
   Bperox(6,6)=1980.
   Bperox(7,7)=1000.
   Bperox(8,8)=1000.
   Bperox(9,9)=1000.
   Bperox(10,10)=1000.
   DO I=1,nperox
   DO J=1,nperox
     IF(I.NE.J) THEN
       Aperox(I,J)=2.0*SQRT(Aperox(I,I)*Aperox(J,J))
       Bperox(I,J)=0.5*(Bperox(I,I)+Bperox(J,J))
     ENDIF
   ENDDO
   ENDDO
   Aperox(3,1)=1.3e-12
   Aperox(1,3)=1.3e-12
   Bperox(3,1)=640.0
   Bperox(1,3)=640.0
!
   DO I=1,nperox
     RK_PARAM(I)=0.0
   ENDDO
   DO I=1,nperox
   DO J=1,nperox
     RK_PEROX(I,J)=ARR3(Aperox(I,J),-Bperox(I,J),TEMP)
     RK_PARAM(I)=RK_PARAM(I)+RK_PEROX(I,J)*SPEROX(J)
   ENDDO
   ENDDO
   peroxy=RK_PARAM(K)
!
   END FUNCTION peroxy
!__________________________________________________
!__________________________________________________
   REAL(KIND=dp) FUNCTION peroxy_OLD(K,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10, &
     TEMP,C_M)
   REAL(KIND=dp), INTENT(IN) :: X1,X2,X3,X4,X5,X6,X7,X8,X9,X10
   REAL(KIND=dp), INTENT(IN) :: TEMP,C_M
   INTEGER :: nperox, I, J, K
   PARAMETER(nperox=10)
   REAL(KIND=dp) :: Aperox(nperox,nperox),Bperox(nperox,nperox)
   REAL(KIND=dp) :: RK_PEROX(nperox,nperox)
   REAL(KIND=dp) :: RK_PARAM(nperox),SPEROX(nperox)
!
   SPEROX(1)=X1
   SPEROX(2)=X2
   SPEROX(3)=X3
   SPEROX(4)=X4
   SPEROX(5)=X5
   SPEROX(6)=X6
   SPEROX(7)=X7
   SPEROX(8)=X8
   SPEROX(9)=X9
   SPEROX(10)=X10
!
   Aperox(1,1)=2.5e-13
   Aperox(2,2)=6.8e-14
   Aperox(3,3)=2.9e-12
   Aperox(4,4)=8.0e-12
   Aperox(5,5)=1.0e-12
   Aperox(6,6)=5.3e-16
   Aperox(7,7)=3.1e-14
   Aperox(8,8)=3.1e-14
   Aperox(9,9)=3.1e-14
   Aperox(10,10)=3.1e-14
   Bperox(1,1)=190.
   Bperox(2,2)=0.0
   Bperox(3,3)=500.
   Bperox(4,4)=0.0
   Bperox(5,5)=0.0
   Bperox(6,6)=1980.
   Bperox(7,7)=1000.
   Bperox(8,8)=1000.
   Bperox(9,9)=1000.
   Bperox(10,10)=1000.
   DO I=1,nperox
   DO J=1,nperox
     IF(I.NE.J) THEN
       Aperox(I,J)=2.0*SQRT(Aperox(I,I)*Aperox(J,J))
       Bperox(I,J)=0.5*(Bperox(I,I)+Bperox(J,J))
     ENDIF
   ENDDO
   ENDDO
   Aperox(3,1)=1.3e-12
   Aperox(1,3)=1.3e-12
   Bperox(3,1)=640.0
   Bperox(1,3)=640.0
!
   DO I=1,nperox
     RK_PARAM(I)=0.0
   ENDDO
   DO I=1,nperox
   DO J=1,nperox
     RK_PEROX(I,J)=ARR3(Aperox(I,J),-Bperox(I,J),TEMP)
     RK_PARAM(I)=RK_PARAM(I)+RK_PEROX(I,J)*SPEROX(J)
   ENDDO
   ENDDO
   peroxy_OLD=RK_PARAM(K)
!
   END FUNCTION peroxy_OLD
!__________________________________________________

#ENDINLINE