  The 2D Discrete Gaussian model gives each height function $\eta :
\mathbb{Z}^2\to\mathbb{Z}$ a probability proportional to $\exp(-\beta
\mathcal{H}(\eta))$, where $\beta$ is the inverse-temperature and
$\mathcal{H}(\eta) = \sum_{x\sim y}(\eta_x-\eta_y)^2$ sums over
nearest-neighbor bonds. We consider the model at large fixed $\beta$, where it
is flat unlike its continuous analog (the Gaussian Free Field).
  We first establish that the maximum height in an $L\times L$ box with 0
boundary conditions concentrates on two integers $M,M+1$ with $M\sim
\sqrt{(1/2\pi\beta)\log L\log\log L}$. The key is a large deviation estimate
for the height at the origin in $\mathbb{Z}^2$, dominated by "harmonic
pinnacles", integer approximations of a harmonic variational problem. Second,
in this model conditioned on $\eta\geq 0$ (a floor), the average height rises,
and in fact the height of almost all sites concentrates on levels $H,H+1$ where
$H\sim M/\sqrt{2}$.
  This in particular pins down the asymptotics, and corrects the order, in
results of Bricmont, El-Mellouki and Fr\"ohlich (1986), where it was argued
that the maximum and the height of the surface above a floor are both of order
$\sqrt{\log L}$.
  Finally, our methods extend to other classical surface models (e.g.,
restricted SOS), featuring connections to $p$-harmonic analysis and alternating
sign matrices.
