From 0f786d42c04e955859c43c5a75723dc8e1dbfce6 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sat, 7 Jun 2025 10:59:05 -0500 Subject: [PATCH 01/28] adding scaffolding for modal run --- src/run_modal.F90 | 143 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 src/run_modal.F90 diff --git a/src/run_modal.F90 b/src/run_modal.F90 new file mode 100644 index 000000000..6804e351b --- /dev/null +++ b/src/run_modal.F90 @@ -0,0 +1,143 @@ +!> \file +!> The pmc_run_modal module. + +!> 1D modal simulation. + +module pmc_run_modal + + use pmc_util + use pmc_aero_dist + use pmc_scenario + use pmc_env_state + use pmc_aero_data + use pmc_output + use pmc_bin_grid + use pmc_aero_binned + use pmc_gas_data + use pmc_gas_state + use pmc_constants + + !> Options controlling the operation of run_modal() + type run_modal_opt_t + !> Final time (s). + real(kind=dp) :: t_max + !> Timestep (s). + real(kind=dp) :: del_t + !> Output interval (0 disables) (s). + real(kind=dp) :: t_output + !> Progress interval (0 disables) (s). + real(kind=dp) :: t_progress + !> Whether to do emissions and dilution. + logical :: do_emissions + !> Output prefix. + character(len=300) :: prefix + !> UUID of the simulation. + character(len=PMC_UUID_LEN) :: uuid + end type run_modal_opt_t + +contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Run a modal simulation. + subroutine run_modal(aero_data, aero_dist, scenario, env_state, & + gas_data, bin_grid, run_modal_opt) + + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Aerosol distribution. + type(aero_dist_t), intent(inout) :: aero_dist + !> Environment data. + type(scenario_t), intent(inout) :: scenario + !> Environment state. + type(env_state_t), intent(inout) :: env_state + !> Gas data. + type(gas_data_t), intent(in) :: gas_data + !> Bin grid. + type(bin_grid_t), intent(in) :: bin_grid + !> Modal options. + type(run_modal_opt_t), intent(in) :: run_modal_opt + + type(env_state_t) :: old_env_state + type(gas_state_t) :: gas_state + type(aero_binned_t) :: aero_binned + + real(kind=dp) time, last_output_time, last_progress_time, removed + + integer i, i_time, n_time, i_summary, i_mode, mode + logical do_output, do_progress + + call check_time_multiple("t_max", run_modal_opt%t_max, & + "del_t", run_modal_opt%del_t) + call check_time_multiple("t_output", run_modal_opt%t_output, & + "del_t", run_modal_opt%del_t) + call check_time_multiple("t_progress", run_modal_opt%t_progress, & + "del_t", run_modal_opt%del_t) + + if (aero_data_n_spec(aero_data) /= 1) then + call die_msg(927384615, & + 'run_modal() can only use one aerosol species') + end if + + ! output data structure + call gas_state_set_size(gas_state, gas_data_n_spec(gas_data)) + + ! Initialize time. + time = 0d0 + last_output_time = 0d0 + last_progress_time = 0d0 + i_summary = 1 + removed = 0d0 + + ! initial output + call check_event(time, run_modal_opt%del_t, run_modal_opt%t_output, & + last_output_time, do_output) + if (do_output) then + call aero_binned_add_aero_dist(aero_binned, bin_grid, aero_data, & + aero_dist) + call output_modal(run_modal_opt%prefix, aero_binned, aero_dist, aero_data, & + env_state, gas_data, gas_state, bin_grid, scenario, i_summary, time, & + run_modal_opt%del_t, run_modal_opt%uuid) + call aero_binned_zero(aero_binned) + end if + + ! Main time-stepping loop + n_time = nint(run_modal_opt%t_max / run_modal_opt%del_t) + do i_time = 1,n_time + + time = run_modal_opt%t_max * real(i_time, kind=dp) & + / real(n_time, kind=dp) + + old_env_state = env_state + + call scenario_update_env_state(scenario, env_state, time) + call scenario_update_gas_state(scenario, run_modal_opt%del_t, & + env_state, old_env_state, gas_data, gas_state) + call scenario_update_aero_modes(aero_dist, run_modal_opt%del_t, & + env_state, aero_data%density(1), scenario, removed) + + call check_event(time, run_modal_opt%del_t, run_modal_opt%t_output, & + last_output, do_output) + if (do_output) then + i_summary = i_summary + 1 + call aero_binned_add_aero_dist(aero_binned, bin_grid, aero_data, & + aero_dist) + call output_modal(run_modal_opt%prefix, aero_binned, aero_dist, aero_data, & + env_state, gas_data, gas_state, bin_grid, scenario,i_summary, time, & + run_modal_opt%del_t, run_modal_opt%uuid) + call aero_binned_zero(aero_binned) + end if + + call check_event(time, run_modal_opt%del_t, run_modal_opt%t_progress, & + last_progress, do_progress) + if (do_progress) then + write(*, '(a6,a8)') 'step', 'time' + write(*, '(i6,f8.1)') i_time, time + end if + end do + + end subroutine run_modal + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end module pmc_run_modal \ No newline at end of file From 9e603335f3f16d0756a5ece49613758ea9bf0e40 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sat, 7 Jun 2025 11:07:45 -0500 Subject: [PATCH 02/28] add modal run to driver --- src/partmc.F90 | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/partmc.F90 b/src/partmc.F90 index b873970f2..24360d6d8 100644 --- a/src/partmc.F90 +++ b/src/partmc.F90 @@ -300,6 +300,8 @@ subroutine partmc_run(spec_name) call partmc_exact(file) elseif (trim(run_type) == 'sectional') then call partmc_sect(file) + elseif (trim(run_type) == 'modal') then + call partmc_modal(file) else call die_msg(719261940, "unknown run_type: " // trim(run_type)) end if @@ -782,4 +784,73 @@ end subroutine partmc_sect !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!> Run a modal code simulation. + subroutine partmc_modal(file) + + !> Spec file + type(spec_file_t), intent(inout) :: file + + type(run_modal_opt_t) :: run_modal_opt + type(aero_data_t) :: aero_data + type(aero_dist_t) :: aero_dist_init + type(scenario_t) :: scenario + type(env_state_t) :: env_state + type(bin_grid_t) :: bin_grid + type(gas_data_t) :: gas_data + character(len=PMC_MAX_FILENAME_LEN) :: sub_filename + type(spec_file_t) :: sub_file + + ! only serial code here + if (pmc_mpi_rank() /= 0) then + return + end if + + call spec_file_read_string(file, 'output_prefix', run_modal_opt%prefix) + + call spec_file_read_real(file, 't_max', run_modal_opt%t_max) + call spec_file_read_real(file, 'del_t', run_modal_opt%del_t) + call spec_file_read_real(file, 't_output', run_modal_opt%t_output) + call spec_file_read_real(file, 't_progress', run_modal_opt%t_progress) + + call spec_file_read_radius_bin_grid(file, bin_grid) + + call spec_file_read_string(file, 'gas_data', sub_filename) + call spec_file_open(sub_filename, sub_file) + call spec_file_read_gas_data(sub_file, gas_data) + call spec_file_close(sub_file) + + call spec_file_read_string(file, 'aerosol_data', sub_filename) + call spec_file_open(sub_filename, sub_file) + call spec_file_read_aero_data(sub_file, aero_data) + call spec_file_close(sub_file) + + call spec_file_read_fractal(file, aero_data%fractal) + + call spec_file_read_string(file,'aerosol_init', sub_filename) + call spec_file_open(sub_filename, sub_file) + call spec_file_read_aero_dist(sub_file, aero_data, aero_dist_init) + call spec_file_close(sub_file) + + call spec_file_read_scenario(file, gas_data, aero_data, scenario) + call spec_file_read_env_state(file, env_state) + + call spec_file_close(file) + + ! All data from spec file read. Do the modal run. + + call pmc_srand(0,0) + + call uuid4_str(run_modal_opt%uuid) + + call scenario_init_env_state(scenario, env_state, 0d0) + + call run_modal(aero_data, aero_dist_init, scenario, & + env_state, gas_data, bin_grid, run_modal_opt) + + call pmc_rand_finalize() + + end subroutine partmc_modal + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + end program partmc From 3d70db27f0c504e9059ad883b10b14e37a9e1678 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 8 Jun 2025 12:27:49 -0500 Subject: [PATCH 03/28] add run_modal.F90 to CMakeLists --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 788097ae3..a0791bc94 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -233,7 +233,7 @@ add_library(partmclib src/aero_state.F90 src/integer_varray.F90 src/aero_particle.F90 src/aero_particle_array.F90 src/mpi.F90 src/netcdf.F90 src/aero_info.F90 src/aero_info_array.F90 src/nucleate.F90 src/condense.F90 src/fractal.F90 src/chamber.F90 - src/camp_interface.F90 src/photolysis.F90 src/sys.F90 + src/camp_interface.F90 src/photolysis.F90 src/sys.F90 src/run_modal.F90 src/tchem_interface.F90 src/aero_component.F90 ${SUNDIALS_SRC} ${GSL_SRC} ${TCHEM_SRC} From 8ac4cd2a281d1f6c10d974af006097c5389c957e Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 8 Jun 2025 13:33:53 -0500 Subject: [PATCH 04/28] add input and output routines for modal run --- src/output.F90 | 145 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) diff --git a/src/output.F90 b/src/output.F90 index b102f0275..31d416401 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -782,6 +782,151 @@ subroutine input_sectional(filename, index, time, del_t, uuid, bin_grid, & end subroutine input_sectional +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Write the current modal data. + subroutine output_modal(prefix, aero_binned, aero_dist, aero_data, & + env_state, gas_data, gas_state, bin_grid, & + scenario, index, time, del_t, uuid) + + !> Prefix of filename to write. + character(len=*), intent(in) :: prefix + !> Binned aerosol distribution. + type(aero_binned_t), intent(in) :: aero_binned + !> Aerosol distribution. + type(aero_dist_t), intent(in) :: aero_dist + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(in) :: env_state + !> Gas data. + type(gas_data_t), intent(in) :: gas_data + !> Gas state. + type(gas_state_t), intent(in) :: gas_state + !> Bin grid. + type(bin_grid_t), intent(in) :: bin_grid + !> Scenario data. + type(scenario_t), intent(in) :: scenario + !> Filename index. + integer, intent(in) :: index + !> Current time (s). + real(kind=dp), intent(in) :: time + !> Current output timestep (s). + real(kind=dp), intent(in) :: del_t + !> UUID of the simulation. + character(len=PMC_UUID_LEN), intent(in) :: uuid + + integer :: ncid + character(len=len(prefix)+100) :: filename + + write(filename, '(a,a,i8.8,a)') trim(prefix), & + '_', index, '.nc' + call pmc_nc_open_write(filename, ncid) + call pmc_nc_write_info(ncid, uuid, & + "PartMC version " // trim(PARTMC_VERSION)) + call write_time(ncid, time, del_t, index) + + ! Write data. + call aero_binned_output_netcdf(aero_binned, ncid, bin_grid, aero_data) + call aero_dist_output_netcdf(aero_dist, ncid) + call env_state_output_netcdf(env_state, ncid) + call gas_data_output_netcdf(gas_data, ncid) + call gas_state_output_netcdf(gas_state, ncid, gas_data) + call aero_data_output_netcdf(aero_data, ncid) + + if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then + call scenario_output_drydep_param(scenario, ncid) + end if + + call pmc_nc_check(nf90_close(ncid)) + + end subroutine output_modal + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Input modal data. + subroutine input_modal(filename, index, time, del_t, uuid, aero_dist, & + aero_binned, aero_data, env_state, gas_data, & + gas_state, bin_grid, scenario) + + !> Filename to read. + character(len=*), intent(in) :: filename + !> Filename index. + integer, intent(out) :: index + !> Current time (s). + real(kind=dp), intent(out) :: time + !> Current output timestep. + real(kind=dp), intent(out) :: del_t + !> UUID of the simulation. + character(len=PMC_UUID_LEN), intent(out) :: uuid + !> Aerosol distribution + type(aero_dist_t), optional, intent(inout) :: aero_dist + !> Binned aerosol distribution. + type(aero_binned_t), optional, intent(inout) :: aero_binned + !> Aerosol data. + type(aero_data_t), optional, intent(inout) :: aero_data + !> Environment state. + type(env_state_t), optional, intent(inout) :: env_state + !> Gas data. + type(gas_data_t), optional, intent(inout) :: gas_data + !> Gas state. + type(gas_state_t), optional, intent(inout) :: gas_state + !> Bin grid. + type(bin_grid_t), optional, intent(inout) :: bin_grid + !> Scenario data. + type(scenario_t), optional, intent(inout) :: scenario + + integer :: ncid + + call assert_msg(348927561, pmc_mpi_rank() == 0, & + "can only call from process 0") + + call pmc_nc_open_read(filename, ncid) + + call pmc_nc_check(nf90_get_att(ncid, NF90_GLOBAL, "UUID", uuid)) + + call pmc_nc_read_real(ncid, time, "time") + call pmc_nc_read_real(ncid, del_t, "timestep") + call pmc_nc_read_integer(ncid, index, "timestep_index") + + if (present(aero_data)) then + call aero_data_input_netcdf(aero_data, ncid) + end if + + if (present(aero_dist)) then + call aero_dist_input_netcdf(aero_dist, ncid) + end if + + if (present(env_state)) then + call env_state_input_netcdf(env_state, ncid) + end if + + if (present(gas_data)) then + call gas_data_input_netcdf(gas_data, ncid) + if (present(gas_state)) then + call gas_state_input_netcdf(gas_state, ncid, gas_data) + end if + else + call assert_msg(739182654, present(gas_state) .eqv. .false., & + "cannot input gas_state without gas_data") + end if + + if (present(bin_grid)) then + call bin_grid_input_netcdf(bin_grid, ncid, "diam") + if (present(aero_binned)) then + call aero_binned_input_netcdf(aero_binned, ncid, bin_grid, aero_data) + end if + else + call assert_msg(582491376, present(aero_binned) .eqv. .false., & + "cannot input aero_binned without bin_grid") + end if + + if (present(scenario)) then + call scenario_input_drydep_param(scenario, ncid) + end if + + end subroutine input_modal + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Read the specification for an output type from a spec file and From 6526f22c1fa180c68728124719527d9b36b69918 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 8 Jun 2025 15:56:51 -0500 Subject: [PATCH 05/28] add input and output routines for modal parameters --- src/aero_dist.F90 | 92 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) diff --git a/src/aero_dist.F90 b/src/aero_dist.F90 index 78a696918..80b733b1f 100644 --- a/src/aero_dist.F90 +++ b/src/aero_dist.F90 @@ -271,6 +271,98 @@ subroutine spec_file_read_aero_dist(file, aero_data, & end subroutine spec_file_read_aero_dist +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Write distribution parameters for each mode in the aerosol distribution. + subroutine aero_dist_output_netcdf(aero_dist, ncid) + + !> Aerosol distribution to write. + type(aero_dist_t), intent(in) :: aero_dist + !> NetCDF file ID, in data mode. + integer, intent(in) :: ncid + + integer :: i_mode + type(aero_mode_t) :: aero_mode + character(len=1) :: num + + do i_mode = 1,aero_dist_n_mode(aero_dist) + aero_mode = aero_dist%mode(i_mode) + write(num, '(I1.1)') i_mode + call pmc_nc_write_real(ncid, aero_mode%char_radius, "char_radius_mode_" & + // trim(num), unit="m") + call pmc_nc_write_real(ncid, aero_mode%log10_std_dev_radius, "std_dev_mode_" & + // trim(num), unit="1", description="log base 10 of geometric standard" & + // " deviation of radius") + call pmc_nc_write_real(ncid, aero_mode%num_conc, "num_conc_mode_" & + // trim(num), unit="m^-3", description="total number concentration for mode") + end do + + call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "num_modes", & + description="total number of modes") + + ! integer :: i_mode, n_mode + ! real(kind=dp), allocatable :: char_rads(:), log10_std_dev_rads(:), num_concs(:) + + ! n_mode = aero_dist_n_mode(aero_dist) + + ! allocate(char_rads(n_mode), log10_std_dev_rads(n_mode), num_concs(n_mode)) + + ! do i_mode = 1,n_mode + ! char_rads(i_mode) = aero_dist%mode(i_mode)%char_radius + ! log10_std_dev_rads(i_mode) = aero_dist%mode(i_mode)%log10_std_dev_radius + ! num_concs(i_mode) = aero_dist%mode(i_mode)%num_conc + ! end do + + ! call pmc_nc_write_real_1d(ncid, char_rads, "char_radius_modes", dim_name="num_modes", & + ! long_name="characteristic radius of each mode", & + ! unit="m") + ! call pmc_nc_write_real_1d(ncid, log10_std_dev_rads, "std_dev_radius_modes", dim_name="num_modes", & + ! long_name="log_10 of geometric std dev for each mode ", & + ! unit="1") + ! call pmc_nc_write_real_1d(ncid, num_concs, "num_conc_modes", dim_name="num_modes", & + ! long_name="number concentration for each mode", & + ! unit="m^-3") + ! call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "num_modes", & + ! dimdescription="total number of modes") + + ! deallocate(char_rads, log10_std_dev_rads, num_concs) + + end subroutine aero_dist_output_netcdf + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Read distribution parameters for each mode in the aerosol distribution. + subroutine aero_dist_input_netcdf(aero_dist, ncid) + + !> Aerosol distribution to read. + type(aero_dist_t), intent(inout) :: aero_dist + !> NetCDF file ID, in data mode + integer, intent(in) :: ncid + + integer :: i_mode, n_mode + character(len=1) :: num + type(aero_mode_t), allocatable :: modes(:) + + call pmc_nc_read_integer(ncid, n_mode, "num_modes") + call assert(836472910, n_mode < 100) + + if (allocated(modes)) deallocate(modes) + allocate(modes(n_mode)) + + do i_mode=1,n_mode + write(num, '(I1.1)') i_mode + call pmc_nc_read_real(ncid, modes(i_mode)%char_radius, "char_radius_mode_" & + // trim(num)) + call pmc_nc_read_real(ncid, modes(i_mode)%log10_std_dev_radius, "std_dev_mode_" & + // trim(num)) + call pmc_nc_read_real(ncid, modes(i_mode)%num_conc, "num_conc_mode_" & + // trim(num)) + end do + + aero_dist%mode = modes + + end subroutine aero_dist_input_netcdf + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Read an array of aero_dists with associated times and rates from From f3dd4538629e4331a708763f77a8ea3951fffc59 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Mon, 9 Jun 2025 15:06:17 -0500 Subject: [PATCH 06/28] fix use statements and typos --- src/output.F90 | 1 + src/partmc.F90 | 5 +++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 31d416401..70d9dfa76 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -76,6 +76,7 @@ module pmc_output use pmc_env_state use pmc_util use pmc_gas_data + use pmc_scenario use pmc_mpi #ifdef PMC_USE_MPI use mpi diff --git a/src/partmc.F90 b/src/partmc.F90 index 24360d6d8..96d3bef89 100644 --- a/src/partmc.F90 +++ b/src/partmc.F90 @@ -226,6 +226,7 @@ program partmc use pmc_run_part use pmc_run_exact use pmc_run_sect + use pmc_run_modal use pmc_spec_file use pmc_gas_data use pmc_gas_state @@ -828,10 +829,10 @@ subroutine partmc_modal(file) call spec_file_read_string(file,'aerosol_init', sub_filename) call spec_file_open(sub_filename, sub_file) - call spec_file_read_aero_dist(sub_file, aero_data, aero_dist_init) + call spec_file_read_aero_dist(sub_file, aero_data, .false., aero_dist_init) call spec_file_close(sub_file) - call spec_file_read_scenario(file, gas_data, aero_data, scenario) + call spec_file_read_scenario(file, gas_data, aero_data, .false., scenario) call spec_file_read_env_state(file, env_state) call spec_file_close(file) From 2872bdc136c83d30e954b8b4772cb95d4d49f7a0 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Mon, 9 Jun 2025 15:06:45 -0500 Subject: [PATCH 07/28] fix typos --- src/run_modal.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/run_modal.F90 b/src/run_modal.F90 index 6804e351b..e75f8b7db 100644 --- a/src/run_modal.F90 +++ b/src/run_modal.F90 @@ -117,7 +117,7 @@ subroutine run_modal(aero_data, aero_dist, scenario, env_state, & env_state, aero_data%density(1), scenario, removed) call check_event(time, run_modal_opt%del_t, run_modal_opt%t_output, & - last_output, do_output) + last_output_time, do_output) if (do_output) then i_summary = i_summary + 1 call aero_binned_add_aero_dist(aero_binned, bin_grid, aero_data, & @@ -129,7 +129,7 @@ subroutine run_modal(aero_data, aero_dist, scenario, env_state, & end if call check_event(time, run_modal_opt%del_t, run_modal_opt%t_progress, & - last_progress, do_progress) + last_progress_time, do_progress) if (do_progress) then write(*, '(a6,a8)') 'step', 'time' write(*, '(i6,f8.1)') i_time, time From fbe3daf6837447f4dbcb9dd9e53b5c587faac6fd Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Mon, 9 Jun 2025 15:08:18 -0500 Subject: [PATCH 08/28] add subroutines for computing integrated deposition velocities for the modal run; handle input and output for processing --- src/scenario.F90 | 272 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) diff --git a/src/scenario.F90 b/src/scenario.F90 index ead729f94..9c5a16253 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -35,6 +35,11 @@ module pmc_scenario !> Type code for a loss rate function for chamber experiments. integer, parameter :: SCENARIO_LOSS_FUNCTION_CHAMBER = 5 + !> Type code for Zhang et al., 2001 dry deposition parameterization. + integer, parameter :: SCENARIO_DRYDEP_ZHANG = 0 + !> Type code for Emeroson et al., 2020 dry deposition parameterization. + integer, parameter :: SCENARIO_DRYDEP_EMERSON = 1 + !> Parameter to switch between algorithms for particle loss. !! A value of 0 will always use the naive algorithm, and !! a value of 1 will always use the accept-reject algorithm. @@ -97,6 +102,8 @@ module pmc_scenario !> Type of loss rate function. integer :: loss_function_type + !> Parameterization for dry deposition. + integer :: drydep_param !> Chamber parameters for wall loss and sedimentation. type(chamber_t) :: chamber end type scenario_t @@ -392,6 +399,81 @@ subroutine scenario_update_aero_binned(scenario, delta_t, env_state, & end subroutine scenario_update_aero_binned +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, density, & + scenario, removed) + + !> Aerosol distribution. + type(aero_dist_t), intent(inout) :: aero_dist + !> Timestep. + real(kind=dp), intent(in) :: del_t + !> Environment state. + type(env_state_t), intent(in) :: env_state + !> Density for each mode. + real(kind=dp), intent(in) :: density + !> Scenario + type(scenario_t), intent(inout) :: scenario + !> Removed (third moment) + real(kind=dp), intent(inout) :: removed + + !> Current mode + type(aero_mode_t) :: aero_mode + + real(kind=dp) :: N, d_pg, ln_sigma_g + real(kind=dp) :: m_0_rate, m_3_rate + real(kind=dp) :: M, new_M + real(kind=dp) :: new_N, new_d_pg + integer :: i_mode + + if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_INVALID) then + return + else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_NONE) then + return + else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_CONSTANT) then + return + else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then + ! loss + do i_mode = 1,aero_dist_n_mode(aero_dist) + aero_mode = aero_dist%mode(i_mode) + N = aero_mode%num_conc + + if (N == 0d0) then + cycle + end if + + d_pg = aero_mode%char_radius * 2.0d0 + ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) + + ! Integrated deposition rate for the 0-th moment (exactly equal to number conc.) + m_0_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, & + 0.0d0, density, env_state) + new_N = N * exp(m_0_rate * del_t) + + if (new_N < 1d0) then + aero_dist%mode(i_mode)%num_conc = 0d0 + aero_dist%mode(i_mode)%char_radius = 0d0 + cycle + end if + + aero_dist%mode(i_mode)%num_conc = new_N + + ! Integrated deposition rate for the 3-rd moment (proportional to volume conc.) + m_3_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, 3.0d0, & + density, env_state) + M = N * d_pg**3.0d0 * exp((3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)) + new_M = M * exp(m_3_rate * del_t) + ! Added variables for debugging + removed = removed + (M - new_M) + + ! New geometric mean diameter + new_d_pg = (new_M / new_N * exp(-(3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)))**(1.0d0/3.0d0) + aero_dist%mode(i_mode)%char_radius = new_d_pg / 2.0d0 + end do + end if + + end subroutine scenario_update_aero_modes + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Evaluate a loss rate function. @@ -537,6 +619,158 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & end function scenario_loss_rate_dry_dep +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Compute and return the integrated dry deposition rate for a given + !> lognormal aerosol mode. + real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & + aero_mode, moment, density, env_state) + + !> Scenario data. + type(scenario_t), intent(in) :: scenario + !> Aerosol mode. + type(aero_mode_t), intent(in) :: aero_mode + !> Moment to calculate loss rate for. + real(kind=dp), intent(in) :: moment + !> Particle density assumed for entire mode (kg m^-3). + real(kind=dp), intent(in) :: density + !> Environment state. + type(env_state_t), intent(in) :: env_state + + real(kind=dp) :: V_d_hat, V_g_hat + real(kind=dp) :: V_g_bar + real(kind=dp) :: d_pg, ln_sigma_g + real(kind=dp) :: density_air + real(kind=dp) :: visc_d, visc_k + real(kind=dp) :: gas_speed, gas_mean_free_path + real(kind=dp) :: knud + real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu + real(kind=dp) :: D_bar, D_hat + real(kind=dp) :: St, Sc + real(kind=dp) :: u_mean, u_star, z_ref, z_rough + real(kind=dp) :: R_a, R_s + real(kind=dp) :: E_B, E_IN, E_IM, R1 + real(kind=dp) :: C_B, C_IN, C_IM + + ! Hardcoded meteorological variables and + ! parameterization-dependent LUC values. + z_ref = 20.0d0 ! Reference height + u_mean = 5.0d0 ! Mean wind speed at reference height + ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 + z_rough = 1.05d0 + A = 5.0d0 / 1000.0d0 + alpha = .8d0 + eps_0 = 3.0d0 + + if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then + gamma = 2.0d0 / 3.0d0 ! Not dependent on LUC + C_B = .2d0 + C_IN = 2.5d0 + C_IM = .4d0 + nu = .8d0 + beta = 1.7d0 + else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then + gamma = .56d0 ! LUC-dependent for Zhang. + C_B = 1.0d0 + C_IN = .5d0 + C_IM = 1.0d0 + nu = 2.0d0 + beta = 2.0d0 + end if + + ! particle diameter equal to geometric mean diameter + d_pg = aero_mode%char_radius * 2.0d0 + ! natural log of geometric standard deviation + ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) + ! density of air + density_air = (const%air_molec_weight * env_state%pressure) & + / (const%univ_gas_const * env_state%temp) + ! dynamic viscosity + visc_d = 1.8325d-5 * (416.16 / (env_state%temp + 120.0d0)) & + * (env_state%temp / 296.16)**1.5d0 + ! kinematic viscosity + visc_k = visc_d / density_air + ! gas speed + gas_speed = & + sqrt((8.0d0 * const%boltzmann * env_state%temp * const%avagadro) / & + (const%pi * const%air_molec_weight)) + ! gas mean free path + gas_mean_free_path = (2.0d0 * visc_d) / (density_air * gas_speed) + ! Knudsen number + knud = (2.0d0 * gas_mean_free_path) / d_pg + ! Settling velcoity + V_g_bar = (density * d_pg**2.0d0 * const%std_grav) / (18.0d0 * visc_d) + ! Compute integrated settling velocity + V_g_hat = V_g_bar * (exp((4.0d0 * moment + 4.0d0) / 2.0d0 * ln_sigma_g**2.0d0) + 1.246d0 * & + knud * exp((2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0)) + + ! Aerodynamic resistance (assuming neutral stability) + u_star = 0.4d0 * u_mean / log(z_ref / z_rough) + R_a = (1.0d0 / (0.4d0 * u_star)) * log(z_ref / z_rough) + + ! Brownian diffusivity + D_bar = (const%boltzmann * env_state%temp) / & + (3.0d0 * const%pi * visc_d * d_pg) + ! Compute integrated Brownian diffusivity + D_hat = D_bar * ((exp((-2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0) + 1.246d0 * & + knud * exp((-4.0d0 * moment + 4.0d0) / 2.0d0 * ln_sigma_g**2.0d0))) + ! Schmidt number based on integrated diffusivity + Sc = visc_k / D_hat + ! Collection efficiency due to Brownian diffusion + E_B = C_B * Sc**(-gamma) + + ! Collection efficiency due to interception + E_IN = C_IN * (d_pg / A)**nu + + ! Stokes number based on integrated settling velocity + St = (V_g_hat * u_star) / (const%std_grav * A) + ! Collection efficiency due to impaction + E_IM = C_IM * (St / (alpha + St))**beta + + ! Rebound correction + R1 = exp(-St**0.5d0) + + ! Surface resistance + R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) + + ! Integrated deposition velocity + V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s + R_a * R_s * V_g_hat)) + + ! Loss rate + scenario_integrated_loss_rate_dry_dep = V_d_hat / env_state%height + + end function scenario_integrated_loss_rate_dry_dep + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Updates a given array to contain the integrated deposition rates for the + !> given moment of each aerosol mode in the distribution. + subroutine scenario_modal_dry_dep_rates(scenario, aero_dist, moment, & + density, env_state, rates) + + !> Scenario data. + type(scenario_t), intent(in) :: scenario + !> Aerosol distribution. + type(aero_dist_t), intent(in) :: aero_dist + !> Moment to compute rate for. + real(kind=dp), intent(in) :: moment + !> Density for the mode. + real(kind=dp), intent(in) :: density + !> Environment state. + type(env_state_t), intent(in) :: env_state + !> Rates. + real(kind=dp), intent(inout) :: rates(:) + + integer :: i_mode, n_mode + + do i_mode = 1,size(rates) + rates(i_mode) = scenario_integrated_loss_rate_dry_dep( & + scenario, aero_dist%mode(i_mode), moment, density, env_state) & + * env_state%height + end do + + end subroutine + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute and return the max loss rate function for a given volume. @@ -889,6 +1123,15 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_VOLUME else if (trim(function_name) == 'drydep') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_DRYDEP + call spec_file_read_string(file, "dry_dep param", function_name) + if (trim(function_name) == 'Zhang_2001') then + scenario%drydep_param = SCENARIO_DRYDEP_ZHANG + else if (trim(function_name) == 'Emerson_2020') then + scenario%drydep_param = SCENARIO_DRYDEP_EMERSON + else + call die_msg(395827406, "Unknown dry deposition parameterization: " & + // trim(function_name)) + end if else if (trim(function_name) == 'chamber') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_CHAMBER call spec_file_read_chamber(file, scenario%chamber) @@ -1000,6 +1243,35 @@ end subroutine spec_file_read_scenario !! - \ref input_format_scenario --- the environment data !! containing the mixing layer height profile +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Output the dry deposition parameterization to NetCDF file. + subroutine scenario_output_drydep_param(scenario, ncid) + + !> Scenario data. + type(scenario_t), intent(in) :: scenario + !> NetCDF file ID, in data mode. + integer, intent(in) :: ncid + + call pmc_nc_write_integer(ncid, scenario%drydep_param, & + "drydep_param") + + end subroutine scenario_output_drydep_param + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Input the dry deposition parameterization to NetCDF file. + subroutine scenario_input_drydep_param(scenario, ncid) + + !> Scenario data. + type(scenario_t), intent(inout) :: scenario + !> NetCDF file ID, in data mode. + integer, intent(in) :: ncid + + call pmc_nc_read_integer(ncid, scenario%drydep_param, "drydep_param") + + end subroutine scenario_input_drydep_param + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Determines the number of bytes required to pack the given value. From b27627feca2b2c89b6055f1f4d01050b5089ede1 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 09:39:20 -0500 Subject: [PATCH 09/28] add support for QUADPACK; comment out dry dep scenario for now --- CMakeLists.txt | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a0791bc94..ccc6a2315 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,7 @@ option(ENABLE_TCHEM "Enable TChem chemistry support" OFF) option(ENABLE_MPI "Enable MPI parallel support" OFF) option(ENABLE_SUNDIALS "Enable SUNDIALS solver for condensation support" OFF) option(ENABLE_C_SORT "Enable C sorting routines" OFF) +option(ENABLE_QUADPACK "Enable QUADPACK for numerical integration" OFF) ###################################################################### # CPack @@ -179,6 +180,25 @@ if(ENABLE_SUNDIALS) add_definitions(-DPMC_USE_SUNDIALS) endif() +###################################################################### +# QUADPACK + +if (ENABLE_QUADPACK) + find_path(QUADPACK_INCLUDE_DIR quadpack_double.mod + DOC "QUADPACK include directory" + PATHS $ENV{QUADPACK_HOME}/include) + find_library(QUADPACK_LIB quadpack + DOC "QUADPACK library" + PATHS $ENV{QUADPACK_HOME}/lib) + if(EXISTS "${QUADPACK_INCLUDE_DIR}/quadpack.mod" AND EXISTS "${QUADPACK_LIB}") + include_directories(${QUADPACK_INCLUDE_DIR}) + add_definitions(-DPMC_USE_QUADPACK) + else() + message(WARNING "QUADPACK library not found - disabling QUADPACK support") + set(ENABLE_QUADPACK OFF) + endif() +endif() + ###################################################################### # tests @@ -241,7 +261,7 @@ add_library(partmclib src/aero_state.F90 src/integer_varray.F90 target_link_libraries(partmclib ${NETCDF_LIBS} ${SUNDIALS_LIBS} ${MOSAIC_LIB} ${GSL_LIBS} ${CAMP_LIB} ${TCHEM_LIB} ${YAML_LIB} - ${KOKKOS_LIB} ${KOKKOSKERNEL_LIB} ${TINES_LIB} ${CPP_LIB} ${LAPACK_LIB}) + ${KOKKOS_LIB} ${KOKKOSKERNEL_LIB} ${TINES_LIB} ${CPP_LIB} ${LAPACK_LIB} ${QUADPACK_LIB}) if (ENABLE_TCHEM) find_package(OpenMP REQUIRED) @@ -413,3 +433,10 @@ add_executable(chamber_process target_link_libraries(chamber_process partmclib) ###################################################################### +# scenarios/7_dry_dep/dry_dep_modal_process + +add_executable(dry_dep_modal_process + scenarios/7_dry_dep/dry_dep_modal_process.F90) +target_link_libraries(dry_dep_modal_process partmclib) + +###################################################################### \ No newline at end of file From e530f6fe283824294ad8554a363f2774e4f7cf02 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 09:40:30 -0500 Subject: [PATCH 10/28] remove debugging code for modal run --- src/output.F90 | 1 + src/run_modal.F90 | 5 ++--- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 70d9dfa76..3101c9ee4 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -834,6 +834,7 @@ subroutine output_modal(prefix, aero_binned, aero_dist, aero_data, & call gas_data_output_netcdf(gas_data, ncid) call gas_state_output_netcdf(gas_state, ncid, gas_data) call aero_data_output_netcdf(aero_data, ncid) + call bin_grid_output_netcdf(bin_grid, ncid, "diam", unit="m") if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then call scenario_output_drydep_param(scenario, ncid) diff --git a/src/run_modal.F90 b/src/run_modal.F90 index e75f8b7db..bd4f1b9d3 100644 --- a/src/run_modal.F90 +++ b/src/run_modal.F90 @@ -62,7 +62,7 @@ subroutine run_modal(aero_data, aero_dist, scenario, env_state, & type(gas_state_t) :: gas_state type(aero_binned_t) :: aero_binned - real(kind=dp) time, last_output_time, last_progress_time, removed + real(kind=dp) time, last_output_time, last_progress_time integer i, i_time, n_time, i_summary, i_mode, mode logical do_output, do_progress @@ -87,7 +87,6 @@ subroutine run_modal(aero_data, aero_dist, scenario, env_state, & last_output_time = 0d0 last_progress_time = 0d0 i_summary = 1 - removed = 0d0 ! initial output call check_event(time, run_modal_opt%del_t, run_modal_opt%t_output, & @@ -114,7 +113,7 @@ subroutine run_modal(aero_data, aero_dist, scenario, env_state, & call scenario_update_gas_state(scenario, run_modal_opt%del_t, & env_state, old_env_state, gas_data, gas_state) call scenario_update_aero_modes(aero_dist, run_modal_opt%del_t, & - env_state, aero_data%density(1), scenario, removed) + env_state, aero_data%density(1), scenario) call check_event(time, run_modal_opt%del_t, run_modal_opt%t_output, & last_output_time, do_output) From 562c066679fd6d4cda68c64fe7d207e944abc78f Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 09:42:18 -0500 Subject: [PATCH 11/28] clean up implementation of numerical integration with QUADPACK --- src/scenario.F90 | 235 ++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 201 insertions(+), 34 deletions(-) diff --git a/src/scenario.F90 b/src/scenario.F90 index 9c5a16253..e63d03255 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -21,6 +21,9 @@ module pmc_scenario #ifdef PMC_USE_MPI use mpi #endif +#ifdef PMC_USE_QUADPACK + use quadpack_double, only: dqagse +#endif !> Type code for an undefined or invalid loss function. integer, parameter :: SCENARIO_LOSS_FUNCTION_INVALID = 0 @@ -401,8 +404,8 @@ end subroutine scenario_update_aero_binned !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, density, & - scenario, removed) + subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, & + density, scenario) !> Aerosol distribution. type(aero_dist_t), intent(inout) :: aero_dist @@ -414,8 +417,6 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, density, & real(kind=dp), intent(in) :: density !> Scenario type(scenario_t), intent(inout) :: scenario - !> Removed (third moment) - real(kind=dp), intent(inout) :: removed !> Current mode type(aero_mode_t) :: aero_mode @@ -445,9 +446,9 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, density, & d_pg = aero_mode%char_radius * 2.0d0 ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) - ! Integrated deposition rate for the 0-th moment (exactly equal to number conc.) + ! Integrated deposition rate for the 0-th moment (num. conc.) m_0_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, & - 0.0d0, density, env_state) + 0.0d0, density, env_state) new_N = N * exp(m_0_rate * del_t) if (new_N < 1d0) then @@ -463,8 +464,6 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, density, & density, env_state) M = N * d_pg**3.0d0 * exp((3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)) new_M = M * exp(m_3_rate * del_t) - ! Added variables for debugging - removed = removed + (M - new_M) ! New geometric mean diameter new_d_pg = (new_M / new_N * exp(-(3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)))**(1.0d0/3.0d0) @@ -622,7 +621,7 @@ end function scenario_loss_rate_dry_dep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute and return the integrated dry deposition rate for a given - !> lognormal aerosol mode. + !> lognormal aerosol mode (modal approximation). real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & aero_mode, moment, density, env_state) @@ -652,30 +651,37 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & real(kind=dp) :: E_B, E_IN, E_IM, R1 real(kind=dp) :: C_B, C_IN, C_IM +#ifdef PMC_USE_QUADPACK + ! Do numerical integration when QUADPACK is available + scenario_integrated_loss_rate_dry_dep = scenario_integrated_loss_rate_dry_dep_quadpack( & + scenario, aero_mode, moment, density, env_state) + return +#endif + ! Hardcoded meteorological variables and ! parameterization-dependent LUC values. z_ref = 20.0d0 ! Reference height u_mean = 5.0d0 ! Mean wind speed at reference height ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 - z_rough = 1.05d0 - A = 5.0d0 / 1000.0d0 - alpha = .8d0 + z_rough = .8d0 + A = 2.0d0 / 1000.0d0 + alpha = 1.0d0 eps_0 = 3.0d0 if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then - gamma = 2.0d0 / 3.0d0 ! Not dependent on LUC - C_B = .2d0 - C_IN = 2.5d0 - C_IM = .4d0 - nu = .8d0 - beta = 1.7d0 - else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then gamma = .56d0 ! LUC-dependent for Zhang. - C_B = 1.0d0 - C_IN = .5d0 - C_IM = 1.0d0 - nu = 2.0d0 - beta = 2.0d0 + C_B = .2d0 + C_IN = 2.5d0 + C_IM = .4d0 + nu = .8d0 + beta = 1.7d0 + else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then + gamma = .56d0 ! LUC-dependent for Zhang. + C_B = 1.0d0 + C_IN = .5d0 + C_IM = 1.0d0 + nu = 2.0d0 + beta = 2.0d0 end if ! particle diameter equal to geometric mean diameter @@ -684,16 +690,15 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) ! density of air density_air = (const%air_molec_weight * env_state%pressure) & - / (const%univ_gas_const * env_state%temp) + / (const%univ_gas_const * env_state%temp) ! dynamic viscosity visc_d = 1.8325d-5 * (416.16 / (env_state%temp + 120.0d0)) & - * (env_state%temp / 296.16)**1.5d0 + * (env_state%temp / 296.16)**1.5d0 ! kinematic viscosity visc_k = visc_d / density_air ! gas speed - gas_speed = & - sqrt((8.0d0 * const%boltzmann * env_state%temp * const%avagadro) / & - (const%pi * const%air_molec_weight)) + gas_speed = sqrt((8.0d0 * const%boltzmann * env_state%temp * const%avagadro) / & + (const%pi * const%air_molec_weight)) ! gas mean free path gas_mean_free_path = (2.0d0 * visc_d) / (density_air * gas_speed) ! Knudsen number @@ -702,7 +707,7 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & V_g_bar = (density * d_pg**2.0d0 * const%std_grav) / (18.0d0 * visc_d) ! Compute integrated settling velocity V_g_hat = V_g_bar * (exp((4.0d0 * moment + 4.0d0) / 2.0d0 * ln_sigma_g**2.0d0) + 1.246d0 * & - knud * exp((2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0)) + knud * exp((2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0)) ! Aerodynamic resistance (assuming neutral stability) u_star = 0.4d0 * u_mean / log(z_ref / z_rough) @@ -710,10 +715,10 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & ! Brownian diffusivity D_bar = (const%boltzmann * env_state%temp) / & - (3.0d0 * const%pi * visc_d * d_pg) + (3.0d0 * const%pi * visc_d * d_pg) ! Compute integrated Brownian diffusivity D_hat = D_bar * ((exp((-2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0) + 1.246d0 * & - knud * exp((-4.0d0 * moment + 4.0d0) / 2.0d0 * ln_sigma_g**2.0d0))) + knud * exp((-4.0d0 * moment + 4.0d0) / 2.0d0 * ln_sigma_g**2.0d0))) ! Schmidt number based on integrated diffusivity Sc = visc_k / D_hat ! Collection efficiency due to Brownian diffusion @@ -734,13 +739,175 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) ! Integrated deposition velocity - V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s + R_a * R_s * V_g_hat)) + V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s)) ! Loss rate scenario_integrated_loss_rate_dry_dep = V_d_hat / env_state%height end function scenario_integrated_loss_rate_dry_dep +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#ifdef PMC_USE_QUADPACK + real(kind=dp) function scenario_integrated_loss_rate_dry_dep_quadpack(scenario, & + aero_mode, moment, density, env_state) + + !> Scenario data. + type(scenario_t), intent(in) :: scenario + !> Aerosol mode. + type(aero_mode_t), intent(in) :: aero_mode + !> Moment to calculate loss rate for. + real(kind=dp), intent(in) :: moment + !> Particle density assumed for entire mode (kg m^-3). + real(kind=dp), intent(in) :: density + !> Environment state. + type(env_state_t), intent(in) :: env_state + + real(kind=dp) :: d_pg, ln_sigma_g + real(kind=dp) :: density_air + real(kind=dp) :: visc_d, visc_k + real(kind=dp) :: gas_speed, gas_mean_free_path + real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu + real(kind=dp) :: u_mean, u_star, z_ref, z_rough + real(kind=dp) :: C_B, C_IN, C_IM, M_k + real(kind=dp) :: lower, upper, epsabs, epsrel, result, abserr, ref + integer :: limit, neval, ier, last + real(kind=dp), allocatable :: alist(:), blist(:), rlist(:), elist(:) + integer, allocatable :: iord(:) + + ! Hardcoded meteorological variables and + ! parameterization-dependent LUC values. + z_ref = 20.0d0 ! Reference height + u_mean = 5.0d0 ! Mean wind speed at reference height + ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 + z_rough = .8d0 + A = 2.0d0 / 1000.0d0 + alpha = 1.0d0 + eps_0 = 3.0d0 + + if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then + gamma = .56d0 ! LUC-dependent for Zhang. + C_B = .2d0 + C_IN = 2.5d0 + C_IM = .4d0 + nu = .8d0 + beta = 1.7d0 + else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then + gamma = .56d0 ! LUC-dependent for Zhang. + C_B = 1.0d0 + C_IN = .5d0 + C_IM = 1.0d0 + nu = 2.0d0 + beta = 2.0d0 + end if + + ! particle diameter equal to geometric mean diameter + d_pg = aero_mode%char_radius * 2.0d0 + ! natural log of geometric standard deviation + ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) + ! density of air + density_air = (const%air_molec_weight * env_state%pressure) & + / (const%univ_gas_const * env_state%temp) + ! dynamic viscosity + visc_d = 1.8325d-5 * (416.16 / (env_state%temp + 120.0d0)) & + * (env_state%temp / 296.16)**1.5d0 + ! kinematic viscosity + visc_k = visc_d / density_air + ! gas speed + gas_speed = sqrt((8.0d0 * const%boltzmann * env_state%temp * const%avagadro) / & + (const%pi * const%air_molec_weight)) + ! gas mean free path + gas_mean_free_path = (2.0d0 * visc_d) / (density_air * gas_speed) + + ref = exp(log(d_pg) + moment*ln_sigma_g**2) + lower = ref / (exp(ln_sigma_g)*10.0) + upper = ref * exp(ln_sigma_g)*10.0 + + ! Set integration parameters + limit = 1000 + epsabs = 1.0d-10 + epsrel = 1.0d-6 + + ! Allocate arrays for QUADPACK integration + allocate(alist(limit), blist(limit), rlist(limit), elist(limit), iord(limit)) + + ! Call QUADPACK integration routine + call dqagse(dep_vel_integrand, lower, upper, epsabs, epsrel, limit, & + result, abserr, neval, ier, alist, blist, rlist, elist, iord, last) + + if (ier > 0) then + call die_msg(837465, "QUADPACK integration failed (error code: " & + // trim(integer_to_string(ier))) + endif + + M_k = d_pg**moment * exp(moment**2 * ln_sigma_g**2 / 2.0d0) + + ! Integration result + scenario_integrated_loss_rate_dry_dep_quadpack = 1.0d0 / M_k * result / env_state%height + + ! Clean up + deallocate(alist, blist, rlist, elist, iord) + + contains + + real(kind=dp) function dep_vel_integrand(d_p) + real(kind=dp), intent(in) :: d_p + + ! Local variables + real(kind=dp) :: V_d, V_s + real(kind=dp) :: knud_local, cunning + real(kind=dp) :: diff_p, Sc, St + real(kind=dp) :: u_star, R_a, R_s + real(kind=dp) :: E_B, E_IN, E_IM, R1 + real(kind=dp) :: ln_dp, ln_dp_g, n_ddp + + ! Knudsen number for this diameter + knud_local = (2.0d0 * gas_mean_free_path) / d_p + + ! Cunningham correction factor + cunning = 1.0d0 + knud_local * (1.257d0 + 0.4d0 * exp(-1.1d0 / knud_local)) + + ! Settling velocity + V_s = (density * d_p**2.0d0 * const%std_grav * cunning) / (18.0d0 * visc_d) + + ! Friction velocity and aerodynamic resistance + u_star = 0.4d0 * u_mean / log(z_ref / z_rough) + R_a = (1.0d0 / (0.4d0 * u_star)) * log(z_ref / z_rough) + + ! Particle diffusivity and Schmidt number + diff_p = (const%boltzmann * env_state%temp * cunning) / & + (3.0d0 * const%pi * visc_d * d_p) + Sc = visc_k / diff_p + + ! Collection efficiencies + E_B = C_B * Sc**(-gamma) + E_IN = C_IN * (d_p / A)**nu + + ! Stokes number and impaction efficiency + St = (V_s * u_star) / (const%std_grav * A) + E_IM = C_IM * (St / (alpha + St))**beta + + ! Rebound correction + R1 = exp(-St**0.5d0) + + ! Surface resistance + R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) + + ! Deposition velocity + V_d = V_s + (1.0d0 / (R_a + R_s)) + + ! Log-normal size distribution + ln_dp = log(d_p) + ln_dp_g = log(d_pg) + n_ddp = (1.0d0/(sqrt(2.0d0 * const%pi) * d_p * ln_sigma_g)) * & + exp(-((ln_dp - ln_dp_g)**2) / (2.0d0 * ln_sigma_g**2)) + + ! Final integrand + dep_vel_integrand = d_p**moment * V_d * n_ddp + end function dep_vel_integrand + end function scenario_integrated_loss_rate_dry_dep_quadpack +#endif + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Updates a given array to contain the integrated deposition rates for the @@ -1123,7 +1290,7 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_VOLUME else if (trim(function_name) == 'drydep') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_DRYDEP - call spec_file_read_string(file, "dry_dep param", function_name) + call spec_file_read_string(file, "dry_dep_param", function_name) if (trim(function_name) == 'Zhang_2001') then scenario%drydep_param = SCENARIO_DRYDEP_ZHANG else if (trim(function_name) == 'Emerson_2020') then From 374d5714eb7514831566813ccd9434fb5fa11397 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 10:25:20 -0500 Subject: [PATCH 12/28] formatting --- CMakeLists.txt | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ccc6a2315..e62d1a3f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,7 +31,7 @@ option(ENABLE_TCHEM "Enable TChem chemistry support" OFF) option(ENABLE_MPI "Enable MPI parallel support" OFF) option(ENABLE_SUNDIALS "Enable SUNDIALS solver for condensation support" OFF) option(ENABLE_C_SORT "Enable C sorting routines" OFF) -option(ENABLE_QUADPACK "Enable QUADPACK for numerical integration" OFF) +option(ENABLE_QUADPACK "Enable QUADPACK for numerical integration (dry dep)" OFF) ###################################################################### # CPack @@ -190,13 +190,13 @@ if (ENABLE_QUADPACK) find_library(QUADPACK_LIB quadpack DOC "QUADPACK library" PATHS $ENV{QUADPACK_HOME}/lib) - if(EXISTS "${QUADPACK_INCLUDE_DIR}/quadpack.mod" AND EXISTS "${QUADPACK_LIB}") - include_directories(${QUADPACK_INCLUDE_DIR}) - add_definitions(-DPMC_USE_QUADPACK) - else() - message(WARNING "QUADPACK library not found - disabling QUADPACK support") - set(ENABLE_QUADPACK OFF) - endif() + if(EXISTS "${QUADPACK_INCLUDE_DIR}/quadpack.mod" AND EXISTS "${QUADPACK_LIB}") + include_directories(${QUADPACK_INCLUDE_DIR}) + add_definitions(-DPMC_USE_QUADPACK) + else() + message(WARNING "QUADPACK library not found - disabling QUADPACK support") + set(ENABLE_QUADPACK OFF) + endif() endif() ###################################################################### @@ -435,8 +435,8 @@ target_link_libraries(chamber_process partmclib) ###################################################################### # scenarios/7_dry_dep/dry_dep_modal_process -add_executable(dry_dep_modal_process - scenarios/7_dry_dep/dry_dep_modal_process.F90) -target_link_libraries(dry_dep_modal_process partmclib) +# add_executable(dry_dep_modal_process +# scenarios/7_dry_dep/dry_dep_modal_process.F90) +# target_link_libraries(dry_dep_modal_process partmclib) ###################################################################### \ No newline at end of file From 537b2c32a1f1b2842d52ab08b81576e46faf85e6 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 10:32:20 -0500 Subject: [PATCH 13/28] remove unused flag for emissions/dilution --- src/run_modal.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/run_modal.F90 b/src/run_modal.F90 index bd4f1b9d3..ee22ebd16 100644 --- a/src/run_modal.F90 +++ b/src/run_modal.F90 @@ -27,8 +27,6 @@ module pmc_run_modal real(kind=dp) :: t_output !> Progress interval (0 disables) (s). real(kind=dp) :: t_progress - !> Whether to do emissions and dilution. - logical :: do_emissions !> Output prefix. character(len=300) :: prefix !> UUID of the simulation. From 6da5eeb8ff4ecad488195143d9e32b7764e77ea2 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 13 Jul 2025 10:39:36 -0500 Subject: [PATCH 14/28] fix conditional for inputting dry deposition parameterization --- src/output.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/output.F90 b/src/output.F90 index 3101c9ee4..7f888e1a5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -923,7 +923,7 @@ subroutine input_modal(filename, index, time, del_t, uuid, aero_dist, & "cannot input aero_binned without bin_grid") end if - if (present(scenario)) then + if (present(scenario) .and. scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then call scenario_input_drydep_param(scenario, ncid) end if From bdc2de6166a7e02aef3624505233423cd99e7c4f Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 7 Sep 2025 13:57:11 -0500 Subject: [PATCH 15/28] allow for reading drydep_params from text file otherwise fall back on default parameters --- src/output.F90 | 8 -- src/scenario.F90 | 342 ++++++++++++++++++++++++++++------------------ src/spec_file.F90 | 19 +++ 3 files changed, 230 insertions(+), 139 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 7f888e1a5..4cba0a9fe 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -836,10 +836,6 @@ subroutine output_modal(prefix, aero_binned, aero_dist, aero_data, & call aero_data_output_netcdf(aero_data, ncid) call bin_grid_output_netcdf(bin_grid, ncid, "diam", unit="m") - if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then - call scenario_output_drydep_param(scenario, ncid) - end if - call pmc_nc_check(nf90_close(ncid)) end subroutine output_modal @@ -923,10 +919,6 @@ subroutine input_modal(filename, index, time, del_t, uuid, aero_dist, & "cannot input aero_binned without bin_grid") end if - if (present(scenario) .and. scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then - call scenario_input_drydep_param(scenario, ncid) - end if - end subroutine input_modal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/scenario.F90 b/src/scenario.F90 index e63d03255..6f130ac8a 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -38,16 +38,39 @@ module pmc_scenario !> Type code for a loss rate function for chamber experiments. integer, parameter :: SCENARIO_LOSS_FUNCTION_CHAMBER = 5 - !> Type code for Zhang et al., 2001 dry deposition parameterization. - integer, parameter :: SCENARIO_DRYDEP_ZHANG = 0 - !> Type code for Emeroson et al., 2020 dry deposition parameterization. - integer, parameter :: SCENARIO_DRYDEP_EMERSON = 1 - !> Parameter to switch between algorithms for particle loss. !! A value of 0 will always use the naive algorithm, and !! a value of 1 will always use the accept-reject algorithm. real(kind=dp), parameter :: SCENARIO_LOSS_ALG_THRESHOLD = 1.0d0 + !> Parameters for simulating dry deposition + type drydep_params_t + !> Reference height + real(kind=dp) :: z_ref = 20.0d0 + !> Mean wind speed at reference height + real(kind=dp) :: u_mean = 5.0d0 + !> Roughness length (land-use category dependent) + real(kind=dp) :: z_rough = 0.8d0 ! crops, mixed farming (LUC 7 from Zhang et al., 2001) + !> Characteristic radius of collectors (land-use category dependent) + real(kind=dp) :: A = 2.0d0 / 1000.0d0 + !> Impaction parameter + real(kind=dp) :: alpha = 1.0d0 + !> Surface resistance parameter + real(kind=dp) :: eps_0 = 3.0d0 + !> Diffusivity paramter + real(kind=dp) :: gamma = 0.56d0 ! LUC-dependent for Zhang et al., 2001 + !> Brownian diffusion coefficient + real(kind=dp) :: C_B = 1.0d0 + !> Interception coefficient + real(kind=dp) :: C_IN = 0.5d0 + !> Impaction coefficient + real(kind=dp) :: C_IM = 1.0d0 + !> Interception exponent + real(kind=dp) :: nu = 2.0d0 + !> Impaction exponent + real(kind=dp) :: beta = 2.0d0 + end type drydep_params_t + !> Scenario data. !! !! This is everything needed to drive the scenario being simulated. @@ -105,8 +128,8 @@ module pmc_scenario !> Type of loss rate function. integer :: loss_function_type - !> Parameterization for dry deposition. - integer :: drydep_param + !> Dry deposition parameters + type(drydep_params_t) :: drydep !> Chamber parameters for wall loss and sedimentation. type(chamber_t) :: chamber end type scenario_t @@ -503,7 +526,7 @@ real(kind=dp) function scenario_loss_rate(scenario, vol, density, & else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) & then scenario_loss_rate = scenario_loss_rate_dry_dep(vol, density, & - aero_data, env_state) + aero_data, env_state, scenario) else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_CHAMBER) & then scenario_loss_rate = chamber_loss_rate_wall(scenario%chamber, vol, & @@ -523,7 +546,7 @@ end function scenario_loss_rate !! All equations used here are written in detail in the file !! \c doc/deposition/deposition.tex. real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & - env_state) + env_state, scenario) !> Particle volume (m^3). real(kind=dp), intent(in) :: vol @@ -533,6 +556,8 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & type(aero_data_t), intent(in) :: aero_data !> Environment state. type(env_state_t), intent(in) :: env_state + !> Scenario data. + type(scenario_t), intent(in) :: scenario real(kind=dp) :: V_d real(kind=dp) :: V_s @@ -543,22 +568,25 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & real(kind=dp) :: knud, cunning real(kind=dp) :: grav real(kind=dp) :: R_s, R_a - real(kind=dp) :: alpha, beta, gamma, A, eps_0 + ! real(kind=dp) :: alpha, beta, gamma, A, eps_0 real(kind=dp) :: diff_p real(kind=dp) :: von_karman real(kind=dp) :: St, Sc, u_star real(kind=dp) :: E_B, E_IM, E_IN, R1 - real(kind=dp) :: u_mean, z_ref, z_rough - - ! User set variables - u_mean = 5.0d0 ! Mean wind speed at reference height - z_ref = 20.0d0 ! Reference height - ! Setting for LUC = 7, SC = 1 - See Table 3 - z_rough = .1d0 ! According to land category - A = 2.0d0 / 1000.0d0 ! Dependent on land type - alpha = 1.2d0 ! From table - beta = 2.0d0 ! From text - gamma = .54d0 ! From table + ! real(kind=dp) :: u_mean, z_ref, z_rough + type(drydep_params_t) :: drydep_params + + drydep_params = scenario%drydep + + ! ! User set variables + ! u_mean = 5.0d0 ! Mean wind speed at reference height + ! z_ref = 20.0d0 ! Reference height + ! ! Setting for LUC = 7, SC = 1 - See Table 3 + ! z_rough = .1d0 ! According to land category + ! A = 2.0d0 / 1000.0d0 ! Dependent on land type + ! alpha = 1.2d0 ! From table + ! beta = 2.0d0 ! From text + ! gamma = .54d0 ! From table ! particle diameter d_p = aero_data_vol2diam(aero_data, vol) @@ -587,28 +615,27 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & ! Aerodynamic resistance ! For neutral stability - u_star = .4d0 * u_mean / log(z_ref / z_rough) - R_a = (1.0d0 / (.4d0 * u_star)) * log(z_ref / z_rough) + u_star = .4d0 * drydep_params%u_mean / log(drydep_params%z_ref / drydep_params%z_rough) + R_a = (1.0d0 / (.4d0 * u_star)) * log(drydep_params%z_ref / drydep_params%z_rough) ! Brownian diffusion efficiency diff_p = (const%boltzmann * env_state%temp * cunning) / & (3.d0 * const%pi * visc_d * d_p) Sc = visc_k / diff_p - E_B = Sc**(-gamma) + E_B = drydep_params%C_B * Sc**(-drydep_params%gamma) ! Interception efficiency ! Characteristic radius of large collectors - E_IN = .5d0 * (d_p / A)**2.0d0 + E_IN = drydep_params%C_IN * (d_p / drydep_params%A)**2.0d0 ! Impaction efficiency - St = (V_s * u_star) / (grav * A) - E_IM = (St / (alpha + St))**beta + St = (V_s * u_star) / (grav * drydep_params%A) + E_IM = drydep_params%C_IM * (St / (drydep_params%alpha + St))**drydep_params%beta ! Rebound correction R1 = exp(-St**.5d0) ! Surface resistance - eps_0 = 3.0d0 ! Taken to be 3 - R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) + R_s = 1.0d0 / (drydep_params%eps_0 * u_star * (E_B + E_IN + E_IM) * R1) ! Dry deposition V_d = V_s + (1.0d0 / (R_a + R_s + R_a * R_s * V_s)) @@ -643,13 +670,14 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & real(kind=dp) :: visc_d, visc_k real(kind=dp) :: gas_speed, gas_mean_free_path real(kind=dp) :: knud - real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu + ! real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu real(kind=dp) :: D_bar, D_hat real(kind=dp) :: St, Sc - real(kind=dp) :: u_mean, u_star, z_ref, z_rough + real(kind=dp) :: u_star real(kind=dp) :: R_a, R_s real(kind=dp) :: E_B, E_IN, E_IM, R1 - real(kind=dp) :: C_B, C_IN, C_IM + ! real(kind=dp) :: C_B, C_IN, C_IM + type(drydep_params_t) :: drydep_params #ifdef PMC_USE_QUADPACK ! Do numerical integration when QUADPACK is available @@ -658,31 +686,37 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & return #endif - ! Hardcoded meteorological variables and - ! parameterization-dependent LUC values. - z_ref = 20.0d0 ! Reference height - u_mean = 5.0d0 ! Mean wind speed at reference height - ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 - z_rough = .8d0 - A = 2.0d0 / 1000.0d0 - alpha = 1.0d0 - eps_0 = 3.0d0 - - if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then - gamma = .56d0 ! LUC-dependent for Zhang. - C_B = .2d0 - C_IN = 2.5d0 - C_IM = .4d0 - nu = .8d0 - beta = 1.7d0 - else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then - gamma = .56d0 ! LUC-dependent for Zhang. - C_B = 1.0d0 - C_IN = .5d0 - C_IM = 1.0d0 - nu = 2.0d0 - beta = 2.0d0 - end if + ! ! Hardcoded meteorological variables and + ! ! parameterization-dependent LUC values. + ! z_ref = 20.0d0 ! Reference height + ! u_mean = 5.0d0 ! Mean wind speed at reference height + ! ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 + ! z_rough = .8d0 + ! A = 2.0d0 / 1000.0d0 + ! alpha = 1.0d0 + ! eps_0 = 3.0d0 + + ! if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then + ! gamma = .56d0 ! LUC-dependent for Zhang. + ! C_B = .2d0 + ! C_IN = 2.5d0 + ! C_IM = .4d0 + ! nu = .8d0 + ! beta = 1.7d0 + ! else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then + ! gamma = .56d0 ! LUC-dependent for Zhang. + ! C_B = 1.0d0 + ! C_IN = .5d0 + ! C_IM = 1.0d0 + ! nu = 2.0d0 + ! beta = 2.0d0 + ! end if + drydep_params = scenario%drydep + + print*, drydep_params%z_ref, drydep_params%u_mean, drydep_params%z_rough, & + drydep_params%A, drydep_params%alpha, drydep_params%eps_0, & + drydep_params%gamma, drydep_params%C_B, drydep_params%C_IN, & + drydep_params%C_IM, drydep_params%nu, drydep_params%beta ! particle diameter equal to geometric mean diameter d_pg = aero_mode%char_radius * 2.0d0 @@ -710,8 +744,8 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & knud * exp((2.0d0 * moment + 1.0d0) / 2.0d0 * ln_sigma_g**2.0d0)) ! Aerodynamic resistance (assuming neutral stability) - u_star = 0.4d0 * u_mean / log(z_ref / z_rough) - R_a = (1.0d0 / (0.4d0 * u_star)) * log(z_ref / z_rough) + u_star = 0.4d0 * drydep_params%u_mean / log(drydep_params%z_ref / drydep_params%z_rough) + R_a = (1.0d0 / (0.4d0 * u_star)) * log(drydep_params%z_ref / drydep_params%z_rough) ! Brownian diffusivity D_bar = (const%boltzmann * env_state%temp) / & @@ -722,21 +756,21 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & ! Schmidt number based on integrated diffusivity Sc = visc_k / D_hat ! Collection efficiency due to Brownian diffusion - E_B = C_B * Sc**(-gamma) + E_B = drydep_params%C_B * Sc**(-drydep_params%gamma) ! Collection efficiency due to interception - E_IN = C_IN * (d_pg / A)**nu + E_IN = drydep_params%C_IN * (d_pg / drydep_params%A)**drydep_params%nu ! Stokes number based on integrated settling velocity - St = (V_g_hat * u_star) / (const%std_grav * A) + St = (V_g_hat * u_star) / (const%std_grav * drydep_params%A) ! Collection efficiency due to impaction - E_IM = C_IM * (St / (alpha + St))**beta + E_IM = drydep_params%C_IM * (St / (drydep_params%alpha + St))**drydep_params%beta ! Rebound correction R1 = exp(-St**0.5d0) ! Surface resistance - R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) + R_s = 1.0d0 / (drydep_params%eps_0 * u_star * (E_B + E_IN + E_IM) * R1) ! Integrated deposition velocity V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s)) @@ -767,39 +801,42 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep_quadpack(scenario, real(kind=dp) :: density_air real(kind=dp) :: visc_d, visc_k real(kind=dp) :: gas_speed, gas_mean_free_path - real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu - real(kind=dp) :: u_mean, u_star, z_ref, z_rough - real(kind=dp) :: C_B, C_IN, C_IM, M_k + real(kind=dp) :: u_star + real(kind=dp) :: M_k real(kind=dp) :: lower, upper, epsabs, epsrel, result, abserr, ref integer :: limit, neval, ier, last real(kind=dp), allocatable :: alist(:), blist(:), rlist(:), elist(:) integer, allocatable :: iord(:) - ! Hardcoded meteorological variables and - ! parameterization-dependent LUC values. - z_ref = 20.0d0 ! Reference height - u_mean = 5.0d0 ! Mean wind speed at reference height - ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 - z_rough = .8d0 - A = 2.0d0 / 1000.0d0 - alpha = 1.0d0 - eps_0 = 3.0d0 - - if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then - gamma = .56d0 ! LUC-dependent for Zhang. - C_B = .2d0 - C_IN = 2.5d0 - C_IM = .4d0 - nu = .8d0 - beta = 1.7d0 - else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then - gamma = .56d0 ! LUC-dependent for Zhang. - C_B = 1.0d0 - C_IN = .5d0 - C_IM = 1.0d0 - nu = 2.0d0 - beta = 2.0d0 - end if + type(drydep_params_t) :: drydep_params + + drydep_params = scenario%drydep + + ! ! Hardcoded meteorological variables and + ! ! parameterization-dependent LUC values. + ! z_ref = 20.0d0 ! Reference height + ! u_mean = 5.0d0 ! Mean wind speed at reference height + ! ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 + ! z_rough = .8d0 + ! A = 2.0d0 / 1000.0d0 + ! alpha = 1.0d0 + ! eps_0 = 3.0d0 + + ! if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then + ! gamma = .56d0 ! LUC-dependent for Zhang. + ! C_B = .2d0 + ! C_IN = 2.5d0 + ! C_IM = .4d0 + ! nu = .8d0 + ! beta = 1.7d0 + ! else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then + ! gamma = .56d0 ! LUC-dependent for Zhang. + ! C_B = 1.0d0 + ! C_IN = .5d0 + ! C_IM = 1.0d0 + ! nu = 2.0d0 + ! beta = 2.0d0 + ! end if ! particle diameter equal to geometric mean diameter d_pg = aero_mode%char_radius * 2.0d0 @@ -871,8 +908,8 @@ real(kind=dp) function dep_vel_integrand(d_p) V_s = (density * d_p**2.0d0 * const%std_grav * cunning) / (18.0d0 * visc_d) ! Friction velocity and aerodynamic resistance - u_star = 0.4d0 * u_mean / log(z_ref / z_rough) - R_a = (1.0d0 / (0.4d0 * u_star)) * log(z_ref / z_rough) + u_star = 0.4d0 * drydep_params%u_mean / log(drydep_params%z_ref / drydep_params%z_rough) + R_a = (1.0d0 / (0.4d0 * u_star)) * log(drydep_params%z_ref / drydep_params%z_rough) ! Particle diffusivity and Schmidt number diff_p = (const%boltzmann * env_state%temp * cunning) / & @@ -880,18 +917,18 @@ real(kind=dp) function dep_vel_integrand(d_p) Sc = visc_k / diff_p ! Collection efficiencies - E_B = C_B * Sc**(-gamma) - E_IN = C_IN * (d_p / A)**nu + E_B = drydep_params%C_B * Sc**(-drydep_params%gamma) + E_IN = drydep_params%C_IN * (d_p / drydep_params%A)**drydep_params%nu ! Stokes number and impaction efficiency - St = (V_s * u_star) / (const%std_grav * A) - E_IM = C_IM * (St / (alpha + St))**beta + St = (V_s * u_star) / (const%std_grav * drydep_params%A) + E_IM = C_IM * (St / (drydep_params%alpha + St))**drydep_params%beta ! Rebound correction R1 = exp(-St**0.5d0) ! Surface resistance - R_s = 1.0d0 / (eps_0 * u_star * (E_B + E_IN + E_IM) * R1) + R_s = 1.0d0 / (drydep_params%eps_0 * u_star * (E_B + E_IN + E_IM) * R1) ! Deposition velocity V_d = V_s + (1.0d0 / (R_a + R_s)) @@ -1186,6 +1223,7 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & character(len=PMC_MAX_FILENAME_LEN) :: sub_filename type(spec_file_t) :: sub_file character(len=SPEC_LINE_MAX_VAR_LEN) :: function_name + type(spec_line_t) :: line ! note that we have to hard-code the list for doxygen below @@ -1290,15 +1328,27 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_VOLUME else if (trim(function_name) == 'drydep') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_DRYDEP - call spec_file_read_string(file, "dry_dep_param", function_name) - if (trim(function_name) == 'Zhang_2001') then - scenario%drydep_param = SCENARIO_DRYDEP_ZHANG - else if (trim(function_name) == 'Emerson_2020') then - scenario%drydep_param = SCENARIO_DRYDEP_EMERSON - else - call die_msg(395827406, "Unknown dry deposition parameterization: " & - // trim(function_name)) + + call spec_file_read_line_no_eof(file, line) + call spec_file_unread_line(file) + if (line%name /= 'drydep_params') then + call warn_msg(735291468, "using default dry deposition parameters") + else + call spec_file_read_string(file, 'drydep_params', sub_filename) + call spec_file_open(sub_filename, sub_file) + call spec_file_read_drydep_params(sub_file, scenario%drydep) + call spec_file_close(sub_file) end if + + ! call spec_file_read_string(file, "dry_dep_param", function_name) + ! if (trim(function_name) == 'Zhang_2001') then + ! scenario%drydep_param = SCENARIO_DRYDEP_ZHANG + ! else if (trim(function_name) == 'Emerson_2020') then + ! scenario%drydep_param = SCENARIO_DRYDEP_EMERSON + ! else + ! call die_msg(395827406, "Unknown dry deposition parameterization: " & + ! // trim(function_name)) + ! end if else if (trim(function_name) == 'chamber') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_CHAMBER call spec_file_read_chamber(file, scenario%chamber) @@ -1412,32 +1462,62 @@ end subroutine spec_file_read_scenario !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Output the dry deposition parameterization to NetCDF file. - subroutine scenario_output_drydep_param(scenario, ncid) - - !> Scenario data. - type(scenario_t), intent(in) :: scenario - !> NetCDF file ID, in data mode. - integer, intent(in) :: ncid - - call pmc_nc_write_integer(ncid, scenario%drydep_param, & - "drydep_param") - - end subroutine scenario_output_drydep_param - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !> Input the dry deposition parameterization to NetCDF file. - subroutine scenario_input_drydep_param(scenario, ncid) - - !> Scenario data. - type(scenario_t), intent(inout) :: scenario - !> NetCDF file ID, in data mode. - integer, intent(in) :: ncid + !> Read dry deposition parameters from a spec file. + subroutine spec_file_read_drydep_params(file, dry_dep_params) - call pmc_nc_read_integer(ncid, scenario%drydep_param, "drydep_param") + !> Spec file. + type(spec_file_t), intent(inout) :: file + !> Dry deposition parameters. + type(drydep_params_t), intent(inout) :: dry_dep_params - end subroutine scenario_input_drydep_param + !> \page input_format_chamber Input File Format: Dry Deposition Parameters + !! + !! Dry deposition is simulatied using the specified parameters: + !! - \b drydep_param (integer): flag for the dry deposition parameterization + !! - \b z_ref (real, unit m): the reference height \f$z_{\rm ref}\f$ used + !! in the calculation of aerodynamic resistance \f$R_a\f$ + !! - \b u_mean (real, unit m s^{-1}): the wind speed at the reference + !! height + !! - \b z_rough (real, unit m): the roughness length associated with + !! the surface + !! - \b A (real, unit m): the characteristic radius of collectors + !! associated with the surface + !! - \b alpha (real, dimensionless): the parameter \f$\alpha\f$ used in the + !! calculation of impaction efficiency \f$E_{\rm IM}\f$ + !! - \b eps_0 (real, dimensionless): the empirical constant used in the + !! calculation of surface resistance \f$R_s\f$ + !! - \b gamma (real, dimensionless): the exponent \f$\gamma\f$ used in the + !! calculation of Brownian diffusion collection efficiency \f$E_{\rm B}\f$ + !! - \b C_B (real, dimensionless): the coefficient for Brownian diffusion + !! collection efficiency \f$E_{\rm B}\f$ + !! - \b C_IN (real, dimensionless): the coefficient for interception + !! collection efficiency \f$E_{\rm IN}\f$ + !! - \b C_IM (real, dimensionless): the coefficient for impaction collection + !! efficiency \f$E_{\rm IM}\f$ + !! - \b nu (real, dimensionless): the exponent \f$\nu\f$ used in the calculation + !! of interception collection efficiency \f$E_{\rm IN}\f$ + !! - \b beta (real, dimensionless): the exponent \f$\beta\f$ used in the + !! calculation of impaction collection efficiency \f$E_{\rm IM}\f$ + !! + !! See also: + !! - \ref spec_file_format --- the input file text format + !! - \ref input_format_scenario --- the prescribed profiles of + !! other environment data + + call spec_file_read_real(file, "z_ref", dry_dep_params%z_ref) + call spec_file_read_real(file, "u_mean", dry_dep_params%u_mean) + call spec_file_read_real(file, "z_rough", dry_dep_params%z_rough) + call spec_file_read_real(file, "A", dry_dep_params%A) + call spec_file_read_real(file, "alpha", dry_dep_params%alpha) + call spec_file_read_real(file, "eps_0", dry_dep_params%eps_0) + call spec_file_read_real(file, "gamma", dry_dep_params%gamma) + call spec_file_read_real(file, "C_B", dry_dep_params%C_B) + call spec_file_read_real(file, "C_IN", dry_dep_params%C_IN) + call spec_file_read_real(file, "C_IM", dry_dep_params%C_IM) + call spec_file_read_real(file, "nu", dry_dep_params%nu) + call spec_file_read_real(file, "beta", dry_dep_params%beta) + + end subroutine spec_file_read_drydep_params !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/spec_file.F90 b/src/spec_file.F90 index 746621fdc..7422c3701 100644 --- a/src/spec_file.F90 +++ b/src/spec_file.F90 @@ -463,6 +463,25 @@ subroutine spec_file_check_read_iostat(file, ios, type) end subroutine spec_file_check_read_iostat +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Set the line number and file read pointer back one entry. + subroutine spec_file_unread_line(file) + + !> Spec file. + type(spec_file_t), intent(inout) :: file + + integer :: ios + + file%line_num = file%line_num - 1 + backspace(file%unit, iostat=ios) + if (ios /= 0) then + call spec_file_die_msg(362985714, file, & + 'error backspacing: IOSTAT = ' // trim(integer_to_string(ios))) + end if + + end subroutine spec_file_unread_line + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Convert a string to an integer. From ca4142bce542de3072576df172597d6ca21acb2d Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 7 Sep 2025 14:19:02 -0500 Subject: [PATCH 16/28] remove debugging and commented code --- src/scenario.F90 | 95 ++++-------------------------------------------- 1 file changed, 7 insertions(+), 88 deletions(-) diff --git a/src/scenario.F90 b/src/scenario.F90 index 6f130ac8a..0747188b9 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -459,12 +459,12 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, & else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then ! loss do i_mode = 1,aero_dist_n_mode(aero_dist) - aero_mode = aero_dist%mode(i_mode) - N = aero_mode%num_conc + aero_mode = aero_dist%mode(i_mode) + N = aero_mode%num_conc - if (N == 0d0) then - cycle - end if + if (N == 0d0) then + cycle + end if d_pg = aero_mode%char_radius * 2.0d0 ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) @@ -484,7 +484,7 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, & ! Integrated deposition rate for the 3-rd moment (proportional to volume conc.) m_3_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, 3.0d0, & - density, env_state) + density, env_state) M = N * d_pg**3.0d0 * exp((3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)) new_M = M * exp(m_3_rate * del_t) @@ -568,26 +568,14 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & real(kind=dp) :: knud, cunning real(kind=dp) :: grav real(kind=dp) :: R_s, R_a - ! real(kind=dp) :: alpha, beta, gamma, A, eps_0 real(kind=dp) :: diff_p real(kind=dp) :: von_karman real(kind=dp) :: St, Sc, u_star real(kind=dp) :: E_B, E_IM, E_IN, R1 - ! real(kind=dp) :: u_mean, z_ref, z_rough type(drydep_params_t) :: drydep_params drydep_params = scenario%drydep - ! ! User set variables - ! u_mean = 5.0d0 ! Mean wind speed at reference height - ! z_ref = 20.0d0 ! Reference height - ! ! Setting for LUC = 7, SC = 1 - See Table 3 - ! z_rough = .1d0 ! According to land category - ! A = 2.0d0 / 1000.0d0 ! Dependent on land type - ! alpha = 1.2d0 ! From table - ! beta = 2.0d0 ! From text - ! gamma = .54d0 ! From table - ! particle diameter d_p = aero_data_vol2diam(aero_data, vol) ! density of air @@ -670,13 +658,11 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & real(kind=dp) :: visc_d, visc_k real(kind=dp) :: gas_speed, gas_mean_free_path real(kind=dp) :: knud - ! real(kind=dp) :: alpha, beta, gamma, A, eps_0, nu real(kind=dp) :: D_bar, D_hat real(kind=dp) :: St, Sc real(kind=dp) :: u_star real(kind=dp) :: R_a, R_s real(kind=dp) :: E_B, E_IN, E_IM, R1 - ! real(kind=dp) :: C_B, C_IN, C_IM type(drydep_params_t) :: drydep_params #ifdef PMC_USE_QUADPACK @@ -686,38 +672,8 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & return #endif - ! ! Hardcoded meteorological variables and - ! ! parameterization-dependent LUC values. - ! z_ref = 20.0d0 ! Reference height - ! u_mean = 5.0d0 ! Mean wind speed at reference height - ! ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 - ! z_rough = .8d0 - ! A = 2.0d0 / 1000.0d0 - ! alpha = 1.0d0 - ! eps_0 = 3.0d0 - - ! if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then - ! gamma = .56d0 ! LUC-dependent for Zhang. - ! C_B = .2d0 - ! C_IN = 2.5d0 - ! C_IM = .4d0 - ! nu = .8d0 - ! beta = 1.7d0 - ! else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then - ! gamma = .56d0 ! LUC-dependent for Zhang. - ! C_B = 1.0d0 - ! C_IN = .5d0 - ! C_IM = 1.0d0 - ! nu = 2.0d0 - ! beta = 2.0d0 - ! end if drydep_params = scenario%drydep - print*, drydep_params%z_ref, drydep_params%u_mean, drydep_params%z_rough, & - drydep_params%A, drydep_params%alpha, drydep_params%eps_0, & - drydep_params%gamma, drydep_params%C_B, drydep_params%C_IN, & - drydep_params%C_IM, drydep_params%nu, drydep_params%beta - ! particle diameter equal to geometric mean diameter d_pg = aero_mode%char_radius * 2.0d0 ! natural log of geometric standard deviation @@ -812,32 +768,6 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep_quadpack(scenario, drydep_params = scenario%drydep - ! ! Hardcoded meteorological variables and - ! ! parameterization-dependent LUC values. - ! z_ref = 20.0d0 ! Reference height - ! u_mean = 5.0d0 ! Mean wind speed at reference height - ! ! LUC 7 (crops, mixed farming) from Zhang et al., 2001 - ! z_rough = .8d0 - ! A = 2.0d0 / 1000.0d0 - ! alpha = 1.0d0 - ! eps_0 = 3.0d0 - - ! if (scenario%drydep_param == SCENARIO_DRYDEP_EMERSON) then - ! gamma = .56d0 ! LUC-dependent for Zhang. - ! C_B = .2d0 - ! C_IN = 2.5d0 - ! C_IM = .4d0 - ! nu = .8d0 - ! beta = 1.7d0 - ! else if (scenario%drydep_param == SCENARIO_DRYDEP_ZHANG) then - ! gamma = .56d0 ! LUC-dependent for Zhang. - ! C_B = 1.0d0 - ! C_IN = .5d0 - ! C_IM = 1.0d0 - ! nu = 2.0d0 - ! beta = 2.0d0 - ! end if - ! particle diameter equal to geometric mean diameter d_pg = aero_mode%char_radius * 2.0d0 ! natural log of geometric standard deviation @@ -947,7 +877,7 @@ end function scenario_integrated_loss_rate_dry_dep_quadpack !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Updates a given array to contain the integrated deposition rates for the + !> Updates an array to contain the integrated deposition rates for the !> given moment of each aerosol mode in the distribution. subroutine scenario_modal_dry_dep_rates(scenario, aero_dist, moment, & density, env_state, rates) @@ -1328,7 +1258,6 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_VOLUME else if (trim(function_name) == 'drydep') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_DRYDEP - call spec_file_read_line_no_eof(file, line) call spec_file_unread_line(file) if (line%name /= 'drydep_params') then @@ -1339,16 +1268,6 @@ subroutine spec_file_read_scenario(file, gas_data, aero_data, & call spec_file_read_drydep_params(sub_file, scenario%drydep) call spec_file_close(sub_file) end if - - ! call spec_file_read_string(file, "dry_dep_param", function_name) - ! if (trim(function_name) == 'Zhang_2001') then - ! scenario%drydep_param = SCENARIO_DRYDEP_ZHANG - ! else if (trim(function_name) == 'Emerson_2020') then - ! scenario%drydep_param = SCENARIO_DRYDEP_EMERSON - ! else - ! call die_msg(395827406, "Unknown dry deposition parameterization: " & - ! // trim(function_name)) - ! end if else if (trim(function_name) == 'chamber') then scenario%loss_function_type = SCENARIO_LOSS_FUNCTION_CHAMBER call spec_file_read_chamber(file, scenario%chamber) From a855f694a40bb45dbb00637f1a83b214396e62b4 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sat, 13 Sep 2025 09:47:04 -0500 Subject: [PATCH 17/28] update MPI routines to handle drydep parameters --- src/scenario.F90 | 97 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/src/scenario.F90 b/src/scenario.F90 index 0747188b9..6894ad884 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -1464,6 +1464,7 @@ integer function pmc_mpi_pack_size_scenario(val) + pmc_mpi_pack_size_real_array(val%aero_dilution_time) & + pmc_mpi_pack_size_real_array(val%aero_dilution_rate) & + pmc_mpi_pack_size_integer(val%loss_function_type) & + + pmc_mpi_pack_size_drydep(val%drydep) & + pmc_mpi_pack_size_chamber(val%chamber) if (allocated(val%gas_emission_time)) then do i = 1,size(val%gas_emission) @@ -1526,6 +1527,7 @@ subroutine pmc_mpi_pack_scenario(buffer, position, val) call pmc_mpi_pack_real_array(buffer, position, val%aero_dilution_time) call pmc_mpi_pack_real_array(buffer, position, val%aero_dilution_rate) call pmc_mpi_pack_integer(buffer, position, val%loss_function_type) + call pmc_mpi_pack_drydep(buffer, position, val%drydep) call pmc_mpi_pack_chamber(buffer, position, val%chamber) if (allocated(val%gas_emission_time)) then do i = 1,size(val%gas_emission) @@ -1586,6 +1588,7 @@ subroutine pmc_mpi_unpack_scenario(buffer, position, val) call pmc_mpi_unpack_real_array(buffer, position, val%aero_dilution_time) call pmc_mpi_unpack_real_array(buffer, position, val%aero_dilution_rate) call pmc_mpi_unpack_integer(buffer, position, val%loss_function_type) + call pmc_mpi_unpack_drydep(buffer, position, val%drydep) call pmc_mpi_unpack_chamber(buffer, position, val%chamber) if (allocated(val%gas_emission)) deallocate(val%gas_emission) if (allocated(val%gas_background)) deallocate(val%gas_background) @@ -1623,6 +1626,100 @@ subroutine pmc_mpi_unpack_scenario(buffer, position, val) end subroutine pmc_mpi_unpack_scenario +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Determines the number of bytes required to pack the given value. + integer function pmc_mpi_pack_size_drydep(val) + + !> Value to pack. + type(drydep_params_t), intent(in) :: val + + integer :: total_size, i, N + + pmc_mpi_pack_size_drydep = & + pmc_mpi_pack_size_real(val%z_ref) & + + pmc_mpi_pack_size_real(val%u_mean) & + + pmc_mpi_pack_size_real(val%z_rough) & + + pmc_mpi_pack_size_real(val%A) & + + pmc_mpi_pack_size_real(val%alpha) & + + pmc_mpi_pack_size_real(val%eps_0) & + + pmc_mpi_pack_size_real(val%gamma) & + + pmc_mpi_pack_size_real(val%C_B) & + + pmc_mpi_pack_size_real(val%C_IN) & + + pmc_mpi_pack_size_real(val%C_IM) & + + pmc_mpi_pack_size_real(val%nu) & + + pmc_mpi_pack_size_real(val%beta) + + end function pmc_mpi_pack_size_drydep + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Packs the given value into the buffer, advancing position. + subroutine pmc_mpi_pack_drydep(buffer, position, val) + + !> Memory buffer. + character, intent(inout) :: buffer(:) + !> Current buffer position. + integer, intent(inout) :: position + !> Value to pack. + type(drydep_params_t), intent(in) :: val + +#ifdef PMC_USE_MPI + integer :: prev_position + + prev_position = position + call pmc_mpi_pack_real(buffer, position, val%z_ref) + call pmc_mpi_pack_real(buffer, position, val%u_mean) + call pmc_mpi_pack_real(buffer, position, val%z_rough) + call pmc_mpi_pack_real(buffer, position, val%A) + call pmc_mpi_pack_real(buffer, position, val%alpha) + call pmc_mpi_pack_real(buffer, position, val%eps_0) + call pmc_mpi_pack_real(buffer, position, val%gamma) + call pmc_mpi_pack_real(buffer, position, val%C_B) + call pmc_mpi_pack_real(buffer, position, val%C_IN) + call pmc_mpi_pack_real(buffer, position, val%C_IM) + call pmc_mpi_pack_real(buffer, position, val%nu) + call pmc_mpi_pack_real(buffer, position, val%beta) + call assert(371592468, & + position - prev_position <= pmc_mpi_pack_size_drydep(val)) +#endif + + end subroutine pmc_mpi_pack_drydep + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Unpacks the given value from the buffer, advancing position. + subroutine pmc_mpi_unpack_drydep(buffer, position, val) + + !> Memory buffer. + character, intent(inout) :: buffer(:) + !> Current buffer position. + integer, intent(inout) :: position + !> Value to pack. + type(drydep_params_t), intent(inout) :: val + +#ifdef PMC_USE_MPI + integer :: prev_position + + prev_position = position + call pmc_mpi_unpack_real(buffer, position, val%z_ref) + call pmc_mpi_unpack_real(buffer, position, val%u_mean) + call pmc_mpi_unpack_real(buffer, position, val%z_rough) + call pmc_mpi_unpack_real(buffer, position, val%A) + call pmc_mpi_unpack_real(buffer, position, val%alpha) + call pmc_mpi_unpack_real(buffer, position, val%eps_0) + call pmc_mpi_unpack_real(buffer, position, val%gamma) + call pmc_mpi_unpack_real(buffer, position, val%C_B) + call pmc_mpi_unpack_real(buffer, position, val%C_IN) + call pmc_mpi_unpack_real(buffer, position, val%C_IM) + call pmc_mpi_unpack_real(buffer, position, val%nu) + call pmc_mpi_unpack_real(buffer, position, val%beta) + call assert(582741936, & + position - prev_position <= pmc_mpi_pack_size_drydep(val)) +#endif + + end subroutine pmc_mpi_unpack_drydep + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module pmc_scenario From 093b3349485d0e5deeb5ca73a09e9731168eda3a Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 12:41:20 -0500 Subject: [PATCH 18/28] update processing script for modal runs and add exectuable --- CMakeLists.txt | 9 +- scenarios/7_drydep/drydep_modal_process.F90 | 106 ++++++++++++++++++++ 2 files changed, 111 insertions(+), 4 deletions(-) create mode 100644 scenarios/7_drydep/drydep_modal_process.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index e62d1a3f4..f8c6153b9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -433,10 +433,11 @@ add_executable(chamber_process target_link_libraries(chamber_process partmclib) ###################################################################### -# scenarios/7_dry_dep/dry_dep_modal_process +# scenarios/7_drydep/drydep_modal_process -# add_executable(dry_dep_modal_process -# scenarios/7_dry_dep/dry_dep_modal_process.F90) -# target_link_libraries(dry_dep_modal_process partmclib) +add_executable(drydep_modal_process + scenarios/7_drydep/drydep_modal_process.F90 + scenarios/7_drydep/drydep_modal_process.F90) +target_link_libraries(drydep_modal_process partmclib) ###################################################################### \ No newline at end of file diff --git a/scenarios/7_drydep/drydep_modal_process.F90 b/scenarios/7_drydep/drydep_modal_process.F90 new file mode 100644 index 000000000..174aa4181 --- /dev/null +++ b/scenarios/7_drydep/drydep_modal_process.F90 @@ -0,0 +1,106 @@ +!> Read NetCDF output files from a modal runn and process them. + +program process + + use pmc_output + use pmc_stats + use pmc_aero_dist + use pmc_scenario + + character(len=PMC_MAX_FILENAME_LEN), parameter :: prefix & += "out/modal_test" + + character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename + type(bin_grid_t) :: bin_grid + type(aero_dist_t) :: aero_dist + type(aero_binned_t) :: aero_binned + type(aero_data_t) :: aero_data + type(env_state_t) :: env_state + type(gas_data_t) :: gas_data + type(gas_state_t) :: gas_state + type(scenario_t) :: scenario + integer :: ncid, index, i_mode, i_index, dum, n_index + real(kind=dp) :: time, del_t, tot_num_conc, density, tot_mass_conc + character(len=PMC_UUID_LEN) :: uuid + real(kind=dp), allocatable :: times(:), num_conc(:) + type(stats_1d_t) :: stats_tot_num_conc, stats_num_conc, stats_rates_0, & + stats_rates_3, stats_tot_mass_conc + real(kind=dp), allocatable :: rates_0(:), rates_3(:) + character(len=PMC_MAX_FILENAME_LEN), allocatable :: file_list(:) + + call pmc_mpi_init() + + call input_filename_list(prefix, file_list) + n_index = size(file_list) + + allocate(times(n_index)) + + do i_index = 1,n_index + call make_filename(in_filename, prefix, ".nc", i_index) + write(*,*) "Processing " // trim(in_filename) + call input_modal(in_filename, index, time, del_t, uuid, aero_dist=aero_dist, & + aero_binned=aero_binned, aero_data=aero_data, env_state=env_state, & + gas_data=gas_data, gas_state=gas_state, bin_grid=bin_grid, scenario=scenario) + times(i_index) = time + + density = aero_data%density(1) + + num_conc = aero_binned%num_conc * bin_grid%widths + tot_num_conc = sum(num_conc) + + tot_mass_conc = sum(aero_binned%vol_conc(:,1) * bin_grid%widths) * aero_data%density(1) + + if (.not. allocated(rates_0)) allocate(rates_0(aero_dist_n_mode(aero_dist))) + if (.not. allocated(rates_3)) allocate(rates_3(aero_dist_n_mode(aero_dist))) + + call scenario_modal_dry_dep_rates(scenario, aero_dist, 0.0d0, density, & + env_state, rates_0) + call stats_1d_add(stats_rates_0, rates_0) + call scenario_modal_dry_dep_rates(scenario, aero_dist, 3.0d0, density, & + env_state, rates_3) + call stats_1d_add(stats_rates_3, rates_3) + + call stats_1d_add(stats_num_conc, num_conc) + call stats_1d_add_entry(stats_tot_num_conc, tot_num_conc, i_index) + + call stats_1d_add_entry(stats_tot_mass_conc, tot_mass_conc, i_index) + + call make_filename(out_filename, prefix, "_process.nc", index) + write(*,*) "Writing " // trim(out_filename) + call pmc_nc_open_write(out_filename, ncid) + call pmc_nc_write_info(ncid, uuid, "1_urban_plume modal process") + call env_state_output_netcdf(env_state, ncid) + call aero_data_output_netcdf(aero_data, ncid) + call aero_dist_output_netcdf(aero_dist, ncid) + call aero_binned_output_netcdf(aero_binned, ncid, bin_grid, aero_data) + + call stats_1d_output_netcdf(stats_num_conc, ncid, "num concs. per bin", & + dim_name="diam", unit="m^{-3}") + call stats_1d_output_netcdf(stats_rates_0, ncid, "loss_rates_0", & + dim_name="modes", unit="m s^{-1}") + call stats_1d_output_netcdf(stats_rates_3, ncid, "loss_rates_3", & + dim_name="modes", unit="m s^{-1}") + + call aero_binned_zero(aero_binned) + + call stats_1d_clear(stats_num_conc) + call stats_1d_clear(stats_rates_0) + call stats_1d_clear(stats_rates_3) + + call pmc_nc_close(ncid) + end do + + call make_filename(out_filename, prefix, "_process.nc") + write(*,*) "Writing " // trim(out_filename) + call pmc_nc_open_write(out_filename, ncid) + call pmc_nc_write_info(ncid, uuid, "1_urban_plume modal process") + call pmc_nc_write_real_1d(ncid, times, "time", dim_name="time", unit="s") + call stats_1d_output_netcdf(stats_tot_num_conc, ncid, "tot_num_conc", & + dim_name="time", unit="m^{-3}") + call stats_1d_output_netcdf(stats_tot_mass_conc, ncid, "tot_mass_conc", & + dim_name="time", unit="ug m^{-3}") + call pmc_nc_close(ncid) + + call pmc_mpi_finalize() + +end program process \ No newline at end of file From cc7ff7e50c032dbb11bb40c1b18d1a6134f3318b Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 12:42:44 -0500 Subject: [PATCH 19/28] update aero_dist netcdf input and output for modal runs to operate with arrays --- src/aero_dist.F90 | 82 +++++++++++++++++++---------------------------- 1 file changed, 33 insertions(+), 49 deletions(-) diff --git a/src/aero_dist.F90 b/src/aero_dist.F90 index 80b733b1f..552a868d5 100644 --- a/src/aero_dist.F90 +++ b/src/aero_dist.F90 @@ -281,51 +281,32 @@ subroutine aero_dist_output_netcdf(aero_dist, ncid) !> NetCDF file ID, in data mode. integer, intent(in) :: ncid - integer :: i_mode - type(aero_mode_t) :: aero_mode - character(len=1) :: num + integer :: i_mode, n_mode + real(kind=dp), allocatable :: char_rads(:), log10_std_dev_rads(:), num_concs(:) - do i_mode = 1,aero_dist_n_mode(aero_dist) - aero_mode = aero_dist%mode(i_mode) - write(num, '(I1.1)') i_mode - call pmc_nc_write_real(ncid, aero_mode%char_radius, "char_radius_mode_" & - // trim(num), unit="m") - call pmc_nc_write_real(ncid, aero_mode%log10_std_dev_radius, "std_dev_mode_" & - // trim(num), unit="1", description="log base 10 of geometric standard" & - // " deviation of radius") - call pmc_nc_write_real(ncid, aero_mode%num_conc, "num_conc_mode_" & - // trim(num), unit="m^-3", description="total number concentration for mode") + n_mode = aero_dist_n_mode(aero_dist) + + allocate(char_rads(n_mode), log10_std_dev_rads(n_mode), num_concs(n_mode)) + + do i_mode = 1,n_mode + char_rads(i_mode) = aero_dist%mode(i_mode)%char_radius + log10_std_dev_rads(i_mode) = aero_dist%mode(i_mode)%log10_std_dev_radius + num_concs(i_mode) = aero_dist%mode(i_mode)%num_conc end do + call pmc_nc_write_real_1d(ncid, char_rads, "char_radius_modes", dim_name="num_modes", & + long_name="characteristic radius of each mode", & + unit="m") + call pmc_nc_write_real_1d(ncid, log10_std_dev_rads, "std_dev_radius_modes", dim_name="num_modes", & + long_name="log_10 of geometric std dev for each mode", & + unit="1") + call pmc_nc_write_real_1d(ncid, num_concs, "num_conc_modes", dim_name="num_modes", & + long_name="number concentration for each mode", & + unit="m^-3") call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "num_modes", & description="total number of modes") - ! integer :: i_mode, n_mode - ! real(kind=dp), allocatable :: char_rads(:), log10_std_dev_rads(:), num_concs(:) - - ! n_mode = aero_dist_n_mode(aero_dist) - - ! allocate(char_rads(n_mode), log10_std_dev_rads(n_mode), num_concs(n_mode)) - - ! do i_mode = 1,n_mode - ! char_rads(i_mode) = aero_dist%mode(i_mode)%char_radius - ! log10_std_dev_rads(i_mode) = aero_dist%mode(i_mode)%log10_std_dev_radius - ! num_concs(i_mode) = aero_dist%mode(i_mode)%num_conc - ! end do - - ! call pmc_nc_write_real_1d(ncid, char_rads, "char_radius_modes", dim_name="num_modes", & - ! long_name="characteristic radius of each mode", & - ! unit="m") - ! call pmc_nc_write_real_1d(ncid, log10_std_dev_rads, "std_dev_radius_modes", dim_name="num_modes", & - ! long_name="log_10 of geometric std dev for each mode ", & - ! unit="1") - ! call pmc_nc_write_real_1d(ncid, num_concs, "num_conc_modes", dim_name="num_modes", & - ! long_name="number concentration for each mode", & - ! unit="m^-3") - ! call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "num_modes", & - ! dimdescription="total number of modes") - - ! deallocate(char_rads, log10_std_dev_rads, num_concs) + deallocate(char_rads, log10_std_dev_rads, num_concs) end subroutine aero_dist_output_netcdf @@ -340,27 +321,30 @@ subroutine aero_dist_input_netcdf(aero_dist, ncid) integer, intent(in) :: ncid integer :: i_mode, n_mode - character(len=1) :: num type(aero_mode_t), allocatable :: modes(:) + real(kind=dp), allocatable :: char_rads(:), log10_std_dev_rads(:), num_concs(:) call pmc_nc_read_integer(ncid, n_mode, "num_modes") - call assert(836472910, n_mode < 100) if (allocated(modes)) deallocate(modes) allocate(modes(n_mode)) - do i_mode=1,n_mode - write(num, '(I1.1)') i_mode - call pmc_nc_read_real(ncid, modes(i_mode)%char_radius, "char_radius_mode_" & - // trim(num)) - call pmc_nc_read_real(ncid, modes(i_mode)%log10_std_dev_radius, "std_dev_mode_" & - // trim(num)) - call pmc_nc_read_real(ncid, modes(i_mode)%num_conc, "num_conc_mode_" & - // trim(num)) + allocate(char_rads(n_mode), log10_std_dev_rads(n_mode), num_concs(n_mode)) + + call pmc_nc_read_real_1d(ncid, char_rads, "char_radius_modes") + call pmc_nc_read_real_1d(ncid, log10_std_dev_rads, "std_dev_radius_modes") + call pmc_nc_read_real_1d(ncid, num_concs, "num_conc_modes") + + do i_mode = 1, n_mode + modes(i_mode)%char_radius = char_rads(i_mode) + modes(i_mode)%log10_std_dev_radius = log10_std_dev_rads(i_mode) + modes(i_mode)%num_conc = num_concs(i_mode) end do aero_dist%mode = modes + deallocate(char_rads, log10_std_dev_rads, num_concs) + end subroutine aero_dist_input_netcdf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! From 6f186e912cedc96c8310b67732a5bf47e2921404 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 12:46:34 -0500 Subject: [PATCH 20/28] remove duplicated line --- CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f8c6153b9..1a0c54df7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -436,7 +436,6 @@ target_link_libraries(chamber_process partmclib) # scenarios/7_drydep/drydep_modal_process add_executable(drydep_modal_process - scenarios/7_drydep/drydep_modal_process.F90 scenarios/7_drydep/drydep_modal_process.F90) target_link_libraries(drydep_modal_process partmclib) From cd9d90e8644e80c69894d1352ae6a46ea2044636 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 13:08:24 -0500 Subject: [PATCH 21/28] update dockerignore --- .dockerignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.dockerignore b/.dockerignore index e725e276f..ba64e72f6 100644 --- a/.dockerignore +++ b/.dockerignore @@ -25,3 +25,5 @@ scenarios/4_chamber/out scenarios/5_coag_brownian/out !scenarios/6_camp scenarios/6_camp/out +!scenarios/7_drydep/ +scenarios/7_drydep/out From ed7219703485ff7dc69a6edc4920008d7622a132 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 13:16:17 -0500 Subject: [PATCH 22/28] dry_dep --> drydep --- scenarios/7_drydep/drydep_modal_process.F90 | 4 +- src/scenario.F90 | 58 ++++++++++----------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/scenarios/7_drydep/drydep_modal_process.F90 b/scenarios/7_drydep/drydep_modal_process.F90 index 174aa4181..1648542e3 100644 --- a/scenarios/7_drydep/drydep_modal_process.F90 +++ b/scenarios/7_drydep/drydep_modal_process.F90 @@ -53,10 +53,10 @@ program process if (.not. allocated(rates_0)) allocate(rates_0(aero_dist_n_mode(aero_dist))) if (.not. allocated(rates_3)) allocate(rates_3(aero_dist_n_mode(aero_dist))) - call scenario_modal_dry_dep_rates(scenario, aero_dist, 0.0d0, density, & + call scenario_modal_drydep_rates(scenario, aero_dist, 0.0d0, density, & env_state, rates_0) call stats_1d_add(stats_rates_0, rates_0) - call scenario_modal_dry_dep_rates(scenario, aero_dist, 3.0d0, density, & + call scenario_modal_drydep_rates(scenario, aero_dist, 3.0d0, density, & env_state, rates_3) call stats_1d_add(stats_rates_3, rates_3) diff --git a/src/scenario.F90 b/src/scenario.F90 index 6894ad884..0934af908 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -470,7 +470,7 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, & ln_sigma_g = aero_mode%log10_std_dev_radius / log10(exp(1.0d0)) ! Integrated deposition rate for the 0-th moment (num. conc.) - m_0_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, & + m_0_rate = -1.0d0 * scenario_integrated_loss_rate_drydep(scenario, aero_mode, & 0.0d0, density, env_state) new_N = N * exp(m_0_rate * del_t) @@ -483,7 +483,7 @@ subroutine scenario_update_aero_modes(aero_dist, del_t, env_state, & aero_dist%mode(i_mode)%num_conc = new_N ! Integrated deposition rate for the 3-rd moment (proportional to volume conc.) - m_3_rate = -1.0d0 * scenario_integrated_loss_rate_dry_dep(scenario, aero_mode, 3.0d0, & + m_3_rate = -1.0d0 * scenario_integrated_loss_rate_drydep(scenario, aero_mode, 3.0d0, & density, env_state) M = N * d_pg**3.0d0 * exp((3.0d0**2.0d0)/2.0d0 * (ln_sigma_g**2.0d0)) new_M = M * exp(m_3_rate * del_t) @@ -525,7 +525,7 @@ real(kind=dp) function scenario_loss_rate(scenario, vol, density, & scenario_loss_rate = 1d15*vol else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) & then - scenario_loss_rate = scenario_loss_rate_dry_dep(vol, density, & + scenario_loss_rate = scenario_loss_rate_drydep(vol, density, & aero_data, env_state, scenario) else if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_CHAMBER) & then @@ -545,7 +545,7 @@ end function scenario_loss_rate !> Compute and return the dry deposition rate for a given particle. !! All equations used here are written in detail in the file !! \c doc/deposition/deposition.tex. - real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & + real(kind=dp) function scenario_loss_rate_drydep(vol, density, aero_data, & env_state, scenario) !> Particle volume (m^3). @@ -629,15 +629,15 @@ real(kind=dp) function scenario_loss_rate_dry_dep(vol, density, aero_data, & V_d = V_s + (1.0d0 / (R_a + R_s + R_a * R_s * V_s)) ! The loss rate - scenario_loss_rate_dry_dep = V_d / env_state%height + scenario_loss_rate_drydep = V_d / env_state%height - end function scenario_loss_rate_dry_dep + end function scenario_loss_rate_drydep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Compute and return the integrated dry deposition rate for a given !> lognormal aerosol mode (modal approximation). - real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & + real(kind=dp) function scenario_integrated_loss_rate_drydep(scenario, & aero_mode, moment, density, env_state) !> Scenario data. @@ -667,7 +667,7 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & #ifdef PMC_USE_QUADPACK ! Do numerical integration when QUADPACK is available - scenario_integrated_loss_rate_dry_dep = scenario_integrated_loss_rate_dry_dep_quadpack( & + scenario_integrated_loss_rate_drydep = scenario_integrated_loss_rate_drydep_quadpack( & scenario, aero_mode, moment, density, env_state) return #endif @@ -732,14 +732,14 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep(scenario, & V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s)) ! Loss rate - scenario_integrated_loss_rate_dry_dep = V_d_hat / env_state%height + scenario_integrated_loss_rate_drydep = V_d_hat / env_state%height - end function scenario_integrated_loss_rate_dry_dep + end function scenario_integrated_loss_rate_drydep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifdef PMC_USE_QUADPACK - real(kind=dp) function scenario_integrated_loss_rate_dry_dep_quadpack(scenario, & + real(kind=dp) function scenario_integrated_loss_rate_drydep_quadpack(scenario, & aero_mode, moment, density, env_state) !> Scenario data. @@ -810,7 +810,7 @@ real(kind=dp) function scenario_integrated_loss_rate_dry_dep_quadpack(scenario, M_k = d_pg**moment * exp(moment**2 * ln_sigma_g**2 / 2.0d0) ! Integration result - scenario_integrated_loss_rate_dry_dep_quadpack = 1.0d0 / M_k * result / env_state%height + scenario_integrated_loss_rate_drydep_quadpack = 1.0d0 / M_k * result / env_state%height ! Clean up deallocate(alist, blist, rlist, elist, iord) @@ -872,14 +872,14 @@ real(kind=dp) function dep_vel_integrand(d_p) ! Final integrand dep_vel_integrand = d_p**moment * V_d * n_ddp end function dep_vel_integrand - end function scenario_integrated_loss_rate_dry_dep_quadpack + end function scenario_integrated_loss_rate_drydep_quadpack #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Updates an array to contain the integrated deposition rates for the !> given moment of each aerosol mode in the distribution. - subroutine scenario_modal_dry_dep_rates(scenario, aero_dist, moment, & + subroutine scenario_modal_drydep_rates(scenario, aero_dist, moment, & density, env_state, rates) !> Scenario data. @@ -898,7 +898,7 @@ subroutine scenario_modal_dry_dep_rates(scenario, aero_dist, moment, & integer :: i_mode, n_mode do i_mode = 1,size(rates) - rates(i_mode) = scenario_integrated_loss_rate_dry_dep( & + rates(i_mode) = scenario_integrated_loss_rate_drydep( & scenario, aero_dist%mode(i_mode), moment, density, env_state) & * env_state%height end do @@ -1382,12 +1382,12 @@ end subroutine spec_file_read_scenario !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Read dry deposition parameters from a spec file. - subroutine spec_file_read_drydep_params(file, dry_dep_params) + subroutine spec_file_read_drydep_params(file, drydep_params) !> Spec file. type(spec_file_t), intent(inout) :: file !> Dry deposition parameters. - type(drydep_params_t), intent(inout) :: dry_dep_params + type(drydep_params_t), intent(inout) :: drydep_params !> \page input_format_chamber Input File Format: Dry Deposition Parameters !! @@ -1423,18 +1423,18 @@ subroutine spec_file_read_drydep_params(file, dry_dep_params) !! - \ref input_format_scenario --- the prescribed profiles of !! other environment data - call spec_file_read_real(file, "z_ref", dry_dep_params%z_ref) - call spec_file_read_real(file, "u_mean", dry_dep_params%u_mean) - call spec_file_read_real(file, "z_rough", dry_dep_params%z_rough) - call spec_file_read_real(file, "A", dry_dep_params%A) - call spec_file_read_real(file, "alpha", dry_dep_params%alpha) - call spec_file_read_real(file, "eps_0", dry_dep_params%eps_0) - call spec_file_read_real(file, "gamma", dry_dep_params%gamma) - call spec_file_read_real(file, "C_B", dry_dep_params%C_B) - call spec_file_read_real(file, "C_IN", dry_dep_params%C_IN) - call spec_file_read_real(file, "C_IM", dry_dep_params%C_IM) - call spec_file_read_real(file, "nu", dry_dep_params%nu) - call spec_file_read_real(file, "beta", dry_dep_params%beta) + call spec_file_read_real(file, "z_ref", drydep_params%z_ref) + call spec_file_read_real(file, "u_mean", drydep_params%u_mean) + call spec_file_read_real(file, "z_rough", drydep_params%z_rough) + call spec_file_read_real(file, "A", drydep_params%A) + call spec_file_read_real(file, "alpha", drydep_params%alpha) + call spec_file_read_real(file, "eps_0", drydep_params%eps_0) + call spec_file_read_real(file, "gamma", drydep_params%gamma) + call spec_file_read_real(file, "C_B", drydep_params%C_B) + call spec_file_read_real(file, "C_IN", drydep_params%C_IN) + call spec_file_read_real(file, "C_IM", drydep_params%C_IM) + call spec_file_read_real(file, "nu", drydep_params%nu) + call spec_file_read_real(file, "beta", drydep_params%beta) end subroutine spec_file_read_drydep_params From d2191d547082025ec107be79d3b39b7f82980f53 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sun, 14 Sep 2025 13:39:44 -0500 Subject: [PATCH 23/28] add input files and scripts for modal drydep runs --- scenarios/7_drydep/1_run_modal.sh | 12 ++++ scenarios/7_drydep/2_process_modal.sh | 10 ++++ scenarios/7_drydep/aero_back.dat | 6 ++ scenarios/7_drydep/aero_back_dist.dat | 1 + scenarios/7_drydep/aero_data.dat | 2 + scenarios/7_drydep/aero_emit.dat | 6 ++ scenarios/7_drydep/aero_emit_dist.dat | 1 + scenarios/7_drydep/aero_init_comp.dat | 2 + scenarios/7_drydep/aero_init_dist.dat | 7 +++ scenarios/7_drydep/drydep_modal.spec | 33 +++++++++++ scenarios/7_drydep/drydep_params.dat | 12 ++++ scenarios/7_drydep/gas_back.dat | 79 +++++++++++++++++++++++++++ scenarios/7_drydep/gas_data.dat | 78 ++++++++++++++++++++++++++ scenarios/7_drydep/gas_emit.dat | 26 +++++++++ scenarios/7_drydep/height.dat | 5 ++ scenarios/7_drydep/pres.dat | 4 ++ scenarios/7_drydep/temp.dat | 4 ++ 17 files changed, 288 insertions(+) create mode 100755 scenarios/7_drydep/1_run_modal.sh create mode 100755 scenarios/7_drydep/2_process_modal.sh create mode 100644 scenarios/7_drydep/aero_back.dat create mode 100644 scenarios/7_drydep/aero_back_dist.dat create mode 100644 scenarios/7_drydep/aero_data.dat create mode 100644 scenarios/7_drydep/aero_emit.dat create mode 100644 scenarios/7_drydep/aero_emit_dist.dat create mode 100644 scenarios/7_drydep/aero_init_comp.dat create mode 100644 scenarios/7_drydep/aero_init_dist.dat create mode 100644 scenarios/7_drydep/drydep_modal.spec create mode 100644 scenarios/7_drydep/drydep_params.dat create mode 100644 scenarios/7_drydep/gas_back.dat create mode 100644 scenarios/7_drydep/gas_data.dat create mode 100644 scenarios/7_drydep/gas_emit.dat create mode 100644 scenarios/7_drydep/height.dat create mode 100644 scenarios/7_drydep/pres.dat create mode 100644 scenarios/7_drydep/temp.dat diff --git a/scenarios/7_drydep/1_run_modal.sh b/scenarios/7_drydep/1_run_modal.sh new file mode 100755 index 000000000..9354c9fc5 --- /dev/null +++ b/scenarios/7_drydep/1_run_modal.sh @@ -0,0 +1,12 @@ +#!/bin/sh + +# exit on error +set -e +# turn on command echoing +set -v + +mkdir -p out + +../../build/partmc drydep_modal.spec + +# Now run ./2_process_drydep_modal.sh to process the data \ No newline at end of file diff --git a/scenarios/7_drydep/2_process_modal.sh b/scenarios/7_drydep/2_process_modal.sh new file mode 100755 index 000000000..6c92ff47e --- /dev/null +++ b/scenarios/7_drydep/2_process_modal.sh @@ -0,0 +1,10 @@ +#!/bin/sh + +# exit on error +set -e +# turn on command echoing +set -v + +# The data should have already been generated by ./1_run_dry_dep_modal.sh + +../../build/dry_dep_modal_process \ No newline at end of file diff --git a/scenarios/7_drydep/aero_back.dat b/scenarios/7_drydep/aero_back.dat new file mode 100644 index 000000000..72f0e739a --- /dev/null +++ b/scenarios/7_drydep/aero_back.dat @@ -0,0 +1,6 @@ +# time (s) +# rate (s^{-1}) +# aerosol distribution filename +time 0 +rate 0 +dist aero_back_dist.dat \ No newline at end of file diff --git a/scenarios/7_drydep/aero_back_dist.dat b/scenarios/7_drydep/aero_back_dist.dat new file mode 100644 index 000000000..ee6d35bde --- /dev/null +++ b/scenarios/7_drydep/aero_back_dist.dat @@ -0,0 +1 @@ +# no background distribution \ No newline at end of file diff --git a/scenarios/7_drydep/aero_data.dat b/scenarios/7_drydep/aero_data.dat new file mode 100644 index 000000000..c6559d053 --- /dev/null +++ b/scenarios/7_drydep/aero_data.dat @@ -0,0 +1,2 @@ +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) +SO4 1800 0 96d-3 0.65 diff --git a/scenarios/7_drydep/aero_emit.dat b/scenarios/7_drydep/aero_emit.dat new file mode 100644 index 000000000..2cb60a378 --- /dev/null +++ b/scenarios/7_drydep/aero_emit.dat @@ -0,0 +1,6 @@ +# time (s) +# rate (s^{-1}) +# aerosol distribution filename +time 0 +rate 0 +dist aero_emit_dist.dat diff --git a/scenarios/7_drydep/aero_emit_dist.dat b/scenarios/7_drydep/aero_emit_dist.dat new file mode 100644 index 000000000..57995ab14 --- /dev/null +++ b/scenarios/7_drydep/aero_emit_dist.dat @@ -0,0 +1 @@ +# no emissions \ No newline at end of file diff --git a/scenarios/7_drydep/aero_init_comp.dat b/scenarios/7_drydep/aero_init_comp.dat new file mode 100644 index 000000000..cbe61c6c9 --- /dev/null +++ b/scenarios/7_drydep/aero_init_comp.dat @@ -0,0 +1,2 @@ +# mass fractions +SO4 1 diff --git a/scenarios/7_drydep/aero_init_dist.dat b/scenarios/7_drydep/aero_init_dist.dat new file mode 100644 index 000000000..721acd06a --- /dev/null +++ b/scenarios/7_drydep/aero_init_dist.dat @@ -0,0 +1,7 @@ +mode_name init_mode +mass_frac aero_init_comp.dat # composition proportions of species +diam_type geometric # type of diameter specified +mode_type log_normal # type of distribution +num_conc 3.2e9 # particle number concentration (#/m^3) +geom_mean_diam .00000001000000000000 +log10_geom_std_dev .04139268515822504074 \ No newline at end of file diff --git a/scenarios/7_drydep/drydep_modal.spec b/scenarios/7_drydep/drydep_modal.spec new file mode 100644 index 000000000..193865999 --- /dev/null +++ b/scenarios/7_drydep/drydep_modal.spec @@ -0,0 +1,33 @@ +run_type modal # modal run +output_prefix out/drydep_modal # prefix for output files + +t_max 28800 # total simulation time (s) +del_t 60 # timestep (s) +t_output 3600 # output interval (0 disables) (s) +t_progress 0 # progress printing interval (0 disables) (s) + +n_bin 1000 # number of bins (for processing purposes) +d_min .00000040000000000000e-1 +d_max .00025000000000000000e1 + +gas_data gas_data.dat # file containing gas data +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pres.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas concentrations file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function drydep # loss function specification +drydep_params drydep_params.dat # dry deposition parameters + +rel_humidity 0.95 # initial relative humidity (1) +latitude 0 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 21600 # start time (s since 00:00 UTC) +start_day 200 # start day of year (UTC) \ No newline at end of file diff --git a/scenarios/7_drydep/drydep_params.dat b/scenarios/7_drydep/drydep_params.dat new file mode 100644 index 000000000..b97f6c4ce --- /dev/null +++ b/scenarios/7_drydep/drydep_params.dat @@ -0,0 +1,12 @@ +z_ref 10.0 # reference height (m) +u_mean 2.0 # mean wind speed at the reference height (m) +z_rough 3.8 # surface roughness length (m) +A 0.0024 # characteristic radius of collectors (m) +alpha 5.0 +eps_0 6.0 +gamma .756 +C_B 0.82 +C_IN 92.5 +C_IM 10.4 +nu 11.8 +beta 12.7 \ No newline at end of file diff --git a/scenarios/7_drydep/gas_back.dat b/scenarios/7_drydep/gas_back.dat new file mode 100644 index 000000000..8c540b3b6 --- /dev/null +++ b/scenarios/7_drydep/gas_back.dat @@ -0,0 +1,79 @@ +# time (s) +# rate (s^{-1}) +# concentrations (ppb) +time 0 +rate 1.5e-5 +NO 0.1E+00 +NO2 1.0E+00 +NO3 0.0E+00 +N2O5 0.0E+00 +HONO 0.0E+00 +HNO3 1.0E+00 +HNO4 0.0E+00 +O3 5.0E+01 +O1D 0.0E+00 +O3P 0.0E+00 +OH 0.0E+00 +HO2 0.0E+00 +H2O2 1.1E+00 +CO 2.1E+02 +SO2 0.8E+00 +H2SO4 0.0E+00 +NH3 0.5E+00 +HCl 0.7E+00 +CH4 2.2E+03 +C2H6 1.0E+00 +CH3O2 0.0E+00 +ETHP 0.0E+00 +HCHO 1.2E+00 +CH3OH 1.2E-01 +CH3OOH 0.5E+00 +ETHOOH 0.0E+00 +ALD2 1.0E+00 +HCOOH 0.0E+00 +PAR 2.0E+00 +AONE 1.0E+00 +MGLY 0.0E+00 +ETH 0.2E+00 +OLET 2.3E-02 +OLEI 3.1E-04 +TOL 0.1E+00 +XYL 0.1E+00 +CRES 0.0E+00 +TO2 0.0E+00 +CRO 0.0E+00 +OPEN 0.0E+00 +ONIT 0.1E+00 +PAN 0.8E+00 +RCOOH 0.2E+00 +ROOH 2.5E-02 +C2O3 0.0E+00 +RO2 0.0E+00 +ANO2 0.0E+00 +NAP 0.0E+00 +ARO1 0.0E+00 +ARO2 0.0E+00 +ALK1 0.0E+00 +OLE1 0.0E+00 +XO2 0.0E+00 +XPAR 0.0E+00 +ISOP 0.5E+00 +API 0.0E+00 +LIM 0.0E+00 +API1 0.0E+00 +API2 0.0E+00 +LIM1 0.0E+00 +LIM2 0.0E+00 +ISOPRD 0.0E+00 +ISOPP 0.0E+00 +ISOPN 0.0E+00 +ISOPO2 0.0E+00 +DMS 0.0E+00 +MSA 0.0E+00 +DMSO 0.0E+00 +DMSO2 0.0E+00 +CH3SO2H 0.0E+00 +CH3SCH2OO 0.0E+00 +CH3SO2 0.0E+00 +CH3SO3 0.0E+00 +CH3SO2OO 0.0E+00 diff --git a/scenarios/7_drydep/gas_data.dat b/scenarios/7_drydep/gas_data.dat new file mode 100644 index 000000000..0345bcb1f --- /dev/null +++ b/scenarios/7_drydep/gas_data.dat @@ -0,0 +1,78 @@ +# list of gas species +H2SO4 +HNO3 +HCl +NH3 +NO +NO2 +NO3 +N2O5 +HONO +HNO4 +O3 +O1D +O3P +OH +HO2 +H2O2 +CO +SO2 +CH4 +C2H6 +CH3O2 +ETHP +HCHO +CH3OH +ANOL +CH3OOH +ETHOOH +ALD2 +HCOOH +RCOOH +C2O3 +PAN +ARO1 +ARO2 +ALK1 +OLE1 +API1 +API2 +LIM1 +LIM2 +PAR +AONE +MGLY +ETH +OLET +OLEI +TOL +XYL +CRES +TO2 +CRO +OPEN +ONIT +ROOH +RO2 +ANO2 +NAP +XO2 +XPAR +ISOP +ISOPRD +ISOPP +ISOPN +ISOPO2 +API +LIM +DMS +MSA +DMSO +DMSO2 +CH3SO2H +CH3SCH2OO +CH3SO2 +CH3SO3 +CH3SO2OO +CH3SO2CH2OO +SULFHOX diff --git a/scenarios/7_drydep/gas_emit.dat b/scenarios/7_drydep/gas_emit.dat new file mode 100644 index 000000000..1473855a4 --- /dev/null +++ b/scenarios/7_drydep/gas_emit.dat @@ -0,0 +1,26 @@ +# time (s) +# rate = scaling parameter +# emissions (mol m^{-2} s^{-1}) +time 0 3600 7200 10800 14400 18000 21600 25200 28800 32400 36000 39600 43200 46800 50400 54000 57600 61200 64800 68400 72000 75600 79200 82800 90000 93600 97200 100800 104400 108000 +rate 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +#rate 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0 0 0 0 0 0 0 0 0 0 0 0 +#rate 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +SO2 4.234E-09 5.481E-09 5.089E-09 5.199E-09 5.221E-09 5.284E-09 5.244E-09 5.280E-09 5.560E-09 5.343E-09 4.480E-09 3.858E-09 3.823E-09 3.607E-09 3.533E-09 3.438E-09 2.866E-09 2.667E-09 2.636E-09 2.573E-09 2.558E-09 2.573E-09 2.715E-09 3.170E-09 4.234E-09 5.481E-09 5.089E-09 5.199E-09 5.221E-09 5.284E-09 +NO2 3.024E-09 3.334E-09 3.063E-09 3.281E-09 3.372E-09 3.523E-09 3.402E-09 3.551E-09 3.413E-09 3.985E-09 3.308E-09 2.933E-09 2.380E-09 1.935E-09 1.798E-09 1.537E-09 9.633E-10 8.873E-10 7.968E-10 6.156E-10 5.920E-10 6.320E-10 9.871E-10 1.901E-09 .024E-09 3.334E-09 3.063E-09 3.281E-09 3.372E-09 3.523E-09 +NO 5.749E-08 6.338E-08 5.825E-08 6.237E-08 6.411E-08 6.699E-08 6.468E-08 6.753E-08 6.488E-08 7.575E-08 6.291E-08 5.576E-08 4.524E-08 3.679E-08 3.419E-08 2.924E-08 1.832E-08 1.687E-08 1.515E-08 1.171E-08 1.125E-08 1.202E-08 1.877E-08 3.615E-08 5.749E-08 6.338E-08 5.825E-08 6.237E-08 6.411E-08 6.699E-08 +#NH3 1.786E-08 1.741E-08 3.278E-08 2.932E-08 3.281E-08 3.761E-08 3.300E-08 3.609E-08 2.694E-08 1.349E-08 1.083E-08 5.106E-09 4.174E-09 4.577E-09 5.453E-09 5.476E-09 1.992E-09 5.414E-09 1.968E-09 1.935E-09 1.981E-09 2.069E-09 2.165E-09 5.493E-09 +NH3 8.93E-09 8.705E-09 1.639E-08 1.466E-08 1.6405E-08 1.8805E-08 1.65E-08 1.8045E-08 1.347E-08 6.745E-09 5.415E-09 2.553E-09 2.087E-09 2.2885E-09 2.7265E-09 2.738E-09 9.96E-10 2.707E-09 9.84E-10 9.675E-10 9.905E-10 1.0345E-09 1.0825E-09 2.7465E-09 8.93E-09 8.705E-09 1.639E-08 1.466E-08 1.6405E-08 1.8805E-08 +CO 7.839E-07 5.837E-07 4.154E-07 4.458E-07 4.657E-07 4.912E-07 4.651E-07 4.907E-07 6.938E-07 8.850E-07 8.135E-07 4.573E-07 3.349E-07 2.437E-07 2.148E-07 1.662E-07 8.037E-08 7.841E-08 6.411E-08 2.551E-08 2.056E-08 3.058E-08 1.083E-07 3.938E-07 7.839E-07 5.837E-07 4.154E-07 4.458E-07 4.657E-07 4.912E-07 +ALD2 1.702E-09 1.283E-09 9.397E-10 1.024E-09 1.076E-09 1.132E-09 1.068E-09 1.130E-09 1.651E-09 2.132E-09 1.985E-09 1.081E-09 7.847E-10 5.676E-10 5.003E-10 3.838E-10 1.784E-10 1.766E-10 1.430E-10 5.173E-11 4.028E-11 6.349E-11 2.428E-10 8.716E-10 1.702E-09 1.283E-09 9.397E-10 1.024E-09 1.076E-09 1.132E-09 +HCHO 4.061E-09 3.225E-09 2.440E-09 2.639E-09 2.754E-09 2.888E-09 2.741E-09 2.885E-09 4.088E-09 5.186E-09 4.702E-09 2.601E-09 1.923E-09 1.412E-09 1.252E-09 9.776E-10 4.687E-10 4.657E-10 3.836E-10 1.717E-10 1.448E-10 1.976E-10 6.193E-10 2.090E-09 4.061E-09 3.225E-09 2.440E-09 2.639E-09 2.754E-09 2.888E-09 +ETH 1.849E-08 1.391E-08 1.010E-08 1.095E-08 1.148E-08 1.209E-08 1.142E-08 1.205E-08 1.806E-08 2.320E-08 2.149E-08 1.146E-08 8.384E-09 6.124E-09 5.414E-09 4.119E-09 1.953E-09 1.927E-09 1.575E-09 6.164E-10 4.973E-10 7.420E-10 2.653E-09 9.477E-09 1.849E-08 1.391E-08 1.010E-08 1.095E-08 1.148E-08 1.209E-08 +OLEI 5.948E-09 4.573E-09 3.374E-09 3.668E-09 3.851E-09 4.050E-09 3.841E-09 4.052E-09 6.094E-09 7.795E-09 7.215E-09 3.738E-09 2.718E-09 1.973E-09 1.729E-09 1.338E-09 6.333E-10 6.394E-10 5.126E-10 2.089E-10 1.708E-10 2.480E-10 8.947E-10 3.057E-09 5.948E-09 4.573E-09 3.374E-09 3.668E-09 3.851E-09 4.050E-09 +OLET 5.948E-09 4.573E-09 3.374E-09 3.668E-09 3.851E-09 4.050E-09 3.841E-09 4.052E-09 6.094E-09 7.795E-09 7.215E-09 3.738E-09 2.718E-09 1.973E-09 1.729E-09 1.338E-09 6.333E-10 6.394E-10 5.126E-10 2.089E-10 1.708E-10 2.480E-10 8.947E-10 3.057E-09 5.948E-09 4.573E-09 3.374E-09 3.668E-09 3.851E-09 4.050E-09 +TOL 6.101E-09 8.706E-09 7.755E-09 8.024E-09 8.202E-09 8.410E-09 8.218E-09 8.407E-09 1.020E-08 1.139E-08 7.338E-09 4.184E-09 3.078E-09 2.283E-09 2.010E-09 1.575E-09 8.966E-10 6.705E-10 5.395E-10 2.462E-10 2.106E-10 2.852E-10 9.300E-10 3.144E-09 6.101E-09 8.706E-09 7.755E-09 8.024E-09 8.202E-09 8.410E-09 +XYL 5.599E-09 4.774E-09 3.660E-09 3.909E-09 4.060E-09 4.239E-09 4.060E-09 4.257E-09 6.036E-09 7.448E-09 6.452E-09 3.435E-09 2.525E-09 1.859E-09 1.650E-09 1.302E-09 6.852E-10 6.773E-10 5.437E-10 2.697E-10 2.358E-10 3.059E-10 8.552E-10 2.861E-10 5.599E-09 4.774E-09 3.660E-09 3.909E-09 4.060E-09 4.239E-09 +AONE 7.825E-10 2.858E-09 2.938E-09 2.947E-09 2.948E-09 2.951E-09 2.947E-09 2.954E-09 3.032E-09 2.766E-09 1.313E-09 1.015E-09 8.363E-10 7.040E-10 6.404E-10 6.264E-10 5.661E-10 1.538E-10 1.500E-10 1.395E-10 1.476E-10 1.503E-10 2.256E-10 4.244E-10 7.825E-10 2.858E-09 2.938E-09 2.947E-09 2.948E-09 2.951E-09 +PAR 1.709E-07 1.953E-07 1.698E-07 1.761E-07 1.808E-07 1.865E-07 1.822E-07 1.859E-07 2.412E-07 2.728E-07 2.174E-07 1.243E-07 9.741E-08 7.744E-08 6.931E-08 5.805E-08 3.900E-08 3.317E-08 2.956E-08 2.306E-08 2.231E-08 2.395E-08 4.284E-08 9.655E-08 1.709E-07 1.953E-07 1.698E-07 1.761E-07 1.808E-07 1.865E-07 +ISOP 2.412E-10 2.814E-10 3.147E-10 4.358E-10 5.907E-10 6.766E-10 6.594E-10 5.879E-10 5.435E-10 6.402E-10 5.097E-10 9.990E-11 7.691E-11 5.939E-11 5.198E-11 4.498E-11 3.358E-11 2.946E-11 2.728E-11 2.183E-11 1.953E-11 1.890E-11 2.948E-11 1.635E-10 2.412E-10 2.814E-10 3.147E-10 4.358E-10 5.907E-10 6.766E-10 +CH3OH 2.368E-10 6.107E-10 6.890E-10 6.890E-10 6.890E-10 6.889E-10 6.886E-10 6.890E-10 6.890E-10 5.414E-10 3.701E-10 2.554E-10 1.423E-10 6.699E-11 2.912E-11 2.877E-11 2.825E-11 2.056E-12 2.056E-12 2.056E-12 2.435E-12 2.435E-12 4.030E-11 1.168E-10 2.368E-10 6.107E-10 6.890E-10 6.890E-10 6.890E-10 6.889E-10 +ANOL 5.304E-09 7.960E-09 7.649E-09 7.649E-09 7.432E-09 7.428E-09 7.431E-09 7.434E-09 7.434E-09 6.979E-09 5.666E-09 4.361E-09 4.148E-09 3.289E-09 2.858E-09 2.856E-09 1.127E-09 9.615E-10 9.616E-10 9.616E-10 9.654E-10 9.654E-10 1.397E-09 2.264E-09 5.304E-09 7.960E-09 7.649E-09 7.649E-09 7.432E-09 7.428E-09 + diff --git a/scenarios/7_drydep/height.dat b/scenarios/7_drydep/height.dat new file mode 100644 index 000000000..efcf74e95 --- /dev/null +++ b/scenarios/7_drydep/height.dat @@ -0,0 +1,5 @@ +# time (s) +# height (m) +time 0 3600 7200 10800 14400 18000 21600 25200 28800 32400 36000 39600 43200 46800 50400 54000 57600 61200 64800 68400 72000 75600 79200 82800 86400 90000 93600 97200 100800 104400 108000 +height 171.045 228.210 296.987 366.002 410.868 414.272 417.807 414.138 397.465 376.864 364.257 352.119 338.660 322.028 305.246 258.497 240.478 187.229 145.851 128.072 110.679 97.628 93.034 93.034 93.034 93.034 93.034 93.034 93.034 93.034 93.034 + diff --git a/scenarios/7_drydep/pres.dat b/scenarios/7_drydep/pres.dat new file mode 100644 index 000000000..07215617a --- /dev/null +++ b/scenarios/7_drydep/pres.dat @@ -0,0 +1,4 @@ +# time (s) +# pressure (Pa) +time 0 +pressure 1e5 diff --git a/scenarios/7_drydep/temp.dat b/scenarios/7_drydep/temp.dat new file mode 100644 index 000000000..791422ba5 --- /dev/null +++ b/scenarios/7_drydep/temp.dat @@ -0,0 +1,4 @@ +# time (s) +# temp (K) +time 0 3600 7200 10800 14400 18000 21600 25200 28800 32400 36000 39600 43200 46800 50400 54000 57600 61200 64800 68400 72000 75600 79200 82800 86400 90000 93600 97200 100800 104400 108000 111600 115200 118800 122400 126000 129600 133200 136800 140400 144000 147600 151200 154800 158400 162000 165600 169200 172800 +temp 290.016 292.5 294.5 296.112 297.649 299.049 299.684 299.509 299.002 298.432 296.943 295.153 293.475 292.466 291.972 291.96 291.512 291.481 290.5 290.313 290.317 290.362 290.245 290.228 291.466 292.5 294.5 296.112 297.649 299.049 299.684 299.509 299.002 298.432 296.943 295.153 293.475 292.466 291.972 291.96 291.512 291.481 290.5 290.313 290.317 290.362 290.245 290.228 291.466 From 02c7e11ef94c5862d26b87b7f4e8ea6f6063085d Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Sat, 20 Sep 2025 08:54:14 -0500 Subject: [PATCH 24/28] update doc for spec_file_read_drydep_params --- src/scenario.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/scenario.F90 b/src/scenario.F90 index 0934af908..7b59a6017 100644 --- a/src/scenario.F90 +++ b/src/scenario.F90 @@ -1392,7 +1392,6 @@ subroutine spec_file_read_drydep_params(file, drydep_params) !> \page input_format_chamber Input File Format: Dry Deposition Parameters !! !! Dry deposition is simulatied using the specified parameters: - !! - \b drydep_param (integer): flag for the dry deposition parameterization !! - \b z_ref (real, unit m): the reference height \f$z_{\rm ref}\f$ used !! in the calculation of aerodynamic resistance \f$R_a\f$ !! - \b u_mean (real, unit m s^{-1}): the wind speed at the reference From 1c3610c2fb5e59072c1b640ecdc5d50cc7a2a876 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Wed, 8 Oct 2025 14:38:14 -0500 Subject: [PATCH 25/28] fix aero_dist naming convention --- src/aero_dist.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aero_dist.F90 b/src/aero_dist.F90 index 552a868d5..2a485a8b5 100644 --- a/src/aero_dist.F90 +++ b/src/aero_dist.F90 @@ -303,7 +303,7 @@ subroutine aero_dist_output_netcdf(aero_dist, ncid) call pmc_nc_write_real_1d(ncid, num_concs, "num_conc_modes", dim_name="num_modes", & long_name="number concentration for each mode", & unit="m^-3") - call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "num_modes", & + call pmc_nc_write_integer(ncid, aero_dist_n_mode(aero_dist), "n_mode", & description="total number of modes") deallocate(char_rads, log10_std_dev_rads, num_concs) @@ -324,7 +324,7 @@ subroutine aero_dist_input_netcdf(aero_dist, ncid) type(aero_mode_t), allocatable :: modes(:) real(kind=dp), allocatable :: char_rads(:), log10_std_dev_rads(:), num_concs(:) - call pmc_nc_read_integer(ncid, n_mode, "num_modes") + call pmc_nc_read_integer(ncid, n_mode, "n_mode") if (allocated(modes)) deallocate(modes) allocate(modes(n_mode)) From eeb1c472a13c014bc24c6599069f9a2428aea9ef Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Wed, 8 Oct 2025 14:38:48 -0500 Subject: [PATCH 26/28] drydep_params.dat comments update --- scenarios/7_drydep/drydep_params.dat | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/scenarios/7_drydep/drydep_params.dat b/scenarios/7_drydep/drydep_params.dat index b97f6c4ce..5765ddeb8 100644 --- a/scenarios/7_drydep/drydep_params.dat +++ b/scenarios/7_drydep/drydep_params.dat @@ -2,11 +2,11 @@ z_ref 10.0 # reference height (m) u_mean 2.0 # mean wind speed at the reference height (m) z_rough 3.8 # surface roughness length (m) A 0.0024 # characteristic radius of collectors (m) -alpha 5.0 -eps_0 6.0 -gamma .756 -C_B 0.82 -C_IN 92.5 -C_IM 10.4 -nu 11.8 -beta 12.7 \ No newline at end of file +alpha 5.0 # parameter used in impaction efficiency +eps_0 6.0 # emprical constant used in surface resistance +gamma .756 # exponent for Brownian diffusion collection efficiency +C_B 0.82 # coefficient for Brownian diffusion +C_IN 92.5 # coefficient for interception +C_IM 10.4 # coefficient for impaction +nu 11.8 # exponent for interception collection efficiency +beta 12.7 # exponent for impaction collection efficiency From e7e23ea7adcf3e98f3ef0bb83080858f29b84b98 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Wed, 8 Oct 2025 14:44:49 -0500 Subject: [PATCH 27/28] text file clean up --- scenarios/7_drydep/aero_init_dist.dat | 4 ++-- scenarios/7_drydep/drydep_modal.spec | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/scenarios/7_drydep/aero_init_dist.dat b/scenarios/7_drydep/aero_init_dist.dat index 721acd06a..09d97302b 100644 --- a/scenarios/7_drydep/aero_init_dist.dat +++ b/scenarios/7_drydep/aero_init_dist.dat @@ -3,5 +3,5 @@ mass_frac aero_init_comp.dat # composition proportions of species diam_type geometric # type of diameter specified mode_type log_normal # type of distribution num_conc 3.2e9 # particle number concentration (#/m^3) -geom_mean_diam .00000001000000000000 -log10_geom_std_dev .04139268515822504074 \ No newline at end of file +geom_mean_diam 1e-8 # geometric mean diameter (m) +log10_geom_std_dev 0.041 # log_10 of geometric std dev of diameter \ No newline at end of file diff --git a/scenarios/7_drydep/drydep_modal.spec b/scenarios/7_drydep/drydep_modal.spec index 193865999..062458c1a 100644 --- a/scenarios/7_drydep/drydep_modal.spec +++ b/scenarios/7_drydep/drydep_modal.spec @@ -1,5 +1,5 @@ run_type modal # modal run -output_prefix out/drydep_modal # prefix for output files +output_prefix data/modal_dpg_00001000000000000000_sig_2_5_emerson_grass t_max 28800 # total simulation time (s) del_t 60 # timestep (s) @@ -7,8 +7,8 @@ t_output 3600 # output interval (0 disables) (s) t_progress 0 # progress printing interval (0 disables) (s) n_bin 1000 # number of bins (for processing purposes) -d_min .00000040000000000000e-1 -d_max .00025000000000000000e1 +d_min 4e-8 # minimum diameter (m) +d_max 2.5e-3 # maximum diameter (m) gas_data gas_data.dat # file containing gas data aerosol_data aero_data.dat # file containing aerosol data @@ -23,7 +23,7 @@ gas_background gas_back.dat # background gas concentrations file aero_emissions aero_emit.dat # aerosol emissions file aero_background aero_back.dat # aerosol background file loss_function drydep # loss function specification -drydep_params drydep_params.dat # dry deposition parameters +drydep_params drydep_grass_emerson.dat rel_humidity 0.95 # initial relative humidity (1) latitude 0 # latitude (degrees, -90 to 90) From bd9f931865947f75c8c8a0a507cfbf8044588852 Mon Sep 17 00:00:00 2001 From: zdaq12 Date: Wed, 8 Oct 2025 14:51:22 -0500 Subject: [PATCH 28/28] newlines at end of files --- CMakeLists.txt | 2 +- scenarios/7_drydep/1_run_modal.sh | 4 ++-- scenarios/7_drydep/2_process_modal.sh | 2 +- scenarios/7_drydep/aero_back.dat | 2 +- scenarios/7_drydep/aero_back_dist.dat | 2 +- scenarios/7_drydep/aero_emit_dist.dat | 2 +- scenarios/7_drydep/aero_init_dist.dat | 2 +- scenarios/7_drydep/drydep_modal.spec | 2 +- scenarios/7_drydep/drydep_modal_process.F90 | 6 +++--- src/run_modal.F90 | 2 +- 10 files changed, 13 insertions(+), 13 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1a0c54df7..9f5f165cc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -439,4 +439,4 @@ add_executable(drydep_modal_process scenarios/7_drydep/drydep_modal_process.F90) target_link_libraries(drydep_modal_process partmclib) -###################################################################### \ No newline at end of file +###################################################################### diff --git a/scenarios/7_drydep/1_run_modal.sh b/scenarios/7_drydep/1_run_modal.sh index 9354c9fc5..f6ead5b7d 100755 --- a/scenarios/7_drydep/1_run_modal.sh +++ b/scenarios/7_drydep/1_run_modal.sh @@ -5,8 +5,8 @@ set -e # turn on command echoing set -v -mkdir -p out +mkdir -p data ../../build/partmc drydep_modal.spec -# Now run ./2_process_drydep_modal.sh to process the data \ No newline at end of file +# Now run ./2_process_drydep_modal.sh to process the data diff --git a/scenarios/7_drydep/2_process_modal.sh b/scenarios/7_drydep/2_process_modal.sh index 6c92ff47e..1cebe6af5 100755 --- a/scenarios/7_drydep/2_process_modal.sh +++ b/scenarios/7_drydep/2_process_modal.sh @@ -7,4 +7,4 @@ set -v # The data should have already been generated by ./1_run_dry_dep_modal.sh -../../build/dry_dep_modal_process \ No newline at end of file +../../build/drydep_modal_process diff --git a/scenarios/7_drydep/aero_back.dat b/scenarios/7_drydep/aero_back.dat index 72f0e739a..3879b8f6b 100644 --- a/scenarios/7_drydep/aero_back.dat +++ b/scenarios/7_drydep/aero_back.dat @@ -3,4 +3,4 @@ # aerosol distribution filename time 0 rate 0 -dist aero_back_dist.dat \ No newline at end of file +dist aero_back_dist.dat diff --git a/scenarios/7_drydep/aero_back_dist.dat b/scenarios/7_drydep/aero_back_dist.dat index ee6d35bde..ee194a062 100644 --- a/scenarios/7_drydep/aero_back_dist.dat +++ b/scenarios/7_drydep/aero_back_dist.dat @@ -1 +1 @@ -# no background distribution \ No newline at end of file +# no background distribution diff --git a/scenarios/7_drydep/aero_emit_dist.dat b/scenarios/7_drydep/aero_emit_dist.dat index 57995ab14..8d50b7844 100644 --- a/scenarios/7_drydep/aero_emit_dist.dat +++ b/scenarios/7_drydep/aero_emit_dist.dat @@ -1 +1 @@ -# no emissions \ No newline at end of file +# no emissions diff --git a/scenarios/7_drydep/aero_init_dist.dat b/scenarios/7_drydep/aero_init_dist.dat index 09d97302b..7f510934c 100644 --- a/scenarios/7_drydep/aero_init_dist.dat +++ b/scenarios/7_drydep/aero_init_dist.dat @@ -4,4 +4,4 @@ diam_type geometric # type of diameter specified mode_type log_normal # type of distribution num_conc 3.2e9 # particle number concentration (#/m^3) geom_mean_diam 1e-8 # geometric mean diameter (m) -log10_geom_std_dev 0.041 # log_10 of geometric std dev of diameter \ No newline at end of file +log10_geom_std_dev 0.041 # log_10 of geometric std dev of diameter diff --git a/scenarios/7_drydep/drydep_modal.spec b/scenarios/7_drydep/drydep_modal.spec index 062458c1a..679ad83dc 100644 --- a/scenarios/7_drydep/drydep_modal.spec +++ b/scenarios/7_drydep/drydep_modal.spec @@ -30,4 +30,4 @@ latitude 0 # latitude (degrees, -90 to 90) longitude 0 # longitude (degrees, -180 to 180) altitude 0 # altitude (m) start_time 21600 # start time (s since 00:00 UTC) -start_day 200 # start day of year (UTC) \ No newline at end of file +start_day 200 # start day of year (UTC) diff --git a/scenarios/7_drydep/drydep_modal_process.F90 b/scenarios/7_drydep/drydep_modal_process.F90 index 1648542e3..72b409e1f 100644 --- a/scenarios/7_drydep/drydep_modal_process.F90 +++ b/scenarios/7_drydep/drydep_modal_process.F90 @@ -8,7 +8,7 @@ program process use pmc_scenario character(len=PMC_MAX_FILENAME_LEN), parameter :: prefix & -= "out/modal_test" += "data/modal_dpg_00001000000000000000_sig_2_5_emerson_grass" character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename type(bin_grid_t) :: bin_grid @@ -93,7 +93,7 @@ program process call make_filename(out_filename, prefix, "_process.nc") write(*,*) "Writing " // trim(out_filename) call pmc_nc_open_write(out_filename, ncid) - call pmc_nc_write_info(ncid, uuid, "1_urban_plume modal process") + call pmc_nc_write_info(ncid, uuid, "7_drydep modal process") call pmc_nc_write_real_1d(ncid, times, "time", dim_name="time", unit="s") call stats_1d_output_netcdf(stats_tot_num_conc, ncid, "tot_num_conc", & dim_name="time", unit="m^{-3}") @@ -103,4 +103,4 @@ program process call pmc_mpi_finalize() -end program process \ No newline at end of file +end program process diff --git a/src/run_modal.F90 b/src/run_modal.F90 index ee22ebd16..06d6a0d59 100644 --- a/src/run_modal.F90 +++ b/src/run_modal.F90 @@ -137,4 +137,4 @@ end subroutine run_modal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module pmc_run_modal \ No newline at end of file +end module pmc_run_modal