CoMorph KGO-breaking refactoring

[MichaelWhitall]

Are there any linked Issues or Pull Requests?

No response

Brief description

The changes under Issue #178 / PR #292 and Issue #713 / PR #717 left various lines of "temporary code" in place so-as to avoid changing known-good-output under the wholesale refactoring of the CoMorph convection scheme. This included many lines that do nothing scientifically, but are just hacks to force calculations to yield the same floating-point rounding-errors as before. Now we remove these, with no expected science impact, but a KGO-change in CoMorph due to changing rounding-errors in various calculations.

lfric_apps branch:
vn3.2_comorph_refact2 -> vn3.2_comorph_refact3
Note the base-line for this diff is the branch for #717, which needs to be lodged before this.

Equivalent changes in the UM:
vn14.2_comorph_refact2 -> vn14.2_comorph_refact3
Note the base-line code for this diff is a branch which already includes the KGO-preserving refactoring from PRs #292 and #717 (which are not being lodged in the UM trunk).

Further details of the issue.

Main design changes:

1. The code to impose a height-varying minimum limit on the turbulence length-scale (which then scales the parcel radius), using ass_min_radius and min_radius_fac, has been moved from set_par_fields.F90 (within comorph itself) to calc_turb_len.F90 (where the turbulence length-scale is calculated for input to comorph). ass_min_radius and min_radius_fac are then no-longer used within comorph itself so are deleted from comorph_constants_mod. Also the variable par_radius_amp_um (the "precip ramp" factor scaling the parcel radius, calculated in calc_turb_len and passed into comorph) now includes the scaling by par_radius_knob, so that par_radius_knob no-longer needs to be applied inside comorph itself (hence the scaling of par_gen_radius_fac by par_radius_knob in set_constants_from_um.F90 is deleted). These changes simplify the code's handling of parcel radius scaling without changing the science (but changes answers at the level of floating-point rounding error, as the order of calculations has changed). Note the code in the upgrade macro (versions.py) to divide the existing value of ass_min_radius by par_radius_knob; this is needed to preserve existing behaviour after moving the application of the min radius limit to before the point where the scaling by par_radius_knob comes in (now via par_radius_amp). We also now apply the updated CoMA9 standard values for the "min_radius" in the CoMorph stand-alone test, in set_test_profiles (with par_radius_knob=0.45, upper value of the limit = 500.0/0.45 = 1111.1, and height-dependence factor = 0.25/0.45 = 0.55555).

2. Related to the above, we also add the tuning factors used to compute the turbulence length-scale (and the height-dependence of its minimum limit) to the namelist, as turb_len_fac (which replaces the confusingly named and confusingly rescaled par_gen_radius_fac) and min_radius_fac. See the changes in the rose meta-data, upgrade macro, um_physics_init_mod, comorph_um_namelist_mod, set_constants_from_um.

3. The handling of the model's lower and upper boundaries is simplified in the calculation of the subsidence Tv increment for a given mass-flux (delta_tv in conv_sweep_ctl, used in the implicit detrainment calculation). Instead of having some complicated code in calc_delta_tv to handle the special case where the "next" model-level doesn't exist, we just skip the call to calc_delta_tv entirely at the last model-level (k=1 for downdrafts) and just leave delta_tv set to zero in this case. This basically makes no difference, since downdrafts that reach k=1 are forced to fully detrain (in set_det), regardless of delta_tv. But this change doesn't quite preserve KGO (it may alter how much detrainment of downdrafts occurs in the 1st versus 2nd half-level step within k=1, introducing rounding differences in the sum of resolved-scale source terms over the two half-levels. Note this change was already implemented in conv_sweep_ctl but commented-out, with a TEMPORARY CODE comment which we can now remove.

4. In env_half_mod, subroutine env_half_interp is responsible for interpolating the environment virtual temperature and vertical velocity onto model-level interfaces (aka rho-levels). Here, some code intended to treat the special case of the model's lower and upper boundaries (where either k-1 or k+1 do not exist to interpolate from) was wrongly active for the inwards interpolation (between k and k+1 at k=1), not just the outwards interpolation that would go out-of-bounds. This is fixed by replacing the if-test on l_last_level with an explicit test on whether the "next" k index will go out-of-bounds. We can then remove l_last_level from the argument lists in the calls from conv_sweep_ctl -> conv_sweep_compress -> env_half_interp since it is no-longer used in these routines. As with the previous change, this will have no science impact, since downdrafts reaching k=1 are forced to fully detrain anyway (and updrafts don't reach the model-top unless something has gone wrong).

5. The vertical interpolation of virtual temperature to model-level interfaces (aka "half-levels", "rho-levels") is refactored. Instead of doing this on-the-fly in env_half_mod subroutine env_half_interp, we add a new bespoke subroutine interp_virt_temp which does the interpolation for Tv in a full 3D array at a higher level in the code. For now this just does linear interpolation similar to the existing calculation in env_half_interp, but will facilitate a not-yet-included change to add an option to use a higher-order vertical interpolation, giving more accurate treatment at inversions. Subroutine interp_virt_temp converts Tv to potential temperature, then interpolates, then converts back to Tv. This improves the accuracy and means we assume neutral stratification at the top and bottom of the model (where we assume the interpolated variable has same value as the adjacent full level) instead of assuming constant Tv which implies strong stable stratification. As noted in the above points, this makes no real difference since downdrafts must fully detrain in the bottom model-level anyway. The new 3D array "virt_temp_half" is declared in comorph_ctl where interp_virt_temp is called. It is then passed down through the call tree argument lists comorph_main -> conv_sweep_ctl, -> conv_sweep_compress -> env_half_interp. Subroutine env_half_interp is then refactored to replace its existing interpolation of env Tv with just compressing the value from the 3D array "virt_temp_half" instead (the existing vertical interpolation of env w in env_half_interp is retained).

6. The handling of parcel properties stored at the BL-top (used for homogenising convective source-terms within the BL) is simplified. This affects subroutines save_parcel_bl_top, homog_conv_bl_ctl, homog_conv_bl, which handle the saved parcel properties in par_bl_top. The original version of save_parcel_bl_top tests whether the current theta-level k (which the parcel has now crossed) is within the BL, and whether the previous theta-level (k-1 for updrafts, k+1 for downdrafts) is. And then it also tested whether the mass-flux at the end of level k is still > zero. However, the code is simpler and more intuitive if we just test whether theta-level k is the level spanning the BL-top and always save the parcel properties at this point (i.e. if level k-1 is fully within the BL and level k isn't). This removes the "TEMPORARY CODE" block which tests on mass-flux > zero in save_parcel_bl_top. We then adjust the k-indexing in homog_conv_bl_ctl, homog_conv_bl to account for this change, allowing us to remove a "TEMPORARY CODE" block there too.

7. The above paved the way for new functionality to homogenise partially / smoothly across the BL-top-spanning model-level in homog_conv_bl, instead of having a nasty jump in behaviour of the homogenisation when the BL-top height crosses a model-level boundary. To do this, the saved properties of the parcel at the BL-top are first adiabatically adjusted from the rho-level where they're defined to the BL-top-spanning theta-level (the added call to fields_k_pressure_adjust). The accurate BL-top height z_bl_top is now passed in from homog_conv_bl_ctl. Then we interpolate to estimate conserved variable properties at the accurate BL-top height, by adding a fraction of the convective source-terms at the spanning theta-level onto the parcel (from conservation, d/dz of mass-flux * parcel properties = -env source-terms). This is in the added calculation of "homog_frac" etc. Finally, the part of the existing convective source-terms in the BL-top-spanning theta-level that is above the accurate BL-top height is protected from the homogenisation (only scaled down by 1-homog_frac at this point instead of reset to zero). NOTE: the new vertically-smooth homogenisation is not actually active yet because on the trunk the BL-top height passed into comorph from LFRic conv_comorph_kernel_mod (UM comorph_interface_um), zh_homog, is rounded up to the next whole rho-level height so-as to preserve existing behaviour. A later PR will add the option to let this vary smoothly.

8. In homog_conv_bl, mass-weighted vertical means of various fields are computed (used to estimate a vertically-uniform perturbation to apply in a vertically-uniform initiating mass-source within the boundary-layer). The existing code first computes a mass-weighted vertical integral of each field within the boundary-layer, then normalises by the vertical integral of mass afterwards to retrieve a mean. The new code is simplified by dividing the mass of each model-level (stored in layer_mass_cmpr) by its vertical integral a-priori, removing the need for the various extra normalisation calculations. We also safely move the rescaling of layer_mass_cmpr to convert it to the initiating mass on each level (such that its vertical integral over the BL equals the mass-flux at the BL-top) from old L435 to later at new L841. Before, the vertical integrals were confusingly mass-weighted using original layer-mass for some variables, but the proportional initiating mass for others.

9. Calls to dry_adiabat now update temperature directly, instead of updating virtual temperature or adjusting an "exner_ratio" (initialised to 1.0) which the temperature is then multiplied by afterwards (instances in fields_type_mod, homog_conv_bl, test_unstable). This completes a similar change already implemented in most of comorph in Issue CoMorph convection scheme refactoring #178 / PR CoMorph convection scheme refactoring #292, changing those instances which change answers at rounding-error levels and so were excluded from the earlier PR. Note the changes in init_test, to pass "temperature" into test_unstable; this is for the call to dry_adiabat inside test_unstable. Ensuring all dry_adiabat calls update T directly instead of scaling other things paves the way for a new option (to be lodged in another PR) to use dry static energy instead of potential temperature (in-which case the adiabatic adjustment is additive rather than multiplicative, so the scaling approach would be wrong).

10. Refactoring of averaging of parcel initial properties over sub-grid regions in init_mass_moist_frac... Each grid-cell is subdivided into liquid-cloud, mixed-phase-cloud, ice/rain and clear sub-grid regions, with initiating parcel properties calculated separately from each region and then averaged together. Mass-flux-weighted contributions from each region are added in the calls to add_region_parcel. Previously, the calls to normalise_init_parcel divided the summed parcel properties by the total mass-flux to retrieve mass-flux-weighted means. We have now removed the need for this final normalisation step by changing add_region_parcel so that it computes appropriate mass-flux-weights and accumulates a running mean of the parcel properties, instead of summing them scaled by mass-flux (this is algebraically identical but changes rounding-errors). Consistent with this, the normalisations done in normalise_init_parcel have been deleted. Note we still need that subroutine because it also applies safety-checks to the parcel initial moisture perturbation, but we have changed its name to finalise_init_parcel (since it no-longer normalises anything, its old name would be misleading).

11. Refactoring of sub-grid region diagnostics copied from calc_env_regions... There is some existing code here to set a fields super-array "fields_reg" with the properties of the current sub-grid region (populating the condensed water mixing-ratios using set_region_cond_fields) and then copy any requested diagnostics for the region properties into output arrays using fields_diags_copy. The existing code does this for all grid-points in the list, even where the current region has zero fraction. The field values at grid-cells with zero fraction were set to zero before calling fields_diags_copy, so we'd just copy zeros into the output. There were 2 problems with this:
    a) We did redundant data-copying and fields diags calculations on zeros.
    b) For the virtual temperature excess (i.e. buoyancy) diagnostic calculated in fields_diags_copy, the code would subtract the env Tv from a region Tv of zero, giving a spurious huge negative buoyancy at points with no fraction. It was annoying having to filter these out when plotting this diagnostic.
    These issues are avoided by compressing all inputs to fields_diags_copy onto points where current region fraction is nonzero, then scattering the diagnostics back to the full list afterwards. There is a comorph job in the UM rose-stem that outputs the region Tv excess diagnostics, so this change couldn't be included in PR CoMorph convection scheme refactoring #292 without changing answers.

Miscellaneous minor changes

1. The value of pi is specified to full precision in comorph_constants_mod (previously rounded to 10^-6).

2. The scaling of the subsidence virtual temperature increment delta_tv by the timestep and implicitness weight alpha_detrain is moved from calc_delta_tv (at control-level) to init_sublevs (within the plume-model). This simplifies the code / reduces the amount of "science" being confusingly done in the "control" routines.

3. In calc_delta_tv, dry-adiabatically-adjust the virtual temperature at the next level (virt_temp_2) directly in the call to dry_adiabat instead of passing out exner_ratio and scaling it in the calculation of delta_tv_cmpr.

4. A safety-check is added in homog_conv_bl at L568, to avoid rare div-by-zero where calculating the initiating parcel enthalpy perturbation (I think this could only happen when other changes not included here were added, but seems safest to include this check here).

5. In parcel_type_mod, in subroutine parcel_combine (which calculates mass-flux weighted mean properties where initiating parcels from different heights are combined together), the calculation of the mass-flux weight is modified to reduce rounding-error when "parcel_m" has a mass-flux much smaller than "parcel_a". We now do weight_m = mass_m / (mass_m + mass_a) instead of weight_m = 1 - weight_a = 1 - mass_a / (mass_m + mass_a). The former code is algebraically identical, but had larger rounding errors (it could spuriously fall to exactly zero when mass_m << mass_a).

6. Refactoring in calc_qvl_supersat... This subroutine calculates (among other things) the virtual temperature the parcel would have if all liquid-cloud were evaporated (used to estimate the buoyancy at the saturation boundary). The original calculation made a copy of the entire condensate super-array q_cond_l, only to add the liquid onto the vapour, set the liquid to zero, then pass into calc_virt_temp to compute the virtual temperature. The redundant copy of q_cond is now avoided by calculating the total-water mixing-ratio q_tot from the original q_vap and q_cond (q_tot is unaffected by conversion of water between species), and applying the 1/q_tot factor in the formula for virtual temperature in calc_qvl_supersat instead of doing it inside calc_vit_temp. This makes no difference to science but changes rounding-errors.

7. Compiler INLINE directives around the call to calc_env_region_tq_nb from calc_env_regions have been removed. These were only added in PR CoMorph convection scheme refactoring #292 to force the code to reproduce KGO (that PR moved the contained code inside its own subroutine, and that caused change-in-answers at bit-level due to a change in compiler optimisation behaviour). Since this PR is changing answers anyway, the compiler directives are now redundant and we should just let the compiler make its own choice.

8. In calc_init_par_fields, redundant TEMPORARY CODE lines that were added in PR CoMorph convection scheme refactoring #292 to force it to preserve answers exactly are removed. The affected code limits the parcel initial RH perturbation to ensure the parcel-core is not supersaturated. The temporary lines had been added to force the ratio of the actual core RH perturbation over the namelist value par_gen_rhpert to be exactly 1, which preserved KGO (the code was generalised in CoMorph convection scheme refactoring #292 to allow the RH perturbation to vary with convective type, though this functionality is not yet used). At present, removing the temporary lines makes the ratio go to very close to but not exactly 1, due to rounding errors. But the difference is negligible.

9. In calc_turb_parcel, the vertical interpolation of turbulent parcel initial perturbations from rho-levels to theta-levels has been tweaked slightly. This code sets the interpolation weight towards the neighbouring rho-level which has theta_l, q_t etc closer to the current theta-level value (i.e. interpolating in scalar-variable-space instead of height-space). The slightly rearranged code avoids duplicate calculation of the difference terms by storing them in "diffm", "diffp" and reusing. We also now don't let these terms go to zero, so that when there is exactly zero vertical gradient the interpolation weight goes to 0.5 instead of 0.0 (this hardly ever happens so negligible impact).

10. In cfl_limit_init_mass, we delete the code for special treatment when there is only 1 convection type (the multi-type implementation we now default to is algebraically identical to it when n_conv_types == 1, but changes answers at bit-level, hence we retained the single-type implementation under an IF test in CoMorph convection scheme refactoring #292 to preserve answers).

11. In calc_cape, the calculation of the height-interval across the current sub-level step has been simplified / made more succinct (in-line max,min limits on the k-index at the top and bottom replace if-tests on i_lev > i_prev and i_lev < i_next). The new calculation should be a slightly cheaper way to get the same answer, but changes answers at bit-level. This change also avoids a spurious compiler warning from gfortran rigorous compiles; the compiler complained that "i_lev" may reference "sublevs" out-of-bounds, which wasn't actually the case. But the refactored code both avoids the warning and is tidier :)

12. calc_mean_q_cl_with_core; rearrange calculation inside if-test; we have just multiplied both sides of the comparison by "core_mean_ratio" (should be slightly cheaper as now uses multiply instead of divide). Changes answers at bit-level.

13. 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 code detected grid-levels spanning the saturation height based on the in-plume liquid-water content "q_cl" being zero at the start and non-zero at the end. It then interpolated to find the saturation height using the variables "prev_ss", "next_ss" (parcel supersaturation w.r.t. liquid water at start and end of the level-step). However, occasionally due to rounding errors, it was possible to meet the first condition for a cloud-base-spanning grid-level based on "q_cl", but not have "ss" change sign between prev and next. This could make the interpolation to find "sat_height" yield nonsense. The existing code trapped these instances using additional checks (Beware Earth radius in single precision #326 recently added the trap to reset "interp" to zero when this happens, to avoid a very rare div-by-zero which didn't propagate as was overwritten by other checks, but could cause rigorous compile runs to fail). We now remove the need for this complexity by just detecting the cloud-base-spanning grid-level based on change of sign of "ss" rather than "q_cl" (then the interpolation is guaranteed to yield sensible answers, and the subsequent checks on prev_ss(ic) * next_ss(ic) >= zero are redundant and are deleted).

14. calc_sat_height; redundant TEMPORARY CODE lines removed. In the case of sublevs(:,j_delta_tv,i_sat) (env subsidence virtual temperature increment at the accurate saturation height of the updraft), this has already been interpolated between the values at the start and end of the level-step in the preceding loop. At present the value at start of level-step is set equal to the value at the end in init_sublevs (other changes not included in this PR will allow delta_tv to vary linearly over the level-step). So the interpolation should just yield the same value at "prev" and "next", but introduced inconsequential rounding-error differences. Hence the temporary line to force sublevs(:,j_delta_tv,i_sat) to exactly equal sublevs(:,j_delta_tv,i_prev) was added in CoMorph convection scheme refactoring #292 to preserve KGO, but we now safely remove it. The 2nd TEMPORARY CODE block removed (L331) redundantly repeated the vertical interpolation of the env virtual temperature to the saturation height (whereas now we just use the value already stored in sublevs at L303). This temporary code was only needed in CoMorph convection scheme refactoring #292 to avoid a change in compiler optimisation behaviour which changed answers for CCE high-optimisation jobs.

15. Safety-check added to avoid div-by-zero in ice_rain_to_graupel, in calculation which updates the graupel heat / vapour exchange coefficients (I think this is a bug that only got exposed when using other changes from comorph B, so fixing it is unlikely to impact any existing configurations; best to fix it just in case!)

16. Brackets added in calculations in check_negatives (part of the comorph microphysics implicit phase-change solver which modifies the solution to avoid creating negative condensate mass). The added brackets avert an occasional problem caused by rounding errors, in the case where condensation / deposition has been exactly cancelled by freezing / melting.

17. Safety-check added in set_dqsatdt to avoid div-by-zero. We're only likely to hit this check in the case where a massless grid-point has been included in the compression list (such that parcel fields including temperature might be zero). In this case, the result of the div-by-zero would not have propagated anyway, but the div-by-zero might trigger run-time errors in rigorous compile jobs so its best to avoid it.

18. Refactored calculation of non-liquid-cloud region temperature in calc_env_region_tq_nb. The existing code adds on "qc_tot" in the calculations of "qc_tot_noliq", then redundantly subtracts it off again where that is used a couple lines down. We simplify and remove the redundant add / subtract by instead calculating the excess of "qc_tot" above the grid-mean value in the first place, replacing "qc_tot_noliq" with "qc_excess_noliq". This makes no difference algebraically, but changes rounding errors so changes answers at bit-level.

19. in momentum_eqn (calculates drag on the horizontal and vertical wind components of the parcel), we move where the scaling by timestep "delta_t" is done, so-as to avoid a redundant divide (algebraically identical but changes answers).

20. In set_det (detrainment calculation), the condition for using a taylor expansion when the full formula loses precision is tightened-up (full formula still needs frac < 1 by a numerical tolerance even if x_edge is). Note this will have negligible impact on CoMA9 because it has min_cmr=2.0, which restricts the power in the assumed power-law PDF to be > zero. I think we only get frac > x_edge when the PDF power is negative.

21. Refactored heterogeneous nucleation (subroutine ice_nucleation) to explicitly set q_ice to a negligibly-small non-zero seed value, instead of using a redundant rate coefficient ("coef_hetnuc") that we hardwired to be negligibly small anyway. Note the subroutine has been reordered so we now do homogeneous freezing first, which has thoroughly confused the diff viewer. Note I've also moved the copying of "dq_frz" into output diagnostic arrays to before the code which recomputes the compression indices (another change that confounds the diff viewer); this is because for simplicity the rearranged code now defines "dq_frz" on the full list of points, combining the increments from heterogeneous and homogeneous nucleation, and we need to use the compression indices "index_ic_liq" (defined before we remove the fully frozen points) in the loop to copy to the diagnostics.

    Consistent with this change, the timestep "delta_t" is removed from the argument list where ice_nucleation is called from microphysics_1, the comments explaining the issues with performing heterogeneous freezing at this point in the timestep are updated, and the now unused "coef_hetnuc" is deleted from comorph_constants_mod.

22. In set_test_profiles (which sets made-up profiles used in the comorph stand-alone test), I've modified the profiles to make them look more realistic. They now have a pseudo well-mixed vapour mixing-ratio profile near the surface instead of having RH constant with height. I also now apply a 1-2-1 smoother in height to Tv and qv, to make the inversion slightly less sharp. These changes don't affect any model-runs, only the stand-alone test.

23. In the lfric_apps comorph namelist meta-data (interfaces/physics_schemes_interface/rose-meta/um-convection/HEAD/rose-meta.conf), the precision of the namelist input "par_radius_knob" is changed from "double" to "default" (i.e. from 64-bit to 32-bit). It had only been set to "double" in CoMorph further KGO-preserving refactoring #717 so-as to preserve KGO with the previous trunk code which hardwired the comorph_um_namelist_mod variable (which is 64-bit) in um_physics_init. Reading the value from the namelist to a 32-bit variable and then converting to 64-bit introduced rounding-error which slightly changed the value used in comorph. However this is inconsequential; since the present branch changes answers anyway, we may as well set the precision of "par_radius_knob" the same as all the other comorph namelist inputs (don't want to mislead people into thinking this variable had to be extremely accurate for the scheme to work!)