From de213c9200ec819d019cb890be2ca25e2c1abfec Mon Sep 17 00:00:00 2001 From: Robert Conrick Date: Wed, 22 Jan 2025 18:27:54 -0800 Subject: [PATCH] Add RCON Microphysics (#2144) TYPE: new feature KEYWORDS: microphysics, physics, rainfall, warm rain, cloud water SOURCE: Robert Conrick (U. of Washington); robert.conrick@gmail.com DESCRIPTION OF CHANGES: This is the release of the RCON microphysics scheme, the intent of which is to improve warm rain representation within the Thompson-Eidhammer scheme. LIST OF MODIFIED FILES: modified: Registry/Registry.EM_COMMON modified: Registry/registry.var modified: dyn_em/solve_em.F modified: phys/Makefile modified: phys/module_microphysics_driver.F new file: phys/module_mp_rcon.F modified: phys/module_physics_init.F The code has passed the regression tests. RELEASE NOTE: Release of the RCON Microphysics package into WRF, which improves upon the warm rain representation of 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. 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 accommodate 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. Accompanying published reference: 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. --- Registry/Registry.EM_COMMON | 9 +- Registry/registry.var | 2 + dyn_em/module_initialize_real.F | 36 +- dyn_em/solve_em.F | 5 +- main/depend.common | 10 + main/real_em.F | 16 +- phys/Makefile | 1 + phys/module_microphysics_driver.F | 89 +- phys/module_mp_rcon.F | 6221 +++++++++++++++++++++++++++++ phys/module_physics_init.F | 20 + run/README.namelist | 24 +- share/module_check_a_mundo.F | 22 +- 12 files changed, 6410 insertions(+), 45 deletions(-) create mode 100644 phys/module_mp_rcon.F 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 ) )