calculating zostoga

[paolap]

There's few functions involved in this in original app and want to make sure before I convert them to xarray that there's not some function we can import from elsewhere.
From APP4 app_functions.py roughly:

tmp=((1. -
   rho_from_theta(T[t,:],35.00,press)/rho_from_theta(4.00,35.00,press))*dz[t,:]).sum(0)
zostoga[t]=(tmp*areacello).sum(0).sum(0)/areacello.sum()

where T is pot_temp
rho_from_theta() is a function to calculate density from temperature, salinity (constant=35.0), pressure
and pressure itself is calculated from depth and latitude

I should be able to replace these two functions with gsw package functions
gsw.conversions.p_from_z
and
one of the gsw.density.sigma functions from https://teos-10.github.io/GSW-Python/density.html

Not sure how different is conservative temperature form potential temperature in the model and which one of the functions to use!!

[anton-seaice]

@taimoorsohail - we thought/ hoped you might have a suggestion for this ?

[taimoorsohail]

Thanks for the ping @anton-seaice.

In most cases, the temperature variable in models (i.e., potential temperature) and the salinity variable (often practical salinity) can actually be considered as the same as "conservative temperature" and "absolute salinity" to be used in TEOS-10 functions, from McDougall et al., 2021:

"... the most appropriate interpretation of the model's temperature variable is actually Conservative Temperature. This perhaps unexpected interpretation is needed to ensure that the air–sea heat flux that leaves and arrives in atmosphere and sea ice models is the same as that which arrives in and leaves the ocean model."

and "... the only interpretation for the model's salinity variable is Preformed Salinity, and the use of the TEOS-10 equation of state will then yield the correct specific volume.".

There are slight differences on whether the model equation of state is based on EOS-80 or TEOS-10, but either way you can assume that the native model T and S can be used in place of CT and SA in the gsw functions. See Section 2.7 here for more information on the difference between EOS-80 and TEOS-10 models.

So, the right functions in gsw would be gsw.conversions.p_from_z.

You are calculating the fractional change in in-situ density in the APP4 code snippet, so ideally you'd want to use gsw.density.rho in place of rho_from_theta above.

Hope this helps!

[taimoorsohail]

@rmholmes was a co-author on the McDougall et al., 2021 paper and may have more insight.

[paolap]

Thanks this all very useful, just generally I would like to make sure all the calculations we are adding are somehow endorsed by the community and model experts. In fact, ideally we should just be able to access them from a repository that provides the core of the calculations and have only to create wrappers around them, so they can work in this tool framework and with different versions of the model.
So far for the ocean I added meridional overturning calculation, I was able for that to refer to cosima cookbook but as it exists as a notebook I couldn't import it.

[rmholmes]

@rmholmes was a co-author on the McDougall et al., 2021 paper and may have more insight.

Sorry for the slow reply. What @taimoorsohail said looks right to me. I think you can take this interpretation of the model temperature as CT (despite its labelling, in EOS-80 models, as potential temperature) as community endorsed - it will be included as such in an upcoming "recommendations for CMIP7" paper (basically an update to https://www.geosci-model-dev.net/9/3231/2016/).

[paolap]

I managed to get a working function together, here are some test results using data from CMIP6:

[0.51817106 0.52007097 0.52156734 0.52226752 0.52083167 0.51813325] the new function using gsw functions
[0.5203042 0.5222081 0.5237081 0.5244111 0.52297163 0.5202664 ] the APP4
[0.52065015 0.52255523 0.5240563 0.5247597 0.5233193 0.5206123 what’s in the published data

The last two should be the same, but they might be differing because these runs usually have some issues with time axis and I didn't took that into account.
There are still some difference but they are of the order of 0.002 so negligible considering this is a sum over all ocean points.
overall is more likely that the gsw functions are "more" correct (or at least more consistent with others calculating the same) than the APP4 ones.
Changes are in branch issue199 and will be push soon to a new release