CoMorph minor bug-fixes

[MichaelWhitall]

Version

main

Are there any linked Issues or Pull Requests?

No response

Brief description

This issue mops-up various minor bug-fixes to the CoMorph convection scheme, which have so-far been implemented in the UM comorph_dev branch UM comorph_dev branch for CoMorph B, but ought to be lodged so they can be included in UP.

Equivalent changes in the UM:
vn14.2_comorph_refact3 -> vn14.2_comorph_fixes1
Note the base-line for this diff is a branch which already includes the refactoring from Issues #178, #713 and #251 (which are not being lodged in the UM trunk).

Further details of the issue.

1. In comorph_constants_mod, we have increased the hardwired ad-hoc minimum limit on liquid-cloud droplet size r_min from 1 micron to 2 micron. This is used in comorph's in-plume microphysics scheme. The min limit is meant to be a crude representation of the CCN size (if the mass of liquid in the parcel is so small that the consistent droplet radius is smaller than this, the particle radius used for estimating condensation / evaporation rates, accretion etc defaults to the size of the dry particles that the water has condensed onto). To recap, comorph's microphysics scheme does not explicitly impose full adjustment to liquid-saturation; it is able to represent the small supersaturation occurring in rapid updrafts. Observationally, you get a peak in supersaturation slightly above cloud-base, due to the cloud-droplets being very small so that condensation is less efficient. As the droplets grow during ascent, the efficiency of condensation rapidly increases so that supersaturation quickly declines with height. What I found was that comorph tries to represent the peak in supersaturation just above cloud-base, but the model vertical grid is too coarse to resolve it properly. The updraft buoyancy at cloud-base (which largely controls the mass-flux at cloud-base via selective detrainment just below) is calculated by an interpolation which uses the parcel virtual temperature at the first grid-level above cloud-base. That in-turn was found to be sensitive to the in-plume supersaturation (higher supersaturation -> less condensation -> less buoyancy). And the vertical profile of supersaturation is inadequately resolved in the vertical, with a spurious dependency on whether the saturation height falls just above a grid-level or just beneath one. This caused occasional noisy / intermittent behaviour of comorph in shallow-Cu regimes when the saturation-height crossed a model-level boundary. Doubling r_min was sufficient to reduce the size of the poorly-resolved near-cloud-base supersaturation peak enough to avoid spurious noise in the mass-flux evolution.

2. In comorph_ctl, there is an initial basic convective triggering test (call to init_test) to determine which grid-points to do the full convective triggering calculations at (call to conv_genesis_ctl). In the existing code, init_test uses the fields at start-of-timestep (fields_n), whereas conv_genesis_ctl uses the "latest" fields already updated by other processes in the timestep (fields_np1). This inconsistency can occasionally cause convective triggering to spuriously fail, e.g. if the latest fields have liquid-cloud in a dry-statically-stable but moist-unstable environment, but the cloud hadn't appeared yet at start-of-timestep. We now correct this by passing fields_np1 into init_test, consistent with conv_genesis_ctl.

3. In conv_genesis_ctl, the surface boundary conditions passed into init_mass_moist_frac for the convective triggering calculations have been improved. The existing code passes in values of the primary fields at the current level (k), and also the levels above (k+1) and below (k-1) (inside init_mass_moist_frac -> calc_turb_parcel, the vertical gradients of various fields are used in the interpolation of turbulence fields from rho-levels to theta-levels). The values compressed onto potential convective triggering grid-points are held in a "multi-level compression array" fields_cmpr (and similar for turbulence fields and virtual temperature), which had vertical indices k_c = 0 for k, kp1_c = 1 for k+1, and km1_c = -1 for k-1. If at the top or bottom model-level, the existing code just set kp1_c or km1_c equal to k_c, so that the data passed in for the out-of-bounds grid-level just pointed to the values at level k. i.e. all variables were assumed to have zero vertical gradient between k=k_bot_conv and the surface (and between k=k_top_init and k_top_init+1). There were 2 problems with this:

   a) Assuming temperature and Tv have zero vertical gradient actually implies strong stable stratification.

   b) Applying a special condition at k=k_top_init was unnecessary, as we've already imposed elsewhere that k_top_init (the highest model-level where convection may trigger) must be less than k_top_conv (the highest model-level where convection maybe active).

   We address these problems. First, the confusing "slight-of-hand" with kp1_c and km1_c indices that point to index 0 (for level k) at the boundaries is removed. We now explicitly set the data in index -1 (for level k-1) to specify the lower boundary when k=k_bot_conv, under a new if-block. For now we assume winds go to zero at the surface (no-slip boundary), temperature and vapour extrapolate using the lapse-rate between k and k+1, and other fields have zero gradient as before. There is now no special treatment at k=k_top_init, since the fields are available at k_top_init+1 anyway. The special indices k_c, kp1_c, km1_c, kph_c, kmh_c for subscripting the multi-layer compression arrays are now deleted, and we just reference them using indices -1, 0 1 (as was already effectively hardwired in the declaration of the arrays).

4. The existing code in calc_turb_parcel (which interpolates the turbulent perturbations from rho-levels to theta-levels based on the local vertical gradients of scalars) had a trap to avoid using the gradients at the top and bottom model-levels (it defaulted to just using linear interpolation in height-space). Now that the above change sets the neighbouring model-level values of the fields appropriately even at the top and bottom model-levels, the if test on k not being at the top or bottom and the fall-back linear interpolation code are safely removed.

   I've also changed the method used to do the vertical interpolation of the turbulent w-perturbation in calc_turb_parcel; it is now based on a wind-sheer weighting rather than linear in height. This is designed to down-weight the contribution from the neighbouring level with the stronger wind-change relative to the current level. i.e. if the winds at the current level are similar to those above but very different to those below, then the current model-level is better-characterised by the properties above and so we assign its TKE closer to the value above. This often acts to down-weight the contribution from small surface values of the w perturbation when at k=1, when the near-surface wind profile is poorly-resolved in the vertical. This was needed once the changes in interp_turb (described further down this page) were added, yielding realistic smaller values of w_var at the surface. Too-small values of the w perturbation led to excessive parcel initial perturbations to T,q,u,v (since they are set as the turbulent flux divided by sqrt(w_var) ).

5. The existing routine par_gen_distinct_layers groups the found convective triggering parcels from all model-levels into a smaller number of distinct "layers". The existing code follows the simple rule that parcels triggering from adjacent model-levels must belong to the same layer, whereas if there is a gap with no convective triggering and then more parcels above it, those must belong to a separate layer. All parcels that have been grouped in the same layer are then integrated together in conv_sweep_ctl, but each layer is integrated independently. These grouping rules usually put convective parcels triggering from the surface dry-statically-unstable layer together in one layer, and any triggering from cloud at the mixed-layer top in a separate layer. However, sometimes (when the mixed-layer is shallow and poorly resolved in the vertical), there would be no dry-stable but cloud-free model-level between the surface unstable-layer and the mixed-layer-top cloud. When this happened, the surface-triggered and cloud-triggered convective parcels would be combined into a single layer. They often have quite different properties (parcel radius, MSE, etc), so taking a mass-flux-weighted mean over both and integrating them upwards as one can give a very different mass-flux profile compared to integrating the two separately. I spotted that in SCM runs of the TOGA-COARE case, the model would flip-flop between having a gap between the surface unstable-layer and the mixed-layer-top cloud and not having a gap. This led to spurious noisy jumps in the mass-flux profile.

   This issue has been addressed by tweaking the grouping rules in par_gen_distinct_layers; if one model-level contains liquid-cloud and another does not, their convective parcels are grouped into separate layers, even if they are vertically-adjacent. This forces the surface-triggered and cloud-triggered parcels to always be integrated separately, removing the noisy behaviour when the intervening gap model-level comes and goes. To facilitate this, the 3D liquid-cloud mixing-ratio fields % q_cl has been passed into par_gen_distinct_layers from conv_genesis_ctl.

6. In conv_incr_ctl, a check has been added to remove tiny negative values (for any positive-only fields), after adding the resolved-scale source terms from all convective plumes. It was found that, even though the source terms are limited elsewhere to avoid creating negative values, rounding-errors could still create tiny negatives where algebraically the result should have been zero. An existing check to remove negative values of the cloud-fractions is then removed (as that is now done by the new check on all positive-only fields). We also add a check to ensure the bulk cloud-fraction is within the plausible bounds consistent with the liquid and ice cloud-fractions (similarly the bounds could be violated by a tiny margin due to rounding-errors).

7. In mass_rearrange, a bug was found that caused very-rare but large violations of conservation for all transported scalars. This subroutine performs comorph's "compensating subsidence" calculation, by redistributing all the transported fields in the vertical so-as to get back the same dry-mass on each model-level as we started with, after entrainment and detrainment by the convection. It does this by searching up the column and trying to place the mass left after entrainment and detrainment into the next level, checking whether we have more or less mass than is required to fill the current model-level to its previous value, and moving to the next level when we have a surplus. Sometimes, due to rounding errors, we still have a tiny amount of surplus mass leftover after filling the uppermost convecting level. In this case, we just ditch the surplus mass (reset layer_mass_k to zero) and move on (there might still be another layer of convection higher-up in the column, and we don't want to spuriously move tiny bits of mass between the different layers). HOWEVER, when this happens the existing code wrongly sets the variable layer_mass_k2_added to zero as well. Occasionally, this can prompt later calculations in mass_rearrange_calc to spuriously add mass from a higher-up convection layer to model-level k2, and very-rarely it then fails to fill other model-levels, leaving spurious zeros in the fields (especially damaging for temperature!)

   This has been fixed by setting layer_mass_k2_added equal to the expected full mass of the layer (indicating there is nothing more to add here), and slightly adjusting the logic in mass_rearrange_calc to account for this.

8. The existing code in parcel_type_mod stores a field edge_virt_temp in the parcel properties array. This holds the virtual temperature at the outer edge of the in-plume assumed PDF. This was set equal to the environment virtual temperature after compensating subsidence, which is solved in the implicit detrainment calculation. This ensures the detrainment occurs so-as to keep the edge of the PDF neutrally buoyant. In the case where the whole plume remains buoyant (so that there is no detrainment), the imposed equality Tv_edge = Tv_env implies the PDF immediately stretches to keep its edge neutrally-buoyant (so that there is no discontinuity in Tv at the updraft edge). However, the detrainment calculation also imposes safety-limits on the positive or negative skewness of the assumed PDF, which sometimes makes it impossible to consistently match Tv_edge with Tv_env. I also found that sudden changes in the assumed PDF-shape when the detrainment rate falls to zero could create spurious noisy behaviour of the mass-flux profile.

   To avoid these problems, I've relaxed the assumption that Tv_edge = Tv_env. This is still the case when detrainment occurs, but in the absence of detrainment we now allow Tv_edge > Tv_env. When this happens, Tv_edge is relaxed smoothly towards Tv_env over a height-scale proportional to the parcel radius, instead of being forced to adjust immediately to Tv_env. This is done by the modifications in subroutine update_edge_virt_temp.

   To facilitate this, we store the new field "env_virt_temp" in the parcel super-array defined in parcel_type_mod alongside the existing parcel field "edge_virt_temp", and update both consistently throughout the plume-model. Some calculations that were using "edge_virt_temp" assuming it was equal to the env Tv are changed to use "env_virt_temp" instead (parcel_diags_type_mod).

9. In the parcel initiation calculations in init_mass_moist_frac, in set_par_fields we now store the grid-mean virtual temperature in the parcel field env_virt_temp instead of in the edge_virt_temp field as discussed above. Meanwhile edge_virt_temp can now be set independently in the different sub-grid regions (liquid-cloud, mixed-phase-cloud, ice / rain / graupel, clear-sky). This facilitates new options (to be added in later PR) to relax the assumption that all the sub-grid regions are neutrally-buoyant. In add_region_parcel, we now calculate the virtual temperature of the current sub-grid region (in the added call to calc_virt_temp using the unperturbed parcel fields "fields_par"). The values from the different regions are combined together and stored in the parcel "edge_virt_temp" field via the new call to core_combine from add_region_parcel. This reuses the existing code currently used to combine parcel core and edge properties from different model-levels, now to combine them from different sub-grid regions. This reduces needless code-duplication, but results in a change in the method used to combine the parcels from different regions (the old code took a mass-flux-weighted mean over the different region core properties, whereas the new code picks the region with the most buoyant core properties). I see the greater consistency here as an improvement (same method whatever the context for combining core properties), and it facilitates future options (to be added in a later PR) to combine the core properties using alternative methods. Note that subroutine core_combine needs as input the flag "l_down", indicating whether it is being called from a downdraft or an updraft, so this is now passed into add_region_parcel from init_mass_moist_frac.

10. In subroutine core_combine, we have fixed a bug that could spuriously set the combined parcel's core properties to those of a source parcel which had no mass-flux (e.g. if "parcel_m" had fully detrained before reaching level k to be combined with "parcel_a"). The amended code (which checks the mass-flux weights to determine whether each parcel has non-zero mass) should be safer. To do this, the mass-flux weights have been added to the argument list where core_combine is called from parcel_combine in parcel_type_mod.

11. In interp_turb (which interpolates the BL scheme turbulent fluxes and w-variance to rho-levels for input to comorph), we have changed the method used to estimate the turbulent w-variance at the surface. The existing code estimated it using the similarity theory:

    w_var = ( 1/4 w*^3 + u*^3 )^(2/3)

    where u* is the friction velocity, and w*^3 = boundary-layer-depth times surface buoyancy-flux is a convective velocity-scale. However, this was found to give much larger values at the surface than at k=1, k=2, etc, when in reality we expect w_var to increase with height in convective boundary-layers. I think this formula gives a scaling for the turbulent w-variance in the boundary-layer as a whole, not the reduced value expected near-surface.

    The spuriously enhanced value of w_var at the surface led to comorph's turbulent T,q perturbations being strangely reduced at the lowest model-level (when they generally decrease with height in the mixed-layer), since the perturbations scale with the fluxes divided by sqrt(w_var). This has been avoided by ditching the above equation for w_var and extrapolating from the values of w_var at k=1 and k=2 (with plausible limits). We still use the friction velocity^2 as a min limit, but the code rarely hits this limit.

    This change means the surface buoyancy flux "fb_surf" is no-longer used in comorph; we therefore delete it from the argument lists going from atmos_physics2 -> other_conv_ctl -> comorph_interface_um -> interp_turb.

12. In comorph_interface_um, the boundary-layer turbulent fluxes of heat, moisture and momentum are divided by density to get <w'Tl'>, <w'qt'>, <w'u'>, <w'v'>, as used in comorph. The existing code normalises the heat-flux "ftl" by wet density "rho_wet". However, it turns out then when the model is using mixing-ratios, the heat capacity of moisture should be neglected for consistency, so we should normalise the heat-flux by dry-density "rho_dry" instead. This has been corrected.

13. In the comorph microphysics, in ice_nucleation we adjust the homogeneous freezing of liquid to ice at the homogeneous freezing threshold (-40oC) to act more smoothly when the temperature is very close to the threshold. I found there were noisy jumps in the updraft buoyancy in the upper troposphere when the height at-which the parcel first falls to -40oC crosses a model-level. If the parcel is at -39.99 oC at level k at one timestep, no homog freezing occurs, then at the next timestep the parcel is at -40.01 oC at level k, all liquid is suddenly frozen at that level. This noisy behaviour is avoided by applying a correction to the homog freezing increment in ice_nucleation, so that instead of automatically freezing all liquid at -40oC, the freezing is limited to keep the temperature at or below -40oC after the latent heat of freezing is added on. This is done by first computing the temperature increment from freezing all liquid as before (dT_frz), then calculating a negative correction to the freezing increment dq_frz such that we scale it down to avoid heating beyond the homog freezing threshold:

    dq_frz -> dq_frz (T_hom - T_b4) / dT_frz
    = dq_frz - dq_frz (T_b4 - T_hom + dT_frz) / dT_frz

    (where T_hom is the homogeneous freezing threshold, and T_b4 is the temperature before homogeneous freezing was applied).

14. Subroutine calc_env_region_tq_nb calculates the temperature T and vapour content qv of the liquid-cloud, mixed-phase cloud, ice/rain and clear-sky sub-grid regions of the grid-box (these are used to attempt to trigger separate convective parcels from each region). For the liquid-cloud and mixed-phase cloud regions, T, qv are set such that the region is neutrally-buoyant and saturated w.r.t. liquid-water. The existing code estimates them by linearising the equations for virtual temperature Tv and saturation vapour mixing-ratio qsat about the grid-mean T and qv. However, I found that the linearisation of qsat sometimes introduced noticeable errors, so that the calculated T, qv of the liquid-cloud were slightly subsaturated. This occasionally spuriously / noisily suppresses convective triggering from liquid-cloud.

    This has been fixed by improving the accuracy of the saturation calculation for the liquid and mixed-phase cloud regions in calc_env_region_tq_nb. After the initial linear estimate, we call set_qsat_liq again to refine the saturated vapour mixing-ratio estimate at the sub-grid region's temperature. If the new qsat estimate differs sufficiently from the value at the grid-mean T, we recalculate dqsat/dT used in the linearisation based on a finite difference (this time the linearisation will exactly cross the actual qsat curve at the current-guess value of the region's T). Then we update our linearised estimate of saturated neutrally-buoyant qv using the new dqsat/dT (note the in-region temperature is updated consistent with this qv later in the subroutine, and we haven't changed that).

    Note the added calls to set_qsat_liq required passing "pressure" in through the argument list from calc_env_region_tq.

15. Subroutine calc_env_region_tq_nb also estimates a T, qv difference between the sub-grid region containing rain / graupel / ice but no liquid-cloud (labelled "icr") versus the region with no condensate (labelled "dry"). The existing code sets the "icr" region to be closer to saturation than the "dry" region by an amount equal to its total hydrometeor mixing-ratio (when supersaturated, we expect vapour-deposition onto the hydrometeors to've reduced the vapour-content of the "icr" region, whereas when subsaturated, we expect evaporation of precip to've increased the vapour-content of the "icr" region). The parameterised difference in vapour in the "icr" region was limited to avoid crossing saturation. But if it hit this limit (due to very high precip mixing-ratio), it made the "icr" region fully saturated. This led to a noisy behaviour where high rain water-contents occurred in grid-mean subsaturated, dry-statically-stable, moist-statically-unstable conditions. If the rain water content was just high enough, the "icr" region would fully saturate and so trigger convection (since the test lifting done in region_parcel_calcs would follow a moist adiabat). But if the rain water content was just a tiny but lower, the "icr" region would be slightly subsaturated and convection would not trigger. Further, rain evaporation alone should not be able to fully saturate the air (since the evaporation rate declines as you approach saturation), so parameterising saturated conditions in the "icr" region doesn't seem physical. I also thought this problem might be contributing to grid-point-storm activity, by allowing excessive convective triggering when the falling rain-mass gets large.

    To address these problems, the calculation of the "icr" versus "dry" region vapour difference in calc_env_region_tq_nb has been rewritten so that the "icr" region gets closer to saturation with increasing hydrometeor mass, but never quite reaches it (see the new comments added from L527 for the derivation of the new approach).

16. By historical accident, the calculations in cor_init_mass_liq_1 (which does an implicit correction to the initiating mass-flux from liquid-cloud) inconsistently use the gradient of qsat with temperature dqsat/dT calculated at the dry temperature T in some places, and calculated at the liquid-water temperature Tl in others (before CoMorph convection scheme refactoring #292, the values used in the calls to calc_qss_forcing_init were calculated on-the-fly inside that routine where Tl was not available; only when the calculation of dqsat/dT terms was rationalised in CoMorph convection scheme refactoring #292 did I spot the inconsistency). These calculations should most-correctly all consistently use dqsat/dT calculated at Tl. Therefore the calculation of dqsat/dT at T in init_mass_moist_frac has been deleted, it is no-longer passed into cor_init_mass_liq_1, and we now pass the existing value calculated at Tl into calc_qss_forcing_init in its place. The impact of this change is likely very small, but bigger than rounding-error.

17. In set_par_fields, the existing code applies the parcel initial RH perturbation (based on the namelist input "par_gen_rhpert") to both updrafts and downdrafts. But it seemed ill-justified to me to preferentially moisten downdraft initiating parcels, so the updated code only applies the RH perturbation to updrafts and sets it to zero for downdrafts. In practice this makes little difference, since in comorph downdrafts usually trigger from liquid-cloud in a moist-unstable environment (precip-driven downdrafts aren't properly represented yet). The RH perturbation gets limited to keep the initial parcel at or below liquid-saturation, so the RH perturbation already gets removed for downdrafts triggering from liquid-cloud.

18. The estimation of sub-level-step mass-fluxes based on the sub-level buoyancy profile in the CAPE calculation (calc_cape) has been deleted. We now use the more-accurate values of sub-level-step mass-fluxes calculated in the detrainment calculation, in set_det.

19. The calculation of the ratio of parcel core buoyancy over mean buoyancy in calc_core_mean_ratio (which sets the shape of the assumed PDF used in the detrainment) has been refactored to avoid noisy behaviour / spurious changes in detrainment when the buoyancies are extremely small or inconsistent. The added use of abs() and max() functions ensures that:

    • In the limit that the buoyancies are small enough to be comparable with numerical error of the virtual temperature, the ratio smoothly converges to a fall-back value of "par_gen_core_fac" (the ratio set for initiating parcels).
    • If the buoyancies go inconsistent such that really the ratio is ill-defined (e.g. core buoy smaller than mean buoy, or opposite signs), we extrapolate a smooth continuum response for the ratio, not a sudden noisy jump when a sign flips.

20. Subroutine calc_sat_height (called from parcel_dyn, the main routine that encapsulates moist and dynamic processes in the plume model) detects points where the plume has crossed from subsaturated to saturated (or vice-versa) and performs an interpolation to find the accurate cloud-base height. The existing version attempted to find the saturation height independently for both the parcel core and parcel mean properties. These then each had separate height addresses in the sublevs array, declared in conv_level_step ("i_sat" for the mean saturation height, and "i_core_sat" for the core saturation height). However, given that there's a continuous PDF of RH within the plume, the height where the mean happens to cross saturation didn't really match any discontinuity relevant to the detrainment calculation (further, it was ill-defined at the point where it was calculated in calc_sat_height, since the selective detrainment subsequently changes the parcel mean RH and therefore the mean saturation height). In recognition of this, there was already some slightly convoluted code in calc_sat_height to reset the mean saturation height properties to those of the core (by setting "i_sat" = "i_core_sat") if the core hit saturation before the mean did (which was usually the case!)

    Therefore, we have simplified this area by deleting the separate variable "i_core_sat" from conv_level_step, parcel_dyn, and now just store the core saturation height in the address for "i_sat". We now only do the interpolation to find the accurate saturation height in the call to calc_sat_height for the parcel core, and skip it in the call for the parcel mean (hence the loop to find cloud-base-straddling points is now protected by if (.not. l_mean_with_core) then). This allowed the complicated bit of code which modified the saturation height interpolation for the mean parcel properties to be deleted (along with the associated compression list indices "index_ic_new"; some later calculations now need to be defined on the existing "index_ic_sat" indices instead).

    Then towards the end of calc_sat_height there was a block of code to interpolate the core saturation height properties to the mean saturation height and vice-versa; this has been replaced by a new calculation of the parcel mean buoyancy at the single saturation height (which now corresponds to the core saturation height). This has been improved; in the usual case where the parcel core has reached saturation first (hence no condensation has yet occurred in the parcel-mean properties), we find the parcel mean buoyancy at the core saturation height by interpolating the parcel virtual temperature in the absence of condensation ("prev_tvl", "next_tvl") to that height. The existing code simply interpolated the mean buoyancy profile from the sub-level heights above and below, allowing increased buoyancy above the saturation height to contribute to the buoyancy at the saturation height; this was a bug.

21. At the end of calc_sat_height, we have added a new block of code to calculate a new output "next_sat_buoy"; the buoyancy at the PDF-position corresponding to the saturation boundary at the end-of-level-step height (the diff viewer has confusingly spliced this with parts of the deleted core sat height / mean sat height interpolation discussed above). The new code handles the situation where, at end-of-level-step, the parcel core properties are saturated but part of the in-plume PDF is still subsaturated. In this case, we expect there to be a buoyancy minimum at the saturation boundary within the plume. The existing detrainment code underestimates the width of the in-plume buoyancy PDF when this happens, since it does not account for the PDF having a minimum somewhere in the middle, rather than at the edge). See the added in-line comments for the full derivation of the calculation of "next_sat_buoy". This variable is then passed out through parcel_dyn to conv_level_step and then into set_det, where it is used in new code to set a minimum limit on the width of the assumed PDF used for detrainment...

22. In set_det (the main detrainment routine), we temporarily modify the core buoyancies held in the "sublevs" super-array before they are used to compute the detrainment rate, then reset them back again (the original core buoyancy is held in "tmp_buoy" so we can copy it back into "sublevs" afterwards). This is to impose a minimum limit on the width of the assumed PDF used in the detrainment calculation. This avoids problems with noisy behaviour / erratic detrainment rates when the mean and core buoyancies are nearly equal, so that the assumed in-plume buoyancy PDF would otherwise be extremely narrow. For simplicity, comorph assumes all in-plume variables are perfectly correlated with eachother in the PDF; the lack of decorrelated scatter means the spread of buoyancies spuriously collapses to zero if the core and mean parcel buoyancies cross over. Imposing a minimum limit on the buoyancy width is an attempt to represent the residual decorrelated scatter, giving smoother and more robust behaviour when the core versus mean parcel properties collapse in buoyancy space.

    The minimum buoyancy width is parameterised as:

    a) A dimensionless constant "min_width_fac" times the width of the in-plume supersaturation PDF (in terms of vapour mixing-ratio converted to its virtual effect on buoyancy, by scaling by dTv/dqv = Tv (Rv/Rd - 1).

    b) The difference between the parcel mean buoyancy and the buoyancy at the saturation boundary within the PDF ("next_sat_buoy", calculated in calc_sat_height as discussed above).

    c) A fall-back tiny numerical tolerance "safety_thresh" * Tv, which forces the width to be at least a few times greater than the expected rounding error of the virtual temperatures differenced to calculate buoyancy.

    If the core minus mean buoyancy width is smaller than the minimum limit "min_width", we reset the core buoyancy to mean buoyancy + "min_width" to broaden the PDF, and store the fraction by-which we have broadened the PDF in "cmm_frac_lev" / "cmm_frac" (where the latter stores the fraction at the sub-level step where the most detrainment occurred). "cmm_frac" is then used to scale down the difference between the detrained versus non-detrained fields due to selective detrainment (so that in the limit that the original buoyancy PDF width had collapsed to zero, the difference vanishes, since the buoyancy PDF is no-longer represented at all by the difference between the core versus mean properties of the parcel).

    The PDF-width broadening factor "cmm_frac" is also passed out of set_det to conv_level_step and then into the aforementioned update_edge_virt_temp, where the rate of relaxation of the parcel edge Tv towards the environment Tv is scaled down if we have broadened the PDF. This reduces the tendency to get noisy behaviour with very narrow PDFs due to small changes in edge Tv yielding very rapid changes in "core_mean_ratio" (which sets PDF shape / skewness) when the core and mean buoyancies are too close together. i.e. if the PDF is broken, relaxing the edge Tv is likely to break it even worse, so leave it alone.

23. Where the entrained mass is used in the calculation of the update to parcel radius (subroutine update_par_radius), we now use the entrained mass calculated before the application of the CFL limit, not the final limited entrained mass (the CFL limit is a purely numerical constraint on the entrainment and wasn't meant to also limit the radius). Growth of the mass-flux with height is assumed to lead to growth of the thermals and hence increase of radius with height. The change here avoids a potential problem where, when the mass-flux gets huge, the updraft radius fails to increase with height following the stated equation, so that entrainment rates (which scale with 1/R) remain too high aloft. To implement this change, a new variable "ent_mass_orig" is declared in conv_level_step, passed out from set_ent (where it is set equal to "ent_mass_d" before the CFL-limit is applied), and passed into update_par_radius. Here, those options that use the ratio of next over prev mass-flux now estimate the next mass-flux we would've had without the CFL-limiting of the entrainment (actual next mass-flux + the missing entrained mass ent_mass_orig-ent_mass, scaled down by the non-detrained fraction of the total mass-flux).