[MITgcm-support] How to get correct buoyancy from output variables

[Lo Piccolo, Anna]

Hi all,

I am running some hydrostatic simulations, and I'm trying to close the
buoyancy budget. I've found that if I compute it in two ways:

1) As the vertical derivative of the hydrostatic pressure (using a
centered difference), using the diagnostic output PHIHYD

2) As -g*RHOANOMA/RHOCONST

I end up with discrepancies even at t=0  (0.1% relative difference
around the mixed layer). It is possible that this comes from the
numerical derivative in PHIHYD, but I want to make sure I am doing
this correctly. Given the difference, what would be the best way to
compute the buoyancy given how it is used to force momentum?
Thank you.

Best regards,
Anna Lo Piccolo