Skip to content

Commit

Permalink
Merge branch 'max_hail_mp' into mynnthomp_oct2023
Browse files Browse the repository at this point in the history
  • Loading branch information
@grantfirl
grantfirl committed Oct 30, 2023
2 parents ac108da + ddf6a5c commit 5351c16
Show file tree
Hide file tree
Showing 4 changed files with 134 additions and 35 deletions.
150 changes: 118 additions & 32 deletions physics/module_mp_thompson.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

!+---+-----------------------------------------------------------------+
Expand Down Expand Up @@ -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

!+---+-----------------------------------------------------------------+
Expand Down Expand Up @@ -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

!+---+-----------------------------------------------------------------+
Expand Down Expand Up @@ -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
Expand Down
7 changes: 4 additions & 3 deletions physics/module_mp_thompson_make_number_concentrations.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down
4 changes: 4 additions & 0 deletions physics/mp_thompson.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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, &
Expand Down
8 changes: 8 additions & 0 deletions physics/mp_thompson.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

4 comments on commit 5351c16

@RuiyuSun
Copy link
Collaborator

Choose a reason for hiding this comment

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

Interesting! @joeolson42. I conducted an experiment with rew_min = 7um and compared with previous experiments, p8 and HR1. Please see the a summary here. . Based on the limited test results (from 6 cases) increasing the minimum REW led to a large global averaged downward SW bias.

@RuiyuSun
Copy link
Collaborator

Choose a reason for hiding this comment

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

@joeolson42 I am sorry this is the where I meant to add my comment.

@joeolson42
Copy link
Collaborator

@joeolson42 joeolson42 commented on 5351c16 Oct 31, 2023 via email

Choose a reason for hiding this comment

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

@RuiyuSun
Copy link
Collaborator

Choose a reason for hiding this comment

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

@joeolson42 The plot has been update.

Please sign in to comment.