The numerical simulation of blood flows in the human body with a certain level of clinical accuracy is important for the understanding of the human physiology. The success of the modeling relies on a robust numerical method with the corresponding software that can handle the complex geometry, the complex fluid flows and run efficiently on a supercomputer. In this work, we introduce a highly parallel domain decomposition method to solve the three-dimensional incompressible Navier-Stokes equations on a patient-specific artery at the full-body scale from neck to feet with 222 outlets and a minimum diameter around 1.0 mm. A locally refined, unstructured mesh is used to resolve the complex fluid flow. Moreover, a two-level method is introduced to determine the model parameters in the Windkessel outlet boundary condition to guarantee clinically correct flow distributions to 14 major regions. A fully implicit Newton-Krylov-Schwarz method is used to solve the nonlinear algebraic system at each time step and numerical experiments show that the proposed method is robust with respect to the complex geometry, the graph-based partition of the complex mesh, the ill-conditioned sparse systems with locally dense blocks, and different model parameters and is scalable with up to 15,360 processor cores. With the proposed method, one simulation of the blood flow in a full-body arterial network can be obtained in about 8 hours per cardiac cycle, which enables its potential use in a wide range of clinical scenarios.