Undefined behaviour in `calculate_Tin_from_qin`

[kieranricardo]

Hi there,

I've come across some undefined behaviour in calculate_Tin_from_qin. Depending on the compiler and compilation options, when the sqrt argument is negative the following will either yield a NaN or Tmltk. Replacing min with ieee_min_num fixes this for my use case, and will consistently return Tmltk when the sqrt argument is negative, independent of the compiler.

https://github.com/CICE-Consortium/Icepack/blob/f6ff8f7c4d4cb6feabe3651b13204cf43fc948e3/columnphysics/icepack_therm_shared.F90#L80C1-L81

Tin =  min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) /  &
                         (c2*aa1),Tmltk)

Happy to submit a PR from my fork https://github.com/kieranricardo/Icepack/tree/bugfix-calculate_Tin_from_qin-undefined-behaviour :)

[eclare108213]

The quantity under the square root probably should not ever be negative. What are the values of aa1, bb1, cc1 when this fails? Is it a roundoff problem or a bigger issue that needs to be looked into?

[dabail10]

Any updates on this one?

[anton-seaice]

Kieran is on leave this week, we resolved why the quantity was negative.

We could just ignore this, or if we think returning Tmltk is better than returning Nan when the first argument is Nan, then make the change?

[dabail10]

I am happy to put in a fix for this. Maybe we should just put in a check for the argument in the square root and set it to zero when negative. Then if the result is greater than Tmltk, then you just get that. I suspect this is a roundoff issue as highlighted by @eclare108213 .

[anton-seaice]

@apcraig - Do you have any experience with ieee_min_num and if it is better than min? I am not sure about compiler support for it?

[dabail10]

I did some looking into this and ieee_min_num is a fortran 2018 feature, and so we need to be careful here. I think the typical approach here is to do the following instead:

Tin =  min((-bb1 - sqrt(max(bb1*bb1 - c4*aa1*cc1),puny)) /  &
                         (c2*aa1),Tmltk)

[anton-seaice]

@kieranricardo resolved why the quantity was negative some time ago, we think it was an issue with how the coupling was getting set-up.

My view is we shouldn't change anything, because we don't think its valid for the quantity under the square root to be negative, and giving an error helps with debugging.

I guess a more graceful handling would be to check if bb1*bb1 - c4*aa1*cc1 is negative, and abort with a meaningful error if it is.