From 3305acf7ce60870e7dc093858add3b5d8c49a0fc Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Tue, 4 Aug 2026 10:29:09 +0100 Subject: [PATCH 1/7] adding in the hail_diagnostics from fcm vn13.9_hail_diag --- src/generic_diagnostic_variables.F90 | 50 ++++++++++++++++++++++++++++ src/micro_main.F90 | 23 +++++++++++++ 2 files changed, 73 insertions(+) diff --git a/src/generic_diagnostic_variables.F90 b/src/generic_diagnostic_variables.F90 index cd2b16a..b259ecc 100644 --- a/src/generic_diagnostic_variables.F90 +++ b/src/generic_diagnostic_variables.F90 @@ -323,6 +323,12 @@ MODULE generic_diagnostic_variables LOGICAL :: l_dqs = .FALSE. LOGICAL :: l_dqg = .FALSE. + !--------------------------------- + ! hail diagnostics + !--------------------------------- + LOGICAL :: l_hail = .FALSE. + + !-------------------------------- ! 2D variable arrays @@ -354,6 +360,12 @@ MODULE generic_diagnostic_variables REAL, ALLOCATABLE :: dbz_l(:,:,:) REAL, ALLOCATABLE :: dbz_r(:,:,:) + REAL, ALLOCATABLE :: hail_d_max_sfc(:,:) + REAL, ALLOCATABLE :: hail_d_crit(:,:) + REAL, ALLOCATABLE :: hail_d0_thresh(:,:) + REAL, ALLOCATABLE :: hail_flag_sfc(:,:) + REAL, ALLOCATABLE :: hail_precip_rate(:,:) + ! Process rate diagnostics REAL, ALLOCATABLE :: phomc(:,:,:) REAL, ALLOCATABLE :: pinuc(:,:,:) @@ -631,6 +643,23 @@ SUBROUTINE allocate_diagnostic_space(i_start, i_end, j_start, j_end, k_start, k_ END IF ! casdiags % l_radar +IF ( casdiags % l_hail ) THEN + + ALLOCATE ( casdiags % hail_d_max_sfc(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_d_crit(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_flag_sfc(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_d0_thresh(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_precip_rate(i_start:i_end, j_start:j_end) ) + + casdiags % hail_d_max_sfc(:,:) = zero_real_wp + casdiags % hail_d_crit(:,:) = zero_real_wp + casdiags % hail_flag_sfc(:,:) = zero_real_wp + casdiags % hail_d0_thresh(:,:) = zero_real_wp + casdiags % hail_precip_rate(:,:) = zero_real_wp + + +END IF ! casdiags % l_hail + IF (casdiags % l_phomc) THEN ALLOCATE ( casdiags % phomc(i_start:i_end, j_start:j_end, k_start:k_end) ) !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(k) & @@ -2589,6 +2618,26 @@ SUBROUTINE deallocate_diagnostic_space() DEALLOCATE ( casdiags % phomc ) END IF +IF (casdiags % l_hail) THEN + + IF ( ALLOCATED ( casdiags % hail_precip_rate ) ) THEN + DEALLOCATE( casdiags % hail_precip_rate ) + END IF + IF ( ALLOCATED ( casdiags % hail_d0_thresh ) ) THEN + DEALLOCATE( casdiags % hail_d0_thresh ) + END IF + IF ( ALLOCATED ( casdiags % hail_flag_sfc ) ) THEN + DEALLOCATE( casdiags % hail_flag_sfc ) + END IF + IF ( ALLOCATED ( casdiags % hail_d_crit ) ) THEN + DEALLOCATE( casdiags % hail_d_crit ) + END IF + IF ( ALLOCATED ( casdiags % hail_d_max_sfc ) ) THEN + DEALLOCATE( casdiags % hail_d_max_sfc ) + END IF + +ENDIF + IF (casdiags % l_radar) THEN IF ( ALLOCATED ( casdiags % dbz_r ) ) THEN DEALLOCATE( casdiags % dbz_r ) @@ -2640,6 +2689,7 @@ SUBROUTINE deallocate_diagnostic_space() casdiags % l_process_rates = .FALSE. casdiags % l_tendency_dg = .FALSE. casdiags % l_radar = .FALSE. +casdiags % l_hail = .FALSE. IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) diff --git a/src/micro_main.F90 b/src/micro_main.F90 index c9e2a37..3dc12a2 100644 --- a/src/micro_main.F90 +++ b/src/micro_main.F90 @@ -419,6 +419,8 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & USE yomhook, ONLY: lhook, dr_hook USE parkind1, ONLY: jprb, jpim + + use hail_diagnostic_fast_mod, ONLY: diagnose_hail_fast implicit none @@ -504,6 +506,10 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & real(wp) :: precip_g1d(nz, nxy_inner) real(wp) :: waterpath + + real(wp) :: D_crit, D0_thresh, D_max_sfc, precip_rate_hail + real(wp) :: hail_flag ! .true. if any hail reaches the surface + integer :: ierr ! 0=ok, 1=no graupel, 2=no melting level real(wp) :: dbz_tot_c(nz), dbz_g_c(nz), dbz_i_c(nz), & dbz_s_c(nz), dbz_l_c(nz), dbz_r_c(nz) @@ -783,6 +789,23 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & casdiags % dbz_r(i,j, k_start:k_end) = dbz_r_c(:) end if ! casdiags % l_radar + + + if ( casdiags % l_hail ) then + + call diagnose_hail_fast( ixy_inner, nz, nq, qfields(:,:,ixy_inner), cffields(:,:,ixy_inner), & + D_crit, D0_thresh, D_max_sfc, precip_rate_hail, hail_flag, ierr ) + + + casdiags % hail_d_max_sfc(i,j)=D_max_sfc + casdiags % hail_d_crit(i,j)=D_crit + casdiags % hail_flag_sfc(i,j)=hail_flag + casdiags % hail_d0_thresh(i,j)=D0_thresh + casdiags % hail_precip_rate(i,j)=precip_rate_hail + + + + end if ! casdiags % l_hail if ( casdiags % l_tendency_dg ) then DO k = k_start, k_end From 4eb9eedb189490ed44a7b46dd203e4a1f885614b Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Tue, 4 Aug 2026 13:43:03 +0100 Subject: [PATCH 2/7] adding in new subroutine --- src/hail_diag_fast_mod.F90 | 206 +++++++++++++++++++++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 src/hail_diag_fast_mod.F90 diff --git a/src/hail_diag_fast_mod.F90 b/src/hail_diag_fast_mod.F90 new file mode 100644 index 0000000..41b9da3 --- /dev/null +++ b/src/hail_diag_fast_mod.F90 @@ -0,0 +1,206 @@ +!=============================================================================== +! MODULE: hail_diagnostic_fast_mod + +!=============================================================================== +module hail_diagnostic_fast_mod + + USE precision, ONLY: wp + + USE umPrintMgr, ONLY: umPrint, umMessage + + implicit none + private + public :: diagnose_hail_fast + + ! ---- physical constants ---- + real(wp), parameter :: pi = 3.14159265358979_wp + real(wp), parameter :: Rd = 287.05_wp ! J/kg/K, dry air + real(wp), parameter :: Rv = 461.5_wp ! J/kg/K, water vapor + real(wp), parameter :: Lf = 3.34e5_wp ! J/kg, latent heat of fusion + real(wp), parameter :: Lv = 2.5e6_wp ! J/kg, latent heat of vaporization + real(wp), parameter :: Kw0 = 0.56_wp ! W/m/K, thermal conductivity of water (~0C) + real(wp), parameter :: Ka0 = 2.4e-2_wp ! W/m/K, thermal conductivity of air (~0C) + real(wp), parameter :: Dv0 = 2.21e-5_wp ! m2/s, vapor diffusivity in air (~0C, 1000hPa) + real(wp), parameter :: T0 = 273.15_wp ! K, melting point + real(wp), parameter :: es0 = 611.2_wp ! Pa, saturation vapor pressure at 0C + real(wp), parameter :: Pr_no = 0.71_wp ! Prandtl number of air + real(wp), parameter :: Sc_no = 0.60_wp ! Schmidt number of air + real(wp), parameter :: rho0_ref = 1.2_wp ! kg/m3, reference air density for V(D) + real(wp), parameter :: visc = 1.7e-5_wp ! visc kg/m/s + + ! ---- terminal velocity power law: V = a_v * (D_cm)^b_v * sqrt(rho0_ref/rho_air) ---- + real(wp), parameter :: a_v = 14.0_wp ! for cm; MASON 1971 4.41(Dmm)0.5 + real(wp), parameter :: b_v = 0.5_wp + + ! ---- minimum thresholds ---- + real(wp), parameter :: qg_min = 1.0e-6_wp ! kg/kg, minimum graupel to bother + +contains + + subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, & + D_crit, D0_thresh, D_max_sfc, precip_rate, & + hail_flag, ierr ) + + USE distributions, ONLY: query_distributions, dist_lambda, dist_mu, dist_n0 + USE mphys_parameters, ONLY: graupel_params + USE passive_fields, ONLY: Tdegk, dz, pressure, rho + use mphys_switches, only: i_qv, i_qg + + real(wp), parameter :: Nthresh=1e-6 + + integer, intent(in) :: nz, nq + integer, intent(in) :: ixy_inner + + real(wp), intent(in) :: qfields(nz,nq) + real(wp), intent(in) :: cffields(nz,nq) + + real(wp), intent(out) :: D_crit ! melting-level critical diameter [m] + real(wp), intent(out) :: D0_thresh ! melting-level size at the concentration threshold [m] + real(wp), intent(out) :: D_max_sfc ! largest hailstone diameter at surface [m] + real(wp), intent(out) :: precip_rate ! surface hail mass flux [kg/m2/s] + ! (x3600 for mm/hr liquid-equivalent) + real(wp), intent(out) :: hail_flag ! .true. if any hail reaches the surface + integer, intent(out) :: ierr ! 0=ok, 1=no graupel, 2=no melting level + + real(wp) :: rhoa(nz), qv(nz), t(nz), p(nz), qg(nz), dz_in(nz) + real(wp) :: qg_ml, ng_ml, t_ml, p_ml, qv_ml, z_ml, rho_air_ml, rho_air_sfc + real(wp) :: N_t, M_t, lam, N0, Mtotal, n_exp + real(wp) :: rho_g, mu_g + integer :: k_ml, nn, k + real(wp) :: a, a0, v, da + + integer, parameter :: nh = 7 + real(wp), parameter :: rhoi = 900.0 + real(wp) :: haild(nh) + real(wp) :: haildf(nh) + real(wp) :: haildedge(nh+1) + real(wp) :: hailprecipf(nh) + 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 + + rho_g=graupel_params%density + qv=qfields(:, i_qv) + qg=qfields(:, i_qg) + t(:)=TdegK(:,ixy_inner) + p(:)=pressure(:,ixy_inner) + rhoa(:)=rho(:,ixy_inner) + dz_in(:)=dz(:,ixy_inner) + + D_crit = 0.0_wp + D0_thresh = 0.0_wp + D_max_sfc = 0.0_wp + precip_rate = 0.0_wp + hail_flag = 0.0_wp + ierr = 0 + + if (maxval(qg) < qg_min) then + ierr = 1 + return + end if + + + call find_melting_level(nz, dz_in, t, k_ml, z_ml, ierr) + if (ierr /= 0) return + mu_g=dist_mu(k_ml,graupel_params%id) + lam=dist_lambda(k_ml, graupel_params%id) + N_t=dist_n0(k_ml, graupel_params%id) + + do nn=1,nh + haild(nn)=(haildedge(nn+1)+haildedge(nn))/2.0 !initial diameters + haildf(nn)=0.0_wp !final diameter + hailprecipf(nn)=0.0_wp ! precip in bins + hailbin(nn)=(haildedge(nn+1)-haildedge(nn)) ! bin width + 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 + ierr = 2 + return + end if + + do nn=1,nh + 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 !diam in cm + da=mason_melt(a,a0,dz_in(k),t(k)-273.15,v,rhoi) + a=a-da + end do + if ( a .gt. 0.0 .and. hailconc(nn) .gt. Nthresh) then + haildf(nn)=2.0*a + hailprecipf(nn)=hailconc(nn)*v*pi/6.0*haildf(nn)**3*rhoi*rhoa(1) + D_max_sfc=2.0*a + hail_flag=1.0 + else + haildf(nn)=0.0 + hailprecipf(nn)=0.0 + end if + end do + precip_rate=sum(hailprecipf(:)) !kg m-2 s-1 + + + + + end subroutine diagnose_hail_fast + + subroutine find_melting_level(nz, dz, t, k_ml, z_ml, ierr) + + integer, intent(in) :: nz + real(wp), intent(in) :: dz(nz), t(nz) + integer, intent(out) :: k_ml + real(wp), intent(out) :: z_ml + integer, intent(out) :: ierr + integer :: k + real(wp) :: frac + + ierr = 3 + k_ml = -1 + z_ml = -999.0_wp + + if (t(1) >= T0) then + do k = nz-1, 1, -1 + + if (t(k+1) < T0 .and. t(k) >= T0) then + frac = (T0 - t(k)) / (t(k+1) - t(k)) + z_ml = SUM(dz(1:k-1)) + frac*dz(k) + k_ml = k + ierr = 0 + + return + end if + end do + end if + + if (t(1) < T0) then + ! surface already below freezing -- no melting layer above the surface + z_ml = 0.0_wp + k_ml = 1 + ierr = 2 + end if + + end subroutine find_melting_level + + !================== + ! Mason melting + !================== + function mason_melt(a,a0,dz,t,v,rhoi) result(da) + real(wp), intent(in) :: a,a0,dz,t,rhoi,v + real(wp) :: da, C, Re, beta + beta=0.0 ! ignore condensation/evap + + Re=v*(2*a)/visc + C=1.6+0.3*Re**0.5 + da=(Kw0*t*dz/Lf/rhoi/v)/ & + ((a0-a)*a/a0+(Kw0/(C*(Ka0+Lv*Dv0*beta)))*a**2/a0) + + + end function mason_melt + + + +end module hail_diagnostic_fast_mod + + + From db68903d2df172966abcedba42cf6a94b35dcc0e Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Tue, 4 Aug 2026 10:29:09 +0100 Subject: [PATCH 3/7] adding in the hail_diagnostics from fcm vn13.9_hail_diag --- src/generic_diagnostic_variables.F90 | 50 ++++++++++++++++++++++++++++ src/micro_main.F90 | 23 +++++++++++++ 2 files changed, 73 insertions(+) diff --git a/src/generic_diagnostic_variables.F90 b/src/generic_diagnostic_variables.F90 index cd2b16a..b259ecc 100644 --- a/src/generic_diagnostic_variables.F90 +++ b/src/generic_diagnostic_variables.F90 @@ -323,6 +323,12 @@ MODULE generic_diagnostic_variables LOGICAL :: l_dqs = .FALSE. LOGICAL :: l_dqg = .FALSE. + !--------------------------------- + ! hail diagnostics + !--------------------------------- + LOGICAL :: l_hail = .FALSE. + + !-------------------------------- ! 2D variable arrays @@ -354,6 +360,12 @@ MODULE generic_diagnostic_variables REAL, ALLOCATABLE :: dbz_l(:,:,:) REAL, ALLOCATABLE :: dbz_r(:,:,:) + REAL, ALLOCATABLE :: hail_d_max_sfc(:,:) + REAL, ALLOCATABLE :: hail_d_crit(:,:) + REAL, ALLOCATABLE :: hail_d0_thresh(:,:) + REAL, ALLOCATABLE :: hail_flag_sfc(:,:) + REAL, ALLOCATABLE :: hail_precip_rate(:,:) + ! Process rate diagnostics REAL, ALLOCATABLE :: phomc(:,:,:) REAL, ALLOCATABLE :: pinuc(:,:,:) @@ -631,6 +643,23 @@ SUBROUTINE allocate_diagnostic_space(i_start, i_end, j_start, j_end, k_start, k_ END IF ! casdiags % l_radar +IF ( casdiags % l_hail ) THEN + + ALLOCATE ( casdiags % hail_d_max_sfc(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_d_crit(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_flag_sfc(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_d0_thresh(i_start:i_end, j_start:j_end) ) + ALLOCATE ( casdiags % hail_precip_rate(i_start:i_end, j_start:j_end) ) + + casdiags % hail_d_max_sfc(:,:) = zero_real_wp + casdiags % hail_d_crit(:,:) = zero_real_wp + casdiags % hail_flag_sfc(:,:) = zero_real_wp + casdiags % hail_d0_thresh(:,:) = zero_real_wp + casdiags % hail_precip_rate(:,:) = zero_real_wp + + +END IF ! casdiags % l_hail + IF (casdiags % l_phomc) THEN ALLOCATE ( casdiags % phomc(i_start:i_end, j_start:j_end, k_start:k_end) ) !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(k) & @@ -2589,6 +2618,26 @@ SUBROUTINE deallocate_diagnostic_space() DEALLOCATE ( casdiags % phomc ) END IF +IF (casdiags % l_hail) THEN + + IF ( ALLOCATED ( casdiags % hail_precip_rate ) ) THEN + DEALLOCATE( casdiags % hail_precip_rate ) + END IF + IF ( ALLOCATED ( casdiags % hail_d0_thresh ) ) THEN + DEALLOCATE( casdiags % hail_d0_thresh ) + END IF + IF ( ALLOCATED ( casdiags % hail_flag_sfc ) ) THEN + DEALLOCATE( casdiags % hail_flag_sfc ) + END IF + IF ( ALLOCATED ( casdiags % hail_d_crit ) ) THEN + DEALLOCATE( casdiags % hail_d_crit ) + END IF + IF ( ALLOCATED ( casdiags % hail_d_max_sfc ) ) THEN + DEALLOCATE( casdiags % hail_d_max_sfc ) + END IF + +ENDIF + IF (casdiags % l_radar) THEN IF ( ALLOCATED ( casdiags % dbz_r ) ) THEN DEALLOCATE( casdiags % dbz_r ) @@ -2640,6 +2689,7 @@ SUBROUTINE deallocate_diagnostic_space() casdiags % l_process_rates = .FALSE. casdiags % l_tendency_dg = .FALSE. casdiags % l_radar = .FALSE. +casdiags % l_hail = .FALSE. IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) diff --git a/src/micro_main.F90 b/src/micro_main.F90 index c9e2a37..3dc12a2 100644 --- a/src/micro_main.F90 +++ b/src/micro_main.F90 @@ -419,6 +419,8 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & USE yomhook, ONLY: lhook, dr_hook USE parkind1, ONLY: jprb, jpim + + use hail_diagnostic_fast_mod, ONLY: diagnose_hail_fast implicit none @@ -504,6 +506,10 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & real(wp) :: precip_g1d(nz, nxy_inner) real(wp) :: waterpath + + real(wp) :: D_crit, D0_thresh, D_max_sfc, precip_rate_hail + real(wp) :: hail_flag ! .true. if any hail reaches the surface + integer :: ierr ! 0=ok, 1=no graupel, 2=no melting level real(wp) :: dbz_tot_c(nz), dbz_g_c(nz), dbz_i_c(nz), & dbz_s_c(nz), dbz_l_c(nz), dbz_r_c(nz) @@ -783,6 +789,23 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & casdiags % dbz_r(i,j, k_start:k_end) = dbz_r_c(:) end if ! casdiags % l_radar + + + if ( casdiags % l_hail ) then + + call diagnose_hail_fast( ixy_inner, nz, nq, qfields(:,:,ixy_inner), cffields(:,:,ixy_inner), & + D_crit, D0_thresh, D_max_sfc, precip_rate_hail, hail_flag, ierr ) + + + casdiags % hail_d_max_sfc(i,j)=D_max_sfc + casdiags % hail_d_crit(i,j)=D_crit + casdiags % hail_flag_sfc(i,j)=hail_flag + casdiags % hail_d0_thresh(i,j)=D0_thresh + casdiags % hail_precip_rate(i,j)=precip_rate_hail + + + + end if ! casdiags % l_hail if ( casdiags % l_tendency_dg ) then DO k = k_start, k_end From d8c9d6135ac3f882134ef3bcc21c0168d599dae7 Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Tue, 4 Aug 2026 13:43:03 +0100 Subject: [PATCH 4/7] adding in new subroutine --- src/hail_diag_fast_mod.F90 | 206 +++++++++++++++++++++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 src/hail_diag_fast_mod.F90 diff --git a/src/hail_diag_fast_mod.F90 b/src/hail_diag_fast_mod.F90 new file mode 100644 index 0000000..41b9da3 --- /dev/null +++ b/src/hail_diag_fast_mod.F90 @@ -0,0 +1,206 @@ +!=============================================================================== +! MODULE: hail_diagnostic_fast_mod + +!=============================================================================== +module hail_diagnostic_fast_mod + + USE precision, ONLY: wp + + USE umPrintMgr, ONLY: umPrint, umMessage + + implicit none + private + public :: diagnose_hail_fast + + ! ---- physical constants ---- + real(wp), parameter :: pi = 3.14159265358979_wp + real(wp), parameter :: Rd = 287.05_wp ! J/kg/K, dry air + real(wp), parameter :: Rv = 461.5_wp ! J/kg/K, water vapor + real(wp), parameter :: Lf = 3.34e5_wp ! J/kg, latent heat of fusion + real(wp), parameter :: Lv = 2.5e6_wp ! J/kg, latent heat of vaporization + real(wp), parameter :: Kw0 = 0.56_wp ! W/m/K, thermal conductivity of water (~0C) + real(wp), parameter :: Ka0 = 2.4e-2_wp ! W/m/K, thermal conductivity of air (~0C) + real(wp), parameter :: Dv0 = 2.21e-5_wp ! m2/s, vapor diffusivity in air (~0C, 1000hPa) + real(wp), parameter :: T0 = 273.15_wp ! K, melting point + real(wp), parameter :: es0 = 611.2_wp ! Pa, saturation vapor pressure at 0C + real(wp), parameter :: Pr_no = 0.71_wp ! Prandtl number of air + real(wp), parameter :: Sc_no = 0.60_wp ! Schmidt number of air + real(wp), parameter :: rho0_ref = 1.2_wp ! kg/m3, reference air density for V(D) + real(wp), parameter :: visc = 1.7e-5_wp ! visc kg/m/s + + ! ---- terminal velocity power law: V = a_v * (D_cm)^b_v * sqrt(rho0_ref/rho_air) ---- + real(wp), parameter :: a_v = 14.0_wp ! for cm; MASON 1971 4.41(Dmm)0.5 + real(wp), parameter :: b_v = 0.5_wp + + ! ---- minimum thresholds ---- + real(wp), parameter :: qg_min = 1.0e-6_wp ! kg/kg, minimum graupel to bother + +contains + + subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, & + D_crit, D0_thresh, D_max_sfc, precip_rate, & + hail_flag, ierr ) + + USE distributions, ONLY: query_distributions, dist_lambda, dist_mu, dist_n0 + USE mphys_parameters, ONLY: graupel_params + USE passive_fields, ONLY: Tdegk, dz, pressure, rho + use mphys_switches, only: i_qv, i_qg + + real(wp), parameter :: Nthresh=1e-6 + + integer, intent(in) :: nz, nq + integer, intent(in) :: ixy_inner + + real(wp), intent(in) :: qfields(nz,nq) + real(wp), intent(in) :: cffields(nz,nq) + + real(wp), intent(out) :: D_crit ! melting-level critical diameter [m] + real(wp), intent(out) :: D0_thresh ! melting-level size at the concentration threshold [m] + real(wp), intent(out) :: D_max_sfc ! largest hailstone diameter at surface [m] + real(wp), intent(out) :: precip_rate ! surface hail mass flux [kg/m2/s] + ! (x3600 for mm/hr liquid-equivalent) + real(wp), intent(out) :: hail_flag ! .true. if any hail reaches the surface + integer, intent(out) :: ierr ! 0=ok, 1=no graupel, 2=no melting level + + real(wp) :: rhoa(nz), qv(nz), t(nz), p(nz), qg(nz), dz_in(nz) + real(wp) :: qg_ml, ng_ml, t_ml, p_ml, qv_ml, z_ml, rho_air_ml, rho_air_sfc + real(wp) :: N_t, M_t, lam, N0, Mtotal, n_exp + real(wp) :: rho_g, mu_g + integer :: k_ml, nn, k + real(wp) :: a, a0, v, da + + integer, parameter :: nh = 7 + real(wp), parameter :: rhoi = 900.0 + real(wp) :: haild(nh) + real(wp) :: haildf(nh) + real(wp) :: haildedge(nh+1) + real(wp) :: hailprecipf(nh) + 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 + + rho_g=graupel_params%density + qv=qfields(:, i_qv) + qg=qfields(:, i_qg) + t(:)=TdegK(:,ixy_inner) + p(:)=pressure(:,ixy_inner) + rhoa(:)=rho(:,ixy_inner) + dz_in(:)=dz(:,ixy_inner) + + D_crit = 0.0_wp + D0_thresh = 0.0_wp + D_max_sfc = 0.0_wp + precip_rate = 0.0_wp + hail_flag = 0.0_wp + ierr = 0 + + if (maxval(qg) < qg_min) then + ierr = 1 + return + end if + + + call find_melting_level(nz, dz_in, t, k_ml, z_ml, ierr) + if (ierr /= 0) return + mu_g=dist_mu(k_ml,graupel_params%id) + lam=dist_lambda(k_ml, graupel_params%id) + N_t=dist_n0(k_ml, graupel_params%id) + + do nn=1,nh + haild(nn)=(haildedge(nn+1)+haildedge(nn))/2.0 !initial diameters + haildf(nn)=0.0_wp !final diameter + hailprecipf(nn)=0.0_wp ! precip in bins + hailbin(nn)=(haildedge(nn+1)-haildedge(nn)) ! bin width + 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 + ierr = 2 + return + end if + + do nn=1,nh + 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 !diam in cm + da=mason_melt(a,a0,dz_in(k),t(k)-273.15,v,rhoi) + a=a-da + end do + if ( a .gt. 0.0 .and. hailconc(nn) .gt. Nthresh) then + haildf(nn)=2.0*a + hailprecipf(nn)=hailconc(nn)*v*pi/6.0*haildf(nn)**3*rhoi*rhoa(1) + D_max_sfc=2.0*a + hail_flag=1.0 + else + haildf(nn)=0.0 + hailprecipf(nn)=0.0 + end if + end do + precip_rate=sum(hailprecipf(:)) !kg m-2 s-1 + + + + + end subroutine diagnose_hail_fast + + subroutine find_melting_level(nz, dz, t, k_ml, z_ml, ierr) + + integer, intent(in) :: nz + real(wp), intent(in) :: dz(nz), t(nz) + integer, intent(out) :: k_ml + real(wp), intent(out) :: z_ml + integer, intent(out) :: ierr + integer :: k + real(wp) :: frac + + ierr = 3 + k_ml = -1 + z_ml = -999.0_wp + + if (t(1) >= T0) then + do k = nz-1, 1, -1 + + if (t(k+1) < T0 .and. t(k) >= T0) then + frac = (T0 - t(k)) / (t(k+1) - t(k)) + z_ml = SUM(dz(1:k-1)) + frac*dz(k) + k_ml = k + ierr = 0 + + return + end if + end do + end if + + if (t(1) < T0) then + ! surface already below freezing -- no melting layer above the surface + z_ml = 0.0_wp + k_ml = 1 + ierr = 2 + end if + + end subroutine find_melting_level + + !================== + ! Mason melting + !================== + function mason_melt(a,a0,dz,t,v,rhoi) result(da) + real(wp), intent(in) :: a,a0,dz,t,rhoi,v + real(wp) :: da, C, Re, beta + beta=0.0 ! ignore condensation/evap + + Re=v*(2*a)/visc + C=1.6+0.3*Re**0.5 + da=(Kw0*t*dz/Lf/rhoi/v)/ & + ((a0-a)*a/a0+(Kw0/(C*(Ka0+Lv*Dv0*beta)))*a**2/a0) + + + end function mason_melt + + + +end module hail_diagnostic_fast_mod + + + From dda84c28fbbaf71f779a072150d9d786a3abd55a Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Mon, 10 Aug 2026 10:23:05 +0100 Subject: [PATCH 5/7] removed unused variables --- src/generic_diagnostic_variables.F90 | 12 ---------- src/hail_diag_fast_mod.F90 | 34 +++++++++++----------------- src/micro_main.F90 | 7 +++--- 3 files changed, 16 insertions(+), 37 deletions(-) diff --git a/src/generic_diagnostic_variables.F90 b/src/generic_diagnostic_variables.F90 index b259ecc..a859a0f 100644 --- a/src/generic_diagnostic_variables.F90 +++ b/src/generic_diagnostic_variables.F90 @@ -361,8 +361,6 @@ MODULE generic_diagnostic_variables REAL, ALLOCATABLE :: dbz_r(:,:,:) REAL, ALLOCATABLE :: hail_d_max_sfc(:,:) - REAL, ALLOCATABLE :: hail_d_crit(:,:) - REAL, ALLOCATABLE :: hail_d0_thresh(:,:) REAL, ALLOCATABLE :: hail_flag_sfc(:,:) REAL, ALLOCATABLE :: hail_precip_rate(:,:) @@ -646,15 +644,11 @@ SUBROUTINE allocate_diagnostic_space(i_start, i_end, j_start, j_end, k_start, k_ IF ( casdiags % l_hail ) THEN ALLOCATE ( casdiags % hail_d_max_sfc(i_start:i_end, j_start:j_end) ) - ALLOCATE ( casdiags % hail_d_crit(i_start:i_end, j_start:j_end) ) ALLOCATE ( casdiags % hail_flag_sfc(i_start:i_end, j_start:j_end) ) - ALLOCATE ( casdiags % hail_d0_thresh(i_start:i_end, j_start:j_end) ) ALLOCATE ( casdiags % hail_precip_rate(i_start:i_end, j_start:j_end) ) casdiags % hail_d_max_sfc(:,:) = zero_real_wp - casdiags % hail_d_crit(:,:) = zero_real_wp casdiags % hail_flag_sfc(:,:) = zero_real_wp - casdiags % hail_d0_thresh(:,:) = zero_real_wp casdiags % hail_precip_rate(:,:) = zero_real_wp @@ -2623,15 +2617,9 @@ SUBROUTINE deallocate_diagnostic_space() IF ( ALLOCATED ( casdiags % hail_precip_rate ) ) THEN DEALLOCATE( casdiags % hail_precip_rate ) END IF - IF ( ALLOCATED ( casdiags % hail_d0_thresh ) ) THEN - DEALLOCATE( casdiags % hail_d0_thresh ) - END IF IF ( ALLOCATED ( casdiags % hail_flag_sfc ) ) THEN DEALLOCATE( casdiags % hail_flag_sfc ) END IF - IF ( ALLOCATED ( casdiags % hail_d_crit ) ) THEN - DEALLOCATE( casdiags % hail_d_crit ) - END IF IF ( ALLOCATED ( casdiags % hail_d_max_sfc ) ) THEN DEALLOCATE( casdiags % hail_d_max_sfc ) END IF diff --git a/src/hail_diag_fast_mod.F90 b/src/hail_diag_fast_mod.F90 index 41b9da3..8ad4753 100644 --- a/src/hail_diag_fast_mod.F90 +++ b/src/hail_diag_fast_mod.F90 @@ -8,24 +8,19 @@ module hail_diagnostic_fast_mod USE umPrintMgr, ONLY: umPrint, umMessage + USE mphys_constants, ONLY: pi, Rd, Rv, Lf, Lv, ka, Dv, rho0 + implicit none private public :: diagnose_hail_fast ! ---- physical constants ---- - real(wp), parameter :: pi = 3.14159265358979_wp - real(wp), parameter :: Rd = 287.05_wp ! J/kg/K, dry air - real(wp), parameter :: Rv = 461.5_wp ! J/kg/K, water vapor - real(wp), parameter :: Lf = 3.34e5_wp ! J/kg, latent heat of fusion - real(wp), parameter :: Lv = 2.5e6_wp ! J/kg, latent heat of vaporization + real(wp), parameter :: Kw0 = 0.56_wp ! W/m/K, thermal conductivity of water (~0C) - real(wp), parameter :: Ka0 = 2.4e-2_wp ! W/m/K, thermal conductivity of air (~0C) - real(wp), parameter :: Dv0 = 2.21e-5_wp ! m2/s, vapor diffusivity in air (~0C, 1000hPa) real(wp), parameter :: T0 = 273.15_wp ! K, melting point real(wp), parameter :: es0 = 611.2_wp ! Pa, saturation vapor pressure at 0C real(wp), parameter :: Pr_no = 0.71_wp ! Prandtl number of air real(wp), parameter :: Sc_no = 0.60_wp ! Schmidt number of air - real(wp), parameter :: rho0_ref = 1.2_wp ! kg/m3, reference air density for V(D) real(wp), parameter :: visc = 1.7e-5_wp ! visc kg/m/s ! ---- terminal velocity power law: V = a_v * (D_cm)^b_v * sqrt(rho0_ref/rho_air) ---- @@ -37,11 +32,12 @@ module hail_diagnostic_fast_mod contains - subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, & - D_crit, D0_thresh, D_max_sfc, precip_rate, & + subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, & + D_max_sfc, precip_rate, & hail_flag, ierr ) - USE distributions, ONLY: query_distributions, dist_lambda, dist_mu, dist_n0 + USE distributions, ONLY: query_distributions, dist_lambda, & + dist_mu, dist_n0 USE mphys_parameters, ONLY: graupel_params USE passive_fields, ONLY: Tdegk, dz, pressure, rho use mphys_switches, only: i_qv, i_qg @@ -52,19 +48,17 @@ subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, integer, intent(in) :: ixy_inner real(wp), intent(in) :: qfields(nz,nq) - real(wp), intent(in) :: cffields(nz,nq) +! real(wp), intent(in) :: cffields(nz,nq) !not used at the moment. - real(wp), intent(out) :: D_crit ! melting-level critical diameter [m] - real(wp), intent(out) :: D0_thresh ! melting-level size at the concentration threshold [m] real(wp), intent(out) :: D_max_sfc ! largest hailstone diameter at surface [m] real(wp), intent(out) :: precip_rate ! surface hail mass flux [kg/m2/s] ! (x3600 for mm/hr liquid-equivalent) real(wp), intent(out) :: hail_flag ! .true. if any hail reaches the surface integer, intent(out) :: ierr ! 0=ok, 1=no graupel, 2=no melting level - real(wp) :: rhoa(nz), qv(nz), t(nz), p(nz), qg(nz), dz_in(nz) - real(wp) :: qg_ml, ng_ml, t_ml, p_ml, qv_ml, z_ml, rho_air_ml, rho_air_sfc - real(wp) :: N_t, M_t, lam, N0, Mtotal, n_exp + real(wp) :: rhoa(nz), qv(nz), t(nz), qg(nz), dz_in(nz) + real(wp) :: z_ml + real(wp) :: N_t, lam real(wp) :: rho_g, mu_g integer :: k_ml, nn, k real(wp) :: a, a0, v, da @@ -88,8 +82,6 @@ subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, rhoa(:)=rho(:,ixy_inner) dz_in(:)=dz(:,ixy_inner) - D_crit = 0.0_wp - D0_thresh = 0.0_wp D_max_sfc = 0.0_wp precip_rate = 0.0_wp hail_flag = 0.0_wp @@ -124,7 +116,7 @@ subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, cffields, 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 !diam in cm + 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) a=a-da end do @@ -193,7 +185,7 @@ function mason_melt(a,a0,dz,t,v,rhoi) result(da) Re=v*(2*a)/visc C=1.6+0.3*Re**0.5 da=(Kw0*t*dz/Lf/rhoi/v)/ & - ((a0-a)*a/a0+(Kw0/(C*(Ka0+Lv*Dv0*beta)))*a**2/a0) + ((a0-a)*a/a0+(Kw0/(C*(Ka+Lv*Dv*beta)))*a**2/a0) end function mason_melt diff --git a/src/micro_main.F90 b/src/micro_main.F90 index 3dc12a2..a9da331 100644 --- a/src/micro_main.F90 +++ b/src/micro_main.F90 @@ -793,14 +793,13 @@ subroutine shipway_microphysics(il, iu, jl, ju, kl, ku, dt, & if ( casdiags % l_hail ) then - call diagnose_hail_fast( ixy_inner, nz, nq, qfields(:,:,ixy_inner), cffields(:,:,ixy_inner), & - D_crit, D0_thresh, D_max_sfc, precip_rate_hail, hail_flag, ierr ) + call diagnose_hail_fast( ixy_inner, nz, nq, & + qfields(:,:,ixy_inner), & + D_max_sfc, precip_rate_hail, hail_flag, ierr ) casdiags % hail_d_max_sfc(i,j)=D_max_sfc - casdiags % hail_d_crit(i,j)=D_crit casdiags % hail_flag_sfc(i,j)=hail_flag - casdiags % hail_d0_thresh(i,j)=D0_thresh casdiags % hail_precip_rate(i,j)=precip_rate_hail From e1912708bcff6145a616c4ab284291a3dddaf2c5 Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Tue, 11 Aug 2026 09:02:09 +0100 Subject: [PATCH 6/7] simplified hail_diag --- src/hail_diag_fast_mod.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hail_diag_fast_mod.F90 b/src/hail_diag_fast_mod.F90 index 8ad4753..44e8c9f 100644 --- a/src/hail_diag_fast_mod.F90 +++ b/src/hail_diag_fast_mod.F90 @@ -78,7 +78,6 @@ subroutine diagnose_hail_fast( ixy_inner, nz, nq, qfields, & qv=qfields(:, i_qv) qg=qfields(:, i_qg) t(:)=TdegK(:,ixy_inner) - p(:)=pressure(:,ixy_inner) rhoa(:)=rho(:,ixy_inner) dz_in(:)=dz(:,ixy_inner) From 59aa5de0a227aa2c72bdba8d506d442733faa4af Mon Sep 17 00:00:00 2001 From: paulfield2024 Date: Fri, 14 Aug 2026 14:21:46 +0100 Subject: [PATCH 7/7] Add paulfield2024 to CONTRIBUTORS.md --- CONTRIBUTORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index c1d9657..a5524b0 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -5,3 +5,4 @@ | james-bruten-mo | James Bruten | Met Office | 2025-12-09 | | t00sa | Sam Clarke-Green | Met Office | 2026-03-02 | | Pierre-siddall| Pierre Siddall | Met Office | 2026-03-11 | +| paulfield2024 | Paul Field | Met Office | 2026-08-14 |