Io's sublimation-driven atmosphere is modeled using the direct simulation Monte Carlo (DSMC) method. These rarefied gas dynamics simulations improve upon earlier models by using a three-dimensional domain encompassing the entire planet computed in parallel. The effects of plasma heating, planetary rotation, inhomogeneous surface frost, molecular residence time of SO 2 on the exposed (non-volatile) rocky surface, and surface temperature distribution are investigated. Circumplanetary flow is predicted to develop from the warm dayside toward the cooler nightside. Io's rotation leads to a highly asymmetric frost surface temperature distribution (due to the frost's high thermal inertia) which results in circumplanetary flow that is not axi-symmetric about the subsolar point. The non-equilibrium thermal structure of the atmosphere, specifically vibrational and rotational temperatures, is also examined. Plasma heating is found to significantly inflate the atmosphere on both the dayside and nightside. The plasma energy flux causes high temperatures at high altitudes but plasma energy depletion through the dense gas column above the warmest frost permits gas temperatures cooler than the surface at low altitudes. A frost map (Douté, S., Schmitt, B., Lopes-Gautier, R., Carlson, R., Soderblom, L., Shirley, J., and the Galileo NIMS Team . Icarus 149, 107-132) is used to control the sublimated flux of SO 2 which can result in inhomogeneous column densities that vary by nearly a factor of four for the same surface temperature. A short residence time for SO 2 molecules on the "rock" component is found to smooth lateral atmospheric inhomogeneities caused by variations in the surface frost distribution, creating an atmosphere that looks nearly identical to one with uniform frost coverage. A longer residence time is found to agree better with mid-infrared observations (Spencer, J.R., Lellouch, E., Richter, M.J., López-Valverde, M.A., Jessup, K.L, Greathouse, T.K., Flaud, J. . Icarus 176, 283-304) and reproduce the observed anti-jovian/sub-jovian column density asymmetry. The computed peak dayside column density for Io assuming a surface frost temperature of 115 K agrees with those suggested by Lyman-α observations (Feaga, L.M., McGrath, M., Feldman, P.D. . Icarus 201, 570-584). On the other hand, the peak dayside column density at 120 K is a factor of five larger and is higher than the upper range of observations (Jessup, K.L., Spencer, J.R., Ballester, G.E., Howell, R.R., Roesler, F., Vigel, M., Yelle, R. . Icarus 169, 197-215; Spencer et al., 2005).