In this work, we present a significantly improved numerical algorithm that has been developed in order to perform direct numerical simulations (DNS) of viscoelastic turbulent channel flow. As in previous work [Housiadas and Beris, Phys. Fluids, 15 (2003) 2369], a spectral approximation is used for every dependent variable that involves Fourier modes along the streamwise and spanwise periodic directions, and Chebyshev orthogonal polynomials along the shear direction. However, we introduced several improvements, the most important of which are: (a) a fully implicit second order time integration scheme, implemented using two nested iteration loops, applied to the constitutive and the momentum equation, respectively, (b) dealiasing along the periodic directions for all the non-linear terms and (c) machine accurate satisfaction of the divergence-free velocity incompressibility condition. The new code, described in this work, is much more stable than the previous one and it has been validated in simulations performed with the FENE-P constitutive equation against results obtained with the previous code. In addition, the new code allowed us to perform successful simulations with the Giesekus and the Oldroyd-B models at simulation conditions that were previously unattainable. A comparison between the results of the Oldroyd-B, Giesekus (alpha = 1/900) and the FENE-P (L = 30) constitutive models offered here shows, for the same Weissenberg number, significant increases to the drag reduction as we move from FENE-P to Giesekus and to Oldroyd-B. (C) 2004 Elsevier B.V. All rights reserved.