palph2 Subroutine

  subroutine palph2(bt,bp,dene,deni,dnitot,falpe,falpi,palpnb, &
       ifalphap,pchargepv,pneutpv,ten,tin,vol,palpmw,pneutmw,pchargemw, &
       betaft,palppv,palpipv,palpepv,pfuscmw,powfmw)

    !! (Concluding part of) fusion power and fast alpha pressure
    !! calculations
    !! author: P J Knight, CCFE, Culham Science Centre
    !! bp       : input real :  poloidal field (T)
    !! bt       : input real :  toroidal field on axis (T)
    !! dene     : input real :  electron density (/m3)
    !! deni     : input real :  fuel ion density (/m3)
    !! dnitot   : input real :  total ion density (/m3)
    !! falpe    : input real :  fraction of alpha energy to electrons
    !! falpi    : input real :  fraction of alpha energy to ions
    !! ifalphap : input integer :  switch for fast alpha pressure method
    !! palpnb   : input real :  alpha power from hot neutral beam ions (MW)
    !! pchargepv : input real : other charged particle fusion power/volume (MW/m3)
    !! pneutpv  : input/output real : neutron fusion power per volume (MW/m3)
    !! ten      : input real :  density-weighted electron temperature (keV)
    !! tin      : input real :  density-weighted ion temperature (keV)
    !! vol      : input real :  plasma volume (m3)
    !! palpmw   : output real : alpha power (MW)
    !! pneutmw  : output real : neutron fusion power (MW)
    !! pchargemw : output real : other charged particle fusion power (MW)
    !! betaft   : output real : fast alpha beta component
    !! palppv   : input/output real : alpha power per volume (MW/m3)
    !! palpepv  : output real : alpha power per volume to electrons (MW/m3)
    !! palpipv  : output real : alpha power per volume to ions (MW/m3)
    !! pfuscmw  : output real : charged particle fusion power (MW)
    !! powfmw   : output real : fusion power (MW)
    !! This subroutine completes the calculation of the fusion power
    !! fast alpha pressure, and determines other alpha particle quantities.
    !! ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al,
    !! ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
    !! D J Ward, UKAEA Fusion: F/PL/PJK/PROCESS/CODE/050
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: echarge, rmu0
    use physics_variables, only: falpha, fdeut
    implicit none

    !  Arguments

    integer, intent(in) :: ifalphap
    real(dp), intent(in) :: bp, bt, dene, deni, dnitot, falpe, &
         falpi, palpnb, pchargepv, ten, tin, vol
    real(dp), intent(inout) :: palppv, pneutpv
    real(dp), intent(out) :: palpmw, pneutmw, pchargemw, betaft, palpepv, &
         palpipv, pfuscmw, powfmw

    !  Local variables

    real(dp) :: betath, fact, fact2, palppv_no_nb

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

    ! Store the palppv value without the NB alpha power
    palppv_no_nb = palppv

    !  Add neutral beam alpha power / volume
    palppv = palppv + palpnb/vol

    !  Add extra neutron power
    pneutpv = pneutpv + 4.0D0*palpnb/vol

    !  Total alpha power
    palpmw = palppv*vol

    !  Total non-alpha charged particle power
    pchargemw = pchargepv*vol

    !  Total neutron power
    pneutmw = pneutpv*vol

    !  Total fusion power
    powfmw = palpmw + pneutmw + pchargemw

    !  Charged particle fusion power
    pfuscmw = palpmw + pchargemw

    !  Alpha power to electrons and ions (used with electron
    !  and ion power balance equations only)
    !  No consideration of pchargepv here...
    palpipv = falpha * palppv*falpi
    palpepv = falpha * palppv*falpe

    !  Determine average fast alpha density
    if (fdeut < 1.0D0) then

       betath = 2.0D3*rmu0*echarge * (dene*ten + dnitot*tin)/(bt**2 + bp**2)

       ! jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T)
       ! IPDG89 fast alpha scaling
       if (ifalphap == 0) then
          fact = min( 0.30D0, &
               0.29D0*(deni/dene)**2 * ( (ten+tin)/20.0D0 - 0.37D0) )

       ! Modified scaling, D J Ward
       else
          fact = min( 0.30D0, &
               0.26D0*(deni/dene)**2 * &
               sqrt( max(0.0D0, ((ten+tin)/20.0D0 - 0.65D0)) ) )
       end if

       fact = max(fact,0.0D0)
       fact2 = palppv / palppv_no_nb
       betaft = betath * fact*fact2

    else  !  negligible alpha production, palppv = palpnb = 0
       betaft = 0.0D0
    end if

  end subroutine palph2