Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MP38 Thompson 2-mom graupel/hail (prev #1667) #1728

Closed
wants to merge 11 commits into from
3 changes: 3 additions & 0 deletions Registry/Registry.EM_COMMON
Original file line number Diff line number Diff line change
Expand Up @@ -2371,6 +2371,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 @@ -2989,6 +2990,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 @@ -3030,6 +3032,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
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))
Copy link
Collaborator

Choose a reason for hiding this comment

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

The NSSL schemes also have qvolg. Would they give a reasonable diagnostic or should this be limited to Thompson?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It should be limited to Thompson, @dudhia. That's because the diagnostic assumes there are 9 possible densities and it's unlikely NSSL makes the same assumptions. Should I add a Case (THOMPSONGH) statement here?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@gthompsnWRF would you be able to confirm if this is necessary, or should this work with NSSL too?

! !$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
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