! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! The Reaction Rates File
! 
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! 
! File                 : HETCHEM_Rates.f90
! Time                 : Thu Oct  8 08:22:02 2020
! Working directory    : /Users/lechriso/git/boxmox/boxmox/tmp_HETCHEM_12963_BOXMOX
! Equation file        : HETCHEM.kpp
! Output root filename : HETCHEM
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



MODULE HETCHEM_Rates

  USE HETCHEM_Parameters
  USE HETCHEM_Global
  USE HETCHEM_Util
  IMPLICIT NONE

CONTAINS



! Begin Rate Law Functions from KPP_HOME/util/UserRateLaws

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  User-defined Rate Law functions
!  Note: the default argument type for rate laws, as read from the equations file, is single precision
!        but all the internal calculations are performed in double precision
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!~~~>  Arrhenius
!   REAL(kind=dp) FUNCTION ARR( A0,B0,C0 )
!      REAL A0,B0,C0
!      ARR =  DBLE(A0) * EXP(-DBLE(B0)/TEMP) * (TEMP/300.0_dp)**DBLE(C0)
!   END FUNCTION ARR

! Arrhenius (added by psaide 15-07-2009)
    REAL(kind=dp) FUNCTION ARR( A0,B0,C0, TEMP )
      REAL(kind=dp) :: TEMP
      REAL(kind=dp) A0,B0,C0
      ARR =  A0 * EXP( -B0 /TEMP ) * (TEMP/300._dp)**C0
    END FUNCTION ARR

!~~~> Simplified Arrhenius, with two arguments
!~~~> Note: The argument B0 has a changed sign when compared to ARR
!   REAL(kind=dp) FUNCTION ARR2( A0,B0 )
!      REAL A0,B0
!      ARR2 =  DBLE(A0) * EXP( DBLE(B0)/TEMP )
!   END FUNCTION ARR2
! CK 20140313 WRF-chem uses different kind of Arrhenius function
!             with TEMP as argument...
   REAL(kind=dp) FUNCTION ARR2( A0,B0, TEMP )
      REAL(kind=dp) :: TEMP
      REAL(kind=dp) A0,B0
      ARR2 = A0 * EXP( -B0 /TEMP )
   END FUNCTION ARR2
! CK 20140313 end

!   REAL(kind=dp) FUNCTION EP2(A0,C0,A2,C2,A3,C3)
!      REAL A0,C0,A2,C2,A3,C3
!      REAL(kind=dp) K0,K2,K3
!      K0 = DBLE(A0) * EXP(-DBLE(C0)/TEMP)
!      K2 = DBLE(A2) * EXP(-DBLE(C2)/TEMP)
!      K3 = DBLE(A3) * EXP(-DBLE(C3)/TEMP)
!      K3 = K3*CFACTOR*1.0E6_dp
!      EP2 = K0 + K3/(1.0_dp+K3/K2 )
!   END FUNCTION EP2

!   REAL(kind=dp) FUNCTION EP3(A1,C1,A2,C2)
!      REAL A1, C1, A2, C2
!      REAL(kind=dp) K1, K2
!      K1 = DBLE(A1) * EXP(-DBLE(C1)/TEMP)
!      K2 = DBLE(A2) * EXP(-DBLE(C2)/TEMP)
!      EP3 = K1 + K2*(1.0E6_dp*CFACTOR)
!   END FUNCTION EP3

! EP2 (added by psaide 22-07-2009)
   REAL(kind=dp) FUNCTION EP2( A0,C0,A2,C2,A3,C3,TEMP,cair)
      REAL(kind=dp) :: TEMP
      REAL(kind=dp) :: cair
      REAL(kind=dp) A0,C0,A2,C2,A3,C3
      REAL(kind=dp) K0,K2,K3

      K0 = A0 * EXP(-C0 /TEMP)
      K2 = A2 * EXP(-C2 /TEMP)
      K3 = A3 * EXP(-C3 /TEMP)
!      K3 = K3 * 2.45E13_dp * 1.0E6_dp
      K3 = K3 * cair
      EP2 = K0 + K3/(1._dp+K3/K2 )
   END FUNCTION EP2

! EP3 (added by psaide 15-07-2009)
    REAL(kind=dp) FUNCTION EP3(A1,C1,A2,C2,TEMP,cair)
      REAL(kind=dp) :: TEMP
      REAL(kind=dp) :: cair
      REAL(kind=dp) A1, C1, A2, C2
      REAL(kind=dp) K1, K2

      K1 = A1 * EXP(-C1 /TEMP)
      K2 = A2 * EXP(-C2 /TEMP)
!      EP3 = K1 + K2*(1.0E6_dp * 2.45E13_dp)
      EP3 = K1 + K2*cair
    END FUNCTION EP3


!   REAL(kind=dp) FUNCTION FALL ( A0,B0,C0,A1,B1,C1,CF)
!      REAL A0,B0,C0,A1,B1,C1,CF
!      REAL(kind=dp) K0, K1
!      K0 = DBLE(A0) * EXP(-DBLE(B0)/TEMP)* (TEMP/300.0_dp)**DBLE(C0)
!      K1 = DBLE(A1) * EXP(-DBLE(B1)/TEMP)* (TEMP/300.0_dp)**DBLE(C1)
!      K0 = K0*CFACTOR*1.0E6_dp
!      K1 = K0/K1
!      FALL = (K0/(1.0_dp+K1))*   &
!           DBLE(CF)**(1.0_dp/(1.0_dp+(LOG10(K1))**2))
!   END FUNCTION FALL

! FALL (added by psaide 15-07-2009)
    REAL(kind=dp) FUNCTION FALL( A0,B0,C0,A1,B1,C1,CF,TEMP,cair)

      INTRINSIC LOG10

      REAL(kind=dp) :: TEMP
      REAL(kind=dp) :: cair
      REAL(kind=dp) A0,B0,C0,A1,B1,C1,CF
      REAL(kind=dp) K0, K1

      K0 = A0 * EXP(-B0 /TEMP)* (TEMP/300._dp)**C0
      K1 = A1 * EXP(-B1 /TEMP)* (TEMP/300._dp)**C1
!      K0 = K0 * 2.45E13_dp * 1.0E6_dp
      K0 = K0 * cair
      K1 = K0/K1
      FALL = (K0/(1._dp+K1))*CF**(1._dp/(1._dp+(LOG10(K1))**2))

    END FUNCTION FALL

! FALL2 (added by psaide 28-10-2009)
    REAL(KIND=dp) FUNCTION F2( A0,B0,C0,A1,B1,C1,CF,CN,TEMP,cair)

      INTRINSIC LOG10

      REAL(KIND=dp) :: TEMP
      REAL(KIND=dp) :: cair
      REAL(KIND=dp) A0,B0,C0,A1,B1,C1,CF,CN
      REAL(KIND=dp) K0, K1

      K0 = A0 * EXP(-B0 /TEMP)* (TEMP/300._dp)**C0
      K1 = A1 * EXP(-B1 /TEMP)* (TEMP/300._dp)**C1
!      K0 = K0 * 2.45E13_dp * 1.0E6_dp
      K0 = K0 * cair
      K1 = K0/K1
      F2 = (K0/(1._dp+K1))*CF**(1._dp/(1._dp+(LOG10(K1)/CN)**2))

    END FUNCTION F2


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

  ELEMENTAL REAL(kind=dp) FUNCTION k_3rd(temp,cair,k0_300K,n,kinf_300K,m,fc)

    INTRINSIC LOG10

    REAL(kind=dp), INTENT(IN) :: temp      ! temperature [K]
    REAL(kind=dp), INTENT(IN) :: cair      ! air concentration [molecules/cm3]
    REAL, INTENT(IN) :: k0_300K   ! low pressure limit at 300 K
    REAL, INTENT(IN) :: n         ! exponent for low pressure limit
    REAL, INTENT(IN) :: kinf_300K ! high pressure limit at 300 K
    REAL, INTENT(IN) :: m         ! exponent for high pressure limit
    REAL, INTENT(IN) :: fc        ! broadening factor (usually fc=0.6)
    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
    k_3rd   = k0_T/(1._dp+k_ratio)*fc**(1._dp/(1._dp+LOG10(k_ratio)**2))

  END FUNCTION k_3rd

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

  ELEMENTAL REAL(kind=dp) FUNCTION k_arr (k_298,tdep,temp)
    ! Arrhenius function

    REAL,     INTENT(IN) :: k_298 ! k at T = 298.15K
    REAL,     INTENT(IN) :: tdep  ! temperature dependence
    REAL(kind=dp), INTENT(IN) :: temp  ! temperature

    INTRINSIC EXP

    k_arr = k_298 * EXP(tdep*(1._dp/temp-3.3540E-3_dp)) ! 1/298.15=3.3540e-3

  END FUNCTION k_arr

! CK 20140313 in the following some more rate functions needed by
!             WRF-chem packages

!------------------------------------
! Troe reactions (as in Stockwell et al, 1997)

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

    INTRINSIC LOG10

    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
    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
    TROE   = k0_T/(1._dp+k_ratio)*0.6_dp**(1._dp/(1._dp+LOG10(k_ratio)**2))

   END FUNCTION TROE



!-------------------------------------------
! Troe equilibrium reactions (as in Stockwell et al, 1997)

    REAL(kind=dp) FUNCTION TROEE(A, B, k0_300K,n,kinf_300K,m,temp,cair)

    INTRINSIC LOG10

    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
    REAL(kind=dp),     INTENT(IN) :: A, B
    REAL(kind=dp)             :: zt_help, k0_T, kinf_T, k_ratio, troe


    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
    troe   = k0_T/(1._dp+k_ratio)*0.6_dp**(1._dp/(1._dp+LOG10(k_ratio)**2))

    TROEE = A * EXP( - B / temp) * troe



  END FUNCTION TROEE

!------------------------
! k=T^2 C exp (-D/T) reactions

   REAL(kind=dp) FUNCTION THERMAL_T2(c, d ,temp)
    REAL(kind=dp), INTENT(IN) :: temp      ! temperature [K]
    REAL(kind=dp),     INTENT(IN) :: c, d


     THERMAL_T2= temp**2._dp * c * EXP(- d  / temp)

   END FUNCTION THERMAL_T2

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  End of User-defined Rate Law functions
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! End Rate Law Functions from KPP_HOME/util/UserRateLaws


! Begin INLINED Rate Law Functions



! ****************************************************************************
! J - the photolysis rate retrieval function
! ****************************************************************************
      REAL(kind=dp) FUNCTION J(name)

      CHARACTER(LEN=*), INTENT(IN)     :: name

      CHARACTER(LEN=15) :: xname

      INTEGER                          :: i

      xname(:)           = " "
      xname(1:LEN(name)) = name

      J                  = jvalues(hash(xname))

      RETURN

      END FUNCTION J
! ****************************************************************************



REAL(KIND=dp) FUNCTION SFC_REACTION(TEMP, SAD, MW, GAMMA)

  REAL(KIND=dp), INTENT(IN)  :: TEMP, & ! temperature [K]
                                SAD, &  ! aerosol surface area [m2/m3]
                                MW, &   ! molecular weight [kg/mol]
                                GAMMA   ! uptake coeff. [-]

  REAL(KIND=dp)              :: OMEGA

  ! mean molecular velocity (m/s)
  OMEGA = 1.455e2 * SQRT(TEMP / (MW * 1e3))

  ! first order uptake, Fuchs and Sutugin, 1971
  ! dCg = 1/4 * gamma * |v_mol| * A
  SFC_REACTION = 0.25 * GAMMA * OMEGA * SAD

END FUNCTION SFC_REACTION

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

!------------------------------------------------------------
!	... dummy arguments
!------------------------------------------------------------
    REAL(kind=dp), INTENT(IN) :: base      ! base exponent
    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_N2O5_M( temp, c_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


! End INLINED Rate Law Functions

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Update_SUN - update SUN light using TIME
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  SUBROUTINE Update_SUN()
      !USE HETCHEM_Parameters
      !USE HETCHEM_Global

    IMPLICIT NONE

    REAL(kind=dp) :: SunRise, SunSet
    REAL(kind=dp) :: Thour, Tlocal, Ttmp 
    ! PI - Value of pi
    REAL(kind=dp), PARAMETER :: PI = 3.14159265358979d0
    
    SunRise = 4.5_dp 
    SunSet  = 19.5_dp 
    Thour = TIME/3600.0_dp 
    Tlocal = Thour - (INT(Thour)/24)*24

    IF ((Tlocal>=SunRise).AND.(Tlocal<=SunSet)) THEN
       Ttmp = (2.0*Tlocal-SunRise-SunSet)/(SunSet-SunRise)
       IF (Ttmp.GT.0) THEN
          Ttmp =  Ttmp*Ttmp
       ELSE
          Ttmp = -Ttmp*Ttmp
       END IF
       SUN = ( 1.0_dp + COS(PI*Ttmp) )/2.0_dp 
    ELSE
       SUN = 0.0_dp 
    END IF

 END SUBROUTINE Update_SUN

! End of Update_SUN function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Update_RCONST - function to update rate constants
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Update_RCONST ( )




! Begin INLINED RCONST



  c_m = C(ind_M)


! End INLINED RCONST

  RCONST(1) = (ARR2(1.200000e-13_dp,2450.00_dp,TEMP))
  RCONST(2) = (JPL_TROE(2.400000e-30_dp,3.00_dp,1.600000e-12_dp,-0.10_dp,0.60_dp,TEMP,C_M))
  RCONST(3) = (usr_N2O5_M(temp,c_m))
  RCONST(4) = (SFC_REACTION(TEMP,AER_SAD,0.108009_dp,1.0e-3))
      
END SUBROUTINE Update_RCONST

! End of Update_RCONST function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Update_PHOTO - function to update photolytical rate constants
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Update_PHOTO ( )


   USE HETCHEM_Global

      
END SUBROUTINE Update_PHOTO

! End of Update_PHOTO function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



END MODULE HETCHEM_Rates

