Slopes may have sign error, but thickness parameterization is fine (double sign error)

[dhruvbalwada]

I was looking through a few of the MOM6 F90 files where slopes are defined and used:

• MOM_isopycnal_slopes.F90 (e.g line 325)
• MOM_thickness_diffuse.F90 (e.g. line 991).

I noticed that the slope is defined as (for example):

slope = (e(i,j,K)-e(i+1,j,K)) * G%IdxCu(I,j)

I am wondering why is this definition of slope, for example in x direction, $= - d \eta /dx$ used, which also corresponds to ($\partial_x \rho /\partial_z \rho$)? Rather than the more standard definition, which would have a minus in front of what is used ($d \eta/dx$ or $- \partial_x \rho / \partial_z \rho$)? Is this an error?.

So, I checked what the slope diagnostic looks like in an example where the slope should be physically clear. Assuming that MOM6 uses the usual physical definition of slopes : isopycnal shoaling towards the surface as we move north should be positive. This is the setup in Phillips_2Layer (also drawn in Figure 2 of Hallberg 2013).

Here is a figure of the diagnostics from model:

The slope is negative, even though it should be positive.

The concern then was whether the parameterized transport resulting from the thickness package was right? The answer is yes. The model manages to get the right parameterized transport, even with a sign error in slope. This happens because there is seemingly an error in the definition of stream function as well, at least relative to most of the literature, eg Ferrari et al 2010. From this paper equation 9 for $\gamma$ should be considered and the corresponding transport is in equation 5, which matches the definition of $uh^*$ here (the doc denotes $\gamma$ from Ferrari et al 2010 as $\Psi$).
This streamfunction $\gamma = - K S$ (called $\psi$ in MOM6 doc).
However, in the code (eg line 994 in MOM_thickness_diffuse.F90) the stream function is defined as $KS$. Since the slope has a sign error and then there is error in sign of stream function too, the parameterized streamfunction and transport end up with the right signs.

The doc here is also confusing. First it defines $\psi$ as KS, which does match what is in the code but not what we want. Then N2 is defined as $-g \partial_z \rho/ \rho_0$, but then M2 is defined as $g \nabla \rho/ \rho_0$ (Is this sign difference between M2 and N2 conventional?).
I checked the definitions of streamfunctions in Hallberg 2013 and I believe they are consistent with what is in in the rest of the literature (eg Ferrari et al 2010), but not with the code or the doc in MOM6.

[Hallberg-NOAA]

Thank you, @dhruvbalwada, for pointing out the sign error in the neutral_slope_x and neutral_slope_y and streamfunction diagnostics, and also for doing the extra work to determine that the internal sign conventions in the MOM6 code are at least self-consistent, if unconventional, and that the solutions themselves are not problematic.

A pull request to correct the sign error in the two slope diagnostics could be small (about 4-lines) and self-contained, but it would impact people who are depending on the diagnostic output not changing. This change should be discussed broadly within the community to see if anyone objects (e.g., if it would somehow disrupt operational uses of MOM6) or whether we should add a run-time flag to control the sign convention for the slope diagnostics.

A similar set of considerations would arise with the streamfunction diagnostics, and they could be handled separately, although I am not so certain that everyone always adopts the same sign convention, and I would not, for example, want to be accused of discrimination against left-handed people or residents of the Southern Hemisphere! Assuming that there is a consensus on the sign convention to use for streamfunctions, this should probably be handled separately from the correction to the slope signs to facilitate the review process.

As for the sign conventions for the internal variables in the code, those changes would be much more pervasive, but they should not lead to any answer changes, so they would not need as much deliberate coordination across the community. However, they should be dealt with in a separate pull request from the ones that address the signs of the diagnostics, in case we miss something and need to go back later to debug any problems that this conversion introduces (e.g., when synchronizing with code that is not yet on the main code branch).

If you would be willing to contribute one or more of the series of 3 pull requests to make the changes you are suggesting, that would be greatly appreciated.

[Hallberg-NOAA]

Upon further investigation (while finally getting around to putting together a PR to fix this issue), it turns out that the strange sign convention for slopes was only being used when the model was being run in pure layered mode and no equation of state was being used. Note that there is a sign difference between the case that uses the equation of state

MOM6/src/parameterizations/lateral/MOM_thickness_diffuse.F90

Line 1018 in e253883

Sfn_unlim_u(I,K) = -(KH_u(I,j,K)*G%dy_Cu(I,j))*Slope

and the one that does not.

MOM6/src/parameterizations/lateral/MOM_thickness_diffuse.F90

Line 1047 in e253883

Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope)

Both expressions ultimately give the right flattening effects on isopycnal slopes, but there is a bug in the sign of the diagnostics of the slope when no equation of state is being used. A similar issue exists with the slopes returned from calc_isonuetral_slopes(). Fortunately, it appears that apart from the interface height diffusion code in thickness_diffuse() and the diagnostics where this problem was originally noted, every other instance where these slopes are being used, it is only the absolute value or the square of the slopes that are being used.

Because this is ultimately a diagnostic sign error in the case that is not using an equation of state, I think that it is safe enough to just fix this problem without using a run-time flag to reproduce the previous bugs with the slope.

As for the sign convention for the overturning stream function (which is the same for all cases), I am not planning to change that now, although it would be easy enough to add a run-time flag flipping the sign of the conversion factor where the streamfunction diagnostics are registered with the lines

MOM6/src/parameterizations/lateral/MOM_thickness_diffuse.F90

Lines 2377 to 2388 in e253883

	CS%id_sfn_x = register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, Time, &
	'Parameterized Zonal Overturning Streamfunction', &
	'm3 s-1', conversion=GV%H_to_m*US%L_to_m**2*US%s_to_T)
	CS%id_sfn_y = register_diag_field('ocean_model', 'GM_sfn_y', diag%axesCvi, Time, &
	'Parameterized Meridional Overturning Streamfunction', &
	'm3 s-1', conversion=GV%H_to_m*US%L_to_m**2*US%s_to_T)
	CS%id_sfn_unlim_x = register_diag_field('ocean_model', 'GM_sfn_unlim_x', diag%axesCui, Time, &
	'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
	'm3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T)
	CS%id_sfn_unlim_y = register_diag_field('ocean_model', 'GM_sfn_unlim_y', diag%axesCvi, Time, &
	'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
	'm3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T)

at the end of thickness_diffuse_init().

[Hallberg-NOAA]

This issue has been addressed by MOM6 PR #633 .