We use the method of least squares with Lagrange multipliers to fit an ocean general circulation model to the Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) estimate of near sea surface temperature (NSST) at the Last Glacial Maximum (LGM; circa 23–19 thousand years ago). Compared to a modern simulation, the resulting global, last-glacial ocean state estimate, which fits the MARGO data within uncertainties in a free-running coupled ocean–sea ice simulation, has global-mean NSSTs that are 2°C lower and greater sea ice extent in all seasons in both the Northern and Southern Hemispheres. Increased brine rejection by sea ice formation in the Southern Ocean contributes to a stronger abyssal stratification set principally by salinity, qualitatively consistent with pore fluid measurements. The upper cell of the glacial Atlantic overturning circulation is deeper and stronger. Dye release experiments show similar distributions of Southern Ocean source waters in the glacial and modern western Atlantic, suggesting that LGM NSST data do not require a major reorganization of abyssal water masses. Outstanding challenges in reconstructing LGM ocean conditions include reducing effects from model biases and finding computationally efficient ways to incorporate abyssal tracers in global circulation inversions. Progress will be aided by the development of coupled ocean–atmosphere–ice inverse models, by improving high-latitude model processes that connect the upper and abyssal oceans, and by the collection of additional paleoclimate observations.