diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 4c92f156b3..9cffd556e8 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -2989,6 +2989,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 @@ -3030,6 +3031,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 diff --git a/phys/module_diag_nwp.F b/phys/module_diag_nwp.F index 65b4805620..9879b496a7 100644 --- a/phys/module_diag_nwp.F +++ b/phys/module_diag_nwp.F @@ -38,7 +38,7 @@ 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 & ) @@ -46,11 +46,12 @@ SUBROUTINE diagnostic_output_nwp( & 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 !====================================================================== @@ -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) :: & @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/phys/module_diagnostics_driver.F b/phys/module_diagnostics_driver.F index b2ec19fb59..8baf9dfb06 100644 --- a/phys/module_diagnostics_driver.F +++ b/phys/module_diagnostics_driver.F @@ -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, & @@ -494,6 +494,57 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & ,QR_CURR=moist(ims,kms,jms,P_QR) & ! gthompsn ,NR_CURR=scalar(ims,kms,jms,P_QNR) & ! gthompsn ,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 (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 & diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index a276d8228e..4d4bbb59f9 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -165,7 +165,7 @@ SUBROUTINE microphysics_driver( & ! Framework USE module_state_description, ONLY : & KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME & - ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & + ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, THOMPSONGH, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & ,GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, NSSL_2MOM, NSSL_2MOMCCN, NSSL_2MOMG, MADWRF_MP & ,NSSL_1MOM,NSSL_1MOMLFO, FER_MP_HIRES_ADVECT & ! ,NSSL_3MOM & ,WSM7SCHEME, WDM7SCHEME & @@ -1122,6 +1122,79 @@ SUBROUTINE microphysics_driver( & ELSE CALL wrf_error_fatal ( 'arguments not present for calling thompson_et_al' ) ENDIF +! + CASE (THOMPSONGH) + CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson' ) + IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. & + PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. & + PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. & + PRESENT( QNR_CURR) .AND. PRESENT ( QNI_CURR) .AND. & + PRESENT( QNG_CURR) .AND. PRESENT ( QVOLG_CURR) .AND. & + PRESENT( QNC_CURR) .AND. PRESENT ( QNWFA_CURR) .AND. & + PRESENT( QNIFA_CURR).AND.PRESENT ( QNWFA2D) .AND. & + PRESENT( SNOWNC) .AND. PRESENT ( SNOWNCV) .AND. & + PRESENT( GRAUPELNC).AND. PRESENT ( GRAUPELNCV) .AND. & + PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN +#if ( WRF_CHEM == 1 ) + qv_b4mp(its:ite,kts:kte,jts:jte) = qv_curr(its:ite,kts:kte,jts:jte) + qc_b4mp(its:ite,kts:kte,jts:jte) = qc_curr(its:ite,kts:kte,jts:jte) + qi_b4mp(its:ite,kts:kte,jts:jte) = qi_curr(its:ite,kts:kte,jts:jte) + qs_b4mp(its:ite,kts:kte,jts:jte) = qs_curr(its:ite,kts:kte,jts:jte) +#endif + CALL mp_gt_driver( & + QV=qv_curr, & + QC=qc_curr, & + QR=qr_curr, & + QI=qi_curr, & + QS=qs_curr, & + QG=qg_curr, & + QB=qvolg_curr, & + NI=qni_curr, & + NR=qnr_curr, & + NC=qnc_curr, & + NG=qng_curr, & + NWFA=qnwfa_curr, & + NIFA=qnifa_curr, & + NBCA=qnbca_curr, & + NWFA2D=qnwfa2d, & + NIFA2D=qnifa2d, & + NBCA2D=qnbca2d, & + aer_init_opt=config_flags%aer_init_opt, & + wif_input_opt=config_flags%wif_input_opt, & + TH=th, & + PII=pi_phy, & + P=p, & + W=w, & + DZ=dz8w, & + DT_IN=dt, & + ITIMESTEP=itimestep, & + RAINNC=RAINNC, & + RAINNCV=RAINNCV, & + SNOWNC=SNOWNC, & + SNOWNCV=SNOWNCV, & + GRAUPELNC=GRAUPELNC, & + GRAUPELNCV=GRAUPELNCV, & + SR=SR, & +#if ( WRF_CHEM == 1 ) + WETSCAV_ON=config_flags%wetscav_onoff == 1, & + RAINPROD=rainprod, & + EVAPPROD=evapprod, & +#endif + REFL_10CM=refl_10cm, & + diagflag=diagflag, & + do_radar_ref=do_radar_ref, & + re_cloud=re_cloud, & + re_ice=re_ice, & + re_snow=re_snow, & + has_reqc=has_reqc, & ! G. Thompson + has_reqi=has_reqi, & ! G. Thompson + has_reqs=has_reqs, & ! G. Thompson + IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, & + IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, & + ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte) + ELSE + CALL wrf_error_fatal ( 'arguments not present for calling thompson_et_al' ) + ENDIF ! CASE (THOMPSON) CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson' ) diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index 8c24f3fd42..570d80d99c 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -38,10 +38,13 @@ !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805 !..Last modified: 24 Jan 2018 Aerosol additions to v3.5.1 code 9/2013 !.. Cloud fraction additions 11/2014 part of pre-v3.7 +!.. 2-moment graupel/hail additions Aug. 2018 (released in v4.4) !+---+-----------------------------------------------------------------+ !wrft:model_layer:physics !+---+-----------------------------------------------------------------+ ! + + MODULE module_mp_thompson USE module_wrf_error @@ -57,16 +60,23 @@ MODULE module_mp_thompson LOGICAL, PRIVATE:: is_aerosol_aware = .false. LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true. LOGICAL, PARAMETER, PRIVATE:: homogIce = .true. + LOGICAL, PRIVATE:: is_hail_aware = .false. INTEGER, PARAMETER, PRIVATE:: IFDRY = 0 REAL, PARAMETER, PRIVATE:: T_0 = 273.15 REAL, PARAMETER, PRIVATE:: PI = 3.1415926536 + INTEGER, PRIVATE:: NRHG, idx_bgmax + !..Densities of rain, snow, graupel, and cloud ice. REAL, PARAMETER, PRIVATE:: rho_w = 1000.0 - REAL, PARAMETER, PRIVATE:: rho_s = 100.0 - REAL, PARAMETER, PRIVATE:: rho_g = 500.0 + REAL, PRIVATE:: rho_s = 100.0 REAL, PARAMETER, PRIVATE:: rho_i = 890.0 + INTEGER, PARAMETER, PRIVATE:: NRHG9 = 9 + INTEGER, PARAMETER, PRIVATE:: NRHG1 = 1 ! NRHG will be set to 1 for mp8,28 or 9 for mp38 + REAL, DIMENSION(NRHG9), PARAMETER, PRIVATE:: & + rho_g = (/50., 100., 200., 300., 400., 500., 600., 700., 800./) + INTEGER, PARAMETER :: idx_bg1 = 5 ! representative rho_g index when NRHG=1 !..Prescribed number of cloud droplets. Set according to known data or !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and @@ -119,7 +129,10 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: bm_r = 3.0 REAL, PARAMETER, PRIVATE:: am_s = 0.069 REAL, PARAMETER, PRIVATE:: bm_s = 2.0 - REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0 + REAL, DIMENSION(NRHG9), PARAMETER,PRIVATE:: am_g = (/ & + & PI*rho_g(1)/6.0, PI*rho_g(2)/6.0, PI*rho_g(3)/6.0, & + & PI*rho_g(4)/6.0, PI*rho_g(5)/6.0, PI*rho_g(6)/6.0, & + & PI*rho_g(7)/6.0, PI*rho_g(8)/6.0, PI*rho_g(9)/6.0 /) REAL, PARAMETER, PRIVATE:: bm_g = 3.0 REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0 REAL, PARAMETER, PRIVATE:: bm_i = 3.0 @@ -133,9 +146,19 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: av_s = 40.0 REAL, PARAMETER, PRIVATE:: bv_s = 0.55 REAL, PARAMETER, PRIVATE:: fv_s = 100.0 - REAL, PARAMETER, PRIVATE:: av_g = 442.0 - REAL, PARAMETER, PRIVATE:: bv_g = 0.89 - REAL, PARAMETER, PRIVATE:: av_i = 1847.5 + REAL, PARAMETER, PRIVATE:: av_g_old = 442.0 + REAL, PARAMETER, PRIVATE:: bv_g_old = 0.89 + REAL, DIMENSION(NRHG9), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation + & av_g = (/ 45.9173813, 67.0867386, 98.0158463, 122.353378, & + & 143.204224, 161.794724, 178.762115, 194.488785, & + & 209.225876/) + REAL, DIMENSION(NRHG9), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation + & bv_g = (/0.640961647, 0.640961647, 0.640961647, 0.640961647, & + & 0.640961647, 0.640961647, 0.640961647, 0.640961647, & + & 0.640961647/) + REAL, PARAMETER, PRIVATE:: a_coeff = 0.47244157 + REAL, PARAMETER, PRIVATE:: b_coeff = 0.54698726 + REAL, PARAMETER, PRIVATE:: av_i = 1493.9 REAL, PARAMETER, PRIVATE:: bv_i = 1.0 REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8 REAL, PARAMETER, PRIVATE:: bv_c = 2.0 @@ -349,7 +372,7 @@ MODULE module_mp_thompson !.. allocatable (2009Jun12, J. Michalakes). INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4 - REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: & tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, & tnr_racg, tnr_gacr REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & @@ -377,10 +400,11 @@ MODULE module_mp_thompson REAL, PRIVATE:: oig1, oig2, obmi REAL, DIMENSION(13), PRIVATE:: cre, crg REAL, PRIVATE:: ore1, org1, org2, org3, obmr - REAL, DIMENSION(18), PRIVATE:: cse, csg + REAL, DIMENSION(17), PRIVATE:: cse, csg REAL, PRIVATE:: oams, obms, ocms - REAL, DIMENSION(12), PRIVATE:: cge, cgg - REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg + REAL, DIMENSION(12,NRHG9), PRIVATE:: cge, cgg + REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, obmg + REAL, DIMENSION(NRHG9), PRIVATE:: oamg, ocmg !..Declaration of precomputed constants in various rate eqns. REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi @@ -397,7 +421,7 @@ MODULE module_mp_thompson CONTAINS - SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & + SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, ng, & nwfa, nifa, nbca, wif_input_opt, frc_urb2d, & dx, dy, is_start, & ids, ide, jds, jde, kds, kde, & @@ -411,6 +435,9 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & its,ite, jts,jte, kts,kte REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: hgt +!..OPTIONAL variables that control application of hail-aware scheme + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: ng + !..OPTIONAL variables that control application of aerosol-aware scheme REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa, nbca @@ -427,6 +454,18 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1 LOGICAL:: micro_init, has_CCN, has_IN +!+---+ + + if (PRESENT(ng)) then + is_hail_aware = .TRUE. + NRHG = NRHG9 + else + NRHG = NRHG1 + av_g(idx_bg1) = av_g_old + bv_g(idx_bg1) = bv_g_old + endif + idx_bgmax = MAX(NRHG, idx_bg1) + is_aerosol_aware = .FALSE. micro_init = .FALSE. has_CCN = .FALSE. @@ -566,15 +605,15 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & !..Allocate space for lookup tables (J. Michalakes 2009Jun08). if (.NOT. ALLOCATED(tcg_racg) ) then - ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) + ALLOCATE(tcg_racg(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) micro_init = .TRUE. endif - if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) - if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) - if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) - if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) - if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,NRHG,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) @@ -627,7 +666,7 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & !..Compute min ice diam from mass, min snow/graupel mass from diam. D0i = (xm0i/am_i)**(1./bm_i) xm0s = am_s * D0s**bm_s - xm0g = am_g * D0g**bm_g + xm0g = am_g(idx_bgmax) * D0g**bm_g !..These constants various exponents and gamma() assoc with cloud, !.. rain, snow, and graupel. @@ -702,37 +741,43 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & cse(14) = bm_s + bv_s cse(15) = mu_s + 1. cse(16) = 1.0 + (1.0 + bv_s)/2. - cse(17) = cse(16) + mu_s + 1. - cse(18) = bv_s + mu_s + 3. - do n = 1, 18 + cse(17) = bm_s + bv_s + 2. + do n = 1, 17 csg(n) = WGAMMA(cse(n)) enddo oams = 1./am_s obms = 1./bm_s ocms = oams**obms - cge(1) = bm_g + 1. - cge(2) = mu_g + 1. - cge(3) = bm_g + mu_g + 1. - cge(4) = bm_g*2. + mu_g + 1. - cge(5) = bm_g*2. + mu_g + bv_g + 1. - cge(6) = bm_g + mu_g + bv_g + 1. - cge(7) = bm_g + mu_g + bv_g + 2. - cge(8) = bm_g + mu_g + bv_g + 3. - cge(9) = mu_g + bv_g + 3. - cge(10) = mu_g + 2. - cge(11) = 0.5*(bv_g + 5. + 2.*mu_g) - cge(12) = 0.5*(bv_g + 5.) + mu_g + cge(1,:) = bm_g + 1. + cge(2,:) = mu_g + 1. + cge(3,:) = bm_g + mu_g + 1. + cge(4,:) = bm_g*2. + mu_g + 1. + cge(10,:) = mu_g + 2. + cge(12,:) = bm_g*0.5 + mu_g + 1. + do m = 1, NRHG + cge(5,m) = bm_g*2. + mu_g + bv_g(m) + 1. + cge(6,m) = bm_g + mu_g + bv_g(m) + 1. + cge(7,m) = bm_g*0.5 + mu_g + bv_g(m) + 1. + cge(8,m) = mu_g + bv_g(m) + 1. ! not used + cge(9,m) = mu_g + bv_g(m) + 3. + cge(11,m) = 0.5*(bv_g(m) + 5. + 2.*mu_g) + enddo + do m = 1, NRHG do n = 1, 12 - cgg(n) = WGAMMA(cge(n)) + cgg(n,m) = WGAMMA(cge(n,m)) + enddo enddo oamg = 1./am_g obmg = 1./bm_g - ocmg = oamg**obmg - oge1 = 1./cge(1) - ogg1 = 1./cgg(1) - ogg2 = 1./cgg(2) - ogg3 = 1./cgg(3) + do m = 1, NRHG + oamg(m) = 1./am_g(m) + ocmg(m) = oamg(m)**obmg + enddo + oge1 = 1./cge(1,1) + ogg1 = 1./cgg(1,1) + ogg2 = 1./cgg(2,1) + ogg3 = 1./cgg(3,1) !+---+-----------------------------------------------------------------+ !..Simplify various rate eqns the best we can now. @@ -744,7 +789,7 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & t2_qr_qi = PI*.25*am_r*av_r * crg(8) !..Graupel collecting cloud water - t1_qg_qc = PI*.25*av_g * cgg(9) +! t1_qg_qc = PI*.25*av_g * cgg(9) !..Snow collecting cloud water t1_qs_qc = PI*.25*av_s @@ -765,12 +810,12 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s) !..Sublimation/depositional growth of graupel - t1_qg_sd = 0.86 * cgg(10) - t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11) + t1_qg_sd = 0.86 * cgg(10,1) +! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11) !..Melting of graupel - t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10) - t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11) + t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10,1) +! t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11) !..Constants for helping find lookup table indexes. nic2 = NINT(ALOG10(r_c(1))) @@ -857,14 +902,16 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & do m = 1, ntb_r do k = 1, ntb_r1 + do n = 1, NRHG do j = 1, ntb_g do i = 1, ntb_g1 - tcg_racg(i,j,k,m) = 0.0d0 - tmr_racg(i,j,k,m) = 0.0d0 - tcr_gacr(i,j,k,m) = 0.0d0 - tmg_gacr(i,j,k,m) = 0.0d0 - tnr_racg(i,j,k,m) = 0.0d0 - tnr_gacr(i,j,k,m) = 0.0d0 + tcg_racg(i,j,n,k,m) = 0.0d0 + tmr_racg(i,j,n,k,m) = 0.0d0 + tcr_gacr(i,j,n,k,m) = 0.0d0 + tmg_gacr(i,j,n,k,m) = 0.0d0 + tnr_racg(i,j,n,k,m) = 0.0d0 + tnr_gacr(i,j,n,k,m) = 0.0d0 + enddo enddo enddo enddo @@ -985,7 +1032,7 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & xam_s = am_s xbm_s = bm_s xmu_s = mu_s - xam_g = am_g + xam_g = am_g(idx_bg1) xbm_g = bm_g xmu_g = mu_g call radar_init @@ -1020,7 +1067,7 @@ END SUBROUTINE thompson_init !+---+-----------------------------------------------------------------+ !..This is a wrapper routine designed to transfer values from 3D to 1D. !+---+-----------------------------------------------------------------+ - SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & + SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, qb, ni, nr, nc, ng,& nwfa, nifa, nbca, nwfa2d, nifa2d, nbca2d, & aer_init_opt, wif_input_opt, & th, pii, p, w, dz, dt_in, itimestep, & @@ -1046,7 +1093,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & qv, qc, qr, qi, qs, qg, ni, nr, th REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - nc, nwfa, nifa, nbca + nc, nwfa, nifa, nbca, qb, ng REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, & nbca2d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & @@ -1073,8 +1120,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !..Local variables REAL, DIMENSION(kts:kte):: & - qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, nbca1d, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & + ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d,& t1d, p1d, w1d, dz1d, rho, dBZ REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d #if ( WRF_CHEM == 1 ) @@ -1086,6 +1133,9 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max REAL:: nwfa1 INTEGER:: i, j, k + REAL:: ygra1, zans1, mvd_r + DOUBLE PRECISION:: lamg, lam_exp, lamr, N0_min, N0_exp + INTEGER:: k_0, k_inj INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr @@ -1202,12 +1252,39 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nifa1d(k) = naIN1*0.01/rho(k) nbca1d(k) = 5.55E6/rho(k) enddo - nwfa1 = 11.1E6 + endif + +!..If not the variable-density graupel-hail hybrid, then set the vol mixing +!.. ratio to mass mixing ratio divided by constant density (500kg/m3) value. + + if (is_hail_aware) then + do k = kts, kte + ng1d(k) = ng(i,k,j) + qb1d(k) = qb(i,k,j) + enddo + else + + do k = kte, kts, -1 + if (qg1d(k).gt.R1) then + ygra1 = alog10(max(1.E-9, qg1d(k)*rho(k))) + zans1 = 3.0 + 2./7.*(ygra1+8.) + zans1 = MAX(2., MIN(zans1, 6.)) + N0_exp = 10.**(zans1) + lam_exp = (N0_exp*am_g(idx_bg1)*cgg(1,1)/(rho(k)*qg1d(k)))**oge1 + lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg + ng1d(k) = cgg(2,1)*ogg3*rho(k)*qg1d(k)*lamg**bm_g / am_g(idx_bg1) + ng1d(k) = MAX(R2, ng1d(k)/rho(k)) + qb1d(k) = qg1d(k)/rho_g(idx_bg1) + else + ng1d(k) = 0 + qb1d(k) = 0 + endif + enddo endif call mp_thompson(aer_init_opt, wif_input_opt, & - qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, & + nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, & pptrain, pptsnow, pptgraul, pptice, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod1d, evapprod1d, & @@ -1262,6 +1339,12 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & endif enddo endif + if (is_hail_aware) then + do k = kts, kte + ng(i,k,j) = ng1d(k) + qb(i,k,j) = qb1d(k) + enddo + endif do k = kts, kte qv(i,k,j) = qv1d(k) @@ -1373,7 +1456,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ENDIF call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, kts, kte, i, j, kediagloc) + ng1d, qb1d, t1d, p1d, dBZ, kts, kte, i, j, kediagloc) do k = kts, kte refl_10cm(i,k,j) = MAX(-35., dBZ(k)) enddo @@ -1428,8 +1511,9 @@ END SUBROUTINE mp_gt_driver !+---+-----------------------------------------------------------------+ ! subroutine mp_thompson (aer_init_opt, wif_input_opt, & - qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dzq, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & + ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, & + t1d, p1d, w1d, dzq, & pptrain, pptsnow, pptgraul, pptice, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod, evapprod, & @@ -1442,8 +1526,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt INTEGER, INTENT(IN):: kts, kte, ii, jj REAL, DIMENSION(kts:kte), INTENT(INOUT):: & - qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & + ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN):: dt @@ -1454,8 +1538,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & #endif !..Local variables - REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & - qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten, & + REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, qrten, & + qsten, qgten, qbten, niten, nrten, ncten, ngten, nwfaten, nifaten, & nbcaten DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd @@ -1484,13 +1568,16 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & prs_ide DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, & - prg_gcw, prg_rci, prg_rcs, & - prg_rcg, prg_ihm + prg_gcw, prg_rci, prg_rcs, prg_rcg, prg_ihm, & + png_rcs, png_rcg, png_scw, png_gde, & + pbg_scw, pbg_rfz, pbg_gcw, pbg_rci, pbg_rcs, pbg_rcg, & + pbg_sml, pbg_gml DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 - REAL, DIMENSION(kts:kte):: temp, pres, qv - REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa, nbca + REAL, DIMENSION(kts:kte):: temp, twet, pres, qv + REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, rb + REAL, DIMENSION(kts:kte):: ni, nr, nc, ns, ng, nwfa, nifa, nbca REAL, DIMENSION(kts:kte):: rho, rhof, rhof2 REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati @@ -1498,11 +1585,12 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & tcond, lvap, ocp, lvt2 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g - REAL, DIMENSION(kts:kte):: mvd_r, mvd_c + DOUBLE PRECISION:: N0_melt + REAL, DIMENSION(kts:kte):: mvd_r, mvd_c, mvd_g REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, & - smoc, smod, smoe, smof + smoc, smod, smoe, smof, smog - REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c + REAL, DIMENSION(kts:kte):: sed_r,sed_s,sed_g,sed_i,sed_n,sed_c,sed_b REAL:: rgvm, delta_tp, orho, lfus2 REAL, DIMENSION(5):: onstep @@ -1513,26 +1601,32 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & REAL:: zeta1, zeta, taud, tau REAL:: stoke_r, stoke_s, stoke_g, stoke_i REAL:: vti, vtr, vts, vtg, vtc + REAL:: xrho_g, afall, vtg1, vtg2 + REAL:: bfall = 3*b_coeff - 1 REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, & - vtck, vtnck + vtngk, vtck, vtnck REAL, DIMENSION(kts:kte):: vts_boost + REAL:: M0, slam1, slam2 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow REAL:: a_, b_, loga_, A1, A2, tf REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat - REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr + REAL:: dew_t, Tlcl, The + REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr, xrg, xng, xrb REAL:: xsat, rate_max, sump, ratio REAL:: clap, fcd, dfcd REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl - REAL:: r_frac, g_frac + REAL:: r_frac, g_frac, const_Ri, rime_dens REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr REAL:: Ef_ra, Ef_sa, Ef_ga - REAL:: dtsave, odts, odt, odzq, hgt_agl, SR - REAL:: xslw1, ygra1, zans1, eva_factor - INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq + REAL:: dtsave, odts, odt, odzq, hgt_agl + REAL:: ygra1, zans1, eva_factor + REAL:: SR, melt_f + INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq, k_melting INTEGER, DIMENSION(5):: ksed1 INTEGER:: nir, nis, nig, nii, nic, niin INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, & - idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in + idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in, idx_bgmax + INTEGER, DIMENSION(kts:kte):: idx_bg LOGICAL:: melti, no_micro LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg @@ -1578,6 +1672,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & qrten(k) = 0. qsten(k) = 0. qgten(k) = 0. + ngten(k) = 0. + qbten(k) = 0. niten(k) = 0. nrten(k) = 0. ncten(k) = 0. @@ -1644,6 +1740,20 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & prg_rcs(k) = 0. prg_rcg(k) = 0. prg_ihm(k) = 0. + ! new source/sink terms for 3-moment graupel + png_scw(k) = 0. + png_rcs(k) = 0. + png_rcg(k) = 0. + png_gde(k) = 0. + + pbg_scw(k) = 0. + pbg_rfz(k) = 0. + pbg_gcw(k) = 0. + pbg_rci(k) = 0. + pbg_rcs(k) = 0. + pbg_rcg(k) = 0. + pbg_sml(k) = 0. + pbg_gml(k) = 0. pna_rca(k) = 0. pna_sca(k) = 0. @@ -1666,7 +1776,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & endif #endif -!..Bug fix (2016Jun15), prevent use of uninitialized value(s). +!..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments. do k = kts, kte smo0(k) = 0. smo1(k) = 0. @@ -1676,6 +1786,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & smod(k) = 0. smoe(k) = 0. smof(k) = 0. + smog(k) = 0. + ns(k) = 0. mvd_r(k) = 0. mvd_c(k) = 0. enddo @@ -1802,13 +1914,40 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & endif if (qg1d(k) .gt. R1) then no_micro = .false. - rg(k) = qg1d(k)*rho(k) L_qg(k) = .true. + rg(k) = qg1d(k)*rho(k) + ng(k) = MAX(R2, ng1d(k)*rho(k)) + rb(k) = MAX(qg1d(k)/rho_g(idx_bgmax), qb1d(k)) + rb(k) = MIN(qg1d(k)/rho_g(1), rb(k)) + qb1d(k) = rb(k) + idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,idx_bgmax)) + if (ng(k).le. R2) then + mvd_g(k) = 1.5E-3 + lamg = (3.0 + mu_g + 0.672) / mvd_g(k) + ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) + endif + lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg + mvd_g(k) = (3.0 + mu_g + 0.672) / lamg + if (mvd_g(k) .gt. 25.4E-3) then + mvd_g(k) = 25.4E-3 + lamg = (3.0 + mu_g + 0.672) / mvd_g(k) + ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) + elseif (mvd_g(k) .lt. D0r) then + mvd_g(k) = D0r + lamg = (3.0 + mu_g + 0.672) / mvd_g(k) + ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) + endif else qg1d(k) = 0.0 + ng1d(k) = 0.0 + qb1d(k) = 0.0 + idx_bg(k) = idx_bg1 rg(k) = R1 + ng(k) = R2 + rb(k) = R1/rho(k)/rho_g(idx_bgmax) L_qg(k) = .false. endif + if (.not. is_hail_aware) idx_bg(k) = idx_bg1 enddo !+---+-----------------------------------------------------------------+ @@ -1816,8 +1955,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & ! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj ! CALL wrf_debug(550, mp_debug) ! do k = kts, kte -! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 11(1x,e13.6))') & -! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), nwfa(k), nifa(k) +! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 13(1x,e13.6))') & +! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), ng(k), rb(k), nwfa(k), nifa(k) ! CALL wrf_debug(550, mp_debug) ! enddo ! endif @@ -1829,6 +1968,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978. !+---+-----------------------------------------------------------------+ + k_melting = 0 do k = kts, kte tempc = temp(k) - 273.15 rhof(k) = SQRT(RHO_NOT/rho(k)) @@ -1839,6 +1979,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & qvsi(k) = rsif(pres(k), temp(k)) else qvsi(k) = qvs(k) + k_melting = MAX(k, k_melting) endif satw(k) = qv(k)/qvs(k) sati(k) = qv(k)/qvsi(k) @@ -1857,8 +1998,20 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & vsc2(k) = SQRT(rho(k)/visco(k)) lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936 + twet(k) = temp(k) enddo + if (k_melting .gt. 0) then + do k = kts, k_melting + if (satw(k) .lt. 0.999) then + dew_t = MIN(temp(k)-0.001, t_dew(pres(k), qv(k))) + Tlcl = t_lcl(temp(k), dew_t) + The = theta_e(pres(k), temp(k), qv(k), Tlcl) + twet(k) = MIN(temp(k), compT_fr_The(The, pres(k))) + endif + enddo + endif + !+---+-----------------------------------------------------------------+ !..If no existing hydrometeor species and no chance to initiate ice or !.. condense cloud water, just exit quickly! @@ -1926,6 +2079,14 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc(k) = a_ * smo2(k)**b_ +!..Calculate snow number concentration (explicit integral, not smo0) + M0 = smob(k)/smoc(k) + Mrat = smob(k)*M0*M0*M0 + slam1 = M0 * Lam0 + slam2 = M0 * Lam1 + ns(k) = Mrat*Kap0/slam1 & + + Mrat*Kap1*M0**mu_s*csg(15)/slam2**cse(15) + !..Calculate bv_s+2 (th) moment. Useful for riming. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) & + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 & @@ -1952,6 +2113,19 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16) smof(k) = a_ * smo2(k)**b_ +!..Calculate bm_s + bv_s+2 (th) moment. Useful for riming into graupel. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(17) & + + sa(4)*tc0*cse(17) + sa(5)*tc0*tc0 & + + sa(6)*cse(17)*cse(17) + sa(7)*tc0*tc0*cse(17) & + + sa(8)*tc0*cse(17)*cse(17) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(17)*cse(17)*cse(17) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(17) + sb(4)*tc0*cse(17) & + + sb(5)*tc0*tc0 + sb(6)*cse(17)*cse(17) & + + sb(7)*tc0*tc0*cse(17) + sb(8)*tc0*cse(17)*cse(17) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(17)*cse(17)*cse(17) + smog(k) = a_ * smo2(k)**b_ + enddo !+---+-----------------------------------------------------------------+ @@ -1959,14 +2133,9 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) - zans1 = MAX(2., MIN(zans1, 6.)) - N0_exp = 10.**(zans1) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + N0_g(k) = ng(k)*ogg2*lamg**cge(2,1) enddo endif @@ -2003,6 +2172,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.E6) lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr mvd_c(k) = (3.0+nu_c+0.672) / lamc + mvd_c(k) = MAX(D0c, MIN(mvd_c(k), D0r)) endif !..Autoconversion follows Berry & Reinhardt (1974) with characteristic @@ -2018,7 +2188,7 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & tau = 3.72/(rc(k)*taud) prr_wau(k) = zeta/tau prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) - pnr_wau(k) = prr_wau(k) / (am_r*nu_c*200.*D0r*D0r*D0r) ! RAIN2M + pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*D0r*D0r*D0r) ! RAIN2M pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M endif @@ -2071,8 +2241,15 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & !..fallspeed due to riming of snow do k = kts, kte vts_boost(k) = 1.0 - xDs = 0.0 - if (L_qs(k)) xDs = smoc(k) / smob(k) + orho = 1./rho(k) + + if (L_qs(k)) then + xDs = smoc(k) / smob(k) + rho_s = MAX(rho_g(1), MIN(0.13/xDs, rho_i-100.)) + else + xDs = 0. + rho_s = 100. + endif !..Temperature lookup table indexes. tempc = temp(k) - 273.15 @@ -2187,8 +2364,8 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & idx_g = MAX(1, MIN(idx_g, ntb_g)) lamg = 1./ilamg(k) - lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g - N0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1) + lam_exp = lamg * (cgg(3,1)*ogg2*ogg1)**bm_g + N0_exp = ogg1*rg(k)/am_g(idx_bg(k)) * lam_exp**cge(1,1) nig = NINT(DLOG10(N0_exp)) do nn = nig-1, nig+1 n = nn @@ -2237,9 +2414,12 @@ subroutine mp_thompson (aer_init_opt, wif_input_opt, & !..Graupel collecting cloud water. In CE, assume Dc<