Non-reproducible trigonometric functions in viscous BBL

[marshallward]

There are trigonometric functions in the calculation of the bottom boundary layer viscosity, set_viscous_BBL(), of the form cos(acos(x)/3 - 2*pi/3).

These functions may to be responsible for an answer change when transitioning from an Intel to AMD processor:

@@ -298,11 +298,11 @@
 h-point: mean=   2.5651029259221975E+01 min=   0.0000000000000000E+00 max=   4.0983512515433851E+01 Start set_viscous_BBL S
 h-point: c= 937092751 sw= 934896663 se= 934896663 nw= 939288839 ne= 939288839 Start set_viscous_BBL S
 u-point: mean=   4.7372083769301478E-07 min=   0.0000000000000000E+00 max=   5.1746049326678756E-02 u Ray [uv]
-u-point: c= 575868012 u Ray [uv]
+u-point: c= 575868009 u Ray [uv]
 v-point: mean=   4.5227642461442761E-07 min=   0.0000000000000000E+00 max=   1.1190647973926611E-02 v Ray [uv]
-v-point: c= 574147674 v Ray [uv]
+v-point: c= 574147672 v Ray [uv]
 u-point: mean=   9.2146868327139069E-04 min=   0.0000000000000000E+00 max=   1.5000000000000007E-03 u kv_bbl_[uv]

Compilers were identical across machines. AFAIK, Intel does not use the C math library to compute these functions.

CHANNEL_DRAG = False restores answers, hinting that these functions may be responsible for the change in answers. In any case, they are recognized as a potential source of answer changes.

We may want to consider replacing these expressions with a bit-reproducible approximation.

Snippet containing the lines is shown below.

MOM6/src/parameterizations/vertical/MOM_set_viscosity.F90

Lines 939 to 963 in 25b57f4

	elseif (crv > 0) then
	! There may be a minimum depth, and there are
	! analytic expressions for L for all cases.
	if (vol < Vol_2_reg) then
	! In this case, there is a contiguous open region and
	! vol = 0.5*L^2*(slope + crv/3*(3-4L)).
	if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7?
	! There is a very good approximation here for massless layers.
	L0 = sqrt(2.0*vol*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0)
	else
	L(K) = apb_4a * (1.0 - &
	2.0 * cos(C1_3*acos(a2x48_apb3*vol - 1.0) - C2pi_3))
	endif
	! To check the answers.
	! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol
	else ! There are two separate open regions.
	! vol = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L)
	! At the deepest volume, L = slope/crv, at the top L = 1.
	!L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol)) - C2pi_3)
	tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol)
	tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
	L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3)
	! To check the answers.
	! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol
	endif

[Hallberg-NOAA]

This code is directly solving for the root of a cubic equation using trigonometric functions. These roots are well bracketed and well-behaved, so we could replace this code snippet with an iterative Newton's method solver that should of comparable cost and accuracy, but we would need to preserve this code wrapped in an answer-date flag because it has been used for several years in published solutions.

[Hallberg-NOAA]

The new options for solving the cubic equations for the open face lengths that were introduced with #543 allow the user to avoid using these trigonometric functions by setting TRIG_CHANNEL_DRAG_WIDTHS. The fact that the new options are both more accurate and more efficient should be motivation for them to be used routinely, but the existing solver that uses trigonometric functions does need to be retained to avoid changing answers for existing solutions. This issue has been as thoroughly addressed as it can be, so it should now be closed.