An efficient multigrid approach for solving a discretized elliptic equation whose boundary values are determined in part by integral relations is developed, analyzed and tested. The algorithm is implemented on a problem that is solved during integration of the three-dimensional Quasigeostrophic equations, which model large-scale rotating stratified flows, where the integral constraints represent mass conservation.