Generalize cloud chemistry

src/chemistry/aerosol/aero_wetdep_cam.F90

@@ -309,12 +309,12 @@ subroutine add_hist_fields(name,baseunits)
          call add_default (trim(name)//'SFWET', 1, ' ')
       endif
       if ( history_aerosol ) then
-         call add_default (trim(name)//'SFSEC', 1, ' ')
          call add_default (trim(name)//'SFSIC', 1, ' ')
          call add_default (trim(name)//'SFSIS', 1, ' ')
          call add_default (trim(name)//'SFSBC', 1, ' ')
          call add_default (trim(name)//'SFSBS', 1, ' ')
          if (convproc_do_aer) then
+            call add_default (trim(name)//'SFSEC', 1, ' ')
             call add_default (trim(name)//'SFSES', 1, ' ')
             call add_default (trim(name)//'SFSBD', 1, ' ')
          end if

src/chemistry/aerosol/aerosol_properties_mod.F90

@@ -109,7 +109,7 @@ end function aero_number_transported
      !  long wave species refractive indices
      !  species morphology
      !------------------------------------------------------------------------
-     subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, &
+     subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, spec_mw, &
           spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, &
           dryrad)
        import :: aerosol_properties, r8
@@ -118,6 +118,7 @@ subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, &
        integer, intent(in) :: species_ndx         ! species index
        real(r8), optional, intent(out) :: density ! density (kg/m3)
        real(r8), optional, intent(out) :: hygro   ! hygroscopicity
+       real(r8), optional, intent(out) :: spec_mw ! species molecular weight
        character(len=*), optional, intent(out) :: spectype  ! species type
        character(len=*), optional, intent(out) :: specname  ! species name
        character(len=*), optional, intent(out) :: specmorph ! species morphology

src/chemistry/aerosol/aerosol_state_mod.F90

@@ -61,6 +61,8 @@ module aerosol_state_mod
      procedure(aero_wet_diam), deferred :: wet_diameter
      procedure :: convcld_actfrac
      procedure :: sol_factb_interstitial
+     procedure(aero_aqu_gain_binfraction), deferred :: aqu_gain_binfraction
+
   end type aerosol_state
 
   ! for state fields
@@ -279,6 +281,21 @@ function aero_wet_diam(self, bin_idx, ncol, nlev) result(diam)
 
      end function aero_wet_diam
 
+     !------------------------------------------------------------------------------
+     ! aqueous chemistry partitioning -- used in sox_cldaero_update
+     !------------------------------------------------------------------------------
+     subroutine aero_aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqgain)
+       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) :: type                ! aerosol species type
+       real(r8), intent(in) :: qcw(:,:,:)                  ! cloud-borne aerosol volume mixing ratio
+       real(r8), intent(in) :: delso4_o3rxn(:,:)           ! sulfate concentration change due to oxidation
+       real(r8), intent(out) :: faqgain(:,:,:)             ! fraction gain in each mode / bin
+
+     end subroutine aero_aqu_gain_binfraction
+
   end interface
 
 contains

src/chemistry/aerosol/bulk_aerosol_properties_mod.F90

@@ -152,7 +152,7 @@ end function number_transported
   !  long wave species refractive indices
   !  species morphology
   !------------------------------------------------------------------------
-  subroutine get(self, bin_ndx, species_ndx, density, hygro, &
+  subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, &
                  spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, &
                  dryrad)
 
@@ -161,6 +161,7 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, &
     integer, intent(in) :: species_ndx         ! species index
     real(r8), optional, intent(out) :: density ! density (kg/m3)
     real(r8), optional, intent(out) :: hygro   ! hygroscopicity
+    real(r8), optional, intent(out) :: spec_mw ! species molecular weight
     character(len=*), optional, intent(out) :: spectype  ! species type
     character(len=*), optional, intent(out) :: specname  ! species name
     character(len=*), optional, intent(out) :: specmorph ! species morphology
@@ -217,6 +218,23 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, &
     if (present(dryrad)) then
        call rad_aer_get_props(self%list_idx_, bin_ndx,  dryrad_aer=dryrad)
     end if
+    if (present(spec_mw)) then
+       call rad_aer_get_props(self%list_idx_, bin_ndx,  aername=aername)
+
+       select case ( to_lower( aername(:4) ) )
+       case('sulf','volc')
+          spec_mw = 96._r8
+       case('bcar','bcph','ocar','ocph')
+          spec_mw = 12._r8
+       case('dust')
+          spec_mw = 12._r8 !!! ????
+       case('sslt','seas','ssam','sscm')
+          spec_mw = 57._r8
+       case default
+          spec_mw = nan
+          call endrun('ERROR: bulk_aerosol_properties_mod%get aername not recognized : '//aername)
+       end select
+    end if
 
   end subroutine get
 

src/chemistry/aerosol/bulk_aerosol_state_mod.F90

@@ -45,6 +45,7 @@ module bulk_aerosol_state_mod
      procedure :: wet_diameter
      procedure :: convcld_actfrac
      procedure :: wgtpct
+     procedure :: aqu_gain_binfraction
 
      final :: destructor
 
@@ -419,4 +420,20 @@ function wgtpct(self, ncol, nlev) result(wtp)
 
   end function wgtpct
 
+  !------------------------------------------------------------------------------
+  ! aqueous chemistry partitioning -- used in sox_cldaero_update
+  !------------------------------------------------------------------------------
+  subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqgain)
+
+    class(bulk_aerosol_state), intent(in) :: self
+    class(aerosol_properties), intent(in) :: aero_props
+    character(len=*), intent(in) :: type
+    real(r8), intent(in) :: qcw(:,:,:)
+    real(r8), intent(in) :: delso4_o3rxn(:,:)
+    real(r8), intent(out) :: faqgain(:,:,:) ! fraction gain in each mode / bin
+
+    faqgain(:,:,:) = 1._r8
+
+  end subroutine aqu_gain_binfraction
+
 end module bulk_aerosol_state_mod

src/chemistry/aerosol/carma_aerosol_properties_mod.F90

@@ -264,7 +264,7 @@ end function number_transported
   !  long wave species refractive indices
   !  species morphology
   !------------------------------------------------------------------------
-  subroutine get(self, bin_ndx, species_ndx, density, hygro, &
+  subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, &
                  spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, &
                  dryrad)
     use cam_abortutils, only: endrun
@@ -274,6 +274,7 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, &
     integer, intent(in) :: species_ndx         ! species index
     real(r8), optional, intent(out) :: density ! density (kg/m3)
     real(r8), optional, intent(out) :: hygro   ! hygroscopicity
+    real(r8), optional, intent(out) :: spec_mw ! species molecular weight
     character(len=*), optional, intent(out) :: spectype  ! species type
     character(len=*), optional, intent(out) :: specname  ! species name
     character(len=*), optional, intent(out) :: specmorph ! species morphology
@@ -282,6 +283,8 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, &
     real(r8), optional, intent(out) :: num_to_mass_aer ! ratio of number to mass concentration
     real(r8), optional, intent(out) :: dryrad  ! dry radius (m)
 
+    character(len=32) :: type
+
     if (present(density)) then
        call rad_aer_get_bin_props_by_idx(self%list_idx_, bin_ndx, species_ndx, density_aer=density)
     end if
@@ -307,6 +310,25 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, &
           call rad_aer_get_info_by_bin_spec(self%list_idx_, bin_ndx, species_ndx, spec_name=specname)
        end if
     end if
+    if (present(spec_mw)) then
+       call rad_aer_get_bin_props_by_idx(self%list_idx_, bin_ndx, species_ndx, spectype=type)
+       select case (trim(type))
+       case('sulfate')
+          spec_mw = 96._r8
+       case('black-c')
+          spec_mw = 12._r8
+       case('s-organic')
+          spec_mw = 250._r8
+       case('p-organic')
+          spec_mw = 12._r8
+       case('dust')
+          spec_mw = 12._r8
+       case('seasalt')
+          spec_mw = 57._r8
+       case default
+          spec_mw = nan
+       end select
+    end if
 
     if (present(num_to_mass_aer)) then
        ! num_to_mass_aer for sectional aerosols should not be read from file

src/chemistry/aerosol/carma_aerosol_state_mod.F90

@@ -52,6 +52,7 @@ module carma_aerosol_state_mod
      procedure :: wet_volume
      procedure :: water_volume
      procedure :: wet_diameter
+     procedure :: aqu_gain_binfraction
 
      final :: destructor
 
@@ -605,4 +606,78 @@ function wet_diameter(self, bin_idx, ncol, nlev) result(diam)
 
   end function wet_diameter
 
+  !------------------------------------------------------------------------------
+  ! aqueous chemistry partitioning -- used in sox_cldaero_update
+  !------------------------------------------------------------------------------
+  subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqgain)
+
+    class(carma_aerosol_state), intent(in) :: self
+    class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object
+    character(len=*), intent(in) :: type                ! aerosol species type
+    real(r8), intent(in) :: qcw(:,:,:)                  ! cloud-borne aerosol volume mixing ratio
+    real(r8), intent(in) :: delso4_o3rxn(:,:)           ! sulfate concentration change due to oxidation
+    real(r8), intent(out) :: faqgain(:,:,:)             ! fraction gain in each mode / bin
+
+    real(r8) :: raddry(pcols,pver)   ! dry radius (m)
+    real(r8) :: rhodry(pcols,pver)   ! dry density (kg/m3)
+    character(len=aero_name_len) :: bin_name, shortname
+    integer :: igroup, ibin, jbin, rc, nchr
+    integer :: icol, klev, nbins, ncol
+    real(r8), allocatable :: rad_cm(:,:,:)
+    real(r8), allocatable :: wt_mass(:)
+    real(r8) :: wt_sum
+
+    ! To calculate the fraction of sulfate mass produced by aq.chemistry that will be
+    ! distributed across different bins, we need to conserve particle number. In CARMA,
+    ! this is done following the formula for mass transfer for, following
+    ! Yu et al., 2015, A3.2:  https://doi.org/10.1002/2014MS000421. The mass fraction
+    ! for each bin can be calculated by normalizing the mass transfer rate in each bin:
+    ! M/D^2 as: fra = Mi/D^2/sum(Mi/D^2).
+
+    ncol = self%state%ncol
+    nbins = aero_props%nbins()
+    faqgain(:,:,:) = 0._r8
+
+    allocate(wt_mass(nbins))
+    allocate(rad_cm(nbins,pcols,pver))
+    rad_cm(:,:,:) = 0._r8
+
+    do ibin = 1, nbins
+       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:),*) jbin
+       call carma_get_dry_radius(self%state, igroup, jbin, raddry, rhodry, rc)
+       if (index(bin_name,'MXAER')>0) then
+          rad_cm(ibin,:ncol,:) = raddry(:ncol,:)*1.e2_r8 ! m -> cm
+       end if
+    end do
+
+    lev_loop: do klev = 1,pver
+       col_loop: do icol = 1,ncol
+
+          !faqgain = fraction of total so4_c gain going to bin n
+          wt_sum = 0._r8
+          wt_mass(:) = 0._r8
+
+          do ibin = 1, nbins
+             if (rad_cm(ibin,icol,klev) > 0._r8) then
+                wt_mass(ibin) = delso4_o3rxn(icol,klev) / rad_cm(ibin,icol,klev) / rad_cm(ibin,icol,klev)
+                wt_sum = wt_sum + wt_mass(ibin)
+             end if
+          end do
+          do ibin = 1, nbins
+             if (wt_mass(ibin) > 0._r8) then
+                faqgain(ibin,icol,klev) = wt_mass(ibin)/wt_sum
+             end if
+          end do
+
+       end do col_loop
+    end do lev_loop
+
+    deallocate(rad_cm, wt_mass)
+
+  end subroutine aqu_gain_binfraction
+
 end module carma_aerosol_state_mod