  A density-functional approach is used to calculate the inhomogeneous vortex
density distribution in the flux liquid phase at the planar surface of a
layered superconductor, where the external magnetic field is perpendicular to
the superconducting layers and parallel to the surface. The interactions with
image vortices are treated within a mean field approximation as a functional of
the vortex density. Near the freezing transition strong vortex density
fluctuations are found to persist far into the bulk liquid. We also calculate
the height of the Bean-Livingston surface barrier.
