A technique to determine the large-scale time-averaged ocean circulation from hydrographic observations and surface flux estimates is described. It is based on an inversion of the Bryan-Cox ocean general circulation model. We have constructed an approximate adjoint to that model which is computationally simpler and more economic than the exact adjoint. The optimization algorithm, although not optimal in a statistical sense, allows calculation of all state variables such that they are consistent with the equilibrium dynamics of the circulation model and agree as closely as possible with the observed data. To verify the technique, we have performed identical twin experiments with the circulation model in an idealized geometry. It is found that in principle the true model state including surface fluxes can be recovered with acceptable accuracy, even if no information on the surface fluxes is available. Under ideal conditions, the resulting rms errors of the surface fluxes were as low as 3 W/m2 and 0.2 m/year for heat and freshwater, respectively. Regions of deep-water formation due to convection show however larger errors on a small spatial scale, depending on the nonlinear, threshold-like nature of convective adjustment. The optimized solutions are distinctly sensitive to the integration time interval, with optimal values around three years. The results suggest that the procedure is suitable to obtain a consistent description of the oceanic state, and in particular more accurate estimates of the air-sea heat and freshwater fluxes.