[1682] | 1 | !> @file lpm_droplet_collision.f90 |
---|
[2000] | 2 | !------------------------------------------------------------------------------! |
---|
[2696] | 3 | ! This file is part of the PALM model system. |
---|
[1036] | 4 | ! |
---|
[2000] | 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. |
---|
[1036] | 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 | ! |
---|
[2101] | 17 | ! Copyright 1997-2017 Leibniz Universitaet Hannover |
---|
[2000] | 18 | !------------------------------------------------------------------------------! |
---|
[1036] | 19 | ! |
---|
[849] | 20 | ! Current revisions: |
---|
| 21 | ! ------------------ |
---|
[2375] | 22 | ! |
---|
| 23 | ! |
---|
[1321] | 24 | ! Former revisions: |
---|
| 25 | ! ----------------- |
---|
| 26 | ! $Id: lpm_droplet_collision.f90 2696 2017-12-14 17:12:51Z kanani $ |
---|
[2375] | 27 | ! Changed ONLY-dependencies |
---|
| 28 | ! |
---|
| 29 | ! 2312 2017-07-14 20:26:51Z hoffmann |
---|
[2312] | 30 | ! Consideration of aerosol mass during collision. Average impact algorithm has |
---|
| 31 | ! been removed. |
---|
| 32 | ! |
---|
| 33 | ! 2274 2017-06-09 13:27:48Z Giersch |
---|
[2274] | 34 | ! Changed error messages |
---|
[1321] | 35 | ! |
---|
[2274] | 36 | ! 2123 2017-01-18 12:34:59Z hoffmann |
---|
| 37 | ! |
---|
[2123] | 38 | ! 2122 2017-01-18 12:22:54Z hoffmann |
---|
| 39 | ! Some reformatting of the code. |
---|
| 40 | ! |
---|
[2001] | 41 | ! 2000 2016-08-20 18:09:15Z knoop |
---|
| 42 | ! Forced header and separation lines into 80 columns |
---|
| 43 | ! |
---|
[1885] | 44 | ! 1884 2016-04-21 11:11:40Z hoffmann |
---|
| 45 | ! Conservation of mass should only be checked if collisions took place. |
---|
| 46 | ! |
---|
[1861] | 47 | ! 1860 2016-04-13 13:21:28Z hoffmann |
---|
[2312] | 48 | ! Interpolation of dissipation rate adjusted to more reasonable values. |
---|
[1861] | 49 | ! |
---|
[1823] | 50 | ! 1822 2016-04-07 07:49:42Z hoffmann |
---|
| 51 | ! Integration of a new collision algortithm based on Shima et al. (2009) and |
---|
| 52 | ! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented |
---|
| 53 | ! collision algorithm is called average_impact. Moreover, both algorithms are |
---|
| 54 | ! now positive definit due to their construction, i.e., no negative weighting |
---|
| 55 | ! factors should occur. |
---|
| 56 | ! |
---|
[1683] | 57 | ! 1682 2015-10-07 23:56:08Z knoop |
---|
[2312] | 58 | ! Code annotations made doxygen readable |
---|
| 59 | ! |
---|
[1360] | 60 | ! 1359 2014-04-11 17:15:14Z hoffmann |
---|
[2312] | 61 | ! New particle structure integrated. |
---|
[1360] | 62 | ! Kind definition added to all floating point numbers. |
---|
| 63 | ! |
---|
[1323] | 64 | ! 1322 2014-03-20 16:38:49Z raasch |
---|
| 65 | ! REAL constants defined as wp_kind |
---|
| 66 | ! |
---|
[1321] | 67 | ! 1320 2014-03-20 08:40:49Z raasch |
---|
[1320] | 68 | ! ONLY-attribute added to USE-statements, |
---|
| 69 | ! kind-parameters added to all INTEGER and REAL declaration statements, |
---|
| 70 | ! kinds are defined in new module kinds, |
---|
| 71 | ! revision history before 2012 removed, |
---|
| 72 | ! comment fields (!:) to be used for variable explanations added to |
---|
| 73 | ! all variable declaration statements |
---|
[1072] | 74 | ! |
---|
[1093] | 75 | ! 1092 2013-02-02 11:24:22Z raasch |
---|
| 76 | ! unused variables removed |
---|
| 77 | ! |
---|
[1072] | 78 | ! 1071 2012-11-29 16:54:55Z franke |
---|
[1071] | 79 | ! Calculation of Hall and Wang kernel now uses collision-coalescence formulation |
---|
| 80 | ! proposed by Wang instead of the continuous collection equation (for more |
---|
| 81 | ! information about new method see PALM documentation) |
---|
| 82 | ! Bugfix: message identifiers added |
---|
[849] | 83 | ! |
---|
[1037] | 84 | ! 1036 2012-10-22 13:43:42Z raasch |
---|
| 85 | ! code put under GPL (PALM 3.9) |
---|
| 86 | ! |
---|
[850] | 87 | ! 849 2012-03-15 10:35:09Z raasch |
---|
| 88 | ! initial revision (former part of advec_particles) |
---|
[849] | 89 | ! |
---|
[850] | 90 | ! |
---|
[849] | 91 | ! Description: |
---|
| 92 | ! ------------ |
---|
[1682] | 93 | !> Calculates change in droplet radius by collision. Droplet collision is |
---|
| 94 | !> calculated for each grid box seperately. Collision is parameterized by |
---|
[1822] | 95 | !> using collision kernels. Two different kernels are available: |
---|
[1682] | 96 | !> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which |
---|
| 97 | !> considers collision due to pure gravitational effects. |
---|
| 98 | !> Wang kernel: Beside gravitational effects (treated with the Hall-kernel) also |
---|
| 99 | !> the effects of turbulence on the collision are considered using |
---|
| 100 | !> parameterizations of Ayala et al. (2008, New J. Phys., 10, |
---|
| 101 | !> 075015) and Wang and Grabowski (2009, Atmos. Sci. Lett., 10, |
---|
| 102 | !> 1-8). This kernel includes three possible effects of turbulence: |
---|
| 103 | !> the modification of the relative velocity between the droplets, |
---|
| 104 | !> the effect of preferential concentration, and the enhancement of |
---|
[2312] | 105 | !> collision efficiencies. |
---|
[849] | 106 | !------------------------------------------------------------------------------! |
---|
[1682] | 107 | SUBROUTINE lpm_droplet_collision (i,j,k) |
---|
[849] | 108 | |
---|
[1320] | 109 | USE arrays_3d, & |
---|
[1822] | 110 | ONLY: diss, ql_v, ql_vp |
---|
[849] | 111 | |
---|
[1320] | 112 | USE cloud_parameters, & |
---|
[2375] | 113 | ONLY: rho_l, rho_s |
---|
[1320] | 114 | |
---|
| 115 | USE constants, & |
---|
| 116 | ONLY: pi |
---|
| 117 | |
---|
| 118 | USE control_parameters, & |
---|
[1822] | 119 | ONLY: dt_3d, message_string, dz |
---|
[1320] | 120 | |
---|
| 121 | USE cpulog, & |
---|
| 122 | ONLY: cpu_log, log_point_s |
---|
| 123 | |
---|
| 124 | USE grid_variables, & |
---|
[1822] | 125 | ONLY: dx, dy |
---|
[1320] | 126 | |
---|
| 127 | USE kinds |
---|
| 128 | |
---|
| 129 | USE lpm_collision_kernels_mod, & |
---|
[1822] | 130 | ONLY: ckernel, recalculate_kernel |
---|
[1320] | 131 | |
---|
| 132 | USE particle_attributes, & |
---|
[2312] | 133 | ONLY: curvature_solution_effects, dissipation_classes, hall_kernel, & |
---|
| 134 | iran_part, number_of_particles, particles, particle_type, & |
---|
[2375] | 135 | prt_count, use_kernel_tables, wang_kernel |
---|
[1320] | 136 | |
---|
[1822] | 137 | USE random_function_mod, & |
---|
| 138 | ONLY: random_function |
---|
| 139 | |
---|
[1359] | 140 | USE pegrid |
---|
| 141 | |
---|
[849] | 142 | IMPLICIT NONE |
---|
| 143 | |
---|
[1682] | 144 | INTEGER(iwp) :: eclass !< |
---|
| 145 | INTEGER(iwp) :: i !< |
---|
| 146 | INTEGER(iwp) :: j !< |
---|
| 147 | INTEGER(iwp) :: k !< |
---|
| 148 | INTEGER(iwp) :: n !< |
---|
[1822] | 149 | INTEGER(iwp) :: m !< |
---|
[1682] | 150 | INTEGER(iwp) :: rclass_l !< |
---|
| 151 | INTEGER(iwp) :: rclass_s !< |
---|
[849] | 152 | |
---|
[1822] | 153 | REAL(wp) :: collection_probability !< probability for collection |
---|
| 154 | REAL(wp) :: ddV !< inverse grid box volume |
---|
| 155 | REAL(wp) :: epsilon !< dissipation rate |
---|
| 156 | REAL(wp) :: factor_volume_to_mass !< 4.0 / 3.0 * pi * rho_l |
---|
[2312] | 157 | REAL(wp) :: xm !< droplet mass of super-droplet m |
---|
| 158 | REAL(wp) :: xn !< droplet mass of super-droplet n |
---|
| 159 | REAL(wp) :: xsm !< aerosol mass of super-droplet m |
---|
| 160 | REAL(wp) :: xsn !< aerosol mass of super-droplet n |
---|
[1359] | 161 | |
---|
[2312] | 162 | REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !< weighting factor |
---|
| 163 | REAL(wp), DIMENSION(:), ALLOCATABLE :: mass !< total mass of super droplet |
---|
| 164 | REAL(wp), DIMENSION(:), ALLOCATABLE :: aero_mass !< total aerosol mass of super droplet |
---|
[1359] | 165 | |
---|
[849] | 166 | CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' ) |
---|
| 167 | |
---|
[1822] | 168 | number_of_particles = prt_count(k,j,i) |
---|
[2312] | 169 | factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l |
---|
| 170 | ddV = 1.0_wp / ( dx * dy * dz ) |
---|
[849] | 171 | ! |
---|
[1822] | 172 | !-- Collision requires at least one super droplet inside the box |
---|
| 173 | IF ( number_of_particles > 0 ) THEN |
---|
[849] | 174 | |
---|
[1822] | 175 | IF ( use_kernel_tables ) THEN |
---|
[849] | 176 | ! |
---|
[1822] | 177 | !-- Fast method with pre-calculated collection kernels for |
---|
[1359] | 178 | !-- discrete radius- and dissipation-classes. |
---|
| 179 | IF ( wang_kernel ) THEN |
---|
[1860] | 180 | eclass = INT( diss(k,j,i) * 1.0E4_wp / 600.0_wp * & |
---|
[1359] | 181 | dissipation_classes ) + 1 |
---|
| 182 | epsilon = diss(k,j,i) |
---|
| 183 | ELSE |
---|
| 184 | epsilon = 0.0_wp |
---|
| 185 | ENDIF |
---|
[2312] | 186 | |
---|
[1359] | 187 | IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001_wp ) THEN |
---|
| 188 | eclass = 0 ! Hall kernel is used |
---|
| 189 | ELSE |
---|
| 190 | eclass = MIN( dissipation_classes, eclass ) |
---|
| 191 | ENDIF |
---|
| 192 | |
---|
[2312] | 193 | ELSE |
---|
[849] | 194 | ! |
---|
[2312] | 195 | !-- Collection kernels are re-calculated for every new |
---|
[1822] | 196 | !-- grid box. First, allocate memory for kernel table. |
---|
| 197 | !-- Third dimension is 1, because table is re-calculated for |
---|
| 198 | !-- every new dissipation value. |
---|
| 199 | ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) ) |
---|
| 200 | ! |
---|
| 201 | !-- Now calculate collection kernel for this box. Note that |
---|
| 202 | !-- the kernel is based on the previous time step |
---|
| 203 | CALL recalculate_kernel( i, j, k ) |
---|
[2312] | 204 | |
---|
| 205 | ENDIF |
---|
[1822] | 206 | ! |
---|
[2312] | 207 | !-- Temporary fields for total mass of super-droplet, aerosol mass, and |
---|
| 208 | !-- weighting factor are allocated. |
---|
| 209 | ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles)) |
---|
| 210 | IF ( curvature_solution_effects ) ALLOCATE(aero_mass(1:number_of_particles)) |
---|
[849] | 211 | |
---|
[2312] | 212 | mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * & |
---|
| 213 | particles(1:number_of_particles)%radius**3 * & |
---|
| 214 | factor_volume_to_mass |
---|
[849] | 215 | |
---|
[2312] | 216 | weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor |
---|
[849] | 217 | |
---|
[2312] | 218 | IF ( curvature_solution_effects ) THEN |
---|
| 219 | aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * & |
---|
| 220 | particles(1:number_of_particles)%aux1**3 * & |
---|
| 221 | 4.0 / 3.0 * pi * rho_s |
---|
| 222 | ENDIF |
---|
| 223 | ! |
---|
| 224 | !-- Calculate collision/coalescence |
---|
| 225 | DO n = 1, number_of_particles |
---|
[849] | 226 | |
---|
[2312] | 227 | DO m = n, number_of_particles |
---|
| 228 | ! |
---|
| 229 | !-- For collisions, the weighting factor of at least one super-droplet |
---|
| 230 | !-- needs to be larger or equal to one. |
---|
| 231 | IF ( MIN( weight(n), weight(m) ) .LT. 1.0 ) CYCLE |
---|
| 232 | ! |
---|
| 233 | !-- Get mass of individual droplets (aerosols) |
---|
| 234 | xn = mass(n) / weight(n) |
---|
| 235 | xm = mass(m) / weight(m) |
---|
| 236 | IF ( curvature_solution_effects ) THEN |
---|
| 237 | xsn = aero_mass(n) / weight(n) |
---|
| 238 | xsm = aero_mass(m) / weight(m) |
---|
| 239 | ENDIF |
---|
| 240 | ! |
---|
| 241 | !-- Probability that the necessary collisions take place |
---|
| 242 | IF ( use_kernel_tables ) THEN |
---|
| 243 | rclass_l = particles(n)%class |
---|
| 244 | rclass_s = particles(m)%class |
---|
[849] | 245 | |
---|
[2312] | 246 | collection_probability = MAX( weight(n), weight(m) ) * & |
---|
| 247 | ckernel(rclass_l,rclass_s,eclass) * ddV * dt_3d |
---|
| 248 | ELSE |
---|
| 249 | collection_probability = MAX( weight(n), weight(m) ) * & |
---|
| 250 | ckernel(n,m,1) * ddV * dt_3d |
---|
| 251 | ENDIF |
---|
[849] | 252 | ! |
---|
[2312] | 253 | !-- Calculate the number of collections and consider multiple collections. |
---|
| 254 | !-- (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...) |
---|
| 255 | IF ( collection_probability - FLOOR(collection_probability) & |
---|
| 256 | .GT. random_function( iran_part ) ) THEN |
---|
| 257 | collection_probability = FLOOR(collection_probability) + 1.0_wp |
---|
| 258 | ELSE |
---|
| 259 | collection_probability = FLOOR(collection_probability) |
---|
| 260 | ENDIF |
---|
[1822] | 261 | |
---|
[2312] | 262 | IF ( collection_probability .GT. 0.0 ) THEN |
---|
[849] | 263 | ! |
---|
[2312] | 264 | !-- Super-droplet n collects droplets of super-droplet m |
---|
| 265 | IF ( weight(n) .LT. weight(m) ) THEN |
---|
[849] | 266 | |
---|
[2312] | 267 | mass(n) = mass(n) + weight(n) * xm * collection_probability |
---|
| 268 | weight(m) = weight(m) - weight(n) * collection_probability |
---|
| 269 | mass(m) = mass(m) - weight(n) * xm * collection_probability |
---|
| 270 | IF ( curvature_solution_effects ) THEN |
---|
| 271 | aero_mass(n) = aero_mass(n) + weight(n) * xsm * collection_probability |
---|
| 272 | aero_mass(m) = aero_mass(m) - weight(n) * xsm * collection_probability |
---|
[1822] | 273 | ENDIF |
---|
[849] | 274 | |
---|
[2312] | 275 | ELSEIF ( weight(m) .LT. weight(n) ) THEN |
---|
[849] | 276 | |
---|
[2312] | 277 | mass(m) = mass(m) + weight(m) * xn * collection_probability |
---|
| 278 | weight(n) = weight(n) - weight(m) * collection_probability |
---|
| 279 | mass(n) = mass(n) - weight(m) * xn * collection_probability |
---|
| 280 | IF ( curvature_solution_effects ) THEN |
---|
| 281 | aero_mass(m) = aero_mass(m) + weight(m) * xsn * collection_probability |
---|
| 282 | aero_mass(n) = aero_mass(n) - weight(m) * xsn * collection_probability |
---|
| 283 | ENDIF |
---|
[849] | 284 | |
---|
[2312] | 285 | ELSE |
---|
[1822] | 286 | ! |
---|
[2312] | 287 | !-- Collisions of particles of the same weighting factor. |
---|
| 288 | !-- Particle n collects 1/2 weight(n) droplets of particle m, |
---|
| 289 | !-- particle m collects 1/2 weight(m) droplets of particle n. |
---|
| 290 | !-- The total mass mass changes accordingly. |
---|
| 291 | !-- If n = m, the first half of the droplets coalesces with the |
---|
| 292 | !-- second half of the droplets; mass is unchanged because |
---|
| 293 | !-- xm = xn for n = m. |
---|
[1822] | 294 | !-- |
---|
[2312] | 295 | !-- Note: For m = n this equation is an approximation only |
---|
| 296 | !-- valid for weight >> 1 (which is usually the case). The |
---|
| 297 | !-- approximation is weight(n)-1 = weight(n). |
---|
| 298 | mass(n) = mass(n) + 0.5_wp * weight(n) * ( xm - xn ) |
---|
| 299 | mass(m) = mass(m) + 0.5_wp * weight(m) * ( xn - xm ) |
---|
| 300 | IF ( curvature_solution_effects ) THEN |
---|
| 301 | aero_mass(n) = aero_mass(n) + 0.5_wp * weight(n) * ( xsm - xsn ) |
---|
| 302 | aero_mass(m) = aero_mass(m) + 0.5_wp * weight(m) * ( xsn - xsm ) |
---|
[1822] | 303 | ENDIF |
---|
[2312] | 304 | weight(n) = weight(n) - 0.5_wp * weight(m) |
---|
| 305 | weight(m) = weight(n) |
---|
[1071] | 306 | |
---|
[2312] | 307 | ENDIF |
---|
[1071] | 308 | |
---|
[2312] | 309 | ENDIF |
---|
[1071] | 310 | |
---|
[2312] | 311 | ENDDO |
---|
[1071] | 312 | |
---|
[2312] | 313 | ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass |
---|
[1071] | 314 | |
---|
[2312] | 315 | ENDDO |
---|
[1071] | 316 | |
---|
[2312] | 317 | IF ( ANY(weight < 0.0_wp) ) THEN |
---|
| 318 | WRITE( message_string, * ) 'negative weighting factor' |
---|
| 319 | CALL message( 'lpm_droplet_collision', 'PA0028', & |
---|
| 320 | 2, 2, -1, 6, 1 ) |
---|
| 321 | ENDIF |
---|
[1822] | 322 | |
---|
[2312] | 323 | particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) / & |
---|
| 324 | ( weight(1:number_of_particles) & |
---|
| 325 | * factor_volume_to_mass & |
---|
| 326 | ) & |
---|
| 327 | )**0.33333333333333_wp |
---|
[1822] | 328 | |
---|
[2312] | 329 | IF ( curvature_solution_effects ) THEN |
---|
| 330 | particles(1:number_of_particles)%aux1 = ( aero_mass(1:number_of_particles) / & |
---|
| 331 | ( weight(1:number_of_particles) & |
---|
| 332 | * 4.0_wp / 3.0_wp * pi * rho_s & |
---|
| 333 | ) & |
---|
| 334 | )**0.33333333333333_wp |
---|
| 335 | ENDIF |
---|
[1822] | 336 | |
---|
[2312] | 337 | particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles) |
---|
[1822] | 338 | |
---|
[2312] | 339 | DEALLOCATE( weight, mass ) |
---|
| 340 | IF ( curvature_solution_effects ) DEALLOCATE( aero_mass ) |
---|
| 341 | IF ( .NOT. use_kernel_tables ) DEALLOCATE( ckernel ) |
---|
| 342 | |
---|
[849] | 343 | ! |
---|
[1884] | 344 | !-- Check if LWC is conserved during collision process |
---|
| 345 | IF ( ql_v(k,j,i) /= 0.0_wp ) THEN |
---|
| 346 | IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp .OR. & |
---|
| 347 | ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp ) THEN |
---|
| 348 | WRITE( message_string, * ) ' LWC is not conserved during', & |
---|
| 349 | ' collision! ', & |
---|
| 350 | ' LWC after condensation: ', ql_v(k,j,i), & |
---|
| 351 | ' LWC after collision: ', ql_vp(k,j,i) |
---|
| 352 | CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 ) |
---|
| 353 | ENDIF |
---|
[1359] | 354 | ENDIF |
---|
[1884] | 355 | |
---|
[1359] | 356 | ENDIF |
---|
[2312] | 357 | |
---|
[849] | 358 | CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' ) |
---|
| 359 | |
---|
[2312] | 360 | END SUBROUTINE lpm_droplet_collision |
---|