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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 788097ae3..9f5f165cc 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 (dry dep)" 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 @@ -233,7 +253,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} @@ -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_drydep/drydep_modal_process + +add_executable(drydep_modal_process + scenarios/7_drydep/drydep_modal_process.F90) +target_link_libraries(drydep_modal_process partmclib) + +###################################################################### diff --git a/scenarios/7_drydep/1_run_modal.sh b/scenarios/7_drydep/1_run_modal.sh new file mode 100755 index 000000000..f6ead5b7d --- /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 data + +../../build/partmc drydep_modal.spec + +# 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 new file mode 100755 index 000000000..1cebe6af5 --- /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/drydep_modal_process diff --git a/scenarios/7_drydep/aero_back.dat b/scenarios/7_drydep/aero_back.dat new file mode 100644 index 000000000..3879b8f6b --- /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 diff --git a/scenarios/7_drydep/aero_back_dist.dat b/scenarios/7_drydep/aero_back_dist.dat new file mode 100644 index 000000000..ee194a062 --- /dev/null +++ b/scenarios/7_drydep/aero_back_dist.dat @@ -0,0 +1 @@ +# no background distribution 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..8d50b7844 --- /dev/null +++ b/scenarios/7_drydep/aero_emit_dist.dat @@ -0,0 +1 @@ +# no emissions 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..7f510934c --- /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 1e-8 # geometric mean diameter (m) +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 new file mode 100644 index 000000000..679ad83dc --- /dev/null +++ b/scenarios/7_drydep/drydep_modal.spec @@ -0,0 +1,33 @@ +run_type modal # modal run +output_prefix data/modal_dpg_00001000000000000000_sig_2_5_emerson_grass + +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 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 +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_grass_emerson.dat + +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) diff --git a/scenarios/7_drydep/drydep_modal_process.F90 b/scenarios/7_drydep/drydep_modal_process.F90 new file mode 100644 index 000000000..72b409e1f --- /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 & += "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 + 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_drydep_rates(scenario, aero_dist, 0.0d0, density, & + env_state, rates_0) + call stats_1d_add(stats_rates_0, rates_0) + call scenario_modal_drydep_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, "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}") + 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 diff --git a/scenarios/7_drydep/drydep_params.dat b/scenarios/7_drydep/drydep_params.dat new file mode 100644 index 000000000..5765ddeb8 --- /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 # 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 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 diff --git a/src/aero_dist.F90 b/src/aero_dist.F90 index 78a696918..2a485a8b5 100644 --- a/src/aero_dist.F90 +++ b/src/aero_dist.F90 @@ -271,6 +271,82 @@ 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, 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), "n_mode", & + description="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 + 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, "n_mode") + + if (allocated(modes)) deallocate(modes) + allocate(modes(n_mode)) + + 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 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Read an array of aero_dists with associated times and rates from diff --git a/src/output.F90 b/src/output.F90 index b102f0275..4cba0a9fe 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 @@ -782,6 +783,144 @@ 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) + call bin_grid_output_netcdf(bin_grid, ncid, "diam", unit="m") + + 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 + + end subroutine input_modal + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Read the specification for an output type from a spec file and diff --git a/src/partmc.F90 b/src/partmc.F90 index b873970f2..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 @@ -300,6 +301,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 +785,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, .false., aero_dist_init) + call spec_file_close(sub_file) + + 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) + + ! 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 diff --git a/src/run_modal.F90 b/src/run_modal.F90 new file mode 100644 index 000000000..06d6a0d59 --- /dev/null +++ b/src/run_modal.F90 @@ -0,0 +1,140 @@ +!> \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 + !> 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 + + 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 + + ! 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) + + call check_event(time, run_modal_opt%del_t, run_modal_opt%t_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, & + 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_time, 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 diff --git a/src/scenario.F90 b/src/scenario.F90 index ead729f94..7b59a6017 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 @@ -40,6 +43,34 @@ module pmc_scenario !! 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. @@ -97,6 +128,8 @@ module pmc_scenario !> Type of loss rate function. integer :: loss_function_type + !> Dry deposition parameters + type(drydep_params_t) :: drydep !> Chamber parameters for wall loss and sedimentation. type(chamber_t) :: chamber end type scenario_t @@ -392,6 +425,77 @@ 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) + + !> 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 + + !> 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 (num. conc.) + 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) + + 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_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) + + ! 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. @@ -421,8 +525,8 @@ 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, & - aero_data, env_state) + 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 scenario_loss_rate = chamber_loss_rate_wall(scenario%chamber, vol, & @@ -441,8 +545,8 @@ 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, & - env_state) + real(kind=dp) function scenario_loss_rate_drydep(vol, density, aero_data, & + env_state, scenario) !> Particle volume (m^3). real(kind=dp), intent(in) :: vol @@ -452,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 @@ -462,22 +568,13 @@ 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 - - ! 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 + type(drydep_params_t) :: drydep_params + + drydep_params = scenario%drydep ! particle diameter d_p = aero_data_vol2diam(aero_data, vol) @@ -506,36 +603,307 @@ 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)) ! 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_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_drydep(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) :: 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 + type(drydep_params_t) :: drydep_params + +#ifdef PMC_USE_QUADPACK + ! Do numerical integration when QUADPACK is available + scenario_integrated_loss_rate_drydep = scenario_integrated_loss_rate_drydep_quadpack( & + scenario, aero_mode, moment, density, env_state) + return +#endif + + drydep_params = scenario%drydep + + ! 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 * 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) / & + (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 = drydep_params%C_B * Sc**(-drydep_params%gamma) + + ! Collection efficiency due to interception + 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 * drydep_params%A) + ! Collection efficiency due to impaction + 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 / (drydep_params%eps_0 * u_star * (E_B + E_IN + E_IM) * R1) - end function scenario_loss_rate_dry_dep + ! Integrated deposition velocity + V_d_hat = V_g_hat + (1.0d0 / (R_a + R_s)) + + ! Loss rate + scenario_integrated_loss_rate_drydep = V_d_hat / env_state%height + + end function scenario_integrated_loss_rate_drydep + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#ifdef PMC_USE_QUADPACK + real(kind=dp) function scenario_integrated_loss_rate_drydep_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) :: 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(:) + + type(drydep_params_t) :: drydep_params + + drydep_params = scenario%drydep + + ! 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_drydep_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 * 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) / & + (3.0d0 * const%pi * visc_d * d_p) + Sc = visc_k / diff_p + + ! Collection efficiencies + 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 * 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 / (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)) + + ! 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_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_drydep_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_drydep( & + scenario, aero_dist%mode(i_mode), moment, density, env_state) & + * env_state%height + end do + + end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -785,6 +1153,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 @@ -889,6 +1258,16 @@ 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 + 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 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 +1379,64 @@ end subroutine spec_file_read_scenario !! - \ref input_format_scenario --- the environment data !! containing the mixing layer height profile +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Read dry deposition parameters from a spec file. + 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) :: drydep_params + + !> \page input_format_chamber Input File Format: Dry Deposition Parameters + !! + !! Dry deposition is simulatied using the specified parameters: + !! - \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", 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 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Determines the number of bytes required to pack the given value. @@ -1026,6 +1463,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) @@ -1088,6 +1526,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) @@ -1148,6 +1587,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) @@ -1185,6 +1625,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 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.