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..db2f740226 100644 --- a/src/chemistry/aerosol/aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/aerosol_properties_mod.F90 @@ -109,7 +109,7 @@ end function aero_number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, & + subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, spec_mw, & spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) import :: aerosol_properties, r8 @@ -118,6 +118,7 @@ subroutine aero_props_get(self, bin_ndx, species_ndx, density, hygro, & integer, intent(in) :: species_ndx ! species index real(r8), optional, intent(out) :: density ! density (kg/m3) real(r8), optional, intent(out) :: hygro ! hygroscopicity + real(r8), optional, intent(out) :: spec_mw ! species molecular weight character(len=*), optional, intent(out) :: spectype ! species type character(len=*), optional, intent(out) :: specname ! species name character(len=*), optional, intent(out) :: specmorph ! species morphology diff --git a/src/chemistry/aerosol/aerosol_state_mod.F90 b/src/chemistry/aerosol/aerosol_state_mod.F90 index 7fb1689b1a..c1bc31f73a 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 ! aerosol properties object + character(len=*), intent(in) :: type ! aerosol species type + real(r8), intent(in) :: qcw(:,:,:) ! cloud-borne aerosol volume mixing ratio + real(r8), intent(in) :: delso4_o3rxn(:,:) ! sulfate concentration change due to oxidation + real(r8), intent(out) :: faqgain(:,:,:) ! fraction gain in each mode / bin + + end subroutine aero_aqu_gain_binfraction + end interface contains diff --git a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 index 1ca5f3d801..4b5e3552c7 100644 --- a/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/bulk_aerosol_properties_mod.F90 @@ -152,7 +152,7 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & + subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, & spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) @@ -161,6 +161,7 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & integer, intent(in) :: species_ndx ! species index real(r8), optional, intent(out) :: density ! density (kg/m3) real(r8), optional, intent(out) :: hygro ! hygroscopicity + real(r8), optional, intent(out) :: spec_mw ! species molecular weight character(len=*), optional, intent(out) :: spectype ! species type character(len=*), optional, intent(out) :: specname ! species name character(len=*), optional, intent(out) :: specmorph ! species morphology @@ -217,6 +218,23 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & if (present(dryrad)) then call rad_aer_get_props(self%list_idx_, bin_ndx, dryrad_aer=dryrad) end if + if (present(spec_mw)) then + call rad_aer_get_props(self%list_idx_, bin_ndx, aername=aername) + + select case ( to_lower( aername(:4) ) ) + case('sulf','volc') + spec_mw = 96._r8 + case('bcar','bcph','ocar','ocph') + spec_mw = 12._r8 + case('dust') + spec_mw = 12._r8 !!! ???? + case('sslt','seas','ssam','sscm') + spec_mw = 57._r8 + case default + spec_mw = nan + call endrun('ERROR: bulk_aerosol_properties_mod%get aername not recognized : '//aername) + end select + end if end subroutine get 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..718adfda60 100644 --- a/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 +++ b/src/chemistry/aerosol/carma_aerosol_properties_mod.F90 @@ -264,7 +264,7 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & + subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, & spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) use cam_abortutils, only: endrun @@ -274,6 +274,7 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & integer, intent(in) :: species_ndx ! species index real(r8), optional, intent(out) :: density ! density (kg/m3) real(r8), optional, intent(out) :: hygro ! hygroscopicity + real(r8), optional, intent(out) :: spec_mw ! species molecular weight character(len=*), optional, intent(out) :: spectype ! species type character(len=*), optional, intent(out) :: specname ! species name character(len=*), optional, intent(out) :: specmorph ! species morphology @@ -282,6 +283,8 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & real(r8), optional, intent(out) :: num_to_mass_aer ! ratio of number to mass concentration real(r8), optional, intent(out) :: dryrad ! dry radius (m) + character(len=32) :: type + if (present(density)) then call rad_aer_get_bin_props_by_idx(self%list_idx_, bin_ndx, species_ndx, density_aer=density) end if @@ -307,6 +310,25 @@ subroutine get(self, bin_ndx, species_ndx, density, hygro, & call rad_aer_get_info_by_bin_spec(self%list_idx_, bin_ndx, species_ndx, spec_name=specname) end if end if + if (present(spec_mw)) then + call rad_aer_get_bin_props_by_idx(self%list_idx_, bin_ndx, species_ndx, spectype=type) + select case (trim(type)) + case('sulfate') + spec_mw = 96._r8 + case('black-c') + spec_mw = 12._r8 + case('s-organic') + spec_mw = 250._r8 + case('p-organic') + spec_mw = 12._r8 + case('dust') + spec_mw = 12._r8 + case('seasalt') + spec_mw = 57._r8 + case default + spec_mw = nan + end select + end if if (present(num_to_mass_aer)) then ! num_to_mass_aer for sectional aerosols should not be read from file diff --git a/src/chemistry/aerosol/carma_aerosol_state_mod.F90 b/src/chemistry/aerosol/carma_aerosol_state_mod.F90 index e2b52f14ba..d8b197bd17 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 ! aerosol properties object + character(len=*), intent(in) :: type ! aerosol species type + real(r8), intent(in) :: qcw(:,:,:) ! cloud-borne aerosol volume mixing ratio + real(r8), intent(in) :: delso4_o3rxn(:,:) ! sulfate concentration change due to oxidation + real(r8), intent(out) :: faqgain(:,:,:) ! fraction gain in each mode / bin + + real(r8) :: raddry(pcols,pver) ! dry radius (m) + real(r8) :: rhodry(pcols,pver) ! dry density (kg/m3) + character(len=aero_name_len) :: bin_name, shortname + integer :: igroup, ibin, jbin, rc, nchr + integer :: icol, klev, nbins, ncol + real(r8), allocatable :: rad_cm(:,:,:) + real(r8), allocatable :: wt_mass(:) + real(r8) :: wt_sum + + ! To calculate the fraction of sulfate mass produced by aq.chemistry that will be + ! distributed across different bins, we need to conserve particle number. In CARMA, + ! this is done following the formula for mass transfer for, following + ! Yu et al., 2015, A3.2: https://doi.org/10.1002/2014MS000421. The mass fraction + ! for each bin can be calculated by normalizing the mass transfer rate in each bin: + ! M/D^2 as: fra = Mi/D^2/sum(Mi/D^2). + + ncol = self%state%ncol + nbins = aero_props%nbins() + faqgain(:,:,:) = 0._r8 + + allocate(wt_mass(nbins)) + allocate(rad_cm(nbins,pcols,pver)) + rad_cm(:,:,:) = 0._r8 + + do ibin = 1, nbins + call rad_aer_get_info_by_bin(0, ibin, bin_name=bin_name) + nchr = len_trim(bin_name)-2 + shortname = bin_name(:nchr) + call carma_get_group_by_name(shortname, igroup, rc) + read(bin_name(nchr+1:),*) jbin + call carma_get_dry_radius(self%state, igroup, jbin, raddry, rhodry, rc) + if (index(bin_name,'MXAER')>0) then + rad_cm(ibin,:ncol,:) = raddry(:ncol,:)*1.e2_r8 ! m -> cm + end if + end do + + lev_loop: do klev = 1,pver + col_loop: do icol = 1,ncol + + !faqgain = fraction of total so4_c gain going to bin n + wt_sum = 0._r8 + wt_mass(:) = 0._r8 + + do ibin = 1, nbins + if (rad_cm(ibin,icol,klev) > 0._r8) then + wt_mass(ibin) = delso4_o3rxn(icol,klev) / rad_cm(ibin,icol,klev) / rad_cm(ibin,icol,klev) + wt_sum = wt_sum + wt_mass(ibin) + end if + end do + do ibin = 1, nbins + if (wt_mass(ibin) > 0._r8) then + faqgain(ibin,icol,klev) = wt_mass(ibin)/wt_sum + end if + end do + + end do col_loop + end do lev_loop + + deallocate(rad_cm, wt_mass) + + end subroutine aqu_gain_binfraction + end module carma_aerosol_state_mod diff --git a/src/chemistry/aerosol/mo_setsox.F90 b/src/chemistry/aerosol/mo_setsox.F90 index 7e0077adc6..1b1a83fa63 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 @@ -17,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. @@ -28,7 +29,7 @@ module mo_setsox !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine sox_inti + subroutine sox_inti(aero_props) !----------------------------------------------------------------------- ! ... initialize the hetero sox routine !----------------------------------------------------------------------- @@ -38,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), target, intent(in) :: aero_props logical :: modal_aerosols @@ -142,17 +146,15 @@ subroutine sox_inti return end if - call sox_cldaero_init() + call sox_cldaero_init(aero_props) end subroutine sox_inti !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine setsox( state, & + subroutine setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset,& dtime, & press, & pdel, & @@ -198,7 +200,6 @@ subroutine setsox( 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 @@ -211,11 +212,10 @@ 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 - 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) @@ -268,20 +268,19 @@ subroutine setsox( 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) 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= @@ -354,7 +353,7 @@ subroutine setsox( 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 @@ -362,6 +361,7 @@ 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 @@ -681,8 +681,8 @@ subroutine setsox( 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 @@ -757,9 +757,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,10 +877,26 @@ subroutine setsox( state, & end do col_loop1 end do ver_loop1 - call sox_cldaero_update( & - 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, & + 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 + 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 @@ -900,7 +916,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..191e62118f 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 @@ -407,7 +407,7 @@ end function number_transported ! long wave species refractive indices ! species morphology !------------------------------------------------------------------------ - subroutine get(self, bin_ndx, species_ndx, density, hygro, & + subroutine get(self, bin_ndx, species_ndx, density, hygro, spec_mw, & spectype, specname, specmorph, refindex_sw, refindex_lw, num_to_mass_aer, & dryrad) use cam_abortutils, only: endrun @@ -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..daab068375 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 ! 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 + 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/sox_cldaero_mod.F90 b/src/chemistry/aerosol/sox_cldaero_mod.F90 similarity index 52% rename from src/chemistry/carma_aero/sox_cldaero_mod.F90 rename to src/chemistry/aerosol/sox_cldaero_mod.F90 index 319c177b26..dcc6a64374 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 @@ -8,13 +8,12 @@ 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 aerosol_properties_mod, only: aerosol_properties + use aerosol_state_mod, only: aerosol_state implicit none private @@ -24,110 +23,79 @@ 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 + class(aerosol_properties), pointer :: aero_props =>null() + + logical :: has_msa = .false. contains !---------------------------------------------------------------------------------- !---------------------------------------------------------------------------------- - subroutine sox_cldaero_init - - integer :: l, m, ii - logical :: history_aerosol ! Output the MAM aerosol tendencies + 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 ') + if ( id_so2<1 ) then + call endrun('sox_cldaero_init: SO2 is not included in chemistry -- 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 => 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 ) + 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 - - real(r8) :: so4mmr(pcols,pver) 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 - ! local indexing for bins - !integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr + 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 - conc_obj => cldaero_allocate() + if (aero_props%model_is('MAM')) then + ntot_amode = aero_props%nbins() + else + ntot_amode = 0 + end if + + mode7 = ntot_amode == 7 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 @@ -140,44 +108,48 @@ 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( & - 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, & + subroutine sox_cldaero_update( aero_state, & + 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 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 - type(physics_state), intent(in) :: state ! Physics state variables + class(aerosol_state), intent(in) :: aero_state integer, intent(in) :: ncol - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset real(r8), intent(in) :: dtime ! time step (sec) @@ -196,16 +168,13 @@ subroutine sox_cldaero_update( & 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) 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,76 +185,59 @@ 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_ox(ncol,pver) - real(r8) :: delso4_o3rxn, & - dso4dt_aqrxn, dso4dt_hprxn, & - dso4dt_gasuptk, dmsadt_gasuptk_toso4, & + real(r8) :: delnh3, delnh4 + real(r8) :: dso4dt_aqrxn, dso4dt_hprxn, & + dso4dt_gasuptk, dmsadt_gasuptk, & + dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, & dqdt_aq, dqdt_wr, dqdt - real(r8) :: delnh3 real(r8) :: fwetrem, uptkrate - integer :: l, n, mm - integer :: ntot_msa_c - + integer :: l, m, n, mm integer :: i,k 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_ox = 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_ox(: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_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 @@ -294,45 +246,35 @@ 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) - if (id_nh3>0) then delnh3 = nh3g(i,k) - xnh3(i,k) + delnh4 = - delnh3 endif - ! 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 + ! 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 @@ -340,36 +282,49 @@ subroutine sox_cldaero_update( & ! 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 - !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 @@ -380,6 +335,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)))) ) @@ -410,7 +366,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 present + endif !! when cloud is present endif cloud enddo col_loop enddo lev_loop @@ -420,44 +376,51 @@ subroutine sox_cldaero_update( & !============================================================== do k = 1,pver - 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) + 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(:,k,mm) = max(qcw(:,k,mm), small_value ) + 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_nh3 > 0 ) qin(:ncol,k,id_nh3) = max( qin(:ncol,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 + mw_so4 = -huge(1._r8) - specmw_so4_amode = 96.0_r8 - do n = 1, nbins + do n = 1, aero_props%nbins() ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero - do l = 1, nspec(n) - call rad_aer_get_bin_props_by_idx(0, n, l,spectype=spectype) + 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 - mm = bin_idx(n, l) + 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)*specmw_so4_amode/mbar(i,k) & + 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)*specmw_so4_amode/mbar(i,k) & + 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 @@ -468,37 +431,37 @@ subroutine sox_cldaero_update( & 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/bulk_aero/aero_model.F90 b/src/chemistry/bulk_aero/aero_model.F90 index 798a59e7f3..cd96f61d21 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,15 +1076,17 @@ 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, & - loffset, & delt, & pmid, & pdel, & 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 0c5a7cc923..0000000000 --- a/src/chemistry/bulk_aero/sox_cldaero_mod.F90 +++ /dev/null @@ -1,143 +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 - - 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 - - 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( & - 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 - - 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/carma_aero/aero_model.F90 b/src/chemistry/carma_aero/aero_model.F90 index 778e8093f6..dd0c96a8ac 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 @@ -251,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') @@ -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,11 +652,9 @@ 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, & - loffset, & delt, & pmid, & pdel, & @@ -705,10 +675,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 +725,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 +890,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 +926,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 8c2210c1f0..bde88507ce 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 @@ -693,7 +680,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 @@ -712,6 +699,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/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 ) diff --git a/src/chemistry/modal_aero/aero_model.F90 b/src/chemistry/modal_aero/aero_model.F90 index 6e450ab2b4..a091798087 100644 --- a/src/chemistry/modal_aero/aero_model.F90 +++ b/src/chemistry/modal_aero/aero_model.F90 @@ -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 implicit none private @@ -103,6 +105,8 @@ 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 @@ -215,9 +219,26 @@ 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() ! aqueous chem initialization - call sox_inti() + 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') @@ -250,7 +271,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 @@ -985,11 +1005,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), 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 @@ -1040,6 +1061,13 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re real(r8), pointer :: fldcw(:,:) real(r8), pointer :: sulfeq(:,:,:) + real(r8) :: qqcw(ncol,pver,ncnst_tot) + + integer :: mm + character(len=32) :: specname + class(aerosol_state), pointer :: aero_state + + aero_state => modal_aerosol_state(state, pbuf) ! ! ... initialize nh3 ! @@ -1084,11 +1112,20 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re ! aqueous chemistry ... if( has_sox ) then - call setsox( state, & + + ! 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) + qqcw(:,:,mm) = vmrcw(:,:,chem_map_ndx(mm)) + end do + end do + + call setsox( aero_state, state, & pbuf, & ncol, & - lchnk, & - loffset, & delt, & pmid, & pdel, & @@ -1098,7 +1135,7 @@ subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, re cldfr, & cldnum, & invariants, & - vmrcw, & + qqcw, & vmr, & xphlwc, & aqso4, & @@ -1107,6 +1144,14 @@ 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) + vmrcw(:,:,chem_map_ndx(mm)) = qqcw(:,:,mm) + end do + end do + do n = 1, ntot_amode l = lptr_so4_cw_amode(n) if (l > 0) then 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 1765e51c34..0000000000 --- a/src/chemistry/modal_aero/sox_cldaero_mod.F90 +++ /dev/null @@ -1,553 +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 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 cldaero_mod, only : cldaero_uptakerate - use chem_mods, only : gas_pcnst - - implicit none - private - - public :: sox_cldaero_init - public :: sox_cldaero_create_obj - public :: sox_cldaero_update - public :: sox_cldaero_destroy_obj - - integer :: id_msa, id_h2so4, id_so2, id_h2o2, id_nh3 - - real(r8), parameter :: small_value = 1.e-20_r8 - -contains - -!---------------------------------------------------------------------------------- -!---------------------------------------------------------------------------------- - - 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' ) - - 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 - ! - - 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 - - - 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 - - logical :: mode7 - logical :: mode5 - - mode7 = ntot_amode == 7 - mode5 = ntot_amode == 5 - - 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 - - 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 if (mode5) then -#if ( defined MODAL_AERO_5MODE ) - 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 - id_so4_5a = lptr_so4_cw_amode(5) - loffset -#endif - conc_obj%so4c(:ncol,:) & - = qcw(:,:,id_so4_1a) & - + qcw(:,:,id_so4_2a) & - + qcw(:,:,id_so4_3a) & - + qcw(:,:,id_so4_5a) - - ! 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 - - 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 - - endif - - end function sox_cldaero_create_obj - -!---------------------------------------------------------------------------------- -! Update the mixing ratios -!---------------------------------------------------------------------------------- - subroutine sox_cldaero_update( & - 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 - - 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,gas_pcnst), & - dqdt_aqh2so4(ncol,pver,gas_pcnst), & - dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), & - sflx(1:ncol) - - real(r8) :: faqgain_msa(ntot_amode), faqgain_so4(ntot_amode), qnum_c(ntot_amode) - - 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 - - integer :: i,k - real(r8) :: xl - - ! 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 - - ! 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 - - 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) - - IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED - - delso4_o3rxn = 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 - 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 (id_msa > 0) 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 - 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 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 - 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 (id_msa > 0) 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 PRESENTED - endif cloud - enddo col_loop - enddo lev_loop - - !============================================================== - ! ... Update the mixing ratios - !============================================================== - 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 - - 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 ) - - 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 - 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 - enddo - enddo - endif - 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 - 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 - 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 - 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 - 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