This study demonstrates how the adjoint-based framework traditionally used to compute gradients in PDE optimization problems can be extended to handle general constraints on the state variables. This is accomplished by constructing a projection of the gradient of the objective functional onto a subspace tangent to the manifold defined by the constraint. This projection is realized by solving an adjoint problem defined in terms of the same adjoint operator as used in the system employed to determine the gradient, but with a different forcing. We focus on the “optimize-then-discretize” paradigm in the infinite-dimensional setting where the required regularity of both the gradient and of the projection is ensured. The proposed approach is illustrated with two examples: a simple test problem describing optimization of heat transfer in one direction and a more involved problem where an optimal closure is found for a turbulent flow described by the Navier-Stokes system in two dimensions, both considered subject to different state constraints. The accuracy of the gradients and projections computed by solving suitable adjoint systems is carefully verified and the presented computational results show that the solutions of the optimization problems obtained with the proposed approach satisfy the state constraints with a good accuracy, although not exactly.