Historical profiles of ocean temperature and salinity have been mapped to a regular grid using a comprehensive interpolation system that integrates the four-dimensional weighted least squares, the minimum-water-mass-change, and vertical stabilization. The amalgamated interpolation method has multiple advantages. It takes into account not only the three-dimensional spatial data projection but also the temporal signals in the dataset, and also minimizes the unrealistic creation or destruction in salinity-temperature structure when stabilizing and interpolating ocean profiles. Additionally, it avoids the flattened nature of linear interpolation and the over-shooting that is common with cubic splines. The resultant mapping has been used to estimate change of global ocean heat content, and we find that the heat content of the ocean in the middle of last century is smaller than when estimated using linear interpolation, with this difference in heat content being equivalent to about five years worth of global warming.