  Based on the original combination of analytical methods, computer algebra
tools and numerical calculations, proposed recently in Refs. [1]-[3], the
non-perturbative vacuum polarization effects in the 2+1 D supercritical
Dirac-Coulomb system with $Z > Z_{cr,1}$ are explored. Both the vacuum charge
density $\rho_{VP}(\vec{r})$ and vacuum energy $\mathcal{E}_{VP}$ are
considered. The main result of the work is that in the overcritical region
$\mathcal{E}_{VP}$ turns out to be a rapidly decreasing function $\sim -
\eta_{eff}\, Z^3/R\, $ with $\eta_{eff}>0 $ and $R$ being the size of the
external Coulomb source. Due to a lot of details of calculation the whole work
is divided into two parts I and II. In the present part I we consider the
evaluation and behavior of the vacuum density $\rho_{VP}$, which further is
used in the part II for evaluation of the vacuum energy, with emphasis on the
renormalization, convergence of the partial expansion for $\rho_{VP}$ and
behavior of the integral induced charge $Q_{VP}$ in the overcritical region.
