get_neutron_rate Subroutine

subroutine get_neutron_rate(plasma, eb, rate)
    !+ Gets neutron 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 :: err_status, neb, nt, ebi, tii
    real(Float64) :: dlogE, dlogT, logEmin, logTmin
    real(Float64) :: logeb, logti, lograte, denp
    type(InterpolCoeffs2D) :: c
    real(Float64) :: b11, b12, b21, b22

    real(Float64) :: mullog10 = log(10.0d0)

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

    !!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_NEUTRON_RATE: Eb or Ti out of range of D_D table. Setting D_D rates to zero"
            write(*,'("eb = ",ES10.3," [keV]")') eb
            write(*,'("ti = ",ES10.3," [keV]")') plasma%ti
        endif
        denp = 0.d0
    endif

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

    if (lograte.lt.tables%D_D%minlog_rate) then
        rate = 0.d0
    else
        rate = denp * exp(lograte*mullog10)
    endif

end subroutine get_neutron_rate