From 8b5bfe593e171715aab487eebee249e9874d65e6 Mon Sep 17 00:00:00 2001 From: Timothy Juliano <46979614+twjuliano@users.noreply.github.com> Date: Thu, 20 Jan 2022 13:13:48 -0700 Subject: [PATCH] =?UTF-8?q?Thompson=20AA=20enhancements:=20BC=20aerosol,?= =?UTF-8?q?=20biomass=20burning=20emissions,=20and=20=E2=80=A6=20(#1616)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TYPE: enhancement KEYWORDS: real, Thompson aerosol-aware, microphysics SOURCE: Timothy W. Juliano and Pedro A. Jimenez (NCAR/RAL) DESCRIPTION OF CHANGES: Problem: Several enhancements have been made to the Thompson Aerosol-Aware microphysics scheme (mp_physics=28) related to the addition of: 1. black carbon (BC) aerosol category from either climatology (e.g., GOCART) or first guess aerosol source (e.g., GEOS-5); 2. biomass burning aerosol (organic carbon (OC) and BC) emissions from a first guess aerosol source (e.g., GEOS-5); 3. time-varying surface emissions from a first guess aerosol source (e.g., GEOS-5) Solution: 1. Addition of a new category (3D scalar variable, QNBCA), that represents BC aerosol for the Thompson Aerosol Aware (AA) microphysics scheme. This variable is handled in a similar fashion as QNWFA and QNIFA (water-friendly and ice-friendly aerosols, respectively), including the surface emissions. In terms of the microphysical activity of QNBCA, at present we consider only its removal due to wet scavenging by rain, snow, and graupel. We recognize that this is a limitation of the current implementation; however, more detailed investigation into the water/ice nucleating abilities of BC are needed and left for future work. As a result of this limitation, activation of BC aerosol is made through the &domains namelist option, wif_input_opt: =1 retains the current meaning (water- and ice-friendly aerosols only) and =2 adds the BC aerosol on top of the water- and ice-friendly aerosols. This allows a user to active the Thompson AA scheme as in its original implementation if desired. To account for the radiative contribution of BC aerosols (which are strongly absorbing compared to the QNWFA and QNIFA aerosols that are strongly scattering) to diagnosed AOD at 550nm, we extend the look-up table in subroutine gt_aod. The new look-up table values are computed using Mie code provided by Trude Eidhammer (NCAR/RAL). For BC, we use the following parameter values: * Modal radius = 11.8 nm following Chin et al. (2002) to maintain consistency with the original implementation of the Thompson AA scheme * Geometric standard deviation = 2.0 following Chin et al. (2002) to maintain consistency with the original implementation of the Thompson AA scheme * Real and imaginary indices of refraction = 1.85 and 0.71, respectively, following Bond et al. (2006) * Hygroscopicity = 0.2 following Engelhart et al. (2012) Citations: Chin, M., Ginoux, P., Kinne, S., Torres, O., Holben, B. N., Duncan, B. N., Martin, R. V., Logan, J. A., Higurashi, A., & Nakajima, T. (2002). Tropospheric Aerosol Optical Thickness from the GOCART Model and Comparisons with Satellite and Sun Photometer Measurements, Journal of the Atmospheric Sciences, 59(3), 461-483. Retrieved Dec 22, 2021, from https://journals.ametsoc.org/view/journals/atsc/59/3/1520-0469_2002_059_0461_taotft_2.0.co_2.xml Bond, T. C., G. Habib, and R. W. Bergstrom (2006), Limitations in the enhancement of visible light absorption due to mixing state, J. Geophys. Res., 111, D20211, doi:10.1029/2006JD007315. Engelhart, G. J., Hennigan, C. J., Miracolo, M. A., Robinson, A. L., and Pandis, S. N.: Cloud condensation nuclei activity of fresh primary and aged biomass burning aerosol, Atmos. Chem. Phys., 12, 7285–7293, https://doi.org/10.5194/acp-12-7285-2012, 2012. 2. Addition of OC and BC biomass burning aerosol emissions. These two aerosols are important during periods of active wildfires. Therefore, only when using a first guess aerosol source that has information about biomass burning emissions (e.g., GEOS-5), a user may include these effects through a new &physics namelist option: wif_fire_emit (logical). We note that if wif_fire_emit=.true. and wif_input_opt=1 (i.e., water- and ice-friendly aerosols only), then only OC biomass burning aerosols are emitted, as OC is included in the water-friendly category. However, if wif_fire_emit=.true. and wif_input_opt=2 (i.e., water- and ice-friendly plus BC aerosols), then both OC and BC biomass burning aerosols are emitted. This logic is handled using a surrogate integer variable (aer_fire_emit_opt) that is invisible to the user and thus defined on the Registry as a derived variable. Using the integer variable allows us to properly handle the new emission variables using packages while maintaining simplicity for the user to set a single, logical namelist option. To complement this enhancement, by default we distribute the OC/BC fire aerosols evenly throughout the PBL column (&physics namelist option: wif_fire_inj=1) as a simple plume rise parameterization. This is controllable through the namelist, and we set a warning message in check_a_mundo if the user is not running a PBL scheme with wif_fire_inj=1, as resolved vertical motions may be sufficient such that the parameterization should be turned off. 3. Modifications have been made to the METGRID.TBL file to handle the processing of the new QNBCA aerosol, as well as the surface emissions of anthropogenic QNBCA (QNBCA2D – similar to QNWFA2D and QNIFA2D) and surface emissions of OC and BC biomass burning aerosols (QNOCBB2D and QNBCBB2D, respectively). The modifications to METGRID.TBL may be found in the WPS#190 (https://github.com/wrf-model/WPS/pull/190): We add entries to METGRID.TBL to handle the black carbon aerosol category in addition to biomass burning emissions. Specifically, we: * add monthly climatology entries for black carbon aerosol (B_WIF_*) which generate FLAG_QNBCA_CL for processing in real * add first guess entry for black carbon aerosol (QNBCA) which generates FLAG_QNBCA for processing in real * add first guess entry for anthropogenic emission of black carbon aerosol (QNBCA2D) which generates FLAG_QNBCA2D * add first guess entries for biomass burning emissions of organic carbon (QNOCBB2D) and black carbon (QNBCBB2D) which generate FLAG_QNOCBB2D and FLAG_QNBCBB2D, respectively Note that the new .dat file for monthly GOCART climatology is hosted on Google Drive: https://drive.google.com/file/d/1BYflyu65kP5giRYbTzKo6y4iSnTfb1Fw/view?usp=sharing 4. The ability to have time-varying aerosol emissions has been added. This capability is handled through &physics namelist option qna_update=1, similar to sst_update=1. The I/O for qna_update is done through auxinput17 and the file generated during real is called wrfqnainp_d0* (individual files for each domain, again similar to wrflowinp for sst_update=1). LIST OF MODIFIED FILES: M Registry/Registry.EM_COMMON M Registry/registry.new3d_wif M dyn_em/module_em.F M dyn_em/module_first_rk_step_part1.F M dyn_em/module_initialize_real.F M dyn_em/solve_em.F M dyn_em/start_em.F M main/real_em.F M phys/module_bl_mynn.F M phys/module_microphysics_driver.F M phys/module_mp_thompson.F M phys/module_pbl_driver.F M phys/module_physics_init.F M phys/module_radiation_driver.F M share/input_wrf.F M share/mediation_integrate.F M share/module_check_a_mundo.F M share/module_optional_input.F TESTS CONDUCTED: 1. We have conducted numerous simulations with the new code for a 21-day period in 2016 and a 9-day period in 2020 during active wildfire events. AOD/irradiance quantities from the model output have been compared to observations. For the 2020 wildfire event, our findings are summarized in a manuscript that is currently under review (Juliano, T. W., P. A. Jiménez, B. Kosović, T. Eidhammer, G. Thompson, J. Fast, L. Berg, A. Motley, and A. Polidori, 2021: Smoke from 2020 United States wildfires responsible for substantial solar energy forecast errors, in review at Environmental Research Letters). 2. Jenkins tests are all passing. RELEASE NOTE: A black carbon aerosol category has been added to the Thompson Aerosol-Aware microphysics scheme. Moreover, code enhancements are introduced to allow for time-varying surface aerosol emissions, in addition to consideration of biomass burning organic and black carbon aerosols when using a first guess aerosol source (e.g., GEOS-5). --- Registry/Registry.EM_COMMON | 11 +- Registry/registry.new3d_wif | 24 +++ dyn_em/module_em.F | 33 ++++ dyn_em/module_first_rk_step_part1.F | 3 + dyn_em/module_initialize_real.F | 202 ++++++++++++++++++++++- dyn_em/solve_em.F | 34 +++- dyn_em/start_em.F | 3 +- main/real_em.F | 83 ++++++++-- phys/module_bl_mynn.F | 128 ++++++++++++--- phys/module_microphysics_driver.F | 29 +++- phys/module_mp_thompson.F | 246 +++++++++++++++++++++++++--- phys/module_pbl_driver.F | 27 ++- phys/module_physics_init.F | 19 ++- phys/module_radiation_driver.F | 85 ++++++++-- share/input_wrf.F | 2 +- share/mediation_integrate.F | 76 +++++++++ share/module_check_a_mundo.F | 96 ++++++++++- share/module_optional_input.F | 36 +++- 18 files changed, 1016 insertions(+), 121 deletions(-) diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 9ddd2be5f7..c1e6d5685f 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -488,8 +488,11 @@ state real dfi_qh ikjftb dfi_moist 1 - \ rusdf=(bdy_interp:dt) "DFI_QHAIL" "Hail mixing ratio" "kg kg-1" state real qvold ikj misc 1 - rdu "QVOLD" "Water vapor mixing ratio, old time step" "kg kg-1" state real rimi ikj misc 1 - irh "RIMI" "riming intensity" "fraction" -state real qnwfa2d ij misc 1 - i014rhdu "QNWFA2D" "Surface aerosol number conc emission" "kg-1 s-1" -state real qnifa2d ij misc 1 - i014rhdu "QNIFA2D" "Surface dust number conc emission" "kg-1 s-1" +state real qnwfa2d ij misc 1 - i014{17}rhdu "QNWFA2D" "Surface aerosol number conc emission" "kg-1 s-1" +state real qnifa2d ij misc 1 - i014{17}rhdu "QNIFA2D" "Surface dust number conc emission" "kg-1 s-1" +state real qnbca2d ij misc 1 - i014{17}rhdu "QNBCA2D" "Surface black carbon number conc emission" "kg-1 s-1" +state real qnocbb2d ij misc 1 - i014{17}rhdu "QNOCBB2D" "Surface organic carbon biomass burning number conc emission" "kg-1 s-1" +state real qnbcbb2d ij misc 1 - i014{17}rhdu "QNBCBB2D" "Surface black carbon biomass burning number conc emission" "kg-1 s-1" state real re_cloud ikj misc 1 - r "RE_CLOUD" "Effective radius cloud water" "m" state real re_ice ikj misc 1 - r "RE_ICE" "Effective radius cloud ice" "m" state real re_snow ikj misc 1 - r "RE_SNOW" "Effective radius snow" "m" @@ -2506,6 +2509,7 @@ rconfig real mp_zero_out_thresh namelist,physics 1 1.e-8 rconfig real seaice_threshold namelist,physics 1 100 h "seaice_threshold" "tsk below which which water points are set to sea ice for slab scheme" "K" rconfig logical bmj_rad_feedback namelist,physics max_domains .false. - "if true include radiative effects of bmj clouds" "" rconfig integer sst_update namelist,physics 1 0 h "sst_update" "update sst from wrflowinp file 0=no, 1=yes" "" +rconfig integer qna_update namelist,physics 1 0 h "qna_update" "update aerosol in Thompson-MP-Aero from wrfqnainp file 0=no, 1=yes" "" rconfig integer sst_skin namelist,physics 1 0 h "sst_skin" "calculate sst skin temperature 0=no, 1=yes" "" rconfig integer tmn_update namelist,physics 1 0 h "tmn_update" "update tmn from calculation 0=no, 1=yes" "" @@ -2602,6 +2606,9 @@ rconfig integer tracer_pblmix namelist,physics max_domains 1 rconfig logical use_aero_icbc namelist,physics 1 .false. rh "use_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero" "logical flag" rconfig logical use_rap_aero_icbc namelist,physics 1 .false. r "use_rap_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero from RAP" "logical flag" rconfig integer aer_init_opt derived 1 0 irh "aer_init_opt" "surrogate to handle aerosol IC/BC data in Thompson-MP-Aero: 0=no IC/BC aerosol, 1=climo, 2=first guess" +rconfig logical wif_fire_emit namelist,physics 1 .false. irh "wif_fire_emit" "Activate biomass burning emissionsu in Thompson-MP-Aero" " " +rconfig integer aer_fire_emit_opt derived 1 0 irh "aer_fire_emit_opt" "surrogate to handle aerosol fire emissions in Thompson-MP-Aero: 0=no fire emissions, 1=OC only, 2=OC+BC" " " +rconfig integer wif_fire_inj namelist,physics max_domains 1 irh "wif_fire_inj" "Vertically distribute biomass burning emissions in Thompson-MP-Aero" " " rconfig integer use_mp_re namelist,physics 1 1 h "use_mp_re" "use effective radii computed in some mp schemes in RRTMG" "flag" rconfig logical insert_init_cloud namelist,physics 1 .false. irh "insert_init_cloud" "Insert diagnostic initial cloud using icloud=3 method (Thompson)" "logical flag" diff --git a/Registry/registry.new3d_wif b/Registry/registry.new3d_wif index 1779dfe169..1ecc06d80f 100644 --- a/Registry/registry.new3d_wif +++ b/Registry/registry.new3d_wif @@ -18,6 +18,7 @@ rconfig integer wif_input_opt namelist,domains 1 1 state real qnwfa_gc i{wif}j dyn_em 1 Z i1 "QNWFA" "water-friendly aerosol num concentration" "# kg-1" state real qnifa_gc i{wif}j dyn_em 1 Z i1 "QNIFA" "water-friendly aerosol num concentration" "# kg-1" +state real qnbca_gc i{wif}j dyn_em 1 Z i1 "QNBCA" "black carbon aerosol num concentration" "# kg-1" state real p_wif_gc i{wif}j dyn_em 1 Z i1 "P_WIF" "Pressure for first guess aerosol fields from metgrid" "Pa" state real p_wif_now i{wif}j dyn_em 1 Z - "P_WIF_NOW" "Pressure for Water Ice Friendly Aerosols Now" "Pa" @@ -62,15 +63,38 @@ state real i_wif_oct i{wif}j dyn_em 1 Z i1 state real i_wif_nov i{wif}j dyn_em 1 Z i1 "I_WIF_NOV" "Background Ice Friendly Aerosol option Nov" "# kg-1" state real i_wif_dec i{wif}j dyn_em 1 Z i1 "I_WIF_DEC" "Background Ice Friendly Aerosol option Dec" "# kg-1" +state real b_wif_now i{wif}j dyn_em 1 Z - "B_WIF_NOW" "Background Black Carbon Aerosol option Now" "# kg-1" +state real b_wif_jan i{wif}j dyn_em 1 Z i1 "B_WIF_JAN" "Background Black Carbon Aerosol option Jan" "# kg-1" +state real b_wif_feb i{wif}j dyn_em 1 Z i1 "B_WIF_FEB" "Background Black Carbon Aerosol option Feb" "# kg-1" +state real b_wif_mar i{wif}j dyn_em 1 Z i1 "B_WIF_MAR" "Background Black Carbon Aerosol option Mar" "# kg-1" +state real b_wif_apr i{wif}j dyn_em 1 Z i1 "B_WIF_APR" "Background Black Carbon Aerosol option Apr" "# kg-1" +state real b_wif_may i{wif}j dyn_em 1 Z i1 "B_WIF_MAY" "Background Black Carbon Aerosol option May" "# kg-1" +state real b_wif_jun i{wif}j dyn_em 1 Z i1 "B_WIF_JUN" "Background Black Carbon Aerosol option Jun" "# kg-1" +state real b_wif_jul i{wif}j dyn_em 1 Z i1 "B_WIF_JUL" "Background Black Carbon Aerosol option Jul" "# kg-1" +state real b_wif_aug i{wif}j dyn_em 1 Z i1 "B_WIF_AUG" "Background Black Carbon Aerosol option Aug" "# kg-1" +state real b_wif_sep i{wif}j dyn_em 1 Z i1 "B_WIF_SEP" "Background Black Carbon Aerosol option Sep" "# kg-1" +state real b_wif_oct i{wif}j dyn_em 1 Z i1 "B_WIF_OCT" "Background Black Carbon Aerosol option Oct" "# kg-1" +state real b_wif_nov i{wif}j dyn_em 1 Z i1 "B_WIF_NOV" "Background Black Carbon Aerosol option Nov" "# kg-1" +state real b_wif_dec i{wif}j dyn_em 1 Z i1 "B_WIF_DEC" "Background Black Carbon Aerosol option Dec" "# kg-1" + package use_wif_input wif_input_opt==1 - state:p_wif_gc,p_wif_now,p_wif_jan,p_wif_feb,p_wif_mar,p_wif_apr,p_wif_may,p_wif_jun,p_wif_jul,p_wif_aug,p_wif_sep,p_wif_oct,p_wif_nov,p_wif_dec,qnwfa_gc,w_wif_now,w_wif_jan,w_wif_feb,w_wif_mar,w_wif_apr,w_wif_may,w_wif_jun,w_wif_jul,w_wif_aug,w_wif_sep,w_wif_oct,w_wif_nov,w_wif_dec,qnifa_gc,i_wif_now,i_wif_jan,i_wif_feb,i_wif_mar,i_wif_apr,i_wif_may,i_wif_jun,i_wif_jul,i_wif_aug,i_wif_sep,i_wif_oct,i_wif_nov,i_wif_dec +package use_wif_input_bc wif_input_opt==2 - scalar:qnbca;state:qnbca2d,p_wif_gc,p_wif_now,p_wif_jan,p_wif_feb,p_wif_mar,p_wif_apr,p_wif_may,p_wif_jun,p_wif_jul,p_wif_aug,p_wif_sep,p_wif_oct,p_wif_nov,p_wif_dec,qnwfa_gc,w_wif_now,w_wif_jan,w_wif_feb,w_wif_mar,w_wif_apr,w_wif_may,w_wif_jun,w_wif_jul,w_wif_aug,w_wif_sep,w_wif_oct,w_wif_nov,w_wif_dec,qnifa_gc,i_wif_now,i_wif_jan,i_wif_feb,i_wif_mar,i_wif_apr,i_wif_may,i_wif_jun,i_wif_jul,i_wif_aug,i_wif_sep,i_wif_oct,i_wif_nov,i_wif_dec,b_wif_now,b_wif_jan,b_wif_feb,b_wif_mar,b_wif_apr,b_wif_may,b_wif_jun,b_wif_jul,b_wif_aug,b_wif_sep,b_wif_oct,b_wif_nov,b_wif_dec + +package emit_oc_bb aer_fire_emit_opt==1 - state:qnocbb2d +package emit_ocbc_bb aer_fire_emit_opt==2 - state:qnocbb2d,qnbcbb2d + state real qnwfa ikjftb scalar 1 - \ i0rhusdf=(bdy_interp:dt) "QNWFA" "water-friendly aerosol number con" "# kg(-1)" state real qnifa ikjftb scalar 1 - \ i0rhusdf=(bdy_interp:dt) "QNIFA" "ice-friendly aerosol number con" "# kg(-1)" +state real qnbca ikjftb scalar 1 - \ + i0rhusdf=(bdy_interp:dt) "QNBCA" "black carbon aerosol number con" "# kg(-1)" state real dfi_qnwfa ikjftb dfi_scalar 1 - \ i0rhusdf=(bdy_interp:dt) "DFI_QNWFA" "DFI water-friendly aerosol number con" "# kg(-1)" state real dfi_qnifa ikjftb dfi_scalar 1 - \ i0rhusdf=(bdy_interp:dt) "DFI_QNIFA" "DFI ice-friendly aerosol number con" "# kg(-1)" +state real dfi_qnbca ikjftb dfi_scalar 1 - \ + i0rhusdf=(bdy_interp:dt) "DFI_QNBCA" "DFI black carbon aerosol number con" "# kg(-1)" diff --git a/dyn_em/module_em.F b/dyn_em/module_em.F index c10189e2ee..b71934b641 100644 --- a/dyn_em/module_em.F +++ b/dyn_em/module_em.F @@ -2519,6 +2519,39 @@ SUBROUTINE bound_tke ( tke, tke_upper_bound, & END SUBROUTINE bound_tke +!----------------------------------------------------------------------- + +SUBROUTINE bound_qna ( qna, & + 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(INOUT) :: qna + + INTEGER :: i,k,j + +! +! +! bounds aerosol between zero and current value +! +! + + DO j=jts,min(jte,jde-1) + DO k=kts,kte-1 + DO i=its,min(ite,ide-1) + qna(i,k,j) = max(qna(i,k,j),0.) + ENDDO + ENDDO + ENDDO + + END SUBROUTINE bound_qna + !---------------------------------------------------------------------------------- !cyl: Implement the forward Lagrangian trajectory calculation in WRF !Chiaying Lee RSMAS/UM diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 21e16e6c87..7a5f0e7dfc 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -390,6 +390,8 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & , QNDROP=scalar(ims,kms,jms,P_QNDROP), F_QNDROP=F_QNDROP & & ,QNIFA=scalar(ims,kms,jms,P_QNIFA),F_QNIFA=F_QNIFA & !Trude & ,QNWFA=scalar(ims,kms,jms,P_QNWFA),F_QNWFA=F_QNWFA & !Trude + & ,QNBCA=scalar(ims,kms,jms,P_QNBCA) & + & ,wif_input_opt=config_flags%wif_input_opt & & ,qc_tot=grid%qc_tot, qi_tot=grid%qi_tot & ! Solar diag & ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc & & ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc & @@ -1100,6 +1102,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,QG_CURR=moist(ims,kms,jms,P_QG), F_QG=F_QG & & ,QNIFA_CURR=scalar(ims,kms,jms,P_QNIFA),F_QNIFA=F_QNIFA & & ,QNWFA_CURR=scalar(ims,kms,jms,P_QNWFA),F_QNWFA=F_QNWFA & + & ,QNBCA_CURR=scalar(ims,kms,jms,P_QNBCA),F_QNBCA=F_QNBCA & & ,HOL=HOL, MOL=grid%mol, REGIME=grid%REGIME & !mynn mp@ & ,QKE=grid%qke, Sh3d=grid%sh3d & diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index 6dbf43b809..96629232bf 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -430,6 +430,7 @@ SUBROUTINE init_domain_rk ( grid & hold_ups = ( internal_time_loop .EQ. 1 ) .OR. & ( config_flags%grid_fdda .NE. 0 ) .OR. & ( config_flags%sst_update .EQ. 1 ) .OR. & + ( config_flags%qna_update .EQ. 1 ) .OR. & ( config_flags%all_ic_times ) .OR. & ( config_flags%polar ) @@ -450,6 +451,8 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_message ( a_message ) WRITE ( a_message,* ) ' ( config_flags%sst_update .EQ. 1 ) ', ( config_flags%sst_update .EQ. 1 ) CALL wrf_message ( a_message ) + WRITE ( a_message,* ) ' ( config_flags%qna_update .EQ. 1 ) ', ( config_flags%qna_update .EQ. 1 ) + CALL wrf_message ( a_message ) WRITE ( a_message,* ) ' ( config_flags%all_ic_times ) ', ( config_flags%all_ic_times ) CALL wrf_message ( a_message ) WRITE ( a_message,* ) ' ( config_flags%smooth_cg_topo ) ', ( config_flags%smooth_cg_topo ) @@ -2266,15 +2269,17 @@ SUBROUTINE init_domain_rk ( grid & ! Number of vertical levels: num_wif_levels ! Variable names (assumed to end with _jan, _feb, _mar, ..., _nov, _dec): w_wif, i_wif ! Option to interpolate data: wif_input_opt = 1 (water and ice friendly aerosols) + ! = 2 (water and ice friendly + black carbon aerosols) ! Stored in scalar arrays, tested and assumed to be upside down. ! There are two data fields plus pressure - they are 3d, so no easy way to loop over them. ! QNWFA - Number concentration water-friendly aerosols ! QNIFA - Number concentration ice-friendly aerosols + ! QNBCA - Number concentration black carbon aerosols aer_init_opt = config_flags%aer_init_opt if_thompsonaero_3d: IF (config_flags%mp_physics .EQ. THOMPSONAERO .AND. & - config_flags%wif_input_opt .EQ. 1) THEN + config_flags%wif_input_opt .GT. 0) THEN select_aer_init_opt_3d: select case (aer_init_opt) @@ -2351,6 +2356,7 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_error_fatal ('mp_physics=28 and use_aero_icbc=.true. but aerosol climatology field(s) missing' ) end if do_pres_cl + ! Water-friendly aerosol do_qnwfa_cl: if (flag_qnwfa_cl .EQ. 1) then DO k = 1, config_flags%num_wif_levels DO j = jts, MIN(jte,jde-1) @@ -2406,6 +2412,7 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_error_fatal ('mp_physics=28 but there is not QNWFA aerosol information from climatology' ) end if do_qnwfa_cl + ! Ice-friendly aerosol do_qnifa_cl: if (flag_qnifa_cl .EQ. 1) then DO k = 1, config_flags%num_wif_levels WRITE(a_message,*) ' transferring each K-level ', k, ' to QNIFA, sample Jan data, ', grid %i_wif_jan(its,k,jts) @@ -2471,10 +2478,79 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_error_fatal ('mp_physics=28 but there is not QNIFA aerosol information from climatology' ) end if do_qnifa_cl + ! Black carbon aerosol + if (config_flags%wif_input_opt .EQ. 2) then + do_qnbca_cl: if (flag_qnbca_cl .EQ. 1) then + DO k = 1, config_flags%num_wif_levels + WRITE(a_message,*) ' transferring each K-level ', k, ' to QNBCA, sample Jan data, ', grid %b_wif_jan(its,k,jts) + CALL wrf_debug ( 1 , a_message) + DO j = jts, MIN(jte,jde-1) + DO i = its, MIN(ite,ide-1) + grid%qntemp(i, 1, j) = grid %b_wif_jan(i,k,j) + grid%qntemp(i, 2, j) = grid %b_wif_feb(i,k,j) + grid%qntemp(i, 3, j) = grid %b_wif_mar(i,k,j) + grid%qntemp(i, 4, j) = grid %b_wif_apr(i,k,j) + grid%qntemp(i, 5, j) = grid %b_wif_may(i,k,j) + grid%qntemp(i, 6, j) = grid %b_wif_jun(i,k,j) + grid%qntemp(i, 7, j) = grid %b_wif_jul(i,k,j) + grid%qntemp(i, 8, j) = grid %b_wif_aug(i,k,j) + grid%qntemp(i, 9, j) = grid %b_wif_sep(i,k,j) + grid%qntemp(i,10, j) = grid %b_wif_oct(i,k,j) + grid%qntemp(i,11, j) = grid %b_wif_nov(i,k,j) + grid%qntemp(i,12, j) = grid %b_wif_dec(i,k,j) + END DO + END DO + IF ( k .EQ. 1 ) THEN + WRITE(a_message,*) ' QNBCA (jan and feb) ', grid%qntemp(its,1,jts),grid%qntemp(its,2,jts) + CALL wrf_debug ( 1 , a_message) + END IF + CALL monthly_interp_to_date ( grid%qntemp , current_date , grid%qntemp2 , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + IF ( k .eq. 1 ) THEN + write(a_message,*) ' QNBCA (now) ', grid%qntemp2(its,jts) + CALL wrf_debug ( 1 , a_message) + END IF + IF ( wif_upside_down ) THEN + DO j = jts, MIN(jte,jde-1) + DO i = its, MIN(ite,ide-1) + grid %b_wif_now(i,config_flags%num_wif_levels+1-k,j) = grid%qntemp2(i,j) + END DO + END DO + ELSE IF ( .NOT. wif_upside_down ) THEN + DO j = jts, MIN(jte,jde-1) + DO i = its, MIN(ite,ide-1) + grid %b_wif_now(i, k,j) = grid%qntemp2(i,j) + END DO + END DO + END IF + END DO + + CALL wrf_debug (0 , 'Vertically-interpolating QNBCA climatology from WPS data to fill scalar') + CALL vert_interp ( grid %b_wif_now , grid%p_wif_now , scalar(:,:,:,P_qnbca) , grid%pb , & + grid%hgtmaxw , grid%hgttrop , grid%pmaxw , grid%ptrop , & + grid%pmaxwnn , grid%ptropnn , & + 0 , 0 , & + config_flags%maxw_horiz_pres_diff , config_flags%trop_horiz_pres_diff , & + config_flags%maxw_above_this_level , & + config_flags%num_wif_levels , 'Q' , & + interp_type , linear_interp , extrap_type , & + .false. , use_levels_below_ground , use_surface , & + zap_close_levels , force_sfc_in_vinterp , grid%id , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + else + CALL wrf_error_fatal ('mp_physics=28 and wif_input_opt=2 but there is not QNBCA aerosol information from climatology' ) + end if do_qnbca_cl + end if + case (2) ! First guess aerosol (GEOS-5, etc.) CALL wrf_debug (0 , 'COMMENT: Using first guess aerosols') + ! Water-friendly aerosol do_qnwfa: if (flag_qnwfa .EQ. 1) then if (flag_p_wif .EQ. 1 ) then ! Interpolate according to native pressure field from aerosol forcing model CALL wrf_debug (0 , 'Vertically-interpolating QNWFA first guess from WPS data to fill scalar using native pressure field') @@ -2516,6 +2592,7 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_error_fatal ('mp_physics=28 but there is not QNWFA aerosol information from first guess' ) end if do_qnwfa + ! Ice-friendly aerosol do_qnifa: if (flag_qnifa .EQ. 1) then if (flag_p_wif .EQ. 1) then CALL wrf_debug (0 , 'Vertically-interpolating QNIFA first guess from WPS data to fill scalar using native pressure field') @@ -2557,6 +2634,50 @@ SUBROUTINE init_domain_rk ( grid & CALL wrf_error_fatal ('mp_physics=28 but there is not QNIFA aerosol information from first guess' ) end if do_qnifa + ! Black carbon aerosol + if (config_flags%wif_input_opt .EQ. 2) then + do_qnbca: if (flag_qnbca .EQ. 1) then + if (flag_p_wif .EQ. 1) then + CALL wrf_debug (0 , 'Vertically-interpolating QNBCA first guess from WPS data to fill scalar using native pressure field') + CALL wrf_debug (0 , 'COMMENT: BE SURE num_wif_levels IN NAMELIST EQUALS num_wif_levels IN MET_EM FILES!') + CALL vert_interp ( grid%qnbca_gc , grid%p_wif_gc , scalar(:,:,:,P_QNBCA) , grid%pb , & + grid%hgtmaxw , grid%hgttrop , grid%pmaxw , grid%ptrop , & + grid%pmaxwnn , grid%ptropnn , & + 0 , 0 , & + config_flags%maxw_horiz_pres_diff , config_flags%trop_horiz_pres_diff , & + config_flags%maxw_above_this_level , & + config_flags%num_wif_levels , 'Q' , & + interp_type , linear_interp , extrap_type , & + .false. , use_levels_below_ground , use_surface , & + zap_close_levels , force_sfc_in_vinterp , grid%id , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + else ! Interpolate according to metgrid pressure field + if (config_flags%num_wif_levels .EQ. num_metgrid_levels) then ! Check to make sure that the number of aerosol levels is consistent with the metgrid pressure levels + CALL wrf_debug (0 , 'Vertically-interpolating QNBCA first guess from WPS data to fill scalar using metgrid pressure field') + CALL vert_interp ( grid%qnbca_gc , grid%pd_gc , scalar(:,:,:,P_QNBCA) , grid%pb , & + grid%hgtmaxw , grid%hgttrop , grid%pmaxw , grid%ptrop , & + grid%pmaxwnn , grid%ptropnn , & + 0 , 0 , & + config_flags%maxw_horiz_pres_diff , config_flags%trop_horiz_pres_diff , & + config_flags%maxw_above_this_level , & + config_flags%num_wif_levels , 'Q' , & + interp_type , linear_interp , extrap_type , & + .false. , use_levels_below_ground , use_surface , & + zap_close_levels , force_sfc_in_vinterp , grid%id , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + else + CALL wrf_error_fatal ('num_qnifa_levels not equal to num_metgrid_levels') + end if + end if + else + CALL wrf_error_fatal ('mp_physics=28 and wif_input_opt=2 but there is not QNBCA aerosol information from first guess' ) + end if do_qnbca + end if + case default CALL wrf_debug (0 , 'aer_init_opt = ', aer_init_opt) @@ -4312,7 +4433,7 @@ SUBROUTINE init_domain_rk ( grid & !+---+-----------------------------------------------------------------+ if_thompsonaero_2d: IF (config_flags%mp_physics .EQ. THOMPSONAERO .AND. & - config_flags%wif_input_opt .EQ. 1) THEN + config_flags%wif_input_opt .GT. 0) THEN select_aer_init_opt_2d: select case (aer_init_opt) @@ -4346,6 +4467,17 @@ SUBROUTINE init_domain_rk ( grid & enddo enddo + ! Black carbon aerosol + if (config_flags%wif_input_opt .EQ. 2) then + CALL wrf_debug (0 , 'Calculating anthropogenic surface emissions of QNBCA using climatology') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + z1 = (grid%phb(i,2,j)-grid%phb(i,1,j))/g + grid%qnbca2d(i,j) = grid%b_wif_now(i,1,j) * 0.000098 * (50./z1) * (1. + grid%frc_urb2d(i,j)) + enddo + enddo + end if + case (2) ! First guess aerosol (GEOS-5, etc.) ! Water-friendly aerosol @@ -4386,6 +4518,72 @@ SUBROUTINE init_domain_rk ( grid & enddo end if + ! Black carbon aerosol + if (config_flags%wif_input_opt .EQ. 2) then + if (flag_qnbca2d .EQ. 1) then + CALL wrf_debug (0 , 'Calculating surface emissions of QNBCA using first guess') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + z1 = (grid%phb(i,2,j)-grid%phb(i,1,j))/g + grid%qnbca2d(i,j) = grid%qnbca2d(i,j) * grid%alt(i,1,j) / z1 + enddo + enddo + else + CALL wrf_debug (0 , 'Using first guess aerosol option, but no surface emissions of QNBCA found') + CALL wrf_debug (0 , 'Setting surface emissions of QNBCA to zero') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + grid%qnbca2d(i,j) = 0.0 + enddo + enddo + end if + end if + + ! Biomass burning aerosol + if (config_flags%aer_fire_emit_opt .GT. 0) then + ! Organic carbon first + if (flag_qnocbb2d .EQ. 1) then + CALL wrf_debug (0 , 'Calculating biomass burning surface emissions of organic carbon aerosol using first guess') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + z1 = (grid%phb(i,2,j)-grid%phb(i,1,j))/g + grid%qnocbb2d(i,j) = grid%qnocbb2d(i,j) * grid%alt(i,1,j) / z1 + enddo + enddo + else + CALL wrf_debug (0 , 'Using first guess aerosol option, but no biomass burning surface emissions of organic carbon aerosol found') + CALL wrf_debug (0 , 'Setting biomass burning surface emissions of organic carbon aerosol to zero') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + grid%qnocbb2d(i,j) = 0.0 + enddo + enddo + end if + + ! Black carbon second + if (config_flags%aer_fire_emit_opt .EQ. 2) then + if (flag_qnbcbb2d .EQ. 1) then + CALL wrf_debug (0 , 'Calculating biomass burning surface emissions of black carbon aerosol using first guess') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + z1 = (grid%phb(i,2,j)-grid%phb(i,1,j))/g + grid%qnbcbb2d(i,j) = grid%qnbcbb2d(i,j) * grid%alt(i,1,j) / z1 + enddo + enddo + else + CALL wrf_debug (0 , 'Using first guess aerosol option, but no biomass burning surface emissions of black carbon aerosol found') + CALL wrf_debug (0 , 'Setting biomass burning surface emissions of black carbon aerosol to zero') + do j = jts, min(jde-1,jte) + do i = its, min(ide-1,ite) + grid%qnbcbb2d(i,j) = 0.0 + enddo + enddo + end if + end if + else + CALL wrf_debug (0 , 'Skipping biomass burning surface emissions') + end if + case default CALL wrf_debug (0 , 'aer_init_opt = ', aer_init_opt) diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 5af1767de8..be73d1e304 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -2749,7 +2749,7 @@ SUBROUTINE solve_em ( grid , config_flags & IF( rk_step == 1 ) THEN IF ( config_flags%nested .OR. & ( config_flags%specified .AND. config_flags%have_bcs_scalar ) .OR. & - ( ( is .EQ. P_QNWFA .OR. is .EQ. P_QNIFA) .AND. config_flags%use_aero_icbc ) ) THEN + ( ( is .EQ. P_QNWFA .OR. is .EQ. P_QNIFA .OR. is .EQ. P_QNBCA) .AND. config_flags%aer_init_opt .GT. 0) ) THEN CALL relax_bdy_scalar ( scalar_tend(ims,kms,jms,is), & scalar(ims,kms,jms,is), grid%mut, & @@ -2813,6 +2813,19 @@ SUBROUTINE solve_em ( grid , config_flags & jts=grid%j_start(ij), jte=grid%j_end(ij), & kts=k_start , kte=k_end ) +! bound the aerosol fields (greater than 0) when using first guess aerosol +! as fields may be highly heterogeneous compared to climatology + + IF ( ( ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) .OR. ( is .EQ. P_QNBCA ) ) .AND. & + ( config_flags%aer_init_opt .EQ. 2 ) ) THEN + CALL bound_qna( scalar(ims,kms,jms,is), & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + grid%i_start(ij), grid%i_end(ij), & + grid%j_start(ij), grid%j_end(ij), & + k_start , k_end ) + END IF + IF ( config_flags%specified ) THEN IF (is.EQ.P_QDCN.OR.is.EQ.P_QTCN.OR.is.EQ.P_QNIN) THEN ! for ntu3m CALL flow_dep_bdy_fixed_inflow(scalar(ims,kms,jms,is), & @@ -2835,14 +2848,14 @@ SUBROUTINE solve_em ( grid , config_flags & grid%i_start(ij), grid%i_end(ij), & grid%j_start(ij), grid%j_end(ij), & k_start, k_end ) - ELSE IF ( ( ( ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) ) .AND. & - ( .NOT. config_flags%use_aero_icbc ) ) & + ELSE IF ( ( ( ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) .OR. ( is .EQ. P_QNBCA ) ) .AND. & + ( config_flags%aer_init_opt .EQ. 0 ) ) & .OR. & - ( ( .NOT. ( ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) ) ) .AND. & + ( ( .NOT. ( ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) .OR. ( is .EQ. P_QNBCA ) ) ) .AND. & ( .NOT. config_flags%have_bcs_scalar ) ) ) THEN -! A = ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) -! B = config_flags%use_aero_icbc +! A = ( is .EQ. P_QNWFA ) .OR. ( is .EQ. P_QNIFA ) .OR. ( is .EQ. P_QNBCA ) +! B = config_flags%aer_init_opt .GT. 0 ! C = config_glags%have_bcs_scalar ! Test| A | B | C | ( A AND NOT B ) OR ( NOT A AND NOT C ) @@ -3720,8 +3733,9 @@ SUBROUTINE solve_em ( grid , config_flags & !====================== ! Variables required for CAMMGMP Scheme & ,DLF=grid%dlf,DLF2=grid%dlf2,T_PHY=grid%t_phy,P_HYD=grid%p_hyd & - & ,P8W_HYD=grid%p_hyd_w,TKE_PBL=grid%tke_pbl & - & ,Z_AT_W=grid%z_at_w,QFX=grid%qfx,RLIQ=grid%rliq & + & ,P8W_HYD=grid%p_hyd_w,TKE_PBL=grid%tke_pbl,PBLH=grid%PBLH & + & ,Z_AT_MASS=grid%z,Z_AT_W=grid%z_at_w & + & ,QFX=grid%qfx,RLIQ=grid%rliq & & ,TURBTYPE3D=grid%turbtype3d,SMAW3D=grid%smaw3d & & ,WSEDL3D=grid%wsedl3d,CLDFRA_OLD_MP=grid%cldfra_old_mp & & ,CLDFRA_MP=grid%cldfra_mp,CLDFRA_MP_ALL=grid%cldfra_mp_ALL & @@ -3773,6 +3787,7 @@ SUBROUTINE solve_em ( grid , config_flags & & , QNG_CURR=scalar(ims,kms,jms,P_QNG), F_QNG=F_QNG & & , QNWFA_CURR=scalar(ims,kms,jms,P_QNWFA), F_QNWFA=F_QNWFA & ! for Thompson water-friendly aerosol & , QNIFA_CURR=scalar(ims,kms,jms,P_QNIFA), F_QNIFA=F_QNIFA & ! for Thompson ice-friendly aerosol + & , QNBCA_CURR=scalar(ims,kms,jms,P_QNBCA), F_QNBCA=F_QNBCA & ! for Thompson black carbon aerosol & , QNH_CURR=scalar(ims,kms,jms,P_QNH), F_QNH=F_QNH & ! for milbrandt2mom and nssl_2mom & , QNIC_CURR=scalar(ims,kms,jms,P_QNIC), F_QNIC=F_QNIC & & , QNIP_CURR=scalar(ims,kms,jms,P_QNIP), F_QNIP=F_QNIP & @@ -3854,7 +3869,8 @@ SUBROUTINE solve_em ( grid , config_flags & & , RI_CURR=grid%rimi & & , re_cloud=grid%re_cloud, re_ice=grid%re_ice, re_snow=grid%re_snow & ! G. Thompson & , has_reqc=grid%has_reqc, has_reqi=grid%has_reqi, has_reqs=grid%has_reqs & ! G. Thompson - & , qnwfa2d=grid%qnwfa2d, qnifa2d=grid%qnifa2d & ! G. Thompson + & , qnwfa2d=grid%qnwfa2d, qnifa2d=grid%qnifa2d, qnbca2d=grid%qnbca2d & ! G. Thompson + & , qnocbb2d=grid%qnocbb2d, qnbcbb2d=grid%qnbcbb2d & ! for biomass burning emissions & , diagflag=diag_flag, do_radar_ref=config_flags%do_radar_ref & & , ke_diag=ke_diag & & ,u=grid%u_phy,v=grid%v_phy & diff --git a/dyn_em/start_em.F b/dyn_em/start_em.F index 90a53b9ab9..019208789b 100644 --- a/dyn_em/start_em.F +++ b/dyn_em/start_em.F @@ -1035,7 +1035,8 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & grid%stepbl,grid%stepra,grid%stepcu, & grid%w0avg, grid%rainnc, grid%rainc, grid%raincv, grid%rainncv, & grid%snownc, grid%snowncv, grid%graupelnc, grid%graupelncv, & - z_at_q, grid%alt, grid%qnwfa2d, scalar(ims,kms,jms,1), num_scalar, & ! G. Thompson + z_at_q, grid%alt, grid%qnwfa2d, grid%qnbca2d, & ! G. Thompson + scalar(ims,kms,jms,1), num_scalar, & ! G. Thompson grid%re_cloud, grid%re_ice, grid%re_snow, & ! G. Thompson grid%has_reqc, grid%has_reqi, grid%has_reqs, & ! G. Thompson grid%phb,grid%ph_2,grid%p,grid%pb, & ! for ntu3m diff --git a/main/real_em.F b/main/real_em.F index 704e372bb5..796a0ac219 100644 --- a/main/real_em.F +++ b/main/real_em.F @@ -397,7 +397,8 @@ END SUBROUTINE start_domain internal_time_loop = loop IF ( ( grid%id .GT. 1 ) .AND. ( loop .GT. 1 ) .AND. & ( model_config_rec%grid_fdda(grid%id) .EQ. 0 ) .AND. & - ( model_config_rec%sst_update .EQ. 0 ) ) EXIT + ( model_config_rec%sst_update .EQ. 0 ) .AND. & + ( model_config_rec%qna_update .EQ. 0 ) ) EXIT print *,' ' print *,'-----------------------------------------------------------------------------' @@ -696,8 +697,8 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren INTEGER :: ijds , ijde , spec_bdy_width INTEGER :: i , j , k , idts - INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, grid_fdda - INTEGER , SAVE :: id, id2, id4 + INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, qna_update, grid_fdda + INTEGER , SAVE :: id, id2, id4, id17 CHARACTER (LEN=256) :: inpname , bdyname CHARACTER(LEN= 4) :: loop_char CHARACTER (LEN=256) :: message @@ -709,7 +710,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren REAL , DIMENSION(:,:,:) , ALLOCATABLE , SAVE :: mbdy2dtemp1 REAL , DIMENSION(:,:,:) , ALLOCATABLE , SAVE :: ubdy3dtemp2 , vbdy3dtemp2 , tbdy3dtemp2 , pbdy3dtemp2 , qbdy3dtemp2 REAL , DIMENSION(:,:,:) , ALLOCATABLE , SAVE :: mbdy2dtemp2 - REAL , DIMENSION(:,:,:) , ALLOCATABLE , SAVE :: qn1bdy3dtemp1, qn1bdy3dtemp2, qn2bdy3dtemp1, qn2bdy3dtemp2 + REAL , DIMENSION(:,:,:) , ALLOCATABLE , SAVE :: qn1bdy3dtemp1, qn1bdy3dtemp2, qn2bdy3dtemp1, qn2bdy3dtemp2, qn3bdy3dtemp1, qn3bdy3dtemp2 real::t1,t2 INTEGER :: open_status @@ -747,6 +748,7 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren spec_bdy_width = model_config_rec%spec_bdy_width interval_seconds = model_config_rec%interval_seconds sst_update = model_config_rec%sst_update + qna_update = model_config_rec%qna_update grid_fdda = model_config_rec%grid_fdda(grid%id) ! Domain check. We cannot decompose the domain into pieces that are @@ -800,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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN IF ( ALLOCATED ( qn1bdy3dtemp1 ) ) DEALLOCATE ( qn1bdy3dtemp1 ) IF ( ALLOCATED ( qn2bdy3dtemp1 ) ) DEALLOCATE ( qn2bdy3dtemp1 ) IF ( ALLOCATED ( qn1bdy3dtemp2 ) ) DEALLOCATE ( qn1bdy3dtemp2 ) @@ -809,6 +811,12 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ALLOCATE ( qn2bdy3dtemp1(ims:ime,kms:kme,jms:jme) ) ALLOCATE ( qn1bdy3dtemp2(ims:ime,kms:kme,jms:jme) ) ALLOCATE ( qn2bdy3dtemp2(ims:ime,kms:kme,jms:jme) ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + IF ( ALLOCATED ( qn3bdy3dtemp1 ) ) DEALLOCATE ( qn3bdy3dtemp1 ) + IF ( ALLOCATED ( qn3bdy3dtemp2 ) ) DEALLOCATE ( qn3bdy3dtemp2 ) + ALLOCATE ( qn3bdy3dtemp1(ims:ime,kms:kme,jms:jme) ) + ALLOCATE ( qn3bdy3dtemp2(ims:ime,kms:kme,jms:jme) ) + END IF END IF END IF @@ -834,6 +842,15 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren END IF CALL output_auxinput4 ( id4, grid , config_flags , ierr ) END IF + + IF(qna_update .EQ. 1)THEN + CALL construct_filename1( inpname , 'wrfqnainp' , grid%id , 2 ) + CALL open_w_dataset ( id17, TRIM(inpname) , grid , config_flags , output_auxinput17 , "DATASET=AUXINPUT17", ierr ) + IF ( ierr .NE. 0 ) THEN + CALL wrf_error_fatal( 'real: error opening wrfqnainp for writing' ) + END IF + CALL output_auxinput17 ( id17, grid , config_flags , ierr ) + END IF END IF IF ( ( time_loop_max .EQ. 1 ) .OR. ( config_flags%polar ) ) THEN @@ -868,13 +885,18 @@ 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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .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 ) CALL couple ( grid%mu_2 , grid%mub , qn2bdy3dtemp1 , grid%scalar(:,:,:,P_QNIFA) , '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 ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + CALL couple ( grid%mu_2 , grid%mub , qn3bdy3dtemp1 , grid%scalar(:,:,:,P_QNBCA) , '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 ) + END IF END IF END IF @@ -943,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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .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 , & @@ -956,6 +978,14 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & ips , ipe , jps , jpe , kps , kpe ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + CALL stuff_bdy ( qn3bdy3dtemp1 , grid%scalar_bxs(:,:,:,P_QNBCA), grid%scalar_bxe(:,:,:,P_QNBCA), & + grid%scalar_bys(:,:,:,P_QNBCA), grid%scalar_bye(:,:,:,P_QNBCA), & + 'T' , spec_bdy_width , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + ips , ipe , jps , jpe , kps , kpe ) + END IF END IF END IF @@ -966,6 +996,10 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren CALL output_auxinput4 ( id4, grid , config_flags , ierr ) END IF + IF(qna_update .EQ. 1)THEN + CALL output_auxinput17 ( id17, grid , config_flags , ierr ) + END IF + ! Open the boundary and the fdda file. IF ( loop .eq. 2 ) THEN @@ -1037,13 +1071,18 @@ 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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .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 ) CALL couple ( grid%mu_2 , grid%mub , qn2bdy3dtemp2 , grid%scalar(:,:,:,P_QNIFA) , '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 ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + CALL couple ( grid%mu_2 , grid%mub , qn3bdy3dtemp2 , grid%scalar(:,:,:,P_QNBCA) , '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 ) + END IF END IF END IF @@ -1129,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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .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), & @@ -1146,6 +1185,16 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & ips , ipe , jps , jpe , kps , kpe ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + CALL stuff_bdytend ( qn3bdy3dtemp2 , qn3bdy3dtemp1 , REAL(interval_seconds) , & + grid%scalar_btxs(:,:,:,P_QNBCA), grid%scalar_btxe(:,:,:,P_QNBCA), & + grid%scalar_btys(:,:,:,P_QNBCA), grid%scalar_btye(:,:,:,P_QNBCA), & + 'T' , & + spec_bdy_width , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + ips , ipe , jps , jpe , kps , kpe ) + END IF END IF END IF @@ -1257,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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .AND. config_flags%aer_init_opt .gt. 0) THEN DO j = jps , jpe DO k = kps , kpe DO i = ips , ipe @@ -1334,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%use_aero_icbc) THEN + IF (config_flags%mp_physics.eq.THOMPSONAERO .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 , & @@ -1347,6 +1396,14 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & ips , ipe , jps , jpe , kps , kpe ) + IF ( ( config_flags%wif_input_opt .EQ. 2 ) .AND. ( f_qnbca ) ) THEN + CALL stuff_bdy ( qn3bdy3dtemp1 , grid%scalar_bxs(:,:,:,P_QNBCA), grid%scalar_bxe(:,:,:,P_QNBCA), & + grid%scalar_bys(:,:,:,P_QNBCA), grid%scalar_bye(:,:,:,P_QNBCA), & + 'T' , spec_bdy_width , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + ips , ipe , jps , jpe , kps , kpe ) + END IF END IF END IF @@ -1373,6 +1430,10 @@ SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max , curren CALL close_dataset ( id4 , config_flags , "DATASET=AUXINPUT4" ) END IF + IF(qna_update .EQ. 1)THEN + CALL close_dataset ( id17 , config_flags , "DATASET=AUXINPUT17" ) + END IF + END IF END IF diff --git a/phys/module_bl_mynn.F b/phys/module_bl_mynn.F index aef56bf5ce..6f7fc3fa68 100644 --- a/phys/module_bl_mynn.F +++ b/phys/module_bl_mynn.F @@ -2766,25 +2766,25 @@ SUBROUTINE mynn_tendencies(kts,kte, & &u,v,th,tk,qv,qc,qi,qnc,qni, & &p,exner, & &thl,sqv,sqc,sqi,sqw, & - &qnwfa,qnifa, & + &qnwfa,qnifa,qnbca, & &ust,flt,flq,flqv,flqc,wspd,qcg, & &uoce,voce, & &tsq,qsq,cov, & &tcd,qcd, & &dfm,dfh,dfq, & &Du,Dv,Dth,Dqv,Dqc,Dqi,Dqnc,Dqni, & - &Dqnwfa,Dqnifa, & + &Dqnwfa,Dqnifa,Dqnbca, & &vdfg1,diss_heat, & &s_aw,s_awthl,s_awqt,s_awqv,s_awqc, & &s_awu,s_awv, & &s_awqnc,s_awqni, & - &s_awqnwfa,s_awqnifa, & + &s_awqnwfa,s_awqnifa,s_awqnbca, & &sub_thl,sub_sqv, & &sub_u,sub_v, & &det_thl,det_sqv,det_sqc, & &det_u,det_v, & &FLAG_QC,FLAG_QI,FLAG_QNC,FLAG_QNI, & - &FLAG_QNWFA,FLAG_QNIFA, & + &FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, & &cldfra_bl1d, & &bl_mynn_cloudmix, & &bl_mynn_mixqt, & @@ -2805,7 +2805,7 @@ SUBROUTINE mynn_tendencies(kts,kte, & bl_mynn_edmf,bl_mynn_edmf_mom, & bl_mynn_mixscalars LOGICAL, INTENT(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,& - FLAG_QNWFA,FLAG_QNIFA + FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA !! grav_settling = 1 or 2 for gravitational settling of droplets !! grav_settling = 0 otherwise @@ -2817,16 +2817,17 @@ SUBROUTINE mynn_tendencies(kts,kte, & ! mass-flux plumes REAL, DIMENSION(kts:kte+1), INTENT(in) :: s_aw,s_awthl,s_awqt,& - &s_awqnc,s_awqni,s_awqv,s_awqc,s_awu,s_awv,s_awqnwfa,s_awqnifa + &s_awqnc,s_awqni,s_awqv,s_awqc,s_awu,s_awv,s_awqnwfa,s_awqnifa,& + s_awqnbca ! tendencies from mass-flux environmental subsidence and detrainment REAL, DIMENSION(kts:kte), INTENT(in) :: sub_thl,sub_sqv, & &sub_u,sub_v,det_thl,det_sqv,det_sqc,det_u,det_v REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,tk,qv,qc,qi,qni,qnc,& &rho,p,exner,dfq,dz,tsq,qsq,cov,tcd,qcd,cldfra_bl1d,diss_heat REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc,sqi,& - &qnwfa,qnifa,dfm,dfh + &qnwfa,qnifa,qnbca,dfm,dfh REAL, DIMENSION(kts:kte), INTENT(inout) :: du,dv,dth,dqv,dqc,dqi,& - &dqni,dqnc,dqnwfa,dqnifa + &dqni,dqnc,dqnwfa,dqnifa,dqnbca REAL, INTENT(IN) :: delt,ust,flt,flq,flqv,flqc,wspd,uoce,voce,qcg ! REAL, INTENT(IN) :: delt,ust,flt,flq,qcg,& @@ -2836,7 +2837,7 @@ SUBROUTINE mynn_tendencies(kts,kte, & REAL, DIMENSION(kts:kte) :: dtz,vt,vq,dfhc,dfmc !Kh for clouds (Pr < 2) REAL, DIMENSION(kts:kte) :: sqv2,sqc2,sqi2,sqw2,qni2,qnc2, & !AFTER MIXING - qnwfa2,qnifa2 + qnwfa2,qnifa2,qnbca2 REAL, DIMENSION(kts:kte) :: zfac,plumeKh REAL, DIMENSION(kts:kte) :: a,b,c,d,x REAL, DIMENSION(kts:kte+1) :: rhoz, & !rho on model interface @@ -3477,6 +3478,48 @@ SUBROUTINE mynn_tendencies(kts,kte, & qnifa2=qnifa ENDIF +!============================================ +! Black carbon aerosols ( qnbca ). +!============================================ +IF (bl_mynn_cloudmix > 0 .AND. FLAG_QNBCA .AND. & + bl_mynn_mixscalars > 0) THEN + + k=kts + + a(k)= -dtz(k)*khdz(k)/rho(k) + b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))/rho(k) - & + & 0.5*dtz(k)*s_aw(k+1)*nonloc + c(k)= -dtz(k)*khdz(k+1)/rho(k) - 0.5*dtz(k)*s_aw(k+1)*nonloc + d(k)=qnbca(k) - dtz(k)*s_awqnbca(k+1)*nonloc + + DO k=kts+1,kte-1 + a(k)= -dtz(k)*khdz(k)/rho(k) + 0.5*dtz(k)*s_aw(k)*nonloc + b(k)=1.+dtz(k)*(khdz(k) + khdz(k+1))/rho(k) + & + & 0.5*dtz(k)*(s_aw(k)-s_aw(k+1))*nonloc + c(k)= -dtz(k)*khdz(k+1)/rho(k) - 0.5*dtz(k)*s_aw(k+1)*nonloc + d(k)=qnbca(k) + dtz(k)*(s_awqnbca(k)-s_awqnbca(k+1))*nonloc + ENDDO + +! prescribed value + a(kte)=0. + b(kte)=1. + c(kte)=0. + d(kte)=qnbca(kte) + +! CALL tridiag(kte,a,b,c,d) +! CALL tridiag2(kte,a,b,c,d,x) + CALL tridiag3(kte,a,b,c,d,x) + + DO k=kts,kte + !qnbca2(k)=d(k-kts+1) + qnbca2(k)=x(k) + ENDDO + +ELSE + !If not mixing aerosols, set "updated" array equal to original array + qnbca2=qnbca +ENDIF + !!============================================ !! Compute tendencies and convert to mixing ratios for WRF. @@ -3634,11 +3677,18 @@ SUBROUTINE mynn_tendencies(kts,kte, & ! Ice-friendly aerosols !===================== Dqnifa(k)=(qnifa2(k) - qnifa(k))/delt + IF (FLAG_QNBCA) THEN + !===================== + ! Black carbon aerosols + !===================== + Dqnbca(k)=(qnbca2(k) - qnbca(k))/delt + ENDIF ENDDO ELSE DO k=kts,kte Dqnwfa(k)=0. Dqnifa(k)=0. + Dqnbca(k)=0. ENDDO ENDIF @@ -3877,7 +3927,7 @@ SUBROUTINE mynn_bl_driver( & &grav_settling, & &delt,dz,dx,znt, & &u,v,w,th,qv,qc,qi,qnc,qni, & - &qnwfa,qnifa, & + &qnwfa,qnifa,qnbca, & &p,exner,rho,T3D, & &xland,ts,qsfc,qcg,ps, & &ust,ch,hfx,qfx,rmol,wspd, & @@ -3894,6 +3944,7 @@ SUBROUTINE mynn_bl_driver( & &RQVBLTEN,RQCBLTEN,RQIBLTEN, & &RQNCBLTEN,RQNIBLTEN, & &RQNWFABLTEN,RQNIFABLTEN, & + &RQNBCABLTEN, & &exch_h,exch_m, & &Pblh,kpbl, & &el_pbl, & @@ -3916,7 +3967,8 @@ SUBROUTINE mynn_bl_driver( & &spp_pbl,pattern_spp_pbl, & &RTHRATEN, & &FLAG_QC,FLAG_QI,FLAG_QNC, & - &FLAG_QNI,FLAG_QNWFA,FLAG_QNIFA & + &FLAG_QNI,FLAG_QNWFA,FLAG_QNIFA, & + &FLAG_QNBCA & &,IDS,IDE,JDS,JDE,KDS,KDE & &,IMS,IME,JMS,JME,KMS,KME & &,ITS,ITE,JTS,JTE,KTS,KTE) @@ -3941,7 +3993,7 @@ SUBROUTINE mynn_bl_driver( & INTEGER, INTENT(in) :: icloud_bl LOGICAL, INTENT(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,& - FLAG_QNWFA,FLAG_QNIFA + FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA INTEGER,INTENT(IN) :: & & IDS,IDE,JDS,JDE,KDS,KDE & @@ -3970,7 +4022,7 @@ SUBROUTINE mynn_bl_driver( & REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,& &u,v,w,th,qv,p,exner,rho,T3D REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), OPTIONAL, INTENT(in)::& - &qc,qi,qni,qnc,qnwfa,qnifa + &qc,qi,qni,qnc,qnwfa,qnifa,qnbca REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,& &ch,rmol,ts,qsfc,qcg,ps,hfx,qfx,wspd,uoce,voce,vdfg,znt @@ -3982,7 +4034,7 @@ SUBROUTINE mynn_bl_driver( & REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: & &RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN,& &RQIBLTEN,RQNIBLTEN,RTHRATEN,RQNCBLTEN, & - &RQNWFABLTEN,RQNIFABLTEN + &RQNWFABLTEN,RQNIFABLTEN,RQNBCABLTEN REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: & &exch_h,exch_m @@ -4039,7 +4091,8 @@ SUBROUTINE mynn_bl_driver( & REAL, DIMENSION(KTS:KTE) :: thetav,sh,u1,v1,w1,p1,ex1,dz1,th1,tk1,rho1,& & qke1,tsq1,qsq1,cov1,qv1,qi1,qc1,du1,dv1,dth1,dqv1,dqc1,dqi1, & - & k_m1,k_h1,qni1,dqni1,qnc1,dqnc1,qnwfa1,qnifa1,dqnwfa1,dqnifa1 + & k_m1,k_h1,qni1,dqni1,qnc1,dqnc1,qnwfa1,qnifa1,qnbca1, & + dqnwfa1,dqnifa1,dqnbca1 !JOE: mass-flux variables REAL, DIMENSION(KTS:KTE) :: dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf @@ -4049,7 +4102,7 @@ SUBROUTINE mynn_bl_driver( & det_thl,det_sqv,det_sqc,det_u,det_v REAL,DIMENSION(KTS:KTE+1) :: s_aw1,s_awthl1,s_awqt1,& s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1,& - s_awqnc1,s_awqni1,s_awqnwfa1,s_awqnifa1 + s_awqnc1,s_awqni1,s_awqnwfa1,s_awqnifa1,s_awqnbca1 REAL, DIMENSION(KTS:KTE+1) :: zw REAL :: cpm,sqcg,flt,flq,flqv,flqc,pmz,phh,exnerg,zet,& @@ -4152,6 +4205,7 @@ SUBROUTINE mynn_bl_driver( & dqnc1(kts:kte)=0.0 dqnwfa1(kts:kte)=0.0 dqnifa1(kts:kte)=0.0 + dqnbca1(kts:kte)=0.0 qc_bl1D(kts:kte)=0.0 qi_bl1D(kts:kte)=0.0 cldfra_bl1D(kts:kte)=0.0 @@ -4362,6 +4416,7 @@ SUBROUTINE mynn_bl_driver( & dqnc1(k)=0.0 dqnwfa1(k)=0.0 dqnifa1(k)=0.0 + dqnbca1(k)=0.0 IF(PRESENT(qi) .AND. FLAG_QI)THEN qi1(k)= qi(i,k,j) sqi(k)= qi(i,k,j)/(1.+qv(i,k,j)) @@ -4424,6 +4479,11 @@ SUBROUTINE mynn_bl_driver( & ELSE qnifa1(k)=0.0 ENDIF + IF (PRESENT(qnbca) .AND. FLAG_QNBCA ) THEN + qnbca1(k)=qnbca(i,k,j) + ELSE + qnbca1(k)=0.0 + ENDIF p1(k) = p(i,k,j) ex1(k)= exner(i,k,j) el(k) = el_pbl(i,k,j) @@ -4455,6 +4515,7 @@ SUBROUTINE mynn_bl_driver( & s_awqni1(k)=0. s_awqnwfa1(k)=0. s_awqnifa1(k)=0. + s_awqnbca1(k)=0. sub_thl(k)=0. sub_sqv(k)=0. sub_u(k)=0. @@ -4513,6 +4574,7 @@ SUBROUTINE mynn_bl_driver( & s_awqni1(kte+1)=0. s_awqnwfa1(kte+1)=0. s_awqnifa1(kte+1)=0. + s_awqnbca1(kte+1)=0. #if (WRF_CHEM == 1) DO ic = 1,nchem s_awchem1(kte+1,ic)=0. @@ -4690,7 +4752,7 @@ SUBROUTINE mynn_bl_driver( & &bl_mynn_mixscalars, & &u1,v1,w1,th1,thl,thetav,tk1, & &sqw,sqv,sqc,qke1, & - &qnc1,qni1,qnwfa1,qnifa1, & + &qnc1,qni1,qnwfa1,qnifa1,qnbca1, & &ex1,Vt,Vq,sgm, & &ust(i,j),flt,flq,flqv,flqc, & &PBLH(i,j),KPBL(i,j),DX, & @@ -4706,6 +4768,7 @@ SUBROUTINE mynn_bl_driver( & & s_awu1,s_awv1,s_awqke1, & & s_awqnc1,s_awqni1, & & s_awqnwfa1,s_awqnifa1, & + & s_awqnbca1, & & sub_thl,sub_sqv, & & sub_u,sub_v, & & det_thl,det_sqv,det_sqc, & @@ -4717,7 +4780,7 @@ SUBROUTINE mynn_bl_driver( & & qc_bl1D_old,cldfra_bl1D_old, & & FLAG_QC,FLAG_QI, & & FLAG_QNC,FLAG_QNI, & - & FLAG_QNWFA,FLAG_QNIFA, & + & FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA,& & Psig_shcu(i,j), & & nupdraft(i,j),ktop_plume(i,j), & & maxmf(i,j),ztop_plume, & @@ -4764,7 +4827,7 @@ SUBROUTINE mynn_bl_driver( & &u1, v1, th1, tk1, qv1, & &qc1, qi1, qnc1, qni1, & &p1, ex1, thl, sqv, sqc, sqi, sqw,& - &qnwfa1, qnifa1, & + &qnwfa1, qnifa1, qnbca1, & &ust(i,j),flt,flq,flqv,flqc, & &wspd(i,j),qcg(i,j), & &uoce(i,j),voce(i,j), & @@ -4773,19 +4836,21 @@ SUBROUTINE mynn_bl_driver( & &dfm, dfh, dfq, & &Du1, Dv1, Dth1, Dqv1, & &Dqc1, Dqi1, Dqnc1, Dqni1, & - &Dqnwfa1, Dqnifa1, & + &Dqnwfa1, Dqnifa1, Dqnbca1, & &vdfg(i,j), diss_heat, & ! mass flux components &s_aw1,s_awthl1,s_awqt1, & &s_awqv1,s_awqc1,s_awu1,s_awv1, & &s_awqnc1,s_awqni1, & &s_awqnwfa1,s_awqnifa1, & + &s_awqnbca1, & &sub_thl,sub_sqv, & &sub_u,sub_v, & &det_thl,det_sqv,det_sqc, & &det_u,det_v, & &FLAG_QC,FLAG_QI,FLAG_QNC, & &FLAG_QNI,FLAG_QNWFA,FLAG_QNIFA, & + &FLAG_QNBCA, & &cldfra_bl1d, & &bl_mynn_cloudmix, & &bl_mynn_mixqt, & @@ -4839,11 +4904,13 @@ SUBROUTINE mynn_bl_driver( & IF (PRESENT(qni) .AND. FLAG_QNI) RQNIBLTEN(i,k,j)=dqni1(k) IF (PRESENT(qnwfa) .AND. FLAG_QNWFA) RQNWFABLTEN(i,k,j)=dqnwfa1(k) IF (PRESENT(qnifa) .AND. FLAG_QNIFA) RQNIFABLTEN(i,k,j)=dqnifa1(k) + IF (PRESENT(qnbca) .AND. FLAG_QNBCA) RQNBCABLTEN(i,k,j)=dqnbca1(k) ELSE IF (PRESENT(qnc) .AND. FLAG_QNC) RQNCBLTEN(i,k,j)=0. IF (PRESENT(qni) .AND. FLAG_QNI) RQNIBLTEN(i,k,j)=0. IF (PRESENT(qnwfa) .AND. FLAG_QNWFA) RQNWFABLTEN(i,k,j)=0. IF (PRESENT(qnifa) .AND. FLAG_QNIFA) RQNIFABLTEN(i,k,j)=0. + IF (PRESENT(qnbca) .AND. FLAG_QNBCA) RQNBCABLTEN(i,k,j)=0. ENDIF IF(icloud_bl > 0)THEN @@ -5196,7 +5263,7 @@ SUBROUTINE DMP_mf( & & scalar_opt, & & u,v,w,th,thl,thv,tk, & & qt,qv,qc,qke, & - qnc,qni,qnwfa,qnifa, & + qnc,qni,qnwfa,qnifa,qnbca, & & exner,vt,vq,sgm, & & ust,flt,flq,flqv,flqc, & & pblh,kpbl,DX,landsea,ts, & @@ -5210,6 +5277,7 @@ SUBROUTINE DMP_mf( & & s_awu,s_awv,s_awqke, & & s_awqnc,s_awqni, & & s_awqnwfa,s_awqnifa, & + & s_awqnbca, & & sub_thl,sub_sqv, & & sub_u,sub_v, & & det_thl,det_sqv,det_sqc, & @@ -5223,7 +5291,7 @@ SUBROUTINE DMP_mf( & ! inputs - flags for moist arrays & F_QC,F_QI, & F_QNC,F_QNI, & - & F_QNWFA,F_QNIFA, & + & F_QNWFA,F_QNIFA,F_QNBCA, & & Psig_shcu, & ! output info &nup2,ktop,maxmf,ztop, & @@ -5243,11 +5311,11 @@ SUBROUTINE DMP_mf( & REAL, DIMENSION(KTS:KTE) :: rstoch_col REAL,DIMENSION(KTS:KTE), INTENT(IN) :: U,V,W,TH,THL,TK,QT,QV,QC,& - exner,dz,THV,P,qke,qnc,qni,qnwfa,qnifa + exner,dz,THV,P,qke,qnc,qni,qnwfa,qnifa,qnbca REAL,DIMENSION(KTS:KTE+1), INTENT(IN) :: ZW !height at full-sigma REAL, INTENT(IN) :: DT,UST,FLT,FLQ,FLQV,FLQC,PBLH,& DX,Psig_shcu,landsea,ts - LOGICAL, OPTIONAL :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA + LOGICAL, OPTIONAL :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA,F_QNBCA ! outputs - updraft properties REAL,DIMENSION(KTS:KTE), INTENT(OUT) :: edmf_a,edmf_w, & @@ -5267,6 +5335,7 @@ SUBROUTINE DMP_mf( & s_awqni, & s_awqnwfa, & s_awqnifa, & + s_awqnbca, & s_awu, & s_awv, & s_awqke, s_aw2 @@ -5281,7 +5350,7 @@ SUBROUTINE DMP_mf( & ! first model layer REAL,DIMENSION(KTS:KTE+1,1:NUP) :: UPW,UPTHL,UPQT,UPQC,UPQV, & UPA,UPU,UPV,UPTHV,UPQKE,UPQNC, & - UPQNI,UPQNWFA,UPQNIFA + UPQNI,UPQNWFA,UPQNIFA,UPQNBCA ! entrainment variables REAL,DIMENSION(KTS:KTE,1:NUP) :: ENT,ENTf INTEGER,DIMENSION(KTS:KTE,1:NUP) :: ENTi @@ -5289,7 +5358,7 @@ SUBROUTINE DMP_mf( & INTEGER :: K,I,k50 REAL :: fltv,wstar,qstar,thstar,sigmaW,sigmaQT,sigmaTH,z0, & pwmin,pwmax,wmin,wmax,wlv,Psig_w,maxw,maxqc,wpbl - REAL :: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,QNCn,QNIn,QNWFAn,QNIFAn, & + REAL :: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,QNCn,QNIn,QNWFAn,QNIFAn,QNBCAn, & Wn2,Wn,EntEXP,EntW,BCOEFF,THVkm1,THVk,Pk ! w parameters @@ -5394,6 +5463,7 @@ SUBROUTINE DMP_mf( & UPQNI=0. UPQNWFA=0. UPQNIFA=0. + UPQNBCA=0. #if (WRF_CHEM == 1) IF (bl_mynn_mixchem == 1) THEN UPCHEM(KTS:KTE+1,1:NUP,1:nchem)=0.0 @@ -5425,6 +5495,7 @@ SUBROUTINE DMP_mf( & s_awqni=0. s_awqnwfa=0. s_awqnifa=0. + s_awqnbca=0. #if (WRF_CHEM == 1) IF (bl_mynn_mixchem == 1) THEN s_awchem(kts:kte+1,1:nchem) = 0.0 @@ -5626,6 +5697,7 @@ SUBROUTINE DMP_mf( & UPQNI(1,I)=(QNI(KTS)*DZ(KTS+1)+QNI(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPQNWFA(1,I)=(QNWFA(KTS)*DZ(KTS+1)+QNWFA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPQNIFA(1,I)=(QNIFA(KTS)*DZ(KTS+1)+QNIFA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) + UPQNBCA(1,I)=(QNBCA(KTS)*DZ(KTS+1)+QNBCA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) ENDDO #if (WRF_CHEM == 1) @@ -5686,6 +5758,7 @@ SUBROUTINE DMP_mf( & QNIn=UPQNI(k-1,I)*(1.-EntExp) + QNI(k)*EntExp QNWFAn=UPQNWFA(k-1,I)*(1.-EntExp) + QNWFA(k)*EntExp QNIFAn=UPQNIFA(k-1,I)*(1.-EntExp) + QNIFA(k)*EntExp + QNBCAn=UPQNBCA(k-1,I)*(1.-EntExp) + QNBCA(k)*EntExp !capture the updated qc, qt & thl modified by entranment alone, !since they will be modified later if condensation occurs. @@ -5834,6 +5907,7 @@ SUBROUTINE DMP_mf( & UPQNI(K,I)=QNIn UPQNWFA(K,I)=QNWFAn UPQNIFA(K,I)=QNIFAn + UPQNBCA(K,I)=QNBCAn UPA(K,I)=UPA(K-1,I) #if (WRF_CHEM == 1) IF (bl_mynn_mixchem == 1) THEN @@ -5923,6 +5997,7 @@ SUBROUTINE DMP_mf( & s_awqni(k+1)= s_awqni(K+1) + UPA(K,i)*UPW(K,i)*UPQNI(K,i)*Psig_w s_awqnwfa(k+1)= s_awqnwfa(K+1) + UPA(K,i)*UPW(K,i)*UPQNWFA(K,i)*Psig_w s_awqnifa(k+1)= s_awqnifa(K+1) + UPA(K,i)*UPW(K,i)*UPQNIFA(K,i)*Psig_w + s_awqnbca(k+1)= s_awqnbca(K+1) + UPA(K,i)*UPW(K,i)*UPQNBCA(K,i)*Psig_w ENDDO ENDDO ENDIF @@ -5952,6 +6027,7 @@ SUBROUTINE DMP_mf( & s_awqni= s_awqni*adjustment s_awqnwfa= s_awqnwfa*adjustment s_awqnifa= s_awqnifa*adjustment + s_awqnbca= s_awqnbca*adjustment IF (momentum_opt > 0) THEN s_awu = s_awu*adjustment s_awv = s_awv*adjustment diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 55a0cfbb1c..6bdda02bbc 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -29,7 +29,8 @@ SUBROUTINE microphysics_driver( & ,gmt,xtime & !====================== !Variables required for CAMMGMP Scheme - ,dlf,dlf2,t_phy,p_hyd,p8w_hyd,tke_pbl,z_at_w,qfx & + ,dlf,dlf2,t_phy,p_hyd,p8w_hyd,tke_pbl,pblh & + ,z_at_mass,z_at_w,qfx & ,rliq,turbtype3d,smaw3d,wsedl3d,cldfra_old_mp & ,cldfra_mp,cldfra_mp_all,lradius,iradius & ,cldfrai,cldfral,cldfra_conv & @@ -49,8 +50,8 @@ SUBROUTINE microphysics_driver( & ,qndrop_curr,qni_curr,qh_curr,qnh_curr & ,qzr_curr,qzi_curr,qzs_curr,qzg_curr,qzh_curr & ,qns_curr,qnr_curr,qng_curr,qnn_curr,qnc_curr & - ,qnwfa_curr,qnifa_curr & ! for water/ice-friendly aerosols - ,f_qnwfa,f_qnifa & ! for water/ice-friendly aerosols + ,qnwfa_curr,qnifa_curr,qnbca_curr & ! for water/ice-friendly/black carbon aerosols + ,f_qnwfa,f_qnifa,f_qnbca & ! for water/ice-friendly/black carbon aerosols ,qvolg_curr,qvolh_curr & ,qdcn_curr,qtcn_curr,qccn_curr,qrcn_curr & ! for ntu3m ,qnin_curr,fi_curr,fs_curr,vi_curr,vs_curr & ! for ntu3m @@ -105,7 +106,8 @@ SUBROUTINE microphysics_driver( & ,rainprod, evapprod & ,qv_b4mp, qc_b4mp, qi_b4mp, qs_b4mp & #endif - ,qnwfa2d, qnifa2d & ! for water/ice-friendly aerosols + ,qnwfa2d, qnifa2d, qnbca2d & ! for water/ice-friendly/black carbon aerosols + ,qnocbb2d, qnbcbb2d & ! for biomass burning aerosols ,refl_10cm & ! HM, 9/22/09, add for refl ,vmi3d & ! for P3 ,di3d & ! for P3 @@ -494,8 +496,10 @@ SUBROUTINE microphysics_driver( & t_phy, & !Temprature at the mid points (K) p_hyd, & !Hydrostatic pressure(Pa) p8w_hyd, & !Hydrostatic Pressure at level interface (Pa) + z_at_mass, & !Height above sea level at grid cell center (m) z_at_w, & !Height above sea level at layer interfaces (m) tke_pbl, & !Turbulence kinetic energy + pblh, & !Planetary boundary layer height (m) turbtype3d, & !Turbulent interface types [ no unit ] smaw3d, & !Normalized Galperin instability function for momentum [no units] alt, & !inverse density(m3/kg) @@ -590,7 +594,8 @@ SUBROUTINE microphysics_driver( & itype,itype_2,itype_3 ! for Jensen ISHMAEL LOGICAL, OPTIONAL, INTENT(IN ) :: channel_switch REAL, OPTIONAL, INTENT(INOUT ) :: naer ! aerosol number concentration (/kg) - REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) , OPTIONAL :: qnwfa2d, qnifa2d + REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) , OPTIONAL :: qnwfa2d, qnifa2d, qnbca2d, & + qnocbb2d, qnbcbb2d REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(INOUT ) :: & @@ -614,7 +619,7 @@ SUBROUTINE microphysics_driver( & ,kext_qic,kext_qip,kext_qid,tempc,height & ,kext_ft_qic,kext_ft_qip,kext_ft_qid & ,kext_ft_qs,kext_ft_qg & - ,qnwfa_curr,qnifa_curr & ! Added by G. Thompson + ,qnwfa_curr,qnifa_curr,qnbca_curr & ! Added by G. Thompson ,qvolg_curr,qvolh_curr, qrimef_curr REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & @@ -690,7 +695,7 @@ SUBROUTINE microphysics_driver( & ,f_qvoli,f_qaoli & ! for Jensen ISHMAEL ,f_qvoli2,f_qaoli2 & ! for Jensen ISHMAEL ,f_qi3,f_qni3,f_qvoli3,f_qaoli3 & ! for Jensen ISHMAEL - ,f_qnwfa, f_qnifa ! Added by G. Thompson + ,f_qnwfa, f_qnifa, f_qnbca ! Added by G. Thompson LOGICAL, OPTIONAL, INTENT(IN) :: diagflag @@ -1018,12 +1023,22 @@ SUBROUTINE microphysics_driver( & NC=qnc_curr, & NWFA=qnwfa_curr, & NIFA=qnifa_curr, & + NBCA=qnbca_curr, & NWFA2D=qnwfa2d, & NIFA2D=qnifa2d, & + NBCA2D=qnbca2d, & + NOCBB2D=qnocbb2d, & + NBCBB2D=qnbcbb2d, & + aer_init_opt=config_flags%aer_init_opt, & + wif_input_opt=config_flags%wif_input_opt, & + aer_fire_emit_opt=config_flags%aer_fire_emit_opt, & + wif_fire_inj=config_flags%wif_fire_inj, & TH=th, & PII=pi_phy, & P=p, & W=w, & + PBLH=PBLH, & + z_at_mass=z_at_mass, & DZ=dz8w, & DT_IN=dt, & ITIMESTEP=itimestep, & diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index ecc4055a1d..c1e91da495 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -85,6 +85,8 @@ MODULE module_mp_thompson 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. @@ -395,7 +397,8 @@ MODULE module_mp_thompson CONTAINS - SUBROUTINE thompson_init(hgt, orho, nwfa2d, nwfa, nifa, & + SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, & + nwfa, nifa, nbca, wif_input_opt, frc_urb2d, & dx, dy, is_start, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & @@ -410,16 +413,18 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nwfa, nifa, & !..OPTIONAL variables that control application of aerosol-aware scheme - REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa - REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d + 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, max_test, z1 + REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1 LOGICAL:: micro_init, has_CCN, has_IN is_aerosol_aware = .FALSE. @@ -514,6 +519,47 @@ SUBROUTINE thompson_init(hgt, orho, nwfa2d, nwfa, nifa, & 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 @@ -974,9 +1020,13 @@ END SUBROUTINE thompson_init !+---+-----------------------------------------------------------------+ !..This is a wrapper routine designed to transfer values from 3D to 1D. !+---+-----------------------------------------------------------------+ - SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & - nwfa, nifa, nwfa2d, nifa2d, & - th, pii, p, w, dz, dt_in, itimestep, & + SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & + nwfa, nifa, nbca, nwfa2d, nifa2d, nbca2d, & + nocbb2d, nbcbb2d, & + aer_init_opt, wif_input_opt, & + aer_fire_emit_opt, wif_fire_inj, & + th, pii, p, w, z_at_mass, dz, & + dt_in, itimestep, PBLH, & RAINNC, RAINNCV, & SNOWNC, SNOWNCV, & GRAUPELNC, GRAUPELNCV, SR, & @@ -999,8 +1049,9 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & qv, qc, qr, qi, qs, qg, ni, nr, th REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - nc, nwfa, nifa - REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d + nc, nwfa, nifa, nbca + REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, & + nbca2d, nocbb2d, nbcbb2d 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 @@ -1010,6 +1061,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & #endif REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & pii, p, w, dz + REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & + z_at_mass + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(IN):: & + PBLH REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & RAINNC, RAINNCV, SR REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: & @@ -1018,6 +1073,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & refl_10cm REAL, INTENT(IN):: dt_in INTEGER, INTENT(IN):: itimestep + INTEGER, OPTIONAL, INTENT(IN):: & + aer_init_opt, wif_input_opt, aer_fire_emit_opt, wif_fire_inj #if ( WRF_CHEM == 1 ) LOGICAL, INTENT(in) :: wetscav_on #endif @@ -1025,7 +1082,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !..Local variables REAL, DIMENSION(kts:kte):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, & + nr1d, nc1d, nwfa1d, nifa1d, nbca1d, & t1d, p1d, w1d, dz1d, rho, dBZ REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d #if ( WRF_CHEM == 1 ) @@ -1035,8 +1092,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic REAL:: dt, pptrain, pptsnow, pptgraul, pptice REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max - REAL:: nwfa1 - INTEGER:: i, j, k + REAL:: nwfa1, noc_emit, nbc_emit + INTEGER:: i, j, k, 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 @@ -1139,6 +1196,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & 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 @@ -1146,12 +1208,14 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & 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 nwfa1 = 11.1E6 endif - call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, & + call mp_thompson(aer_init_opt, wif_input_opt, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, & pptrain, pptsnow, pptgraul, pptice, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod1d, evapprod1d, & @@ -1175,19 +1239,60 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & 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 + +!..Add organic carbon (OC) and black carbon (BC) biomass burning emissions + if (aer_fire_emit_opt .gt. 0) then + if (wif_fire_inj .eq. 1) then ! distribute aerosol through PBL + k_inj = kts + do while (z_at_mass(i,k_inj,j) .lt. PBLH(i,j)) + k_inj = k_inj + 1 + end do + noc_emit = (nocbb2d(i,j)*dt_in)/k_inj ! num oc aer per grid cell + if (aer_fire_emit_opt .eq. 2) then + nbc_emit = (nbcbb2d(i,j)*dt_in)/k_inj ! num bc aer per grid cell + end if + do k = kts, k_inj + nwfa1d(k) = nwfa1d(k) + noc_emit + if (aer_fire_emit_opt .eq. 2) then + nbca1d(k) = nbca1d(k) + nbc_emit + end if + end do + else ! emit at sfc + nwfa1d(kts) = nwfa1d(kts) + nocbb2d(i,j)*dt_in + if (aer_fire_emit_opt .eq. 2) then + nbca1d(kts) = nbca1d(kts) + nbcbb2d(i,j)*dt_in + end if + endif + end if +!.. +!..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 @@ -1355,8 +1460,9 @@ END SUBROUTINE mp_gt_driver !.. Thompson et al. (2004, 2008). !+---+-----------------------------------------------------------------+ ! - subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, & + subroutine mp_thompson (aer_init_opt, wif_input_opt, & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dzq, & pptrain, pptsnow, pptgraul, pptice, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod, evapprod, & @@ -1366,10 +1472,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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, ni1d, & - nr1d, nc1d, nwfa1d, nifa1d, t1d + nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN):: dt @@ -1381,7 +1488,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !..Local variables REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & - qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten + qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten, & + nbcaten DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd @@ -1389,7 +1497,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnc_scw, pnc_gcw DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, & - pnd_rcd, pnd_scd, pnd_gcd + pnd_rcd, pnd_scd, pnd_gcd, pnb_rcb, pnb_scb, pnb_gcb DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, & prr_rcg, prr_sml, prr_gml, & @@ -1415,7 +1523,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 REAL, DIMENSION(kts:kte):: temp, pres, qv - REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa + REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, 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 @@ -1508,6 +1616,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ncten(k) = 0. nwfaten(k) = 0. nifaten(k) = 0. + nbcaten(k) = 0. prw_vcd(k) = 0. @@ -1576,6 +1685,10 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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 @@ -1608,8 +1721,29 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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)) - 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 (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 if (qc1d(k) .gt. R1) then no_micro = .false. @@ -1948,6 +2082,16 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnd_rcd(k) = rhof(k)*t1_qr_qc*Ef_ra*nifa(k)*N0_r(k) & *((lamr+fv_r)**(-cre(9))) pnd_rcd(k) = MIN(DBLE(nifa(k)*odts), pnd_rcd(k)) + + if (present(wif_input_opt)) then + if (wif_input_opt .eq. 2) then + Ef_ra = Eff_aero(mvd_r(k),0.0236E-6, & + visco(k),rho(k),temp(k),'r') + pnb_rcb(k) = rhof(k)*t1_qr_qc*Ef_ra*nbca(k)*N0_r(k) & + *((lamr+fv_r)**(-cre(9))) + pnb_rcb(k) = MIN(DBLE(nbca(k)*odts), pnb_rcb(k)) + endif + endif endif enddo @@ -2154,6 +2298,14 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & Ef_sa = Eff_aero(xDs,0.8E-6,visco(k),rho(k),temp(k),'s') pnd_scd(k) = rhof(k)*t1_qs_qc*Ef_sa*nifa(k)*smoe(k) pnd_scd(k) = MIN(DBLE(nifa(k)*odts), pnd_scd(k)) + + if (present(wif_input_opt)) then + if (wif_input_opt .eq. 2) then + Ef_sa = Eff_aero(xDs,0.0236E-6,visco(k),rho(k),temp(k),'s') + pnb_scb(k) = rhof(k)*t1_qs_qc*Ef_sa*nbca(k)*smoe(k) + pnb_scb(k) = MIN(DBLE(nbca(k)*odts), pnb_scb(k)) + endif + endif endif if (rg(k) .gt. r_g(1)) then xDg = (bm_g + mu_g + 1.) * ilamg(k) @@ -2166,6 +2318,15 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnd_gcd(k) = rhof(k)*t1_qg_qc*Ef_ga*nifa(k)*N0_g(k) & *ilamg(k)**cge(9) pnd_gcd(k) = MIN(DBLE(nifa(k)*odts), pnd_gcd(k)) + + if (present(wif_input_opt)) then + if (wif_input_opt .eq. 2) then + Ef_ga = Eff_aero(xDg,0.0236E-6,visco(k),rho(k),temp(k),'g') + pnb_gcb(k) = rhof(k)*t1_qg_qc*Ef_ga*nbca(k)*N0_g(k) & + *ilamg(k)**cge(9) + pnb_gcb(k) = MIN(DBLE(nbca(k)*odts), pnb_gcb(k)) + endif + endif endif !..Rain collecting snow. Cannot assume Wisner (1972) approximation @@ -2609,6 +2770,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + pna_gca(k) + pni_iha(k)) * orho nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) & + pnd_gcd(k)) * orho + if (wif_input_opt .eq. 2) then + nbcaten(k) = nbcaten(k) - (pnb_rcb(k) + pnb_scb(k) & + + pnb_gcb(k)) * orho + else + nbcaten(k) = 0.0 + endif if (dustyIce) then nifaten(k) = nifaten(k) - pni_inu(k)*orho else @@ -3499,10 +3666,35 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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)) - 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 (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 diff --git a/phys/module_pbl_driver.F b/phys/module_pbl_driver.F index 9cf044a0f4..7e6213b77c 100644 --- a/phys/module_pbl_driver.F +++ b/phys/module_pbl_driver.F @@ -87,6 +87,7 @@ SUBROUTINE pbl_driver( & ! Optional aerosol ,qnwfa_curr,f_qnwfa & ,qnifa_curr,f_qnifa & + ,qnbca_curr,f_qnbca & ! Optional moisture tracer flags ,f_qv,f_qc,f_qr & ,f_qi,f_qs,f_qg & @@ -142,7 +143,7 @@ SUBROUTINE pbl_driver( & FITCHSCHEME,SHINHONGSCHEME, & TEMFPBLSCHEME,GBMPBLSCHEME,EEPSSCHEME, & CAMMGMPSCHEME,p_qi,p_qni,p_qnc,param_first_scalar,& !CAMMGMPSCHEME, p_qni,p_qnc is used for camuwpbl scheme - p_qnwfa,p_qnifa + p_qnwfa,p_qnifa,p_qnbca #if ( WRFPLUS == 1 ) USE module_state_description, ONLY : SURFDRAGSCHEME #endif @@ -564,7 +565,7 @@ SUBROUTINE pbl_driver( & & INTENT(OUT) :: maxMF REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & - INTENT(INOUT) :: qnwfa_curr,qnifa_curr + INTENT(INOUT) :: qnwfa_curr,qnifa_curr,qnbca_curr REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & & INTENT(INOUT) :: exch_tke ! for GBM PBL scheme @@ -656,7 +657,8 @@ SUBROUTINE pbl_driver( & ,f_qnc & !used in CAMUWPBL ,f_qni & !used in CAMUWPBL ,f_qnwfa & - ,f_qnifa + ,f_qnifa & + ,f_qnbca REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, INTENT(INOUT) :: & @@ -669,7 +671,7 @@ SUBROUTINE pbl_driver( & ,rqiblten,rqsblten,rqgblten,rqniblten !rqniblten used in CAMUWPBL REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: & !local: added for MYNN - rqncblten,rqnwfablten,rqnifablten + rqncblten,rqnwfablten,rqnifablten,rqnbcablten REAL, DIMENSION( ims:ime, jms:jme ) , & OPTIONAL , & @@ -783,7 +785,7 @@ SUBROUTINE pbl_driver( & LOGICAL :: flag_bep LOGICAL :: flag_myjsfc LOGICAL :: flag_qv, flag_qc, flag_qr, flag_qi, flag_qs, flag_qg, & - flag_qnc,flag_qni, flag_qnwfa, flag_qnifa + flag_qnc,flag_qni, flag_qnwfa, flag_qnifa, flag_qnbca CHARACTER*256 :: message REAL :: next_bl_time LOGICAL :: run_param , doing_adapt_dt , decided @@ -829,6 +831,7 @@ SUBROUTINE pbl_driver( & flag_qni = .FALSE. ; IF ( PRESENT( F_QNI ) ) flag_qni = F_QNI !Used in CAMUWPBL flag_qnwfa = .FALSE. ; IF ( PRESENT( F_QNWFA ) ) flag_qnwfa = F_QNWFA flag_qnifa = .FALSE. ; IF ( PRESENT( F_QNIFA ) ) flag_qnifa = F_QNIFA + flag_qnbca = .FALSE. ; IF ( PRESENT( F_QNBCA ) ) flag_qnbca = F_QNBCA if (bl_pbl_physics .eq. 0) return @@ -1605,7 +1608,7 @@ SUBROUTINE pbl_driver( & &u=u_phy,v=v_phy,w=w,th=th_phy,qv=qv_curr, & &qc=qc_curr,qi=qi_curr, & &qnc=qnc_curr,qni=qni_curr, & - &QNWFA=qnwfa_curr,QNIFA=qnifa_curr, & + &QNWFA=qnwfa_curr,QNIFA=qnifa_curr,QNBCA=qnbca_curr, & &p=p_phy,exner=pi_phy,rho=rho,T3D=t_phy, & &xland=xland,ts=tsk,qsfc=qsfc,qcg=qcg,ps=psfc, & &ust=ust,ch=ch,hfx=hfx,qfx=qfx,rmol=rmol,wspd=wspd, & @@ -1622,6 +1625,7 @@ SUBROUTINE pbl_driver( & &RQVBLTEN=rqvblten,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten,& &RQNCBLTEN=rqncblten,RQNIBLTEN=rqniblten, & &RQNWFABLTEN=rqnwfablten,RQNIFABLTEN=rqnifablten, & + &RQNBCABLTEN=rqnbcablten, & &EXCH_H=exch_h,EXCH_M=exch_m, & &pblh=pblh,KPBL=KPBL & &,el_pbl=el_pbl & @@ -1651,6 +1655,7 @@ SUBROUTINE pbl_driver( & &,FLAG_QC=flag_qc,FLAG_QI=flag_qi & &,FLAG_QNC=flag_qnc,FLAG_QNI=flag_qni & &,FLAG_QNWFA=flag_qnwfa,FLAG_QNIFA=flag_qnifa & + &,FLAG_QNBCA=flag_qnbca & ,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 & @@ -1745,6 +1750,16 @@ SUBROUTINE pbl_driver( & enddo enddo ENDIF + IF (PRESENT( qnbca_curr ) .AND. Flag_qnbca) THEN + !print*,"Updating qnbca after mynn-edmf",P_QNBCA + DO j=j_start(ij),j_end(ij) + DO k=kts,kte + DO i=i_start(ij),i_end(ij) + scalar_tend(i,k,j,P_QNBCA)=RQNBCABLTEN(i,k,j) + enddo + enddo + enddo + ENDIF ENDIF CASE (EEPSSCHEME) diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index d6e65c01ec..8725908b3a 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -51,7 +51,8 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & STEPBL,STEPRA,STEPCU, & W0AVG, RAINNC, RAINC, RAINCV, RAINNCV, & SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV, & - z_at_q, inv_dens, qnwfa2d, scalar, num_sc, & ! G. Thompson + z_at_q, inv_dens, qnwfa2d, qnbca2d, & ! G. Thompson + scalar, num_sc, & ! G. Thompson re_cloud, re_ice, re_snow, & ! G. Thompson has_reqc, has_reqi, has_reqs, & ! G. Thompson #if ( EM_CORE == 1 ) @@ -461,7 +462,7 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: z_at_q ! G. Thompson REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: inv_dens ! G. Thompson - REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qnwfa2d ! G. Thompson + REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qnwfa2d, qnbca2d ! G. Thompson INTEGER, INTENT(IN) :: num_sc ! G. Thompson REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_sc), INTENT(INOUT) :: scalar ! G. Thompson @@ -1664,7 +1665,8 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & nssl_rho_qh, nssl_rho_qhl, & nssl_rho_qs, & ccn_conc, & ! RAS - z_at_q, inv_dens, qnwfa2d, scalar, num_sc, & ! G. Thompson + z_at_q, inv_dens, qnwfa2d, qnbca2d, & ! G. Thompson + frc_urb2d, scalar, num_sc, & ! G. Thompson ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) @@ -4408,7 +4410,8 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, nssl_rho_qh, nssl_rho_qhl, & nssl_rho_qs, & ccn_conc, & ! RAS - z_at_q, inv_dens, qnwfa2d, scalar, num_sc, & ! G. Thompson + z_at_q, inv_dens, qnwfa2d, qnbca2d, & ! G. Thompson + frc_urb2d, scalar, num_sc, & ! G. Thompson ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) @@ -4476,7 +4479,8 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, REAL , DIMENSION(:) ,INTENT(INOUT) :: mp_restart_state,tbpvs_state,tbpvs0_state LOGICAL , INTENT(IN) :: allowed_to_read REAL, INTENT(INOUT) :: ccn_conc ! RAS - REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT):: qnwfa2d ! G. Thompson + REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT):: qnwfa2d, qnbca2d ! G. Thompson + REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN):: frc_urb2d REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: z_at_q ! G. Thompson REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: inv_dens ! G. Thompson INTEGER, INTENT(IN) :: num_sc ! G. Thompson @@ -4584,9 +4588,12 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, IF(start_of_simulation.or.restart.or.config_flags%cycling) & CALL thompson_init(HGT=z_at_q, & ORHO=inv_dens, & - NWFA2D=qnwfa2d, & + 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, & diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index c63dbee151..d076fe7f85 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -101,6 +101,8 @@ SUBROUTINE radiation_driver ( & , QNDROP, F_QNDROP & ,QNIFA,F_QNIFA & ,QNWFA,F_QNWFA & + ,QNBCA & + ,wif_input_opt & ,qc_tot, qi_tot & ,ACSWUPT ,ACSWUPTC & ,ACSWDNT ,ACSWDNTC & @@ -798,12 +800,13 @@ SUBROUTINE radiation_driver ( & INTENT(INOUT ) :: & pb & ,qv,qc,qr,qi,qs,qg,qh,qndrop, & - qnifa,qnwfa, & ! Trude + qnifa,qnwfa,qnbca, & ! Trude qi2,qi3 ! for P3 LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qh,f_qndrop,& f_qnifa,f_qnwfa, & ! trude f_qi2,f_qi3 ! for P3 + INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt ! Solar diag REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT), OPTIONAL :: qc_tot, qi_tot !shbaek @@ -1508,7 +1511,7 @@ SUBROUTINE radiation_driver ( & ENDIF ENDIF -! Allocate aerosol arrays used by aer_opt = 3 option, and explicit AOD from QNWFA+QNIFA (Trude) +! Allocate aerosol arrays used by aer_opt = 3 option, and explicit AOD from QNWFA+QNIFA(+QNBCA) (Trude) IF (PRESENT(f_qnwfa) .AND. PRESENT(f_qnifa) .AND. PRESENT(taod5503d) .AND. PRESENT(taod5502d)) THEN IF (F_QNWFA .AND. aer_opt.eq.3 .AND. & (sw_physics.eq.RRTMG_SWSCHEME & @@ -1519,7 +1522,7 @@ SUBROUTINE radiation_driver ( & .OR. sw_physics.eq.RRTMK_SWSCHEME & #endif )) THEN - CALL wrf_debug (150, 'DEBUG-GT: computing 3D AOD from QNWFA+QNIFA') + CALL wrf_debug (150, 'DEBUG-GT: computing 3D AOD from QNWFA+QNIFA(+QNBCA)') allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) @@ -1539,8 +1542,8 @@ SUBROUTINE radiation_driver ( & end do end do - call gt_aod (p, DZ8W, t, qv, qnwfa, qnifa, taod5503d, & - ims,ime, jms,jme, kms,kme,its,ite, jts,jte, kts,kte) + call gt_aod (p, DZ8W, t, qv, qnwfa, qnifa, qnbca, taod5503d, & + wif_input_opt, ims,ime, jms,jme, kms,kme,its,ite, jts,jte, kts,kte) do j=jts,jte do i=its,ite @@ -5063,8 +5066,8 @@ END SUBROUTINE aer_p_int !+---+-----------------------------------------------------------------+ - SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & - & ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte) + SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa,nbca, taod5503d, & + & wif_input_opt, ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte) USE module_mp_thompson, only: RSLF @@ -5075,97 +5078,133 @@ SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: & & t_phy,p_phy, DZ8W, & - & qvapor, nifa, nwfa + & qvapor, nifa, nwfa, nbca REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT):: taod5503d - + INTEGER, INTENT(IN):: wif_input_opt + !..Local variables. - REAL, DIMENSION(its:ite,kts:kte,jts:jte):: AOD_wfa, AOD_ifa - REAL:: RH, a_RH,b_RH, rh_d,rh_f, rhoa,qvsat, unit_bext1,unit_bext3 + REAL, DIMENSION(its:ite,kts:kte,jts:jte):: AOD_wfa, AOD_ifa, AOD_bca + REAL:: RH, a_RH,b_RH, rh_d,rh_f, rhoa,qvsat, unit_bext1,unit_bext3,unit_bext4 REAL:: ntemp INTEGER :: i, k, j, RH_idx, RH_idx1, RH_idx2, t_idx INTEGER, PARAMETER:: rind=8 REAL, DIMENSION(rind), PARAMETER:: rh_arr = & & (/10., 60., 70., 80., 85., 90., 95., 99.8/) - REAL, DIMENSION(rind,4,2) :: lookup_tabl ! RH, temp, water-friendly, ice-friendly + REAL, DIMENSION(rind,4,3) :: lookup_tabl ! RH, temp, water-friendly, ice-friendly lookup_tabl(1,1,1) = 5.73936E-15 lookup_tabl(1,1,2) = 2.63577E-12 + lookup_tabl(1,1,3) = 7.83852E-16 lookup_tabl(1,2,1) = 5.73936E-15 lookup_tabl(1,2,2) = 2.63577E-12 + lookup_tabl(1,2,3) = 7.83803E-16 lookup_tabl(1,3,1) = 5.73936E-15 lookup_tabl(1,3,2) = 2.63577E-12 + lookup_tabl(1,3,3) = 7.83770E-16 lookup_tabl(1,4,1) = 5.73936E-15 lookup_tabl(1,4,2) = 2.63577E-12 + lookup_tabl(1,4,3) = 7.83692E-16 lookup_tabl(2,1,1) = 6.93515E-15 lookup_tabl(2,1,2) = 2.72095E-12 + lookup_tabl(2,1,3) = 9.96291E-16 lookup_tabl(2,2,1) = 6.93168E-15 - lookup_tabl(2,2,2) = 2.72092E-12 + lookup_tabl(2,2,2) = 2.72092E-12 + lookup_tabl(2,2,3) = 9.94950E-16 lookup_tabl(2,3,1) = 6.92570E-15 - lookup_tabl(2,3,2) = 2.72091E-12 + lookup_tabl(2,3,2) = 2.72091E-12 + lookup_tabl(2,3,3) = 9.93238E-16 lookup_tabl(2,4,1) = 6.91833E-15 lookup_tabl(2,4,2) = 2.72087E-12 + lookup_tabl(2,4,3) = 9.91346E-16 lookup_tabl(3,1,1) = 7.24707E-15 lookup_tabl(3,1,2) = 2.77219E-12 + lookup_tabl(3,1,3) = 1.11936E-15 lookup_tabl(3,2,1) = 7.23809E-15 lookup_tabl(3,2,2) = 2.77222E-12 + lookup_tabl(3,2,3) = 1.11687E-15 lookup_tabl(3,3,1) = 7.23108E-15 lookup_tabl(3,3,2) = 2.77201E-12 + lookup_tabl(3,3,3) = 1.11447E-15 lookup_tabl(3,4,1) = 7.21800E-15 lookup_tabl(3,4,2) = 2.77111E-12 + lookup_tabl(3,4,3) = 1.11061E-15 lookup_tabl(4,1,1) = 8.95130E-15 lookup_tabl(4,1,2) = 2.87263E-12 + lookup_tabl(4,1,3) = 1.36116E-15 lookup_tabl(4,2,1) = 9.01582E-15 lookup_tabl(4,2,2) = 2.87252E-12 + lookup_tabl(4,2,3) = 1.35479E-15 lookup_tabl(4,3,1) = 9.13216E-15 lookup_tabl(4,3,2) = 2.87241E-12 + lookup_tabl(4,3,3) = 1.34787E-15 lookup_tabl(4,4,1) = 9.16219E-15 lookup_tabl(4,4,2) = 2.87211E-12 + lookup_tabl(4,4,3) = 1.33910E-15 lookup_tabl(5,1,1) = 1.06695E-14 lookup_tabl(5,1,2) = 2.96752E-12 + lookup_tabl(5,1,3) = 1.58848E-15 lookup_tabl(5,2,1) = 1.06370E-14 lookup_tabl(5,2,2) = 2.96726E-12 + lookup_tabl(5,2,3) = 1.57854E-15 lookup_tabl(5,3,1) = 1.05999E-14 lookup_tabl(5,3,2) = 2.96702E-12 + lookup_tabl(5,3,3) = 1.56648E-15 lookup_tabl(5,4,1) = 1.05443E-14 lookup_tabl(5,4,2) = 2.96603E-12 + lookup_tabl(5,4,3) = 1.55057E-15 lookup_tabl(6,1,1) = 1.37908E-14 lookup_tabl(6,1,2) = 3.15081E-12 + lookup_tabl(6,1,3) = 2.02033E-15 lookup_tabl(6,2,1) = 1.37172E-14 lookup_tabl(6,2,2) = 3.15020E-12 + lookup_tabl(6,2,3) = 1.99948E-15 lookup_tabl(6,3,1) = 1.36362E-14 lookup_tabl(6,3,2) = 3.14927E-12 + lookup_tabl(6,3,3) = 1.97488E-15 lookup_tabl(6,4,1) = 1.35287E-14 lookup_tabl(6,4,2) = 3.14817E-12 + lookup_tabl(6,4,3) = 1.94523E-15 lookup_tabl(7,1,1) = 2.26019E-14 lookup_tabl(7,1,2) = 3.66798E-12 + lookup_tabl(7,1,3) = 3.14156E-15 lookup_tabl(7,2,1) = 2.24435E-14 lookup_tabl(7,2,2) = 3.66540E-12 + lookup_tabl(7,2,3) = 3.08114E-15 lookup_tabl(7,3,1) = 2.23254E-14 lookup_tabl(7,3,2) = 3.66173E-12 + lookup_tabl(7,3,3) = 3.01021E-15 lookup_tabl(7,4,1) = 2.20496E-14 lookup_tabl(7,4,2) = 3.65796E-12 + lookup_tabl(7,4,3) = 2.92456E-15 lookup_tabl(8,1,1) = 4.41983E-13 lookup_tabl(8,1,2) = 7.50091E-11 + lookup_tabl(8,1,3) = 1.95503E-14 lookup_tabl(8,2,1) = 3.93335E-13 lookup_tabl(8,2,2) = 6.79097E-11 + lookup_tabl(8,2,3) = 1.74308E-14 lookup_tabl(8,3,1) = 3.45569E-13 lookup_tabl(8,3,2) = 6.07845E-11 + lookup_tabl(8,3,3) = 1.53194E-14 lookup_tabl(8,4,1) = 2.96971E-13 - lookup_tabl(8,4,2) = 5.36085E-11 + lookup_tabl(8,4,2) = 5.36085E-11 + lookup_tabl(8,4,3) = 1.32479E-14 DO j=jts,jte DO k=kts,kte DO i=its,ite AOD_wfa(i,k,j) = 0. AOD_ifa(i,k,j) = 0. + if ( wif_input_opt .eq. 2 ) then + AOD_bca(i,k,j) = 0. + end if END DO END DO END DO @@ -5231,6 +5270,11 @@ SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & unit_bext3 = lookup_tabl(RH_idx1,t_idx,2) & & + (lookup_tabl(RH_idx2,t_idx,2) & & - lookup_tabl(RH_idx1,t_idx,2))*rh_f + if ( wif_input_opt .eq. 2 ) then + unit_bext4 = lookup_tabl(RH_idx1,t_idx,3) & + & + (lookup_tabl(RH_idx2,t_idx,3) & + & - lookup_tabl(RH_idx1,t_idx,3))*rh_f + end if ntemp = MAX(1., MIN(99999.E6, nwfa(i,k,j))) AOD_wfa(i,k,j) = unit_bext1*ntemp*dz8w(i,k,j)*rhoa @@ -5238,6 +5282,11 @@ SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & ntemp = MAX(0.01, MIN(9999.E6, nifa(i,k,j))) AOD_ifa(i,k,j) = unit_bext3*ntemp*dz8w(i,k,j)*rhoa + if ( wif_input_opt .eq. 2 ) then + ntemp = MAX(1., MIN(9999.E9, nbca(i,k,j))) + AOD_bca(i,k,j) = unit_bext4*ntemp*dz8w(i,k,j)*rhoa + end if + END DO END DO END DO @@ -5245,7 +5294,11 @@ SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & DO j=jts,jte DO k=kts,kte DO i=its,ite - taod5503d(i,k,j) = MAX(1.E-3, aod_wfa(i,k,j) + aod_ifa(i,k,j)) + if ( wif_input_opt .eq. 2 ) then + taod5503d(i,k,j) = MAX(1.E-3, aod_wfa(i,k,j) + aod_ifa(i,k,j) + aod_bca(i,k,j)) + else + taod5503d(i,k,j) = MAX(1.E-3, aod_wfa(i,k,j) + aod_ifa(i,k,j)) + end if END DO END DO END DO diff --git a/share/input_wrf.F b/share/input_wrf.F index ae1aff98c6..76631def2c 100644 --- a/share/input_wrf.F +++ b/share/input_wrf.F @@ -1124,7 +1124,7 @@ SUBROUTINE input_wrf ( fid , grid , config_flags , switch , ierr ) CASE ( input_only, auxinput1_only, auxinput2_only, & auxinput3_only, auxinput4_only, auxinput5_only, & auxinput6_only, auxinput7_only, auxinput8_only, & - auxinput9_only, auxinput10_only ) + auxinput9_only, auxinput10_only, auxinput17_only ) #if ( WRF_CHEM == 1 ) IF( (config_flags%io_style_emissions .eq. 1) .and. & ((switch.eq.auxinput5_only) .or. (switch.eq.auxinput6_only) .or. & diff --git a/share/mediation_integrate.F b/share/mediation_integrate.F index e37fa37e92..98bc99643a 100644 --- a/share/mediation_integrate.F +++ b/share/mediation_integrate.F @@ -214,6 +214,20 @@ SUBROUTINE med_before_solve_io ( grid , config_flags ) ENDIF ENDIF #endif + ELSE IF( ialarm .EQ. AUXINPUT17_ALARM ) THEN + IF( config_flags%qna_update .EQ. 1) THEN + IF( WRFU_AlarmIsRinging( grid%alarms( ialarm ), rc=rc ) ) THEN + call wrf_debug(00,' CALL med_read_qna_emissions ') + CALL med_read_qna_emissions ( grid , config_flags ) + CALL WRFU_AlarmRingerOff( grid%alarms( ialarm ), rc=rc ) + call wrf_debug(00,' Back from CALL med_read_qna_emissions ') + ENDIF + ELSE + IF( WRFU_AlarmIsRinging( grid%alarms( ialarm ), rc=rc ) ) THEN + CALL med_auxinput_in ( grid, ialarm, config_flags ) + CALL WRFU_AlarmRingerOff( grid%alarms( ialarm ), rc=rc ) + ENDIF + ENDIF #if ( EM_CORE == 1 ) ELSE IF( ialarm .EQ. AUXINPUT11_ALARM ) THEN IF( config_flags%obs_nudge_opt .EQ. 1) THEN @@ -2433,3 +2447,65 @@ END SUBROUTINE med_read_wrf_chem_emissopt3 #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +SUBROUTINE med_read_qna_emissions ( grid , config_flags ) + ! Driver layer + USE module_domain , ONLY : domain , domain_clock_get + USE module_io_domain + USE module_timing + USE module_configure , ONLY : grid_config_rec_type + ! Model layer + USE module_bc_time_utilities +#ifdef DM_PARALLEL + USE module_dm +#endif + USE module_date_time + USE module_utility + + IMPLICIT NONE + + ! Arguments + TYPE(domain) :: grid + + TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags + + ! Local data + LOGICAL, EXTERNAL :: wrf_dm_on_monitor + + INTEGER :: ierr, efid + REAL :: time, tupdate + real, allocatable :: dumc0(:,:,:) + CHARACTER (LEN=256) :: message, current_date_char, date_string + CHARACTER (LEN=256) :: inpname + +#include "wrf_io_flags.h" + + CALL domain_clock_get( grid, current_timestr=current_date_char ) + + CALL construct_filename1 ( inpname , 'wrfqnainp' , grid%id , 2 ) + WRITE(message,*)'mediation_integrate: med_read_qna_emissions: Open file ',TRIM(inpname) + CALL wrf_message( TRIM(message) ) + + if( grid%auxinput17_oid .NE. 0 ) then + CALL close_dataset ( grid%auxinput17_oid , config_flags , "DATASET=AUXINPUT17" ) + endif + + CALL open_r_dataset ( grid%auxinput17_oid, TRIM(inpname) , grid , config_flags, & + "DATASET=AUXINPUT17", ierr ) + IF ( ierr .NE. 0 ) THEN + WRITE( message , * ) 'med_read_qna_emissions: error opening ', TRIM( inpname ) + CALL wrf_error_fatal( TRIM( message ) ) + ENDIF + + WRITE(message,*)'mediation_integrate: med_read_qna_emissions: Read surface aerosol emissions at time ',& + TRIM(current_date_char) + CALL wrf_message( TRIM(message) ) + + CALL wrf_debug (00 , 'mediation_integrate: calling input_auxinput17' ) + CALL input_auxinput17 ( grid%auxinput17_oid, grid , config_flags , ierr ) + + CALL close_dataset ( grid%auxinput17_oid , config_flags , "DATASET=AUXINPUT17" ) + + CALL wrf_debug (00 , 'mediation_integrate: med_read_qna_emissions: exit' ) + +END SUBROUTINE med_read_qna_emissions diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index caa2feeec4..6e4b1d3146 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -1408,6 +1408,63 @@ END FUNCTION bep_bem_ngr_u END IF END IF +!----------------------------------------------------------------------- +! If qna_update = 0, set io_form_auxinput17 to 0 so WRF will not try to +! input the data; auxinput_interval must also be 0 +!----------------------------------------------------------------------- + + IF ( model_config_rec%qna_update .EQ. 0 ) THEN + model_config_rec%io_form_auxinput17 = 0 + DO i = 1, model_config_rec % max_dom + wrf_err_message = '--- NOTE: qna_update is 0, ' // & + 'setting io_form_auxinput17 = 0 and auxinput17_interval = 0 for all domains' + CALL wrf_debug ( 1, TRIM( wrf_err_message ) ) + model_config_rec%auxinput17_interval(i) = 0 + model_config_rec%auxinput17_interval_y(i) = 0 + model_config_rec%auxinput17_interval_d(i) = 0 + model_config_rec%auxinput17_interval_h(i) = 0 + model_config_rec%auxinput17_interval_m(i) = 0 + model_config_rec%auxinput17_interval_s(i) = 0 + ENDDO + ELSE + IF ( model_config_rec%io_form_auxinput17 .EQ. 0 ) THEN + wrf_err_message = '--- ERROR: If qna_update /= 0, io_form_auxinput17 must be /= 0' + CALL wrf_message ( wrf_err_message ) + wrf_err_message = '--- Set io_form_auxinput17 in the time_control namelist (probably to 2).' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + count_fatal_error = count_fatal_error + 1 + END IF + END IF + +!----------------------------------------------------------------------- +! If qna_update = 1, we need to make sure that two nml items are set: +! 1. io_form_auxinput17 = 2 (only for one domain) +! 2. auxinput17_interval = NON-ZERO (just check most coarse domain) +!----------------------------------------------------------------------- + + IF ( model_config_rec%qna_update .EQ. 1 ) THEN + IF ( model_config_rec%io_form_auxinput17 .EQ. 0 ) THEN + wrf_err_message = '--- ERROR: If qna_update /= 0, io_form_auxinput17 must be /= 0' + CALL wrf_debug ( 0, TRIM(wrf_err_message) ) + wrf_err_message = '--- Set io_form_auxinput17 in the time_control namelist (probably to 2).' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + count_fatal_error = count_fatal_error + 1 + END IF + + IF ( ( model_config_rec%auxinput17_interval(1) .EQ. 0 ) .AND. & + ( model_config_rec%auxinput17_interval_y(1) .EQ. 0 ) .AND. & + ( model_config_rec%auxinput17_interval_d(1) .EQ. 0 ) .AND. & + ( model_config_rec%auxinput17_interval_h(1) .EQ. 0 ) .AND. & + ( model_config_rec%auxinput17_interval_m(1) .EQ. 0 ) .AND. & + ( model_config_rec%auxinput17_interval_s(1) .EQ. 0 ) ) THEN + wrf_err_message = '--- ERROR: If qna_update /= 0, one of the auxinput17_interval settings must be /= 0' + CALL wrf_debug ( 0, TRIM(wrf_err_message) ) + wrf_err_message = '--- Set auxinput17_interval_s to the same value as interval_seconds (usually a pretty good guess).' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + count_fatal_error = count_fatal_error + 1 + END IF + END IF + !----------------------------------------------------------------------- ! The qndropsource relies on the flag PROGN (when not running chemistry) ! and is always allocated when running WRF Chem. @@ -2353,7 +2410,8 @@ END FUNCTION bep_bem_ngr_u 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%use_aero_icbc .AND. model_config_rec%scalar_pblmix(i) .NE. 1 ) 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 oops = oops + 1 END IF @@ -2396,9 +2454,39 @@ END FUNCTION bep_bem_ngr_u END IF END DO - IF ( model_config_rec%aer_init_opt > 0 ) THEN - model_config_rec%wif_input_opt = 1 - END IF +!----------------------------------------------------------------------- +! Set aer_fire_emit_opt for Thompson-MP-Aero (mp_physics=28) +!----------------------------------------------------------------------- + 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%aer_init_opt .EQ. 2) THEN + IF ( model_config_rec%wif_input_opt .EQ. 1 ) THEN + model_config_rec%aer_fire_emit_opt = 1 + ELSE IF ( model_config_rec%wif_input_opt .EQ. 2 ) THEN + model_config_rec%aer_fire_emit_opt = 2 + END IF + ELSE IF ( model_config_rec%aer_init_opt .EQ. 0 .OR. model_config_rec%aer_init_opt .EQ. 1) THEN + wrf_err_message = '--- ERROR: wif_fire_emit=.true. but selected aerosol source does not contain fire emissions ' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + wrf_err_message = '--- ERROR: Please use first guess aerosol source with fire emissions and set use_rap_aero_icbc=.true. ' + CALL wrf_debug ( 0, TRIM( wrf_err_message ) ) + count_fatal_error = count_fatal_error + 1 + END IF + END IF + END DO + +!----------------------------------------------------------------------- +! Set warning message if wif_fire_inj for Thompson-MP-Aero (mp_physics=28) +! 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%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 ) ) + END IF + END IF + END DO !----------------------------------------------------------------------- ! DJW Check that we're not using ndown and vertical nesting. diff --git a/share/module_optional_input.F b/share/module_optional_input.F index 13dd13b90b..5d44626d0f 100644 --- a/share/module_optional_input.F +++ b/share/module_optional_input.F @@ -7,9 +7,9 @@ MODULE module_optional_input flag_qg , flag_qh , & flag_qni , flag_qnc , flag_qnr , & flag_qns , flag_qng , flag_qnh , & - flag_qnwfa2d , flag_qnifa2d , & - flag_qnwfa , flag_qnifa , & - flag_qnwfa_cl , flag_qnifa_cl , & + flag_qnwfa2d , flag_qnifa2d , flag_qnbca2d , flag_qnocbb2d , flag_qnbcbb2d , & + flag_qnwfa , flag_qnifa , flag_qnbca , & + flag_qnwfa_cl , flag_qnifa_cl , flag_qnbca_cl , & flag_p_wif , & flag_sh , flag_speccldl , flag_speccldf @@ -305,10 +305,15 @@ SUBROUTINE optional_moist ( grid , fid , & flag_qnh = 0 flag_qnwfa2d = 0 flag_qnifa2d = 0 + flag_qnbca2d = 0 + flag_qnocbb2d = 0 + flag_qnbcbb2d = 0 flag_qnwfa = 0 flag_qnifa = 0 + flag_qnbca = 0 flag_qnwfa_cl = 0 flag_qnifa_cl = 0 + flag_qnbca_cl = 0 flag_p_wif = 0 flag_sh = 0 flag_speccldl = 0 @@ -393,6 +398,21 @@ SUBROUTINE optional_moist ( grid , fid , & IF ( ierr .EQ. 0 ) THEN flag_qnifa2d = itmp END IF + flag_name(1:8) = 'QNBCA2D ' + CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) + IF ( ierr .EQ. 0 ) THEN + flag_qnbca2d = itmp + END IF + flag_name(1:8) = 'QNOCBB2D' + CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) + IF ( ierr .EQ. 0 ) THEN + flag_qnocbb2d = itmp + END IF + flag_name(1:8) = 'QNBCBB2D' + CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) + IF ( ierr .EQ. 0 ) THEN + flag_qnbcbb2d = itmp + END IF flag_name(1:8) = 'QNWFA ' CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN @@ -403,6 +423,11 @@ SUBROUTINE optional_moist ( grid , fid , & IF ( ierr .EQ. 0 ) THEN flag_qnifa = itmp END IF + flag_name(1:8) = 'QNBCA ' + CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) + IF ( ierr .EQ. 0 ) THEN + flag_qnbca = itmp + END IF flag_name(1:8) = 'QNWFA_CL' CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN @@ -413,6 +438,11 @@ SUBROUTINE optional_moist ( grid , fid , & IF ( ierr .EQ. 0 ) THEN flag_qnifa_cl = itmp END IF + flag_name(1:8) = 'QNBCA_CL' + CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) + IF ( ierr .EQ. 0 ) THEN + flag_qnbca_cl = itmp + END IF flag_name(1:8) = 'P_WIF ' CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN