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

1) 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!!

[bolding]

I have merged the gsw branch into master and updated gotm.yaml files in gotm-cases.

Github actions succeeds :-)

A few things needs to be sorted out:

1) Somebody must check const_NNS.F90, const_NNT.F90, convert_fluxes.F90.

2) stratification.F90 has been cleaned up - and in addition to use gsw_Nsquared() the use of a prognostic equation for buoyancy has been moved out. Right now calculations uses a constant latitude of 0.

3) prognostic_buoyancy.F 90 must be tested

4) Seems boundary conditions for NN and SS are set wrongly (all just set to 0 at index 0 and nlev)- sometime ago I talked to Hans about and we should update those.

Copy of message from Hans: Shear is easy, just divide the respective u_star by (kappa z_0). NN at the bottom should be zero (unless you consider geothermal heating or release of heat stored in the sediment in the Wadden Sea). Most tricky is NN at the surface. For stable stratification, I would divide the surface buoyancy flux by the eddy diffusivity (Prt kappa z0 u_star), where Prt is the turbulent Prandtl number (would be about 0.7 for neutral stratification, and quite tricky for strong stratification). For unstable stratification (e.g. night-time cooling) I would have no clue. I think a good check after implementation would be to plot profiles of NN and SS. The boundary values should be continuous to the adjacent values. Try best with GOTM at high resolution. Greetings from Stege!

[lumlauf]

Hi Karsten,

I just had quick look - but already the first test case I tried (our standard "entrainment" case) did not work properly. It seems that NN is not properly initalized (NN contains some arbitrary numbers or equals zero).

I checked for the reasons, and the implementation of the new TEOS-10 stuff in const_NNS.F90 and const_NNT.F90 is completely missing at the moment. There is just a commented line you apparently inserted.  Also the infrastructure for the different versions of the equation of state (full TEOS-10, linear TEOS-10, linear custom) needs to be generated. In the old version, these decisions were made via eos_alpha() etc., which, however, does not exist any more. So this needs to be fixed before I can do any meaningful testing.

Some comments on your other points below:

On 8/12/2022 11:33, Karsten Bolding wrote:

I have merged the gsw branch into master and updated gotm.yaml files in gotm-cases.

Github actions succeeds :-)

A few things needs to be sorted out:

1.

Somebody must check const_NNS.F90, const_NNT.F90, convert_fluxes.F90.

Also in convert_fluxes, it seems the infrastructure to specify the type of the equation of state (full TEOS-10, linear TEOS-10, linear custom) is lacking. So at the moment, only the full version can be used. Or am I overlooking something here? If this is fixed, I'll be happy to do some testing.

1.

2.

stratification.F90 has been cleaned up - and in addition to use
gsw_Nsquared() the use of a prognostic equation for buoyancy has
been moved out. Right now calculations uses a constant latitude of 0.

OK, but TEOS-10 really requires that lat and lon are known to assume the correct ion composition for salinity. But don't we have lat and lon from the yaml file?

1.

2.

prognostic_buoyancy.F 90 must be tested

Yes, I can do this. But I just realized that the "buoyancy" group in gotm.yaml only contains "NN_ini" in the standard version. I needed to use "--detail full" to also see "method" and "surf_ini". This is very confusing for inexperienced users as they don't even see that a prognostic equation can be solved. Can you please fix this?

1.

2.

Seems boundary conditions for NN and SS are set wrongly (all just
set to 0 at index 0 and nlev)- sometime ago I talked to Hans about
and we should update those.

Copy of message from Hans: Shear is easy, just divide the respective u_star by (kappa z_0). NN at the bottom should be zero (unless you consider geothermal heating or release of heat stored in the sediment in the Wadden Sea). Most tricky is NN at the surface. For stable stratification, I would divide the surface buoyancy flux by the eddy diffusivity (Prt kappa z0 u_star), where Prt is the turbulent Prandtl number (would be about 0.7 for neutral stratification, and quite tricky for strong stratification). For unstable stratification (e.g. night-time cooling) I would have no clue. I think a good check after implementation would be to plot profiles of NN and SS. The boundary values should be continuous to the adjacent values. Try best with GOTM at high resolution. Greetings from Stege!

Well. Regarding NN, I think that zero would be the more correct value. We assume a standard log-layer expression near the surface which implies that NN is zero. Only if we went for some kind of Monin-Obukhov-type boundary conditions, it would make sense to specify a non-zero NN. Regarding SS, I agree with Hans that this parameter could be computed from the log-layer expression he mentions. BUT we also have the option to inject TKE, which would result in a different near-surface velocity structure (not log any more) and thus a different value of SS at the surface. As these boundary values are of purely cosmetic nature (for output) anyway, I suggest to stay with SS=0 (or even insert "missing_values"?) so users can see immediately that they should be careful.

Cheers,

Lars

— Reply to this email directly, view it on GitHub https://github.com/gotm-model/code/issues/25#issuecomment-1212917480, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEVOCRMLGLBAASSR2LLN2ALVYYK4ZANCNFSM5FOBRRWA. You are receiving this because you commented.Web Bug from https://github.com/notifications/beacon/AEVOCRIIH4Y4DIKZVHHL3J3VYYK4ZA5CNFSM5FOBRRWKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOJBF2N2A.gifMessage ID: @.***>

[ { @.": "http://schema.org", @.": "EmailMessage", "potentialAction": { @.": "ViewAction", "target": "https://github.com/gotm-model/code/issues/25#issuecomment-1212917480", "url": "https://github.com/gotm-model/code/issues/25#issuecomment-1212917480", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { @.": "Organization", "name": "GitHub", "url": "https://github.com" } } ]

--

Lars Umlauf

Physical Oceanography and Instrumentation Leibniz-Institute for Baltic Sea Research

phone : ++49 381 5197 223 fax : ++49 381 5197 114 web :www.io-warnemuende.de/lars-umlauf-en.html

address: Leibniz-Institute for Baltic Sea Research Seestrasse 15 D-18119 Rostock-Warnemuende Germany

[bolding]

Hi Lars

My first point was that somebody needs to look at const_NNS.F90, const_NNT.F90, convert_fluxes.F90 :-) So it was known configurations using those would not work.

TEOS-10 comes with 3 different routines relating to the calculation of alpha and beta - toolbox/gsw_rho_alpha_beta_bsq.f90 toolbox/gsw_rho_alpha_beta.f90 toolbox/gsw_specvol_alpha_beta.f90

where the difference between the versions are if rho or specvol is calculated - and the unit of alpha and beta. Which one do you prefer?

Don't you have something like this?

convert_fluxes() - the same logic as before - as far as I can see. Again the right gsw routine should be chosen. And then the interface can be be cleaned up.

Concerning what to include in minimum, default and maximum configuration - I think that is for later - let us make it run first.

Karsten

[bolding]

added the calls for alpha and beta calculations - must be checked which version to use

[lumlauf]

Hi Karsten

On 8/15/2022 08:55, Karsten Bolding wrote:

Hi Lars

My first point was that somebody needs to look at const_NNS.F90, const_NNT.F90, convert_fluxes.F90 :-) So it was known configurations using those would not work.

TEOS-10 comes with 3 different routines relating to the calculation of alpha and beta - toolbox/gsw_rho_alpha_beta_bsq.f90 toolbox/gsw_rho_alpha_beta.f90 toolbox/gsw_specvol_alpha_beta.f90

where the difference between the versions are if rho or specvol is calculated - and the unit of alpha and beta. Which one do you prefer?

All these functions compute exactly the same alpha and beta. However, they also compute a lot of other stuff that usually is not needed. So I guess it would be smartest from a computational efficiency point of view to always use the "slimest" function for a particular purpose. For example, in const_NNT(), only alpha is needed. So why not simply use gsw_alpha() instead of the more complicated functions? Same for const_NNS().

Regarding the "unit" for the fixed alpha and beta (specified in the yaml-file), I suggest to go for the official TEOS-10 conventions (1/K for alpha, and kg/g for beta) rather than what we had until know in GOTM. Everything else would be confusing.

Don't you have something like this? image https://user-images.githubusercontent.com/11144809/184589172-7b0cebe5-d0c2-4bc0-afbb-0f8e619369e2.png

What a tried to say in my last email was the following: at the moment, the functions const_NNT() and const_NNS() do not "know" if the user has selected the full teos-10 version, or one of the two simpler linearlized versions in the yaml file. At the moment, always the full version is used. So the infrastructure (passing of the fixed alpha, beta, decision parameters, ...) must be implemented before I can do the testing. Or am I overlooking something important here?

convert_fluxes() - the same logic as before - as far as I can see. Again the right gsw routine should be chosen. And then the interface can be be cleaned up.

Same problem with the lacking infrastructure for chosing the linearized/full versions.

Concerning what to include in minimum, default and maximum configuration - I think that is for later - let us make it run first.

Agreed.

Lars

Karsten

— Reply to this email directly, view it on GitHub https://github.com/gotm-model/code/issues/25#issuecomment-1214681204, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEVOCRJ45EWWPLHKILR2N2TVZHSWBANCNFSM5FOBRRWA. You are receiving this because you commented.Web Bug from https://github.com/notifications/beacon/AEVOCRJBV3VM3BSQKT342YTVZHSWBA5CNFSM5FOBRRWKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOJBTJA5A.gifMessage ID: @.***>

[ { @.": "http://schema.org", @.": "EmailMessage", "potentialAction": { @.": "ViewAction", "target": "https://github.com/gotm-model/code/issues/25#issuecomment-1214681204", "url": "https://github.com/gotm-model/code/issues/25#issuecomment-1214681204", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { @.": "Organization", "name": "GitHub", "url": "https://github.com" } } ]

--

Lars Umlauf

Physical Oceanography and Instrumentation Leibniz-Institute for Baltic Sea Research

phone : ++49 381 5197 223 fax : ++49 381 5197 114 web :www.io-warnemuende.de/lars-umlauf-en.html

address: Leibniz-Institute for Baltic Sea Research Seestrasse 15 D-18119 Rostock-Warnemuende Germany