This work outlines a second order accurate, coupled, conservative new numerical scheme for solving a two dimensional incompressible turbulent flow filed. Mean vorticity, ω, and mean stream function, ψ, are used as the mean flow dependent variables. The turbulent kinetic energy k and the turbulent energy decay rate, ϵ, are used to define the turbulent state. In the present computational scheme two systems of equations and variables are considered: the mean flow system, ψ-ω, and the turbulent state system, k-ϵ. Every system is solved implicity in a coupled double loop manner, and all the flow equations are solved iteratively in the global sense. Since the turbulence boundary conditions have a non-regular variation near a solid wall, they are coupled to the equations implicitly in both systems. In this way the numerical instabilities due to the irregular form of the equations near the solid walls are suppressed. The rate of convergence of the new numerical scheme of the coupled systems ψ-ω and k-ϵ is twice that realized when solving these equations separately. The necessary conditions for convergence of the numerical equations are investigated as well as the rate of convergence features. The detailed stability conditions are derived. As an example, the axisymmetric mixing of two confined jets with an internal heat source is considered with this numerical scheme.