In this paper, we present a parallel domain decomposition algorithm for the simulation of blood flows in patient-specific artery. The flow may be turbulent in certain situations such as when there is stenosis or aneurysm. An accurate simulation of the turbulent effect is important for the understanding of the hemodynamics. Direct numerical simulation is computationally expensive in practical applications. As a result, most researchers choose to focus on a portion of the artery or use a low-dimensional approximation of the artery. In this paper, we focus on the large eddy simulation (LES) of blood flows in the abdominal aorta. To make the model more physiologically accurate, we consider a resistive outflow boundary condition which is more accurate than the traditional traction free condition. Different from the common decoupled approach where the resistive boundary condition is pre-calculated and then imposed as a Neumann condition, we prescribe it implicitly as an integral term on the LES equations and solve the coupled system monolithically. The governing system of equations is discretized by a stabilized finite element method in space and an implicit second-order backward differentiation scheme in time. A parallel Newton–Krylov–Schwarz algorithm is applied for solving the resulting nonlinear system with analytic Jacobian. Due to the integral nature of the resistive boundary condition, the Jacobian matrix has a dense block corresponding to all the variables on the outlet boundaries. Impacts of the resistive boundary condition with different parameters on the simulation results and the performance of the algorithm are investigated in detail. Numerical experiments show that the algorithm is stable with large time step size and is robust with respect to other parameters of the solution algorithm. We also report the parallel scalability of the algorithm on a supercomputer with over one thousand processor cores.