diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 151a8e9..f587870 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -1,8 +1,9 @@ # Contributors -| GitHub user | Real Name | Affiliation | Date | -| --------------- | ---------------- | ----------- | ---------- | -| 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 | +| GitHub user | Real Name | Affiliation | Date | +| ----------- | --------- | ----------- | ---- | +| 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 | | yaswant | Yaswant Pradhan | Met Office | 2027-07-17 | +| paulfield2024 | Paul Field | Met Office | 2026-08-14 | \ No newline at end of file diff --git a/src/generic_diagnostic_variables.F90 b/src/generic_diagnostic_variables.F90 index cd2b16a..a859a0f 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,10 @@ MODULE generic_diagnostic_variables REAL, ALLOCATABLE :: dbz_l(:,:,:) REAL, ALLOCATABLE :: dbz_r(:,:,:) + REAL, ALLOCATABLE :: hail_d_max_sfc(:,:) + REAL, ALLOCATABLE :: hail_flag_sfc(:,:) + REAL, ALLOCATABLE :: hail_precip_rate(:,:) + ! Process rate diagnostics REAL, ALLOCATABLE :: phomc(:,:,:) REAL, ALLOCATABLE :: pinuc(:,:,:) @@ -631,6 +641,19 @@ 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_flag_sfc(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_flag_sfc(:,:) = 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 +2612,20 @@ 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_flag_sfc ) ) THEN + DEALLOCATE( casdiags % hail_flag_sfc ) + 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 +2677,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/hail_diag_fast_mod.F90 b/src/hail_diag_fast_mod.F90 new file mode 100644 index 0000000..44e8c9f --- /dev/null +++ b/src/hail_diag_fast_mod.F90 @@ -0,0 +1,197 @@ +!=============================================================================== +! MODULE: hail_diagnostic_fast_mod + +!=============================================================================== +module hail_diagnostic_fast_mod + + USE precision, ONLY: wp + + 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 :: Kw0 = 0.56_wp ! W/m/K, thermal conductivity of water (~0C) + 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 :: 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, & + 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) !not used at the moment. + + 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), 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 + + 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) + rhoa(:)=rho(:,ixy_inner) + dz_in(:)=dz(:,ixy_inner) + + 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 * sqrt(rho0) !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*(Ka+Lv*Dv*beta)))*a**2/a0) + + + end function mason_melt + + + +end module hail_diagnostic_fast_mod + + + diff --git a/src/micro_main.F90 b/src/micro_main.F90 index c9e2a37..a9da331 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,22 @@ 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), & + D_max_sfc, precip_rate_hail, hail_flag, ierr ) + + + casdiags % hail_d_max_sfc(i,j)=D_max_sfc + casdiags % hail_flag_sfc(i,j)=hail_flag + 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