34 real(kind=dp) :: t_max
36 real(kind=dp) :: del_t
38 real(kind=dp) :: t_output
40 real(kind=dp) :: t_progress
42 logical :: do_coagulation
44 character(len=300) :: prefix
46 integer :: coag_kernel_type
48 character(len=PMC_UUID_LEN) :: uuid
56 subroutine run_sect(bin_grid, gas_data, aero_data, aero_dist, &
57 scenario, env_state, run_sect_opt)
84 real(kind=dp) time, last_output_time, last_progress_time
89 integer i, j, i_time, num_t, i_summary
90 logical do_output, do_progress
93 "del_t", run_sect_opt%del_t)
95 "del_t", run_sect_opt%del_t)
97 "del_t", run_sect_opt%del_t)
106 'run_sect() can only use one aerosol species')
114 r(i) = bin_grid%centers(i) * 1d6
116 * aero_data%density(1) * 1d6
126 last_progress_time = 0d0
132 run_sect_opt%coag_kernel_type, env_state, k_bin)
136 ck(i,j) = ck(i,j) * 1d6
143 ck(i,j) = ck(i,j) * run_sect_opt%del_t * bin_grid%widths(i)
148 call check_event(time, run_sect_opt%del_t, run_sect_opt%t_output, &
149 last_output_time, do_output)
152 aero_binned, gas_data, gas_state, env_state, i_summary, &
153 time, run_sect_opt%t_output, run_sect_opt%uuid)
157 num_t = nint(run_sect_opt%t_max / run_sect_opt%del_t)
160 if (run_sect_opt%do_coagulation)
then 161 g = aero_binned%vol_conc(:,1) * aero_data%density(1)
163 taul, tauu, prod, ploss, c, ima, g, r, e, ck, ec)
164 aero_binned%vol_conc(:,1) = g / aero_data%density(1)
165 aero_binned%num_conc = aero_binned%vol_conc(:,1) &
169 time = run_sect_opt%t_max *
real(i_time, kind=dp) &
170 /
real(num_t, kind=
dp)
172 old_env_state = env_state
175 env_state, old_env_state, gas_data, gas_state)
177 env_state, old_env_state, bin_grid, aero_data, aero_binned)
180 call check_event(time, run_sect_opt%del_t, run_sect_opt%t_output, &
181 last_output_time, do_output)
183 i_summary = i_summary + 1
185 aero_binned, gas_data, gas_state, env_state, i_summary, &
186 time, run_sect_opt%t_output, run_sect_opt%uuid)
190 call check_event(time, run_sect_opt%del_t, run_sect_opt%t_progress, &
191 last_progress_time, do_progress)
192 if (do_progress)
then 193 write(*,
'(a6,a8)')
'step',
'time' 194 write(*,
'(i6,f8.1)') i_time, time
203 subroutine coad(n_bin, dt, taug, taup, taul, tauu, prod, ploss, &
204 c, ima, g, r, e, ck, ec)
208 real(kind=dp) taug(n_bin)
209 real(kind=dp) taup(n_bin)
210 real(kind=dp) taul(n_bin)
211 real(kind=dp) tauu(n_bin)
212 real(kind=dp) prod(n_bin)
213 real(kind=dp) ploss(n_bin)
214 real(kind=dp) c(n_bin,n_bin)
215 integer ima(n_bin,n_bin)
216 real(kind=dp) g(n_bin)
217 real(kind=dp) r(n_bin)
218 real(kind=dp) e(n_bin)
219 real(kind=dp) ck(n_bin,n_bin)
220 real(kind=dp) ec(n_bin,n_bin)
222 real(kind=dp),
parameter :: gmin = 1d-60
224 integer i, i0, i1, j, k, kp
225 real(kind=dp) x0, gsi, gsj, gsk, gk, x1, flux
233 do i0 = 1,(n_bin - 1)
234 if (g(i0) .gt. gmin)
goto 2000
237 do i1 = (n_bin - 1),1,-1
238 if (g(i1) .gt. gmin)
goto 2010
247 x0 = ck(i,j) * g(i) * g(j)
248 x0 = min(x0, g(i) * e(j))
250 if (j .ne. k) x0 = min(x0, g(j) * e(i))
256 ploss(i) = ploss(i) + gsi
257 ploss(j) = ploss(j) + gsj
264 if (gk .gt. gmin)
then 265 x1 = log(g(kp) / gk + 1d-60)
266 flux = gsk / x1 * (exp(0.5d0 * x1) &
267 - exp(x1 * (0.5d0 - c(i,j))))
272 prod(k) = prod(k) + gsk - flux
273 prod(kp) = prod(kp) + flux
284 subroutine courant(n_bin, log_width, e, ima, c)
287 integer,
intent(in) :: n_bin
289 real(kind=dp),
intent(in) :: log_width
291 real(kind=dp),
intent(in) :: e(n_bin)
293 integer,
intent(out) :: ima(n_bin,n_bin)
295 real(kind=dp),
intent(out) :: c(n_bin,n_bin)
315 if (c(i,j) .lt. 1d0 - 1d-08)
then 317 c(i,j) = log(x0 / e(k-1)) / (3d0 * log_width)
322 ima(i,j) = min(n_bin - 1, kk)
337 integer,
intent(in) :: n_bin
339 real(kind=dp),
intent(in) :: k(n_bin,n_bin)
341 real(kind=dp),
intent(out) :: k_smooth(n_bin,n_bin)
343 integer i, j, im, ip, jm, jp
349 jp = min0(j + 1, n_bin)
350 ip = min0(i + 1, n_bin)
351 k_smooth(i,j) = 0.125d0 * (k(i,jm) + k(im,j) &
352 + k(ip,j) + k(i,jp)) &
355 k_smooth(i,j) = 0.5d0 * k_smooth(i,j)
subroutine check_time_multiple(first_name, first_time, second_name, second_time)
Check that the first time interval is close to an integer multiple of the second, and warn if it is n...
subroutine scenario_update_aero_binned(scenario, delta_t, env_state, old_env_state, bin_grid, aero_data, aero_binned)
Do emissions and background dilution from the environment for a binned aerosol distribution.
subroutine coad(n_bin, dt, taug, taup, taul, tauu, prod, ploss, c, ima, g, r, e, ck, ec)
Collision subroutine, exponential approach.
The bin_grid_t structure and associated subroutines.
real(kind=dp) elemental function aero_data_rad2vol(aero_data, r)
Convert geometric radius (m) to mass-equivalent volume (m^3).
subroutine aero_binned_add_aero_dist(aero_binned, bin_grid, aero_data, aero_dist)
Add an aero_dist_t to an aero_binned_t.
subroutine run_sect(bin_grid, gas_data, aero_data, aero_dist, scenario, env_state, run_sect_opt)
Run a sectional simulation.
The env_state_t structure and associated subroutines.
subroutine scenario_update_gas_state(scenario, delta_t, env_state, old_env_state, gas_data, gas_state)
Do gas emissions and background dilution.
The aero_dist_t structure and associated subroutines.
Generic coagulation kernel.
The gas_data_t structure and associated subroutines.
subroutine scenario_update_env_state(scenario, env_state, time)
Update time-dependent contents of the environment. scenario_init_env_state() should have been called ...
subroutine output_sectional(prefix, bin_grid, aero_data, aero_binned, gas_data, gas_state, env_state, index, time, del_t, uuid)
Write the current sectional data.
A complete aerosol distribution, consisting of several modes.
Current environment state.
subroutine check_event(time, timestep, interval, last_time, do_event)
Computes whether an event is scheduled to take place.
subroutine die_msg(code, error_msg)
Error immediately.
subroutine courant(n_bin, log_width, e, ima, c)
Determines the Courant number for each bin pair.
Write data in NetCDF format.
elemental integer function bin_grid_size(bin_grid)
Return the number of bins in the grid, or -1 if the bin grid is not allocated.
Options controlling the operation of run_sect().
elemental integer function gas_data_n_spec(gas_data)
Return the number of gas species.
The aero_data_t structure and associated subroutines.
subroutine gas_state_set_size(gas_state, n_spec)
Sets the sizes of the gas state.
1D grid, either logarithmic or linear.
The gas_state_t structure and associated subroutines.
integer function find_1d(n, x_vals, x)
Find the position of a real number in an arbitrary 1D array.
integer, parameter dp
Kind of a double precision real number.
Current state of the gas mixing ratios in the system.
elemental integer function aero_data_n_spec(aero_data)
Return the number of aerosol species, or -1 if uninitialized.
Aerosol material properties and associated data.
Common utility subroutines.
The aero_binned_t structure and associated subroutines.
The scenario_t structure and associated subroutines.
Aerosol number and volume distributions stored per bin.
subroutine bin_kernel(n_bin, bin_r, aero_data, coag_kernel_type, env_state, k)
Computes an array of kernel values for each bin pair. k(i,j) is the kernel value at the centers of bi...
subroutine smooth_bin_kernel(n_bin, k, k_smooth)
Smooths kernel values for bin pairs, and halves the self-rate.