||This paper describes in detail a numerical scheme designed for direct numerical simulation (DNS) of turbulent drag reduction. The hybrid spatial scheme includes Fourier spectral accuracy in two directions and 6th-order compact finite differences for first and second-order wall-normal derivatives, while time marching can be up to 4-th order accurate. High-resolution and high-drag reduction viscoelastic DNS are made possible through domain decomposition with a two-dimensional MPI Cartesian grid alternatively splitting two directions of space (’pencil’ decomposition). The resulting algorithm has been shown to scale properly up to 16384 cores on the Blue Gene/P at IDRIS-CNRS, France. Drag reduction is modeled for the three-dimensional wall-bounded channel flow of a FENE-P dilute polymer solution which mimics injection of heavy-weight flexible polymers in a Newtonian solvent. We present results for 4 high-drag reduction viscoelastic flows with friction Reynolds numbers Reτ0=180, 395, 590 and 1000, all of them sharing the same friction Weissenberg number Weτ0=115 and the same rheological parameters. A primary analysis of the DNS database indicates that turbulence modification by the presence of polymers is Reynolds-number dependent. This translates into a smaller percent drag reduction with increasing Reynolds number, from 64% at Reτ0=180 down to 59% at Reτ0=1000, and a steeper mean current at small Reynolds number. The Reynolds number dependence is also visible in second-order statistics and in the vortex structures visualized with iso-surfaces of the Q-criterion.