Skip to content

Hail diagnostic vn14.2 - #21

Merged
James Bruten (james-bruten-mo) merged 13 commits into
MetOffice:mainfrom
paulfield2024:hail_diagnostic_vn14.2
Oct 8, 2026
Merged

James Bruten (james-bruten-mo) merged 13 commits into
MetOffice:mainfrom
paulfield2024:hail_diagnostic_vn14.2

Conversation

@paulfield2024

@paulfield2024 paulfield2024 commented Aug 14, 2026 •

Copy link
Copy Markdown
Contributor

PR Summary

Sci/Tech Reviewer: mo-sabel
Code Reviewer: Pierre Siddall (@Pierre-siddall)

Hail diagnostic based on Mason 1956 "Melting of Hailstones"

Code Quality Checklist

(Some checks are automatically carried out via the CI pipeline)

  • I have performed a self-review of my own code
  • My code follows the project's style guidelines
  • Comments have been included that aid understanding and enhance the readability of the code
  • My changes generate no new warnings

Testing

  • If shared files have been modified, I have run the UM and LFRic Apps rose
    stem suites
  • If any tests fail (rose-stem or CI) the reason is understood and
    acceptable (eg. kgo changes)
  • I have added tests to cover new functionality as appropriate (eg. system
    tests, unit tests, etc.)

trac.log

Security Considerations

  • I have reviewed my changes for potential security issues
  • Sensitive data is properly handled (if applicable)
  • Authentication and authorisation are properly implemented (if applicable)

Performance Impact

  • Performance of the code has been considered and, if applicable, suitable performance measurements have been conducted

AI Assistance and Attribution

  • Some of the content of this change has been produced with the assistance of Generative AI tool name (e.g., Met Office Github Copilot Enterprise, Github Copilot Personal, ChatGPT GPT-4, etc) and I have followed the Simulation Systems AI policy (including attribution labels)

Documentation

  • Where appropriate I have updated documentation related to this change and confirmed that it builds correctly

Sci/Tech Review

  • I understand this area of code and the changes being added
  • The proposed changes correspond to the pull request description
  • Documentation is sufficient (do documentation papers need updating)
  • Sufficient testing has been completed

Please alert the code reviewer via a tag when you have approved the SR

Code Review

  • All dependencies have been resolved
  • Related Issues have been properly linked and addressed
  • CLA compliance has been confirmed
  • Code quality standards have been met
  • Tests are adequate and have passed
  • Documentation is complete and accurate
  • Security considerations have been addressed
  • Performance impact is acceptable

@github-actions github-actions Bot added the cla-signed The CLA has been signed as part of this PR - added by GA label Aug 14, 2026
@paulfield2024

paulfield2024 commented Aug 14, 2026 •

Copy link
Copy Markdown
Contributor Author

This is linked to UM PR to plumb through the hail diagnostic requests
UM#118

Purpose: Implements a new physically-based hail-size/melting diagnostic in
CASIM microphysics.

Files changed (3):

src/hail_diag_fast_mod.F90 (new, ~206 lines)
New module hail_diagnostic_fast_mod with diagnose_hail_fast: builds a
hail size distribution from the graupel PSD at the melting level, then
integrates hailstone melting downward to the surface using Mason (1971)
melting physics. Returns the largest surviving hailstone diameter at the
surface, a critical/threshold diameter, a hail flag, and the surface hail
precipitation rate.
src/generic_diagnostic_variables.F90
Adds an l_hail flag and five new 2D diagnostic arrays
(hail_d_max_sfc, hail_d_crit, hail_d0_thresh, hail_flag_sfc,
hail_precip_rate) with allocate/deallocate/reset logic, following the
existing pattern used for other diagnostics (e.g. radar).
src/micro_main.F90
Calls diagnose_hail_fast per grid column when l_hail is set, and stores
the results into the new diagnostic arrays.
Impact on results/performance: Additive and opt-in (gated by l_hail,
only set when the UM requests the corresponding STASH items). Adds one
physics-based diagnostic call per column; negligible cost when hail
diagnostics are disabled. No existing microphysics behavior is changed.

Detail: src/hail_diag_fast_mod.F90 (module hail_diagnostic_fast_mod)

Algorithm, per grid column:

Find the melting level (find_melting_level) — scans from the top
of the column down for the 0°C crossing (T0 = 273.15 K) and
linearly interpolates its height. If the whole column (including the
surface) is already below freezing, there's no melting layer aloft and
it returns ierr=2.

Build a hail size distribution at the melting level — reads the
graupel particle size distribution parameters (dist_mu, dist_lambda,
dist_n0) from CASIM's distributions module at the melting level, then
discretizes hail into 7 fixed diameter bins (edges hardcoded from
1 mm to 300 mm) and computes the number concentration in each bin from
the gamma-distribution graupel PSD. If total concentration across all
bins is below a threshold, returns ierr=2 (no significant hail).
Melt each bin down to the surface — for each size bin, steps
level-by-level from the melting level down to the surface, shrinking the
stone's radius each level via mason_melt (a Mason 1971 melting-rate
formulation using latent heat of fusion, thermal conductivities of water/
air, and a ventilation coefficient derived from the Reynolds number of
the falling stone). Terminal velocity uses a simple power law
(V = a_v * D_cm^b_v).
Accumulate surface results — if a stone in a given bin survives
(radius > 0) with concentration above threshold, it contributes to the
surface hail precipitation rate (concentration × velocity × mass × air density), updates the flag (hail_flag = 1), and updates
D_max_sfc to the largest surviving diameter. precip_rate is the sum
of all bins' contributions (kg m⁻² s⁻¹).
Notes:

beta is set to 0 to explicitly ignore condensation evaporation effects. This allows hailstones to fall further.
density of hail is set to 900kg/m3 whereas the density of the graupel for the distribution from which the hail are spawned is lower. Again this allows the hailstones to fall further before melting completely.

The 7 hail-size bins and their edges (1 mm–300 mm) are hardcoded rather
than derived from the actual distribution — a coarse, fixed
discretization.

How it fits together
The two branches are complementary: casim/hail_diagnostic adds the physics
and storage for the new hail diagnostics, and um/hail_diagnostic wires up
STASH so a UM run can request them as output. Together they let a run output:
a hail flag, max hail diameter at the surface (mm), hail precipitation rate,
and accumulated hail precipitation amount.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this looks good. The method is based on literature sources and assumes a "maximum" scenario or worst case scenario for graupel getting to the ground. This is a sensible choice for the initial implementation.

real(wp) :: hailbin(nh)
real(wp) :: hailconc(nh)

data haildedge /1e-3, 3e-3,6e-3,1e-2,3e-2,6e-2,1e-1, 3e-1/ ! initial hail size edges at melting layer diameter in m

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

bins centres not identical to description in um PR, 20 mm (here) not 15mm (in PR explanation)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i couldnt see this. Its mentioned as bin edges in the PR description above.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment thread src/hail_diag_fast_mod.F90 Outdated
integer :: k_ml, nn, k
real(wp) :: a, a0, v, da

integer, parameter :: nh = 7

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can this be derived from hailedge, in case bins change in future?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i will just add a comment to mention that this is linked to haildedge

Comment thread src/hail_diag_fast_mod.F90 Outdated
a=a0
do k = k_ml,1,-1
v=a_v*(2.0*a*100.0)**b_v * sqrt(rho0) !diam in cm
da=mason_melt(a,a0,dz_in(k),t(k)-273.15,v,rhoi)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use T0 here, as it's defined anyway?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

@paul-barrett

Copy link
Copy Markdown

References Mason 1971 - did AI pick this and alter the reference? Was initially using the paper from 1956 rather than the book of 1971

@mo-sabel mo-sabel left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Overall, the approach looks physically reasonable and provides a pragmatic way of diagnosing hail from the existing CASIM graupel distribution. One underlying assumption underpinning the diagnostic is that the graupel PSD at the melting level is representative of potential hail. It may be worth considering how the diagnostic behaves in cloud regimes that can generate graupel without necessarily being regarded as hail-producing environments (e.g. mixed-phase stratocumulus or shallow convection). The existing concentration thresholds and subsequent melting calculations may already filter out many such cases, but testing the diagnostic in a benign mixed-phase cloud case could help assess the sensitivity of the results to this assumption.

Comment thread src/hail_diag_fast_mod.F90 Outdated
real(wp) :: da, C, Re, beta
beta=0.0 ! ignore condensation/evap

Re=v*(2*a)/visc

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Double check the units.
visc = 1.7e-5 kg/m/s (dynamic viscosity of air?)
would give units of Re m3/kg
Do you need to multiply by air density?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep - done. As its close to 1 and sqrt'd there should be little impact.

Comment thread src/hail_diag_fast_mod.F90 Outdated
a0=haild(nn)/2.0 ! convert to radius at melting level
a=a0
do k = k_ml,1,-1
v=a_v*(2.0*a*100.0)**b_v * sqrt(rho0) !diam in cm

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this have air density as per the comment?
! ---- terminal velocity power law: V = a_v * (D_cm)^b_v * sqrt(rho0_ref/rho_air) ----

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep done. Again it is close to 1 and sqrt'd so there will be little impact.

hailconc(nn)=N_t*lam**(mu_g+1)/gamma(mu_g+1)*haild(nn)**mu_g*exp(-lam*haild(nn))*hailbin(nn) !conc from graupel dist
end do

if ( maxval(hailconc) < Nthresh) then

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Description in PR states "If total concentration across all
bins is below a threshold, returns ierr=2 (no significant hail)". Perhaps reword to "If the concentration in every bin is below a threshold, returns ierr=2 (no significant hail)". With the current wording I thought it was the total conc in the PSD

@paulfield2024
paulfield2024 marked this pull request as ready for review September 17, 2026 07:44

@mo-sabel mo-sabel left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding air density to fallspeed and Re looks good. Now also includes air density argument in mason_melt()

@mo-sabel mo-sabel left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just needs air density added in arguments in mason_melt()
Already captured in next commit c942480

@mo-sabel mo-sabel left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Approved but noting that rhoa needs added as an argument in mason_melt. That has been implemented in the next commit

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi paulfield2024, thanks for doing this just a trivial suggestion regarding some seemingly redundant variable initialisation. I think removing these would be helpful to ensure failures don't pass silently .

Comment thread src/hail_diag_fast_mod.F90 Outdated
@james-bruten-mo
James Bruten (james-bruten-mo) merged commit 996db30 into MetOffice:main Oct 8, 2026
5 of 7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

cla-signed The CLA has been signed as part of this PR - added by GA

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants