module pfunc
implicit none

contains

!=======================================================================
! Reaction rate function
!
! version: April 2019 - Yassin Jaziri
!=======================================================================


function pfunc1(nlayer,t,dens,param)

use types_asis

implicit none

! input/output

integer                 :: nlayer       ! number of atmospheric layers
real, dimension(nlayer) :: pfunc1       ! output values
type(rtype1)            :: param        ! parameters
real, dimension(nlayer) :: t            ! temperature (K)
real, dimension(nlayer) :: dens         ! total number density (molecule.cm-3)


pfunc1(:) = param%a*exp(-param%b/t(:))*((t(:)/param%t0)**param%c)*(dens(:)**param%d)

return
end function pfunc1


function pfunc2(nlayer,t,dens,param)

use types_asis

implicit none

! input/output

integer                 :: nlayer       ! number of atmospheric layers
real, dimension(nlayer) :: pfunc2       ! output values
type(rtype2)            :: param        ! parameters
real, dimension(nlayer) :: t            ! temperature (K)
real, dimension(nlayer) :: dens         ! total number density (molecule.cm-3)

! local

real, dimension(nlayer) :: ak0, ak1, ak2, rate, xpo

ak0(:)  = param%k0*exp(-param%a/t(:))*((t(:)/param%t0)**param%n)
ak1(:)  = param%kinf*exp(-param%b/t(:))*((t(:)/param%t0)**param%m)
ak2(:)  = param%g*exp(-param%h/t(:))
rate(:) = (ak0(:)*dens(:)**param%dup)/(1. + ak0(:)*(dens(:)**param%ddown)/ak1(:))
xpo(:)  = 1./(1. + alog10((ak0(:)*dens(:))/ak1(:))**2)

pfunc2(:) = ak2(:) + rate(:)*(param%fc**xpo(:))

return
end function pfunc2

function pfunc3(nlayer,t,dens,param)

use types_asis

implicit none

! input/output

integer                 :: nlayer       ! number of atmospheric layers
real, dimension(nlayer) :: pfunc3       ! output values
type(rtype3)            :: param        ! parameters
real, dimension(nlayer) :: t            ! temperature (K)
real, dimension(nlayer) :: dens         ! total number density (molecule.cm-3)

! local

real, dimension(nlayer) :: ak0, ak1, rate, fc, c, N, d, xpo

ak0(:)  = param%k0*exp(-param%a/t(:))*((t(:)/param%t0)**param%n)
ak1(:)  = param%kinf*exp(-param%b/t(:))*((t(:)/param%t0)**param%m)
rate(:) = (ak0(:)*dens(:)**param%dup)/(1. + ak0(:)*(dens(:)**param%ddown)/ak1(:))
fc(:)   = (1-param%atroe)*exp(-t(:)/param%btroe) + param%atroe*exp(-t(:)/param%ctroe) + exp(-param%dtroe/t(:))
c(:)    = -0.4-0.67*alog10(fc(:))
N(:)    = 0.75-1.27*alog10(fc(:))
d(:)    = 0.14
xpo(:)  = 1./(1. + ( (alog10((ak0(:)*dens(:))/ak1(:))+c(:))/(N(:)-d(:)*(alog10((ak0(:)*dens(:))/ak1(:))+c(:))) )**2)

pfunc3(:) = rate(:)*(fc(:)**xpo(:))

return
end function pfunc3

end module pfunc