From 5d85b707e71fe69a9a88d8f44c86af7305f47d00 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 15 Apr 2022 22:13:25 -0500 Subject: [PATCH 01/33] add subroutine for single particle --- src/output.F90 | 49 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index b3e33b2d0..627dc958e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1,4 +1,4 @@ -! Copyright (C) 2005-2022 Nicole Riemer and Matthew West +! Copyright (C) 2005-2018 Nicole Riemer and Matthew West ! Licensed under the GNU General Public License version 2 or (at your ! option) any later version. See the file COPYING for details. @@ -82,7 +82,7 @@ module pmc_output #endif !> PartMC verson number. - character(len=100), parameter :: PARTMC_VERSION = "2.6.1" + character(len=100), parameter :: PARTMC_VERSION = "2.5.0" !> Type code for undefined or invalid output. integer, parameter :: OUTPUT_TYPE_INVALID = 0 @@ -247,6 +247,51 @@ subroutine make_filename(filename, prefix, suffix, index, i_repeat, & end subroutine make_filename +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Modified by xx24, Feb/12/2022 + !> Make a filename from a given prefix and other information. + subroutine for_single_particle(filename, prefix_single_particle, suffix, index, i_repeat, & + write_rank, write_n_proc) + +!> Filename to create. +character(len=*), intent(out) :: filename +!> Filename prefix. +character(len=*), intent(in) :: prefix_single_particle +!> Filename suffix. +character(len=*), intent(in) :: suffix +!> Filename index. +integer, intent(in), optional :: index +!> Current repeat number. +integer, intent(in), optional :: i_repeat +!> Rank to write into file. +integer, intent(in), optional :: write_rank +!> Number of processes to write into file. +integer, intent(in), optional :: write_n_proc + +integer :: ncid, use_rank, use_n_proc +character(len=100) :: proc_string, index_string, repeat_string + +if (present(write_rank)) then + use_rank = write_rank +else + use_rank = pmc_mpi_rank() +end if +if (present(write_n_proc)) then + use_n_proc = write_n_proc +else + use_n_proc = pmc_mpi_size() +end if + +repeat_string = "" +proc_string = "" +index_string = "" +if (present(i_repeat)) write(repeat_string, '(a,i4.4)') '_', i_repeat +if (use_n_proc > 1) write(proc_string, '(a,i4.4)') '_', (use_rank + 1) +if (present(index)) write(index_string, '(a,i8.8)') '_', index +write(filename, '(a,a,a,a,a)') trim(prefix_single_particle), trim(repeat_string), & + trim(proc_string), trim(index_string), trim(suffix) + +end subroutine for_single_particle !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Helper routine to write time variables. Do not call directly. From d04a6df396a45a1a550cdb4345626f9c531f8a92 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 15 Apr 2022 22:15:14 -0500 Subject: [PATCH 02/33] add detailed species --- .../1_urban_plume/urban_plume_process.F90 | 518 +++++++++++++++++- 1 file changed, 509 insertions(+), 9 deletions(-) diff --git a/scenarios/1_urban_plume/urban_plume_process.F90 b/scenarios/1_urban_plume/urban_plume_process.F90 index 6b5e9f03e..818844673 100644 --- a/scenarios/1_urban_plume/urban_plume_process.F90 +++ b/scenarios/1_urban_plume/urban_plume_process.F90 @@ -10,26 +10,93 @@ program process use pmc_output use pmc_stats + use pmc_aero_state + use pmc_aero_data + use pmc_env_state + use pmc_aero_particle + use pmc_aero_particle_array character(len=PMC_MAX_FILENAME_LEN), parameter :: prefix & = "out/urban_plume" - - character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename - type(bin_grid_t) :: diam_grid, bc_grid, sc_grid + character(len=PMC_MAX_FILENAME_LEN), parameter :: prefix_new & + = "out/single_particle" + + character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename, & + single_particle_output + type(bin_grid_t) :: diam_grid, bc_grid, oc_grid, sc_grid, & + so4_grid, no3_grid, nh4_grid, soa_grid, & + cl_grid, msa_grid, aro1_grid, aro2_grid, & + alk1_grid, ole1_grid, api1_grid, api2_grid, & + lim1_grid, lim2_grid, co3_grid, na_grid, & + ca_grid, oin_grid, h2o_grid type(aero_data_t) :: aero_data + type(aero_particle_t) :: aero_particle type(aero_state_t) :: aero_state type(env_state_t) :: env_state integer :: ncid, index, repeat, i_index, i_repeat, n_index, n_repeat - real(kind=dp) :: time, del_t, tot_num_conc, tot_mass_conc + integer :: number + real(kind=dp) :: crit_diam + real(kind=dp) :: time, del_t, tot_num_conc, tot_mass_conc, & + tot_dry_mass_conc, tot_bc_mass_conc, tot_oc_mass_conc, & + tot_so4_mass_conc, tot_no3_mass_conc, tot_nh4_mass_conc, & + tot_soa_mass_conc, tot_cl_mass_conc, tot_msa_mass_conc, & + tot_aro1_mass_conc, tot_aro2_mass_conc, tot_alk1_mass_conc, & + tot_ole1_mass_conc, tot_api1_mass_conc, tot_api2_mass_conc, & + tot_lim1_mass_conc, tot_lim2_mass_conc, tot_co3_mass_conc, & + tot_na_mass_conc, tot_ca_mass_conc, tot_oin_mass_conc, & + tot_h2o_mass_conc real(kind=dp) :: d_alpha, d_gamma, chi character(len=PMC_UUID_LEN) :: uuid 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(:,:), oc_masses(:), oc_fracs(:), diam_oc_dist(:,:), & + crit_rhs(:), scs(:), num_dist(:), diam_sc_dist(:,:), & + so4_masses(:), so4_fracs(:), diam_so4_dist(:,:), & + no3_masses(:), no3_fracs(:), diam_no3_dist(:,:), & + nh4_masses(:), nh4_fracs(:), diam_nh4_dist(:,:), & + soa_masses(:), soa_fracs(:), diam_soa_dist(:,:), & + cl_masses(:), cl_fracs(:), diam_cl_dist(:,:), & + msa_masses(:), msa_fracs(:), diam_msa_dist(:,:), & + aro1_masses(:), aro1_fracs(:), diam_aro1_dist(:,:), & + aro2_masses(:), aro2_fracs(:), diam_aro2_dist(:,:), & + alk1_masses(:), alk1_fracs(:), diam_alk1_dist(:,:), & + ole1_masses(:), ole1_fracs(:), diam_ole1_dist(:,:), & + api1_masses(:), api1_fracs(:), diam_api1_dist(:,:), & + api2_masses(:), api2_fracs(:), diam_api2_dist(:,:), & + lim1_masses(:), lim1_fracs(:), diam_lim1_dist(:,:), & + lim2_masses(:), lim2_fracs(:), diam_lim2_dist(:,:), & + co3_masses(:), co3_fracs(:), diam_co3_dist(:,:), & + na_masses(:), na_fracs(:), diam_na_dist(:,:), & + ca_masses(:), ca_fracs(:), diam_ca_dist(:,:), & + oin_masses(:), oin_fracs(:), diam_oin_dist(:,:), & + h2o_masses(:), h2o_fracs(:), diam_h2o_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 + stats_tot_mass_conc, stats_d_gamma, stats_chi, & + stats_tot_bc_mass_conc, stats_tot_oc_mass_conc, & + stats_tot_so4_mass_conc, stats_tot_no3_mass_conc, & + stats_tot_nh4_mass_conc, stats_tot_dry_mass_conc, & + stats_tot_soa_mass_conc, stats_tot_cl_mass_conc, & + stats_tot_msa_mass_conc, stats_tot_aro1_mass_conc, & + stats_tot_aro2_mass_conc, stats_tot_alk1_mass_conc, & + stats_tot_ole1_mass_conc, stats_tot_api1_mass_conc, & + stats_tot_api2_mass_conc, stats_tot_lim1_mass_conc, & + stats_tot_lim2_mass_conc, stats_tot_co3_mass_conc, & + stats_tot_na_mass_conc, stats_tot_ca_mass_conc, & + stats_tot_oin_mass_conc, stats_tot_h2o_mass_conc + + type(stats_2d_t) :: stats_diam_sc_dist, stats_diam_bc_dist, & + stats_diam_oc_dist, stats_diam_so4_dist, & + stats_diam_no3_dist, stats_diam_nh4_dist, & + stats_diam_soa_dist, stats_diam_cl_dist, & + stats_diam_msa_dist, stats_diam_aro1_dist, & + stats_diam_aro2_dist, stats_diam_alk1_dist, & + stats_diam_ole1_dist, stats_diam_api1_dist, & + stats_diam_api2_dist, stats_diam_lim1_dist, & + stats_diam_lim2_dist, stats_diam_co3_dist, & + stats_diam_na_dist, stats_diam_oin_dist, & + stats_diam_ca_dist, stats_diam_h2o_dist + character(len=AERO_NAME_LEN), allocatable :: mixing_state_groups(:,:) call pmc_mpi_init() @@ -38,6 +105,26 @@ program process call bin_grid_make(diam_grid, BIN_GRID_TYPE_LOG, 180, 1d-9, 1d-3) call bin_grid_make(bc_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(oc_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(so4_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(no3_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(nh4_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(soa_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(cl_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(msa_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(aro1_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(aro2_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(alk1_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(ole1_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(api1_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(api2_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(lim1_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(lim2_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(co3_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(na_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(ca_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(oin_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) + call bin_grid_make(h2o_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) call bin_grid_make(sc_grid, BIN_GRID_TYPE_LOG, 50, 1d-4, 1d0) allocate(times(n_index)) @@ -47,8 +134,28 @@ program process mixing_state_groups(2,:) = ["API1 ", "API2 ", "LIM1 ", "LIM2 "] mixing_state_groups(3,:) = ["SO4 ", "NO3 ", "NH4 ", " "] + scs = [ real(kind=dp) :: ] ! silence compiler warnings bc_fracs = [ real(kind=dp) :: ] + oc_fracs = [ real(kind=dp) :: ] + so4_fracs = [ real(kind=dp) :: ] + no3_fracs = [ real(kind=dp) :: ] + nh4_fracs = [ real(kind=dp) :: ] + cl_fracs = [ real(kind=dp) :: ] + msa_fracs = [ real(kind=dp) :: ] + aro1_fracs = [ real(kind=dp) :: ] + aro2_fracs = [ real(kind=dp) :: ] + alk1_fracs = [ real(kind=dp) :: ] + ole1_fracs = [ real(kind=dp) :: ] + api1_fracs = [ real(kind=dp) :: ] + api2_fracs = [ real(kind=dp) :: ] + lim1_fracs = [ real(kind=dp) :: ] + lim2_fracs = [ real(kind=dp) :: ] + co3_fracs = [ real(kind=dp) :: ] + na_fracs = [ real(kind=dp) :: ] + ca_fracs = [ real(kind=dp) :: ] + oin_fracs = [ real(kind=dp) :: ] + h2o_fracs = [ real(kind=dp) :: ] do i_index = 1,n_index do i_repeat = 1,n_repeat @@ -73,19 +180,268 @@ program process masses = aero_state_masses(aero_state, aero_data) tot_mass_conc = sum(masses * num_concs) call stats_1d_add_entry(stats_tot_mass_conc, tot_mass_conc, i_index) - +!for dry mass dry_masses = aero_state_masses(aero_state, aero_data, & exclude=(/"H2O"/)) + tot_dry_mass_conc = sum(dry_masses * num_concs) + call stats_1d_add_entry(stats_tot_dry_mass_conc, tot_dry_mass_conc, i_index) +!for bc bc_masses = aero_state_masses(aero_state, aero_data, & include=(/"BC"/)) + tot_bc_mass_conc = sum(bc_masses * num_concs) + call stats_1d_add_entry(stats_tot_bc_mass_conc, tot_bc_mass_conc, i_index) + bc_fracs = bc_masses / dry_masses diam_bc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & bc_grid, bc_fracs, num_concs) call stats_2d_add(stats_diam_bc_dist, diam_bc_dist) +!for OC + oc_masses = aero_state_masses(aero_state, aero_data, & + include=(/"OC"/)) + tot_oc_mass_conc = sum(oc_masses * num_concs) + call stats_1d_add_entry(stats_tot_oc_mass_conc, tot_oc_mass_conc, i_index) + + oc_fracs = oc_masses / dry_masses + diam_oc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + oc_grid, oc_fracs, num_concs) + call stats_2d_add(stats_diam_oc_dist, diam_oc_dist) +!for SO4 + so4_masses = aero_state_masses(aero_state, aero_data, & + include=(/"SO4"/)) + tot_so4_mass_conc = sum(so4_masses * num_concs) + call stats_1d_add_entry(stats_tot_so4_mass_conc, tot_so4_mass_conc, i_index) + + so4_fracs = so4_masses / dry_masses + diam_so4_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + so4_grid, so4_fracs, num_concs) + call stats_2d_add(stats_diam_so4_dist, diam_so4_dist) +!for NO3 + no3_masses = aero_state_masses(aero_state, aero_data, & + include=(/"NO3"/)) + tot_no3_mass_conc = sum(no3_masses * num_concs) + call stats_1d_add_entry(stats_tot_no3_mass_conc, tot_no3_mass_conc, i_index) + + no3_fracs = no3_masses / dry_masses + diam_no3_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + no3_grid, no3_fracs, num_concs) + call stats_2d_add(stats_diam_no3_dist, diam_no3_dist) +!for NH4 + nh4_masses = aero_state_masses(aero_state, aero_data, & + include=(/"NH4"/)) + tot_nh4_mass_conc = sum(nh4_masses * num_concs) + call stats_1d_add_entry(stats_tot_nh4_mass_conc, tot_nh4_mass_conc, i_index) + + nh4_fracs = nh4_masses / dry_masses + diam_nh4_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + nh4_grid, nh4_fracs, num_concs) + call stats_2d_add(stats_diam_nh4_dist, diam_nh4_dist) +!for SOA + soa_masses = aero_state_masses(aero_state, aero_data, & + include=(/"ARO1", "ARO2", "ALK1", "OLE1", "API1", "API2", "LIM1", "LIM2"/)) + tot_soa_mass_conc = sum(soa_masses * num_concs) + call stats_1d_add_entry(stats_tot_soa_mass_conc, tot_soa_mass_conc, i_index) + + soa_fracs = soa_masses / dry_masses + diam_soa_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + soa_grid, soa_fracs, num_concs) + call stats_2d_add(stats_diam_soa_dist, diam_soa_dist) +!for Cl + cl_masses = aero_state_masses(aero_state, aero_data, & + include=(/"Cl"/)) + tot_cl_mass_conc = sum(cl_masses * num_concs) + call stats_1d_add_entry(stats_tot_cl_mass_conc, tot_cl_mass_conc, i_index) + + cl_fracs = cl_masses / dry_masses + diam_cl_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + cl_grid, cl_fracs, num_concs) + call stats_2d_add(stats_diam_cl_dist, diam_cl_dist) +!for MSA + msa_masses = aero_state_masses(aero_state, aero_data, & + include=(/"MSA"/)) + tot_msa_mass_conc = sum(msa_masses * num_concs) + call stats_1d_add_entry(stats_tot_msa_mass_conc, tot_msa_mass_conc, i_index) + + msa_fracs = msa_masses / dry_masses + diam_msa_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + msa_grid, msa_fracs, num_concs) + call stats_2d_add(stats_diam_msa_dist, diam_msa_dist) +!for aro1 + aro1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"ARO1"/)) + tot_aro1_mass_conc = sum(aro1_masses * num_concs) + call stats_1d_add_entry(stats_tot_aro1_mass_conc, tot_aro1_mass_conc, i_index) + + aro1_fracs = aro1_masses / dry_masses + diam_aro1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + aro1_grid, aro1_fracs, num_concs) + call stats_2d_add(stats_diam_aro1_dist, diam_aro1_dist) +!for aro2 + aro2_masses = aero_state_masses(aero_state, aero_data, & + include=(/"ARO2"/)) + tot_aro2_mass_conc = sum(aro2_masses * num_concs) + call stats_1d_add_entry(stats_tot_aro2_mass_conc, tot_aro2_mass_conc, i_index) + + aro2_fracs = aro2_masses / dry_masses + diam_aro2_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + aro2_grid, aro2_fracs, num_concs) + call stats_2d_add(stats_diam_aro2_dist, diam_aro2_dist) +!for alk1 + alk1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"ALK1"/)) + tot_alk1_mass_conc = sum(alk1_masses * num_concs) + call stats_1d_add_entry(stats_tot_alk1_mass_conc, tot_alk1_mass_conc, i_index) + + alk1_fracs = alk1_masses / dry_masses + diam_alk1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + alk1_grid, alk1_fracs, num_concs) + call stats_2d_add(stats_diam_alk1_dist, diam_alk1_dist) +!for ole1 + ole1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"OLE1"/)) + tot_ole1_mass_conc = sum(ole1_masses * num_concs) + call stats_1d_add_entry(stats_tot_ole1_mass_conc, tot_ole1_mass_conc, i_index) + + ole1_fracs = ole1_masses / dry_masses + diam_ole1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + ole1_grid, ole1_fracs, num_concs) + call stats_2d_add(stats_diam_ole1_dist, diam_ole1_dist) +!for ole1 + ole1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"OLE1"/)) + tot_ole1_mass_conc = sum(ole1_masses * num_concs) + call stats_1d_add_entry(stats_tot_ole1_mass_conc, tot_ole1_mass_conc, i_index) + + ole1_fracs = ole1_masses / dry_masses + diam_ole1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + ole1_grid, ole1_fracs, num_concs) + call stats_2d_add(stats_diam_ole1_dist, diam_ole1_dist) +!for api1 + api1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"API1"/)) + tot_api1_mass_conc = sum(api1_masses * num_concs) + call stats_1d_add_entry(stats_tot_api1_mass_conc, tot_api1_mass_conc, i_index) + + api1_fracs = api1_masses / dry_masses + diam_api1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + api1_grid, api1_fracs, num_concs) + call stats_2d_add(stats_diam_api1_dist, diam_api1_dist) +!for api2 + api2_masses = aero_state_masses(aero_state, aero_data, & + include=(/"API2"/)) + tot_api2_mass_conc = sum(api2_masses * num_concs) + call stats_1d_add_entry(stats_tot_api2_mass_conc, tot_api2_mass_conc, i_index) + + api2_fracs = api2_masses / dry_masses + diam_api2_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + api2_grid, api2_fracs, num_concs) + call stats_2d_add(stats_diam_api2_dist, diam_api2_dist) +!for lim1 + lim1_masses = aero_state_masses(aero_state, aero_data, & + include=(/"LIM1"/)) + tot_lim1_mass_conc = sum(lim1_masses * num_concs) + call stats_1d_add_entry(stats_tot_lim1_mass_conc, tot_lim1_mass_conc, i_index) + + lim1_fracs = lim1_masses / dry_masses + diam_lim1_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + lim1_grid, lim1_fracs, num_concs) + call stats_2d_add(stats_diam_lim1_dist, diam_lim1_dist) +!for lim2 + lim2_masses = aero_state_masses(aero_state, aero_data, & + include=(/"LIM2"/)) + tot_lim2_mass_conc = sum(lim2_masses * num_concs) + call stats_1d_add_entry(stats_tot_lim2_mass_conc, tot_lim2_mass_conc, i_index) + + lim2_fracs = lim2_masses / dry_masses + diam_lim2_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + lim2_grid, lim2_fracs, num_concs) + call stats_2d_add(stats_diam_lim2_dist, diam_lim2_dist) +!for co3 + co3_masses = aero_state_masses(aero_state, aero_data, & + include=(/"CO3"/)) + tot_co3_mass_conc = sum(co3_masses * num_concs) + call stats_1d_add_entry(stats_tot_co3_mass_conc, tot_co3_mass_conc, i_index) + + co3_fracs = co3_masses / dry_masses + diam_co3_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + co3_grid, co3_fracs, num_concs) + call stats_2d_add(stats_diam_co3_dist, diam_co3_dist) +!for na + na_masses = aero_state_masses(aero_state, aero_data, & + include=(/"Na"/)) + tot_na_mass_conc = sum(na_masses * num_concs) + call stats_1d_add_entry(stats_tot_na_mass_conc, tot_na_mass_conc, i_index) + + na_fracs = na_masses / dry_masses + diam_na_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + na_grid, na_fracs, num_concs) + call stats_2d_add(stats_diam_na_dist, diam_na_dist) +!for ca + ca_masses = aero_state_masses(aero_state, aero_data, & + include=(/"Ca"/)) + tot_ca_mass_conc = sum(ca_masses * num_concs) + call stats_1d_add_entry(stats_tot_ca_mass_conc, tot_ca_mass_conc, i_index) + + ca_fracs = ca_masses / dry_masses + diam_ca_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + ca_grid, ca_fracs, num_concs) + call stats_2d_add(stats_diam_ca_dist, diam_ca_dist) +!for oin + oin_masses = aero_state_masses(aero_state, aero_data, & + include=(/"OIN"/)) + tot_oin_mass_conc = sum(oin_masses * num_concs) + call stats_1d_add_entry(stats_tot_oin_mass_conc, tot_oin_mass_conc, i_index) + + oin_fracs = oin_masses / dry_masses + diam_oin_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + oin_grid, oin_fracs, num_concs) + call stats_2d_add(stats_diam_oin_dist, diam_oin_dist) +!for h2o + h2o_masses = aero_state_masses(aero_state, aero_data, & + include=(/"H2O"/)) + tot_h2o_mass_conc = sum(h2o_masses * num_concs) + call stats_1d_add_entry(stats_tot_h2o_mass_conc, tot_h2o_mass_conc, i_index) + + h2o_fracs = h2o_masses / dry_masses + diam_h2o_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + h2o_grid, h2o_fracs, num_concs) + call stats_2d_add(stats_diam_h2o_dist, diam_h2o_dist) +!for critical RH supersaturation crit_rhs = aero_state_crit_rel_humids(aero_state, aero_data, & env_state) scs = crit_rhs - 1d0 + + call for_single_particle(single_particle_output, prefix_new, ".csv", index) + open(15, file= trim(single_particle_output)) + write(15, *) 'dry_diameters', dry_diameters + write(15, *) 'masses', masses + write(15, *) 'dry_masses', dry_masses + write(15, *) 'bc_masses', bc_masses + write(15, *) 'oc_masses', oc_masses + write(15, *) 'so4_masses', so4_masses + write(15, *) 'no3_masses', no3_masses + write(15, *) 'nh4_masses', nh4_masses + write(15, *) 'soa_masses', soa_masses + write(15, *) 'cl_masses', cl_masses + write(15, *) 'msa_masses', msa_masses + write(15, *) 'aro1_masses', aro1_masses + write(15, *) 'aro2_masses', aro2_masses + write(15, *) 'alk1_masses', alk1_masses + write(15, *) 'ole1_masses', ole1_masses + write(15, *) 'api1_masses', api1_masses + write(15, *) 'api2_masses', api2_masses + write(15, *) 'lim1_masses', lim1_masses + write(15, *) 'lim2_masses', lim2_masses + write(15, *) 'co3_masses', co3_masses + write(15, *) 'na_masses', na_masses + write(15, *) 'ca_masses', ca_masses + write(15, *) 'oin_masses', oin_masses + write(15, *) 'h2o_masses', h2o_masses + write(15, *) 'crit_rhs', crit_rhs + write(15, *) 'scs', scs + write(15, *) 'num_con', num_concs + close(15) + diam_sc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & sc_grid, scs, num_concs) call stats_2d_add(stats_diam_sc_dist, diam_sc_dist) @@ -105,6 +461,26 @@ program process call pmc_nc_write_info(ncid, uuid, "1_urban_plume process") call bin_grid_output_netcdf(diam_grid, ncid, "diam", unit="m") call bin_grid_output_netcdf(bc_grid, ncid, "bc_frac", unit="1") + call bin_grid_output_netcdf(oc_grid, ncid, "oc_frac", unit="1") + call bin_grid_output_netcdf(so4_grid, ncid, "so4_frac", unit="1") + call bin_grid_output_netcdf(no3_grid, ncid, "no3_frac", unit="1") + call bin_grid_output_netcdf(nh4_grid, ncid, "nh4_frac", unit="1") + call bin_grid_output_netcdf(soa_grid, ncid, "soa_frac", unit="1") + call bin_grid_output_netcdf(cl_grid, ncid, "cl_frac", unit="1") + call bin_grid_output_netcdf(msa_grid, ncid, "msa_frac", unit="1") + call bin_grid_output_netcdf(aro1_grid, ncid, "aro1_frac", unit="1") + call bin_grid_output_netcdf(aro2_grid, ncid, "aro2_frac", unit="1") + call bin_grid_output_netcdf(alk1_grid, ncid, "alk1_frac", unit="1") + call bin_grid_output_netcdf(ole1_grid, ncid, "ole1_frac", unit="1") + call bin_grid_output_netcdf(api1_grid, ncid, "api1_frac", unit="1") + call bin_grid_output_netcdf(api2_grid, ncid, "api2_frac", unit="1") + call bin_grid_output_netcdf(lim1_grid, ncid, "lim1_frac", unit="1") + call bin_grid_output_netcdf(lim2_grid, ncid, "lim2_frac", unit="1") + call bin_grid_output_netcdf(co3_grid, ncid, "co3_frac", unit="1") + call bin_grid_output_netcdf(na_grid, ncid, "na_frac", unit="1") + call bin_grid_output_netcdf(ca_grid, ncid, "ca_frac", unit="1") + call bin_grid_output_netcdf(oin_grid, ncid, "oin_frac", unit="1") + call bin_grid_output_netcdf(h2o_grid, ncid, "h2o_frac", unit="1") call bin_grid_output_netcdf(sc_grid, ncid, "sc", unit="1") call stats_1d_output_netcdf(stats_num_dist, ncid, "num_dist", & @@ -115,6 +491,86 @@ program process dim_name_1="diam", dim_name_2="bc_frac", unit="m^{-3}") call stats_2d_clear(stats_diam_bc_dist) + call stats_2d_output_netcdf(stats_diam_oc_dist, ncid, "diam_oc_dist", & + dim_name_1="diam", dim_name_2="oc_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_oc_dist) + + call stats_2d_output_netcdf(stats_diam_so4_dist, ncid, "diam_so4_dist", & + dim_name_1="diam", dim_name_2="so4_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_so4_dist) + + call stats_2d_output_netcdf(stats_diam_no3_dist, ncid, "diam_no3_dist", & + dim_name_1="diam", dim_name_2="no3_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_no3_dist) + + call stats_2d_output_netcdf(stats_diam_nh4_dist, ncid, "diam_nh4_dist", & + dim_name_1="diam", dim_name_2="nh4_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_nh4_dist) + + call stats_2d_output_netcdf(stats_diam_soa_dist, ncid, "diam_soa_dist", & + dim_name_1="diam", dim_name_2="soa_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_soa_dist) + + call stats_2d_output_netcdf(stats_diam_cl_dist, ncid, "diam_cl_dist", & + dim_name_1="diam", dim_name_2="cl_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_cl_dist) + + call stats_2d_output_netcdf(stats_diam_msa_dist, ncid, "diam_msa_dist", & + dim_name_1="diam", dim_name_2="msa_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_cl_dist) + + call stats_2d_output_netcdf(stats_diam_aro1_dist, ncid, "diam_aro1_dist", & + dim_name_1="diam", dim_name_2="aro1_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_aro1_dist) + + call stats_2d_output_netcdf(stats_diam_aro2_dist, ncid, "diam_aro2_dist", & + dim_name_1="diam", dim_name_2="aro2_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_aro2_dist) + + call stats_2d_output_netcdf(stats_diam_alk1_dist, ncid, "diam_alk1_dist", & + dim_name_1="diam", dim_name_2="alk1_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_alk1_dist) + + call stats_2d_output_netcdf(stats_diam_ole1_dist, ncid, "diam_ole1_dist", & + dim_name_1="diam", dim_name_2="ole1_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_ole1_dist) + + call stats_2d_output_netcdf(stats_diam_api1_dist, ncid, "diam_api1_dist", & + dim_name_1="diam", dim_name_2="api1_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_api1_dist) + + call stats_2d_output_netcdf(stats_diam_api2_dist, ncid, "diam_api2_dist", & + dim_name_1="diam", dim_name_2="api2_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_api2_dist) + + call stats_2d_output_netcdf(stats_diam_lim1_dist, ncid, "diam_lim1_dist", & + dim_name_1="diam", dim_name_2="lim1_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_lim1_dist) + + call stats_2d_output_netcdf(stats_diam_lim2_dist, ncid, "diam_lim2_dist", & + dim_name_1="diam", dim_name_2="lim2_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_lim2_dist) + + call stats_2d_output_netcdf(stats_diam_co3_dist, ncid, "diam_co3_dist", & + dim_name_1="diam", dim_name_2="co3_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_co3_dist) + + call stats_2d_output_netcdf(stats_diam_na_dist, ncid, "diam_na_dist", & + dim_name_1="diam", dim_name_2="na_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_na_dist) + + call stats_2d_output_netcdf(stats_diam_ca_dist, ncid, "diam_ca_dist", & + dim_name_1="diam", dim_name_2="ca_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_ca_dist) + + call stats_2d_output_netcdf(stats_diam_oin_dist, ncid, "diam_oin_dist", & + dim_name_1="diam", dim_name_2="oin_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_oin_dist) + + call stats_2d_output_netcdf(stats_diam_h2o_dist, ncid, "diam_h2o_dist", & + dim_name_1="diam", dim_name_2="h2o_frac", unit="m^{-3}") + call stats_2d_clear(stats_diam_h2o_dist) + call stats_2d_output_netcdf(stats_diam_sc_dist, ncid, "diam_sc_dist", & dim_name_1="diam", dim_name_2="sc", unit="m^{-3}") call stats_2d_clear(stats_diam_sc_dist) @@ -131,6 +587,50 @@ program process dim_name="time", unit="m^{-3}") call stats_1d_output_netcdf(stats_tot_mass_conc, ncid, "tot_mass_conc", & dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_dry_mass_conc, ncid, "tot_dry_mass_conc",& + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_bc_mass_conc, ncid, "tot_bc_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_oc_mass_conc, ncid, "tot_oc_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_so4_mass_conc, ncid, "tot_so4_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_no3_mass_conc, ncid, "tot_no3_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_nh4_mass_conc, ncid, "tot_nh4_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_soa_mass_conc, ncid, "tot_soa_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_cl_mass_conc, ncid, "tot_cl_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_msa_mass_conc, ncid, "tot_msa_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_aro1_mass_conc, ncid, "tot_aro1_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_aro2_mass_conc, ncid, "tot_aro2_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_alk1_mass_conc, ncid, "tot_alk1_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_ole1_mass_conc, ncid, "tot_ole1_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_api1_mass_conc, ncid, "tot_api1_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_api2_mass_conc, ncid, "tot_api2_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_lim1_mass_conc, ncid, "tot_lim1_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_lim2_mass_conc, ncid, "tot_lim2_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_co3_mass_conc, ncid, "tot_co3_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_na_mass_conc, ncid, "tot_na_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_ca_mass_conc, ncid, "tot_ca_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_oin_mass_conc, ncid, "tot_oin_mass_conc", & + dim_name="time", unit="kg m^{-3}") + call stats_1d_output_netcdf(stats_tot_h2o_mass_conc, ncid, "tot_h2o_mass_conc", & + dim_name="time", unit="kg m^{-3}") call stats_1d_output_netcdf(stats_d_alpha, ncid, "d_alpha", & dim_name="time", unit="1") call stats_1d_output_netcdf(stats_d_gamma, ncid, & From cd6d346484d5ae5968cd7bc222af6f1df69e99bb Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 28 Oct 2022 16:25:28 -0500 Subject: [PATCH 03/33] first commit --- src/aero_particle.F90 | 92 +++++++++++++++++++++++++++++++++++++++++++ src/aero_state.F90 | 70 +++++++++++++++++++++++++++++++- src/env_state.F90 | 22 ++++++++++- 3 files changed, 181 insertions(+), 3 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 41618c9a6..6aadd1ee5 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -285,7 +285,20 @@ elemental real(kind=dp) function aero_particle_species_volume( & aero_particle_species_volume = aero_particle%vol(i_spec) end function aero_particle_species_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Organic volume of a single species in the particle (m^3). xtxu, Oct/19/2022 + elemental real(kind=dp) function aero_particle_organic_species_volume( & + aero_particle, aero_data) + !> Particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Species number to find volume of. + integer, intent(in) :: i_spec_organics + + aero_particle_organic_species_volume = aero_particle%vol(i_spec_organics) + + end function aero_particle_organic_species_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Total dry volume of the particle (m^3). @@ -756,6 +769,85 @@ end function aero_particle_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical relative humidity (1). + real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(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) :: kappa, crit_diam, dry_diam, A_varying_sigma, sigma + !pass this sigma to env_state + + sigma = aero_particle_varying_sigma(aero_particle, aero_data) + A_varying_sigma = env_state_A_varying_sigma( sigma, env_state) + + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + crit_diam = aero_particle_crit_diameter(aero_particle, aero_data, & + env_state) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + if (kappa < 1d-30) then + aero_particle_crit_rel_humid_varying_sigma = exp(A_varying_sigma & + / crit_diam) + else + aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & + / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & + exp(A_varying_sigma / crit_diam) + end if + + end function aero_particle_crit_rel_humid_varying_sigma + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the varying sigma. + real(kind=dp) function aero_particle_varying_sigma(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) :: v_delta, coverage, frac_core_water + real(kind=dp) :: sigma_core + + !> Minimum shell thickness + real(kind=dp) :: delta_min = 1.6d-10 + !> sigma_core + real(kind=dp) :: sigma_inorganic = 0.058d0 + !> sigma_shell + real(kind=dp) :: sigma_shell = 0.03d0 + + + !> minimum shell volume, v_delta + v_delta = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & + / 3d0) * ( aero_particle_radius(aero_particle, aero_data) & + - delta_min)**3) + + !> coverage parameter + coverage = min(aero_particle_organic_species_volume(aero_particle, & + i_spec_organics) / v_delta, 1d0) + + !> fraction of water for inorganic core + frac_core_water = ((aero_particle_diameter(aero_particle, aero_data) & + - 2d0 * delta_min)**3 - aero_particle_dry_volume(aero_particle, & + aero_data)**3) / (aero_particle_diameter(aero_particle, aero_data) & + - 2d0 * delta_min)**3 + + !> surface tension of inorganic core and water + sigma_core = f_core_water * const%water_surf_eng + (1d0 - & + frac_core_water) * sigma_inorganic + + !> calculate effective surface tension value + aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & + coverage * sigma_shell + + end function aero_particle_varying_sigma + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical diameter (m). !! !! The method is as follows. We need to solve the polynomial diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 7c45482e6..3bb4c4e62 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1120,6 +1120,49 @@ function aero_state_volumes(aero_state, aero_data, include, exclude) end function aero_state_volumes +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> calculate organic species volume fraction, xtxu, Oct/17/2022 + + function aero_state_organic_volumes(aero_state, aero_data) + + !> Aerosol state. + type(aero_state_t), intent(in) :: aero_state + !> Aerosol data. + type(aero_data_t), optional, intent(in) :: aero_data + !> Species names to include in the mass. + character(len=*) :: include(:) + + !> Return volumes array (m^3). + real(kind=dp) :: aero_state_organic_volumes(aero_state_n_part(aero_state)) + + logical, allocatable :: use_species(:) + integer :: i_name_organics, i_spec_organics + + include = ["MSA ", "ARO1", "ARO2", "ALK1", "OLE1", & + "API1", "API2", "LIM1", "LIM2", "OC "] + + allocate(use_species(aero_data_n_spec(aero_data))) + + use_species = .false. + do i_name_organics = 1, size(include) + i_spec_organics = aero_data_spec_by_name(aero_data, include(i_name_organics)) + call assert_msg(111852070, i_spec_organics > 0, & + "unknown species: " // trim(include(i_name_organics))) + use_species(i_spec_organics) = .true. + end do + + aero_state_organic_volumes = 0d0 ! v_beta + do i_spec_organics = 1,aero_data_n_spec(aero_data) + if (use_species(i_spec_organics)) then + aero_state_organic_volumes = aero_state_organic_volumes & + + aero_particle_organic_species_volume( & + aero_state%apa%particle(1:aero_state_n_part(aero_state)), & + i_spec_organics) + end if + end do + + end function aero_state_organic_volumes + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the masses of all particles. @@ -1478,7 +1521,7 @@ end function aero_state_approx_crit_rel_humids !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical relative humidity for all particles (1). +! !> Returns the critical relative humidity for all particles (1). function aero_state_crit_rel_humids(aero_state, aero_data, env_state) !> Aerosol state. @@ -1500,6 +1543,31 @@ function aero_state_crit_rel_humids(aero_state, aero_data, env_state) end function aero_state_crit_rel_humids +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the critical relative humidity for all particles (1). + function aero_state_crit_rel_humids_varying_sigma(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_varying_sigma(aero_state_n_part(aero_state)) + + integer :: i_part + + do i_part = 1,aero_state_n_part(aero_state) + aero_state_crit_rel_humids_varying_sigma(i_part) = & + aero_particle_crit_rel_humid_varying_sigma( & + aero_state%apa%particle(i_part), aero_data, env_state) + end do + + end function aero_state_crit_rel_humids_varying_sigma + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Does the same thing as aero_state_to_bin() but based on dry radius. diff --git a/src/env_state.F90 b/src/env_state.F90 index 682ed32ae..57ccb581c 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -174,8 +174,8 @@ real(kind=dp) function env_state_air_molar_den(env_state) end function env_state_air_molar_den !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Condensation \f$A\f$ parameter. + ! constant sigma + ! !> Condensation \f$A\f$ parameter. real(kind=dp) function env_state_A(env_state) !> Environment state. @@ -187,7 +187,25 @@ real(kind=dp) function env_state_A(env_state) end function env_state_A !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! varying sigma + !> Condensation \f$A\f$ parameter. + real(kind=dp) function env_state_A_varying_sigma(sigma, env_state) + !> Environment state. + type(env_state_t), intent(in) :: env_state + real(kind=dp), intent(in) :: sigma + !only pass surface tension here + !> Aerosol particle. + !type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + !type(aero_data_t), intent(in) :: aero_data + + env_state_A_varying_sigma = 4d0 * sigma * const%water_molec_weight / & + (const%univ_gas_const * env_state%temp * const%water_density) + + end function env_state_A_varying_sigma + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Convert (ppb) to (molecules m^{-3}). real(kind=dp) function env_state_ppb_to_conc(env_state, ppb) From 7ecd973c70d0c1d4b3fa434c1ec8db48213eb5db Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 30 Oct 2022 20:23:26 -0500 Subject: [PATCH 04/33] add function for organic volume --- src/aero_particle.F90 | 34 +++++++++++++------ src/aero_state.F90 | 77 +++++++++++++++++++++--------------------- src/pmc_constants.mod | Bin 0 -> 1523 bytes 3 files changed, 63 insertions(+), 48 deletions(-) create mode 100644 src/pmc_constants.mod diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 6aadd1ee5..243e5af03 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -287,18 +287,32 @@ elemental real(kind=dp) function aero_particle_species_volume( & end function aero_particle_species_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Organic volume of a single species in the particle (m^3). xtxu, Oct/19/2022 - elemental real(kind=dp) function aero_particle_organic_species_volume( & + !> Organic volume of a single species 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 - !> Species number to find volume of. - integer, intent(in) :: i_spec_organics - - aero_particle_organic_species_volume = aero_particle%vol(i_spec_organics) + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Species names to include in the volume. + character(len=6), allocatable :: org_spec(:) + ! character(len=*), intent(in) :: include(:) + integer :: i_org_spec + integer :: n_org_spec + + 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_species_volume + end function aero_particle_organic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Total dry volume of the particle (m^3). @@ -306,7 +320,7 @@ elemental real(kind=dp) function aero_particle_dry_volume(aero_particle, & aero_data) !> Particle. - type(aero_particle_t), intent(in) :: aero_particle + type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. type(aero_data_t), intent(in) :: aero_data @@ -828,8 +842,8 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) - delta_min)**3) !> coverage parameter - coverage = min(aero_particle_organic_species_volume(aero_particle, & - i_spec_organics) / v_delta, 1d0) + coverage = min(aero_particle_organic_volume(aero_particle, & + aero_data) / v_delta, 1d0) !> fraction of water for inorganic core frac_core_water = ((aero_particle_diameter(aero_particle, aero_data) & diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 3bb4c4e62..92de3c309 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1121,47 +1121,48 @@ function aero_state_volumes(aero_state, aero_data, include, exclude) end function aero_state_volumes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> calculate organic species volume fraction, xtxu, Oct/17/2022 - - function aero_state_organic_volumes(aero_state, aero_data) - - !> Aerosol state. - type(aero_state_t), intent(in) :: aero_state - !> Aerosol data. - type(aero_data_t), optional, intent(in) :: aero_data - !> Species names to include in the mass. - character(len=*) :: include(:) - - !> Return volumes array (m^3). - real(kind=dp) :: aero_state_organic_volumes(aero_state_n_part(aero_state)) + !> calculate organic species volume fraction + !> might be useful later. + +! function aero_state_organic_volumes(aero_state, aero_data) + +! !> Aerosol state. +! type(aero_state_t), intent(in) :: aero_state +! !> Aerosol data. +! type(aero_data_t), optional, intent(in) :: aero_data +! !> Species names to include in the mass. +! character(len=*) :: include(:) + +! !> Return volumes array (m^3). +! real(kind=dp) :: aero_state_organic_volumes(aero_state_n_part(aero_state)) - logical, allocatable :: use_species(:) - integer :: i_name_organics, i_spec_organics +! logical, allocatable :: use_species(:) +! integer :: i_name_organics, i_spec_organics - include = ["MSA ", "ARO1", "ARO2", "ALK1", "OLE1", & - "API1", "API2", "LIM1", "LIM2", "OC "] +! include = ["MSA ", "ARO1", "ARO2", "ALK1", "OLE1", & +! "API1", "API2", "LIM1", "LIM2", "OC "] - allocate(use_species(aero_data_n_spec(aero_data))) - - use_species = .false. - do i_name_organics = 1, size(include) - i_spec_organics = aero_data_spec_by_name(aero_data, include(i_name_organics)) - call assert_msg(111852070, i_spec_organics > 0, & - "unknown species: " // trim(include(i_name_organics))) - use_species(i_spec_organics) = .true. - end do - - aero_state_organic_volumes = 0d0 ! v_beta - do i_spec_organics = 1,aero_data_n_spec(aero_data) - if (use_species(i_spec_organics)) then - aero_state_organic_volumes = aero_state_organic_volumes & - + aero_particle_organic_species_volume( & - aero_state%apa%particle(1:aero_state_n_part(aero_state)), & - i_spec_organics) - end if - end do - - end function aero_state_organic_volumes +! allocate(use_species(aero_data_n_spec(aero_data))) + +! use_species = .false. +! do i_name_organics = 1, size(include) +! i_spec_organics = aero_data_spec_by_name(aero_data, include(i_name_organics)) +! call assert_msg(111852070, i_spec_organics > 0, & +! "unknown species: " // trim(include(i_name_organics))) +! use_species(i_spec_organics) = .true. +! end do + +! aero_state_organic_volumes = 0d0 ! v_beta +! do i_spec_organics = 1,aero_data_n_spec(aero_data) +! if (use_species(i_spec_organics)) then +! aero_state_organic_volumes = aero_state_organic_volumes & +! + aero_particle_organic_species_volume( & +! aero_state%apa%particle(1:aero_state_n_part(aero_state)), & +! i_spec_organics) +! end if +! end do + +! end function aero_state_organic_volumes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/pmc_constants.mod b/src/pmc_constants.mod new file mode 100644 index 0000000000000000000000000000000000000000..d1af6b07f2dc42179ee4de88f99c145a051cf23f GIT binary patch literal 1523 zcmV-X}}NE8)kQd91lm;loBIk)e<=ic_|aqP?-!^X>1_Eh4$Mf{+P;1J#&A)m} zPu(pDSj3&qfI8MQ9bnh~YLCC!ca}Y)_N+RcIAar)C;TENm;j`X49|=O+;ONe#9DS0 ze)1E5TO+~8-#re`O@o`Jsm2ybADM<0RrBdsTDdS3~N=4^(yQ7~#Uo$%b z^qNXEQ4ajMCPuo({RMaWSV*^__qy?S7mbdtg#8I_t_R>NdZ_5Cs7c!1%{46eL9kka z!C4gDsH4|kdyDHx_#1NS#$$hZJAdu1H!MwWl7ZBd>kMShKqQ<(v|;QA^R{fczDWk6 zC`!HRRbwD>2I==r29o5@YcABQ z7|4+9V=(`k1-FY{W?^#xDBS+b)dC3o&EnU$+oVwXjw(m}fMs%9aZovjkR@rneZNWg zGNk{;73dZIbMyA^_k3B6*gVph1O#gK@0Sq(>U+BqAZqhq9|QMlgi7JX+BLQ6`M0Q(${(DIj$no>ccqn_S3cq+A0kI(PzGhc4D4&0=djjT{3C-F? z3A$15%}?LJ(NH$}i8@#*(ev}qYhxw=*O<|IDSa?0s!s|F%l(yyu>=jSm&QMsupg>5 z;zgJ6w&(y#yxF|arh4ojs9nM=@10%gPwBGpWSrgUd;p;W2p+M*`FQyp5c@zy91pJH zfkqO%a-Avhew5Td%1oaFX6imTV|QlR6wx6a!5GLN6W5{G7=ANerXKyHh*6w9!aD@K zcZDJ;qwzC!9BYuJ06Z}uOHU`%%;@6D#LPmBf&)v0|<UaY>q4t2<=1U_>)Wouj0pkJa{;}oiLkkkw|^lRI{&tzn`9UhXH;sXL+r&oS(;$9@aTEzIAy z%u6BFg{&F;FO;6o#?Ym8*ynn!%mU5ouYh7gaK#WmXN-JpN#?!T3V;=aw+2*IL`F6=zOqsYbP9&Qp5FEI(O>doS_ti;PIsBG`1s6XIH*{e64 Zsa*R=u7d{&3vTFl(Z6_0gNijD004sw1|t9f literal 0 HcmV?d00001 From 8ac8937b424baf9f231d25a8aaf4d6c23e407982 Mon Sep 17 00:00:00 2001 From: Xu Date: Thu, 3 Nov 2022 16:59:02 -0500 Subject: [PATCH 05/33] fix frac_core_water error --- src/aero_particle.F90 | 42 +++++++++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 243e5af03..a1927e23b 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -297,7 +297,6 @@ real(kind=dp) function aero_particle_organic_volume( & type(aero_data_t), intent(in) :: aero_data !> Species names to include in the volume. character(len=6), allocatable :: org_spec(:) - ! character(len=*), intent(in) :: include(:) integer :: i_org_spec integer :: n_org_spec @@ -315,6 +314,33 @@ real(kind=dp) function aero_particle_organic_volume( & end function aero_particle_organic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> inrganic volume of a single species in the particle (m^3). + real(kind=dp) function aero_particle_inorganic_volume( & + aero_particle, aero_data) + + !> Particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Species names to include in the volume. + character(len=6), allocatable :: inorg_spec(:) + integer :: i_inorg_spec + integer :: n_inorg_spec + + inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + "Na ", "Ca ", "OIN ", "BC "] + + aero_particle_inorganic_volume = 0d0 + + do n_inorg_spec = 1, size(inorg_spec) + i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + aero_particle_inorganic_volume = aero_particle_inorganic_volume & + + aero_particle%vol(i_inorg_spec) + end do + + end function aero_particle_inorganic_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Total dry volume of the particle (m^3). elemental real(kind=dp) function aero_particle_dry_volume(aero_particle, & aero_data) @@ -835,10 +861,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> sigma_shell real(kind=dp) :: sigma_shell = 0.03d0 - !> minimum shell volume, v_delta v_delta = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & - / 3d0) * ( aero_particle_radius(aero_particle, aero_data) & + / 3d0) * (aero_particle_radius(aero_particle, aero_data) & - delta_min)**3) !> coverage parameter @@ -846,10 +871,13 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) aero_data) / v_delta, 1d0) !> fraction of water for inorganic core - frac_core_water = ((aero_particle_diameter(aero_particle, aero_data) & - - 2d0 * delta_min)**3 - aero_particle_dry_volume(aero_particle, & - aero_data)**3) / (aero_particle_diameter(aero_particle, aero_data) & - - 2d0 * delta_min)**3 + !> v_core = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 + v_core = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & + / 3d0) * (aero_particle_radius(aero_particle, aero_data) & + - delta_min)**3) + + frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & + aero_data) / v_core !> surface tension of inorganic core and water sigma_core = f_core_water * const%water_surf_eng + (1d0 - & From 728ebbe9cfa22d8f25f84223b8815e4dbc7268af Mon Sep 17 00:00:00 2001 From: Xu Date: Thu, 3 Nov 2022 17:03:33 -0500 Subject: [PATCH 06/33] fix v_delta error --- src/aero_particle.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index a1927e23b..74f6c101f 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -852,7 +852,7 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) type(aero_data_t), intent(in) :: aero_data real(kind=dp) :: v_delta, coverage, frac_core_water - real(kind=dp) :: sigma_core + real(kind=dp) :: v_core, sigma_core !> Minimum shell thickness real(kind=dp) :: delta_min = 1.6d-10 @@ -862,9 +862,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) real(kind=dp) :: sigma_shell = 0.03d0 !> minimum shell volume, v_delta - v_delta = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & - / 3d0) * (aero_particle_radius(aero_particle, aero_data) & - - delta_min)**3) + !> v_delta = (4*pi/3)*(D/2 - delta_min)^3 + v_delta = ((4d0 * const%pi) / 3d0) * (aero_particle_radius(aero_particle, & + aero_data) - delta_min)**3 !> coverage parameter coverage = min(aero_particle_organic_volume(aero_particle, & @@ -873,8 +873,8 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> fraction of water for inorganic core !> v_core = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 v_core = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & - / 3d0) * (aero_particle_radius(aero_particle, aero_data) & - - delta_min)**3) + / 3d0) * (aero_particle_radius(aero_particle, aero_data) & + - delta_min)**3) frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & aero_data) / v_core From 8a1a4ffa162f8ebe794532e38ae32b28f2e713a5 Mon Sep 17 00:00:00 2001 From: Xu Date: Thu, 3 Nov 2022 20:43:57 -0500 Subject: [PATCH 07/33] add critical diameter for varying sigma --- src/aero_particle.F90 | 57 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 3 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 74f6c101f..ef2785f2a 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -824,7 +824,7 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, !pass this sigma to env_state sigma = aero_particle_varying_sigma(aero_particle, aero_data) - A_varying_sigma = env_state_A_varying_sigma( sigma, env_state) + A_varying_sigma = env_state_A_varying_sigma(sigma, env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) crit_diam = aero_particle_crit_diameter(aero_particle, aero_data, & @@ -856,9 +856,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> Minimum shell thickness real(kind=dp) :: delta_min = 1.6d-10 - !> sigma_core + !> sigma_core real(kind=dp) :: sigma_inorganic = 0.058d0 - !> sigma_shell + !> sigma_shell real(kind=dp) :: sigma_shell = 0.03d0 !> minimum shell volume, v_delta @@ -976,6 +976,57 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & end function aero_particle_crit_diameter + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the critical diameter (m) for varying simga. + real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& + 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 + + integer, parameter :: CRIT_DIAM_MAX_ITER = 100 + + real(kind=dp) :: kappa, dry_diam, A_varying_sigma, c4, c3, & + c0, d, f, df, dd + integer :: i_newton + + A_varying_sigma = env_state_A_varying_sigma(sigma, env_state) + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + if (kappa < 1d-30) then + ! bail out early for hydrophobic particles + aero_particle_crit_diameter_varying_sigma = dry_diam + return + end if + + c4 = - 3d0 * dry_diam**3 * kappa / A_varying_sigma + c3 = - dry_diam**3 * (2d0 - kappa) + c0 = dry_diam**6 * (1d0 - kappa) + + ! Newton's method for f(d) = d^6 + c4 d^4 + c3 d^3 + c0 + d = max(sqrt(-4d0 / 3d0 * c4), (-c3)**(1d0/3d0)) + do i_newton = 1,CRIT_DIAM_MAX_ITER + f = d**6 + c4 * d**4 + c3 * d**3 + c0 + df = 6 * d**5 + 4 * c4 * d**3 + 3 * c3 * d**2 + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-14) then + exit + end if + end do + call warn_assert_msg(408545686, i_newton < CRIT_DIAM_MAX_ITER, & + "critical diameter Newton loop failed to converge") + call warn_assert_msg(353290871, d >= dry_diam, & + "critical diameter Newton loop converged to invalid solution") + aero_particle_crit_diameter_varying_sigma = d + + end function aero_particle_crit_diameter_varying_sigma + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Coagulate two particles together to make a new one. The new From 008902d2c3474056fe402fffcbd69699e9776267 Mon Sep 17 00:00:00 2001 From: Xu Date: Thu, 3 Nov 2022 20:45:23 -0500 Subject: [PATCH 08/33] fix env_state_A_varying sigma error --- src/env_state.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/env_state.F90 b/src/env_state.F90 index 57ccb581c..276527613 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -199,7 +199,7 @@ real(kind=dp) function env_state_A_varying_sigma(sigma, env_state) !type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. !type(aero_data_t), intent(in) :: aero_data - + env_state_A_varying_sigma = 4d0 * sigma * const%water_molec_weight / & (const%univ_gas_const * env_state%temp * const%water_density) From 718bb07a7124b5d367a0b0e71a9592468f278678 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 4 Nov 2022 00:22:40 -0500 Subject: [PATCH 09/33] add env_state_A_varying_sigma --- src/env_state.F90 | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/env_state.F90 b/src/env_state.F90 index 276527613..e31c411d0 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -189,18 +189,13 @@ end function env_state_A !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! varying sigma !> Condensation \f$A\f$ parameter. - real(kind=dp) function env_state_A_varying_sigma(sigma, env_state) + real(kind=dp) function env_state_A_varying_sigma(varying_sigma, env_state) !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp), intent(in) :: sigma - !only pass surface tension here - !> Aerosol particle. - !type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. - !type(aero_data_t), intent(in) :: aero_data - - env_state_A_varying_sigma = 4d0 * sigma * const%water_molec_weight / & + real(kind=dp), intent(in) :: varying_sigma + + env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & (const%univ_gas_const * env_state%temp * const%water_density) end function env_state_A_varying_sigma From b54ca506ad218d1aacfd5735f1cb8afda1e4ec50 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 4 Nov 2022 01:31:19 -0500 Subject: [PATCH 10/33] rename sigma to vary_sigma --- src/aero_particle.F90 | 10 +++++----- src/env_state.F90 | 11 +++-------- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index ef2785f2a..90218eba9 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -820,11 +820,11 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp) :: kappa, crit_diam, dry_diam, A_varying_sigma, sigma + real(kind=dp) :: kappa, crit_diam, dry_diam, A_varying_sigma, varying_sigma !pass this sigma to env_state - sigma = aero_particle_varying_sigma(aero_particle, aero_data) - A_varying_sigma = env_state_A_varying_sigma(sigma, env_state) + varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + A_varying_sigma = env_state_A_varying_sigma(varying_sigma, env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) crit_diam = aero_particle_crit_diameter(aero_particle, aero_data, & @@ -992,10 +992,10 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& integer, parameter :: CRIT_DIAM_MAX_ITER = 100 real(kind=dp) :: kappa, dry_diam, A_varying_sigma, c4, c3, & - c0, d, f, df, dd + c0, d, f, df, dd, varying_sigma integer :: i_newton - A_varying_sigma = env_state_A_varying_sigma(sigma, env_state) + A_varying_sigma = env_state_A_varying_sigma(varying_sigma, env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) kappa = aero_particle_solute_kappa(aero_particle, aero_data) if (kappa < 1d-30) then diff --git a/src/env_state.F90 b/src/env_state.F90 index 57ccb581c..5fbee9c55 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -189,18 +189,13 @@ end function env_state_A !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! varying sigma !> Condensation \f$A\f$ parameter. - real(kind=dp) function env_state_A_varying_sigma(sigma, env_state) + real(kind=dp) function env_state_A_varying_sigma(varying_sigma, env_state) !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp), intent(in) :: sigma - !only pass surface tension here - !> Aerosol particle. - !type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. - !type(aero_data_t), intent(in) :: aero_data + real(kind=dp), intent(in) :: varying_sigma - env_state_A_varying_sigma = 4d0 * sigma * const%water_molec_weight / & + env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & (const%univ_gas_const * env_state%temp * const%water_density) end function env_state_A_varying_sigma From 74d12246739d7b8b926c92bbec9c19e30609fe2c Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 5 Nov 2022 16:15:29 -0500 Subject: [PATCH 11/33] correct v_delta and v_core --- src/aero_particle.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 90218eba9..5766dd4a1 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -862,8 +862,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) real(kind=dp) :: sigma_shell = 0.03d0 !> minimum shell volume, v_delta - !> v_delta = (4*pi/3)*(D/2 - delta_min)^3 - v_delta = ((4d0 * const%pi) / 3d0) * (aero_particle_radius(aero_particle, & + !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 + v_delta = aero_particle_volume(aero_particle) - ((4d0 * const%pi) & + / 3d0) * (aero_particle_radius(aero_particle, & aero_data) - delta_min)**3 !> coverage parameter @@ -871,10 +872,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) aero_data) / v_delta, 1d0) !> fraction of water for inorganic core - !> v_core = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 - v_core = aero_particle_volume(aero_particle) - (((4d0 * const%pi) & - / 3d0) * (aero_particle_radius(aero_particle, aero_data) & - - delta_min)**3) + !> v_core = (4*pi/3)*(D/2 - delta_min)^3 + v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_radius( & + aero_particle, aero_data) - delta_min)**3 frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & aero_data) / v_core From 4d9868e304c128c6e8943912cc535e9067a9b4c2 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 5 Nov 2022 18:05:15 -0500 Subject: [PATCH 12/33] add sc_varying_sigma related variables --- .../1_urban_plume/urban_plume_process.F90 | 45 +++++++++++++++---- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/scenarios/1_urban_plume/urban_plume_process.F90 b/scenarios/1_urban_plume/urban_plume_process.F90 index 818844673..a05dd1896 100644 --- a/scenarios/1_urban_plume/urban_plume_process.F90 +++ b/scenarios/1_urban_plume/urban_plume_process.F90 @@ -24,6 +24,7 @@ program process character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename, & single_particle_output type(bin_grid_t) :: diam_grid, bc_grid, oc_grid, sc_grid, & + sc_varying_sigma_grid, & so4_grid, no3_grid, nh4_grid, soa_grid, & cl_grid, msa_grid, aro1_grid, aro2_grid, & alk1_grid, ole1_grid, api1_grid, api2_grid, & @@ -51,6 +52,9 @@ program process dry_masses(:), masses(:), bc_masses(:), bc_fracs(:), & diam_bc_dist(:,:), oc_masses(:), oc_fracs(:), diam_oc_dist(:,:), & crit_rhs(:), scs(:), num_dist(:), diam_sc_dist(:,:), & + sc_dist(:), sc_varying_sigma_dist(:), & + crit_rhs_varying_sigma(:), scs_varying_sigma(:), & + diam_sc_varying_sigma_dist(:,:), & so4_masses(:), so4_fracs(:), diam_so4_dist(:,:), & no3_masses(:), no3_fracs(:), diam_no3_dist(:,:), & nh4_masses(:), nh4_fracs(:), diam_nh4_dist(:,:), & @@ -83,9 +87,11 @@ program process stats_tot_api2_mass_conc, stats_tot_lim1_mass_conc, & stats_tot_lim2_mass_conc, stats_tot_co3_mass_conc, & stats_tot_na_mass_conc, stats_tot_ca_mass_conc, & - stats_tot_oin_mass_conc, stats_tot_h2o_mass_conc + stats_tot_oin_mass_conc, stats_tot_h2o_mass_conc, & + stats_sc_dist, stats_sc_varying_sigma_dist type(stats_2d_t) :: stats_diam_sc_dist, stats_diam_bc_dist, & + stats_diam_sc_varying_sigma_dist, & stats_diam_oc_dist, stats_diam_so4_dist, & stats_diam_no3_dist, stats_diam_nh4_dist, & stats_diam_soa_dist, stats_diam_cl_dist, & @@ -126,6 +132,7 @@ program process call bin_grid_make(oin_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) call bin_grid_make(h2o_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) call bin_grid_make(sc_grid, BIN_GRID_TYPE_LOG, 50, 1d-4, 1d0) + call bin_grid_make(sc_varying_sigma_grid, BIN_GRID_TYPE_LOG, 50, 1d-4, 1d0) allocate(times(n_index)) @@ -136,6 +143,7 @@ program process scs = [ real(kind=dp) :: ] ! silence compiler warnings + scs_varying_sigma = [ real(kind=dp) :: ] bc_fracs = [ real(kind=dp) :: ] oc_fracs = [ real(kind=dp) :: ] so4_fracs = [ real(kind=dp) :: ] @@ -406,10 +414,26 @@ program process h2o_grid, h2o_fracs, num_concs) call stats_2d_add(stats_diam_h2o_dist, diam_h2o_dist) -!for critical RH supersaturation +!for critical RH supersaturation crit_rhs = aero_state_crit_rel_humids(aero_state, aero_data, & env_state) scs = crit_rhs - 1d0 + diam_sc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + sc_grid, scs, num_concs) + sc_dist = bin_grid_histogram_1d(sc_grid, scs, num_concs) + call stats_1d_add(stats_sc_dist, sc_dist) + call stats_2d_add(stats_diam_sc_dist, diam_sc_dist) + +!for critical RH supersaturation with varying sigma + crit_rhs_varying_sigma = aero_state_crit_rel_humids_varying_sigma( & + aero_state, aero_data, env_state) + scs_varying_sigma = crit_rhs_varying_sigma - 1d0 + diam_sc_varying_sigma_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & + sc_varying_sigma_grid, scs_varying_sigma, num_concs) + sc_varying_sigma_dist = bin_grid_histogram_1d(sc_varying_sigma_grid, & + scs_varying_sigma, num_concs) + call stats_1d_add(stats_sc_varying_sigma_dist, sc_varying_sigma_dist) + call stats_2d_add(stats_diam_sc_varying_sigma_dist, diam_sc_varying_sigma_dist) call for_single_particle(single_particle_output, prefix_new, ".csv", index) open(15, file= trim(single_particle_output)) @@ -440,13 +464,11 @@ program process write(15, *) 'crit_rhs', crit_rhs write(15, *) 'scs', scs write(15, *) 'num_con', num_concs + write(15, *) 'crit_rhs_varying_sigma', crit_rhs_varying_sigma + write(15, *) 'scs_varying_sigma', scs_varying_sigma close(15) - diam_sc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & - sc_grid, scs, num_concs) - call stats_2d_add(stats_diam_sc_dist, diam_sc_dist) - - call aero_state_mixing_state_metrics(aero_state, aero_data, & + call aero_state_mixing_state_metrics(aero_state, aero_data, & d_alpha, d_gamma, chi, groups=mixing_state_groups) call stats_1d_add_entry(stats_d_alpha, d_alpha, i_index) @@ -482,6 +504,8 @@ program process call bin_grid_output_netcdf(oin_grid, ncid, "oin_frac", unit="1") call bin_grid_output_netcdf(h2o_grid, ncid, "h2o_frac", unit="1") call bin_grid_output_netcdf(sc_grid, ncid, "sc", unit="1") + call bin_grid_output_netcdf(sc_varying_sigma_grid, ncid, & + "sc_varying_sigma", unit="1") call stats_1d_output_netcdf(stats_num_dist, ncid, "num_dist", & dim_name="diam", unit="m^{-3}") @@ -572,9 +596,14 @@ program process call stats_2d_clear(stats_diam_h2o_dist) call stats_2d_output_netcdf(stats_diam_sc_dist, ncid, "diam_sc_dist", & - dim_name_1="diam", dim_name_2="sc", unit="m^{-3}") + 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_varying_sigma_dist, ncid, & + "diam_sc_varying_sigma_dist", & + dim_name_1="diam", dim_name_2="sc_varying_sigma", unit="m^{-3}") + call stats_2d_clear(stats_diam_sc_varying_sigma_dist) + call pmc_nc_close(ncid) end do From cb2125d9c531084303328552f20c2c3f6de3e441 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 5 Nov 2022 22:14:19 -0500 Subject: [PATCH 13/33] add sigma value to aero_data.dat --- scenarios/1_urban_plume/aero_data.dat | 42 +++++++++++++-------------- 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/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 +-- 14 files changed, 110 insertions(+), 110 deletions(-) diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index be5c1cb19..8cd99073b 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) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 0.53 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 0.53 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +SO4 1800 0 96d-3 0.65 0.05d0 +NO3 1800 0 62d-3 0.65 0.05d0 +Cl 2200 0 35.5d-3 0.53 0.05d0 +NH4 1800 0 18d-3 0.65 0.05d0 +MSA 1800 0 95d-3 0.53 0.03d0 +ARO1 1400 0 150d-3 0.1 0.03d0 +ARO2 1400 0 150d-3 0.1 0.03d0 +ALK1 1400 0 140d-3 0.1 0.03d0 +OLE1 1400 0 140d-3 0.1 0.03d0 +API1 1400 0 184d-3 0.1 0.03d0 +API2 1400 0 184d-3 0.1 0.03d0 +LIM1 1400 0 200d-3 0.1 0.03d0 +LIM2 1400 0 200d-3 0.1 0.03d0 +CO3 2600 0 60d-3 0.53 0.05d0 +Na 2200 0 23d-3 0.53 0.05d0 +Ca 2600 0 40d-3 0.53 0.05d0 +OIN 2600 0 1d-3 0.1 0.05d0 +OC 1000 0 1d-3 0.001 0.03d0 +BC 1800 0 1d-3 0 0.05d0 +H2O 1000 0 18d-3 0 0.073d0 diff --git a/test/additive/aero_data.dat b/test/additive/aero_data.dat index 41cdc37fc..ec050a51b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 diff --git a/test/average/aero_data.dat b/test/average/aero_data.dat index 5984b0e75..2b44b297f 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 -H2O 1000 0 18d-3 0 -SO4 1800 0 96d-3 0.53 -NO3 1800 0 62d-3 0.53 -Cl 2200 0 35.5d-3 0.53 -NH4 1800 0 18d-3 0.53 -CO3 2600 0 60d-3 0.53 -MSA 1800 0 95d-3 0.53 -Na 2200 0 23d-3 0.53 -Ca 2600 0 40d-3 0.53 -OC 1000 0 1d-3 0.1 -BC 1800 0 1d-3 0 -OIN 2600 0 1d-3 0.1 -ARO1 1000 0 150d-3 0.1 -ARO2 1000 0 150d-3 0.1 -ALK1 1000 0 140d-3 0.1 -OLE1 1000 0 140d-3 0.1 -API1 1000 0 184d-3 0.1 -API2 1000 0 184d-3 0.1 -LIM1 1000 0 200d-3 0.1 -LIM2 1000 0 200d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 +SO4 1800 0 96d-3 0.65 0.05d0 +NO3 1800 0 62d-3 0.65 0.05d0 +Cl 2200 0 35.5d-3 0.53 0.05d0 +NH4 1800 0 18d-3 0.65 0.05d0 +MSA 1800 0 95d-3 0.53 0.03d0 +ARO1 1400 0 150d-3 0.1 0.03d0 +ARO2 1400 0 150d-3 0.1 0.03d0 +ALK1 1400 0 140d-3 0.1 0.03d0 +OLE1 1400 0 140d-3 0.1 0.03d0 +API1 1400 0 184d-3 0.1 0.03d0 +API2 1400 0 184d-3 0.1 0.03d0 +LIM1 1400 0 200d-3 0.1 0.03d0 +LIM2 1400 0 200d-3 0.1 0.03d0 +CO3 2600 0 60d-3 0.53 0.05d0 +Na 2200 0 23d-3 0.53 0.05d0 +Ca 2600 0 40d-3 0.53 0.05d0 +OIN 2600 0 1d-3 0.1 0.05d0 +OC 1000 0 1d-3 0.001 0.03d0 +BC 1800 0 1d-3 0 0.05d0 diff --git a/test/bidisperse/aero_data.dat b/test/bidisperse/aero_data.dat index 41cdc37fc..ec050a51b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 diff --git a/test/brownian/aero_data.dat b/test/brownian/aero_data.dat index 41cdc37fc..ec050a51b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 diff --git a/test/condense/aero_data.dat b/test/condense/aero_data.dat index af22101b2..36c4c9b1b 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) -SO4 1800 1 96d-3 0 -NO3 1800 1 62d-3 0 -Cl 2200 1 35.5d-3 0 -NH4 1800 1 18d-3 0 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 1 60d-3 0 -Na 2200 1 23d-3 0 -Ca 2600 1 40d-3 0 -OIN 2600 0 1d-3 0.1 -OC 1400 0 1d-3 0.1 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +SO4 1800 0 96d-3 0.65 0.05d0 +NO3 1800 0 62d-3 0.65 0.05d0 +Cl 2200 0 35.5d-3 0.53 0.05d0 +NH4 1800 0 18d-3 0.65 0.05d0 +MSA 1800 0 95d-3 0.53 0.03d0 +ARO1 1400 0 150d-3 0.1 0.03d0 +ARO2 1400 0 150d-3 0.1 0.03d0 +ALK1 1400 0 140d-3 0.1 0.03d0 +OLE1 1400 0 140d-3 0.1 0.03d0 +API1 1400 0 184d-3 0.1 0.03d0 +API2 1400 0 184d-3 0.1 0.03d0 +LIM1 1400 0 200d-3 0.1 0.03d0 +LIM2 1400 0 200d-3 0.1 0.03d0 +CO3 2600 0 60d-3 0.53 0.05d0 +Na 2200 0 23d-3 0.53 0.05d0 +Ca 2600 0 40d-3 0.53 0.05d0 +OIN 2600 0 1d-3 0.1 0.05d0 +OC 1000 0 1d-3 0.001 0.03d0 +BC 1800 0 1d-3 0 0.05d0 +H2O 1000 0 18d-3 0 0.073d0 diff --git a/test/emission/aero_data.dat b/test/emission/aero_data.dat index 83d7a7091..3671fc4ac 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 -SB 100 0 18d-3 0 -SE 1500 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +SI 10000 0 18d-3 0 0.073d0 +SB 100 0 18d-3 0 0.073d0 +SE 1500 0 18d-3 0 0.073d0 diff --git a/test/fractal/aero_data.dat b/test/fractal/aero_data.dat index 21e82154f..2c227b908 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) -BC 4200 0 1d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +BC 4200 0 1d-3 0 0.073d0 diff --git a/test/loss/aero_data.dat b/test/loss/aero_data.dat index 41cdc37fc..b7e15530b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 \ No newline at end of file diff --git a/test/mixing_state/aero_data.dat b/test/mixing_state/aero_data.dat index 740f9933b..fabcf9815 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) -A 1000 0 100d-3 0.1 -B 1500 0 100d-3 0.1 -C 2000 0 100d-3 0.1 -D 2500 0 100d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +A 1000 0 100d-3 0.1 0.073d0 +B 1500 0 100d-3 0.1 0.073d0 +C 2000 0 100d-3 0.1 0.073d0 +D 2500 0 100d-3 0.1 0.073d0 diff --git a/test/mosaic/aero_data.dat b/test/mosaic/aero_data.dat index 5984b0e75..388fd1690 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 -H2O 1000 0 18d-3 0 -SO4 1800 0 96d-3 0.53 -NO3 1800 0 62d-3 0.53 -Cl 2200 0 35.5d-3 0.53 -NH4 1800 0 18d-3 0.53 -CO3 2600 0 60d-3 0.53 -MSA 1800 0 95d-3 0.53 -Na 2200 0 23d-3 0.53 -Ca 2600 0 40d-3 0.53 -OC 1000 0 1d-3 0.1 -BC 1800 0 1d-3 0 -OIN 2600 0 1d-3 0.1 -ARO1 1000 0 150d-3 0.1 -ARO2 1000 0 150d-3 0.1 -ALK1 1000 0 140d-3 0.1 -OLE1 1000 0 140d-3 0.1 -API1 1000 0 184d-3 0.1 -API2 1000 0 184d-3 0.1 -LIM1 1000 0 200d-3 0.1 -LIM2 1000 0 200d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 +SO4 1800 0 96d-3 0.65 0.05d0 +NO3 1800 0 62d-3 0.65 0.05d0 +Cl 2200 0 35.5d-3 0.53 0.05d0 +NH4 1800 0 18d-3 0.65 0.05d0 +CO3 2600 0 60d-3 0.53 0.05d0 +MSA 1800 0 95d-3 0.53 0.03d0 +Na 2200 0 23d-3 0.53 0.05d0 +Ca 2600 0 40d-3 0.53 0.05d0 +OC 1000 0 1d-3 0.001 0.03d0 +BC 1800 0 1d-3 0 0.05d0 +OIN 2600 0 1d-3 0.1 0.05d0 +ARO1 1400 0 150d-3 0.1 0.03d0 +ARO2 1400 0 150d-3 0.1 0.03d0 +ALK1 1400 0 140d-3 0.1 0.03d0 +OLE1 1400 0 140d-3 0.1 0.03d0 +API1 1400 0 184d-3 0.1 0.03d0 +API2 1400 0 184d-3 0.1 0.03d0 +LIM1 1400 0 200d-3 0.1 0.03d0 +LIM2 1400 0 200d-3 0.1 0.03d0 diff --git a/test/nucleate/aero_data.dat b/test/nucleate/aero_data.dat index 1f0f968fa..dee42a8e6 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) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +SO4 1800 0 96d-3 0.65 0.05d0 +NO3 1800 0 62d-3 0.65 0.05d0 diff --git a/test/parallel/aero_data.dat b/test/parallel/aero_data.dat index 41cdc37fc..b7e15530b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 \ No newline at end of file diff --git a/test/sedi/aero_data.dat b/test/sedi/aero_data.dat index 41cdc37fc..b7e15530b 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) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +H2O 1000 0 18d-3 0 0.073d0 \ No newline at end of file From c04e719ab2097ba8f220ee1e38e33db4a9538a15 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 5 Nov 2022 22:16:27 -0500 Subject: [PATCH 14/33] add aero_data%sigma --- src/aero_data.F90 | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/aero_data.F90 b/src/aero_data.F90 index db4131b12..eba341f6b 100644 --- a/src/aero_data.F90 +++ b/src/aero_data.F90 @@ -62,6 +62,8 @@ module pmc_aero_data real(kind=dp), allocatable :: molec_weight(:) !> Length \c aero_data_n_spec(aero_data), kappas (1). real(kind=dp), allocatable :: kappa(:) + !> 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(:) !> Fractal particle parameters. @@ -430,9 +432,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) == 4) .or. (n_species == 0))) then + if (.not. ((size(species_data, 2) == 5) .or. (n_species == 0))) then call die_msg(428926381, 'each line in ' // trim(file%name) & - // ' should contain exactly 5 values') + // ' should contain exactly 6 values') end if ! allocate and copy over the data @@ -442,12 +444,14 @@ subroutine spec_file_read_aero_data(file, aero_data) call ensure_integer_array_size(aero_data%num_ions, n_species) call ensure_real_array_size(aero_data%molec_weight, n_species) call ensure_real_array_size(aero_data%kappa, 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) aero_data%num_ions(i) = nint(species_data(i,2)) aero_data%molec_weight(i) = species_data(i,3) aero_data%kappa(i) = species_data(i,4) + aero_data%sigma(i) = species_data(i,5) call assert_msg(232362742, & (aero_data%num_ions(i) == 0) .or. (aero_data%kappa(i) == 0d0), & "ions and kappa both non-zero for species " & @@ -519,6 +523,7 @@ integer function pmc_mpi_pack_size_aero_data(val) + pmc_mpi_pack_size_integer_array(val%num_ions) & + pmc_mpi_pack_size_real_array(val%molec_weight) & + pmc_mpi_pack_size_real_array(val%kappa) & + + pmc_mpi_pack_size_real_array(val%sigma) & + pmc_mpi_pack_size_string_array(val%source_name) & + pmc_mpi_pack_size_fractal(val%fractal) @@ -547,6 +552,7 @@ subroutine pmc_mpi_pack_aero_data(buffer, position, val) call pmc_mpi_pack_integer_array(buffer, position, val%num_ions) call pmc_mpi_pack_real_array(buffer, position, val%molec_weight) call pmc_mpi_pack_real_array(buffer, position, val%kappa) + 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_fractal(buffer, position, val%fractal) call assert(183834856, & @@ -578,6 +584,7 @@ subroutine pmc_mpi_unpack_aero_data(buffer, position, val) call pmc_mpi_unpack_integer_array(buffer, position, val%num_ions) call pmc_mpi_unpack_real_array(buffer, position, val%molec_weight) call pmc_mpi_unpack_real_array(buffer, position, val%kappa) + 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_fractal(buffer, position, val%fractal) call assert(188522823, & @@ -760,6 +767,9 @@ subroutine aero_data_output_netcdf(aero_data, ncid) call pmc_nc_write_real_1d(ncid, aero_data%kappa, & "aero_kappa", (/ dimid_aero_species /), unit="1", & long_name="hygroscopicity parameters (kappas) of aerosol species") + 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 fractal_output_netcdf(aero_data%fractal, ncid) end subroutine aero_data_output_netcdf @@ -801,6 +811,7 @@ subroutine aero_data_input_netcdf(aero_data, ncid) call pmc_nc_read_integer_1d(ncid, aero_data%num_ions, "aero_num_ions") call pmc_nc_read_real_1d(ncid, aero_data%molec_weight, "aero_molec_weight") call pmc_nc_read_real_1d(ncid, aero_data%kappa, "aero_kappa") + 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)) @@ -899,6 +910,7 @@ subroutine aero_data_initialize(aero_data, camp_core) allocate(aero_data%num_ions(num_spec)) allocate(aero_data%molec_weight(num_spec)) allocate(aero_data%kappa(num_spec)) + allocate(aero_data%sigma(num_spec)) allocate(aero_data%camp_particle_spec_id(num_spec)) ! Assume no aerosol water @@ -936,6 +948,12 @@ subroutine aero_data_initialize(aero_data, camp_core) call die_msg(944207343, "Missing kappa 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(944207343, "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 From 74c86923f1d34aaad2c6132723175353dcc71d85 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 5 Nov 2022 22:20:11 -0500 Subject: [PATCH 15/33] apply sigma value of each species into effective surface tension calculation --- src/aero_particle.F90 | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 5766dd4a1..5e6009233 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -861,21 +861,20 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> sigma_shell real(kind=dp) :: sigma_shell = 0.03d0 - !> minimum shell volume, v_delta - !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 - v_delta = aero_particle_volume(aero_particle) - ((4d0 * const%pi) & - / 3d0) * (aero_particle_radius(aero_particle, & - aero_data) - delta_min)**3 - !> coverage parameter - coverage = min(aero_particle_organic_volume(aero_particle, & - aero_data) / v_delta, 1d0) - - !> fraction of water for inorganic core !> v_core = (4*pi/3)*(D/2 - delta_min)^3 v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_radius( & aero_particle, aero_data) - delta_min)**3 + + !> minimum shell volume, v_delta + !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 = D^3/6 - v_core + v_delta = aero_particle_volume(aero_particle) - v_core + !> coverage parameter + coverage = min(aero_particle_organic_volume(aero_particle, & + aero_data) / v_delta, 1d0) + + !> fraction of water for inorganic core frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & aero_data) / v_core From 4440c3a675337d085fef6defd33a8572486358f5 Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 6 Nov 2022 14:48:03 -0600 Subject: [PATCH 16/33] apply sigma value of each species to effective surface tension method --- src/aero_particle.F90 | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 5766dd4a1..d868ddd30 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -311,6 +311,12 @@ real(kind=dp) function aero_particle_organic_volume( & + aero_particle%vol(i_org_spec) end do + sigma_shell = 0d0 + do n_org_spec = 1, size(org_spec) + i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) + sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & + * aero_data%sigma(i_org_spec) / aero_particle_organic_volume + end do end function aero_particle_organic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -338,6 +344,13 @@ real(kind=dp) function aero_particle_inorganic_volume( & + aero_particle%vol(i_inorg_spec) end do + sigma_core_without_water = 0d0 + do n_inorg_spec = 1, size(inorg_spec) + i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & + i_inorg_spec) * aero_data%sigma(i_inorg_spec) / & + aero_particle_inorganic_volume + end do end function aero_particle_inorganic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -856,32 +869,25 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> Minimum shell thickness real(kind=dp) :: delta_min = 1.6d-10 - !> sigma_core - real(kind=dp) :: sigma_inorganic = 0.058d0 - !> sigma_shell - real(kind=dp) :: sigma_shell = 0.03d0 - + + !> v_core = (4*pi/3)*(D/2 - delta_min)^3 + v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_radius( & + aero_particle, aero_data) - delta_min)**3 + !> minimum shell volume, v_delta !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 - v_delta = aero_particle_volume(aero_particle) - ((4d0 * const%pi) & - / 3d0) * (aero_particle_radius(aero_particle, & - aero_data) - delta_min)**3 + v_delta = aero_particle_volume(aero_particle) - v_core !> coverage parameter coverage = min(aero_particle_organic_volume(aero_particle, & aero_data) / v_delta, 1d0) !> fraction of water for inorganic core - !> v_core = (4*pi/3)*(D/2 - delta_min)^3 - v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_radius( & - aero_particle, aero_data) - delta_min)**3 - frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & aero_data) / v_core !> surface tension of inorganic core and water - sigma_core = f_core_water * const%water_surf_eng + (1d0 - & - frac_core_water) * sigma_inorganic + sigma_core = f_core_water * const%water_surf_eng + sigma_core_without_water !> calculate effective surface tension value aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & @@ -1208,4 +1214,4 @@ end subroutine aero_particle_check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module pmc_aero_particle +end module pmc_aero_particle \ No newline at end of file From ae27428e2ab2ce780c875eb102bf4dfe8ff843e7 Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 6 Nov 2022 15:02:13 -0600 Subject: [PATCH 17/33] correct crit_diam --- src/aero_particle.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 105e86829..56cee9ea9 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -840,7 +840,7 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, A_varying_sigma = env_state_A_varying_sigma(varying_sigma, env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - crit_diam = aero_particle_crit_diameter(aero_particle, aero_data, & + crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & env_state) kappa = aero_particle_solute_kappa(aero_particle, aero_data) if (kappa < 1d-30) then @@ -1214,4 +1214,4 @@ end subroutine aero_particle_check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module pmc_aero_particle +end module pmc_aero_particle \ No newline at end of file From 7edba1dc062e4d4add8f1b6f9cd58e83586d30c4 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 11 Nov 2022 13:31:16 -0600 Subject: [PATCH 18/33] use equation following files --- src/aero_particle.F90 | 115 ++++++++++++++++++++++++++---------------- 1 file changed, 71 insertions(+), 44 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 56cee9ea9..860f82318 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -311,47 +311,46 @@ real(kind=dp) function aero_particle_organic_volume( & + aero_particle%vol(i_org_spec) end do - sigma_shell = 0d0 - do n_org_spec = 1, size(org_spec) - i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) - sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & - * aero_data%sigma(i_org_spec) / aero_particle_organic_volume - end do + ! sigma_shell = 0d0 + ! do n_org_spec = 1, size(org_spec) + ! i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) + ! sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & + ! * aero_data%sigma(i_org_spec) / v_delta + ! end do end function aero_particle_organic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> inrganic volume of a single species in the particle (m^3). - real(kind=dp) function aero_particle_inorganic_volume( & - aero_particle, aero_data) - - !> Particle. - type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. - type(aero_data_t), intent(in) :: aero_data - !> Species names to include in the volume. - character(len=6), allocatable :: inorg_spec(:) - integer :: i_inorg_spec - integer :: n_inorg_spec + ! real(kind=dp) function aero_particle_inorganic_volume( & + ! aero_particle, aero_data) + + ! !> Particle. + ! type(aero_particle_t), intent(in) :: aero_particle + ! !> Aerosol data. + ! type(aero_data_t), intent(in) :: aero_data + ! !> Species names to include in the volume. + ! character(len=6), allocatable :: inorg_spec(:) + ! integer :: i_inorg_spec + ! integer :: n_inorg_spec - inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - "Na ", "Ca ", "OIN ", "BC "] - - aero_particle_inorganic_volume = 0d0 - - do n_inorg_spec = 1, size(inorg_spec) - i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - aero_particle_inorganic_volume = aero_particle_inorganic_volume & - + aero_particle%vol(i_inorg_spec) - end do - - sigma_core_without_water = 0d0 - do n_inorg_spec = 1, size(inorg_spec) - i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & - i_inorg_spec) * aero_data%sigma(i_inorg_spec) / & - aero_particle_inorganic_volume - end do - end function aero_particle_inorganic_volume + ! inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + ! "Na ", "Ca ", "OIN ", "BC "] + + ! aero_particle_inorganic_volume = 0d0 + + ! do n_inorg_spec = 1, size(inorg_spec) + ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + ! aero_particle_inorganic_volume = aero_particle_inorganic_volume & + ! + aero_particle%vol(i_inorg_spec) + ! end do + + ! sigma_core_without_water = 0d0 + ! do n_inorg_spec = 1, size(inorg_spec) + ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + ! sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & + ! i_inorg_spec) * aero_data%sigma(i_inorg_spec) / v_core + ! end do + ! end function aero_particle_inorganic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Total dry volume of the particle (m^3). @@ -864,30 +863,58 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - real(kind=dp) :: v_delta, coverage, frac_core_water - real(kind=dp) :: v_core, sigma_core + ! real(kind=dp) :: v_delta, coverage, frac_core_water + ! real(kind=dp) :: v_core, sigma_core + real(kind=dp) :: v_core, v_delta, coverage + + character(len=6), allocatable :: shell_spec(:) + character(len=6), allocatable :: core_spec(:) + integer :: i_shell_spec + integer :: n_shell_spec + integer :: i_core_spec + integer :: n_core_spec !> Minimum shell thickness real(kind=dp) :: delta_min = 1.6d-10 !> v_core = (4*pi/3)*(D/2 - delta_min)^3 - v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_radius( & - aero_particle, aero_data) - delta_min)**3 + v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_diameter(aero_particle, & + aero_data) / 2d0 - delta_min)**3 !> minimum shell volume, v_delta - !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 + !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 = D^3/6 - v_core v_delta = aero_particle_volume(aero_particle) - v_core !> coverage parameter coverage = min(aero_particle_organic_volume(aero_particle, & aero_data) / v_delta, 1d0) + shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + + sigma_shell = 0d0 + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & + * aero_data%sigma(i_shell_spec) / v_delta + end do + + core_spec = ["H2O ", "SO4 ", "NO3 ", "Cl ", "NH4 ", & + "CO3 ", "Na ", "Ca ", "OIN ", "BC "] + + sigma_core = 0d0 + do n_core_spec = 1, size(core_spec) + i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) + sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & + aero_data%sigma(i_core_spec) / v_core + end do + !> fraction of water for inorganic core - frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & - aero_data) / v_core + !frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & + ! aero_data) / v_core !> surface tension of inorganic core and water - sigma_core = f_core_water * const%water_surf_eng + sigma_core_without_water + !sigma_core = f_core_water * const%water_surf_eng + sigma_core_without_water !> calculate effective surface tension value aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & From 7568c36929845a4dc7f01fe278809062a73f96fb Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 11 Nov 2022 13:50:01 -0600 Subject: [PATCH 19/33] add varying sigma defination --- src/env_state.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/env_state.F90 b/src/env_state.F90 index 5fbee9c55..240cb8ba0 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -193,8 +193,9 @@ real(kind=dp) function env_state_A_varying_sigma(varying_sigma, env_state) !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp), intent(in) :: varying_sigma + real(kind=dp):: varying_sigma + varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & (const%univ_gas_const * env_state%temp * const%water_density) From efde0198178d8a0e93ad2cff17e996b329c90f2e Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 12 Nov 2022 20:56:37 -0600 Subject: [PATCH 20/33] upload before new test --- src/aero_particle.F90 | 176 ++++++++++++++++++++++++------------------ src/env_state.F90 | 16 ++-- 2 files changed, 107 insertions(+), 85 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 860f82318..347172d23 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -13,6 +13,7 @@ module pmc_aero_particle use pmc_spec_file use pmc_env_state use pmc_mpi + use pmc_constants #ifdef PMC_USE_MPI use mpi #endif @@ -287,72 +288,6 @@ elemental real(kind=dp) function aero_particle_species_volume( & end function aero_particle_species_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Organic volume of a single species 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 - !> Species names to include in the volume. - character(len=6), allocatable :: org_spec(:) - integer :: i_org_spec - integer :: n_org_spec - - 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 - - ! sigma_shell = 0d0 - ! do n_org_spec = 1, size(org_spec) - ! i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) - ! sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & - ! * aero_data%sigma(i_org_spec) / v_delta - ! end do - end function aero_particle_organic_volume -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> inrganic volume of a single species in the particle (m^3). - ! real(kind=dp) function aero_particle_inorganic_volume( & - ! aero_particle, aero_data) - - ! !> Particle. - ! type(aero_particle_t), intent(in) :: aero_particle - ! !> Aerosol data. - ! type(aero_data_t), intent(in) :: aero_data - ! !> Species names to include in the volume. - ! character(len=6), allocatable :: inorg_spec(:) - ! integer :: i_inorg_spec - ! integer :: n_inorg_spec - - ! inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - ! "Na ", "Ca ", "OIN ", "BC "] - - ! aero_particle_inorganic_volume = 0d0 - - ! do n_inorg_spec = 1, size(inorg_spec) - ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - ! aero_particle_inorganic_volume = aero_particle_inorganic_volume & - ! + aero_particle%vol(i_inorg_spec) - ! end do - - ! sigma_core_without_water = 0d0 - ! do n_inorg_spec = 1, size(inorg_spec) - ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - ! sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & - ! i_inorg_spec) * aero_data%sigma(i_inorg_spec) / v_core - ! end do - ! end function aero_particle_inorganic_volume -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Total dry volume of the particle (m^3). elemental real(kind=dp) function aero_particle_dry_volume(aero_particle, & aero_data) @@ -821,6 +756,71 @@ end function aero_particle_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Organic volume of a single species 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 + !> Species names to include in the volume. + character(len=6), allocatable :: org_spec(:) + integer :: i_org_spec + integer :: n_org_spec + + 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 + + ! sigma_shell = 0d0 + ! do n_org_spec = 1, size(org_spec) + ! i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) + ! sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & + ! * aero_data%sigma(i_org_spec) / v_delta + ! end do + end function aero_particle_organic_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> inrganic volume of a single species in the particle (m^3). + ! real(kind=dp) function aero_particle_inorganic_volume( & + ! aero_particle, aero_data) + + ! !> Particle. + ! type(aero_particle_t), intent(in) :: aero_particle + ! !> Aerosol data. + ! type(aero_data_t), intent(in) :: aero_data + ! !> Species names to include in the volume. + ! character(len=6), allocatable :: inorg_spec(:) + ! integer :: i_inorg_spec + ! integer :: n_inorg_spec + + ! inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + ! "Na ", "Ca ", "OIN ", "BC "] + + ! aero_particle_inorganic_volume = 0d0 + + ! do n_inorg_spec = 1, size(inorg_spec) + ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + ! aero_particle_inorganic_volume = aero_particle_inorganic_volume & + ! + aero_particle%vol(i_inorg_spec) + ! end do + + ! sigma_core_without_water = 0d0 + ! do n_inorg_spec = 1, size(inorg_spec) + ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + ! sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & + ! i_inorg_spec) * aero_data%sigma(i_inorg_spec) / v_core + ! end do + ! end function aero_particle_inorganic_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the critical relative humidity (1). real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, & aero_data, env_state) @@ -832,11 +832,15 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp) :: kappa, crit_diam, dry_diam, A_varying_sigma, varying_sigma + real(kind=dp) :: kappa, crit_diam, dry_diam, A_varying_sigma!, varying_sigma !pass this sigma to env_state - varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) - A_varying_sigma = env_state_A_varying_sigma(varying_sigma, env_state) + ! varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + ! A_varying_sigma = env_state_A_varying_sigma(aero_particle, aero_data, env_state) + + A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & + * const%water_molec_weight / (const%univ_gas_const * env_state%temp & + * const%water_density) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & @@ -863,9 +867,7 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - ! real(kind=dp) :: v_delta, coverage, frac_core_water - ! real(kind=dp) :: v_core, sigma_core - real(kind=dp) :: v_core, v_delta, coverage + real(kind=dp) :: v_core, v_delta, coverage, shell character(len=6), allocatable :: shell_spec(:) character(len=6), allocatable :: core_spec(:) @@ -889,15 +891,32 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) coverage = min(aero_particle_organic_volume(aero_particle, & aero_data) / v_delta, 1d0) - shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & - "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + write(*,*) aero_particle_volume(aero_particle) / & + (v_core + aero_particle_organic_volume(aero_particle, aero_data)) + ! shell = aero_particle_diameter(aero_particle, aero_data) / 2d0 - & + ! ((aero_particle_diameter(aero_particle, aero_data))**3 - & + ! 3*aero_particle_organic_volume(aero_particle, aero_data)/(4*const%pi))**(1d0/3d0) + ! write(*,*) shell + + ! write(*, *) (aero_particle_diameter(aero_particle, aero_data) / 2d0)/ delta_min + ! write(*,*) aero_particle_organic_volume(aero_particle, aero_data) / v_delta + + shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + sigma_shell = 0d0 do n_shell_spec = 1, size(shell_spec) i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & * aero_data%sigma(i_shell_spec) / v_delta + ! write(*,"(E24.3)") "sigma_shell = ", sigma_shell + ! write(*,*) "sigma_species = ", aero_data%sigma(i_shell_spec) + ! write(*,*) "species", "species_sequence", "species_volume", "species_sigma" + ! write(*,*) shell_spec(n_shell_spec), i_shell_spec, aero_particle%vol(i_shell_spec), & + ! aero_data%sigma(i_shell_spec) end do + ! write(*,*) "sigma_shell", sigma_shell core_spec = ["H2O ", "SO4 ", "NO3 ", "Cl ", "NH4 ", & "CO3 ", "Na ", "Ca ", "OIN ", "BC "] @@ -908,7 +927,7 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & aero_data%sigma(i_core_spec) / v_core end do - + ! write(*,*) "sigma_core", sigma_core !> fraction of water for inorganic core !frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & ! aero_data) / v_core @@ -1025,10 +1044,13 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& integer, parameter :: CRIT_DIAM_MAX_ITER = 100 real(kind=dp) :: kappa, dry_diam, A_varying_sigma, c4, c3, & - c0, d, f, df, dd, varying_sigma + c0, d, f, df, dd!, varying_sigma integer :: i_newton - A_varying_sigma = env_state_A_varying_sigma(varying_sigma, env_state) + ! A_varying_sigma = env_state_A_varying_sigma(aero_particle, aero_data, env_state) + A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & + * 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) if (kappa < 1d-30) then diff --git a/src/env_state.F90 b/src/env_state.F90 index 240cb8ba0..16b5e4ab1 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -189,17 +189,17 @@ end function env_state_A !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! varying sigma !> Condensation \f$A\f$ parameter. - real(kind=dp) function env_state_A_varying_sigma(varying_sigma, env_state) + ! real(kind=dp) function env_state_A_varying_sigma(aero_particle, aero_data, env_state) - !> Environment state. - type(env_state_t), intent(in) :: env_state - real(kind=dp):: varying_sigma + ! !> Environment state. + ! type(env_state_t), intent(in) :: env_state + ! real(kind=dp):: varying_sigma - varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) - env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & - (const%univ_gas_const * env_state%temp * const%water_density) + ! varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + ! env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & + ! (const%univ_gas_const * env_state%temp * const%water_density) - end function env_state_A_varying_sigma + ! end function env_state_A_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Convert (ppb) to (molecules m^{-3}). From 7f1e083e2d7fd8110f4d5ac7127c9d4f6c2df972 Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 13 Nov 2022 15:02:17 -0600 Subject: [PATCH 21/33] modified effective surface tension method --- src/aero_particle.F90 | 257 +++++++++++++++++++++++------------------- 1 file changed, 144 insertions(+), 113 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 347172d23..34e493a7f 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -756,46 +756,38 @@ end function aero_particle_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Organic volume of a single species in the particle (m^3). - real(kind=dp) function aero_particle_organic_volume( & - aero_particle, aero_data) + ! !> Organic volume of a single species 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 - !> Species names to include in the volume. - character(len=6), allocatable :: org_spec(:) - integer :: i_org_spec - integer :: n_org_spec + ! !> Particle. + ! type(aero_particle_t), intent(in) :: aero_particle + ! !> Aerosol data. + ! type(aero_data_t), intent(in) :: aero_data + ! !> Species names to include in the volume. + ! character(len=6), allocatable :: org_spec(:) + ! integer :: i_org_spec + ! integer :: n_org_spec - org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & - "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + ! org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + ! "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] - aero_particle_organic_volume = 0d0 + ! 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 - - ! sigma_shell = 0d0 - ! do n_org_spec = 1, size(org_spec) - ! i_org_spec = aero_data_spec_by_name(aero_data, org_spec(n_org_spec)) - ! sigmal_shell = sigma_shell + aero_particle%vol(i_org_spec) & - ! * aero_data%sigma(i_org_spec) / v_delta - ! end do - end function aero_particle_organic_volume + ! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> inrganic volume of a single species in the particle (m^3). - ! real(kind=dp) function aero_particle_inorganic_volume( & - ! aero_particle, aero_data) + ! real(kind=dp) function aero_particle_inorganic_volume(aero_particle, aero_data) ! !> Particle. ! type(aero_particle_t), intent(in) :: aero_particle - ! !> Aerosol data. + ! !> Aerosol data.ls ! type(aero_data_t), intent(in) :: aero_data ! !> Species names to include in the volume. ! character(len=6), allocatable :: inorg_spec(:) @@ -813,52 +805,42 @@ end function aero_particle_organic_volume ! + aero_particle%vol(i_inorg_spec) ! end do - ! sigma_core_without_water = 0d0 - ! do n_inorg_spec = 1, size(inorg_spec) - ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - ! sigma_core_without_water = sigma_core_without_water + aero_particle%vol( & - ! i_inorg_spec) * aero_data%sigma(i_inorg_spec) / v_core - ! end do ! end function aero_particle_inorganic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical relative humidity (1). - real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(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) :: kappa, crit_diam, dry_diam, A_varying_sigma!, varying_sigma - !pass this sigma to env_state + ! !> core volume of the particle (m^3). + ! real(kind=dp) function aero_particle_core_volume(aero_particle, aero_data) - ! varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) - ! A_varying_sigma = env_state_A_varying_sigma(aero_particle, aero_data, env_state) - - A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & - * const%water_molec_weight / (const%univ_gas_const * env_state%temp & - * const%water_density) + ! !> Particle. + ! type(aero_particle_t), intent(in) :: aero_particle + ! !> Aerosol data.ls + ! type(aero_data_t), intent(in) :: aero_data + ! !> Species names to include in the volume. + ! character(len=6), allocatable :: core_spec(:) + ! integer :: i_core_spec + ! integer :: n_core_spec + + ! core_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + ! "Na ", "Ca ", "OIN ", "BC ", "H2O "] - dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & - env_state) - kappa = aero_particle_solute_kappa(aero_particle, aero_data) - if (kappa < 1d-30) then - aero_particle_crit_rel_humid_varying_sigma = exp(A_varying_sigma & - / crit_diam) - else - aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & - / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & - exp(A_varying_sigma / crit_diam) - end if + ! aero_particle_core_volume = 0d0 - end function aero_particle_crit_rel_humid_varying_sigma + ! do n_core_spec = 1, size(core_spec) + ! i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) + ! aero_particle_core_volume = aero_particle_core_volume & + ! + aero_particle%vol(i_core_spec) + ! end do - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! sigma_core = 0d0 + ! do n_core_spec = 1, size(core_spec) + ! i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) + ! sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & + ! aero_data%sigma(i_core_spec) / v_core + ! end do + ! end function aero_particle_core_volume +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the varying sigma. real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) @@ -866,81 +848,131 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - - real(kind=dp) :: v_core, v_delta, coverage, shell - + !> inorganic species names to include in the volume. + character(len=6), allocatable :: inorg_spec(:) + integer :: i_inorg_spec + integer :: n_inorg_spec + !> shell species names to include in the volume. character(len=6), allocatable :: shell_spec(:) - character(len=6), allocatable :: core_spec(:) integer :: i_shell_spec integer :: n_shell_spec + !> core/organic species names to include in the volume. + character(len=6), allocatable :: core_spec(:) integer :: i_core_spec integer :: n_core_spec + + real(kind=dp) :: v_core, r_core, r1, r2, delta_min, v_delta, coverage + + !> calculate inorganic volume + inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + "Na ", "Ca ", "OIN ", "BC "] - !> Minimum shell thickness - real(kind=dp) :: delta_min = 1.6d-10 + aero_particle_inorganic_volume = 0d0 - !> v_core = (4*pi/3)*(D/2 - delta_min)^3 - v_core = ((4d0 * const%pi) / 3d0) * (aero_particle_diameter(aero_particle, & - aero_data) / 2d0 - delta_min)**3 + do n_inorg_spec = 1, size(inorg_spec) + i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) + aero_particle_inorganic_volume = aero_particle_inorganic_volume & + + aero_particle%vol(i_inorg_spec) + end do - !> minimum shell volume, v_delta - !> v_delta = D^3/6 - (4*pi/3)*(D/2 - delta_min)^3 = D^3/6 - v_core - v_delta = aero_particle_volume(aero_particle) - v_core + !> calculate organic volume + shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] - !> coverage parameter - coverage = min(aero_particle_organic_volume(aero_particle, & - aero_data) / v_delta, 1d0) + aero_particle_organic_volume = 0d0 - write(*,*) aero_particle_volume(aero_particle) / & - (v_core + aero_particle_organic_volume(aero_particle, aero_data)) + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + aero_particle_organic_volume = aero_particle_organic_volume & + + aero_particle%vol(i_shell_spec) + end do - ! shell = aero_particle_diameter(aero_particle, aero_data) / 2d0 - & - ! ((aero_particle_diameter(aero_particle, aero_data))**3 - & - ! 3*aero_particle_organic_volume(aero_particle, aero_data)/(4*const%pi))**(1d0/3d0) - ! write(*,*) shell + !> calculate core volume + core_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & + "Na ", "Ca ", "OIN ", "BC ", "H2O "] - ! write(*, *) (aero_particle_diameter(aero_particle, aero_data) / 2d0)/ delta_min - ! write(*,*) aero_particle_organic_volume(aero_particle, aero_data) / v_delta + aero_particle_core_volume = 0d0 - shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & - "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + do n_core_spec = 1, size(core_spec) + i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) + aero_particle_core_volume = aero_particle_core_volume & + + aero_particle%vol(i_core_spec) + end do + + v_core = aero_particle_core_volume + r_core = (3d0 * v_core / 4d0 * const%pi)**(1d0/3d0) + + !> r1 is the radius for dry inorganic core + !> r2 is the radius for dry inorganic + orgnaic core + r1 = (3d0 * aero_particle_inorganic_volume / 4d0 * const%pi)**(1d0/3d0) + r2 = (3d0 * (aero_particle_organic_volume + aero_particle_inorganic_volume) / & + 4d0 * const%pi)**(1d0/3d0) + !> Minimum shell thickness + delta_min = r2 - r1 + + !> minimum shell volume, v_delta + v_delta = (4d0 * const%pi / 3d0) * (r_core + delta_min)**3 - v_core + + !> coverage parameter + coverage = min(aero_particle_organic_volume / v_delta, 1d0) + !> calculate sigma for shell sigma_shell = 0d0 do n_shell_spec = 1, size(shell_spec) i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & * aero_data%sigma(i_shell_spec) / v_delta - ! write(*,"(E24.3)") "sigma_shell = ", sigma_shell - ! write(*,*) "sigma_species = ", aero_data%sigma(i_shell_spec) - ! write(*,*) "species", "species_sequence", "species_volume", "species_sigma" - ! write(*,*) shell_spec(n_shell_spec), i_shell_spec, aero_particle%vol(i_shell_spec), & - ! aero_data%sigma(i_shell_spec) end do - ! write(*,*) "sigma_shell", sigma_shell - - core_spec = ["H2O ", "SO4 ", "NO3 ", "Cl ", "NH4 ", & - "CO3 ", "Na ", "Ca ", "OIN ", "BC "] + !> calculate sigma for core sigma_core = 0d0 do n_core_spec = 1, size(core_spec) i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & aero_data%sigma(i_core_spec) / v_core end do - ! write(*,*) "sigma_core", sigma_core - !> fraction of water for inorganic core - !frac_core_water = 1 - aero_particle_inorganic_volume(aero_particle, & - ! aero_data) / v_core - - !> surface tension of inorganic core and water - !sigma_core = f_core_water * const%water_surf_eng + sigma_core_without_water !> calculate effective surface tension value aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & coverage * sigma_shell + write(*, *) aero_particle_varying_sigma + end function aero_particle_varying_sigma + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical relative humidity (1). + real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, & + aero_data, env_state) - end function aero_particle_varying_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 + + real(kind=dp) :: kappa, crit_diam, dry_diam, A_varying_sigma + !pass this sigma to env_state + + A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & + * const%water_molec_weight / (const%univ_gas_const * env_state%temp & + * const%water_density) + + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & + env_state) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + if (kappa < 1d-30) then + aero_particle_crit_rel_humid_varying_sigma = exp(A_varying_sigma & + / crit_diam) + else + aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & + / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & + exp(A_varying_sigma / crit_diam) + end if + + end function aero_particle_crit_rel_humid_varying_sigma + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the critical diameter (m). !! @@ -1044,10 +1076,9 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& integer, parameter :: CRIT_DIAM_MAX_ITER = 100 real(kind=dp) :: kappa, dry_diam, A_varying_sigma, c4, c3, & - c0, d, f, df, dd!, varying_sigma + c0, d, f, df, dd integer :: i_newton - ! A_varying_sigma = env_state_A_varying_sigma(aero_particle, aero_data, env_state) A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & * const%water_molec_weight / (const%univ_gas_const * env_state%temp & * const%water_density) From 857f62886aedeb5a2791a01275f1a71864be83cf Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 13 Nov 2022 17:03:28 -0600 Subject: [PATCH 22/33] add IF for particles with large organic fraction --- src/aero_particle.F90 | 82 ++++++++++++++++++++++++++++++------------- 1 file changed, 57 insertions(+), 25 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 34e493a7f..7e65f572b 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -861,7 +861,9 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) integer :: i_core_spec integer :: n_core_spec - real(kind=dp) :: v_core, r_core, r1, r2, delta_min, v_delta, coverage + real(kind=dp) :: v_core, r_core, r1, r2, delta_r, v_delta, coverage + !> Minimum shell thickness + real(kind=dp) :: delta_min = 1.6d-10 !> calculate inorganic volume inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & @@ -901,30 +903,10 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) v_core = aero_particle_core_volume r_core = (3d0 * v_core / 4d0 * const%pi)**(1d0/3d0) - - !> r1 is the radius for dry inorganic core - !> r2 is the radius for dry inorganic + orgnaic core - r1 = (3d0 * aero_particle_inorganic_volume / 4d0 * const%pi)**(1d0/3d0) - r2 = (3d0 * (aero_particle_organic_volume + aero_particle_inorganic_volume) / & - 4d0 * const%pi)**(1d0/3d0) - !> Minimum shell thickness - delta_min = r2 - r1 - !> minimum shell volume, v_delta - v_delta = (4d0 * const%pi / 3d0) * (r_core + delta_min)**3 - v_core + v_delta = (4d0 * const%pi / 3d0) * (r_core + delta_min)**3 - v_core - !> coverage parameter - coverage = min(aero_particle_organic_volume / v_delta, 1d0) - - !> calculate sigma for shell - sigma_shell = 0d0 - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & - * aero_data%sigma(i_shell_spec) / v_delta - end do - - !> calculate sigma for core + !> calculate sigma for core, sigma_core may not need to use, but it is fixed. sigma_core = 0d0 do n_core_spec = 1, size(core_spec) i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) @@ -932,9 +914,59 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) aero_data%sigma(i_core_spec) / v_core end do - !> calculate effective surface tension value - aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & + + !> r1 is the radius for dry inorganic core + !> r2 is the radius for dry inorganic + orgnaic core + r1 = (3d0 * aero_particle_inorganic_volume / 4d0 * const%pi)**(1d0/3d0) + r2 = (3d0 * (aero_particle_organic_volume + aero_particle_inorganic_volume) / & + 4d0 * const%pi)**(1d0/3d0) + !> shell thickness + delta_r = r2 - r1 + + if (delta_r/delta_min <= 1d0) then + ! use delta_min = 0.16 nm + coverage = aero_particle_organic_volume / v_delta + + !> calculate sigma for shell + sigma_shell = 0d0 + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & + * aero_data%sigma(i_shell_spec) / v_delta + end do + + aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & coverage * sigma_shell + + else if (aero_particle_organic_volume < v_delta) then + ! use delta_min = 0.16 nm + coverage = aero_particle_organic_volume / v_delta + + !> calculate sigma for shell + sigma_shell = 0d0 + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & + * aero_data%sigma(i_shell_spec) / v_delta + end do + + aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & + coverage * sigma_shell + else + ! v_delta = v_organic + ! coverage = 1 + !> calculate sigma for shell + sigma_shell = 0d0 + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & + * aero_data%sigma(i_shell_spec) / aero_particle_organic_volume + end do + + aero_particle_varying_sigma = sigma_shell + end if + + write(*, *) aero_particle_varying_sigma end function aero_particle_varying_sigma From c47115da8be682c57e52430b0f52514064ad3ac8 Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 13 Nov 2022 17:04:16 -0600 Subject: [PATCH 23/33] print sc_dist --- scenarios/1_urban_plume/urban_plume_process.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/scenarios/1_urban_plume/urban_plume_process.F90 b/scenarios/1_urban_plume/urban_plume_process.F90 index a05dd1896..223b1330f 100644 --- a/scenarios/1_urban_plume/urban_plume_process.F90 +++ b/scenarios/1_urban_plume/urban_plume_process.F90 @@ -604,6 +604,14 @@ program process dim_name_1="diam", dim_name_2="sc_varying_sigma", unit="m^{-3}") call stats_2d_clear(stats_diam_sc_varying_sigma_dist) + call stats_1d_output_netcdf(stats_sc_dist, ncid, "sc_dist", & + dim_name="sc", unit="C") + call stats_1d_clear(stats_sc_dist) + + call stats_1d_output_netcdf(stats_sc_varying_sigma_dist, ncid, "sc_varying_sigma_dist", & + dim_name="sc", unit="C") + call stats_1d_clear(stats_sc_varying_sigma_dist) + call pmc_nc_close(ncid) end do From 415caff753fbefa17a46b85b3a3786479d503d20 Mon Sep 17 00:00:00 2001 From: Xu Date: Wed, 30 Nov 2022 17:06:47 -0600 Subject: [PATCH 24/33] update surface tension volume --- scenarios/1_urban_plume/aero_data.dat | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index 8cd99073b..fdc00fb4d 100644 --- a/scenarios/1_urban_plume/aero_data.dat +++ b/scenarios/1_urban_plume/aero_data.dat @@ -1,8 +1,8 @@ # dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma -SO4 1800 0 96d-3 0.65 0.05d0 -NO3 1800 0 62d-3 0.65 0.05d0 -Cl 2200 0 35.5d-3 0.53 0.05d0 -NH4 1800 0 18d-3 0.65 0.05d0 +SO4 1800 0 96d-3 0.65 0.073d0 +NO3 1800 0 62d-3 0.65 0.073d0 +Cl 2200 0 35.5d-3 0.53 0.073d0 +NH4 1800 0 18d-3 0.65 0.073d0 MSA 1800 0 95d-3 0.53 0.03d0 ARO1 1400 0 150d-3 0.1 0.03d0 ARO2 1400 0 150d-3 0.1 0.03d0 @@ -12,10 +12,10 @@ API1 1400 0 184d-3 0.1 API2 1400 0 184d-3 0.1 0.03d0 LIM1 1400 0 200d-3 0.1 0.03d0 LIM2 1400 0 200d-3 0.1 0.03d0 -CO3 2600 0 60d-3 0.53 0.05d0 -Na 2200 0 23d-3 0.53 0.05d0 -Ca 2600 0 40d-3 0.53 0.05d0 -OIN 2600 0 1d-3 0.1 0.05d0 +CO3 2600 0 60d-3 0.53 0.073d0 +Na 2200 0 23d-3 0.53 0.073d0 +Ca 2600 0 40d-3 0.53 0.073d0 +OIN 2600 0 1d-3 0.1 0.00d0 OC 1000 0 1d-3 0.001 0.03d0 -BC 1800 0 1d-3 0 0.05d0 +BC 1800 0 1d-3 0 0.00d0 H2O 1000 0 18d-3 0 0.073d0 From b7cc3272489b662f6fb2e1e74fd2ee6f63bff458 Mon Sep 17 00:00:00 2001 From: Xu Date: Wed, 30 Nov 2022 18:27:10 -0600 Subject: [PATCH 25/33] update sigma value for each species --- scenarios/1_urban_plume/aero_data.dat | 42 +++++++++++++-------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index 43422f091..f86a8afc4 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) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +SO4 1800 0 96d-3 0.65 0.073 +NO3 1800 0 62d-3 0.65 0.073 +Cl 2200 0 35.5d-3 1.28 0.073 +NH4 1800 0 18d-3 0.65 0.073 +MSA 1800 0 95d-3 0.53 0.03 +ARO1 1400 0 150d-3 0.1 0.03 +ARO2 1400 0 150d-3 0.1 0.03 +ALK1 1400 0 140d-3 0.1 0.03 +OLE1 1400 0 140d-3 0.1 0.03 +API1 1400 0 184d-3 0.1 0.03 +API2 1400 0 184d-3 0.1 0.03 +LIM1 1400 0 200d-3 0.1 0.03 +LIM2 1400 0 200d-3 0.1 0.03 +CO3 2600 0 60d-3 0.53 0.073 +Na 2200 0 23d-3 1.28 0.073 +Ca 2600 0 40d-3 0.53 0.073 +OIN 2600 0 1d-3 0.1 0 +OC 1000 0 1d-3 0.001 0.03 +BC 1800 0 1d-3 0 0 +H2O 1000 0 18d-3 0 0.073 From 3aa7cb7b344b5150d0e3a0739d31e717c08c568c Mon Sep 17 00:00:00 2001 From: Xu Date: Sun, 4 Dec 2022 00:44:29 -0600 Subject: [PATCH 26/33] seperate sigma from env_state_A --- src/aero_particle.F90 | 81 +++++++++++++++++++++---------------------- src/env_state.F90 | 17 +-------- 2 files changed, 40 insertions(+), 58 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 7e65f572b..0ecaac3ce 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -746,14 +746,45 @@ real(kind=dp) function aero_particle_crit_rel_humid(aero_particle, & env_state) kappa = aero_particle_solute_kappa(aero_particle, aero_data) if (kappa < 1d-30) then - aero_particle_crit_rel_humid = exp(A / crit_diam) + aero_particle_crit_rel_humid = exp(A * const%water_surf_eng / crit_diam) else aero_particle_crit_rel_humid = (crit_diam**3 - dry_diam**3) & - / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * exp(A / crit_diam) + / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * exp(A * & + const%water_surf_eng / crit_diam) end if end function aero_particle_crit_rel_humid +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Returns the critical relative humidity (1). + real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(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) :: kappa, crit_diam, dry_diam, A, varying_sigma + + A = env_state_A(env_state) + varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & + env_state) + kappa = aero_particle_solute_kappa(aero_particle, aero_data) + if (kappa < 1d-30) then + aero_particle_crit_rel_humid_varying_sigma = exp(A * varying_sigma & + / crit_diam) + else + aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & + / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & + exp(A * varying_sigma / crit_diam) + end if + + end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !> Organic volume of a single species in the particle (m^3). @@ -969,40 +1000,6 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) write(*, *) aero_particle_varying_sigma end function aero_particle_varying_sigma - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical relative humidity (1). - real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(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) :: kappa, crit_diam, dry_diam, A_varying_sigma - !pass this sigma to env_state - - A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & - * const%water_molec_weight / (const%univ_gas_const * env_state%temp & - * const%water_density) - - dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & - env_state) - kappa = aero_particle_solute_kappa(aero_particle, aero_data) - if (kappa < 1d-30) then - aero_particle_crit_rel_humid_varying_sigma = exp(A_varying_sigma & - / crit_diam) - else - aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & - / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & - exp(A_varying_sigma / crit_diam) - end if - - end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1069,7 +1066,7 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & return end if - c4 = - 3d0 * dry_diam**3 * kappa / A + c4 = - 3d0 * dry_diam**3 * kappa / (A * const%water_surf_eng) c3 = - dry_diam**3 * (2d0 - kappa) c0 = dry_diam**6 * (1d0 - kappa) @@ -1107,13 +1104,13 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& integer, parameter :: CRIT_DIAM_MAX_ITER = 100 - real(kind=dp) :: kappa, dry_diam, A_varying_sigma, c4, c3, & + real(kind=dp) :: kappa, dry_diam, A, c4, c3, & c0, d, f, df, dd integer :: i_newton - A_varying_sigma = 4d0 * aero_particle_varying_sigma(aero_particle, aero_data) & - * const%water_molec_weight / (const%univ_gas_const * env_state%temp & - * const%water_density) + A = env_state_A(env_state) + varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) kappa = aero_particle_solute_kappa(aero_particle, aero_data) if (kappa < 1d-30) then @@ -1122,7 +1119,7 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& return end if - c4 = - 3d0 * dry_diam**3 * kappa / A_varying_sigma + c4 = - 3d0 * dry_diam**3 * kappa / (A * varying_sigma) c3 = - dry_diam**3 * (2d0 - kappa) c0 = dry_diam**6 * (1d0 - kappa) diff --git a/src/env_state.F90 b/src/env_state.F90 index 16b5e4ab1..983fdd823 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -181,26 +181,11 @@ real(kind=dp) function env_state_A(env_state) !> Environment state. type(env_state_t), intent(in) :: env_state - env_state_A = 4d0 * const%water_surf_eng * const%water_molec_weight & + env_state_A = 4d0 * const%water_molec_weight & / (const%univ_gas_const * env_state%temp * const%water_density) end function env_state_A -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! varying sigma - !> Condensation \f$A\f$ parameter. - ! real(kind=dp) function env_state_A_varying_sigma(aero_particle, aero_data, env_state) - - ! !> Environment state. - ! type(env_state_t), intent(in) :: env_state - ! real(kind=dp):: varying_sigma - - ! varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) - ! env_state_A_varying_sigma = 4d0 * varying_sigma * const%water_molec_weight / & - ! (const%univ_gas_const * env_state%temp * const%water_density) - - ! end function env_state_A_varying_sigma - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Convert (ppb) to (molecules m^{-3}). real(kind=dp) function env_state_ppb_to_conc(env_state, ppb) From bc725c332e55874d1f0fd75bd915060c74d525a0 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 25 Feb 2023 16:01:48 -0600 Subject: [PATCH 27/33] update Newton's loop --- src/aero_particle.F90 | 401 +++++++++++++++++++----------------------- 1 file changed, 185 insertions(+), 216 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 0ecaac3ce..94c07b06f 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -718,8 +718,9 @@ real(kind=dp) function aero_particle_approx_crit_rel_humid(aero_particle, & real(kind=dp) :: kappa, diam, C, A + A = env_state_A(env_state) kappa = aero_particle_solute_kappa(aero_particle, aero_data) - C = sqrt(4d0 * env_state_A(env_state)**3 / 27d0) + C = sqrt(4d0 * (A * const%water_surf_eng)**3 / 27d0) diam = aero_particle_diameter(aero_particle, aero_data) aero_particle_approx_crit_rel_humid = C / sqrt(kappa * diam**3) + 1d0 @@ -756,7 +757,8 @@ real(kind=dp) function aero_particle_crit_rel_humid(aero_particle, & end function aero_particle_crit_rel_humid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical relative humidity (1). + + ! > Returns the critical relative humidity (1). real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, & aero_data, env_state) @@ -767,136 +769,84 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, !> Environment state. type(env_state_t), intent(in) :: env_state - real(kind=dp) :: kappa, crit_diam, dry_diam, A, varying_sigma - + real(kind=dp) :: kappa, dry_diam, A, sigma, crit_diam_varying_sigma + real(kind=dp) :: v_delta, v_beta !, v_inorg, r_inor, v_dry_shell + real(kind=dp) :: sigma_core, sigma_shell, c_beta + ! real(kind=dp) :: sigma_core = 0.073d0 + ! real(kind=dp) :: sigma_shell = 0.030d0 + real(kind=dp) :: delta_min = 1.6d-10 + + sigma_core = aero_particle_sigma_core(aero_particle, aero_data) + sigma_shell = aero_particle_sigma_shell(aero_particle, aero_data) + A = env_state_A(env_state) - varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - crit_diam = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, & - env_state) + crit_diam_varying_sigma = aero_particle_crit_diameter_varying_sigma(& + aero_particle, aero_data, env_state) + v_beta = aero_particle_organic_volume(aero_particle, aero_data) + v_delta = 4d0 * const%pi * ((crit_diam_varying_sigma / 2d0)**3 - & + (crit_diam_varying_sigma / 2d0 - delta_min)**3) / 3d0 + kappa = aero_particle_solute_kappa(aero_particle, aero_data) - if (kappa < 1d-30) then - aero_particle_crit_rel_humid_varying_sigma = exp(A * varying_sigma & - / crit_diam) - else - aero_particle_crit_rel_humid_varying_sigma = (crit_diam**3 - dry_diam**3) & - / (crit_diam**3 - dry_diam**3 * (1 - kappa)) * & - exp(A * varying_sigma / crit_diam) - end if + + if (kappa < 1d-30) then ! try to catch pure BC particles + aero_particle_crit_rel_humid_varying_sigma = exp(A * const%water_surf_eng / dry_diam) + else + if (v_beta > v_delta) then + sigma = sigma_shell + else if (v_beta < 1d-30) then + sigma = sigma_core + else + sigma = sigma_core + (v_beta/v_delta) * (sigma_shell - sigma_core) + end if + + aero_particle_crit_rel_humid_varying_sigma = (crit_diam_varying_sigma**3 & + - dry_diam**3) / (crit_diam_varying_sigma**3 - dry_diam**3 * (1 - kappa)) * & + exp(A * sigma / crit_diam_varying_sigma) + end if end function aero_particle_crit_rel_humid_varying_sigma -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! !> Organic volume of a single species 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 - ! !> Species names to include in the volume. - ! character(len=6), allocatable :: org_spec(:) - ! integer :: i_org_spec - ! integer :: n_org_spec - - ! 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> inrganic volume of a single species in the particle (m^3). - ! real(kind=dp) function aero_particle_inorganic_volume(aero_particle, aero_data) - - ! !> Particle. - ! type(aero_particle_t), intent(in) :: aero_particle - ! !> Aerosol data.ls - ! type(aero_data_t), intent(in) :: aero_data - ! !> Species names to include in the volume. - ! character(len=6), allocatable :: inorg_spec(:) - ! integer :: i_inorg_spec - ! integer :: n_inorg_spec - - ! inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - ! "Na ", "Ca ", "OIN ", "BC "] - - ! aero_particle_inorganic_volume = 0d0 - - ! do n_inorg_spec = 1, size(inorg_spec) - ! i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - ! aero_particle_inorganic_volume = aero_particle_inorganic_volume & - ! + aero_particle%vol(i_inorg_spec) - ! end do - - ! end function aero_particle_inorganic_volume -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Organic volume of a single species in the particle (m^3). + real(kind=dp) function aero_particle_organic_volume(aero_particle, aero_data) - ! !> core volume of the particle (m^3). - ! real(kind=dp) function aero_particle_core_volume(aero_particle, aero_data) - - ! !> Particle. - ! type(aero_particle_t), intent(in) :: aero_particle - ! !> Aerosol data.ls - ! type(aero_data_t), intent(in) :: aero_data - ! !> Species names to include in the volume. - ! character(len=6), allocatable :: core_spec(:) - ! integer :: i_core_spec - ! integer :: n_core_spec + !> Particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Species names to include in the volume. + character(len=6), allocatable :: org_spec(:) + integer :: i_org_spec + integer :: n_org_spec - ! core_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - ! "Na ", "Ca ", "OIN ", "BC ", "H2O "] + org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] - ! aero_particle_core_volume = 0d0 - - ! do n_core_spec = 1, size(core_spec) - ! i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) - ! aero_particle_core_volume = aero_particle_core_volume & - ! + aero_particle%vol(i_core_spec) - ! end do - - ! sigma_core = 0d0 - ! do n_core_spec = 1, size(core_spec) - ! i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) - ! sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & - ! aero_data%sigma(i_core_spec) / v_core - ! end do + 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_core_volume + end function aero_particle_organic_volume !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Returns the varying sigma. - real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) - !> Aerosol particle. + ! > inrganic volume of a single species in the particle (m^3). + real(kind=dp) function aero_particle_inorganic_volume(aero_particle, aero_data) + + !> Particle. type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. + !> Aerosol data.ls type(aero_data_t), intent(in) :: aero_data - !> inorganic species names to include in the volume. + !> Species names to include in the volume. character(len=6), allocatable :: inorg_spec(:) integer :: i_inorg_spec integer :: n_inorg_spec - !> shell species names to include in the volume. - character(len=6), allocatable :: shell_spec(:) - integer :: i_shell_spec - integer :: n_shell_spec - !> core/organic species names to include in the volume. - character(len=6), allocatable :: core_spec(:) - integer :: i_core_spec - integer :: n_core_spec - - real(kind=dp) :: v_core, r_core, r1, r2, delta_r, v_delta, coverage - !> Minimum shell thickness - real(kind=dp) :: delta_min = 1.6d-10 - - !> calculate inorganic volume + inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & "Na ", "Ca ", "OIN ", "BC "] @@ -908,17 +858,22 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) + aero_particle%vol(i_inorg_spec) end do - !> calculate organic volume - shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & - "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + end function aero_particle_inorganic_volume - aero_particle_organic_volume = 0d0 - - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - aero_particle_organic_volume = aero_particle_organic_volume & - + aero_particle%vol(i_shell_spec) - end do +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the sigma of core. + real(kind=dp) function aero_particle_sigma_core(aero_particle, aero_data) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + + !> core/organic species names to include in the volume. + character(len=6), allocatable :: core_spec(:) + integer :: i_core_spec + integer :: n_core_spec !> calculate core volume core_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & @@ -928,78 +883,63 @@ real(kind=dp) function aero_particle_varying_sigma(aero_particle, aero_data) do n_core_spec = 1, size(core_spec) i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) + ! write(*, *) i_core_spec, ',', aero_particle%vol(i_core_spec), ',' aero_particle_core_volume = aero_particle_core_volume & + aero_particle%vol(i_core_spec) end do - v_core = aero_particle_core_volume - r_core = (3d0 * v_core / 4d0 * const%pi)**(1d0/3d0) - !> minimum shell volume, v_delta - v_delta = (4d0 * const%pi / 3d0) * (r_core + delta_min)**3 - v_core - - !> calculate sigma for core, sigma_core may not need to use, but it is fixed. - sigma_core = 0d0 + aero_particle_sigma_core = 0d0 do n_core_spec = 1, size(core_spec) i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) - sigma_core = sigma_core + aero_particle%vol(i_core_spec) * & - aero_data%sigma(i_core_spec) / v_core + aero_particle_sigma_core = aero_particle_sigma_core + aero_particle%vol(i_core_spec) * & + aero_data%sigma(i_core_spec) / aero_particle_core_volume end do + end function aero_particle_sigma_core + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the sigma of shell + real(kind=dp) function aero_particle_sigma_shell(aero_particle, aero_data) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + + !> shell/org species names to include in the volume. + character(len=6), allocatable :: shell_spec(:) + integer :: i_shell_spec + integer :: n_shell_spec + + !> calculate organic volume + shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + + aero_particle_shell_volume = 0d0 - !> r1 is the radius for dry inorganic core - !> r2 is the radius for dry inorganic + orgnaic core - r1 = (3d0 * aero_particle_inorganic_volume / 4d0 * const%pi)**(1d0/3d0) - r2 = (3d0 * (aero_particle_organic_volume + aero_particle_inorganic_volume) / & - 4d0 * const%pi)**(1d0/3d0) - !> shell thickness - delta_r = r2 - r1 - - if (delta_r/delta_min <= 1d0) then - ! use delta_min = 0.16 nm - coverage = aero_particle_organic_volume / v_delta - - !> calculate sigma for shell - sigma_shell = 0d0 - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & - * aero_data%sigma(i_shell_spec) / v_delta - end do - - aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & - coverage * sigma_shell - - else if (aero_particle_organic_volume < v_delta) then - ! use delta_min = 0.16 nm - coverage = aero_particle_organic_volume / v_delta - - !> calculate sigma for shell - sigma_shell = 0d0 - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & - * aero_data%sigma(i_shell_spec) / v_delta - end do - - aero_particle_varying_sigma = (1d0 - coverage) * sigma_core + & - coverage * sigma_shell + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + ! write(*,*) i_shell_spec, ' ,' , aero_particle%vol(i_shell_spec) , ' ,' + aero_particle_shell_volume = aero_particle_shell_volume & + + aero_particle%vol(i_shell_spec) + end do + + aero_particle_sigma_shell = 0d0 + + if (aero_particle_shell_volume < 1d-30) then + aero_particle_sigma_shell = 0d0 else - ! v_delta = v_organic - ! coverage = 1 - !> calculate sigma for shell - sigma_shell = 0d0 - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - sigmal_shell = sigma_shell + aero_particle%vol(i_shell_spec) & - * aero_data%sigma(i_shell_spec) / aero_particle_organic_volume - end do - - aero_particle_varying_sigma = sigma_shell + do n_shell_spec = 1, size(shell_spec) + i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) + aero_particle_sigma_shell = aero_particle_sigma_shell + aero_particle%vol(i_shell_spec) * & + aero_data%sigma(i_shell_spec) / aero_particle_shell_volume + end do end if + ! write(*,*) aero_particle_sigma_shell - write(*, *) aero_particle_varying_sigma - end function aero_particle_varying_sigma + end function aero_particle_sigma_shell !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1052,7 +992,7 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & !> Environment state. type(env_state_t), intent(in) :: env_state - integer, parameter :: CRIT_DIAM_MAX_ITER = 100 + integer, parameter :: CRIT_DIAM_MAX_ITER = 200 real(kind=dp) :: kappa, dry_diam, A, c4, c3, c0, d, f, df, dd integer :: i_newton @@ -1082,16 +1022,16 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & end if end do call warn_assert_msg(408545686, i_newton < CRIT_DIAM_MAX_ITER, & - "critical diameter Newton loop failed to converge") + "critical diameter of old Newton loop failed to converge") call warn_assert_msg(353290871, d >= dry_diam, & - "critical diameter Newton loop converged to invalid solution") + "critical diameter of old Newton loop converged to invalid solution") aero_particle_crit_diameter = d end function aero_particle_crit_diameter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Returns the critical diameter (m) for varying simga. + ! > Returns the critical diameter (m) for varying simga. real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& aero_particle, aero_data, env_state) @@ -1102,43 +1042,72 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& !> Environment state. type(env_state_t), intent(in) :: env_state - integer, parameter :: CRIT_DIAM_MAX_ITER = 100 + integer, parameter :: CRIT_DIAM_MAX_ITER = 200 - real(kind=dp) :: kappa, dry_diam, A, c4, c3, & - c0, d, f, df, dd + real(kind=dp) :: kappa, dry_diam, A, v_beta, c_1, c_b, c_3, c_4, c_5 + real(kind=dp) :: crit_diameter, d, v_delta, sigma, d_sigma + real(kind=dp) :: d_v_delta, R, d_R, f, df, dd, T integer :: i_newton + real(kind=dp) :: sigma_core, sigma_shell + ! real(kind=dp) :: sigma_core = 0.073d0 + ! real(kind=dp) :: sigma_shell = 0.030d0 + real(kind=dp) :: delta_min = 1.6d-10 - A = env_state_A(env_state) - varying_sigma = aero_particle_varying_sigma(aero_particle, aero_data) + sigma_core = aero_particle_sigma_core(aero_particle, aero_data) + sigma_shell = aero_particle_sigma_shell(aero_particle, aero_data) - dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) kappa = aero_particle_solute_kappa(aero_particle, aero_data) - if (kappa < 1d-30) then - ! bail out early for hydrophobic particles - aero_particle_crit_diameter_varying_sigma = dry_diam - return - end if - - c4 = - 3d0 * dry_diam**3 * kappa / (A * varying_sigma) - c3 = - dry_diam**3 * (2d0 - kappa) - c0 = dry_diam**6 * (1d0 - kappa) - - ! Newton's method for f(d) = d^6 + c4 d^4 + c3 d^3 + c0 - d = max(sqrt(-4d0 / 3d0 * c4), (-c3)**(1d0/3d0)) - do i_newton = 1,CRIT_DIAM_MAX_ITER - f = d**6 + c4 * d**4 + c3 * d**3 + c0 - df = 6 * d**5 + 4 * c4 * d**3 + 3 * c3 * d**2 - dd = f / df - d = d - dd - if (abs(dd / d) < 1d-14) then - exit - end if + dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) + A = env_state_A(env_state) + v_beta = aero_particle_organic_volume(aero_particle, aero_data) + + T = env_state%temp + + c_1 = 3 * dry_diam**3 * kappa / A + c_2= (2 - kappa) * dry_diam**3 + c_3 = (1 - kappa) * dry_diam**6 + c_4 = 2 * const%pi * delta_min ! d_2_v_delta_d + c_5 = v_beta * (sigma_shell - sigma_core) + + ! crit_diameter = aero_particle_crit_diameter(aero_particle, aero_data, env_state) + ! d = 5 * crit_diameter + d = 30 * dry_diam + + do i_newton = 1, CRIT_DIAM_MAX_ITER + v_delta = 4 * const%pi * ((d/2)**3 - (d/2 - delta_min)**3) / 3 + d_v_delta = 2 * const%pi * delta_min * (d - delta_min) + + if (v_beta > v_delta) then + sigma = sigma_shell + d_sigma = 0 + else if (v_beta < 1d-30) then + sigma = sigma_core + d_sigma = 0 + else + sigma = sigma_core + c_5 / v_delta + d_sigma = - c_5 * d_v_delta / v_delta**2 + end if + + R = sigma - d * d_sigma + d_R = d * c_5 * (v_delta * c_4 - 2 * d_v_delta**2) / v_delta**3 + + f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 + df = 6 * d**5 - c_1 * (4 * d**3 * R - d**4 * d_R) / R**2 - 3 * c_2 * d**2 + + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-12) then + exit + end if end do + + ! write(*,*) T, ' ,', dry_diam, ' ,', 30*dry_diam, ' ,', v_beta, ' ,', kappa, ' ,',& + ! i_newton, ' ,', sigma_shell, ' ,', sigma_core, ' ,', d, ' ,' ,abs(dd/d) call warn_assert_msg(408545686, i_newton < CRIT_DIAM_MAX_ITER, & - "critical diameter Newton loop failed to converge") + "critical diameter for new Newton loop failed to converge") call warn_assert_msg(353290871, d >= dry_diam, & - "critical diameter Newton loop converged to invalid solution") - aero_particle_crit_diameter_varying_sigma = d + "critical diameter for new Newton loop converged to invalid solution") + aero_particle_crit_diameter_sigma = d end function aero_particle_crit_diameter_varying_sigma @@ -1323,4 +1292,4 @@ end subroutine aero_particle_check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module pmc_aero_particle \ No newline at end of file +end module pmc_aero_particle From 0d518a21ca9f0ace75513171a6be8f7da1f76bdd Mon Sep 17 00:00:00 2001 From: Xu Date: Wed, 8 Mar 2023 10:34:58 -0600 Subject: [PATCH 28/33] update effective surface tension method in newton's loop and SSc --- src/aero_particle.F90 | 446 ++++++++++++++++++++++-------------------- 1 file changed, 234 insertions(+), 212 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 94c07b06f..8072245c4 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -768,179 +768,64 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, type(aero_data_t), intent(in) :: aero_data !> Environment state. type(env_state_t), intent(in) :: env_state - - real(kind=dp) :: kappa, dry_diam, A, sigma, crit_diam_varying_sigma - real(kind=dp) :: v_delta, v_beta !, v_inorg, r_inor, v_dry_shell - real(kind=dp) :: sigma_core, sigma_shell, c_beta - ! real(kind=dp) :: sigma_core = 0.073d0 - ! real(kind=dp) :: sigma_shell = 0.030d0 + + real(kind=dp) :: A, dry_diam, kappa, d, v_solid, r_solid + real(kind=dp) :: v_org, v_sol, r_core, v_delta + real(kind=dp) :: sigma_soluble, sigma_organic, delta_sigma real(kind=dp) :: delta_min = 1.6d-10 - sigma_core = aero_particle_sigma_core(aero_particle, aero_data) - sigma_shell = aero_particle_sigma_shell(aero_particle, aero_data) - A = env_state_A(env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - crit_diam_varying_sigma = aero_particle_crit_diameter_varying_sigma(& - aero_particle, aero_data, env_state) - v_beta = aero_particle_organic_volume(aero_particle, aero_data) - v_delta = 4d0 * const%pi * ((crit_diam_varying_sigma / 2d0)**3 - & - (crit_diam_varying_sigma / 2d0 - delta_min)**3) / 3d0 - kappa = aero_particle_solute_kappa(aero_particle, aero_data) - - if (kappa < 1d-30) then ! try to catch pure BC particles - aero_particle_crit_rel_humid_varying_sigma = exp(A * const%water_surf_eng / dry_diam) - else - if (v_beta > v_delta) then - sigma = sigma_shell - else if (v_beta < 1d-30) then - sigma = sigma_core - else - sigma = sigma_core + (v_beta/v_delta) * (sigma_shell - sigma_core) - end if - - aero_particle_crit_rel_humid_varying_sigma = (crit_diam_varying_sigma**3 & - - dry_diam**3) / (crit_diam_varying_sigma**3 - dry_diam**3 * (1 - kappa)) * & - exp(A * sigma / crit_diam_varying_sigma) - end if - - end function aero_particle_crit_rel_humid_varying_sigma - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Organic volume of a single species 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 - !> Species names to include in the volume. - character(len=6), allocatable :: org_spec(:) - integer :: i_org_spec - integer :: n_org_spec - - 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 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! > inrganic volume of a single species in the particle (m^3). - real(kind=dp) function aero_particle_inorganic_volume(aero_particle, aero_data) - - !> Particle. - type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data.ls - type(aero_data_t), intent(in) :: aero_data - !> Species names to include in the volume. - character(len=6), allocatable :: inorg_spec(:) - integer :: i_inorg_spec - integer :: n_inorg_spec - - inorg_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - "Na ", "Ca ", "OIN ", "BC "] - - aero_particle_inorganic_volume = 0d0 - - do n_inorg_spec = 1, size(inorg_spec) - i_inorg_spec = aero_data_spec_by_name(aero_data, inorg_spec(n_inorg_spec)) - aero_particle_inorganic_volume = aero_particle_inorganic_volume & - + aero_particle%vol(i_inorg_spec) - end do - - end function aero_particle_inorganic_volume - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Returns the sigma of core. - real(kind=dp) function aero_particle_sigma_core(aero_particle, aero_data) - - !> Aerosol particle. - type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. - type(aero_data_t), intent(in) :: aero_data - - !> core/organic species names to include in the volume. - character(len=6), allocatable :: core_spec(:) - integer :: i_core_spec - integer :: n_core_spec - - !> calculate core volume - core_spec = ["SO4 ", "NO3 ", "Cl ", "NH4 ", "CO3 ", & - "Na ", "Ca ", "OIN ", "BC ", "H2O "] - - aero_particle_core_volume = 0d0 - - do n_core_spec = 1, size(core_spec) - i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) - ! write(*, *) i_core_spec, ',', aero_particle%vol(i_core_spec), ',' - aero_particle_core_volume = aero_particle_core_volume & - + aero_particle%vol(i_core_spec) - end do - - aero_particle_sigma_core = 0d0 - do n_core_spec = 1, size(core_spec) - i_core_spec = aero_data_spec_by_name(aero_data, core_spec(n_core_spec)) - aero_particle_sigma_core = aero_particle_sigma_core + aero_particle%vol(i_core_spec) * & - aero_data%sigma(i_core_spec) / aero_particle_core_volume - end do - - end function aero_particle_sigma_core - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Returns the sigma of shell - real(kind=dp) function aero_particle_sigma_shell(aero_particle, aero_data) - - !> Aerosol particle. - type(aero_particle_t), intent(in) :: aero_particle - !> Aerosol data. - type(aero_data_t), intent(in) :: aero_data - - !> shell/org species names to include in the volume. - character(len=6), allocatable :: shell_spec(:) - integer :: i_shell_spec - integer :: n_shell_spec - - !> calculate organic volume - shell_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & - "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + d = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, env_state) + v_solid = aero_particle_solid_volume(aero_particle, aero_data) + r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) + v_org = aero_particle_organic_volume(aero_particle, aero_data) + v_sol = const%pi * d**3 / 6d0 - v_org - v_solid + + if (v_solid > 0d0) then + v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) + else + v_delta_solid = 0d0 + end if - aero_particle_shell_volume = 0d0 + r_core = ((3d0 * (v_sol + v_solid)) / (4d0 * const%pi))**(1d0 / 3d0) + v_delta = (4d0 * const%pi / 3) * ((r_core + delta_min)**3 - (r_core)**3) - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - ! write(*,*) i_shell_spec, ' ,' , aero_particle%vol(i_shell_spec) , ' ,' - aero_particle_shell_volume = aero_particle_shell_volume & - + aero_particle%vol(i_shell_spec) - end do + sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) + sigma_organic = aero_particle_sigma_organic(aero_particle, aero_data) - aero_particle_sigma_shell = 0d0 + delta_sigma = sigma_organic - sigma_soluble - if (aero_particle_shell_volume < 1d-30) then - aero_particle_sigma_shell = 0d0 + if (d == dry_diam) then + if (v_delta_solid > (v_org + v_sol)) then + sigma = v_org * sigma_organic / v_delta_solid + & + v_sol * sigma_soluble / v_delta_solid + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + else + if (v_org > v_delta_solid) then + sigma = sigma_organic + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + else + sigma = sigma_soluble + v_org * delta_sigma / v_delta_solid + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + end if + end if else - do n_shell_spec = 1, size(shell_spec) - i_shell_spec = aero_data_spec_by_name(aero_data, shell_spec(n_shell_spec)) - aero_particle_sigma_shell = aero_particle_sigma_shell + aero_particle%vol(i_shell_spec) * & - aero_data%sigma(i_shell_spec) / aero_particle_shell_volume - end do - end if - - ! write(*,*) aero_particle_sigma_shell - - end function aero_particle_sigma_shell + if (v_sol > v_delta_solid) then + if (v_org > v_delta) then + sigma = sigma_organic + aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & + (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + else + sigma = sigma_soluble + v_org * delta_sigma / v_delta + aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & + (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + end if + end if + end if + end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the critical diameter (m). @@ -992,7 +877,7 @@ real(kind=dp) function aero_particle_crit_diameter(aero_particle, & !> Environment state. type(env_state_t), intent(in) :: env_state - integer, parameter :: CRIT_DIAM_MAX_ITER = 200 + integer, parameter :: CRIT_DIAM_MAX_ITER = 100 real(kind=dp) :: kappa, dry_diam, A, c4, c3, c0, d, f, df, dd integer :: i_newton @@ -1028,7 +913,6 @@ 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) for varying simga. @@ -1042,75 +926,213 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& !> Environment state. type(env_state_t), intent(in) :: env_state - integer, parameter :: CRIT_DIAM_MAX_ITER = 200 - - real(kind=dp) :: kappa, dry_diam, A, v_beta, c_1, c_b, c_3, c_4, c_5 - real(kind=dp) :: crit_diameter, d, v_delta, sigma, d_sigma - real(kind=dp) :: d_v_delta, R, d_R, f, df, dd, T + integer, parameter :: CRIT_DIAM_MAX_ITER = 100 integer :: i_newton - real(kind=dp) :: sigma_core, sigma_shell - ! real(kind=dp) :: sigma_core = 0.073d0 - ! real(kind=dp) :: sigma_shell = 0.030d0 real(kind=dp) :: delta_min = 1.6d-10 + real(kind=dp) :: kappa, dry_diam, A, T + real(kind=dp) :: sigma_soluble, sigma_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) :: sigma, d_sigma, d_2_sigma, R, d_R, f, df, dd + real(kind=dp) :: c_1, c_2, c_3, c_4, c_5 - sigma_core = aero_particle_sigma_core(aero_particle, aero_data) - sigma_shell = aero_particle_sigma_shell(aero_particle, aero_data) - + A = env_state_A(env_state) + T = env_state%temp kappa = aero_particle_solute_kappa(aero_particle, aero_data) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) - A = env_state_A(env_state) - v_beta = aero_particle_organic_volume(aero_particle, aero_data) + sigma_organic = aero_particle_sigma_organic(aero_particle, aero_data) + v_solid = aero_particle_solid_volume(aero_particle, aero_data) - T = env_state%temp + v_org = aero_particle_organic_volume(aero_particle, aero_data) + r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) - c_1 = 3 * dry_diam**3 * kappa / A - c_2= (2 - kappa) * dry_diam**3 - c_3 = (1 - kappa) * dry_diam**6 - c_4 = 2 * const%pi * delta_min ! d_2_v_delta_d - c_5 = v_beta * (sigma_shell - sigma_core) + if (v_solid > 0d0) then + v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) + else + v_delta_insol = 0d0 + end if - ! crit_diameter = aero_particle_crit_diameter(aero_particle, aero_data, env_state) - ! d = 5 * crit_diameter - d = 30 * dry_diam + 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 = 30d0 * dry_diam do i_newton = 1, CRIT_DIAM_MAX_ITER - v_delta = 4 * const%pi * ((d/2)**3 - (d/2 - delta_min)**3) / 3 - d_v_delta = 2 * const%pi * delta_min * (d - delta_min) - - if (v_beta > v_delta) then - sigma = sigma_shell - d_sigma = 0 - else if (v_beta < 1d-30) then - sigma = sigma_core - d_sigma = 0 - else - sigma = sigma_core + c_5 / v_delta - d_sigma = - c_5 * d_v_delta / v_delta**2 - end if - R = sigma - d * d_sigma - d_R = d * c_5 * (v_delta * c_4 - 2 * d_v_delta**2) / v_delta**3 - - f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 - df = 6 * d**5 - c_1 * (4 * d**3 * R - d**4 * d_R) / R**2 - 3 * c_2 * d**2 - - dd = f / df - d = d - dd - if (abs(dd / d) < 1d-12) then - exit + sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) + c_5 = v_org * (sigma_organic - sigma_soluble) + v_sol = const%pi * d**3 / 6 - v_solid - v_org ! volume of soluble inorganics + water + + if (v_sol > v_delta_solid) then ! fully cover BC/OIN + r_core = ((3d0 * (v_sol + v_solid))/(4 * const%pi))**(1d0 / 3d0) ! BC/OIN + soluble inorganics + v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) + if (v_org > v_delta) then ! soluble inorganics and BC/OIN as core, and organics can fully cover + sigma = sigma_organic + d_sigma = 0 + f = d**6 - c_1 * d**4 / sigma - c_2 * d**3 + c_3 + df = 6 * d**5 - 4 * c_1 * d**3 / sigma - 3 * c_2 * d**2 + else ! soluble inorganics and BC/OIN as core, but organics cannot fully cover + sigma = sigma_soluble + c_5 / v_delta + d_sigma = - c_5 * d_v_delta / v_delta**2 + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_sigma = - c_5 * d_v_delta / v_delta**2 + R = sigma - d * d_sigma + d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_sigma + f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 + df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 + end if + + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-12) then + exit + end if + + else ! cannot fully cover BC/OIN + d = dry_diam + exit end if + end do - ! write(*,*) T, ' ,', dry_diam, ' ,', 30*dry_diam, ' ,', v_beta, ' ,', kappa, ' ,',& - ! i_newton, ' ,', sigma_shell, ' ,', sigma_core, ' ,', d, ' ,' ,abs(dd/d) 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_sigma = d + aero_particle_crit_diameter_varying_sigma = d end function aero_particle_crit_diameter_varying_sigma +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Organic volume of a single species 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 + !> Species names to include in the volume. + character(len=6), allocatable :: org_spec(:) + integer :: i_org_spec + integer :: n_org_spec + + 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 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! > soluble inrganic volume of a single species 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.ls + type(aero_data_t), intent(in) :: aero_data + !> Species names to include in the volume. + character(len=6), allocatable :: solid_spec(:) + integer :: i_solid_spec + integer :: n_solid_spec + + 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 sigma of sol. + real(kind=dp) function aero_particle_sigma_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 + + character(len=6), allocatable :: sol_spec(:) + integer :: i_sol_spec + integer :: n_sol_spec + real(kind=dp) :: d, aero_particle_soluble_volume + + !> calculate sol volume + 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_sigma_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_sigma_soluble = aero_particle_sigma_soluble & + + aero_particle%vol(i_sol_spec) * & + aero_data%sigma(i_sol_spec) / v_sol + end do + aero_particle_sigma_soluble = aero_particle_sigma_soluble + v_water & + * const%water_surf_eng / v_sol + + end function aero_particle_sigma_soluble + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the sigma of shell + real(kind=dp) function aero_particle_sigma_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 + + !> org/org species names to include in the volume. + character(len=6), allocatable :: org_spec(:) + integer :: i_org_spec + integer :: n_org_spec + + !> calculate organic volume + org_spec = ["MSA ", "ARO1 ", "ARO2 ", "ALK1 ", "OLE1 ", & + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] + + org_volume = aero_particle_organic_volume(aero_particle, aero_data) + aero_particle_sigma_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_sigma_organic = aero_particle_sigma_organic + aero_particle%vol(i_org_spec) * & + aero_data%sigma(i_org_spec) / org_volume + end do + end function aero_particle_sigma_organic !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Coagulate two particles together to make a new one. The new From 84fdc25c4348c32a4953ce8d82e92b4ec9691a96 Mon Sep 17 00:00:00 2001 From: Xu Date: Mon, 27 Mar 2023 12:26:30 -0500 Subject: [PATCH 29/33] fix error for v_org = 0 --- src/aero_particle.F90 | 34 +++++++++++++++++++++++----------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 8072245c4..e473e9a95 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -788,42 +788,54 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, else v_delta_solid = 0d0 end if - +· r_core = ((3d0 * (v_sol + v_solid)) / (4d0 * const%pi))**(1d0 / 3d0) v_delta = (4d0 * const%pi / 3) * ((r_core + delta_min)**3 - (r_core)**3) sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) sigma_organic = aero_particle_sigma_organic(aero_particle, aero_data) - delta_sigma = sigma_organic - sigma_soluble + delta_sigma = sigma_organic - sigma_soluble· if (d == dry_diam) then if (v_delta_solid > (v_org + v_sol)) then sigma = v_org * sigma_organic / v_delta_solid + & v_sol * sigma_soluble / v_delta_solid - aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + ! write(*,*) "a", aero_particle_crit_rel_humid_varying_sigma else if (v_org > v_delta_solid) then sigma = sigma_organic aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) + ! write(*,*) "b", aero_particle_crit_rel_humid_varying_sigma else - sigma = sigma_soluble + v_org * delta_sigma / v_delta_solid + if (v_org < 1d-30) then + ! sigma_organic = 0d0 + sigma = sigma_soluble + else + sigma = sigma_soluble + v_org * delta_sigma / v_delta_solid + end if aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) end if end if else - if (v_sol > v_delta_solid) then if (v_org > v_delta) then sigma = sigma_organic aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) else - sigma = sigma_soluble + v_org * delta_sigma / v_delta + if (v_org < 1d-30) then + ! sigma_organic = 0d0 + sigma = sigma_soluble + else + sigma = sigma_soluble + v_org * delta_sigma / v_delta + end if aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) end if end if - end if + + write(*,*) aero_particle_crit_rel_humid_varying_sigma end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -965,7 +977,7 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& c_5 = v_org * (sigma_organic - sigma_soluble) v_sol = const%pi * d**3 / 6 - v_solid - v_org ! volume of soluble inorganics + water - if (v_sol > v_delta_solid) then ! fully cover BC/OIN + if (v_sol + v_org > v_delta_solid) then ! fully cover BC/OIN r_core = ((3d0 * (v_sol + v_solid))/(4 * const%pi))**(1d0 / 3d0) ! BC/OIN + soluble inorganics v_delta = (4d0 * const%pi / 3d0) * ((r_core + delta_min)**3 - (r_core)**3) if (v_org > v_delta) then ! soluble inorganics and BC/OIN as core, and organics can fully cover @@ -987,7 +999,7 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& dd = f / df d = d - dd - if (abs(dd / d) < 1d-12) then + if (abs(dd / d) < 1d-11) then exit end if @@ -997,13 +1009,13 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& end if end do - + ! write(*,*) abs(dd), abs(dd / d) 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_varying_sigma = d - + end function aero_particle_crit_diameter_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 5cc4777ca2cb5e86e6b4f6bb3f7cf5b3ebb0a611 Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 14 Apr 2023 19:43:03 -0500 Subject: [PATCH 30/33] fix defination warning --- src/aero_particle.F90 | 94 +++++++++++++++---------------------------- 1 file changed, 32 insertions(+), 62 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index e473e9a95..b5f31021b 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -770,8 +770,8 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, type(env_state_t), intent(in) :: env_state real(kind=dp) :: A, dry_diam, kappa, d, v_solid, r_solid - real(kind=dp) :: v_org, v_sol, r_core, v_delta - real(kind=dp) :: sigma_soluble, sigma_organic, delta_sigma + real(kind=dp) :: v_org, v_sol, r_core, v_delta, v_delta_solid + real(kind=dp) :: sigma_soluble, sigma_organic, delta_sigma, sigma real(kind=dp) :: delta_min = 1.6d-10 A = env_state_A(env_state) @@ -788,54 +788,35 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, else v_delta_solid = 0d0 end if -· + r_core = ((3d0 * (v_sol + v_solid)) / (4d0 * const%pi))**(1d0 / 3d0) v_delta = (4d0 * const%pi / 3) * ((r_core + delta_min)**3 - (r_core)**3) sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) sigma_organic = aero_particle_sigma_organic(aero_particle, aero_data) - delta_sigma = sigma_organic - sigma_soluble· + delta_sigma = sigma_organic - sigma_soluble if (d == dry_diam) then - if (v_delta_solid > (v_org + v_sol)) then - sigma = v_org * sigma_organic / v_delta_solid + & - v_sol * sigma_soluble / v_delta_solid - aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) - ! write(*,*) "a", aero_particle_crit_rel_humid_varying_sigma - else - if (v_org > v_delta_solid) then - sigma = sigma_organic - aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) - ! write(*,*) "b", aero_particle_crit_rel_humid_varying_sigma - else - if (v_org < 1d-30) then - ! sigma_organic = 0d0 - sigma = sigma_soluble - else - sigma = sigma_soluble + v_org * delta_sigma / v_delta_solid - end if - aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) - end if - end if + sigma = v_org * sigma_organic / v_delta_solid + & + v_sol * sigma_soluble / v_delta_solid + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) else - if (v_org > v_delta) then - sigma = sigma_organic - aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & - (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + if (v_org > v_delta) then + sigma = sigma_organic + else + if (v_org < 1d-30) then + ! sigma_organic = 0d0 + sigma = sigma_soluble else - if (v_org < 1d-30) then - ! sigma_organic = 0d0 - sigma = sigma_soluble - else - sigma = sigma_soluble + v_org * delta_sigma / v_delta - end if - aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & - (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + sigma = sigma_soluble + v_org * delta_sigma / v_delta end if end if + aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & + (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + end if - write(*,*) aero_particle_crit_rel_humid_varying_sigma + ! write(*,*) aero_particle_crit_rel_humid_varying_sigma end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -961,7 +942,7 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& if (v_solid > 0d0) then v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) else - v_delta_insol = 0d0 + v_delta_solid = 0d0 end if c_1 = 3d0 * dry_diam**3 * kappa / A @@ -1026,11 +1007,8 @@ 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 - !> Species names to include in the volume. - character(len=6), allocatable :: org_spec(:) - integer :: i_org_spec - integer :: n_org_spec - + 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 "] @@ -1052,13 +1030,10 @@ real(kind=dp) function aero_particle_solid_volume(aero_particle, aero_data) type(aero_particle_t), intent(in) :: aero_particle !> Aerosol data.ls type(aero_data_t), intent(in) :: aero_data - !> Species names to include in the volume. - character(len=6), allocatable :: solid_spec(:) - integer :: i_solid_spec - integer :: n_solid_spec - + 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) @@ -1081,14 +1056,12 @@ real(kind=dp) function aero_particle_sigma_soluble(aero_particle, & !> Environment state. type(env_state_t), intent(in) :: env_state - character(len=6), allocatable :: sol_spec(:) - integer :: i_sol_spec - integer :: n_sol_spec + integer :: i_sol_spec, n_sol_spec real(kind=dp) :: d, aero_particle_soluble_volume - - !> calculate sol 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 "] + "Na ", "Ca "] v_solid = aero_particle_solid_volume(aero_particle, aero_data) v_org = aero_particle_organic_volume(aero_particle, aero_data) @@ -1126,14 +1099,11 @@ real(kind=dp) function aero_particle_sigma_organic(aero_particle, aero_data) !> Aerosol data. type(aero_data_t), intent(in) :: aero_data - !> org/org species names to include in the volume. - character(len=6), allocatable :: org_spec(:) - integer :: i_org_spec - integer :: n_org_spec - - !> calculate organic volume + 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 "] + "API1 ", "API2 ", "LIM1 ", "LIM2 ", "OC "] org_volume = aero_particle_organic_volume(aero_particle, aero_data) aero_particle_sigma_organic = 0d0 From b65199fe6201aa78b35c963a5872c62d8dbe6fe2 Mon Sep 17 00:00:00 2001 From: Xu Date: Sat, 8 Jul 2023 22:26:29 -0500 Subject: [PATCH 31/33] final effective surface tension update --- scenarios/1_urban_plume/aero_data.dat | 4 +- src/aero_particle.F90 | 173 ++++++++++++++------------ src/aero_state.F90 | 46 +------ 3 files changed, 93 insertions(+), 130 deletions(-) diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index f86a8afc4..4a771ead8 100644 --- a/scenarios/1_urban_plume/aero_data.dat +++ b/scenarios/1_urban_plume/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) sigma +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) surface tension (J m^{-2}) SO4 1800 0 96d-3 0.65 0.073 NO3 1800 0 62d-3 0.65 0.073 Cl 2200 0 35.5d-3 1.28 0.073 @@ -18,4 +18,4 @@ Ca 2600 0 40d-3 0.53 OIN 2600 0 1d-3 0.1 0 OC 1000 0 1d-3 0.001 0.03 BC 1800 0 1d-3 0 0 -H2O 1000 0 18d-3 0 0.073 +H2O 1000 0 18d-3 0 0.073 \ No newline at end of file diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index b5f31021b..bd8615bea 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -768,56 +768,24 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, 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, v_solid, r_solid - real(kind=dp) :: v_org, v_sol, r_core, v_delta, v_delta_solid - real(kind=dp) :: sigma_soluble, sigma_organic, delta_sigma, sigma - real(kind=dp) :: delta_min = 1.6d-10 + real(kind=dp) :: A, dry_diam, kappa, d, sigma + A = env_state_A(env_state) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) kappa = aero_particle_solute_kappa(aero_particle, aero_data) - d = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, env_state) - v_solid = aero_particle_solid_volume(aero_particle, aero_data) - r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) - v_org = aero_particle_organic_volume(aero_particle, aero_data) - v_sol = const%pi * d**3 / 6d0 - v_org - v_solid + d = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, env_state, sigma) - if (v_solid > 0d0) then - v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) - else - v_delta_solid = 0d0 - end if - - r_core = ((3d0 * (v_sol + v_solid)) / (4d0 * const%pi))**(1d0 / 3d0) - v_delta = (4d0 * const%pi / 3) * ((r_core + delta_min)**3 - (r_core)**3) - - sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) - sigma_organic = aero_particle_sigma_organic(aero_particle, aero_data) + crit_rhs = aero_particle_crit_rel_humid(aero_particle, & + aero_data, env_state) - delta_sigma = sigma_organic - sigma_soluble - if (d == dry_diam) then - sigma = v_org * sigma_organic / v_delta_solid + & - v_sol * sigma_soluble / v_delta_solid - aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / d) - else - if (v_org > v_delta) then - sigma = sigma_organic - else - if (v_org < 1d-30) then - ! sigma_organic = 0d0 - sigma = sigma_soluble - else - sigma = sigma_soluble + v_org * delta_sigma / v_delta - end if - end if + aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / dry_diam) + else aero_particle_crit_rel_humid_varying_sigma = (d**3 - dry_diam**3) / & - (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) - end if - - ! write(*,*) aero_particle_crit_rel_humid_varying_sigma - + (d**3 - dry_diam**3 * (1d0 - kappa)) * exp(A * sigma / d) + end if + end function aero_particle_crit_rel_humid_varying_sigma !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -910,7 +878,7 @@ end function aero_particle_crit_diameter ! > Returns the critical diameter (m) for varying simga. real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& - aero_particle, aero_data, env_state) + aero_particle, aero_data, env_state, sigma) !> Aerosol particle. type(aero_particle_t), intent(in) :: aero_particle @@ -919,6 +887,8 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& !> Environment state. type(env_state_t), intent(in) :: env_state + real(kind=dp), intent(out) :: sigma + integer, parameter :: CRIT_DIAM_MAX_ITER = 100 integer :: i_newton real(kind=dp) :: delta_min = 1.6d-10 @@ -926,23 +896,21 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& real(kind=dp) :: sigma_soluble, sigma_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) :: sigma, d_sigma, d_2_sigma, R, d_R, f, df, dd + real(kind=dp) :: d_sigma, d_2_sigma, R, d_R, f, df, dd real(kind=dp) :: c_1, c_2, c_3, c_4, c_5 - + A = env_state_A(env_state) T = env_state%temp kappa = aero_particle_solute_kappa(aero_particle, aero_data) dry_diam = aero_particle_dry_diameter(aero_particle, aero_data) sigma_organic = aero_particle_sigma_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) - r_solid = ((3d0 * v_solid) / (4d0 * const%pi))**(1d0 / 3d0) - if (v_solid > 0d0) then - v_delta_solid = (4d0 * const%pi / 3d0) * ((r_solid + delta_min)**3 - (r_solid)**3) - else - v_delta_solid = 0d0 + if (kappa < 1d-30) then + ! bail out early for hydrophobic particles + aero_particle_crit_diameter_varying_sigma = dry_diam + return end if c_1 = 3d0 * dry_diam**3 * kappa / A @@ -950,47 +918,86 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& c_3 = (1d0 - kappa) * dry_diam**6 c_4 = 2d0 * const%pi * delta_min ! d_2_v_delta_d - d = 30d0 * dry_diam - - do i_newton = 1, CRIT_DIAM_MAX_ITER - - sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) - c_5 = v_org * (sigma_organic - sigma_soluble) - v_sol = const%pi * d**3 / 6 - v_solid - v_org ! volume of soluble inorganics + water - - if (v_sol + v_org > v_delta_solid) then ! fully cover BC/OIN - r_core = ((3d0 * (v_sol + v_solid))/(4 * const%pi))**(1d0 / 3d0) ! BC/OIN + soluble inorganics + d = 50*sqrt(4d0 / 3d0 * c_1) + + if (v_solid == 0d0) then + do i_newton = 1, CRIT_DIAM_MAX_ITER + sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) + c_5 = v_org * (sigma_organic - sigma_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 ! soluble inorganics and BC/OIN as core, and organics can fully cover + if (v_org > v_delta) then sigma = sigma_organic - d_sigma = 0 f = d**6 - c_1 * d**4 / sigma - c_2 * d**3 + c_3 - df = 6 * d**5 - 4 * c_1 * d**3 / sigma - 3 * c_2 * d**2 - else ! soluble inorganics and BC/OIN as core, but organics cannot fully cover - sigma = sigma_soluble + c_5 / v_delta - d_sigma = - c_5 * d_v_delta / v_delta**2 - d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) - d_sigma = - c_5 * d_v_delta / v_delta**2 - R = sigma - d * d_sigma - d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) - d_R = - d * d_2_sigma - f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 - df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 - end if - + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / sigma - 3d0 * c_2 * d**2 + else + if (v_org == 0) then + sigma = sigma_soluble + f = d**6 - c_1 * d**4 / sigma - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / sigma - 3d0 * c_2 * d**2 + else + sigma = sigma_soluble + c_5 / v_delta + d_sigma = - c_5 * d_v_delta / v_delta**2 + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_sigma = - c_5 * d_v_delta / v_delta**2 + R = sigma - d * d_sigma + d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_sigma + f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 + df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 + 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 + sigma_soluble = aero_particle_sigma_soluble(aero_particle, aero_data, env_state, d) + c_5 = v_org * (sigma_organic - sigma_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 + sigma = sigma_organic + f = d**6 - c_1 * d**4 / sigma - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / sigma - 3d0 * c_2 * d**2 + else + if (v_org == 0) then + sigma = sigma_soluble + f = d**6 - c_1 * d**4 / sigma - c_2 * d**3 + c_3 + df = 6d0 * d**5 - 4d0 * c_1 * d**3 / sigma - 3d0 * c_2 * d**2 + else + sigma = sigma_soluble + c_5 / v_delta + d_sigma = - c_5 * d_v_delta / v_delta**2 + d_v_delta = 2d0 * const%pi * delta_min * (d - delta_min) + d_sigma = - c_5 * d_v_delta / v_delta**2 + R = sigma - d * d_sigma + d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) + d_R = - d * d_2_sigma + f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 + df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 + end if + end if + dd = f / df + d = d - dd + if (abs(dd / d) < 1d-11) then + exit + end if + else + sigma = (v_org * sigma_organic + v_sol * sigma_soluble) / v_delta_solid + d = dry_diam + exit + end if + end do + end if - else ! cannot fully cover BC/OIN - d = dry_diam - exit - end if - - end do - ! write(*,*) abs(dd), abs(dd / d) 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, & diff --git a/src/aero_state.F90 b/src/aero_state.F90 index b926a5e29..6b6e3cae4 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1120,50 +1120,6 @@ function aero_state_volumes(aero_state, aero_data, include, exclude) end function aero_state_volumes -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> calculate organic species volume fraction - !> might be useful later. - -! function aero_state_organic_volumes(aero_state, aero_data) - -! !> Aerosol state. -! type(aero_state_t), intent(in) :: aero_state -! !> Aerosol data. -! type(aero_data_t), optional, intent(in) :: aero_data -! !> Species names to include in the mass. -! character(len=*) :: include(:) - -! !> Return volumes array (m^3). -! real(kind=dp) :: aero_state_organic_volumes(aero_state_n_part(aero_state)) - -! logical, allocatable :: use_species(:) -! integer :: i_name_organics, i_spec_organics - -! include = ["MSA ", "ARO1", "ARO2", "ALK1", "OLE1", & -! "API1", "API2", "LIM1", "LIM2", "OC "] - -! allocate(use_species(aero_data_n_spec(aero_data))) - -! use_species = .false. -! do i_name_organics = 1, size(include) -! i_spec_organics = aero_data_spec_by_name(aero_data, include(i_name_organics)) -! call assert_msg(111852070, i_spec_organics > 0, & -! "unknown species: " // trim(include(i_name_organics))) -! use_species(i_spec_organics) = .true. -! end do - -! aero_state_organic_volumes = 0d0 ! v_beta -! do i_spec_organics = 1,aero_data_n_spec(aero_data) -! if (use_species(i_spec_organics)) then -! aero_state_organic_volumes = aero_state_organic_volumes & -! + aero_particle_organic_species_volume( & -! aero_state%apa%particle(1:aero_state_n_part(aero_state)), & -! i_spec_organics) -! end if -! end do - -! end function aero_state_organic_volumes - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the masses of all particles. @@ -1558,7 +1514,7 @@ function aero_state_crit_rel_humids_varying_sigma(aero_state, aero_data, env_sta !> Return value. real(kind=dp) :: aero_state_crit_rel_humids_varying_sigma(aero_state_n_part(aero_state)) - + integer :: i_part do i_part = 1,aero_state_n_part(aero_state) From 210d246bd6f282c86e7bed448305dcaff207fc3a Mon Sep 17 00:00:00 2001 From: Xu Date: Fri, 25 Aug 2023 14:56:00 -0500 Subject: [PATCH 32/33] delete uncessary variable --- src/aero_particle.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index bd8615bea..76d566530 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -776,9 +776,6 @@ real(kind=dp) function aero_particle_crit_rel_humid_varying_sigma(aero_particle, kappa = aero_particle_solute_kappa(aero_particle, aero_data) d = aero_particle_crit_diameter_varying_sigma(aero_particle, aero_data, env_state, sigma) - crit_rhs = aero_particle_crit_rel_humid(aero_particle, & - aero_data, env_state) - if (d == dry_diam) then aero_particle_crit_rel_humid_varying_sigma = exp(A * sigma / dry_diam) else From 7dde0c55b25ed56b0ad533068bb8948f6b09d993 Mon Sep 17 00:00:00 2001 From: Xu Date: Tue, 12 Dec 2023 16:55:10 -0600 Subject: [PATCH 33/33] change f(d) --- src/aero_particle.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index 76d566530..51cec694d 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -941,8 +941,9 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& R = sigma - d * d_sigma d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) d_R = - d * d_2_sigma - f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 - df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 + 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 @@ -978,8 +979,9 @@ real(kind=dp) function aero_particle_crit_diameter_varying_sigma(& R = sigma - d * d_sigma d_2_sigma = (c_5 / v_delta**3) * (2d0 * d_v_delta**2 - v_delta * c_4) d_R = - d * d_2_sigma - f = d**6 - c_1 * d**4 / R - c_2 * d**3 + c_3 - df = 6d0 * d**5 - c_1 * (4d0 * d**3 * R - d**4 * d_R) / R**2 - 3d0 * c_2 * d**2 + 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