!=============================================================================== ! Modal Aerosol Model ! JUly 2015 B.Singh Added unified convection code !=============================================================================== 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 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, spec_class_gas use spmd_utils, only: masterproc 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 use modal_aero_data,only: ntot_amode, modename_amode 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 ! aerosol surface area for chemistry ! These are made public to be used by MMF w/ ECPP public :: calc_1_impact_rate public :: dlndg_nimptblgrow public :: nimptblgrow_mind public :: nimptblgrow_maxd public :: scavimptblnum public :: scavimptblvol ! 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 :: nevapr_idx = 0 integer :: rprddp_idx = 0 integer :: rprdsh_idx = 0 integer :: nevapr_shcu_idx = 0 integer :: nevapr_dpcu_idx = 0 integer :: icwmrdp_idx = 0 integer :: icwmrsh_idx = 0 integer :: sh_frac_idx = 0 integer :: dp_frac_idx = 0 integer :: imozart = -1 logical :: history_aero_prevap_resusp = .false. ! controls output of prevap resusp tendencies to history ! variables for table lookup of aerosol impaction/interception scavenging rates integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12 real(r8) :: dlndg_nimptblgrow real(r8) :: scavimptblnum(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode) real(r8) :: scavimptblvol(nimptblgrow_mind:nimptblgrow_maxd, ntot_amode) ! for aero_model_surfarea called from mo_usrrxt integer :: aitken_idx = -1 integer, dimension(ntot_amode) :: num_idx = -1 integer :: index_tot_mass(ntot_amode,10) = -1 integer :: index_chm_mass(ntot_amode,10) = -1 integer :: ndx_h2so4 character(len=fieldname_len) :: dgnum_name(ntot_amode) character(len=fieldname_len) :: scavcoef_name(ntot_amode) !For aero_model_wetdep subroutine integer :: strt_loop, end_loop, stride_loop !loop indices for the lphase loop ! Namelist variables integer :: mam_amicphys_optaa logical :: sscav_tuning, convproc_do_aer, convproc_do_gas, resus_fix logical :: get_presc_aero_data 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) 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, seasalt_emis_scale, sscav_tuning, & sol_factb_interstitial, sol_factic_interstitial !----------------------------------------------------------------------------- ! 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(sscav_tuning, 1, mpilog, 0, mpicom) call mpibcast(seasalt_emis_scale, 1, mpir8, 0, mpicom) #endif wetdep_list = aer_wetdep_list drydep_list = aer_drydep_list end subroutine aero_model_readnl !============================================================================= !============================================================================= subroutine aero_model_register(imozart_in, species_class) use modal_aero_initialize_data, only : modal_aero_register integer, intent(in) :: imozart_in integer, intent(inout) :: species_class(:) imozart = imozart_in call modal_aero_register(species_class) end subroutine aero_model_register !============================================================================= !============================================================================= subroutine aero_model_init( pbuf2d, species_class, iflagaa ) use mo_chem_utls, only: get_inv_ndx use cam_history, only: addfld, horiz_only, add_default use phys_control, only: phys_getopts use mo_chem_utls, only: get_rxt_ndx, get_spc_ndx use modal_aero_data, only: cnst_name_cw, rain_evap_to_coarse_aero, mam_prevap_resusp_optaa use modal_aero_initialize_data, only: modal_aero_initialize use modal_aero_convproc, only: deepconv_wetdep_history 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 mo_chem_utls, only: get_het_ndx use gas_wetdep_opts, only: gas_wetdep_cnt, gas_wetdep_list, gas_wetdep_method ! REASTER 08/04/2015 ! args type(physics_buffer_desc), pointer :: pbuf2d(:,:) integer, intent(inout) :: species_class(:) integer, intent(in) :: iflagaa ! local vars integer :: id, l, m, n, nspc logical :: history_aerosol ! Output MAM or SECT aerosol tendencies logical :: history_verbose ! produce verbose history output character(len=*), parameter :: subrname = 'aero_model_init' character(len=20) :: dummy character(len=fieldname_len) :: wetdep_name, depflx_name character(len=6) :: test_name character(len=100) :: errmes character(len=2) :: unit_basename ! Units 'kg' or '1' if ( masterproc ) write(iulog,'(a,i5)') 'aero_model_init iflagaa=', iflagaa ! REASTER 08/04/2015 call phys_getopts( history_aerosol_out=history_aerosol, & history_verbose_out=history_verbose, & convproc_do_aer_out = convproc_do_aer, & convproc_do_gas_out = convproc_do_gas, & resus_fix_out = resus_fix, & get_presc_aero_data_out = get_presc_aero_data, & mam_amicphys_optaa_out = mam_amicphys_optaa ) ! REASTER 08/04/2015 ! REASTER 08/04/2015 BEGIN ! This section cannot execute until chemini, ..., chm_diags_inti have been called if ( iflagaa == 2 ) then if ( masterproc ) then write(iulog,'(a,i5,2x,a)') 'gas_wetdep_cnt,meth', gas_wetdep_cnt, gas_wetdep_method do m = 1, gas_wetdep_cnt write(iulog,'(a,i5,2x,a)') 'gas_wetdep_list ', m, trim(gas_wetdep_list(m)) end do end if ! These WD_ and DF_ fields should always been in a MAM history file, ! but for now they are conditional on convproc_do_aer if ( convproc_do_aer ) then do m = 1,gas_pcnst call cnst_get_ind( solsym(m), l, abrtf=.false. ) if ( ( history_aerosol ) .and. (l > 0) ) then if ( species_class(l) == spec_class_gas ) then !RCE - only output WD_xxx and DF_xxx for gases wetdep_name = 'WD_'//trim(solsym(m)) depflx_name = 'DF_'//trim(solsym(m)) nspc = get_het_ndx(solsym(m)) if (nspc > 0) call add_default( wetdep_name, 1, ' ' ) call add_default( depflx_name, 1, ' ' ) endif endif end do ! m = 1,gas_pcnst endif return endif ! ( iflagaa == 2 ) m = mam_prevap_resusp_optaa if ( rain_evap_to_coarse_aero ) then if ( mam_prevap_resusp_optaa /= 30 ) then mam_prevap_resusp_optaa = 30 if ( masterproc ) write(iulog,'(2a,i4,a)') 'aero_model_init - ', & 'mam_prevap_resusp_optaa changed from ', m, ' to 30 because rain_evap_to_coarse_aero = .true.' end if else if ( mam_prevap_resusp_optaa == 10 .and. resus_fix ) then ! this case is mam_prevap_resusp_optaa = 20, so change the default value mam_prevap_resusp_optaa = 20 if ( masterproc ) write(iulog,'(2a)') 'aero_model_init - ', & 'mam_prevap_resusp_optaa changed from 10 to 20 because resus_fix = .true.' else if ( mam_prevap_resusp_optaa == 10 .and. convproc_do_aer ) then ! when convproc_do_aer, the prevap_resusp tendencies are always written ! to history, so change mam_prevap_resusp_optaa to 11 mam_prevap_resusp_optaa = 11 if ( masterproc ) write(iulog,'(2a)') 'aero_model_init - ', & 'mam_prevap_resusp_optaa changed from 10 to 11 because convproc_do_aer = .true.' endif endif ! *** activate this to override mam_prevap_resusp_optaa value for testing *** ! mam_prevap_resusp_optaa = 20 ! if ( masterproc ) write(iulog,'(2a)') 'aero_model_init - ', & ! 'mam_prevap_resusp_optaa changed from ?? to 20 for special test run' if ( masterproc ) then write(iulog,'(2a,4l5,2i5)') 'aero_model_init - ', & 'convproc_do_aer & _gas, resus_fix, rain_evap_to_coarse, mam_prevap_resusp_optaa_v1/v2', & convproc_do_aer, convproc_do_gas, resus_fix, & rain_evap_to_coarse_aero, m, mam_prevap_resusp_optaa endif ! The unified convective transport/removal for aerosols does not ! do gases yet, and convproc_do_gas is just a place holder. For that reason, ! (1) All of the "if ( convproc_do_aer .or. convproc_do_gas ) then" statements ! in aero_model.F90 have been changed to "if ( convproc_do_aer ) then" ! (2) convproc_do_aer=.false. and convproc_do_gas=.true. is no longer allowed. if ( ( .not. convproc_do_aer ) .and. convproc_do_gas ) then errmes = 'aero_model_init - ' // & 'convproc_do_aer MUST BE .true. when convproc_do_gas is .true.' call endrun( errmes ) endif if (masterproc) then if ( convproc_do_aer .and. ( .not. resus_fix ) ) then write(iulog,*)'WARNING: resus_fix=.false. and convproc_do_aer=.true.' //& ' is not a well tested configuration, may produce incorrect results!!' endif endif ! mam_prevap_resusp_optaa values ! 0 = no resuspension ! 10 = original mam method with resus_fix=.false. (so4_a1 --> so4_a1, so4_c1 --> so4_c1) ! 20 = original mam method with resus_fix=.true. (so4_a1 & so4_c1 --> so4_a1) ! 30 = resuspend to coarse mode, full non-linear method (so4_a1 & so4_c1 --> so4_a3) ! 11 = like 10 but output column resuspension tendencies (rcscavt & rsscavt) to history ! 21 = like 20 but a with a few xxx = max( 0.0, xxx) added in werdepa_v2 m = 0 if ( mam_prevap_resusp_optaa == 0 ) m = 1 if ( mam_prevap_resusp_optaa == 10 ) m = 1 if ( mam_prevap_resusp_optaa == 11 ) m = 1 if ( mam_prevap_resusp_optaa == 20 ) m = 1 if ( mam_prevap_resusp_optaa == 21 ) m = 1 if ( rain_evap_to_coarse_aero ) then if ( mam_prevap_resusp_optaa == 30 ) m = 1 endif if (m <= 0) then write(errmes,'(2a,l5,i10)') 'aero_model_init - ', & 'illegal rain_evap_to_coarse_aero, bad mam_prevap_resusp_optaa = ', & rain_evap_to_coarse_aero, mam_prevap_resusp_optaa call endrun( errmes ) endif history_aero_prevap_resusp = .false. if ( mam_prevap_resusp_optaa /= 10 ) history_aero_prevap_resusp = .true. ! REASTER 08/04/2015 END dgnum_idx = pbuf_get_index('DGNUM') dgnumwet_idx = pbuf_get_index('DGNUMWET') !BSINGH: Decide the loop counters for the lphase loop in aero_model_wetdep subroutine !for cases with and without the unified convective transport !Counters for "without" unified convective treatment (i.e. default case) strt_loop = 1 end_loop = 2 stride_loop = 1 if (convproc_do_aer) then !BSINGH (09/12/2014):Do cloudborne first for unified convection scheme so that the resuspension of cloudborne !can be saved then applied to interstitial (RCE) strt_loop = 2 end_loop = 1 stride_loop = -1 endif call rad_cnst_get_info(0, nmodes=nmodes) call modal_aero_initialize(pbuf2d, imozart, species_class) call modal_aero_bcscavcoef_init() call mam_prevap_resusp_init( ) ! REASTER 08/04/2015 call dust_init() call seasalt_init() call wetdep_init() fracis_idx = pbuf_get_index('FRACIS') prain_idx = pbuf_get_index('PRAIN') nevapr_idx = pbuf_get_index('NEVAPR') 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') icwmrdp_idx = pbuf_get_index('ICWMRDP') icwmrsh_idx = pbuf_get_index('ICWMRSH') sh_frac_idx = pbuf_get_index('SH_FRAC') dp_frac_idx = pbuf_get_index('DP_FRAC') 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, abrtf=.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, abrtf=.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','1','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 #if (defined MODAL_AERO_9MODE || MODAL_AERO_4MODE_MOM || defined MODAL_AERO_5MODE_AGEDCARBON) dummy = 'SSTSFMBL_OM' call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux of marine organic matter at surface') if (history_aerosol) then call add_default (dummy, 1, ' ') endif dummy = 'F_eff' call addfld (dummy,horiz_only, 'A','1','Effective enrichment factor of marine organic matter') if (history_aerosol) then call add_default (dummy, 1, ' ') endif #endif 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, ' ') if ( history_verbose ) then call add_default (trim(drydep_list(m))//'TBF', 1, ' ') call add_default (trim(drydep_list(m))//'GVF', 1, ' ') endif endif enddo do m = 1,nwetdep if ( masterproc ) write(iulog,'(a,i3,2x,a)') 'm, wetdep_list', m, trim(wetdep_list(m)) ! REASTER 08/04/2015 ! 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 ( history_aero_prevap_resusp ) then call addfld (trim(wetdep_list(m))//'SFSEC', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, convective) at surface') !RCE call addfld (trim(wetdep_list(m))//'SFSES', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, stratiform) at surface') !RCE if (convproc_do_aer .and. deepconv_wetdep_history) then call addfld (trim(wetdep_list(m))//'SFSED', & horiz_only, 'A','kg/m2/s','Wet deposition flux (precip evap, deep convective) at surface') !RCE endif endif if (convproc_do_aer .and. deepconv_wetdep_history) then call addfld (trim(wetdep_list(m))//'SFSID', & horiz_only, 'A','kg/m2/s','Wet deposition flux (incloud, deep convective) at surface') !RCE call addfld (trim(wetdep_list(m))//'SFSBD', & horiz_only, 'A','kg/m2/s','Wet deposition flux (belowcloud, deep convective) at surface') !RCE endif 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, ' ') if ( history_verbose ) then 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 ( history_aero_prevap_resusp ) then call add_default (trim(wetdep_list(m))//'SFSEC', 1, ' ') call add_default (trim(wetdep_list(m))//'SFSES', 1, ' ') endif endif endif enddo ! m = 1,nwetdep 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 if ( history_verbose ) then call add_default( 'GS_'//trim(solsym(m)), 1, ' ') call add_default( 'AQ_'//trim(solsym(m)), 1, ' ') else select case (trim(solsym(m))) case ('O3','H2O2','H2SO4','SO2','DMS','SOAG') call add_default( 'AQ_'//trim(solsym(m)), 1, ' ') end select end if endif call cnst_get_ind(trim(solsym(m)), nspc, abrtf=.false. ) ! REASTER 08/04/2015 ! if(nspc > 0 .and. .not.cnst_name_cw(nspc) == ' ') then ! REASTER 08/04/2015 if( nspc > 0 ) then ! REASTER 08/04/2015 if ( .not. cnst_name_cw(nspc) == ' ') then ! REASTER 08/04/2015 if ( history_aero_prevap_resusp ) then call addfld (trim(cnst_name_cw(nspc))//'SFSEC',horiz_only, 'A','kg/m2/s', & trim(cnst_name_cw(nspc))//' wet deposition flux (precip evap, convective) at surface') !RCE call addfld (trim(cnst_name_cw(nspc))//'SFSES',horiz_only, 'A','kg/m2/s', & trim(cnst_name_cw(nspc))//' wet deposition flux (precip evap, stratiform) at surface') !RCE if(history_aerosol .and. history_verbose) then call add_default (trim(cnst_name_cw(nspc))//'SFSEC', 1, ' ') !RCE call add_default (trim(cnst_name_cw(nspc))//'SFSES', 1, ' ') !RCE endif endif endif 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 ( history_aerosol ) then if (history_verbose) then call add_default (trim(cnst_name_cw(n))//'GVF', 1, ' ') call add_default (trim(cnst_name_cw(n))//'TBF', 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, ' ') endif if (get_presc_aero_data .or. history_verbose) then call add_default( cnst_name_cw(n), 1, ' ' ) endif call add_default (trim(cnst_name_cw(n))//'SFWET', 1, ' ') call add_default (trim(cnst_name_cw(n))//'DDF', 1, ' ') endif endif enddo do n=1,ntot_amode dgnum_name(n) = ' ' write(dgnum_name(n),fmt='(a,i1)') 'dgnumwet',n call addfld( dgnum_name(n), (/ 'lev' /), 'I', 'm', 'Aerosol mode wet diameter' ) if ( history_aerosol .and. history_verbose ) then call add_default( dgnum_name(n), 1, ' ' ) endif scavcoef_name(n) = ' ' write(scavcoef_name(n),fmt='(a,i1)') 'scavcoef',n call addfld( scavcoef_name(n), (/ 'lev' /), 'I', '-', 'Aerosol mode scavenging coefficient' ) end do ndx_h2so4 = get_spc_ndx('H2SO4') ! for aero_model_surfarea called from mo_usrrxt do l=1,ntot_amode if ( trim(modename_amode(l)) == 'aitken' ) then aitken_idx = l end if 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 dgnumwet_idx = pbuf_get_index('DGNUMWET') if ( aitken_idx < 0 ) then errmes = 'usrrxt_inti: cannot find aitken_idx' call endrun(errmes) end if ! ! define indeces associated with the various names (defined in ! chemistry/modal_aero/modal_aero_initialize_data.F90) ! #if ( defined MODAL_AERO_3MODE ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('pom_a1') index_tot_mass(1,3) = get_spc_ndx('soa_a1') index_tot_mass(1,4) = get_spc_ndx('bc_a1' ) index_tot_mass(1,5) = get_spc_ndx('dst_a1') index_tot_mass(1,6) = get_spc_ndx('ncl_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('soa_a1') index_chm_mass(1,3) = get_spc_ndx('bc_a1' ) ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('soa_a2') index_tot_mass(2,3) = get_spc_ndx('ncl_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('soa_a2') ! ! coarse mode ! index_tot_mass(3,1) = get_spc_ndx('dst_a3') index_tot_mass(3,2) = get_spc_ndx('ncl_a3') index_tot_mass(3,3) = get_spc_ndx('so4_a3') index_chm_mass(3,1) = get_spc_ndx('so4_a3') ! #elif ( defined MODAL_AERO_4MODE_MOM ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('pom_a1') index_tot_mass(1,3) = get_spc_ndx('soa_a1') index_tot_mass(1,4) = get_spc_ndx('bc_a1' ) index_tot_mass(1,5) = get_spc_ndx('dst_a1') index_tot_mass(1,6) = get_spc_ndx('ncl_a1') index_tot_mass(1,7) = get_spc_ndx('mom_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('soa_a1') index_chm_mass(1,3) = get_spc_ndx('bc_a1' ) ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('soa_a2') index_tot_mass(2,3) = get_spc_ndx('ncl_a2') index_tot_mass(2,4) = get_spc_ndx('mom_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('soa_a2') ! ! coarse mode ! index_tot_mass(3,1) = get_spc_ndx('dst_a3') index_tot_mass(3,2) = get_spc_ndx('ncl_a3') index_tot_mass(3,3) = get_spc_ndx('so4_a3') index_chm_mass(3,1) = get_spc_ndx('so4_a3') ! ! POM mode ! index_tot_mass(4,1) = get_spc_ndx('pom_a4') index_tot_mass(4,2) = get_spc_ndx('bc_a4') index_tot_mass(4,3) = get_spc_ndx('mom_a4') index_chm_mass(4,1) = get_spc_ndx('bc_a1' ) ! #elif ( defined MODAL_AERO_4MODE ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('pom_a1') index_tot_mass(1,3) = get_spc_ndx('soa_a1') index_tot_mass(1,4) = get_spc_ndx('bc_a1' ) index_tot_mass(1,5) = get_spc_ndx('dst_a1') index_tot_mass(1,6) = get_spc_ndx('ncl_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('soa_a1') index_chm_mass(1,3) = get_spc_ndx('bc_a1' ) ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('soa_a2') index_tot_mass(2,3) = get_spc_ndx('ncl_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('soa_a2') ! ! coarse mode ! index_tot_mass(3,1) = get_spc_ndx('dst_a3') index_tot_mass(3,2) = get_spc_ndx('ncl_a3') index_tot_mass(3,3) = get_spc_ndx('so4_a3') index_chm_mass(3,1) = get_spc_ndx('so4_a3') ! ! POM mode ! index_tot_mass(4,1) = get_spc_ndx('pom_a4') index_tot_mass(4,2) = get_spc_ndx('bc_a4') index_chm_mass(4,1) = get_spc_ndx('bc_a1' ) ! #elif ( defined MODAL_AERO_5MODE_AGEDCARBON ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('soa_a1') index_tot_mass(1,3) = get_spc_ndx('dst_a1') index_tot_mass(1,4) = get_spc_ndx('ncl_a1') index_tot_mass(1,5) = get_spc_ndx('mom_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('soa_a1') ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('soa_a2') index_tot_mass(2,3) = get_spc_ndx('ncl_a2') index_tot_mass(2,4) = get_spc_ndx('mom_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('soa_a2') ! ! coarse mode ! index_tot_mass(3,1) = get_spc_ndx('dst_a3') index_tot_mass(3,2) = get_spc_ndx('ncl_a3') index_tot_mass(3,3) = get_spc_ndx('so4_a3') index_tot_mass(3,4) = get_spc_ndx('bc_a3') index_tot_mass(3,5) = get_spc_ndx('pom_a3') index_tot_mass(3,6) = get_spc_ndx('soa_a3') index_tot_mass(3,7) = get_spc_ndx('mom_a3') index_chm_mass(3,1) = get_spc_ndx('so4_a3') index_chm_mass(3,2) = get_spc_ndx('soa_a3') index_chm_mass(3,3) = get_spc_ndx('bc_a3' ) ! ! POM mode ! index_tot_mass(4,1) = get_spc_ndx('pom_a4') index_tot_mass(4,2) = get_spc_ndx('bc_a4') index_tot_mass(4,3) = get_spc_ndx('mom_a4') index_chm_mass(4,1) = get_spc_ndx('bc_a4' ) ! ! aged carbon mode ! index_tot_mass(5,1) = get_spc_ndx('so4_a5') index_tot_mass(5,2) = get_spc_ndx('pom_a5') index_tot_mass(5,3) = get_spc_ndx('soa_a5') index_tot_mass(5,4) = get_spc_ndx('bc_a5' ) index_tot_mass(5,5) = get_spc_ndx('ncl_a5') index_tot_mass(5,6) = get_spc_ndx('mom_a5') index_chm_mass(5,1) = get_spc_ndx('so4_a5') index_chm_mass(5,2) = get_spc_ndx('soa_a5') index_chm_mass(5,3) = get_spc_ndx('bc_a5' ) ! #elif ( defined MODAL_AERO_7MODE ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('nh4_a1') index_tot_mass(1,3) = get_spc_ndx('pom_a1') index_tot_mass(1,4) = get_spc_ndx('soa_a1') index_tot_mass(1,5) = get_spc_ndx('bc_a1' ) index_tot_mass(1,6) = get_spc_ndx('ncl_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('nh4_a1') index_chm_mass(1,3) = get_spc_ndx('soa_a1') index_chm_mass(1,4) = get_spc_ndx('bc_a1' ) ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('nh4_a2') index_tot_mass(2,3) = get_spc_ndx('soa_a2') index_tot_mass(2,4) = get_spc_ndx('ncl_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('nh4_a2') index_chm_mass(2,3) = get_spc_ndx('soa_a2') ! ! primary carbon mode not added ! ! fine sea salt ! index_tot_mass(4,1) = get_spc_ndx('so4_a4') index_tot_mass(4,2) = get_spc_ndx('nh4_a4') index_tot_mass(4,3) = get_spc_ndx('ncl_a4') index_chm_mass(4,1) = get_spc_ndx('so4_a4') index_chm_mass(4,2) = get_spc_ndx('nh4_a4') ! ! fine soil dust ! index_tot_mass(5,1) = get_spc_ndx('so4_a5') index_tot_mass(5,2) = get_spc_ndx('nh4_a5') index_tot_mass(5,3) = get_spc_ndx('dst_a5') index_chm_mass(5,1) = get_spc_ndx('so4_a5') index_chm_mass(5,2) = get_spc_ndx('nh4_a5') ! ! coarse sea salt ! index_tot_mass(6,1) = get_spc_ndx('so4_a6') index_tot_mass(6,2) = get_spc_ndx('nh4_a6') index_tot_mass(6,3) = get_spc_ndx('ncl_a6') index_chm_mass(6,1) = get_spc_ndx('so4_a6') index_chm_mass(6,2) = get_spc_ndx('nh4_a6') ! ! coarse soil dust ! index_tot_mass(7,1) = get_spc_ndx('so4_a7') index_tot_mass(7,2) = get_spc_ndx('nh4_a7') index_tot_mass(7,3) = get_spc_ndx('dst_a7') index_chm_mass(7,1) = get_spc_ndx('so4_a7') index_chm_mass(7,2) = get_spc_ndx('nh4_a7') ! #elif ( defined MODAL_AERO_9MODE ) ! ! accumulation mode #1 ! index_tot_mass(1,1) = get_spc_ndx('so4_a1') index_tot_mass(1,2) = get_spc_ndx('nh4_a1') index_tot_mass(1,3) = get_spc_ndx('pom_a1') index_tot_mass(1,4) = get_spc_ndx('soa_a1') index_tot_mass(1,5) = get_spc_ndx('bc_a1' ) index_tot_mass(1,6) = get_spc_ndx('ncl_a1') index_tot_mass(1,7) = get_spc_ndx('mpoly_a1') index_tot_mass(1,8) = get_spc_ndx('mprot_a1') index_tot_mass(1,9) = get_spc_ndx('mlip_a1') index_chm_mass(1,1) = get_spc_ndx('so4_a1') index_chm_mass(1,2) = get_spc_ndx('nh4_a1') index_chm_mass(1,3) = get_spc_ndx('soa_a1') index_chm_mass(1,4) = get_spc_ndx('bc_a1' ) ! ! aitken mode ! index_tot_mass(2,1) = get_spc_ndx('so4_a2') index_tot_mass(2,2) = get_spc_ndx('nh4_a2') index_tot_mass(2,3) = get_spc_ndx('soa_a2') index_tot_mass(2,4) = get_spc_ndx('ncl_a2') index_tot_mass(2,5) = get_spc_ndx('mpoly_a2') index_tot_mass(2,6) = get_spc_ndx('mprot_a2') index_tot_mass(2,7) = get_spc_ndx('mlip_a2') index_chm_mass(2,1) = get_spc_ndx('so4_a2') index_chm_mass(2,2) = get_spc_ndx('nh4_a2') index_chm_mass(2,3) = get_spc_ndx('soa_a2') ! ! primary carbon mode not added ! ! fine sea salt ! index_tot_mass(4,1) = get_spc_ndx('so4_a4') index_tot_mass(4,2) = get_spc_ndx('nh4_a4') index_tot_mass(4,3) = get_spc_ndx('ncl_a4') index_tot_mass(4,4) = get_spc_ndx('mpoly_a4') index_tot_mass(4,5) = get_spc_ndx('mprot_a4') index_tot_mass(4,6) = get_spc_ndx('mlip_a4') index_chm_mass(4,1) = get_spc_ndx('so4_a4') index_chm_mass(4,2) = get_spc_ndx('nh4_a4') ! ! fine soil dust ! index_tot_mass(5,1) = get_spc_ndx('so4_a5') index_tot_mass(5,2) = get_spc_ndx('nh4_a5') index_tot_mass(5,3) = get_spc_ndx('dst_a5') index_chm_mass(5,1) = get_spc_ndx('so4_a5') index_chm_mass(5,2) = get_spc_ndx('nh4_a5') ! ! coarse sea salt ! index_tot_mass(6,1) = get_spc_ndx('so4_a6') index_tot_mass(6,2) = get_spc_ndx('nh4_a6') index_tot_mass(6,3) = get_spc_ndx('ncl_a6') index_chm_mass(6,1) = get_spc_ndx('so4_a6') index_chm_mass(6,2) = get_spc_ndx('nh4_a6') ! ! coarse soil dust ! index_tot_mass(7,1) = get_spc_ndx('so4_a7') index_tot_mass(7,2) = get_spc_ndx('nh4_a7') index_tot_mass(7,3) = get_spc_ndx('dst_a7') index_chm_mass(7,1) = get_spc_ndx('so4_a7') index_chm_mass(7,2) = get_spc_ndx('nh4_a7') ! ! marine organics - accumulation marine ! index_tot_mass(8,1) = get_spc_ndx('mpoly_a8') index_tot_mass(8,2) = get_spc_ndx('mprot_a8') index_tot_mass(8,3) = get_spc_ndx('mlip_a8') ! ! marine organics - Aitken marine ! index_tot_mass(9,1) = get_spc_ndx('mpoly_a9') index_tot_mass(9,2) = get_spc_ndx('mprot_a9') index_tot_mass(9,3) = get_spc_ndx('mlip_a9') ! #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_ma', 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 if (vlc_dry(i,pver,jvlc) .ne. 0._r8) then 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) endif 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 if (vlc_dry(i,pver,jvlc) .ne. 0._r8) then 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) endif 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 if (vlc_dry(i,pver,jvlc) .ne. 0._r8) then 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) endif 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 ! REASTER 08/04/2015 BEGIN !============================================================================= !============================================================================= subroutine mam_prevap_resusp_init( ) use modal_aero_data, only: & lmassptr_amode, lspectype_amode, & modeptr_coarse, & nspec_amode, ntot_amode, numptr_amode, & mam_prevap_resusp_optaa, mmtoo_prevap_resusp, ntoo_prevap_resusp use phys_control, only: phys_getopts integer :: lspec, lspec2 integer :: mm, mmtoo, mm2 integer :: n, ntoo, nch character(len=100) :: msg if ( masterproc ) then write(iulog,'(/a)') 'mam_prevap_resusp_init' write(iulog,'(a,i10)') 'mam_prevap_resusp_optaa', mam_prevap_resusp_optaa end if ! calculate pointers for resuspension ! when mam_prevap_resusp_optaa = 30, mmtoo_prevap_resusp values are ! >0 for aerosol mass species with coarse mode counterpart ! -1 for aerosol mass species WITHOUT coarse mode counterpart ! -2 for aerosol number species ! 0 for other species ! when mam_prevap_resusp_optaa = 0, 10, 11, 20, 21, mmtoo_prevap_resusp values are 0 mmtoo_prevap_resusp(:) = 0 ntoo_prevap_resusp(:) = 0 if ( mam_prevap_resusp_optaa == 30 ) then #if ( defined MODAL_AERO_3MODE ) || ( defined MODAL_AERO_4MODE ) || ( defined MODAL_AERO_4MODE_MOM ) || ( defined MODAL_AERO_5MODE_AGEDCARBON ) ntoo = modeptr_coarse #else call endrun( 'modal_aero_wetscav_init: new resuspension not implemented for 7-mode or 9-mode MAM.') #endif do n = 1, ntot_amode ! loop over aerosol modes that was wet-removed do lspec = 1, nspec_amode(n) ! loop over chem constituents that was wet-removed mm = lmassptr_amode(lspec,n) ! q-array index of the species that was wet-removed nch = len( trim( cnst_name(mm) ) ) - 1 if (n >= 10) nch = nch - 1 if (n >= 100) nch = nch - 1 ! lspectype = lspectype_amode(lspec,n) mmtoo = -1 ! q-array index of the coarse mode species that gets the resuspension do lspec2 = 1, nspec_amode(ntoo) ! match based on the cnst_name (except for the last 1-2 characters) mm2 = lmassptr_amode(lspec2,ntoo) if ( cnst_name(mm)(1:nch) == cnst_name(mm2)(1:nch) ) then mmtoo = mm2 exit end if end do if (mmtoo <= 0) then msg = 'modal_aero_wetscav_init: cannot find mmtoo for resuspension of: ' // cnst_name(mm) write(iulog,'(a)') msg call endrun( msg ) end if if ( masterproc ) then write(iulog,'(a,3(2x,a))') 'modal_aero_wetscav_init mmfrm/too: ', & cnst_name(mm), cnst_name(mmtoo), cnst_name(numptr_amode(ntoo)) end if mmtoo_prevap_resusp(mm) = mmtoo ntoo_prevap_resusp(mm) = ntoo end do ! lspec mm = numptr_amode(n) mmtoo_prevap_resusp(mm) = -2 ntoo_prevap_resusp(mm) = ntoo end do ! n end if if ( masterproc ) then do mm = 1, pcnst mmtoo = mmtoo_prevap_resusp(mm) ntoo = ntoo_prevap_resusp(mm) msg = ' ' if (mmtoo > 0) msg = cnst_name(mmtoo) write(iulog,'(2a,3(1x,i9),2x,a)') 'name, mm, mmtoo, ntoo = ', & cnst_name(mm), mm, mmtoo, ntoo, trim(msg) end do end if end subroutine mam_prevap_resusp_init ! REASTER 08/04/2015 END !============================================================================= !============================================================================= subroutine aero_model_wetdep(dt, dlf, dlf2, cmfmc2, state, &!Intent-ins sh_e_ed_ratio, mu, md, du, eu, ed, dp, dsubcld, jt, maxg, ideep, lengath,& species_class, & cam_out, & !Intent-inout pbuf, & !Pointer ptend, & !Intent-out clear_rh ) !optional use modal_aero_deposition, only: set_srf_wetdep use wetdep, only: wetdepa_v2, wetdep_inputs_set, & wetdep_inputs_unset, 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 use mo_constants, only: pi use infnan, only: nan, assignment(=) ! 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] real(r8), intent(in) :: dlf2(:,:) ! Shal conv cldwtr detrainment (kg/kg/s - grid avg) real(r8), intent(in) :: cmfmc2(pcols,pverp) ! Shal conv mass flux (kg/m2/s) real(r8), intent(in) :: sh_e_ed_ratio(pcols,pver) ! shallow conv [ent/(ent+det)] ratio ! mu, md, ..., ideep, lengath are all deep conv variables ! *** AND ARE GATHERED *** real(r8), intent(in) :: mu(pcols,pver) ! Updraft mass flux (positive) real(r8), intent(in) :: md(pcols,pver) ! Downdraft mass flux (negative) real(r8), intent(in) :: du(pcols,pver) ! Mass detrain rate from updraft real(r8), intent(in) :: eu(pcols,pver) ! Mass entrain rate into updraft real(r8), intent(in) :: ed(pcols,pver) ! Mass entrain rate into downdraft ! eu, ed, du are "d(massflux)/dp" and are all positive real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc integer, intent(in) :: jt(pcols) ! Index of cloud top for each column integer, intent(in) :: maxg(pcols) ! Index of cloud top for each column integer, intent(in) :: ideep(pcols) ! Gathering array integer, intent(in) :: lengath ! Gathered min lon indices over which to operate integer, intent(in) :: species_class(:) type(cam_out_t), intent(inout) :: cam_out ! export state type(physics_buffer_desc), pointer :: pbuf(:) type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies real(r8), optional, intent(in) :: clear_rh(pcols,pver) ! optional clear air relative humidity ! that gets passed to modal_aero_wateruptake_dr ! local vars integer :: i integer :: jnv ! index for scavcoefnv 3rd dimension integer :: jnummaswtr ! indicates current aerosol species type (0 = number, 1 = dry mass, 2 = water) integer, parameter :: jaeronumb=0, jaeromass=1, jaerowater=2 integer :: k integer :: lchnk ! chunk identifier integer :: lphase ! index for interstitial / cloudborne aerosol integer :: lspec ! index for aerosol number / chem-mass / water-mass integer :: lspectype integer :: lcoardust, lcoarnacl ! indices for coarse mode dust and seasalt masses integer :: m, mtmp ! mode index integer :: mm, mmai, mmtoo ! tracer (q-array) index integer :: ncol ! number of atmospheric columns integer :: mam_prevap_resusp_optaa10, mam_prevap_resusp_optcc real(r8) :: iscavt(pcols, pver) 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) :: sol_factbi, sol_factii, sol_factiic real(r8) :: sflx(pcols) ! deposition flux real(r8) :: d1p_prevap_resusp, v1p_prevap_resusp real(r8) :: dqdt_tmp(pcols,pver) ! temporary array to hold tendency for the "current" aerosol species real(r8) :: dqdt_sv(pcols,pver,pcnst) ! temporary array to hold tendency for all interstitial aerosol 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) :: qqcw_tmp(pcols,pver) ! temporary array to hold qqcw ! rce 2010/05/01 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 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) :: tmp_evapdp, tmp_evapsh !RCE real(r8) :: tmp_precdp, tmp_precsh !RCE real(r8) :: tmp_resudp, tmp_resush !RCE real(r8) :: sflxec(pcols), sflxecdp(pcols) ! deposition flux !RCE real(r8) :: sflxic(pcols), sflxicdp(pcols) ! deposition flux !RCE real(r8) :: sflxbc(pcols), sflxbcdp(pcols) ! deposition flux !RCE real(r8) :: rcscavt(pcols, pver) !RCE real(r8) :: rsscavt(pcols, pver) !RCE real(r8) :: qqcw_in(pcols,pver), qqcw_sav(pcols,pver,0:maxd_aspectype) ! temporary array to hold qqcw for the current mode !RCE ! real(r8) :: rtscavt_sv(pcols, pver, 0:maxd_aspectype) !RCE real(r8) :: rtscavt_sv(pcols, pver, pcnst) ! REASTER 08/12/2015 real(r8) :: rcscavt_cn_sv(pcols, pver) ! REASTER 08/12/2015 real(r8) :: rsscavt_cn_sv(pcols, pver) ! REASTER 08/12/2015 real(r8), pointer :: fldcw(:,:) real(r8), pointer :: dgnumwet(:,:,:) real(r8), pointer :: qaerwat(:,:,:) ! aerosol water real(r8), pointer :: rate1ord_cw2pr_st(:,:) real(r8), pointer :: fracis(:,:,:) ! fraction of transported species that are insoluble integer, parameter:: nsrflx_mzaer2cnvpr = 2 !RCE 2012/01/12 bgn real(r8) :: aerdepwetis(pcols,pcnst) ! aerosol wet deposition (interstitial) real(r8) :: aerdepwetcw(pcols,pcnst) ! aerosol wet deposition (cloud water) real(r8) :: qsrflx_mzaer2cnvpr(pcols,pcnst,nsrflx_mzaer2cnvpr) real(r8) :: rprddpsum(pcols), rprdshsum(pcols) ! RCE 2012/01/12 real(r8) :: evapcdpsum(pcols), evapcshsum(pcols) ! RCE 2012/01/12 real(r8), pointer :: rprddp(:,:) ! rain production, deep convection real(r8), pointer :: rprdsh(:,:) ! rain production, deep convection real(r8), pointer :: evapcsh(:,:) ! Evaporation rate of shallow convective precipitation >=0. real(r8), pointer :: evapcdp(:,:) ! Evaporation rate of deep convective precipitation >=0. real(r8), pointer :: icwmrdp(:,:) ! in cloud water mixing ratio, deep convection real(r8), pointer :: icwmrsh(:,:) ! in cloud water mixing ratio, deep convection real(r8), pointer :: sh_frac(:,:) ! Shallow convective cloud fraction real(r8), pointer :: dp_frac(:,:) ! Deep convective cloud fraction character(len=100) :: msg type(wetdep_inputs_t) :: dep_inputs lchnk = state%lchnk ncol = state%ncol call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep_ma', 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 ! If not using prognostic aerosol call the diagnostic version ! Calculate aerosol size distribution parameters ! for prognostic modal aerosols the transfer of mass between aitken and ! accumulation modes is done in conjunction with the dry radius calculation call t_startf('calcsize') call modal_aero_calcsize_sub(state, dt, pbuf, ptend) call t_stopf('calcsize') ! Aerosol water uptake call t_startf('wateruptake') if (present(clear_rh)) then ! clear_rh allows us to provide alternate calculation of clear air RH call modal_aero_wateruptake_dr(state, pbuf, clear_rh_in=clear_rh) else call modal_aero_wateruptake_dr(state, pbuf) endif 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, rate1_cw2pr_st_idx, rate1ord_cw2pr_st) call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) ) !Compute variables needed for convproc unified convective transport call pbuf_get_field(pbuf, rprddp_idx, rprddp ) call pbuf_get_field(pbuf, rprdsh_idx, rprdsh ) call pbuf_get_field(pbuf, nevapr_shcu_idx, evapcsh ) call pbuf_get_field(pbuf, nevapr_dpcu_idx, evapcdp ) evapcdpsum(:) = 0.0_r8 rprddpsum(:) = 0.0_r8 !RCE 2012/01/12 bgn evapcshsum(:) = 0.0_r8 rprdshsum(:) = 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 enddo !RCE 2012/01/12 end 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 !RCE aerdepwetis(:,:) = 0.0_r8 !RCE aerdepwetcw(:,:) = 0.0_r8 !RCE qqcw_tmp(:,:) = 0.0_r8 !RCE else !Fill in nans for the default case 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 if ( mam_prevap_resusp_optaa >= 20 ) then ! resuspension goes to a different phase or mode rtscavt_sv(:,:,:) = 0.0_r8 rcscavt_cn_sv(:,:) = 0.0_r8 rsscavt_cn_sv(:,:) = 0.0_r8 endif mmode_loop_aa: & ! REASTER 08/11/2015 BEGIN do mtmp = 1, ntot_amode ! main loop over aerosol modes m = mtmp if (ntot_amode == 4) then ! for mam4, do accum, aitken, pcarbon, then coarse if (mtmp == modeptr_coarse) then m = ntot_amode else if (mtmp > modeptr_coarse) then m = mtmp - 1 endif endif ! REASTER 08/11/2015 END !BSINGH: loop counters (strt_loop,end_loop and stride_loop) are selected based on whether ! convproc_do_aer is true or false in the aero_model_init subroutine above lphase_loop_aa: & 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) ) if (sscav_tuning) then sol_factb = 0.03_r8 ! all below-cloud scav ON (0.1 "tuning factor") ! tuned 1/6 else sol_factb = 0.1_r8 ! all below-cloud scav ON (0.1 "tuning factor") endif ! sol_factb = 0.03_r8 ! all below-cloud scav ON (0.1 "tuning factor") ! tuned 1/6 sol_facti = 0.0_r8 ! strat in-cloud scav totally OFF for institial sol_factic = 0.4_r8 ! xl 2010/05/20 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 if (sscav_tuning) then f_act_conv = 0.4_r8 ! rce 2010/05/02 else f_act_conv = 0.8_r8 ! rce 2010/05/02 endif 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") if (sscav_tuning) then sol_facti = min(0.6_r8, sol_facti_cloud_borne) ! strat in-cloud scav totally ON for cloud-borne ! tuned 1/6 else sol_facti = sol_facti_cloud_borne ! strat in-cloud scav cloud-borne tuning factor endif 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 ! RCE 2012/01/12 ! 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 6 "sol_fact" parameters: ! sol_facti, sol_factic, sol_factb for liquid cloud ! sol_factii, sol_factiic, sol_factbi for ice cloud ! the ice cloud parameters are optional, and if not provided, they default to ! one of the other sol_fact parameters (see subr. wetdepa about this) ! for now, we set the ice cloud parameters equal ! to their liquid cloud counterparts ! currently the ice parameters are not used in wetdepa as ! wetdepa sets "weight" (the ice cloud fraction) to 0.0 ! if this changes, we will have to give more thought to ! the ice cloud parameter values ! sol_factbi = sol_factb sol_factii = sol_facti sol_factiic = sol_factic(1,1) ! REASTER 08/12/2015 - changed ordering (mass then number) for prevap resuspend to coarse lspec_loop_aa: & do lspec = 1, nspec_amode(m)+2 ! loop over number + chem constituents + water mmai = 0 if (lspec <= nspec_amode(m)) then ! non-water mass jnummaswtr = jaeromass if (lphase == 1) then mm = lmassptr_amode(lspec,m) jnv = 2 else mm = lmassptrcw_amode(lspec,m) mmai = lmassptr_amode(lspec,m) jnv = 0 endif else if (lspec == nspec_amode(m)+1) then ! number jnummaswtr = jaeronumb if (lphase == 1) then mm = numptr_amode(m) jnv = 1 else mm = numptrcw_amode(m) mmai = numptr_amode(m) jnv = 0 endif else ! water mass ! bypass wet removal of aerosol water jnummaswtr = jaerowater 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 mam_prevap_resusp_optaa10 = mam_prevap_resusp_optaa/10 ! mam_prevap_resusp_optcc values control the prevap_resusp calculations in wetdepa_v2: ! 0 = no resuspension ! 1 = linear resuspension of aerosol mass or number following original mam coding ! and history_aero_prevap_resusp = .false. ! 2 = same as 1 but history_aero_prevap_resusp = .true. ! 3 = same as 2 but with some added "xxx = max( 0, xxx)" lines ! ! 130 = non-linear resuspension of aerosol mass based on scavenged aerosol mass ! 230 = non-linear resuspension of aerosol number based on raindrop number ! the 130 thru 230 all use the new prevap_resusp code block in subr wetdepa_v2 ! and all have history_aero_prevap_resusp = .true. ! mam_prevap_resusp_optcc = 0 if ( mam_prevap_resusp_optaa == 0 ) then mam_prevap_resusp_optcc = 0 else if ( mam_prevap_resusp_optaa == 10 ) then mam_prevap_resusp_optcc = 1 else if ( mam_prevap_resusp_optaa == 11) then mam_prevap_resusp_optcc = 2 else if ( mam_prevap_resusp_optaa == 20 ) then mam_prevap_resusp_optcc = 2 else if ( mam_prevap_resusp_optaa == 21 ) then mam_prevap_resusp_optcc = 3 else if ( mam_prevap_resusp_optaa == 30 ) then if ( jnummaswtr == jaeromass ) then mam_prevap_resusp_optcc = 130 else if ( jnummaswtr == jaeronumb .and. & lphase == 1 .and. & m == modeptr_coarse ) then mam_prevap_resusp_optcc = 230 endif endif if ( mam_prevap_resusp_optcc /= 0 .and. & mam_prevap_resusp_optcc /= 1 .and. & mam_prevap_resusp_optcc /= 2 .and. & mam_prevap_resusp_optcc /= 3 .and. & mam_prevap_resusp_optcc /= 130 .and. & mam_prevap_resusp_optcc /= 230 ) then write(msg,'(a,2(1x,i10))') & 'aero_model_wetdep - bad mam_prevap_resusp_optaa & cc =', & mam_prevap_resusp_optaa, mam_prevap_resusp_optcc call endrun( msg ) endif ! 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 (jnummaswtr == jaeromass) 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 lphase_jnmw_conditional: & if ((lphase == 1) .and. (jnummaswtr /= jaerowater)) then ! 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 !PMA ! if (m .ne. 3) scavcoefnv(:,:,jnv) = scavcoefnv(:,:,jnv) * 10._r8 call wetdepa_v2( & ncol, dt, & state%t, state%pmid, state%q(:,:,1), state%pdel, & dep_inputs%cmfdqr, dep_inputs%evapc, dlf, dep_inputs%conicw, & dep_inputs%prain, dep_inputs%qme, dep_inputs%evapr, dep_inputs%totcond, & dep_inputs%cldt, dep_inputs%cldcu, & dep_inputs%cldv, dep_inputs%cldvcu, dep_inputs%cldvst, & sol_factb, sol_factbi, sol_facti, sol_factii, sol_factic, sol_factiic, & mam_prevap_resusp_optcc, .false., scavcoefnv(:,:,jnv), rate1ord_cw2pr_st, f_act_conv, & q_tmp, qqcw_in(:,:), & fracis(:,:,mm), dqdt_tmp, iscavt, & icscavt, isscavt, bcscavt, bsscavt, rcscavt, rsscavt ) ! REASTER 08/12/2015 BEGIN if ( mam_prevap_resusp_optaa10 == 3 ) then ! resuspension goes to coarse mode ! first deduct the current resuspension from the dqdt_tmp of the current species dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - ( rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ) ! then add the current resuspension to the rtscavt_sv of the appropriate coarse mode species mmtoo = mmtoo_prevap_resusp(mm) if (mmtoo > 0) rtscavt_sv(1:ncol,:,mmtoo) = rtscavt_sv(1:ncol,:,mmtoo) & + ( rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ) ! then add the rtscavt_sv of the current species to the dqdt_tmp of the current species ! note that for so4_a3 and mam3, the rtscavt_sv at this point will have resuspension contributions ! from so4_a1/2/3 and so4c1/2/3 dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt_sv(1:ncol,:,mm) endif ! REASTER 08/12/2015 END ! if (convproc_do_aer) then if ( mam_prevap_resusp_optaa10 == 2 ) then ! add resuspension of cloudborne species to dqdt of interstitial species ! dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt_sv(1:ncol,:,lspec) ! RCE 2012/01/12 dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) + rtscavt_sv(1:ncol,:,mm) ! REASTER 08/12/2015 endif do_hygro_sum_del = .false. ! if ( lspec > 0 ) do_hygro_sum_del = .true. if ( jnummaswtr == jaeromass ) do_hygro_sum_del = .true. if (convproc_do_aer) then do_hygro_sum_del = .false. ! if ( (lspec > 0) .and. do_aero_water_removal ) then !RCE 2012/01/12 if ( (jnummaswtr == jaeromass) .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) if (mm == 24) then call outfld(scavcoef_name(1),scavcoefnv(1:ncol,1:pver,jnv), ncol, lchnk ) endif if (mm == 29) then call outfld(scavcoef_name(2),scavcoefnv(1:ncol,1:pver,jnv), ncol, lchnk ) endif if (mm == 37) then call outfld(scavcoef_name(3),scavcoefnv(1:ncol,1:pver,jnv), ncol, lchnk ) endif if (mm == 41) then call outfld(scavcoef_name(4),scavcoefnv(1:ncol,1:pver,jnv), ncol, lchnk ) endif 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 ( history_aero_prevap_resusp ) then ! here the prevap resuspension is in rcscavt & rsscavt and column integral is written to history !BSINGH(09/15/2014):Following two nested do-loops are new additions for unified convection !BSINGH(09/15/2014):After these do-loops, code was added by RCE, the comments by RCE are kept as it is 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 if (convproc_do_aer) 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) else if (convproc_do_aer) sflxec(:) = 0._r8 end if ! ( history_aero_prevap_resusp ) !RCE 2012/01/12 bgn - next ~40 lines are new ! 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 ( convproc_do_aer ) then if ( deepconv_wetdep_history ) then do i = 1, ncol tmp_precdp = max( rprddpsum(i), 1.0e-35_r8 ) tmp_precsh = max( rprdshsum(i), 1.0e-35_r8 ) tmp_evapdp = max( evapcdpsum(i), 0.1e-35_r8 ) tmp_evapsh = max( evapcshsum(i), 0.1e-35_r8 ) ! assume that in- and below-cloud removal are proportional to column precip production tmpa = tmp_precdp / (tmp_precdp + tmp_precsh) 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)] if ( history_aero_prevap_resusp ) then tmp_resudp = tmpa * min( (tmp_evapdp/tmp_precdp), 1.0_r8 ) tmp_resush = (1.0_r8 - tmpa) * min( (tmp_evapsh/tmp_precsh), 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 else sflxecdp(i) = 0.0_r8 end if end do call outfld( trim(cnst_name(mm))//'SFSBD', sflxbcdp, pcols, lchnk) else 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) end if ! ( convproc_do_aer ) !RCE 2012/01/12 end - prev ~40 lines are new if (do_hygro_sum_del) then tmpa = spechygro(lspectype_amode(lspec,m))/ & specdens_amode(lspectype_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 else if ((lphase == 1) .and. (jnummaswtr == jaerowater)) then lphase_jnmw_conditional 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 ! (do_lphase1) elseif (lphase == 2) then lphase_jnmw_conditional ! lphase == 2 do_lphase2 = .true. ! REASTER 08/11/2015 start ! There is no cloud-borne aerosol water in the model, so the do_lphase2 code block ! should NEVER execute for lspec = nspec_amode(m)+1 (i.e., jnummaswtr = jaerowater). ! The code only worked because the "do lspec" loop cycles when lspec = nspec_amode(m)+1, ! but that does not make the code correct. ! if (convproc_do_aer) then ! do_lphase2 = .false. ! if (lspec <= nspec_amode(m)) do_lphase2 = .true. ! endif if (jnummaswtr == jaerowater) do_lphase2 = .false. ! REASTER 08/11/2015 end do_lphase2_conditional: & 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,:) !RCE 2012/01/12 else qqcw_tmp(:,:) = 0.0_r8 ! rce 2010/05/01 fldcw => qqcw_get_field(pbuf, mm,lchnk) endif call wetdepa_v2( & ncol, dt, & state%t, state%pmid, state%q(:,:,1), state%pdel, & dep_inputs%cmfdqr, dep_inputs%evapc, dlf, dep_inputs%conicw, & dep_inputs%prain, dep_inputs%qme, dep_inputs%evapr, dep_inputs%totcond, & dep_inputs%cldt, dep_inputs%cldcu, & dep_inputs%cldv, dep_inputs%cldvcu, dep_inputs%cldvst, & sol_factb, sol_factbi, sol_facti, sol_factii, sol_factic, sol_factiic, & mam_prevap_resusp_optcc, .true., scavcoefnv(:,:,jnv), rate1ord_cw2pr_st, f_act_conv, & fldcw, qqcw_tmp, & fracis_cw, dqdt_tmp, iscavt, & icscavt, isscavt, bcscavt, bsscavt, rcscavt, rsscavt ) ! REASTER 08/12/2015 BEGIN if ( mam_prevap_resusp_optaa10 == 3 ) then ! resuspension goes to coarse mode ! first deduct the current resuspension from the dqdt_tmp of the current species dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - ( rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ) ! then add the current resuspension to the rtscavt_sv of the appropriate coarse mode species mmtoo = mmtoo_prevap_resusp(mm) if (mmtoo > 0) rtscavt_sv(1:ncol,:,mmtoo) = rtscavt_sv(1:ncol,:,mmtoo) & + ( rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ) endif ! REASTER 08/12/2015 END ! if (convproc_do_aer) then if ( mam_prevap_resusp_optaa10 == 2 ) then ! the original wetdepa_v2 adds the resuspension of cloudborne aerosol species ! to the dqdt of cloudborne species (as a source) ! when resus_fix=.true. and/or mam_prevap_resusp_optaa=20,21 ) then ! > save resuspension of cloudborne species to rtscavt_sv ! > deduct it from the tendency (dqdt_tmp) of the cloudborne species ! > add it to the tendency (dqdt_tmp) of the interstitial species (which is done above) ! rtscavt_sv(1:ncol,:,lspec) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ! RCE 2012/01/12 if ( 0 < mmai .and. mmai <= pcnst ) then rtscavt_sv(1:ncol,:,mmai) = rcscavt(1:ncol,:) + rsscavt(1:ncol,:) ! REASTER 08/11/2015 else write(msg,'(a,3(1x,i5))') 'aero_model_wetdep - bad mmai - m, mm, mmai =', m, mm, mmai call endrun( msg ) endif ! dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt_sv(1:ncol,:,lspec) ! RCE 2012/01/12 dqdt_tmp(1:ncol,:) = dqdt_tmp(1:ncol,:) - rtscavt_sv(1:ncol,:,mmai) ! REASTER 08/11/2015 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 ! !BSINGH(09/15/2014):Following two nested do-loops are new additions for unified convection ! REASTER 08/10/2015 - NO. They are new additions for resus_fix=.true. if ( history_aero_prevap_resusp ) 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) !RCE 2012/01/12 end - prev ~40 lines are changed endif endif do_lphase2_conditional endif lphase_jnmw_conditional enddo lspec_loop_aa ! lspec = 1, nspec_amode(m)+2 enddo lphase_loop_aa ! lphase = 1, 2 enddo mmode_loop_aa ! 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_wetdep(aerdepwetis, aerdepwetcw, cam_out) endif if (convproc_do_aer) then call pbuf_get_field(pbuf, icwmrdp_idx, icwmrdp ) call pbuf_get_field(pbuf, icwmrsh_idx, icwmrsh ) call pbuf_get_field(pbuf, sh_frac_idx, sh_frac ) call pbuf_get_field(pbuf, dp_frac_idx, dp_frac ) call t_startf('ma_convproc') call ma_convproc_intr( state, ptend, pbuf, dt, & dp_frac, icwmrdp, rprddp, evapcdp, & sh_frac, icwmrsh, rprdsh, evapcsh, & dlf, dlf2, cmfmc2, sh_e_ed_ratio, & nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr, aerdepwetis, & mu, md, du, eu, ed, dp, dsubcld, jt, maxg, ideep, lengath, & species_class, mam_prevap_resusp_optaa, & history_aero_prevap_resusp ) call t_stopf('ma_convproc') endif call wetdep_inputs_unset(dep_inputs) end subroutine aero_model_wetdep !============================================================================= !============================================================================= !------------------------------------------------------------------------- ! provides 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, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total ) use ref_pres, only : top_lev=>trop_cloud_top_lev use modal_aero_data, only : sigmag_amode use mo_constants, only : pi use modal_aero_data, only : nspec_amode use modal_aero_data, only : alnsg_amode ! 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(:) 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), target :: sad_mode(pcols,pver,ntot_amode) real(r8) :: rho_air real(r8), pointer, dimension(:,:,:) :: dgnumwet integer :: l integer :: i,k,m real(r8) :: chm_mass,tot_mass call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) ! ! compute surface aero for each mode; however, at this point we only use Aitken mode (mode 2 in MAM3; how ! can we move from hard-wiring this?) as the surface area for chemical reactions. ! sad_mode = 0._r8 sad_total = 0._r8 dm_aer = 0._r8 do i = 1,ncol do k = ltrop(i),pver 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*dgnumwet(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_total(i,k) = sum(sad_mode(i,k,:)) dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm enddo enddo sfc(:,:,:) = sad_mode(:,:,:) end subroutine aero_model_surfarea !============================================================================= !============================================================================= subroutine aero_model_gasaerexch( loffset, ncol, lchnk, delt, & latndx, lonndx, 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_amicphys, only : modal_aero_amicphys_intr 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 mo_setsox, only : setsox, has_sox use modal_aero_data, only : cnst_name_cw, qqcw_get_field use phys_control, only : phys_getopts !----------------------------------------------------------------------- ! ... 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) :: latndx(pcols) ! latitude indices integer, intent(in) :: lonndx(pcols) ! longitude indices 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 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), pointer :: fldcw(:,:) logical :: use_ECPP call phys_getopts( use_ECPP_out = use_ECPP ) call pbuf_get_field(pbuf, dgnum_idx, dgnum, start=(/1,1,1/), kount=(/pcols,pver,ntot_amode/) ) 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),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. use_ECPP) then ! non-MMF Model 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 & ) endif ! Tendency due to aqueous chemistry ! When mam_amicphys_optaa > 0, dvmrdt & dvmrcwdt to hold vmr & vmrcw ! before aqueous chemistry, and cannot be used to hold aq. chem. tendencies ! ***Note - should calc & output tendencies for cloud-borne aerosol species ! rather than interstitial here if (mam_amicphys_optaa <= 0) then dvmrdt = (vmr - dvmrdt) / delt dvmrcwdt = (vmrcw - dvmrcwdt) / delt endif do m = 1, gas_pcnst wrk(:) = 0._r8 do k = 1,pver if (mam_amicphys_optaa <= 0) then ! here dvmrdt is (delta vmr from aqueous chemistry)/(delt) wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit else ! here dvmrdt is vmr before aqueous chemistry, so need to calculate (delta vmr)/(delt) wrk(:ncol) = wrk(:ncol) + ((vmr(:ncol,k,m)-dvmrdt(:ncol,k,m))/delt) & * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit endif end do name = 'AQ_'//trim(solsym(m)) call outfld( name, wrk(:ncol), ncol, lchnk ) enddo !------------------------------------------------------ else if (use_ECPP) then ! MMF w/ ECPP ! when ECPP is used, aqueous chemistry is done in ECPP, and not updated here. ! Minghuai Wang, 2010-02 (Minghuai.Wang@pnl.gov) if (mam_amicphys_optaa <= 0) then dvmrdt = 0.0_r8 dvmrcwdt = 0.0_r8 else dvmrdt(:ncol,:,:) = vmr(:ncol,:,:) dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:) end if end if ! use_ECPP !------------------------------------------------------ if (mam_amicphys_optaa <= 0) then ! do gas-aerosol exchange, nucleation, and coagulation using old routines ! 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') call modal_aero_gasaerexch_sub( & lchnk, ncol, nstep, & loffset, delt, & tfld, pmid, pdel, & vmr, vmrcw, & dvmrdt, dvmrcwdt, & dgnum, dgnumwet ) 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') ! do aerosol nucleation (new particle formation) call t_startf('modal_nucl') 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') ! do aerosol coagulation call t_startf('modal_coag') call modal_aero_coag_sub( & lchnk, ncol, nstep, & loffset, delt, & tfld, pmid, pdel, & vmr, & dgnum, dgnumwet, & wetdens ) call t_stopf('modal_coag') else ! (mam_amicphys_optaa > 0) ! do gas-aerosol exchange, nucleation, and coagulation using new routines call t_startf('modal_aero_amicphys') ! note that: ! vmr0 holds vmr before gas-phase chemistry ! dvmrdt and dvmrcwdt hold vmr and vmrcw before aqueous chemistry call modal_aero_amicphys_intr( & 1, 1, & 1, 1, & lchnk, ncol, nstep, & loffset, delt, & latndx, lonndx, & tfld, pmid, pdel, & zm, pblh, & qh2o, cldfr, & vmr, vmrcw, & vmr0, & dvmrdt, dvmrcwdt, & dgnum, dgnumwet, & #if ( defined NPF ) wetdens, invariants ) #else wetdens ) #endif ! subroutine modal_aero_amicphys_intr( & ! mdo_gasaerexch, mdo_rename, & ! mdo_newnuc, mdo_coag, & ! lchnk, ncol, nstep, & ! loffset, deltat, & ! latndx, lonndx, & ! t, pmid, pdel, & ! zm, pblh, & ! qv, cld, & ! q, qqcw, & ! q_pregaschem, & ! q_precldchem, qqcw_precldchem, & ! dgncur_a, dgncur_awet, & ! wetdens_host, & ! qaerwat ) call t_stopf('modal_aero_amicphys') endif ! (mam_amicphys_optaa <= 0 OR > 0) 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, & has_mam_mom, F_eff_out, nslt_om 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) :: u10(pcols) ! Needed in Gantt et al. calculation of organic mass fraction real(r8) :: F_eff(pcols) ! optional diagnostic output -- organic enrichment ratio 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 u10(: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)=u10(: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 F_eff(:)=0._r8 call seasalt_emis(u10, u10cubed, lchnk, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx, seasalt_emis_scale, F_eff) ! Write out salt mass fluxes to history files do m=1,seasalt_nbin-nslt_om 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 ! accumulated flux call outfld('SSTSFMBL',sflx(:),pcols,lchnk) ! Write out marine organic mass fluxes to history files if ( has_mam_mom ) then sflx(:)=0._r8 do m=seasalt_nbin-nslt_om+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) end do ! accumulated flux call outfld('SSTSFMBL_OM',sflx(:),pcols,lchnk) if ( F_eff_out ) then call outfld('F_eff',F_eff(:),pcols,lchnk) endif end if endif end subroutine aero_model_emissions !=============================================================================== ! private methods !=============================================================================== !=============================================================================== 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, ll, 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) lunerr = 6 dlndg_nimptblgrow = log( 1.25_r8 ) modeloop: do mode = 1, ntot_amode sigmag = sigmag_amode(mode) ll = lspectype_amode(1,mode) rhodryaero = specdens_amode(ll) 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 use ieee_arithmetic, only: ieee_is_nan use phys_control, only: phys_getopts ! !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 logical :: use_MMF ! 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 call phys_getopts( use_MMF_out = use_MMF ) !------------------------------------------------------------------------ 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 ! in the MMF NaN's were occurring here but the root cause was not ! identified, so this check was added to work around the issue if (use_MMF) then if ( ieee_is_nan(vlc_grv(i,k)) ) vlc_grv(i,k) = 0.0_r8 end if 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 end module aero_model