diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 76f485293d..f21c49f18e 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -3045,7 +3045,8 @@ package ntu mp_physics==56 - moist:qv,qc package etampnew mp_physics==95 - moist:qv,qc,qr,qs;scalar:qt;state:f_ice_phy,f_rain_phy,f_rimef_phy package gsfcgcescheme mp_physics==97 - moist:qv,qc,qr,qi,qs,qg package madwrf_mp mp_physics==96 - moist:qv,qc,qi,qs - +package rcon_mp_scheme mp_physics==29 - 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,cloudnc + package nssl2mconc nssl_2moment_on==1 - scalar:qndrop,qnr,qni,qns,qng;state:re_cloud,re_ice,re_snow package nssl3mg nssl_3moment==1 - scalar:qzr,qzg package nssl3m nssl_3moment==2 - scalar:qzr,qzg,qzh @@ -3097,6 +3098,7 @@ package jensen_ishmael_dfi mp_physics_dfi==55 - dfi_moist:dfi package ntu_dfi mp_physics_dfi==56 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg,dfi_qh;dfi_scalar:dfi_qnc,dfi_qnr,dfi_qni,dfi_qns,dfi_qng,dfi_qnh,dfi_qdcn,dfi_qtcn,dfi_qccn,dfi_qrcn,dfi_qnin,dfi_fi,dfi_fs,dfi_vi,dfi_vs,dfi_vg,dfi_ai,dfi_as,dfi_ag,dfi_ah,dfi_i3m package etampnew_dfi mp_physics_dfi==95 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qs;dfi_scalar:dfi_qt package gsfcgcescheme_dfi mp_physics_dfi==97 - dfi_moist:dfi_qv,dfi_qc,dfi_qr,dfi_qi,dfi_qs,dfi_qg +package rcon_dfi mp_physics_dfi==29 - 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 noprogn progn==0 - - package progndrop progn==1 - scalar:qndrop;dfi_scalar:dfi_qndrop;state:qndropsource @@ -3654,3 +3656,8 @@ rconfig integer windfarm_wake_model namelist,physics max_domai # wake overlap method, M1, M2, M3, M4 [1, 2, 3, 4] rconfig integer windfarm_overlap_method namelist,physics max_domains 4 rh "windfarm_overlap_method" "" "" rconfig real windfarm_deg namelist,physics max_domains 0 - "windfarm_deg" "for windfarm ideal case" "degree" + + +# outputs for RCON model. +state real CLOUDNC ij misc 1 - rh "CLOUDNC" "ACCUMULATED TOTAL GRID SCALE CLOUD PRECIPITATION" "mm" + diff --git a/Registry/registry.var b/Registry/registry.var index 32cc1471db..930d5a59d1 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -596,6 +596,7 @@ package wdm6scheme mp_physics==16 - moist:qv,qc # Note: Options 17, 19, 21, 22 are deprecated but still reserved for compatibility -- for now package nssl_2mom mp_physics==18 - moist:qv,qc,qr,qi,qs,qg package thompsonaero mp_physics==28 - moist:qv,qc,qr,qi,qs,qg +package rcon_mp_scheme mp_physics==29 - moist:qv,qc,qr,qi,qs,qg package p3_1category mp_physics==50 - moist:qv,qc,qr,qi package p3_1category_nc mp_physics==51 - moist:qv,qc,qr,qi package p3_2category mp_physics==52 - moist:qv,qc,qr,qi,qi2 @@ -627,6 +628,7 @@ package wdm6_4dvar mp_physics_4dvar==16 - g_moist:g_q # Note: Options 17, 19, 21, 22 are deprecated but still reserved for compatibility -- for now package nssl_2mom_4dvar mp_physics_4dvar==18 - g_moist:g_qv,g_qc,g_qr,g_qi,g_qs,g_qg,g_qh;a_moist:a_qv,a_qc,a_qr,a_qi,a_qs,a_qg,a_qh package thompsonaero_4dvar mp_physics_4dvar==28 - g_moist:g_qv,g_qc,g_qr,g_qi,g_qs,g_qg;a_moist:a_qv,a_qc,a_qr,a_qi,a_qs,a_qg +package rcon_mp_scheme_4dvar mp_physics_4dvar==29 - g_moist:g_qv,g_qc,g_qr,g_qi,g_qs,g_qg;a_moist:a_qv,a_qc,a_qr,a_qi,a_qs,a_qg package p3_1category_4dvar mp_physics_4dvar==50 - g_moist:g_qv,g_qc,g_qr,g_qi;a_moist:a_qv,a_qc,a_qr,a_qi package p3_1category_nc_4dvar mp_physics_4dvar==51 - g_moist:g_qv,g_qc,g_qr,g_qi;a_moist:a_qv,a_qc,a_qr,a_qi package p3_2category_4dvar mp_physics_4dvar==52 - g_moist:g_qv,g_qc,g_qr,g_qi,g_qi2;a_moist:a_qv,a_qc,a_qr,a_qi,a_qi2 diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index 07d9173e02..914ea3c6a2 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -237,21 +237,27 @@ SUBROUTINE init_domain_rk ( grid & geogrid_flag_error = geogrid_flag_error + 1 END IF - IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. & + IF ( ( (config_flags%mp_physics .EQ. thompsonaero .OR. & + config_flags%mp_physics .EQ. rcon_mp_scheme & + ) ) .AND. & ( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_erod .EQ. 0 ) ) THEN - CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_erod = 0 ' ) + CALL wrf_message ( '----- ERROR: mp=28 or mp=29 AND dust_emis= 1 AND flag_erod = 0 ' ) geogrid_flag_error = geogrid_flag_error + 1 END IF - IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. & + IF ( ( (config_flags%mp_physics .EQ. thompsonaero .OR. & + config_flags%mp_physics .EQ. rcon_mp_scheme & + ) ) .AND. & ( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_clayfrac .EQ. 0 ) ) THEN - CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_clayfrac = 0 ' ) + CALL wrf_message ( '----- ERROR: mp=28 or mp=29 AND dust_emis= 1 AND flag_clayfrac = 0 ' ) geogrid_flag_error = geogrid_flag_error + 1 END IF - IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. & - ( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_sandfrac .EQ. 0 ) ) THEN - CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_sandfrac = 0 ' ) + IF ( ( (config_flags%mp_physics .EQ. thompsonaero .OR. & + config_flags%mp_physics .EQ. rcon_mp_scheme & + ) ) .AND. & + ( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_sandfrac .EQ. 0 ) ) THEN + CALL wrf_message ( '----- ERROR: mp=28 or mp=29 AND dust_emis= 1 AND flag_sandfrac = 0 ' ) geogrid_flag_error = geogrid_flag_error + 1 END IF @@ -2321,10 +2327,12 @@ SUBROUTINE init_domain_rk ( grid & ! QNWFA - Number concentration water-friendly aerosols ! QNIFA - Number concentration ice-friendly aerosols ! QNBCA - Number concentration black carbon aerosols + ! Also used in RCON Microphysics, mp=29 aer_init_opt = config_flags%aer_init_opt - if_thompsonaero_3d: IF (config_flags%mp_physics .EQ. THOMPSONAERO .AND. & + if_thompsonaero_3d: IF ((config_flags%mp_physics .EQ. THOMPSONAERO & + .OR. config_flags%mp_physics .EQ. RCON_MP_SCHEME) .AND. & config_flags%wif_input_opt .GT. 0) THEN select_aer_init_opt_3d: select case (aer_init_opt) @@ -2731,9 +2739,10 @@ SUBROUTINE init_domain_rk ( grid & end select select_aer_init_opt_3d - ELSE IF (config_flags%mp_physics .EQ. THOMPSONAERO .and. & + ELSE IF ((config_flags%mp_physics .EQ. THOMPSONAERO & + .OR. config_flags%mp_physics .EQ. RCON_MP_SCHEME) .and. & config_flags%wif_input_opt .EQ. 0 ) THEN - CALL wrf_error_fatal ('wif_input_opt=0 but mp_physics=28' ) + CALL wrf_error_fatal ('wif_input_opt=0 but mp_physics=28 or mp_physics=29' ) END IF if_thompsonaero_3d !========================================================================================= @@ -4493,7 +4502,8 @@ SUBROUTINE init_domain_rk ( grid & !.. to read biomass burning aerosol emissions !+---+-----------------------------------------------------------------+ - if_thompsonaero_2d: IF (config_flags%mp_physics .EQ. THOMPSONAERO .AND. & + if_thompsonaero_2d: IF ((config_flags%mp_physics .EQ. THOMPSONAERO & + .OR. config_flags%mp_physics .EQ. RCON_MP_SCHEME) .AND. & config_flags%wif_input_opt .GT. 0) THEN select_aer_init_opt_2d: select case (aer_init_opt) @@ -4782,7 +4792,9 @@ SUBROUTINE init_domain_rk ( grid & !+---+-----------------------------------------------------------------+ IF ( config_flags%mp_physics .EQ. THOMPSON .OR. & - config_flags%mp_physics .EQ. THOMPSONAERO ) THEN + config_flags%mp_physics .EQ. THOMPSONAERO .OR. & + config_flags%mp_physics .EQ. RCON_MP_SCHEME & + ) THEN !..As it occurs up above, temporarily utilizing the v_1 variable, !.. to hold temperature, which it does when time_loop=0. diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index c5f47a50a6..251f1ab98a 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -3977,7 +3977,10 @@ END SUBROUTINE CMAQ_DRIVER & ,pert_thom_qc=config_flags%pert_thom_qc & & ,pert_thom_qi=config_flags%pert_thom_qi & & ,pert_thom_qs=config_flags%pert_thom_qs & - & ,pert_thom_ni=config_flags%pert_thom_ni ) + & ,pert_thom_ni=config_flags%pert_thom_ni & + & ,cloudnc=grid%cloudnc & + ) + BENCH_END(micro_driver_tim) diff --git a/main/depend.common b/main/depend.common index 9ca1327f5f..2fa3fa8d02 100644 --- a/main/depend.common +++ b/main/depend.common @@ -971,6 +971,14 @@ module_mp_thompson.o: \ module_mp_radar.o +module_mp_rcon.o: \ + ../frame/module_domain.o \ + ../share/module_model_constants.o \ + ../frame/module_timing.o \ + ../frame/module_wrf_error.o \ + module_mp_radar.o + + module_mp_nssl_2mom.o: \ ../frame/module_wrf_error.o \ ../share/module_model_constants.o @@ -1335,6 +1343,7 @@ module_physics_init.o: \ module_fdda_spnudging.o \ module_fddaobs_rtfdda.o \ module_mp_thompson.o \ + module_mp_rcon.o \ module_mp_gsfcgce.o \ module_mp_gsfcgce_4ice_nuwrf.o \ module_mp_morr_two_moment.o \ @@ -1381,6 +1390,7 @@ module_microphysics_driver.o: \ module_mp_wsm6r.o \ module_mp_fer_hires.o \ module_mp_thompson.o \ + module_mp_rcon.o \ module_mp_gsfcgce.o \ module_mp_gsfcgce_4ice_nuwrf.o \ module_mp_morr_two_moment.o \ diff --git a/main/real_em.F b/main/real_em.F index 796a0ac219..0a0075c100 100644 --- a/main/real_em.F +++ b/main/real_em.F @@ -15,7 +15,7 @@ PROGRAM real_data USE module_configure, ONLY : grid_config_rec_type, model_config_rec, & initial_config, get_config_as_buffer, set_config_as_buffer USE module_timing - USE module_state_description, ONLY : realonly, THOMPSONAERO + USE module_state_description, ONLY : realonly, THOMPSONAERO, RCON_MP_SCHEME #ifdef NO_LEAP_CALENDAR USE module_symbols_util, ONLY: wrfu_cal_noleap #else @@ -802,7 +802,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ALLOCATE ( qbdy3dtemp2(ims:ime,kms:kme,jms:jme) ) ALLOCATE ( mbdy2dtemp2(ims:ime,1:1, jms:jme) ) - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN IF ( ALLOCATED ( qn1bdy3dtemp1 ) ) DEALLOCATE ( qn1bdy3dtemp1 ) IF ( ALLOCATED ( qn2bdy3dtemp1 ) ) DEALLOCATE ( qn2bdy3dtemp1 ) IF ( ALLOCATED ( qn1bdy3dtemp2 ) ) DEALLOCATE ( qn1bdy3dtemp2 ) @@ -885,7 +885,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren END DO END DO - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME).AND. config_flags%aer_init_opt .gt. 0) THEN CALL couple ( grid%mu_2 , grid%mub , qn1bdy3dtemp1 , grid%scalar(:,:,:,P_QNWFA) , 't' , grid%msfty , & grid%c1h, grid%c2h, grid%c1h, grid%c2h, & ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe ) @@ -965,7 +965,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ims , ime , jms , jme , 1 , 1 , & ips , ipe , jps , jpe , 1 , 1 ) - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN CALL stuff_bdy ( qn1bdy3dtemp1 , grid%scalar_bxs(:,:,:,P_QNWFA), grid%scalar_bxe(:,:,:,P_QNWFA), & grid%scalar_bys(:,:,:,P_QNWFA), grid%scalar_bye(:,:,:,P_QNWFA), & 'T' , spec_bdy_width , & @@ -1071,7 +1071,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren END DO END DO - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN CALL couple ( grid%mu_2 , grid%mub , qn1bdy3dtemp2 , grid%scalar(:,:,:,P_QNWFA) , 't' , grid%msfty , & grid%c1h, grid%c2h, grid%c1h, grid%c2h, & ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe ) @@ -1168,7 +1168,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ids , ide , jds , jde , 1 , 1 , & ims , ime , jms , jme , 1 , 1 , & ips , ipe , jps , jpe , 1 , 1 ) - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN CALL stuff_bdytend ( qn1bdy3dtemp2 , qn1bdy3dtemp1 , REAL(interval_seconds) , & grid%scalar_btxs(:,:,:,P_QNWFA), grid%scalar_btxe(:,:,:,P_QNWFA), & grid%scalar_btys(:,:,:,P_QNWFA), grid%scalar_btye(:,:,:,P_QNWFA), & @@ -1306,7 +1306,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren END DO END DO - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN DO j = jps , jpe DO k = kps , kpe DO i = ips , ipe @@ -1383,7 +1383,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ims , ime , jms , jme , 1 , 1 , & ips , ipe , jps , jpe , 1 , 1 ) - IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN + IF ((config_flags%mp_physics.eq.THOMPSONAERO .OR. config_flags%mp_physics.eq.RCON_MP_SCHEME) .AND. config_flags%aer_init_opt .gt. 0) THEN CALL stuff_bdy ( qn1bdy3dtemp1 , grid%scalar_bxs(:,:,:,P_QNWFA), grid%scalar_bxe(:,:,:,P_QNWFA), & grid%scalar_bys(:,:,:,P_QNWFA), grid%scalar_bye(:,:,:,P_QNWFA), & 'T' , spec_bdy_width , & diff --git a/phys/Makefile b/phys/Makefile index 193a02a9c2..1fdac24766 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -97,6 +97,7 @@ MODULES = \ module_mp_etanew.o \ module_mp_fer_hires.o \ module_mp_thompson.o \ + module_mp_rcon.o \ module_fire_emis.o \ module_mp_SBM_polar_radar.o \ module_mp_full_sbm.o \ diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 53514d346e..f07d1de73a 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -161,12 +161,13 @@ SUBROUTINE microphysics_driver( & ,perts_qsnow, perts_ni & ,pert_thom_qv,pert_thom_qc,pert_thom_qi & ,pert_thom_qs,pert_thom_ni & + ,cloudnc & ) ! Framework USE module_state_description, ONLY : & KESSLERSCHEME, LINSCHEME, SBU_YLINSCHEME, WSM3SCHEME, WSM5SCHEME & - ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, THOMPSONGH, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & + ,WSM6SCHEME, ETAMPNEW, FER_MP_HIRES, THOMPSON, THOMPSONAERO, RCON_MP_SCHEME, THOMPSONGH, FAST_KHAIN_LYNN_SHPUND, MORR_TWO_MOMENT & ,GSFCGCESCHEME, WDM5SCHEME, WDM6SCHEME, NSSL_2MOM, MADWRF_MP & ,FER_MP_HIRES_ADVECT & ,WSM7SCHEME, WDM7SCHEME & @@ -220,6 +221,7 @@ SUBROUTINE microphysics_driver( & USE module_mp_wsm7 USE module_mp_etanew USE module_mp_fer_hires + USE module_mp_rcon USE module_mp_thompson USE module_mp_full_sbm #if ( BUILD_SBM_FAST == 1 ) @@ -684,6 +686,7 @@ SUBROUTINE microphysics_driver( & ,GRAUPELNCV & ,HAILNC & ,HAILNCV & + ,CLOUDNC & ,hail_maxk1, hail_max2d #if ( WRF_CHEM == 1) @@ -1025,6 +1028,90 @@ SUBROUTINE microphysics_driver( & CALL wrf_error_fatal ( 'arguments not present for calling mkessler' ) ENDIF #endif +! +! + CASE (RCON_MP_SCHEME) + + CALL wrf_debug ( 100 , 'microphysics_driver: calling RCON' ) + 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( QNC_CURR) .AND. PRESENT ( QNWFA_CURR) .AND. & + PRESENT( QNIFA_CURR).AND.PRESENT ( QNWFA2D) .AND. & + PRESENT( QNIFA2D) .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_rcon_driver( & + QV=qv_curr, & + QC=qc_curr, & + QR=qr_curr, & + QI=qi_curr, & + QS=qs_curr, & + QG=qg_curr, & + NI=qni_curr, & + NR=qnr_curr, & + NC=qnc_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, & + CLOUDNC=CLOUDNC, & + 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, & + ke_diag = ke_diag, & + do_radar_ref=do_radar_ref, & + re_cloud=re_cloud, & + re_ice=re_ice, & + re_snow=re_snow, & + has_reqc=has_reqc, & + has_reqi=has_reqi, & + has_reqs=has_reqs, & + 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) + + + IF (config_flags%aer_fire_emit_opt.gt.0) then + CALL wrf_debug ( 200 , ' call fire_emis_simple_plumerise' ) + CALL fire_emis_simple_plumerise (config_flags%wif_fire_inj, config_flags%aer_fire_emit_opt & + ,z_at_mass, pblh, qnwfa_curr, qnbca_curr & + ,qnocbb2d, qnbcbb2d, dt, ids, ide, jds, jde, kds, kde & + ,ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte ) + ENDIF + ELSE + CALL wrf_error_fatal ( 'arguments not present for calling rcon' ) + ENDIF ! CASE (THOMPSONAERO) if (pert_thom .and. multi_perturb == 1) then diff --git a/phys/module_mp_rcon.F b/phys/module_mp_rcon.F new file mode 100644 index 0000000000..de22441244 --- /dev/null +++ b/phys/module_mp_rcon.F @@ -0,0 +1,6221 @@ +!+---+-----------------------------------------------------------------+ +! +! RELEASED JANUARY 2025 +! +! This is the WRF module for the RCON microphysics scheme, the intent +! of which is to improve warm rain representation within the Thompson-Eidhammer scheme. +! +! RCON is based heavily on the Thompson-Eidhammer scheme with a couple significant +! changes that improve upon the code in module_mp_rcon.F to generate more realistic +! rainfall during warm rain events with additional benefits for cold rain, especially +! warm processes during cold rain events. +! +! Most deviations from the Thompson-Eidhammer scheme are marked with one of the following comments: +! !RC, !RC !\RC (for blocks) +! +! +! Among the most significant changes for rain productions are 1) the use of a wider cloud water DSD +! of lognormal shape instead of the gamma DSD used by the Thompson–Eidhammer parameterization and +! 2) enhancement of the cloud-to-rain autoconversion parameterization to accomodate the new shape. +! The changes here also allow for sedimentation of cloud water within the lowest model layer, which +! effectively creates a drizzle mode in the scheme. One new diagnostic is provided in wrfout, CLOUDNC, +! which is the accumulated cloud water at the surface, i.e. drizzle. +! +! A future version of this scheme will introduce an expanded aerosol table to increase accuracy of cloud representation. +! +! +! Further details and extensive evaluation of the rain/cloud changes is found in the following paper: +! Conrick, R., C. F. Mass, and L. McMurdie, 2023: Improving Simulations of Warm Rain in a Bulk Microphysics Scheme. +! Mon. Wea. Rev., 152, 169–185, https://doi.org/10.1175/MWR-D-23-0035.1. +! +! +! +! NOTE 1: The intent of this scheme is for the Thompson aerosol options to be used. +! Either the climotological aerosol data or the 'default' profile may be used. +! NOTE 2: The Thompson hail scheme is NOT supported at this time. +! +! +! Questions/comments/bugs? +! Contact Robert Conrick at robert.conrick@gmail.com or Clifford Mass at cmass@uw.edu +! +! +! +!+---+-----------------------------------------------------------------+ +!wrft:model_layer:physics +!+---+-----------------------------------------------------------------+ +! + + + MODULE module_mp_rcon + + USE module_wrf_error + USE module_mp_radar +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + USE module_dm, ONLY : wrf_dm_max_real +#endif + USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG + + IMPLICIT NONE + + LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. + 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 + +!..Densities of rain, snow, graupel, and cloud ice. + REAL, PARAMETER, PRIVATE:: rho_w = 1000.0 + REAL, PRIVATE:: rho_s = 100.0 + REAL, PARAMETER, PRIVATE:: rho_i = 890.0 + INTEGER, PARAMETER, PRIVATE:: NRHG = 9 + INTEGER, PARAMETER, PRIVATE:: NRHG1 = 1 + INTEGER :: dimNRHG + + REAL, DIMENSION(NRHG), PARAMETER, PRIVATE:: & + rho_g = (/50., 100., 200., 300., 400., 500., 600., 700., 800./) + INTEGER, PARAMETER :: idx_bg1 = 5 ! index for rhog when mp=8 or 28 + +!..Prescribed number of cloud droplets. Set according to known data or +!.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and +!.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter, +!.. mu_c, calculated based on Nt_c is important in autoconversion +!.. scheme. In 2-moment cloud water, Nt_c represents a maximum of +!.. droplet concentration and nu_c is also variable depending on local +!.. droplet number concentration. + REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6 + REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6 + +!..Declaration of constants for assumed CCN/IN aerosols when none in +!.. the input data. Look inside the init routine for modifications +!.. due to surface land-sea points or vegetation characteristics. + REAL, PARAMETER, PRIVATE:: naIN0 = 1.5E6 + REAL, PARAMETER, PRIVATE:: naIN1 = 0.5E6 + REAL, PARAMETER, PRIVATE:: naCCN0 = 300.0E6 + REAL, PARAMETER, PRIVATE:: naCCN1 = 50.0E6 + REAL, PARAMETER, PRIVATE:: naBC0 = 150.0E6 + REAL, PARAMETER, PRIVATE:: naBC1 = 25.0E6 + +!..Generalized gamma distributions for rain, graupel and cloud ice. +!.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential. + REAL, PARAMETER, PRIVATE:: mu_r = 0.0 + REAL, PARAMETER, PRIVATE:: mu_g = 0.0 + REAL, PARAMETER, PRIVATE:: mu_i = 0.0 + REAL, PRIVATE:: mu_c + +!..Sum of two gamma distrib for snow (Field et al. 2005). +!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3) +!.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)] +!.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively +!.. calculated as function of ice water content and temperature. + REAL, PARAMETER, PRIVATE:: mu_s = 0.6357 + REAL, PARAMETER, PRIVATE:: Kap0 = 490.6 + REAL, PARAMETER, PRIVATE:: Kap1 = 17.46 + REAL, PARAMETER, PRIVATE:: Lam0 = 20.78 + REAL, PARAMETER, PRIVATE:: Lam1 = 3.29 + +!..Y-intercept parameter for graupel is not constant and depends on +!.. mixing ratio. Also, when mu_g is non-zero, these become equiv +!.. y-intercept for an exponential distrib and proper values are +!.. computed based on same mixing ratio and total number concentration. + REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2 + REAL, PARAMETER, PRIVATE:: gonv_max = 1.E6 + +!..Mass power law relations: mass = am*D**bm +!.. Snow from Field et al. (2005), others assume spherical form. + REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0 + REAL, PARAMETER, PRIVATE:: bm_r = 3.0 + REAL, PARAMETER, PRIVATE:: am_s = 0.069 + REAL, PARAMETER, PRIVATE:: bm_s = 2.0 + REAL, DIMENSION(NRHG), 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 + +!..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D) +!.. Rain from Ferrier (1994), ice, snow, and graupel from +!.. Thompson et al (2008). Coefficient fv is zero for graupel/ice. + REAL, PARAMETER, PRIVATE:: av_r = 4854.0 + REAL, PARAMETER, PRIVATE:: bv_r = 1.0 + REAL, PARAMETER, PRIVATE:: fv_r = 195.0 + 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_old = 442.0 + REAL, PARAMETER, PRIVATE:: bv_g_old = 0.89 + REAL, DIMENSION(NRHG), 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(NRHG), 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 + +!..Capacitance of sphere and plates/aggregates: D**3, D**2 + REAL, PARAMETER, PRIVATE:: C_cube = 0.5 + REAL, PARAMETER, PRIVATE:: C_sqrd = 0.15 + +!..Collection efficiencies. Rain/snow/graupel collection of cloud +!.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and +!.. get computed elsewhere because they are dependent on stokes +!.. number. + REAL, PARAMETER, PRIVATE:: Ef_si = 0.05 + REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95 + REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75 + REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95 + +!..Minimum microphys values +!.. R1 value, 1.E-12, cannot be set lower because of numerical +!.. problems with Paul Field's moments and should not be set larger +!.. because of truncation problems in snow/ice growth. + REAL, PARAMETER, PRIVATE:: R1 = 1.E-12 + REAL, PARAMETER, PRIVATE:: R2 = 1.E-6 + REAL, PARAMETER, PRIVATE:: eps = 1.E-15 + +!..Constants in Cooper curve relation for cloud ice number. + REAL, PARAMETER, PRIVATE:: TNO = 5.0 + REAL, PARAMETER, PRIVATE:: ATO = 0.304 + +!..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment. + REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0) + +!..Schmidt number + REAL, PARAMETER, PRIVATE:: Sc = 0.632 + REAL, PRIVATE:: Sc3 + +!..Homogeneous freezing temperature + REAL, PARAMETER, PRIVATE:: HGFR = 235.16 + +!..Water vapor and air gas constants at constant pressure + REAL, PARAMETER, PRIVATE:: Rv = 461.5 + REAL, PARAMETER, PRIVATE:: oRv = 1./Rv + REAL, PARAMETER, PRIVATE:: R = 287.04 + REAL, PARAMETER, PRIVATE:: Cp = 1004.0 + REAL, PARAMETER, PRIVATE:: R_uni = 8.314 ! J (mol K)-1 + + DOUBLE PRECISION, PARAMETER, PRIVATE:: k_b = 1.38065E-23 ! Boltzmann constant [J/K] + DOUBLE PRECISION, PARAMETER, PRIVATE:: M_w = 18.01528E-3 ! molecular mass of water [kg/mol] + DOUBLE PRECISION, PARAMETER, PRIVATE:: M_a = 28.96E-3 ! molecular mass of air [kg/mol] + DOUBLE PRECISION, PARAMETER, PRIVATE:: N_avo = 6.022E23 ! Avogadro number [1/mol] + DOUBLE PRECISION, PARAMETER, PRIVATE:: ma_w = M_w / N_avo ! mass of water molecule [kg] + REAL, PARAMETER, PRIVATE:: ar_volume = 4./3.*PI*(2.5e-6)**3 ! assume radius of 0.025 micrometer, 2.5e-6 cm + +!..Enthalpy of sublimation, vaporization, and fusion at 0C. + REAL, PARAMETER, PRIVATE:: lsub = 2.834E6 + REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6 + REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0 + REAL, PARAMETER, PRIVATE:: olfus = 1./lfus + +!..Ice initiates with this mass (kg), corresponding diameter calc. +!..Min diameters and mass of cloud, rain, snow, and graupel (m, kg). + REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12 + REAL, PARAMETER, PRIVATE:: D0c = 1.E-6 + REAL, PARAMETER, PRIVATE:: D0r = 50.E-6 + REAL, PARAMETER, PRIVATE:: D0s = 300.E-6 + REAL, PARAMETER, PRIVATE:: D0g = 350.E-6 + REAL, PRIVATE:: D0i, xm0s, xm0g + +!..Lookup table dimensions + INTEGER, PARAMETER, PRIVATE:: nbins = 100 + INTEGER, PARAMETER, PRIVATE:: nbc = nbins + INTEGER, PARAMETER, PRIVATE:: nbi = nbins + INTEGER, PARAMETER, PRIVATE:: nbr = nbins + INTEGER, PARAMETER, PRIVATE:: nbs = nbins + INTEGER, PARAMETER, PRIVATE:: nbg = nbins + INTEGER, PARAMETER, PRIVATE:: ntb_c = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_i = 64 + INTEGER, PARAMETER, PRIVATE:: ntb_r = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_s = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_g = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55 + INTEGER, PARAMETER, PRIVATE:: ntb_t = 9 + INTEGER, PRIVATE:: nic1, nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3 + INTEGER, PARAMETER, PRIVATE:: ntb_arc = 7 + INTEGER, PARAMETER, PRIVATE:: ntb_arw = 9 + INTEGER, PARAMETER, PRIVATE:: ntb_art = 7 + INTEGER, PARAMETER, PRIVATE:: ntb_arr = 5 + INTEGER, PARAMETER, PRIVATE:: ntb_ark = 4 + INTEGER, PARAMETER, PRIVATE:: ntb_IN = 55 + INTEGER, PRIVATE:: niIN2 + + DOUBLE PRECISION, DIMENSION(nbins+1):: xDx + DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc + DOUBLE PRECISION, DIMENSION(nbi):: Di, dti + DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr + DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts + DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg + DOUBLE PRECISION, DIMENSION(nbc):: t_Nc + +!..Lookup tables for cloud water content (kg/m**3). + REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: & + r_c = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for cloud ice content (kg/m**3). + REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: & + r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, & + 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, & + 1.e-9,2.e-9,3.e-9,4.e-9,5.e-9,6.e-9,7.e-9,8.e-9,9.e-9, & + 1.e-8,2.e-8,3.e-8,4.e-8,5.e-8,6.e-8,7.e-8,8.e-8,9.e-8, & + 1.e-7,2.e-7,3.e-7,4.e-7,5.e-7,6.e-7,7.e-7,8.e-7,9.e-7, & + 1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3/) + +!..Lookup tables for rain content (kg/m**3). + REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: & + r_r = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for graupel content (kg/m**3). + REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: & + r_g = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for snow content (kg/m**3). + REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: & + r_s = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for rain y-intercept parameter (/m**4). + REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: & + N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & + 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, & + 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, & + 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, & + 1.e10/) + +!..Lookup tables for graupel y-intercept parameter (/m**4). + REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: & + N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & + 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & + 1.e6/) + +!..Lookup tables for ice number concentration (/m**3). + REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: & + Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, & + 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, & + 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & + 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & + 1.e6/) + +!..Aerosol table parameter: Number of available aerosols, vertical +!.. velocity, temperature, aerosol mean radius, and hygroscopicity. + REAL, DIMENSION(ntb_arc), PARAMETER, PRIVATE:: & + ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) + REAL, DIMENSION(ntb_arw), PARAMETER, PRIVATE:: & + ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) + REAL, DIMENSION(ntb_art), PARAMETER, PRIVATE:: & + ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) + REAL, DIMENSION(ntb_arr), PARAMETER, PRIVATE:: & + ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) + REAL, DIMENSION(ntb_ark), PARAMETER, PRIVATE:: & + ta_Ka = (/0.2, 0.4, 0.6, 0.8/) + +!..Lookup tables for IN concentration (/m**3) from 0.001 to 1000/Liter. + REAL, DIMENSION(ntb_IN), PARAMETER, PRIVATE:: & + Nt_IN = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, & + 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, & + 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & + 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & + 1.e6/) + +!..For snow moments conversions (from Field et al. 2005) + REAL, DIMENSION(10), PARAMETER, PRIVATE:: & + sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, & + 0.31255, 0.000204, 0.003199, 0.0, -0.015952/) + REAL, DIMENSION(10), PARAMETER, PRIVATE:: & + sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, & + 0.060366, 0.000079, 0.000594, 0.0, -0.003577/) + +!..Temperatures (5 C interval 0 to -40) used in lookup tables. + REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: & + Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./) + +!..Lookup tables for various accretion/collection terms. +!.. ntb_x refers to the number of elements for rain, snow, graupel, +!.. and temperature array indices. Variables beginning with t-p/c/m/n +!.. represent lookup tables. Save compile-time memory by making +!.. allocatable (2009Jun12, J. Michalakes). + INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8 + INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4 + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: & + tcg_racg, tmr_racg, tcr_gacr, & ! tmg_gacr + tnr_racg, tnr_gacr + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, & + tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, & + tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2 + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + tpi_qcfz, tni_qcfz + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: & + tps_iaus, tni_iaus, tpi_ide + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: & + tpc_wev, tnc_wev + REAL (KIND=R4SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: tnccn_act + +!..Variables holding a bunch of exponents and gamma values (cloud water, +!.. cloud ice, rain, snow, then graupel). + REAL, DIMENSION(5,15), PRIVATE:: cce, ccg + REAL, DIMENSION(15), PRIVATE:: ocg1, ocg2 + REAL, DIMENSION(7), PRIVATE:: cie, cig + REAL, PRIVATE:: oig1, oig2, obmi + REAL, DIMENSION(13), PRIVATE:: cre, crg + REAL, PRIVATE:: ore1, org1, org2, org3, obmr + REAL, DIMENSION(17), PRIVATE:: cse, csg + REAL, PRIVATE:: oams, obms, ocms + REAL, DIMENSION(12,NRHG), PRIVATE:: cge, cgg + REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, obmg + REAL, DIMENSION(NRHG), 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 + REAL:: t1_qr_ev, t2_qr_ev + REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd + REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me + +!+---+ +!+---+-----------------------------------------------------------------+ +!..END DECLARATIONS +!+---+-----------------------------------------------------------------+ +!+---+ +!ctrlL + + CONTAINS + + SUBROUTINE rcon_init(hgt, orho, nwfa2d, nbca2d, ng, & + nwfa, nifa, nbca, wif_input_opt, frc_urb2d, & + dx, dy, is_start, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + its, ite, jts, jte, kts, kte) + + IMPLICIT NONE + + INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + 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 + REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d, nbca2d + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: orho + REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN) :: frc_urb2d + REAL, OPTIONAL, INTENT(IN) :: DX, DY + LOGICAL, OPTIONAL, INTENT(IN) :: is_start + INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt + CHARACTER*256:: mp_debug + + + INTEGER:: i, j, k, l, m, n + REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1 + LOGICAL:: micro_init, has_CCN, has_IN + +!+---+ + + if (PRESENT(ng)) then + is_hail_aware = .TRUE. + dimNRHG = NRHG + else + av_g(idx_bg1) = av_g_old + bv_g(idx_bg1) = bv_g_old + dimNRHG = NRHG1 + endif + + is_aerosol_aware = .FALSE. + micro_init = .FALSE. + has_CCN = .FALSE. + has_IN = .FALSE. + + write(mp_debug,*) ' DEBUG checking column of hgt ', its+1,jts+1 + CALL wrf_debug(250, mp_debug) + do k = kts, kte + write(mp_debug,*) ' DEBUGT k, hgt = ', k, hgt(its+1,k,jts+1) + CALL wrf_debug(250, mp_debug) + enddo + + if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE. + + if (is_aerosol_aware) then + +!..Check for existing aerosol data, both CCN and IN aerosols. If missing +!.. fill in just a basic vertical profile, somewhat boundary-layer following. + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + max_test = wrf_dm_max_real ( MAXVAL(nwfa(its:ite-1,:,jts:jte-1)) ) +#else + max_test = MAXVAL ( nwfa(its:ite-1,:,jts:jte-1) ) +#endif + + if (max_test .lt. eps) then + write(mp_debug,*) ' Apparently there are no initial CCN aerosols.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' checked column at point (i,j) = ', its,jts + CALL wrf_debug(100, mp_debug) + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + if (hgt(i,1,j).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1,j).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) + endif + niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 + nwfa(i,1,j) = naCCN1+naCCN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niCCN3) + z1=hgt(i,2,j)-hgt(i,1,j) + nwfa2d(i,j) = nwfa(i,1,j) * 0.000196 * (50./z1) + do k = 2, kte + nwfa(i,k,j) = naCCN1+naCCN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niCCN3) + enddo + enddo + enddo + else + has_CCN = .TRUE. + write(mp_debug,*) ' Apparently initial CCN aerosols are present.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts)) + CALL wrf_debug(100, mp_debug) + endif + + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + max_test = wrf_dm_max_real ( MAXVAL(nifa(its:ite-1,:,jts:jte-1)) ) +#else + max_test = MAXVAL ( nifa(its:ite-1,:,jts:jte-1) ) +#endif + + if (max_test .lt. eps) then + write(mp_debug,*) ' Apparently there are no initial IN aerosols.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' checked column at point (i,j) = ', its,jts + CALL wrf_debug(100, mp_debug) + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + if (hgt(i,1,j).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1,j).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) + endif + niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 + nifa(i,1,j) = naIN1+naIN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niIN3) + do k = 2, kte + nifa(i,k,j) = naIN1+naIN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niIN3) + enddo + enddo + enddo + else + has_IN = .TRUE. + write(mp_debug,*) ' Apparently initial IN aerosols are present.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts)) + CALL wrf_debug(100, mp_debug) + endif + + + if ( wif_input_opt .eq. 2 ) then + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + max_test = wrf_dm_max_real ( MAXVAL(nbca(its:ite-1,:,jts:jte-1)) ) +#else + max_test = MAXVAL ( nbca(its:ite-1,:,jts:jte-1) ) +#endif + + if (max_test .lt. eps) then + write(mp_debug,*) ' Apparently there are no initial BC aerosols.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' checked column at point (i,j) = ', its,jts + CALL wrf_debug(100, mp_debug) + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + if (hgt(i,1,j).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1,j).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) + endif + niBC3 = -1.0*ALOG(naBC1/naBC0)/h_01 + nbca(i,1,j) = naBC1+naBC0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niBC3) + z1=hgt(i,2,j)-hgt(i,1,j) + nbca2d(i,j) = nbca(i,1,j) * 0.000098 * (50./z1) * (1. + frc_urb2d(i,j)) + do k = 2, kte + nbca(i,k,j) = naBC1+naBC0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niBC3) + enddo + enddo + enddo + else + write(mp_debug,*) ' Apparently initial BC aerosols are present.' + CALL wrf_debug(100, mp_debug) + write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nbca(its,:,jts)) + CALL wrf_debug(100, mp_debug) + endif + + endif + + endif + + +!..Allocate space for lookup tables (J. Michalakes 2009Jun08). + + if (.NOT. ALLOCATED(tcg_racg) ) then + ALLOCATE(tcg_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + micro_init = .TRUE. + endif + + if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + ! if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,dimNRHG,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)) + if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) + + if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,nbc,45,ntb_IN)) + if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,nbc,45,ntb_IN)) + + if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45,ntb_IN)) + if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45,ntb_IN)) + if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45,ntb_IN)) + if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45,ntb_IN)) + + if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1)) + if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1)) + if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1)) + + if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc)) + if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc)) + + if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r)) + if (.NOT. ALLOCATED(tpc_wev)) ALLOCATE(tpc_wev(nbc,ntb_c,nbc)) + if (.NOT. ALLOCATED(tnc_wev)) ALLOCATE(tnc_wev(nbc,ntb_c,nbc)) + + if (.NOT. ALLOCATED(tnccn_act)) & + ALLOCATE(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark)) + + if (micro_init) then + +!..From Martin et al. (1994), assign gamma shape parameter mu for cloud +!.. drops according to general dispersion characteristics (disp=~0.25 +!.. for Maritime and 0.45 for Continental). +!.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime +!.. to 2 for really dirty air. This not used in 2-moment cloud water +!.. scheme and nu_c used instead and varies from 2 to 15 (integer-only). + mu_c = MIN(15., (1000.E6/Nt_c + 2.)) + +!..Schmidt number to one-third used numerous times. + Sc3 = Sc**(1./3.) + +!..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(NRHG) * D0g**bm_g + +!..These constants various exponents and gamma() assoc with cloud, +!.. rain, snow, and graupel. + do n = 1, 15 + cce(1,n) = n + 1. + cce(2,n) = bm_r + n + 1. + cce(3,n) = bm_r + n + 4. + cce(4,n) = n + bv_c + 1. + cce(5,n) = bm_r + n + bv_c + 1. + ccg(1,n) = WGAMMA(cce(1,n)) + ccg(2,n) = WGAMMA(cce(2,n)) + ccg(3,n) = WGAMMA(cce(3,n)) + ccg(4,n) = WGAMMA(cce(4,n)) + ccg(5,n) = WGAMMA(cce(5,n)) + ocg1(n) = 1./ccg(1,n) + ocg2(n) = 1./ccg(2,n) + enddo + + cie(1) = mu_i + 1. + cie(2) = bm_i + mu_i + 1. + cie(3) = bm_i + mu_i + bv_i + 1. + cie(4) = mu_i + bv_i + 1. + cie(5) = mu_i + 2. + cie(6) = bm_i*0.5 + mu_i + bv_i + 1. + cie(7) = bm_i*0.5 + mu_i + 1. + cig(1) = WGAMMA(cie(1)) + cig(2) = WGAMMA(cie(2)) + cig(3) = WGAMMA(cie(3)) + cig(4) = WGAMMA(cie(4)) + cig(5) = WGAMMA(cie(5)) + cig(6) = WGAMMA(cie(6)) + cig(7) = WGAMMA(cie(7)) + oig1 = 1./cig(1) + oig2 = 1./cig(2) + obmi = 1./bm_i + + cre(1) = bm_r + 1. + cre(2) = mu_r + 1. + cre(3) = bm_r + mu_r + 1. + cre(4) = bm_r*2. + mu_r + 1. + cre(5) = mu_r + bv_r + 1. + cre(6) = bm_r + mu_r + bv_r + 1. + cre(7) = bm_r*0.5 + mu_r + bv_r + 1. + cre(8) = bm_r + mu_r + bv_r + 3. + cre(9) = mu_r + bv_r + 3. + cre(10) = mu_r + 2. + cre(11) = 0.5*(bv_r + 5. + 2.*mu_r) + cre(12) = bm_r*0.5 + mu_r + 1. + cre(13) = bm_r*2. + mu_r + bv_r + 1. + do n = 1, 13 + crg(n) = WGAMMA(cre(n)) + enddo + obmr = 1./bm_r + ore1 = 1./cre(1) + org1 = 1./crg(1) + org2 = 1./crg(2) + org3 = 1./crg(3) + + cse(1) = bm_s + 1. + cse(2) = bm_s + 2. + cse(3) = bm_s*2. + cse(4) = bm_s + bv_s + 1. + cse(5) = bm_s*2. + bv_s + 1. + cse(6) = bm_s*2. + 1. + cse(7) = bm_s + mu_s + 1. + cse(8) = bm_s + mu_s + 2. + cse(9) = bm_s + mu_s + 3. + cse(10) = bm_s + mu_s + bv_s + 1. + cse(11) = bm_s*2. + mu_s + bv_s + 1. + cse(12) = bm_s*2. + mu_s + 1. + cse(13) = bv_s + 2. + cse(14) = bm_s + bv_s + cse(15) = mu_s + 1. + cse(16) = 1.0 + (1.0 + bv_s)/2. + 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(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,m) = WGAMMA(cge(n,m)) + enddo + enddo + oamg = 1./am_g + obmg = 1./bm_g + 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. +!+---+-----------------------------------------------------------------+ + +!..Rain collecting cloud water and cloud ice + t1_qr_qc = PI*.25*av_r * crg(9) + t1_qr_qi = PI*.25*av_r * crg(9) + t2_qr_qi = PI*.25*am_r*av_r * crg(8) + +!..Graupel collecting cloud water +! t1_qg_qc = PI*.25*av_g * cgg(9) + +!..Snow collecting cloud water + t1_qs_qc = PI*.25*av_s + +!..Snow collecting cloud ice + t1_qs_qi = PI*.25*av_s + +!..Evaporation of rain; ignore depositional growth of rain. + t1_qr_ev = 0.78 * crg(10) + t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11) + +!..Sublimation/depositional growth of snow + t1_qs_sd = 0.86 + t2_qs_sd = 0.28*Sc3*SQRT(av_s) + +!..Melting of snow + t1_qs_me = PI*4.*C_sqrd*olfus * 0.86 + 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,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,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))) + nii2 = NINT(ALOG10(r_i(1))) + nii3 = NINT(ALOG10(Nt_i(1))) + nir2 = NINT(ALOG10(r_r(1))) + nir3 = NINT(ALOG10(N0r_exp(1))) + nis2 = NINT(ALOG10(r_s(1))) + nig2 = NINT(ALOG10(r_g(1))) + nig3 = NINT(ALOG10(N0g_exp(1))) + niIN2 = NINT(ALOG10(Nt_IN(1))) + +!..Create bins of cloud water (from min diameter up to 100 microns). + Dc(1) = D0c*1.0d0 + dtc(1) = D0c*1.0d0 + do n = 2, nbc + Dc(n) = Dc(n-1) + 1.0D-6 + dtc(n) = (Dc(n) - Dc(n-1)) + enddo + +!..Create bins of cloud ice (from min diameter up to 5x min snow size). + xDx(1) = D0i*1.0d0 + xDx(nbi+1) = 2.0d0*D0s + do n = 2, nbi + xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) & + *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1))) + enddo + do n = 1, nbi + Di(n) = DSQRT(xDx(n)*xDx(n+1)) + dti(n) = xDx(n+1) - xDx(n) + enddo + +!..Create bins of rain (from min diameter up to 5 mm). + xDx(1) = D0r*1.0d0 + xDx(nbr+1) = 0.005d0 + do n = 2, nbr + xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) & + *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1))) + enddo + do n = 1, nbr + Dr(n) = DSQRT(xDx(n)*xDx(n+1)) + dtr(n) = xDx(n+1) - xDx(n) + enddo + +!..Create bins of snow (from min diameter up to 2 cm). + xDx(1) = D0s*1.0d0 + xDx(nbs+1) = 0.02d0 + do n = 2, nbs + xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) & + *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1))) + enddo + do n = 1, nbs + Ds(n) = DSQRT(xDx(n)*xDx(n+1)) + dts(n) = xDx(n+1) - xDx(n) + enddo + +!..Create bins of graupel (from min diameter up to 5 cm). + xDx(1) = D0g*1.0d0 + xDx(nbg+1) = 0.05d0 + do n = 2, nbg + xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) & + *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1))) + enddo + do n = 1, nbg + Dg(n) = DSQRT(xDx(n)*xDx(n+1)) + dtg(n) = xDx(n+1) - xDx(n) + enddo + +!..Create bins of cloud droplet number concentration (1 to 3000 per cc). + xDx(1) = 1.0d0 + xDx(nbc+1) = 3000.0d0 + do n = 2, nbc + xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbc) & + *DLOG(xDx(nbc+1)/xDx(1)) +DLOG(xDx(1))) + enddo + do n = 1, nbc + t_Nc(n) = DSQRT(xDx(n)*xDx(n+1)) * 1.D6 + enddo + nic1 = DLOG(t_Nc(nbc)/t_Nc(1)) + +!+---+-----------------------------------------------------------------+ +!..Create lookup tables for most costly calculations. +!+---+-----------------------------------------------------------------+ + + do m = 1, ntb_r + do k = 1, ntb_r1 + do n = 1, dimNRHG + do j = 1, ntb_g + do i = 1, ntb_g1 + 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 + enddo + + do m = 1, ntb_r + do k = 1, ntb_r1 + do j = 1, ntb_t + do i = 1, ntb_s + tcs_racs1(i,j,k,m) = 0.0d0 + tmr_racs1(i,j,k,m) = 0.0d0 + tcs_racs2(i,j,k,m) = 0.0d0 + tmr_racs2(i,j,k,m) = 0.0d0 + tcr_sacr1(i,j,k,m) = 0.0d0 + tms_sacr1(i,j,k,m) = 0.0d0 + tcr_sacr2(i,j,k,m) = 0.0d0 + tms_sacr2(i,j,k,m) = 0.0d0 + tnr_racs1(i,j,k,m) = 0.0d0 + tnr_racs2(i,j,k,m) = 0.0d0 + tnr_sacr1(i,j,k,m) = 0.0d0 + tnr_sacr2(i,j,k,m) = 0.0d0 + enddo + enddo + enddo + enddo + + do m = 1, ntb_IN + do k = 1, 45 + do j = 1, ntb_r1 + do i = 1, ntb_r + tpi_qrfz(i,j,k,m) = 0.0d0 + tni_qrfz(i,j,k,m) = 0.0d0 + tpg_qrfz(i,j,k,m) = 0.0d0 + tnr_qrfz(i,j,k,m) = 0.0d0 + enddo + enddo + do j = 1, nbc + do i = 1, ntb_c + tpi_qcfz(i,j,k,m) = 0.0d0 + tni_qcfz(i,j,k,m) = 0.0d0 + enddo + enddo + enddo + enddo + + do j = 1, ntb_i1 + do i = 1, ntb_i + tps_iaus(i,j) = 0.0d0 + tni_iaus(i,j) = 0.0d0 + tpi_ide(i,j) = 0.0d0 + enddo + enddo + + do j = 1, nbc + do i = 1, nbr + t_Efrw(i,j) = 0.0 + enddo + do i = 1, nbs + t_Efsw(i,j) = 0.0 + enddo + enddo + + do k = 1, ntb_r + do j = 1, ntb_r1 + do i = 1, nbr + tnr_rev(i,j,k) = 0.0d0 + enddo + enddo + enddo + + do k = 1, nbc + do j = 1, ntb_c + do i = 1, nbc + tpc_wev(i,j,k) = 0.0d0 + tnc_wev(i,j,k) = 0.0d0 + enddo + enddo + enddo + + do m = 1, ntb_ark + do l = 1, ntb_arr + do k = 1, ntb_art + do j = 1, ntb_arw + do i = 1, ntb_arc + tnccn_act(i,j,k,l,m) = 1.0 + enddo + enddo + enddo + enddo + enddo + + CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ') + WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & + ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g + CALL wrf_debug(150, wrf_err_message) + +!..Read a static file containing CCN activation of aerosols. The +!.. data were created from a parcel model by Feingold & Heymsfield with +!.. further changes by Eidhammer and Kriedenweis. + if (is_aerosol_aware) then + CALL wrf_debug(200, ' calling table_ccnAct routine') + call table_ccnAct + endif + +!..Collision efficiency between rain/snow and cloud water. + CALL wrf_debug(200, ' creating qc collision eff tables') + call table_Efrw + call table_Efsw + +!..Drop evaporation. + CALL wrf_debug(200, ' creating rain evap table') + call table_dropEvap + +!..Initialize various constants for computing radar reflectivity. + xam_r = am_r + xbm_r = bm_r + xmu_r = mu_r + xam_s = am_s + xbm_s = bm_s + xmu_s = mu_s + xam_g = am_g(idx_bg1) + xbm_g = bm_g + xmu_g = mu_g + call radar_init + + if (.not. iiwarm) then + +!..Rain collecting graupel & graupel collecting rain. + CALL wrf_debug(200, ' creating rain collecting graupel table') + + call qr_acr_qg(dimNRHG) + +!..Rain collecting snow & snow collecting rain. + CALL wrf_debug(200, ' creating rain collecting snow table') + call qr_acr_qs + +!..Cloud water and rain freezing (Bigg, 1953). + CALL wrf_debug(200, ' creating freezing of water drops table') + call freezeH2O + +!..Conversion of some ice mass into snow category. + CALL wrf_debug(200, ' creating ice converting to snow table') + call qi_aut_qs + + endif + + CALL wrf_debug(150, ' ... DONE microphysical lookup tables') + + endif + + END SUBROUTINE rcon_init +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..This is a wrapper routine designed to transfer values from 3D to 1D. +!+---+-----------------------------------------------------------------+ + SUBROUTINE mp_rcon_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, & + RAINNC, RAINNCV, & + CLOUDNC, & !RC + SNOWNC, SNOWNCV, & + GRAUPELNC, GRAUPELNCV, SR, & +#if ( WRF_CHEM == 1 ) + wetscav_on, rainprod, evapprod, & +#endif + refl_10cm, diagflag, ke_diag, do_radar_ref, & + re_cloud, re_ice, re_snow, & + has_reqc, has_reqi, has_reqs, & + ids,ide, jds,jde, kds,kde, & ! domain dims + ims,ime, jms,jme, kms,kme, & ! memory dims + its,ite, jts,jte, kts,kte & ! tile dims + ) !RC + + + implicit none + +!..Subroutine arguments + INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + 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, qb, ng + REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, & + nbca2d + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & + re_cloud, re_ice, re_snow + INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs +#if ( WRF_CHEM == 1 ) + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & + rainprod, evapprod +#endif + + !RC + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: CLOUDNC + + !\RC + + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & + pii, p, w, dz + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & + RAINNC, RAINNCV, SR + REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: & + SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & + refl_10cm + REAL, INTENT(IN):: dt_in + INTEGER, INTENT(IN):: itimestep + INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt +#if ( WRF_CHEM == 1 ) + LOGICAL, INTENT(in) :: wetscav_on +#endif + +!..Local variables + REAL, DIMENSION(kts:kte):: & + 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 ) + REAL, DIMENSION(kts:kte):: & + rainprod1d, evapprod1d +#endif + + REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic + REAL:: dt, pptrain, pptsnow, pptgraul, pptice, pptcloud !RC; added pptcloud + REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max + REAL:: nwfa1 + REAL:: ygra1, zans1 + DOUBLE PRECISION:: lamg, lam_exp, lamr, N0_min, N0_exp + INTEGER:: i, j, k, 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 + INTEGER:: i_start, j_start, i_end, j_end + LOGICAL, OPTIONAL, INTENT(IN) :: diagflag + INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref, ke_diag + CHARACTER*256:: mp_debug + + integer :: kediagloc + +!+---+ + + i_start = its + j_start = jts + i_end = MIN(ite, ide-1) + j_end = MIN(jte, jde-1) + +!..For idealized testing by developer. +! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. & +! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then +! i_start = its + 2 +! i_end = ite +! j_start = jts +! j_end = jte +! endif + + dt = dt_in + + qc_max = 0. + qr_max = 0. + qs_max = 0. + qi_max = 0. + qg_max = 0 + ni_max = 0. + nr_max = 0. + imax_qc = 0 + imax_qr = 0 + imax_qi = 0 + imax_qs = 0 + imax_qg = 0 + imax_ni = 0 + imax_nr = 0 + jmax_qc = 0 + jmax_qr = 0 + jmax_qi = 0 + jmax_qs = 0 + jmax_qg = 0 + jmax_ni = 0 + jmax_nr = 0 + kmax_qc = 0 + kmax_qr = 0 + kmax_qi = 0 + kmax_qs = 0 + kmax_qg = 0 + kmax_ni = 0 + kmax_nr = 0 + do i = 1, 256 + mp_debug(i:i) = char(0) + enddo + + if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) & + .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then + write(mp_debug,*) 'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE' + CALL wrf_debug(0, mp_debug) + endif + + j_loop: do j = j_start, j_end + i_loop: do i = i_start, i_end + + pptrain = 0. + pptcloud = 0. !RC + pptsnow = 0. + pptgraul = 0. + pptice = 0. + RAINNCV(i,j) = 0. + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = 0. + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = 0. + ENDIF + SR(i,j) = 0. + + do k = kts, kte + t1d(k) = th(i,k,j)*pii(i,k,j) + p1d(k) = p(i,k,j) + w1d(k) = w(i,k,j) + dz1d(k) = dz(i,k,j) + qv1d(k) = qv(i,k,j) + qc1d(k) = qc(i,k,j) + qi1d(k) = qi(i,k,j) + qr1d(k) = qr(i,k,j) + qs1d(k) = qs(i,k,j) + qg1d(k) = qg(i,k,j) + ni1d(k) = ni(i,k,j) + nr1d(k) = nr(i,k,j) + rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) + enddo + if (is_aerosol_aware) then + do k = kts, kte + nc1d(k) = nc(i,k,j) + nwfa1d(k) = nwfa(i,k,j) + nifa1d(k) = nifa(i,k,j) + if (wif_input_opt .eq. 2) then + nbca1d(k) = nbca(i,k,j) + else + nbca1d(k) = 0.0 + endif + enddo + nwfa1 = nwfa2d(i,j) + else + do k = kts, kte + nc1d(k) = Nt_c/rho(k) + nwfa1d(k) = 11.1E6/rho(k) + nifa1d(k) = naIN1*0.01/rho(k) + nbca1d(k) = 5.55E6/rho(k) + enddo + 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_rcon(aer_init_opt, wif_input_opt, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, & + nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, & + pptrain, pptcloud, pptsnow, pptgraul, pptice, & +#if ( WRF_CHEM == 1 ) + wetscav_on, rainprod1d, evapprod1d, & +#endif + kts, kte, dt, i, j) !RC + + pcp_ra(i,j) = pptrain + pcp_sn(i,j) = pptsnow + pcp_gr(i,j) = pptgraul + pcp_ic(i,j) = pptice + RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice + pptcloud !RC + RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + pptcloud !RC + CLOUDNC(i,j) = CLOUDNC(i,j) + pptcloud !RC + IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN + SNOWNCV(i,j) = pptsnow + pptice + SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice + ENDIF + IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN + GRAUPELNCV(i,j) = pptgraul + GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul + ENDIF + SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) + + +!.. +!..BEGIN AEROSOL EMISSIONS +!.. +!..Reset lowest model level to initial state aerosols (fake sfc source). +!.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol +!.. number tendency (number per kg per second). + if (is_aerosol_aware) then +!..Add anthropogenic emissions +!-GT nwfa1d(kts) = nwfa1 + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt_in + if (wif_input_opt .eq. 2) then + nbca1d(kts) = nbca1d(kts) + nbca2d(i,j)*dt_in + else + nbca1d(kts) = 0.0 + endif +!.. +!..END AEROSOL EMISSIONS +!.. + + do k = kts, kte + nc(i,k,j) = nc1d(k) + nwfa(i,k,j) = nwfa1d(k) + nifa(i,k,j) = nifa1d(k) + if (wif_input_opt .eq. 2) then + nbca(i,k,j) = nbca1d(k) + else + nbca(i,k,j) = 0.0 + 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) + qc(i,k,j) = qc1d(k) + qi(i,k,j) = qi1d(k) + qr(i,k,j) = qr1d(k) + qs(i,k,j) = qs1d(k) + qg(i,k,j) = qg1d(k) + ni(i,k,j) = ni1d(k) + nr(i,k,j) = nr1d(k) + th(i,k,j) = t1d(k)/pii(i,k,j) +#if ( WRF_CHEM == 1 ) + IF ( wetscav_on ) THEN + rainprod(i,k,j) = rainprod1d(k) + evapprod(i,k,j) = evapprod1d(k) + ENDIF +#endif + if (qc1d(k) .gt. qc_max) then + imax_qc = i + jmax_qc = j + kmax_qc = k + qc_max = qc1d(k) + elseif (qc1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (qr1d(k) .gt. qr_max) then + imax_qr = i + jmax_qr = j + kmax_qr = k + qr_max = qr1d(k) + elseif (qr1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (nr1d(k) .gt. nr_max) then + imax_nr = i + jmax_nr = j + kmax_nr = k + nr_max = nr1d(k) + elseif (nr1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (qs1d(k) .gt. qs_max) then + imax_qs = i + jmax_qs = j + kmax_qs = k + qs_max = qs1d(k) + elseif (qs1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (qi1d(k) .gt. qi_max) then + imax_qi = i + jmax_qi = j + kmax_qi = k + qi_max = qi1d(k) + elseif (qi1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (qg1d(k) .gt. qg_max) then + imax_qg = i + jmax_qg = j + kmax_qg = k + qg_max = qg1d(k) + elseif (qg1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (ni1d(k) .gt. ni_max) then + imax_ni = i + jmax_ni = j + kmax_ni = k + ni_max = ni1d(k) + elseif (ni1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + endif + if (qv1d(k) .lt. 0.0) then + write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), & + ' at i,j,k=', i,j,k + CALL wrf_debug(150, mp_debug) + if (k.lt.kte-2 .and. k.gt.kts+1) then + write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) + CALL wrf_debug(150, mp_debug) + qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) + else + qv(i,k,j) = 1.E-7 + endif + endif + enddo + + IF ( PRESENT (diagflag) ) THEN + if (diagflag .and. do_radar_ref == 1) then + + IF ( present(ke_diag) ) THEN + kediagloc = ke_diag + ELSE + kediagloc = kte + ENDIF + + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + ng1d, qb1d, t1d, p1d, dBZ, kts, kte, i, j, kediagloc) + do k = kts, kte + refl_10cm(i,k,j) = MAX(-35., dBZ(k)) + enddo + endif + ENDIF + + IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN + do k = kts, kte + re_qc1d(k) = RE_QC_BG + re_qi1d(k) = RE_QI_BG + re_qs1d(k) = RE_QS_BG + enddo + call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & + re_qc1d, re_qi1d, re_qs1d, kts, kte) + do k = kts, kte + re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc1d(k), 75.E-6)) !RC + re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi1d(k), 125.E-6)) + re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs1d(k), 999.E-6)) + enddo + ENDIF + + enddo i_loop + enddo j_loop + +! DEBUG - GT + write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & + 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & + 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & + 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & + 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & + 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & + 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & + 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' + CALL wrf_debug(150, mp_debug) +! END DEBUG - GT + + do i = 1, 256 + mp_debug(i:i) = char(0) + enddo + + END SUBROUTINE mp_rcon_driver + +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +!.. This subroutine computes the moisture tendencies of water vapor, +!.. cloud droplets, rain, cloud ice (pristine), snow, and graupel. +!.. Previously this code was based on Reisner et al (1998), but few of +!.. those pieces remain. A complete description is now found in +!.. Thompson et al. (2004, 2008). +!+---+-----------------------------------------------------------------+ +! + subroutine mp_rcon (aer_init_opt, wif_input_opt, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & + ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, & + t1d, p1d, w1d, dzq, & + pptrain, pptcloud, pptsnow, pptgraul, pptice, & !RC +#if ( WRF_CHEM == 1 ) + wetscav_on, rainprod, evapprod, & +#endif + kts, kte, dt, ii, jj) + + implicit none + +!..Sub arguments + 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, 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, pptcloud !RC + REAL, INTENT(IN):: dt +#if ( WRF_CHEM == 1 ) + REAL, DIMENSION(kts:kte), INTENT(INOUT):: & + rainprod, evapprod + LOGICAL, INTENT(IN) :: wetscav_on +#endif + +!..Local variables + 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 + + DOUBLE PRECISION, DIMENSION(kts:kte):: pnc_wcd, & + pnc_scw, pnc_gcw + + DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, & + pnd_rcd, pnd_scd, pnd_gcd, pnb_rcb, pnb_scb, pnb_gcb + + DOUBLE PRECISION, DIMENSION(kts:kte):: prr_rcs, & + prr_rcg, prr_sml, prr_gml, & + prr_rci, prv_rev, & + pnr_rcs, pnr_rcg, & + pnr_rci, pnr_sml, pnr_gml, & + pnr_rev, pnr_rcr, pnr_rfz + + DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, & + pni_ihm, pri_wfz, pni_wfz, & + pri_rfz, pni_rfz, pri_ide, & + pni_ide, pri_rci, pni_rci, & + pni_sci, pni_iau, pri_iha, pni_iha + + DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, & + prs_scw, prs_sde, prs_ihm, & + prs_ide + + DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, & + 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, 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 + REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, & + tcond, lvap, ocp, lvt2 + + DOUBLE PRECISION, DIMENSION(kts:kte):: ilamg, N0_r, N0_g + 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, smog + + 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 + DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg + DOUBLE PRECISION:: lami, ilami + REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m + DOUBLE PRECISION:: Dr_star, Dc_star + 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, & + 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:: 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, 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 + REAL:: xslw1, 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 + INTEGER, DIMENSION(kts:kte):: idx_bg + + LOGICAL:: melti, no_micro + LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg + LOGICAL:: debug_flag + CHARACTER*256:: mp_debug + INTEGER:: nu_c + + + !RC; processes and microphysical quantities necessary for this scheme. + !RC; included below in case someone wants to output easier. + DOUBLE PRECISION, DIMENSION(kts:kte):: & + prr_wau,pnr_wau,pnc_wau, & + prr_rcw,pnc_rcw,ilamr,ilamc + !/RC + + !RC; lognormal variables + REAL:: ln_varx,ln_tmp1,ln_tmp2 + REAL, DIMENSION(kts:kte):: ln_dn + REAL, DIMENSION(kts:kte):: ln_sigma + REAL:: sigma_d + !/RC + +!+---+ + + debug_flag = .false. +! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true. + if(debug_flag) then + write(mp_debug, *) 'DEBUG INFO, mp_rcon at (i,j) ', ii, ', ', jj + CALL wrf_debug(550, mp_debug) + endif + + no_micro = .true. + dtsave = dt + odt = 1./dt + odts = 1./dtsave + iexfrq = 1 + +!+---+-----------------------------------------------------------------+ +!.. Source/sink terms. First 2 chars: "pr" represents source/sink of +!.. mass while "pn" represents source/sink of number. Next char is one +!.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for +!.. cloud water, "s" for snow, and "g" for graupel. Next chars +!.. represent processes: "de" for sublimation/deposition, "ev" for +!.. evaporation, "fz" for freezing, "ml" for melting, "au" for +!.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop +!.. secondary ice production, and "c" for collection followed by the +!.. character for the species being collected. ALL of these terms are +!.. positive (except for deposition/sublimation terms which can switch +!.. signs based on super/subsaturation) and are treated as negatives +!.. where necessary in the tendency equations. +!+---+-----------------------------------------------------------------+ + + do k = kts, kte + tten(k) = 0. + qvten(k) = 0. + qcten(k) = 0. + qiten(k) = 0. + qrten(k) = 0. + qsten(k) = 0. + qgten(k) = 0. + ngten(k) = 0. + qbten(k) = 0. + niten(k) = 0. + nrten(k) = 0. + ncten(k) = 0. + nwfaten(k) = 0. + nifaten(k) = 0. + nbcaten(k) = 0. + + prw_vcd(k) = 0. + + pnc_wcd(k) = 0. + pnc_wau(k) = 0. + pnc_rcw(k) = 0. + pnc_scw(k) = 0. + pnc_gcw(k) = 0. + + prv_rev(k) = 0. + prr_wau(k) = 0. + prr_rcw(k) = 0. + prr_rcs(k) = 0. + prr_rcg(k) = 0. + prr_sml(k) = 0. + prr_gml(k) = 0. + prr_rci(k) = 0. + pnr_wau(k) = 0. + pnr_rcs(k) = 0. + pnr_rcg(k) = 0. + pnr_rci(k) = 0. + pnr_sml(k) = 0. + pnr_gml(k) = 0. + pnr_rev(k) = 0. + pnr_rcr(k) = 0. + pnr_rfz(k) = 0. + + pri_inu(k) = 0. + pni_inu(k) = 0. + pri_ihm(k) = 0. + pni_ihm(k) = 0. + pri_wfz(k) = 0. + pni_wfz(k) = 0. + pri_rfz(k) = 0. + pni_rfz(k) = 0. + pri_ide(k) = 0. + pni_ide(k) = 0. + pri_rci(k) = 0. + pni_rci(k) = 0. + pni_sci(k) = 0. + pni_iau(k) = 0. + pri_iha(k) = 0. + pni_iha(k) = 0. + + prs_iau(k) = 0. + prs_sci(k) = 0. + prs_rcs(k) = 0. + prs_scw(k) = 0. + prs_sde(k) = 0. + prs_ihm(k) = 0. + prs_ide(k) = 0. + + prg_scw(k) = 0. + prg_rfz(k) = 0. + prg_gde(k) = 0. + prg_gcw(k) = 0. + prg_rci(k) = 0. + 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. + pna_gca(k) = 0. + + pnd_rcd(k) = 0. + pnd_scd(k) = 0. + pnd_gcd(k) = 0. + + pnb_rcb(k) = 0. + pnb_scb(k) = 0. + pnb_gcb(k) = 0. + enddo +#if ( WRF_CHEM == 1 ) + if ( wetscav_on ) then + do k = kts, kte + rainprod(k) = 0. + evapprod(k) = 0. + enddo + endif +#endif + +!..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments. + do k = kts, kte + smo0(k) = 0. + smo1(k) = 0. + smo2(k) = 0. + smob(k) = 0. + smoc(k) = 0. + smod(k) = 0. + smoe(k) = 0. + smof(k) = 0. + smog(k) = 0. + ns(k) = 0. + mvd_r(k) = 0. + mvd_c(k) = 0. + enddo + +!+---+-----------------------------------------------------------------+ +!..Put column of data into local arrays. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + temp(k) = t1d(k) + qv(k) = MAX(1.E-10, qv1d(k)) + pres(k) = p1d(k) + rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) + if (is_aerosol_aware) then + if (aer_init_opt .lt. 2) then ! Constant or climatology (e.g., GOCART) + nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) + nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) + if (wif_input_opt .eq. 2) then ! Considering BC aerosol + nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k))) + else + nbca(k) = 0.0 + endif + else ! First guess (e.g., GEOS-5) + nwfa(k) = MAX(0.0, nwfa1d(k)*rho(k)) + nifa(k) = MAX(0.0, nifa1d(k)*rho(k)) + if (wif_input_opt .eq. 2) then ! Considering BC aerosol + nbca(k) = MAX(0.0, nbca1d(k)*rho(k)) + else + nbca(k) = 0.0 + endif + endif + else + nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) + nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) + nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k))) + endif + + + ln_sigma(k) = 0.20 !RC + + !RC; This whole if block is written to accomodate lognormal cloud water. + if (qc1d(k) .gt. R1) then + no_micro = .false. + rc(k) = qc1d(k)*rho(k) + nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) + L_qc(k) = .true. + + ! RC; set variable sigma + sigma_d = (rc(k)/nc(k))**(0.3333) + ln_sigma(k) = sigma_d*(-1.185E3) + 0.815 + ln_sigma(k) = MAX(ln_sigma(k) , 0.2) + ln_sigma(k) = MIN(ln_sigma(k) , 0.7) + !/RC + + xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) + if (xDc.lt. D0c) then + xDc = D0c + elseif (xDc.gt. D0r*2.) then + xDc = D0r*2. + endif + + ln_dn(k) = xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)) !RC + nc(k) = MAX( 2. , MIN( DBLE(Nt_c_max) , (1./am_r)*rc(k)*EXP(-4.5*ln_sigma(k)*ln_sigma(k))*(1. / ( xDc*EXP(-1.5*ln_sigma(k)*ln_sigma(k)) ))**3. ) ) !RC + mvd_c(k) = MIN( D0r*2., ln_dn(k)*EXP(3.*ln_sigma(k)*ln_sigma(k)) ) !RC + + if (.NOT. is_aerosol_aware) nc(k) = Nt_c + + else + qc1d(k) = 0.0 + nc1d(k) = 0.0 + rc(k) = R1 + nc(k) = 2. + L_qc(k) = .false. + endif + !\RC + + + if (qi1d(k) .gt. R1) then + no_micro = .false. + ri(k) = qi1d(k)*rho(k) + ni(k) = MAX(R2, ni1d(k)*rho(k)) + if (ni(k).le. R2) then + lami = cie(2)/5.E-6 + ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + endif + L_qi(k) = .true. + lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi + ilami = 1./lami + xDi = (bm_i + mu_i + 1.) * ilami + if (xDi.lt. 5.E-6) then + lami = cie(2)/5.E-6 + ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + elseif (xDi.gt. 300.E-6) then + lami = cie(2)/300.E-6 + ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i + endif + else + qi1d(k) = 0.0 + ni1d(k) = 0.0 + ri(k) = R1 + ni(k) = R2 + L_qi(k) = .false. + endif + + if (qr1d(k) .gt. R1) then + no_micro = .false. + rr(k) = qr1d(k)*rho(k) + nr(k) = MAX(R2, nr1d(k)*rho(k)) + if (nr(k).le. R2) then + mvd_r(k) = 1.0E-3 + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r + endif + L_qr(k) = .true. + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + if (mvd_r(k) .gt. 2.5E-3) then + mvd_r(k) = 2.5E-3 + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r + elseif (mvd_r(k) .lt. D0r*0.75) then + mvd_r(k) = D0r*0.75 + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r + endif + else + qr1d(k) = 0.0 + nr1d(k) = 0.0 + rr(k) = R1 + nr(k) = R2 + L_qr(k) = .false. + endif + if (qs1d(k) .gt. R1) then + no_micro = .false. + rs(k) = qs1d(k)*rho(k) + L_qs(k) = .true. + else + qs1d(k) = 0.0 + rs(k) = R1 + L_qs(k) = .false. + endif + if (qg1d(k) .gt. R1) then + no_micro = .false. + 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(NRHG), 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,NRHG)) + 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(NRHG) + L_qg(k) = .false. + endif + if (.not. is_hail_aware) idx_bg(k) = idx_bg1 + enddo + +!+---+-----------------------------------------------------------------+ +! if (debug_flag) then +! 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, 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 +!+---+-----------------------------------------------------------------+ + +!+---+-----------------------------------------------------------------+ +!..Derive various thermodynamic variables frequently used. +!.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from +!.. 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)) + rhof2(k) = SQRT(rhof(k)) + qvs(k) = rslf(pres(k), temp(k)) + delQvs(k) = MAX(0.0, rslf(pres(k), 273.15)-qv(k)) + if (tempc .le. 0.0) then + 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) + ssatw(k) = satw(k) - 1. + ssati(k) = sati(k) - 1. + if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0 + if (abs(ssati(k)).lt. eps) ssati(k) = 0.0 + if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false. + diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k)) + if (tempc .ge. 0.0) then + visco(k) = (1.718+0.0049*tempc)*1.0E-5 + else + visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5 + endif + ocp(k) = 1./(Cp*(1.+0.887*qv(k))) + 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! +!+---+-----------------------------------------------------------------+ + + if (no_micro) return + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope, and useful moments for snow. +!+---+-----------------------------------------------------------------+ + if (.not. iiwarm) then + do k = kts, kte + if (.not. L_qs(k)) CYCLE + tc0 = MIN(-0.1, temp(k)-273.15) + smob(k) = rs(k)*oams + +!..All other moments based on reference, 2nd moment. If bm_s.ne.2, +!.. then we must compute actual 2nd moment and use as reference. + if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then + smo2(k) = smob(k) + else + loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & + + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & + + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & + + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & + + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & + + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & + + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & + + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & + + sb(10)*bm_s*bm_s*bm_s + smo2(k) = (smob(k)/a_)**(1./b_) + endif + +!..Calculate 0th moment. Represents snow number concentration. + loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0 + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0 + smo0(k) = a_ * smo2(k)**b_ + +!..Calculate 1st moment. Useful for depositional growth and melting. + loga_ = sa(1) + sa(2)*tc0 + sa(3) & + + sa(4)*tc0 + sa(5)*tc0*tc0 & + + sa(6) + sa(7)*tc0*tc0 & + + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 & + + sa(10) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 & + + sb(5)*tc0*tc0 + sb(6) & + + sb(7)*tc0*tc0 + sb(8)*tc0 & + + sb(9)*tc0*tc0*tc0 + sb(10) + smo1(k) = a_ * smo2(k)**b_ + +!..Calculate bm_s+1 (th) moment. Useful for diameter calcs. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + + 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 & + + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) & + + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(13)*cse(13)*cse(13) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) & + + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) & + + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13) + smoe(k) = a_ * smo2(k)**b_ + +!..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) & + + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 & + + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) & + + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(16)*cse(16)*cse(16) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) & + + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) & + + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) & + + 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 + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for graupel. +!+---+-----------------------------------------------------------------+ + + do k = kte, kts, -1 + lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg + ilamg(k) = 1./lamg + N0_g(k) = ng(k)*ogg2*lamg**cge(2,1) + enddo + + endif + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for rain. +!+---+-----------------------------------------------------------------+ + do k = kte, kts, -1 + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + ilamr(k) = 1./lamr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + N0_r(k) = nr(k)*org2*lamr**cre(2) + enddo + +!+---+-----------------------------------------------------------------+ +!..Compute warm-rain process terms (except evap done later). +!+---+-----------------------------------------------------------------+ + + do k = kts, kte + +!..Rain self-collection follows Seifert, 1994 and drop break-up +!.. follows Verlinde and Cotton, 1993. + if (L_qr(k) .and. mvd_r(k).gt. D0r) then + Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) + pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) + endif + + + !RC + mvd_c(k) = D0c + ln_dn(k) = 0.0 + if (L_qc(k)) then + xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) ! dv_bar in paper + ln_dn(k) = xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)) + mvd_c(k) = MIN( D0r*2., ln_dn(k)*EXP(3.*ln_sigma(k)*ln_sigma(k)) ) + endif + !/RC + + !RC; Autoconversion follows Nickerson et al. 1986 + if (rc(k) .gt. 0.01e-3) then + ln_varx = EXP(9*ln_sigma(k)*ln_sigma(k))-1. + ln_tmp1 = 0.067*rc(k)*rc(k)*(1.0E16*((rc(k)/nc(k))**1.3333)*((ln_varx)**0.5) - 2.7 ) + ln_tmp2 = 1.0E4*( (rc(k)*(ln_varx**0.5)/nc(k) )**obmr ) - 1.2 + + prr_wau(k) = MAX( 0.d0,DBLE(ln_tmp1)*DBLE(ln_tmp2)) + prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) + pnr_wau(k) = prr_wau(k) / (am_r*200.*D0r*D0r*D0r) + pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & + / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) + endif + !RC; End new autoconversion + +!..Rain collecting cloud water. In CE, assume Dc<1). Either way, only bother to do sedimentation below +!.. 1st level that contains any sedimenting particles (k=ksed1 on down). +!.. New in v3.0+ is computing separate for rain, ice, snow, and +!.. graupel species thus making code faster with credit to J. Schmidt. +!+---+-----------------------------------------------------------------+ + nstep = 0 + onstep(:) = 1.0 + ksed1(:) = 1 + do k = kte+1, kts, -1 + vtrk(k) = 0. + vtnrk(k) = 0. + vtik(k) = 0. + vtnik(k) = 0. + vtsk(k) = 0. + vtgk(k) = 0. + vtngk(k) = 0. + vtck(k) = 0. + vtnck(k) = 0. + enddo + + if (ANY(L_qr .eqv. .true.)) then + do k = kte, kts, -1 + vtr = 0. + rhof(k) = SQRT(RHO_NOT/rho(k)) + + if (rr(k).gt. R1) then + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) & + *((lamr+fv_r)**(-cre(6))) + vtrk(k) = vtr +! First below is technically correct: +! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) & +! *((lamr+fv_r)**(-cre(5))) +! Test: make number fall faster (but still slower than mass) +! Goal: less prominent size sorting + vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) & + *((lamr+fv_r)**(-cre(7))) + vtnrk(k) = vtr + else + vtrk(k) = vtrk(k+1) + vtnrk(k) = vtnrk(k+1) + endif + + if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then + ksed1(1) = MAX(ksed1(1), k) + delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k))) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(1) .eq. kte) ksed1(1) = kte-1 + if (nstep .gt. 0) onstep(1) = 1./REAL(nstep) + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qc .eqv. .true.)) then + hgt_agl = 0. + do k = kts, kte-1 + if (rc(k) .gt. R2) ksed1(5) = k + hgt_agl = hgt_agl + dzq(k) + if (hgt_agl .gt. 500.0) goto 151 + enddo + 151 continue + + do k = ksed1(5), kts, -1 + vtc = 0. + !RC; Adjustments to cloud fallspeed integrals. Allows drizzle mode. + if (rc(k) .gt. R1 .and. w1d(k) .lt. 0.1E-1) then !RC; allow sedimentation only when w<0.01m/s + xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) + vtc = (av_c*(xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)))**bv_c)*EXP(5.*ln_sigma(k)*ln_sigma(k)) + vtck(k) = MIN(0.1d0,vtc) !RC; These get big, so we need to limit it to 0.1 + vtc = (av_c*(xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)))**bv_c)*EXP(1.3*ln_sigma(k)*ln_sigma(k)) + vtnck(k) = MIN(0.1d0,vtc) + !\RC + endif + enddo + endif + +!+---+-----------------------------------------------------------------+ + + if (.not. iiwarm) then + + if (ANY(L_qi .eqv. .true.)) then + nstep = 0 + do k = kte, kts, -1 + vti = 0. + + if (ri(k).gt. R1) then + lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi + ilami = 1./lami + vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i + vtik(k) = vti +! First below is technically correct: +! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i +! Goal: less prominent size sorting + vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i + vtnik(k) = vti + else + vtik(k) = vtik(k+1) + vtnik(k) = vtnik(k+1) + endif + + if (vtik(k) .gt. 1.E-3) then + ksed1(2) = MAX(ksed1(2), k) + delta_tp = dzq(k)/vtik(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(2) .eq. kte) ksed1(2) = kte-1 + if (nstep .gt. 0) onstep(2) = 1./REAL(nstep) + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qs .eqv. .true.)) then + nstep = 0 + do k = kte, kts, -1 + vts = 0. + + if (rs(k).gt. R1) then + xDs = smoc(k) / smob(k) + Mrat = 1./xDs + ils1 = 1./(Mrat*Lam0 + fv_s) + ils2 = 1./(Mrat*Lam1 + fv_s) + t1_vts = Kap0*csg(4)*ils1**cse(4) + t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10) + ils1 = 1./(Mrat*Lam0) + ils2 = 1./(Mrat*Lam1) + t3_vts = Kap0*csg(1)*ils1**cse(1) + t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7) + vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) + if (prr_sml(k) .gt. 0.0) then + SR = rs(k)/(rs(k)+rr(k)) + vtsk(k) = vts*SR + (1.-SR)*vtrk(k) + else + vtsk(k) = vts*vts_boost(k) + endif + else + vtsk(k) = vtsk(k+1) + endif + + if (vtsk(k) .gt. 1.E-3) then + ksed1(3) = MAX(ksed1(3), k) + delta_tp = dzq(k)/vtsk(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(3) .eq. kte) ksed1(3) = kte-1 + if (nstep .gt. 0) onstep(3) = 1./REAL(nstep) + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qg .eqv. .true.)) then + nstep = 0 + do k = kte, kts, -1 + vtg = 0. + + if (rg(k).gt. R1) then + if (is_hail_aware) then + xrho_g = MAX(rho_g(1),MIN(rg(k)/rho(k)/rb(k),rho_g(NRHG))) + afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff + afall = afall * visco(k)**(1.0-2.0*b_coeff) + else + afall = av_g_old + bfall = bv_g_old + endif + vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bfall + vtgk(k) = vtg +! Goal: less prominent size sorting +! the ELSE section below is technically (mathematically) correct: + if (mu_g .eq. 0) then + vtg = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) * ilamg(k)**bfall + else + vtg = rhof(k)*afall*cgg(8,idx_bg(k))*ogg2 * ilamg(k)**bfall + endif + vtngk(k) = vtg + else + vtgk(k) = vtgk(k+1) + vtngk(k) = vtngk(k+1) + endif + + if (vtgk(k) .gt. 1.E-3) then + ksed1(4) = MAX(ksed1(4), k) + delta_tp = dzq(k)/vtgk(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(4) .eq. kte) ksed1(4) = kte-1 + if (nstep .gt. 0) onstep(4) = 1./REAL(nstep) + endif + endif + +!+---+-----------------------------------------------------------------+ +!..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD, +!.. whereas neglect m(D) term for number concentration. Therefore, +!.. cloud ice has proper differential sedimentation. +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qr .eqv. .true.)) then + nstep = NINT(1./onstep(1)) + do n = 1, nstep + do k = kte, kts, -1 + sed_r(k) = vtrk(k)*rr(k) + sed_n(k) = vtnrk(k)*nr(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho + nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho + rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1)) + nr(k) = MAX(R2, nr(k) - sed_n(k)*odzq*DT*onstep(1)) + do k = ksed1(1), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) & + *odzq*onstep(1)*orho + nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*onstep(1)*orho + rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) & + *odzq*DT*onstep(1)) + nr(k) = MAX(R2, nr(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*DT*onstep(1)) + enddo + + if (rr(kts).gt.R1*1000.) & + pptrain = pptrain + sed_r(kts)*DT*onstep(1) + enddo + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qc .eqv. .true.)) then + do k = kte, kts, -1 + sed_c(k) = vtck(k)*rc(k) + sed_n(k) = vtnck(k)*nc(k) + enddo + do k = ksed1(5), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho + ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho + rc(k) = MAX(R1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*DT) + nc(k) = MAX(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*DT) + enddo + !RC; next two lines we allow cloud water to contribute to surface precip. + !RC; No need for the onstep() variable since cloud water falls so slowly. + if (rc(kts) .gt. R1*1000.) then + pptcloud = pptcloud + sed_c(kts)*DT + endif + + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qi .eqv. .true.)) then + nstep = NINT(1./onstep(2)) + do n = 1, nstep + do k = kte, kts, -1 + sed_i(k) = vtik(k)*ri(k) + sed_n(k) = vtnik(k)*ni(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho + niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho + ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2)) + ni(k) = MAX(R2, ni(k) - sed_n(k)*odzq*DT*onstep(2)) + do k = ksed1(2), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) & + *odzq*onstep(2)*orho + niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*onstep(2)*orho + ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) & + *odzq*DT*onstep(2)) + ni(k) = MAX(R2, ni(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*DT*onstep(2)) + enddo + + if (ri(kts).gt.R1*1000.) & + pptice = pptice + sed_i(kts)*DT*onstep(2) + enddo + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qs .eqv. .true.)) then + nstep = NINT(1./onstep(3)) + do n = 1, nstep + do k = kte, kts, -1 + sed_s(k) = vtsk(k)*rs(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho + rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3)) + do k = ksed1(3), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) & + *odzq*onstep(3)*orho + rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) & + *odzq*DT*onstep(3)) + enddo + + if (rs(kts).gt.R1*1000.) & + pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) + enddo + endif + +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qg .eqv. .true.)) then + nstep = NINT(1./onstep(4)) + do n = 1, nstep + do k = kte, kts, -1 + sed_g(k) = vtgk(k)*rg(k) + sed_n(k) = vtngk(k)*ng(k) + sed_b(k) = vtgk(k)*rb(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho + ngten(k) = ngten(k) - sed_n(k)*odzq*onstep(4)*orho + qbten(k) = qbten(k) - sed_b(k)*odzq*onstep(4) + rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4)) + ng(k) = MAX(R2, ng(k) - sed_n(k)*odzq*DT*onstep(4)) + rb(k) = MAX(R1/rho(k)/rho_g(NRHG), rb(k) - sed_b(k)*odzq*DT*onstep(4)) + do k = ksed1(4), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) & + *odzq*onstep(4)*orho + ngten(k) = ngten(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*onstep(4)*orho + qbten(k) = qbten(k) + (sed_b(k+1)-sed_b(k)) & + *odzq*onstep(4) + rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) & + *odzq*DT*onstep(4)) + ng(k) = MAX(R2, ng(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*DT*onstep(4)) + rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), rb(k) + (sed_b(k+1)-sed_b(k)) & + *odzq*DT*onstep(4)) + enddo + + if (rg(kts).gt.R1*1000.) & + pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) + enddo + endif + +!+---+-----------------------------------------------------------------+ +!.. Instantly melt any cloud ice into cloud water if above 0C and +!.. instantly freeze any cloud water found below HGFR. +!+---+-----------------------------------------------------------------+ + if (.not. iiwarm) then + do k = kts, kte + xri = MAX(0.0, qi1d(k) + qiten(k)*DT) + if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then + qcten(k) = qcten(k) + xri*odt + ncten(k) = ncten(k) + ni1d(k)*odt + qiten(k) = qiten(k) - xri*odt + niten(k) = -ni1d(k)*odt + tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY) + endif + + xrc = MAX(0.0, qc1d(k) + qcten(k)*DT) + if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then + lfus2 = lsub - lvap(k) + xnc = nc1d(k) + ncten(k)*DT + qiten(k) = qiten(k) + xrc*odt + niten(k) = niten(k) + xnc*odt + qcten(k) = qcten(k) - xrc*odt + ncten(k) = ncten(k) - xnc*odt + tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY) + endif + enddo + endif + +!+---+-----------------------------------------------------------------+ +!.. All tendencies computed, apply and pass back final values to parent. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + t1d(k) = t1d(k) + tten(k)*DT + qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) + qc1d(k) = qc1d(k) + qcten(k)*DT + nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max)) + if (is_aerosol_aware) then + if (aer_init_opt .lt. 2) then + nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & + (nwfa1d(k)+nwfaten(k)*DT))) + nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & + (nifa1d(k)+nifaten(k)*DT))) + if (wif_input_opt .eq. 2) then + nbca1d(k) = MAX(5.55E6, MIN(9999.E6, & + (nbca1d(k)+nbcaten(k)*DT))) + else + nbca1d(k) = 0.0 + endif + else + nwfa1d(k) = MAX(0.0, (nwfa1d(k)+nwfaten(k)*DT)) + nifa1d(k) = MAX(0.0, (nifa1d(k)+nifaten(k)*DT)) + if (wif_input_opt .eq. 2) then + nbca1d(k) = MAX(0.0, (nbca1d(k)+nbcaten(k)*DT)) + else + nbca1d(k) = 0.0 + endif + endif + else + nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & + (nwfa1d(k)+nwfaten(k)*DT))) + nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & + (nifa1d(k)+nifaten(k)*DT))) + nbca1d(k) = MAX(5.55E6, MIN(9999.E6, & + (nbca1d(k)+nbcaten(k)*DT))) + endif + + if (qc1d(k) .le. R1) then + qc1d(k) = 0.0 + nc1d(k) = 0.0 + else + !RC; Lognormal size adjustments + xDc = ( ( (1./am_r)*(qc1d(k)/nc1d(k)) )**obmr ) + if (xDc.lt. D0c) then + xDc = D0c + elseif (xDc.gt. D0r*2.) then + xDc = D0r*2. + endif + nc1d(k) = MAX(2., MIN( DBLE(Nt_c_max/rho(k)), (1./am_r)*qc1d(k)*EXP(-4.5*ln_sigma(k)*ln_sigma(k))*(1. / ( xDc*EXP(-1.5*ln_sigma(k)*ln_sigma(k)) ))**3. )) + !RC; end adjustments + endif + + qi1d(k) = qi1d(k) + qiten(k)*DT + ni1d(k) = MAX(R2/rho(k), ni1d(k) + niten(k)*DT) + if (qi1d(k) .le. R1) then + qi1d(k) = 0.0 + ni1d(k) = 0.0 + else + lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi + ilami = 1./lami + xDi = (bm_i + mu_i + 1.) * ilami + if (xDi.lt. 5.E-6) then + lami = cie(2)/5.E-6 + elseif (xDi.gt. 300.E-6) then + lami = cie(2)/300.E-6 + endif + ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & + 999.D3/rho(k)) + endif + qr1d(k) = qr1d(k) + qrten(k)*DT + nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) + if (qr1d(k) .le. R1) then + qr1d(k) = 0.0 + nr1d(k) = 0.0 + else + lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + if (mvd_r(k) .gt. 2.5E-3) then + mvd_r(k) = 2.5E-3 + elseif (mvd_r(k) .lt. D0r*0.75) then + mvd_r(k) = D0r*0.75 + endif + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r + endif + qs1d(k) = qs1d(k) + qsten(k)*DT + if (qs1d(k) .le. R1) qs1d(k) = 0.0 + qg1d(k) = qg1d(k) + qgten(k)*DT + ng1d(k) = MAX(R2/rho(k), ng1d(k) + ngten(k)*DT) + if (qg1d(k) .le. R1) then + qg1d(k) = 0.0 + ng1d(k) = 0.0 + qb1d(k) = 0.0 + else + qb1d(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT) + qb1d(k) = MIN(qg1d(k)/rho_g(1), qb1d(k)) + idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/qb1d(k) *0.01)+1,NRHG)) + if (.not. is_hail_aware) idx_bg(k) = idx_bg1 + lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng1d(k)/qg1d(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 + elseif (mvd_g(k) .lt. D0r) then + mvd_g(k) = D0r + endif + lamg = (3.0 + mu_g + 0.672) / mvd_g(k) + ng1d(k) = cgg(2,1)*ogg3*qg1d(k)*lamg**bm_g / am_g(idx_bg(k)) + endif + + !RC; We set ilamc to 0 if qc 0. + IMPLICIT NONE + REAL, INTENT(IN):: XX + DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 + DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & + COF = (/76.18009172947146D0, -86.50532032941677D0, & + 24.01409824083091D0, -1.231739572450155D0, & + .1208650973866179D-2, -.5395239384953D-5/) + DOUBLE PRECISION:: SER,TMP,X,Y + INTEGER:: J + + X=XX + Y=X + TMP=X+5.5D0 + TMP=(X+0.5D0)*LOG(TMP)-TMP + SER=1.000000000190015D0 + DO 11 J=1,6 + Y=Y+1.D0 + SER=SER+COF(J)/Y +11 CONTINUE + GAMMLN=TMP+LOG(STP*SER/X) + END FUNCTION GAMMLN +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMP(A,X) +! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X) +! --- SEE ABRAMOWITZ AND STEGUN 6.5.1 +! --- USES GCF,GSER + IMPLICIT NONE + REAL, INTENT(IN):: A,X + REAL:: GAMMCF,GAMSER,GLN + GAMMP = 0. + IF((X.LT.0.) .OR. (A.LE.0.)) THEN + PRINT *, 'BAD ARGUMENTS IN GAMMP' + RETURN + ELSEIF(X.LT.A+1.)THEN + CALL GSER(GAMSER,A,X,GLN) + GAMMP=GAMSER + ELSE + CALL GCF(GAMMCF,A,X,GLN) + GAMMP=1.-GAMMCF + ENDIF + END FUNCTION GAMMP +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION WGAMMA(y) + + IMPLICIT NONE + REAL, INTENT(IN):: y + + WGAMMA = EXP(GAMMLN(y)) + + END FUNCTION WGAMMA +!+---+-----------------------------------------------------------------+ +! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS +! A FUNCTION OF TEMPERATURE AND PRESSURE +! + REAL FUNCTION RSLF(P,T) + + IMPLICIT NONE + REAL, INTENT(IN):: P, T + REAL:: ESL,X + REAL, PARAMETER:: C0= .611583699E03 + REAL, PARAMETER:: C1= .444606896E02 + REAL, PARAMETER:: C2= .143177157E01 + REAL, PARAMETER:: C3= .264224321E-1 + REAL, PARAMETER:: C4= .299291081E-3 + REAL, PARAMETER:: C5= .203154182E-5 + REAL, PARAMETER:: C6= .702620698E-8 + REAL, PARAMETER:: C7= .379534310E-11 + REAL, PARAMETER:: C8=-.321582393E-13 + + X=MAX(-80.,T-273.16) + +! ESL=612.2*EXP(17.67*X/(T-29.65)) + ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) + ESL=MIN(ESL, P*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres. + RSLF=.622*ESL/(P-ESL) + +! ALTERNATIVE +! ; Source: Murphy and Koop, Review of the vapour pressure of ice and +! supercooled water for atmospheric applications, Q. J. R. +! Meteorol. Soc (2005), 131, pp. 1539-1565. +! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T +! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22 +! / T - 9.44523 * ALOG(T) + 0.014025 * T)) + + END FUNCTION RSLF +!+---+-----------------------------------------------------------------+ +! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A +! FUNCTION OF TEMPERATURE AND PRESSURE +! + REAL FUNCTION RSIF(P,T) + + IMPLICIT NONE + REAL, INTENT(IN):: P, T + REAL:: ESI,X + REAL, PARAMETER:: C0= .609868993E03 + REAL, PARAMETER:: C1= .499320233E02 + REAL, PARAMETER:: C2= .184672631E01 + REAL, PARAMETER:: C3= .402737184E-1 + REAL, PARAMETER:: C4= .565392987E-3 + REAL, PARAMETER:: C5= .521693933E-5 + REAL, PARAMETER:: C6= .307839583E-7 + REAL, PARAMETER:: C7= .105785160E-9 + REAL, PARAMETER:: C8= .161444444E-12 + + X=MAX(-80.,T-273.16) + ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) + ESI=MIN(ESI, P*0.15) + RSIF=.622*ESI/max(1.e-4,(P-ESI)) + +! ALTERNATIVE +! ; Source: Murphy and Koop, Review of the vapour pressure of ice and +! supercooled water for atmospheric applications, Q. J. R. +! Meteorol. Soc (2005), 131, pp. 1539-1565. +! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T) + + END FUNCTION RSIF + +!+---+-----------------------------------------------------------------+ + real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa) + implicit none + + REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa + +!..Local vars + REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx + REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc + REAL, PARAMETER:: p_c1 = 1000. + REAL, PARAMETER:: p_rho_c = 0.76 + REAL, PARAMETER:: p_alpha = 1.0 + REAL, PARAMETER:: p_gam = 2. + REAL, PARAMETER:: delT = 5. + REAL, PARAMETER:: T0x = -40. + REAL, PARAMETER:: Sw0x = 0.97 + REAL, PARAMETER:: delSi = 0.1 + REAL, PARAMETER:: hdm = 0.15 + REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c + REAL, PARAMETER:: aap = 1. + REAL, PARAMETER:: bbp = 0. + REAL, PARAMETER:: y1p = -35. + REAL, PARAMETER:: y2p = -25. + REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15) + +!+---+ + + xni = 0.0 +! satw = qv/qvs +! sati = qv/qvsi +! siw = qvs/qvsi +! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) & +! + (8.2584e-6*(tempc*tempc*tempc)) +! si0x = 1.+(10.**p_x) +! if (sati.ge.si0x .and. satw.lt.0.985) then +! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm) +! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.) +! dsw = delta_p (satw, Sw0x, 1., 0., 1.) +! fc = dtt*dsi*0.5 +! hx = min(fc+((1.-fc)*dsw), 1.) +! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c +! if (tempc .le. y1p) then +! n_in = ntilde +! elseif (tempc .ge. y2p) then +! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639) +! else +! if (tempc .le. -30.) then +! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c +! else +! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639) +! endif +! ntilde = MIN(ntilde, nmax) +! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax) +! dab = delta_p (tempc, y1p, y2p, aap, bbp) +! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax) +! endif +! mux = hx*p_alpha*n_in*rho +! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.) +! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then + nifa_cc = MAX(0.5, nifa*RHO_NOT0*1.E-6/rho) +! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015] + xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010] + * (nifa_cc**((-0.0264*(tempc))+0.0033)) + xni = xni*rho/RHO_NOT0 * 1000. +! endif + + iceDeMott = MAX(0., xni) + + end FUNCTION iceDeMott + +!+---+-----------------------------------------------------------------+ +!..Newer research since Koop et al (2001) suggests that the freezing +!.. rate should be lower than original paper, so J_rate is reduced +!.. by two orders of magnitude. + + real function iceKoop(temp, qv, qvs, naero, dt) + implicit none + + REAL, INTENT(IN):: temp, qv, qvs, naero, DT + REAL:: mu_diff, a_w_i, delta_aw, log_J_rate, J_rate, prob_h, satw + REAL:: xni + + xni = 0.0 + satw = qv/qvs + mu_diff = 210368.0 + (131.438*temp) - (3.32373E6/temp) & + & - (41729.1*alog(temp)) + a_w_i = exp(mu_diff/(R_uni*temp)) + delta_aw = satw - a_w_i + log_J_rate = -906.7 + (8502.0*delta_aw) & + & - (26924.0*delta_aw*delta_aw) & + & + (29180.0*delta_aw*delta_aw*delta_aw) + log_J_rate = MIN(20.0, log_J_rate) + J_rate = 10.**log_J_rate ! cm-3 s-1 + prob_h = MIN(1.-exp(-J_rate*ar_volume*DT), 1.) + if (prob_h .gt. 0.) then + xni = MIN(prob_h*naero, 1000.E3) + endif + + iceKoop = MAX(0.0, xni) + + end FUNCTION iceKoop + +!+---+-----------------------------------------------------------------+ +!.. Helper routine for Phillips et al (2008) ice nucleation. Trude + + REAL FUNCTION delta_p (yy, y1, y2, aa, bb) + IMPLICIT NONE + + REAL, INTENT(IN):: yy, y1, y2, aa, bb + REAL:: dab, A, B, a0, a1, a2, a3 + + A = 6.*(aa-bb)/((y2-y1)*(y2-y1)*(y2-y1)) + B = aa+(A*y1*y1*y1/6.)-(A*y1*y1*y2*0.5) + a0 = B + a1 = A*y1*y2 + a2 = -A*(y1+y2)*0.5 + a3 = A/3. + + if (yy.le.y1) then + dab = aa + else if (yy.ge.y2) then + dab = bb + else + dab = a0+(a1*yy)+(a2*yy*yy)+(a3*yy*yy*yy) + endif + + if (dab.lt.aa) then + dab = aa + endif + if (dab.gt.bb) then + dab = bb + endif + delta_p = dab + + END FUNCTION delta_p + +!+---+-----------------------------------------------------------------+ +!ctrlL + +!+---+-----------------------------------------------------------------+ +!..Compute _radiation_ effective radii of cloud water, ice, and snow. +!.. These are entirely consistent with microphysics assumptions, not +!.. constant or otherwise ad hoc as is internal to most radiation +!.. schemes. Since only the smallest snowflakes should impact +!.. radiation, compute from first portion of complicated Field number +!.. distribution, not the second part, which is the larger sizes. +!+---+-----------------------------------------------------------------+ + + subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & + & re_qc1d, re_qi1d, re_qs1d, kts, kte) + + IMPLICIT NONE + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte + REAL, DIMENSION(kts:kte), INTENT(IN):: & + & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d + REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d +!..Local variables + INTEGER:: k + REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs + REAL:: smo2, smob, smoc + REAL:: tc0, loga_, a_, b_ + DOUBLE PRECISION:: lamc, lami + LOGICAL:: has_qc, has_qi, has_qs + INTEGER:: inu_c + real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, & + & 504,720,990,1320,1716,2184,2730,3360,4080,4896/) + + DOUBLE PRECISION:: ln_sigma,ln_dn,sigma_d !RC + + has_qc = .false. + has_qi = .false. + has_qs = .false. + + re_qc1d(:) = RE_QC_BG + re_qi1d(:) = RE_QI_BG + re_qs1d(:) = RE_QS_BG + + do k = kts, kte + rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) + rc(k) = MAX(R1, qc1d(k)*rho(k)) + nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) + if (.NOT. is_aerosol_aware) nc(k) = Nt_c + if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. + ri(k) = MAX(R1, qi1d(k)*rho(k)) + ni(k) = MAX(R2, ni1d(k)*rho(k)) + if (ri(k).gt.R1 .and. ni(k).gt.R2) has_qi = .true. + rs(k) = MAX(R1, qs1d(k)*rho(k)) + if (rs(k).gt.R1) has_qs = .true. + enddo + + !RC; New effective radii calculation for cloud. + if (has_qc) then + do k = kts, kte + re_qc1d(k) = 2.49E-6 + if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE + sigma_d = (rc(k)/nc(k))**(0.3333) + ln_sigma = sigma_d*(-1.185E3) + 0.815 + ln_sigma = MAX(ln_sigma , 0.2) + ln_sigma = MIN(ln_sigma , 0.7) + + ln_dn = EXP(-3.0*0.5*ln_sigma**2)*(6.*rc(k)/(PI*nc(k)*1000.))**0.3333 + re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * (ln_dn*EXP(2.5*ln_sigma*ln_sigma) ) ), 75.E-6)) + enddo + endif + !/RC + + if (has_qi) then + do k = kts, kte + if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE + lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi + re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) + enddo + endif + + if (has_qs) then + do k = kts, kte + if (rs(k).le.R1) CYCLE + tc0 = MIN(-0.1, t1d(k)-273.15) + smob = rs(k)*oams + +!..All other moments based on reference, 2nd moment. If bm_s.ne.2, +!.. then we must compute actual 2nd moment and use as reference. + if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then + smo2 = smob + else + loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & + & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & + & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & + & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & + & + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & + & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & + & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & + & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & + & + sb(10)*bm_s*bm_s*bm_s + smo2 = (smob/a_)**(1./b_) + endif +!..Calculate bm_s+1 (th) moment. Useful for diameter calcs. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + & + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) + smoc = a_ * smo2**b_ + re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) + enddo + endif + + end subroutine calc_effectRad + +!+---+-----------------------------------------------------------------+ +!..Compute radar reflectivity assuming 10 cm wavelength radar and using +!.. Rayleigh approximation. Only complication is melted snow/graupel +!.. which we treat as water-coated ice spheres and use Uli Blahak's +!.. library of routines. The meltwater fraction is simply the amount +!.. of frozen species remaining from what initially existed at the +!.. melting level interface. +!+---+-----------------------------------------------------------------+ + + subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + ng1d, qb1d, t1d, p1d, dBZ, kts, kte, ii, jj, ke_diag) + + IMPLICIT NONE + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte, ii, jj, ke_diag + REAL, DIMENSION(kts:kte), INTENT(IN):: & + qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, ng1d, qb1d, t1d, p1d + REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ +! REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ + +!..Local variables + REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof + REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg, ng, rb + + DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g + REAL, DIMENSION(kts:kte):: mvd_r + REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz + REAL:: oM3, M0, Mrat, slam1, slam2, xDs + REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts + REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt + + REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel + INTEGER, DIMENSION(kts:kte):: idx_bg + + DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg + REAL:: a_, b_, loga_, tc0 + DOUBLE PRECISION:: fmelt_s, fmelt_g + + INTEGER:: i, k, k_0, kbot, n, ktop + LOGICAL:: melti + LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg + + DOUBLE PRECISION:: cback, x, eta, f_d + REAL:: xslw1, ygra1, zans1 + INTEGER :: k_0loop + +!+---+ + + do k = kts, kte + dBZ(k) = -35.0 + enddo + +!+---+-----------------------------------------------------------------+ +!..Put column of data into local arrays. +!+---+-----------------------------------------------------------------+ + + do k = kts, kte + temp(k) = t1d(k) + qv(k) = MAX(1.E-10, qv1d(k)) + pres(k) = p1d(k) + rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) + rhof(k) = SQRT(RHO_NOT/rho(k)) + rc(k) = MAX(R1, qc1d(k)*rho(k)) + if (qr1d(k) .gt. R1) then + rr(k) = qr1d(k)*rho(k) + nr(k) = MAX(R2, nr1d(k)*rho(k)) + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + ilamr(k) = 1./lamr + N0_r(k) = nr(k)*org2*lamr**cre(2) + mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k) + L_qr(k) = .true. + else + rr(k) = R1 + nr(k) = R1 + mvd_r(k) = 50.E-6 + L_qr(k) = .false. + endif + if (qs1d(k) .gt. R2) then + rs(k) = qs1d(k)*rho(k) + L_qs(k) = .true. + else + rs(k) = R1 + L_qs(k) = .false. + endif + if (qg1d(k) .gt. R2) then + rg(k) = qg1d(k)*rho(k) + ng(k) = MAX(R2, ng1d(k)*rho(k)) + rb(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k)) + rb(k) = MIN(qg1d(k)/rho_g(1), rb(k)) + idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG)) + if (.not. is_hail_aware) idx_bg(k) = idx_bg1 + L_qg(k) = .true. + else + rg(k) = R1 + ng(k) = R2 + idx_bg(k) = idx_bg1 + L_qg(k) = .false. + endif + enddo + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope, and useful moments for snow. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + smo2(k) = 0. + smob(k) = 0. + smoc(k) = 0. + smoz(k) = 0. + enddo + if ( ( ke_diag > kts .and. ANY(L_qs .eqv. .true.) ) .or. & + (ke_diag == kts .and. L_qs(kts) .eqv. .true. ) ) then + do k = kts, ke_diag ! kte + if (.not. L_qs(k)) CYCLE + tc0 = MIN(-0.1, temp(k)-273.15) + smob(k) = rs(k)*oams + +!..All other moments based on reference, 2nd moment. If bm_s.ne.2, +!.. then we must compute actual 2nd moment and use as reference. + if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then + smo2(k) = smob(k) + else + loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & + & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & + & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & + & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & + & + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & + & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & + & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & + & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & + & + sb(10)*bm_s*bm_s*bm_s + smo2(k) = (smob(k)/a_)**(1./b_) + endif + +!..Calculate bm_s+1 (th) moment. Useful for diameter calcs. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + & + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) + smoc(k) = a_ * smo2(k)**b_ + +!..Calculate bm_s*2 (th) moment. Useful for reflectivity. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) & + & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 & + & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) & + & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 & + & + sa(10)*cse(3)*cse(3)*cse(3) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) & + & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) & + & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) & + & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3) + smoz(k) = a_ * smo2(k)**b_ + enddo + endif + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for graupel. +!+---+-----------------------------------------------------------------+ + + if (ANY(L_qg .eqv. .true.)) then + do k = kte, kts, -1 + lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg + ilamg(k) = 1./lamg + N0_g(k) = ng(k)*ogg2*lamg**cge(2,1) + enddo + else + ilamg(:) = 0. + N0_g(:) = 0. + endif + +!+---+-----------------------------------------------------------------+ +!..Locate K-level of start of melting (k_0 is level above). +!+---+-----------------------------------------------------------------+ + melti = .false. + k_0 = kts + do k = kte-1, kts, -1 + if ( (temp(k).gt.273.15) .and. L_qr(k) & + & .and. (L_qs(k+1).or.L_qg(k+1)) ) then + k_0 = MAX(k+1, k_0) + melti=.true. + goto 195 + endif + enddo + 195 continue + +! Set loop limit for wet ice according to whether the full 3D field is needed or just k=1 + k_0loop = Min(k_0, ke_diag+1) + +!+---+-----------------------------------------------------------------+ +!..Do not do the so-called bright band if using the variable density +!.. graupel-hail category since the density increases during melting and +!.. will account for bright band behavior explicitly. +!+---+-----------------------------------------------------------------+ +! if (is_hail_aware) melti = .false. + +!+---+-----------------------------------------------------------------+ +!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps) +!.. and non-water-coated snow and graupel when below freezing are +!.. simple. Integrations of m(D)*m(D)*N(D)*dD. +!+---+-----------------------------------------------------------------+ + + do k = kts, ke_diag ! kte + ze_rain(k) = 1.e-22 + ze_snow(k) = 1.e-22 + ze_graupel(k) = 1.e-22 + if (L_qr(k)) ze_rain(k) = N0_r(k)*crg(4)*ilamr(k)**cre(4) + if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & + & * (am_s/900.0)*(am_s/900.0)*smoz(k) + if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & + & * (am_g(idx_bg(k))/900.0)*(am_g(idx_bg(k))/900.0) & + & * N0_g(k)*cgg(4,1)*ilamg(k)**cge(4,1) + enddo + +!+---+-----------------------------------------------------------------+ +!..Special case of melting ice (snow/graupel) particles. Assume the +!.. ice is surrounded by the liquid water. Fraction of meltwater is +!.. extremely simple based on amount found above the melting level. +!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting +!.. routines). +!+---+-----------------------------------------------------------------+ + + if (.not. iiwarm .and. melti .and. k_0.ge.2) then + do k = k_0loop-1, kts, -1 + +!..Reflectivity contributed by melting snow + if (L_qs(k) .and. L_qs(k_0) ) then + fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) + eta = 0.d0 + oM3 = 1./smoc(k) + M0 = (smob(k)*oM3) + Mrat = smob(k)*M0*M0*M0 + slam1 = M0 * Lam0 + slam2 = M0 * Lam1 + do n = 1, nrbins + x = am_s * xxDs(n)**bm_s + call rayleigh_soak_wetgraupel (x, DBLE(ocms), DBLE(obms), & + & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, & + & CBACK, mixingrulestring_s, matrixstring_s, & + & inclusionstring_s, hoststring_s, & + & hostmatrixstring_s, hostinclusionstring_s) + f_d = Mrat*(Kap0*DEXP(-slam1*xxDs(n)) & + & + Kap1*(M0*xxDs(n))**mu_s * DEXP(-slam2*xxDs(n))) + eta = eta + f_d * CBACK * simpson(n) * xdts(n) + enddo + ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) + endif + +!..Reflectivity contributed by melting graupel + +! if (L_qg(k) .and. L_qg(k_0) ) then +! fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) +! eta = 0.d0 +! lamg = 1./ilamg(k) +! do n = 1, nrbins +! x = am_g(idx_bg(k)) * xxDg(n)**bm_g +! call rayleigh_soak_wetgraupel (x, DBLE(ocmg(idx_bg(k))), DBLE(obmg), & +! & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, & +! & CBACK, mixingrulestring_g, matrixstring_g, & +! & inclusionstring_g, hoststring_g, & +! & hostmatrixstring_g, hostinclusionstring_g) +! f_d = N0_g(k)*xxDg(n)**mu_g * DEXP(-lamg*xxDg(n)) +! eta = eta + f_d * CBACK * simpson(n) * xdtg(n) +! enddo +! ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta) +! endif + + enddo + endif + + do k = ke_diag, kts, -1 + dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18) + enddo + + +!..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix). +! do k = kte, kts, -1 +! vt_dBZ(k) = 1.E-3 +! if (rs(k).gt.R2) then +! Mrat = smob(k) / smoc(k) +! ils1 = 1./(Mrat*Lam0 + fv_s) +! ils2 = 1./(Mrat*Lam1 + fv_s) +! t1_vts = Kap0*csg(5)*ils1**cse(5) +! t2_vts = Kap1*Mrat**mu_s*csg(11)*ils2**cse(11) +! ils1 = 1./(Mrat*Lam0) +! ils2 = 1./(Mrat*Lam1) +! t3_vts = Kap0*csg(6)*ils1**cse(6) +! t4_vts = Kap1*Mrat**mu_s*csg(12)*ils2**cse(12) +! vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) +! if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then +! vts_dbz_wt = vts_dbz_wt*1.5 +! elseif (temp(k).ge.275.15) then +! vts_dbz_wt = vts_dbz_wt*2.0 +! endif +! else +! vts_dbz_wt = 1.E-3 +! endif + +! if (rr(k).gt.R1) then +! lamr = 1./ilamr(k) +! vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) & +! & / (crg(4)*lamr**(-cre(4))) +! else +! vtr_dbz_wt = 1.E-3 +! endif + +! if (rg(k).gt.R2) then +! lamg = 1./ilamg(k) +! vtg_dbz_wt = rhof(k)*av_g(idx_bg(k))*cgg(5,idx_bg(k))*lamg**(-cge(5,idx_bg(k))) & +! & / (cgg(4,1)*lamg**(-cge(4,1))) +! else +! vtg_dbz_wt = 1.E-3 +! endif + +! vt_dBZ(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) & +! & + vtg_dbz_wt*ze_graupel(k)) & +! & / (ze_rain(k)+ze_snow(k)+ze_graupel(k)) +! enddo + + end subroutine calc_refl10cm +! +!+---+-----------------------------------------------------------------+ +! + real function theta_e(pres_Pa,temp_K,w_non,tlcl_K) +!.. +!.. The following code was based on Bolton (1980) eqn #43 +!.. and claims to have 0.3 K maximum error within -35 < T < 35 C +!.. pres_Pa = Pressure in Pascals +!.. temp_K = Temperature in Kelvin +!.. w_non = mixing ratio (non-dimensional = kg/kg) +!.. tlcl_K = Temperature at Lifting Condensation Level (K) +!.. + IMPLICIT NONE + + real:: pres_Pa, temp_K, w_non, tlcl_K + real:: pp, tt, rr, tlc, power, xx, p1, p2 + +!+---+ + + pp = pres_Pa + tt = temp_K + rr = w_non + 1.e-8 + tlc = tlcl_K + + power=(0.2854*(1.0 - (0.28*rr) )) + xx = tt * (100000.0/pp)**power + + p1 = (3.376/tlc) - 0.00254 + p2 = (rr*1000.0) * (1.0 + 0.81*rr) + + theta_e = xx * exp(p1*p2) + + return + end +! +!+---+-----------------------------------------------------------------+ +! + real function t_lcl(temp_K,tdew_K) +!.. +!.. The following code was based on Bolton (1980) eqn #15 +!.. and claims to have 0.1 K maximum error within -35 < T < 35 C +!.. temp_K = Temperature in Kelvin +!.. tdew_K = Dewpoint T at Lifting Condensation Level (K) +!.. + IMPLICIT NONE + + real:: temp_K, tdew_K + real:: tt, tttd, denom + +!+---+ + + tt = temp_K + tttd= tdew_K + denom= ( 1.0/(tttd-56.0) ) + (log(tt/tttd)/800.) + t_lcl = ( 1.0 / denom ) + 56.0 + return + end +! +!+---+-----------------------------------------------------------------+ +! + real function t_dew(pres_Pa,w_non) +!.. +!.. pres_Pa = Pressure in Pascals +!.. w_non = mixing ratio (non-dimensional = kg/kg) +!.. + IMPLICIT NONE + + real:: pres_Pa, w_non + real:: p, RR, ES, ESLN + +!+---+ + + p = pres_Pa + RR=w_non+1e-8 + ES=P*RR/(.622+RR) + ESLN=LOG(ES) + T_Dew=(35.86*ESLN-4947.2325)/(ESLN-23.6837) + return + end +! +!+---+-----------------------------------------------------------------+ +! + real function theta_wetb(thetae_K) +!.. +!.. Eqn below was gotten from polynomial fit to data in +!.. Smithsonian Meteorological Tables showing Theta-e +!.. and Theta-w +!.. + IMPLICIT NONE + + real:: thetae_K + real:: x, answer + +!+---+ + + real*8 c(0:6), d(0:6) + data c/-1.00922292e-10, -1.47945344e-8, -1.7303757e-6 & + & ,-0.00012709, 1.15849867e-6, -3.518296861e-9 & + & ,3.5741522e-12/ + data d/0.00000000, -3.5223513e-10, -5.7250807e-8 & + & ,-5.83975422e-6, 4.72445163e-8, -1.13402845e-10 & + & ,8.729580402e-14/ + + x=min(475.0,thetae_K) + + if( x .le. 335.5 ) then + answer = c(0)+x*(c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*(c(5)+ & + & x*c(6) ))))) + else + answer = d(0)+x*(d(1)+x*(d(2)+x*(d(3)+x*(d(4)+x*(d(5)+ & + & x*d(6) ))))) + endif + + theta_wetb = answer + 273.15 + + return + end +! +!+---+-----------------------------------------------------------------+ +! + real function compT_fr_The(thelcl_K,pres_Pa) +!.. +!.. pres_Pa = Pressure in Pascals +!.. thelcl = Theta-e at LCL (units in Kelvin) +!.. +!.. Temperature (K) is returned given Theta-e at LCL +!.. and a pressure. This describes a moist-adiabat. +!.. This temperature is the parcel temp at level Pres +!.. along moist adiabat described by theta-e. +!.. + IMPLICIT NONE + + real:: thelcl_K, pres_Pa + real:: guess, epsilon, w1, w2, tenu, tenup, cor, thwlcl_K + integer:: iter + +!+---+ + + guess= (thelcl_K - 0.5 * ( max(thelcl_K-270., 0.))**1.05) & + & * (pres_Pa/100000.0)**.2 + epsilon=0.01 + do iter=1,100 + w1 = rslf(pres_Pa,guess) + w2 = rslf(pres_Pa,guess+1.) + tenu = theta_e(pres_Pa,guess,w1,guess) + tenup = theta_e(pres_Pa,guess+1,w2,guess+1.) + cor = (thelcl_K - tenu) / (tenup - tenu) + guess = guess + cor + if( (cor.lt.epsilon) .and. (-cor.lt.epsilon) ) then + compT_fr_The=guess + return + endif + enddo +! print*, ' convergence not reached ' + thwlcl_K=theta_wetb(thelcl_K) + compT_fr_The = thwlcl_K*((pres_Pa/100000.0)**0.286) + + return + end + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +END MODULE module_mp_rcon +!+---+-----------------------------------------------------------------+ diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 9d419edf7d..67ab230518 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -1004,6 +1004,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & (config_flags%ra_sw_physics .eq. goddardswscheme ) ) .and. & (config_flags%mp_physics .eq. THOMPSON .or. & config_flags%mp_physics .eq. THOMPSONAERO .or. & + config_flags%mp_physics .eq. RCON_MP_SCHEME .or. & (config_flags%mp_physics .eq. NSSL_2MOM .and. config_flags%nssl_2moment_on == 1) .or. & config_flags%mp_physics .eq. WSM3SCHEME .or. & config_flags%mp_physics .eq. WSM5SCHEME .or. & @@ -1037,6 +1038,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & IF ( config_flags%swint_opt .eq. 2 ) THEN IF (( config_flags%mp_physics == THOMPSON .OR. & + config_flags%mp_physics == RCON_MP_SCHEME .OR. & config_flags%mp_physics == THOMPSONAERO .OR. & config_flags%mp_physics == WSM3SCHEME .OR. & config_flags%mp_physics == WSM5SCHEME .OR. & @@ -4379,6 +4381,7 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !------------------------------------------------------------------ + USE module_mp_rcon USE module_mp_wsm3 USE module_mp_wsm5 USE mp_wsm6 @@ -4537,6 +4540,23 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, IMS=ims, IME=ime, JMS=jms, JME=jme, KMS=kms, KME=kme, & ITS=its, ITE=ite, JTS=jts, JTE=jte, KTS=kts, KTE=kte) + CASE (RCON_MP_SCHEME) + IF(start_of_simulation.or.restart.or.config_flags%cycling) & + CALL rcon_init(HGT=z_at_q, & + ORHO=inv_dens, & + NWFA2D=qnwfa2d, NBCA2D=qnbca2d, & + NWFA=scalar(ims,kms,jms,P_QNWFA), & + NIFA=scalar(ims,kms,jms,P_QNIFA), & + NBCA=scalar(ims,kms,jms,P_QNBCA), & + wif_input_opt=config_flags%wif_input_opt, & + FRC_URB2D=frc_urb2d, & + DX=DX, DY=DY, & + is_start=start_of_simulation, & + 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) + + CASE (THOMPSONGH) ! Cycling the WRF forecast with moving nests will cause this initialization to be ! called for each nest move. This is potentially very computationally expensive. diff --git a/run/README.namelist b/run/README.namelist index 4f1d31f71b..59f69e1ec1 100644 --- a/run/README.namelist +++ b/run/README.namelist @@ -255,10 +255,10 @@ Namelist variables specifically for the WPS input for real: = 2: it uses an alternative way (less biased when compared against input data) to compute height in program real and pressure in model. - wif_input_opt = 0 ! = 1: option to process the Water Ice Friendly Aerosol input from metgrid for use with mp_physics=28 - = 2: since V4.4, option to use black carbon aerosol category with mp_physics=28, as well as its radiative effect. Must include + wif_input_opt = 0 ! = 1: option to process the Water Ice Friendly Aerosol input from metgrid for use with mp_physics=28,29 + = 2: since V4.4, option to use black carbon aerosol category with mp_physics=28,29, as well as its radiative effect. Must include file QNWFA_QNIFA_QNBCA_SIGMA_MONTHLY.dat during WPS - num_wif_levels = 30 ! number of levels in the Thompson Water Ice Friendly aerosols (mp_physic=28) + num_wif_levels = 30 ! number of levels in the Thompson Water Ice Friendly aerosols (mp_physic=28,29) p_top_requested = 5000 ! p_top (Pa) to use in the model vert_refine_fact = 1 ! vertical refinement factor for ndown, not used for concurrent vertical grid refinement vert_refine_method (max_dom) = 0 ! vertical refinement method @@ -498,15 +498,17 @@ Namelist variables for controlling the adaptive time step option: This option has two climatological aerosol input options: use_aero_icbc = .F. : use constant values use_aero_icbc = .T. : use climatological aerosol input from WPS - use_rap_aero_icbc = .false. ! Set to .true. to ingest real-time data containing aerosols (new in 4.4) +use_rap_aero_icbc = .false. ! Set to .true. to ingest real-time data containing aerosols (new in 4.4) qna_update = 0 ! set to 1 to update time-varying sfc aerosol emission from climatology or real-time data - with mp_physics = 28. Use with input file ‘wrfqnainp_d0*’ + with mp_physics = 28,29. Use with input file ‘wrfqnainp_d0*’ (must set auxinput17_interval and io_form_auxinput17; new in 4.4) wif_fire_emit = .false. ! set to .true. to include biomass burning organic and black carbon - aerosols with mp_physics = 28 (new in 4.4) + aerosols with mp_physics = 28,29 (new in 4.4) wif_fire_inj = 1 ! (default) vertically distribute biomass burning emissions - in mp_physics = 28 (new in 4.4) + in mp_physics = 28,29 (new in 4.4) + = 29, RCON scheme, Thompson aerosol-aware with liquid-phase and cloud water modifications. + This option supports the thompson aerosol options above. = 30, HUJI (Hebrew University of Jerusalem, Israel) spectral bin microphysics, fast version = 32, HUJI spectral bin microphysics, full version @@ -567,8 +569,8 @@ Namelist variables for controlling the adaptive time step option: 0: do not use; 1: use effective radii (The mp schemes that compute effective radii are 3,4,6,7,8,10,14,16,18,24,26,28,50-53,55) - force_read_thompson = .false. ! whether to read tables for mp_physics = 8,28 - write_thompson_tables = .true. ! whether to read or compute tables for mp_phyiscs = 8,28 + force_read_thompson = .false. ! whether to read tables for mp_physics = 8,28,29 + write_thompson_tables = .true. ! whether to read or compute tables for mp_phyiscs = 8,28,29 write_thompson_mp38table = .false. ! whether to read table (qr_acr_qg_mp38V1.dat) for mp_physics = 38 ra_lw_physics (max_dom) longwave radiation option @@ -993,7 +995,7 @@ Namelist variables for controlling the adaptive time step option: rdmaxalb = .true. ! use snow albedo from geogrid; false means using values from table rdlai2d = .false. ! use LAI from input; false means using values from table if sst_update=1, LAI will also be in wrflowinp file - dust_emis = 0 ! Enable (0=no, 1=yes) surface dust emission scheme to enter mp_physics=28 QNIFA (ice-friendly aerosol variable) + dust_emis = 0 ! Enable (0=no, 1=yes) surface dust emission scheme to enter mp_physics=28,29 QNIFA (ice-friendly aerosol variable) erosion_dim = 3 ! In conjunction with dust_emis=1, this value can only be set equal to 3 (erodibility information) bucket_mm = -1. ! bucket reset value for water accumulations (value in mm, -1.=inactive) bucket_J = -1. ! bucket reset value for energy accumulations (value in J, -1.=inactive) @@ -1829,7 +1831,7 @@ data. To introduce new data sets, mods are required in the Registry and in module_initialize_real.F. There is a space-holder/practice set-up for "GCA". The actual data set for Thompson mp=28 (WIF) that utilizes QNWFA and QNIFA (water and ice friendly aerosols) has -been tested. +been tested. Also compatible with mp=29 (RCON scheme). &domains num_wif_levels = 30 diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 0bc9c2f37e..c655474244 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -2456,12 +2456,12 @@ END FUNCTION bep_bem_ngr_u #endif !----------------------------------------------------------------------- -! grav_settling = 1 must be turned off for mp_physics=28. +! grav_settling = 1 must be turned off for mp_physics=28 or mp_physics=29 !----------------------------------------------------------------------- oops = 0 DO i = 1, model_config_rec % max_dom IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO ) THEN + IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME ) THEN IF ( model_config_rec%grav_settling(i) .NE. FOGSETTLING0 ) THEN model_config_rec%grav_settling(i) = 0 oops = oops + 1 @@ -2474,12 +2474,12 @@ END FUNCTION bep_bem_ngr_u END IF !----------------------------------------------------------------------- -! scalar_pblmix = 1 should be turned on for mp_physics=28. But can be off for MYNN (when bl_mynn_mixscalars = 1) +! scalar_pblmix = 1 should be turned on for mp_physics=28 and mp_physics=29. But can be off for MYNN (when bl_mynn_mixscalars = 1) !----------------------------------------------------------------------- oops = 0 DO i = 1, model_config_rec % max_dom IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO ) THEN + IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME ) THEN IF ( (model_config_rec%use_aero_icbc .OR. model_config_rec%use_rap_aero_icbc) & .AND. model_config_rec%scalar_pblmix(i) .NE. 1 ) THEN model_config_rec%scalar_pblmix(i) = 1 @@ -2511,10 +2511,10 @@ END FUNCTION bep_bem_ngr_u END IF !----------------------------------------------------------------------- -! Set aer_init_opt for Thompson-MP-Aero (mp_physics=28) +! Set aer_init_opt for Thompson-MP-Aero (mp_physics=28) AND rcon mp (mp_physics=29) !----------------------------------------------------------------------- DO i = 1, model_config_rec % max_dom - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO ) THEN + IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME ) THEN IF ( model_config_rec%use_aero_icbc ) THEN model_config_rec%aer_init_opt = 1 ELSE IF ( model_config_rec%use_rap_aero_icbc ) THEN @@ -2524,10 +2524,10 @@ END FUNCTION bep_bem_ngr_u END DO !----------------------------------------------------------------------- -! Check if qna_update=0 when aer_init_opt>1 for Thompson-MP-Aero (mp_physics=28) +! Check if qna_update=0 when aer_init_opt>1 for Thompson-MP-Aero (mp_physics=28) AND rcon mp (mp_physics=29 !----------------------------------------------------------------------- DO i = 1, model_config_rec % max_dom - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO ) THEN + IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME ) THEN IF ( model_config_rec%aer_init_opt .GT. 1 .and. model_config_rec%qna_update .EQ. 0 ) THEN wrf_err_message = '--- ERROR: Time-varying sfc aerosol emissions not set for mp_physics=28 ' CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) @@ -2539,10 +2539,10 @@ END FUNCTION bep_bem_ngr_u END DO !----------------------------------------------------------------------- -! Set aer_fire_emit_opt for Thompson-MP-Aero (mp_physics=28) +! Set aer_fire_emit_opt for Thompson-MP-Aero (mp_physics=28) AND rcon mp (mp_physics=29) !----------------------------------------------------------------------- DO i = 1, model_config_rec % max_dom - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .AND. model_config_rec%wif_fire_emit) THEN + IF ( (model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME) .AND. model_config_rec%wif_fire_emit) THEN IF ( model_config_rec%aer_init_opt .EQ. 2) THEN IF ( model_config_rec%wif_input_opt .EQ. 1 ) THEN model_config_rec%aer_fire_emit_opt = 1 @@ -2564,7 +2564,7 @@ END FUNCTION bep_bem_ngr_u ! is turned on when no PBL scheme is active !----------------------------------------------------------------------- DO i = 1, model_config_rec % max_dom - IF ( model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .AND. model_config_rec%wif_fire_emit) THEN + IF ( (model_config_rec%mp_physics(i) .EQ. THOMPSONAERO .OR. model_config_rec%mp_physics(i) .EQ. RCON_MP_SCHEME) .AND. model_config_rec%wif_fire_emit) THEN IF ( model_config_rec%bl_pbl_physics(i) .EQ. 0 ) THEN wrf_err_message = '--- WARNING: PBL scheme not active but wif_fire_inj=1 for mp_physics=28 ' CALL wrf_debug ( 0, TRIM( wrf_err_message ) )