Baroclinic sea level calculation

[c2xu]

This replaces PR #984 and is based on algorithms in a manuscript that we plan to submit to JPO. Reference to this work will be added as soon as it's available.

[Hallberg-NOAA]

The comment here should make it clear that this is the anomaly in the pressure change across a layer relative to what would be found if the density of the layer were rho_s. Without this clarification, it is very hard to see how the code in this routine would generate the baroclinic sea level anomaly. This comment also pertains to the descriptions of dg, dp_int, dg_int and p_int in the next few lines, all of which are also anomalies relative to the reference state.

[Hallberg-NOAA]

This call to find_eta is only needed in Boussinesq mode. Because this can be an expensive call in non-Boussinesq mode, it should be avoided when it is not needed. In fact, given the substantial differences in the algorithms and the meanings of the variables between Boussinesq and non-Boussinesq modes, I would be tempted to go so far as to split the body of this routine into two separate routines that are called depending on whether the Boussinsesq approximation is being made.

[Hallberg-NOAA]

I have worked through all of the calculations in this revised PR, and I am convinced that all of the expressions are correct (or at least plausible). The one big thing that would have helped in working through the code is for there to be more explicit comments describing some of the variables in find_bsl() to make it clear that they are anomalies relative to what would be found if the density were RHO_BSL, and that they can be negative. I think that this probably applies to dp, dp_int, and pt in Boussinesq mode while in non-Boussinesq mode it applies to gz, dg, and dg_int and p_int, but not to pb, pt or dp. (Without this reminder in the comments I spent much more time than I should looking for the subtraction without which the calculations do not make sense!) Given the very different meanings of these variables, it might even make sense to refactor find_bsl() to have separate internal routines (with a common top-level interfaces) for the Boussinesq and non-Boussinesq cases.

There also needs to be a more explicit description in the MOM_parameter_doc description for RHO_BSL explaining that this is the reference density for determining the surface height anomalies from on the baroclinic density structure based on integrals of the differences between the actual density and this value.

With these relatively straightforward changes to the documentation in the code and perhaps some strategic refactoring greater clarity (but without changing any of the actual calculations), I think that this PR would be ready to be accepted.

[Hallberg-NOAA]

I suspect that there are cases where dividing by bathyT here will not work because in some cases (such as the Great Lakes) bathyT can be equal to or above the reference geoid and still have water. Would it make more sense to use the vertical sum of the layer thicknesses (perhaps in a max with GV%H_subroundoff) in the denominator instead?

[Hallberg-NOAA]

I think that the denominator here should be the sum of the layer thicknesses.

Also, there should be spaces around the semicolon on the next line to follow the pattern used elsewhere.

[chengz2]

Because whether the Boussinesq mode is activated must be determined during the run time, I feel it's difficult to split this subroutine into multiple subroutines under the same interface.

[chengz2]

For BSL calculation, what matters is not whether bathyT is above or below the reference geoid, but where the mean sea level is. The bsl calculation will be incorrect if the mean sea level is non-zero, even if bathyT is below the geoid. bathyT is in fact defined as bathyT(i,j) = G%bathyT(i,j) + dZ_ref. Am I correct that this dZ_ref is a correction to the mean sea level if it's non-zero? If not, is there another variable available in MOM6 that I can use? Or do I have to calculate it manually? Thanks, @Hallberg-NOAA.

[chengz2]

For BSL calculation, what matters is not whether bathyT is above or below the reference geoid, but where the mean sea level is. The bsl calculation will be incorrect if the mean sea level is non-zero, even if bathyT is below the geoid. bathyT is in fact defined as bathyT(i,j) = G%bathyT(i,j) + dZ_ref. Am I correct that this dZ_ref is a correction to the mean sea level if it's non-zero? If not, is there another variable available in MOM6 that I can use? Or do I have to calculate it manually? Thanks, @Hallberg-NOAA.

Also, in theory, the denominator should be the distance between the bottom and the mean sea level, instead of the instantaneous sea level, i.e., it should be time-independent. That's another reason I'm trying to avoid using the sum of layer thickness, because if the SSH variation is large relative to the total depth, this will amplify the error.

[Hallberg-NOAA]

Splitting this subroutine into 3 (one that is called from outside, and a Boussinesq version or a non-Boussinesq version that are called from the first one, depending on whether or not a particular run is Boussinesq) should be straightforward. This is the pattern that was adopted in MOM_PressureForce.F90, although in this case I would keep all 3 subroutines in the same module.

[Hallberg-NOAA]

You can accommodate spatially varying time-mean sea levels by using (G%meanSL(i,j) + G%bathyT(i,j)) in place of G%bathyT(i,j).

However, in the cases here it seems to me that we are dividing an integral by the total depth or pressure to give what amounts to a vertical average. If this can be interpreted as a vertical average, I would think that because the bounds of integration in the numerator vary with the rapid time-scales, the normalization in the denominator should also exhibit the same rapid variations.

[c2xu]

Splitting this subroutine into 3 (one that is called from outside, and a Boussinesq version or a non-Boussinesq version that are called from the first one, depending on whether or not a particular run is Boussinesq) should be straightforward. This is the pattern that was adopted in MOM_PressureForce.F90, although in this case I would keep all 3 subroutines in the same module.

find_bsl() is split into three subroutines as requested.

[c2xu]

You can accommodate spatially varying time-mean sea levels by using (G%meanSL(i,j) + G%bathyT(i,j)) in place of G%bathyT(i,j).

However, in the cases here it seems to me that we are dividing an integral by the total depth or pressure to give what amounts to a vertical average. If this can be interpreted as a vertical average, I would think that because the bounds of integration in the numerator vary with the rapid time-scales, the normalization in the denominator should also exhibit the same rapid variations.

I thought about it and decided not to make this change for the following reasons:

1. Barotropic-baroclinic sea level decomposition is likely to be useful only in the deep ocean where baroclinic internal waves are "large", yielding sea level fluctuations of a few cm. I don't think it should ever be calculated at locations such as the Great Lakes.

2. In the deep ocean, mean sea level fluctuations relative to the total depth is on the order of $10^{-3}$, so I don't think it is necessary to replace G%bathyT(i,j) by (G%meanSL(i,j) + G%bathyT(i,j)) in the denominator.

3. This is perhaps the most confusing part: In theory, the denominator, as well as the upper limit of integration, should not consist of temporal variation of the sea level, because this calculation is based on the linear decomposition of the sea level, and the baroclinic sea level is an anomaly to the mean state around which the dynamics are linearized. Introducing temporal dependency to the denominator will introduce a nonlinearity that breaks this linear decomposition. But again, because the temporal variation of the sea level relative to the total depth is on the order of $10^{-3}$, whether or not including this temporal variation in the denominator should not change the answer in the deep ocean in any significant way.

4. That said, the temporal varying, instantaneous sea level, instead of the mean sea level, is used as the upper limit of integration in this algorithm. This is not a problem in general because the density difference between the mean sea level and the instantaneous sea level, if any, should be negligible. Because we use surface density as the reference density, the density anomaly at the surface is essentially zero. Hence, the vertical integral of this density anomaly between the mean sea level and the instantaneous sea level should also be negligible.

I've added some comments explaining some of these aspects in the code. A more complete, mathematically rigours explanation of this algorithm, and the physics behind it, will be available once our paper is published.