We perform molecular dynamics simulations of a coarse-grained model of a polymer-solvent mixture to study solvent evaporation from supported and freestanding polymer films near the bulk glass transition temperature T(g). We find that the evaporation process is characterized by three time (t) regimes: An early regime where the initially large surplus of solvent at the film-vapor interface evaporates and the film thickness h varies little with t, an intermediate regime where h decreases strongly, and a final regime where h slowly converges toward the asymptotic value of the dry film. In the intermediate regime the decrease of h goes along with an increase of the monomer density at the retracting interface. This polymer-rich "crust" is a nonequilibrium effect caused by the fast evaporation rate in our simulation. The interfacial excess of polymer gradually vanishes as the film approaches the dry state. In the intermediate and final time regimes it is possible to describe the simulation data for h(t) and the solvent density profile phi(L)(y,t) by the numerical solution of a one-dimensional diffusion model depending only on the y direction perpendicular to the interface. The key parameter of this model is the mutual diffusion coefficient D(L) of the solvent in the film. Above T(g) we find that a constant D(L) allows to describe the simulation data, whereas near T(g) agreement between simulation and modeling can only be obtained if the diffusion coefficient depends on y through two factors: A factor describing the slowing down of the dynamics with decreasing solvent concentration phi(L)(y,t) and a factor parametrizing the smooth gradient toward enhanced dynamics as the film-vapor interface is approached.