Obukhov length calculation fails

[fjansson]

I'm trying to run a simulation with zero-flux boundary conditions at the ground. This seems to be possible, by choosing flux boundary conditions (isurf=4) and setting the temperature and qt fluxes to 0. However, this causes the Obukhov length calculation to fail, before the first time step has been completed.

To reproduce:
Fresh clone of dales master branch, compiled in DEBUG mode.
cases/example
in namoptions.001, set wtsurf=0 and wqsurf=0.
isurf=4 is already set, meaning flux boundary conditions at the ground.
Running this results in:

.....
&NAMCAPE
 LCAPE=F,
 DTAV=  60.000000000000000     ,
 /
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
STOP Obukhov length calculation does not converge!

The stop occurs in subroutine getobl(), in modsurface.f90, line 1307 (i.e. the second convergence check).
The reason may be that on line 1280, Rib gets set to 0 since thv = thvs, which seems to happen when the boundary fluxes are set to 0.
The case Rib = 0 is then not handled in the if:s that follow on 1285, should it be?

[fjansson]

Some more details on the Obukhov length.
Rib in the code is the bulk Richardson number, eq (29) in Heus et al 2010. If the theta_l gradient at the ground is 0, Rib is also 0.

Rib = 0 corresponds to an infinite Obukhov length L. According to Wikipedia, an infinite L is expected at dawn and dusk, when the ground heat flux passes through 0.

There is already a cap for L in the code: if (abs(L)>1e6) L = sign(1.0e6,L) (modsurface.f90, line 1316).
I propose to set L=1e6 also if Rib == 0 (or perhaps if abs(Rib) is smaller than some threshold)

As a workaround, I'm specifying a tiny but non-zero flux for now.

[fjansson]

Fixed in version 4.2, commit cd92f73