Using TEOS-10 routines for density

[bolding]

On gotm-devel we started a discussion on using TEOS-10 instead of the present density routines in GOTM - https://groups.google.com/u/1/g/gotm-devel/c/njWS8GBx0gw.

Based on experience the discussion becomes easier to follow if we use the 'Issues' feature of GitHub.

So that is why I've moved it here.

[bolding]

Comments to Lars's mail.

As TEOS-10 provides routines for conversion between in-situ. potential and conservative temperature + practical and absolute salinity I think this part is easy. Basically our T ans S must be conservative temperature and absolute salinity. When we read in profiles we (assume) in-situ temperature and practical salinity. The conversion requires specification of pressure in dbar - for now I use depths in meters as a proxy. Somebody can at a later stage investigate the difference using real pressure :-)

I looked at the GSW routine - gsw_nsquared() - it is not easy to convert to user provided pressure points. I know Jorn has worked on it in a Python version of gsw_nsquared(). It involves manually interpolating S and T to - and calculating alpha and beta - for each pressure level. It would be good to see what the error is by using a zoomed grid instead as if it is equidistant.

I have some preliminary timings for using gsw_nsquared() instead of old GOTM routines - and it is much faster.

https://github.com/gotm-model/code/blob/gsw/src/util/equation_of_state.F90 contains the present version - documentation not updated - but I know Hans will love to do that.

Density is used by some models in FABM and will be calculated - and saved.

Makes it easier to use the same rho0 throughout - but maybe also more confusing - it is not the variable changed most often.

Concerning horizontal pressure gradients we assume gradients in S and T are given in conservative and absolute terms. And to keep it potential p should be p0 (likely 0) in the calculations - and not an array.

I've updated some of the code - notably equation_of_state.F90 and will commit other changes as I progress.

Note - the code in the repository will not compile and run right now.

I'll also commit a gotm.yaml.gsw to the nns_annual test case to be used for actual testing.

[lumlauf]

Regarding the computation of N2, two options are on the table at the moment:

• We compute it based on the local gradient of potential density at local reference pressure, as before.

• We use the TEOS-10 implementation in gsw_Nsquared(). This has the problem that, in GOTM, the pressure of the interfaces is not exactly the average of the pressure of the adjacent T/S-points for non-equidistant grids. This, however, is assumed in gsw_Nsquared().

My suggestion would be to stick with the second solution because the TEOS-10 guys know what they're doing, and their code might also be more efficient.

Instead of modifying their code (which would be difficult to keep in sync with future TEOS-10 releases, as Hans noted already), we could simply assimilate the content of gsw_Nsquared() in our own GOTM code, including the modified pressure points. I just checked the Matlab version of the code in gsw_Nsquared(). The code is not very difficult. The interface to TEOS-10 would only consist of two function calls:

gsw_grav() ,
gsw_specvol_alpha_beta() ,

which compute gravity and the expansion coefficients. These interfaces are very unlikely to change, and they are otherwise easy to modify because their function is clear. Any other alternative would also have an interface, e.g. via gsw_rho(), that we would need to maintain. So not much would be gained.

Lars

[bolding]

On the Fortran side there are two problems - https://github.com/TEOS-10/GSW-Fortran/blob/master/toolbox/gsw_nsquared.f90

p_mid is intent(out) so we can't be sure to pass an array in.

2. the 0.5 coefficient would have to be changed as a weighted average of distances instead

forall (k = 1: nz-1) -
grav_local(k) = 0.5_r8*(grav(k) + grav(k+1))
dsa(k) = (sa(k+1) - sa(k))
sa_mid(k) = 0.5_r8*(sa(k) + sa(k+1))
dct(k) = (ct(k+1) - ct(k))
ct_mid(k) = 0.5_r8*(ct(k) + ct(k+1))
dp(k) = (p(k+1) - p(k))
p_mid(k) = 0.5_r8*(p(k) + p(k+1))
end forall

likely modifying this routine and have it adopted by the TEOS-10 people would be the right thing to do - but it will take time.
p_mid should be intent(inout) and the scale factor would default to 0.5 but could be re-calculated based on values in a provide p_mid.

[lumlauf]

What I meant with "assimilating" in my previous comment was that we simply copy their code in gsw_Nsquared() into a new routine that will become part of GOTM, and that we can modify as we wish. So we cold easily change things like intent() and the 0.5 factor. This new routine would then contain the interface to the two TEOS-10 functions mentioned above, and we would only need to make sure that these interfaces are kept up to date. I don't know if that was clear?

We could in parallel (or better before we start) talk to the TEOS-10 people and see if there is any hope they will modify their code accordingly.

[bolding]

The code is now compilable and runable for the nns_annual with the included gotm.yaml.gsw

Here are the steps:
cd /tmp
git clone -b gsw --recursive --single-branch git@github.com:gotm-model/code.git gotm-gsw
cd gotm-gsw/ && cmake -B build && cmake --build build --target
./build/gotm -v

should show (gsw branch) and TEOS-10 version

now cd to the folder of the nns_annual setup:
/tmp/gotm-gsw/build/gotm /tmp/gotm-gsw/gotm.yaml.gsw

Capital GSW in the code will show the places with significant changes

[bolding]

Does this make sense? Lars suggest rho0 should be the same used in meanflow as the Boussinesq approximation.

I can change the order of the different methods. I also think about chaning to alpha and beta to be consistent with TEOS-10 naming conversion.

[lumlauf]

This looks great, Karsten. It is also good to have the rho0 factor defined here in the context of the EOS. And I am very much in favor of finally getting rid of the dtr0 and dstr0 factors, and replacing them by the standard alpha/beta expansion factors. I ALWAYS forget where to put the rho0 factor to convert them into each other...

[bolding]

So we should move rho_0 from meanflow to here? That might be a good idea and with less confusion. Then it would not be below the linear: section - but just below the method. Then it would be defined only one place.

[lumlauf]

Yes, this sounds like the best solution, at least to me.

[bolding]

I think gravity will go the same way.

[lumlauf]

Hm, this is a more difficult one. I'm aware that TEOS-10 includes a routine to compute gravity as a function of latitude for internal use (e.g., to convert pressure to depth, etc.). Still, however, it should be clear that gravity is definitely not a thermodynamic quantity, and should therefore not appear in the EOS section. So I would leave it where it is, and simply ignore that TEOS-10 internally uses a (minimally) different version.

[bolding]

Back to this - after merging the plume code - and hope to close it soon.

As discussed by Lars above we can use the method from present GOTM or use gsw_Nsquared for the calculation of NN. Note also the issues mentioned by Lars regarding non-equidistant grids.

I've done some testing between the two for equidistant grids - i.e. ddu and ddl are both 0 - for the nns_annual case.

This is the difference of NN

with the following impact in temperature

Concerning timings - gsw_Nsquared wins straight up.

TIMINGS 9.0000000003698233E-006 4.9999999999883471E-005
TIMINGS 9.9999999996214228E-006 5.6000000004274852E-005
TIMINGS 9.9999999996214228E-006 6.1999999999784450E-005
TIMINGS 1.0000000001397780E-005 9.8999999998738986E-005

being multiple times faster.

Next question is to deal with the non-equidistant case - and if it matters at all with the small dz we use.

Note - that the method for testing the gsw branch I gave above is still valid (after I've committed the code I have on my workstation)

[bolding]

To view the difference between the gsw and master branch - try

git fetch --all
git difftool -y ..origin/master

To compile old, new and both line 89, 90 allows that. Comment in/out.

In stratification.F90 3 files are written for looking at differences the files are written from line 190 and forward

I've renamed some variables in the NetCDF file to keep e.g. all temperatures together in pyncview.

I've added a yaml file that can be used with the nns_annual case

gotm.yaml.zip

[bolding]

I would really like to have this merged in as I have to merge with changed master again and again.

Before merging into master I need somebody to at least have a look at it - try it out.

Update gotm.yaml can be generated by running the compiled executable - some sections have been renamed, a few new variables have been added - but relative minor changes to the outside.

I cant remember the versioning scheme we use - but most others seems to have changed to just a sequential scheme - where e.g. 7 follows 6 - and there are no stable/unstable version. Comments?

If there are no feedback within a month I merge the gsw branch in!!