diff --git a/physics/GFS_PBL_generic_common.F90 b/physics/GFS_PBL_generic_common.F90 index 9b3f83b57..5f09f5347 100644 --- a/physics/GFS_PBL_generic_common.F90 +++ b/physics/GFS_PBL_generic_common.F90 @@ -12,7 +12,7 @@ module GFS_PBL_generic_common contains subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & - imp_physics_thompson, ltaerosol, & + imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & imp_physics_zhao_carr, imp_physics_nssl,& nssl_hail_on, nssl_ccn_on, kk, & @@ -23,7 +23,7 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & imp_physics_thompson, & imp_physics_mg, ntgl, imp_physics_gfdl, & imp_physics_zhao_carr,imp_physics_nssl - logical, intent(in ) :: ltaerosol, nssl_hail_on, nssl_ccn_on + logical, intent(in ) :: ltaerosol, mraerosol, nssl_hail_on, nssl_ccn_on integer, intent(out) :: kk character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -38,6 +38,8 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & ! Thompson if(ltaerosol) then kk = 12 + else if(mraerosol) then + kk = 10 else kk = 9 endif @@ -70,4 +72,4 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & end subroutine set_aerosol_tracer_index - end module GFS_PBL_generic_common \ No newline at end of file + end module GFS_PBL_generic_common diff --git a/physics/GFS_PBL_generic_post.F90 b/physics/GFS_PBL_generic_post.F90 index 67b5443cd..0d13dc5d8 100644 --- a/physics/GFS_PBL_generic_post.F90 +++ b/physics/GFS_PBL_generic_post.F90 @@ -12,7 +12,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev,nqrimef, & trans_aero, ntchs, ntchm, ntccn, nthl, nthnc, ntgv, nthv, & imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, imp_physics_mg, & - imp_physics_fer_hires, imp_physics_nssl, nssl_ccn_on, ltaerosol, nssl_hail_on, & + imp_physics_fer_hires, imp_physics_nssl, nssl_ccn_on, ltaerosol, mraerosol, nssl_hail_on, & cplflx, cplaqm, cplchm, lssav, flag_for_pbl_generic_tend, ldiag3d, lsidea, hybedmf, do_shoc, satmedmf, & shinhong, do_ysu, dvdftra, dusfc1, dvsfc1, dtsfc1, dqsfc1, dtf, dudt, dvdt, dtdt, htrsw, htrlw, xmu, & dqdt, dusfc_cpl, dvsfc_cpl, dtsfc_cpl, dtend, dtidx, index_of_temperature, index_of_x_wind, index_of_y_wind, & @@ -36,7 +36,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires integer, intent(in) :: imp_physics_nssl logical, intent(in) :: nssl_ccn_on, nssl_hail_on - logical, intent(in) :: ltaerosol, cplflx, cplaqm, cplchm, lssav, ldiag3d, lsidea, use_med_flux + logical, intent(in) :: ltaerosol, cplflx, cplaqm, cplchm, lssav, ldiag3d, lsidea, use_med_flux, mraerosol logical, intent(in) :: hybedmf, do_shoc, satmedmf, shinhong, do_ysu logical, intent(in) :: flag_for_pbl_generic_tend @@ -104,7 +104,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, if (trans_aero) then ! Set kk if chemistry-aerosol tracers are diffused call set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & - imp_physics_thompson, ltaerosol, & + imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & imp_physics_zhao_carr, imp_physics_nssl,& nssl_hail_on, nssl_ccn_on, kk, & @@ -165,6 +165,21 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqdt(i,k,ntia) = dvdftra(i,k,12) enddo enddo + else if(mraerosol) then + do k=1,levs + do i=1,im + dqdt(i,k,ntqv) = dvdftra(i,k,1) + dqdt(i,k,ntcw) = dvdftra(i,k,2) + dqdt(i,k,ntiw) = dvdftra(i,k,3) + dqdt(i,k,ntrw) = dvdftra(i,k,4) + dqdt(i,k,ntsw) = dvdftra(i,k,5) + dqdt(i,k,ntgl) = dvdftra(i,k,6) + dqdt(i,k,ntlnc) = dvdftra(i,k,7) + dqdt(i,k,ntinc) = dvdftra(i,k,8) + dqdt(i,k,ntrnc) = dvdftra(i,k,9) + dqdt(i,k,ntoz) = dvdftra(i,k,10) + enddo + enddo else do k=1,levs do i=1,im diff --git a/physics/GFS_PBL_generic_post.meta b/physics/GFS_PBL_generic_post.meta index 11b8ac514..b20142991 100644 --- a/physics/GFS_PBL_generic_post.meta +++ b/physics/GFS_PBL_generic_post.meta @@ -274,6 +274,13 @@ dimensions = () type = logical intent = in +[mraerosol] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [nssl_ccn_on] standard_name = nssl_ccn_on long_name = CCN activation flag in NSSL micro diff --git a/physics/GFS_PBL_generic_pre.F90 b/physics/GFS_PBL_generic_pre.F90 index 0dbdf7225..b9f7bb880 100644 --- a/physics/GFS_PBL_generic_pre.F90 +++ b/physics/GFS_PBL_generic_pre.F90 @@ -14,8 +14,8 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm, & ntccn, nthl, nthnc, ntgv, nthv, & imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & - imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, imp_physics_nssl, & - ltaerosol, nssl_ccn_on, nssl_hail_on, & + imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, imp_physics_nssl, & + ltaerosol, mraerosol, nssl_ccn_on, nssl_hail_on, & hybedmf, do_shoc, satmedmf, qgrs, vdftra, save_u, save_v, save_t, save_q, & flag_for_pbl_generic_tend, ldiag3d, qdiag3d, lssav, ugrs, vgrs, tgrs, errmsg, errflg) @@ -33,7 +33,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, logical, intent(in) :: trans_aero, ldiag3d, qdiag3d, lssav integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires - logical, intent(in) :: ltaerosol, hybedmf, do_shoc, satmedmf, flag_for_pbl_generic_tend + logical, intent(in) :: ltaerosol, hybedmf, do_shoc, satmedmf, flag_for_pbl_generic_tend, mraerosol integer, intent(in) :: imp_physics_nssl logical, intent(in) :: nssl_hail_on, nssl_ccn_on @@ -108,6 +108,22 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, enddo enddo rtg_ozone_index = 10 + elseif(mraerosol) then + do k=1,levs + do i=1,im + vdftra(i,k,1) = qgrs(i,k,ntqv) + vdftra(i,k,2) = qgrs(i,k,ntcw) + vdftra(i,k,3) = qgrs(i,k,ntiw) + vdftra(i,k,4) = qgrs(i,k,ntrw) + vdftra(i,k,5) = qgrs(i,k,ntsw) + vdftra(i,k,6) = qgrs(i,k,ntgl) + vdftra(i,k,7) = qgrs(i,k,ntlnc) + vdftra(i,k,8) = qgrs(i,k,ntinc) + vdftra(i,k,9) = qgrs(i,k,ntrnc) + vdftra(i,k,10) = qgrs(i,k,ntoz) + enddo + enddo + rtg_ozone_index = 10 else do k=1,levs do i=1,im @@ -242,7 +258,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, ! if (trans_aero) then call set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & - imp_physics_thompson, ltaerosol, & + imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & imp_physics_zhao_carr, imp_physics_nssl,& nssl_hail_on, nssl_ccn_on, kk, & @@ -297,4 +313,4 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, end subroutine GFS_PBL_generic_pre_run - end module GFS_PBL_generic_pre \ No newline at end of file + end module GFS_PBL_generic_pre diff --git a/physics/GFS_PBL_generic_pre.meta b/physics/GFS_PBL_generic_pre.meta index 5f765d508..a09b34b48 100644 --- a/physics/GFS_PBL_generic_pre.meta +++ b/physics/GFS_PBL_generic_pre.meta @@ -280,6 +280,13 @@ dimensions = () type = logical intent = in +[mraerosol] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [nssl_ccn_on] standard_name = nssl_ccn_on long_name = CCN activation flag in NSSL micro @@ -429,4 +436,4 @@ units = 1 dimensions = () type = integer - intent = out \ No newline at end of file + intent = out diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 835725eee..b442866d9 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -30,7 +30,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, & imp_physics_fer_hires, iovr_rand, iovr_maxrand, iovr_max, iovr_dcorr, & iovr_exp, iovr_exprand, idcor_con, idcor_hogan, idcor_oreopoulos, & julian, yearlen, lndp_var_list, lsswr, lslwr, & - ltaerosol, lgfdlmprad, uni_cld, effr_in, do_mynnedmf, lmfshal, & + ltaerosol, mraerosol, lgfdlmprad, uni_cld, effr_in, do_mynnedmf, lmfshal, & lmfdeep2, fhswr, fhlwr, solhr, sup, con_eps, epsm1, fvirt, & rog, rocp, con_rd, xlat_d, xlat, xlon, coslat, sinlat, tsfc, slmsk, & prsi, prsl, prslk, tgrs, sfc_wts, mg_cld, effrr_in, pert_clds, & @@ -121,8 +121,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, & character(len=3), dimension(:), intent(in) :: lndp_var_list logical, intent(in) :: lextop, lsswr, lslwr, ltaerosol, lgfdlmprad, & - uni_cld, effr_in, do_mynnedmf, & - lmfshal, lmfdeep2, pert_clds + uni_cld, effr_in, do_mynnedmf, & + lmfshal, lmfdeep2, pert_clds, mraerosol logical, intent(in) :: aero_dir_fdb real(kind=kind_phys), dimension(:,:), intent(in) :: smoke_ext, dust_ext @@ -722,7 +722,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, & enddo enddo ! for Thompson MP - prepare variables for calc_effr - if_thompson: if (imp_physics == imp_physics_thompson .and. ltaerosol) then + if_thompson: if (imp_physics == imp_physics_thompson .and. (ltaerosol .or. mraerosol)) then do k=1,LMK do i=1,IM qvs = qlyr(i,k) @@ -866,7 +866,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp, & ! Update number concentration, consistent with sub-grid clouds (GF, MYNN) or without (all others) do k=1,lm do i=1,im - if (ltaerosol .and. qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then + if ((ltaerosol .or. mraerosol) .and. qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) endif if (qi_mp(i,k)>1.e-12 .and. ni_mp(i,k)<100.) then diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 25102d078..63ab11d3e 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -557,6 +557,13 @@ dimensions = () type = logical intent = in +[mraerosol] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [lgfdlmprad] standard_name = flag_for_GFDL_microphysics_radiation_interaction long_name = flag for GFDL microphysics-radiation interaction diff --git a/physics/GFS_rrtmgp_cloud_mp.F90 b/physics/GFS_rrtmgp_cloud_mp.F90 index ca9457b4c..2acf8b4da 100644 --- a/physics/GFS_rrtmgp_cloud_mp.F90 +++ b/physics/GFS_rrtmgp_cloud_mp.F90 @@ -45,7 +45,7 @@ module GFS_rrtmgp_cloud_mp subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, & i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, i_cldliq_nc, i_cldice_nc, i_twa, kdt, & imfdeepcnv, imfdeepcnv_gf, imfdeepcnv_samf, doSWrad, doLWrad, effr_in, lmfshal, & - ltaerosol, icloud, imp_physics, imp_physics_thompson, imp_physics_gfdl, & + ltaerosol,mraerosol, icloud, imp_physics, imp_physics_thompson, imp_physics_gfdl, & lgfdlmprad, do_mynnedmf, uni_cld, lmfdeep2, p_lev, p_lay, t_lay, qs_lay, q_lay, & relhum, lsmask, xlon, xlat, dx, tv_lay, effrin_cldliq, effrin_cldice, & effrin_cldrain, effrin_cldsnow, tracer, cnv_mixratio, cld_cnv_frac, qci_conv, & @@ -86,6 +86,7 @@ subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldic effr_in, & ! Provide hydrometeor radii from macrophysics? lmfshal, & ! Flag for mass-flux shallow convection scheme used by Xu-Randall ltaerosol, & ! Flag for aerosol option + mraerosol, & ! Flag for aerosol option lgfdlmprad, & ! Flag for GFDLMP radiation interaction do_mynnedmf, & ! Flag to activate MYNN-EDMF uni_cld, & ! Flag for unified cloud scheme @@ -253,7 +254,7 @@ subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldic ! Update particle size using modified mixing-ratios from Thompson. call cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, & i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol,& - effrin_cldliq, effrin_cldice, effrin_cldsnow) + mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow) cld_reliq = effrin_cldliq cld_reice = effrin_cldice cld_resnow = effrin_cldsnow @@ -819,13 +820,13 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) !! \section cmp_reff_Thompson_gen General Algorithm subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, & i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol, & - effrin_cldliq, effrin_cldice, effrin_cldsnow) + mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow) implicit none ! Inputs integer, intent(in) :: nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, & i_cldliq_nc, i_twa - logical, intent(in) :: ltaerosol + logical, intent(in) :: ltaerosol, mraerosol real(kind_phys), intent(in) :: con_eps,con_rd real(kind_phys), dimension(:,:),intent(in) :: q_lay, p_lay, t_lay real(kind_phys), dimension(:,:,:),intent(in) :: tracer @@ -856,6 +857,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho endif + elseif (mraerosol) then + nc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq_nc) / (1.-q_lay(iCol,iLay)) + if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then + nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho + endif else nc_mp(iCol,iLay) = nt_c*orho endif diff --git a/physics/GFS_rrtmgp_cloud_mp.meta b/physics/GFS_rrtmgp_cloud_mp.meta index 88530d84c..4f48d53ef 100644 --- a/physics/GFS_rrtmgp_cloud_mp.meta +++ b/physics/GFS_rrtmgp_cloud_mp.meta @@ -189,6 +189,13 @@ dimensions = () type = logical intent = in +[mraerosol] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [imfdeepcnv] standard_name = control_for_deep_convection_scheme long_name = flag for mass-flux deep convection scheme diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 17d09a62e..b71c8f1f7 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -71,6 +71,7 @@ MODULE module_mp_thompson LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. LOGICAL, PRIVATE:: is_aerosol_aware = .false. + LOGICAL, PRIVATE:: merra2_aerosol_aware = .false. LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true. LOGICAL, PARAMETER, PRIVATE:: homogIce = .true. @@ -438,12 +439,14 @@ MODULE module_mp_thompson !>\section gen_thompson_init thompson_init General Algorithm !> @{ SUBROUTINE thompson_init(is_aerosol_aware_in, & + merra2_aerosol_aware_in, & mpicomm, mpirank, mpiroot, & threads, errmsg, errflg) IMPLICIT NONE LOGICAL, INTENT(IN) :: is_aerosol_aware_in + LOGICAL, INTENT(IN) :: merra2_aerosol_aware_in INTEGER, INTENT(IN) :: mpicomm, mpirank, mpiroot INTEGER, INTENT(IN) :: threads CHARACTER(len=*), INTENT(INOUT) :: errmsg @@ -454,14 +457,23 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & real :: stime, etime LOGICAL, PARAMETER :: precomputed_tables = .FALSE. -! Set module variable is_aerosol_aware +! Set module variable is_aerosol_aware/merra2_aerosol_aware is_aerosol_aware = is_aerosol_aware_in + merra2_aerosol_aware = merra2_aerosol_aware_in + if (is_aerosol_aware .and. merra2_aerosol_aware) then + errmsg = 'Logic error in thompson_init: only one of the two options can be true, ' // & + 'not both: is_aerosol_aware or merra2_aerosol_aware' + errflg = 1 + return + end if if (mpirank==mpiroot) then - if (is_aerosol_aware) then - write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' - else - write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' - end if + if (is_aerosol_aware) then + write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' + else if(merra2_aerosol_aware) then + write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' + else + write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' + end if end if micro_init = .FALSE. @@ -1130,15 +1142,12 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then - if (present(errmsg)) then + if (present(errmsg) .and. present(errflg)) then write(errmsg, '(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' - else - write(*,'(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' - end if - if (present(errflg)) then errflg = 1 return else + write(*,'(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' stop end if end if @@ -1148,24 +1157,32 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & .not.present(nifa) .or. & .not.present(nwfa2d) .or. & .not.present(nifa2d) )) then - if (present(errmsg)) then + if (present(errmsg) .and. present(errflg)) then write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', & ' and nifa2d for aerosol-aware version of Thompson microphysics' + errflg = 1 + return else write(*, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', & ' and nifa2d for aerosol-aware version of Thompson microphysics' + stop end if - if (present(errflg)) then + else if (merra2_aerosol_aware .and. (.not.present(nc) .or. & + .not.present(nwfa) .or. & + .not.present(nifa) )) then + if (present(errmsg) .and. present(errflg)) then + write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' errflg = 1 return else + write(*, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & + ' for merra2 aerosol-aware version of Thompson microphysics' stop end if - else if (.not.is_aerosol_aware .and. (present(nwfa) .or. & - present(nifa) .or. & - present(nwfa2d) .or. & - present(nifa2d) )) then - write(*,*) 'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware is FALSE' + else if (.not.is_aerosol_aware .and. .not.merra2_aerosol_aware .and. & + (present(nwfa) .or. present(nifa) .or. present(nwfa2d) .or. present(nifa2d))) then + write(*,*) 'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE' end if end if test_only_once @@ -1395,7 +1412,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & qcten1(k) = 0. endif initialize_extended_diagnostics enddo - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then do k = kts, kte nc1d(k) = nc(i,k,j) nwfa1d(k) = nwfa(i,k,j) @@ -2188,7 +2205,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) & / am_r*lamc**bm_r) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c else qc1d(k) = 0.0 nc1d(k) = 0.0 @@ -2854,7 +2871,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !! Implemented by T. Eidhammer and G. Thompson 2012Dec18 !+---+-----------------------------------------------------------------+ - if (dustyIce .AND. is_aerosol_aware) then + if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then xni = iceDeMott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k)) else xni = 1.0 *1000. ! Default is 1.0 per Liter @@ -2902,7 +2919,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !! we may need to relax the temperature and ssati constraints. if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps & .and. temp(k).lt.253.15) ) then - if (dustyIce .AND. is_aerosol_aware) then + if (dustyIce .AND. (is_aerosol_aware .or. merra2_aerosol_aware)) then xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k)) xnc = xnc*(1.0 + 50.*rand3) else @@ -2916,7 +2933,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Freezing of aqueous aerosols based on Koop et al (2001, Nature) xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave - if (is_aerosol_aware .AND. homogIce .AND. (xni.le.499.E3) & + if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.499.E3) & & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave) pni_iha(k) = xnc*odts @@ -3370,7 +3387,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & if ((qc1d(k) + qcten(k)*DT) .gt. R1) then rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k) nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c L_qc(k) = .true. else rc(k) = R1 @@ -3538,7 +3555,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & prw_vcd(k) = clap*odt !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION if (clap .gt. eps) then - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k))) else xnc = Nt_c @@ -3546,7 +3563,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho !+---+-----------------------------------------------------------------+ ! EVAPORATION - elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.E-6 .AND. is_aerosol_aware) then + elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.E-6 .AND. & + (is_aerosol_aware .or. merra2_aerosol_aware)) then tempc = temp(k) - 273.15 otemp = 1./temp(k) rvs = rho(k)*qvs(k) @@ -3610,7 +3628,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k)) if (rc(k).eq.R1) L_qc(k) = .false. nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k)) temp(k) = t1d(k) + DT*tten(k) rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) @@ -5760,7 +5778,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) - if (.NOT. is_aerosol_aware) nc(k) = Nt_c + if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware)) nc(k) = Nt_c if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. ri(k) = MAX(R1, qi1d(k)*rho(k)) ni(k) = MAX(R2, ni1d(k)*rho(k)) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 99ab601f9..80762a5f1 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -30,14 +30,15 @@ module mp_thompson !! \section arg_table_mp_thompson_init Argument Table !! \htmlinclude mp_thompson_init.html !! - subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & - restart, imp_physics, & - imp_physics_thompson, convert_dry_rho,& - spechum, qc, qr, qi, qs, qg, ni, nr, & - is_aerosol_aware, nc, nwfa2d, nifa2d, & - nwfa, nifa, tgrs, prsl, phil, area, & - mpicomm, mpirank, mpiroot, & - threads, ext_diag, diag3d, & + subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & + restart, imp_physics, & + imp_physics_thompson, convert_dry_rho, & + spechum, qc, qr, qi, qs, qg, ni, nr, & + is_aerosol_aware, merra2_aerosol_aware, & + nc, nwfa2d, nifa2d, & + nwfa, nifa, tgrs, prsl, phil, area, & + aerfld, mpicomm, mpirank, mpiroot, & + threads, ext_diag, diag3d, & errmsg, errflg) implicit none @@ -61,11 +62,13 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & real(kind_phys), intent(inout) :: nr(:,:) ! Aerosols logical, intent(in ) :: is_aerosol_aware + logical, intent(in ) :: merra2_aerosol_aware real(kind_phys), intent(inout) :: nc(:,:) real(kind_phys), intent(inout) :: nwfa(:,:) real(kind_phys), intent(inout) :: nifa(:,:) real(kind_phys), intent(inout) :: nwfa2d(:) real(kind_phys), intent(inout) :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) ! State variables real(kind_phys), intent(in ) :: tgrs(:,:) real(kind_phys), intent(in ) :: prsl(:,:) @@ -115,10 +118,17 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & end if end if + if (is_aerosol_aware .and. merra2_aerosol_aware) then + write(errmsg,'(*(a))') "Logic error: Only one Thompson aerosol option can be true, either is_aerosol_aware or merra2_aerosol_aware)" + errflg = 1 + return + end if + ! Call Thompson init - call thompson_init(is_aerosol_aware_in=is_aerosol_aware, mpicomm=mpicomm, & - mpirank=mpirank, mpiroot=mpiroot, threads=threads, & - errmsg=errmsg, errflg=errflg) + call thompson_init(is_aerosol_aware_in=is_aerosol_aware, & + merra2_aerosol_aware_in=merra2_aerosol_aware, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) if (errflg /= 0) return ! For restart runs, the init is done here @@ -198,6 +208,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) enddo enddo + else if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) else if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosols are present.' if (MAXVAL(nwfa2d) .lt. eps) then @@ -260,6 +272,13 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ! Copy to local array for calculating cloud effective radii below nc_local = nc + + else if (merra2_aerosol_aware) then + + ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. + where(qc .LE. 0.0) nc=0.0 + where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho + where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 else @@ -278,7 +297,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ni = ni/(1.0_kind_phys+qv) nr = nr/(1.0_kind_phys+qv) - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then nc = nc/(1.0_kind_phys+qv) nwfa = nwfa/(1.0_kind_phys+qv) nifa = nifa/(1.0_kind_phys+qv) @@ -299,13 +318,15 @@ end subroutine mp_thompson_init subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & con_eps, convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & - is_aerosol_aware, nc, nwfa, nifa, & + is_aerosol_aware, & + merra2_aerosol_aware, nc, nwfa, nifa,& nwfa2d, nifa2d, aero_ind_fdb, & tgrs, prsl, phii, omega, & sedi_semi, decfl, dtp, dt_inner, & first_time_step, istep, nsteps, & prcp, rain, graupel, ice, snow, sr, & refl_10cm, reset_dBZ, do_radar_ref, & + aerfld, & mpicomm, mpirank, mpiroot, blkno, & ext_diag, diag3d, reset_diag3d, & spp_wts_mp, spp_mp, n_var_spp, & @@ -336,12 +357,13 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & real(kind_phys), intent(inout) :: nr(:,:) ! Aerosols logical, intent(in) :: is_aerosol_aware, reset_dBZ - ! The following arrays are not allocated if is_aerosol_aware is false + logical, intent(in) :: merra2_aerosol_aware real(kind_phys), optional, intent(inout) :: nc(:,:) real(kind_phys), optional, intent(inout) :: nwfa(:,:) real(kind_phys), optional, intent(inout) :: nifa(:,:) real(kind_phys), optional, intent(in ) :: nwfa2d(:) real(kind_phys), optional, intent(in ) :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) logical, optional, intent(in ) :: aero_ind_fdb ! State variables and timestep information real(kind_phys), intent(inout) :: tgrs(:,:) @@ -494,6 +516,14 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ' nc, nwfa, nifa, nwfa2d, nifa2d' errflg = 1 return + else if (merra2_aerosol_aware .and. .not. (present(nc) .and. & + present(nwfa) .and. & + present(nifa) )) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & + ' merra2 aerosol-aware microphysics require the', & + ' following optional arguments: nc, nwfa, nifa' + errflg = 1 + return end if ! Consistency cheecks - subcycling and inner loop at the same time are not supported if (nsteps>1 .and. dt_inner < dtp) then @@ -649,9 +679,11 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ncten3 => diag3d(:,:,36:36) qcten3 => diag3d(:,:,37:37) end if set_extended_diagnostic_pointers - - !> - Call mp_gt_driver() with or without aerosols - if (is_aerosol_aware) then + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + !> - Call mp_gt_driver() with or without aerosols, with or without effective radii, ... + if (is_aerosol_aware .or. merra2_aerosol_aware) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & aero_ind_fdb=aero_ind_fdb, & @@ -750,7 +782,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ni = ni/(1.0_kind_phys+qv) nr = nr/(1.0_kind_phys+qv) - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then nc = nc/(1.0_kind_phys+qv) nwfa = nwfa/(1.0_kind_phys+qv) nifa = nifa/(1.0_kind_phys+qv) @@ -846,4 +878,26 @@ subroutine mp_thompson_finalize(errmsg, errflg) end subroutine mp_thompson_finalize + subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + ! To calculate nifa and nwfa from bins of aerosols. + ! In GOCART and MERRA2, aerosols are given as mixing ratio (kg/kg). To + ! convert from kg/kg to #/kg, the "unit mass" (mass of one particle) + ! within the mass bins is calculated. A lognormal size distribution + ! within aerosol bins is used to find the size based upon the median + ! mass. NIFA is mainly summarized over five dust bins and NWFA over the + ! other 10 bins. The parameters besides each bins are carefully tuned + ! for a good performance of the scheme. + implicit none + integer, intent(in)::ncol, nlev + real (kind=kind_phys), dimension(:,:,:), intent(in) :: aerfld + real (kind=kind_phys), dimension(:,:), intent(out ):: nifa, nwfa + + nifa=(aerfld(:,:,1)/4.0737762+aerfld(:,:,2)/30.459203+aerfld(:,:,3)/153.45048+ & + aerfld(:,:,4)/1011.5142+ aerfld(:,:,5)/5683.3501)*1.e15 + + nwfa=((aerfld(:,:,6)/0.0045435214+aerfld(:,:,7)/0.2907854+aerfld(:,:,8)/12.91224+ & + aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*1.+aerfld(:,:,11)/0.3053104*5+ & + aerfld(:,:,15)/0.3232698*1)*1.e15 + end subroutine get_niwfa + end module mp_thompson diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index e376ef233..93a3ae7de 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -144,6 +144,13 @@ dimensions = () type = logical intent = in +[merra2_aerosol_aware] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [nc] standard_name = mass_number_concentration_of_cloud_liquid_water_particles_in_air long_name = cloud droplet number concentration @@ -216,6 +223,14 @@ type = real kind = kind_phys intent = in +[aerfld] + standard_name = mass_mixing_ratio_of_aerosol_from_gocart_or_merra2 + long_name = mass mixing ratio of aerosol from gocart or merra2 + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_layer_dimension,number_of_aerosol_tracers_MG) + type = real + kind = kind_phys + intent = in [mpicomm] standard_name = mpi_communicator long_name = MPI communicator @@ -395,6 +410,13 @@ dimensions = () type = logical intent = in +[merra2_aerosol_aware] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [nc] standard_name = mass_number_concentration_of_cloud_liquid_water_particles_in_air_of_new_state long_name = cloud droplet number concentration @@ -595,6 +617,14 @@ dimensions = () type = logical intent = in +[aerfld] + standard_name = mass_mixing_ratio_of_aerosol_from_gocart_or_merra2 + long_name = mass mixing ratio of aerosol from gocart or merra2 + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_aerosol_tracers_MG) + type = real + kind = kind_phys + intent = in [mpicomm] standard_name = mpi_communicator long_name = MPI communicator diff --git a/physics/mynnedmf_wrapper.F90 b/physics/mynnedmf_wrapper.F90 index a7bbdf5a1..9ad6f4aa3 100644 --- a/physics/mynnedmf_wrapper.F90 +++ b/physics/mynnedmf_wrapper.F90 @@ -160,7 +160,7 @@ SUBROUTINE mynnedmf_wrapper_run( & & imp_physics_thompson, imp_physics_wsm6, & & chem3d, frp, mix_chem, rrfs_smoke, fire_turb, nchem, ndvel, & & imp_physics_nssl, nssl_ccn_on, & - & ltaerosol, spp_wts_pbl, spp_pbl, lprnt, huge, errmsg, errflg ) + & ltaerosol, mraerosol, spp_wts_pbl, spp_pbl, lprnt, huge, errmsg, errflg ) ! should be moved to inside the mynn: use machine, only: kind_phys @@ -187,7 +187,7 @@ SUBROUTINE mynnedmf_wrapper_run( & logical, intent(in) :: & & bl_mynn_tkeadvect, & & bl_mynn_tkebudget, & - & ltaerosol, & + & ltaerosol, mraerosol, & & lprnt, & & do_mynnsfclay, & & flag_for_pbl_generic_tend, & @@ -475,6 +475,32 @@ SUBROUTINE mynnedmf_wrapper_run( & qnifa(i,k) = qgrs_ice_aer_num_conc(i,k) enddo enddo + else if(mraerosol) then + FLAG_QI = .true. + FLAG_QNI= .true. + FLAG_QC = .true. + FLAG_QNC= .true. + FLAG_QNWFA= .false. + FLAG_QNIFA= .false. + p_qc = 2 + p_qr = 0 + p_qi = 2 + p_qs = 0 + p_qg = 0 + p_qnc= 0 + p_qni= 0 + do k=1,levs + do i=1,im + sqv(i,k) = qgrs_water_vapor(i,k) + sqc(i,k) = qgrs_liquid_cloud(i,k) + sqi(i,k) = qgrs_ice_cloud(i,k) + qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) + qni(i,k) = qgrs_cloud_ice_num_conc(i,k) + ozone(i,k) = qgrs_ozone(i,k) + qnwfa(i,k) = 0. + qnifa(i,k) = 0. + enddo + enddo else FLAG_QI = .true. FLAG_QNI= .true. @@ -860,6 +886,23 @@ SUBROUTINE mynnedmf_wrapper_run( & ! !qgrs_ice_aer_num_conc(i,k) = qgrs_ice_aer_num_conc(i,k) + RQNIFABLTEN(i,k)*delt ! enddo !enddo + else if(mraerosol) then + do k=1,levs + do i=1,im + dqdt_water_vapor(i,k) = RQVBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_liquid_cloud(i,k) = RQCBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) + dqdt_ice_cloud(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) + enddo + enddo + if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then + call dtend_helper(100+ntqv,RQVBLTEN) + call dtend_helper(100+ntcw,RQCBLTEN) + call dtend_helper(100+ntlnc,RQNCBLTEN) + call dtend_helper(100+ntiw,RQIBLTEN) + call dtend_helper(100+ntinc,RQNIBLTEN) + endif else !Thompson (2008) do k=1,levs diff --git a/physics/mynnedmf_wrapper.meta b/physics/mynnedmf_wrapper.meta index 33f97113f..a44a13f1b 100644 --- a/physics/mynnedmf_wrapper.meta +++ b/physics/mynnedmf_wrapper.meta @@ -1401,6 +1401,13 @@ dimensions = () type = logical intent = in +[mraerosol] + standard_name = do_merra2_aerosol_awareness + long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics + units = flag + dimensions = () + type = logical + intent = in [spp_wts_pbl] standard_name = spp_weights_for_pbl_scheme long_name = spp weights for pbl scheme