Prevent numerical instability in pt_downscale_radiations

[mmazzolini]

Hi everyone,

I have experience some numerical instability when downscaling the Shortwave Radiation. I am downscaling ERA5 to a set of representative points in a basin. Some of them have large slope and I think that's why I encountered this problem.
I noticed the problem as in some points I got nonphysical values of SW, caused by a large SW_direct component. Here's an example of a west-facing and 30 degrees steep point:

I get a spike of SW radiation in the evening (sorry for the wrong time axis, times are in UTC but the basin is at UTC-7). That is caused by the fact that the ratio between cos_illumination and mu0 explodes, likely because there is a very small value of mu0 just before the sunset.

I thought of solving it by setting a minimum value of mu0 as 0.25 mu0_stable = mu0.where(mu0 > 0.25, 0.25). That value corresponds to a sun zenith of 75 degrees. I called this mu0_stable and substituted it in the formulas where mu0 appears in the denominator.

Moreover, while I was exploring these terms I also discovered that, for some approximation errors I guess, the calculation of surf_interp['SW_direct'] leads in some cases to values slightly larger than the SW radiation on the top of the atmosphere (SWtoa). When that happens, I got that ka was at times negative, as you can see in the next figure (representing another west-facing point but steeper). Negative ka also leads to exploding SW_direct, this time in the formula for down_pt['SW_direct_tmp'].

So I removed the values of surf_interp['SW_direct'] that are larger than SWtoa by substituting them with SWtoa:

surf_interp['SW_direct'][surf_interp['SW_direct'] > SWtoa] = SWtoa[surf_interp['SW_direct'] > SWtoa]

I've already discussed about the first issue with @krisaalstad. However, I'm not sure at all that these fixes work well in other cases. I guess there needs to be some testing in other basins before accepting these?
If needed I can provide you with the dataset I'm working with.

Cheers,
Marco.

[ArcticSnow]

Thank you @mmazzolini ! Good catch.

I won't have time soon to verify this non physical spikes occur in for instance the the example watersheds. Or in case you have time on your hand, I'll let you check too with the example cases: https://github.com/ArcticSnow/TopoPyScale_examples/tree/main

Also, we now have another downscaling workflow where this would need to be corrected too: https://github.com/ArcticSnow/TopoPyScale/blob/main/TopoPyScale/topo_scale_zarr.py I put this for when we follow through the PR to not forget to correct the bug here as well, if consistent.

This points to the need to put in place some thorough testing scripts, looking at hourly behavior of the output variables.