E3SM-Project
SHOC: unintentional discontinuous functions
In discussing EKAT issue 117, we came across two instances where divide-by-0 avoidance and clipping code lead to a discontinuous function:
Example 1:
relvar(:,:) = 1.0_r8
relvarmax = 10.0_r8
where (rcm(:ncol,:pver) /= 0.0 .and. rcm2(:ncol,:pver) /= 0.0) &
relvar(:ncol,:pver) = min(relvarmax,max(0.001_r8,rcm(:ncol,:pver)**2.0/rcm2(:ncol,:pver)))
For simplicity, assume these are all scalars, not arrays. Suppose rcm2 = 0 and rcm = 1. Then rcm/rcm2 = +Inf. In that case, one would think relvar should be set to relvarmax = 10. But, instead, the above code makes relvar = 1. Yet if rcm2 were just a bit above 0, say rcm2 = 1e-6, the above code would indeed make it 10. So that makes the overall function discontinuous. A similar situation occurs when rcm = 0 but rcm2 != 0. relvar is set to 1 but if rcm were just a bit above 0, it would instead be 0.001.
Example 2:
testvar=(a*sqrtqw2_1*sqrtthl2_1+(1._rtype-a)*sqrtqw2_2*sqrtthl2_2)
if (testvar .eq. 0._rtype) then
r_qwthl_1=0._rtype
else
r_qwthl_1=max(-1.0_rtype,min(1.0_rtype,(qwthlsec-a*(qw1_1-qw_first) &
*(thl1_1-thl_first)-(1._rtype-a)*(qw1_2-qw_first) &
*(thl1_2-thl_first))/testvar))
endif
If testvar = 0, then r_qwthl_1 = 0. But if it deviates from 0 by a tiny bit, then r_qwthl_1 = 1 or -1.
A related question: if, say in example 1, rcm2=0, will the result of relvar end up being actually used later on, or is whatever value we put in going to be ignored? If it's the latter, then the not-so-pretty discontinuity is benign (though it might quickly turn malign if for some reason we change the code and end up using the result, without eliminating the discontinuity).
Luca, not sure about example 1. The output in example 2 appears to be used without qualification elsewhere in SHOC.
Either way, I agree that a continuous treatment of the clipping is preferable. If the clipping is one-sided, then it should be as easy as init-ing the var to the clip value (no additional cost). If two sided, it might get trickier, I'm not sure.
@ambrad and/or @bogensch - naive question, but could this discontinuity in the limiter contribute to cold T biases? I'm guessing it is an edge case, but worth checking anyway.
Seems unlikely. But we want to remove discontinuities caused by the limiters regardless.
