Thanks @taimoorsohail.
Just to clarify for others - the problem is with the pre-processing that computes the ESMF/CMIP output variable hfds from the raw MOM5 output diagnostics. The old calculation (“abhi-hfds” above) did not include frazil formation below the surface layer (frazil-2d vs. frazil-3d) or some components of the heat flux associated with surface volume fluxes (a few terms to do with ice-ocean fluxes). The “ryan-hfds” calculation above correctly includes all of these terms.
Technical note: To close the heat budget exactly you should technically also be comparing the time integral of area-integrated hfds to snapshots of the global OHC, not time-averages of OHC. This point is minor compared to the difference that Taimoor is showing above in the first plot, as shown by the second plot.