Found a bug in routine back_to_cell_average and a bug in routine p3_main_part3 in P3 microphysics SCREAM F90 version

[linlintamu]

Here we report two bugs found in the P3 microphysics SCREAM F90 version code.
1- in routine back_to_cell_average, the one line code that converts in-cloud immersion freezing of cloud drops mass into grid cell mean value: qc2qi_hetero_freeze_tend = qc2qi_hetero_freeze_tend***il_cldm** ! Immersion freezing of cloud drops.

The cloud fraction that is multiplied by should be liquid cloud fraction "cld_frac_l" rather than the intersection of liquid and ice cloud "il_cldm".

The corresponding conversion of in-cloud immersion freezing of cloud drop number into grid cell mean value: nc2ni_immers_freeze_tend = nc2ni_immers_freeze_tend*cld_frac_l ! Number change associated with freexzing of cld drops use the correct cloud fraction to multiply by.

2- in routine p3_main_part3, when diagnosing the rain effective radius: diag_eff_radius_qr(k) = 1.5_rtype/lamr(k).
Now it is hard-coded using constant 1.5 divided by lamr(k). This is correct only when the spectral width parameter (mu_r) in the rain gamma size distribution is set to zero. However, the model now use mu_r = mu_r_constant set in routine get_rain_dsd2. And mu_r_constant = 1 set in micro_p3_utils.F90. As a result, the rain effective radius is no longer 1.5/lamr(k). To be highly flexible in accepting any given mu_r values and automate the diagnosis of rain effective radius, I suggest the equation should be diag_eff_radius_qr(k) = 0.5_rtype*(mu_r(k)+3._rtype)/lamr(k). This equation gives you diag_eff_radius_qr(k) = 2_rtype/lamr(k) when mu_r(k) = 1.

[PeterCaldwell]

@linlintamu and @hassanbeydoun - item 1 above seems to say that immersion freezing should be applied anywhere there's liquid cloud even if there's no ice cloud but the fix Hassan put in #2901 makes immersion freezing occur only where there's both liquid and ice cloud. My understanding of immersion freezing is that it's caused by supercooled liquid coming in contact with aerosols conducive to freezing... in which case it can happen outside of ice cloud. This might be important to enable ice initiation in a supercooled liquid cloud. What am I missing? In any case, I think we need a comment in the code changed by the PR explaining why we make the choice we do.

[hassanbeydoun]

@PeterCaldwell that's a good point. In the current code, qc2qi_hetero_freeze_tend is multiplied by il_cldm while nc2ni_immers_freeze_tend is multiplied by cld_frac_l so that's inconsistent. They should be multiplied by the same thing since they are the same process. The question is should they be multiplied by cld_frac_l or il_cldm? I think it should be cld_frac_l as you say. I will make the fix so qc2qi_hetero_freeze_tend is multiplied by cld_frac_l.

[PeterCaldwell]

Yup, agreed there really is a bug because we were inconsistent between ni and qi treatments. And yeah, I think cld_frac_l is the way to go for both cases (but I could be wrong).

[PeterCaldwell]

I agree with Lin's 2nd bug but add a few notes here because it took me a while to confirm... Eq 5 of MG1 (https://journals.ametsoc.org/view/journals/clim/21/15/2008jcli2105.1.xml ) says effective radius = gamma(mu+4)/gamma(mu+3) * 1/(2*lambda) where mu and lambda are parameters of the gamma PDF. And barely remembered rules for the gamma function yield that gamma(mu+4)/gamma(mu+3) = mu+3. So yes, it all works out.