Direct numerical simulation (DNS) of incompressible flows in a zero-pressure gradient turbulent boundary layer (TBL) is conducted by a finite difference method in which a fourth order upwind scheme is applied to discretize the convective terms while a re-scaling approach is used to set inlet flow conditions. The Reynolds numbers based on free flow velocity and momentum thickness at the recycle section are respectively 687, 1074, and 1430. The DNS has obtained favorable results indicating that the turbulence statistics are quite satisfactory as compared with the existing numerical and experimental results. The three dimensional turbulent structures at the momentum thickness Reynolds number of 1430 in several different instants are illustrated by the iso-surface of swirl strength square (the square of imaginary part of the complex eigenvalue of velocity gradient tensor) together with velocity vectors in three different cross sections. It is found that there are three kinds of vortical structures: quasi-symmetrical and asymmetrical hairpin vortices, and worm-like vortices. The DNS based on the numerical method can certainly reveal the main characteristics of the TBL flows at the given Reynolds numbers.