ncore Function

  function ncore(rhopedn, nped, nsep, nav, alphan)

    !! Central density of 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
    !! rhopedn : input real : normalised minor radius pedestal position
    !! nped : input real : pedestal density (/m3)
    !! nsep : input real : separatrix density (/m3)
    !! nav : input real : volume average density (/m3)
    !! alphan : input real : density peaking parameter
    !! This routine calculates the central density
    !! of a pedestalised profile.
    !! J.Johner, Fusion Science and Technology 59 (2011), pp 308-349
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		use error_handling, only: report_error
    implicit none

    real(dp) :: ncore

    !  Arguments

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

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

    ncore = 1.0D0 / (3.0D0*rhopedn**2) * ( &
         3.0D0*nav*(1.0D0 + alphan) &
         + nsep*(1.0D0 + alphan)*(-2.0D0 + rhopedn + rhopedn**2) &
         - nped*( (1.0D0 + alphan)*(1.0D0 + rhopedn) + &
         (alphan - 2.0D0)*rhopedn**2 ) )

    if (ncore < 1.0D-6) then
      ! Prevent ncore from going negative (and terminating the optimisation) by
      ! kludging to small positive value. Allows solver to continue and
      ! hopefully be constrained away from this point (e.g. constraint 81, ne0 > neped)
      ncore = 1.0D-6
      ! write(*,*) 'Error in ncore: negative central density'
      ! write(*,*) 'nped =', nped, ' nsep =', nsep
      ! write(*,*) 'nav =', nav, ' ncore =', ncore
      ! call report_error(212)
    end if

  end function ncore