diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 40d66eec0..b1e5875da 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -993,6 +993,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rainprod, evapprod, & #endif refl_10cm, diagflag, do_radar_ref, & + max_hail_diam_sfc, & vt_dbz_wt, first_time_step, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & @@ -1062,6 +1063,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & GRAUPELNC, GRAUPELNCV REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & refl_10cm + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & + max_hail_diam_sfc REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & vt_dbz_wt LOGICAL, INTENT(IN) :: first_time_step @@ -1679,6 +1682,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & (nsteps>1 .and. istep==nsteps) .or. & (nsteps==1 .and. ndt==1)) THEN + max_hail_diam_sfc(i,j) = hail_mass_99th_percentile(kts, kte, qg1d, t1d, p1d, qv1d) + !> - Call calc_refl10cm() diagflag_present: IF ( PRESENT (diagflag) ) THEN @@ -2464,17 +2469,17 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo - + ! do k = kte, kts, -1 + ! ygra1 = alog10(max(1.E-9, rg(k))) + ! zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 + ! N0_exp = 10.**(zans1) + ! N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + ! lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + ! lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ! ilamg(k) = 1./lamg + ! N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + ! enddo + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -3541,17 +3546,17 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo - + ! do k = kte, kts, -1 + ! ygra1 = alog10(max(1.E-9, rg(k))) + ! zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 + ! N0_exp = 10.**(zans1) + ! N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + ! lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + ! lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ! ilamg(k) = 1./lamg + ! N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + ! enddo + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -6096,16 +6101,17 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo + ! do k = kte, kts, -1 + ! ygra1 = alog10(max(1.E-9, rg(k))) + ! zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 + ! N0_exp = 10.**(zans1) + ! N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + ! lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + ! lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ! ilamg(k) = 1./lamg + ! N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + ! enddo + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -6482,6 +6488,86 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) END SUBROUTINE semi_lagrange_sedim +!>\ingroup aathompson +!! @brief Calculates graupel size distribution parameters +!! +!! Calculates graupel intercept and slope parameters for +!! for a vertical column +!! +!! @param[in] kts integer start index for vertical column +!! @param[in] kte integer end index for vertical column +!! @param[in] rand1 real random number for stochastic physics +!! @param[in] rg real array, size(kts:kte) for graupel mass concentration [kg m^3] +!! @param[out] ilamg double array, size(kts:kte) for inverse graupel slope parameter [m] +!! @param[out] N0_g double array, size(kts:kte) for graupel intercept paramter [m-4] +subroutine graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) + + implicit none + + integer, intent(in) :: kts, kte + real, intent(in) :: rand1 + real, intent(in) :: rg(:) + double precision, intent(out) :: ilamg(:), N0_g(:) + + integer :: k + real :: ygra1, zans1 + double precision :: N0_exp, lam_exp, lamg + + do k = kte, kts, -1 + ygra1 = alog10(max(1.e-9, rg(k))) + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 + N0_exp = 10.**(zans1) + N0_exp = max1(dble(gonv_min), min(N0_exp, dble(gonv_max))) + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ilamg(k) = 1./lamg + N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + enddo + +end subroutine graupel_psd_parameters + +!>\ingroup aathompson +!! @brief Calculates graupel/hail maximum diameter +!! +!! Calculates graupel/hail maximum diameter (currently the 99th percentile of mass distribtuion) +!! for a vertical column +!! +!! @param[in] kts integer start index for vertical column +!! @param[in] kte integer end index for vertical column +!! @param[in] qg real array, size(kts:kte) for graupel mass mixing ratio [kg kg^-1] +!! @param[in] temperature double array, size(kts:kte) temperature [K] +!! @param[in] pressure double array, size(kts:kte) pressure [Pa] +!! @param[in] qv real array, size(kts:kte) water vapor mixing ratio [kg kg^-1] +!! @param[out] max_hail_diam real maximum hail diameter [m] +function hail_mass_99th_percentile(kts, kte, qg, temperature, pressure, qv) result(max_hail_diam) + + implicit none + + integer, intent(in) :: kts, kte + real, intent(in) :: qg(:), temperature(:), pressure(:), qv(:) + real :: max_hail_diam + + integer :: k + real :: rho(kts:kte), rg(kts:kte), max_hail_column(kts:kte) + double precision :: ilamg(kts:kte), N0_g(kts:kte) + real, parameter :: random_number = 0. + + max_hail_column = 0. + rg = 0. + do k = kts, kte + rho(k) = 0.622*pressure(k)/(R*temperature(k)*(max(1.e-10, qv(k))+0.622)) + if (qg(k) .gt. R1) then + rg(k) = qg(k)*rho(k) + endif + enddo + + call graupel_psd_parameters(kts, kte, random_number, rg, ilamg, N0_g) + + where(rg .gt. 1.e-9) max_hail_column = 10.05 * ilamg + max_hail_diam = max_hail_column(kts) + +end function hail_mass_99th_percentile + !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson diff --git a/physics/module_mp_thompson_make_number_concentrations.F90 b/physics/module_mp_thompson_make_number_concentrations.F90 index 72a1055dd..496942f21 100644 --- a/physics/module_mp_thompson_make_number_concentrations.F90 +++ b/physics/module_mp_thompson_make_number_concentrations.F90 @@ -4,6 +4,7 @@ !>\ingroup aathompson module module_mp_thompson_make_number_concentrations + use machine, only : kind_phys use physcons, only: PI => con_pi implicit none @@ -35,7 +36,7 @@ elemental real function make_IceNumber (Q_ice, temp) !IMPLICIT NONE REAL, PARAMETER:: Ice_density = 890.0 !REAL, PARAMETER:: PI = 3.1415926536 - real, intent(in):: Q_ice, temp + real(kind_phys), intent(in):: Q_ice, temp integer idx_rei real corr, reice, deice double precision lambda @@ -134,7 +135,7 @@ elemental real function make_DropletNumber (Q_cloud, qnwfa) !IMPLICIT NONE - real, intent(in):: Q_cloud, qnwfa + real(kind_phys), intent(in):: Q_cloud, qnwfa !real, parameter:: PI = 3.1415926536 real, parameter:: am_r = PI*1000./6. @@ -173,7 +174,7 @@ elemental real function make_RainNumber (Q_rain, temp) IMPLICIT NONE - real, intent(in):: Q_rain, temp + real(kind_phys), intent(in):: Q_rain, temp double precision:: lambda, N0, qnr !real, parameter:: PI = 3.1415926536 real, parameter:: am_r = PI*1000./6. diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index c456e87cd..7b5b83b37 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -329,6 +329,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step, istep, nsteps, & prcp, rain, graupel, ice, snow, sr, & refl_10cm, fullradar_diag, & + max_hail_diam_sfc, & do_radar_ref, aerfld, & mpicomm, mpirank, mpiroot, blkno, & ext_diag, diag3d, reset_diag3d, & @@ -387,6 +388,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & real(kind_phys), intent( out) :: sr(:) ! Radar reflectivity real(kind_phys), intent(inout) :: refl_10cm(:,:) + real(kind_phys), intent(inout) :: max_hail_diam_sfc(:) logical, intent(in ) :: do_radar_ref logical, intent(in) :: sedi_semi integer, intent(in) :: decfl @@ -698,6 +700,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & kme_stoch=kme_stoch, & @@ -738,6 +741,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 5918e4dd9..293dd7625 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -610,6 +610,14 @@ type = real kind = kind_phys intent = out +[max_hail_diam_sfc] + standard_name = max_hail_diameter_sfc + long_name = instantaneous maximum hail diameter at lowest model level + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [fullradar_diag] standard_name = do_full_radar_reflectivity long_name = flag for computing full radar reflectivity