diff --git a/bld/build-namelist b/bld/build-namelist index faad098c53..e670c6ca3c 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -3898,6 +3898,34 @@ 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'"); +} +# ...for prescribed BAM with chemistry (WACCM-4) sad_chem_spec_types should +# be empty as sulfate is from supplemental SAD from sulf_file, not here. + # 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 3eeaaf3499..5ff967b2c4 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -6784,6 +6784,43 @@ 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). + +In bulk (BAM) configurations, types not in rad_climate are provided by +a supplemental path in aero_model_surfarea. + +Default: 'sulfate','s-organic','p-organic','black-c','ammonium' (MAM, CARMA); +'sulfate','black-c','p-organic','s-organic','nitrate' (prognostic bulk); +UNSET (prescribed bulk) + + + +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..2242717cd8 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,32 @@ 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, pi, 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(in) :: pi ! pi mathematical constant + + 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 +994,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..2c884c521e 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,257 @@ 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, pi, sad, reff, sfc, dm_aer) + use aerosol_spec_utils, only : spec_type_in_list + + 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(in) :: pi ! pi mathematical constant + + 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 + 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 /) + + ! 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) + character(len=32), allocatable :: spectypes(:,:) + real(r8), pointer :: mmr(:,:) ! interstitial aerosol mass, number mixing ratios + integer :: ispc, ibin, ndx, astat + + mmr => null() + + nspec_max = 1 ! 1 species per bin in bulk representation + nbins = aero_props%nbins() + allocate(mmr_ptr(nbins,nspec_max), spectypes(nbins,nspec_max), stat=astat) + if( astat/= 0 ) call endrun('bulk_aerosol_state_mod%surf_area_dens: mmr_ptr,spectypes allocate error') + + spectypes(:,:) = ' ' + + do ibin = 1, aero_props%nbins() + do ispc = 1, aero_props%nspecies(ibin) + call aero_props%get(bin_ndx=ibin, species_ndx=ispc, spectype=spectypes(ibin,ispc)) + if (spec_type_in_list(spectypes(ibin,ispc),types_list)) then + call self%get_ambient_mmr(species_ndx=ispc, bin_ndx=ibin, mmr=mmr_ptr(ibin,ispc)%fld) + else + mmr_ptr(ibin,ispc)%fld => null() + end if + end do + end do + + 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 + + ndx = 0 + + ver_loop: do k = 1,nlev + col_loop: do i = 1,ncol + if (k < beglev(i) .or. k > endlev(i)) cycle col_loop + !------------------------------------------------------------------------- + ! ... 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) + + if (.not. associated(mmr_ptr(ibin,ispc)%fld)) cycle + + mmr => mmr_ptr(ibin,ispc)%fld + + select case ( trim(spectypes(ibin,ispc)) ) + 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) + deallocate(mmr_ptr, spectypes) + + 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..0f2ec81c57 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,145 @@ 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, pi, 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(in) :: pi ! pi mathematical constant + + 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..aecbaf8f41 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, pi, sad, reff, sfc, dm_aer) + 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(in) :: pi ! pi mathematical constant + + 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..4d05083a68 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 @@ -40,8 +41,8 @@ module aero_model ! Misc private data - integer :: so4_ndx, cb2_ndx, oc2_ndx, nit_ndx - integer :: soa_ndx, soai_ndx, soam_ndx, soab_ndx, soat_ndx, soax_ndx + integer :: so4_ndx, nit_idx + integer :: soa_idx, soai_idx, soam_idx, soab_idx, soat_idx, soax_idx ! Namelist variables character(len=16) :: wetdep_list(pcnst) = ' ' @@ -63,6 +64,14 @@ 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) = ' ' + + ! sfc/dm_aer slots reserved beyond the bulk aerosol bins for + ! supplemental_surf_area_dens (offline sulfate, nitrate, secondary + ! organics); mo_usrrxt sizes the arrays with this many extra slots + integer, parameter, public :: n_supplemental_sad = 3 + contains !============================================================================= @@ -87,6 +96,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 +122,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) @@ -392,16 +404,18 @@ subroutine aero_model_init( pbuf2d ) endif endif + ! SO4 is only tested for presence -- when it is absent the offline sulfate + ! climatology supplies the sulfate SAD so4_ndx = get_spc_ndx( 'SO4' ) - soa_ndx = get_spc_ndx( 'SOA' ) - soai_ndx = get_spc_ndx( 'SOAI' ) - soam_ndx = get_spc_ndx( 'SOAM' ) - soab_ndx = get_spc_ndx( 'SOAB' ) - soat_ndx = get_spc_ndx( 'SOAT' ) - soax_ndx = get_spc_ndx( 'SOAX' ) - cb2_ndx = get_spc_ndx( 'CB2' ) - oc2_ndx = get_spc_ndx( 'OC2' ) - nit_ndx = get_spc_ndx( 'NH4NO3' ) + + ! constituent indices used to get mmr for supplemental SADs when available. + call cnst_get_ind( 'SOA', soa_idx, abort=.false. ) + call cnst_get_ind( 'SOAI', soai_idx, abort=.false. ) + call cnst_get_ind( 'SOAM', soam_idx, abort=.false. ) + call cnst_get_ind( 'SOAB', soab_idx, abort=.false. ) + call cnst_get_ind( 'SOAT', soat_idx, abort=.false. ) + call cnst_get_ind( 'SOAX', soax_idx, abort=.false. ) + call cnst_get_ind( 'NH4NO3', nit_idx, abort=.false. ) end subroutine aero_model_init @@ -714,26 +728,17 @@ 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 + use mo_constants, only : pi, r2d ! 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(:,:,:) @@ -742,16 +747,95 @@ subroutine aero_model_surfarea( & real(r8), intent(out) :: sad_ssa(:,:) ! local vars + integer :: i,k, lchnk, ncol + integer :: beglev(pcols) + integer :: endlev(pcols) + + class(aerosol_state), pointer :: aero_state + + sad_ssa = -huge(1._r8) + + lchnk = state%lchnk + ncol = state%ncol + + aero_state => aerosol_instances_get_state(iaermod_, 0, lchnk) + + !------------------------------------------------------------------------- + ! Levels over which tropospheric SADs are computed. + ! Zero the trop SADs above 300 hPa cutoff poleward of 50 deg and above the + ! tropopause elsewhere: + !------------------------------------------------------------------------- + endlev(:ncol)=pver + do i = 1,ncol + if ( abs( state%lat(i)*r2d ) > 50._r8 ) then + beglev(i) = pver+1 + do k = 1,pver + if ( pmid(i,k) >= 30000._r8 ) then + beglev(i) = k + exit + end if + end do + else + beglev(i) = ltrop(i) + end if + end do + + call aero_state%surf_area_dens(aero_props, sad_chem_spec_types, ncol, pver, beglev, endlev, & + relhum, pmid, temp, pi, sad_total, reff_trop, sfc, dm_aer) + + ! surfaces of aerosols the rad_climate list cannot provide + ! (offline sulfate, ammonium nitrate, secondary organics) + call supplemental_surf_area_dens( state, ncol, ltrop, beglev, endlev, & + relhum, pmid, temp, sfc, dm_aer, sad_total ) + end subroutine aero_model_surfarea + + !============================================================================= + ! Supplemental surface area densities for heterogeneous chemistry, for + ! aerosols that are not represented within rad_climate right now. + ! and therefore cannot be provided by bulk_aerosol_state%surf_area_dens: + ! - WACCM4: sulfate from the offline sulf_file climatology since SO4 is not + ! in the chemical mechanism. + ! - FMOZ: ammonium nitrate and secondary organics + ! A switch here verifies 'nitrate', 's-organic' is in sad_chem_spec_types + ! for compatibility. + ! + ! After nitrate, s-organic are added to rad_climate we can retire this code + ! gradually and represent natively in the abstract aerosol interface + ! surf_area_dens: + !============================================================================= + subroutine supplemental_surf_area_dens( state, ncol, ltrop, beglev, endlev, & + relhum, pmid, temp, sfc, dm_aer, sad_total ) + + use mo_constants, only : pi, avo => avogadro, boltz_cgs + use aerosol_spec_utils, only : spec_type_in_list + use mo_sulf, only : sulf_interp + + ! dummy args + type(physics_state), intent(in) :: state ! Physics state variables + integer, intent(in) :: ncol + integer, intent(in) :: ltrop(:) ! chemical tropopause level index + integer, intent(in) :: beglev(:) ! beginning model level index + integer, intent(in) :: endlev(:) ! ending model level index + real(r8), intent(in) :: relhum(:,:) + real(r8), intent(in) :: pmid(:,:) + real(r8), intent(in) :: temp(:,:) + + real(r8), intent(inout) :: sfc(:,:,:) + real(r8), intent(inout) :: dm_aer(:,:,:) + real(r8), intent(inout) :: sad_total(:,:) + + ! local vars integer :: i,k + integer :: ndx_sulf, ndx_nit, ndx_soa + logical :: do_sulf, do_nit, do_soa + real(r8) :: sulfate(pcols,pver) ! offline sulfate vmr (mol/mol) + real(r8) :: m_air ! total atm density (/cm^3) real(r8) :: rho_air - real(r8) :: v, n, n_exp, r_rd, r_sd + real(r8) :: v, n, n_exp 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) :: dm_orgc, dm_orgc_wet, log_sd_orgc, 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 !----------------------------------------------------------------- @@ -769,29 +853,40 @@ subroutine aero_model_surfarea( & 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) + ! Pascals to dyne/cm^2, for the total density paired with the offline + ! sulfate vmr; matches the indexm invariant built by mo_setinv + real(r8), parameter :: Pa_xfac = 10._r8 + integer :: irh, rh_l, rh_u - real(r8) :: factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss - logical :: zero_aerosols + real(r8) :: factor, rfac_sulf, rfac_oc !----------------------------------------------------------------- ! ... 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 + 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 /) - sad_ssa = -huge(1._r8) + do_sulf = ( so4_ndx <= 0 ) .and. & + ( .not. spec_type_in_list('sulfate', sad_chem_spec_types) ) + do_nit = ( nit_idx > 0 ) .and. spec_type_in_list('nitrate', sad_chem_spec_types) + do_soa = ( soa_idx > 0 .or. soai_idx > 0 .or. soam_idx > 0 .or. & + soab_idx > 0 .or. soat_idx > 0 .or. soax_idx > 0 ) .and. & + spec_type_in_list('s-organic', sad_chem_spec_types) + + if (.not. (do_sulf .or. do_nit .or. do_soa)) return - 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/ + if ( do_sulf ) then + call sulf_interp( ncol, state%lchnk, sulfate ) + end if + + ! trailing slots reserved for the supplemental surfaces; slots for + ! inactive branches stay at the zero set by surf_area_dens + ndx_sulf = size(sfc,3) - n_supplemental_sad + 1 + ndx_nit = size(sfc,3) - n_supplemental_sad + 2 + ndx_soa = size(sfc,3) - n_supplemental_sad + 3 !----------------------------------------------------------------- ! ... exponent for calculating number density @@ -800,16 +895,13 @@ subroutine aero_model_surfarea( & 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 + if (k < beglev(i) .or. k > endlev(i)) cycle col_loop !------------------------------------------------------------------------- ! ... air density (kg/m3) !------------------------------------------------------------------------- @@ -820,7 +912,6 @@ subroutine aero_model_surfarea( & 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 @@ -834,159 +925,88 @@ subroutine aero_model_surfarea( & 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 - + ! offline sulfate climatology only within the chemical troposphere: + if ( do_sulf .and. k >= ltrop(i) ) then !------------------------------------------------------------------------- - ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate) + ! 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) !------------------------------------------------------------------------- - 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 + m_air = Pa_xfac * pmid(i,k) / (boltz_cgs*temp(i,k)) + v = sulfate(i,k) * m_air * 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 + + sfc(i,k,ndx_sulf) = sfc_sulf + dm_aer(i,k,ndx_sulf) = dm_sulf_wet + end if + if ( do_nit ) then !------------------------------------------------------------------------- - ! hydrophylic organic carbon (follow same procedure as sulfate) + ! ammonium nitrate (follow same procedure as sulfate, using size and density of 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 + v = state%q(i,k,nit_idx) * 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 + sfc(i,k,ndx_nit) = sfc_nit + dm_aer(i,k,ndx_nit) = dm_sulf_wet + end if + + if ( do_soa ) then !------------------------------------------------------------------------- ! secondary organic carbon (follow same procedure as sulfate) !------------------------------------------------------------------------- - if( soa_ndx > 0 ) then - v = mmr(i,k,soa_ndx) * rho_air/rho_orgc + if( soa_idx > 0 ) then + v = state%q(i,k,soa_idx) * 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 + if( soai_idx > 0 ) then + v = state%q(i,k,soai_idx) * 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 + if( soam_idx > 0 ) then + v = state%q(i,k,soam_idx) * 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 + if( soab_idx > 0 ) then + v = state%q(i,k,soab_idx) * 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 + if( soat_idx > 0 ) then + v = state%q(i,k,soat_idx) * 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 + if( soax_idx > 0 ) then + v = state%q(i,k,soax_idx) * 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 @@ -995,34 +1015,31 @@ subroutine aero_model_surfarea( & end if sfc_soa = sfc_soa + sfc_soai + sfc_soam + sfc_soab + sfc_soat + sfc_soax + sfc(i,k,ndx_soa) = sfc_soa + dm_aer(i,k,ndx_soa) = dm_orgc_wet end if - 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 /) - !------------------------------------------------------------------------- - ! ... add up total surface area density for output + ! ... add supplemental surface area densities to the total !------------------------------------------------------------------------- - sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc + sad_total(i,k) = sad_total(i,k) + sfc(i,k,ndx_sulf) + sfc(i,k,ndx_nit) + sfc(i,k,ndx_soa) enddo col_loop enddo ver_loop - end subroutine aero_model_surfarea + end subroutine supplemental_surf_area_dens + !------------------------------------------------------------------------- ! 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..dc2f57dd05 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,14 @@ 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) = ' ' + + ! sfc/dm_aer slots mo_usrrxt must reserve beyond the aerosol bins; all CARMA + ! surfaces come from the aerosol representation, so no extra slots are needed + integer, parameter, public :: n_supplemental_sad = 0 + logical :: convproc_do_aer class(aerosol_properties), pointer :: aero_props =>null() @@ -100,7 +111,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 +136,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 +370,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 +427,17 @@ 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 ) + + use mo_constants, only : pi ! 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 +446,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 :: 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, pi, sad_trop, reff_trop, sfc, dm_aer ) end subroutine aero_model_surfarea @@ -439,29 +469,41 @@ 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 + use mo_constants, only: pi ! 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, pi, strato_sad, reff_strat) + end if end subroutine aero_model_strat_surfarea @@ -758,117 +800,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..662dcca297 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,20 @@ 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) = ' ' + + ! sfc/dm_aer slots mo_usrrxt must reserve beyond the aerosol bins; all modal + ! surfaces come from the aerosol representation, so no extra slots are needed + integer, parameter, public :: n_supplemental_sad = 0 + + ! 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 +143,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 +169,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 +212,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 +533,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 +928,17 @@ 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 ) + + use mo_constants, only : pi ! 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 +947,30 @@ 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 :: lchnk, ncol + real(r8) :: reff_ssa(pcols,pver) + + class(aerosol_state), pointer :: aero_state - call pbuf_get_field(pbuf, dgnumwet_idx, dgnumwet ) + sad_ssa = -huge(1._r8) + + 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, pi, 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, pi, sad_trop, reff_trop, sfc, dm_aer ) end subroutine aero_model_surfarea @@ -966,34 +978,42 @@ 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 ) + + use mo_constants, only : pi ! 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, pi, strato_sad, reff_strat) end subroutine aero_model_strat_surfarea @@ -1327,112 +1347,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..51dd46ff73 100644 --- a/src/chemistry/mozart/mo_gas_phase_chemdr.F90 +++ b/src/chemistry/mozart/mo_gas_phase_chemdr.F90 @@ -298,7 +298,7 @@ subroutine gas_phase_chemdr(state, lchnk, ncol, imozart, q, & use mo_usrrxt, only : usrrxt, has_ice_trp_rxts use mo_setinv, only : setinv use mo_negtrc, only : negtrc - use mo_sulf, only : sulf_interp + use mo_sulf, only : sulf_interp, read_sulf use mo_setext, only : setext use fire_emissions, only : fire_emissions_vrt use mo_sethet, only : sethet @@ -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 @@ -811,6 +811,7 @@ subroutine gas_phase_chemdr(state, lchnk, ncol, imozart, q, & if ( .not. carma_hetchem_feedback ) then if( so4_ndx < 1 ) then ! get offline so4 field if not prognostic call sulf_interp( ncol, lchnk, sulfate ) + if (read_sulf) call outfld( 'SULFATE', sulfate(:ncol,:), ncol, lchnk ) else sulfate(:,:) = vmr(:,:,so4_ndx) endif diff --git a/src/chemistry/mozart/mo_sulf.F90 b/src/chemistry/mozart/mo_sulf.F90 index cb3bdcf7ba..64ebdadf61 100644 --- a/src/chemistry/mozart/mo_sulf.F90 +++ b/src/chemistry/mozart/mo_sulf.F90 @@ -21,13 +21,14 @@ module mo_sulf private public :: sulf_inti, set_sulf_time, sulf_interp, sulf_readnl + public :: read_sulf ! TRUE when offline sulfate data are available save type(trfld), pointer :: fields(:) => null() type(trfile) :: file - logical :: read_sulf = .false. + logical, protected :: read_sulf = .false. character(len=16) :: fld_name = 'SULFATE' character(len=256) :: filename = 'NONE' @@ -208,9 +209,6 @@ subroutine sulf_interp( ncol, lchnk, ccm_sulf ) !----------------------------------------------------------------------- ! ... Time interpolate sulfatei to current time !----------------------------------------------------------------------- - use cam_history, only : outfld - - implicit none !----------------------------------------------------------------------- ! ... Dummy arguments @@ -229,8 +227,6 @@ subroutine sulf_interp( ncol, lchnk, ccm_sulf ) ccm_sulf(:ncol,:) = fields(1)%data(:ncol,:,lchnk) - call outfld( 'SULFATE', ccm_sulf(:ncol,:), ncol, lchnk ) - end subroutine sulf_interp end module mo_sulf diff --git a/src/chemistry/mozart/mo_usrrxt.F90 b/src/chemistry/mozart/mo_usrrxt.F90 index 432c4038f1..01092243ea 100644 --- a/src/chemistry/mozart/mo_usrrxt.F90 +++ b/src/chemistry/mozart/mo_usrrxt.F90 @@ -796,7 +796,7 @@ subroutine usrrxt( state, rxt, temp, tempi, tempe, invariants, h2ovmr, use physics_buffer, only : physics_buffer_desc use physics_types, only : physics_state use carma_flags_mod, only : carma_hetchem_feedback - use aero_model, only : aero_model_surfarea + use aero_model, only : aero_model_surfarea, n_supplemental_sad use radiative_aerosol,only : rad_aer_get_info use time_manager, only : get_curr_calday use infnan, only : nan @@ -1049,7 +1049,7 @@ subroutine usrrxt( state, rxt, temp, tempi, tempe, invariants, h2ovmr, real(r8), parameter :: pH = 4.5e+00_r8 real(r8), pointer :: sfc(:), dm_aer(:) - integer :: ntot_amode, nbins + integer :: ntot_amode, nbins, naero real(r8), pointer :: sfc_array(:,:,:), dm_array(:,:,:) !TS2 @@ -1087,12 +1087,18 @@ subroutine usrrxt( state, rxt, temp, tempi, tempe, invariants, h2ovmr, calday = get_curr_calday() + ! one slot per aerosol mode/bin, plus the slots for supplemental SADs from species not in rad_climate + ! (bulk only: offline sulfate, ammonium nitrate, secondary organics) if (ntot_amode>0) then - allocate(sfc_array(pcols,pver,ntot_amode), dm_array(pcols,pver,ntot_amode) ) + allocate(sfc_array(pcols,pver,ntot_amode+n_supplemental_sad), & + dm_array (pcols,pver,ntot_amode+n_supplemental_sad) ) else if (nbins>0) then - allocate(sfc_array(pcols,pver,nbins), dm_array(pcols,pver,nbins) ) + allocate(sfc_array(pcols,pver,nbins+n_supplemental_sad), & + dm_array (pcols,pver,nbins+n_supplemental_sad) ) else - allocate(sfc_array(pcols,pver,5), dm_array(pcols,pver,5) ) + call rad_aer_get_info(0, naero=naero) + allocate(sfc_array(pcols,pver,naero+n_supplemental_sad), & + dm_array (pcols,pver,naero+n_supplemental_sad) ) endif sfc_array(:,:,:) = 0._r8 @@ -1109,8 +1115,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