Skip to content

Commit

Permalink
MP38 Thompson 2-mom graupel/hail (prev #1667, #1728) (#1808)
Browse files Browse the repository at this point in the history 
Option to compute two-moment prognostics for graupel/hail

TYPE: enhancement, new feature

KEYWORDS: microphysics, Thompson microphysics, graupel, hail, double-moment

SOURCE: Code developed by Greg Thompson (JCSDA, UCAR) and Anders Jensen (RAL, NCAR).
Implemented in WRF v4.4 by Maria Frediani (RAL, NCAR)

DESCRIPTION OF CHANGES:

This code update includes

a package to compute two-moment prognostics for graupel/hail and a predicted density graupel category (mp_physics=38); an update to the Y-intercept relationship for one-moment graupel; and it replaces air temperature for wet-bulb temperature in riming and mixed phase processes.

Problem:

One-dimensional graupel/hail growth does not couple to the storm dynamics and is insufficient for predicting more detailed microphysical storm characteristics and hazards such as hail size, density, and fall speed, which can be used to provide guidance on the timing and spatial extent of damaging hail.

Sensitivity studies have shown that using a constant intercept parameter and a constant density can significantly constrain predicted hail size: simulated storms either produced only pea-size or baseball-size hail (Gilmore et al. 2004).

Improving the representation of riming and mixed phase processes leads to improvements in predicted storm dynamics and propagation speed through microphysical feedbacks and also improves the spatial distribution and type of precipitation at the surface.

Solution:

Changes related with the new package mp_physics=38 include:

Variable density for graupel (rho_g)
Parameters become a function of rho_g (am_g, av_g, bv_g, cge, cgg, oamg)
Extra dimension in lookup tables to account for graupel variable density (rho_g)
New source/sink terms for 3-moment graupel
Computation of radar reflectivity and nwp diagnostics using graupel volume mixing ratio
Additional modifications affecting mp_physics=8 and mp_physics=28 include:

Fall speed power law relations (av_i from 1847.5 to 1493.9)
Reduced dimension of cse, csg (from 18 to 17)
Use of wet-bulb temperature for riming and mixed-phase process
Modified relationship for the Y-intercept of one-moment graupel to shift the properties of the graupel category to become more hail-like, resulting in a category that represents both graupel and hail.

LIST OF MODIFIED FILES:
M Registry/Registry.EM_COMMON
M phys/module_diag_nwp.F
M phys/module_diagnostics_driver.F
M phys/module_microphysics_driver.F
M phys/module_mp_thompson.F
M phys/module_physics_init.F

TESTS CONDUCTED:
1. The modifications were initially demonstrated using the original development made for WRF v4.0 (Jensen et al 2021, under review, MWR-D-21-0319). The operational mp28, mp28 with modified graupel Y-intercept, and mp38 were evaluated for a case study during the PECAN campaign using observed hail sizes from storm reports and estimated from radar. The evaluation showed clear improvement of the simulated reflectivity values in the upper-levels of discrete storms, coinciding with a significant reduction in the areal extent of graupel aloft, also seen when using the updated one-moment scheme. The two-moment and predicted density graupel scheme was also better able to predict a wide variety of hail sizes at the surface, including large (>2-inch in diameter) hail that was observed during this case.

The implementation for this develop branch (aiming at the release v4.4) was tested using a case study from the Relampago campaign and results from mp28, mp38-v4.0, mp38-develop-v4.4 were compared. This comparison indicates that the implementation was successful.

2. It passed regression tests.

RELEASE NOTES: A package to compute two-moment prognostics for graupel/hail and a predicted density graupel category is added in the Thompson scheme (mp_physics=38); Other updates to the scheme include a change to the Y-intercept relationship for one-moment graupel; and replacement of air temperature for wet-bulb temperature in riming and mixed phase processes.

The code requires a data file to run. This data file: qr_acr_qg_mp38V1.dat can be found on NCAR's computer under /glade/work/wrfhelp/WRF_files/, and online at http://www2.mmm.ucar.edu/wrf/src/wrf_files/. If you prefer to compute this file, set namelist write_thompson_mp38table = true. Note that it can take up to 18 min to compute this table using a 12-CPU job, 4 min on 128-CPU, and several hours if computed on a single CPU.
  • Loading branch information
@mefrediani
mefrediani authored Jan 26, 2023
1 parent 57f4a63 commit 0afe119
Show file tree
Hide file tree
Showing 6 changed files with 1,045 additions and 240 deletions.
3 changes: 3 additions & 0 deletions Registry/Registry.EM_COMMON
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2374,6 +2374,7 @@ rconfig real rankine_lid namelist,tc 1 -999. i
rconfig character physics_suite namelist,physics 1 "none" rh "Physics suite selection" "physics suite to use for all domains: CONUS, tropical, or none" "character string"
rconfig logical force_read_thompson namelist,physics 1 .false.
rconfig logical write_thompson_tables namelist,physics 1 .true.
rconfig logical write_thompson_mp38table namelist,physics 1 .false.
rconfig integer mp_physics namelist,physics max_domains -1 irh "mp_physics" "" ""
#rconfig integer milbrandt_ccntype namelist,physics max_domains 0 rh "milbrandt select maritime(1)/continental(2)" "" ""
rconfig real nssl_cccn namelist,physics max_domains 0.5e9 rh "Base CCN concentration for NSSL microphysics" "" ""
Expand Down Expand Up @@ -2995,6 +2996,7 @@ package nssl_2momg mp_physics==22 - moist:qv,qc
package wsm7scheme mp_physics==24 - moist:qv,qc,qr,qi,qs,qg,qh;state:re_cloud,re_ice,re_snow
package wdm7scheme mp_physics==26 - moist:qv,qc,qr,qi,qs,qg,qh;scalar:qnn,qnc,qnr;state:re_cloud,re_ice,re_snow
package thompsonaero mp_physics==28 - moist:qv,qc,qr,qi,qs,qg;scalar:qni,qnr,qnc,qnwfa,qnifa,qnbca;state:re_cloud,re_ice,re_snow,qnwfa2d,qnifa2d,taod5503d,taod5502d
package thompsongh mp_physics==38 - moist:qv,qc,qr,qi,qs,qg;scalar:qni,qnr,qnc,qng,qvolg,qnwfa,qnifa,qnbca;state:re_cloud,re_ice,re_snow,qnwfa2d,qnifa2d,taod5503d,taod5502d
package p3_1category mp_physics==50 - moist:qv,qc,qr,qi;scalar:qni,qnr,qir,qib;state:re_cloud,re_ice,vmi3d,rhopo3d,di3d,refl_10cm,th_old,qv_old
package p3_1category_nc mp_physics==51 - moist:qv,qc,qr,qi;scalar:qnc,qni,qnr,qir,qib;state:re_cloud,re_ice,vmi3d,rhopo3d,di3d,refl_10cm,th_old,qv_old
package p3_2category mp_physics==52 - moist:qv,qc,qr,qi,qi2;scalar:qnc,qni,qnr,qir,qib,qni2,qir2,qib2;state:re_cloud,re_ice,vmi3d,rhopo3d,di3d,vmi3d_2,rhopo3d_2,di3d_2,refl_10cm,th_old,qv_old
Expand Down Expand Up @@ -3036,6 +3038,7 @@ package nssl_1momlfo_dfi mp_physics_dfi==21 - dfi_moist:dfi
package wsm7scheme_dfi mp_physics_dfi==24 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg,dfi_qh;state:dfi_re_cloud,dfi_re_ice,dfi_re_snow
package wdm7scheme_dfi mp_physics_dfi==26 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg,dfi_qh;dfi_scalar:dfi_qnn,dfi_qnc,dfi_qnr;state:dfi_re_cloud,dfi_re_ice,dfi_re_snow
package thompsonaero_dfi mp_physics_dfi==28 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg;dfi_scalar:dfi_qni,dfi_qnr,dfi_qnc,dfi_qnwfa,dfi_qnifa,dfi_qnbca;state:dfi_re_cloud,dfi_re_ice,dfi_re_snow
package thompsongh_dfi mp_physics_dfi==38 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg;dfi_scalar:dfi_qni,dfi_qnr,dfi_qng,dfi_qvolg,dfi_qnc,dfi_qnwfa,dfi_qnifa,dfi_qnbca;state:dfi_re_cloud,dfi_re_ice,dfi_re_snow
package p3_1category_dfi mp_physics_dfi==50 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi;dfi_scalar:dfi_qni,dfi_qnr,dfi_qir,dfi_qib;state:dfi_re_cloud,dfi_re_ice
package p3_1category_nc_dfi mp_physics_dfi==51 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi;dfi_scalar:dfi_qnc,dfi_qni,dfi_qnr,dfi_qir,dfi_qib;state:dfi_re_cloud,dfi_re_ice
package p3_2category_dfi mp_physics_dfi==52 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qi2;dfi_scalar:dfi_qnc,dfi_qni,dfi_qnr,dfi_qir,dfi_qib,dfi_qni2,dfi_qir2,dfi_qib2;state:dfi_re_cloud,dfi_re_ice
Expand Down
112 changes: 96 additions & 16 deletions phys/module_diag_nwp.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,19 +38,20 @@ SUBROUTINE diagnostic_output_nwp( &
,grpl_max,grpl_colint,refd_max,refl_10cm &
,hail_maxk1,hail_max2d &
,qg_curr &
,ng_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn)
,ng_curr,qvolg_curr,qh_curr,nh_curr,qr_curr,nr_curr & ! Optional (gthompsn)
,rho,ph,phb,g &
,max_time_step,adaptive_ts &
)
!----------------------------------------------------------------------

USE module_state_description, ONLY : &
KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, &
WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, &
WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, THOMPSONGH, &
MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, &
NSSL_2MOM, NSSL_2MOMG, NSSL_2MOMCCN, NSSL_1MOM, NSSL_1MOMLFO, &
MILBRANDT2MOM , CAMMGMPSCHEME, FULL_KHAIN_LYNN, MORR_TM_AERO, &
FAST_KHAIN_LYNN_SHPUND !,MILBRANDT3MOM, NSSL_3MOM
USE MODULE_MP_THOMPSON, ONLY: idx_bg1

IMPLICIT NONE
!======================================================================
Expand Down Expand Up @@ -159,7 +160,7 @@ SUBROUTINE diagnostic_output_nwp( &
,ph,phb

REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: &
ng_curr, qh_curr, nh_curr &
ng_curr, qvolg_curr, qh_curr, nh_curr &
,qr_curr, nr_curr

REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
Expand All @@ -177,7 +178,10 @@ SUBROUTINE diagnostic_output_nwp( &
,hail_maxk1,hail_max2d &
,refd_max

REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: temp_qg, temp_ng
REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: &
temp_qg, temp_ng, temp_qb &
,temp_qr, temp_nr
INTEGER, DIMENSION(ims:ime,kms:kme,jms:jme):: idx_bg

INTEGER :: idump

Expand All @@ -191,12 +195,13 @@ SUBROUTINE diagnostic_output_nwp( &
INTEGER, PARAMETER:: ngbins=50
DOUBLE PRECISION, DIMENSION(ngbins+1):: xxDx
DOUBLE PRECISION, DIMENSION(ngbins):: xxDg, xdtg
REAL:: xrho_g, xam_g, xbm_g, xmu_g
REAL:: xrho_g, xam_g, xxam_g, xbm_g, xmu_g
REAL, DIMENSION(9), PARAMETER:: xrho_gh = (/ 50,100,200,300,400,500,600,700,800 /)
REAL, DIMENSION(3):: cge, cgg
DOUBLE PRECISION:: f_d, sum_ng, sum_t, lamg, ilamg, N0_g, lam_exp, N0exp
DOUBLE PRECISION:: lamr, N0min
REAL:: xslw1, ygra1, zans1
INTEGER:: ng, n
REAL:: mvd_r, xslw1, ygra1, zans1
INTEGER:: ng, n, k_0

REAL :: time_from_output
INTEGER :: max_time_step
Expand Down Expand Up @@ -455,6 +460,39 @@ SUBROUTINE diagnostic_output_nwp( &
! !$OMP END PARALLEL DO
ENDIF

IF (PRESENT(qvolg_curr)) THEN
WRITE(outstring,*) 'GT-Debug, this mp scheme, ', mphysics_opt, ' has graupel volume mixing ratio'
CALL wrf_debug (150, TRIM(outstring))
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO k=kms,kme-1
DO i=i_start(ij),i_end(ij)
temp_qb(i,k,j) = MAX(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(9), qvolg_curr(i,k,j))
temp_qb(i,k,j) = MIN(temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(1), temp_qb(i,k,j))
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
ELSE
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO k=kms,kme-1
DO i=i_start(ij),i_end(ij)
idx_bg(i,k,j) = idx_bg1
temp_qb(i,k,j) = temp_qg(i,k,j)/rho(i,k,j)/xrho_gh(idx_bg(i,k,j))
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO
ENDIF


! Calculate the max hail size from graupel/hail parameters in microphysics scheme (gthompsn 12Mar2015)
! First, we do know multiple schemes have assumed inverse-exponential distribution with constant
! intercept parameter and particle density. Please leave next 4 settings alone for common
Expand Down Expand Up @@ -599,20 +637,56 @@ SUBROUTINE diagnostic_output_nwp( &
cgg(n) = WGAMMA(cge(n))
enddo

if (.not.(PRESENT(ng_curr))) then
! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO k=kms,kme-1
DO i=i_start(ij),i_end(ij)
temp_qr(i,k,j) = MAX(1.E-10, qr_curr(i,k,j)*rho(i,k,j))
temp_nr(i,k,j) = MAX(1.E-8, nr_curr(i,k,j)*rho(i,k,j))
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO

! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO i=i_start(ij),i_end(ij)
DO k=kme-1, kms, -1
if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
zans1 = (3.0 + 2./7. * (ALOG10(temp_qg(i,k,j))+8.))
ygra1 = alog10(max(1.E-9, temp_qg(i,k,j)))
zans1 = 3.0 + 2./7.*(ygra1+8.)
zans1 = MAX(2., MIN(zans1, 6.))
N0exp = 10.**zans1
N0exp = 10.**(zans1)
lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g)
N0_g = N0exp/(cgg(2)*lam_exp) * lamg**cge(2)
temp_ng(i,k,j) = N0_g*cgg(2)*lamg**(-cge(2))
temp_ng(i,k,j) = MAX(1.E-8, N0_g*cgg(2)*lamg**(-cge(2)))
ENDDO
ENDDO
ENDDO
ENDDO
! !$OMP END PARALLEL DO

endif


CASE (THOMPSONGH)

scheme_has_graupel = .true.

! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
DO ij = 1 , num_tiles
DO j=j_start(ij),j_end(ij)
DO k=kme-1, kms, -1
DO i=i_start(ij),i_end(ij)
idx_bg(i,k,j) = NINT(temp_qg(i,k,j)/temp_qb(i,k,j) *0.01)+1
idx_bg(i,k,j) = MAX(1, MIN(idx_bg(i,k,j), 9))
ENDDO
ENDDO
ENDDO
Expand Down Expand Up @@ -724,7 +798,13 @@ SUBROUTINE diagnostic_output_nwp( &
DO k=kms,kme-1
DO i=i_start(ij),i_end(ij)
if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
lamg = (xam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g)
IF (PRESENT(qvolg_curr)) THEN
xxam_g = 3.1415926536/6.0*xrho_gh(idx_bg(i,k,j))
IF (xrho_gh(idx_bg(i,k,j)).lt.350.0) CYCLE
ELSE
xxam_g = xam_g
ENDIF
lamg = (xxam_g*cgg(3)/cgg(2)*temp_ng(i,k,j)/temp_qg(i,k,j))**(1./xbm_g)
N0_g = temp_ng(i,k,j)/cgg(2)*lamg**cge(2)
sum_ng = 0.0d0
sum_t = 0.0d0
Expand All @@ -744,10 +824,10 @@ SUBROUTINE diagnostic_output_nwp( &
else
hail_max = 1.E-4
endif
if (hail_max .gt. 1E-2) then
WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000.
CALL wrf_debug (350, TRIM(outstring))
endif
! if (hail_max .gt. 1E-2) then
! WRITE(outstring,*) 'GT-Debug-Hail, ', hail_max*1000.
! CALL wrf_debug (350, TRIM(outstring))
! endif
hail_max_sp = hail_max
if (k.eq.kms) then
hail_maxk1(i,j) = MAX(hail_maxk1(i,j), hail_max_sp)
Expand Down
55 changes: 53 additions & 2 deletions phys/module_diagnostics_driver.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,10 +34,10 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, &
DO_SOLAR_OUTPUT, &
PARAM_FIRST_SCALAR, &
P_QG, P_QH, P_QV, P_QC, P_QI, P_QS, &
P_QNG, P_QH, P_QNH, P_QR, P_QNR, &
P_QNG, P_QH, P_QNH, P_QR, P_QNR, P_QVOLG, &
F_QV, F_QC, F_QI, F_QS, &
KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME, &
WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, &
WSM6SCHEME, ETAMPNEW, THOMPSON, THOMPSONAERO, THOMPSONGH, &
MORR_TWO_MOMENT, GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, &
NSSL_2MOM, NSSL_2MOMCCN, NSSL_1MOM, NSSL_1MOMLFO, &
MILBRANDT2MOM , CAMMGMPSCHEME, FAST_KHAIN_LYNN_SHPUND, FULL_KHAIN_LYNN, &
Expand Down Expand Up @@ -506,6 +506,57 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, &
,ADAPTIVE_TS=config_flags%use_adaptive_time_step &
)

CASE (THOMPSONGH)

CALL diagnostic_output_nwp( &
U=grid%u_2 ,V=grid%v_2 &
,TEMP=grid%t_phy ,P8W=p8w &
,DT=grid%dt ,SBW=config_flags%spec_bdy_width &
,XTIME=grid%xtime &
! Selection flag
,MPHYSICS_OPT=config_flags%mp_physics & ! gthompsn
,GSFCGCE_HAIL=config_flags%gsfcgce_hail & ! gthompsn
,GSFCGCE_2ICE=config_flags%gsfcgce_2ice & ! gthompsn
,MPUSE_HAIL=config_flags%hail_opt & ! gthompsn
,NSSL_ALPHAH=config_flags%nssl_alphah & ! gthompsn
,NSSL_ALPHAHL=config_flags%nssl_alphahl & ! gthompsn
,NSSL_CNOH=config_flags%nssl_cnoh & ! gthompsn
,NSSL_CNOHL=config_flags%nssl_cnohl & ! gthompsn
,NSSL_RHO_QH=config_flags%nssl_rho_qh & ! gthompsn
,NSSL_RHO_QHL=config_flags%nssl_rho_qhl & ! gthompsn
,CURR_SECS2=curr_secs2 &
,NWP_DIAGNOSTICS=config_flags%nwp_diagnostics &
,DIAGFLAG=diag_flag &
,HISTORY_INTERVAL=grid%history_interval &
,ITIMESTEP=grid%itimestep &
,U10=grid%u10,V10=grid%v10,W=grid%w_2 &
,WSPD10MAX=grid%wspd10max &
,UP_HELI_MAX=grid%up_heli_max &
,W_UP_MAX=grid%w_up_max,W_DN_MAX=grid%w_dn_max &
,ZNW=grid%znw,W_COLMEAN=grid%w_colmean &
,NUMCOLPTS=grid%numcolpts,W_MEAN=grid%w_mean &
,GRPL_MAX=grid%grpl_max,GRPL_COLINT=grid%grpl_colint &
,REFD_MAX=grid%refd_max &
,refl_10cm=grid%refl_10cm &
,HAIL_MAXK1=grid%hail_maxk1,HAIL_MAX2D=grid%hail_max2d & ! gthompsn
,QG_CURR=moist(ims,kms,jms,P_QG) &
,QR_CURR=moist(ims,kms,jms,P_QR) & ! gthompsn
,NR_CURR=scalar(ims,kms,jms,P_QNR) & ! gthompsn
,NG_CURR=scalar(ims,kms,jms,P_QNG) & ! THOMPSONGH
,QVOLG_CURR=scalar(ims,kms,jms,P_QVOLG) & ! THOMPSONGH
,RHO=grid%rho,PH=grid%ph_2,PHB=grid%phb,G=g &
! Dimension arguments
,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe &
,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) &
,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) &
,KTS=k_start, KTE=min(k_end,kde-1) &
,NUM_TILES=grid%num_tiles &
,MAX_TIME_STEP=grid%max_time_step &
,ADAPTIVE_TS=config_flags%use_adaptive_time_step &
)

CASE (MORR_TWO_MOMENT, MORR_TM_AERO) ! TWG add

CALL diagnostic_output_nwp( &
Expand Down
Loading

0 comments on commit 0afe119

Please sign in to comment.