1 | !> @file lpm_collision_kernels.f90 |
---|
2 | !------------------------------------------------------------------------------! |
---|
3 | ! This file is part of the PALM model system. |
---|
4 | ! |
---|
5 | ! PALM is free software: you can redistribute it and/or modify it under the |
---|
6 | ! terms of the GNU General Public License as published by the Free Software |
---|
7 | ! Foundation, either version 3 of the License, or (at your option) any later |
---|
8 | ! version. |
---|
9 | ! |
---|
10 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
11 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
12 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
13 | ! |
---|
14 | ! You should have received a copy of the GNU General Public License along with |
---|
15 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
16 | ! |
---|
17 | ! Copyright 1997-2019 Leibniz Universitaet Hannover |
---|
18 | !------------------------------------------------------------------------------! |
---|
19 | ! |
---|
20 | ! Current revisions: |
---|
21 | ! ----------------- |
---|
22 | ! |
---|
23 | ! |
---|
24 | ! Former revisions: |
---|
25 | ! ----------------- |
---|
26 | ! $Id: lpm_collision_kernels.f90 3655 2019-01-07 16:51:22Z suehring $ |
---|
27 | ! Modularization of all bulk cloud physics code components |
---|
28 | ! |
---|
29 | ! 2718 2018-01-02 08:49:38Z maronga |
---|
30 | ! Corrected "Former revisions" section |
---|
31 | ! |
---|
32 | ! 2696 2017-12-14 17:12:51Z kanani |
---|
33 | ! Change in file header (GPL part) |
---|
34 | ! |
---|
35 | ! 2101 2017-01-05 16:42:31Z suehring |
---|
36 | ! |
---|
37 | ! 2000 2016-08-20 18:09:15Z knoop |
---|
38 | ! Forced header and separation lines into 80 columns |
---|
39 | ! |
---|
40 | ! 1880 2016-04-20 09:36:50Z hoffmann |
---|
41 | ! Bugfix: The index of the larger particle has to be chosen for interpolation. |
---|
42 | ! |
---|
43 | ! 1873 2016-04-18 14:50:06Z maronga |
---|
44 | ! Module renamed (removed _mod) |
---|
45 | ! |
---|
46 | ! 1858 2016-04-13 13:12:11Z hoffmann |
---|
47 | ! Interpolation of collision kernels adjusted to more reasonable values. |
---|
48 | ! Reformatting of the code. |
---|
49 | ! |
---|
50 | ! 1850 2016-04-08 13:29:27Z maronga |
---|
51 | ! Module renamed |
---|
52 | ! |
---|
53 | ! 1822 2016-04-07 07:49:42Z hoffmann |
---|
54 | ! PALM kernel has been deleted. |
---|
55 | ! Bugfix in the calculation of the turbulent enhancement factor of the |
---|
56 | ! collection efficiency. |
---|
57 | ! |
---|
58 | ! Unused variables removed. |
---|
59 | ! |
---|
60 | ! 1776 2016-03-02 17:54:58Z hoffmann |
---|
61 | ! Bugfix: Collection efficiencies must be calculated for the larger droplet. |
---|
62 | ! |
---|
63 | ! 1682 2015-10-07 23:56:08Z knoop |
---|
64 | ! Code annotations made doxygen readable |
---|
65 | ! |
---|
66 | ! 1519 2015-01-08 10:20:42Z hoffmann |
---|
67 | ! Bugfix: Using the new particle structure, particles are not sorted by size. |
---|
68 | ! Hence, computation of collision efficiencies must ensure that the ratio of |
---|
69 | ! two colliding droplets is < 1. |
---|
70 | ! |
---|
71 | ! 1359 2014-04-11 17:15:14Z hoffmann |
---|
72 | ! New particle structure integrated. |
---|
73 | ! Kind definition added to all floating point numbers. |
---|
74 | ! |
---|
75 | ! 1346 2014-03-27 13:18:20Z heinze |
---|
76 | ! Bugfix: REAL constants provided with KIND-attribute especially in call of |
---|
77 | ! intrinsic function like MAX, MIN, SIGN |
---|
78 | ! |
---|
79 | ! 1322 2014-03-20 16:38:49Z raasch |
---|
80 | ! REAL constants defined as wp_kind |
---|
81 | ! |
---|
82 | ! 1320 2014-03-20 08:40:49Z |
---|
83 | ! ONLY-attribute added to USE-statements, |
---|
84 | ! kind-parameters added to all INTEGER and REAL declaration statements, |
---|
85 | ! kinds are defined in new module kinds, |
---|
86 | ! revision history before 2012 removed, |
---|
87 | ! comment fields (!:) to be used for variable explanations added to |
---|
88 | ! all variable declaration statements |
---|
89 | ! |
---|
90 | ! 1092 2013-02-02 11:24:22Z raasch |
---|
91 | ! unused variables removed |
---|
92 | ! |
---|
93 | ! 1071 2012-11-29 16:54:55Z franke |
---|
94 | ! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20 |
---|
95 | ! |
---|
96 | ! 1036 2012-10-22 13:43:42Z raasch |
---|
97 | ! code put under GPL (PALM 3.9) |
---|
98 | ! |
---|
99 | ! 1007 2012-09-19 14:30:36Z franke |
---|
100 | ! converted all units to SI units and replaced some parameters by corresponding |
---|
101 | ! PALM parameters |
---|
102 | ! Bugfix: factor in calculation of enhancement factor for collision efficencies |
---|
103 | ! changed from 10. to 1.0 |
---|
104 | ! |
---|
105 | ! 849 2012-03-15 10:35:09Z raasch |
---|
106 | ! routine collision_efficiency_rogers added (moved from former advec_particles |
---|
107 | ! to here) |
---|
108 | ! |
---|
109 | ! 835 2012-02-22 11:21:19Z raasch $ |
---|
110 | ! Bugfix: array diss can be used only in case of Wang kernel |
---|
111 | ! |
---|
112 | ! 828 2012-02-21 12:00:36Z raasch |
---|
113 | ! code has been completely reformatted, routine colker renamed |
---|
114 | ! recalculate_kernel, |
---|
115 | ! routine init_kernels added, radius is now communicated to the collision |
---|
116 | ! routines by array radclass |
---|
117 | ! |
---|
118 | ! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4 |
---|
119 | ! |
---|
120 | ! 825 2012-02-19 03:03:44Z raasch |
---|
121 | ! routine renamed from wang_kernel to lpm_collision_kernels, |
---|
122 | ! turbulence_effects on collision replaced by wang_kernel |
---|
123 | ! |
---|
124 | ! 790 2011-11-29 03:11:20Z raasch |
---|
125 | ! initial revision |
---|
126 | ! |
---|
127 | ! Description: |
---|
128 | ! ------------ |
---|
129 | !> This module calculates collision efficiencies either due to pure gravitational |
---|
130 | !> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or |
---|
131 | !> including the effects of turbulence (Wang kernel, see Wang and |
---|
132 | !> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8, and Ayala et al., 2008: |
---|
133 | !> New J. Phys., 10, 075016). The original code has been |
---|
134 | !> provided by L.-P. Wang but is substantially reformatted and speed optimized |
---|
135 | !> here. |
---|
136 | !------------------------------------------------------------------------------! |
---|
137 | MODULE lpm_collision_kernels_mod |
---|
138 | |
---|
139 | |
---|
140 | USE basic_constants_and_equations_mod, & |
---|
141 | ONLY: g, pi |
---|
142 | |
---|
143 | USE kinds |
---|
144 | |
---|
145 | USE particle_attributes, & |
---|
146 | ONLY: collision_kernel, dissipation_classes, particles, & |
---|
147 | radius_classes |
---|
148 | |
---|
149 | USE pegrid |
---|
150 | |
---|
151 | |
---|
152 | IMPLICIT NONE |
---|
153 | |
---|
154 | PRIVATE |
---|
155 | |
---|
156 | PUBLIC ckernel, init_kernels, rclass_lbound, rclass_ubound, & |
---|
157 | recalculate_kernel |
---|
158 | |
---|
159 | REAL(wp) :: epsilon !< |
---|
160 | REAL(wp) :: rclass_lbound !< |
---|
161 | REAL(wp) :: rclass_ubound !< |
---|
162 | REAL(wp) :: urms !< |
---|
163 | |
---|
164 | REAL(wp), DIMENSION(:), ALLOCATABLE :: epsclass !< dissipation rate class |
---|
165 | REAL(wp), DIMENSION(:), ALLOCATABLE :: radclass !< radius class |
---|
166 | REAL(wp), DIMENSION(:), ALLOCATABLE :: winf !< |
---|
167 | |
---|
168 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ec !< |
---|
169 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ecf !< |
---|
170 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gck !< |
---|
171 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hkernel !< |
---|
172 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hwratio !< |
---|
173 | |
---|
174 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel !< |
---|
175 | |
---|
176 | SAVE |
---|
177 | |
---|
178 | ! |
---|
179 | !-- Public interfaces |
---|
180 | INTERFACE init_kernels |
---|
181 | MODULE PROCEDURE init_kernels |
---|
182 | END INTERFACE init_kernels |
---|
183 | |
---|
184 | INTERFACE recalculate_kernel |
---|
185 | MODULE PROCEDURE recalculate_kernel |
---|
186 | END INTERFACE recalculate_kernel |
---|
187 | |
---|
188 | |
---|
189 | CONTAINS |
---|
190 | |
---|
191 | |
---|
192 | !------------------------------------------------------------------------------! |
---|
193 | ! Description: |
---|
194 | ! ------------ |
---|
195 | !> Initialization of the collision efficiency matrix with fixed radius and |
---|
196 | !> dissipation classes, calculated at simulation start only. |
---|
197 | !------------------------------------------------------------------------------! |
---|
198 | |
---|
199 | SUBROUTINE init_kernels |
---|
200 | |
---|
201 | IMPLICIT NONE |
---|
202 | |
---|
203 | INTEGER(iwp) :: i !< |
---|
204 | INTEGER(iwp) :: j !< |
---|
205 | INTEGER(iwp) :: k !< |
---|
206 | |
---|
207 | |
---|
208 | ! |
---|
209 | !-- Calculate collision efficiencies for fixed radius- and dissipation |
---|
210 | !-- classes |
---|
211 | IF ( collision_kernel(6:9) == 'fast' ) THEN |
---|
212 | |
---|
213 | ALLOCATE( ckernel(1:radius_classes,1:radius_classes, & |
---|
214 | 0:dissipation_classes), epsclass(1:dissipation_classes), & |
---|
215 | radclass(1:radius_classes) ) |
---|
216 | |
---|
217 | ! |
---|
218 | !-- Calculate the radius class bounds with logarithmic distances |
---|
219 | !-- in the interval [1.0E-6, 1000.0E-6] m |
---|
220 | rclass_lbound = LOG( 1.0E-6_wp ) |
---|
221 | rclass_ubound = LOG( 1000.0E-6_wp ) |
---|
222 | radclass(1) = EXP( rclass_lbound ) |
---|
223 | DO i = 2, radius_classes |
---|
224 | radclass(i) = EXP( rclass_lbound + & |
---|
225 | ( rclass_ubound - rclass_lbound ) * & |
---|
226 | ( i - 1.0_wp ) / ( radius_classes - 1.0_wp ) ) |
---|
227 | ENDDO |
---|
228 | |
---|
229 | ! |
---|
230 | !-- Set the class bounds for dissipation in interval [0.0, 600.0] cm**2/s**3 |
---|
231 | DO i = 1, dissipation_classes |
---|
232 | epsclass(i) = 0.06_wp * REAL( i, KIND=wp ) / dissipation_classes |
---|
233 | ENDDO |
---|
234 | ! |
---|
235 | !-- Calculate collision efficiencies of the Wang/ayala kernel |
---|
236 | ALLOCATE( ec(1:radius_classes,1:radius_classes), & |
---|
237 | ecf(1:radius_classes,1:radius_classes), & |
---|
238 | gck(1:radius_classes,1:radius_classes), & |
---|
239 | winf(1:radius_classes) ) |
---|
240 | |
---|
241 | DO k = 1, dissipation_classes |
---|
242 | |
---|
243 | epsilon = epsclass(k) |
---|
244 | urms = 2.02_wp * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp ) |
---|
245 | |
---|
246 | CALL turbsd |
---|
247 | CALL turb_enhance_eff |
---|
248 | CALL effic |
---|
249 | |
---|
250 | DO j = 1, radius_classes |
---|
251 | DO i = 1, radius_classes |
---|
252 | ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j) |
---|
253 | ENDDO |
---|
254 | ENDDO |
---|
255 | |
---|
256 | ENDDO |
---|
257 | |
---|
258 | ! |
---|
259 | !-- Calculate collision efficiencies of the Hall kernel |
---|
260 | ALLOCATE( hkernel(1:radius_classes,1:radius_classes), & |
---|
261 | hwratio(1:radius_classes,1:radius_classes) ) |
---|
262 | |
---|
263 | CALL fallg |
---|
264 | CALL effic |
---|
265 | |
---|
266 | DO j = 1, radius_classes |
---|
267 | DO i = 1, radius_classes |
---|
268 | hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 & |
---|
269 | * ec(i,j) * ABS( winf(j) - winf(i) ) |
---|
270 | ckernel(i,j,0) = hkernel(i,j) ! hall kernel stored on index 0 |
---|
271 | ENDDO |
---|
272 | ENDDO |
---|
273 | |
---|
274 | ! |
---|
275 | !-- Test output of efficiencies |
---|
276 | IF ( j == -1 ) THEN |
---|
277 | |
---|
278 | PRINT*, '*** Hall kernel' |
---|
279 | WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6_wp, & |
---|
280 | i = 1,radius_classes ) |
---|
281 | DO j = 1, radius_classes |
---|
282 | WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j), & |
---|
283 | ( hkernel(i,j), i = 1,radius_classes ) |
---|
284 | ENDDO |
---|
285 | |
---|
286 | DO k = 1, dissipation_classes |
---|
287 | DO i = 1, radius_classes |
---|
288 | DO j = 1, radius_classes |
---|
289 | IF ( hkernel(i,j) == 0.0_wp ) THEN |
---|
290 | hwratio(i,j) = 9999999.9_wp |
---|
291 | ELSE |
---|
292 | hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j) |
---|
293 | ENDIF |
---|
294 | ENDDO |
---|
295 | ENDDO |
---|
296 | |
---|
297 | PRINT*, '*** epsilon = ', epsclass(k) |
---|
298 | WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i) * 1.0E6_wp, & |
---|
299 | i = 1,radius_classes ) |
---|
300 | DO j = 1, radius_classes |
---|
301 | WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j) * 1.0E6_wp, & |
---|
302 | ( hwratio(i,j), i = 1,radius_classes ) |
---|
303 | ENDDO |
---|
304 | ENDDO |
---|
305 | |
---|
306 | ENDIF |
---|
307 | |
---|
308 | DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf ) |
---|
309 | |
---|
310 | ENDIF |
---|
311 | |
---|
312 | END SUBROUTINE init_kernels |
---|
313 | |
---|
314 | |
---|
315 | !------------------------------------------------------------------------------! |
---|
316 | ! Description: |
---|
317 | ! ------------ |
---|
318 | !> Calculation of collision kernels during each timestep and for each grid box |
---|
319 | !------------------------------------------------------------------------------! |
---|
320 | SUBROUTINE recalculate_kernel( i1, j1, k1 ) |
---|
321 | |
---|
322 | USE arrays_3d, & |
---|
323 | ONLY: diss |
---|
324 | |
---|
325 | USE particle_attributes, & |
---|
326 | ONLY: number_of_particles, prt_count, radius_classes, wang_kernel |
---|
327 | |
---|
328 | IMPLICIT NONE |
---|
329 | |
---|
330 | INTEGER(iwp) :: i !< |
---|
331 | INTEGER(iwp) :: i1 !< |
---|
332 | INTEGER(iwp) :: j !< |
---|
333 | INTEGER(iwp) :: j1 !< |
---|
334 | INTEGER(iwp) :: k1 !< |
---|
335 | |
---|
336 | |
---|
337 | number_of_particles = prt_count(k1,j1,i1) |
---|
338 | radius_classes = number_of_particles ! necessary to use the same |
---|
339 | ! subroutines as for |
---|
340 | ! precalculated kernels |
---|
341 | |
---|
342 | ALLOCATE( ec(1:number_of_particles,1:number_of_particles), & |
---|
343 | radclass(1:number_of_particles), winf(1:number_of_particles) ) |
---|
344 | |
---|
345 | ! |
---|
346 | !-- Store particle radii on the radclass array |
---|
347 | radclass(1:number_of_particles) = particles(1:number_of_particles)%radius |
---|
348 | |
---|
349 | IF ( wang_kernel ) THEN |
---|
350 | epsilon = diss(k1,j1,i1) ! dissipation rate in m**2/s**3 |
---|
351 | ELSE |
---|
352 | epsilon = 0.0_wp |
---|
353 | ENDIF |
---|
354 | urms = 2.02_wp * ( epsilon / 0.04_wp )**( 0.33333333333_wp ) |
---|
355 | |
---|
356 | IF ( wang_kernel .AND. epsilon > 1.0E-7_wp ) THEN |
---|
357 | ! |
---|
358 | !-- Call routines to calculate efficiencies for the Wang kernel |
---|
359 | ALLOCATE( gck(1:number_of_particles,1:number_of_particles), & |
---|
360 | ecf(1:number_of_particles,1:number_of_particles) ) |
---|
361 | |
---|
362 | CALL turbsd |
---|
363 | CALL turb_enhance_eff |
---|
364 | CALL effic |
---|
365 | |
---|
366 | DO j = 1, number_of_particles |
---|
367 | DO i = 1, number_of_particles |
---|
368 | ckernel(1+i-1,1+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j) |
---|
369 | ENDDO |
---|
370 | ENDDO |
---|
371 | |
---|
372 | DEALLOCATE( gck, ecf ) |
---|
373 | |
---|
374 | ELSE |
---|
375 | ! |
---|
376 | !-- Call routines to calculate efficiencies for the Hall kernel |
---|
377 | CALL fallg |
---|
378 | CALL effic |
---|
379 | |
---|
380 | DO j = 1, number_of_particles |
---|
381 | DO i = 1, number_of_particles |
---|
382 | ckernel(i,j,1) = pi * ( radclass(j) + radclass(i) )**2 & |
---|
383 | * ec(i,j) * ABS( winf(j) - winf(i) ) |
---|
384 | ENDDO |
---|
385 | ENDDO |
---|
386 | |
---|
387 | ENDIF |
---|
388 | |
---|
389 | DEALLOCATE( ec, radclass, winf ) |
---|
390 | |
---|
391 | END SUBROUTINE recalculate_kernel |
---|
392 | |
---|
393 | |
---|
394 | !------------------------------------------------------------------------------! |
---|
395 | ! Description: |
---|
396 | ! ------------ |
---|
397 | !> Calculation of effects of turbulence on the geometric collision kernel |
---|
398 | !> (by including the droplets' average radial relative velocities and their |
---|
399 | !> radial distribution function) following the analytic model by Aayala et al. |
---|
400 | !> (2008, New J. Phys.). For details check the second part 2 of the publication, |
---|
401 | !> page 37ff. |
---|
402 | !> |
---|
403 | !> Input parameters, which need to be replaced by PALM parameters: |
---|
404 | !> water density, air density |
---|
405 | !------------------------------------------------------------------------------! |
---|
406 | SUBROUTINE turbsd |
---|
407 | |
---|
408 | USE control_parameters, & |
---|
409 | ONLY: molecular_viscosity |
---|
410 | |
---|
411 | USE particle_attributes, & |
---|
412 | ONLY: radius_classes |
---|
413 | |
---|
414 | IMPLICIT NONE |
---|
415 | |
---|
416 | INTEGER(iwp) :: i !< |
---|
417 | INTEGER(iwp) :: j !< |
---|
418 | |
---|
419 | REAL(wp) :: ao !< |
---|
420 | REAL(wp) :: ao_gr !< |
---|
421 | REAL(wp) :: bbb !< |
---|
422 | REAL(wp) :: be !< |
---|
423 | REAL(wp) :: b1 !< |
---|
424 | REAL(wp) :: b2 !< |
---|
425 | REAL(wp) :: ccc !< |
---|
426 | REAL(wp) :: c1 !< |
---|
427 | REAL(wp) :: c1_gr !< |
---|
428 | REAL(wp) :: c2 !< |
---|
429 | REAL(wp) :: d1 !< |
---|
430 | REAL(wp) :: d2 !< |
---|
431 | REAL(wp) :: eta !< |
---|
432 | REAL(wp) :: e1 !< |
---|
433 | REAL(wp) :: e2 !< |
---|
434 | REAL(wp) :: fao_gr !< |
---|
435 | REAL(wp) :: fr !< |
---|
436 | REAL(wp) :: grfin !< |
---|
437 | REAL(wp) :: lambda !< |
---|
438 | REAL(wp) :: lambda_re !< |
---|
439 | REAL(wp) :: lf !< |
---|
440 | REAL(wp) :: rc !< |
---|
441 | REAL(wp) :: rrp !< |
---|
442 | REAL(wp) :: sst !< |
---|
443 | REAL(wp) :: tauk !< |
---|
444 | REAL(wp) :: tl !< |
---|
445 | REAL(wp) :: t2 !< |
---|
446 | REAL(wp) :: tt !< |
---|
447 | REAL(wp) :: t1 !< |
---|
448 | REAL(wp) :: vk !< |
---|
449 | REAL(wp) :: vrms1xy !< |
---|
450 | REAL(wp) :: vrms2xy !< |
---|
451 | REAL(wp) :: v1 !< |
---|
452 | REAL(wp) :: v1v2xy !< |
---|
453 | REAL(wp) :: v1xysq !< |
---|
454 | REAL(wp) :: v2 !< |
---|
455 | REAL(wp) :: v2xysq !< |
---|
456 | REAL(wp) :: wrfin !< |
---|
457 | REAL(wp) :: wrgrav2 !< |
---|
458 | REAL(wp) :: wrtur2xy !< |
---|
459 | REAL(wp) :: xx !< |
---|
460 | REAL(wp) :: yy !< |
---|
461 | REAL(wp) :: z !< |
---|
462 | |
---|
463 | REAL(wp), DIMENSION(1:radius_classes) :: st !< Stokes number |
---|
464 | REAL(wp), DIMENSION(1:radius_classes) :: tau !< inertial time scale |
---|
465 | |
---|
466 | lambda = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon ) |
---|
467 | lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity ) |
---|
468 | tl = urms**2 / epsilon |
---|
469 | lf = 0.5_wp * urms**3 / epsilon |
---|
470 | tauk = SQRT( molecular_viscosity / epsilon ) |
---|
471 | eta = ( molecular_viscosity**3 / epsilon )**0.25_wp |
---|
472 | vk = eta / tauk |
---|
473 | |
---|
474 | ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re ) |
---|
475 | tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk |
---|
476 | |
---|
477 | ! |
---|
478 | !-- Get terminal velocity of droplets |
---|
479 | CALL fallg |
---|
480 | |
---|
481 | DO i = 1, radius_classes |
---|
482 | tau(i) = winf(i) / g ! inertial time scale |
---|
483 | st(i) = tau(i) / tauk ! Stokes number |
---|
484 | ENDDO |
---|
485 | |
---|
486 | ! |
---|
487 | !-- Calculate average radial relative velocity at contact (wrfin) |
---|
488 | z = tt / tl |
---|
489 | be = SQRT( 2.0_wp ) * lambda / lf |
---|
490 | bbb = SQRT( 1.0_wp - 2.0_wp * be**2 ) |
---|
491 | d1 = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb ) |
---|
492 | e1 = lf * ( 1.0_wp + bbb ) * 0.5_wp |
---|
493 | d2 = ( 1.0_wp - bbb ) * 0.5_wp / bbb |
---|
494 | e2 = lf * ( 1.0_wp - bbb ) * 0.5_wp |
---|
495 | ccc = SQRT( 1.0_wp - 2.0_wp * z**2 ) |
---|
496 | b1 = ( 1.0_wp + ccc ) * 0.5_wp / ccc |
---|
497 | c1 = tl * ( 1.0_wp + ccc ) * 0.5_wp |
---|
498 | b2 = ( 1.0_wp - ccc ) * 0.5_wp / ccc |
---|
499 | c2 = tl * ( 1.0_wp - ccc ) * 0.5_wp |
---|
500 | |
---|
501 | DO i = 1, radius_classes |
---|
502 | |
---|
503 | v1 = winf(i) |
---|
504 | t1 = tau(i) |
---|
505 | |
---|
506 | DO j = 1, i |
---|
507 | rrp = radclass(i) + radclass(j) |
---|
508 | v2 = winf(j) |
---|
509 | t2 = tau(j) |
---|
510 | |
---|
511 | v1xysq = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) & |
---|
512 | - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1) |
---|
513 | v1xysq = v1xysq * urms**2 / t1 |
---|
514 | vrms1xy = SQRT( v1xysq ) |
---|
515 | |
---|
516 | v2xysq = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) & |
---|
517 | - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2) |
---|
518 | v2xysq = v2xysq * urms**2 / t2 |
---|
519 | vrms2xy = SQRT( v2xysq ) |
---|
520 | |
---|
521 | IF ( winf(i) >= winf(j) ) THEN |
---|
522 | v1 = winf(i) |
---|
523 | t1 = tau(i) |
---|
524 | v2 = winf(j) |
---|
525 | t2 = tau(j) |
---|
526 | ELSE |
---|
527 | v1 = winf(j) |
---|
528 | t1 = tau(j) |
---|
529 | v2 = winf(i) |
---|
530 | t2 = tau(i) |
---|
531 | ENDIF |
---|
532 | |
---|
533 | v1v2xy = b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - & |
---|
534 | b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - & |
---|
535 | b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + & |
---|
536 | b2 * d2* zhi(c2,e2,v1,t1,v2,t2) |
---|
537 | fr = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 ) |
---|
538 | v1v2xy = v1v2xy * fr * urms**2 / tau(i) / tau(j) |
---|
539 | wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy |
---|
540 | IF ( wrtur2xy < 0.0_wp ) wrtur2xy = 0.0_wp |
---|
541 | wrgrav2 = pi / 8.0_wp * ( winf(j) - winf(i) )**2 |
---|
542 | wrfin = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) |
---|
543 | |
---|
544 | ! |
---|
545 | !-- Calculate radial distribution function (grfin) |
---|
546 | IF ( st(j) > st(i) ) THEN |
---|
547 | sst = st(j) |
---|
548 | ELSE |
---|
549 | sst = st(i) |
---|
550 | ENDIF |
---|
551 | |
---|
552 | xx = -0.1988_wp * sst**4 + 1.5275_wp * sst**3 - 4.2942_wp * & |
---|
553 | sst**2 + 5.3406_wp * sst |
---|
554 | IF ( xx < 0.0_wp ) xx = 0.0_wp |
---|
555 | yy = 0.1886_wp * EXP( 20.306_wp / lambda_re ) |
---|
556 | |
---|
557 | c1_gr = xx / ( g / vk * tauk )**yy |
---|
558 | |
---|
559 | ao_gr = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2 |
---|
560 | fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re ) |
---|
561 | rc = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta |
---|
562 | |
---|
563 | grfin = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5_wp ) |
---|
564 | IF ( grfin < 1.0_wp ) grfin = 1.0_wp |
---|
565 | |
---|
566 | ! |
---|
567 | !-- Calculate general collection kernel (without the consideration of |
---|
568 | !-- collection efficiencies) |
---|
569 | gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin |
---|
570 | gck(j,i) = gck(i,j) |
---|
571 | |
---|
572 | ENDDO |
---|
573 | ENDDO |
---|
574 | |
---|
575 | END SUBROUTINE turbsd |
---|
576 | |
---|
577 | REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 ) |
---|
578 | ! |
---|
579 | !-- Function used in the Ayala et al. (2008) analytical model for turbulent |
---|
580 | !-- effects on the collision kernel |
---|
581 | IMPLICIT NONE |
---|
582 | |
---|
583 | REAL(wp) :: a !< |
---|
584 | REAL(wp) :: aa1 !< |
---|
585 | REAL(wp) :: b !< |
---|
586 | REAL(wp) :: tau0 !< |
---|
587 | REAL(wp) :: vsett !< |
---|
588 | |
---|
589 | aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b |
---|
590 | phi_w = 1.0_wp / aa1 - 0.5_wp * vsett / b / aa1**2 |
---|
591 | |
---|
592 | END FUNCTION phi_w |
---|
593 | |
---|
594 | REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 ) |
---|
595 | ! |
---|
596 | !-- Function used in the Ayala et al. (2008) analytical model for turbulent |
---|
597 | !-- effects on the collision kernel |
---|
598 | IMPLICIT NONE |
---|
599 | |
---|
600 | REAL(wp) :: a !< |
---|
601 | REAL(wp) :: aa1 !< |
---|
602 | REAL(wp) :: aa2 !< |
---|
603 | REAL(wp) :: aa3 !< |
---|
604 | REAL(wp) :: aa4 !< |
---|
605 | REAL(wp) :: aa5 !< |
---|
606 | REAL(wp) :: aa6 !< |
---|
607 | REAL(wp) :: b !< |
---|
608 | REAL(wp) :: tau1 !< |
---|
609 | REAL(wp) :: tau2 !< |
---|
610 | REAL(wp) :: vsett1 !< |
---|
611 | REAL(wp) :: vsett2 !< |
---|
612 | |
---|
613 | aa1 = vsett2 / b - 1.0_wp / tau2 - 1.0_wp / a |
---|
614 | aa2 = vsett1 / b + 1.0_wp / tau1 + 1.0_wp / a |
---|
615 | aa3 = ( vsett1 - vsett2 ) / b + 1.0_wp / tau1 + 1.0_wp / tau2 |
---|
616 | aa4 = ( vsett2 / b )**2 - ( 1.0_wp / tau2 + 1.0_wp / a )**2 |
---|
617 | aa5 = vsett2 / b + 1.0_wp / tau2 + 1.0_wp / a |
---|
618 | aa6 = 1.0_wp / tau1 - 1.0_wp / a + ( 1.0_wp / tau2 + 1.0_wp / a) * & |
---|
619 | vsett1 / vsett2 |
---|
620 | zhi = (1.0_wp / aa1 - 1.0_wp / aa2 ) * ( vsett1 - vsett2 ) * 0.5_wp / & |
---|
621 | b / aa3**2 + ( 4.0_wp / aa4 - 1.0_wp / aa5**2 - 1.0_wp / aa1**2 ) & |
---|
622 | * vsett2 * 0.5_wp / b /aa6 + ( 2.0_wp * ( b / aa2 - b / aa1 ) - & |
---|
623 | vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3 |
---|
624 | |
---|
625 | END FUNCTION zhi |
---|
626 | |
---|
627 | |
---|
628 | !------------------------------------------------------------------------------! |
---|
629 | ! Description: |
---|
630 | ! ------------ |
---|
631 | !> Parameterization of terminal velocity following Rogers et al. (1993, J. Appl. |
---|
632 | !> Meteorol.) |
---|
633 | !------------------------------------------------------------------------------! |
---|
634 | SUBROUTINE fallg |
---|
635 | |
---|
636 | USE particle_attributes, & |
---|
637 | ONLY: radius_classes |
---|
638 | |
---|
639 | IMPLICIT NONE |
---|
640 | |
---|
641 | INTEGER(iwp) :: j !< |
---|
642 | |
---|
643 | REAL(wp), PARAMETER :: k_cap_rog = 4.0_wp !< parameter |
---|
644 | REAL(wp), PARAMETER :: k_low_rog = 12.0_wp !< parameter |
---|
645 | REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter |
---|
646 | REAL(wp), PARAMETER :: b_rog = 10.43_wp !< parameter |
---|
647 | REAL(wp), PARAMETER :: c_rog = 0.6_wp !< parameter |
---|
648 | REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< seperation diameter |
---|
649 | |
---|
650 | REAL(wp) :: diameter !< droplet diameter in mm |
---|
651 | |
---|
652 | |
---|
653 | DO j = 1, radius_classes |
---|
654 | |
---|
655 | diameter = radclass(j) * 2000.0_wp |
---|
656 | |
---|
657 | IF ( diameter <= d0_rog ) THEN |
---|
658 | winf(j) = k_cap_rog * diameter * ( 1.0_wp - & |
---|
659 | EXP( -k_low_rog * diameter ) ) |
---|
660 | ELSE |
---|
661 | winf(j) = a_rog - b_rog * EXP( -c_rog * diameter ) |
---|
662 | ENDIF |
---|
663 | |
---|
664 | ENDDO |
---|
665 | |
---|
666 | END SUBROUTINE fallg |
---|
667 | |
---|
668 | |
---|
669 | !------------------------------------------------------------------------------! |
---|
670 | ! Description: |
---|
671 | ! ------------ |
---|
672 | !> Interpolation of collision efficiencies (Hall, 1980, J. Atmos. Sci.) |
---|
673 | !------------------------------------------------------------------------------! |
---|
674 | SUBROUTINE effic |
---|
675 | |
---|
676 | USE particle_attributes, & |
---|
677 | ONLY: radius_classes |
---|
678 | |
---|
679 | IMPLICIT NONE |
---|
680 | |
---|
681 | INTEGER(iwp) :: i !< |
---|
682 | INTEGER(iwp) :: iq !< |
---|
683 | INTEGER(iwp) :: ir !< |
---|
684 | INTEGER(iwp) :: j !< |
---|
685 | INTEGER(iwp) :: k !< |
---|
686 | |
---|
687 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< |
---|
688 | |
---|
689 | LOGICAL, SAVE :: first = .TRUE. !< |
---|
690 | |
---|
691 | REAL(wp) :: ek !< |
---|
692 | REAL(wp) :: particle_radius !< |
---|
693 | REAL(wp) :: pp !< |
---|
694 | REAL(wp) :: qq !< |
---|
695 | REAL(wp) :: rq !< |
---|
696 | |
---|
697 | REAL(wp), DIMENSION(1:21), SAVE :: rat !< |
---|
698 | |
---|
699 | REAL(wp), DIMENSION(1:15), SAVE :: r0 !< |
---|
700 | |
---|
701 | REAL(wp), DIMENSION(1:15,1:21), SAVE :: ecoll !< |
---|
702 | |
---|
703 | ! |
---|
704 | !-- Initial assignment of constants |
---|
705 | IF ( first ) THEN |
---|
706 | |
---|
707 | first = .FALSE. |
---|
708 | r0 = (/ 6.0_wp, 8.0_wp, 10.0_wp, 15.0_wp, 20.0_wp, 25.0_wp, & |
---|
709 | 30.0_wp, 40.0_wp, 50.0_wp, 60.0_wp, 70.0_wp, 100.0_wp, & |
---|
710 | 150.0_wp, 200.0_wp, 300.0_wp /) |
---|
711 | |
---|
712 | rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp, & |
---|
713 | 0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp, & |
---|
714 | 0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp, & |
---|
715 | 0.90_wp, 0.95_wp, 1.00_wp /) |
---|
716 | |
---|
717 | ecoll(:,1) = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, & |
---|
718 | 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, & |
---|
719 | 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /) |
---|
720 | ecoll(:,2) = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp, & |
---|
721 | 0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp, & |
---|
722 | 0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /) |
---|
723 | ecoll(:,3) = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp, & |
---|
724 | 0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp, & |
---|
725 | 0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /) |
---|
726 | ecoll(:,4) = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp, & |
---|
727 | 0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp, & |
---|
728 | 0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /) |
---|
729 | ecoll(:,5) = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp, & |
---|
730 | 0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp, & |
---|
731 | 0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /) |
---|
732 | ecoll(:,6) = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp, & |
---|
733 | 0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp, & |
---|
734 | 0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
735 | ecoll(:,7) = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp, & |
---|
736 | 0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp, & |
---|
737 | 0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
738 | ecoll(:,8) = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp, & |
---|
739 | 0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp, & |
---|
740 | 0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
741 | ecoll(:,9) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp, & |
---|
742 | 0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp, & |
---|
743 | 0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
744 | ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp, & |
---|
745 | 0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp, & |
---|
746 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
747 | ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp, & |
---|
748 | 0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp, & |
---|
749 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
750 | ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp, & |
---|
751 | 0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp, & |
---|
752 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
753 | ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp, & |
---|
754 | 0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp, & |
---|
755 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
756 | ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp, & |
---|
757 | 0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp, & |
---|
758 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
759 | ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp, & |
---|
760 | 0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp, & |
---|
761 | 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
762 | ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp, & |
---|
763 | 0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp, & |
---|
764 | 0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
765 | ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp, & |
---|
766 | 0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp, & |
---|
767 | 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /) |
---|
768 | ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp, & |
---|
769 | 0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp, & |
---|
770 | 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /) |
---|
771 | ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp, & |
---|
772 | 0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp, & |
---|
773 | 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /) |
---|
774 | ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, & |
---|
775 | 0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp, & |
---|
776 | 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /) |
---|
777 | ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, & |
---|
778 | 0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp, & |
---|
779 | 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /) |
---|
780 | ENDIF |
---|
781 | |
---|
782 | ! |
---|
783 | !-- Calculate the radius class index of particles with respect to array r |
---|
784 | !-- Radius has to be in microns |
---|
785 | ALLOCATE( ira(1:radius_classes) ) |
---|
786 | DO j = 1, radius_classes |
---|
787 | particle_radius = radclass(j) * 1.0E6_wp |
---|
788 | DO k = 1, 15 |
---|
789 | IF ( particle_radius < r0(k) ) THEN |
---|
790 | ira(j) = k |
---|
791 | EXIT |
---|
792 | ENDIF |
---|
793 | ENDDO |
---|
794 | IF ( particle_radius >= r0(15) ) ira(j) = 16 |
---|
795 | ENDDO |
---|
796 | |
---|
797 | ! |
---|
798 | !-- Two-dimensional linear interpolation of the collision efficiency. |
---|
799 | !-- Radius has to be in microns |
---|
800 | DO j = 1, radius_classes |
---|
801 | DO i = 1, j |
---|
802 | |
---|
803 | ir = MAX( ira(i), ira(j) ) |
---|
804 | rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) ) |
---|
805 | iq = INT( rq * 20 ) + 1 |
---|
806 | iq = MAX( iq , 2) |
---|
807 | |
---|
808 | IF ( ir < 16 ) THEN |
---|
809 | IF ( ir >= 2 ) THEN |
---|
810 | pp = ( ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp ) - & |
---|
811 | r0(ir-1) ) / ( r0(ir) - r0(ir-1) ) |
---|
812 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
813 | ec(j,i) = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) & |
---|
814 | * ecoll(ir-1,iq-1) & |
---|
815 | + pp * ( 1.0_wp - qq ) * ecoll(ir,iq-1) & |
---|
816 | + qq * ( 1.0_wp - pp ) * ecoll(ir-1,iq) & |
---|
817 | + pp * qq * ecoll(ir,iq) |
---|
818 | ELSE |
---|
819 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
820 | ec(j,i) = ( 1.0_wp - qq ) * ecoll(1,iq-1) + qq * ecoll(1,iq) |
---|
821 | ENDIF |
---|
822 | ELSE |
---|
823 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
824 | ek = ( 1.0_wp - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) |
---|
825 | ec(j,i) = MIN( ek, 1.0_wp ) |
---|
826 | ENDIF |
---|
827 | |
---|
828 | IF ( ec(j,i) < 1.0E-20_wp ) ec(j,i) = 0.0_wp |
---|
829 | |
---|
830 | ec(i,j) = ec(j,i) |
---|
831 | |
---|
832 | ENDDO |
---|
833 | ENDDO |
---|
834 | |
---|
835 | DEALLOCATE( ira ) |
---|
836 | |
---|
837 | END SUBROUTINE effic |
---|
838 | |
---|
839 | |
---|
840 | !------------------------------------------------------------------------------! |
---|
841 | ! Description: |
---|
842 | ! ------------ |
---|
843 | !> Interpolation of turbulent enhancement factor for collision efficencies |
---|
844 | !> following Wang and Grabowski (2009, Atmos. Sci. Let.) |
---|
845 | !------------------------------------------------------------------------------! |
---|
846 | SUBROUTINE turb_enhance_eff |
---|
847 | |
---|
848 | USE particle_attributes, & |
---|
849 | ONLY: radius_classes |
---|
850 | |
---|
851 | IMPLICIT NONE |
---|
852 | |
---|
853 | INTEGER(iwp) :: i !< |
---|
854 | INTEGER(iwp) :: iq !< |
---|
855 | INTEGER(iwp) :: ir !< |
---|
856 | INTEGER(iwp) :: j !< |
---|
857 | INTEGER(iwp) :: k !< |
---|
858 | INTEGER(iwp) :: kk !< |
---|
859 | |
---|
860 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< |
---|
861 | |
---|
862 | LOGICAL, SAVE :: first = .TRUE. !< |
---|
863 | |
---|
864 | REAL(wp) :: particle_radius !< |
---|
865 | REAL(wp) :: pp !< |
---|
866 | REAL(wp) :: qq !< |
---|
867 | REAL(wp) :: rq !< |
---|
868 | REAL(wp) :: y1 !< |
---|
869 | REAL(wp) :: y2 !< |
---|
870 | REAL(wp) :: y3 !< |
---|
871 | |
---|
872 | REAL(wp), DIMENSION(1:11), SAVE :: rat !< |
---|
873 | REAL(wp), DIMENSION(1:7), SAVE :: r0 !< |
---|
874 | |
---|
875 | REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_100 !< |
---|
876 | REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_400 !< |
---|
877 | |
---|
878 | ! |
---|
879 | !-- Initial assignment of constants |
---|
880 | IF ( first ) THEN |
---|
881 | |
---|
882 | first = .FALSE. |
---|
883 | |
---|
884 | r0 = (/ 10.0_wp, 20.0_wp, 30.0_wp, 40.0_wp, 50.0_wp, 60.0_wp, & |
---|
885 | 100.0_wp /) |
---|
886 | |
---|
887 | rat = (/ 0.0_wp, 0.1_wp, 0.2_wp, 0.3_wp, 0.4_wp, 0.5_wp, 0.6_wp, & |
---|
888 | 0.7_wp, 0.8_wp, 0.9_wp, 1.0_wp /) |
---|
889 | ! |
---|
890 | !-- Tabulated turbulent enhancement factor at 100 cm**2/s**3 |
---|
891 | ecoll_100(:,1) = (/ 1.74_wp, 1.74_wp, 1.773_wp, 1.49_wp, & |
---|
892 | 1.207_wp, 1.207_wp, 1.0_wp /) |
---|
893 | ecoll_100(:,2) = (/ 1.46_wp, 1.46_wp, 1.421_wp, 1.245_wp, & |
---|
894 | 1.069_wp, 1.069_wp, 1.0_wp /) |
---|
895 | ecoll_100(:,3) = (/ 1.32_wp, 1.32_wp, 1.245_wp, 1.123_wp, & |
---|
896 | 1.000_wp, 1.000_wp, 1.0_wp /) |
---|
897 | ecoll_100(:,4) = (/ 1.250_wp, 1.250_wp, 1.148_wp, 1.087_wp, & |
---|
898 | 1.025_wp, 1.025_wp, 1.0_wp /) |
---|
899 | ecoll_100(:,5) = (/ 1.186_wp, 1.186_wp, 1.066_wp, 1.060_wp, & |
---|
900 | 1.056_wp, 1.056_wp, 1.0_wp /) |
---|
901 | ecoll_100(:,6) = (/ 1.045_wp, 1.045_wp, 1.000_wp, 1.014_wp, & |
---|
902 | 1.028_wp, 1.028_wp, 1.0_wp /) |
---|
903 | ecoll_100(:,7) = (/ 1.070_wp, 1.070_wp, 1.030_wp, 1.038_wp, & |
---|
904 | 1.046_wp, 1.046_wp, 1.0_wp /) |
---|
905 | ecoll_100(:,8) = (/ 1.000_wp, 1.000_wp, 1.054_wp, 1.042_wp, & |
---|
906 | 1.029_wp, 1.029_wp, 1.0_wp /) |
---|
907 | ecoll_100(:,9) = (/ 1.223_wp, 1.223_wp, 1.117_wp, 1.069_wp, & |
---|
908 | 1.021_wp, 1.021_wp, 1.0_wp /) |
---|
909 | ecoll_100(:,10) = (/ 1.570_wp, 1.570_wp, 1.244_wp, 1.166_wp, & |
---|
910 | 1.088_wp, 1.088_wp, 1.0_wp /) |
---|
911 | ecoll_100(:,11) = (/ 20.3_wp, 20.3_wp, 14.6_wp, 8.61_wp, & |
---|
912 | 2.60_wp, 2.60_wp, 1.0_wp /) |
---|
913 | ! |
---|
914 | !-- Tabulated turbulent enhancement factor at 400 cm**2/s**3 |
---|
915 | ecoll_400(:,1) = (/ 4.976_wp, 4.976_wp, 3.593_wp, 2.519_wp, & |
---|
916 | 1.445_wp, 1.445_wp, 1.0_wp /) |
---|
917 | ecoll_400(:,2) = (/ 2.984_wp, 2.984_wp, 2.181_wp, 1.691_wp, & |
---|
918 | 1.201_wp, 1.201_wp, 1.0_wp /) |
---|
919 | ecoll_400(:,3) = (/ 1.988_wp, 1.988_wp, 1.475_wp, 1.313_wp, & |
---|
920 | 1.150_wp, 1.150_wp, 1.0_wp /) |
---|
921 | ecoll_400(:,4) = (/ 1.490_wp, 1.490_wp, 1.187_wp, 1.156_wp, & |
---|
922 | 1.126_wp, 1.126_wp, 1.0_wp /) |
---|
923 | ecoll_400(:,5) = (/ 1.249_wp, 1.249_wp, 1.088_wp, 1.090_wp, & |
---|
924 | 1.092_wp, 1.092_wp, 1.0_wp /) |
---|
925 | ecoll_400(:,6) = (/ 1.139_wp, 1.139_wp, 1.130_wp, 1.091_wp, & |
---|
926 | 1.051_wp, 1.051_wp, 1.0_wp /) |
---|
927 | ecoll_400(:,7) = (/ 1.220_wp, 1.220_wp, 1.190_wp, 1.138_wp, & |
---|
928 | 1.086_wp, 1.086_wp, 1.0_wp /) |
---|
929 | ecoll_400(:,8) = (/ 1.325_wp, 1.325_wp, 1.267_wp, 1.165_wp, & |
---|
930 | 1.063_wp, 1.063_wp, 1.0_wp /) |
---|
931 | ecoll_400(:,9) = (/ 1.716_wp, 1.716_wp, 1.345_wp, 1.223_wp, & |
---|
932 | 1.100_wp, 1.100_wp, 1.0_wp /) |
---|
933 | ecoll_400(:,10) = (/ 3.788_wp, 3.788_wp, 1.501_wp, 1.311_wp, & |
---|
934 | 1.120_wp, 1.120_wp, 1.0_wp /) |
---|
935 | ecoll_400(:,11) = (/ 36.52_wp, 36.52_wp, 19.16_wp, 22.80_wp, & |
---|
936 | 26.0_wp, 26.0_wp, 1.0_wp /) |
---|
937 | |
---|
938 | ENDIF |
---|
939 | |
---|
940 | ! |
---|
941 | !-- Calculate the radius class index of particles with respect to array r0 |
---|
942 | !-- The droplet radius has to be given in microns. |
---|
943 | ALLOCATE( ira(1:radius_classes) ) |
---|
944 | |
---|
945 | DO j = 1, radius_classes |
---|
946 | particle_radius = radclass(j) * 1.0E6_wp |
---|
947 | DO k = 1, 7 |
---|
948 | IF ( particle_radius < r0(k) ) THEN |
---|
949 | ira(j) = k |
---|
950 | EXIT |
---|
951 | ENDIF |
---|
952 | ENDDO |
---|
953 | IF ( particle_radius >= r0(7) ) ira(j) = 8 |
---|
954 | ENDDO |
---|
955 | |
---|
956 | ! |
---|
957 | !-- Two-dimensional linear interpolation of the turbulent enhancement factor. |
---|
958 | !-- The droplet radius has to be given in microns. |
---|
959 | DO j = 1, radius_classes |
---|
960 | DO i = 1, j |
---|
961 | |
---|
962 | ir = MAX( ira(i), ira(j) ) |
---|
963 | rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) ) |
---|
964 | |
---|
965 | DO kk = 2, 11 |
---|
966 | IF ( rq <= rat(kk) ) THEN |
---|
967 | iq = kk |
---|
968 | EXIT |
---|
969 | ENDIF |
---|
970 | ENDDO |
---|
971 | |
---|
972 | y1 = 1.0_wp ! turbulent enhancement factor at 0 m**2/s**3 |
---|
973 | |
---|
974 | IF ( ir < 8 ) THEN |
---|
975 | IF ( ir >= 2 ) THEN |
---|
976 | pp = ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp - & |
---|
977 | r0(ir-1) ) / ( r0(ir) - r0(ir-1) ) |
---|
978 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
979 | y2 = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) * ecoll_100(ir-1,iq-1) + & |
---|
980 | pp * ( 1.0_wp - qq ) * ecoll_100(ir,iq-1) + & |
---|
981 | qq * ( 1.0_wp - pp ) * ecoll_100(ir-1,iq) + & |
---|
982 | pp * qq * ecoll_100(ir,iq) |
---|
983 | y3 = ( 1.0-pp ) * ( 1.0_wp - qq ) * ecoll_400(ir-1,iq-1) + & |
---|
984 | pp * ( 1.0_wp - qq ) * ecoll_400(ir,iq-1) + & |
---|
985 | qq * ( 1.0_wp - pp ) * ecoll_400(ir-1,iq) + & |
---|
986 | pp * qq * ecoll_400(ir,iq) |
---|
987 | ELSE |
---|
988 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
989 | y2 = ( 1.0_wp - qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq) |
---|
990 | y3 = ( 1.0_wp - qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq) |
---|
991 | ENDIF |
---|
992 | ELSE |
---|
993 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
994 | y2 = ( 1.0_wp - qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq) |
---|
995 | y3 = ( 1.0_wp - qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq) |
---|
996 | ENDIF |
---|
997 | ! |
---|
998 | !-- Linear interpolation of turbulent enhancement factor |
---|
999 | IF ( epsilon <= 0.01_wp ) THEN |
---|
1000 | ecf(j,i) = ( epsilon - 0.01_wp ) / ( 0.0_wp - 0.01_wp ) * y1 & |
---|
1001 | + ( epsilon - 0.0_wp ) / ( 0.01_wp - 0.0_wp ) * y2 |
---|
1002 | ELSEIF ( epsilon <= 0.06_wp ) THEN |
---|
1003 | ecf(j,i) = ( epsilon - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 & |
---|
1004 | + ( epsilon - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3 |
---|
1005 | ELSE |
---|
1006 | ecf(j,i) = ( 0.06_wp - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 & |
---|
1007 | + ( 0.06_wp - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3 |
---|
1008 | ENDIF |
---|
1009 | |
---|
1010 | IF ( ecf(j,i) < 1.0_wp ) ecf(j,i) = 1.0_wp |
---|
1011 | |
---|
1012 | ecf(i,j) = ecf(j,i) |
---|
1013 | |
---|
1014 | ENDDO |
---|
1015 | ENDDO |
---|
1016 | |
---|
1017 | END SUBROUTINE turb_enhance_eff |
---|
1018 | |
---|
1019 | END MODULE lpm_collision_kernels_mod |
---|