  We construct atmosphere models for strongly magnetized neutron stars with
surface fields $B\sim 10^{12}-10^{15}$ G and effective temperatures $T_{\rm
eff}\sim 10^6-10^7$ K. The atmospheres directly determine the characteristics
of thermal emission from isolated neutron stars, including radio pulsars, soft
gamma-ray repeaters, and anomalous X-ray pulsars. In our models, the atmosphere
is composed of pure hydrogen or helium and is assumed to be fully ionized. The
radiative opacities include free-free absorption and scattering by both
electrons and ions computed for the two photon polarization modes in the
magnetized electron-ion plasma. Since the radiation emerges from deep layers in
the atmosphere with $\rho\ga 10^2$ g/cm$^3$, plasma effects can significantly
modify the photon opacities by changing the properties of the polarization
modes. In the case where the magnetic field and the surface normal are
parallel, we solve the full, angle-dependent, coupled radiative transfer
equations for both polarization modes. We also construct atmosphere models for
general field orientations based on the diffusion approximation of the
transport equations and compare the results with models based on full radiative
transport. In general, the emergent thermal radiation exhibits significant
deviation from blackbody, with harder spectra at high energies. The spectra
also show a broad feature ($\Delta E/\Ebi\sim 1$) around the ion cyclotron
resonance $\Ebi=0.63 (Z/A)(B/10^{14}{G})$ keV, where $Z$ and $A$ are the atomic
charge and atomic mass of the ion, respectively; this feature is particularly
pronounced when $\Ebi\ga 3k\Teff$. Detection of the resonance feature would
provide a direct measurement of the surface magnetic fields on magnetars.
