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)
75 integer ima(bin_grid_size(bin_grid),bin_grid_size(bin_grid))
76 real(kind=
dp) g(bin_grid_size(bin_grid)), r(bin_grid_size(bin_grid))
77 real(kind=
dp) e(bin_grid_size(bin_grid))
78 real(kind=
dp) k_bin(bin_grid_size(bin_grid),bin_grid_size(bin_grid))
79 real(kind=
dp) ck(bin_grid_size(bin_grid),bin_grid_size(bin_grid))
80 real(kind=
dp) ec(bin_grid_size(bin_grid),bin_grid_size(bin_grid))
81 real(kind=
dp) taug(bin_grid_size(bin_grid)), taup(bin_grid_size(bin_grid))
82 real(kind=
dp) taul(bin_grid_size(bin_grid)), tauu(bin_grid_size(bin_grid))
83 real(kind=
dp) prod(bin_grid_size(bin_grid)), ploss(bin_grid_size(bin_grid))
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')
113 do i = 1,bin_grid_size(bin_grid)
114 r(i) = bin_grid%centers(i) * 1d6
116 * aero_data%density(1) * 1d6
123 call courant(bin_grid_size(bin_grid), bin_grid%widths(1), e, ima, c)
126 last_progress_time = 0d0
131 call bin_kernel(bin_grid_size(bin_grid), bin_grid%centers, aero_data, &
132 run_sect_opt%coag_kernel_type, env_state, k_bin)
134 do i = 1,bin_grid_size(bin_grid)
135 do j = 1,bin_grid_size(bin_grid)
136 ck(i,j) = ck(i,j) * 1d6
141 do i = 1,bin_grid_size(bin_grid)
142 do j = 1,bin_grid_size(bin_grid)
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)
162 call coad(bin_grid_size(bin_grid), run_sect_opt%del_t, taug, taup, &
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)