Possible zero in the denominator in MOM_mixed_layer_restrat

[gustavo-marques]

We have a global test case where the corner point is on the Equator. Once we set GUST_CONST=0.0, following the discussion about default parameter changes in one of the dev/calls, this test started failing. The reason being a zero in the denominator of the following calculation:

mom_mixrate = (0.41*9.8696)*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)

See relevant lines of code here

There is no issue when GUST_CONST > 0 because that enforces u_star to be > 0. However, since the corner point is on the Equator absf is problematic and perhaps must be bounded:

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))

@alperaltuntas has come up with a possible way of avoiding this issue, but it would be good to get @Hallberg-NOAA's option before we send a PR.

[Hallberg-NOAA]

At first glance this is a good candidate for using Adcroft's rule of reciprocals, namely if (x==0) then Inv_x = 0.0 else Inv_x = 1./x. Alternately we could just trap the special case where absf = 0 in the mom_mixrate expression on this line, in which a mathematical cancellation of terms gives if (absf == 0) mom_mixrate = (0.41*9.8696) * u_star / (4.0*(h_vel+dz_neglect)), but this would just lead to problems a with an infinite mixing growth timescale a few lines later. The code subsequently applies limits on the restratifying streamfunctions to avoid any CFL violations, so a very large (but not infinite) mixing growth timescale might give plausible answers.

Another approach might be to take the viscous limit on the mixing rate as floor, namely something like mom_mixrate_min = pi2 * Kv / (h_vel + dz_neglect)**2, where pi2 ~= 9.8696 and Kv is a (molecular or turbulent?) viscosity. This would probably still give such long timescales that the CFL limits would still be what really determined the streamfunctions in the end, but at least there would be a physical explanation for this timescale.

[gustavo-marques]

@alperaltuntas, can you please post here your initial fix to this issue so we can compare all the alternatives?

[Hallberg-NOAA]

@gustavo-marques and @alperaltuntas, do you have an update on this issue?

[alperaltuntas]

@Hallberg-NOAA, we haven't implemented a fix for this yet. IIRC, the solution I temporarily tested was to skip the computations on land cells, although I am not sure which mask2d is an appropriate test here. Again, this issue occurs only when ustar==0.0, so our current temporary fix is to run with small GUST_CONST when we have to( i.e., when any grid cell corners coincide with the Equator). Also, I am happy to test any of the solutions you proposed. if I interpret your first solution proposal correctly, the fix is to (effectively) set mom_mixrate and timescale to zero when ustar==0.0 and absf==0. If so, I can implement and test it.

[Hallberg-NOAA]

This issue has been addressed by https://github.com/NOAA-GFDL/MOM6/pull/266, which has now been merged onto the main branch of MOM6, and hence this issue can be closed.