Hi Sougata,
The model uses the multi-dimensional piece-wise parabolic method (MDPPM, Colella & Woodward, 1984) for both horizontal and vertical tracer advection with a monotonicity-preserving flux limiter following Suresh and Huynh (1997). But I’m not convinced that will be useful to you. @adele157 mentions here that she has reproduced this calculation offline, so maybe you could ask her.
Another inaccuracy you will have is that the time-averaged budget diagnostics (like temp_xflux_adv) will relate to the difference in temperature at the two instants in time at the beginning and ending of the averaging time period (this is what I mean by “snapshots”), as opposed to the differences in time-averages that I presume you are using there. That might be a minor difference for your situation, I’m not sure.
What depths are you doing this analysis over? Are you trying to do a mixed layer heat budget? If so, there are another few issues (I’m working on a robust way of doing this now).
Please also keep in mind that the individual components Xadv, Yadv, Zadv do not have a robust physical meaning as they depend on the temperature scale that you are using (i.e. they could be very different if you use Celsius vs. Kelvin). Their sum should be robust to this (if the volume of the region isn’t changing).