We present and study a parallel subgrid stabilized algorithm for simulating the incompressible Navier-Stokes equations with nonlinear slip boundary conditions. The algorithm uses finite element discretizations and an approach of completely overlapping region division for the parallelization, where an elliptic projection operator is applied to define the stabilization term. It has the following appealing features: 1) each subproblem used to calculate a local solution in a subdomain is actually a global problem defined on a global mesh that is locally refined around the subdomain; 2) it can re-use existing sequential software in coding, where both parallelization and stabilization are easy to implement based on existing sequential codes without extensive effort; 3) it inherits the advantages of the subgrid stabilization method and the completely overlapping region division approach; 4) with suitable algorithmic parameters, it is able to yield an optimal convergence rate for the approximate solutions with a comparable accuracy to that of the solutions from the global subgrid stabilized method, with considerable reduction in computational time. With the help of a local a priori estimate, we estimate error bounds of the obtained solutions from the algorithm, and perform some numerical tests to validate the theoretical prediction and illustrate the applicability of the proposed algorithm.