We present a Monte Carlo model to compute cometary dust tail brightnesses. The model, anologous in many ways to that developed by Fulle (1989, Astron. Astrophys. 217, 283), computes the keplerian orbit of each grain released from the nucleus surface in the Sun gravity field reduced by the radiation pressure. It can be used to derive the time dependence of the particle size distribution functions and the terminal ejection velocities as a function of size and time. The ejection anisotropy, which is an important parameter regarding the nucleus surface structure, can also be estimated from the model. The Monte Carlo procedure allows us also to compute the total dust grain mass which will remain in bound orbits, contributing to the zodiacal dust cloud. We will present some applications of the model to observed cometary dust tails.