tcore Function

  function tcore(rhopedt, tped, tsep, tav, alphat, tbeta)

    !! Central temperature for a pedestal profile
    !! author: R Kemp, CCFE, Culham Science Centre
    !! author: H Lux, CCFE, Culham Science Centre
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rhopedt : input real : normalised minor radius pedestal position
    !! tped : input real : pedestal temperature (keV)
    !! tsep : input real : separatrix temperature (keV)
    !! tav : input real : volume average temperature (keV)
    !! alphat : input real : temperature peaking parameter
    !! tbeta : input real : second temperature exponent
    !! This routine calculates the core temperature (keV)
    !! of a pedestalised profile.
    !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		use constants, only: pi
		use maths_library, only: gamfun
    implicit none

    real(dp) :: tcore

    !  Arguments

    real(dp), intent(in) :: rhopedt, tped, tsep, tav, alphat, tbeta

    !  Local variables

    real(dp), parameter :: numacc = 1.0D-7
    real(dp) :: gamfac

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

    !  For integer values of alphat, the limit of
    !  gamfun(-alphat)*sin(pi*alphat) needs to be calculated directly

    gamfac = gamfun(1.0D0 + alphat + 2.0D0/tbeta) / &
         gamfun((2.0D0 + tbeta)/tbeta) / rhopedT**2

    if (abs(alphat-nint(alphat)) <= numacc) then
       gamfac = -gamfac / gamfun(1.0D0 + alphat)
    else
       gamfac = gamfac * gamfun(-alphat)*sin(pi*alphat)/pi
    end if

    !  Calculate core temperature

    tcore = tped + gamfac * ( tped*rhopedT**2 - tav + &
         (1.0D0 - rhopedT)/3.0D0 * &
         ( (1.0D0 + 2.0D0*rhopedT)*tped + (2.0D0 + rhopedT)*tsep ) )

  end function tcore