The ground-state properties of superfluid nuclear systems with 1S0 pairing are studied within a local energy-density functional (LEDF) approach. A new form of the LEDF is proposed with a volume part which fits the Friedman-Pandharipande and Wiringa-Fiks-Fabrocini equation of state at low and moderate densities and allows an extrapolation to higher densities which preserves causality. For inhomogeneous systems, a surface term is added, with two free parameters, which has a fractional form like a Padé approximant containing the square of the density gradient in both the numerator and denominator. In addition to the direct and exchange Coulomb interaction energy, an effective density-dependent Coulomb-nuclear correlation term is included with one more free parameter. A three-parameter fit to the masses and radii of about 100 spherical nuclei has shown that the latter term gives a contribution of the same order of magnitude as the Nolen-Schiffer anomaly in the Coulomb displacement energy. The root-mean-square deviations from experimental masses and radii with the proposed LEDF come out about a factor of two smaller than those obtained with the conventional functionals based on the Skyrme or finite-range Gogny force, or on relativistic mean-field theory. The generalized variational principle is formulated leading to the self-consistent Gor'kov equations which are sovled exactly, with physical boundary conditions both for the bound and scattering states. The method is used to calculate the differential observables such as odd-even mass differences and staggering in charge radii. With a zero-range density-dependent cutoff pairing interaction incorporating a density-gradient term, the evolution of these observables is reproduced reasonably well, including the kinks at magic neutron numbers and the sizes of the associated staggering. An extrapolation from the pairing properties of finite nuclei to pairing in infinite nuclear matter is discussed. A "reference" value of the pairing gap ΔF≈ 3.3 MeV is found for subsaturated nuclear matter at about 0.65 of the equilibrium density. With the formulated LEDF approach, we study also the dilute limit in both the weak and strong coupling regimes. Within the sum rules approach it is shown that the density-dependent pairing may also induce sizeable staggering and kinks in the evolution of the mean energies of multipole excitations.