A finite-difference time integration method for the calculation of seismic ray intensity is developed. Discontinuities in the depth derivatives of the velocity-depth function at layer boundaries cause anomalies in the intensity distance curves calculated using the standard integral formulation. The time integration method overcomes these difficulties. Calculations for a simple analytic case and a Gutenberg earth model demonstrate the difficulties with the standard integral method and the superior performance of the time integration scheme. The method may also be applied to laterally inhomogeneous earth models.