diff --git a/bld/build-namelist b/bld/build-namelist index faad098c53..469b883b6a 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -3898,6 +3898,31 @@ if ($chem =~ /_mam(\d)/) { } } +# SAD process participation spec_types +if ($chem =~ /_mam/) { + add_default($nl, 'sad_chem_spec_types', + 'val'=>"'sulfate','s-organic','p-organic','black-c','ammonium'"); + # sad_seasalt_spec_types is MAM-only: CARMA does not compute separate + # sea-salt SAD (its aero_model_surfarea sets sad_ssa to a sentinel value). + add_default($nl, 'sad_seasalt_spec_types', + 'val'=>"'seasalt'"); + # b4b default matches sad_chem_spec_types; mo_sad.F90 currently assumes + # all strat SAD is sulfate aerosol (see sad_strat_calc/sulfate_sad_calc). + add_default($nl, 'sad_strat_spec_types', + 'val'=>"'sulfate','s-organic','p-organic','black-c','ammonium'"); +} +if ($carma ne 'none' and $chem !~ /_mam/) { + add_default($nl, 'sad_chem_spec_types', + 'val'=>"'sulfate','s-organic','p-organic','black-c','ammonium'"); + # b4b default matches sad_chem_spec_types; mo_sad.F90 currently assumes + # all strat SAD is sulfate aerosol (see sad_strat_calc/sulfate_sad_calc). + add_default($nl, 'sad_strat_spec_types', + 'val'=>"'sulfate','s-organic','p-organic','black-c','ammonium'"); +} +if ($aer_model eq 'bam' and $prescribed_aero_model eq 'none') { + add_default($nl, 'sad_chem_spec_types', + 'val'=>"'sulfate','black-c','p-organic','s-organic','nitrate'"); +} # Gravity wave drag settings # By default, orographic waves are always on diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index ca66966f56..e5ea7cb8a7 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -6708,6 +6708,37 @@ Tuning for seasalt_emis Default: set by build-namelist. + +Species types contributing to chemistry-relevant surface area density (SAD) +for heterogeneous reaction rates. The SAD determines how much aerosol surface +is available for gas-phase chemical reactions (computed in aero_model, used +by mo_usrrxt). +Default: 'sulfate','s-organic','p-organic','black-c','ammonium' + + + +Species types contributing to sea-salt surface area density (SAD). +Computed separately from chemistry SAD and used for sea-salt-specific +heterogeneous reactions. +Default: 'seasalt' + + + +Species types contributing to stratospheric surface area density (SAD). +Used by aero_model_strat_surfarea to compute strato_sad, which is passed +to sad_strat_calc (mo_sad.F90) for STS/NAT/ICE thermodynamic partitioning. +Note: sad_strat_calc currently assumes all incoming SAD is sulfate aerosol +(via sulfate_sad_calc and calc_radius_lbs). Supporting non-sulfate species +in the stratospheric SAD pathway requires refactoring mo_sad.F90 to handle +multiple aerosol types independently. The physically correct default would +be 'sulfate' only, but the current default matches sad_chem_spec_types for +bit-for-bit consistency with prior hardcoded behavior. +Default: 'sulfate','s-organic','p-organic','black-c','ammonium' + + In-cloud solubility factor used in BULK aerosol wet removal diff --git a/src/chemistry/aerosol/aerosol_spec_utils.F90 b/src/chemistry/aerosol/aerosol_spec_utils.F90 new file mode 100644 index 0000000000..e12cb7b741 --- /dev/null +++ b/src/chemistry/aerosol/aerosol_spec_utils.F90 @@ -0,0 +1,25 @@ +! Utility module for aerosol species types +module aerosol_spec_utils + + implicit none + private + public :: spec_type_in_list + +contains + + ! returns TRUE if species type is found in type_list + logical function spec_type_in_list(spec_type, type_list) + character(len=*), intent(in) :: spec_type + character(len=*), intent(in) :: type_list(:) + integer :: i + spec_type_in_list = .false. + do i = 1, size(type_list) + if (len_trim(type_list(i)) == 0) cycle + if (trim(spec_type) == trim(type_list(i))) then + spec_type_in_list = .true. + return + end if + end do + end function spec_type_in_list + +end module aerosol_spec_utils diff --git a/src/chemistry/aerosol/aerosol_state_mod.F90 b/src/chemistry/aerosol/aerosol_state_mod.F90 index a39e3653d5..e3adcd04f9 100644 --- a/src/chemistry/aerosol/aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/aerosol_state_mod.F90 @@ -65,6 +65,7 @@ module aerosol_state_mod procedure :: convcld_actfrac procedure :: sol_factb_interstitial procedure(aero_aqu_gain_binfraction), deferred :: aqu_gain_binfraction + procedure(aero_surf_area_dens), deferred :: surf_area_dens end type aerosol_state @@ -299,6 +300,31 @@ subroutine aero_aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, end subroutine aero_aqu_gain_binfraction + !------------------------------------------------------------------------ + ! aerosol surface area density + !------------------------------------------------------------------------ + subroutine aero_surf_area_dens(self, aero_props, types_list, ncol, nlev, beglev, endlev, & + relhum, pmid, temp, sad, reff, sfc, dm_aer) + import :: aerosol_state, aerosol_properties, r8 + + class(aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object + character(len=*), intent(in) :: types_list(:) ! list of aerosol types to include + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + integer, intent(in) :: beglev(:) ! beginning model level index + integer, intent(in) :: endlev(:) ! ending model level index + real(r8), intent(in) :: relhum(:,:) ! relative humidity + real(r8), intent(in) :: pmid(:,:) ! mid-level pressure (Pa) + real(r8), intent(in) :: temp(:,:) ! temperature (K) + + real(r8), intent(out) :: sad(:,:) ! surface area density (cm2/cm3) + real(r8), intent(out) :: reff(:,:) ! effective radius (units cm) + real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin (cm2/cm3) + real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin (cm) + + end subroutine aero_surf_area_dens + end interface contains @@ -967,5 +993,4 @@ function sol_factb_interstitial(self, bin_ndx, ncol, nlev, aero_props) result(so end function sol_factb_interstitial - end module aerosol_state_mod diff --git a/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 index ae2348ec34..fbc798b366 100644 --- a/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 @@ -59,6 +59,7 @@ module bulk_aerosol_state_mod procedure :: convcld_actfrac procedure :: wgtpct procedure :: aqu_gain_binfraction + procedure :: surf_area_dens procedure :: get_bulk_num_and_mass ! for bit-for-bit procedure :: nuclice_get_numdens => nuclice_get_numdens_bam @@ -512,6 +513,232 @@ subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqga end subroutine aqu_gain_binfraction + !------------------------------------------------------------------------ + ! aerosol surface area density + !------------------------------------------------------------------------ + subroutine surf_area_dens(self, aero_props, types_list, ncol, nlev, beglev, endlev, & + relhum, pmid, temp, sad, reff, sfc, dm_aer) + use mo_constants, only : pi, avo => avogadro + use aerosol_spec_utils, only : spec_type_in_list + use ppgrid, only: pcols, pver + + class(bulk_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object + character(len=*), intent(in) :: types_list(:) ! list of aerosol types to include + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + integer, intent(in) :: beglev(:) ! beginning model level index + integer, intent(in) :: endlev(:) ! ending model level index + real(r8), intent(in) :: relhum(:,:) ! relative humidity + real(r8), intent(in) :: pmid(:,:) ! mid-level pressure (Pa) + real(r8), intent(in) :: temp(:,:) ! temperature (K) + + real(r8), intent(out) :: sad(:,:) ! surface area density (cm2/cm3) + real(r8), intent(out) :: reff(:,:) ! effective radius (units cm) + real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin (cm2/cm3) + real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin (cm) + + ! local vars + integer :: i,k + integer :: irh, rh_l, rh_u + real(r8) :: rho_air + real(r8) :: factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss + real(r8) :: dm_sulf_wet + real(r8) :: dm_orgc_wet + real(r8) :: dm_bc_wet + real(r8) :: num, vol + real(r8) :: s_exp + + real(r8), allocatable :: sadbins(:,:,:) + real(r8), allocatable :: diabins(:,:,:) + + !----------------------------------------------------------------- + ! ... parameters for log-normal distribution by number + ! references: + ! Chin et al., JAS, 59, 461, 2003 + ! Liao et al., JGR, 108(D1), 4001, 2003 + ! Martin et al., JGR, 108(D3), 4097, 2003 + !----------------------------------------------------------------- + real(r8), parameter :: rm_sulf = 6.95e-6_r8 ! mean radius of sulfate particles (cm) (Chin) + real(r8), parameter :: sd_sulf = 2.03_r8 ! standard deviation of radius for sulfate (Chin) + real(r8), parameter :: rho_sulf = 1.7e3_r8 ! density of sulfate aerosols (kg/m3) (Chin) + + real(r8), parameter :: rm_orgc = 2.12e-6_r8 ! mean radius of organic carbon particles (cm) (Chin) + real(r8), parameter :: sd_orgc = 2.20_r8 ! standard deviation of radius for OC (Chin) + real(r8), parameter :: rho_orgc = 1.8e3_r8 ! density of OC aerosols (kg/m3) (Chin) + + real(r8), parameter :: rm_bc = 1.18e-6_r8 ! mean radius of soot/BC particles (cm) (Chin) + real(r8), parameter :: sd_bc = 2.00_r8 ! standard deviation of radius for BC (Chin) + real(r8), parameter :: rho_bc = 1.0e3_r8 ! density of BC aerosols (kg/m3) (Chin) + + real(r8), parameter :: mw_so4 = 98.e-3_r8 ! so4 molecular wt (kg/mole) + + !----------------------------------------------------------------- + ! ... exponent for calculating number density + !----------------------------------------------------------------- + real(r8), parameter :: n_exp = exp( -4.5_r8*log(sd_sulf)*log(sd_sulf) ) + + real(r8), parameter :: dm_sulf = 2._r8 * rm_sulf + real(r8), parameter :: dm_orgc = 2._r8 * rm_orgc + real(r8), parameter :: dm_bc = 2._r8 * rm_bc + + real(r8), parameter :: log_sd_sulf = log(sd_sulf) + real(r8), parameter :: log_sd_orgc = log(sd_orgc) + real(r8), parameter :: log_sd_bc = log(sd_bc) + + !----------------------------------------------------------------- + ! ... table for hygroscopic growth effect on radius (Chin et al) + ! (no growth effect for mineral dust) + !----------------------------------------------------------------- + real(r8), parameter :: table_rh(7) = (/ 0.0_r8, 0.5_r8, 0.7_r8, 0.8_r8, 0.9_r8, 0.95_r8, 0.99_r8 /) + real(r8), parameter :: table_rfac_sulf(7) = (/ 1.0_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 1.9_r8, 2.2_r8 /) + real(r8), parameter :: table_rfac_oc(7) = (/ 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 2.2_r8 /) + real(r8), parameter :: table_rfac_bc(7) = (/ 1.0_r8, 1.0_r8, 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.9_r8 /) + real(r8), parameter :: table_rfac_ss(7) = (/ 1.0_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.4_r8, 2.9_r8, 4.8_r8 /) + + real(r8), pointer :: mmr(:,:) ! interstitial aerosol mass, number mixing ratios + integer :: ispc, ibin, ndx, astat + character(len=32) :: spectype + + sad = 0.0_r8 + reff = 0.0_r8 + if (present(sfc)) sfc = 0.0_r8 + if (present(dm_aer)) dm_aer = 0.0_r8 + + allocate(sadbins(ncol,nlev,aero_props%nbins()),stat=astat) + if( astat/= 0 ) call endrun('bulk_aerosol_state_mod%surf_area_dens: sadbins allocate error') + allocate(diabins(ncol,nlev,aero_props%nbins()),stat=astat) + if( astat/= 0 ) call endrun('bulk_aerosol_state_mod%surf_area_dens: diabins allocate error') + + sadbins = 0._r8 + diabins = 0._r8 + + ver_loop: do k = 1,nlev + col_loop: do i = 1,ncol + !------------------------------------------------------------------------- + ! ... air density (kg/m3) + !------------------------------------------------------------------------- + rho_air = pmid(i,k)/(temp(i,k)*287.04_r8) + !------------------------------------------------------------------------- + ! ... aerosol growth interpolated from M.Chin's table + !------------------------------------------------------------------------- + if (relhum(i,k) >= table_rh(7)) then + rfac_sulf = table_rfac_sulf(7) + rfac_oc = table_rfac_oc(7) + rfac_bc = table_rfac_bc(7) + else + do irh = 2,7 + if (relhum(i,k) <= table_rh(irh)) then + exit + end if + end do + rh_l = irh-1 + rh_u = irh + + factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l)) + + rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l)) + rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l)) + rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l)) + end if + + dm_sulf_wet = dm_sulf * rfac_sulf + dm_orgc_wet = dm_orgc * rfac_oc + dm_bc_wet = dm_bc * rfac_bc + + dm_bc_wet = min(dm_bc_wet ,50.e-6_r8) ! maximum size is 0.5 micron (Chin) + dm_orgc_wet = min(dm_orgc_wet,50.e-6_r8) ! maximum size is 0.5 micron (Chin) + + ndx=0 + do ibin = 1, aero_props%nbins() + do ispc = 1, aero_props%nspecies(ibin) + call aero_props%get(bin_ndx=ibin, species_ndx=ispc, spectype=spectype) + call self%get_ambient_mmr(species_ndx=ispc, bin_ndx=ibin, mmr=mmr) + + if (.not. spec_type_in_list(spectype, types_list)) cycle + + select case ( trim(spectype) ) + case('sulfate') + ndx = ndx+1 + !------------------------------------------------------------------------- + ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air) + ! vol=volume density (m^3/m^3) + ! rho_aer=density of aerosol (kg/m^3) + ! vol=m*rho_air/rho_aer [kg/kg * (kg/m3)_air/(kg/m3)_aer] + !------------------------------------------------------------------------- + vol = mmr(i,k) * rho_air/rho_sulf + !------------------------------------------------------------------------- + ! calculate the number density of aerosol (aerosols/cm3) + ! assuming a lognormal distribution + ! num = (aerosols/cm3) + ! dm = geometric mean diameter + ! + ! because only the dry mass of the aerosols is known, we + ! use the mean dry radius + !------------------------------------------------------------------------- + num = vol * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp + !------------------------------------------------------------------------- + ! find surface area of aerosols using dm_wet, log_sd + ! (increase of sd due to RH is negligible) + ! and number density calculated above as distribution + ! parameters + ! sfc = surface area of wet aerosols (cm^2/cm^3) + !------------------------------------------------------------------------- + s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf) + sadbins(i,k,ndx) = num * pi * (dm_sulf_wet**2._r8) * s_exp + diabins(i,k,ndx) = dm_sulf_wet + case('black-c') + if ( aero_props%hydrophilic(ibin) ) then + ndx = ndx+1 + vol = mmr(i,k) * rho_air/rho_bc + num = vol * (6._r8/pi)*(1._r8/(dm_bc**3._r8))*n_exp + s_exp = exp(2._r8*log_sd_bc*log_sd_bc) + sadbins(i,k,ndx) = num * pi * (dm_bc_wet**2._r8) * s_exp + diabins(i,k,ndx) = dm_bc_wet + end if + case('p-organic') + if ( aero_props%hydrophilic(ibin) ) then + ndx = ndx+1 + vol = mmr(i,k) * rho_air/rho_orgc + num = vol * (6._r8/pi)*(1._r8/(dm_orgc**3))*n_exp + s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) + sadbins(i,k,ndx) = num * pi * (dm_orgc_wet**2._r8) * s_exp + diabins(i,k,ndx) = dm_orgc_wet + end if + case('s-organic') + ndx = ndx+1 + vol = mmr(i,k) * rho_air/rho_orgc + num = vol * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp + s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) + sadbins(i,k,ndx) = num * pi * (dm_orgc_wet**2._r8) * s_exp + diabins(i,k,ndx) = dm_orgc_wet + case('nitrate') + ndx = ndx+1 + vol = mmr(i,k) * rho_air/rho_sulf + num = vol * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp + s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf) + sadbins(i,k,ndx) = num * pi * (dm_sulf_wet**2._r8) * s_exp + diabins(i,k,ndx) = dm_sulf_wet + end select + end do + end do + + !------------------------------------------------------------------------- + ! ... add up total surface area density for output + !------------------------------------------------------------------------- + sad(i,k) = sum(sadbins(i,k,:)) + + enddo col_loop + enddo ver_loop + + if (present(sfc)) sfc(:ncol,:nlev,:ndx) = sadbins(:ncol,:nlev,:ndx) + if (present(dm_aer)) dm_aer(:ncol,:nlev,:ndx) = diabins(:ncol,:nlev,:ndx) + + deallocate(sadbins) + deallocate(diabins) + + end subroutine surf_area_dens + !------------------------------------------------------------------------------ ! Compute BAM number concentration (#/m3) and mass concentration (kg/m3) ! for a single bin. Applies bam_sulfate_scale only to SULFATE (not volcanic). diff --git a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 index 21efbb5e42..ee9086b5b8 100644 --- a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 @@ -13,7 +13,7 @@ module carma_aerosol_state_mod use physconst, only: pi use carma_intr, only: carma_get_total_mmr, carma_get_dry_radius, carma_get_number, carma_get_number_cld use carma_intr, only: carma_get_group_by_name, carma_get_kappa, carma_get_dry_radius, carma_get_wet_radius - use carma_intr, only: carma_get_wght_pct + use carma_intr, only: carma_get_wght_pct, carma_effecitive_radius, carma_get_sad use ppgrid, only: begchunk, endchunk, pcols, pver implicit none @@ -50,6 +50,7 @@ module carma_aerosol_state_mod procedure :: water_volume procedure :: wet_diameter procedure :: aqu_gain_binfraction + procedure :: surf_area_dens final :: destructor @@ -664,4 +665,144 @@ subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqga end subroutine aqu_gain_binfraction + !------------------------------------------------------------------------ + ! aerosol surface area density + !------------------------------------------------------------------------ + subroutine surf_area_dens(self, aero_props, types_list, ncol, nlev, beglev, endlev, & + relhum, pmid, temp, sad, reff, sfc, dm_aer) + use aerosol_spec_utils, only : spec_type_in_list + + class(carma_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object + character(len=*), intent(in) :: types_list(:) ! list of aerosol types to include + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + integer, intent(in) :: beglev(:) ! beginning model level index + integer, intent(in) :: endlev(:) ! ending model level index + real(r8), intent(in) :: relhum(:,:) ! relative humidity + real(r8), intent(in) :: pmid(:,:) ! mid-level pressure (Pa) + real(r8), intent(in) :: temp(:,:) ! temperature (K) + + real(r8), intent(out) :: sad(:,:) ! surface area density (cm2/cm3) + real(r8), intent(out) :: reff(:,:) ! effective radius (units cm) + real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin (cm2/cm3) + real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin (cm) + + ! local vars + real(r8) :: reffaer(pcols,pver) ! bulk effective radius in cm + + integer :: icol, ilev, ibin, ispec, ierr + real(r8) :: chm_mass, tot_mass + character(len=32) :: spectype + real(r8) :: wetr(pcols,pver) ! CARMA bin wet radius in cm + real(r8) :: wetrho(pcols,pver) ! CARMA bin wet density + real(r8) :: sad_carma(pcols,pver) ! CARMA bin wet surface area density in cm2/cm3 + real(r8), pointer :: aer_bin_mmr(:,:) + + character(len=aero_name_len) :: bin_name, shortname + integer :: igroup, indxbin, rc, nchr + + real(r8), allocatable :: sad_bins(:,:,:) + + integer :: nbins, nspec_max + type(ptr2d_t), allocatable :: mmr_ptr(:,:) ! interstitial mmr field per (bin,species) + logical, allocatable :: in_list(:,:) ! species type is in types_list + + nbins = aero_props%nbins() + + allocate(sad_bins(ncol,nlev,nbins), stat=ierr) + if (ierr/=0) then + call endrun('carma_aerosol_state::surf_area_dens: not able to allocate sad_bins') + end if + + nspec_max = 0 + do ibin = 1,nbins + nspec_max = max(nspec_max, aero_props%nspecies(ibin)) + end do + + allocate(mmr_ptr(nbins,nspec_max), in_list(nbins,nspec_max), stat=ierr) + if (ierr/=0) then + call endrun('carma_aerosol_state::surf_area_dens: not able to allocate bin lookup arrays') + end if + + sad = 0._r8 + reff = 0._r8 + sad_bins = 0._r8 + if (present(dm_aer)) dm_aer = 0._r8 + + ! + ! Compute surface aero for each bin. + ! Total over all bins as the surface area for chemical reactions. + ! + + reffaer = carma_effecitive_radius(self%host_%state) + + in_list(:,:) = .false. + + do ibin=1,nbins ! loop over aerosol bins + call rad_aer_get_info_by_bin(self%list_idx_, ibin, bin_name=bin_name) + + nchr = len_trim(bin_name)-2 + shortname = bin_name(:nchr) + + call carma_get_group_by_name(shortname, igroup, rc) + + read(bin_name(nchr+1:),*) indxbin + + call carma_get_wet_radius(self%host_%state, igroup, indxbin, wetr, wetrho, rc) ! m + wetr(:ncol,:) = wetr(:ncol,:) * 1.e2_r8 ! cm + call carma_get_sad(self%host_%state, igroup, indxbin, sad_carma, rc) + + if (present(dm_aer)) dm_aer(:ncol,:,ibin) = 2._r8 * wetr(:ncol,:) ! convert wet radius (cm) to wet diameter (cm) + sad_bins(:ncol,:,ibin) = sad_carma(:ncol,:) ! cm^2/cm^3 + + ! Resolve species types and mmr field pointers that do not vary over columns or levels here: + do ispec = 1,aero_props%nspecies(ibin) + call aero_props%get(bin_ndx=ibin, species_ndx=ispec, spectype=spectype) + call self%get_ambient_mmr(species_ndx=ispec, bin_ndx=ibin, mmr=mmr_ptr(ibin,ispec)%fld) + in_list(ibin,ispec) = spec_type_in_list(spectype, types_list) + end do + end do + + do icol = 1, ncol + do ilev = beglev(icol),endlev(icol) + do ibin = 1, nbins ! loop over aerosol bins + ! + ! compute a mass weighting of the number + ! + tot_mass = 0._r8 + chm_mass = 0._r8 + do ispec=1,aero_props%nspecies(ibin) + + aer_bin_mmr => mmr_ptr(ibin,ispec)%fld + + tot_mass = tot_mass + aer_bin_mmr(icol,ilev) + + if (in_list(ibin,ispec)) then + chm_mass = chm_mass + aer_bin_mmr(icol,ilev) + end if + + end do + if ( tot_mass > 0._r8 ) then + ! surface area density + sad_bins(icol,ilev,ibin) = chm_mass / tot_mass * sad_bins(icol,ilev,ibin) ! cm^2/cm^3 + else + sad_bins(icol,ilev,ibin) = 0._r8 + end if + end do + sad(icol,ilev) = sum(sad_bins(icol,ilev,:)) + reff(icol,ilev) = reffaer(icol,ilev) + + end do + end do + + if (present(sfc)) then + sfc(:ncol,:,:) = sad_bins(:ncol,:,:) + end if + + deallocate(sad_bins) + deallocate(mmr_ptr, in_list) + + end subroutine surf_area_dens + end module carma_aerosol_state_mod diff --git a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 index 17f3ad5e60..3342f1f7cf 100644 --- a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 @@ -48,6 +48,7 @@ module modal_aerosol_state_mod procedure :: convcld_actfrac procedure :: wgtpct procedure :: aqu_gain_binfraction + procedure :: surf_area_dens final :: destructor @@ -799,4 +800,171 @@ subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqga end subroutine aqu_gain_binfraction + !------------------------------------------------------------------------ + ! aerosol surface area density + !------------------------------------------------------------------------ + subroutine surf_area_dens(self, aero_props, types_list, ncol, nlev, beglev, endlev, & + relhum, pmid, temp, sad, reff, sfc, dm_aer) + use mo_constants, only : pi + use aerosol_spec_utils, only : spec_type_in_list + + class(modal_aerosol_state), intent(in) :: self + class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object + character(len=*), intent(in) :: types_list(:) ! list of aerosol types to include + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: nlev ! number of levels + integer, intent(in) :: beglev(:) ! beginning model level index + integer, intent(in) :: endlev(:) ! ending model level index + real(r8), intent(in) :: relhum(:,:) ! relative humidity + real(r8), intent(in) :: pmid(:,:) ! mid-level pressure (Pa) + real(r8), intent(in) :: temp(:,:) ! temperature (K) + + real(r8), intent(out) :: sad(:,:) ! surface area density (cm2/cm3) + real(r8), intent(out) :: reff(:,:) ! effective radius (units cm) + real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin (cm2/cm3) + real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin (cm) + + ! local vars + real(r8), allocatable ::sad_mode(:,:,:) + real(r8), allocatable ::vol_mode(:,:,:) + + real(r8) :: radeff(ncol,nlev) + real(r8) :: vol(ncol,nlev) + real(r8) :: rho_air + integer :: i,k,l,m + real(r8) :: chm_mass, tot_mass + integer :: idx_chm_val, ierr + real(r8), pointer :: mmr(:,:) ! mass mixing ratios (ncol,nlev) + real(r8), pointer :: num(:,:) + real(r8), pointer :: diam(:,:,:) + + character(len=32) :: modetype + character(len=32) :: spectype + + ! for per-mode / per-species lookups outside the column loop: + integer :: nbins, nspec_max + type(ptr2d_t), allocatable :: mmr_ptr(:,:) ! interstitial mmr field per (bin,species) + type(ptr2d_t), allocatable :: num_ptr(:) ! interstitial number field per bin + logical, allocatable :: in_list(:,:) ! species type is in types_list + logical, allocatable :: use_mode(:) ! mode contributes to the mass weighting + real(r8), allocatable :: alogsig(:) ! ln(geometric standard deviation) per mode + + character(len=*), parameter :: subname = 'modal_aerosol_state_mod::surf_area_dens' + + call pbuf_get_field(self%host_%pbuf, pbuf_get_index('DGNUMWET'), diam) + + ! + ! Compute surface aero for each mode. + ! Total over all modes as the surface area for chemical reactions. + ! + nbins = aero_props%nbins() + + allocate(sad_mode(ncol,nlev,nbins),stat=ierr) + if (ierr/=0) then + call endrun(subname//': sad_mode array allocate error') + end if + + allocate(vol_mode(ncol,nlev,nbins),stat=ierr) + if (ierr/=0) then + call endrun(subname//': vol_mode array allocate error') + end if + + ! resolve mode geometry, species types and field pointers that do not vary over columns or levels here: + nspec_max = 0 + do l = 1,nbins + nspec_max = max(nspec_max, aero_props%nspecies(l)) + end do + + allocate(mmr_ptr(nbins,nspec_max), num_ptr(nbins), in_list(nbins,nspec_max), & + use_mode(nbins), alogsig(nbins), stat=ierr) + if (ierr/=0) then + call endrun(subname//': mode lookup arrays allocate error') + end if + + in_list(:,:) = .false. + + do l = 1,nbins + call self%get_ambient_num(bin_ndx=l, num=num_ptr(l)%fld) + alogsig(l) = aero_props%alogsig(l) + + ! ignore primary carbon mode + call rad_aer_get_info(self%list_idx_, l, mode_type=modetype) + use_mode(l) = trim(modetype) /= 'primary_carbon' ! ignore the primary_carbon mode + + if (use_mode(l)) then + do m = 1,aero_props%nspecies(l) + call aero_props%get(bin_ndx=l, species_ndx=m, spectype=spectype) + call self%get_ambient_mmr(species_ndx=m, bin_ndx=l, mmr=mmr_ptr(l,m)%fld) + in_list(l,m) = spec_type_in_list(spectype, types_list) + end do + end if + end do + + sad = 0._r8 + sad_mode = 0._r8 + vol = 0._r8 + vol_mode = 0._r8 + reff = 0._r8 + + do i = 1,ncol + do k = beglev(i), endlev(i) + rho_air = pmid(i,k)/(temp(i,k)*287.04_r8) + do l=1,nbins + ! + ! compute a mass weighting of the number + ! + tot_mass = 0._r8 + chm_mass = 0._r8 + + if ( use_mode(l) ) then + + do m=1,aero_props%nspecies(l) + + mmr => mmr_ptr(l,m)%fld + + tot_mass = tot_mass + mmr(i,k) + + if (in_list(l,m)) then + chm_mass = chm_mass + mmr(i,k) + end if + end do + + end if + + if ( tot_mass > 0._r8 ) then + + num => num_ptr(l)%fld + + ! surface area density + sad_mode(i,k,l) = chm_mass /tot_mass & + * num(i,k)*rho_air*pi*diam(i,k,l)**2._r8 & + * exp(2._r8*alogsig(l)**2._r8) ! m^2/m^3 + + sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3 + + ! volume calculation, for use in effective radius calculation + vol_mode(i,k,l) = chm_mass/tot_mass & + * num(i,k)*rho_air*pi/6._r8*diam(i,k,l)**3._r8 & + * exp(4.5_r8*alogsig(l)**2._r8) ! m^3/m^3 = cm^3/cm^3 + else + sad_mode(i,k,l) = 0._r8 + vol_mode(i,k,l) = 0._r8 + end if + end do + sad(i,k) = sum(sad_mode(i,k,:)) + vol(i,k) = sum(vol_mode(i,k,:)) + reff(i,k) = 3._r8*vol(i,k)/sad(i,k) + + enddo + enddo + + if (present(sfc)) sfc(:ncol,:,:) = sad_mode(:ncol,:,:) + if (present(dm_aer)) dm_aer(:ncol,:,:) = diam(:ncol,:,:) * 1.e2_r8 ! convert m to cm + + deallocate(sad_mode) + deallocate(vol_mode) + deallocate(mmr_ptr, num_ptr, in_list, use_mode, alogsig) + + end subroutine surf_area_dens + end module modal_aerosol_state_mod diff --git a/src/chemistry/bulk_aero/aero_model.F90 b/src/chemistry/bulk_aero/aero_model.F90 index 834079dfc6..12a396e403 100644 --- a/src/chemistry/bulk_aero/aero_model.F90 +++ b/src/chemistry/bulk_aero/aero_model.F90 @@ -22,6 +22,7 @@ module aero_model use aerosol_properties_mod, only: aerosol_properties use aerosol_instances_mod, only: aerosol_instances_get_props, & aerosol_instances_get_state, aerosol_instances_get_num_models + use aerosol_state_mod, only: aerosol_state implicit none private @@ -63,6 +64,9 @@ module aero_model class(aerosol_properties), pointer :: aero_props =>null() integer :: iaermod_ = -1 + integer, parameter :: max_sad_spec = 16 + character(len=32) :: sad_chem_spec_types(max_sad_spec) = ' ' + contains !============================================================================= @@ -87,6 +91,8 @@ subroutine aero_model_readnl(nlfile) namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list namelist /aerosol_nl/ aer_sol_facti, aer_sol_factb, aer_scav_coef + namelist /aerosol_nl/ sad_chem_spec_types + !----------------------------------------------------------------------------- aer_sol_facti = nan aer_sol_factb = nan @@ -111,6 +117,7 @@ subroutine aero_model_readnl(nlfile) ! 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(sad_chem_spec_types,len(sad_chem_spec_types(1))*max_sad_spec,mpichar, 0, mpicom) call mpibcast(aer_sol_facti, pcnst, mpir8, 0, mpicom) call mpibcast(aer_sol_factb, pcnst, mpir8, 0, mpicom) call mpibcast(aer_scav_coef, pcnst, mpir8, 0, mpicom) @@ -714,8 +721,8 @@ subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf) ! called from mo_usrrxt !------------------------------------------------------------------------- subroutine aero_model_surfarea( & - state, mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, m, ltrop, & - dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total, reff_trop, sad_ssa ) + state, relhum, pmid, temp, ltrop, & + sfc, dm_aer, sad_total, reff_trop, sad_ssa ) use mo_constants, only : pi, avo => avogadro @@ -723,17 +730,8 @@ subroutine aero_model_surfarea( & type(physics_state), intent(in) :: state ! Physics state variables real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) - real(r8), intent(in) :: mmr(:,:,:) - real(r8), intent(in) :: radmean ! mean radii in cm - real(r8), intent(in) :: strato_sad(:,:) - integer, intent(in) :: ncol integer, intent(in) :: ltrop(:) - real(r8), intent(in) :: dlat(:) ! degrees latitude - integer, intent(in) :: het1_ndx real(r8), intent(in) :: relhum(:,:) - real(r8), intent(in) :: m(:,:) ! 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(:,:,:) @@ -742,287 +740,37 @@ subroutine aero_model_surfarea( & real(r8), intent(out) :: sad_ssa(:,:) ! local vars + integer :: i,k, lchnk, ncol + integer :: beglev(pcols) + integer :: endlev(pcols) - integer :: i,k - real(r8) :: rho_air - real(r8) :: v, n, n_exp, r_rd, r_sd - real(r8) :: dm_sulf, dm_sulf_wet, log_sd_sulf, sfc_sulf, sfc_nit - real(r8) :: dm_orgc, dm_orgc_wet, log_sd_orgc, sfc_oc, sfc_soa - real(r8) :: sfc_soai, sfc_soam, sfc_soab, sfc_soat, sfc_soax - real(r8) :: dm_bc, dm_bc_wet, log_sd_bc, sfc_bc - real(r8) :: rxt_sulf, rxt_nit, rxt_oc, rxt_soa - real(r8) :: c_n2o5, c_ho2, c_no2, c_no3 - real(r8) :: s_exp - - !----------------------------------------------------------------- - ! ... parameters for log-normal distribution by number - ! references: - ! Chin et al., JAS, 59, 461, 2003 - ! Liao et al., JGR, 108(D1), 4001, 2003 - ! Martin et al., JGR, 108(D3), 4097, 2003 - !----------------------------------------------------------------- - real(r8), parameter :: rm_sulf = 6.95e-6_r8 ! mean radius of sulfate particles (cm) (Chin) - real(r8), parameter :: sd_sulf = 2.03_r8 ! standard deviation of radius for sulfate (Chin) - real(r8), parameter :: rho_sulf = 1.7e3_r8 ! density of sulfate aerosols (kg/m3) (Chin) - - real(r8), parameter :: rm_orgc = 2.12e-6_r8 ! mean radius of organic carbon particles (cm) (Chin) - real(r8), parameter :: sd_orgc = 2.20_r8 ! standard deviation of radius for OC (Chin) - real(r8), parameter :: rho_orgc = 1.8e3_r8 ! density of OC aerosols (kg/m3) (Chin) - - real(r8), parameter :: rm_bc = 1.18e-6_r8 ! mean radius of soot/BC particles (cm) (Chin) - real(r8), parameter :: sd_bc = 2.00_r8 ! standard deviation of radius for BC (Chin) - real(r8), parameter :: rho_bc = 1.0e3_r8 ! density of BC aerosols (kg/m3) (Chin) - - real(r8), parameter :: mw_so4 = 98.e-3_r8 ! so4 molecular wt (kg/mole) - - integer :: irh, rh_l, rh_u - real(r8) :: factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss - logical :: zero_aerosols - - !----------------------------------------------------------------- - ! ... table for hygroscopic growth effect on radius (Chin et al) - ! (no growth effect for mineral dust) - !----------------------------------------------------------------- - real(r8), dimension(7) :: table_rh, table_rfac_sulf, table_rfac_bc, table_rfac_oc, table_rfac_ss + class(aerosol_state), pointer :: aero_state sad_ssa = -huge(1._r8) - data table_rh(1:7) / 0.0_r8, 0.5_r8, 0.7_r8, 0.8_r8, 0.9_r8, 0.95_r8, 0.99_r8/ - data table_rfac_sulf(1:7) / 1.0_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 1.9_r8, 2.2_r8/ - data table_rfac_oc(1:7) / 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 2.2_r8/ - data table_rfac_bc(1:7) / 1.0_r8, 1.0_r8, 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.9_r8/ - data table_rfac_ss(1:7) / 1.0_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.4_r8, 2.9_r8, 4.8_r8/ - - !----------------------------------------------------------------- - ! ... exponent for calculating number density - !----------------------------------------------------------------- - n_exp = exp( -4.5_r8*log(sd_sulf)*log(sd_sulf) ) - - dm_sulf = 2._r8 * rm_sulf - dm_orgc = 2._r8 * rm_orgc - dm_bc = 2._r8 * rm_bc - - log_sd_sulf = log(sd_sulf) - log_sd_orgc = log(sd_orgc) - log_sd_bc = log(sd_bc) - - reff_trop(:,:) = 0._r8 - - ver_loop: do k = 1,pver - col_loop: do i = 1,ncol - !------------------------------------------------------------------------- - ! ... air density (kg/m3) - !------------------------------------------------------------------------- - rho_air = pmid(i,k)/(temp(i,k)*287.04_r8) - !------------------------------------------------------------------------- - ! ... aerosol growth interpolated from M.Chin's table - !------------------------------------------------------------------------- - if (relhum(i,k) >= table_rh(7)) then - rfac_sulf = table_rfac_sulf(7) - rfac_oc = table_rfac_oc(7) - rfac_bc = table_rfac_bc(7) - else - do irh = 2,7 - if (relhum(i,k) <= table_rh(irh)) then - exit - end if - end do - rh_l = irh-1 - rh_u = irh - - factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l)) - - rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l)) - rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l)) - rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l)) - end if - - dm_sulf_wet = dm_sulf * rfac_sulf - dm_orgc_wet = dm_orgc * rfac_oc - dm_bc_wet = dm_bc * rfac_bc - - dm_bc_wet = min(dm_bc_wet ,50.e-6_r8) ! maximum size is 0.5 micron (Chin) - dm_orgc_wet = min(dm_orgc_wet,50.e-6_r8) ! maximum size is 0.5 micron (Chin) - - - !------------------------------------------------------------------------- - ! ... sulfate aerosols - !------------------------------------------------------------------------- - zero_aerosols = k < ltrop(i) - if ( abs( dlat(i) ) > 50._r8 ) then - zero_aerosols = pmid(i,k) < 30000._r8 - endif - !------------------------------------------------------------------------- - ! ... use ubvals climatology for stratospheric sulfate surface area density - !------------------------------------------------------------------------- - if( zero_aerosols ) then - sfc_sulf = strato_sad(i,k) - if ( het1_ndx > 0 ) then - sfc_sulf = 0._r8 ! reaction already taken into account in mo_strato_rates.F90 - end if - sfc_nit = 0._r8 - sfc_soa = 0._r8 - sfc_oc = 0._r8 - sfc_bc = 0._r8 - else - - if( so4_ndx > 0 ) then - !------------------------------------------------------------------------- - ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air) - ! v=volume density (m^3/m^3) - ! rho_aer=density of aerosol (kg/m^3) - ! v=m*rho_air/rho_aer [kg/kg * (kg/m3)_air/(kg/m3)_aer] - !------------------------------------------------------------------------- - v = mmr(i,k,so4_ndx) * rho_air/rho_sulf - !------------------------------------------------------------------------- - ! calculate the number density of aerosol (aerosols/cm3) - ! assuming a lognormal distribution - ! n = (aerosols/cm3) - ! dm = geometric mean diameter - ! - ! because only the dry mass of the aerosols is known, we - ! use the mean dry radius - !------------------------------------------------------------------------- - n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp - !------------------------------------------------------------------------- - ! find surface area of aerosols using dm_wet, log_sd - ! (increase of sd due to RH is negligible) - ! and number density calculated above as distribution - ! parameters - ! sfc = surface area of wet aerosols (cm^2/cm^3) - !------------------------------------------------------------------------- - s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf) - sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp - - else - !------------------------------------------------------------------------- - ! if so4 not simulated, use off-line sulfate and calculate as above - ! convert sulfate vmr to volume density of aerosol (cm^3_aerosol/cm^3_air) - !------------------------------------------------------------------------- - v = sulfate(i,k) * m(i,k) * mw_so4 / (avo * rho_sulf) *1.e6_r8 - n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf) - sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp - - end if - - !------------------------------------------------------------------------- - ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate) - !------------------------------------------------------------------------- - if( nit_ndx > 0 ) then - v = mmr(i,k,nit_ndx) * rho_air/rho_sulf - n = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_sulf*log_sd_sulf) - sfc_nit = n * pi * (dm_sulf_wet**2._r8) * s_exp - else - sfc_nit = 0._r8 - end if - - !------------------------------------------------------------------------- - ! hydrophylic organic carbon (follow same procedure as sulfate) - !------------------------------------------------------------------------- - if( oc2_ndx > 0 ) then - v = mmr(i,k,oc2_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_oc = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_oc = 0._r8 - end if - - !------------------------------------------------------------------------- - ! secondary organic carbon (follow same procedure as sulfate) - !------------------------------------------------------------------------- - if( soa_ndx > 0 ) then - v = mmr(i,k,soa_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soa = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soa = 0._r8 - end if - - !------------------------------------------------------------------------- - ! black carbon (follow same procedure as sulfate) - !------------------------------------------------------------------------- - if( cb2_ndx > 0 ) then - v = mmr(i,k,cb2_ndx) * rho_air/rho_bc - n = v * (6._r8/pi)*(1._r8/(dm_bc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_bc*log_sd_bc) - sfc_bc = n * pi * (dm_bc_wet**2._r8) * s_exp - else - sfc_bc = 0._r8 - end if - if( soai_ndx > 0 ) then - v = mmr(i,k,soai_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soai = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soai = 0._r8 - end if - if( soam_ndx > 0 ) then - v = mmr(i,k,soam_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soam = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soam = 0._r8 - end if - if( soab_ndx > 0 ) then - v = mmr(i,k,soab_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soab = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soab = 0._r8 - end if - if( soat_ndx > 0 ) then - v = mmr(i,k,soat_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soat = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soat = 0._r8 - end if - if( soax_ndx > 0 ) then - v = mmr(i,k,soax_ndx) * rho_air/rho_orgc - n = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp - s_exp = exp(2._r8*log_sd_orgc*log_sd_orgc) - sfc_soax = n * pi * (dm_orgc_wet**2._r8) * s_exp - else - sfc_soax = 0._r8 - end if - sfc_soa = sfc_soa + sfc_soai + sfc_soam + sfc_soab + sfc_soat + sfc_soax - - end if + lchnk = state%lchnk + ncol = state%ncol - sfc(i,k,:) = (/ sfc_sulf, sfc_nit, sfc_oc, sfc_soa, sfc_bc /) - dm_aer(i,k,:) = (/ dm_sulf_wet,dm_sulf_wet,dm_orgc_wet,dm_orgc_wet,dm_bc_wet /) + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) - !------------------------------------------------------------------------- - ! ... add up total surface area density for output - !------------------------------------------------------------------------- - sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc + beglev(:ncol)=ltrop(:ncol)+1 + endlev(:ncol)=pver - enddo col_loop - enddo ver_loop + call aero_state%surf_area_dens(aero_props, sad_chem_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, sad_total, reff_trop, sfc, dm_aer) end subroutine aero_model_surfarea !------------------------------------------------------------------------- ! stub !------------------------------------------------------------------------- - subroutine aero_model_strat_surfarea( state, ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat ) + subroutine aero_model_strat_surfarea( state, pmid, temp, ltrop, strato_sad, reff_strat ) ! dummy args type(physics_state), intent(in) :: state ! Physics state variables - integer, intent(in) :: ncol - real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) integer, intent(in) :: ltrop(:) ! tropopause level indices - type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(out) :: strato_sad(:,:) real(r8), intent(out) :: reff_strat(:,:) diff --git a/src/chemistry/carma_aero/aero_model.F90 b/src/chemistry/carma_aero/aero_model.F90 index 6dd994167d..3ebad3d1e7 100644 --- a/src/chemistry/carma_aero/aero_model.F90 +++ b/src/chemistry/carma_aero/aero_model.F90 @@ -31,6 +31,9 @@ module aero_model use carma_intr, only: carma_get_sad use aerosol_properties_mod, only: aero_name_len + use aerosol_state_mod, only: aerosol_state + use aerosol_instances_mod, only: aerosol_instances_get_props, & + aerosol_instances_get_state, aerosol_instances_get_num_models implicit none private @@ -43,7 +46,7 @@ module aero_model public :: aero_model_wetdep ! aerosol wet removal public :: aero_model_emissions ! aerosol emissions public :: aero_model_surfarea ! tropospheric aerosol wet surface area for chemistry - public :: aero_model_strat_surfarea ! stub + public :: aero_model_strat_surfarea ! Misc private data character(len=32), allocatable :: fieldname(:) ! names for interstitial output fields @@ -76,6 +79,10 @@ module aero_model real(r8) :: sol_factb_interstitial = 0.1_r8 real(r8) :: sol_factic_interstitial = 0.4_r8 + integer, parameter :: max_sad_spec = 16 + character(len=32) :: sad_chem_spec_types(max_sad_spec) = ' ' + character(len=32) :: sad_strat_spec_types(max_sad_spec) = ' ' + logical :: convproc_do_aer class(aerosol_properties), pointer :: aero_props =>null() @@ -100,7 +107,8 @@ subroutine aero_model_readnl(nlfile) character(len=*), parameter :: subname = 'aero_model_readnl' ! Namelist variables - namelist /aerosol_nl/ sol_facti_cloud_borne, sol_factb_interstitial, sol_factic_interstitial + namelist /aerosol_nl/ sol_facti_cloud_borne, sol_factb_interstitial, sol_factic_interstitial, & + sad_chem_spec_types, sad_strat_spec_types !----------------------------------------------------------------------------- @@ -124,6 +132,8 @@ subroutine aero_model_readnl(nlfile) 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(sad_chem_spec_types, len(sad_chem_spec_types(1))*max_sad_spec, mpichar, 0, mpicom) + call mpibcast(sad_strat_spec_types, len(sad_strat_spec_types(1))*max_sad_spec, mpichar, 0, mpicom) #endif call aero_wetdep_readnl(nlfile) @@ -356,6 +366,21 @@ subroutine aero_model_init( pbuf2d ) call aero_wetdep_init() + if (masterproc) then + write(iulog,*) 'SAD chemistry spec_types:' + do l = 1, max_sad_spec + if (len_trim(sad_chem_spec_types(l)) > 0) then + write(iulog,*) ' ', trim(sad_chem_spec_types(l)) + end if + end do + write(iulog,*) 'SAD stratospheric spec_types:' + do l = 1, max_sad_spec + if (len_trim(sad_strat_spec_types(l)) > 0) then + write(iulog,*) ' ', trim(sad_strat_spec_types(l)) + end if + end do + end if + end subroutine aero_model_init !============================================================================= @@ -398,24 +423,15 @@ end subroutine aero_model_wetdep ! called from mo_usrrxt !------------------------------------------------------------------------- subroutine aero_model_surfarea( & - state, mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, m, ltrop, & - dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop, sad_ssa ) + state, relhum, pmid, temp, ltrop, & + sfc, dm_aer, sad_trop, reff_trop, sad_ssa ) ! dummy args type(physics_state), intent(in) :: state ! Physics state variables real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) - real(r8), intent(in) :: mmr(:,:,:) - real(r8), intent(in) :: radmean ! mean radii in cm - real(r8), intent(in) :: strato_sad(:,:) - integer, intent(in) :: ncol integer, intent(in) :: ltrop(:) - real(r8), intent(in) :: dlat(:) ! degrees latitude - integer, intent(in) :: het1_ndx real(r8), intent(in) :: relhum(:,:) - real(r8), intent(in) :: m(:,:) ! 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(:,:,:) @@ -424,14 +440,22 @@ subroutine aero_model_surfarea( & real(r8), intent(out) :: sad_ssa(:,:) ! local vars - integer :: beglev(ncol) - integer :: endlev(ncol) + integer :: beglev(pcols) + integer :: endlev(pcols) + + integer :: i,k, lchnk, ncol + + class(aerosol_state), pointer :: aero_state sad_ssa = -huge(1._r8) + lchnk = state%lchnk + ncol = state%ncol beglev(:ncol)=ltrop(:ncol)+1 endlev(:ncol)=pver - call surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, sad_trop, reff_trop, sfc=sfc, dm_aer=dm_aer ) + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) + call aero_state%surf_area_dens(aero_props, sad_chem_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, sad_trop, reff_trop, sfc, dm_aer ) end subroutine aero_model_surfarea @@ -439,29 +463,40 @@ end subroutine aero_model_surfarea ! provides wet stratospheric aerosol surface area info for sectional aerosols ! called from mo_gas_phase_chemdr.F90 !------------------------------------------------------------------------- - subroutine aero_model_strat_surfarea( state, ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat ) + subroutine aero_model_strat_surfarea( state, pmid, temp, ltrop, strato_sad, reff_strat ) use ref_pres, only: clim_modal_aero_top_lev ! dummy args type(physics_state), intent(in) :: state ! Physics state variables - integer, intent(in) :: ncol - real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) integer, intent(in) :: ltrop(:) ! tropopause level indices - type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(out) :: strato_sad(:,:) ! aerosol surface area density (cm2/cm3), zeroed below the tropopause real(r8), intent(out) :: reff_strat(:,:) ! aerosol effective radius (cm), zeroed below the tropopause ! local vars - integer :: beglev(ncol) - integer :: endlev(ncol) + integer :: i,k, lchnk, ncol + + real(r8) :: sfc_tmp(pcols,pver,nbins) + real(r8) :: dm_tmp(pcols,pver,nbins) + real(r8) :: relhum(pcols,pver) + + class(aerosol_state), pointer :: aero_state + integer :: beglev(pcols) + integer :: endlev(pcols) + + lchnk = state%lchnk + ncol = state%ncol beglev(:ncol) = clim_modal_aero_top_lev endlev(:ncol) = ltrop(:ncol) - call surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, strato_sad, reff_strat ) + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) + if (len_trim(sad_strat_spec_types(1)) > 0) then + call aero_state%surf_area_dens(aero_props, sad_strat_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, strato_sad, reff_strat) + end if end subroutine aero_model_strat_surfarea @@ -758,117 +793,8 @@ end subroutine aero_model_emissions !=============================================================================== ! private methods - !============================================================================= !============================================================================= - subroutine surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, sad, reff, sfc, dm_aer ) - use mo_constants, only: pi - use carma_intr, only: carma_effecitive_radius - - ! dummy args - type(physics_state), intent(in) :: state ! Physics state variables - type(physics_buffer_desc), pointer :: pbuf(:) - integer, intent(in) :: ncol - real(r8), intent(in) :: mmr(:,:,:) - integer, intent(in) :: beglev(:) - integer, intent(in) :: endlev(:) - real(r8), intent(out) :: sad(:,:) ! bulk surface area density in cm2/cm3 from beglev to endlev, zero elsewhere - real(r8), intent(out) :: reff(:,:) ! bulk effective radius in cm from beglev to endlev, zero elsewhere - real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin - real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin - - ! local vars - real(r8) :: reffaer(pcols,pver) ! bulk effective radius in cm - - real(r8) :: sad_bin(pcols,pver,nbins) - integer :: icol, ilev, ibin, ispec !!, reff_pbf_ndx - real(r8) :: chm_mass, tot_mass - character(len=32) :: spectype - real(r8) :: wetr(pcols,pver) ! CARMA bin wet radius in cm - real(r8) :: wetrho(pcols,pver) ! CARMA bin wet density - real(r8) :: sad_carma(pcols,pver) ! CARMA bin wet surface area density in cm2/cm3 - real(r8), pointer :: aer_bin_mmr(:,:) - - character(len=aero_name_len) :: bin_name, shortname - integer :: igroup, indxbin, rc, nchr - - sad = 0._r8 - reff = 0._r8 - - ! - ! Compute surface aero for each bin. - ! Total over all bins as the surface area for chemical reactions. - ! - - reffaer = carma_effecitive_radius(state) - - sad = 0._r8 - sad_bin = 0._r8 - reff = 0._r8 - - do ibin=1,nbins ! loop over aerosol bins - call rad_aer_get_info_by_bin(0, ibin, bin_name=bin_name) - - nchr = len_trim(bin_name)-2 - shortname = bin_name(:nchr) - - call carma_get_group_by_name(shortname, igroup, rc) - - read(bin_name(nchr+1:),*) indxbin - - call carma_get_wet_radius(state, igroup, indxbin, wetr, wetrho, rc) ! m - wetr(:ncol,:) = wetr(:ncol,:) * 1.e2_r8 ! cm - call carma_get_sad(state, igroup, indxbin, sad_carma, rc) - - if (present(dm_aer)) then - dm_aer(:ncol,:,ibin) = 2._r8 * wetr(:ncol,:) ! convert wet radius (cm) to wet diameter (cm) - endif - sad_bin(:ncol,:,ibin) = sad_carma(:ncol,:) ! cm^2/cm^3 - end do - - do icol = 1,ncol - do ilev = beglev(icol),endlev(icol) - do ibin=1,nbins ! loop over aerosol bins - ! - ! compute a mass weighting of the number - ! - tot_mass = 0._r8 - chm_mass = 0._r8 - do ispec=1,nspec(ibin) - - call rad_cnst_get_bin_mmr_by_idx(0, ibin, ispec, 'a', state, pbuf, aer_bin_mmr) - - tot_mass = tot_mass + aer_bin_mmr(icol,ilev) - - call rad_aer_get_bin_props_by_idx(0, ibin, ispec, spectype=spectype) - - if ( trim(spectype) == 'sulfate' .or. & - trim(spectype) == 's-organic' .or. & - trim(spectype) == 'p-organic' .or. & - trim(spectype) == 'black-c' .or. & - trim(spectype) == 'ammonium') then - chm_mass = chm_mass + aer_bin_mmr(icol,ilev) - end if - - end do - if ( tot_mass > 0._r8 ) then - ! surface area density - sad_bin(icol,ilev,ibin) = chm_mass / tot_mass * sad_bin(icol,ilev,ibin) ! cm^2/cm^3 - else - sad_bin(icol,ilev,ibin) = 0._r8 - end if - end do - sad(icol,ilev) = sum(sad_bin(icol,ilev,:)) - reff(icol,ilev) = reffaer(icol,ilev) - - end do - end do - - if (present(sfc)) then - sfc(:,:,:) = sad_bin(:,:,:) - endif - - end subroutine surf_area_dens !============================================================================= subroutine mmr2vmr_carma(lchnk, vmr, mbar, mw_carma, ncol, im, rmass) diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index 0768e57c68..bab25f289a 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -85,9 +85,6 @@ module aero_model ! for surf_area_dens integer,allocatable :: num_idx(:) - integer,allocatable :: index_tot_mass(:,:) - integer,allocatable :: index_chm_mass(:,:) - integer,allocatable :: index_ssa_mass(:,:) integer :: ndx_h2so4 character(len=fieldname_len), allocatable :: dgnum_name(:), dgnumwet_name(:) @@ -96,6 +93,16 @@ module aero_model character(len=16) :: drydep_list(pcnst) = ' ' real(r8) :: seasalt_emis_scale + integer, parameter :: max_sad_spec = 16 + character(len=32) :: sad_chem_spec_types(max_sad_spec) = ' ' + character(len=32) :: sad_seasalt_spec_types(max_sad_spec) = ' ' + character(len=32) :: sad_strat_spec_types(max_sad_spec) = ' ' + + ! Mode types excluded from SAD: primary_carbon should not contribute + integer, parameter :: num_sad_exclude_modes = 1 + character(len=32), parameter :: sad_exclude_mode_types(num_sad_exclude_modes) = (/ & + 'primary_carbon ' /) + integer :: ndrydep = 0 integer,allocatable :: drydep_indices(:) logical :: drydep_lq(pcnst) @@ -132,7 +139,8 @@ subroutine aero_model_readnl(nlfile) ! Namelist variables character(len=16) :: aer_drydep_list(pcnst) = ' ' - namelist /aerosol_nl/ aer_drydep_list, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale + namelist /aerosol_nl/ aer_drydep_list, modal_strat_sulfate, modal_accum_coarse_exch, seasalt_emis_scale, & + sad_chem_spec_types, sad_seasalt_spec_types, sad_strat_spec_types !----------------------------------------------------------------------------- @@ -157,6 +165,9 @@ subroutine aero_model_readnl(nlfile) call mpibcast(modal_strat_sulfate, 1, mpilog, 0, mpicom) call mpibcast(seasalt_emis_scale, 1, mpir8, 0, mpicom) call mpibcast(modal_accum_coarse_exch, 1, mpilog, 0, mpicom) + call mpibcast(sad_chem_spec_types, len(sad_chem_spec_types(1))*max_sad_spec, mpichar, 0, mpicom) + call mpibcast(sad_seasalt_spec_types, len(sad_seasalt_spec_types(1))*max_sad_spec, mpichar, 0, mpicom) + call mpibcast(sad_strat_spec_types, len(sad_strat_spec_types(1))*max_sad_spec, mpichar, 0, mpicom) #endif drydep_list = aer_drydep_list @@ -197,6 +208,7 @@ subroutine aero_model_init( pbuf2d ) use modal_aero_gasaerexch, only: modal_aero_gasaerexch_init use modal_aero_newnuc, only: modal_aero_newnuc_init use modal_aero_rename, only: modal_aero_rename_init + use aerosol_spec_utils, only: spec_type_in_list ! args type(physics_buffer_desc), pointer :: pbuf2d(:,:) @@ -517,40 +529,35 @@ subroutine aero_model_init( pbuf2d ) endif end do - allocate(index_tot_mass(nmodes,nspec_max)) - allocate(index_chm_mass(nmodes,nspec_max)) - index_tot_mass = -1 - index_chm_mass = -1 - allocate(index_ssa_mass(nmodes,nspec_max)) - index_ssa_mass = -1 - - ! for surf_area_dens - ! define indices associated with the various aerosol types + ! determine coarse dust mode number do n = 1,nmodes call rad_aer_get_info(0, n, mode_type=mode_type, nspec=nspec) - if ( trim(mode_type) /= 'primary_carbon') then ! ignore the primary_carbon mode - do l = 1, nspec - call rad_aer_get_info(0, n, l, spec_type=spec_type, spec_name=spec_name) - index_tot_mass(n,l) = get_spc_ndx(spec_name) - if ( trim(spec_type) == 'sulfate' .or. & - trim(spec_type) == 's-organic' .or. & - trim(spec_type) == 'p-organic' .or. & - trim(spec_type) == 'black-c' .or. & - trim(spec_type) == 'ammonium') then - index_chm_mass(n,l) = get_spc_ndx(spec_name) - endif - if ( trim(spec_type) == 'seasalt') then - index_ssa_mass(n,l) = get_spc_ndx(spec_name) - endif - enddo - endif - - ! determine coarse dust mode number if (mode_type=='coarse' .or. mode_type=='coarse_dust') then n_coarse_dust = n end if enddo + if (masterproc) then + write(iulog,*) 'SAD chemistry spec_types:' + do l = 1, max_sad_spec + if (len_trim(sad_chem_spec_types(l)) > 0) then + write(iulog,*) ' ', trim(sad_chem_spec_types(l)) + end if + end do + write(iulog,*) 'SAD seasalt spec_types:' + do l = 1, max_sad_spec + if (len_trim(sad_seasalt_spec_types(l)) > 0) then + write(iulog,*) ' ', trim(sad_seasalt_spec_types(l)) + end if + end do + write(iulog,*) 'SAD stratospheric spec_types:' + do l = 1, max_sad_spec + if (len_trim(sad_strat_spec_types(l)) > 0) then + write(iulog,*) ' ', trim(sad_strat_spec_types(l)) + end if + end do + end if + if (has_sox) then do m = 1, ntot_amode @@ -917,24 +924,15 @@ end subroutine aero_model_wetdep ! called from mo_usrrxt !------------------------------------------------------------------------- subroutine aero_model_surfarea( & - state, mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, rho, ltrop, & - dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop, sad_ssa ) + state, relhum, pmid, temp, ltrop, & + sfc, dm_aer, sad_trop, reff_trop, sad_ssa ) ! dummy args type(physics_state), intent(in) :: state ! Physics state variables real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) - real(r8), intent(in) :: mmr(:,:,:) - real(r8), intent(in) :: radmean ! mean radii in cm - real(r8), intent(in) :: strato_sad(:,:) - integer, intent(in) :: ncol integer, intent(in) :: ltrop(:) - real(r8), intent(in) :: dlat(:) ! degrees latitude - integer, intent(in) :: het1_ndx real(r8), intent(in) :: relhum(:,:) - real(r8), intent(in) :: rho(:,:) ! total atm density (/cm^3) - real(r8), intent(in) :: sulfate(:,:) - type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(inout) :: sfc(:,:,:) real(r8), intent(inout) :: dm_aer(:,:,:) @@ -943,22 +941,28 @@ subroutine aero_model_surfarea( & real(r8), intent(out) :: sad_ssa(:,:) ! local vars - real(r8), pointer, dimension(:,:,:) :: dgnumwet - integer :: beglev(ncol) - integer :: endlev(ncol) - integer :: i,k + integer :: beglev(pcols) + integer :: endlev(pcols) + integer :: i,k, lchnk, ncol + real(r8) :: reff_ssa(pcols,pver) + + class(aerosol_state), pointer :: aero_state - call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) + lchnk = state%lchnk + ncol = state%ncol beglev(:ncol)=ltrop(:ncol)+1 endlev(:ncol)=pver - call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, sad_trop, reff_trop, sfc=sfc, sad_ssa=sad_ssa ) - do i = 1,ncol - do k = ltrop(i)+1,pver - dm_aer(i,k,:) = dgnumwet(i,k,:) * 1.e2_r8 ! convert m to cm - enddo - enddo + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) + + if (len_trim(sad_seasalt_spec_types(1))>0) then + call aero_state%surf_area_dens(aero_props, sad_seasalt_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, sad_ssa, reff_ssa ) + end if + + call aero_state%surf_area_dens(aero_props, sad_chem_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, sad_trop, reff_trop, sfc, dm_aer ) end subroutine aero_model_surfarea @@ -966,34 +970,40 @@ end subroutine aero_model_surfarea ! provides WET stratospheric aerosol surface area info for modal aerosols ! if modal_strat_sulfate = TRUE -- called from mo_gas_phase_chemdr !------------------------------------------------------------------------- - subroutine aero_model_strat_surfarea( state, ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat ) + subroutine aero_model_strat_surfarea( state, pmid, temp, ltrop, strato_sad, reff_strat ) ! dummy args type(physics_state), intent(in) :: state ! Physics state variables - integer, intent(in) :: ncol - real(r8), intent(in) :: mmr(:,:,:) real(r8), intent(in) :: pmid(:,:) real(r8), intent(in) :: temp(:,:) integer, intent(in) :: ltrop(:) ! tropopause level indices - type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(out) :: strato_sad(:,:) real(r8), intent(out) :: reff_strat(:,:) ! local vars - real(r8), pointer, dimension(:,:,:) :: dgnumwet - integer :: beglev(ncol) - integer :: endlev(ncol) + integer :: i,k, lchnk, ncol + + real(r8) :: relhum(state%ncol,pver) + + class(aerosol_state), pointer :: aero_state + integer :: beglev(pcols) + integer :: endlev(pcols) reff_strat = 0._r8 strato_sad = 0._r8 if (.not.modal_strat_sulfate) return - call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) + relhum = huge(1._r8) - beglev(:ncol)=top_lev - endlev(:ncol)=ltrop(:ncol) - call surf_area_dens( ncol, mmr, pmid, temp, dgnumwet, beglev, endlev, strato_sad, reff_strat ) + lchnk = state%lchnk + ncol = state%ncol + beglev(:ncol) = top_lev + endlev(:ncol) = ltrop(:ncol) + + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) + call aero_state%surf_area_dens(aero_props, sad_strat_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, strato_sad, reff_strat) end subroutine aero_model_strat_surfarea @@ -1327,112 +1337,6 @@ end subroutine aero_model_emissions !=============================================================================== ! private methods - - !============================================================================= - !============================================================================= - subroutine surf_area_dens( ncol, mmr, pmid, temp, diam, beglev, endlev, sad, reff, sfc, sad_ssa ) - use mo_constants, only : pi - use modal_aero_data, only : nspec_amode, alnsg_amode - - ! dummy args - integer, intent(in) :: ncol - real(r8), intent(in) :: mmr(:,:,:) - real(r8), intent(in) :: pmid(:,:) - real(r8), intent(in) :: temp(:,:) - real(r8), intent(in) :: diam(:,:,:) - integer, intent(in) :: beglev(:) - integer, intent(in) :: endlev(:) - real(r8), intent(out) :: sad(:,:) - real(r8), intent(out) :: reff(:,:) - real(r8),optional, intent(out) :: sfc(:,:,:) - real(r8),optional, intent(out) :: sad_ssa(:,:) - - ! local vars - real(r8) :: sad_mode(pcols,pver,ntot_amode),radeff(pcols,pver) - real(r8) :: vol(pcols,pver),vol_mode(pcols,pver,ntot_amode) - real(r8) :: rho_air - integer :: i,k,l,m - real(r8) :: chm_mass, tot_mass - real(r8) :: ssa_mass - real(r8) :: sad_mode_ssa(pcols,pver,ntot_amode) - - ! - ! Compute surface aero for each mode. - ! Total over all modes as the surface area for chemical reactions. - ! - - sad = 0._r8 - sad_mode = 0._r8 - vol = 0._r8 - vol_mode = 0._r8 - reff = 0._r8 - if (present(sad_ssa)) then - sad_ssa = 0._r8 - sad_mode_ssa = 0._r8 - end if - - do i = 1,ncol - do k = beglev(i),endlev(i) - rho_air = pmid(i,k)/(temp(i,k)*287.04_r8) - do l=1,ntot_amode - ! - ! compute a mass weighting of the number - ! - tot_mass = 0._r8 - chm_mass = 0._r8 - ssa_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)) - if (present(sad_ssa)) then - if ( index_ssa_mass(l,m) > 0 ) & - ssa_mass = ssa_mass + mmr(i,k,index_ssa_mass(l,m)) - end if - end do - if ( tot_mass > 0._r8 ) then - ! surface area density - sad_mode(i,k,l) = chm_mass/tot_mass & - * mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2._r8 & - * exp(2._r8*alnsg_amode(l)**2._r8) ! m^2/m^3 - sad_mode(i,k,l) = 1.e-2_r8 * sad_mode(i,k,l) ! cm^2/cm^3 - - if (present(sad_ssa)) then - sad_mode_ssa(i,k,l) = ssa_mass/tot_mass & - * mmr(i,k,num_idx(l))*rho_air*pi*diam(i,k,l)**2._r8 & - * exp(2._r8*alnsg_amode(l)**2._r8) ! m^2/m^3 - sad_mode_ssa(i,k,l) = 1.e-2_r8 * sad_mode_ssa(i,k,l) ! cm^2/cm^3 - end if - - ! volume calculation, for use in effective radius calculation - vol_mode(i,k,l) = chm_mass/tot_mass & - * mmr(i,k,num_idx(l))*rho_air*pi/6._r8*diam(i,k,l)**3._r8 & - * exp(4.5_r8*alnsg_amode(l)**2._r8) ! m^3/m^3 = cm^3/cm^3 - else - sad_mode(i,k,l) = 0._r8 - vol_mode(i,k,l) = 0._r8 - if (present(sad_ssa)) then - sad_mode_ssa(i,k,l) = 0._r8 - end if - end if - end do - sad(i,k) = sum(sad_mode(i,k,:)) - vol(i,k) = sum(vol_mode(i,k,:)) - reff(i,k) = 3._r8*vol(i,k)/sad(i,k) - if (present(sad_ssa)) then - sad_ssa(i,k) = sum(sad_mode_ssa(i,k,:)) - end if - - enddo - enddo - - if (present(sfc)) then - sfc(:,:,:) = sad_mode(:,:,:) - endif - - end subroutine surf_area_dens - !=============================================================================== !=============================================================================== subroutine modal_aero_bcscavcoef_init diff --git a/src/chemistry/mozart/mo_gas_phase_chemdr.F90 b/src/chemistry/mozart/mo_gas_phase_chemdr.F90 index c071951218..f2a4e531b6 100644 --- a/src/chemistry/mozart/mo_gas_phase_chemdr.F90 +++ b/src/chemistry/mozart/mo_gas_phase_chemdr.F90 @@ -674,7 +674,7 @@ subroutine gas_phase_chemdr(state, lchnk, ncol, imozart, q, & strato_sad(:,:) = 0._r8 ! Prognostic modal stratospheric sulfate: compute dry strato_sad - call aero_model_strat_surfarea( state, ncol, mmr, pmid, tfld, troplevchem, pbuf, strato_sad, reff_strat ) + call aero_model_strat_surfarea( state, pmid, tfld, troplevchem, strato_sad, reff_strat ) endif diff --git a/src/chemistry/mozart/mo_usrrxt.F90 b/src/chemistry/mozart/mo_usrrxt.F90 index 432c4038f1..08b20f2209 100644 --- a/src/chemistry/mozart/mo_usrrxt.F90 +++ b/src/chemistry/mozart/mo_usrrxt.F90 @@ -1109,8 +1109,8 @@ subroutine usrrxt( state, rxt, temp, tempi, tempe, invariants, h2ovmr, sad_trop(:ncol,:pver)=strato_sad(:ncol,:pver) else call aero_model_surfarea( & - state, mmr, rm1, relhum, pmid, temp, strato_sad, sulfate, m, tropchemlev, dlat, & - het1_ndx, pbuf, ncol, sfc_array, dm_array, sad_trop, reff_trop, sad_sslt ) + state, relhum, pmid, temp, tropchemlev, & + sfc_array, dm_array, sad_trop, reff_trop, sad_sslt ) endif endif