get_dd_rate Subroutine

subroutine get_dd_rate(plasma, eb, rate, branch)
    !+ Gets d(d,p)T rate for a beam with energy `eb` interacting with a target plasma
    type(LocalProfiles), intent(in) :: plasma
        !+ Plasma Paramters
    real(Float64), intent(in)       :: eb
        !+ Beam energy [keV]
    real(Float64), intent(out)      :: rate
        !+ Neutron reaction rate [1/s]
    integer, intent(in), optional   :: branch
        !+ Indicates 1 for cfpd rate and 2 for neutron rate

    integer :: err_status, neb, nt, ebi, tii, is, ib
    real(Float64) :: dlogE, dlogT, logEmin, logTmin
    real(Float64) :: logeb, logti, lograte
    type(InterpolCoeffs2D) :: c
    real(Float64) :: b11, b12, b21, b22

    if(present(branch)) then
        ib = branch
    else
        ib = 2
    endif

    logeb = log10(eb)
    logti = log10(plasma%ti)

    !!D_D
    err_status = 1
    logEmin = tables%D_D%logemin
    logTmin = tables%D_D%logtmin
    dlogE = tables%D_D%dlogE
    dlogT = tables%D_D%dlogT
    neb = tables%D_D%nenergy
    nt = tables%D_D%ntemp
    call interpol_coeff(logEmin, dlogE, neb, logTmin, dlogT, nt, &
                        logeb, logti, c, err_status)
    ebi = c%i
    tii = c%j
    b11 = c%b11
    b12 = c%b12
    b21 = c%b21
    b22 = c%b22
    if(err_status.eq.1) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "GET_DD_RATE: Eb or Ti out of range of D_D table. Setting D_D rates to zero"
        endif
        rate = 0.d0
        return
    endif

    lograte = (b11*tables%D_D%log_rate(ebi,tii,ib)   + &
               b12*tables%D_D%log_rate(ebi,tii+1,ib) + &
               b21*tables%D_D%log_rate(ebi+1,tii,ib) + &
               b22*tables%D_D%log_rate(ebi+1,tii+1,ib))

    if (lograte.lt.tables%D_D%minlog_rate) then
        rate = 0.d0
    else
        rate = 0.d0
        do is=1,n_thermal
            if(thermal_mass(is).eq.H2_amu) then
                rate = rate + plasma%deni(is) * exp(lograte*log_10)
            endif
        enddo
    endif

end subroutine get_dd_rate