In this work, we investigate the effect of the introduction of a stress diffusive term into the classical Oldroyd-B constitutive equation on the numerical stability of time-dependent viscoelastic flow calculations. The channel Poiseuille flow at Re much greater than 1 and O(1) We is chosen as a test problem. Through a linear stability analysis, we demonstrate that the introduction of a small amount of (dimensionless) diffusivity, typically of the order 10(-3), does not affect the critical eigenmodes of the viscoelastic Orr-Sommerfeld problem appreciably. However, a diffusive term of that magnitude is shown to have a significant influence on the singular eigenmodes of the classical Oldroyd-B model, associated with the continuum spectra. A finite amplitude perturbation is constructed as a linear superposition of the eigenvectors corresponding to the most unstable eigenvalues of the problem. This is superimposed on the steady Poiseuille flow solution to provide the initial conditions for time-dependent simulations. The numerical algorithm involves a fully spectral spatial discretization and a semi-implicit second order integration in time. For the Oldroyd-B fluid, depending on the magnitude of the initial perturbation, numerical instabilities set in at relatively short times while the components of the conformation tensor increase monotonically in magnitude with time. Introduction of a diffusive term into this model is shown to stabilize the calculations remarkably, and for a three-dimensional simulation with Re = 5000 and We = 1, no instabilities were observed even at very large times. The effect of the magnitude of the diffusivity on the stability and the flow dynamics is addressed through a direct comparison of the results with those obtained for the Oldroyd-B model.