Protect against zero divisions in RH for cold T

[trontrytel]

We are getting NaNs in RH for very low temperatures because our saturation vapor pressure is zero. For example

TD.saturation_vapor_pressure(tps, FT(30), TD.Ice()) = 0
TD.saturation_vapor_pressure(tps, FT(40), TD.Ice()) = 0
TD.saturation_vapor_pressure(tps, FT(50), TD.Ice()) = 4.16f-42

This PR adds a small number in the denominator to prevent that.

eps(Float32(0)) = 1.0f-45

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 90.40%. Comparing base (f6fc4c8) to head (7610265).
Report is 2 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #230   +/-   ##
=======================================
  Coverage   90.40%   90.40%           
=======================================
  Files           9        9           
  Lines        1167     1167           
=======================================
  Hits         1055     1055           
  Misses        112      112           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[szy21]

Could we also limit relative humidity between 0 and 1, and add a test for saturation_vapor_pressure at cold temperature?

[trontrytel]

Could we also limit relative humidity between 0 and 1, and add a test for saturation_vapor_pressure at cold temperature?

I can add the test. I thought that Tapio wanted to avoid actually clipping the RH.

[charleskawczynski]

This looks fine to me. I'm just wondering if we know or can know what inputs could lead to this? Is it dominated by temperature?

[szy21]

Could we also limit relative humidity between 0 and 1, and add a test for saturation_vapor_pressure at cold temperature?

I can add the test. I thought that Tapio wanted to avoid actually clipping the RH.

I think @tapios wants to clip RH to 1 (because this won't prevent that)? Let's double check.

[trontrytel]

This looks fine to me. I'm just wondering if we know or can know what inputs could lead to this? Is it dominated by temperature?

Yes, for T < 40K our parameterization for saturation vapor pressure returns zero Hence the nans

[trontrytel]

Lets wait for me to add the test

[tapios]

Could we also limit relative humidity between 0 and 1, and add a test for saturation_vapor_pressure at cold temperature?

I can add the test. I thought that Tapio wanted to avoid actually clipping the RH.

I think @tapios wants to clip RH to 1 (because this won't prevent that)? Let's double check.

I think clipping to <= 1 makes sense, provided (a) we document this clearly in the doc string, and (b) there is no case where we use this for microphysics and assessing supersaturation. @trontrytel Is there any reason we'd want to diagnose supersaturation here?

[trontrytel]

I think clipping to <= 1 makes sense, provided (a) we document this clearly in the doc string, and (b) there is no case where we use this for microphysics and assessing supersaturation. @trontrytel Is there any reason we'd want to diagnose supersaturation here?

Sounds good. And no, I was just thinking that for debugging Atmos it might be easier to see the actual RH numbers we get from the division than to mask them with clipping. But we can always add saturation_excess to diagnostics if we need it.

[szy21]

Thanks! Just one more thing. Should we add the eps to relative humidity, or should we add it to saturation_vapor_pressure itself? I don't have a strong preference. The latter can prevent issues if we use saturation_vapor_pressure elsewhere, but we can also remember to limit it if we use it in the denominator.

[trontrytel]

I'm voting for adding it to RH only. The only reason why we need to do anything is because we are dividing by it here

[trontrytel]

I'm going to try to merge it in, so that our nightly tiny AMIP can run with it. We can always revisit tomorrow, if we want to change something about the epsilon!

[szy21]

Great. This issue may not show up in the nightly AMIP run now. I'll run one with 55km top later tonight, which will be more likely to have the issue.