nprofile Function

  function nprofile(rho, rhopedn, n0, nped, nsep, alphan)

    !! Implementation of HELIOS-type density pedestal profile
    !! author: R Kemp, CCFE, Culham Science Centre
    !! author: H Lux, CCFE, Culham Science Centre
    !! author: P J Knight, CCFE, Culham Science Centre
    !! rho     : input real : normalised minor radius
    !! rhopedn : input real : normalised minor radius pedestal position
    !! n0      : input real : central density (/m3)
    !! nped    : input real : pedestal density (/m3)
    !! nsep    : input real : separatrix density (/m3)
    !! alphan  : input real : density peaking parameter
    !! This routine calculates the density at a normalised minor
    !! radius position rho for a pedestalised profile.
    !! <P>If <CODE>ipedestal = 0</CODE> the original parabolic
    !! profile form is used instead.
    !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		use error_handling, only: fdiags, report_error
		use physics_variables, only: ipedestal
    implicit none

    real(dp) :: nprofile

    !  Arguments

    real(dp), intent(in) :: rho, rhopedn, n0,  nped, nsep, alphan

    !  Local variables

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

    if (ipedestal == 0) then
       nprofile = n0 * (1.0D0 - rho**2)**alphan
       return
    end if

    !  Error trap; shouldn't happen unless volume-averaged density has
    !  been allowed to drop below nped. This may happen during a HYBRD case,
    !  but should have been prevented for optimisation runs.

    !  Input checks

    if (n0 < nped) then
       fdiags(1) = nped ; fdiags(2) = n0
       call report_error(153)
    end if

    if (rho <= rhopedn) then
       nprofile = nped + (n0 - nped) * (1.0D0 - (rho/rhopedn)**2)**alphan
    else
       nprofile = nsep + (nped - nsep) * (1.0D0 - rho)/(1.0D0 - rhopedn)
    end if

  end function nprofile