From d745b5c2188621e996c62977775f2a816c7282ce Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 5 May 2026 13:36:05 -0600 Subject: [PATCH 1/9] merge aqueous_chem branch --- src/chemistry/aerosol/aero_wetdep_cam.F90 | 2 +- .../aerosol/aerosol_properties_mod.F90 | 3 +- src/chemistry/aerosol/aerosol_state_mod.F90 | 17 + .../aerosol/bulk_aerosol_properties_mod.F90 | 20 +- .../aerosol/bulk_aerosol_state_mod.F90 | 17 + .../aerosol/carma_aerosol_properties_mod.F90 | 24 +- .../aerosol/carma_aerosol_state_mod.F90 | 75 ++++ src/chemistry/aerosol/mo_setsox.F90 | 15 +- .../aerosol/modal_aerosol_properties_mod.F90 | 9 +- .../aerosol/modal_aerosol_state_mod.F90 | 95 +++- src/chemistry/carma_aero/aero_model.F90 | 129 ++---- .../carma_aero/carma_aero_gasaerexch.F90 | 36 +- src/chemistry/carma_aero/sox_cldaero_mod.F90 | 414 +++++++++--------- src/chemistry/modal_aero/aero_model.F90 | 145 +++--- .../modal_aero/modal_aero_gasaerexch.F90 | 50 ++- src/chemistry/modal_aero/sox_cldaero_mod.F90 | 398 ++++++++--------- 16 files changed, 826 insertions(+), 623 deletions(-) diff --git a/src/chemistry/aerosol/aero_wetdep_cam.F90 b/src/chemistry/aerosol/aero_wetdep_cam.F90 index cbb345f8d3..d41abcec94 100644 --- a/src/chemistry/aerosol/aero_wetdep_cam.F90 +++ b/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 diff --git a/src/chemistry/aerosol/aerosol_properties_mod.F90 b/src/chemistry/aerosol/aerosol_properties_mod.F90 index 0c7f70d73a..5e7642a198 100644 --- a/src/chemistry/aerosol/aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/aerosol_properties_mod.F90 @@ -110,7 +110,7 @@ end function aero_number_transported ! species morphology !------------------------------------------------------------------------ subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, & - spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) import :: aerosol_properties, r8 class(aerosol_properties), intent(in) :: self @@ -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 diff --git a/src/chemistry/aerosol/aerosol_state_mod.F90 b/src/chemistry/aerosol/aerosol_state_mod.F90 index 7fb1689b1a..c0fa318c01 100644 --- a/src/chemistry/aerosol/aerosol_state_mod.F90 +++ b/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 + 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 + + end subroutine aero_aqu_gain_binfraction + end interface contains diff --git a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 index 1ca5f3d801..c6e54c64ad 100644 --- a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 @@ -153,12 +153,13 @@ end function number_transported ! species morphology !------------------------------------------------------------------------ subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) class(bulk_aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin index integer, intent(in) :: species_ndx ! species index + real(r8), optional, intent(out) :: spec_mw ! species molecular weight real(r8), optional, intent(out) :: density ! density (kg/m3) real(r8), optional, intent(out) :: hygro ! hygroscopicity character(len=*), optional, intent(out) :: spectype ! species type @@ -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 diff --git a/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 index 38c13f372c..e599f5efc1 100644 --- a/src/chemistry/aerosol/bulk_aerosol_state_mod.F90 +++ b/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 diff --git a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 index bac1c280c4..8ee5a3ba18 100644 --- a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 @@ -265,7 +265,7 @@ end function number_transported ! species morphology !------------------------------------------------------------------------ subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + 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 diff --git a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 index e2b52f14ba..f0119ee007 100644 --- a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 +++ b/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 + 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 + + 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 diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index 7e0077adc6..2d404c5226 100644 --- a/src/chemistry/aerosol/mo_setsox.F90 +++ b/src/chemistry/aerosol/mo_setsox.F90 @@ -3,6 +3,7 @@ module mo_setsox use shr_kind_mod, only : r8 => shr_kind_r8 use cam_logfile, only : iulog use physics_types,only : physics_state + use aerosol_state_mod, only: aerosol_state implicit none @@ -148,7 +149,7 @@ end subroutine sox_inti !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine setsox( state, & + subroutine setsox( aero_state, state, & pbuf, & ncol, & lchnk, & @@ -211,6 +212,7 @@ subroutine setsox( state, & !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- + class(aerosol_state), intent(in) :: aero_state type(physics_state), intent(in) :: state ! Physics state variables type(physics_buffer_desc), pointer, intent(inout) :: pbuf(:) ! Physics buffer integer, intent(in) :: ncol ! num of columns in chunk @@ -362,7 +364,8 @@ subroutine setsox( state, & xso4(:,:) = 0._r8 xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 - + xso4_init = 0._r8 + do k = 1,pver xph(:,k) = xph0 ! initial PH value @@ -757,9 +760,9 @@ subroutine setsox( state, & / xam & ! / (molecule(a)/m3(a)) * AVOGADRO & ! * (molecule(a)/mole(a)) * M3_TO_L ! * (L(a)/m3(a)) = mole(h2o2)/mole(a)/s - + xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production - + !----------------------------------------------- ! ... Partioning !----------------------------------------------- @@ -877,7 +880,7 @@ subroutine setsox( state, & end do col_loop1 end do ver_loop1 - call sox_cldaero_update( & + call sox_cldaero_update( aero_state, & state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d ) @@ -900,7 +903,7 @@ end subroutine setsox !----------------------------------------------------------------- pure integer function get_heff_index(species_name) result(index) use shr_drydep_mod, only: species_name_table - + character(len=*), intent(in) :: species_name do index = 1, size(species_name_table) diff --git a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 index 7fc6e60345..9b2728b9c8 100644 --- a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 @@ -3,7 +3,7 @@ module modal_aerosol_properties_mod use physconst, only: pi use aerosol_properties_mod, only: aerosol_properties, aero_name_len use radiative_aerosol, only: rad_aer_get_info, rad_aer_get_mode_props, rad_aer_get_props - + use modal_aero_data, only: specmw_amode implicit none private @@ -408,7 +408,7 @@ end function number_transported ! species morphology !------------------------------------------------------------------------ subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) use cam_abortutils, only: endrun @@ -417,6 +417,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 @@ -429,6 +430,10 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & density_aer=density, hygro_aer=hygro, spectype=spectype, & refindex_aer_sw=refindex_sw, refindex_aer_lw=refindex_lw) + if (present(spec_mw)) then + spec_mw = specmw_amode(species_ndx,bin_ndx) + end if + if (present(specname)) then call rad_aer_get_info(self%list_idx_, bin_ndx, species_ndx, spec_name=specname) end if diff --git a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 index 7e2b57d130..e1601bd82f 100644 --- a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 @@ -10,9 +10,10 @@ module modal_aerosol_state_mod use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_get_index use physics_types, only: physics_state !REMOVECAM_END - use aerosol_properties_mod, only: aerosol_properties + use aerosol_properties_mod, only: aerosol_properties, aero_name_len use physconst, only: rhoh2o use cam_abortutils, only: endrun + use ppgrid, only: pver implicit none @@ -49,6 +50,7 @@ module modal_aerosol_state_mod procedure :: wet_diameter procedure :: convcld_actfrac procedure :: wgtpct + procedure :: aqu_gain_binfraction final :: destructor @@ -215,7 +217,6 @@ end subroutine get_states ! return aerosol bin size weights for a given bin !------------------------------------------------------------------------------ subroutine icenuc_size_wght_arr(self, bin_ndx, ncol, nlev, species_type, use_preexisting_ice, wght) - use aerosol_properties_mod, only: aero_name_len class(modal_aerosol_state), intent(in) :: self integer, intent(in) :: bin_ndx ! bin number @@ -278,7 +279,6 @@ end subroutine icenuc_size_wght_arr ! return aerosol bin size weights for a given bin, column and vertical layer !------------------------------------------------------------------------------ subroutine icenuc_size_wght_val(self, bin_ndx, col_ndx, lyr_ndx, species_type, use_preexisting_ice, wght) - use aerosol_properties_mod, only: aero_name_len class(modal_aerosol_state), intent(in) :: self integer, intent(in) :: bin_ndx ! bin number @@ -339,9 +339,6 @@ end subroutine icenuc_size_wght_val !------------------------------------------------------------------------------ subroutine icenuc_type_wght(self, bin_ndx, ncol, nlev, species_type, aero_props, rho, wght, cloud_borne) - use aerosol_properties_mod, only: aerosol_properties - use aerosol_properties_mod, only: aero_name_len - class(modal_aerosol_state), intent(in) :: self integer, intent(in) :: bin_ndx ! bin number integer, intent(in) :: ncol ! number of columns @@ -418,7 +415,6 @@ end subroutine update_bin ! as heterogeneous freezing nuclei !------------------------------------------------------------------------------ function hetfrz_size_wght(self, bin_ndx, ncol, nlev) result(wght) - use aerosol_properties_mod, only: aero_name_len class(modal_aerosol_state), intent(in) :: self integer, intent(in) :: bin_ndx ! bin number @@ -634,7 +630,6 @@ end function wet_diameter ! prescribed aerosol activation fraction for convective cloud !------------------------------------------------------------------------------ function convcld_actfrac(self, aero_props, ibin, ispc, ncol, nlev) result(frac) - use aerosol_properties_mod, only: aero_name_len class(modal_aerosol_state), intent(in) :: self class(aerosol_properties), intent(in) :: aero_props ! aerosol properties object @@ -723,4 +718,88 @@ 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(modal_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 + + character(len=aero_name_len) :: modetype, spectype + integer :: i,k,l,m,n,mm, ncol, nbins + integer :: accum_n + real(r8) :: sumf + real(r8), allocatable :: qnum_c(:) + + ncol = self%state%ncol + nbins = aero_props%nbins() + + !------------------------------------------------------------------------- + ! compute factors for partitioning aerosol mass gains among modes. + ! The factors are proportional to the activated particle MR for each + ! mode, which is the MR of cloud drops "associated with" the mode + ! thus we are assuming the cloud drop size is independent of the + ! associated aerosol mode properties (i.e., drops associated with + ! Aitken and coarse sea-salt particles are same size) + ! + ! qnum_c(n) = activated particle number MR for mode n (these are just + ! used for partitioning among modes, so don't need to divide by cldfrc) + !------------------------------------------------------------------------- + + accum_n = -1 + do m = 1, nbins + call rad_aer_get_info(0, m, mode_type=modetype) + if (modetype=='accum') then + accum_n = m + end if + end do + + allocate(qnum_c(nbins)) + + faqgain = 0.0_r8 + + lev_loop: do k = 1,pver + col_loop: do i = 1,ncol + do m = 1, nbins + mm = aero_props%indexer(m,0) + qnum_c(m) = max( 0.0_r8, qcw(i,k,mm) ) + end do + + ! force qnum_c(n) to be positive for n=modeptr_accum or n=1 + n = accum_n + if (n <= 0) n = 1 + qnum_c(n) = max( 1.0e-10_r8, qnum_c(n) ) + + ! faqgain_so4(n) = fraction of total so4_c gain going to mode n + ! these are proportional to the activated particle MR for each mode + sumf = 0.0_r8 + do n = 1, nbins + do l = 1, aero_props%nspecies(n) + call aero_props%get(n,l, spectype=spectype) + if (trim(spectype) == trim(type)) then + faqgain(n,i,k) = qnum_c(n) + sumf = sumf + faqgain(n,i,k) + end if + end do + end do + + if (sumf > 0.0_r8) then + do n = 1, nbins + faqgain(n,i,k) = faqgain(n,i,k) / sumf + end do + end if + ! at this point (sumf <= 0.0) only when all the faqgain_so4 are zero + + end do col_loop + end do lev_loop + + deallocate(qnum_c) + + end subroutine aqu_gain_binfraction + end module modal_aerosol_state_mod diff --git a/src/chemistry/carma_aero/aero_model.F90 b/src/chemistry/carma_aero/aero_model.F90 index 778e8093f6..c53669355e 100644 --- a/src/chemistry/carma_aero/aero_model.F90 +++ b/src/chemistry/carma_aero/aero_model.F90 @@ -65,22 +65,7 @@ module aero_model integer, public, protected :: nspec_max = 0 integer, public, protected :: nbins = 0 integer, public, protected, allocatable :: nspec(:) - - ! local indexing for bins - integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr - integer :: ncnst_tot ! total number of mode number conc + mode species - integer :: ncnst_extd ! twiece total number of mode number conc + mode species - - ! Indices for CARMA species in the ptend%q array. Needed for prognostic aerosol case. - logical, allocatable :: bin_cnst_lq(:,:) - integer, allocatable :: bin_cnst_idx(:,:) - - - ! ptr2d_t is used to create arrays of pointers to 2D fields - type ptr2d_t - real(r8), pointer :: fld(:,:) => null() - end type ptr2d_t - + integer :: ncnst_tot logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies ! in the ptend object @@ -209,7 +194,7 @@ subroutine aero_model_init( pbuf2d ) ! local vars character(len=*), parameter :: subrname = 'aero_model_init' - integer :: m, n, ii, mm + integer :: m, n, mm integer :: idxtmp = -1 logical :: history_aerosol ! Output MAM or SECT aerosol tendencies @@ -273,50 +258,33 @@ subroutine aero_model_init( pbuf2d ) nspec_max = maxval(nspec) - ncnst_tot = nspec(1) - do m = 2, nbins - ncnst_tot = ncnst_tot + nspec(m) - end do - ncnst_extd = 2*ncnst_tot + ncnst_tot = aero_props%ncnst_tot() allocate( & - bin_idx(nbins,nspec_max), & - bin_cnst_lq(nbins,nspec_max), & - bin_cnst_idx(nbins,nspec_max), & fieldname_cw(ncnst_tot), & fieldname(ncnst_tot), stat=ierr ) if (ierr/=0) call endrun(subrname//' : allocate error') - ii = 0 do m = 1, nbins do l = 1, nspec(m) ! loop through species - ii = ii + 1 - bin_idx(m,l) = ii + mm = aero_props%indexer(m,l) if (l <= nspec(m) ) then ! species - call rad_aer_get_info_by_bin_spec(0, m, l, spec_name=fieldname(ii), spec_name_cw=fieldname_cw(ii)) + call rad_aer_get_info_by_bin_spec(0, m, l, spec_name=fieldname(mm), spec_name_cw=fieldname_cw(mm)) else !number - call rad_aer_get_info_by_bin(0, m, num_name=fieldname(ii), num_name_cw=fieldname_cw(ii)) + call rad_aer_get_info_by_bin(0, m, num_name=fieldname(mm), num_name_cw=fieldname_cw(mm)) end if - call cnst_get_ind(fieldname(ii), idxtmp, abort=.false.) + call cnst_get_ind(fieldname(mm), idxtmp, abort=.false.) if (idxtmp.gt.0) then - bin_cnst_lq(m,l) = .true. - bin_cnst_idx(m,l) = idxtmp lq(idxtmp) = .true. call cnst_set_convtran2(idxtmp, .not.convproc_do_aer) - else - bin_cnst_lq(m,l) = .false. - bin_cnst_idx(m,l) = 0 end if - mm = ii - - unit_basename = 'kg' - if (l == nspec(m) + 2) then ! number - unit_basename = ' 1' - end if - + unit_basename = 'kg' + if (l == nspec(m) + 2) then ! number + unit_basename = ' 1' + end if call addfld( fieldname_cw(mm), (/ 'lev' /), 'A', unit_basename//'/kg ', & trim(fieldname_cw(mm))//' in cloud water') @@ -355,7 +323,7 @@ subroutine aero_model_init( pbuf2d ) if (has_sox) then do n = 1, nbins do l = 1, nspec(n) ! not for total mass or number - mm = bin_idx(n, l) + mm = aero_props%indexer(n,l) call addfld (& trim(fieldname_cw(mm))//'AQSO4',horiz_only, 'A','kg/m2/s', & trim(fieldname_cw(mm))//' aqueous phase chemistry') @@ -501,10 +469,13 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re use carma_aero_gasaerexch, only : carma_aero_gasaerexch_sub use time_manager, only : get_nstep + use carma_aerosol_state_mod, only : carma_aerosol_state + use aerosol_state_mod, only: aerosol_state, ptr2d_t + !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- - type(physics_state), intent(in) :: state ! Physics state variables + type(physics_state),target, intent(in) :: state ! Physics state variables integer, intent(in) :: loffset ! offset applied to modal aero "pointers" integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: lchnk ! chunk index @@ -540,8 +511,6 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8) :: del_h2so4_aeruptk(ncol,pver) - real(r8), pointer :: pblh(:) ! pbl height (m) - real(r8), dimension(ncol) :: wrk character(len=32) :: name real(r8) :: dvmrcwdt(ncol,pver,ncnst_tot) @@ -559,14 +528,10 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc real(r8) :: nh3_beg(ncol,pver) real(r8) :: mw_carma(ncnst_tot) - real(r8), pointer :: fldcw(:,:) - real(r8), pointer :: sulfeq(:,:,:) real(r8) :: wetr(pcols,pver) ! CARMA wet radius in cm real(r8) :: wetrho(pcols,pver) ! CARMA wet dens real(r8), allocatable :: rmass(:) ! CARMA rmass - real(r8) :: old_total_mass - real(r8) :: new_total_mass real(r8) :: old_total_number character(len=32) :: spectype @@ -575,6 +540,11 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re integer :: igroup, ibin, rc, nchr, ierr character(len=*), parameter :: subname = 'aero_model_gasaerexch' + class(aerosol_state), pointer :: aero_state + +!---------------------------------------------------------------------- + aero_state => carma_aerosol_state(state, pbuf) + ! ! ... initialize nh3 ! @@ -606,6 +576,10 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re qqcw(ncnst_tot), stat=ierr ) if (ierr /= 0) call endrun(subname//': allocate error') + ! Init pointers to mode number and specie mass mixing ratios in + ! intersitial and cloud borne phases. + call aero_state%get_states( aero_props, raer, qqcw ) + mw_carma(:) = 0.0_r8 do m = 1, nbins ! main loop over aerosol bins ! dryr is the dry bin radius @@ -641,11 +615,9 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! Init pointers to mode number and specie mass mixing ratios in ! intersitial and cloud borne phases. do l = 1, nspec(m) - mm = bin_idx(m, l) + mm = aero_props%indexer(m,l) if (l <= nspec(m)) then - call rad_aer_get_bin_props_by_idx(0, m, l,spectype=spectype) - call rad_cnst_get_bin_mmr_by_idx(0, m, l, 'a', state, pbuf, raer(mm)%fld) - call rad_cnst_get_bin_mmr_by_idx(0, m, l, 'c', state, pbuf, qqcw(mm)%fld) ! cloud-borne aerosol + call aero_props%get(bin_ndx=m, species_ndx=l, spectype=spectype) if (trim(spectype) == 'sulfate') then mw_carma(mm) = 96._r8 end if @@ -680,7 +652,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! aqueous chemistry ... if( has_sox ) then - call setsox( state, & + call setsox( aero_state, state, & pbuf, & ncol, & lchnk, & @@ -705,10 +677,13 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re do n = 1, nbins do l = 1, nspec(n) ! not for total mass or number - mm = bin_idx(n, l) - call outfld( trim(fieldname_cw(mm))//'AQSO4', aqso4(:ncol,mm), ncol, lchnk) - call outfld( trim(fieldname_cw(mm))//'AQH2SO4', aqh2so4(:ncol,mm), ncol, lchnk) - end do + call aero_props%get(bin_ndx=n, species_ndx=l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + mm = aero_props%indexer(n,l) + call outfld( trim(fieldname_cw(mm))//'AQSO4', aqso4(:ncol,n), ncol, lchnk) + call outfld( trim(fieldname_cw(mm))//'AQH2SO4', aqh2so4(:ncol,n), ncol, lchnk) + end if + end do end do call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk) @@ -752,40 +727,16 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! note vmr2qqcw does not change qqcw pointer (different than in MAM) call vmr2mmr_carma ( lchnk, vmrcw, mbar, mw_carma, ncol, loffset, rmass ) - !vmrcw in kg/kg - ! change pointer value for total mmr and number. In order to do this correctly - ! only mass has to be added to each bin (not number). This will require redistributing - ! mass to different bins. Here, we change both mass and number until we have a better - ! solution. - delta_so4mass(:,:,:) = 0.0_r8 - do m = 1, nbins - do l = 1, nspec(m) ! for sulfate only - mm = bin_idx(m, l) - ! sulfate mass that needs to be added to the total mass - call rad_aer_get_bin_props_by_idx(0, m, l,spectype=spectype) - if (trim(spectype) == 'sulfate') then - ! only do loop if vmrcw has changed - do k=1,pver - do i=1,ncol - if (vmrcw(i,k,mm) .gt. mmrcw(i,k,mm) .and. mmrcw(i,k,mm) /= 0.0_r8) then - delta_so4mass(i,k,mm) = ( vmrcw(i,k,mm) - mmrcw(i,k,mm) ) - else - delta_so4mass(i,k,mm) = 0.0_r8 - end if - end do - end do - end if - end do - end do - do m = 1, nbins do l = 1, nspec(m) ! for sulfate only - mm = bin_idx(m, l) + mm = aero_props%indexer(m,l) qqcw(mm)%fld(:ncol,:) = vmrcw(:ncol,:,mm) call outfld( trim(fieldname_cw(mm)), qqcw(mm)%fld(:ncol,:), ncol, lchnk) end do end do + deallocate(aero_state) + nullify(aero_state) end subroutine aero_model_gasaerexch @@ -941,7 +892,7 @@ subroutine mmr2vmr_carma(lchnk, vmr, mbar, mw_carma, ncol, im, rmass) do m = 1, nbins do l = 1, nspec(m) ! for each species, not total mmr or number, information of mw are missing - mm = bin_idx(m, l) + mm = aero_props%indexer(m,l) do k=1,pver vmr(:ncol,k,mm) = mbar(:ncol,k) * vmr(:ncol,k,mm) / mw_carma(mm) end do @@ -977,7 +928,7 @@ subroutine vmr2mmr_carma ( lchnk, vmr, mbar, mw_carma, ncol, im, rmass ) !----------------------------------------------------------------- do m = 1, nbins do l = 1, nspec(m) ! for each species, not total mmr or number, information of mw are missing - mm = bin_idx(m, l) + mm = aero_props%indexer(m,l) do k=1,pver vmr(:ncol,k,mm) = mw_carma(mm) * vmr(:ncol,k,mm) / mbar(:ncol,k) end do diff --git a/src/chemistry/carma_aero/carma_aero_gasaerexch.F90 b/src/chemistry/carma_aero/carma_aero_gasaerexch.F90 index c44fe03704..342692b1f7 100644 --- a/src/chemistry/carma_aero/carma_aero_gasaerexch.F90 +++ b/src/chemistry/carma_aero/carma_aero_gasaerexch.F90 @@ -15,6 +15,7 @@ module carma_aero_gasaerexch use radiative_aerosol, only: rad_aer_get_info, rad_aer_get_info_by_bin, rad_aer_get_bin_props_by_idx, & rad_aer_get_info_by_bin_spec use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field + use carma_aerosol_properties_mod, only: carma_aerosol_properties implicit none private @@ -36,8 +37,8 @@ module carma_aero_gasaerexch character(len=32), allocatable :: fldname_cw(:) ! names for cloud_borne output fields ! local indexing for bins - integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr integer :: ncnst_tot ! total number of mode number conc + mode species + type(carma_aerosol_properties), pointer :: aero_props =>null() real(r8) :: mw_soa = 250._r8 integer :: fracvbs_idx = -1 @@ -121,28 +122,14 @@ subroutine carma_aero_gasaerexch_init nspec_max = maxval(nspec) - ncnst_tot = nspec(1) - do m = 2, nbins - ncnst_tot = ncnst_tot + nspec(m) - end do + aero_props => carma_aerosol_properties() + + ncnst_tot = aero_props%ncnst_tot() - allocate( bin_idx(nbins,nspec_max), & - do_soag_any(nbins), & + allocate( do_soag_any(nbins), & fldname_cw(ncnst_tot), & fldname(ncnst_tot) ) - ! Local indexing compresses the mode and number/mass indicies into one index. - ! This indexing is used by the pointer arrays used to reference state and pbuf - ! fields. - ! for CARMA we add number = 0, total mass = 1, and mass from each constituence into mm. - ii = 0 - do m = 1, nbins - do l = 1, nspec(m) ! do through nspec - ii = ii + 1 - bin_idx(m,l) = ii - end do - end do - ! SOAG / SOA / POM information ! Define number of VBS bins (nsoa) based on number of SOAG chemistry species @@ -236,7 +223,7 @@ subroutine carma_aero_gasaerexch_init ! define history fields for basic gas-aer exchange do m = 1, nbins do l = 1, nspec(m) ! do through nspec - ii = bin_idx(m,l) + ii = aero_props%indexer(m,l) if (l <= nspec(m) ) then ! species call rad_aer_get_info_by_bin_spec(0, m, l, spec_name=fldname(ii) ) ! only write out SOA exchange here @@ -429,7 +416,7 @@ subroutine carma_aero_gasaerexch_sub( state, & n = 0 nn = 0 do l = 1, nspec(m) - mm = bin_idx(m, l) + mm = aero_props%indexer(m,l) call rad_aer_get_bin_props_by_idx(0, m, l, spectype=spectype) if (trim(spectype) == 's-organic') then n = n + 1 @@ -507,6 +494,8 @@ subroutine carma_aero_gasaerexch_sub( state, & num_bin, t, pmid, & wetr_n, uptkrate ) + uptkrate_all = 0.0_r8 + do m = 1, nbins write(fieldname,'("NUMDENS_bin",I2.2)') m @@ -690,7 +679,7 @@ subroutine carma_aero_gasaerexch_sub( state, & if (do_soag_any(m)) then j = 0 do l = 1, nspec(m) - mm = bin_idx(m,l) + mm = aero_props%indexer(m,l) call rad_aer_get_bin_props_by_idx(0, m, l,spectype=spectype) if (trim(spectype) == 's-organic') then j = j + 1 @@ -709,6 +698,9 @@ subroutine carma_aero_gasaerexch_sub( state, & end if end do ! m = ... + deallocate(aero_state) + nullify(aero_state) + end subroutine carma_aero_gasaerexch_sub diff --git a/src/chemistry/carma_aero/sox_cldaero_mod.F90 b/src/chemistry/carma_aero/sox_cldaero_mod.F90 index fead6a2fe8..3d2d0845f8 100644 --- a/src/chemistry/carma_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/carma_aero/sox_cldaero_mod.F90 @@ -8,13 +8,14 @@ module sox_cldaero_mod use ppgrid, only : pcols, pver use mo_chem_utls, only : get_spc_ndx use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate - use cam_logfile, only : iulog - use chem_mods, only : adv_mass use physconst, only : gravit - use phys_control, only : phys_getopts + use phys_control, only : cam_chempkg_is use cldaero_mod, only : cldaero_uptakerate use chem_mods, only : gas_pcnst - use radiative_aerosol, only: rad_aer_get_info, rad_aer_get_info_by_bin, rad_aer_get_bin_props_by_idx + use carma_aerosol_properties_mod, only: carma_aerosol_properties + use aerosol_state_mod, only: aerosol_state + + use modal_aero_data, only : ntot_amode implicit none private @@ -24,18 +25,16 @@ module sox_cldaero_mod public :: sox_cldaero_update public :: sox_cldaero_destroy_obj - integer :: id_msa, id_h2so4, id_so2, id_h2o2, id_nh3 + integer :: id_msa=-1, id_h2so4=-1, id_so2=-1, id_h2o2=-1, id_nh3=-1 real(r8), parameter :: small_value = 1.e-20_r8 - ! description of bin aerosols - integer, public, protected :: nspec_max = 0 + integer :: ncnst_tot = -huge(1) ! total number of mode number conc + mode species integer, public, protected :: nbins = 0 - integer, public, protected, allocatable :: nspec(:) - ! local indexing for bins - integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr - integer :: ncnst_tot ! total number of mode number conc + mode species + type(carma_aerosol_properties), pointer :: aero_props =>null() + + logical :: has_msa = .false. contains @@ -44,58 +43,22 @@ module sox_cldaero_mod subroutine sox_cldaero_init - integer :: l, m, ii - logical :: history_aerosol ! Output the MAM aerosol tendencies - id_msa = get_spc_ndx( 'MSA' ) id_h2so4 = get_spc_ndx( 'H2SO4' ) id_so2 = get_spc_ndx( 'SO2' ) id_h2o2 = get_spc_ndx( 'H2O2' ) id_nh3 = get_spc_ndx( 'NH3' ) + has_msa = id_msa>0 if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then call endrun('sox_cldaero_init:MAM mech does not include necessary species' & //' -- should not invoke sox_cldaero_mod ') endif - call phys_getopts( history_aerosol_out = history_aerosol ) - ! - ! add to history - ! - - ! get info about the modal aerosols - ! get nbins - - call rad_aer_get_info( 0, nbins=nbins) - - allocate( nspec(nbins) ) - - do m = 1, nbins - call rad_aer_get_info_by_bin(0, m, nspec=nspec(m)) - end do - ! add plus one to include number, total mmr and nspec - nspec_max = maxval(nspec) - - ncnst_tot = nspec(1) - do m = 2, nbins - ncnst_tot = ncnst_tot + nspec(m) - end do - - allocate( bin_idx(nbins,nspec_max) ) - - - ! Local indexing compresses the mode and number/mass indicies into one index. - ! This indexing is used by the pointer arrays used to reference state and pbuf - ! fields. - ! for CARMA we add number = 0, total mass = 1, and mass from each constituence into mm. - ii = 0 - do m = 1, nbins - do l = 1, nspec(m) ! loop through species - ii = ii + 1 - bin_idx(m,l) = ii - end do - end do + aero_props => carma_aerosol_properties() + ncnst_tot = aero_props%ncnst_tot() + nbins = aero_props%nbins() end subroutine sox_cldaero_init @@ -110,8 +73,6 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( integer, intent(in) :: ncol integer, intent(in) :: loffset - real(r8) :: so4mmr(pcols,pver) - type(cldaero_conc_t), pointer :: conc_obj character(len=32) :: spectype @@ -119,9 +80,9 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( integer :: l,m integer :: i,k,mm - ! local indexing for bins - !integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr + logical :: mode7 + mode7 = ntot_amode == 7 conc_obj => cldaero_allocate() @@ -140,39 +101,46 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj%nh4c(:,:) = 0._r8 conc_obj%so4c(:,:) = 0._r8 - so4mmr(:,:) = 0._r8 do k = 1,pver do i = 1,ncol - do m = 1, nbins - do l = 1, nspec(m) - mm = bin_idx(m, l) - call rad_aer_get_bin_props_by_idx(0, m, l,spectype=spectype) - if (trim(spectype) == 'sulfate') then - so4mmr(i,k) = so4mmr(i,k) + qcw(i,k,mm) - end if + do m = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + call aero_props%get(m,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + conc_obj%so4c(i,k) = conc_obj%so4c(i,k) + qcw(i,k,mm) + end if + if (trim(spectype) == 'ammonium') then + conc_obj%nh4c(i,k) = conc_obj%nh4c(i,k) + qcw(i,k,mm) + end if end do end do end do end do - conc_obj%so4c = so4mmr - end function sox_cldaero_create_obj + ! *** NOTE *** + ! should refactor this bit after merging to later cam tag -- where aero_props%model_is is avail + if (ntot_amode>0) then + if (.not.mode7) then + conc_obj%so4_fact = 1._r8 + end if + end if + end function sox_cldaero_create_obj !---------------------------------------------------------------------------------- ! Update the mixing ratios !---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( & + subroutine sox_cldaero_update( aero_state, & state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d) - use aerosol_properties_mod, only: aero_name_len use physics_types, only: physics_state - use carma_intr, only: carma_get_group_by_name, carma_get_dry_radius ! args + class(aerosol_state), intent(in) :: aero_state type(physics_state), intent(in) :: state ! Physics state variables integer, intent(in) :: ncol @@ -205,7 +173,7 @@ subroutine sox_cldaero_update( & real(r8), intent(in) :: xh2o2(:,:) real(r8), intent(in) :: xno3c(:,:) - real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) vmrcw(ncol,pver,ncnst_tot) + real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr ) real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry @@ -216,75 +184,60 @@ subroutine sox_cldaero_update( & real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) ! local vars ... - real(r8) :: dryr(pcols,pver) ! CARMA dry radius in cm - real(r8) :: rho(pcols,pver) ! - real(r8) :: dryr_n(nbins,ncol,pver) ! CARMA dry radius in cm + real(r8) :: dqdt_aqso4(ncol,pver,ncnst_tot), & dqdt_aqh2so4(ncol,pver,ncnst_tot), & dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver) - real(r8) :: faqgain_so4(nbins) - real(r8) :: wt_mass(nbins) + real(r8) :: faqgain_msa(nbins,ncol,pver), faqgain_so4(nbins,ncol,pver) + real(r8) :: delso4_3d(ncol,pver) + real(r8) :: delnh3, delnh4 real(r8) :: delso4_o3rxn, & dso4dt_aqrxn, dso4dt_hprxn, & - dso4dt_gasuptk, dmsadt_gasuptk_toso4, & + dso4dt_gasuptk, dmsadt_gasuptk, & + dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, & dqdt_aq, dqdt_wr, dqdt real(r8) :: fwetrem, uptkrate - integer :: l, n, mm - integer :: ntot_msa_c - - integer :: i,k + integer :: l, m, n, mm + integer :: i,k, ndx real(r8) :: xl - real(r8) :: wt_sum - real(r8) :: specmw_so4_amode - + real(r8) :: mw_so4 character(len=32) :: spectype - - character(len=*), parameter :: subname = 'sox_cldaero_update' - character(len=aero_name_len) :: bin_name, shortname - integer :: igroup, ibin, rc, nchr + character(len=32) :: specname ! make sure dqdt is zero initially, for budgets dqdt_aqso4(:,:,:) = 0.0_r8 dqdt_aqh2so4(:,:,:) = 0.0_r8 dqdt_aqhprxn(:,:) = 0.0_r8 dqdt_aqo3rxn(:,:) = 0.0_r8 - dryr_n(:,:,:) = 0.0_r8 - ntot_msa_c = 0.0_r8 aqso4 = 0.0_r8 aqh2so4 = 0.0_r8 aqso4_h2o2 = 0.0_r8 aqso4_o3 = 0.0_r8 - - do n = 1, nbins - call rad_aer_get_info_by_bin(0, n, nspec=nspec(n), bin_name=bin_name) - - - nchr = len_trim(bin_name)-2 - shortname = bin_name(:nchr) - - call carma_get_group_by_name(shortname, igroup, rc) - if (rc/=0) then - call endrun(subname//': ERROR in carma_get_group_by_name') - end if - - read(bin_name(nchr+1:),*) ibin - - call carma_get_dry_radius(state, igroup, ibin, dryr, rho, rc) - if (rc/=0) then - call endrun(subname//': ERROR in carma_get_dry_radius') - end if - - dryr(:ncol,:) = dryr(:ncol,:)*1.e2_r8 ! cm - - if (index(bin_name,'MXAER')>0) then - dryr_n(n,:ncol,:) = dryr(:ncol,:) - end if - end do + delso4_3d = 0.0_r8 + + ! Avoid double counting in-cloud sulfur oxidation when running with + ! GEOS-Chem. If running with GEOS-Chem then sulfur oxidation + ! is performed internally to GEOS-Chem. Here, we just return to the + ! parent routine and thus we do not apply tendencies calculated by MAM. + if ( cam_chempkg_is('geoschem_mam4') ) return + + where (cldfrc(:ncol,:) >= 1.0e-5_r8) + delso4_3d(:ncol,:) = xso4(:ncol,:) - xso4_init(:ncol,:) + end where + + !------------------------------------------------------------------------- + ! Compute factors for partitioning aerosol mass gains among bins / modes. + ! The factors are proportional to the activated particle MR for each + ! bin, which is the MR of cloud drops "associated with" the mode + ! thus we are assuming the cloud drop size is independent of the + ! associated aerosol mode properties + call aero_state%aqu_gain_binfraction(aero_props, 'sulfate', qcw, delso4_3d, faqgain_so4) + if (has_msa) call aero_state%aqu_gain_binfraction(aero_props, 'msa', qcw, delso4_3d, faqgain_msa) lev_loop: do k = 1,pver col_loop: do i = 1,ncol @@ -293,41 +246,37 @@ subroutine sox_cldaero_update( & if (xl .ge. 1.e-8_r8) then !! when cloud is present - delso4_o3rxn = xso4(i,k) - xso4_init(i,k) - - ! the factors are proportional to the activated particle MR for each - ! bin, which is the MR of cloud drops "associated with" the mode - ! thus we are assuming the cloud drop size is independent of the - ! associated aerosol mode properties (i.e., drops associated with - ! Aitken and coarse sea-salt particles are same size) - ! - ! qnum_c(n) = activated particle number MR for mode n (these are just - ! used for partitioning among modes, so don't need to divide by cldfrc) - - !faqgain_so4(n) = fraction of total so4_c gain going to mode n - wt_sum = 0._r8 - wt_mass(:) = 0._r8 - faqgain_so4(:) = 0.0_r8 - do n = 1, nbins - if (dryr_n(n,i,k) > 0._r8) then - wt_mass(n) = delso4_o3rxn / dryr_n(n,i,k) / dryr_n(n,i,k) - wt_sum = wt_sum + wt_mass(n) - end if - end do - do n = 1, nbins - if (wt_mass(n) > 0._r8) then - faqgain_so4(n) = wt_mass(n)/wt_sum - end if - end do + delso4_o3rxn = delso4_3d(i,k) ! xso4(i,k) - xso4_init(i,k) + + if (id_nh3>0) then + delnh3 = nh3g(i,k) - xnh3(i,k) + delnh4 = - delnh3 + endif + + ! faqgain_msa(n) = fraction of total msa_c gain going to mode n uptkrate = cldaero_uptakerate( xl, cldnum(i,k), cfact(i,k), cldfrc(i,k), tfld(i,k), press(i,k) ) ! average uptake rate over dtime uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime + ! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s) + ! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s) dso4dt_gasuptk = xh2so4(i,k) * uptkrate + if (has_msa) then + dmsadt_gasuptk = xmsa(i,k) * uptkrate + else + dmsadt_gasuptk = 0.0_r8 + end if ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4 - dmsadt_gasuptk_toso4 = 0.0_r8 + if (has_msa) then + dmsadt_gasuptk_toso4 = 0.0_r8 + dmsadt_gasuptk_tomsa = dmsadt_gasuptk + else + ! no MSA + dmsadt_gasuptk_tomsa = 0.0_r8 + dmsadt_gasuptk_toso4 = dmsadt_gasuptk + end if !----------------------------------------------------------------------- ! now compute TMR tendencies @@ -337,34 +286,47 @@ subroutine sox_cldaero_update( & dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime dso4dt_hprxn = delso4_hprxn(i,k) / dtime - !write(iulog,*) 'dso4dt_aqrxn ',dso4dt_aqrxn ! fwetrem = fraction of in-cloud-water material that is wet removed ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) ) fwetrem = 0.0_r8 ! don't have so4 & msa wet removal here - ! compute TMR tendencies for so4, not done currently for msa aerosol-in-cloud-water - do n = 1, nbins - do l = 1, nspec(n) - mm = bin_idx(n, l) - call rad_aer_get_bin_props_by_idx(0, n, l,spectype=spectype) - if (trim(spectype) == 'sulfate') then - if (faqgain_so4(n) .gt. 0.0_r8) then - dqdt_aqso4(i,k,mm) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k) - - dqdt_aqh2so4(i,k,mm) = faqgain_so4(n)* & - (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k) - dqdt_aq = dqdt_aqso4(i,k,mm) + dqdt_aqh2so4(i,k,mm) - dqdt_wr = -fwetrem*dqdt_aq - dqdt= dqdt_aq + dqdt_wr - !write(iulog,*) 'qcw(i,k,mm) before ', m, qcw(i,k,mm) - qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime - !write(iulog,*) 'qcw(i,k,mm) after', m, qcw(i,k,mm) - end if - end if + ! compute TMR tendencies for so4 and msa aerosol-in-cloud-water + do m = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + call aero_props%get(m,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + + dqdt_aqso4(i,k,mm) = faqgain_so4(m,i,k)*dso4dt_aqrxn*cldfrc(i,k) + + dqdt_aqh2so4(i,k,mm) = faqgain_so4(m,i,k)* & + (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k) + dqdt_aq = dqdt_aqso4(i,k,mm) + dqdt_aqh2so4(i,k,mm) + dqdt_wr = -fwetrem*dqdt_aq + dqdt = dqdt_aq + dqdt_wr + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + + end if + if (trim(spectype) == 'msa') then + dqdt_aq = faqgain_msa(m,i,k)*dmsadt_gasuptk_tomsa*cldfrc(i,k) + dqdt_wr = -fwetrem*dqdt_aq + dqdt = dqdt_aq + dqdt_wr + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + end if + if (trim(spectype) == 'ammonium') then + if (delnh4 > 0.0_r8) then + dqdt_aq = faqgain_so4(m,i,k)*delnh4/dtime*cldfrc(i,k) + dqdt = dqdt_aq + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + else + dqdt = (qcw(i,k,mm)/max(xnh4c(i,k),1.0e-35_r8)) & + *delnh4/dtime*cldfrc(i,k) + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + endif + end if end do - end do - + end do ! For gas species, tendency includes ! reactive uptake to cloud water that essentially transforms the gas to @@ -375,6 +337,7 @@ subroutine sox_cldaero_update( & ! h2so4 (g) & msa (g) qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k) + if (has_msa) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k) ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k) ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) ) @@ -384,7 +347,6 @@ subroutine sox_cldaero_update( & dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k) dqdt = dqdt_aq + dqdt_wr qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime - qin(i,k,id_so2) = MAX( qin(i,k,id_so2), small_value ) ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k) ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) ) @@ -394,13 +356,19 @@ subroutine sox_cldaero_update( & dqdt_aq = -dso4dt_hprxn*cldfrc(i,k) dqdt = dqdt_aq + dqdt_wr qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime - qin(i,k,id_h2o2) = MAX( qin(i,k,id_h2o2), small_value ) + + ! NH3 + if (id_nh3>0) then + dqdt_aq = delnh3/dtime*cldfrc(i,k) + dqdt = dqdt_aq + qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime + endif ! for SO4 from H2O2/O3 budgets dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k) dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k) - endif !! when cloud is present + endif !! when cloud is present endif cloud enddo col_loop enddo lev_loop @@ -408,66 +376,94 @@ subroutine sox_cldaero_update( & !============================================================== ! ... Update the mixing ratios !============================================================== + do k = 1,pver - ! diagnostics + do n = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(n) + mm = aero_props%indexer(n,l) + call aero_props%get(n,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + if (trim(spectype) == 'msa') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + if (trim(spectype) == 'ammonium') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + end do + end do - specmw_so4_amode = 96.0_r8 - do n = 1, nbins - ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero - do l = 1, nspec(n) - mm = bin_idx(n, l) - aqso4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s - enddo - enddo - - aqh2so4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s - enddo - enddo - end do - end do + qin(:ncol,k,id_so2) = MAX( qin(:ncol,k,id_so2), small_value ) + qin(:ncol,k,id_h2o2) = MAX( qin(:ncol,k,id_h2o2), small_value ) + qin(:ncol,k,id_h2so4) = MAX( qin(:ncol,k,id_h2so4), small_value ) + if ( id_msa > 0 ) qin(:ncol,k,id_msa) = MAX( qin(:ncol,k,id_msa), small_value ) + if ( id_nh3 > 0 ) qin(:ncol,k,id_nh3) = MAX( qin(:ncol,k,id_nh3), small_value ) + + end do + + ! diagnostics + mw_so4 = -huge(1._r8) + + do n = 1, aero_props%nbins() + ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero + do l = 1, aero_props%nspecies(n) + mm = aero_props%indexer(n,l) + call aero_props%get(n,l, spectype=spectype, specname=specname) + if (trim(spectype) == 'sulfate') then + call aero_props%get(n,l, spec_mw=mw_so4) + aqso4(:,n)=0._r8 + do k=1,pver + do i=1,ncol + aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg/m2/s + enddo + enddo + aqh2so4(:,n)=0._r8 + do k=1,pver + do i=1,ncol + aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg/m2/s + enddo + enddo + end if + end do + end do aqso4_h2o2(:) = 0._r8 do k=1,pver do i=1,ncol - aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s + aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s enddo enddo if (present(aqso4_h2o2_3d)) then - aqso4_h2o2_3d(:,:) = 0._r8 - do k=1,pver - do i=1,ncol - aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo + aqso4_h2o2_3d(:,:) = 0._r8 + do k=1,pver + do i=1,ncol + aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s + enddo + enddo end if aqso4_o3(:)=0._r8 do k=1,pver - do i=1,ncol - aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo + do i=1,ncol + aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s + enddo enddo if (present(aqso4_o3_3d)) then - aqso4_o3_3d(:,:)=0._r8 - do k=1,pver - do i=1,ncol - aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo + aqso4_o3_3d(:,:)=0._r8 + do k=1,pver + do i=1,ncol + aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s + enddo + enddo end if end subroutine sox_cldaero_update diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index 6e450ab2b4..be1da86975 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -3,7 +3,7 @@ !=============================================================================== module aero_model use shr_kind_mod, only: r8 => shr_kind_r8 - use constituents, only: pcnst, cnst_name, cnst_get_ind + use constituents, only: pcnst, cnst_name, cnst_get_ind, cnst_mw use ppgrid, only: pcols, pver, pverp use phys_control, only: phys_getopts, cam_physpkg_is use cam_abortutils, only: endrun @@ -30,6 +30,8 @@ module aero_model use modal_aero_wateruptake, only: modal_strat_sulfate use mo_setsox, only: setsox, has_sox use modal_aerosol_properties_mod, only: modal_aerosol_properties + use modal_aerosol_state_mod, only: modal_aerosol_state + use aerosol_state_mod, only: aerosol_state, ptr2d_t implicit none private @@ -100,6 +102,7 @@ module aero_model logical :: modal_accum_coarse_exch = .false. type(modal_aerosol_properties), pointer :: aero_props=>null() + integer :: ncnst_tot integer :: n_coarse_dust=-1 ! dmleung added n_coarse_dust to determine the index for the ! coarse dust mode for different MAM versions. 29 Oct 2025 @@ -203,7 +206,7 @@ subroutine aero_model_init( pbuf2d ) logical :: history_aerosol ! Output MAM or SECT aerosol tendencies logical :: history_chemistry, history_cesm_forcing, history_dust - integer :: l + integer :: l, mm character(len=6) :: test_name character(len=64) :: errmes @@ -215,6 +218,12 @@ subroutine aero_model_init( pbuf2d ) character(len=32) :: spec_type character(len=32) :: mode_type integer :: nspec + character(len=32) :: spectype + character(len=32) :: name_a + character(len=32) :: name_c + + aero_props => modal_aerosol_properties() + ncnst_tot = aero_props%ncnst_tot() ! aqueous chem initialization call sox_inti() @@ -250,7 +259,6 @@ subroutine aero_model_init( pbuf2d ) ! call aero_deposition_cam_init only if the user has not specified ! prescribed aerosol deposition fluxes if (.not.aerodep_flx_prescribed()) then - aero_props => modal_aerosol_properties() call aero_deposition_cam_init(aero_props) endif @@ -989,7 +997,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- - type(physics_state), intent(in) :: state ! Physics state variables + type(physics_state), target, intent(in) :: state ! Physics state variables integer, intent(in) :: loffset ! offset applied to modal aero "pointers" integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: lchnk ! chunk index @@ -1027,9 +1035,9 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8), dimension(ncol) :: wrk character(len=32) :: name - real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst) + real(r8) :: dvmrcwdt(ncol,pver,ncnst_tot) real(r8) :: dvmrdt(ncol,pver,gas_pcnst) - real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr) + real(r8) :: vmrcw(ncol,pver,ncnst_tot) ! cloud-borne aerosol (vmr) real(r8) :: aqso4(ncol,ntot_amode) ! aqueous phase chemistry real(r8) :: aqh2so4(ncol,ntot_amode) ! aqueous phase chemistry @@ -1040,7 +1048,21 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8), pointer :: fldcw(:,:) real(r8), pointer :: sulfeq(:,:,:) -! + character(len=32) :: spectype + character(len=32) :: specname + character(len=32) :: name_a, name_c + real(r8) :: mw(ncnst_tot) + integer :: ndx, ierr, mm + type(ptr2d_t), allocatable :: raer(:) ! aerosol mass, number mixing ratios + type(ptr2d_t), allocatable :: qqcw(:) + class(aerosol_state), pointer :: aero_state + + character(len=*), parameter :: subname = 'aero_model_gasaerexch' + +!---------------------------------------------------------------------- + aero_state => modal_aerosol_state(state, pbuf) + +!! ! ... initialize nh3 ! if ( nh3_ndx > 0 ) then @@ -1076,7 +1098,33 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! ! Aerosol processes ... ! - call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) + allocate( & + raer(ncnst_tot), & + qqcw(ncnst_tot), stat=ierr ) + if (ierr /= 0) call endrun(subname//': allocate error') + + ! Init pointers to mode number and specie mass mixing ratios in + ! intersitial and cloud borne phases. + call aero_state%get_states( aero_props, raer, qqcw ) + + mw(:) = 0.0_r8 + do m = 1, aero_props%nbins() ! main loop over aerosol bins + do l = 0, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + + call cnst_get_ind(name_a,ndx) + mw(mm) = cnst_mw(ndx) + vmrcw(:ncol,:,mm) = qqcw(mm)%fld(:ncol,:) + + end do + end do + + call qqcw2vmr( vmrcw, mw, mbar, ncol ) dvmrdt(:ncol,:,:) = vmr(:ncol,:,:) dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:) @@ -1084,7 +1132,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! aqueous chemistry ... if( has_sox ) then - call setsox( state, & + call setsox( aero_state, state, & pbuf, & ncol, & lchnk, & @@ -1192,7 +1240,14 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re call t_stopf('modal_coag') - call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) + call vmr2qqcw( vmrcw, mw, mbar, ncol ) + + do m = 1, aero_props%nbins() ! main loop over aerosol bins + do l = 0, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + qqcw(mm)%fld(:ncol,:) = vmrcw(:ncol,:,mm) + end do + end do ! diagnostics for cloud-borne aerosols... do n = 1,pcnst @@ -1999,83 +2054,63 @@ end subroutine calc_1_impact_rate !============================================================================= !============================================================================= - subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf) - use modal_aero_data, only : qqcw_get_field - use physics_buffer, only : physics_buffer_desc + subroutine qqcw2vmr(vmr, mw, mbar, ncol) !----------------------------------------------------------------- ! ... Xfrom from mass to volume mixing ratio !----------------------------------------------------------------- - use chem_mods, only : adv_mass, gas_pcnst - - implicit none - !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- - integer, intent(in) :: lchnk, ncol, im + integer, intent(in) :: ncol real(r8), intent(in) :: mbar(ncol,pver) - real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst) - type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), intent(in) :: mw(ncnst_tot) + real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- - integer :: k, m - real(r8), pointer :: fldcw(:,:) + integer :: k,l, m,mm - do m=1,gas_pcnst - if( adv_mass(m) /= 0._r8 ) then - fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) - if(associated(fldcw)) then - do k=1,pver - vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m) - end do - else - vmr(:,:,m) = 0.0_r8 - end if - end if + do m = 1, aero_props%nbins() + do l = 0, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + do k=1,pver + vmr(:ncol,k,mm) = mbar(:ncol,k) * vmr(:ncol,k,mm) / mw(mm) + end do + end do end do + end subroutine qqcw2vmr !============================================================================= !============================================================================= - subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf ) + subroutine vmr2qqcw(vmr, mw, mbar, ncol) !----------------------------------------------------------------- ! ... Xfrom from volume to mass mixing ratio !----------------------------------------------------------------- - use m_spc_id - use chem_mods, only : adv_mass, gas_pcnst - use modal_aero_data, only : qqcw_get_field - use physics_buffer, only : physics_buffer_desc - - implicit none - !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- - integer, intent(in) :: lchnk, ncol, im + integer, intent(in) :: ncol real(r8), intent(in) :: mbar(ncol,pver) - real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst) - type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), intent(in) :: mw(ncnst_tot) + real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- - integer :: k, m - real(r8), pointer :: fldcw(:,:) - !----------------------------------------------------------------- - ! ... The non-group species - !----------------------------------------------------------------- - do m = 1,gas_pcnst - fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) - if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then - do k = 1,pver - fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k) + integer :: k, l, m, mm + + do m = 1, aero_props%nbins() + do l = 0, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + do k=1,pver + vmr(:ncol,k,mm) = mw(mm) * vmr(:ncol,k,mm) / mbar(:ncol,k) end do - end if + end do end do end subroutine vmr2qqcw diff --git a/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 b/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 index 8b36461075..8b5649b79b 100644 --- a/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 +++ b/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 @@ -18,6 +18,7 @@ module modal_aero_gasaerexch use ppgrid, only: pcols, pver use modal_aero_data, only: ntot_amode, numptr_amode, sigmag_amode use modal_aero_data, only: lptr2_soa_g_amode, lptr2_soa_a_amode, lptr2_pom_a_amode + use modal_aerosol_properties_mod, only: modal_aerosol_properties implicit none private @@ -58,6 +59,10 @@ module modal_aero_gasaerexch real (r8), allocatable :: fac_m2v_pcarbon(:) + ! local indexing for bins + integer :: ncnst_tot ! total number of mode number conc + mode species + type(modal_aerosol_properties), pointer :: aero_props =>null() + ! !DESCRIPTION: This module implements ... ! ! !REVISION HISTORY: @@ -88,8 +93,8 @@ subroutine modal_aero_gasaerexch_sub( & loffset, deltat, & t, pmid, pdel, & qh2o, troplev, & - q, qqcw, & - dqdt_other, dqqcwdt_other, & + q, qqcw_in, & + dqdt_other, dqqcwdt_other_in, & dgncur_a, dgncur_awet, & sulfeq ) @@ -109,6 +114,7 @@ subroutine modal_aero_gasaerexch_sub( & use cam_abortutils, only: endrun use spmd_utils, only: iam, masterproc use phys_control, only: cam_chempkg_is +use mo_chem_utls, only: get_spc_ndx implicit none @@ -124,13 +130,13 @@ subroutine modal_aero_gasaerexch_sub( & ! *** MUST BE #/kmol-air for number ! *** MUST BE mol/mol-air for mass ! *** NOTE ncol dimension - real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx) + real(r8), intent(inout) :: qqcw_in(ncol,pver,ncnst_tot) ! like q but for cloud-borner tracers real(r8), intent(in) :: dqdt_other(ncol,pver,pcnstxx) ! TMR tendency from other continuous ! growth processes (aqchem, soa??) ! *** NOTE ncol dimension - real(r8), intent(in) :: dqqcwdt_other(ncol,pver,pcnstxx) + real(r8), intent(in) :: dqqcwdt_other_in(ncol,pver,ncnst_tot) ! like dqdt_other but for cloud-borner tracers real(r8), intent(in) :: t(pcols,pver) ! temperature at model levels (K) real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa) @@ -257,6 +263,25 @@ subroutine modal_aero_gasaerexch_sub( & real(r8) :: dqdtsv1(ncol,pver,pcnstxx) real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx) + real(r8) :: qqcw(ncol,pver,pcnstxx) + real(r8) :: dqqcwdt_other(ncol,pver,pcnstxx) + + integer :: m, mm, ndx + character(len=32) :: name_a, name_c + + do m = 1,aero_props%nbins() + do l = 0,aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + ndx = get_spc_ndx( name_a ) + qqcw(:,:,ndx) = qqcw_in(:,:,mm) + dqqcwdt_other(:,:,ndx) = dqqcwdt_other_in(:,:,mm) + end do + end do !---------------------------------------------------------------------- @@ -944,6 +969,19 @@ subroutine modal_aero_gasaerexch_sub( & end do ! jsrf = ... end do ! l = ... + do m = 1,aero_props%nbins() + do l = 0,aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + ndx = get_spc_ndx( name_a ) + qqcw_in(:,:,mm) = qqcw(:,:,ndx) + end do + end do + return end subroutine modal_aero_gasaerexch_sub @@ -1483,6 +1521,10 @@ subroutine modal_aero_gasaerexch_init logical :: history_aerocom ! Output the aerocom history !----------------------------------------------------------------------- + aero_props => modal_aerosol_properties() + + ncnst_tot = aero_props%ncnst_tot() + call phys_getopts( history_aerosol_out = history_aerosol ) maxspec_pcage = nspec_max diff --git a/src/chemistry/modal_aero/sox_cldaero_mod.F90 b/src/chemistry/modal_aero/sox_cldaero_mod.F90 index ccd5d50381..1701d4a6dd 100644 --- a/src/chemistry/modal_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/modal_aero/sox_cldaero_mod.F90 @@ -8,14 +8,14 @@ module sox_cldaero_mod use ppgrid, only : pcols, pver use mo_chem_utls, only : get_spc_ndx use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate - use modal_aero_data, only : ntot_amode, modeptr_accum, lptr_so4_cw_amode, lptr_msa_cw_amode - use modal_aero_data, only : numptrcw_amode, lptr_nh4_cw_amode - use modal_aero_data, only : cnst_name_cw, specmw_so4_amode - use chem_mods, only : adv_mass use physconst, only : gravit - use phys_control, only : phys_getopts, cam_chempkg_is + use phys_control, only : cam_chempkg_is use cldaero_mod, only : cldaero_uptakerate use chem_mods, only : gas_pcnst + use modal_aerosol_properties_mod, only: modal_aerosol_properties + use aerosol_state_mod, only: aerosol_state + + use modal_aero_data, only : ntot_amode implicit none private @@ -25,10 +25,17 @@ module sox_cldaero_mod public :: sox_cldaero_update public :: sox_cldaero_destroy_obj - integer :: id_msa, id_h2so4, id_so2, id_h2o2, id_nh3 + integer :: id_msa=-1, id_h2so4=-1, id_so2=-1, id_h2o2=-1, id_nh3=-1 real(r8), parameter :: small_value = 1.e-20_r8 + integer :: ncnst_tot = -huge(1) ! total number of mode number conc + mode species + integer, public, protected :: nbins = 0 + + type(modal_aerosol_properties), pointer :: aero_props =>null() + + logical :: has_msa = .false. + contains !---------------------------------------------------------------------------------- @@ -36,24 +43,22 @@ module sox_cldaero_mod subroutine sox_cldaero_init - integer :: l, m - logical :: history_aerosol ! Output the MAM aerosol tendencies - id_msa = get_spc_ndx( 'MSA' ) id_h2so4 = get_spc_ndx( 'H2SO4' ) id_so2 = get_spc_ndx( 'SO2' ) id_h2o2 = get_spc_ndx( 'H2O2' ) id_nh3 = get_spc_ndx( 'NH3' ) + has_msa = id_msa>0 if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then call endrun('sox_cldaero_init:MAM mech does not include necessary species' & //' -- should not invoke sox_cldaero_mod ') endif - call phys_getopts( history_aerosol_out = history_aerosol ) - ! - ! add to history - ! + aero_props => modal_aerosol_properties() + + ncnst_tot = aero_props%ncnst_tot() + nbins = aero_props%nbins() end subroutine sox_cldaero_init @@ -70,11 +75,10 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( type(cldaero_conc_t), pointer :: conc_obj + character(len=32) :: spectype - integer :: id_so4_1a, id_so4_2a, id_so4_3a, id_so4_4a, id_so4_5a, id_so4_6a - integer :: id_nh4_1a, id_nh4_2a, id_nh4_3a, id_nh4_4a, id_nh4_5a, id_nh4_6a - integer :: l,n - integer :: i,k + integer :: l,m + integer :: i,k,mm logical :: mode7 @@ -94,64 +98,40 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( enddo conc_obj%no3c(:,:) = 0._r8 + conc_obj%nh4c(:,:) = 0._r8 + conc_obj%so4c(:,:) = 0._r8 - if (mode7) then -#if ( defined MODAL_AERO_7MODE ) -!put ifdef here so ifort will compile - id_so4_1a = lptr_so4_cw_amode(1) - loffset - id_so4_2a = lptr_so4_cw_amode(2) - loffset - id_so4_3a = lptr_so4_cw_amode(4) - loffset - id_so4_4a = lptr_so4_cw_amode(5) - loffset - id_so4_5a = lptr_so4_cw_amode(6) - loffset - id_so4_6a = lptr_so4_cw_amode(7) - loffset - - id_nh4_1a = lptr_nh4_cw_amode(1) - loffset - id_nh4_2a = lptr_nh4_cw_amode(2) - loffset - id_nh4_3a = lptr_nh4_cw_amode(4) - loffset - id_nh4_4a = lptr_nh4_cw_amode(5) - loffset - id_nh4_5a = lptr_nh4_cw_amode(6) - loffset - id_nh4_6a = lptr_nh4_cw_amode(7) - loffset -#endif - conc_obj%so4c(:ncol,:) & - = qcw(:ncol,:,id_so4_1a) & - + qcw(:ncol,:,id_so4_2a) & - + qcw(:ncol,:,id_so4_3a) & - + qcw(:ncol,:,id_so4_4a) & - + qcw(:ncol,:,id_so4_5a) & - + qcw(:ncol,:,id_so4_6a) - - conc_obj%nh4c(:ncol,:) & - = qcw(:ncol,:,id_nh4_1a) & - + qcw(:ncol,:,id_nh4_2a) & - + qcw(:ncol,:,id_nh4_3a) & - + qcw(:ncol,:,id_nh4_4a) & - + qcw(:ncol,:,id_nh4_5a) & - + qcw(:ncol,:,id_nh4_6a) - else - id_so4_1a = lptr_so4_cw_amode(1) - loffset - id_so4_2a = lptr_so4_cw_amode(2) - loffset - id_so4_3a = lptr_so4_cw_amode(3) - loffset - conc_obj%so4c(:ncol,:) & - = qcw(:,:,id_so4_1a) & - + qcw(:,:,id_so4_2a) & - + qcw(:,:,id_so4_3a) - - ! for 3-mode, so4 is assumed to be nh4hso4 - ! the partial neutralization of so4 is handled by using a - ! -1 charge (instead of -2) in the electro-neutrality equation - conc_obj%nh4c(:ncol,:) = 0._r8 - - ! with 3-mode, assume so4 is nh4hso4, and so half-neutralized - conc_obj%so4_fact = 1._r8 + do k = 1,pver + do i = 1,ncol + do m = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + call aero_props%get(m,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + conc_obj%so4c(i,k) = conc_obj%so4c(i,k) + qcw(i,k,mm) + end if + if (trim(spectype) == 'ammonium') then + conc_obj%nh4c(i,k) = conc_obj%nh4c(i,k) + qcw(i,k,mm) + end if + end do + end do + end do + end do - endif + ! *** NOTE *** + ! should refactor this bit after merging to later cam tag -- where aero_props%model_is is avail + if (ntot_amode>0) then + if (.not.mode7) then + conc_obj%so4_fact = 1._r8 + end if + end if end function sox_cldaero_create_obj !---------------------------------------------------------------------------------- ! Update the mixing ratios !---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( & + subroutine sox_cldaero_update( aero_state, & state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d) @@ -160,6 +140,7 @@ subroutine sox_cldaero_update( & ! args + class(aerosol_state), intent(in) :: aero_state type(physics_state), intent(in) :: state ! Physics state variables integer, intent(in) :: ncol @@ -202,30 +183,30 @@ subroutine sox_cldaero_update( & real(r8), intent(out), optional :: aqso4_h2o2_3d(:,:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - ! local vars ... - real(r8) :: dqdt_aqso4(ncol,pver,gas_pcnst), & - dqdt_aqh2so4(ncol,pver,gas_pcnst), & - dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), & - sflx(1:ncol) + real(r8) :: dqdt_aqso4(ncol,pver,ncnst_tot), & + dqdt_aqh2so4(ncol,pver,ncnst_tot), & + dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver) - real(r8) :: faqgain_msa(ntot_amode), faqgain_so4(ntot_amode), qnum_c(ntot_amode) + real(r8) :: faqgain_msa(nbins,ncol,pver), faqgain_so4(nbins,ncol,pver) + real(r8) :: delso4_3d(ncol,pver) + real(r8) :: delnh3, delnh4 real(r8) :: delso4_o3rxn, & dso4dt_aqrxn, dso4dt_hprxn, & dso4dt_gasuptk, dmsadt_gasuptk, & dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, & dqdt_aq, dqdt_wr, dqdt - real(r8) :: fwetrem, sumf, uptkrate - real(r8) :: delnh3, delnh4 - - integer :: l, n, m - integer :: ntot_msa_c + real(r8) :: fwetrem, uptkrate - integer :: i,k + integer :: l, m, n, mm + integer :: i,k, ndx real(r8) :: xl + real(r8) :: mw_so4 + character(len=32) :: spectype + character(len=32) :: specname ! make sure dqdt is zero initially, for budgets dqdt_aqso4(:,:,:) = 0.0_r8 @@ -233,84 +214,46 @@ subroutine sox_cldaero_update( & dqdt_aqhprxn(:,:) = 0.0_r8 dqdt_aqo3rxn(:,:) = 0.0_r8 + aqso4 = 0.0_r8 + aqh2so4 = 0.0_r8 + aqso4_h2o2 = 0.0_r8 + aqso4_o3 = 0.0_r8 + delso4_3d = 0.0_r8 + ! Avoid double counting in-cloud sulfur oxidation when running with ! GEOS-Chem. If running with GEOS-Chem then sulfur oxidation ! is performed internally to GEOS-Chem. Here, we just return to the ! parent routine and thus we do not apply tendencies calculated by MAM. if ( cam_chempkg_is('geoschem_mam4') ) return + where (cldfrc(:ncol,:) >= 1.0e-5_r8) + delso4_3d(:ncol,:) = xso4(:ncol,:) - xso4_init(:ncol,:) + end where + + !------------------------------------------------------------------------- + ! Compute factors for partitioning aerosol mass gains among bins / modes. + ! The factors are proportional to the activated particle MR for each + ! bin, which is the MR of cloud drops "associated with" the mode + ! thus we are assuming the cloud drop size is independent of the + ! associated aerosol mode properties + call aero_state%aqu_gain_binfraction(aero_props, 'sulfate', qcw, delso4_3d, faqgain_so4) + if (has_msa) call aero_state%aqu_gain_binfraction(aero_props, 'msa', qcw, delso4_3d, faqgain_msa) + lev_loop: do k = 1,pver col_loop: do i = 1,ncol cloud: if (cldfrc(i,k) >= 1.0e-5_r8) then - xl = xlwc(i,k) ! / cldfrc(i,k) + xl = xlwc(i,k) - IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED + if (xl .ge. 1.e-8_r8) then !! when cloud is present - delso4_o3rxn = xso4(i,k) - xso4_init(i,k) + delso4_o3rxn = delso4_3d(i,k) ! xso4(i,k) - xso4_init(i,k) if (id_nh3>0) then delnh3 = nh3g(i,k) - xnh3(i,k) delnh4 = - delnh3 endif - !------------------------------------------------------------------------- - ! compute factors for partitioning aerosol mass gains among modes - ! the factors are proportional to the activated particle MR for each - ! mode, which is the MR of cloud drops "associated with" the mode - ! thus we are assuming the cloud drop size is independent of the - ! associated aerosol mode properties (i.e., drops associated with - ! Aitken and coarse sea-salt particles are same size) - ! - ! qnum_c(n) = activated particle number MR for mode n (these are just - ! used for partitioning among modes, so don't need to divide by cldfrc) - - do n = 1, ntot_amode - qnum_c(n) = 0.0_r8 - l = numptrcw_amode(n) - loffset - if (l > 0) qnum_c(n) = max( 0.0_r8, qcw(i,k,l) ) - end do - - ! force qnum_c(n) to be positive for n=modeptr_accum or n=1 - n = modeptr_accum - if (n <= 0) n = 1 - qnum_c(n) = max( 1.0e-10_r8, qnum_c(n) ) - - ! faqgain_so4(n) = fraction of total so4_c gain going to mode n - ! these are proportional to the activated particle MR for each mode - sumf = 0.0_r8 - do n = 1, ntot_amode - faqgain_so4(n) = 0.0_r8 - if (lptr_so4_cw_amode(n) > 0) then - faqgain_so4(n) = qnum_c(n) - sumf = sumf + faqgain_so4(n) - end if - end do - - if (sumf > 0.0_r8) then - do n = 1, ntot_amode - faqgain_so4(n) = faqgain_so4(n) / sumf - end do - end if - ! at this point (sumf <= 0.0) only when all the faqgain_so4 are zero - ! faqgain_msa(n) = fraction of total msa_c gain going to mode n - ntot_msa_c = 0 - sumf = 0.0_r8 - do n = 1, ntot_amode - faqgain_msa(n) = 0.0_r8 - if (lptr_msa_cw_amode(n) > 0) then - faqgain_msa(n) = qnum_c(n) - ntot_msa_c = ntot_msa_c + 1 - end if - sumf = sumf + faqgain_msa(n) - end do - - if (sumf > 0.0_r8) then - do n = 1, ntot_amode - faqgain_msa(n) = faqgain_msa(n) / sumf - end do - end if - ! at this point (sumf <= 0.0) only when all the faqgain_msa are zero uptkrate = cldaero_uptakerate( xl, cldnum(i,k), cfact(i,k), cldfrc(i,k), tfld(i,k), press(i,k) ) ! average uptake rate over dtime @@ -319,16 +262,18 @@ subroutine sox_cldaero_update( & ! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s) ! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s) dso4dt_gasuptk = xh2so4(i,k) * uptkrate - if (id_msa > 0) then + if (has_msa) then dmsadt_gasuptk = xmsa(i,k) * uptkrate else dmsadt_gasuptk = 0.0_r8 end if ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4 - dmsadt_gasuptk_toso4 = 0.0_r8 - dmsadt_gasuptk_tomsa = dmsadt_gasuptk - if (ntot_msa_c == 0) then + if (has_msa) then + dmsadt_gasuptk_toso4 = 0.0_r8 + dmsadt_gasuptk_tomsa = dmsadt_gasuptk + else + ! no MSA dmsadt_gasuptk_tomsa = 0.0_r8 dmsadt_gasuptk_toso4 = dmsadt_gasuptk end if @@ -347,38 +292,40 @@ subroutine sox_cldaero_update( & fwetrem = 0.0_r8 ! don't have so4 & msa wet removal here ! compute TMR tendencies for so4 and msa aerosol-in-cloud-water - do n = 1, ntot_amode - l = lptr_so4_cw_amode(n) - loffset - if (l > 0) then - dqdt_aqso4(i,k,l) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k) - dqdt_aqh2so4(i,k,l) = faqgain_so4(n)* & - (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k) - dqdt_aq = dqdt_aqso4(i,k,l) + dqdt_aqh2so4(i,k,l) - dqdt_wr = -fwetrem*dqdt_aq - dqdt= dqdt_aq + dqdt_wr - qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime - end if - - l = lptr_msa_cw_amode(n) - loffset - if (l > 0) then - dqdt_aq = faqgain_msa(n)*dmsadt_gasuptk_tomsa*cldfrc(i,k) - dqdt_wr = -fwetrem*dqdt_aq - dqdt = dqdt_aq + dqdt_wr - qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime - end if - - l = lptr_nh4_cw_amode(n) - loffset - if (l > 0) then - if (delnh4 > 0.0_r8) then - dqdt_aq = faqgain_so4(n)*delnh4/dtime*cldfrc(i,k) - dqdt = dqdt_aq - qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime - else - dqdt = (qcw(i,k,l)/max(xnh4c(i,k),1.0e-35_r8)) & - *delnh4/dtime*cldfrc(i,k) - qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime - endif - end if + do m = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + call aero_props%get(m,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + + dqdt_aqso4(i,k,mm) = faqgain_so4(m,i,k)*dso4dt_aqrxn*cldfrc(i,k) + + dqdt_aqh2so4(i,k,mm) = faqgain_so4(m,i,k)* & + (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k) + dqdt_aq = dqdt_aqso4(i,k,mm) + dqdt_aqh2so4(i,k,mm) + dqdt_wr = -fwetrem*dqdt_aq + dqdt = dqdt_aq + dqdt_wr + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + + end if + if (trim(spectype) == 'msa') then + dqdt_aq = faqgain_msa(m,i,k)*dmsadt_gasuptk_tomsa*cldfrc(i,k) + dqdt_wr = -fwetrem*dqdt_aq + dqdt = dqdt_aq + dqdt_wr + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + end if + if (trim(spectype) == 'ammonium') then + if (delnh4 > 0.0_r8) then + dqdt_aq = faqgain_so4(m,i,k)*delnh4/dtime*cldfrc(i,k) + dqdt = dqdt_aq + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + else + dqdt = (qcw(i,k,mm)/max(xnh4c(i,k),1.0e-35_r8)) & + *delnh4/dtime*cldfrc(i,k) + qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime + endif + end if + end do end do ! For gas species, tendency includes @@ -390,7 +337,7 @@ subroutine sox_cldaero_update( & ! h2so4 (g) & msa (g) qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k) - if (id_msa > 0) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k) + if (has_msa) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k) ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k) ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) ) @@ -421,7 +368,7 @@ subroutine sox_cldaero_update( & dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k) dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k) - ENDIF !! WHEN CLOUD IS PRESENTED + endif !! when cloud is present endif cloud enddo col_loop enddo lev_loop @@ -431,60 +378,63 @@ subroutine sox_cldaero_update( & !============================================================== do k = 1,pver - do n = 1, ntot_amode - - l = lptr_so4_cw_amode(n) - loffset - if (l > 0) then - qcw(:,k,l) = MAX(qcw(:,k,l), small_value ) - end if - l = lptr_msa_cw_amode(n) - loffset - if (l > 0) then - qcw(:,k,l) = MAX(qcw(:,k,l), small_value ) - end if - l = lptr_nh4_cw_amode(n) - loffset - if (l > 0) then - qcw(:,k,l) = MAX(qcw(:,k,l), small_value ) - end if - + do n = 1, aero_props%nbins() + do l = 1, aero_props%nspecies(n) + mm = aero_props%indexer(n,l) + call aero_props%get(n,l, spectype=spectype) + if (trim(spectype) == 'sulfate') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + if (trim(spectype) == 'msa') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + if (trim(spectype) == 'ammonium') then + qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) + end if + end do end do - qin(:,k,id_so2) = MAX( qin(:,k,id_so2), small_value ) - qin(:,k,id_h2o2) = MAX( qin(:,k,id_h2o2), small_value ) - qin(:,k,id_h2so4) = MAX( qin(:,k,id_h2so4), small_value ) - if ( id_msa > 0 ) qin(:,k,id_msa) = MAX( qin(:,k,id_msa), small_value ) - if ( id_nh3 > 0 ) qin(:,k,id_nh3) = MAX( qin(:,k,id_nh3), small_value ) + qin(:ncol,k,id_so2) = MAX( qin(:ncol,k,id_so2), small_value ) + qin(:ncol,k,id_h2o2) = MAX( qin(:ncol,k,id_h2o2), small_value ) + qin(:ncol,k,id_h2so4) = MAX( qin(:ncol,k,id_h2so4), small_value ) + if ( id_msa > 0 ) qin(:ncol,k,id_msa) = MAX( qin(:ncol,k,id_msa), small_value ) + if ( id_nh3 > 0 ) qin(:ncol,k,id_nh3) = MAX( qin(:ncol,k,id_nh3), small_value ) end do ! diagnostics - - do n = 1, ntot_amode - m = lptr_so4_cw_amode(n) - l = m - loffset - if (l > 0) then - aqso4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,l)*adv_mass(l)/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s + mw_so4 = -huge(1._r8) + + do n = 1, aero_props%nbins() + ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero + do l = 1, aero_props%nspecies(n) + mm = aero_props%indexer(n,l) + call aero_props%get(n,l, spectype=spectype, specname=specname) + if (trim(spectype) == 'sulfate') then + call aero_props%get(n,l, spec_mw=mw_so4) + aqso4(:,n)=0._r8 + do k=1,pver + do i=1,ncol + aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg/m2/s + enddo enddo - enddo - - aqh2so4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,l)*adv_mass(l)/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s + aqh2so4(:,n)=0._r8 + do k=1,pver + do i=1,ncol + aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg/m2/s + enddo enddo - enddo - endif + end if + end do end do aqso4_h2o2(:) = 0._r8 do k=1,pver do i=1,ncol - aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s + aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s enddo enddo @@ -492,8 +442,8 @@ subroutine sox_cldaero_update( & aqso4_h2o2_3d(:,:) = 0._r8 do k=1,pver do i=1,ncol - aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s + aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s enddo enddo end if @@ -501,8 +451,8 @@ subroutine sox_cldaero_update( & aqso4_o3(:)=0._r8 do k=1,pver do i=1,ncol - aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s + aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s enddo enddo @@ -510,8 +460,8 @@ subroutine sox_cldaero_update( & aqso4_o3_3d(:,:)=0._r8 do k=1,pver do i=1,ncol - aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s + aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & + *pdel(i,k)/gravit ! kg SO4 /m2/s enddo enddo end if From 4cf6fa8332546d911ead4c4d1680e90aebcbdcae Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Tue, 5 May 2026 19:06:38 -0600 Subject: [PATCH 2/9] Cleaner MAM mods modified: src/chemistry/modal_aero/aero_model.F90 modified: src/chemistry/modal_aero/modal_aero_gasaerexch.F90 --- src/chemistry/modal_aero/aero_model.F90 | 168 +++++++++--------- .../modal_aero/modal_aero_gasaerexch.F90 | 50 +----- 2 files changed, 91 insertions(+), 127 deletions(-) diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index be1da86975..23bb03bf8e 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -3,7 +3,7 @@ !=============================================================================== module aero_model use shr_kind_mod, only: r8 => shr_kind_r8 - use constituents, only: pcnst, cnst_name, cnst_get_ind, cnst_mw + use constituents, only: pcnst, cnst_name, cnst_get_ind use ppgrid, only: pcols, pver, pverp use phys_control, only: phys_getopts, cam_physpkg_is use cam_abortutils, only: endrun @@ -31,7 +31,7 @@ module aero_model use mo_setsox, only: setsox, has_sox use modal_aerosol_properties_mod, only: modal_aerosol_properties use modal_aerosol_state_mod, only: modal_aerosol_state - use aerosol_state_mod, only: aerosol_state, ptr2d_t + use aerosol_state_mod, only: aerosol_state implicit none private @@ -102,10 +102,10 @@ module aero_model logical :: modal_accum_coarse_exch = .false. type(modal_aerosol_properties), pointer :: aero_props=>null() - integer :: ncnst_tot integer :: n_coarse_dust=-1 ! dmleung added n_coarse_dust to determine the index for the ! coarse dust mode for different MAM versions. 29 Oct 2025 + integer :: ncnst_tot = -1 contains @@ -206,7 +206,7 @@ subroutine aero_model_init( pbuf2d ) logical :: history_aerosol ! Output MAM or SECT aerosol tendencies logical :: history_chemistry, history_cesm_forcing, history_dust - integer :: l, mm + integer :: l character(len=6) :: test_name character(len=64) :: errmes @@ -218,9 +218,6 @@ subroutine aero_model_init( pbuf2d ) character(len=32) :: spec_type character(len=32) :: mode_type integer :: nspec - character(len=32) :: spectype - character(len=32) :: name_a - character(len=32) :: name_c aero_props => modal_aerosol_properties() ncnst_tot = aero_props%ncnst_tot() @@ -993,11 +990,12 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re use modal_aero_gasaerexch, only : modal_aero_gasaerexch_sub use modal_aero_newnuc, only : modal_aero_newnuc_sub use modal_aero_data, only : cnst_name_cw, qqcw_get_field + use mo_chem_utls, only : get_spc_ndx !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- - type(physics_state), target, intent(in) :: state ! Physics state variables + type(physics_state),target,intent(in) :: state ! Physics state variables integer, intent(in) :: loffset ! offset applied to modal aero "pointers" integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: lchnk ! chunk index @@ -1035,9 +1033,9 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8), dimension(ncol) :: wrk character(len=32) :: name - real(r8) :: dvmrcwdt(ncol,pver,ncnst_tot) + real(r8) :: dvmrcwdt(ncol,pver,gas_pcnst) real(r8) :: dvmrdt(ncol,pver,gas_pcnst) - real(r8) :: vmrcw(ncol,pver,ncnst_tot) ! cloud-borne aerosol (vmr) + real(r8) :: vmrcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr) real(r8) :: aqso4(ncol,ntot_amode) ! aqueous phase chemistry real(r8) :: aqh2so4(ncol,ntot_amode) ! aqueous phase chemistry @@ -1048,21 +1046,15 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8), pointer :: fldcw(:,:) real(r8), pointer :: sulfeq(:,:,:) - character(len=32) :: spectype + real(r8) :: qqcw(ncol,pver,ncnst_tot) + + integer :: ndx, mm character(len=32) :: specname character(len=32) :: name_a, name_c - real(r8) :: mw(ncnst_tot) - integer :: ndx, ierr, mm - type(ptr2d_t), allocatable :: raer(:) ! aerosol mass, number mixing ratios - type(ptr2d_t), allocatable :: qqcw(:) class(aerosol_state), pointer :: aero_state - character(len=*), parameter :: subname = 'aero_model_gasaerexch' - -!---------------------------------------------------------------------- aero_state => modal_aerosol_state(state, pbuf) - -!! +! ! ... initialize nh3 ! if ( nh3_ndx > 0 ) then @@ -1098,33 +1090,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! ! Aerosol processes ... ! - allocate( & - raer(ncnst_tot), & - qqcw(ncnst_tot), stat=ierr ) - if (ierr /= 0) call endrun(subname//': allocate error') - - ! Init pointers to mode number and specie mass mixing ratios in - ! intersitial and cloud borne phases. - call aero_state%get_states( aero_props, raer, qqcw ) - - mw(:) = 0.0_r8 - do m = 1, aero_props%nbins() ! main loop over aerosol bins - do l = 0, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - if (l==0) then - call aero_props%num_names(m, name_a, name_c) - else - call aero_props%mmr_names(m,l, name_a, name_c) - end if - - call cnst_get_ind(name_a,ndx) - mw(mm) = cnst_mw(ndx) - vmrcw(:ncol,:,mm) = qqcw(mm)%fld(:ncol,:) - - end do - end do - - call qqcw2vmr( vmrcw, mw, mbar, ncol ) + call qqcw2vmr( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) dvmrdt(:ncol,:,:) = vmr(:ncol,:,:) dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:) @@ -1132,6 +1098,20 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! aqueous chemistry ... if( has_sox ) then + + do m = 1,aero_props%nbins() + do l = 0,aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + ndx = get_spc_ndx( name_a ) + qqcw(:,:,mm) = vmrcw(:,:,ndx) + end do + end do + call setsox( aero_state, state, & pbuf, & ncol, & @@ -1146,7 +1126,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re cldfr, & cldnum, & invariants, & - vmrcw, & + qqcw, & vmr, & xphlwc, & aqso4, & @@ -1155,6 +1135,19 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re aqso4_o3 & ) + do m = 1,aero_props%nbins() + do l = 0,aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + ndx = get_spc_ndx( name_a ) + vmrcw(:,:,ndx) = qqcw(:,:,mm) + end do + end do + do n = 1, ntot_amode l = lptr_so4_cw_amode(n) if (l > 0) then @@ -1240,14 +1233,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re call t_stopf('modal_coag') - call vmr2qqcw( vmrcw, mw, mbar, ncol ) - - do m = 1, aero_props%nbins() ! main loop over aerosol bins - do l = 0, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - qqcw(mm)%fld(:ncol,:) = vmrcw(:ncol,:,mm) - end do - end do + call vmr2qqcw( lchnk, vmrcw, mbar, ncol, loffset, pbuf ) ! diagnostics for cloud-borne aerosols... do n = 1,pcnst @@ -2054,63 +2040,83 @@ end subroutine calc_1_impact_rate !============================================================================= !============================================================================= - subroutine qqcw2vmr(vmr, mw, mbar, ncol) + subroutine qqcw2vmr(lchnk, vmr, mbar, ncol, im, pbuf) + use modal_aero_data, only : qqcw_get_field + use physics_buffer, only : physics_buffer_desc !----------------------------------------------------------------- ! ... Xfrom from mass to volume mixing ratio !----------------------------------------------------------------- + use chem_mods, only : adv_mass, gas_pcnst + + implicit none + !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- - integer, intent(in) :: ncol + integer, intent(in) :: lchnk, ncol, im real(r8), intent(in) :: mbar(ncol,pver) - real(r8), intent(in) :: mw(ncnst_tot) - real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot) + real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst) + type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- - integer :: k,l, m,mm + integer :: k, m + real(r8), pointer :: fldcw(:,:) - do m = 1, aero_props%nbins() - do l = 0, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - do k=1,pver - vmr(:ncol,k,mm) = mbar(:ncol,k) * vmr(:ncol,k,mm) / mw(mm) - end do - end do + do m=1,gas_pcnst + if( adv_mass(m) /= 0._r8 ) then + fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) + if(associated(fldcw)) then + do k=1,pver + vmr(:ncol,k,m) = mbar(:ncol,k) * fldcw(:ncol,k) / adv_mass(m) + end do + else + vmr(:,:,m) = 0.0_r8 + end if + end if end do - end subroutine qqcw2vmr !============================================================================= !============================================================================= - subroutine vmr2qqcw(vmr, mw, mbar, ncol) + subroutine vmr2qqcw( lchnk, vmr, mbar, ncol, im, pbuf ) !----------------------------------------------------------------- ! ... Xfrom from volume to mass mixing ratio !----------------------------------------------------------------- + use m_spc_id + use chem_mods, only : adv_mass, gas_pcnst + use modal_aero_data, only : qqcw_get_field + use physics_buffer, only : physics_buffer_desc + + implicit none + !----------------------------------------------------------------- ! ... Dummy args !----------------------------------------------------------------- - integer, intent(in) :: ncol + integer, intent(in) :: lchnk, ncol, im real(r8), intent(in) :: mbar(ncol,pver) - real(r8), intent(in) :: mw(ncnst_tot) - real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot) + real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst) + type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------- - integer :: k, l, m, mm - - do m = 1, aero_props%nbins() - do l = 0, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - do k=1,pver - vmr(:ncol,k,mm) = mw(mm) * vmr(:ncol,k,mm) / mbar(:ncol,k) + integer :: k, m + real(r8), pointer :: fldcw(:,:) + !----------------------------------------------------------------- + ! ... The non-group species + !----------------------------------------------------------------- + do m = 1,gas_pcnst + fldcw => qqcw_get_field(pbuf, m+im,lchnk,errorhandle=.true.) + if( adv_mass(m) /= 0._r8 .and. associated(fldcw)) then + do k = 1,pver + fldcw(:ncol,k) = adv_mass(m) * vmr(:ncol,k,m) / mbar(:ncol,k) end do - end do + end if end do end subroutine vmr2qqcw diff --git a/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 b/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 index 8b5649b79b..8b36461075 100644 --- a/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 +++ b/src/chemistry/modal_aero/modal_aero_gasaerexch.F90 @@ -18,7 +18,6 @@ module modal_aero_gasaerexch use ppgrid, only: pcols, pver use modal_aero_data, only: ntot_amode, numptr_amode, sigmag_amode use modal_aero_data, only: lptr2_soa_g_amode, lptr2_soa_a_amode, lptr2_pom_a_amode - use modal_aerosol_properties_mod, only: modal_aerosol_properties implicit none private @@ -59,10 +58,6 @@ module modal_aero_gasaerexch real (r8), allocatable :: fac_m2v_pcarbon(:) - ! local indexing for bins - integer :: ncnst_tot ! total number of mode number conc + mode species - type(modal_aerosol_properties), pointer :: aero_props =>null() - ! !DESCRIPTION: This module implements ... ! ! !REVISION HISTORY: @@ -93,8 +88,8 @@ subroutine modal_aero_gasaerexch_sub( & loffset, deltat, & t, pmid, pdel, & qh2o, troplev, & - q, qqcw_in, & - dqdt_other, dqqcwdt_other_in, & + q, qqcw, & + dqdt_other, dqqcwdt_other, & dgncur_a, dgncur_awet, & sulfeq ) @@ -114,7 +109,6 @@ subroutine modal_aero_gasaerexch_sub( & use cam_abortutils, only: endrun use spmd_utils, only: iam, masterproc use phys_control, only: cam_chempkg_is -use mo_chem_utls, only: get_spc_ndx implicit none @@ -130,13 +124,13 @@ subroutine modal_aero_gasaerexch_sub( & ! *** MUST BE #/kmol-air for number ! *** MUST BE mol/mol-air for mass ! *** NOTE ncol dimension - real(r8), intent(inout) :: qqcw_in(ncol,pver,ncnst_tot) + real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx) ! like q but for cloud-borner tracers real(r8), intent(in) :: dqdt_other(ncol,pver,pcnstxx) ! TMR tendency from other continuous ! growth processes (aqchem, soa??) ! *** NOTE ncol dimension - real(r8), intent(in) :: dqqcwdt_other_in(ncol,pver,ncnst_tot) + real(r8), intent(in) :: dqqcwdt_other(ncol,pver,pcnstxx) ! like dqdt_other but for cloud-borner tracers real(r8), intent(in) :: t(pcols,pver) ! temperature at model levels (K) real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa) @@ -263,25 +257,6 @@ subroutine modal_aero_gasaerexch_sub( & real(r8) :: dqdtsv1(ncol,pver,pcnstxx) real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx) - real(r8) :: qqcw(ncol,pver,pcnstxx) - real(r8) :: dqqcwdt_other(ncol,pver,pcnstxx) - - integer :: m, mm, ndx - character(len=32) :: name_a, name_c - - do m = 1,aero_props%nbins() - do l = 0,aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - if (l==0) then - call aero_props%num_names(m, name_a, name_c) - else - call aero_props%mmr_names(m,l, name_a, name_c) - end if - ndx = get_spc_ndx( name_a ) - qqcw(:,:,ndx) = qqcw_in(:,:,mm) - dqqcwdt_other(:,:,ndx) = dqqcwdt_other_in(:,:,mm) - end do - end do !---------------------------------------------------------------------- @@ -969,19 +944,6 @@ subroutine modal_aero_gasaerexch_sub( & end do ! jsrf = ... end do ! l = ... - do m = 1,aero_props%nbins() - do l = 0,aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - if (l==0) then - call aero_props%num_names(m, name_a, name_c) - else - call aero_props%mmr_names(m,l, name_a, name_c) - end if - ndx = get_spc_ndx( name_a ) - qqcw_in(:,:,mm) = qqcw(:,:,ndx) - end do - end do - return end subroutine modal_aero_gasaerexch_sub @@ -1521,10 +1483,6 @@ subroutine modal_aero_gasaerexch_init logical :: history_aerocom ! Output the aerocom history !----------------------------------------------------------------------- - aero_props => modal_aerosol_properties() - - ncnst_tot = aero_props%ncnst_tot() - call phys_getopts( history_aerosol_out = history_aerosol ) maxspec_pcage = nspec_max From 6f78a51f3550df107c7e2c8e1f82ae80e10b2594 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 6 May 2026 07:43:34 -0600 Subject: [PATCH 3/9] minor changes modified: src/chemistry/aerosol/aerosol_properties_mod.F90 modified: src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 modified: src/chemistry/aerosol/carma_aerosol_properties_mod.F90 modified: src/chemistry/aerosol/modal_aerosol_properties_mod.F90 --- src/chemistry/aerosol/aerosol_properties_mod.F90 | 4 ++-- src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 | 6 +++--- src/chemistry/aerosol/carma_aerosol_properties_mod.F90 | 4 ++-- src/chemistry/aerosol/modal_aerosol_properties_mod.F90 | 4 ++-- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/chemistry/aerosol/aerosol_properties_mod.F90 b/src/chemistry/aerosol/aerosol_properties_mod.F90 index 5e7642a198..db2f740226 100644 --- a/src/chemistry/aerosol/aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/aerosol_properties_mod.F90 @@ -109,8 +109,8 @@ end function aero_number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, & - spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + 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 class(aerosol_properties), intent(in) :: self diff --git a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 index c6e54c64ad..4b5e3552c7 100644 --- a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 @@ -152,16 +152,16 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, & + spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) class(bulk_aerosol_properties), intent(in) :: self integer, intent(in) :: bin_ndx ! bin index integer, intent(in) :: species_ndx ! species index - real(r8), optional, intent(out) :: spec_mw ! species molecular weight 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 diff --git a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 index 8ee5a3ba18..718adfda60 100644 --- a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 @@ -264,8 +264,8 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + 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 diff --git a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 index 9b2728b9c8..191e62118f 100644 --- a/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_properties_mod.F90 @@ -407,8 +407,8 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & - spec_mw, spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & + 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 From 3418224f75a5f66bd2929698b7fbeaebb664ed99 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 6 May 2026 09:00:45 -0600 Subject: [PATCH 4/9] fixes for bulk aerosols, pass in aero_props, use model_is routine for mam7 modified: src/chemistry/aerosol/mo_setsox.F90 modified: src/chemistry/bulk_aero/aero_model.F90 modified: src/chemistry/bulk_aero/sox_cldaero_mod.F90 modified: src/chemistry/carma_aero/aero_model.F90 modified: src/chemistry/carma_aero/sox_cldaero_mod.F90 modified: src/chemistry/modal_aero/aero_model.F90 modified: src/chemistry/modal_aero/sox_cldaero_mod.F90 --- src/chemistry/aerosol/mo_setsox.F90 | 9 ++++++--- src/chemistry/bulk_aero/aero_model.F90 | 17 +++++++++++++--- src/chemistry/bulk_aero/sox_cldaero_mod.F90 | 8 ++++++-- src/chemistry/carma_aero/aero_model.F90 | 2 +- src/chemistry/carma_aero/sox_cldaero_mod.F90 | 21 +++++++++++--------- src/chemistry/modal_aero/aero_model.F90 | 2 +- src/chemistry/modal_aero/sox_cldaero_mod.F90 | 21 +++++++++++--------- 7 files changed, 52 insertions(+), 28 deletions(-) diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index 2d404c5226..adc1169ca4 100644 --- a/src/chemistry/aerosol/mo_setsox.F90 +++ b/src/chemistry/aerosol/mo_setsox.F90 @@ -29,7 +29,7 @@ module mo_setsox !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine sox_inti + subroutine sox_inti(aero_props) !----------------------------------------------------------------------- ! ... initialize the hetero sox routine !----------------------------------------------------------------------- @@ -39,6 +39,9 @@ subroutine sox_inti use phys_control, only : phys_getopts use carma_flags_mod, only : carma_do_cloudborne use sox_cldaero_mod, only : sox_cldaero_init + use aerosol_properties_mod, only : aerosol_properties + + class(aerosol_properties), intent(in) :: aero_props logical :: modal_aerosols @@ -143,7 +146,7 @@ subroutine sox_inti return end if - call sox_cldaero_init() + call sox_cldaero_init(aero_props) end subroutine sox_inti @@ -365,7 +368,7 @@ subroutine setsox( aero_state, state, & xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 xso4_init = 0._r8 - + do k = 1,pver xph(:,k) = xph0 ! initial PH value diff --git a/src/chemistry/bulk_aero/aero_model.F90 b/src/chemistry/bulk_aero/aero_model.F90 index 798a59e7f3..e1d81ba3d0 100644 --- a/src/chemistry/bulk_aero/aero_model.F90 +++ b/src/chemistry/bulk_aero/aero_model.F90 @@ -19,6 +19,7 @@ module aero_model use physics_buffer, only: pbuf_get_field, pbuf_get_index use cam_history, only: outfld use infnan, only: nan, assignment(=) + use bulk_aerosol_properties_mod, only: bulk_aerosol_properties implicit none private @@ -57,6 +58,8 @@ module aero_model real(r8) :: aer_sol_factb(pcnst) ! below-cloud solubility factor real(r8) :: aer_scav_coef(pcnst) + type(bulk_aerosol_properties), pointer :: aero_props =>null() + contains !============================================================================= @@ -151,8 +154,10 @@ subroutine aero_model_init( pbuf2d ) logical :: history_aerosol ! Output MAM or SECT aerosol tendencies logical :: history_dust ! Output dust + aero_props => bulk_aerosol_properties() + ! aqueous chem initialization - call sox_inti() + call sox_inti(aero_props) call phys_getopts( history_aerosol_out = history_aerosol,& history_dust_out = history_dust ) @@ -1030,11 +1035,13 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re use mo_aerosols, only : aerosols_formation, has_aerosols use mo_setsox, only : setsox, has_sox use mo_setsoa, only : setsoa, has_soa + use aerosol_state_mod, only : aerosol_state + use bulk_aerosol_state_mod, only : bulk_aerosol_state !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- - type(physics_state), intent(in) :: state ! Physics state variables + type(physics_state),target, intent(in) :: state ! Physics state variables integer, intent(in) :: loffset ! offset applied to modal aero "pointers" integer, intent(in) :: ncol ! number columns in chunk integer, intent(in) :: lchnk ! chunk index @@ -1069,11 +1076,15 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3 real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc + class(aerosol_state), pointer :: aero_state + +!---------------------------------------------------------------------- + aero_state => bulk_aerosol_state(state, pbuf) ! aqueous chemistry ... if( has_sox ) then - call setsox( state, & + call setsox( aero_state, state, & pbuf, & ncol, & lchnk, & diff --git a/src/chemistry/bulk_aero/sox_cldaero_mod.F90 b/src/chemistry/bulk_aero/sox_cldaero_mod.F90 index 0c5a7cc923..86aa88e6f6 100644 --- a/src/chemistry/bulk_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/bulk_aero/sox_cldaero_mod.F90 @@ -8,6 +8,8 @@ module sox_cldaero_mod use ppgrid, only : pcols, pver use mo_chem_utls, only : get_spc_ndx use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate + use aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state implicit none private @@ -26,7 +28,8 @@ module sox_cldaero_mod !---------------------------------------------------------------------------------- !---------------------------------------------------------------------------------- - subroutine sox_cldaero_init + subroutine sox_cldaero_init(aero_props_in) + class(aerosol_properties), target, intent(in) :: aero_props_in id_so2 = get_spc_ndx( 'SO2' ) id_so4 = get_spc_ndx( 'SO4' ) @@ -60,7 +63,7 @@ end function sox_cldaero_create_obj !---------------------------------------------------------------------------------- ! Update the mixing ratios !---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( & + subroutine sox_cldaero_update( aero_state, & state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d ) @@ -68,6 +71,7 @@ subroutine sox_cldaero_update( & ! args + class(aerosol_state), intent(in) :: aero_state type(physics_state), intent(in) :: state ! Physics state variables integer, intent(in) :: ncol integer, intent(in) :: lchnk ! chunk id diff --git a/src/chemistry/carma_aero/aero_model.F90 b/src/chemistry/carma_aero/aero_model.F90 index c53669355e..a1952acb9d 100644 --- a/src/chemistry/carma_aero/aero_model.F90 +++ b/src/chemistry/carma_aero/aero_model.F90 @@ -236,7 +236,7 @@ subroutine aero_model_init( pbuf2d ) end if ! aqueous chem initialization - call sox_inti() + call sox_inti(aero_props) h2so4_ndx = get_spc_ndx('H2SO4') nh3_ndx = get_spc_ndx('NH3') diff --git a/src/chemistry/carma_aero/sox_cldaero_mod.F90 b/src/chemistry/carma_aero/sox_cldaero_mod.F90 index 3d2d0845f8..ef5bd77a14 100644 --- a/src/chemistry/carma_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/carma_aero/sox_cldaero_mod.F90 @@ -12,11 +12,9 @@ module sox_cldaero_mod use phys_control, only : cam_chempkg_is use cldaero_mod, only : cldaero_uptakerate use chem_mods, only : gas_pcnst - use carma_aerosol_properties_mod, only: carma_aerosol_properties + use aerosol_properties_mod, only: aerosol_properties use aerosol_state_mod, only: aerosol_state - use modal_aero_data, only : ntot_amode - implicit none private @@ -32,7 +30,7 @@ module sox_cldaero_mod integer :: ncnst_tot = -huge(1) ! total number of mode number conc + mode species integer, public, protected :: nbins = 0 - type(carma_aerosol_properties), pointer :: aero_props =>null() + class(aerosol_properties), pointer :: aero_props =>null() logical :: has_msa = .false. @@ -41,7 +39,8 @@ module sox_cldaero_mod !---------------------------------------------------------------------------------- !---------------------------------------------------------------------------------- - subroutine sox_cldaero_init + subroutine sox_cldaero_init(aero_props_in) + class(aerosol_properties), target, intent(in) :: aero_props_in id_msa = get_spc_ndx( 'MSA' ) id_h2so4 = get_spc_ndx( 'H2SO4' ) @@ -55,7 +54,7 @@ subroutine sox_cldaero_init //' -- should not invoke sox_cldaero_mod ') endif - aero_props => carma_aerosol_properties() + aero_props => aero_props_in ncnst_tot = aero_props%ncnst_tot() nbins = aero_props%nbins() @@ -76,12 +75,16 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( type(cldaero_conc_t), pointer :: conc_obj character(len=32) :: spectype - integer :: l,m - integer :: i,k,mm - + integer :: i,k,mm, ntot_amode logical :: mode7 + if (aero_props%model_is('MAM')) then + ntot_amode = aero_props%nbins() + else + ntot_amode = 0 + end if + mode7 = ntot_amode == 7 conc_obj => cldaero_allocate() diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index 23bb03bf8e..a039d7ff71 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -223,7 +223,7 @@ subroutine aero_model_init( pbuf2d ) ncnst_tot = aero_props%ncnst_tot() ! aqueous chem initialization - call sox_inti() + call sox_inti(aero_props) dgnum_idx = pbuf_get_index('DGNUM') dgnumwet_idx = pbuf_get_index('DGNUMWET') diff --git a/src/chemistry/modal_aero/sox_cldaero_mod.F90 b/src/chemistry/modal_aero/sox_cldaero_mod.F90 index 1701d4a6dd..46d3ed839a 100644 --- a/src/chemistry/modal_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/modal_aero/sox_cldaero_mod.F90 @@ -12,11 +12,9 @@ module sox_cldaero_mod use phys_control, only : cam_chempkg_is use cldaero_mod, only : cldaero_uptakerate use chem_mods, only : gas_pcnst - use modal_aerosol_properties_mod, only: modal_aerosol_properties + use aerosol_properties_mod, only: aerosol_properties use aerosol_state_mod, only: aerosol_state - use modal_aero_data, only : ntot_amode - implicit none private @@ -32,7 +30,7 @@ module sox_cldaero_mod integer :: ncnst_tot = -huge(1) ! total number of mode number conc + mode species integer, public, protected :: nbins = 0 - type(modal_aerosol_properties), pointer :: aero_props =>null() + class(aerosol_properties), pointer :: aero_props =>null() logical :: has_msa = .false. @@ -41,7 +39,8 @@ module sox_cldaero_mod !---------------------------------------------------------------------------------- !---------------------------------------------------------------------------------- - subroutine sox_cldaero_init + subroutine sox_cldaero_init(aero_props_in) + class(aerosol_properties), target, intent(in) :: aero_props_in id_msa = get_spc_ndx( 'MSA' ) id_h2so4 = get_spc_ndx( 'H2SO4' ) @@ -55,7 +54,7 @@ subroutine sox_cldaero_init //' -- should not invoke sox_cldaero_mod ') endif - aero_props => modal_aerosol_properties() + aero_props => aero_props_in ncnst_tot = aero_props%ncnst_tot() nbins = aero_props%nbins() @@ -76,12 +75,16 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( type(cldaero_conc_t), pointer :: conc_obj character(len=32) :: spectype - integer :: l,m - integer :: i,k,mm - + integer :: i,k,mm, ntot_amode logical :: mode7 + if (aero_props%model_is('MAM')) then + ntot_amode = aero_props%nbins() + else + ntot_amode = 0 + end if + mode7 = ntot_amode == 7 conc_obj => cldaero_allocate() From 4ae8f26fbc16c64a726454d74544ef2c7c3ac834 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Wed, 6 May 2026 20:12:10 -0600 Subject: [PATCH 5/9] single sox_cldaero_mod module file renamed: src/chemistry/modal_aero/sox_cldaero_mod.F90 -> src/chemistry/aerosol/sox_cldaero_mod.F90 deleted: src/chemistry/bulk_aero/sox_cldaero_mod.F90 deleted: src/chemistry/carma_aero/sox_cldaero_mod.F90 modified: src/chemistry/aerosol/mo_setsox.F90 modified: src/chemistry/aerosol/sox_cldaero_mod.F90 --- src/chemistry/aerosol/mo_setsox.F90 | 25 +- .../sox_cldaero_mod.F90 | 19 +- src/chemistry/bulk_aero/sox_cldaero_mod.F90 | 147 ------ src/chemistry/modal_aero/sox_cldaero_mod.F90 | 483 ------------------ 4 files changed, 33 insertions(+), 641 deletions(-) rename src/chemistry/{carma_aero => aerosol}/sox_cldaero_mod.F90 (97%) delete mode 100644 src/chemistry/bulk_aero/sox_cldaero_mod.F90 delete mode 100644 src/chemistry/modal_aero/sox_cldaero_mod.F90 diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index adc1169ca4..a32b4afaea 100644 --- a/src/chemistry/aerosol/mo_setsox.F90 +++ b/src/chemistry/aerosol/mo_setsox.F90 @@ -273,6 +273,7 @@ subroutine setsox( aero_state, state, & real(r8), parameter :: kh2 = 8.3e5_r8 ! HO2(a) + ho2(a) -> h2o2(a) + o2 Reference: JPL; Bielski et al. 1985 real(r8), parameter :: kh3 = 9.7e7_r8 ! HO2(a) + o2- -> h2o2(a) + o2 Reference: JPL; Bielski et al. 1985 real(r8), parameter :: Ra = GAS_CONSTANT_KMOL / KMOL_TO_MOL * M3_TO_L * PASCAL_TO_ATM ! universal constant (atm)/(M-K) + real(r8), parameter :: small_value = 1.e-20_r8 ! real(r8) :: xdelso4hp(ncol,pver) @@ -883,10 +884,26 @@ subroutine setsox( aero_state, state, & end do col_loop1 end do ver_loop1 - call sox_cldaero_update( aero_state, & - state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & - xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & - aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d ) + aqso4 = 0._r8 + aqh2so4 = 0._r8 + aqso4_h2o2 = 0._r8 + aqso4_o3 = 0._r8 + + if (cloud_borne) then + ! update cloud-borne aerosols + call sox_cldaero_update( aero_state, & + state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & + xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & + aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d ) + else + if (id_so2>0) then + qin(:ncol,:,id_so2) = max( xso2(:ncol,:), small_value ) + endif + if (id_h2o2>0) then + qin(:ncol,:,id_h2o2) = max( xh2o2(:ncol,:), small_value ) + endif + qin(:ncol,:,id_so4) = max( xso4(:ncol,:), small_value ) + endif xphlwc(:,:) = 0._r8 do k = 1, pver diff --git a/src/chemistry/carma_aero/sox_cldaero_mod.F90 b/src/chemistry/aerosol/sox_cldaero_mod.F90 similarity index 97% rename from src/chemistry/carma_aero/sox_cldaero_mod.F90 rename to src/chemistry/aerosol/sox_cldaero_mod.F90 index ef5bd77a14..bd8fbf32d9 100644 --- a/src/chemistry/carma_aero/sox_cldaero_mod.F90 +++ b/src/chemistry/aerosol/sox_cldaero_mod.F90 @@ -1,5 +1,5 @@ !---------------------------------------------------------------------------------- -! CARMA implementation +! Generic aerosol implementation !---------------------------------------------------------------------------------- module sox_cldaero_mod @@ -49,9 +49,8 @@ subroutine sox_cldaero_init(aero_props_in) id_nh3 = get_spc_ndx( 'NH3' ) has_msa = id_msa>0 - if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then - call endrun('sox_cldaero_init:MAM mech does not include necessary species' & - //' -- should not invoke sox_cldaero_mod ') + if ( id_so2<1 ) then + call endrun('sox_cldaero_init: SO2 is not included in chemistry -- should not invoke sox_cldaero_mod...') endif aero_props => aero_props_in @@ -79,6 +78,14 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( integer :: i,k,mm, ntot_amode logical :: mode7 + conc_obj => cldaero_allocate() + + if (aero_props%model_is('BAM')) then + ! no cloud-borne aerosols + conc_obj%xlwc(:ncol,:) = lwc(:ncol,:)*cfact(:ncol,:) ! cloud water L(water)/ + return + end if + if (aero_props%model_is('MAM')) then ntot_amode = aero_props%nbins() else @@ -87,11 +94,9 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( mode7 = ntot_amode == 7 - conc_obj => cldaero_allocate() - do k = 1,pver do i = 1,ncol - if( cldfrc(i,k) >0._r8) then + if(cldfrc(i,k)>0._r8) then conc_obj%xlwc(i,k) = lwc(i,k) *cfact(i,k) ! cloud water L(water)/L(air) conc_obj%xlwc(i,k) = conc_obj%xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell else diff --git a/src/chemistry/bulk_aero/sox_cldaero_mod.F90 b/src/chemistry/bulk_aero/sox_cldaero_mod.F90 deleted file mode 100644 index 86aa88e6f6..0000000000 --- a/src/chemistry/bulk_aero/sox_cldaero_mod.F90 +++ /dev/null @@ -1,147 +0,0 @@ -!---------------------------------------------------------------------------------- -! Bulk aerosol implementation -!---------------------------------------------------------------------------------- -module sox_cldaero_mod - - use shr_kind_mod, only : r8 => shr_kind_r8 - use cam_abortutils, only : endrun - use ppgrid, only : pcols, pver - use mo_chem_utls, only : get_spc_ndx - use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate - use aerosol_properties_mod, only: aerosol_properties - use aerosol_state_mod, only: aerosol_state - - implicit none - private - - public :: sox_cldaero_init - public :: sox_cldaero_create_obj - public :: sox_cldaero_update - public :: sox_cldaero_destroy_obj - - integer :: id_so2, id_so4, id_h2o2 - - real(r8), parameter :: small_value = 1.e-20_r8 - -contains - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - - subroutine sox_cldaero_init(aero_props_in) - class(aerosol_properties), target, intent(in) :: aero_props_in - - id_so2 = get_spc_ndx( 'SO2' ) - id_so4 = get_spc_ndx( 'SO4' ) - id_h2o2 = get_spc_ndx( 'H2O2' ) - - if ( id_so2<1 ) then - call endrun('sox_cldaero_init: SO2 is not included in chemistry -- should not invoke sox_cldaero_mod...') - endif - - end subroutine sox_cldaero_init - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj ) - - real(r8), intent(in) :: cldfrc(:,:) - real(r8), intent(in) :: qcw(:,:,:) - real(r8), intent(in) :: lwc(:,:) - real(r8), intent(in) :: cfact(:,:) - integer, intent(in) :: ncol - integer, intent(in) :: loffset - - type(cldaero_conc_t), pointer :: conc_obj - - conc_obj => cldaero_allocate() - - conc_obj%xlwc(:ncol,:) = lwc(:ncol,:)*cfact(:ncol,:) ! cloud water L(water)/L(air) - - end function sox_cldaero_create_obj - -!---------------------------------------------------------------------------------- -! Update the mixing ratios -!---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( aero_state, & - state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & - delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & - aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d ) - use physics_types, only: physics_state - - ! args - - class(aerosol_state), intent(in) :: aero_state - type(physics_state), intent(in) :: state ! Physics state variables - integer, intent(in) :: ncol - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset - - real(r8), intent(in) :: dtime ! time step (sec) - - real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu ) - real(r8), intent(in) :: pdel(:,:) - real(r8), intent(in) :: press(:,:) - real(r8), intent(in) :: tfld(:,:) - - real(r8), intent(in) :: cldnum(:,:) - real(r8), intent(in) :: cldfrc(:,:) - real(r8), intent(in) :: cfact(:,:) - real(r8), intent(in) :: xlwc(:,:) - - real(r8), intent(in) :: delso4_hprxn(:,:) - real(r8), intent(in) :: xh2so4(:,:) - real(r8), intent(in) :: xso4(:,:) - real(r8), intent(in) :: xso4_init(:,:) - real(r8), intent(in) :: nh3g(:,:) - real(r8), intent(in) :: hno3g(:,:) - real(r8), intent(in) :: xnh3(:,:) - real(r8), intent(in) :: xhno3(:,:) - real(r8), intent(in) :: xnh4c(:,:) - real(r8), intent(in) :: xmsa(:,:) - real(r8), intent(in) :: xso2(:,:) - real(r8), intent(in) :: xh2o2(:,:) - real(r8), intent(in) :: xno3c(:,:) - - real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) - real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr ) - - real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - real(r8), intent(out), optional :: aqso4_h2o2_3d(:,:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - - ! local vars ... - - integer :: k - - !============================================================== - ! ... Update the mixing ratios - !============================================================== - do k = 1,pver - - if (id_so2>0) then - qin(:,k,id_so2) = MAX( xso2(:,k), small_value ) - endif - if (id_h2o2>0) then - qin(:,k,id_h2o2)= MAX( xh2o2(:,k), small_value ) - endif - - qin(:,k,id_so4) = MAX( xso4(:,k), small_value ) - - end do - - end subroutine sox_cldaero_update - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - subroutine sox_cldaero_destroy_obj( conc_obj ) - type(cldaero_conc_t), pointer :: conc_obj - - call cldaero_deallocate( conc_obj ) - - end subroutine sox_cldaero_destroy_obj - -end module sox_cldaero_mod diff --git a/src/chemistry/modal_aero/sox_cldaero_mod.F90 b/src/chemistry/modal_aero/sox_cldaero_mod.F90 deleted file mode 100644 index 46d3ed839a..0000000000 --- a/src/chemistry/modal_aero/sox_cldaero_mod.F90 +++ /dev/null @@ -1,483 +0,0 @@ -!---------------------------------------------------------------------------------- -! Modal aerosol implementation -!---------------------------------------------------------------------------------- -module sox_cldaero_mod - - use shr_kind_mod, only : r8 => shr_kind_r8 - use cam_abortutils, only : endrun - use ppgrid, only : pcols, pver - use mo_chem_utls, only : get_spc_ndx - use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate - use physconst, only : gravit - use phys_control, only : cam_chempkg_is - use cldaero_mod, only : cldaero_uptakerate - use chem_mods, only : gas_pcnst - use aerosol_properties_mod, only: aerosol_properties - use aerosol_state_mod, only: aerosol_state - - implicit none - private - - public :: sox_cldaero_init - public :: sox_cldaero_create_obj - public :: sox_cldaero_update - public :: sox_cldaero_destroy_obj - - integer :: id_msa=-1, id_h2so4=-1, id_so2=-1, id_h2o2=-1, id_nh3=-1 - - real(r8), parameter :: small_value = 1.e-20_r8 - - integer :: ncnst_tot = -huge(1) ! total number of mode number conc + mode species - integer, public, protected :: nbins = 0 - - class(aerosol_properties), pointer :: aero_props =>null() - - logical :: has_msa = .false. - -contains - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - - subroutine sox_cldaero_init(aero_props_in) - class(aerosol_properties), target, intent(in) :: aero_props_in - - id_msa = get_spc_ndx( 'MSA' ) - id_h2so4 = get_spc_ndx( 'H2SO4' ) - id_so2 = get_spc_ndx( 'SO2' ) - id_h2o2 = get_spc_ndx( 'H2O2' ) - id_nh3 = get_spc_ndx( 'NH3' ) - has_msa = id_msa>0 - - if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then - call endrun('sox_cldaero_init:MAM mech does not include necessary species' & - //' -- should not invoke sox_cldaero_mod ') - endif - - aero_props => aero_props_in - - ncnst_tot = aero_props%ncnst_tot() - nbins = aero_props%nbins() - - end subroutine sox_cldaero_init - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj ) - - real(r8), intent(in) :: cldfrc(:,:) - real(r8), intent(in) :: qcw(:,:,:) - real(r8), intent(in) :: lwc(:,:) - real(r8), intent(in) :: cfact(:,:) - integer, intent(in) :: ncol - integer, intent(in) :: loffset - - type(cldaero_conc_t), pointer :: conc_obj - - character(len=32) :: spectype - integer :: l,m - integer :: i,k,mm, ntot_amode - logical :: mode7 - - if (aero_props%model_is('MAM')) then - ntot_amode = aero_props%nbins() - else - ntot_amode = 0 - end if - - mode7 = ntot_amode == 7 - - conc_obj => cldaero_allocate() - - do k = 1,pver - do i = 1,ncol - if( cldfrc(i,k) >0._r8) then - conc_obj%xlwc(i,k) = lwc(i,k) *cfact(i,k) ! cloud water L(water)/L(air) - conc_obj%xlwc(i,k) = conc_obj%xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell - else - conc_obj%xlwc(i,k) = 0._r8 - endif - enddo - enddo - - conc_obj%no3c(:,:) = 0._r8 - conc_obj%nh4c(:,:) = 0._r8 - conc_obj%so4c(:,:) = 0._r8 - - do k = 1,pver - do i = 1,ncol - do m = 1, aero_props%nbins() - do l = 1, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - call aero_props%get(m,l, spectype=spectype) - if (trim(spectype) == 'sulfate') then - conc_obj%so4c(i,k) = conc_obj%so4c(i,k) + qcw(i,k,mm) - end if - if (trim(spectype) == 'ammonium') then - conc_obj%nh4c(i,k) = conc_obj%nh4c(i,k) + qcw(i,k,mm) - end if - end do - end do - end do - end do - - ! *** NOTE *** - ! should refactor this bit after merging to later cam tag -- where aero_props%model_is is avail - if (ntot_amode>0) then - if (.not.mode7) then - conc_obj%so4_fact = 1._r8 - end if - end if - - end function sox_cldaero_create_obj - -!---------------------------------------------------------------------------------- -! Update the mixing ratios -!---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( aero_state, & - state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & - delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & - aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d) - - use physics_types, only: physics_state - - ! args - - class(aerosol_state), intent(in) :: aero_state - type(physics_state), intent(in) :: state ! Physics state variables - - integer, intent(in) :: ncol - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset - - real(r8), intent(in) :: dtime ! time step (sec) - - real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu ) - real(r8), intent(in) :: pdel(:,:) - real(r8), intent(in) :: press(:,:) - real(r8), intent(in) :: tfld(:,:) - - real(r8), intent(in) :: cldnum(:,:) - real(r8), intent(in) :: cldfrc(:,:) - real(r8), intent(in) :: cfact(:,:) - real(r8), intent(in) :: xlwc(:,:) - - real(r8), intent(in) :: delso4_hprxn(:,:) - real(r8), intent(in) :: xh2so4(:,:) - real(r8), intent(in) :: xso4(:,:) - real(r8), intent(in) :: xso4_init(:,:) - real(r8), intent(in) :: nh3g(:,:) - real(r8), intent(in) :: hno3g(:,:) - real(r8), intent(in) :: xnh3(:,:) - real(r8), intent(in) :: xhno3(:,:) - real(r8), intent(in) :: xnh4c(:,:) - real(r8), intent(in) :: xmsa(:,:) - real(r8), intent(in) :: xso2(:,:) - real(r8), intent(in) :: xh2o2(:,:) - real(r8), intent(in) :: xno3c(:,:) - - real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) - real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr ) - - real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - real(r8), intent(out), optional :: aqso4_h2o2_3d(:,:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - - ! local vars ... - - real(r8) :: dqdt_aqso4(ncol,pver,ncnst_tot), & - dqdt_aqh2so4(ncol,pver,ncnst_tot), & - dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver) - - real(r8) :: faqgain_msa(nbins,ncol,pver), faqgain_so4(nbins,ncol,pver) - real(r8) :: delso4_3d(ncol,pver) - - real(r8) :: delnh3, delnh4 - real(r8) :: delso4_o3rxn, & - dso4dt_aqrxn, dso4dt_hprxn, & - dso4dt_gasuptk, dmsadt_gasuptk, & - dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, & - dqdt_aq, dqdt_wr, dqdt - - real(r8) :: fwetrem, uptkrate - - integer :: l, m, n, mm - integer :: i,k, ndx - real(r8) :: xl - real(r8) :: mw_so4 - character(len=32) :: spectype - character(len=32) :: specname - - ! make sure dqdt is zero initially, for budgets - dqdt_aqso4(:,:,:) = 0.0_r8 - dqdt_aqh2so4(:,:,:) = 0.0_r8 - dqdt_aqhprxn(:,:) = 0.0_r8 - dqdt_aqo3rxn(:,:) = 0.0_r8 - - aqso4 = 0.0_r8 - aqh2so4 = 0.0_r8 - aqso4_h2o2 = 0.0_r8 - aqso4_o3 = 0.0_r8 - delso4_3d = 0.0_r8 - - ! Avoid double counting in-cloud sulfur oxidation when running with - ! GEOS-Chem. If running with GEOS-Chem then sulfur oxidation - ! is performed internally to GEOS-Chem. Here, we just return to the - ! parent routine and thus we do not apply tendencies calculated by MAM. - if ( cam_chempkg_is('geoschem_mam4') ) return - - where (cldfrc(:ncol,:) >= 1.0e-5_r8) - delso4_3d(:ncol,:) = xso4(:ncol,:) - xso4_init(:ncol,:) - end where - - !------------------------------------------------------------------------- - ! Compute factors for partitioning aerosol mass gains among bins / modes. - ! The factors are proportional to the activated particle MR for each - ! bin, which is the MR of cloud drops "associated with" the mode - ! thus we are assuming the cloud drop size is independent of the - ! associated aerosol mode properties - call aero_state%aqu_gain_binfraction(aero_props, 'sulfate', qcw, delso4_3d, faqgain_so4) - if (has_msa) call aero_state%aqu_gain_binfraction(aero_props, 'msa', qcw, delso4_3d, faqgain_msa) - - lev_loop: do k = 1,pver - col_loop: do i = 1,ncol - cloud: if (cldfrc(i,k) >= 1.0e-5_r8) then - xl = xlwc(i,k) - - if (xl .ge. 1.e-8_r8) then !! when cloud is present - - delso4_o3rxn = delso4_3d(i,k) ! xso4(i,k) - xso4_init(i,k) - - if (id_nh3>0) then - delnh3 = nh3g(i,k) - xnh3(i,k) - delnh4 = - delnh3 - endif - - ! faqgain_msa(n) = fraction of total msa_c gain going to mode n - - uptkrate = cldaero_uptakerate( xl, cldnum(i,k), cfact(i,k), cldfrc(i,k), tfld(i,k), press(i,k) ) - ! average uptake rate over dtime - uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime - - ! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s) - ! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s) - dso4dt_gasuptk = xh2so4(i,k) * uptkrate - if (has_msa) then - dmsadt_gasuptk = xmsa(i,k) * uptkrate - else - dmsadt_gasuptk = 0.0_r8 - end if - - ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4 - if (has_msa) then - dmsadt_gasuptk_toso4 = 0.0_r8 - dmsadt_gasuptk_tomsa = dmsadt_gasuptk - else - ! no MSA - dmsadt_gasuptk_tomsa = 0.0_r8 - dmsadt_gasuptk_toso4 = dmsadt_gasuptk - end if - - !----------------------------------------------------------------------- - ! now compute TMR tendencies - ! this includes the above aqueous so2 chemistry AND - ! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...) - ! AND the wetremoval of dissolved, unreacted so2 and h2o2 - - dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime - dso4dt_hprxn = delso4_hprxn(i,k) / dtime - - ! fwetrem = fraction of in-cloud-water material that is wet removed - ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) ) - fwetrem = 0.0_r8 ! don't have so4 & msa wet removal here - - ! compute TMR tendencies for so4 and msa aerosol-in-cloud-water - do m = 1, aero_props%nbins() - do l = 1, aero_props%nspecies(m) - mm = aero_props%indexer(m,l) - call aero_props%get(m,l, spectype=spectype) - if (trim(spectype) == 'sulfate') then - - dqdt_aqso4(i,k,mm) = faqgain_so4(m,i,k)*dso4dt_aqrxn*cldfrc(i,k) - - dqdt_aqh2so4(i,k,mm) = faqgain_so4(m,i,k)* & - (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k) - dqdt_aq = dqdt_aqso4(i,k,mm) + dqdt_aqh2so4(i,k,mm) - dqdt_wr = -fwetrem*dqdt_aq - dqdt = dqdt_aq + dqdt_wr - qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime - - end if - if (trim(spectype) == 'msa') then - dqdt_aq = faqgain_msa(m,i,k)*dmsadt_gasuptk_tomsa*cldfrc(i,k) - dqdt_wr = -fwetrem*dqdt_aq - dqdt = dqdt_aq + dqdt_wr - qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime - end if - if (trim(spectype) == 'ammonium') then - if (delnh4 > 0.0_r8) then - dqdt_aq = faqgain_so4(m,i,k)*delnh4/dtime*cldfrc(i,k) - dqdt = dqdt_aq - qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime - else - dqdt = (qcw(i,k,mm)/max(xnh4c(i,k),1.0e-35_r8)) & - *delnh4/dtime*cldfrc(i,k) - qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime - endif - end if - end do - end do - - ! For gas species, tendency includes - ! reactive uptake to cloud water that essentially transforms the gas to - ! a different species. Wet removal associated with this is applied - ! to the "new" species (e.g., so4_c) rather than to the gas. - ! wet removal of the unreacted gas that is dissolved in cloud water. - ! Need to multiply both these parts by cldfrc - - ! h2so4 (g) & msa (g) - qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k) - if (has_msa) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k) - - ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k) - ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) ) - fwetrem = 0.0_r8 ! don't include so2 wet removal here - - dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k) - dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k) - dqdt = dqdt_aq + dqdt_wr - qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime - - ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k) - ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) ) - fwetrem = 0.0_r8 ! don't include h2o2 wet removal here - - dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k) - dqdt_aq = -dso4dt_hprxn*cldfrc(i,k) - dqdt = dqdt_aq + dqdt_wr - qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime - - ! NH3 - if (id_nh3>0) then - dqdt_aq = delnh3/dtime*cldfrc(i,k) - dqdt = dqdt_aq - qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime - endif - - ! for SO4 from H2O2/O3 budgets - dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k) - dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k) - - endif !! when cloud is present - endif cloud - enddo col_loop - enddo lev_loop - - !============================================================== - ! ... Update the mixing ratios - !============================================================== - do k = 1,pver - - do n = 1, aero_props%nbins() - do l = 1, aero_props%nspecies(n) - mm = aero_props%indexer(n,l) - call aero_props%get(n,l, spectype=spectype) - if (trim(spectype) == 'sulfate') then - qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) - end if - if (trim(spectype) == 'msa') then - qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) - end if - if (trim(spectype) == 'ammonium') then - qcw(:ncol,k,mm) = MAX(qcw(:ncol,k,mm), small_value ) - end if - end do - end do - - qin(:ncol,k,id_so2) = MAX( qin(:ncol,k,id_so2), small_value ) - qin(:ncol,k,id_h2o2) = MAX( qin(:ncol,k,id_h2o2), small_value ) - qin(:ncol,k,id_h2so4) = MAX( qin(:ncol,k,id_h2so4), small_value ) - if ( id_msa > 0 ) qin(:ncol,k,id_msa) = MAX( qin(:ncol,k,id_msa), small_value ) - if ( id_nh3 > 0 ) qin(:ncol,k,id_nh3) = MAX( qin(:ncol,k,id_nh3), small_value ) - - end do - - ! diagnostics - mw_so4 = -huge(1._r8) - - do n = 1, aero_props%nbins() - ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero - do l = 1, aero_props%nspecies(n) - mm = aero_props%indexer(n,l) - call aero_props%get(n,l, spectype=spectype, specname=specname) - if (trim(spectype) == 'sulfate') then - call aero_props%get(n,l, spec_mw=mw_so4) - aqso4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s - enddo - enddo - aqh2so4(:,n)=0._r8 - do k=1,pver - do i=1,ncol - aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg/m2/s - enddo - enddo - end if - end do - end do - - aqso4_h2o2(:) = 0._r8 - do k=1,pver - do i=1,ncol - aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo - - if (present(aqso4_h2o2_3d)) then - aqso4_h2o2_3d(:,:) = 0._r8 - do k=1,pver - do i=1,ncol - aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo - end if - - aqso4_o3(:)=0._r8 - do k=1,pver - do i=1,ncol - aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo - - if (present(aqso4_o3_3d)) then - aqso4_o3_3d(:,:)=0._r8 - do k=1,pver - do i=1,ncol - aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*mw_so4/mbar(i,k) & - *pdel(i,k)/gravit ! kg SO4 /m2/s - enddo - enddo - end if - - end subroutine sox_cldaero_update - - !---------------------------------------------------------------------------------- - !---------------------------------------------------------------------------------- - subroutine sox_cldaero_destroy_obj( conc_obj ) - type(cldaero_conc_t), pointer :: conc_obj - - call cldaero_deallocate( conc_obj ) - - end subroutine sox_cldaero_destroy_obj - -end module sox_cldaero_mod From 4cc10c8cfe753869f5225e457cd6b2471e354682 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 7 May 2026 10:51:02 -0600 Subject: [PATCH 6/9] more code clean up modified: src/chemistry/aerosol/mo_setsox.F90 modified: src/chemistry/aerosol/sox_cldaero_mod.F90 modified: src/chemistry/bulk_aero/aero_model.F90 modified: src/chemistry/carma_aero/aero_model.F90 modified: src/chemistry/modal_aero/aero_model.F90 --- src/chemistry/aerosol/mo_setsox.F90 | 25 +++++--------- src/chemistry/aerosol/sox_cldaero_mod.F90 | 15 +++------ src/chemistry/bulk_aero/aero_model.F90 | 2 -- src/chemistry/carma_aero/aero_model.F90 | 2 -- src/chemistry/modal_aero/aero_model.F90 | 40 +++++++++++++---------- 5 files changed, 35 insertions(+), 49 deletions(-) diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index a32b4afaea..6300bc315f 100644 --- a/src/chemistry/aerosol/mo_setsox.F90 +++ b/src/chemistry/aerosol/mo_setsox.F90 @@ -18,7 +18,7 @@ module mo_setsox integer :: id_so4, id_h2so4 logical :: has_sox = .true. - logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2 + logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ho2 logical :: cloud_borne = .false. @@ -155,8 +155,6 @@ end subroutine sox_inti subroutine setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset,& dtime, & press, & pdel, & @@ -202,7 +200,6 @@ subroutine setsox( aero_state, state, & MOLECULAR_WEIGHT_DRY_AIR_G_MOL => mwdry use ppgrid, only : pcols, pver use chem_mods, only : gas_pcnst, nfs - use chem_mods, only : adv_mass use physconst, only : mwdry, gravit use mo_constants, only : pi use sox_cldaero_mod, only : sox_cldaero_update, sox_cldaero_create_obj, sox_cldaero_destroy_obj @@ -219,8 +216,6 @@ subroutine setsox( aero_state, state, & type(physics_state), intent(in) :: state ! Physics state variables type(physics_buffer_desc), pointer, intent(inout) :: pbuf(:) ! Physics buffer integer, intent(in) :: ncol ! num of columns in chunk - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array real(r8), intent(in) :: dtime ! time step (sec) real(r8), intent(in) :: press(:,:) ! midpoint pressure ( Pa ) real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa) @@ -279,15 +274,13 @@ subroutine setsox( aero_state, state, & real(r8) :: xdelso4hp(ncol,pver) real(r8) :: xhnm(ncol,pver) ! air number density (molecules cm-3) - integer :: k, i, iter, file + integer :: k, i, iter real(r8) :: wrk, delta - real(r8) :: xph0, aden, xk, xe, x2 - real(r8) :: tz, xl, px, qz, es, qs, patm + real(r8) :: xph0, xk, xe, x2 + real(r8) :: tz, xl, px, patm real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3 real(r8) :: so2g, h2o2g, o3g - real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a real(r8) :: rah2o2, rao3, pso4, ccc - real(r8) :: cnh3, chno3, com, com1, com2, xra real(r8) :: f_hso3 ! fraction of aqueous S(IV) that's HSO3- real(r8) :: f_so3 ! fraction of aqueous S(IV) that's SO3= @@ -360,7 +353,7 @@ subroutine setsox( aero_state, state, & * 1.e-3_r8 ! Kg(a)/L(a) end do - cldconc => sox_cldaero_create_obj( cldfrc,qcw,lwc, cfact, ncol, loffset ) + cldconc => sox_cldaero_create_obj( cldfrc,qcw,lwc, cfact, ncol ) xso4c => cldconc%so4c xnh4c => cldconc%nh4c xno3c => cldconc%no3c @@ -688,8 +681,8 @@ subroutine setsox( aero_state, state, & end do ! iter if( .not. converged ) then - write(iulog,*) 'setsox: pH failed to converge @ (',i,',',k,'), % change=', & - 100._r8*delta + write(*,*) 'setsox: pH failed to converge @ (',i,',',k,'), % change=', & + 100._r8*delta !!! What should delta be set to ???? end if else xph(i,k) = 1.e-7_r8 @@ -892,8 +885,8 @@ subroutine setsox( aero_state, state, & if (cloud_borne) then ! update cloud-borne aerosols call sox_cldaero_update( aero_state, & - state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & - xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & + ncol, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & + xdelso4hp, xh2so4, xso4, xso4_init, nh3g, xnh3, xnh4c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d ) else if (id_so2>0) then diff --git a/src/chemistry/aerosol/sox_cldaero_mod.F90 b/src/chemistry/aerosol/sox_cldaero_mod.F90 index bd8fbf32d9..e42c341b0f 100644 --- a/src/chemistry/aerosol/sox_cldaero_mod.F90 +++ b/src/chemistry/aerosol/sox_cldaero_mod.F90 @@ -62,14 +62,13 @@ end subroutine sox_cldaero_init !---------------------------------------------------------------------------------- !---------------------------------------------------------------------------------- - function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj ) + function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol) result( conc_obj ) real(r8), intent(in) :: cldfrc(:,:) real(r8), intent(in) :: qcw(:,:,:) real(r8), intent(in) :: lwc(:,:) real(r8), intent(in) :: cfact(:,:) integer, intent(in) :: ncol - integer, intent(in) :: loffset type(cldaero_conc_t), pointer :: conc_obj @@ -140,8 +139,8 @@ end function sox_cldaero_create_obj ! Update the mixing ratios !---------------------------------------------------------------------------------- subroutine sox_cldaero_update( aero_state, & - state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & - delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & + ncol, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, & + delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, xnh3, xnh4c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d) use physics_types, only: physics_state @@ -149,11 +148,8 @@ subroutine sox_cldaero_update( aero_state, & ! args class(aerosol_state), intent(in) :: aero_state - type(physics_state), intent(in) :: state ! Physics state variables integer, intent(in) :: ncol - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset real(r8), intent(in) :: dtime ! time step (sec) @@ -172,14 +168,11 @@ subroutine sox_cldaero_update( aero_state, & real(r8), intent(in) :: xso4(:,:) real(r8), intent(in) :: xso4_init(:,:) real(r8), intent(in) :: nh3g(:,:) - real(r8), intent(in) :: hno3g(:,:) real(r8), intent(in) :: xnh3(:,:) - real(r8), intent(in) :: xhno3(:,:) real(r8), intent(in) :: xnh4c(:,:) real(r8), intent(in) :: xmsa(:,:) real(r8), intent(in) :: xso2(:,:) real(r8), intent(in) :: xh2o2(:,:) - real(r8), intent(in) :: xno3c(:,:) real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr ) @@ -210,7 +203,7 @@ subroutine sox_cldaero_update( aero_state, & real(r8) :: fwetrem, uptkrate integer :: l, m, n, mm - integer :: i,k, ndx + integer :: i,k real(r8) :: xl real(r8) :: mw_so4 character(len=32) :: spectype diff --git a/src/chemistry/bulk_aero/aero_model.F90 b/src/chemistry/bulk_aero/aero_model.F90 index e1d81ba3d0..cd96f61d21 100644 --- a/src/chemistry/bulk_aero/aero_model.F90 +++ b/src/chemistry/bulk_aero/aero_model.F90 @@ -1087,8 +1087,6 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re call setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset, & delt, & pmid, & pdel, & diff --git a/src/chemistry/carma_aero/aero_model.F90 b/src/chemistry/carma_aero/aero_model.F90 index a1952acb9d..dd0c96a8ac 100644 --- a/src/chemistry/carma_aero/aero_model.F90 +++ b/src/chemistry/carma_aero/aero_model.F90 @@ -655,8 +655,6 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re call setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset, & delt, & pmid, & pdel, & diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index a039d7ff71..a091798087 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -106,6 +106,7 @@ module aero_model integer :: n_coarse_dust=-1 ! dmleung added n_coarse_dust to determine the index for the ! coarse dust mode for different MAM versions. 29 Oct 2025 integer :: ncnst_tot = -1 + integer :: chem_map_ndx(gas_pcnst) = -1 contains @@ -218,6 +219,8 @@ subroutine aero_model_init( pbuf2d ) character(len=32) :: spec_type character(len=32) :: mode_type integer :: nspec + integer :: mm + character(len=32) :: name_a, name_c aero_props => modal_aerosol_properties() ncnst_tot = aero_props%ncnst_tot() @@ -225,6 +228,18 @@ subroutine aero_model_init( pbuf2d ) ! aqueous chem initialization call sox_inti(aero_props) + do m = 1,aero_props%nbins() + do l = 0,aero_props%nspecies(m) + mm = aero_props%indexer(m,l) + if (l==0) then + call aero_props%num_names(m, name_a, name_c) + else + call aero_props%mmr_names(m,l, name_a, name_c) + end if + chem_map_ndx(mm) = get_spc_ndx( name_a ) + end do + end do + dgnum_idx = pbuf_get_index('DGNUM') dgnumwet_idx = pbuf_get_index('DGNUMWET') fracis_idx = pbuf_get_index('FRACIS') @@ -1048,9 +1063,8 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8) :: qqcw(ncol,pver,ncnst_tot) - integer :: ndx, mm + integer :: mm character(len=32) :: specname - character(len=32) :: name_a, name_c class(aerosol_state), pointer :: aero_state aero_state => modal_aerosol_state(state, pbuf) @@ -1099,24 +1113,19 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re if( has_sox ) then + ! Temperary code to map cloud-borne aerosol VMRs to aerosol only array (qqcw) + ! needed for setsox interface. When refactoring aero_model_gasaerexch + ! with modal_aero_gasaerexch_sub, this mapping should go away. do m = 1,aero_props%nbins() do l = 0,aero_props%nspecies(m) mm = aero_props%indexer(m,l) - if (l==0) then - call aero_props%num_names(m, name_a, name_c) - else - call aero_props%mmr_names(m,l, name_a, name_c) - end if - ndx = get_spc_ndx( name_a ) - qqcw(:,:,mm) = vmrcw(:,:,ndx) + qqcw(:,:,mm) = vmrcw(:,:,chem_map_ndx(mm)) end do end do call setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset, & delt, & pmid, & pdel, & @@ -1135,16 +1144,11 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re aqso4_o3 & ) + ! Map back to all-species chemistry VMR array do m = 1,aero_props%nbins() do l = 0,aero_props%nspecies(m) mm = aero_props%indexer(m,l) - if (l==0) then - call aero_props%num_names(m, name_a, name_c) - else - call aero_props%mmr_names(m,l, name_a, name_c) - end if - ndx = get_spc_ndx( name_a ) - vmrcw(:,:,ndx) = qqcw(:,:,mm) + vmrcw(:,:,chem_map_ndx(mm)) = qqcw(:,:,mm) end do end do From ed83364bf3c3527a336d72feee88d56312a95b02 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 7 May 2026 17:10:51 -0600 Subject: [PATCH 7/9] fix dangling pointer problem with nag compiler modified: src/chemistry/aerosol/mo_setsox.F90 --- src/chemistry/aerosol/mo_setsox.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index 6300bc315f..1b1a83fa63 100644 --- a/src/chemistry/aerosol/mo_setsox.F90 +++ b/src/chemistry/aerosol/mo_setsox.F90 @@ -41,7 +41,7 @@ subroutine sox_inti(aero_props) use sox_cldaero_mod, only : sox_cldaero_init use aerosol_properties_mod, only : aerosol_properties - class(aerosol_properties), intent(in) :: aero_props + class(aerosol_properties), target, intent(in) :: aero_props logical :: modal_aerosols From 131168cc157789963898f01479a30a455afc9d82 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Thu, 14 May 2026 14:22:00 -0600 Subject: [PATCH 8/9] add comments; minor code clean modified: src/chemistry/aerosol/aerosol_state_mod.F90 modified: src/chemistry/aerosol/carma_aerosol_state_mod.F90 modified: src/chemistry/aerosol/modal_aerosol_state_mod.F90 modified: src/chemistry/aerosol/sox_cldaero_mod.F90 --- src/chemistry/aerosol/aerosol_state_mod.F90 | 10 +++++----- .../aerosol/carma_aerosol_state_mod.F90 | 10 +++++----- .../aerosol/modal_aerosol_state_mod.F90 | 10 +++++----- src/chemistry/aerosol/sox_cldaero_mod.F90 | 17 +++++++---------- 4 files changed, 22 insertions(+), 25 deletions(-) diff --git a/src/chemistry/aerosol/aerosol_state_mod.F90 b/src/chemistry/aerosol/aerosol_state_mod.F90 index c0fa318c01..c1bc31f73a 100644 --- a/src/chemistry/aerosol/aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/aerosol_state_mod.F90 @@ -288,11 +288,11 @@ subroutine aero_aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, import :: aerosol_state, aerosol_properties, r8 class(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 + 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 diff --git a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 index f0119ee007..d8b197bd17 100644 --- a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 @@ -612,11 +612,11 @@ end function wet_diameter 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 - 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 + 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) diff --git a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 index e1601bd82f..daab068375 100644 --- a/src/chemistry/aerosol/modal_aerosol_state_mod.F90 +++ b/src/chemistry/aerosol/modal_aerosol_state_mod.F90 @@ -724,11 +724,11 @@ end function wgtpct subroutine aqu_gain_binfraction(self, aero_props, type, qcw, delso4_o3rxn, faqgain) class(modal_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 + 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 character(len=aero_name_len) :: modetype, spectype integer :: i,k,l,m,n,mm, ncol, nbins diff --git a/src/chemistry/aerosol/sox_cldaero_mod.F90 b/src/chemistry/aerosol/sox_cldaero_mod.F90 index e42c341b0f..dcc6a64374 100644 --- a/src/chemistry/aerosol/sox_cldaero_mod.F90 +++ b/src/chemistry/aerosol/sox_cldaero_mod.F90 @@ -191,11 +191,10 @@ subroutine sox_cldaero_update( aero_state, & dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver) real(r8) :: faqgain_msa(nbins,ncol,pver), faqgain_so4(nbins,ncol,pver) - real(r8) :: delso4_3d(ncol,pver) + real(r8) :: delso4_ox(ncol,pver) real(r8) :: delnh3, delnh4 - real(r8) :: delso4_o3rxn, & - dso4dt_aqrxn, dso4dt_hprxn, & + real(r8) :: dso4dt_aqrxn, dso4dt_hprxn, & dso4dt_gasuptk, dmsadt_gasuptk, & dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, & dqdt_aq, dqdt_wr, dqdt @@ -219,7 +218,7 @@ subroutine sox_cldaero_update( aero_state, & aqh2so4 = 0.0_r8 aqso4_h2o2 = 0.0_r8 aqso4_o3 = 0.0_r8 - delso4_3d = 0.0_r8 + delso4_ox = 0.0_r8 ! Avoid double counting in-cloud sulfur oxidation when running with ! GEOS-Chem. If running with GEOS-Chem then sulfur oxidation @@ -228,7 +227,7 @@ subroutine sox_cldaero_update( aero_state, & if ( cam_chempkg_is('geoschem_mam4') ) return where (cldfrc(:ncol,:) >= 1.0e-5_r8) - delso4_3d(:ncol,:) = xso4(:ncol,:) - xso4_init(:ncol,:) + delso4_ox(:ncol,:) = xso4(:ncol,:) - xso4_init(:ncol,:) end where !------------------------------------------------------------------------- @@ -237,8 +236,8 @@ subroutine sox_cldaero_update( aero_state, & ! bin, which is the MR of cloud drops "associated with" the mode ! thus we are assuming the cloud drop size is independent of the ! associated aerosol mode properties - call aero_state%aqu_gain_binfraction(aero_props, 'sulfate', qcw, delso4_3d, faqgain_so4) - if (has_msa) call aero_state%aqu_gain_binfraction(aero_props, 'msa', qcw, delso4_3d, faqgain_msa) + call aero_state%aqu_gain_binfraction(aero_props, 'sulfate', qcw, delso4_ox, faqgain_so4) + if (has_msa) call aero_state%aqu_gain_binfraction(aero_props, 'msa', qcw, delso4_ox, faqgain_msa) lev_loop: do k = 1,pver col_loop: do i = 1,ncol @@ -247,8 +246,6 @@ subroutine sox_cldaero_update( aero_state, & if (xl .ge. 1.e-8_r8) then !! when cloud is present - delso4_o3rxn = delso4_3d(i,k) ! xso4(i,k) - xso4_init(i,k) - if (id_nh3>0) then delnh3 = nh3g(i,k) - xnh3(i,k) delnh4 = - delnh3 @@ -285,7 +282,7 @@ subroutine sox_cldaero_update( aero_state, & ! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...) ! AND the wetremoval of dissolved, unreacted so2 and h2o2 - dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime + dso4dt_aqrxn = (delso4_ox(i,k) + delso4_hprxn(i,k)) / dtime dso4dt_hprxn = delso4_hprxn(i,k) / dtime ! fwetrem = fraction of in-cloud-water material that is wet removed From 7aee22d72bbb9db71842e9a7185b9b340b743e89 Mon Sep 17 00:00:00 2001 From: Francis Vitt Date: Fri, 15 May 2026 08:16:23 -0600 Subject: [PATCH 9/9] fix geos-chem build modified: src/chemistry/geoschem/chemistry.F90 --- src/chemistry/geoschem/chemistry.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/chemistry/geoschem/chemistry.F90 b/src/chemistry/geoschem/chemistry.F90 index adf71a332a..fd231e78be 100644 --- a/src/chemistry/geoschem/chemistry.F90 +++ b/src/chemistry/geoschem/chemistry.F90 @@ -21,6 +21,7 @@ module chemistry use string_utils, only : to_upper #if defined( MODAL_AERO ) use modal_aero_data, only : ntot_amode + use modal_aerosol_properties_mod, only: modal_aerosol_properties #endif ! GEOS-Chem derived types @@ -169,6 +170,10 @@ module chemistry ! for nitrogen deposition fluxes to surface models logical, parameter :: chem_has_ndep_flx = .false. +#if defined( MODAL_AERO ) + ! MAM aerosol properties object class + type(modal_aerosol_properties), pointer :: aero_props=>null() +#endif contains !================================================================================================ @@ -1573,8 +1578,11 @@ subroutine chem_init(phys_state, pbuf2d) ENDIF #if defined( MODAL_AERO ) + ! instantiate MAM aerosol properties object + aero_props => modal_aerosol_properties() + ! Initialize aqueous chem - CALL SOx_inti() + CALL SOx_inti(aero_props) ! Initialize aerosols CALL aero_model_init( pbuf2d )