Bulkmicro: avoid rain sedimentation when qr ≤ 0

[adoyenne]

Summary

Add an early return in calc_sed_qr_sb to ensure zero rain sedimentation when no rain water is present.

Details

The sedimentation rate of rain water is proportional to the rain water content (qr).
When qr ≤ 0, sedimentation should therefore be zero by definition.

This PR adds a short-circuit at the top of the elemental function to:

• Explicitly return sed_qr = 0 for qr ≤ 0
• Avoid unnecessary and potentially ill-defined calculations of derived rain properties

Impact

• No change for grid points with rain present
• Improved numerical robustness and clarity for dry conditions
• Safe for elemental, vectorized, and OpenACC execution

[fjansson]

Is this for speed or correctness?

There should already be a limiter in place, clipping negative rain values to 0. One thing to watch out for is that early returns may interfere with statistics, not sure if that applies here.

EDIT: I see the early return is safe for statistics (this is in a per-gridpoint function).

[adoyenne]

Is this for speed or correctness?

There should already be a limiter in place, clipping negative rain values to 0. One thing to watch out for is that early returns may interfere with statistics, not sure if that applies here.

EDIT: I see the early return is safe for statistics (this is in a per-gridpoint function).

Good to know. I added this return simply because I was getting a division-by-zero error when running in debug mode. Debug mode is usually more picky, and after adding this procedure here, the error stopped occurring.

[fjansson]

Ok, if the debug mode catches a division by 0 there, we're doing something wrong. Do you happen to know where exactly it occurred? I can imagine that if qr or Nr is 0 (which is allowed), some of those derived quantities D, mu, are not defined. As I remember it many of them also have limiters to deal with this, perhaps not all.

[adoyenne]

Ok, if the debug mode catches a division by 0 there, we're doing something wrong. Do you happen to know where exactly it occurred? I can imagine that if qr or Nr is 0 (which is allowed), some of those derived quantities D, mu, are not defined. As I remember it many of them also have limiters to deal with this, perhaps not all.

I checked again with the up-to-date dev version, and it no longer complains about this issue (at least with my current setup). So it was probably fixed elsewhere, and this short-circuit may no longer be necessary.

However, when running in debug mode I still encountered two crashes that do not appear in the normal run.

The first one occurred in the new lD80R subgrid scheme:

[tcn688:2508765:0:2508765] Caught signal 8 (Floating point exception: floating-point divide by zero)
==== backtrace (tid:2508739) ====
0 0x000000000003ebf0 __GI___sigaction()
1 0x0000000000e78786 __modsubgrid_MOD_closure()
/gpfs/home5/alos/New_DALES/dales_dev_new/src/modsubgrid.f90:418

I addressed it by adding guards against very small values:

             if (dthvdz(i,j,k) > 0 .and. zh(k) > eps1 .and. e120(i,j,k) > eps1) then

                 stab = sqrt(grav/thvf(k)*dthvdz(i,j,k))
                 mixl = (0.4_field_r*zh(k) * cn*e120(i,j,k)/stab) &
                        / (0.4_field_r*zh(k) + cn*e120(i,j,k)/stab)

                 zlt(i,j,k) = min(delta(k), mixl)
              else
                 zlt(i,j,k) = delta(k)
              end if  

The second crash appeared in modlsm.f90:

[tcn790:2231322:0:2231322] Caught signal 8 (Floating point exception: floating-point divide by zero)
==== backtrace (tid:2231322) ====
0 0x000000000003ebf0 __GI___sigaction()
1 0x000000000095e800 __modlsm_MOD_calc_obuk_dirichlet()
/gpfs/home5/alos/New_DALES/dales_dev_new/src/modlsm.f90:3112 ( this line: L = L - fx0/fxdif)

I handled this by protecting the Newton update from a zero derivative:

        if (abs(fxdif) > 1e-12) then
            L = L - fx0/fxdif
        else
            exit   ! derivative unusable → leave iteration
        end if

After these fixes, the debug run completes without crashes.