module m_chemistry

use m_parameters

implicit none

real, private :: coef(10)
real, private :: lambda_onestep
real, private :: lambda1_twostep
real, private :: lambda2_twostep
real, private :: beta1_twostep
real, private :: hrp_twostep

contains

subroutine init_chemistry

    character(len=40) :: nrxn_string

    if (nrxn == 1) then
       reaction_type = "onestep"
    else if (nrxn == 2) then
       reaction_type = "twostep"
    else
       write(nrxn_string, *) nrxn
       reaction_type = trim(nrxn_string) // "-step"
    end if

    if ( reaction_type == "onestep" ) then
        lambda_onestep = pre * exp ( - beta / hrp )
    else if ( reaction_type == "twostep" ) then
        lambda1_twostep = lambda1
        lambda2_twostep = lambda2
        beta1_twostep = beta1
        hrp_twostep = hrp
    else
        WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
        stop
    end if

end subroutine init_chemistry


subroutine update_chemistry (t)

    real :: t

    real :: factor

    real :: relax_duration = 60.

    if ( reaction_type == "onestep" ) then
        lambda_onestep = pre * exp ( - beta / hrp )
    else if ( reaction_type == "twostep" ) then
        if (t < relax_duration) then
            factor = (relax_duration + t) / relax_duration / 2.
        else
            factor = 1.
        end if

        lambda1_twostep = factor * lambda1
        lambda2_twostep = factor * lambda2

    else
        stop
    end if

end subroutine update_chemistry


real function rate_1step (yr, theta)

    real, intent(in)  :: yr
    real, intent(in)  :: theta

    real :: y
    real :: t_reduce

    y=yr
    ! if(yr.lt.0.) y=0.
    ! if(yr.gt.1.) y=1.

    t_reduce=theta
    ! if(theta.lt.0.) t_reduce=0.
    ! if(theta.gt.1.) t_reduce=1.

    if (t_reduce.gt.c_ref) then

      rate_1step = pre*y*exp(-ac/(1.+bc*t_reduce))

    else if (t_reduce.le.c_cut) then

      rate_1step = min_wr

    else

      rate_1step = &
        ((refwr-min_wr)*exp(prof_wr*(t_reduce-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
        / (1.-exp(prof_wr*(c_cut-c_ref)))

   endif

end function rate_1step


real function rate1_2step (ya, yx, theta)

    real, intent(in)  :: ya
    real, intent(in)  :: yx
    real, intent(in)  :: theta

    real :: y1
    real :: y2
    real :: t_reduce

    y1=ya
    if(ya.lt.0.) y1=0.
    if(ya.gt.1.) y1=1.

    y2=yx
    if(yx.lt.0.) y2=0.
    if(yx.gt.1.) y2=1.

    t_reduce=theta
    if(theta.lt.0.) t_reduce=0.
    if(theta.gt.1.) t_reduce=1.

    rate1_2step = lambda1_twostep * y1 * y2 * & 
                  exp (-(beta1_twostep*(1. - t_reduce))/(1. - hrp_twostep*(1. - t_reduce)))

end function rate1_2step


real function rate2_2step (yx, theta)

    real, intent(in)  :: yx
    real, intent(in)  :: theta

    real :: y
    real :: t_reduce

    y=yx
    if(yx.lt.0.) y=0.
    if(yx.gt.1.) y=1.

    t_reduce=theta
    if(theta.lt.0.) t_reduce=0.
    if(theta.gt.1.) t_reduce=1.

    rate2_2step = lambda2_twostep * yx * yx

end function rate2_2step

end module m_chemistry