We construct radiative equilibrium models for strongly magnetized (B greater than or similar to 10(13) G) neutron star atmospheres taking into account magnetic free-free absorption and scattering processes computed for two polarization modes. We include the effects of vacuum polarization in our calculations. We present temperature profiles and the angle-, photon energy-, and polarization-dependent emerging intensity for a range of magnetic field strengths and effective temperatures of the atmospheres. We find that for B less than or similar to 10(14) G, the emerging spectra are bluer than the blackbody corresponding to the effective temperature T-eff with modified Planckian shapes because of the photon energy dependence of the magnetic opacities. However, vacuum polarization resonance significantly modifies the spectra for B similar to 10(15) G, giving rise to power-law tails at high photon energies. The angle dependence (beaming) of the emerging intensity has two maxima: a narrow (pencil) peak at small angles (less than or similar to 5 degrees) with respect to the normal and a broad maximum (fan beam) at intermediate angles (similar to 20 degrees- 60 degrees). The relative importance and the opening angle of the radial beam decreases strongly with increasing magnetic field strength and decreasing photon energy. We finally compute a T-eff -T-c relation for our models, where T-c is the local color temperature of the spectrum emerging from the neutron star surface, and find that T-c/T-eff ranges between 1.1 and 1.8. We discuss the implications of our results for various thermally emitting neutron star models.