P3: `D_gr` is less than `D_th` for some `F_rim`, `ρ_rim`

[haakon-e]

It appears that for some values of F_rim and ρ_rim, the expected relationship D_th < D_gr < D_cr does not hold.

In particular, F_rim=0.93, ρ_rim=975.0 gives rise to D_th > D_gr. Code follows:

1. Using built-in functionality (needs rft: remove P3Distribution struct #582):

import CloudMicrophysics.Parameters as CMP
import CloudMicrophysics.P3Scheme as P3
FT = Float64
params = CMP.ParametersP3(FT)
F_rim = FT(0.93)
ρ_rim = FT(975)
D_th1, D_gr1, = P3.get_thresholds_ρ_g(params, F_rim, ρ_rim)
# 9.73e-5 9.69e-5

2. using previous iterative methods

import RootSolvers as RS

FT = Float64
F_rim = FT(0.93)
ρ_rim = FT(975)

# Constants from CMP.ParametersP3
α_va, β_va = 0.018537721864540644, 1.9
ρ_i = 916.7

function ρ_d_helper(D_cr, D_gr)
   β_m2 = β_va - 2
   6 * α_va * (D_cr^β_m2 - D_gr^β_m2) / π / β_m2 / (D_cr - D_gr)
end
ρ_g_helper(F_rim, ρ_rim, ρ_d) = F_rim * ρ_rim + (1 - F_rim) * ρ_d
D_gr_helper(ρ_g) = (6 * α_va / π / ρ_g)^(1 / (3 - β_va))
D_cr_helper(F_rim, ρ_g) = (1 / (1 - F_rim) * 6 * α_va / π / ρ_g)^(1 / (3 - β_va))

P3_problem(ρ_d) =
   ρ_d - ρ_d_helper(
	   D_cr_helper(F_rim, ρ_g_helper(F_rim, ρ_rim, ρ_d)),
	   D_gr_helper(ρ_g_helper(F_rim, ρ_rim, ρ_d)),
   )

ρ_d =
   RS.find_zero(
	   P3_problem,
	   RS.SecantMethod(FT(0), FT(1000)),
	   RS.CompactSolution(),
   ).root
# 193.4

ρ_g = ρ_g_helper(F_rim, ρ_rim, ρ_d)
# 920.3

D_th = (π * ρ_i / 6 / α_va)^(1 / (β_va - 3))
# 9.73e-5

D_gr_helper(ρ_g)
# 9.69e-5

[trontrytel]

Just to leave some written trail of our conversation: One option would be to assume rime density to be that of the frozen water. I don't know enough about hail etc, so I don't know what processes would make it more dense than fully frozen liquid. If this is physical behavior on the other hand, then we could just swap the order in the integrals, depending on which threshold ends up being smaller.

[haakon-e]

From MM15:

The predicted rime density ρ_r is also limited to values between 50 and 900 kg / m^3. If necessary, B_rim is adjusted to keep ρ_r within this specified range.

So we could think of some reasonable ways to achieve this, possibly in the host code.