pthresh Subroutine

  subroutine pthresh(dene,dnla,bt,rmajor,kappa,sarea,aion,aspect,pthrmw)
    !! author: P J Knight, CCFE, Culham Science Centre
    !! L-mode to H-mode power threshold calculation
    !! dene   : input real :  volume-averaged electron density (/m3)
    !! dnla   : input real :  line-averaged electron density (/m3)
    !! bt     : input real :  toroidal field on axis (T)
    !! rmajor : input real :  plasma major radius (m)
    !! kappa  : input real :  plasma elongation
    !! sarea  : input real :  plasma surface area (m**2)
    !! aion   : input real :  average mass of all ions (amu)
    !! aspect : input real :  aspect ratio
    !! pthrmw(17) : output real array : power threshold (different scalings)
    !! This routine calculates the power threshold for the L-mode to
    !! H-mode transition.
    !! ITER Physics Design Description Document, p.2-2
    !! ITER-FDR Plasma Performance Assessments, p.III-9
    !! Snipes, 24th EPS Conference, Berchtesgaden 1997, p.961
    !! Martin et al, 11th IAEA Tech. Meeting on H-mode Physics and
    !! Transport Barriers, Journal of Physics: Conference Series
    !! 123 (2008) 012033
    !! J A Snipes and the International H-mode Threshold Database
    !! Working Group, 2000, Plasma Phys. Control. Fusion, 42, A299
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use physics_variables, only: rminor, plascur
    implicit none

    !  Arguments

    real(dp), intent(in) :: dene,dnla,bt,rmajor,kappa,sarea,aion,aspect
    real(dp), dimension(21), intent(out) :: pthrmw

    !  Local variables

    real(dp) :: dene20,dnla20,marterr

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

    dene20 = 1.0D-20*dene
    dnla20 = 1.0D-20*dnla

    !  ITER-DDD, D.Boucher
    !  Fit to 1996 H-mode power threshold database: nominal

    pthrmw(1) = 0.45D0 * dene20**0.75D0 * bt * rmajor**2

    !  Fit to 1996 H-mode power threshold database: upper bound

    pthrmw(2) = 0.37D0 * dene20 * bt * rmajor**2.5D0

    !  Fit to 1996 H-mode power threshold database: lower bound

    pthrmw(3) = 0.54D0 * dene20**0.5D0 * bt * rmajor**1.5D0

    !  J. A. Snipes, ITER H-mode Threshold Database Working Group,
    !  Controlled Fusion and Plasma Physics, 24th EPS Conference,
    !  Berchtesgaden, June 1997, vol.21A, part III, p.961

    pthrmw(4) = 0.65D0 * dnla20**0.93D0 * bt**0.86D0 * rmajor**2.15D0

    pthrmw(5) = 0.42D0 * dnla20**0.80D0 * bt**0.90D0 * rmajor**1.99D0 &
         * kappa**0.76D0

    !  Martin et al (2008) for recent ITER scaling, with mass correction
    !  and 95% confidence limits

    pthrmw(6) = 0.0488D0 * dnla20**0.717D0 * bt**0.803D0 &
         * sarea**0.941D0 * (2.0D0/aion)

    marterr = 0.057D0**2 + (0.035D0 * log(dnla20))**2 &
         + (0.032D0 * log(bt))**2 + (0.019D0 * log(sarea))**2
    marterr = sqrt(marterr) * pthrmw(6)

    pthrmw(7) = pthrmw(6) + 2.0D0*marterr
    pthrmw(8) = pthrmw(6) - 2.0D0*marterr

    ! Snipes et al (2000) scaling with mass correction
    ! Nominal, upper and lower

    pthrmw(9) = 1.42D0 * dnla20**0.58D0 * bt**0.82D0 * rmajor &
               * rminor**0.81D0 * (2.0D0/aion)

    pthrmw(10) = 1.547D0 * dnla20**0.615D0 * bt**0.851D0 &
              * rmajor**1.089D0 * rminor**0.876D0 * (2.0D0/aion)

    pthrmw(11) = 1.293D0 * dnla20**0.545D0 * bt**0.789D0 &
              * rmajor**0.911D0 * rminor**0.744D0 * (2.0D0/aion)

    ! Snipes et al (2000) scaling (closed divertor) with mass correction
    ! Nominal, upper and lower

    pthrmw(12) = 0.8D0 * dnla20**0.5D0 * bt**0.53D0 * rmajor**1.51D0 &
               * (2.0D0/aion)

    pthrmw(13) = 0.867D0 * dnla20**0.561D0 * bt**0.588D0 * rmajor**1.587D0 &
               * (2.0D0/aion)

    pthrmw(14) = 0.733D0 * dnla20**0.439D0 * bt**0.472D0 * rmajor**1.433D0 &
               * (2.0D0/aion)

    ! Hubbard et al. 2012 L-I threshold scaling

    ! Nominal
    pthrmw(15) = 2.11 * (plascur/1.0D6)**0.94 * dnla20**0.65

    ! Lower bound
    pthrmw(16) = 2.11 * (plascur/1.0D6)**0.70 * dnla20**0.47

    ! Upper bound
    pthrmw(17) = 2.11 * (plascur/1.0D6)**1.18 * dnla20**0.83

    ! Hubbard et al. 2017 L-I threshold scaling
    pthrmw(18) = 0.162 * dnla20 * sarea * (bt)**0.26

    !  Aspect ratio corrected Martin et al (2008)
    !  Correction: Takizuka 2004, Plasma Phys. Control Fusion 46 A227
    if (aspect.le.2.7D0) then
        pthrmw(19) = pthrmw(6) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0))
        pthrmw(20) = pthrmw(7) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0))
        pthrmw(21) = pthrmw(8) * (0.098D0 * aspect / (1.0D0 - (2.0D0/(1.0D0 + aspect))**0.5D0))
    else
        pthrmw(19) = pthrmw(6)
        pthrmw(20) = pthrmw(7)
        pthrmw(21) = pthrmw(8)
    end if

  end subroutine pthresh