#include atoms_red
#include MOZART_T1.spc
#include MOZART_T1.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
  REAL(KIND=dp)      :: rh

  ! SOA formation
  REAL(KIND=dp)      :: nume, den
  INTEGER, PARAMETER :: vbs_nbin = 4, vbs_nspec = 9,              &
                        vbs_alk4 = 1, vbs_alk5 = 2, vbs_ole1 = 3, &
                        vbs_ole2 = 4, vbs_aro1 = 5, vbs_aro2 = 6, &
                        vbs_isop = 7, vbs_sesq = 8, vbs_terp = 9, &
                        vbs_c1   = 1, vbs_c10  = 2, vbs_c100 = 3, &
                        vbs_c1000 = 4

#ENDINLINE

#INLINE F90_INIT

  ! relative humidity (not really used...)
  rh = 0.5

#ENDINLINE

#INLINE F90_RCONST

  c_m = C(ind_M)
  c_h2o = C(ind_H2O)

#ENDINLINE

#INLINE F90_RATES

REAL(kind=dp) FUNCTION JPL_TROE( k0_300K, n, kinf_300K, m, base, temp, cair )

!------------------------------------------------------------
! ... dummy arguments
!------------------------------------------------------------
    REAL(kind=dp), INTENT(IN) :: base      ! base expononent
    REAL(kind=dp), INTENT(IN) :: temp      ! temperature [K]
    REAL(kind=dp), INTENT(IN) :: cair      ! 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

!------------------------------------------------------------
! ... local variables
!------------------------------------------------------------
    REAL(kind=dp)  :: zt_help, k0_T, kinf_T, k_ratio

    zt_help = 300._dp/temp
    k0_T    = k0_300K   * zt_help**(n) * cair ! k_0   at current T
    kinf_T  = kinf_300K * zt_help**(m)        ! k_inf at current T
    k_ratio = k0_T/kinf_T

    JPL_TROE = k0_T/(1._dp + k_ratio)*base**(1._dp/(1._dp + LOG10(k_ratio)**2))

END FUNCTION JPL_TROE

REAL(KIND=dp) FUNCTION usr_CO_OH_a( temp, c_m )
! for cesm-consistent reaction labels, equivalent to usr8
! CO+OH -> CO2+HO2

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    REAL(KIND=dp), parameter :: boltz = 1.38044e-16_dp

    usr_CO_OH_a = 1.5e-13_dp * (1._dp + 6.e-7_dp*boltz*c_m*temp)

END FUNCTION usr_CO_OH_a

REAL(KIND=dp) FUNCTION usr_O_O2( temp )
! for cesm-consistent reaction labels
! O+O2+M -> O3+M

    REAL(KIND=dp), INTENT(IN) :: temp

    usr_O_O2 = 6.00e-34_dp*(temp/300._dp)**(-2.4_dp)

END FUNCTION usr_O_O2

REAL(KIND=dp) FUNCTION usr_O_O_M( temp )
! for cesm-consistent reaction labels
! O+O+M->O2+M

    REAL(KIND=dp), INTENT(IN) :: temp

    usr_O_O_M = 2.76e-34_dp * exp( 720.0_dp/temp)

END FUNCTION usr_O_O_M

REAL(KIND=dp) FUNCTION usr_HO2_HO2( temp, c_m, c_h2o )
! for cesm-consistent reaction labels, equivalent to usr9
! HO2+HO2 -> H2O2+O2
! H2O included in fc calculation, not as reactant

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m
    REAL(KIND=dp), INTENT(IN) :: c_h2o

    REAL(KIND=dp) :: ko, kinf, fc

    if( c_h2o > 0._dp ) then
       ko   = 2.3e-13_dp * exp( 600._dp/temp )
       kinf = 1.7e-33_dp * c_m * exp( 1000._dp/temp )
       fc = 1._dp + 1.4e-21_dp *c_h2o* exp( 2200._dp/temp )
       usr_HO2_HO2 = (ko + kinf) * fc
    else
       usr_HO2_HO2 = 0._dp
    end if
END FUNCTION usr_HO2_HO2

REAL(KIND=dp) FUNCTION usr_HNO3_OH( temp, c_m )
! for cesm-consistent reaction labels, equivalent to usr5
! HNO3 + OH -> NO3 + H2O

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    REAL(KIND=dp) :: k0, k2

   k0 = c_m * 6.5e-34_dp * exp( 1335._dp/temp )
   k2 = exp( 2199._dp/temp )
   k0 = k0 /(1.0_dp + k0/(2.7e-17_dp*k2))
   k2 = exp( 460._dp/temp )

   usr_HNO3_OH = k0 + 2.4e-14_dp * k2

END FUNCTION usr_HNO3_OH


REAL(KIND=dp) FUNCTION usr_HO2NO2_M( temp, c_m )
! for cesm-consistent reaction labels
! HO2NO2 + M -> HO2+NO2+M
! CESM writes this as dependent on the NO2+HO2 rate
! *** if M is included in reactants, must divide TROEE function by [M]

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_HO2NO2_M = TROEE( 4.76e26_dp,10900._dp, 1.8e-31_dp , 3.2_dp , &
                          4.7e-12_dp , 1.4_dp , TEMP, C_M ) /c_m

END FUNCTION usr_HO2NO2_M

REAL(KIND=dp) FUNCTION usr_PAN_M( temp, c_m )
! for cesm-consistent reaction labels
! PAN+M -> CH3CO3+NO2
! *** if M is included in reactants, must divide TROEE function by [M]

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_PAN_M = TROEE(1.111e28_dp, 14000._dp, 8.5e-29_dp, 6.5_dp, &
                       1.1e-11_dp, 0._dp, TEMP, C_M) /c_m

END FUNCTION usr_PAN_M

REAL(KIND=dp) FUNCTION usr_CH3COCH3_OH( temp )
! for cesm-consistent reaction labels
! CH3COCH3 + OH -> RO2 + H2O

    REAL(KIND=dp), INTENT(IN) :: temp

    usr_CH3COCH3_OH = 3.82e-11_dp*exp( -2000._dp/temp ) + 1.33e-13_dp

END FUNCTION usr_CH3COCH3_OH

REAL(KIND=dp) FUNCTION usr_MCO3_NO2( temp, c_m )
! for cesm-consistent reaction labels
! MCO3+NO2+M -> MPAN+M

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_MCO3_NO2 = 1.1e-11_dp*300._dp/(temp*c_m)

END FUNCTION usr_MCO3_NO2


REAL(KIND=dp) FUNCTION usr_MPAN_M( temp, c_m )
! for cesm-consistent reaction labels
! MPAN+M -> MCO3+NO2+M

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_MPAN_M = 1.2221e17_dp*300._dp*exp( -14000._dp/temp )/(temp*c_m)

END FUNCTION usr_MPAN_M

REAL(KIND=dp) FUNCTION usr_N2O5_aer( rh, temp )
! for cesm-consistent reaction labels, equivalent to usr16
! heterogeneous uptake on aerosols: N2O5 -> 2 HNO3
! **** needs code for running with WRF-Chem!!! ****

    REAL(KIND=dp), INTENT(IN) :: rh                       ! relative humidity
    REAL(KIND=dp), INTENT(IN) :: temp                     ! temperature (K)

    usr_N2O5_aer = 0._dp

END FUNCTION usr_N2O5_aer

REAL(KIND=dp) FUNCTION usr_NO3_aer( rh, temp )
! for cesm-consistent reaction labels, equivalent to usr17
! heterogeneous uptake on aerosols: NO3 -> HNO3
! **** needs code for running with WRF-Chem!!! ****

    REAL(KIND=dp), INTENT(IN) :: rh                       ! relative humidity
    REAL(KIND=dp), INTENT(IN) :: temp                     ! temperature (K)

    usr_NO3_aer = 0._dp

END FUNCTION usr_NO3_aer


REAL(KIND=dp) FUNCTION usr_NO2_aer( rh, temp )
! for cesm-consistent reaction labels, equivalent to usr17a
! heterogeneous uptake on aerosols: NO2 -> 0.5 OH + 0.5 NO + 0.5 HNO3
! **** needs code for running with WRF-Chem!!! ****

    REAL(KIND=dp), INTENT(IN) :: rh                       ! relative humidity
    REAL(KIND=dp), INTENT(IN) :: temp                     ! temperature (K)

    usr_NO2_aer = 0._dp

END FUNCTION usr_NO2_aer

REAL(KIND=dp) FUNCTION usr_DMS_OH( temp, c_m )
! for cesm-consistent reaction labels, equivalent to usr24
! DMS + OH -> 0.5 SO2 + 0.5 HO2

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    REAL(KIND=dp) :: ko, wrk

    wrk   = .21_dp*c_m
    ko    = 1._dp + 5.5e-31_dp*exp( 7460._dp/temp )*wrk
    usr_DMS_OH = 1.7e-42_dp*exp( 7810._dp/temp )*wrk/ko

END FUNCTION usr_DMS_OH


REAL(KIND=dp) FUNCTION usr_HO2_aer( rh, temp )
! for cesm-consistent reaction labels, equivalent to usr26
! heterogeneous uptake on aerosols: HO2 -> 0.5 H2O2
! **** needs code for running with WRF-Chem!!! ****

    REAL(KIND=dp), INTENT(IN) :: rh                       ! relative humidity
    REAL(KIND=dp), INTENT(IN) :: temp                     ! temperature (K)

    usr_HO2_aer = 0._dp

END FUNCTION usr_HO2_aer

REAL(KIND=dp) FUNCTION usr_SO2_OH( temp, c_m )
! for cesm-consistent reaction labels, equivalent to usr23
! SO2 + OH -> SO4

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    REAL(KIND=dp) :: fc, k0
    REAL(KIND=dp) :: wrk

    fc    = 3.e-11_dp * (300._dp/temp) ** 3.3_dp
    wrk   = fc * c_m
    k0    = wrk / (1._dp + wrk/1.5e-12_dp)
    usr_SO2_OH = k0 * .6_dp ** (1._dp/(1._dp + &
                 (log10( wrk/1.5e-12_dp ))**2._dp))

END FUNCTION usr_SO2_OH

REAL(KIND=dp) FUNCTION usr_PBZNIT_M( temp, c_m )
! for cesm-consistent reaction labels
! PBZNIT+M=ACBZO2+NO2+M
! *** if M is included in reactants, must divide TROEE function by [M]

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_PBZNIT_M = TROEE( 1.111e28_dp,14000._dp, 9.7e-29_dp , 5.6_dp , &
                          9.3e-12_dp , 0._dp , TEMP, C_M) /c_m

END FUNCTION usr_PBZNIT_M

REAL(KIND=dp) FUNCTION usr_N2O5_M( temp, c_m )
! for cesm-consistent reaction labels
! N2O5 + M -> NO2 + NO3 + M
! CESM writes this as dependent on the NO2+NO3 rate
! *** if M is included in reactants, must divide TROEE function by [M]

    REAL(KIND=dp), INTENT(IN) :: temp
    REAL(KIND=dp), INTENT(IN) :: c_m

    usr_N2O5_M = TROEE(3.333e26_dp, 10900._dp, 2.2e-30_dp, 4.4_dp, &
                       1.4e-12_dp , .7_dp , TEMP, C_M ) /c_m

END FUNCTION usr_N2O5_M


!__________________________________________________

    REAL(KIND=dp) FUNCTION Keff ( A0,B0,C0, TEMP,X1,X2,y1,y2 )
    REAL(KIND=dp),INTENT(IN) :: X1,X2,y1,y2
    REAL(KIND=dp),INTENT(IN) :: TEMP
    REAL(KIND=dp),INTENT(IN):: A0,B0,C0
    Keff = A0 * EXP(- B0 /TEMP ) &
      *(TEMP/300._dp)**C0*(y1*X1/(X1 + X2 + 1.0e-35) &
       +y2*(1-X1/(X1 + X2 + 1.0e-35)))
    END FUNCTION Keff
!__________________________________________________

    REAL(KIND=dp) FUNCTION Keff2 ( C0,X1,X2,y1,y2 )
    REAL(KIND=dp),INTENT(IN) :: X1,X2,y1,y2
    REAL(KIND=dp),INTENT(IN):: C0
    Keff2 = C0*(y1*X1/(X1 + X2 + 1.0e-35) &
       +y2*(1-X1/(X1 + X2 + 1.0e-35 )))
    END FUNCTION Keff2

!__________________________________________________


    REAL(KIND=dp) FUNCTION vbs_yield ( nume, den, voc_idx, bin_idx )
    REAL(KIND=dp), INTENT(IN) :: nume, den
    INTEGER, INTENT(IN)       :: voc_idx, bin_idx
    INTEGER, PARAMETER        :: vbs_nbin = 4, vbs_nspec = 9

    ! normalized (1 g/m3 density) yield for condensable vapors
    REAL(KIND=dp)             :: vbs_alphlowN(vbs_nbin,vbs_nspec)
    REAL(KIND=dp)             :: vbs_alphhiN(vbs_nbin,vbs_nspec)
    REAL(KIND=dp)             :: vbs_mw_prec(vbs_nspec)
    ! SOA density (g/m3)
    REAL(KIND=dp), PARAMETER  :: dens_aer = 1.5
    ! SOA molecular weight (g/mol)
    REAL(KIND=dp), PARAMETER  :: mw_aer   = 250.0

    ! --------------------------------------------------------------------------
    ! Yields used in Murphy and Pandis, 2009; Tsimpidi et al., 2010;
    ! Ahmadov et al., 2012

    ! low NOx condition
    DATA vbs_alphlowN /   &
    0.0000, 0.0750, 0.0000, 0.0000,   & ! ALK4
    0.0000, 0.3000, 0.0000, 0.0000,   & ! ALK5
    0.0045, 0.0090, 0.0600, 0.2250,   & ! OLE1
    0.0225, 0.0435, 0.1290, 0.3750,   & ! OLE2
    0.0750, 0.2250, 0.3750, 0.5250,   & ! ARO1
    0.0750, 0.3000, 0.3750, 0.5250,   & ! ARO2
    0.0090, 0.0300, 0.0150, 0.0000,   & ! ISOP
    0.0750, 0.1500, 0.7500, 0.9000,   & ! SESQ
    0.1073, 0.0918, 0.3587, 0.6075/     ! TERP

    ! high NOx condition
    DATA vbs_alphhiN /    &
    0.0000, 0.0375, 0.0000, 0.0000,   & ! ALK4
    0.0000, 0.1500, 0.0000, 0.0000,   & ! ALK5
    0.0008, 0.0045, 0.0375, 0.1500,   & ! OLE1
    0.0030, 0.0255, 0.0825, 0.2700,   & ! OLE2
    0.0030, 0.1650, 0.3000, 0.4350,   & ! ARO1
    0.0015, 0.1950, 0.3000, 0.4350,   & ! ARO2
    0.0003, 0.0225, 0.0150, 0.0000,   & ! ISOP
    0.0750, 0.1500, 0.7500, 0.9000,   & ! SESQ
    0.0120, 0.1215, 0.2010, 0.5070/     ! TERP

    DATA vbs_mw_prec /    &
    120.0, & ! ALK4
    150.0, & ! ALK5
    120.0, & ! OLE1
    120.0, & ! OLE2
    150.0, & ! ARO1
    150.0, & ! ARO2
    136.0, & ! ISOP
    250.0, & ! SESQ
    180.0/   ! TERP

    REAL(KIND=dp), PARAMETER  :: yields_dens_aer = 1.5 ! g/m3

    !
    ! --------------------------------------------------------------------------

    REAL(KIND=dp)             :: B, mw_ratio, dens_ratio

    ! Lane et al., ES&T, 2008
    ! B = (RO2 + NO) / ((RO2 + NO) + (RO2 + RO2) + (RO2 + HO2))
    ! with nume = (RO2 + NO) and den = (RO2 + RO2) + (RO2 + HO2)
    B = nume / (nume + den + 1.0e-35_dp)

    ! we need molar yields, not mass yields
    mw_ratio = vbs_mw_prec(voc_idx)/mw_aer

    ! density correction
    dens_ratio = dens_aer / yields_dens_aer

    vbs_yield = (vbs_alphhiN(bin_idx,voc_idx)  * B +             &
                 vbs_alphlowN(bin_idx,voc_idx) * (1.0_dp - B)) * &
                 dens_ratio * mw_ratio

    END FUNCTION vbs_yield


#ENDINLINE