As was shown by the Authors in recent papers, the computational efficiency of the commonly- used drift-diffusion model can be increased significantly by recasting the equations such that the potential is obtained fromthe generalized Ohm’s law rather thanGauss’s law and by adding some source terms to the ion transport equations to ensure that Gauss’s law is satisfied either in quasi- neutral or non-neutral regions. Not only was such approach demonstrated to reduce the stiffness of the systemleading to faster convergence but it was also shown to result in a higher resolution of the converged solution. The combined gains in convergence acceleration and resolution amounted to a one-hundredfold increase in computational efficiency when simulating glow discharges and plasmas with significant quasi-neutral regions. The gains in efficiency were also observed to be unaffected by an external magnetic field. In this paper, it is demonstrated that such a recast of the drift-diffusion model has yet another advantage: its lack of stiffness permits the electron and ion transport equations to be integrated in coupled form along with the Reynolds-averaged Navier- Stokes equations. Test cases relevant to plasma aerodynamics (including magnetic field effects, negative ions, and multiple positive ions) demonstrate that such a coupled system of equations can be converged in more or less the same number of iterations and computing time as that describing non-ionized non-chemically-reacting flows while not sacrificing the generality of the drift-diffusion model.