When high Reynolds turbulent flows are combined with complex and large size geometries, computers are no longer enough powerful to deal with Direct Numerical Simulation (DNS) and with the resolution of all the scales of turbulence motion. Therefore, at a fixed scale suitable for reasonable computations cost, the unresolved scales of the turbulence have to be modeled. These models contain dissipation processes that should not be polluted by the numerical diffusion needed to stabilized approximations for convectiondominated flows. In this paper, we proposed a strongly coupled numerical formulation for the Spalart-Allmaras subscale model, in the framework of Stabilized Finite Element Methods. Computations are performed for compressible Newtonian fluids (2D and 3D) on unstructured grids of high aspect ratio. Results are compared with experimental data and also with solutions obtained by different numerical strategies.