Aerosol surface area densities refactor

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

bld/namelist_files/namelist_definition.xml

@@ -6708,6 +6708,37 @@ Tuning for seasalt_emis
 Default: set by build-namelist.
 </entry>
 
+<entry id="sad_chem_spec_types" type="char*32(16)" category="cam_chem"
+       group="aerosol_nl" valid_values="" >
+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'
+</entry>
+
+<entry id="sad_seasalt_spec_types" type="char*32(16)" category="cam_chem"
+       group="aerosol_nl" valid_values="" >
+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'
+</entry>
+
+<entry id="sad_strat_spec_types" type="char*32(16)" category="cam_chem"
+       group="aerosol_nl" valid_values="" >
+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'
+</entry>
+
 <entry id="aer_sol_facti" type="real(1000)" category="cam_chem"
        group="aerosol_nl" valid_values="" >
 In-cloud solubility factor used in BULK aerosol wet removal

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

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

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).