beamcalc Subroutine

  subroutine beamcalc(nd,nt,ealphadt,ebeam,ecritd,ecritt,tausbme, &
       ftritbm,ibeam,ti,vol,svdt,palfdb,palftb,nhot,ehot)

    !! Neutral beam alpha power and ion energy
    !! author: P J Knight, CCFE, Culham Science Centre
    !! ealphadt : input real :  alpha particle birth energy (D-T) (keV)
    !! ebeam  : input real :  beam energy (keV)
    !! ecritd : input real :  critical energy for electron/ion slowing down of
    !! the beam ion (deuterium neutral beam) (keV)
    !! ecritt : input real :  critical energy for beam slowing down
    !! (tritium neutral beam) (keV)
    !! ftritbm: input real :  beam tritium fraction (0.0 = deuterium beam)
    !! ibeam  : input real :  beam current (A)
    !! nd     : input real :  thermal deuterium density (/m3)
    !! nt     : input real :  thermal tritium density   (/m3)
    !! svdt   : input real :  profile averaged <sigma v> for D-T (m3/s)
    !! tausbme: input real :  beam ion slowing down time on electrons (s)
    !! ti     : input real :  thermal ion temperature (keV)
    !! vol    : input real :  plasma volume (m3) (95% flux surface)
    !! ehot   : output real : average hot beam ion energy (keV)
    !! nhot   : output real : hot beam ion density (/m3)
    !! palfdb : output real : alpha power from deut. beam-background fusion (MW)
    !! palftb : output real : alpha power from trit. beam-background fusion (MW)
    !! This routine calculates the neutral beam alpha power and ion energy.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        use constants, only: echarge, mproton

        implicit none

        ! Arguments
        real(kind(1.0D0)), intent(in) :: ealphadt, ebeam, ecritd, ecritt, &
            ftritbm, ibeam, nd, nt, svdt, tausbme, ti, vol
        real(kind(1.0D0)), intent(out) :: ehot, nhot, palfdb, palftb

        ! Local variables
        integer :: iabm
        real(kind(1.0D0)) :: ebmratd, ebmratt, ehotd, ehott, ifbmd, ifbmt, &
            ndhot, nhotmsd, nhotmst, nthot, presd, prest, s0d, s0t, svdhotn, &
            svthotn, tauseffd, tausefft, vcds, vcritd, vcritt, vcts, xcoefd, &
            xcoeft
        real(kind(1.0D0)) :: atmd, atmt, epsabs, epsrel

        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ! Initialise shared variables
        atmd = 2.0D0   !  atomic mass of deuterium
        atmt = 3.0D0   !  atomic mass of tritium
        epsabs = 1.0D-7  !  absolute error
        epsrel = 1.0D-7  !  relative error

        ! D and T beam current fractions
        ifbmd = ibeam * (1.0D0 - ftritbm)
        ifbmt = ibeam * ftritbm

        ebmratd = ebeam/ecritd
        vcritd = sqrt(2.0D0*echarge*1000.0D0*ecritd/(mproton*atmd))
        tauseffd = tausbme/3.0D0 * log(1.0D0+(ebmratd)**1.5D0)
        nhotmsd = (1.0D0-ftritbm) * ibeam * tauseffd/(echarge * vol)

        ebmratt = ebeam/ecritt
        vcritt = sqrt(2.0D0*echarge*1000.0D0*ecritt/(mproton*atmt))
        tausefft = tausbme/3.0D0 * log(1.0D0+(ebmratt)**1.5D0)
        nhotmst = ftritbm * ibeam * tausefft/(echarge * vol)

        nhot = nhotmsd + nhotmst
        ndhot = nhotmsd
        nthot = nhotmst

        ! Average hot ion energy from Deng & Emmert, UWFDM-718, Jan 87
        vcds = 2.0D0 * ecritd * echarge * 1000.0D0/(2.0D0 * mproton)
        vcts = 2.0D0 * ecritt * echarge * 1000.0D0/(3.0D0 * mproton)

        s0d = ifbmd/(echarge * vol)
        s0t = ifbmt/(echarge * vol)

        xcoefd = atmd * mproton * tausbme * vcds * s0d / &
            (echarge * 1000.0D0 * 3.0D0)
        xcoeft = atmt * mproton * tausbme * vcts * s0t / &
            (echarge * 1000.0D0 * 3.0D0)

        presd = xcoefd * xbrak(ebeam,ecritd)
        prest = xcoeft * xbrak(ebeam,ecritt)

        ehotd = 1.5D0 * presd/ndhot
        ehott = 1.5D0 * prest/nthot
        ehot = (ndhot*ehotd + nthot*ehott)/nhot

        iabm = 2 ; svdhotn = 1.0D-4 * sgvhot(iabm,vcritd,ebeam)
        iabm = 3 ; svthotn = 1.0D-4 * sgvhot(iabm,vcritt,ebeam)

        palfdb = palphabm(ealphadt,ndhot,nt,svdhotn,vol,ti,svdt)
        palftb = palphabm(ealphadt,nthot,nd,svthotn,vol,ti,svdt)

    contains

        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        function xbrak(e0,ec)
            !! Hot ion energy parameter
            !! author: P J Knight, CCFE, Culham Science Centre
            !! e0 : input real :  neutral beam energy (keV)
            !! ec : input real :  critical energy for electron/ion slowing down of
            !! the beam ion (keV)
            !! This routine calculates something to do with the hot ion energy...
            !!             !
            ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            implicit none

            real(dp) :: xbrak

            ! Arguments
            real(dp), intent(in) :: e0, ec

      real(dp) :: ans,t1,t2,t3,t4,xarg,xc,xcs

            ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            xcs = e0/ec
            xc = sqrt(xcs)

            t1 = xcs/2.0D0
            t2 = (log((xcs + 2.0D0*xc + 1.0D0)/(xcs - xc + 1.0D0)))/6.0D0

            xarg = (2.0D0*xc -1.0D0)/sqrt(3.0D0)
            t3 = (atan(xarg))/sqrt(3.0D0)
            t4 = 0.3022999D0

            ans = t1 + t2 - t3 - t4
            xbrak = ans

        end function xbrak

        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function palphabm(ealphadt,nbm,nblk,sigv,vol,ti,svdt)

      !! Alpha power from beam-background fusion
      !! author: P J Knight, CCFE, Culham Science Centre
      !! ealphadt : input real :  alpha particle birth energy (D-T) (keV)
      !! nblk   : input real :  thermal ion density (/m3)
      !! nbm    : input real :  hot beam ion density (/m3)
      !! sigv   : input real :  hot beam fusion reaction rate (m3/s)
      !! svdt   : input real :  profile averaged <sigma v> for D-T (m3/s)
      !! ti     : input real :  thermal ion temperature (keV)
      !! vol    : input real :  plasma volume (m3)
      !! This routine calculates the alpha power from
      !! beam-background fusion.
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: palphabm

      !  Arguments

      real(dp), intent(in) :: ealphadt,nblk,nbm,sigv,svdt,ti,vol

      !  Local variables

      real(dp) :: ratio

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      ratio = svdt / bosch_hale(ti,1)

      palphabm = echarge/1000.0D0 * nbm * nblk * sigv * ealphadt * vol * ratio

    end function palphabm

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function sgvhot(iabm,vcrx,ebeam)

      !! Hot beam fusion reaction rate
      !! author: P J Knight, CCFE, Culham Science Centre
      !! ebeam  : input real :  neutral beam energy (keV)
      !! iabm   : input integer : switch denoting type of ion (2=D,3=T)
      !! vcrx   : input real :  critical velocity for electron/ion slowing down of
      !! the beam ion (m/s)
      !! This routine calculates the hot beam fusion reaction rate in m3/s.
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use error_handling, only: idiags, report_error
      use maths_library, only: quanc8

      implicit none

      real(dp) :: sgvhot

      !  Arguments
      integer, intent(in) :: iabm
      real(dp), intent(in) :: ebeam, vcrx

      !  Local variables

      integer :: nofun
      real(dp) :: abm,abserr,epsabs1,flag,svint,t1,t2, &
           vbeam,vbeams,xv

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      epsabs1 = 1.0D-33

      if (iabm == 2) then
         abm = atmd
      else if (iabm == 3) then
         abm = atmt
      else
         idiags(1) = iabm ; call report_error(84)
      end if

      !  Initialise global variables

      vcritx = vcrx

      !  Beam velocity

      vbeams = ebeam * echarge * 1000.0D0 * 2.0D0/(abm * mproton)
      vbeam = sqrt(vbeams)

      xv = vbeam/vcrx
      t1 = 3.0D0 * vcrx/log(1.0D0+(xv**3))

      call quanc8(fsv,0.0D0,xv,epsabs1,epsrel,svint,abserr,nofun,flag)
      t2 = svint

      sgvhot = t1 * t2

    end function sgvhot

    end subroutine beamcalc