From d0eec18c6d45163ca15269144a7f4160a03a5e40 Mon Sep 17 00:00:00 2001 From: Xiaotian Xu Date: Thu, 16 Jul 2026 14:53:50 -0500 Subject: [PATCH 1/5] Add effective surface tension --- scenarios/1_urban_plume/aero_data.dat | 42 +-- .../1_urban_plume/urban_plume_process.F90 | 20 +- scenarios/2_urban_plume2/aero_data.dat | 42 +-- scenarios/3_condense/aero_data.dat | 42 +-- scenarios/4_chamber/aero_data.dat | 42 +-- scenarios/5_coag_brownian/aero_data.dat | 42 +-- src/aero_data.F90 | 22 +- src/aero_particle.F90 | 300 +++++++++++++++++- src/aero_state.F90 | 27 ++ test/additive/aero_data.dat | 4 +- test/average/aero_data.dat | 42 +-- test/bidisperse/aero_data.dat | 4 +- test/brownian/aero_data.dat | 4 +- test/condense/aero_data.dat | 42 +-- test/emission/aero_data.dat | 8 +- test/fractal/aero_data.dat | 4 +- test/freezing/aero_data.dat | 37 ++- test/loss/aero_data.dat | 4 +- test/mixing_state/aero_data.dat | 10 +- test/mosaic/aero_data.dat | 42 +-- test/nucleate/aero_data.dat | 6 +- test/parallel/aero_data.dat | 4 +- test/sedi/aero_data.dat | 4 +- test/tchem/aero_data.dat | 8 +- test/weighting/aero_data.dat | 8 +- 25 files changed, 584 insertions(+), 226 deletions(-) diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index 26fcd7e32..c741e0256 100644 --- a/scenarios/1_urban_plume/aero_data.dat +++ b/scenarios/1_urban_plume/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 0 96d-3 0.65 0 0 -NO3 1800 0 62d-3 0.65 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.65 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 0 60d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1000 0 1d-3 0.001 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 0 96d-3 0.65 0 0 0.073 +NO3 1800 0 62d-3 0.65 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.65 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 0 60d-3 0.53 0 0 0.073 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1000 0 1d-3 0.001 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/scenarios/1_urban_plume/urban_plume_process.F90 b/scenarios/1_urban_plume/urban_plume_process.F90 index 6b5e9f03e..111e28e81 100644 --- a/scenarios/1_urban_plume/urban_plume_process.F90 +++ b/scenarios/1_urban_plume/urban_plume_process.F90 @@ -26,10 +26,12 @@ program process real(kind=dp), allocatable :: times(:), dry_diameters(:), num_concs(:), & dry_masses(:), masses(:), bc_masses(:), bc_fracs(:), & crit_rhs(:), scs(:), num_dist(:), & - diam_bc_dist(:,:), diam_sc_dist(:,:) + diam_bc_dist(:,:), diam_sc_dist(:,:), & + crit_rhs_est(:), scs_est(:), diam_sc_est_dist(:,:) type(stats_1d_t) :: stats_num_dist, stats_d_alpha, stats_tot_num_conc, & stats_tot_mass_conc, stats_d_gamma, stats_chi - type(stats_2d_t) :: stats_diam_bc_dist, stats_diam_sc_dist + type(stats_2d_t) :: stats_diam_bc_dist, stats_diam_sc_dist, & + stats_diam_sc_est_dist character(len=AERO_NAME_LEN), allocatable :: mixing_state_groups(:,:) call pmc_mpi_init() @@ -48,6 +50,7 @@ program process mixing_state_groups(3,:) = ["SO4 ", "NO3 ", "NH4 ", " "] scs = [ real(kind=dp) :: ] ! silence compiler warnings + scs_est = [ real(kind=dp) :: ] bc_fracs = [ real(kind=dp) :: ] do i_index = 1,n_index @@ -90,6 +93,14 @@ program process sc_grid, scs, num_concs) call stats_2d_add(stats_diam_sc_dist, diam_sc_dist) + ! critical supersaturation using the effective surface tension (EST) + crit_rhs_est = aero_state_crit_rel_humids_est(aero_state, aero_data, & + env_state) + scs_est = crit_rhs_est - 1d0 + diam_sc_est_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + sc_grid, scs_est, num_concs) + call stats_2d_add(stats_diam_sc_est_dist, diam_sc_est_dist) + call aero_state_mixing_state_metrics(aero_state, aero_data, & d_alpha, d_gamma, chi, groups=mixing_state_groups) @@ -119,6 +130,11 @@ program process dim_name_1="diam", dim_name_2="sc", unit="m^{-3}") call stats_2d_clear(stats_diam_sc_dist) + call stats_2d_output_netcdf(stats_diam_sc_est_dist, ncid, & + "diam_sc_est_dist", dim_name_1="diam", dim_name_2="sc", & + unit="m^{-3}") + call stats_2d_clear(stats_diam_sc_est_dist) + call pmc_nc_close(ncid) end do diff --git a/scenarios/2_urban_plume2/aero_data.dat b/scenarios/2_urban_plume2/aero_data.dat index 26fcd7e32..d05878897 100644 --- a/scenarios/2_urban_plume2/aero_data.dat +++ b/scenarios/2_urban_plume2/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 0 96d-3 0.65 0 0 -NO3 1800 0 62d-3 0.65 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.65 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 0 60d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1000 0 1d-3 0.001 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 0 96d-3 0.65 0 0 0.073 +NO3 1800 0 62d-3 0.65 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.65 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 0 60d-3 0.53 0 0 0.073 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1000 0 1d-3 0.001 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/scenarios/3_condense/aero_data.dat b/scenarios/3_condense/aero_data.dat index d0cd57cfa..0e4980b47 100644 --- a/scenarios/3_condense/aero_data.dat +++ b/scenarios/3_condense/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 1 96d-3 0 0 0 -NO3 1800 1 62d-3 0 0 0 -Cl 2200 1 35.5d-3 0 0 0 -NH4 1800 1 18d-3 0 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 1 60d-3 0 0 0 -Na 2200 1 23d-3 0 0 0 -Ca 2600 1 40d-3 0 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1000 0 1d-3 0.1 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 1 96d-3 0 0 0 0.073 +NO3 1800 1 62d-3 0 0 0 0.073 +Cl 2200 1 35.5d-3 0 0 0 0.073 +NH4 1800 1 18d-3 0 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 1 60d-3 0 0 0 0.073 +Na 2200 1 23d-3 0 0 0 0.073 +Ca 2600 1 40d-3 0 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1000 0 1d-3 0.1 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/scenarios/4_chamber/aero_data.dat b/scenarios/4_chamber/aero_data.dat index 26fcd7e32..d05878897 100644 --- a/scenarios/4_chamber/aero_data.dat +++ b/scenarios/4_chamber/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 0 96d-3 0.65 0 0 -NO3 1800 0 62d-3 0.65 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.65 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 0 60d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1000 0 1d-3 0.001 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 0 96d-3 0.65 0 0 0.073 +NO3 1800 0 62d-3 0.65 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.65 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 0 60d-3 0.53 0 0 0.073 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1000 0 1d-3 0.001 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/scenarios/5_coag_brownian/aero_data.dat b/scenarios/5_coag_brownian/aero_data.dat index 26fcd7e32..d05878897 100644 --- a/scenarios/5_coag_brownian/aero_data.dat +++ b/scenarios/5_coag_brownian/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 0 96d-3 0.65 0 0 -NO3 1800 0 62d-3 0.65 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.65 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 0 60d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1000 0 1d-3 0.001 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 0 96d-3 0.65 0 0 0.073 +NO3 1800 0 62d-3 0.65 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.65 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 0 60d-3 0.53 0 0 0.073 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1000 0 1d-3 0.001 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/src/aero_data.F90 b/src/aero_data.F90 index d0fe85384..11b67593f 100644 --- a/src/aero_data.F90 +++ b/src/aero_data.F90 @@ -74,6 +74,8 @@ module pmc_aero_data real(kind=dp), allocatable :: abifm_m(:) !> Length \c aero_data_n_spec(aero_data), abifm_c (1). real(kind=dp), allocatable :: abifm_c(:) + !> Length \c aero_data_n_spec(aero_data), sigma (J m-2). + real(kind=dp), allocatable :: sigma(:) !> Length \c aero_data_n_source(aero_data), source names. character(len=AERO_SOURCE_NAME_LEN), allocatable :: source_name(:) !> Length \c aero_data_n_weight_classes, weight class names. @@ -487,9 +489,9 @@ subroutine spec_file_read_aero_data(file, aero_data) ! check the data size n_species = size(species_data, 1) - if (.not. ((size(species_data, 2) == 6) .or. (n_species == 0))) then + if (.not. ((size(species_data, 2) == 7) .or. (n_species == 0))) then call die_msg(428926381, 'each line in ' // trim(file%name) & - // ' should contain exactly 7 values') + // ' should contain exactly 8 values') end if ! allocate and copy over the data @@ -502,6 +504,7 @@ subroutine spec_file_read_aero_data(file, aero_data) call ensure_real_array_size(aero_data%kappa, n_species) call ensure_real_array_size(aero_data%abifm_m, n_species) call ensure_real_array_size(aero_data%abifm_c, n_species) + call ensure_real_array_size(aero_data%sigma, n_species) do i = 1,n_species aero_data%name(i) = species_name(i)(1:AERO_NAME_LEN) aero_data%density(i) = species_data(i,1) @@ -510,6 +513,7 @@ subroutine spec_file_read_aero_data(file, aero_data) aero_data%kappa(i) = species_data(i,4) aero_data%abifm_m(i) = species_data(i,5) aero_data%abifm_c(i) = species_data(i,6) + aero_data%sigma(i) = species_data(i,7) call assert_msg(232362742, & (aero_data%num_ions(i) == 0) .or. (aero_data%kappa(i) == 0d0), & "ions and kappa both non-zero for species " & @@ -584,6 +588,7 @@ integer function pmc_mpi_pack_size_aero_data(val) + pmc_mpi_pack_size_real_array(val%kappa) & + pmc_mpi_pack_size_real_array(val%abifm_m) & + pmc_mpi_pack_size_real_array(val%abifm_c) & + + pmc_mpi_pack_size_real_array(val%sigma) & + pmc_mpi_pack_size_string_array(val%source_name) & + pmc_mpi_pack_size_string_array(val%weight_class_name) & + pmc_mpi_pack_size_fractal(val%fractal) @@ -616,6 +621,7 @@ subroutine pmc_mpi_pack_aero_data(buffer, position, val) call pmc_mpi_pack_real_array(buffer, position, val%kappa) call pmc_mpi_pack_real_array(buffer, position, val%abifm_m) call pmc_mpi_pack_real_array(buffer, position, val%abifm_c) + call pmc_mpi_pack_real_array(buffer, position, val%sigma) call pmc_mpi_pack_string_array(buffer, position, val%source_name) call pmc_mpi_pack_string_array(buffer, position, val%weight_class_name) call pmc_mpi_pack_fractal(buffer, position, val%fractal) @@ -651,6 +657,7 @@ subroutine pmc_mpi_unpack_aero_data(buffer, position, val) call pmc_mpi_unpack_real_array(buffer, position, val%kappa) call pmc_mpi_unpack_real_array(buffer, position, val%abifm_m) call pmc_mpi_unpack_real_array(buffer, position, val%abifm_c) + call pmc_mpi_unpack_real_array(buffer, position, val%sigma) call pmc_mpi_unpack_string_array(buffer, position, val%source_name) call pmc_mpi_unpack_string_array(buffer, position, val%weight_class_name) call pmc_mpi_unpack_fractal(buffer, position, val%fractal) @@ -942,6 +949,9 @@ subroutine aero_data_output_netcdf(aero_data, ncid) call pmc_nc_write_real_1d(ncid, aero_data%abifm_c, & "aero_abifm_c", (/ dimid_aero_species /), unit="1", & long_name="c parameter of ABIFM") + call pmc_nc_write_real_1d(ncid, aero_data%sigma, & + "aero_sigma", (/ dimid_aero_species /), unit="J/m2", & + long_name="surface tension (sigma) of aerosol species") call pmc_nc_write_integer(ncid, aero_data%i_water, & "aero_i_water", long_name="Index of aerosol water or " & // "0 if water does not exist.") @@ -993,6 +1003,7 @@ subroutine aero_data_input_netcdf(aero_data, ncid) call pmc_nc_read_real_1d(ncid, aero_data%kappa, "aero_kappa") call pmc_nc_read_real_1d(ncid, aero_data%abifm_m, "aero_abifm_m") call pmc_nc_read_real_1d(ncid, aero_data%abifm_c, "aero_abifm_c") + call pmc_nc_read_real_1d(ncid, aero_data%sigma, "aero_sigma") call pmc_nc_check(nf90_inq_varid(ncid, "aero_species", & varid_aero_species)) @@ -1118,6 +1129,7 @@ subroutine aero_data_initialize(aero_data, camp_core) allocate(aero_data%kappa(num_spec)) allocate(aero_data%abifm_m(num_spec)) allocate(aero_data%abifm_c(num_spec)) + allocate(aero_data%sigma(num_spec)) allocate(aero_data%camp_particle_spec_id(num_spec)) ! Assume no aerosol water @@ -1167,6 +1179,12 @@ subroutine aero_data_initialize(aero_data, camp_core) call die_msg(944207346, "Missing abifm_c for aerosol species " & // spec_names(i_spec)%string) end if + prop_name = "sigma" + if (.not. property_set%get_real(prop_name, & + aero_data%sigma(i_spec))) then + call die_msg(944207347, "Missing sigma for aerosol species " & + // spec_names(i_spec)%string) + end if prop_name = "PartMC name" if (property_set%get_string(prop_name, str_val)) then if (str_val == "H2O") then diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index fb6712af2..dded97e95 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -13,6 +13,7 @@ module pmc_aero_particle use pmc_aero_data use pmc_spec_file use pmc_env_state + use pmc_constants use pmc_mpi #ifdef PMC_USE_MPI use mpi @@ -795,7 +796,10 @@ end function aero_particle_approx_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical relative humidity (1). + !> Returns the critical relative humidity (1) using a constant water + !> surface tension. For the effective surface tension (EST) treatment + !> that accounts for organic surface films, use + !> aero_particle_crit_rel_humid_est instead. real(kind=dp) function aero_particle_crit_rel_humid(aero_particle, & aero_data, env_state) @@ -824,6 +828,41 @@ end function aero_particle_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical relative humidity (1) using the effective + !> surface tension (EST), accounting for organic surface films via the + !> per-species sigma in aero_data. For the standard constant surface + !> tension treatment, use aero_particle_crit_rel_humid. + real(kind=dp) function aero_particle_crit_rel_humid_est(aero_particle, & + aero_data, env_state) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(in) :: env_state + + real(kind=dp) :: A, dry_diam, kappa, d, eff_surf_eng + + ! Kelvin A parameter without surface tension; the effective surface + ! tension (eff_surf_eng) is applied explicitly below. + A = 4d0 * const%water_molec_weight & + / (const%univ_gas_const * env_state%temp * const%water_density) + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + d = aero_particle_crit_diameter_est(aero_particle, aero_data, env_state, & + eff_surf_eng) + + if (d == dry_diam) then + aero_particle_crit_rel_humid_est = exp(A * eff_surf_eng / dry_diam) + else + aero_particle_crit_rel_humid_est = (d**3 - dry_diam**3) & + / (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * eff_surf_eng / d) + end if + + end function aero_particle_crit_rel_humid_est + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the critical diameter (m). !! !! The method is as follows. We need to solve the polynomial @@ -909,7 +948,266 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & aero_particle_crit_diameter = d end function aero_particle_crit_diameter + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the critical diameter (m) using the effective surface + !> tension (EST). Also returns the effective surface tension at the + !> critical diameter through eff_surf_eng. + real(kind=dp) function aero_particle_crit_diameter_est(& + aero_particle, aero_data, env_state, eff_surf_eng) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(in) :: env_state + !> Effective surface tension at the critical diameter (J m^-2). + real(kind=dp), intent(out) :: eff_surf_eng + + integer, parameter :: CRIT_DIAM_MAX_ITER = 100 + integer :: i_newton + real(kind=dp) :: delta_min = 1.6d-10 + real(kind=dp) :: kappa, dry_diam, A + real(kind=dp) :: surf_eng_soluble, surf_eng_organic + real(kind=dp) :: v_solid, v_sol, v_org, r_solid + real(kind=dp) :: d, v_delta_solid, r_core, v_delta, d_v_delta + real(kind=dp) :: d_eff_surf_eng, d_2_eff_surf_eng, R, d_R, f, df, dd + real(kind=dp) :: c_1, c_2, c_3, c_4, c_5 + + ! Kelvin A parameter without surface tension; the effective surface + ! tension is applied explicitly through the eff_surf_eng terms below. + A = 4d0 * const%water_molec_weight & + / (const%univ_gas_const * env_state%temp * const%water_density) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + surf_eng_organic = aero_particle_surf_eng_organic(aero_particle, aero_data) + v_solid = aero_particle_solid_volume(aero_particle, aero_data) + v_org = aero_particle_organic_volume(aero_particle, aero_data) + + if (kappa < 1d-30) then + ! bail out early for hydrophobic particles + aero_particle_crit_diameter_est = dry_diam + return + end if + + c_1 = 3d0 * dry_diam**3 * kappa / A + c_2= (2d0 - kappa) * dry_diam**3 + c_3 = (1d0 - kappa) * dry_diam**6 + c_4 = 2d0 * const%pi * delta_min ! d_2_v_delta_d + + d = 50*sqrt(4d0 / 3d0 * c_1) + + if (v_solid == 0d0) then + do i_newton = 1, CRIT_DIAM_MAX_ITER + surf_eng_soluble = aero_particle_surf_eng_soluble(aero_particle, & + aero_data, env_state, d) + c_5 = v_org * (surf_eng_organic - surf_eng_soluble) + v_sol = const%pi * d**3 / 6d0 - v_org ! volume of soluble inorganics + water + r_core = ((3d0 * v_sol)/(4d0 * const%pi))**(1d0 / 3d0) + v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) + if (v_org > v_delta) then + eff_surf_eng = surf_eng_organic + f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 + else + if (v_org == 0) then + eff_surf_eng = surf_eng_soluble + f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 + else + eff_surf_eng = surf_eng_soluble + c_5 / v_delta + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_eff_surf_eng = - c_5 * d_v_delta / v_delta**2 + R = eff_surf_eng - d * d_eff_surf_eng + d_2_eff_surf_eng = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_eff_surf_eng + f = R * (d**6 - c_2 * d**3 + c_3) - c_1 * d**4 + df = d_R * (d**6 - c_2 * d**3 + c_3) + R * (6d0 * d**5 - 3d0 * c_2 * d**2) & + - 4d0 * c_1 * d**3 + end if + end if + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-11) then + exit + end if + end do + else + r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) + v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) + do i_newton = 1, CRIT_DIAM_MAX_ITER + surf_eng_soluble = aero_particle_surf_eng_soluble(aero_particle, & + aero_data, env_state, d) + c_5 = v_org * (surf_eng_organic - surf_eng_soluble) + v_sol = const%pi * d**3 / 6 - v_solid - v_org + if (v_sol + v_org > v_delta_solid) then + r_core = ((3d0 * (v_sol + v_solid))/(4 * const%pi))**(1d0 / 3d0) + v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) + if (v_org > v_delta) then + eff_surf_eng = surf_eng_organic + f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 + else + if (v_org == 0) then + eff_surf_eng = surf_eng_soluble + f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 + else + eff_surf_eng = surf_eng_soluble + c_5 / v_delta + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_eff_surf_eng = - c_5 * d_v_delta / v_delta**2 + R = eff_surf_eng - d * d_eff_surf_eng + d_2_eff_surf_eng = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_eff_surf_eng + f = R * (d**6 - c_2 * d**3 + c_3) - c_1 * d**4 + df = d_R * (d**6 - c_2 * d**3 + c_3) + R * (6d0 * d**5 - 3d0 * c_2 * d**2) & + - 4d0 * c_1 * d**3 + end if + end if + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-11) then + exit + end if + else + eff_surf_eng = (v_org * surf_eng_organic + v_sol * surf_eng_soluble) & + / v_delta_solid + d = dry_diam + exit + end if + end do + end if + + call warn_assert_msg(408545686, i_newton < CRIT_DIAM_MAX_ITER, & + "critical diameter for new Newton loop failed to converge") + call warn_assert_msg(353290871, d >= dry_diam, & + "critical diameter for new Newton loop converged to invalid solution") + aero_particle_crit_diameter_est = d + + end function aero_particle_crit_diameter_est +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Total organic volume in the particle (m^3). + real(kind=dp) function aero_particle_organic_volume(aero_particle, aero_data) + + !> Particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + integer :: i_org_spec, n_org_spec + character(len=AERO_NAME_LEN), parameter, dimension(10) :: & + org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + + aero_particle_organic_volume = 0d0 + + do n_org_spec = 1, size(org_spec) + i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) + aero_particle_organic_volume = aero_particle_organic_volume & + + aero_particle%vol(i_org_spec) + end do + + end function aero_particle_organic_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Total solid (insoluble) volume in the particle (m^3). + real(kind=dp) function aero_particle_solid_volume(aero_particle, aero_data) + + !> Particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + integer :: i_solid_spec, n_solid_spec + character(len=AERO_NAME_LEN), parameter, dimension(2) :: & + solid_spec = ["OIN ", "BC "] + + aero_particle_solid_volume = 0d0 + + do n_solid_spec = 1, size(solid_spec) + i_solid_spec = aero_data_spec_by_name(aero_data, solid_spec(n_solid_spec)) + aero_particle_solid_volume = aero_particle_solid_volume & + + aero_particle%vol(i_solid_spec) + end do + + end function aero_particle_solid_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the volume-weighted surface tension of the soluble core + !> (soluble inorganics plus water) at diameter d (J m^-2). + real(kind=dp) function aero_particle_surf_eng_soluble(aero_particle, & + aero_data, env_state, d) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(in) :: env_state + !> Diameter at which to evaluate the surface tension (m). + real(kind=dp) :: d + + integer :: i_sol_spec, n_sol_spec + real(kind=dp) :: aero_particle_soluble_volume + real(kind=dp) :: v_solid, v_org, v_sol, v_water + character(len=AERO_NAME_LEN), parameter, dimension(7) :: & + sol_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + "Na ", "Ca "] + + v_solid = aero_particle_solid_volume(aero_particle, aero_data) + v_org = aero_particle_organic_volume(aero_particle, aero_data) + + aero_particle_soluble_volume = 0d0 + + do n_sol_spec = 1, size(sol_spec) + i_sol_spec = aero_data_spec_by_name(aero_data, sol_spec(n_sol_spec)) + aero_particle_soluble_volume = aero_particle_soluble_volume & + + aero_particle%vol(i_sol_spec) + end do + + v_sol = d**3 * const%pi / 6 - v_solid - v_org + v_water = v_sol - aero_particle_soluble_volume + + aero_particle_surf_eng_soluble = 0d0 + do n_sol_spec = 1, size(sol_spec) + i_sol_spec = aero_data_spec_by_name(aero_data, sol_spec(n_sol_spec)) + aero_particle_surf_eng_soluble = aero_particle_surf_eng_soluble & + + aero_particle%vol(i_sol_spec) * & + aero_data%sigma(i_sol_spec) / v_sol + end do + aero_particle_surf_eng_soluble = aero_particle_surf_eng_soluble + v_water & + * const%water_surf_eng / v_sol + + end function aero_particle_surf_eng_soluble + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the volume-weighted surface tension of the organic film + !> (J m^-2). + real(kind=dp) function aero_particle_surf_eng_organic(aero_particle, aero_data) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + + real(kind=dp) :: org_volume + integer :: i_org_spec, n_org_spec + character(len=AERO_NAME_LEN), parameter, dimension(10) :: & + org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + + org_volume = aero_particle_organic_volume(aero_particle, aero_data) + aero_particle_surf_eng_organic = 0d0 + + do n_org_spec = 1, size(org_spec) + i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) + aero_particle_surf_eng_organic = aero_particle_surf_eng_organic & + + aero_particle%vol(i_org_spec) * & + aero_data%sigma(i_org_spec) / org_volume + end do + end function aero_particle_surf_eng_organic !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Coagulate two particles together to make a new one. The new diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 528c709aa..0cb5b063b 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1677,6 +1677,33 @@ end function aero_state_crit_rel_humids !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical relative humidity for all particles (1) using + !> the effective surface tension (EST). + function aero_state_crit_rel_humids_est(aero_state, aero_data, env_state) + + !> Aerosol state. + type(aero_state_t), intent(in) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(in) :: env_state + + !> Return value. + real(kind=dp) :: aero_state_crit_rel_humids_est(aero_state_n_part(aero_state)) + + integer :: i_part + + do i_part = 1,aero_state_n_part(aero_state) + aero_state_crit_rel_humids_est(i_part) = & + aero_particle_crit_rel_humid_est( & + aero_state%apa%particle(i_part), aero_data, env_state) + end do + + end function aero_state_crit_rel_humids_est + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Does the same thing as aero_state_to_bin() but based on dry radius. subroutine aero_state_to_binned_dry(bin_grid, aero_data, aero_state, & aero_binned) diff --git a/test/additive/aero_data.dat b/test/additive/aero_data.dat index b38ac540e..11c919d1b 100644 --- a/test/additive/aero_data.dat +++ b/test/additive/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/average/aero_data.dat b/test/average/aero_data.dat index 65d773f46..3dd46c347 100644 --- a/test/average/aero_data.dat +++ b/test/average/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 -SO4 1800 0 96d-3 0.53 0 0 -NO3 1800 0 62d-3 0.53 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.53 0 0 -CO3 2600 0 60d-3 0.53 0 0 -MSA 1800 0 95d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OC 1000 0 1d-3 0.1 0 0 -BC 1800 0 1d-3 0 0 0 -OIN 2600 0 1d-3 0.1 0 0 -ARO1 1000 0 150d-3 0.1 0 0 -ARO2 1000 0 150d-3 0.1 0 0 -ALK1 1000 0 140d-3 0.1 0 0 -OLE1 1000 0 140d-3 0.1 0 0 -API1 1000 0 184d-3 0.1 0 0 -API2 1000 0 184d-3 0.1 0 0 -LIM1 1000 0 200d-3 0.1 0 0 -LIM2 1000 0 200d-3 0.1 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 +SO4 1800 0 96d-3 0.53 0 0 0.073 +NO3 1800 0 62d-3 0.53 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.53 0 0 0.073 +CO3 2600 0 60d-3 0.53 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OC 1000 0 1d-3 0.1 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 0 +ARO1 1000 0 150d-3 0.1 0 0 0.03 +ARO2 1000 0 150d-3 0.1 0 0 0.03 +ALK1 1000 0 140d-3 0.1 0 0 0.03 +OLE1 1000 0 140d-3 0.1 0 0 0.03 +API1 1000 0 184d-3 0.1 0 0 0.03 +API2 1000 0 184d-3 0.1 0 0 0.03 +LIM1 1000 0 200d-3 0.1 0 0 0.03 +LIM2 1000 0 200d-3 0.1 0 0 0.03 diff --git a/test/bidisperse/aero_data.dat b/test/bidisperse/aero_data.dat index cc430aaee..b29e086c6 100644 --- a/test/bidisperse/aero_data.dat +++ b/test/bidisperse/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/brownian/aero_data.dat b/test/brownian/aero_data.dat index cc430aaee..b29e086c6 100644 --- a/test/brownian/aero_data.dat +++ b/test/brownian/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/condense/aero_data.dat b/test/condense/aero_data.dat index 9d466783c..23d3186ad 100644 --- a/test/condense/aero_data.dat +++ b/test/condense/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 1 96d-3 0 0 0 -NO3 1800 1 62d-3 0 0 0 -Cl 2200 1 35.5d-3 0 0 0 -NH4 1800 1 18d-3 0 0 0 -MSA 1800 0 95d-3 0.53 0 0 -ARO1 1400 0 150d-3 0.1 0 0 -ARO2 1400 0 150d-3 0.1 0 0 -ALK1 1400 0 140d-3 0.1 0 0 -OLE1 1400 0 140d-3 0.1 0 0 -API1 1400 0 184d-3 0.1 0 0 -API2 1400 0 184d-3 0.1 0 0 -LIM1 1400 0 200d-3 0.1 0 0 -LIM2 1400 0 200d-3 0.1 0 0 -CO3 2600 1 60d-3 0 0 0 -Na 2200 1 23d-3 0 0 0 -Ca 2600 1 40d-3 0 0 0 -OIN 2600 0 1d-3 0.1 0 0 -OC 1400 0 1d-3 0.1 0 0 -BC 1800 0 1d-3 0 0 0 -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 1 96d-3 0 0 0 0.073 +NO3 1800 1 62d-3 0 0 0 0.073 +Cl 2200 1 35.5d-3 0 0 0 0.073 +NH4 1800 1 18d-3 0 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +ARO1 1400 0 150d-3 0.1 0 0 0.03 +ARO2 1400 0 150d-3 0.1 0 0 0.03 +ALK1 1400 0 140d-3 0.1 0 0 0.03 +OLE1 1400 0 140d-3 0.1 0 0 0.03 +API1 1400 0 184d-3 0.1 0 0 0.03 +API2 1400 0 184d-3 0.1 0 0 0.03 +LIM1 1400 0 200d-3 0.1 0 0 0.03 +LIM2 1400 0 200d-3 0.1 0 0 0.03 +CO3 2600 1 60d-3 0 0 0 0.073 +Na 2200 1 23d-3 0 0 0 0.073 +Ca 2600 1 40d-3 0 0 0 0.073 +OIN 2600 0 1d-3 0.1 0 0 0 +OC 1400 0 1d-3 0.1 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/emission/aero_data.dat b/test/emission/aero_data.dat index d52f270c3..a0b404374 100644 --- a/test/emission/aero_data.dat +++ b/test/emission/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SI 10000 0 18d-3 0 0 0 -SB 100 0 18d-3 0 0 0 -SE 1500 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) surface tension (J m^{-2}) +SI 10000 0 18d-3 0 0 0 0.073 +SB 100 0 18d-3 0 0 0 0.073 +SE 1500 0 18d-3 0 0 0 0.073 diff --git a/test/fractal/aero_data.dat b/test/fractal/aero_data.dat index 9d5459403..7cb98d46f 100644 --- a/test/fractal/aero_data.dat +++ b/test/fractal/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -BC 4200 0 1d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +BC 4200 0 1d-3 0 0 0 0 diff --git a/test/freezing/aero_data.dat b/test/freezing/aero_data.dat index 3735cb147..d1c480b66 100644 --- a/test/freezing/aero_data.dat +++ b/test/freezing/aero_data.dat @@ -1,30 +1,29 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 -#OIN 2600 0 1d-3 0.1 28.13797 -2.92414 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 +#OIN 2600 0 1d-3 0.1 28.13797 -2.92414 0.0 -ILT 2750 0 389.34d-3 0.003 54.48075 -10.66873 +ILT 2750 0 389.34d-3 0.003 54.48075 -10.66873 0.073 # Illite: # dens, molec wght: http://www.webmineral.com/data/Illite.shtml#.ZA0JpezMJ_Q -# kappa: 0.1039/b901585j +# kappa: 0.1039/b901585j sigma: [unknown] -#Al2O3 3970 0 101.96d-3 0.01 14.96639 1.60671 +#Al2O3 3970 0 101.96d-3 0.01 14.96639 1.60671 0.0 # Aluminium oxide: -# dens, molec wght: https://en.wikipedia.org/wiki/Aluminium_oxide -# kappa: [unknown] +# dens, molec wght: https://en.wikipedia.org/wiki/Aluminium_oxide surface tension (J m^{-2}) +# kappa: [unknown] sigma: [unknown] -Fe2O3 5240 0 159.69d-3 0.01 17.62106 1.42411 +Fe2O3 5240 0 159.69d-3 0.01 17.62106 1.42411 0.0 # Iron(III) oxide: -# den, molec wght: https://en.wikipedia.org/wiki/Iron(III)_oxide -# kappa: [unknown] +# den, molec wght: https://en.wikipedia.org/wiki/Iron(III)_oxide surface tension (J m^{-2}) +# kappa: [unknown] sigma: [unknown] -#KLN 2650 0 258.16d-3 0.023 54.58834 -10.54758 +#KLN 2650 0 258.16d-3 0.023 54.58834 -10.54758 0.0 # Kaolinite -# den: https://cameochemicals.noaa.gov/chemical/25036 -# molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q -# kappa: 10.1039/b901585j -# +# den: https://cameochemicals.noaa.gov/chemical/25036 surface tension (J m^{-2}) +# molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q surface tension (J m^{-2}) +# kappa: 10.1039/b901585j sigma: [unknown] -#NVF 5240 0 1d-3 0.01 0.0 0.0 +#NVF 5240 0 1d-3 0.01 0.0 0.0 0.0 # NeVer Freeze (other parameters are same as OIN) -#VOLC 2600 0 1d-3 0.01 28.13797 -2.92414 -# Volcano (Steinke et al., 2011) +#VOLC 2600 0 1d-3 0.01 28.13797 -2.92414 0.0 +# Volcano (Steinke et al., 2011) sigma: [unknown] diff --git a/test/loss/aero_data.dat b/test/loss/aero_data.dat index cc430aaee..b29e086c6 100644 --- a/test/loss/aero_data.dat +++ b/test/loss/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/mixing_state/aero_data.dat b/test/mixing_state/aero_data.dat index 9ee382e85..7d3d6465a 100644 --- a/test/mixing_state/aero_data.dat +++ b/test/mixing_state/aero_data.dat @@ -1,5 +1,5 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -A 1000 0 100d-3 0.1 0 0 -B 1500 0 100d-3 0.1 0 0 -C 2000 0 100d-3 0.1 0 0 -D 2500 0 100d-3 0.1 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +A 1000 0 100d-3 0.1 0 0 0.03 +B 1500 0 100d-3 0.1 0 0 0.03 +C 2000 0 100d-3 0.1 0 0 0.03 +D 2500 0 100d-3 0.1 0 0 0.03 diff --git a/test/mosaic/aero_data.dat b/test/mosaic/aero_data.dat index 65d773f46..3dd46c347 100644 --- a/test/mosaic/aero_data.dat +++ b/test/mosaic/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 -SO4 1800 0 96d-3 0.53 0 0 -NO3 1800 0 62d-3 0.53 0 0 -Cl 2200 0 35.5d-3 1.28 0 0 -NH4 1800 0 18d-3 0.53 0 0 -CO3 2600 0 60d-3 0.53 0 0 -MSA 1800 0 95d-3 0.53 0 0 -Na 2200 0 23d-3 1.28 0 0 -Ca 2600 0 40d-3 0.53 0 0 -OC 1000 0 1d-3 0.1 0 0 -BC 1800 0 1d-3 0 0 0 -OIN 2600 0 1d-3 0.1 0 0 -ARO1 1000 0 150d-3 0.1 0 0 -ARO2 1000 0 150d-3 0.1 0 0 -ALK1 1000 0 140d-3 0.1 0 0 -OLE1 1000 0 140d-3 0.1 0 0 -API1 1000 0 184d-3 0.1 0 0 -API2 1000 0 184d-3 0.1 0 0 -LIM1 1000 0 200d-3 0.1 0 0 -LIM2 1000 0 200d-3 0.1 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 +SO4 1800 0 96d-3 0.53 0 0 0.073 +NO3 1800 0 62d-3 0.53 0 0 0.073 +Cl 2200 0 35.5d-3 1.28 0 0 0.073 +NH4 1800 0 18d-3 0.53 0 0 0.073 +CO3 2600 0 60d-3 0.53 0 0 0.073 +MSA 1800 0 95d-3 0.53 0 0 0.03 +Na 2200 0 23d-3 1.28 0 0 0.073 +Ca 2600 0 40d-3 0.53 0 0 0.073 +OC 1000 0 1d-3 0.1 0 0 0.03 +BC 1800 0 1d-3 0 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 0 +ARO1 1000 0 150d-3 0.1 0 0 0.03 +ARO2 1000 0 150d-3 0.1 0 0 0.03 +ALK1 1000 0 140d-3 0.1 0 0 0.03 +OLE1 1000 0 140d-3 0.1 0 0 0.03 +API1 1000 0 184d-3 0.1 0 0 0.03 +API2 1000 0 184d-3 0.1 0 0 0.03 +LIM1 1000 0 200d-3 0.1 0 0 0.03 +LIM2 1000 0 200d-3 0.1 0 0 0.03 diff --git a/test/nucleate/aero_data.dat b/test/nucleate/aero_data.dat index 3b62e55ae..4410a4a15 100644 --- a/test/nucleate/aero_data.dat +++ b/test/nucleate/aero_data.dat @@ -1,3 +1,3 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SO4 1800 0 96d-3 0.65 0 0 -NO3 1800 0 62d-3 0.65 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SO4 1800 0 96d-3 0.65 0 0 0.073 +NO3 1800 0 62d-3 0.65 0 0 0.073 diff --git a/test/parallel/aero_data.dat b/test/parallel/aero_data.dat index cc430aaee..b29e086c6 100644 --- a/test/parallel/aero_data.dat +++ b/test/parallel/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/sedi/aero_data.dat b/test/sedi/aero_data.dat index 5a703b6c9..8a778915b 100644 --- a/test/sedi/aero_data.dat +++ b/test/sedi/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -H2O 1000 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +H2O 1000 0 18d-3 0 0 0 0.073 diff --git a/test/tchem/aero_data.dat b/test/tchem/aero_data.dat index 892067a0d..6803c3852 100644 --- a/test/tchem/aero_data.dat +++ b/test/tchem/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SI 10000 0 18d-3 0 0 0 -SB 100 0 18d-3 0 0 0 -SE 1500 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SI 10000 0 18d-3 0 0 0 0.0 +SB 100 0 18d-3 0 0 0 0.0 +SE 1500 0 18d-3 0 0 0 0.0 diff --git a/test/weighting/aero_data.dat b/test/weighting/aero_data.dat index 892067a0d..01f80caec 100644 --- a/test/weighting/aero_data.dat +++ b/test/weighting/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) -SI 10000 0 18d-3 0 0 0 -SB 100 0 18d-3 0 0 0 -SE 1500 0 18d-3 0 0 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) surface tension (J m^{-2}) +SI 10000 0 18d-3 0 0 0 0.073 +SB 100 0 18d-3 0 0 0 0.073 +SE 1500 0 18d-3 0 0 0 0.073 From e05a51870f40b052478273a5043852ed579c75ad Mon Sep 17 00:00:00 2001 From: Xiaotian Xu Date: Thu, 16 Jul 2026 19:32:08 -0400 Subject: [PATCH 2/5] fix format --- test/freezing/aero_data.dat | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/test/freezing/aero_data.dat b/test/freezing/aero_data.dat index d1c480b66..3db10fd73 100644 --- a/test/freezing/aero_data.dat +++ b/test/freezing/aero_data.dat @@ -5,25 +5,30 @@ H2O 1000 0 18d-3 0 ILT 2750 0 389.34d-3 0.003 54.48075 -10.66873 0.073 # Illite: # dens, molec wght: http://www.webmineral.com/data/Illite.shtml#.ZA0JpezMJ_Q -# kappa: 0.1039/b901585j sigma: [unknown] +# kappa: 0.1039/b901585j +# sigma: [unknown] #Al2O3 3970 0 101.96d-3 0.01 14.96639 1.60671 0.0 # Aluminium oxide: -# dens, molec wght: https://en.wikipedia.org/wiki/Aluminium_oxide surface tension (J m^{-2}) -# kappa: [unknown] sigma: [unknown] +# dens, molec wght: https://en.wikipedia.org/wiki/Aluminium_oxide +# kappa: [unknown] +# sigma: [unknown] Fe2O3 5240 0 159.69d-3 0.01 17.62106 1.42411 0.0 # Iron(III) oxide: -# den, molec wght: https://en.wikipedia.org/wiki/Iron(III)_oxide surface tension (J m^{-2}) -# kappa: [unknown] sigma: [unknown] +# den, molec wght: https://en.wikipedia.org/wiki/Iron(III)_oxide +# kappa: [unknown] +# sigma: [unknown] #KLN 2650 0 258.16d-3 0.023 54.58834 -10.54758 0.0 # Kaolinite # den: https://cameochemicals.noaa.gov/chemical/25036 surface tension (J m^{-2}) -# molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q surface tension (J m^{-2}) -# kappa: 10.1039/b901585j sigma: [unknown] +# molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q +# kappa: 10.1039/b901585j +# sigma: [unknown] #NVF 5240 0 1d-3 0.01 0.0 0.0 0.0 # NeVer Freeze (other parameters are same as OIN) #VOLC 2600 0 1d-3 0.01 28.13797 -2.92414 0.0 -# Volcano (Steinke et al., 2011) sigma: [unknown] +# Volcano (Steinke et al., 2011) +# sigma: [unknown] From 592fb7da7261d1c7e5d5931ee07ff3a982571b95 Mon Sep 17 00:00:00 2001 From: Xiaotian Xu Date: Thu, 16 Jul 2026 19:34:36 -0400 Subject: [PATCH 3/5] Fix formatting --- test/freezing/aero_data.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/freezing/aero_data.dat b/test/freezing/aero_data.dat index 3db10fd73..bc1c606a4 100644 --- a/test/freezing/aero_data.dat +++ b/test/freezing/aero_data.dat @@ -22,7 +22,7 @@ Fe2O3 5240 0 159.69d-3 0.01 #KLN 2650 0 258.16d-3 0.023 54.58834 -10.54758 0.0 # Kaolinite -# den: https://cameochemicals.noaa.gov/chemical/25036 surface tension (J m^{-2}) +# den: https://cameochemicals.noaa.gov/chemical/25036 # molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q # kappa: 10.1039/b901585j # sigma: [unknown] From 36c06a2a66a541b8e4d3fbd3ec48287660768c4b Mon Sep 17 00:00:00 2001 From: Xiaotian Xu Date: Tue, 28 Jul 2026 08:15:29 -0500 Subject: [PATCH 4/5] add sigma for camp --- .../6_camp/monarch_mod37/custom_species.json | 21 ++++++++++++------- .../species.json | 5 +++++ .../species.json | 9 ++++++++ 3 files changed, 28 insertions(+), 7 deletions(-) diff --git a/scenarios/6_camp/monarch_mod37/custom_species.json b/scenarios/6_camp/monarch_mod37/custom_species.json index 0ef9e68ca..353b502eb 100644 --- a/scenarios/6_camp/monarch_mod37/custom_species.json +++ b/scenarios/6_camp/monarch_mod37/custom_species.json @@ -30,7 +30,8 @@ "num_ions" : 0, "kappa" : 0.1, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0 }, { "monarch name" : "lumped low-density dust species", @@ -43,7 +44,8 @@ "num_ions" : 0, "kappa" : 0.1, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0 }, { "monarch name" : "lumped sea salt species", @@ -56,7 +58,8 @@ "num_ions" : 0, "kappa" : 0.53, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0.073 }, { "name" : "BC_phob", @@ -69,7 +72,8 @@ "num_ions" : 0, "kappa" : 0, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0 }, { "name" : "BC_phil", @@ -82,7 +86,8 @@ "num_ions" : 0, "kappa" : 0.001, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0 }, { "name" : "other_PM", @@ -95,7 +100,8 @@ "num_ions" : 0, "kappa" : 0, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0.073 }, { "name" : "other_other_PM", @@ -108,7 +114,8 @@ "num_ions" : 0, "kappa" : 0, "abifm_m": 0, - "abifm_c": 0 + "abifm_c": 0, + "sigma": 0.073 } ] } diff --git a/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json b/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json index 72bfcc6b0..535b43a6a 100644 --- a/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json +++ b/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json @@ -60,6 +60,7 @@ "kappa" : 0.0, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "note" : "Using C30H54 for molecular weight. TODO find best surrogate" }, { @@ -74,6 +75,7 @@ "kappa" : 0.1, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "note" : [ "Using ketopropanoic acid for molecular weight. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" @@ -91,6 +93,7 @@ "kappa" : 0.1, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "note" : [ "Using oxalic acid for molecular weight. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" @@ -108,6 +111,7 @@ "kappa" : 0.1, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "note" : [ "Using 2-hydroxy-3-isopropyl-6-methyl-cyclohexanone for molecular weight", "TODO find best surrogate", @@ -126,6 +130,7 @@ "kappa" : 0.1, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "note" : [ "Using 2-methyl-5-carboxy-2,4-hexodienoic acid. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" diff --git a/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json b/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json index e0c09ac1c..580b48b83 100644 --- a/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json +++ b/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json @@ -11,6 +11,7 @@ "molecular weight [kg mol-1]" : 0.08806, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "gas-phase product from isoprene oxidation by OH", "notes" : [ "using diffusion coefficient from dry deposition", @@ -24,6 +25,7 @@ "molecular weight [kg mol-1]" : 0.09003, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "gas-phase product from isoprene oxidation by O3", "notes" : [ "using diffusion coefficient from dry deposition", @@ -37,6 +39,7 @@ "molecular weight [kg mol-1]" : 0.17025, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "gas-phase product from monoterpene oxidation by OH", "notes" : [ "using diffusion coefficient from dry deposition", @@ -50,6 +53,7 @@ "molecular weight [kg mol-1]" : 0.202162, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "gas-phase product from monoterpene oxidation by O3", "notes" : [ "using diffusion coefficient from dry deposition", @@ -65,6 +69,7 @@ "molecular weight [kg mol-1]" : 0.41475, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "lumped hydrophobic particulate matter", "num_ions" : 0, "kappa" : 0.0, @@ -79,6 +84,7 @@ "molecular weight [kg mol-1]" : 0.08806, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "First isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -96,6 +102,7 @@ "molecular weight [kg mol-1]" : 0.09003, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "Second isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -113,6 +120,7 @@ "molecular weight [kg mol-1]" : 0.17025, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "First monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -131,6 +139,7 @@ "molecular weight [kg mol-1]" : 0.202162, "abifm_m": 0, "abifm_c": 0, + "sigma" : 0.03, "description" : "Second monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, From efedd713895defb6fd1867ba1171e3cdb3ed42d1 Mon Sep 17 00:00:00 2001 From: Xiaotian Xu Date: Mon, 21 Sep 2026 14:35:34 -0500 Subject: [PATCH 5/5] Fix EST critical-diameter solver to find the global Kohler maximum aero_particle_crit_diameter_est started Newton far from the root and, for organic-film particles whose Kohler curve S(D) is double-peaked, converged to the wrong (outer) root -- a larger critical diameter and lower critical supersaturation -- or diverged to NaN. Replace the far-start Newton with a coarse logarithmic scan for the global maximum of S(D) plus a bracketed Newton (rtsafe) refinement. The EST physics is unchanged; only the root finding. Also address code-review items in the org_film branch: - set eff_surf_eng on the hydrophobic (kappa < 1e-30) early return - guard vol(0) for absent species, 0/0 in surf_eng_organic, and v_sol <= 0 - retire the obsolete d >= dry_diam assert; repurpose the remaining warning to flag a maximum above the scan range - read aero_sigma with must_be_present=.false. for backward compatibility - use util.F90 sphere-geometry helpers (add sphere_diam2vol) - follow the i_spec/n_spec loop-variable convention - precompute per-particle soluble sums (aero_particle_soluble_sums) so the scan avoids repeated species-name lookups (~10x faster post-processing) ctest passes 14/14. --- src/aero_data.F90 | 10 +- src/aero_particle.F90 | 425 +++++++++++++++++++++++++++--------------- src/util.F90 | 13 ++ 3 files changed, 299 insertions(+), 149 deletions(-) diff --git a/src/aero_data.F90 b/src/aero_data.F90 index 11b67593f..4f3fcd9ae 100644 --- a/src/aero_data.F90 +++ b/src/aero_data.F90 @@ -1003,7 +1003,15 @@ subroutine aero_data_input_netcdf(aero_data, ncid) call pmc_nc_read_real_1d(ncid, aero_data%kappa, "aero_kappa") call pmc_nc_read_real_1d(ncid, aero_data%abifm_m, "aero_abifm_m") call pmc_nc_read_real_1d(ncid, aero_data%abifm_c, "aero_abifm_c") - call pmc_nc_read_real_1d(ncid, aero_data%sigma, "aero_sigma") + call pmc_nc_read_real_1d(ncid, aero_data%sigma, "aero_sigma", & + must_be_present=.false.) + if (size(aero_data%sigma) == 0) then + ! Backward compatibility: files written before per-species surface + ! tension was added have no aero_sigma variable. Default to the water + ! value, reproducing the earlier constant-surface-tension behaviour. + call ensure_real_array_size(aero_data%sigma, size(aero_data%density)) + aero_data%sigma = const%water_surf_eng + end if call pmc_nc_check(nf90_inq_varid(ncid, "aero_species", & varid_aero_species)) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index dded97e95..da4087bce 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -953,6 +953,22 @@ end function aero_particle_crit_diameter !> Returns the critical diameter (m) using the effective surface !> tension (EST). Also returns the effective surface tension at the !> critical diameter through eff_surf_eng. + !! + !! With an organic surface film the effective surface tension varies + !! with wet diameter (it transitions from the organic value on a thick + !! film to the soluble-core value once the film is diluted). This can + !! give the equilibrium (Kohler) saturation curve S(D) two local + !! maxima. The critical supersaturation is the largest barrier the + !! droplet must cross, i.e. the *global* maximum of S(D) for D > + !! D_dry, and the critical diameter is the D at which it occurs. + !! + !! We locate the global maximum with a coarse logarithmic scan of S(D) + !! and then refine the corresponding stationary point (a root of the + !! equilibrium condition f(D) = 0, where f = 0 at every extremum of + !! S(D)) with a bracketed Newton iteration (rtsafe). Unlike an + !! unbracketed Newton started far from D_dry, this cannot run away to a + !! negative diameter or converge onto the wrong (outer) root, so it is + !! robust across the whole particle population. real(kind=dp) function aero_particle_crit_diameter_est(& aero_particle, aero_data, env_state, eff_surf_eng) @@ -966,14 +982,16 @@ real(kind=dp) function aero_particle_crit_diameter_est(& real(kind=dp), intent(out) :: eff_surf_eng integer, parameter :: CRIT_DIAM_MAX_ITER = 100 - integer :: i_newton - real(kind=dp) :: delta_min = 1.6d-10 - real(kind=dp) :: kappa, dry_diam, A - real(kind=dp) :: surf_eng_soluble, surf_eng_organic - real(kind=dp) :: v_solid, v_sol, v_org, r_solid - real(kind=dp) :: d, v_delta_solid, r_core, v_delta, d_v_delta - real(kind=dp) :: d_eff_surf_eng, d_2_eff_surf_eng, R, d_R, f, df, dd - real(kind=dp) :: c_1, c_2, c_3, c_4, c_5 + integer, parameter :: CRIT_DIAM_N_SCAN = 300 + real(kind=dp), parameter :: CRIT_DIAM_SCAN_MAX_FACTOR = 1d4 + real(kind=dp), parameter :: delta_min = 1.6d-10 + integer :: i, i_newton, i_best + real(kind=dp) :: kappa, dry_diam, A, surf_eng_organic, v_solid, v_org + real(kind=dp) :: soluble_volume, sum_sol_sigma + real(kind=dp) :: r_solid, v_delta_solid, c_1, c_2, c_3, c_4 + real(kind=dp) :: d, f, df, eff, s, s_best, d_best, log_fac + real(kind=dp) :: x_lo, x_hi, f_lo, f_hi, dd, dxold + logical :: valid ! Kelvin A parameter without surface tension; the effective surface ! tension is applied explicitly through the eff_surf_eng terms below. @@ -984,107 +1002,213 @@ real(kind=dp) function aero_particle_crit_diameter_est(& surf_eng_organic = aero_particle_surf_eng_organic(aero_particle, aero_data) v_solid = aero_particle_solid_volume(aero_particle, aero_data) v_org = aero_particle_organic_volume(aero_particle, aero_data) + ! These soluble sums are constant for the particle, so resolve them (and the + ! species-name lookups they need) once here rather than on every scan point. + call aero_particle_soluble_sums(aero_particle, aero_data, soluble_volume, & + sum_sol_sigma) + + c_1 = 3d0 * dry_diam**3 * kappa / A + c_2 = (2d0 - kappa) * dry_diam**3 + c_3 = (1d0 - kappa) * dry_diam**6 + c_4 = 2d0 * const%pi * delta_min ! d_2_v_delta_d + if (v_solid > 0d0) then + r_solid = sphere_vol2rad(v_solid) + v_delta_solid = sphere_rad2vol(r_solid + delta_min) & + - sphere_rad2vol(r_solid) + end if if (kappa < 1d-30) then - ! bail out early for hydrophobic particles - aero_particle_crit_diameter_est = dry_diam + ! Hydrophobic particle: no soluble material, so there is no Köhler + ! activation barrier. Return the dry diameter; the effective surface + ! tension is undefined without a soluble core, so report it as zero + ! (the film model divides by a soluble volume that is zero at D_dry). + aero_particle_crit_diameter_est = dry_diam + eff_surf_eng = 0d0 return end if - c_1 = 3d0 * dry_diam**3 * kappa / A - c_2= (2d0 - kappa) * dry_diam**3 - c_3 = (1d0 - kappa) * dry_diam**6 - c_4 = 2d0 * const%pi * delta_min ! d_2_v_delta_d + ! Coarse logarithmic scan for the global maximum of the Kohler curve. + log_fac = log(CRIT_DIAM_SCAN_MAX_FACTOR) & + / real(CRIT_DIAM_N_SCAN - 1, kind=dp) + d_best = dry_diam + s_best = - huge(1d0) + i_best = 0 + do i = 1, CRIT_DIAM_N_SCAN + d = dry_diam * exp(log_fac * real(i - 1, kind=dp)) + call est_residual(d, f, df, eff, valid) + if (.not. valid) cycle + s = (d**3 - dry_diam**3) / (d**3 - dry_diam**3 * (1d0 - kappa)) & + * exp(A * eff / d) + if (s > s_best) then + s_best = s + d_best = d + i_best = i + end if + end do - d = 50*sqrt(4d0 / 3d0 * c_1) - - if (v_solid == 0d0) then - do i_newton = 1, CRIT_DIAM_MAX_ITER - surf_eng_soluble = aero_particle_surf_eng_soluble(aero_particle, & - aero_data, env_state, d) - c_5 = v_org * (surf_eng_organic - surf_eng_soluble) - v_sol = const%pi * d**3 / 6d0 - v_org ! volume of soluble inorganics + water - r_core = ((3d0 * v_sol)/(4d0 * const%pi))**(1d0 / 3d0) - v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) - if (v_org > v_delta) then - eff_surf_eng = surf_eng_organic - f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 - df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 - else - if (v_org == 0) then - eff_surf_eng = surf_eng_soluble - f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 - df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 - else - eff_surf_eng = surf_eng_soluble + c_5 / v_delta - d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) - d_eff_surf_eng = - c_5 * d_v_delta / v_delta**2 - R = eff_surf_eng - d * d_eff_surf_eng - d_2_eff_surf_eng = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) - d_R = - d * d_2_eff_surf_eng - f = R * (d**6 - c_2 * d**3 + c_3) - c_1 * d**4 - df = d_R * (d**6 - c_2 * d**3 + c_3) + R * (6d0 * d**5 - 3d0 * c_2 * d**2) & - - 4d0 * c_1 * d**3 - end if - end if - dd = f / df - d = d - dd - if (abs(dd / d) < 1d-11) then - exit - end if - end do - else - r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) - v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) - do i_newton = 1, CRIT_DIAM_MAX_ITER - surf_eng_soluble = aero_particle_surf_eng_soluble(aero_particle, & - aero_data, env_state, d) - c_5 = v_org * (surf_eng_organic - surf_eng_soluble) - v_sol = const%pi * d**3 / 6 - v_solid - v_org - if (v_sol + v_org > v_delta_solid) then - r_core = ((3d0 * (v_sol + v_solid))/(4 * const%pi))**(1d0 / 3d0) - v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) - if (v_org > v_delta) then - eff_surf_eng = surf_eng_organic - f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 - df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 - else - if (v_org == 0) then - eff_surf_eng = surf_eng_soluble - f = d**6 - c_1 * d**4 / eff_surf_eng - c_2 * d**3 + c_3 - df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff_surf_eng - 3d0 * c_2 * d**2 - else - eff_surf_eng = surf_eng_soluble + c_5 / v_delta - d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) - d_eff_surf_eng = - c_5 * d_v_delta / v_delta**2 - R = eff_surf_eng - d * d_eff_surf_eng - d_2_eff_surf_eng = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) - d_R = - d * d_2_eff_surf_eng - f = R * (d**6 - c_2 * d**3 + c_3) - c_1 * d**4 - df = d_R * (d**6 - c_2 * d**3 + c_3) + R * (6d0 * d**5 - 3d0 * c_2 * d**2) & - - 4d0 * c_1 * d**3 - end if - end if + if (i_best <= 1) then + ! no activation barrier resolved above D_dry: treat as non-activating + aero_particle_crit_diameter_est = dry_diam + call est_residual(dry_diam, f, df, eff_surf_eng, valid) + return + end if + + ! The one remaining failure mode: the global maximum lies beyond the top + ! of the scan range (i_best at the last node). The bracket expansion below + ! still tries to recover it, but flag it so it is not silently trusted. + call warn_assert_msg(408545686, i_best < CRIT_DIAM_N_SCAN, & + "critical diameter scan reached its upper bound; " & + // "the Köhler maximum may lie above CRIT_DIAM_SCAN_MAX_FACTOR * D_dry") + + ! Bracket the stationary point around the coarse maximum: f < 0 while + ! S is still rising and f > 0 once S is falling. + x_lo = dry_diam * exp(log_fac * real(i_best - 2, kind=dp)) + x_hi = dry_diam * exp(log_fac & + * real(min(i_best + 1, CRIT_DIAM_N_SCAN) - 1, kind=dp)) + call est_residual(x_lo, f_lo, df, eff, valid) + call est_residual(x_hi, f_hi, df, eff, valid) + do i = 1, CRIT_DIAM_MAX_ITER + if (f_lo <= 0d0) exit + x_lo = x_lo / 1.1d0 + call est_residual(x_lo, f_lo, df, eff, valid) + end do + do i = 1, CRIT_DIAM_MAX_ITER + if (f_hi >= 0d0) exit + x_hi = x_hi * 1.1d0 + call est_residual(x_hi, f_hi, df, eff, valid) + end do + + ! Safeguarded Newton (rtsafe) on f(D) = 0 within [x_lo, x_hi]. + d = d_best + dxold = x_hi - x_lo + dd = dxold + do i_newton = 1, CRIT_DIAM_MAX_ITER + call est_residual(d, f, df, eff, valid) + if (f < 0d0) then + x_lo = d + else + x_hi = d + end if + if ((df == 0d0) & + .or. (((d - x_hi) * df - f) * ((d - x_lo) * df - f) > 0d0) & + .or. (abs(2d0 * f) > abs(dxold * df))) then + dxold = dd + dd = 0.5d0 * (x_hi - x_lo) + d = x_lo + dd + else + dxold = dd dd = f / df d = d - dd - if (abs(dd / d) < 1d-11) then - exit - end if - else - eff_surf_eng = (v_org * surf_eng_organic + v_sol * surf_eng_soluble) & - / v_delta_solid - d = dry_diam - exit - end if - end do - end if + end if + if (abs(dd) < 1d-13 * d) exit + end do - call warn_assert_msg(408545686, i_newton < CRIT_DIAM_MAX_ITER, & - "critical diameter for new Newton loop failed to converge") - call warn_assert_msg(353290871, d >= dry_diam, & - "critical diameter for new Newton loop converged to invalid solution") + ! No convergence assert is needed here: the bracketed solver cannot leave + ! [x_lo, x_hi], and d >= dry_diam holds by construction (the scan starts at + ! D_dry). The only genuine failure mode is caught by the scan-range warning + ! above. (Warning 353290871 is retired for that reason.) + + ! Report the effective surface tension at the critical diameter. + call est_residual(d, f, df, eff_surf_eng, valid) aero_particle_crit_diameter_est = d + contains + + !> Evaluates the equilibrium residual f(d) and its derivative df(d) + !> together with the effective surface tension eff(d) at wet diameter + !> d, for the host particle. valid is .false. in the degenerate + !> regime where an organic film cannot even coat the insoluble core, + !> in which case the particle is treated as non-activating. + subroutine est_residual(d, f, df, eff, valid) + + !> Wet diameter (m). + real(kind=dp), intent(in) :: d + !> Equilibrium residual (zero at every extremum of S(d)). + real(kind=dp), intent(out) :: f + !> Derivative of the residual with respect to d. + real(kind=dp), intent(out) :: df + !> Effective surface tension at d (J m^-2). + real(kind=dp), intent(out) :: eff + !> Whether d is in the activating regime. + logical, intent(out) :: valid + + real(kind=dp) :: surf_eng_soluble, c_5, v_sol, v_water, r_core, v_delta + + v_sol = sphere_diam2vol(d) - v_solid - v_org + if (v_sol <= 0d0) then + ! Diameter below the dry soluble volume: not a physical wet state, and + ! below the activation barrier. Flag invalid (f < 0 so the bracketed + ! solver moves the search up out of it). + eff = const%water_surf_eng + f = -1d0 + df = 0d0 + valid = .false. + return + end if + valid = .true. + ! Volume-weighted surface tension of the soluble core, evaluated from the + ! per-particle sums precomputed in the caller (no per-call species lookup). + v_water = v_sol - soluble_volume + surf_eng_soluble = (sum_sol_sigma + v_water * const%water_surf_eng) / v_sol + c_5 = v_org * (surf_eng_organic - surf_eng_soluble) + if (v_solid == 0d0) then + r_core = sphere_vol2rad(v_sol) + v_delta = sphere_rad2vol(r_core + delta_min) - sphere_rad2vol(r_core) + call est_poly(d, surf_eng_soluble, c_5, v_delta, f, df, eff) + else + if (v_sol + v_org > v_delta_solid) then + r_core = sphere_vol2rad(v_sol + v_solid) + v_delta = sphere_rad2vol(r_core + delta_min) & + - sphere_rad2vol(r_core) + call est_poly(d, surf_eng_soluble, c_5, v_delta, f, df, eff) + else + ! Film cannot form a delta-shell around the solid core: this + ! non-activating regime lies below the activation barrier, so + ! we return f < 0 (df = 0 forces the safeguarded solver to take + ! a bisection step that moves the search up out of it). + eff = (v_org * surf_eng_organic + v_sol * surf_eng_soluble) & + / v_delta_solid + f = -1d0 + df = 0d0 + valid = .false. + end if + end if + + end subroutine est_residual + + !> Evaluates the equilibrium residual, its derivative and the + !> effective surface tension in the shell-forming regime, given the + !> soluble surface tension and delta-shell volume at diameter d. + subroutine est_poly(d, surf_eng_soluble, c_5, v_delta, f, df, eff) + + real(kind=dp), intent(in) :: d, surf_eng_soluble, c_5, v_delta + real(kind=dp), intent(out) :: f, df, eff + + real(kind=dp) :: d_v_delta, d_eff_surf_eng, d_2_eff_surf_eng, R, d_R + + if (v_org > v_delta) then + eff = surf_eng_organic + f = d**6 - c_1 * d**4 / eff - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff - 3d0 * c_2 * d**2 + else if (v_org == 0d0) then + eff = surf_eng_soluble + f = d**6 - c_1 * d**4 / eff - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / eff - 3d0 * c_2 * d**2 + else + eff = surf_eng_soluble + c_5 / v_delta + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_eff_surf_eng = - c_5 * d_v_delta / v_delta**2 + R = eff - d * d_eff_surf_eng + d_2_eff_surf_eng = (c_5 / v_delta**3) & + * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_eff_surf_eng + f = R * (d**6 - c_2 * d**3 + c_3) - c_1 * d**4 + df = d_R * (d**6 - c_2 * d**3 + c_3) & + + R * (6d0 * d**5 - 3d0 * c_2 * d**2) - 4d0 * c_1 * d**3 + end if + + end subroutine est_poly + end function aero_particle_crit_diameter_est !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1095,17 +1219,21 @@ real(kind=dp) function aero_particle_organic_volume(aero_particle, aero_data) type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - integer :: i_org_spec, n_org_spec + integer :: i_org_spec, i_spec character(len=AERO_NAME_LEN), parameter, dimension(10) :: & org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] aero_particle_organic_volume = 0d0 - do n_org_spec = 1, size(org_spec) - i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) - aero_particle_organic_volume = aero_particle_organic_volume & - + aero_particle%vol(i_org_spec) + do i_org_spec = 1, size(org_spec) + i_spec = aero_data_spec_by_name(aero_data, org_spec(i_org_spec)) + ! aero_data_spec_by_name returns 0 for a name absent from aero_data; + ! skip it rather than indexing vol(0) out of bounds. + if (i_spec > 0) then + aero_particle_organic_volume = aero_particle_organic_volume & + + aero_particle%vol(i_spec) + end if end do end function aero_particle_organic_volume @@ -1118,67 +1246,59 @@ real(kind=dp) function aero_particle_solid_volume(aero_particle, aero_data) type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - integer :: i_solid_spec, n_solid_spec + integer :: i_solid_spec, i_spec character(len=AERO_NAME_LEN), parameter, dimension(2) :: & solid_spec = ["OIN ", "BC "] aero_particle_solid_volume = 0d0 - do n_solid_spec = 1, size(solid_spec) - i_solid_spec = aero_data_spec_by_name(aero_data, solid_spec(n_solid_spec)) - aero_particle_solid_volume = aero_particle_solid_volume & - + aero_particle%vol(i_solid_spec) + do i_solid_spec = 1, size(solid_spec) + i_spec = aero_data_spec_by_name(aero_data, solid_spec(i_solid_spec)) + ! skip names absent from aero_data (spec_by_name returns 0) + if (i_spec > 0) then + aero_particle_solid_volume = aero_particle_solid_volume & + + aero_particle%vol(i_spec) + end if end do end function aero_particle_solid_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the volume-weighted surface tension of the soluble core - !> (soluble inorganics plus water) at diameter d (J m^-2). - real(kind=dp) function aero_particle_surf_eng_soluble(aero_particle, & - aero_data, env_state, d) + !> Returns two diameter-independent per-particle sums over the soluble + !> inorganic species: the total soluble volume and the volume-weighted + !> surface-tension sum \f$\sum_i V_i \sigma_i\f$. The caller resolves these + !> once and then forms the soluble-core surface tension at any wet diameter + !> by cheap arithmetic, avoiding per-diameter species-name lookups. + subroutine aero_particle_soluble_sums(aero_particle, aero_data, & + soluble_volume, sum_sol_sigma) !> Aerosol particle. type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - !> Environment state. - type(env_state_t), intent(in) :: env_state - !> Diameter at which to evaluate the surface tension (m). - real(kind=dp) :: d + !> Total soluble inorganic volume (m^3). + real(kind=dp), intent(out) :: soluble_volume + !> Volume-weighted surface-tension sum \f$\sum_i V_i \sigma_i\f$ (m J m^-2). + real(kind=dp), intent(out) :: sum_sol_sigma - integer :: i_sol_spec, n_sol_spec - real(kind=dp) :: aero_particle_soluble_volume - real(kind=dp) :: v_solid, v_org, v_sol, v_water + integer :: i_sol_spec, i_spec character(len=AERO_NAME_LEN), parameter, dimension(7) :: & sol_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & "Na ", "Ca "] - v_solid = aero_particle_solid_volume(aero_particle, aero_data) - v_org = aero_particle_organic_volume(aero_particle, aero_data) - - aero_particle_soluble_volume = 0d0 - - do n_sol_spec = 1, size(sol_spec) - i_sol_spec = aero_data_spec_by_name(aero_data, sol_spec(n_sol_spec)) - aero_particle_soluble_volume = aero_particle_soluble_volume & - + aero_particle%vol(i_sol_spec) - end do - - v_sol = d**3 * const%pi / 6 - v_solid - v_org - v_water = v_sol - aero_particle_soluble_volume - - aero_particle_surf_eng_soluble = 0d0 - do n_sol_spec = 1, size(sol_spec) - i_sol_spec = aero_data_spec_by_name(aero_data, sol_spec(n_sol_spec)) - aero_particle_surf_eng_soluble = aero_particle_surf_eng_soluble & - + aero_particle%vol(i_sol_spec) * & - aero_data%sigma(i_sol_spec) / v_sol + soluble_volume = 0d0 + sum_sol_sigma = 0d0 + do i_sol_spec = 1, size(sol_spec) + i_spec = aero_data_spec_by_name(aero_data, sol_spec(i_sol_spec)) + ! skip names absent from aero_data (spec_by_name returns 0) + if (i_spec > 0) then + soluble_volume = soluble_volume + aero_particle%vol(i_spec) + sum_sol_sigma = sum_sol_sigma & + + aero_particle%vol(i_spec) * aero_data%sigma(i_spec) + end if end do - aero_particle_surf_eng_soluble = aero_particle_surf_eng_soluble + v_water & - * const%water_surf_eng / v_sol - end function aero_particle_surf_eng_soluble + end subroutine aero_particle_soluble_sums !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1192,7 +1312,7 @@ real(kind=dp) function aero_particle_surf_eng_organic(aero_particle, aero_data) type(aero_data_t), intent(in) :: aero_data real(kind=dp) :: org_volume - integer :: i_org_spec, n_org_spec + integer :: i_org_spec, i_spec character(len=AERO_NAME_LEN), parameter, dimension(10) :: & org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] @@ -1200,11 +1320,20 @@ real(kind=dp) function aero_particle_surf_eng_organic(aero_particle, aero_data) org_volume = aero_particle_organic_volume(aero_particle, aero_data) aero_particle_surf_eng_organic = 0d0 - do n_org_spec = 1, size(org_spec) - i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) - aero_particle_surf_eng_organic = aero_particle_surf_eng_organic & - + aero_particle%vol(i_org_spec) * & - aero_data%sigma(i_org_spec) / org_volume + if (org_volume <= 0d0) then + ! No organic material: the film surface tension is undefined. Return + ! zero to avoid a 0/0; the caller only uses this through + ! c_5 = v_org * (sigma_org - sigma_sol), which is zero when v_org = 0. + return + end if + + do i_org_spec = 1, size(org_spec) + i_spec = aero_data_spec_by_name(aero_data, org_spec(i_org_spec)) + if (i_spec > 0) then + aero_particle_surf_eng_organic = aero_particle_surf_eng_organic & + + aero_particle%vol(i_spec) * & + aero_data%sigma(i_spec) / org_volume + end if end do end function aero_particle_surf_eng_organic diff --git a/src/util.F90 b/src/util.F90 index f6ded7fa5..d5ff945db 100644 --- a/src/util.F90 +++ b/src/util.F90 @@ -283,6 +283,19 @@ real(kind=dp) elemental function diam2rad(d) end function diam2rad +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Convert diameter (m) to mass-equivalent volume \f$V\f$ (m^3) for + !> spherical particles. + real(kind=dp) elemental function sphere_diam2vol(d) + + !> Diameter (m). + real(kind=dp), intent(in) :: d + + sphere_diam2vol = const%pi / 6d0 * d**3 + + end function sphere_diam2vol + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Calculate air molecular mean free path \f$l\f$ (m).