Stability of turbulent flux calculations

[RichardEssery]

Some large spikes occur in turbulent heat flux calculations. Although these spikes do not appear to strongly affect snow simulations, numerical stability could be improved by implementation of an iterative solution.

[GiuliaMazzotti]

Review especially:

• Parametrizations of canopy heat capacity
• Computation of bulk Richardson number for CANMOD = 0

[GiuliaMazzotti]

A follow-up comment potentially related to the same issue:

In the subroutine 'THERMAL.F90', see:

! Surface layer
do j = 1, Ny
do i = 1, Nx
Dz1(i,j) = max(Dzsoil(1), Ds(1,i,j))
Ts1(i,j) = Tsoil(1,i,j) + (Tsnow(1,i,j) - Tsoil(1,i,j))Ds(1,i,j)/Dzsoil(1)
ksurf(i,j) = Dzsoil(1) / (2Ds(1,i,j)/ksnow(1,i,j) + (Dzsoil(1) - 2Ds(1,i,j))/ksoil(1,i,j))
if (Ds(1,i,j) > 0.5Dzsoil(1)) ksurf(i,j) = ksnow(1,i,j)
if (Ds(1,i,j) > Dzsoil(1)) Ts1(i,j) = Tsnow(1,i,j)
end do
end do

Three points are unclear to us:

1. Ds(1,i,j) > Dzsoil(1) is never fulfilled if the default parameter settings are used (Dzsoil(1) = 0.1m, Ds(1) >= 0.1m
2. How would a change from default parameters affect these expressions (e.vg. setting Dzsoil(1) = 0.2m)
3. Tsurf is declared in the header of the subroutine but never used

[RichardEssery]

1. The snow layering is too complicated to be fully described in the code comments. Ds(1,i,j) > Dzsoil(1) never occurs if Nsnow(i,j) > 1, but a single snow layer will grow to thickness 2*Dzsnow(1) before dividing.

2. The numerics for shallow snow rely on a trick that requires Dzsnow(1) = Dzsoil(1). This is not entirely satisfactory, as we might well want to use thinner snow layers without changing the soil layers.

3. You are right that Tsurf is not used in THERMAL.f90. I will remove it.

[RichardEssery]

Iterative solution has been added