Summary: A geometric multigrid algorithm is introduced for solving nonsymmetric linear systems resulting from the discretization of the variable density Navier-Stokes equations on nonuniform structured rectilinear grids and high-Reynolds number flows. The restriction operation is defined such that the resulting system on the coarser grids is symmetric, thereby allowing for the use of efficient smoother algorithms. To achieve an optimal rate of convergence, the sequence of interpolation and restriction operations are determined through a dynamic procedure. A parallel partitioning strategy is introduced to minimize communication while maintaining the load balance between all processors. To test the proposed algorithm, we consider two cases: 1) homogeneous isotropic turbulence discretized on uniform grids and 2) turbulent duct flow discretized on stretched grids. Testing the algorithm on systems with up to a billion unknowns shows that the cost varies linearly with the number of unknowns. This \(\mathcal{O}(N)\) behavior confirms the robustness of the proposed multigrid method regarding ill-conditioning of large systems characteristic of multiscale high-Reynolds number turbulent flows. The robustness of our method to density variations is established by considering cases where density varies sharply in space by a factor of up to \(10^{4}\), showing its applicability to two-phase flow problems. Strong and weak scalability studies are carried out, employing up to \(30,000\) processors, to examine the parallel performance of our implementation. Excellent scalability of our solver is shown for a granularity as low as \(10^{4}\) to \(10^{5}\) unknowns per processor. At its tested peak throughput, it solves approximately 4 billion unknowns per second employing over \(16,000\) processors with a parallel efficiency higher than 50%.