palph Subroutine

  subroutine palph(alphan,alphat,deni,fdeut,fhe3,ftrit,ti, &
       palppv,pchargepv,pneutpv,sigvdt,fusionrate,alpharate,protonrate, &
       pdtpv,pdhe3pv,pddpv)

    !! (Initial part of) fusion power and fast alpha pressure calculations
    !! author: P J Knight, CCFE, Culham Science Centre
    !! alphan     : input real :  density profile index
    !! alphat     : input real :  temperature profile index
    !! deni       : input real :  fuel ion density (/m3)
    !! fdeut      : input real :  deuterium fuel fraction
    !! fhe3       : input real :  helium-3 fuel fraction
    !! ftrit      : input real :  tritium fuel fraction
    !! ti         : input real :  ion temperature (keV)
    !! palppv     : output real : alpha particle fusion power per volume (MW/m3)
    !! pchargepv  : output real : other charged particle fusion power/volume (MW/m3)
    !! pneutpv    : output real : neutron fusion power per volume (MW/m3)
    !! sigvdt     : output real : profile averaged <sigma v DT> (m3/s)
    !! fusionrate : output real : fusion reaction rate (reactions/m3/s)
    !! alpharate  : output real : alpha particle production rate (/m3/s)
    !! protonrate : output real : proton production rate (/m3/s)
    !! pdtpv      : output real : D-T fusion power (MW/m3)
    !! pdhe3pv    : output real : D-He3 fusion power (MW/m3)
    !! pddpv      : output real : D-D fusion power (MW/m3)
    !! This subroutine numerically integrates over plasma cross-section to
    !! find the core plasma fusion power.
    !! T&amp;M/PKNIGHT/LOGBOOK24, p.6
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: echarge
    use maths_library, only: quanc8
    implicit none

    !  Arguments

    real(dp), intent(in) :: alphan, alphat, deni, fdeut, &
         fhe3, ftrit, ti
    real(dp), intent(out) :: palppv, pchargepv, pneutpv, sigvdt, &
         fusionrate, alpharate, protonrate, pdtpv, pdhe3pv, pddpv

    !  Local variables

    integer, parameter :: DT=1, DHE3=2, DD1=3, DD2=4
    integer :: ireaction,nofun
    real(dp) :: alow,arate,bhigh,epsq8,errest,etot,flag, &
         fpow,frate,pa,pc,pn,prate,sigmav

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Initialise local quantities
         alow = 0.0D0
         bhigh = 1.0D0
         epsq8 = 1.0D-9

         ! Find fusion power
         ! Integrate over plasma profiles to obtain fusion reaction rate
         palppv = 0.0D0
         pchargepv = 0.0D0
         pneutpv = 0.0D0
         fusionrate = 0.0D0
         alpharate = 0.0D0
         protonrate = 0.0D0
         pddpv = 0.0D0

         do ireaction = 1,4
             ! Fusion reaction rate (m3/s) is calculated in fint for each ireaction
             ! sigmav is the volume-averaged fusion reaction rate (m3/s)
             ! = integral(2 rho sigv(rho).ni(rho)^2 drho) / (deni**2)

             call quanc8(fint,alow,bhigh,epsq8,epsq8,sigmav,errest,nofun,flag)
             if (ireaction == DT) sigvdt = sigmav

             select case (ireaction)

                 case (DT)  ! D + T --> 4He + n reaction

                     etot = 17.59D0 * echarge  ! MJ
                     fpow = 1.0D0 * sigmav * etot * fdeut*ftrit * deni*deni  ! MW/m3
                     pa = 0.2D0 * fpow
                     pc = 0.0D0
                     pn = 0.8D0 * fpow
                     frate = fpow/etot  ! reactions/m3/second
                     arate = frate
                     prate = 0.0D0
                     pdtpv = fpow

                 case (DHE3)  ! D + 3He --> 4He + p reaction

                     etot = 18.35D0 * echarge  ! MJ
                     fpow = 1.0D0 * sigmav * etot * fdeut*fhe3 * deni*deni  ! MW/m3
                     pa = 0.2D0 * fpow
                     pc = 0.8D0 * fpow
                     pn = 0.0D0
                     frate = fpow/etot  ! reactions/m3/second
                     arate = frate
                     prate = frate      ! proton production /m3/second
                     pdhe3pv = fpow

                 case (DD1)  ! D + D --> 3He + n reaction

                     ! The 0.5 branching ratio is assumed to be included in sigmav
                     etot = 3.27D0 * echarge  ! MJ
                     fpow = 1.0D0 * sigmav * etot * 0.5D0*fdeut*fdeut * deni*deni  ! MW/m3
                     pa = 0.0D0
                     pc = 0.25D0 * fpow
                     pn = 0.75D0 * fpow
                     frate = fpow/etot  ! reactions/m3/second
                     arate = 0.0D0
                     prate = 0.0D0      ! Issue #557: No proton production
                     pddpv = pddpv + fpow

                 case (DD2)  !  D + D --> T + p reaction

                     ! The 0.5 branching ratio is assumed to be included in sigmav
                     etot = 4.03D0 * echarge  ! MJ
                     fpow = 1.0D0 * sigmav * etot * 0.5D0*fdeut*fdeut * deni*deni  ! MW/m3
                     pa = 0.0D0
                     pc = fpow
                     pn = 0.0D0
                     frate = fpow/etot  ! reactions/m3/second
                     arate = 0.0D0
                     prate = frate      ! proton production /m3/second
                     pddpv = pddpv + fpow

             end select

             palppv = palppv + pa
             pchargepv = pchargepv + pc
             pneutpv = pneutpv + pn
             fusionrate = fusionrate + frate
             alpharate = alpharate + arate
             protonrate = protonrate + prate

         end do

     contains

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

    function fint(rho)

      !! Integrand for fusion power integration
      !! author: P J Knight, CCFE, Culham Science Centre
      !! rho : input real :  Abscissa of the integration, = normalised
      !! plasma minor radius (0.0 <= rho < 1.0)
      !! This function evaluates the integrand for the fusion power
      !! integration, performed using routine
      !! <A HREF="quanc8.html">QUANC8</A>
      !! in routine <A HREF="palph.html">PALPH</A>.
      !! The fusion reaction assumed is controlled by flag
      !! <CODE>ireaction</CODE> set in <CODE>PALPH</CODE>.
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      use physics_variables, only: te, rhopedt, te0, teped, tesep, tbeta, &
        dene, rhopedn, ne0, neped, nesep
      use profiles_module, only: tprofile, nprofile
      implicit none

      real(dp) :: fint

      !  Arguments

      real(dp), intent(in) :: rho

      !  Local variables

      real(dp) :: nprof, nprofsq, sigv, tiofr

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

      !  Local ion temperature (keV) at r/a = rho

      tiofr = ti/te * tprofile(rho,rhopedt,te0,teped,tesep,alphat,tbeta)

      !  Fusion reaction rate (m3/s)

      sigv = bosch_hale(tiofr,ireaction)

      !  Integrand for the volume averaged fusion reaction rate sigmav:
      !  sigmav = integral(2 rho (sigv(rho) ni(rho)^2) drho),
      !  divided by the square of the volume-averaged ion density
      !  to retain the dimensions m3/s (this is multiplied back in later)

      nprof = 1.0D0/dene * nprofile(rho,rhopedn,ne0,neped,nesep,alphan)
      nprofsq = nprof*nprof

      fint = 2.0D0 * rho * sigv * nprofsq

    end function fint

  end subroutine palph