!=============================================================================== ! Modal Aerosol Model !=============================================================================== module aero_model use shr_kind_mod, only: r8 => shr_kind_r8 use constituents, only: pcnst, cnst_name, cnst_get_ind use ppgrid, only: pcols, pver, pverp use phys_control, only: phys_getopts, cam_physpkg_is use cam_abortutils, only: endrun use cam_logfile, only: iulog use perf_mod, only: t_startf, t_stopf use camsrfexch, only: cam_in_t, cam_out_t use aerodep_flx, only: aerodep_flx_prescribed use physics_types, only: physics_state, physics_ptend, physics_ptend_init use physics_buffer, only: physics_buffer_desc use physics_buffer, only: pbuf_get_field, pbuf_get_index, pbuf_set_field use physconst, only: gravit, rair, rhoh2o use spmd_utils, only: masterproc use infnan, only: nan, assignment(=) use cam_history, only: outfld, fieldname_len use chem_mods, only: gas_pcnst, adv_mass use mo_tracname, only: solsym use modal_aero_data,only: cnst_name_cw, lptr_so4_cw_amode use modal_aero_data,only: ntot_amode, modename_amode, nspec_max use ref_pres, only: top_lev => clim_modal_aero_top_lev use modal_aero_wateruptake, only: modal_strat_sulfate use mo_setsox, only: setsox, has_sox implicit none private public :: aero_model_readnl public :: aero_model_register public :: aero_model_init public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols. public :: aero_model_drydep ! aerosol dry deposition and sediment public :: aero_model_wetdep ! aerosol wet removal public :: aero_model_emissions ! aerosol emissions public :: aero_model_surfarea ! tropopspheric aerosol wet surface area for chemistry public :: aero_model_strat_surfarea ! stratospheric aerosol wet surface area for chemistry public :: calc_1_impact_rate public :: nimptblgrow_mind, nimptblgrow_maxd ! Accessor functions public :: get_scavimptblvol, get_scavimptblnum, get_dlndg_nimptblgrow ! Misc private data ! number of modes integer :: nmodes integer :: pblh_idx = 0 integer :: dgnum_idx = 0 integer :: dgnumwet_idx = 0 integer :: rate1_cw2pr_st_idx = 0 integer :: wetdens_ap_idx = 0 integer :: qaerwat_idx = 0 integer :: fracis_idx = 0 integer :: prain_idx = 0 integer :: rprddp_idx = 0 integer :: rprdsh_idx = 0 integer :: nevapr_shcu_idx = 0 integer :: nevapr_dpcu_idx = 0 integer :: sulfeq_idx = -1 ! variables for table lookup of aerosol impaction/interception scavenging rates integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12 real(r8) :: dlndg_nimptblgrow real(r8),allocatable :: scavimptblnum(:,:) real(r8),allocatable :: scavimptblvol(:,:) ! for surf_area_dens integer,allocatable :: num_idx(:) integer,allocatable :: index_tot_mass(:,:) integer,allocatable :: index_chm_mass(:,:) integer :: ndx_h2so4 character(len=fieldname_len), allocatable :: dgnum_name(:), dgnumwet_name(:) ! Namelist variables character(len=16) :: wetdep_list(pcnst) = ' ' character(len=16) :: drydep_list(pcnst) = ' ' real(r8) :: sol_facti_cloud_borne = 1._r8 real(r8) :: sol_factb_interstitial = 0.1_r8 real(r8) :: sol_factic_interstitial = 0.4_r8 real(r8) :: seasalt_emis_scale integer :: ndrydep = 0 integer,allocatable :: drydep_indices(:) integer :: nwetdep = 0 integer,allocatable :: wetdep_indices(:) logical :: drydep_lq(pcnst) logical :: wetdep_lq(pcnst) logical :: modal_accum_coarse_exch = .false. logical :: convproc_do_aer contains !============================================================================= ! reads aerosol namelist options !============================================================================= subroutine aero_model_readnl(nlfile) use namelist_utils, only: find_group_name use units, only: getunit, freeunit use mpishorthand character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input ! Local variables integer :: unitn, ierr character(len=*), parameter :: subname = 'aero_model_readnl' ! Namelist variables character(len=16) :: aer_wetdep_list(pcnst) = ' ' character(len=16) :: aer_drydep_list(pcnst) = ' ' namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list, sol_facti_cloud_borne, & sol_factb_interstitial, sol_factic_interstitial, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale !----------------------------------------------------------------------------- ! Read namelist if (masterproc) then unitn = getunit() open( unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'aerosol_nl', status=ierr) if (ierr == 0) then read(unitn, aerosol_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if end if close(unitn) call freeunit(unitn) end if #ifdef SPMD ! Broadcast namelist variables call mpibcast(aer_wetdep_list, len(aer_wetdep_list(1))*pcnst, mpichar, 0, mpicom) call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom) call mpibcast(sol_facti_cloud_borne, 1, mpir8, 0, mpicom) call mpibcast(sol_factb_interstitial, 1, mpir8, 0, mpicom) call mpibcast(sol_factic_interstitial, 1, mpir8, 0, mpicom) call mpibcast(modal_strat_sulfate, 1, mpilog, 0, mpicom) call mpibcast(seasalt_emis_scale, 1, mpir8, 0, mpicom) call mpibcast(modal_accum_coarse_exch, 1, mpilog, 0, mpicom) #endif wetdep_list = aer_wetdep_list drydep_list = aer_drydep_list end subroutine aero_model_readnl !============================================================================= !============================================================================= subroutine aero_model_register() use modal_aero_data, only: modal_aero_data_reg call modal_aero_data_reg() end subroutine aero_model_register !============================================================================= !============================================================================= subroutine aero_model_init( pbuf2d ) use mo_chem_utls, only: get_inv_ndx use cam_history, only: addfld, add_default, horiz_only use mo_chem_utls, only: get_rxt_ndx, get_spc_ndx use modal_aero_data, only: cnst_name_cw use modal_aero_data, only: modal_aero_data_init use rad_constituents,only: rad_cnst_get_info use dust_model, only: dust_init, dust_names, dust_active, dust_nbin, dust_nnum use seasalt_model, only: seasalt_init, seasalt_names, seasalt_active,seasalt_nbin use drydep_mod, only: inidrydep use wetdep, only: wetdep_init use modal_aero_calcsize, only: modal_aero_calcsize_init use modal_aero_coag, only: modal_aero_coag_init use modal_aero_deposition, only: modal_aero_deposition_init use modal_aero_gasaerexch, only: modal_aero_gasaerexch_init use modal_aero_newnuc, only: modal_aero_newnuc_init use modal_aero_rename, only: modal_aero_rename_init use modal_aero_convproc, only: ma_convproc_init ! args type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! local vars character(len=*), parameter :: subrname = 'aero_model_init' integer :: m, n, id character(len=20) :: dummy logical :: history_aerosol ! Output MAM or SECT aerosol tendencies integer :: l character(len=6) :: test_name character(len=64) :: errmes character(len=2) :: unit_basename ! Units 'kg' or '1' integer :: errcode character(len=fieldname_len) :: field_name character(len=32) :: spec_name character(len=32) :: spec_type character(len=32) :: mode_type integer :: nspec dgnum_idx = pbuf_get_index('DGNUM') dgnumwet_idx = pbuf_get_index('DGNUMWET') fracis_idx = pbuf_get_index('FRACIS') prain_idx = pbuf_get_index('PRAIN') rprddp_idx = pbuf_get_index('RPRDDP') rprdsh_idx = pbuf_get_index('RPRDSH') nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU') nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU') sulfeq_idx = pbuf_get_index('MAMH2SO4EQ',errcode) call phys_getopts(history_aerosol_out = history_aerosol, & convproc_do_aer_out = convproc_do_aer) call rad_cnst_get_info(0, nmodes=nmodes) call modal_aero_data_init(pbuf2d) call modal_aero_bcscavcoef_init() call modal_aero_rename_init( modal_accum_coarse_exch ) ! calcsize call must follow rename call call modal_aero_calcsize_init( pbuf2d ) call modal_aero_gasaerexch_init ! coag call must follow gasaerexch call call modal_aero_coag_init call modal_aero_newnuc_init ! call modal_aero_deposition_init only if the user has not specified ! prescribed aerosol deposition fluxes if (.not.aerodep_flx_prescribed()) then call modal_aero_deposition_init endif if (convproc_do_aer) then call ma_convproc_init() endif call dust_init() call seasalt_init(seasalt_emis_scale) call wetdep_init() nwetdep = 0 ndrydep = 0 count_species: do m = 1,pcnst if ( len_trim(wetdep_list(m)) /= 0 ) then nwetdep = nwetdep+1 endif if ( len_trim(drydep_list(m)) /= 0 ) then ndrydep = ndrydep+1 endif enddo count_species if (nwetdep>0) & allocate(wetdep_indices(nwetdep)) if (ndrydep>0) & allocate(drydep_indices(ndrydep)) do m = 1,ndrydep call cnst_get_ind ( drydep_list(m), id, abort=.false. ) if (id>0) then drydep_indices(m) = id else call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) ) endif if (masterproc) then write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied' endif enddo do m = 1,nwetdep call cnst_get_ind ( wetdep_list(m), id, abort=.false. ) if (id>0) then wetdep_indices(m) = id else call endrun(subrname//': invalid wetdep species: '//trim(wetdep_list(m)) ) endif if (masterproc) then write(iulog,*) subrname//': '//wetdep_list(m)//' will have wet removal' endif enddo if (ndrydep>0) then call inidrydep(rair, gravit) dummy = 'RAM1' call addfld (dummy,horiz_only, 'A','frac','RAM1') if ( history_aerosol ) then call add_default (dummy, 1, ' ') endif dummy = 'airFV' call addfld (dummy,horiz_only, 'A','frac','FV') if ( history_aerosol ) then call add_default (dummy, 1, ' ') endif endif if (dust_active) then ! emissions diagnostics .... do m = 1, dust_nbin+dust_nnum dummy = trim(dust_names(m)) // 'SF' call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission') if (history_aerosol) then call add_default (dummy, 1, ' ') endif enddo dummy = 'DSTSFMBL' call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface') if (history_aerosol) then call add_default (dummy, 1, ' ') endif dummy = 'LND_MBL' call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor') if (history_aerosol) then call add_default (dummy, 1, ' ') endif endif if (seasalt_active) then dummy = 'SSTSFMBL' call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface') if (history_aerosol) then call add_default (dummy, 1, ' ') endif do m = 1, seasalt_nbin dummy = trim(seasalt_names(m)) // 'SF' call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission') if (history_aerosol) then call add_default (dummy, 1, ' ') endif enddo endif ! set flags for drydep tendencies drydep_lq(:) = .false. do m=1,ndrydep id = drydep_indices(m) drydep_lq(id) = .true. enddo ! set flags for wetdep tendencies wetdep_lq(:) = .false. do m=1,nwetdep id = wetdep_indices(m) wetdep_lq(id) = .true. enddo wetdens_ap_idx = pbuf_get_index('WETDENS_AP') qaerwat_idx = pbuf_get_index('QAERWAT') pblh_idx = pbuf_get_index('pblh') rate1_cw2pr_st_idx = pbuf_get_index('RATE1_CW2PR_ST') call pbuf_set_field(pbuf2d, rate1_cw2pr_st_idx, 0.0_r8) do m = 1,ndrydep ! units if (drydep_list(m)(1:3) == 'num') then unit_basename = ' 1' else unit_basename = 'kg' endif call addfld (trim(drydep_list(m))//'DDF', horiz_only, 'A',unit_basename//'/m2/s ', & trim(drydep_list(m))//' dry deposition flux at bottom (grav + turb)') call addfld (trim(drydep_list(m))//'TBF', horiz_only, 'A',unit_basename//'/m2/s', & trim(drydep_list(m))//' turbulent dry deposition flux') call addfld (trim(drydep_list(m))//'GVF', horiz_only, 'A',unit_basename//'/m2/s ', & trim(drydep_list(m))//' gravitational dry deposition flux') call addfld (trim(drydep_list(m))//'DTQ', (/ 'lev' /), 'A',unit_basename//'/kg/s ', & trim(drydep_list(m))//' dry deposition') call addfld (trim(drydep_list(m))//'DDV', (/ 'lev' /), 'A','m/s', & trim(drydep_list(m))//' deposition velocity') if ( history_aerosol ) then call add_default (trim(drydep_list(m))//'DDF', 1, ' ') call add_default (trim(drydep_list(m))//'TBF', 1, ' ') call add_default (trim(drydep_list(m))//'GVF', 1, ' ') endif enddo do m = 1,nwetdep ! units if (wetdep_list(m)(1:3) == 'num') then unit_basename = ' 1' else unit_basename = 'kg' endif call addfld (trim(wetdep_list(m))//'SFWET', & horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux at surface') call addfld (trim(wetdep_list(m))//'SFSIC', & horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (incloud, convective) at surface') call addfld (trim(wetdep_list(m))//'SFSIS', & horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (incloud, stratiform) at surface') call addfld (trim(wetdep_list(m))//'SFSBC', & horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (belowcloud, convective) at surface') call addfld (trim(wetdep_list(m))//'SFSBS', & horiz_only, 'A',unit_basename//'/m2/s ','Wet deposition flux (belowcloud, stratiform) at surface') if (convproc_do_aer) then call addfld (trim(wetdep_list(m))//'SFSES', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, stratiform) at surface') call addfld (trim(wetdep_list(m))//'SFSBD', & horiz_only, 'A','kg/m2/s','Wet deposition flux (belowcloud, deep convective) at surface') end if call addfld (trim(wetdep_list(m))//'WET',(/ 'lev' /), 'A',unit_basename//'/kg/s ','wet deposition tendency') call addfld (trim(wetdep_list(m))//'SIC',(/ 'lev' /), 'A',unit_basename//'/kg/s ', & trim(wetdep_list(m))//' ic wet deposition') call addfld (trim(wetdep_list(m))//'SIS',(/ 'lev' /), 'A',unit_basename//'/kg/s ', & trim(wetdep_list(m))//' is wet deposition') call addfld (trim(wetdep_list(m))//'SBC',(/ 'lev' /), 'A',unit_basename//'/kg/s ', & trim(wetdep_list(m))//' bc wet deposition') call addfld (trim(wetdep_list(m))//'SBS',(/ 'lev' /), 'A',unit_basename//'/kg/s ', & trim(wetdep_list(m))//' bs wet deposition') if ( history_aerosol ) then call add_default (trim(wetdep_list(m))//'SFWET', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSIC', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSIS', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSBC', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSBS', 1, ' ') if (convproc_do_aer) then call add_default (trim(wetdep_list(m))//'SFSES', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSBD', 1, ' ') end if endif enddo do m = 1,gas_pcnst if ( solsym(m)(1:3) == 'num') then unit_basename = ' 1' ! Units 'kg' or '1' else unit_basename = 'kg' ! Units 'kg' or '1' end if call addfld( 'GS_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', & trim(solsym(m))//' gas chemistry/wet removal (for gas species)') call addfld( 'AQ_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', & trim(solsym(m))//' aqueous chemistry (for gas species)') if ( history_aerosol ) then call add_default( 'GS_'//trim(solsym(m)), 1, ' ') call add_default( 'AQ_'//trim(solsym(m)), 1, ' ') endif enddo do n = 1,pcnst if( .not. (cnst_name_cw(n) == ' ') ) then if (cnst_name_cw(n)(1:3) == 'num') then unit_basename = ' 1' else unit_basename = 'kg' endif call addfld( cnst_name_cw(n), (/ 'lev' /), 'A', unit_basename//'/kg ', & trim(cnst_name_cw(n))//' in cloud water') call addfld (trim(cnst_name_cw(n))//'SFWET', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' wet deposition flux at surface') call addfld (trim(cnst_name_cw(n))//'SFSIC', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' wet deposition flux (incloud, convective) at surface') call addfld (trim(cnst_name_cw(n))//'SFSIS', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' wet deposition flux (incloud, stratiform) at surface') call addfld (trim(cnst_name_cw(n))//'SFSBC', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' wet deposition flux (belowcloud, convective) at surface') call addfld (trim(cnst_name_cw(n))//'SFSBS', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' wet deposition flux (belowcloud, stratiform) at surface') call addfld (trim(cnst_name_cw(n))//'DDF', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' dry deposition flux at bottom (grav + turb)') call addfld (trim(cnst_name_cw(n))//'TBF', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' turbulent dry deposition flux') call addfld (trim(cnst_name_cw(n))//'GVF', horiz_only, 'A', unit_basename//'/m2/s ', & trim(cnst_name_cw(n))//' gravitational dry deposition flux') if (convproc_do_aer) then call addfld (trim(cnst_name_cw(n))//'SFSEC', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, convective) at surface') call addfld (trim(cnst_name_cw(n))//'SFSES', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, stratiform) at surface') call addfld (trim(cnst_name_cw(n))//'SFSBD', & horiz_only, 'A','kg/m2/s','Wet deposition flux (belowcloud, deep convective) at surface') end if if ( history_aerosol ) then call add_default( cnst_name_cw(n), 1, ' ' ) call add_default (trim(cnst_name_cw(n))//'GVF', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFWET', 1, ' ') call add_default (trim(cnst_name_cw(n))//'TBF', 1, ' ') call add_default (trim(cnst_name_cw(n))//'DDF', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSBS', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSIC', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSBC', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSIS', 1, ' ') if (convproc_do_aer) then call add_default (trim(cnst_name_cw(n))//'SFSEC', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSES', 1, ' ') call add_default (trim(cnst_name_cw(n))//'SFSBD', 1, ' ') end if endif endif enddo allocate(dgnum_name(ntot_amode), dgnumwet_name(ntot_amode)) do n=1,ntot_amode dgnum_name(n) = ' ' dgnumwet_name(n) = ' ' write(dgnum_name(n),fmt='(a,i1)') 'dgnum',n write(dgnumwet_name(n),fmt='(a,i1)') 'dgnumwet',n call addfld( dgnum_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode dry diameter' ) call addfld( dgnumwet_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode wet diameter' ) if ( history_aerosol ) then call add_default( dgnum_name(n), 1, ' ' ) call add_default( dgnumwet_name(n), 1, ' ' ) endif if (modal_strat_sulfate) then field_name = ' ' write(field_name,fmt='(a,i1)') 'wtpct_a',n call addfld( field_name, (/ 'lev' /), 'I', '%', 'Aerosol mode weight percent H2SO4' ) if ( history_aerosol ) then call add_default (field_name, 0, 'I') endif field_name = ' ' write(field_name,fmt='(a,i1)') 'sulfeq_a',n call addfld( field_name, (/ 'lev' /), 'I', 'kg/kg', 'H2SO4 equilibrium mixing ratio' ) if ( history_aerosol ) then call add_default (field_name, 0, 'I') endif field_name = ' ' write(field_name,fmt='(a,i1)') 'sulden_a',n call addfld( field_name, (/ 'lev' /), 'I', 'g/cm3', 'Sulfate aerosol particle mass density' ) if ( history_aerosol ) then call add_default (field_name, 0, 'I') endif end if end do ndx_h2so4 = get_spc_ndx('H2SO4') allocate(num_idx(ntot_amode)) num_idx = -1 ! for aero_model_surfarea called from mo_usrrxt do l=1,ntot_amode test_name = ' ' write(test_name,fmt='(a5,i1)') 'num_a',l num_idx(l) = get_spc_ndx( trim(test_name) ) if (num_idx(l) < 0) then write(errmes,fmt='(a,i1)') 'usrrxt_inti: cannot find MAM num_idx ',l write(iulog,*) errmes call endrun(errmes) endif end do allocate(index_tot_mass(nmodes,nspec_max)) allocate(index_chm_mass(nmodes,nspec_max)) index_tot_mass = -1 index_chm_mass = -1 ! for surf_area_dens ! define indeces associated with the various aerosol types do n = 1,nmodes call rad_cnst_get_info(0, n, mode_type=mode_type, nspec=nspec) if ( trim(mode_type) /= 'primary_carbon') then ! ignore the primary_carbon mode do l = 1, nspec call rad_cnst_get_info(0, n, l, spec_type=spec_type, spec_name=spec_name) index_tot_mass(n,l) = get_spc_ndx(spec_name) if ( trim(spec_type) == 'sulfate' .or. & trim(spec_type) == 's-organic' .or. & trim(spec_type) == 'black-c' .or. & trim(spec_type) == 'ammonium') then index_chm_mass(n,l) = get_spc_ndx(spec_name) endif enddo endif enddo if (has_sox) then do m = 1, ntot_amode l = lptr_so4_cw_amode(m) if (l > 0) then call addfld (& trim(cnst_name_cw(l))//'AQSO4',horiz_only, 'A','kg/m2/s', & trim(cnst_name_cw(l))//' aqueous phase chemistry') call addfld (& trim(cnst_name_cw(l))//'AQH2SO4',horiz_only, 'A','kg/m2/s', & trim(cnst_name_cw(l))//' aqueous phase chemistry') if ( history_aerosol ) then call add_default (trim(cnst_name_cw(l))//'AQSO4', 1, ' ') call add_default (trim(cnst_name_cw(l))//'AQH2SO4', 1, ' ') endif end if end do call addfld( 'XPH_LWC', (/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc') call addfld ('AQSO4_H2O2', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to H2O2') call addfld ('AQSO4_O3', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3') if ( history_aerosol ) then call add_default ('XPH_LWC', 1, ' ') call add_default ('AQSO4_H2O2', 1, ' ') call add_default ('AQSO4_O3', 1, ' ') endif endif end subroutine aero_model_init !============================================================================= !============================================================================= subroutine aero_model_drydep ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend ) use dust_sediment_mod, only: dust_sediment_tend use drydep_mod, only: d3ddflux, calcram use modal_aero_data, only: qqcw_get_field use modal_aero_data, only: cnst_name_cw use modal_aero_data, only: alnsg_amode use modal_aero_data, only: sigmag_amode use modal_aero_data, only: nspec_amode use modal_aero_data, only: numptr_amode use modal_aero_data, only: numptrcw_amode use modal_aero_data, only: lmassptr_amode use modal_aero_data, only: lmassptrcw_amode use modal_aero_deposition, only: set_srf_drydep ! args type(physics_state), intent(in) :: state ! Physics state variables real(r8), intent(in) :: obklen(:) real(r8), intent(in) :: ustar(:) ! sfc fric vel type(cam_in_t), target, intent(in) :: cam_in ! import state real(r8), intent(in) :: dt ! time step type(cam_out_t), intent(inout) :: cam_out ! export state type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies type(physics_buffer_desc), pointer :: pbuf(:) ! local vars real(r8), pointer :: landfrac(:) ! land fraction real(r8), pointer :: icefrac(:) ! ice fraction real(r8), pointer :: ocnfrac(:) ! ocean fraction real(r8), pointer :: fvin(:) ! real(r8), pointer :: ram1in(:) ! for dry dep velocities from land model for progseasalts real(r8) :: fv(pcols) ! for dry dep velocities, from land modified over ocean & ice real(r8) :: ram1(pcols) ! for dry dep velocities, from land modified over ocean & ice integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns integer :: jvlc ! index for last dimension of vlc_xxx arrays integer :: lphase ! index for interstitial / cloudborne aerosol integer :: lspec ! index for aerosol number / chem-mass / water-mass integer :: m ! aerosol mode index integer :: mm ! tracer index integer :: i real(r8) :: tvs(pcols,pver) real(r8) :: rho(pcols,pver) ! air density in kg/m3 real(r8) :: sflx(pcols) ! deposition flux real(r8) :: dep_trb(pcols) !kg/m2/s real(r8) :: dep_grv(pcols) !kg/m2/s (total of grav and trb) real(r8) :: pvmzaer(pcols,pverp) ! sedimentation velocity in Pa real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species real(r8) :: rad_drop(pcols,pver) real(r8) :: dens_drop(pcols,pver) real(r8) :: sg_drop(pcols,pver) real(r8) :: rad_aer(pcols,pver) real(r8) :: dens_aer(pcols,pver) real(r8) :: sg_aer(pcols,pver) real(r8) :: vlc_dry(pcols,pver,4) ! dep velocity real(r8) :: vlc_grv(pcols,pver,4) ! dep velocity real(r8):: vlc_trb(pcols,4) ! dep velocity real(r8) :: aerdepdryis(pcols,pcnst) ! aerosol dry deposition (interstitial) real(r8) :: aerdepdrycw(pcols,pcnst) ! aerosol dry deposition (cloud water) real(r8), pointer :: fldcw(:,:) real(r8), pointer :: dgncur_awet(:,:,:) real(r8), pointer :: wetdens(:,:,:) real(r8), pointer :: qaerwat(:,:,:) landfrac => cam_in%landfrac(:) icefrac => cam_in%icefrac(:) ocnfrac => cam_in%ocnfrac(:) fvin => cam_in%fv(:) ram1in => cam_in%ram1(:) lchnk = state%lchnk ncol = state%ncol ! calc ram and fv over ocean and sea ice ... call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,& ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),& state%pdel(:,pver),fvin,fv) call outfld( 'airFV', fv(:), pcols, lchnk ) call outfld( 'RAM1', ram1(:), pcols, lchnk ) ! note that tendencies are not only in sfc layer (because of sedimentation) ! and that ptend is updated within each subroutine for different species call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq) call pbuf_get_field(pbuf, dgnumwet_idx, dgncur_awet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) ) call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) ) call pbuf_get_field(pbuf, qaerwat_idx, qaerwat, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) ) tvs(:ncol,:) = state%t(:ncol,:)!*(1+state%q(:ncol,k) rho(:ncol,:)= state%pmid(:ncol,:)/(rair*state%t(:ncol,:)) ! ! calc settling/deposition velocities for cloud droplets (and cloud-borne aerosols) ! ! *** mean drop radius should eventually be computed from ndrop and qcldwtr rad_drop(:,:) = 5.0e-6_r8 dens_drop(:,:) = rhoh2o sg_drop(:,:) = 1.46_r8 jvlc = 3 call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, & vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), & rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 0, lchnk) jvlc = 4 call modal_aero_depvel_part( ncol,state%t(:,:), state%pmid(:,:), ram1, fv, & vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), & rad_drop(:,:), dens_drop(:,:), sg_drop(:,:), 3, lchnk) do m = 1, ntot_amode ! main loop over aerosol modes do lphase = 1, 2 ! loop over interstitial / cloud-borne forms if (lphase == 1) then ! interstial aerosol - calc settling/dep velocities of mode ! rad_aer = volume mean wet radius (m) ! dgncur_awet = geometric mean wet diameter for number distribution (m) rad_aer(1:ncol,:) = 0.5_r8*dgncur_awet(1:ncol,:,m) & *exp(1.5_r8*(alnsg_amode(m)**2)) ! dens_aer(1:ncol,:) = wet density (kg/m3) dens_aer(1:ncol,:) = wetdens(1:ncol,:,m) sg_aer(1:ncol,:) = sigmag_amode(m) jvlc = 1 call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, & vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), & rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 0, lchnk) jvlc = 2 call modal_aero_depvel_part( ncol, state%t(:,:), state%pmid(:,:), ram1, fv, & vlc_dry(:,:,jvlc), vlc_trb(:,jvlc), vlc_grv(:,:,jvlc), & rad_aer(:,:), dens_aer(:,:), sg_aer(:,:), 3, lchnk) end if do lspec = 0, nspec_amode(m)+1 ! loop over number + constituents + water if (lspec == 0) then ! number if (lphase == 1) then mm = numptr_amode(m) jvlc = 1 else mm = numptrcw_amode(m) jvlc = 3 endif else if (lspec <= nspec_amode(m)) then ! non-water mass if (lphase == 1) then mm = lmassptr_amode(lspec,m) jvlc = 2 else mm = lmassptrcw_amode(lspec,m) jvlc = 4 endif else ! water mass ! bypass dry deposition of aerosol water cycle if (lphase == 1) then mm = 0 ! mm = lwaterptr_amode(m) jvlc = 2 else mm = 0 jvlc = 4 endif endif if (mm <= 0) cycle ! if (lphase == 1) then if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then ptend%lq(mm) = .TRUE. ! use pvprogseasalts instead (means making the top level 0) pvmzaer(:ncol,1)=0._r8 pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc) call outfld( trim(cnst_name(mm))//'DDV', pvmzaer(:,2:pverp), pcols, lchnk ) if(.true.) then ! use phil's method ! convert from meters/sec to pascals/sec ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit ! calculate the tendencies and sfc fluxes from the above velocities call dust_sediment_tend( & ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , & state%q(:,:,mm), pvmzaer, ptend%q(:,:,mm), sflx ) else !use charlie's method call d3ddflux( ncol, vlc_dry(:,:,jvlc), state%q(:,:,mm), state%pmid, & state%pdel, tvs, sflx, ptend%q(:,:,mm), dt ) endif ! apportion dry deposition into turb and gravitational settling for tapes do i=1,ncol dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc) dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc) enddo call outfld( trim(cnst_name(mm))//'DDF', sflx, pcols, lchnk) call outfld( trim(cnst_name(mm))//'TBF', dep_trb, pcols, lchnk ) call outfld( trim(cnst_name(mm))//'GVF', dep_grv, pcols, lchnk ) call outfld( trim(cnst_name(mm))//'DTQ', ptend%q(:,:,mm), pcols, lchnk) aerdepdryis(:ncol,mm) = sflx(:ncol) else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then ! aerosol water ! use pvprogseasalts instead (means making the top level 0) pvmzaer(:ncol,1)=0._r8 pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc) if(.true.) then ! use phil's method ! convert from meters/sec to pascals/sec ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit ! calculate the tendencies and sfc fluxes from the above velocities call dust_sediment_tend( & ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , & qaerwat(:,:,mm), pvmzaer, dqdt_tmp(:,:), sflx ) else !use charlie's method call d3ddflux( ncol, vlc_dry(:,:,jvlc), qaerwat(:,:,mm), state%pmid, & state%pdel, tvs, sflx, dqdt_tmp(:,:), dt ) endif ! apportion dry deposition into turb and gravitational settling for tapes do i=1,ncol dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc) dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc) enddo qaerwat(1:ncol,:,mm) = qaerwat(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) * dt else ! lphase == 2 ! use pvprogseasalts instead (means making the top level 0) pvmzaer(:ncol,1)=0._r8 pvmzaer(:ncol,2:pverp) = vlc_dry(:ncol,:,jvlc) fldcw => qqcw_get_field(pbuf, mm,lchnk) if(.true.) then ! use phil's method ! convert from meters/sec to pascals/sec ! pvprogseasalts(:,1) is assumed zero, use density from layer above in conversion pvmzaer(:ncol,2:pverp) = pvmzaer(:ncol,2:pverp) * rho(:ncol,:)*gravit ! calculate the tendencies and sfc fluxes from the above velocities call dust_sediment_tend( & ncol, dt, state%pint(:,:), state%pmid, state%pdel, state%t , & fldcw(:,:), pvmzaer, dqdt_tmp(:,:), sflx ) else !use charlie's method call d3ddflux( ncol, vlc_dry(:,:,jvlc), fldcw(:,:), state%pmid, & state%pdel, tvs, sflx, dqdt_tmp(:,:), dt ) endif ! apportion dry deposition into turb and gravitational settling for tapes do i=1,ncol dep_trb(i)=sflx(i)*vlc_trb(i,jvlc)/vlc_dry(i,pver,jvlc) dep_grv(i)=sflx(i)*vlc_grv(i,pver,jvlc)/vlc_dry(i,pver,jvlc) enddo fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt call outfld( trim(cnst_name_cw(mm))//'DDF', sflx, pcols, lchnk) call outfld( trim(cnst_name_cw(mm))//'TBF', dep_trb, pcols, lchnk ) call outfld( trim(cnst_name_cw(mm))//'GVF', dep_grv, pcols, lchnk ) aerdepdrycw(:ncol,mm) = sflx(:ncol) endif enddo ! lspec = 0, nspec_amode(m)+1 enddo ! lphase = 1, 2 enddo ! m = 1, ntot_amode ! if the user has specified prescribed aerosol dep fluxes then ! do not set cam_out dep fluxes according to the prognostic aerosols if (.not.aerodep_flx_prescribed()) then call set_srf_drydep(aerdepdryis, aerdepdrycw, cam_out) endif endsubroutine aero_model_drydep !============================================================================= !============================================================================= subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf) use modal_aero_deposition, only: set_srf_wetdep use wetdep, only: wetdepa_v2, wetdep_inputs_set, wetdep_inputs_t use modal_aero_data use modal_aero_calcsize, only: modal_aero_calcsize_sub use modal_aero_wateruptake,only: modal_aero_wateruptake_dr use modal_aero_convproc, only: deepconv_wetdep_history, ma_convproc_intr ! args type(physics_state), intent(in) :: state ! Physics state variables real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s] type(cam_out_t), intent(inout) :: cam_out ! export state type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies type(physics_buffer_desc), pointer :: pbuf(:) ! local vars integer :: m ! tracer index integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns real(r8) :: iscavt(pcols, pver) integer :: mm integer :: i,k real(r8) :: icscavt(pcols, pver) real(r8) :: isscavt(pcols, pver) real(r8) :: bcscavt(pcols, pver) real(r8) :: bsscavt(pcols, pver) real(r8) :: sol_factb, sol_facti real(r8) :: sol_factic(pcols,pver) real(r8) :: sflx(pcols) ! deposition flux integer :: jnv ! index for scavcoefnv 3rd dimension integer :: lphase ! index for interstitial / cloudborne aerosol integer :: strt_loop, end_loop, stride_loop !loop indices for the lphase loop integer :: lspec ! index for aerosol number / chem-mass / water-mass integer :: lcoardust, lcoarnacl ! indices for coarse mode dust and seasalt masses real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for 1 species real(r8) :: f_act_conv(pcols,pver) ! prescribed aerosol activation fraction for convective cloud ! rce 2010/05/01 real(r8) :: f_act_conv_coarse(pcols,pver) ! similar but for coarse mode ! rce 2010/05/02 real(r8) :: f_act_conv_coarse_dust, f_act_conv_coarse_nacl ! rce 2010/05/02 real(r8) :: fracis_cw(pcols,pver) real(r8) :: hygro_sum_old(pcols,pver) ! before removal [sum of (mass*hydro/dens)] real(r8) :: hygro_sum_del(pcols,pver) ! removal change to [sum of (mass*hydro/dens)] real(r8) :: hygro_sum_old_ik, hygro_sum_new_ik real(r8) :: prec(pcols) ! precipitation rate real(r8) :: q_tmp(pcols,pver) ! temporary array to hold "most current" mixing ratio for 1 species real(r8) :: scavcoefnv(pcols,pver,0:2) ! Dana and Hales coefficient (/mm) for ! cloud-borne num & vol (0), ! interstitial num (1), interstitial vol (2) real(r8) :: tmpa, tmpb real(r8) :: tmpdust, tmpnacl real(r8) :: water_old, water_new ! temporary old/new aerosol water mix-rat logical :: isprx(pcols,pver) ! true if precipation real(r8) :: aerdepwetis(pcols,pcnst) ! aerosol wet deposition (interstitial) real(r8) :: aerdepwetcw(pcols,pcnst) ! aerosol wet deposition (cloud water) ! For unified convection scheme logical, parameter :: do_aero_water_removal = .false. ! True if aerosol water reduction by wet removal is to be calculated ! (this has not been fully tested, so best to leave it off) logical :: do_hygro_sum_del, do_lphase1, do_lphase2 real(r8), pointer :: rprddp(:,:) ! rain production, deep convection real(r8), pointer :: rprdsh(:,:) ! rain production, shallow convection real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0. real(r8), pointer :: evapcsh(:,:) ! Evaporation rate of shallow convective precipitation >=0. real(r8) :: rprddpsum(pcols) real(r8) :: rprdshsum(pcols) real(r8) :: evapcdpsum(pcols) real(r8) :: evapcshsum(pcols) real(r8) :: tmp_resudp, tmp_resush real(r8) :: sflxec(pcols), sflxecdp(pcols) ! deposition flux real(r8) :: sflxic(pcols), sflxicdp(pcols) ! deposition flux real(r8) :: sflxbc(pcols), sflxbcdp(pcols) ! deposition flux real(r8) :: rcscavt(pcols, pver) real(r8) :: rsscavt(pcols, pver) real(r8) :: qqcw_in(pcols,pver), qqcw_sav(pcols,pver,0:nspec_max) ! temporary array to hold qqcw for the current mode real(r8) :: rtscavt(pcols, pver, 0:nspec_max) integer, parameter :: nsrflx_mzaer2cnvpr = 2 real(r8) :: qsrflx_mzaer2cnvpr(pcols,pcnst,nsrflx_mzaer2cnvpr) ! End unified convection scheme real(r8), pointer :: fldcw(:,:) real(r8), pointer :: dgnumwet(:,:,:) real(r8), pointer :: qaerwat(:,:,:) ! aerosol water real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble type(wetdep_inputs_t) :: dep_inputs lchnk = state%lchnk ncol = state%ncol call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep', lq=wetdep_lq) ! Do calculations of mode radius and water uptake if: ! 1) modal aerosols are affecting the climate, or ! 2) prognostic modal aerosols are enabled call t_startf('calcsize') ! for prognostic modal aerosols the transfer of mass between aitken and accumulation ! modes is done in conjunction with the dry radius calculation call modal_aero_calcsize_sub(state, ptend, dt, pbuf) call t_stopf('calcsize') call t_startf('wateruptake') call modal_aero_wateruptake_dr(state, pbuf) call t_stopf('wateruptake') if (nwetdep<1) return call wetdep_inputs_set( state, pbuf, dep_inputs ) call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) ) call pbuf_get_field(pbuf, qaerwat_idx, qaerwat, start=(/1,1,1/), kount=(/pcols,pver,nmodes/) ) call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) ) prec(:ncol)=0._r8 do k=1,pver where (prec(:ncol) >= 1.e-7_r8) isprx(:ncol,k) = .true. elsewhere isprx(:ncol,k) = .false. endwhere prec(:ncol) = prec(:ncol) + (dep_inputs%prain(:ncol,k) + dep_inputs%cmfdqr(:ncol,k) - dep_inputs%evapr(:ncol,k)) & *state%pdel(:ncol,k)/gravit end do if(convproc_do_aer) then qsrflx_mzaer2cnvpr(:,:,:) = 0.0_r8 aerdepwetis(:,:) = 0.0_r8 aerdepwetcw(:,:) = 0.0_r8 else qsrflx_mzaer2cnvpr(:,:,:) = nan aerdepwetis(:,:) = nan aerdepwetcw(:,:) = nan endif ! calculate the mass-weighted sol_factic for coarse mode species ! sol_factic_coarse(:,:) = 0.30_r8 ! tuned 1/4 f_act_conv_coarse(:,:) = 0.60_r8 ! rce 2010/05/02 f_act_conv_coarse_dust = 0.40_r8 ! rce 2010/05/02 f_act_conv_coarse_nacl = 0.80_r8 ! rce 2010/05/02 if (modeptr_coarse > 0) then lcoardust = lptr_dust_a_amode(modeptr_coarse) lcoarnacl = lptr_nacl_a_amode(modeptr_coarse) if ((lcoardust > 0) .and. (lcoarnacl > 0)) then do k = 1, pver do i = 1, ncol tmpdust = max( 0.0_r8, state%q(i,k,lcoardust) + ptend%q(i,k,lcoardust)*dt ) tmpnacl = max( 0.0_r8, state%q(i,k,lcoarnacl) + ptend%q(i,k,lcoarnacl)*dt ) if ((tmpdust+tmpnacl) > 1.0e-30_r8) then ! sol_factic_coarse(i,k) = (0.2_r8*tmpdust + 0.4_r8*tmpnacl)/(tmpdust+tmpnacl) ! tuned 1/6 f_act_conv_coarse(i,k) = (f_act_conv_coarse_dust*tmpdust & + f_act_conv_coarse_nacl*tmpnacl)/(tmpdust+tmpnacl) ! rce 2010/05/02 end if end do end do end if end if scavcoefnv(:,:,0) = 0.0_r8 ! below-cloud scavcoef = 0.0 for cloud-borne species ! Counters for "without" unified convective treatment (i.e. default case) strt_loop = 1 end_loop = 2 stride_loop = 1 if (convproc_do_aer) then !Do cloudborne first for unified convection scheme so that the resuspension of cloudborne !can be saved then applied to interstitial strt_loop = 2 end_loop = 1 stride_loop = -1 endif do m = 1, ntot_amode ! main loop over aerosol modes do lphase = strt_loop,end_loop, stride_loop ! loop over interstitial (1) and cloud-borne (2) forms ! sol_factb and sol_facti values ! sol_factb - currently this is basically a tuning factor ! sol_facti & sol_factic - currently has a physical basis, and reflects activation fraction ! ! 2008-mar-07 rce - sol_factb (interstitial) changed from 0.3 to 0.1 ! - sol_factic (interstitial, dust modes) changed from 1.0 to 0.5 ! - sol_factic (cloud-borne, pcarb modes) no need to set it to 0.0 ! because the cloud-borne pcarbon == 0 (no activation) ! ! rce 2010/05/02 ! prior to this date, sol_factic was used for convective in-cloud wet removal, ! and its value reflected a combination of an activation fraction (which varied between modes) ! and a tuning factor ! from this date forward, two parameters are used for convective in-cloud wet removal ! f_act_conv is the activation fraction ! note that "non-activation" of aerosol in air entrained into updrafts should ! be included here ! eventually we might use the activate routine (with w ~= 1 m/s) to calculate ! this, but there is still the entrainment issue ! sol_factic is strictly a tuning factor ! if (lphase == 1) then ! interstial aerosol hygro_sum_old(:,:) = 0.0_r8 hygro_sum_del(:,:) = 0.0_r8 call modal_aero_bcscavcoef_get( m, ncol, isprx, dgnumwet, & scavcoefnv(:,:,1), scavcoefnv(:,:,2) ) sol_factb = sol_factb_interstitial ! all below-cloud scav ON (0.1 "tuning factor") sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial sol_factic = sol_factic_interstitial if (m == modeptr_pcarbon) then ! sol_factic = 0.0_r8 ! conv in-cloud scav OFF (0.0 activation fraction) f_act_conv = 0.0_r8 ! rce 2010/05/02 else if ((m == modeptr_finedust) .or. (m == modeptr_coardust)) then ! sol_factic = 0.2_r8 ! conv in-cloud scav ON (0.5 activation fraction) ! tuned 1/4 f_act_conv = 0.4_r8 ! rce 2010/05/02 else ! sol_factic = 0.4_r8 ! conv in-cloud scav ON (1.0 activation fraction) ! tuned 1/4 f_act_conv = 0.8_r8 ! rce 2010/05/02 end if else ! cloud-borne aerosol (borne by stratiform cloud drops) sol_factb = 0.0_r8 ! all below-cloud scav OFF (anything cloud-borne is located "in-cloud") sol_facti = sol_facti_cloud_borne ! strat in-cloud scav cloud-borne tuning factor sol_factic = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean ! that conv precip collects strat droplets) f_act_conv = 0.0_r8 ! conv in-cloud scav OFF (having this on would mean end if if (convproc_do_aer .and. lphase == 1) then ! if modal aero convproc is turned on for aerosols, then ! turn off the convective in-cloud removal for interstitial aerosols ! (but leave the below-cloud on, as convproc only does in-cloud) ! and turn off the outfld SFWET, SFSIC, SFSID, SFSEC, and SFSED calls ! for (stratiform)-cloudborne aerosols, convective wet removal ! (all forms) is zero, so no action is needed sol_factic = 0.0_r8 endif ! ! rce 2010/05/03 ! wetdepa has "sol_fact" parameters: ! sol_facti, sol_factic, sol_factb for liquid cloud do lspec = 0, nspec_amode(m)+1 ! loop over number + chem constituents + water if (lspec == 0) then ! number if (lphase == 1) then mm = numptr_amode(m) jnv = 1 else mm = numptrcw_amode(m) jnv = 0 endif else if (lspec <= nspec_amode(m)) then ! non-water mass if (lphase == 1) then mm = lmassptr_amode(lspec,m) jnv = 2 else mm = lmassptrcw_amode(lspec,m) jnv = 0 endif else ! water mass ! bypass wet removal of aerosol water if(convproc_do_aer) then if ( .not. do_aero_water_removal ) cycle else cycle endif if (lphase == 1) then mm = 0 ! mm = lwaterptr_amode(m) jnv = 2 else mm = 0 jnv = 0 endif endif if (mm <= 0) cycle ! set f_act_conv for interstitial (lphase=1) coarse mode species ! for the convective in-cloud, we conceptually treat the coarse dust and seasalt ! as being externally mixed, and apply f_act_conv = f_act_conv_coarse_dust/nacl to dust/seasalt ! number and sulfate are conceptually partitioned to the dust and seasalt ! on a mass basis, so the f_act_conv for number and sulfate are ! mass-weighted averages of the values used for dust/seasalt if ((lphase == 1) .and. (m == modeptr_coarse)) then ! sol_factic = sol_factic_coarse f_act_conv = f_act_conv_coarse ! rce 2010/05/02 if (lspec > 0) then if (lmassptr_amode(lspec,m) == lptr_dust_a_amode(m)) then ! sol_factic = 0.2_r8 ! tuned 1/4 f_act_conv = f_act_conv_coarse_dust ! rce 2010/05/02 else if (lmassptr_amode(lspec,m) == lptr_nacl_a_amode(m)) then ! sol_factic = 0.4_r8 ! tuned 1/6 f_act_conv = f_act_conv_coarse_nacl ! rce 2010/05/02 end if end if end if if ((lphase == 1) .and. (lspec <= nspec_amode(m))) then ptend%lq(mm) = .TRUE. dqdt_tmp(:,:) = 0.0_r8 ! q_tmp reflects changes from modal_aero_calcsize and is the "most current" q q_tmp(1:ncol,:) = state%q(1:ncol,:,mm) + ptend%q(1:ncol,:,mm)*dt if(convproc_do_aer) then !Feed in the saved cloudborne mixing ratios from phase 2 qqcw_in(:,:) = qqcw_sav(:,:,lspec) else fldcw => qqcw_get_field(pbuf, mm,lchnk) qqcw_in(:,:) = fldcw(:,:) endif call wetdepa_v2( state%pmid, state%q(:,:,1), state%pdel, & dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, & dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, & dep_inputs%evapr, dep_inputs%totcond, q_tmp, dt, & dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, & dlf, fracis(:,:,mm), sol_factb, ncol, & scavcoefnv(:,:,jnv), & is_strat_cloudborne=.false., & qqcw=qqcw_in(:,:), & f_act_conv=f_act_conv, & icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, & convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt, & sol_facti_in=sol_facti, sol_factic_in=sol_factic ) do_hygro_sum_del = .false. if ( lspec > 0 ) do_hygro_sum_del = .true. if(convproc_do_aer) then do_hygro_sum_del = .false. ! add resuspension of cloudborne species to dqdt of interstitial species dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt(1:ncol,:,lspec) if ( (lspec > 0) .and. do_aero_water_removal ) then do_hygro_sum_del = .true. endif endif ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) call outfld( trim(cnst_name(mm))//'WET', dqdt_tmp(:,:), pcols, lchnk) call outfld( trim(cnst_name(mm))//'SIC', icscavt, pcols, lchnk) call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk) call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk) call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit enddo enddo if (.not.convproc_do_aer) call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk) aerdepwetis(:ncol,mm) = sflx(:ncol) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit enddo enddo if (.not.convproc_do_aer) call outfld( trim(cnst_name(mm))//'SFSIC', sflx, pcols, lchnk) if (convproc_do_aer) sflxic = sflx sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name(mm))//'SFSIS', sflx, pcols, lchnk) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name(mm))//'SFSBC', sflx, pcols, lchnk) if (convproc_do_aer)sflxbc = sflx sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name(mm))//'SFSBS', sflx, pcols, lchnk) if (convproc_do_aer) then sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+rcscavt(i,k)*state%pdel(i,k)/gravit enddo enddo sflxec = sflx sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+rsscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name(mm))//'SFSES', sflx, pcols, lchnk) ! apportion convective surface fluxes to deep and shallow conv ! this could be done more accurately in subr wetdepa ! since deep and shallow rarely occur simultaneously, and these ! fields are just diagnostics, this approximate method is adequate ! only do this for interstitial aerosol, because conv clouds to not ! affect the stratiform-cloudborne aerosol if ( deepconv_wetdep_history) then call pbuf_get_field(pbuf, rprddp_idx, rprddp ) call pbuf_get_field(pbuf, rprdsh_idx, rprdsh ) call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp ) call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh ) rprddpsum(:) = 0.0_r8 rprdshsum(:) = 0.0_r8 evapcdpsum(:) = 0.0_r8 evapcshsum(:) = 0.0_r8 do k = 1, pver rprddpsum(:ncol) = rprddpsum(:ncol) + rprddp(:ncol,k)*state%pdel(:ncol,k)/gravit rprdshsum(:ncol) = rprdshsum(:ncol) + rprdsh(:ncol,k)*state%pdel(:ncol,k)/gravit evapcdpsum(:ncol) = evapcdpsum(:ncol) + evapcdp(:ncol,k)*state%pdel(:ncol,k)/gravit evapcshsum(:ncol) = evapcshsum(:ncol) + evapcsh(:ncol,k)*state%pdel(:ncol,k)/gravit end do do i = 1, ncol rprddpsum(i) = max( rprddpsum(i), 1.0e-35_r8 ) rprdshsum(i) = max( rprdshsum(i), 1.0e-35_r8 ) evapcdpsum(i) = max( evapcdpsum(i), 0.1e-35_r8 ) evapcshsum(i) = max( evapcshsum(i), 0.1e-35_r8 ) ! assume that in- and below-cloud removal are proportional to column precip production tmpa = rprddpsum(i) / (rprddpsum(i) + rprdshsum(i)) tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) ) sflxicdp(i) = sflxic(i)*tmpa sflxbcdp(i) = sflxbc(i)*tmpa ! assume that resuspension is proportional to (wet removal)*[(precip evap)/(precip production)] tmp_resudp = tmpa * min( (evapcdpsum(i)/rprddpsum(i)), 1.0_r8 ) tmp_resush = (1.0_r8 - tmpa) * min( (evapcshsum(i)/rprdshsum(i)), 1.0_r8 ) tmpb = max( tmp_resudp, 1.0e-35_r8 ) / max( (tmp_resudp+tmp_resush), 1.0e-35_r8 ) tmpb = max( 0.0_r8, min( 1.0_r8, tmpb ) ) sflxecdp(i) = sflxec(i)*tmpb end do call outfld( trim(cnst_name(mm))//'SFSBD', sflxbcdp, pcols, lchnk) else sflxec(1:ncol) = 0.0_r8 sflxecdp(1:ncol) = 0.0_r8 end if ! when ma_convproc_intr is used, convective in-cloud wet removal is done there ! the convective (total and deep) precip-evap-resuspension includes in- and below-cloud ! contributions ! so pass the below-cloud contribution to ma_convproc_intr qsrflx_mzaer2cnvpr(1:ncol,mm,1) = sflxec( 1:ncol) qsrflx_mzaer2cnvpr(1:ncol,mm,2) = sflxecdp(1:ncol) endif if (do_hygro_sum_del) then tmpa = spechygro(lspec,m)/ & specdens_amode(lspec,m) tmpb = tmpa*dt hygro_sum_old(1:ncol,:) = hygro_sum_old(1:ncol,:) & + tmpa*q_tmp(1:ncol,:) hygro_sum_del(1:ncol,:) = hygro_sum_del(1:ncol,:) & + tmpb*dqdt_tmp(1:ncol,:) end if else if ((lphase == 1) .and. (lspec == nspec_amode(m)+1)) then do_lphase1 = .true. if(convproc_do_aer) then do_lphase1 = .false. if(do_aero_water_removal)do_lphase1 = .true. endif if(do_lphase1) then ! aerosol water -- because of how wetdepa treats evaporation of stratiform ! precip, it is not appropriate to apply wetdepa to aerosol water ! instead, "hygro_sum" = [sum of (mass*hygro/dens)] is calculated before and ! after wet removal, and new water is calculated using ! new_water = old_water*min(10,(hygro_sum_new/hygro_sum_old)) ! the "min(10,...)" is to avoid potential problems when hygro_sum_old ~= 0 ! also, individual wet removal terms (ic,is,bc,bs) are not output to history ! ptend%lq(mm) = .TRUE. ! dqdt_tmp(:,:) = 0.0_r8 do k = 1, pver do i = 1, ncol ! water_old = max( 0.0_r8, state%q(i,k,mm)+ptend%q(i,k,mm)*dt ) water_old = max( 0.0_r8, qaerwat(i,k,mm) ) hygro_sum_old_ik = max( 0.0_r8, hygro_sum_old(i,k) ) hygro_sum_new_ik = max( 0.0_r8, hygro_sum_old_ik+hygro_sum_del(i,k) ) if (hygro_sum_new_ik >= 10.0_r8*hygro_sum_old_ik) then water_new = 10.0_r8*water_old else water_new = water_old*(hygro_sum_new_ik/hygro_sum_old_ik) end if ! dqdt_tmp(i,k) = (water_new - water_old)/dt qaerwat(i,k,mm) = water_new end do end do ! ptend%q(1:ncol,:,mm) = ptend%q(1:ncol,:,mm) + dqdt_tmp(1:ncol,:) ! call outfld( trim(cnst_name(mm)) ! sflx(:)=0._r8 ! do k=1,pver ! do i=1,ncol ! sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit ! enddo ! enddo ! call outfld( trim(cnst_name(mm)) endif elseif (lphase == 2) then do_lphase2 = .true. if (convproc_do_aer) then do_lphase2 = .false. if (lspec <= nspec_amode(m)) do_lphase2 = .true. endif if (do_lphase2) then dqdt_tmp(:,:) = 0.0_r8 if (convproc_do_aer) then fldcw => qqcw_get_field(pbuf,mm,lchnk) qqcw_sav(1:ncol,:,lspec) = fldcw(1:ncol,:) else fldcw => qqcw_get_field(pbuf, mm,lchnk) endif call wetdepa_v2(state%pmid, state%q(:,:,1), state%pdel, & dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, & dep_inputs%evapc, dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, & dep_inputs%evapr, dep_inputs%totcond, fldcw, dt, & dqdt_tmp, iscavt, dep_inputs%cldvcu, dep_inputs%cldvst, & dlf, fracis_cw, sol_factb, ncol, & scavcoefnv(:,:,jnv), & is_strat_cloudborne=.true., & icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt, & convproc_do_aer=convproc_do_aer, rcscavt=rcscavt, rsscavt=rsscavt, & sol_facti_in=sol_facti, sol_factic_in=sol_factic ) if(convproc_do_aer) then ! save resuspension of cloudborne species rtscavt(1:ncol,:,lspec) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ! wetdepa_v2 adds the resuspension of cloudborne to the dqdt of cloudborne (as a source) ! undo this, so the resuspension of cloudborne can be added to the dqdt of interstitial (above) dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt(1:ncol,:,lspec) endif fldcw(1:ncol,:) = fldcw(1:ncol,:) + dqdt_tmp(1:ncol,:) * dt sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+dqdt_tmp(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name_cw(mm))//'SFWET', sflx, pcols, lchnk) aerdepwetcw(:ncol,mm) = sflx(:ncol) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+icscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name_cw(mm))//'SFSIC', sflx, pcols, lchnk) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+isscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name_cw(mm))//'SFSIS', sflx, pcols, lchnk) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+bcscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name_cw(mm))//'SFSBC', sflx, pcols, lchnk) sflx(:)=0._r8 do k=1,pver do i=1,ncol sflx(i)=sflx(i)+bsscavt(i,k)*state%pdel(i,k)/gravit enddo enddo call outfld( trim(cnst_name_cw(mm))//'SFSBS', sflx, pcols, lchnk) if(convproc_do_aer) then sflx(:)=0.0_r8 do k=1,pver sflx(1:ncol)=sflx(1:ncol)+rcscavt(1:ncol,k)*state%pdel(1:ncol,k)/gravit enddo call outfld( trim(cnst_name_cw(mm))//'SFSEC', sflx, pcols, lchnk) sflx(:)=0.0_r8 do k=1,pver sflx(1:ncol)=sflx(1:ncol)+rsscavt(1:ncol,k)*state%pdel(1:ncol,k)/gravit enddo call outfld( trim(cnst_name_cw(mm))//'SFSES', sflx, pcols, lchnk) endif endif endif enddo ! lspec = 0, nspec_amode(m)+1 enddo ! lphase = 1, 2 enddo ! m = 1, ntot_amode if (convproc_do_aer) then call t_startf('ma_convproc') call ma_convproc_intr( state, ptend, pbuf, dt, & nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis) call t_stopf('ma_convproc') endif ! if the user has specified prescribed aerosol dep fluxes then ! do not set cam_out dep fluxes according to the prognostic aerosols if (.not. aerodep_flx_prescribed()) then call set_srf_wetdep(aerdepwetis, aerdepwetcw, cam_out) endif endsubroutine aero_model_wetdep !------------------------------------------------------------------------- ! provides wet tropospheric aerosol surface area info for modal aerosols ! called from mo_usrrxt !------------------------------------------------------------------------- subroutine aero_model_surfarea( & mmr, radmean, relhum, pmid, temp, strato_sad, & sulfate, rho, ltrop, dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total ) ! dummy args real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: radmean ! mean radii in cm real(r8), intent(in) :: strato_sad(:,:) integer, intent(in) :: ncol integer, intent(in) :: ltrop(:) real(r8), intent(in) :: dlat(:) ! degrees latitude integer, intent(in) :: het1_ndx real(r8), intent(in) :: relhum(:,:) real(r8), intent(in) :: rho(:,:) ! total atm density (/cm^3) real(r8), intent(in) :: sulfate(:,:) type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(inout) :: sfc(:,:,:) real(r8), intent(inout) :: dm_aer(:,:,:) real(r8), intent(inout) :: sad_total(:,:) ! local vars real(r8), pointer, dimension(:,:,:) :: dgnumwet integer :: beglev(ncol) integer :: endlev(ncol) integer :: i,k call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) beglev(:ncol)=ltrop(:ncol)+1 endlev(:ncol)=pver call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, sad_total, sfc=sfc ) do i = 1,ncol do k = ltrop(i)+1,pver dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm enddo enddo end subroutine aero_model_surfarea !------------------------------------------------------------------------- ! provides WET stratospheric aerosol surface area info for modal aerosols ! if modal_strat_sulfate = TRUE -- called from mo_gas_phase_chemdr !------------------------------------------------------------------------- subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad ) ! dummy args integer, intent(in) :: ncol real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) integer, intent(in) :: ltrop(:) ! tropopause level indices type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(out) :: strato_sad(:,:) ! local vars real(r8), pointer, dimension(:,:,:) :: dgnumwet integer :: beglev(ncol) integer :: endlev(ncol) strato_sad = 0._r8 if (.not.modal_strat_sulfate) return call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) beglev(:ncol)=top_lev endlev(:ncol)=ltrop(:ncol) call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, strato_sad ) end subroutine aero_model_strat_surfarea !============================================================================= !============================================================================= subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, & tfld, pmid, pdel, mbar, relhum, & zm, qh2o, cwat, cldfr, cldnum, & airdens, invariants, del_h2so4_gasprod, & vmr0, vmr, pbuf ) use time_manager, only : get_nstep use modal_aero_coag, only : modal_aero_coag_sub use modal_aero_gasaerexch, only : modal_aero_gasaerexch_sub use modal_aero_newnuc, only : modal_aero_newnuc_sub use modal_aero_data, only : cnst_name_cw, qqcw_get_field !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: loffset ! offset applied to modal aero "pointers" integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: lchnk ! chunk index integer, intent(in) :: troplev(pcols) real(r8), intent(in) :: delt ! time step size (sec) real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates real(r8), intent(in) :: tfld(:,:) ! temperature (K) real(r8), intent(in) :: pmid(:,:) ! pressure at model levels (Pa) real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa) real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu ) real(r8), intent(in) :: relhum(:,:) ! relative humidity real(r8), intent(in) :: airdens(:,:) ! total atms density (molec/cm**3) real(r8), intent(in) :: invariants(:,:,:) real(r8), intent(in) :: del_h2so4_gasprod(:,:) real(r8), intent(in) :: zm(:,:) real(r8), intent(in) :: qh2o(:,:) real(r8), intent(in) :: cwat(:,:) ! cloud liquid water content (kg/kg) real(r8), intent(in) :: cldfr(:,:) real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg) real(r8), intent(in) :: vmr0(:,:,:) ! initial mixing ratios (before gas-phase chem changes) real(r8), intent(inout) :: vmr(:,:,:) ! mixing ratios ( vmr ) type(physics_buffer_desc), pointer :: pbuf(:) ! local vars integer :: n, m integer :: i,k,l integer :: nstep real(r8) :: del_h2so4_aeruptk(ncol,pver) real(r8), pointer :: dgnum(:,:,:), dgnumwet(:,:,:), wetdens(:,:,:) real(r8), pointer :: pblh(:) ! pbl height (m) real(r8), dimension(ncol) :: wrk character(len=32) :: name real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst) real(r8) :: dvmrdt(ncol,pver,gas_pcnst) real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr) real(r8) :: aqso4(ncol,ntot_amode) ! aqueous phase chemistry real(r8) :: aqh2so4(ncol,ntot_amode) ! aqueous phase chemistry real(r8) :: aqso4_h2o2(ncol) ! SO4 aqueous phase chemistry due to H2O2 real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3 real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc real(r8), pointer :: fldcw(:,:) real(r8), pointer :: sulfeq(:,:,:) logical :: is_spcam_m2005 is_spcam_m2005 = cam_physpkg_is('spcam_m2005') call pbuf_get_field(pbuf, dgnum_idx, dgnum) call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens ) call pbuf_get_field(pbuf, pblh_idx, pblh) do n=1,ntot_amode call outfld(dgnum_name(n), dgnum(1:ncol,1:pver,n), ncol, lchnk ) call outfld(dgnumwet_name(n), dgnumwet(1:ncol,1:pver,n), ncol, lchnk ) end do ! do gas-aerosol exchange (h2so4, msa, nh3 condensation) nstep = get_nstep() ! calculate tendency due to gas phase chemistry and processes dvmrdt(:ncol,:,:) = (vmr(:ncol,:,:) - vmr0(:ncol,:,:)) / delt do m = 1, gas_pcnst wrk(:) = 0._r8 do k = 1,pver wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m)*adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit end do name = 'GS_'//trim(solsym(m)) call outfld( name, wrk(:ncol), ncol, lchnk ) enddo ! ! Aerosol processes ... ! call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) if (.not. is_spcam_m2005) then ! regular CAM dvmrdt(:ncol,:,:) = vmr(:ncol,:,:) dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:) ! aqueous chemistry ... if( has_sox ) then call setsox( & ncol, & lchnk, & loffset, & delt, & pmid, & pdel, & tfld, & mbar, & cwat, & cldfr, & cldnum, & airdens, & invariants, & vmrcw, & vmr, & xphlwc, & aqso4, & aqh2so4, & aqso4_h2o2, & aqso4_o3 & ) do n = 1, ntot_amode l = lptr_so4_cw_amode(n) if (l > 0) then call outfld( trim(cnst_name_cw(l))//'AQSO4', aqso4(:ncol,n), ncol, lchnk) call outfld( trim(cnst_name_cw(l))//'AQH2SO4', aqh2so4(:ncol,n), ncol, lchnk) end if end do call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk) call outfld( 'AQSO4_O3', aqso4_o3(:ncol), ncol, lchnk) call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol, lchnk ) endif ! Tendency due to aqueous chemistry dvmrdt = (vmr - dvmrdt) / delt dvmrcwdt = (vmrcw - dvmrcwdt) / delt do m = 1, gas_pcnst wrk(:) = 0._r8 do k = 1,pver wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit end do name = 'AQ_'//trim(solsym(m)) call outfld( name, wrk(:ncol), ncol, lchnk ) enddo else if (is_spcam_m2005) then ! SPCAM ECPP ! when ECPP is used, aqueous chemistry is done in ECPP, ! and not updated here. ! Minghuai Wang, 2010-02 (Minghuai.Wang@pnl.gov) ! dvmrdt = 0.0_r8 dvmrcwdt = 0.0_r8 endif ! do gas-aerosol exchange (h2so4, msa, nh3 condensation) if (ndx_h2so4 > 0) then del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) else del_h2so4_aeruptk(:,:) = 0.0_r8 endif call t_startf('modal_gas-aer_exchng') if ( sulfeq_idx>0 ) then call pbuf_get_field( pbuf, sulfeq_idx, sulfeq ) else nullify( sulfeq ) endif call modal_aero_gasaerexch_sub( & lchnk, ncol, nstep, & loffset, delt, & tfld, pmid, pdel, & qh2o, troplev, & vmr, vmrcw, & dvmrdt, dvmrcwdt, & dgnum, dgnumwet, & sulfeq ) if (ndx_h2so4 > 0) then del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,ndx_h2so4) - del_h2so4_aeruptk(1:ncol,:) endif call t_stopf('modal_gas-aer_exchng') call t_startf('modal_nucl') ! do aerosol nucleation (new particle formation) call modal_aero_newnuc_sub( & lchnk, ncol, nstep, & loffset, delt, & tfld, pmid, pdel, & zm, pblh, & qh2o, cldfr, & vmr, & del_h2so4_gasprod, del_h2so4_aeruptk ) call t_stopf('modal_nucl') call t_startf('modal_coag') ! do aerosol coagulation call modal_aero_coag_sub( & lchnk, ncol, nstep, & loffset, delt, & tfld, pmid, pdel, & vmr, & dgnum, dgnumwet, & wetdens ) call t_stopf('modal_coag') call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) ! diagnostics for cloud-borne aerosols... do n = 1,pcnst fldcw => qqcw_get_field(pbuf,n,lchnk,errorhandle=.true.) if(associated(fldcw)) then call outfld( cnst_name_cw(n), fldcw(:,:), pcols, lchnk ) endif end do end subroutine aero_model_gasaerexch !============================================================================= !============================================================================= subroutine aero_model_emissions( state, cam_in ) use seasalt_model, only: seasalt_emis, seasalt_names, seasalt_indices, seasalt_active,seasalt_nbin use dust_model, only: dust_emis, dust_names, dust_indices, dust_active,dust_nbin, dust_nnum use physics_types, only: physics_state ! Arguments: type(physics_state), intent(in) :: state ! Physics state variables type(cam_in_t), intent(inout) :: cam_in ! import state ! local vars integer :: lchnk, ncol integer :: m, mm real(r8) :: soil_erod_tmp(pcols) real(r8) :: sflx(pcols) ! accumulate over all bins for output real(r8) :: u10cubed(pcols) real (r8), parameter :: z0=0.0001_r8 ! m roughness length over oceans--from ocean model lchnk = state%lchnk ncol = state%ncol if (dust_active) then call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp ) ! some dust emis diagnostics ... sflx(:)=0._r8 do m=1,dust_nbin+dust_nnum mm = dust_indices(m) if (m<=dust_nbin) sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm) call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk) enddo call outfld('DSTSFMBL',sflx(:),pcols,lchnk) call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk ) endif if (seasalt_active) then u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2) ! move the winds to 10m high from the midpoint of the gridbox: ! follows Tie and Seinfeld and Pandis, p.859 with math. u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0) ! we need them to the 3.41 power, according to Gong et al., 1997: u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8 sflx(:)=0._r8 call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx ) do m=1,seasalt_nbin mm = seasalt_indices(m) sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm) call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk) enddo call outfld('SSTSFMBL',sflx(:),pcols,lchnk) endif end subroutine aero_model_emissions !=============================================================================== ! private methods !============================================================================= !============================================================================= subroutine surf_area_dens( ncol, mmr, pmid, temp, diam, beglev, endlev, sad, sfc ) use mo_constants, only : pi use modal_aero_data, only : nspec_amode, alnsg_amode ! dummy args integer, intent(in) :: ncol real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) real(r8), intent(in) :: diam(:,:,:) integer, intent(in) :: beglev(:) integer, intent(in) :: endlev(:) real(r8), intent(out) :: sad(:,:) real(r8),optional, intent(out) :: sfc(:,:,:) ! local vars real(r8) :: sad_mode(pcols,pver,ntot_amode) real(r8) :: rho_air integer :: i,k,l,m real(r8) :: chm_mass, tot_mass ! ! Compute surface aero for each mode. ! Total over all modes as the surface area for chemical reactions. ! sad = 0._r8 sad_mode = 0._r8 do i = 1,ncol do k = beglev(i),endlev(i) rho_air = pmid(i,k)/(temp(i,k)*287.04_r8) do l=1,ntot_amode ! ! compute a mass weighting of the number ! tot_mass = 0._r8 chm_mass = 0._r8 do m=1,nspec_amode(l) if ( index_tot_mass(l,m) > 0 ) & tot_mass = tot_mass + mmr(i,k,index_tot_mass(l,m)) if ( index_chm_mass(l,m) > 0 ) & chm_mass = chm_mass + mmr(i,k,index_chm_mass(l,m)) end do if ( tot_mass > 0._r8 ) then sad_mode(i,k,l) = chm_mass/tot_mass * & mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2*& exp(2*alnsg_amode(l)**2) ! m^2/m^3 sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3 else sad_mode(i,k,l) = 0._r8 end if end do sad(i,k) = sum(sad_mode(i,k,:)) enddo enddo if (present(sfc)) then sfc(:,:,:) = sad_mode(:,:,:) endif end subroutine surf_area_dens !=============================================================================== !=============================================================================== subroutine modal_aero_bcscavcoef_init !----------------------------------------------------------------------- ! ! Purpose: ! Computes lookup table for aerosol impaction/interception scavenging rates ! ! Authors: R. Easter ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use modal_aero_data use cam_abortutils, only: endrun implicit none ! local variables integer nnfit_maxd parameter (nnfit_maxd=27) integer i, jgrow, jdens, jpress, jtemp, mode, nnfit integer lunerr real(r8) dg0, dg0_cgs, press, & rhodryaero, rhowetaero, rhowetaero_cgs, rmserr, & scavratenum, scavratevol, sigmag, & temp, wetdiaratio, wetvolratio real(r8) aafitnum(1), xxfitnum(1,nnfit_maxd), yyfitnum(nnfit_maxd) real(r8) aafitvol(1), xxfitvol(1,nnfit_maxd), yyfitvol(nnfit_maxd) allocate(scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)) allocate(scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode)) lunerr = 6 dlndg_nimptblgrow = log( 1.25_r8 ) modeloop: do mode = 1, ntot_amode sigmag = sigmag_amode(mode) rhodryaero = specdens_amode(1,mode) growloop: do jgrow = nimptblgrow_mind, nimptblgrow_maxd wetdiaratio = exp( jgrow*dlndg_nimptblgrow ) dg0 = dgnum_amode(mode)*wetdiaratio wetvolratio = exp( jgrow*dlndg_nimptblgrow*3._r8 ) rhowetaero = 1.0_r8 + (rhodryaero-1.0_r8)/wetvolratio rhowetaero = min( rhowetaero, rhodryaero ) ! ! compute impaction scavenging rates at 1 temp-press pair and save ! nnfit = 0 temp = 273.16_r8 press = 0.75e6_r8 ! dynes/cm2 rhowetaero = rhodryaero dg0_cgs = dg0*1.0e2_r8 ! m to cm rhowetaero_cgs = rhowetaero*1.0e-3_r8 ! kg/m3 to g/cm3 call calc_1_impact_rate( & dg0_cgs, sigmag, rhowetaero_cgs, temp, press, & scavratenum, scavratevol, lunerr ) nnfit = nnfit + 1 if (nnfit .gt. nnfit_maxd) then write(lunerr,9110) call endrun() end if 9110 format( '*** subr. modal_aero_bcscavcoef_init -- nnfit too big' ) xxfitnum(1,nnfit) = 1._r8 yyfitnum(nnfit) = log( scavratenum ) xxfitvol(1,nnfit) = 1._r8 yyfitvol(nnfit) = log( scavratevol ) 5900 continue ! ! skip mlinfit stuff because scav table no longer has dependencies on ! air temp, air press, and particle wet density ! just load the log( scavrate--- ) values ! !! !! do linear regression !! log(scavrate) = a1 + a2*log(wetdens) !! ! call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 1, 1, rmserr ) ! call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 1, 1, rmserr ) ! ! scavimptblnum(jgrow,mode) = aafitnum(1) ! scavimptblvol(jgrow,mode) = aafitvol(1) scavimptblnum(jgrow,mode) = yyfitnum(1) scavimptblvol(jgrow,mode) = yyfitvol(1) enddo growloop enddo modeloop return end subroutine modal_aero_bcscavcoef_init !=============================================================================== !=============================================================================== subroutine modal_aero_depvel_part( ncol, t, pmid, ram1, fv, vlc_dry, vlc_trb, vlc_grv, & radius_part, density_part, sig_part, moment, lchnk ) ! calculates surface deposition velocity of particles ! L. Zhang, S. Gong, J. Padro, and L. Barrie ! A size-seggregated particle dry deposition scheme for an atmospheric aerosol module ! Atmospheric Environment, 35, 549-560, 2001. ! ! Authors: X. Liu ! ! !USES ! use physconst, only: pi,boltz, gravit, rair use mo_drydep, only: n_land_type, fraction_landuse ! !ARGUMENTS: ! implicit none ! real(r8), intent(in) :: t(pcols,pver) !atm temperature (K) real(r8), intent(in) :: pmid(pcols,pver) !atm pressure (Pa) real(r8), intent(in) :: fv(pcols) !friction velocity (m/s) real(r8), intent(in) :: ram1(pcols) !aerodynamical resistance (s/m) real(r8), intent(in) :: radius_part(pcols,pver) ! mean (volume/number) particle radius (m) real(r8), intent(in) :: density_part(pcols,pver) ! density of particle material (kg/m3) real(r8), intent(in) :: sig_part(pcols,pver) ! geometric standard deviation of particles integer, intent(in) :: moment ! moment of size distribution (0 for number, 2 for surface area, 3 for volume) integer, intent(in) :: ncol integer, intent(in) :: lchnk real(r8), intent(out) :: vlc_trb(pcols) !Turbulent deposn velocity (m/s) real(r8), intent(out) :: vlc_grv(pcols,pver) !grav deposn velocity (m/s) real(r8), intent(out) :: vlc_dry(pcols,pver) !dry deposn velocity (m/s) !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! Local Variables integer :: m,i,k,ix !indices real(r8) :: rho !atm density (kg/m**3) real(r8) :: vsc_dyn_atm(pcols,pver) ![kg m-1 s-1] Dynamic viscosity of air real(r8) :: vsc_knm_atm(pcols,pver) ![m2 s-1] Kinematic viscosity of atmosphere real(r8) :: shm_nbr ![frc] Schmidt number real(r8) :: stk_nbr ![frc] Stokes number real(r8) :: mfp_atm(pcols,pver) ![m] Mean free path of air real(r8) :: dff_aer ![m2 s-1] Brownian diffusivity of particle real(r8) :: slp_crc(pcols,pver) ![frc] Slip correction factor real(r8) :: rss_trb ![s m-1] Resistance to turbulent deposition real(r8) :: rss_lmn ![s m-1] Quasi-laminar layer resistance real(r8) :: brownian ! collection efficiency for Browning diffusion real(r8) :: impaction ! collection efficiency for impaction real(r8) :: interception ! collection efficiency for interception real(r8) :: stickfrac ! fraction of particles sticking to surface real(r8) :: radius_moment(pcols,pver) ! median radius (m) for moment real(r8) :: lnsig ! ln(sig_part) real(r8) :: dispersion ! accounts for influence of size dist dispersion on bulk settling velocity ! assuming radius_part is number mode radius * exp(1.5 ln(sigma)) integer :: lt real(r8) :: lnd_frc real(r8) :: wrk1, wrk2, wrk3 ! constants real(r8) gamma(11) ! exponent of schmidt number ! data gamma/0.54d+00, 0.56d+00, 0.57d+00, 0.54d+00, 0.54d+00, & ! 0.56d+00, 0.54d+00, 0.54d+00, 0.54d+00, 0.56d+00, & ! 0.50d+00/ data gamma/0.56e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.56e+00_r8, 0.56e+00_r8, & 0.56e+00_r8, 0.50e+00_r8, 0.54e+00_r8, 0.54e+00_r8, 0.54e+00_r8, & 0.54e+00_r8/ save gamma real(r8) alpha(11) ! parameter for impaction ! data alpha/50.00d+00, 0.95d+00, 0.80d+00, 1.20d+00, 1.30d+00, & ! 0.80d+00, 50.00d+00, 50.00d+00, 2.00d+00, 1.50d+00, & ! 100.00d+00/ data alpha/1.50e+00_r8, 1.20e+00_r8, 1.20e+00_r8, 0.80e+00_r8, 1.00e+00_r8, & 0.80e+00_r8, 100.00e+00_r8, 50.00e+00_r8, 2.00e+00_r8, 1.20e+00_r8, & 50.00e+00_r8/ save alpha real(r8) radius_collector(11) ! radius (m) of surface collectors ! data radius_collector/-1.00d+00, 5.10d-03, 3.50d-03, 3.20d-03, 10.00d-03, & ! 5.00d-03, -1.00d+00, -1.00d+00, 10.00d-03, 10.00d-03, & ! -1.00d+00/ data radius_collector/10.00e-03_r8, 3.50e-03_r8, 3.50e-03_r8, 5.10e-03_r8, 2.00e-03_r8, & 5.00e-03_r8, -1.00e+00_r8, -1.00e+00_r8, 10.00e-03_r8, 3.50e-03_r8, & -1.00e+00_r8/ save radius_collector integer :: iwet(11) ! flag for wet surface = 1, otherwise = -1 ! data iwet/1, -1, -1, -1, -1, & ! -1, -1, -1, 1, -1, & ! 1/ data iwet/-1, -1, -1, -1, -1, & -1, 1, -1, 1, -1, & -1/ save iwet !------------------------------------------------------------------------ do k=1,pver do i=1,ncol lnsig = log(sig_part(i,k)) ! use a maximum radius of 50 microns when calculating deposition velocity radius_moment(i,k) = min(50.0e-6_r8,radius_part(i,k))* & exp((float(moment)-1.5_r8)*lnsig*lnsig) dispersion = exp(2._r8*lnsig*lnsig) rho=pmid(i,k)/rair/t(i,k) ! Quasi-laminar layer resistance: call rss_lmn_get ! Size-independent thermokinetic properties vsc_dyn_atm(i,k) = 1.72e-5_r8 * ((t(i,k)/273.0_r8)**1.5_r8) * 393.0_r8 / & (t(i,k)+120.0_r8) ![kg m-1 s-1] RoY94 p. 102 mfp_atm(i,k) = 2.0_r8 * vsc_dyn_atm(i,k) / & ![m] SeP97 p. 455 (pmid(i,k)*sqrt(8.0_r8/(pi*rair*t(i,k)))) vsc_knm_atm(i,k) = vsc_dyn_atm(i,k) / rho ![m2 s-1] Kinematic viscosity of air slp_crc(i,k) = 1.0_r8 + mfp_atm(i,k) * & (1.257_r8+0.4_r8*exp(-1.1_r8*radius_moment(i,k)/(mfp_atm(i,k)))) / & radius_moment(i,k) ![frc] Slip correction factor SeP97 p. 464 vlc_grv(i,k) = (4.0_r8/18.0_r8) * radius_moment(i,k)*radius_moment(i,k)*density_part(i,k)* & gravit*slp_crc(i,k) / vsc_dyn_atm(i,k) ![m s-1] Stokes' settling velocity SeP97 p. 466 vlc_grv(i,k) = vlc_grv(i,k) * dispersion vlc_dry(i,k)=vlc_grv(i,k) enddo enddo k=pver ! only look at bottom level for next part do i=1,ncol dff_aer = boltz * t(i,k) * slp_crc(i,k) / & ![m2 s-1] (6.0_r8*pi*vsc_dyn_atm(i,k)*radius_moment(i,k)) !SeP97 p.474 shm_nbr = vsc_knm_atm(i,k) / dff_aer ![frc] SeP97 p.972 wrk2 = 0._r8 wrk3 = 0._r8 do lt = 1,n_land_type lnd_frc = fraction_landuse(i,lt,lchnk) if ( lnd_frc /= 0._r8 ) then brownian = shm_nbr**(-gamma(lt)) if (radius_collector(lt) > 0.0_r8) then ! vegetated surface stk_nbr = vlc_grv(i,k) * fv(i) / (gravit*radius_collector(lt)) interception = 2.0_r8*(radius_moment(i,k)/radius_collector(lt))**2.0_r8 else ! non-vegetated surface stk_nbr = vlc_grv(i,k) * fv(i) * fv(i) / (gravit*vsc_knm_atm(i,k)) ![frc] SeP97 p.965 interception = 0.0_r8 endif impaction = (stk_nbr/(alpha(lt)+stk_nbr))**2.0_r8 if (iwet(lt) > 0) then stickfrac = 1.0_r8 else stickfrac = exp(-sqrt(stk_nbr)) if (stickfrac < 1.0e-10_r8) stickfrac = 1.0e-10_r8 endif rss_lmn = 1.0_r8 / (3.0_r8 * fv(i) * stickfrac * (brownian+interception+impaction)) rss_trb = ram1(i) + rss_lmn + ram1(i)*rss_lmn*vlc_grv(i,k) wrk1 = 1.0_r8 / rss_trb wrk2 = wrk2 + lnd_frc*( wrk1 ) wrk3 = wrk3 + lnd_frc*( wrk1 + vlc_grv(i,k) ) endif enddo ! n_land_type vlc_trb(i) = wrk2 vlc_dry(i,k) = wrk3 enddo !ncol return end subroutine modal_aero_depvel_part !=============================================================================== subroutine modal_aero_bcscavcoef_get( m, ncol, isprx, dgn_awet, scavcoefnum, scavcoefvol ) use modal_aero_data !----------------------------------------------------------------------- implicit none integer,intent(in) :: m, ncol logical,intent(in):: isprx(pcols,pver) real(r8), intent(in) :: dgn_awet(pcols,pver,ntot_amode) real(r8), intent(out) :: scavcoefnum(pcols,pver), scavcoefvol(pcols,pver) integer i, k, jgrow real(r8) dumdgratio, xgrow, dumfhi, dumflo, scavimpvol, scavimpnum do k = 1, pver do i = 1, ncol ! do only if no precip if ( isprx(i,k) ) then ! ! interpolate table values using log of (actual-wet-size)/(base-dry-size) dumdgratio = dgn_awet(i,k,m)/dgnum_amode(m) if ((dumdgratio .ge. 0.99_r8) .and. (dumdgratio .le. 1.01_r8)) then scavimpvol = scavimptblvol(0,m) scavimpnum = scavimptblnum(0,m) else xgrow = log( dumdgratio ) / dlndg_nimptblgrow jgrow = int( xgrow ) if (xgrow .lt. 0._r8) jgrow = jgrow - 1 if (jgrow .lt. nimptblgrow_mind) then jgrow = nimptblgrow_mind xgrow = jgrow else jgrow = min( jgrow, nimptblgrow_maxd-1 ) end if dumfhi = xgrow - jgrow dumflo = 1._r8 - dumfhi scavimpvol = dumflo*scavimptblvol(jgrow,m) + & dumfhi*scavimptblvol(jgrow+1,m) scavimpnum = dumflo*scavimptblnum(jgrow,m) + & dumfhi*scavimptblnum(jgrow+1,m) end if ! impaction scavenging removal amount for volume scavcoefvol(i,k) = exp( scavimpvol ) ! impaction scavenging removal amount to number scavcoefnum(i,k) = exp( scavimpnum ) ! scavcoef = impaction scav rate (1/h) for precip = 1 mm/h ! scavcoef = impaction scav rate (1/s) for precip = pfx_inrain ! (scavcoef/3600) = impaction scav rate (1/s) for precip = 1 mm/h ! (pfx_inrain*3600) = in-rain-area precip rate (mm/h) ! impactrate = (scavcoef/3600) * (pfx_inrain*3600) else scavcoefvol(i,k) = 0._r8 scavcoefnum(i,k) = 0._r8 end if end do end do return end subroutine modal_aero_bcscavcoef_get !=============================================================================== subroutine calc_1_impact_rate( & dg0, sigmag, rhoaero, temp, press, & scavratenum, scavratevol, lunerr ) ! ! routine computes a single impaction scavenging rate ! for precipitation rate of 1 mm/h ! ! dg0 = geometric mean diameter of aerosol number size distrib. (cm) ! sigmag = geometric standard deviation of size distrib. ! rhoaero = density of aerosol particles (g/cm^3) ! temp = temperature (K) ! press = pressure (dyne/cm^2) ! scavratenum = number scavenging rate (1/h) ! scavratevol = volume or mass scavenging rate (1/h) ! lunerr = logical unit for error message ! use shr_kind_mod, only: r8 => shr_kind_r8 use mo_constants, only: boltz_cgs, pi, rhowater => rhoh2o_cgs, & gravity => gravity_cgs, rgas => rgas_cgs implicit none ! subr. parameters integer lunerr real(r8) dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol ! local variables integer nrainsvmax parameter (nrainsvmax=50) real(r8) rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),& vfallrainsv(nrainsvmax) integer naerosvmax parameter (naerosvmax=51) real(r8) aaerosv(naerosvmax), & ynumaerosv(naerosvmax), yvolaerosv(naerosvmax) integer i, ja, jr, na, nr real(r8) a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc real(r8) anumsum, avolsum, cair, chi real(r8) d, dr, dum, dumfuchs, dx real(r8) ebrown, eimpact, eintercept, etotal, freepath real(r8) precip, precipmmhr, precipsum real(r8) r, rainsweepout, reynolds, rhi, rhoair, rlo, rnumsum real(r8) scavsumnum, scavsumnumbb real(r8) scavsumvol, scavsumvolbb real(r8) schmidt, sqrtreynolds, sstar, stokes, sx real(r8) taurelax, vfall, vfallstp real(r8) x, xg0, xg3, xhi, xlo, xmuwaterair rlo = .005_r8 rhi = .250_r8 dr = 0.005_r8 nr = 1 + nint( (rhi-rlo)/dr ) if (nr .gt. nrainsvmax) then write(lunerr,9110) call endrun() end if 9110 format( '*** subr. calc_1_impact_rate -- nr > nrainsvmax' ) precipmmhr = 1.0_r8 precip = precipmmhr/36000._r8 ag0 = dg0/2._r8 sx = log( sigmag ) xg0 = log( ag0 ) xg3 = xg0 + 3._r8*sx*sx xlo = xg3 - 4._r8*sx xhi = xg3 + 4._r8*sx dx = 0.2_r8*sx dx = max( 0.2_r8*sx, 0.01_r8 ) xlo = xg3 - max( 4._r8*sx, 2._r8*dx ) xhi = xg3 + max( 4._r8*sx, 2._r8*dx ) na = 1 + nint( (xhi-xlo)/dx ) if (na .gt. naerosvmax) then write(lunerr,9120) call endrun() end if 9120 format( '*** subr. calc_1_impact_rate -- na > naerosvmax' ) ! air molar density cair = press/(rgas*temp) ! air mass density rhoair = 28.966_r8*cair ! molecular freepath freepath = 2.8052e-10_r8/cair ! air dynamic viscosity airdynvisc = 1.8325e-4_r8 * (416.16_r8/(temp+120._r8)) * & ((temp/296.16_r8)**1.5_r8) ! air kinemaic viscosity airkinvisc = airdynvisc/rhoair ! ratio of water viscosity to air viscosity (from Slinn) xmuwaterair = 60.0_r8 ! ! compute rain drop number concentrations ! rrainsv = raindrop radius (cm) ! xnumrainsv = raindrop number concentration (#/cm^3) ! (number in the bin, not number density) ! vfallrainsv = fall velocity (cm/s) ! precipsum = 0._r8 do i = 1, nr r = rlo + (i-1)*dr rrainsv(i) = r xnumrainsv(i) = exp( -r/2.7e-2_r8 ) d = 2._r8*r if (d .le. 0.007_r8) then vfallstp = 2.88e5_r8 * d**2._r8 else if (d .le. 0.025_r8) then vfallstp = 2.8008e4_r8 * d**1.528_r8 else if (d .le. 0.1_r8) then vfallstp = 4104.9_r8 * d**1.008_r8 else if (d .le. 0.25_r8) then vfallstp = 1812.1_r8 * d**0.638_r8 else vfallstp = 1069.8_r8 * d**0.235_r8 end if vfall = vfallstp * sqrt(1.204e-3_r8/rhoair) vfallrainsv(i) = vfall precipsum = precipsum + vfall*(r**3)*xnumrainsv(i) end do precipsum = precipsum*pi*1.333333_r8 rnumsum = 0._r8 do i = 1, nr xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum) rnumsum = rnumsum + xnumrainsv(i) end do ! ! compute aerosol concentrations ! aaerosv = particle radius (cm) ! fnumaerosv = fraction of total number in the bin (--) ! fvolaerosv = fraction of total volume in the bin (--) ! anumsum = 0._r8 avolsum = 0._r8 do i = 1, na x = xlo + (i-1)*dx a = exp( x ) aaerosv(i) = a dum = (x - xg0)/sx ynumaerosv(i) = exp( -0.5_r8*dum*dum ) yvolaerosv(i) = ynumaerosv(i)*1.3333_r8*pi*a*a*a anumsum = anumsum + ynumaerosv(i) avolsum = avolsum + yvolaerosv(i) end do do i = 1, na ynumaerosv(i) = ynumaerosv(i)/anumsum yvolaerosv(i) = yvolaerosv(i)/avolsum end do ! ! compute scavenging ! scavsumnum = 0._r8 scavsumvol = 0._r8 ! ! outer loop for rain drop radius ! jr_loop: do jr = 1, nr r = rrainsv(jr) vfall = vfallrainsv(jr) reynolds = r * vfall / airkinvisc sqrtreynolds = sqrt( reynolds ) ! ! inner loop for aerosol particle radius ! scavsumnumbb = 0._r8 scavsumvolbb = 0._r8 ja_loop: do ja = 1, na a = aaerosv(ja) chi = a/r dum = freepath/a dumfuchs = 1._r8 + 1.246_r8*dum + 0.42_r8*dum*exp(-0.87_r8/dum) taurelax = 2._r8*rhoaero*a*a*dumfuchs/(9._r8*rhoair*airkinvisc) aeromass = 4._r8*pi*a*a*a*rhoaero/3._r8 aerodiffus = boltz_cgs*temp*taurelax/aeromass schmidt = airkinvisc/aerodiffus stokes = vfall*taurelax/r ebrown = 4._r8*(1._r8 + 0.4_r8*sqrtreynolds*(schmidt**0.3333333_r8)) / & (reynolds*schmidt) dum = (1._r8 + 2._r8*xmuwaterair*chi) / & (1._r8 + xmuwaterair/sqrtreynolds) eintercept = 4._r8*chi*(chi + dum) dum = log( 1._r8 + reynolds ) sstar = (1.2_r8 + dum/12._r8) / (1._r8 + dum) eimpact = 0._r8 if (stokes .gt. sstar) then dum = stokes - sstar eimpact = (dum/(dum+0.6666667_r8)) ** 1.5_r8 end if etotal = ebrown + eintercept + eimpact etotal = min( etotal, 1.0_r8 ) rainsweepout = xnumrainsv(jr)*4._r8*pi*r*r*vfall scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja) scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja) enddo ja_loop scavsumnum = scavsumnum + scavsumnumbb scavsumvol = scavsumvol + scavsumvolbb enddo jr_loop scavratenum = scavsumnum*3600._r8 scavratevol = scavsumvol*3600._r8 return end subroutine calc_1_impact_rate !============================================================================= !============================================================================= subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf) use modal_aero_data, only : qqcw_get_field use physics_buffer, only : physics_buffer_desc !----------------------------------------------------------------- ! ... Xfrom from mass to volume mixing ratio !----------------------------------------------------------------- use chem_mods, only : adv_mass, gas_pcnst implicit none !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- integer, intent(in) :: lchnk, ncol, im real(r8), intent(in) :: mbar(ncol,pver) real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst) type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- integer :: k, m real(r8), pointer :: fldcw(:,:) do m=1,gas_pcnst if( adv_mass(m) /= 0._r8 ) then fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) if(associated(fldcw)) then do k=1,pver vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m) end do else vmr(:,:,m) = 0.0_r8 end if end if end do end subroutine qqcw2vmr !============================================================================= !============================================================================= subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf ) !----------------------------------------------------------------- ! ... Xfrom from volume to mass mixing ratio !----------------------------------------------------------------- use m_spc_id use chem_mods, only : adv_mass, gas_pcnst use modal_aero_data, only : qqcw_get_field use physics_buffer, only : physics_buffer_desc implicit none !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- integer, intent(in) :: lchnk, ncol, im real(r8), intent(in) :: mbar(ncol,pver) real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst) type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- integer :: k, m real(r8), pointer :: fldcw(:,:) !----------------------------------------------------------------- ! ... The non-group species !----------------------------------------------------------------- do m = 1,gas_pcnst fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then do k = 1,pver fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k) end do end if end do end subroutine vmr2qqcw function get_dlndg_nimptblgrow() result (dlndg_nimptblgrow_ret) real(r8) :: dlndg_nimptblgrow_ret dlndg_nimptblgrow_ret = dlndg_nimptblgrow end function get_dlndg_nimptblgrow function get_scavimptblvol() result (scavimptblvol_ret) real(r8) :: scavimptblvol_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode) scavimptblvol_ret = scavimptblvol end function get_scavimptblvol function get_scavimptblnum() result (scavimptblnum_ret) real(r8) :: scavimptblnum_ret(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode) scavimptblnum_ret = scavimptblnum end function get_scavimptblnum end module aero_model