!> @file lpm_droplet_collision.f90 !--------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the terms ! of the GNU General Public License as published by the Free Software Foundation, ! either version 3 of the License, or (at your option) any later version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2016 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: lpm_droplet_collision.f90 1885 2016-04-21 11:12:26Z suehring $ ! ! 1884 2016-04-21 11:11:40Z hoffmann ! Conservation of mass should only be checked if collisions took place. ! ! 1860 2016-04-13 13:21:28Z hoffmann ! Interpolation of dissipation rate adjusted to more reasonable values. ! ! 1822 2016-04-07 07:49:42Z hoffmann ! Integration of a new collision algortithm based on Shima et al. (2009) and ! Soelch and Kaercher (2010) called all_or_nothing. The previous implemented ! collision algorithm is called average_impact. Moreover, both algorithms are ! now positive definit due to their construction, i.e., no negative weighting ! factors should occur. ! ! 1682 2015-10-07 23:56:08Z knoop ! Code annotations made doxygen readable ! ! 1359 2014-04-11 17:15:14Z hoffmann ! New particle structure integrated. ! Kind definition added to all floating point numbers. ! ! 1322 2014-03-20 16:38:49Z raasch ! REAL constants defined as wp_kind ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! revision history before 2012 removed, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1092 2013-02-02 11:24:22Z raasch ! unused variables removed ! ! 1071 2012-11-29 16:54:55Z franke ! Calculation of Hall and Wang kernel now uses collision-coalescence formulation ! proposed by Wang instead of the continuous collection equation (for more ! information about new method see PALM documentation) ! Bugfix: message identifiers added ! ! 1036 2012-10-22 13:43:42Z raasch ! code put under GPL (PALM 3.9) ! ! 849 2012-03-15 10:35:09Z raasch ! initial revision (former part of advec_particles) ! ! ! Description: ! ------------ !> Calculates change in droplet radius by collision. Droplet collision is !> calculated for each grid box seperately. Collision is parameterized by !> using collision kernels. Two different kernels are available: !> Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which !> considers collision due to pure gravitational effects. !> Wang kernel: Beside gravitational effects (treated with the Hall-kernel) also !> the effects of turbulence on the collision are considered using !> parameterizations of Ayala et al. (2008, New J. Phys., 10, !> 075015) and Wang and Grabowski (2009, Atmos. Sci. Lett., 10, !> 1-8). This kernel includes three possible effects of turbulence: !> the modification of the relative velocity between the droplets, !> the effect of preferential concentration, and the enhancement of !> collision efficiencies. !------------------------------------------------------------------------------! SUBROUTINE lpm_droplet_collision (i,j,k) USE arrays_3d, & ONLY: diss, ql_v, ql_vp USE cloud_parameters, & ONLY: rho_l USE constants, & ONLY: pi USE control_parameters, & ONLY: dt_3d, message_string, dz USE cpulog, & ONLY: cpu_log, log_point_s USE grid_variables, & ONLY: dx, dy USE kinds USE lpm_collision_kernels_mod, & ONLY: ckernel, recalculate_kernel USE particle_attributes, & ONLY: all_or_nothing, average_impact, dissipation_classes, & hall_kernel, iran_part, number_of_particles, particles, & particle_type, prt_count, use_kernel_tables, wang_kernel USE random_function_mod, & ONLY: random_function USE pegrid IMPLICIT NONE INTEGER(iwp) :: eclass !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: n !< INTEGER(iwp) :: m !< INTEGER(iwp) :: rclass_l !< INTEGER(iwp) :: rclass_s !< REAL(wp) :: collection_probability !< probability for collection REAL(wp) :: ddV !< inverse grid box volume REAL(wp) :: epsilon !< dissipation rate REAL(wp) :: factor_volume_to_mass !< 4.0 / 3.0 * pi * rho_l REAL(wp) :: xm !< mean mass of droplet m REAL(wp) :: xn !< mean mass of droplet n REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !< weighting factor REAL(wp), DIMENSION(:), ALLOCATABLE :: mass !< total mass of super droplet CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' ) number_of_particles = prt_count(k,j,i) factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l ddV = 1 / ( dx * dy * dz ) ! !-- Collision requires at least one super droplet inside the box IF ( number_of_particles > 0 ) THEN ! !-- Now apply the different kernels IF ( use_kernel_tables ) THEN ! !-- Fast method with pre-calculated collection kernels for !-- discrete radius- and dissipation-classes. !-- !-- Determine dissipation class index of this gridbox IF ( wang_kernel ) THEN eclass = INT( diss(k,j,i) * 1.0E4_wp / 600.0_wp * & dissipation_classes ) + 1 epsilon = diss(k,j,i) ELSE epsilon = 0.0_wp ENDIF IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001_wp ) THEN eclass = 0 ! Hall kernel is used ELSE eclass = MIN( dissipation_classes, eclass ) ENDIF ! !-- Droplet collision are calculated using collision-coalescence !-- formulation proposed by Wang (see PALM documentation) !-- Temporary fields for total mass of super-droplet and weighting factors !-- are allocated. ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles)) mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * & particles(1:number_of_particles)%radius**3 * & factor_volume_to_mass weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor IF ( average_impact ) THEN ! select collision algorithm DO n = 1, number_of_particles rclass_l = particles(n)%class xn = mass(n) / weight(n) DO m = n, number_of_particles rclass_s = particles(m)%class xm = mass(m) / weight(m) IF ( xm .LT. xn ) THEN ! !-- Particle n collects smaller particle m collection_probability = ckernel(rclass_l,rclass_s,eclass) * & weight(n) * ddV * dt_3d mass(n) = mass(n) + mass(m) * collection_probability weight(m) = weight(m) - weight(m) * collection_probability mass(m) = mass(m) - mass(m) * collection_probability ELSEIF ( xm .GT. xn ) THEN ! !-- Particle m collects smaller particle n collection_probability = ckernel(rclass_l,rclass_s,eclass) * & weight(m) * ddV * dt_3d mass(m) = mass(m) + mass(n) * collection_probability weight(n) = weight(n) - weight(n) * collection_probability mass(n) = mass(n) - mass(n) * collection_probability ELSE ! !-- Same-size collections. If n = m, weight is reduced by the !-- number of possible same-size collections; the total mass !-- is not changed during same-size collection. !-- Same-size collections of different !-- particles ( n /= m ) are treated as same-size collections !-- of ONE partilce with weight = weight(n) + weight(m) and !-- mass = mass(n) + mass(m). !-- Accordingly, each particle loses the same number of !-- droplets to the other particle, but this has no effect on !-- total mass mass, since the exchanged droplets have the !-- same radius. !-- Note: For m = n this equation is an approximation only !-- valid for weight >> 1 (which is usually the case). The !-- approximation is weight(n)-1 = weight(n). weight(n) = weight(n) - 0.5_wp * weight(n) * & ckernel(rclass_l,rclass_s,eclass) * & weight(m) * ddV * dt_3d IF ( n .NE. m ) THEN weight(m) = weight(m) - 0.5_wp * weight(m) * & ckernel(rclass_l,rclass_s,eclass) * & weight(n) * ddV * dt_3d ENDIF ENDIF ENDDO ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass ENDDO ELSEIF ( all_or_nothing ) THEN ! select collision algorithm DO n = 1, number_of_particles rclass_l = particles(n)%class xn = mass(n) / weight(n) ! mean mass of droplet n DO m = n, number_of_particles rclass_s = particles(m)%class xm = mass(m) / weight(m) ! mean mass of droplet m IF ( weight(n) .LT. weight(m) ) THEN ! !-- Particle n collects weight(n) droplets of particle m collection_probability = ckernel(rclass_l,rclass_s,eclass) * & weight(m) * ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(n) = mass(n) + weight(n) * xm weight(m) = weight(m) - weight(n) mass(m) = mass(m) - weight(n) * xm ENDIF ELSEIF ( weight(m) .LT. weight(n) ) THEN ! !-- Particle m collects weight(m) droplets of particle n collection_probability = ckernel(rclass_l,rclass_s,eclass) * & weight(n) * ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(m) = mass(m) + weight(m) * xn weight(n) = weight(n) - weight(m) mass(n) = mass(n) - weight(m) * xn ENDIF ELSE ! !-- Collisions of particles of the same weighting factor. !-- Particle n collects 1/2 weight(n) droplets of particle m, !-- particle m collects 1/2 weight(m) droplets of particle n. !-- The total mass mass changes accordingly. !-- If n = m, the first half of the droplets coalesces with the !-- second half of the droplets; mass is unchanged because !-- xm = xn for n = m. !-- Note: For m = n this equation is an approximation only !-- valid for weight >> 1 (which is usually the case). The !-- approximation is weight(n)-1 = weight(n). collection_probability = ckernel(rclass_l,rclass_s,eclass) * & weight(n) * ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(n) = mass(n) + 0.5_wp * weight(n) * ( xm - xn ) mass(m) = mass(m) + 0.5_wp * weight(m) * ( xn - xm ) weight(n) = weight(n) - 0.5_wp * weight(m) weight(m) = weight(n) ENDIF ENDIF ENDDO ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass ENDDO ENDIF IF ( ANY(weight < 0.0_wp) ) THEN WRITE( message_string, * ) 'negative weighting' CALL message( 'lpm_droplet_collision', 'PA0028', & 2, 2, -1, 6, 1 ) ENDIF particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) / & ( weight(1:number_of_particles) & * factor_volume_to_mass & ) & )**0.33333333333333_wp particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles) DEALLOCATE(weight, mass) ELSEIF ( .NOT. use_kernel_tables ) THEN ! !-- Collection kernels are calculated for every new !-- grid box. First, allocate memory for kernel table. !-- Third dimension is 1, because table is re-calculated for !-- every new dissipation value. ALLOCATE( ckernel(1:number_of_particles,1:number_of_particles,1:1) ) ! !-- Now calculate collection kernel for this box. Note that !-- the kernel is based on the previous time step CALL recalculate_kernel( i, j, k ) ! !-- Droplet collision are calculated using collision-coalescence !-- formulation proposed by Wang (see PALM documentation) !-- Temporary fields for total mass of super-droplet and weighting factors !-- are allocated. ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles)) mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * & particles(1:number_of_particles)%radius**3 * & factor_volume_to_mass weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor IF ( average_impact ) THEN ! select collision algorithm DO n = 1, number_of_particles xn = mass(n) / weight(n) ! mean mass of droplet n DO m = n, number_of_particles xm = mass(m) / weight(m) !mean mass of droplet m IF ( xm .LT. xn ) THEN ! !-- Particle n collects smaller particle m collection_probability = ckernel(n,m,1) * weight(n) * & ddV * dt_3d mass(n) = mass(n) + mass(m) * collection_probability weight(m) = weight(m) - weight(m) * collection_probability mass(m) = mass(m) - mass(m) * collection_probability ELSEIF ( xm .GT. xn ) THEN ! !-- Particle m collects smaller particle n collection_probability = ckernel(n,m,1) * weight(m) * & ddV * dt_3d mass(m) = mass(m) + mass(n) * collection_probability weight(n) = weight(n) - weight(n) * collection_probability mass(n) = mass(n) - mass(n) * collection_probability ELSE ! !-- Same-size collections. If n = m, weight is reduced by the !-- number of possible same-size collections; the total mass !-- mass is not changed during same-size collection. !-- Same-size collections of different !-- particles ( n /= m ) are treated as same-size collections !-- of ONE partilce with weight = weight(n) + weight(m) and !-- mass = mass(n) + mass(m). !-- Accordingly, each particle loses the same number of !-- droplets to the other particle, but this has no effect on !-- total mass mass, since the exchanged droplets have the !-- same radius. !-- !-- Note: For m = n this equation is an approximation only !-- valid for weight >> 1 (which is usually the case). The !-- approximation is weight(n)-1 = weight(n). weight(n) = weight(n) - 0.5_wp * weight(n) * & ckernel(n,m,1) * weight(m) * & ddV * dt_3d IF ( n .NE. m ) THEN weight(m) = weight(m) - 0.5_wp * weight(m) * & ckernel(n,m,1) * weight(n) * & ddV * dt_3d ENDIF ENDIF ENDDO ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass ENDDO ELSEIF ( all_or_nothing ) THEN ! select collision algorithm DO n = 1, number_of_particles xn = mass(n) / weight(n) ! mean mass of droplet n DO m = n, number_of_particles xm = mass(m) / weight(m) !mean mass of droplet m IF ( weight(n) .LT. weight(m) ) THEN ! !-- Particle n collects smaller particle m collection_probability = ckernel(n,m,1) * weight(m) * & ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(n) = mass(n) + weight(n) * xm weight(m) = weight(m) - weight(n) mass(m) = mass(m) - weight(n) * xm ENDIF ELSEIF ( weight(m) .LT. weight(n) ) THEN ! !-- Particle m collects smaller particle n collection_probability = ckernel(n,m,1) * weight(n) * & ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(m) = mass(m) + weight(m) * xn weight(n) = weight(n) - weight(m) mass(n) = mass(n) - weight(m) * xn ENDIF ELSE ! !-- Collisions of particles of the same weighting factor. !-- Particle n collects 1/2 weight(n) droplets of particle m, !-- particle m collects 1/2 weight(m) droplets of particle n. !-- The total mass mass changes accordingly. !-- If n = m, the first half of the droplets coalesces with the !-- second half of the droplets; mass is unchanged because !-- xm = xn for n = m. !-- !-- Note: For m = n this equation is an approximation only !-- valid for weight >> 1 (which is usually the case). The !-- approximation is weight(n)-1 = weight(n). collection_probability = ckernel(n,m,1) * weight(n) * & ddV * dt_3d IF ( collection_probability .GT. random_function( iran_part ) ) THEN mass(n) = mass(n) + 0.5_wp * weight(n) * ( xm - xn ) mass(m) = mass(m) + 0.5_wp * weight(m) * ( xn - xm ) weight(n) = weight(n) - 0.5_wp * weight(m) weight(m) = weight(n) ENDIF ENDIF ENDDO ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass ENDDO ENDIF IF ( ANY(weight < 0.0_wp) ) THEN WRITE( message_string, * ) 'negative weighting' CALL message( 'lpm_droplet_collision', 'PA0028', & 2, 2, -1, 6, 1 ) ENDIF particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) / & ( weight(1:number_of_particles) & * factor_volume_to_mass & ) & )**0.33333333333333_wp particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles) DEALLOCATE( weight, mass, ckernel ) ENDIF ! !-- Check if LWC is conserved during collision process IF ( ql_v(k,j,i) /= 0.0_wp ) THEN IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp .OR. & ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp ) THEN WRITE( message_string, * ) ' LWC is not conserved during', & ' collision! ', & ' LWC after condensation: ', ql_v(k,j,i), & ' LWC after collision: ', ql_vp(k,j,i) CALL message( 'lpm_droplet_collision', 'PA0040', 2, 2, -1, 6, 1 ) ENDIF ENDIF ENDIF CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' ) END SUBROUTINE lpm_droplet_collision