CoMorph miscellaneous bug-fixes

[MichaelWhitall]

Version

main

Are there any linked Issues or Pull Requests?

No response

Brief description

This issue mops-up numerous 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).

See also the diff between this branch and the comorph_dev branch:
vn14.2_comorph_fixes1 -> vn14.2_comorph_dev
(which should now only contain changes to add new functionality as all refactoring and hardwired bug-fixes should be included in vn14.2_comorph_fixes1).

Note to self: checked this diff as far as parcel_dyn.F90:

• Found some missed code that moves some q_rain back to q_cl in the parcel mean-fields with core, to account for reduced autoconversion.

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 similar for the parcel diagnostics in parcel_diags_type_mod), 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). New code is added towards the end of conv_level_step to consistently set "env_virt_temp" equal to the value after compensating subsidence has been applied (by adding the value of "delta_tv" held in the "sublevs" super-array onto the "virt_temp" value held in the "env_next_super" array, and storing the result in the pacel super-array).

9. We remove an existing assumption that the compensating subsidence increment to env Tv is constant over each level-step; instead the thetav gradient term is allowed to vary linearly with height over each level-step, and the scaling by mass-flux to get "delta_tv" (held in the "sublevs" array) uses the actual sub-level-step mass-fluxes found by the detrainment calculation. These changes improve the accuracy / consistency of the subsidence increment (used in the implicit detrainment and other calculations) at inversions that are poorly-resolved in the vertical, giving a smoother evolution of mass-flux over time. Several changes were required to do this:

   • Subroutine calc_delta_tv (called from conv_sweep_ctl) has been refactored / simplified. The existing version had separate calculations for the 1st and 2nd halves of the current level-step, such that, when ascending from the previous model-level interface to the full-level k it used the Tv gradient between k and k-1, whereas when ascending from full-level k to the next model-level interface it switches discontinuously to using the Tv gradient between k+1 and k (this caused noisy jumps in detrainment at cloud-base when the buoyancy minimum at the saturation height crossed the height of full-level k). The new version now interpolates between the subsidence heating expected at full-level k (which depends on the gradient between k+1 and k), and the subsidence heating expected at the upper model-level interface, half-way between the level k heating rate and that at full-level k+1 (which depends on the gradient between k+2 and k+1). What is set in calc_delta_tv and compressed onto grid-points for the current convection type/layer in conv_sweep_compress is now explicitly the value at the end of the current half-level-step. As such, we have changed the name of the variable in conv_sweep_ctl, conv_sweep_compress, conv_level_step, init_sublevs from the ambiguous "delta_tv" to "delta_tv_next". In conv_sweep_ctl and calc_delta_tv we also carry a variable called "delta_tv_prev", which is just used to interpolate "delta_tv_next" to the next model-level interface.

   • The TEMPORARY CODE in init_sublevs which set "delta_tv" at the start of the level-step ("prev") equal to the value at the end ("next") is deleted. As described above, at the end of each level-step we now set "env_virt_temp" held in the parcel super-array to the current env Tv + "delta_tv" (the subsidence increment). Then at the start of the next timestep, in init_sublevs we now retrieve the previous value of "delta_tv" by exactly reversing that calculation. This fixes a bug in the existing code; due to an inconsistency, the env Tv after compensating subsidence and "delta_tv" values used at the start of one level-step did not exactly match the values at the end of the previous level-step. This led to a spurious drift of the buoyancies and "core_mean_ratio" with height, affecting the detrainment profile and also the entrainment of environment properties into the core (which is parameterised as a function of "core_mean_ratio").

   • In the existing code, the update of the environment Tv increment due to compensating subsidence to account for the fractional detrainment reducing the mass-flux is done in conv_level_step, just after the call to set_det. The original version assumes "delta_tv" on all sub-level steps is scaled down by a constant ratio, based on the fractional detrainment at the end of the level-step. We now delete this calculation (which is under a TEMPORARY CODE comment) and use a more accurate version inside set_det; "delta_tv" on each sub-level is now scaled down by the non-detrained fraction of the mass-flux at that sub-level, instead of assuming the scaling is constant over the level-step. This

   • Towards the end of conv_level_step, some TEMPORARY CODE to update the mean buoyancy and mass-flux at the end of the level-step in the "sublevs" array has been deleted, replaced by more accurate calculations inside the refactored set_det routine.

10. 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.

11. 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.

12. 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.

13. 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.

14. 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).

15. 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.

16. 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).

17. 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.

18. 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.

19. 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.

20. 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.

21. 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, init_sublevs, 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.

22. 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". The calculation is done in the call to parcel_dyn -> calc_sat_height for the parcel-mean properties, but requires as additional inputs the supersaturation and virtual temperature the parcel core would have with all liquid evaporated ("next_core_ss", "next_core_tvl"). Therefore these existing variables defined in conv_level_step are passed into parcel_dyn -> calc_sat_height as new optional arguments (only used in the call for the mean parcel properties). Once "next_sat_buoy" is calculated, it is then passed out through parcel_dyn to conv_level_step as a new optional output argument, and then passed 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...

23. 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 if all liquid-cloud were evaporated (in terms of vapour mixing-ratio converted to its virtual effect on buoyancy, by scaling by dTv/dqv = Tv (Rv/Rd - 1) ). This requires passing the supersaturation of the parcel mean and core properties (if all liquid-cloud were evaporated) from conv_level_step into set_det (the new arguments "par_prev_mean_ss", "par_next_mean_ss", "par_prev_core_ss", "par_next_core_ss").

    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 (note this makes an existing check on core buoy > mean buoy by at least a numerical tolerance redundant, so that block and its associated if-test and fall-back approximate calculations are deleted). We 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.

24. We have added a new subroutine swap_core_edge, called from set_det just before and after the main detrainment calculation. This is to handle the case where the parcel-core is less-buoyant than the parcel-mean". The existing code just ignored the core buoyancy when this happened, and assumed a uniform buoyancy equal to the parce-mean. However, if the parcel mean-fields are still buoyant but the parcel core has become negatively-buoyant, we ought really to selectively-detrain mass from the core end of the PDF (the existing code cannot do this). We now add the ability to do this, by flipping the PDF around. When the core less less buoyant than the mean, subroutine swap_core_edge recalculates the parcel core properties, setting them equal to the extrapolated values at the existing outer edge of the PDF (so that the new extrapolated edge equals the existing core). "core_mean_ratio" is recalculated consistent with this swapped arrangement. The rest of the detrainment calculation then proceeds as before, now detraining mass from the core end of the PDF. Afterwards, swap_core_edge is called again to flip back to having the original (least dilute) properties defined in the core variables. The logical "l_swap_core_edge" flags the point where the PDF has been swapped around, and is used to swap only the same points back again.

25. In the detrainment calculation (set_det), we now compute the detrained versus non-detrained qw = q_vap + q_cl, instead of doing this for q_cl on its own. This turned out to be essential to allow us to have an assumed in-plume PDF where q_cl is positive in the core but falls to zero part-way through the distribution (i.e. have a saturation boundary within the PDF). The existing code checks that any fields which are positive-only do not go negative at the edge of the PDF (which would potentially make them negative in the detrained air). If they do, it makes the PDF narrower (by reducing the core value used in the calculation) so-as to keep the edge value >= zero. When this check was imposed on q_cl, it forced the detrained air to always have positive q_cl (and hence liquid cloud-fraction) if the core and mean did. I think this was particularly a problem when using cloud-shell downdrafts (which are not lodged yet and not included in this PR), but might as well include the fix here.

    We avoid this problem by adding q_vap onto q_cl (in the mean and core properties) in set_det before the detrainment calculations, performing the calculations and positivity checks on their sum, then subtracting q_vap off again afterwards (in the mean, core and detrained properties). In the event that q_cl is then negative after the subtraction, we evaporate the negative value to yield zero q_cl, and conservatively update q_vap and temperature.

26. The loop over sub-level steps to compute the detrainment in set_det has been refactored. The existing code looped over the sub-levels backwards (do i_lev = i_next_max, i_prev+1, -1). This was originally done for computational efficiency; you most-often got the highest detrained fraction at the last sub-level "i_next", so we computed that one first. On other sub-level-steps, it was cheap to pre-estimate whether more of the buoyancy PDF overlapped with zero than at sub-levels that were already calculated, so the full computations at other sub-levels with lower total detrainment than "i_next" could be skipped. However, various problems were found where fields on sub-levels were used (e.g. in the calculation of convective cloud, or in the w-equation) but had not been correctly updated due to the calculations being skipped. Also the cost-saving from skipping those calculations was minimal, since on the majority of model-levels you only have sub-levels "i_prev" and "i_next" and you only evaluated detrainment at "i_next" anyway (the additional intervening sub-level "i_sat" is only used when crossing cloud-base).

    Therefore, we now loop over sub-levels forwards (do i_lev = i_prev+1, i_next_max), and refactor the setting of sub-level fields (the "sublevs" array) to ensure they are all updated correctly at all sub-levels. Note we have moved various things to inside the loop over sub-levels to simplify the code, with just one big loop instead of multiple smaller ones:

    • The code to reverse the sign of the buoyancies for downdrafts before calculating detrainment and then change them back again.
    • Updating of "sublevs" mass-flux and delta_tv fields, by scaling them down by the non-detrained fraction "frac" (this has allowed the array "frac_sublevs", which saved the non-detrained fraction on all sub-levels, to be deleted).
    • The call to wind_w_eqn which updates the parcel vertical velocities using the latest implicitly-solved buoyancies.

    The check_bad_values calls on "frac" and "x_edge" are then retained after all the above, just after the loop over sub-level-steps.

27. Another bug we fix in set_det; while the parcel mean properties get updated to account for selective detrainment, I found that the mean buoyancies held in the sublevs array were accidentally missing this update. It had been added in a commented-out state under a TEMPORARY CODE comment in CoMorph convection scheme refactoring #292 so-as to preserve KGO; we now delete that and add the intended code to update "sublevs(ic,j_mean_buoy,i_lev)" with the selective detrainment increment (now additionally accounting for "cmm_frac" as discussed earlier). The lack of selective detrainment update to mean_buoy meant it was strongly underestimated where detrainment rates were high, sometimes even negative. This field is used to estimate a vertical velocity in the convective cloud-fraction calculation; the bug tended to cause spurious large convective cloud fraction values near cloud-base.

28. 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).

29. In momentum_eqn, we fix a bug in the convective momentum transport. There is a quadratic drag term which relaxes the parcel winds towards the environment winds. This was incorrectly implemented as being proportional to the square of the difference in only the current wind component:

    du(i)/dt ~ -( u_par(i) - u_env(i) )^2

    (where i denotes the zonal, meridional or vertical component).

    This formula is not invariant to the wind orientation (the overall drag is reduced when the wind is diagonal rather than straight along a grid-line); this can't be right. The correct formula for quadratic drag in 3D is:

    du(i)/dt ~ -abs( u_par - u_env ) ( u_par(i) - u_env(i) )

    = sqrt( (u_par_x - u_env_x)^2 + (u_par_y - u_env_y)^2 + (u_par_z - u_env_z)^2 ) ( u_par(i) - u_env(i) )

    The wind vector magnitude term abs( u_par - u_env ) was already calculated and stored in the variable "wind_ex", where it was correctly used already in the drag on the vertical velocity. We just needed to use it in the horizontal drag too, in place of the term ( u_par(i) - u_env(i) ) in "alpha".

30. In precip_res_source (which sets the resolved-scale source terms due to precip falling out of the parcel into the environment), we add some extra code to evaporate liquid-cloud passed from parcel to environment by "droplet settling". In comorph, all condensed water species have a fall-speed; the calculated fall-flux of liquid-cloud in the parcel is generally small but not zero. This leads to a tiny flux of q_cl being passed to the environment without an accompanying cloud-fraction increment, which can cause problems. We now avoid this by assuming that liquid-cloud droplets evaporate immediately when passed to the environment. This change is protected by the new logical "l_evap_q_cl_flux" in comorph_constants_mod, but this is currently hardwired to true.

31. In conv_level_step, we add a missing call to set_par_cloudfrac to ensure the detrained cloud-fractions are updated consistent with the detrained condensate (in the array "det_fields"). "det_fields" are always calculated at the end of the current half-level-step. When stepping from level k to the next interface, "det_fields" have to be adjusted back to level k before being detrained; this is done by the final call to parcel_dyn, acting on "det_fields", and this includes a call to set_par_cloudfrac inside parcel_dyn to set the cloud-fractions consistently. But when stepping from the previous interface to level k, "det_fields" don't need to be adjusted to level k (as they're already defined there) so the call to parcel_dyn is omitted (protected by "if ( .not. l_to_full_level )"). However, this meant the consistent setting of the "det_fields" cloud-fractions within parcel_dyn was also inadvertently omitted. We now add this in via the new call to set_par_cloudfrac, protected by "if ( l_to_full_level )".

32. In set_ent, where we apply the CFL-limit on the entrained mass, we now account for environment mass already removed by the initiation mass-source. i.e. now the sum of init mas + ent mass must be less than the mass held on the model-level, which is safer (the existing code already imposed this constraint via a CFL-limit imposed via the closure rescaling; now that we impose it locally within the plume-model we should be rescaling the profile to reduce the closure less often).

    To implement this, we've added code in conv_sweep_ctl and conv_sweep_compress to compress the initiating mass-source onto the current convection type / layer's compression list, in "sum_massinit_cmpr". This is then passed into conv_level_step -> set_ent where it is used to modify the max allowed entrainment "max_ent" for the CFL constraint.