Add cloud fraction diagnostics

[trontrytel]

This PR adds cloud fraction diagnostics based on diagnosed covariances.

[szy21]

I didn't look at the quad_loop, but the rest looks good, thanks! I may look at it if I have more time, but feel free to merge it.

[Sbozzolo]

Note that since you are changing the type/struct of the cache, this will be a breaking change.

[trontrytel]

Note that since you are changing the type/struct of the cache, this will be a breaking change.

It's because I'm renaming env_thermo_quad int SG_quad?

[Sbozzolo]

Note that since you are changing the type/struct of the cache, this will be a breaking change.

It's because I'm renaming env_thermo_quad int SG_quad?

You are also changing its type to be a subtype of AbstractSGSamplingType.

I don't think anyone downstream uses this, but if the cache is supposed to be a public interface, this will change it.

[akshaysridhar]

Thanks @trontrytel. Some small comments, otherwise looks reasonable!

[trontrytel]

Note that since you are changing the type/struct of the cache, this will be a breaking change.

It's because I'm renaming env_thermo_quad int SG_quad?

You are also changing its type to be a subtype of AbstractSGSamplingType.

I don't think anyone downstream uses this, but if the cache is supposed to be a public interface, this will change it.

I'll consider if I stand behind my name changes enough to do the release, or revert them

[Sbozzolo]

Note that since you are changing the type/struct of the cache, this will be a breaking change.

It's because I'm renaming env_thermo_quad int SG_quad?

You are also changing its type to be a subtype of AbstractSGSamplingType.
I don't think anyone downstream uses this, but if the cache is supposed to be a public interface, this will change it.

I'll consider if I stand behind my name changes enough to do the release, or revert them

You shouldn't feel constrained by this. Also, in the past people changed the cache all the time.

[trontrytel]

For computing the Smagorinsky length scale I changed from

        l_smag =
            c_smag *
            ᶜdz *
            max(0, 1 - N_eff^2 / ᶜPr / (2 * ᶜstrain_rate_norm))^(1 / 4)

to

N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ?
           c_smag * dz * (1 - N_eff^2 / Pr / 2 / ϵ_st)^(1 / 4) : c_smag * dz 

Otherwise the max was giving me zeros almost all the time and was making the whole length scale equal to zero. Does that look reasonable @tapios , @szy21 ?

Additionally I had to change the EnvBuoyGrad struct to pass in also the liquid and ice cloud condensate for the NonEquilibrium moisture option. I think I can do it more elegantly, but maybe in the next PR.

The changes in allocations are small. So if the above two issues are fine, I would bump the version, bump the allocations and merge it in.

[szy21]

For computing the Smagorinsky length scale I changed from

        l_smag =
            c_smag *
            ᶜdz *
            max(0, 1 - N_eff^2 / ᶜPr / (2 * ᶜstrain_rate_norm))^(1 / 4)

to

N_eff > FT(0) && N_eff < sqrt(2 * Pr * ϵ_st) ?
           c_smag * dz * (1 - N_eff^2 / Pr / 2 / ϵ_st)^(1 / 4) : c_smag * dz 

Otherwise the max was giving me zeros almost all the time and was making the whole length scale equal to zero. Does that look reasonable @tapios , @szy21 ?

I think it makes sense that mixing length goes to zero when buoyancy gradient is larger than strain rate, in which case it will reduce to one quadrature point (or we can do many quadrature points with the same value, as the variance is zero). It is common for mixing length to be zero during the spin-up, as we start with zero velocity and the shear is zero. I don't fully understand why this term is still mostly zero after e.g. 6 hours of bomex, which is why I haven't use this as a lower limit for mixing length. But I think we can figure out that later.

So, in short, I think it would be better to keep the check for N_eff > FT(0) but remove N_eff < sqrt(2 * Pr * ϵ_st) for now, and make sure the cloud diagnostics work with zero mixing length.

[trontrytel]

So, in short, I think it would be better to keep the check for N_eff > FT(0) but remove N_eff < sqrt(2 * Pr * ϵ_st) for now, and make sure the cloud diagnostics work with zero mixing length.

It seems weird to have this logic. I would at least want to make things smooth between the different limits

N_eff < 0 => l_smag = c_smag * dz

0 < N_eff < sqrt(2 * Pr * ϵ_st) => l_smag = 0

N_eff > sqrt(2 * Pr * ϵ_st) => c_smag * dz * (1 - N_eff^2 / Pr / 2 / ϵ_st)^(1 / 4)

The couple of cases I tried worked with zero (co)variances, just produced no cloud fraction. Thats why I was looking for ways to make the estimated (co)variances larger. But maybe you are right, that the CI does not provide meaningful conditions to really test it

[tapios]

So, in short, I think it would be better to keep the check for N_eff > FT(0) but remove N_eff < sqrt(2 * Pr * ϵ_st) for now, and make sure the cloud diagnostics work with zero mixing length.

It seems weird to have this logic. I would at least want to make things smooth between the different limits

N_eff < 0 => l_smag = c_smag * dz

0 < N_eff < sqrt(2 * Pr * ϵ_st) => l_smag = 0

N_eff > sqrt(2 * Pr * ϵ_st) => c_smag * dz * (1 - N_eff^2 / Pr / 2 / ϵ_st)^(1 / 4)

The couple of cases I tried worked with zero (co)variances, just produced no cloud fraction. Thats why I was looking for ways to make the estimated (co)variances larger. But maybe you are right, that the CI does not provide meaningful conditions to really test it

The original logic (truncating c_smag * dz by a stratification factor) is fine for vertical mixing. It's possible we need to think about fluctuations generated by horizontal SGS motions (i.e., similar logic with horizontal gradients and a horizontal mixing length).

When you do not get clouds at all, that means the mixing length is also zero near the top of the boundary layer?

[szy21]

It seems weird to have this logic. I would at least want to make things smooth between the different limits

N_eff < 0 => l_smag = c_smag * dz

0 < N_eff < sqrt(2 * Pr * ϵ_st) => l_smag = 0

N_eff > sqrt(2 * Pr * ϵ_st) => c_smag * dz * (1 - N_eff^2 / Pr / 2 / strain_rate)^(1 / 4)

I think the logic is actually:
N_eff < 0 => l_smag = c_smag * dz
0 < N_eff < sqrt(2 * Pr * ϵ_st) => c_smag * dz * (1 - N_eff^2 / Pr / 2 / strain_rate)^(1 / 4)
N_eff > sqrt(2 * Pr * ϵ_st) => l_smag = 0
(i.e. switch the second and third lines). Isn't this smooth?

[trontrytel]

I think the logic is actually: N_eff < 0 => l_smag = c_smag * dz 0 < N_eff < sqrt(2 * Pr * ϵ_st) => c_smag * dz * (1 - N_eff^2 / Pr / 2 / strain_rate)^(1 / 4) N_eff > sqrt(2 * Pr * ϵ_st) => l_smag = 0 (i.e. switch the second and third lines). Isn't this smooth?

You are right. This makes more sense!

[trontrytel]

When you do not get clouds at all, that means the mixing length is also zero near the top of the boundary layer?

I will add mixing length to the output tomorrow and see how it looks in the spherical plots from longer simulations.

Below are some quick plots after 2 days of sphere_held_suarez_rhoe_equilmoist_topography_dcmip from our CI. Thanks to @Sbozzolo new diagnostics and analysis utilities.

[tapios]

Looks pretty reasonable! When you plot relative humidity, could you do it in some such way that RH > 1 on the grid scale is clearly recognizable? I'd like to see where you actually have grid-scale supersaturation.

[trontrytel]

Here is the Smagorinsky mixing length for the same simulation. Empty plot means that its all zero

[szy21]

Thanks!