I followed the instructions in the remark below method 3.

This:

experiment = '025deg_jra55_ryf9091_gadi'

V = cc.querying.getvar(experiment, 'ty_trans', session, frequency = '1 monthly', n = 3, use_cftime = True)

θ = cc.querying.getvar(experiment, 'temp', session, frequency = '1 monthly', n = 3, use_cftime = True)
θ = θ.interp(yt_ocean = V.yu_ocean.values, method = "linear").rename({'yt_ocean': 'yu_ocean'})

Cp = 3992.10322329649

Vm = V.cf.mean('longitude').cf.mean('vertical')
Vnm = V - Vm
mht_method3 = (Cp * Vnm * θ).cf.sum('longitude').cf.sum('vertical').mean('time')

mht_method3.plot()

gave me this

which looks wrong! What did I do wrong?

I vote that we do the "method 3" here the right way and explain nicely. At the moment it's done in some wrong way and we have a caveat at the bottom. But even when I tried to read the caveat and do it the "right" way it wasn't trivial -- I failed apparently.

I think these difficulties here were part of @dhruvbhagtani's difficulties in a recent project of them. So, it'd be nice if we can protect future scientists that want to reproduce these calculations from the troubles and puzzles I was into, e.g., this morning while doing the above (and potentially some of the same puzzles and troubles that @dhruvbhagtani was few months back).

cc @dhruvbhagtani

Does Vm = V.cf.mean('longitude').cf.mean('vertical') properly take into account missing values inside land?

It's easy enough to check whether the answer still depends on the reference temperature, just evaluate the following and see if it gives you the same answer:
mht_method3_difRef = (Cp * Vnm * ( θ + 100.) ).cf.sum('longitude').cf.sum('vertical').mean('time')

I agree it would be nice to illustrate these differences in more detail. However, the interpretation and choice of what is the best approach is not necessarily that simple. A few notes:

1. Methods 1 and 2 are also reference temperature dependent, so you probably also want to retain the original method 3 to compare to these methods. Note that the reference temperature dependence of methods 1 and 2 could be slightly different, as method 2 will only depend on the reference temperature through the net surface mass flux at different latitudes, while methods 1 (and 3 as written) also depend on the reference temperature through tendencies in mass (sea level - likely small).
2. It would be worth referencing the internal vs. external heat transport discussion in Holmes et al. 2019. The reference temperature independent heat transport computed by subtracting the net volume transport should be equivalent to the internal heat transport.
3. There are other approaches to this problem as well. Some have argued in the past for particular choices of reference temperature depending on their application (e.g. Vranes et al. 2002). Others have used divergent-rotational decompositions (e.g. Forget and Ferriera 2019).

Thanks!! I'll give it a go!

But note that I just did whatever I thought I understood from the caveat comment below method 3. So in any case we should update that to be clearer. If all methods are reference-temperature-dependent then perhaps we should mention that because as a non-MHT expert, from reading the notebook I understood that the caveat was only related to method 3.

Does Vm = V.cf.mean('longitude').cf.mean('vertical') properly take into account missing values inside land?

Not sure if here you are referring to cf-xarray, but the code above does exactly what Vm = V.mean(xt_ocean).mean(st_ocean) would do. Not sure if it takes into account land properly. @dhruvbhagtani, do you know?

From what I understand, there is still ambiguity regarding the third method's temperature scale dependence. As of now, yes, the method is temperature scale dependent because there is a small flow at each latitude which we don't remove.

I am a bit torn between providing as much information to someone reading this notebook and being clean and concise. I think the notebook is already quite dense at the moment for someone just starting to understand calculating MHT, and I would vote to keep the code as it is. By this, I vote to keep method 3 temperature scale dependent as it stays consistent with other methods.

While we are on the topic of differences between the three methods, one more is that method 1 and 2 use diagnostics from 025deg_jra55v13_iaf_gmredi6, while method 3 uses diagnostics from 025deg_jra55_ryf9091_gadi. So, it's not an apples to apples comparison anyways.

@navidcy you made a good point about the fact that writing the caveat under method 3 may mislead someone into believing that the other two methods are temperature scale independent. I can edit the wording.

@navidcy I also saw a comment about which method is the best approach. While I agree with @rmholmes that we can't point to one method universally, we can hint the reader towards which method could work best for them. For example (correct me if I'm wrong here) (i) we can use method 2 only in steady state, or (ii) method 1 uses a model diagnostic which might not be present, and method 3 can serve as a good approximation (the accuracy will depend on the frequency of diagnostics and we'll also miss out on correlations between velocity and temperature).

I will go through method 3 and figure out how to continue using pint and not change our results. As of now, I don't know if both methods described above take into account missing values inside land but I can check that too.

The problem with using pint is that (I think) that it enforces doing all computations in Kelvin. Thus, since the method depends on the reference temperature it gives different results!

So, don't worry about pint for now, but let's ensure that the notebook doesn't give the impression that only method 3 is reference-temperature dependent (as I inferred).

If you can do the rewording, I can deal with pint.

Sure, I'll push my changes soon.

From what I understand, there is still ambiguity regarding the third method's temperature scale dependence. As of now, yes, the method is temperature scale dependent because there is a small flow at each latitude which we don't remove.

I am a bit torn between providing as much information to someone reading this notebook and being clean and concise. I think the notebook is already quite dense at the moment for someone just starting to understand calculating MHT, and I would vote to keep the code as it is. By this, I vote to keep method 3 temperature scale dependent as it stays consistent with other methods.

While we are on the topic of differences between the three methods, one more is that method 1 and 2 use diagnostics from 025deg_jra55v13_iaf_gmredi6, while method 3 uses diagnostics from 025deg_jra55_ryf9091_gadi. So, it's not an apples to apples comparison anyways.

@navidcy you made a good point about the fact that writing the caveat under method 3 may mislead someone into believing that the other two methods are temperature scale independent. I can edit the wording.

@navidcy I also saw a comment about which method is the best approach. While I agree with @rmholmes that we can't point to one method universally, we can hint the reader towards which method could work best for them. For example (correct me if I'm wrong here) (i) we can use method 2 only in steady state, or (ii) method 1 uses a model diagnostic which might not be present, and method 3 can serve as a good approximation (the accuracy will depend on the frequency of diagnostics and we'll also miss out on correlations between velocity and temperature).

I will go through method 3 and figure out how to continue using pint and not change our results. As of now, I don't know if both methods described above take into account missing values inside land but I can check that too.

This all sounds good to me @dhruvbhagtani!