The paper deals with the discretization of the optimal control problem \[ \text{minimize } J(v,q)= \tfrac 12\| v-v_d\|^2_{L^2(\Omega)^d}+\tfrac \nu2 \| q\|^2_{L^2(\Omega)^d},\quad \nu>0 \]
subject to the Stokes equations
\[ -\triangle v+\nabla p=f+q,\;\nabla\cdot v=0\;\text{ in}\;\;\Omega,\;v=0\;\text{ on}\;\Gamma \]
and subject to the control constraints \(a\leq q(x)\leq b\) for a.a. \(x\in \Omega\), where \(\Omega\) is a bounded domain in \(\mathbb {R}^d\), \(d=2,3\) and \(\Gamma=\partial\Omega.\) The quantities \(a,b\in \mathbb {R}^d\) are constant vectors, the control constraints are understood pointwise. Both the control and state variables are discretized by finite elements. The optimal control problem is strictly convex and radially unbounded. Hence there exists a uniquely determined optimal solution and the first order necessary conditions are also sufficient for optimality. In order to discretize the optimal control and the adjoint problem the authors consider a \(2\)- or \(3\)-d mesh consisting of open cells which are either triangles, tetrahedra, quadrilaterals or hexahedra. In contrast to previous papers on the subject the control variable is approximated by piecewise constant elements. Using the postprocessing strategy connecting with superconvergence results the authors achieve error estimates of the order \(h^2\) in \(L^\infty\)-norm which is a significant improvement of the accuracy.