!> @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-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: lpm_droplet_collision.f90 1683 2015-10-07 23:57:51Z gronemeier $ ! ! 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. Three different kernels are available: !> PALM kernel: Kernel is approximated using a method from Rogers and !> Yau (1989, A Short Course in Cloud Physics, Pergamon Press). !> All droplets smaller than the treated one are represented by !> one droplet with mean features. Collision efficiencies are taken !> from the respective table in Rogers and Yau. !> 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, ql_v, ql_vp, u, v, w, zu, zw USE cloud_parameters, & ONLY: effective_coll_efficiency USE constants, & ONLY: pi USE control_parameters, & ONLY: dt_3d, message_string, u_gtrans, v_gtrans, dz USE cpulog, & ONLY: cpu_log, log_point_s USE grid_variables, & ONLY: ddx, dx, ddy, dy USE indices, & ONLY: nxl, nxr, nyn, nys, nzb, nzt USE kinds USE lpm_collision_kernels_mod, & ONLY: ckernel, collision_efficiency_rogers, recalculate_kernel USE particle_attributes, & ONLY: deleted_particles, dissipation_classes, hall_kernel, & palm_kernel, particles, particle_type, & prt_count, use_kernel_tables, wang_kernel USE pegrid IMPLICIT NONE INTEGER(iwp) :: eclass !< INTEGER(iwp) :: i !< INTEGER(iwp) :: ii !< INTEGER(iwp) :: inc !< INTEGER(iwp) :: is !< INTEGER(iwp) :: j !< INTEGER(iwp) :: jj !< INTEGER(iwp) :: js !< INTEGER(iwp) :: k !< INTEGER(iwp) :: kk !< INTEGER(iwp) :: n !< INTEGER(iwp) :: pse !< INTEGER(iwp) :: psi !< INTEGER(iwp) :: rclass_l !< INTEGER(iwp) :: rclass_s !< INTEGER(iwp), DIMENSION(prt_count(k,j,i)) :: rclass_v !< LOGICAL, SAVE :: first_flag = .TRUE. !< TYPE(particle_type) :: tmp_particle !< REAL(wp) :: aa !< REAL(wp) :: auxn !< temporary variables REAL(wp) :: auxs !< temporary variables REAL(wp) :: bb !< REAL(wp) :: cc !< REAL(wp) :: dd !< REAL(wp) :: ddV !< REAL(wp) :: delta_r !< REAL(wp) :: delta_v !< REAL(wp) :: epsilon !< REAL(wp) :: gg !< REAL(wp) :: mean_r !< REAL(wp) :: ql_int !< REAL(wp) :: ql_int_l !< REAL(wp) :: ql_int_u !< REAL(wp) :: r3 !< REAL(wp) :: sl_r3 !< REAL(wp) :: sl_r4 !< REAL(wp) :: sum1 !< REAL(wp) :: sum2 !< REAL(wp) :: sum3 !< REAL(wp) :: u_int !< REAL(wp) :: u_int_l !< REAL(wp) :: u_int_u !< REAL(wp) :: v_int !< REAL(wp) :: v_int_l !< REAL(wp) :: v_int_u !< REAL(wp) :: w_int !< REAL(wp) :: w_int_l !< REAL(wp) :: w_int_u !< REAL(wp) :: x !< REAL(wp) :: y !< REAL(wp), DIMENSION(:), ALLOCATABLE :: rad !< REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !< REAL, DIMENSION(prt_count(k,j,i)) :: ck REAL, DIMENSION(prt_count(k,j,i)) :: r3v REAL, DIMENSION(prt_count(k,j,i)) :: sum1v REAL, DIMENSION(prt_count(k,j,i)) :: sum2v CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' ) ! !-- Collision requires at least two particles in the box IF ( prt_count(k,j,i) > 1 ) THEN ! !-- First, sort particles within the gridbox by their size, !-- using Shell's method (see Numerical Recipes) !-- NOTE: In case of using particle tails, the re-sorting of !-- ---- tails would have to be included here! IF ( .NOT. ( ( hall_kernel .OR. wang_kernel ) .AND. & use_kernel_tables ) ) THEN psi = 0 inc = 1 DO WHILE ( inc <= prt_count(k,j,i) ) inc = 3 * inc + 1 ENDDO DO WHILE ( inc > 1 ) inc = inc / 3 DO is = inc+1, prt_count(k,j,i) tmp_particle = particles(psi+is) js = is DO WHILE ( particles(psi+js-inc)%radius > & tmp_particle%radius ) particles(psi+js) = particles(psi+js-inc) js = js - inc IF ( js <= inc ) EXIT ENDDO particles(psi+js) = tmp_particle ENDDO ENDDO ENDIF psi = 1 pse = prt_count(k,j,i) ! !-- Now apply the different kernels IF ( ( hall_kernel .OR. wang_kernel ) .AND. & use_kernel_tables ) THEN ! !-- Fast method with pre-calculated efficiencies 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 / 1000.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) !-- Since new radii after collision are defined by radii of all !-- droplets before collision, temporary fields for new radii and !-- weighting factors are needed ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) rad = 0.0_wp weight = 0.0_wp sum1v(1:prt_count(k,j,i)) = 0.0_wp sum2v(1:prt_count(k,j,i)) = 0.0_wp DO n = 1, prt_count(k,j,i) rclass_l = particles(n)%class ! !-- Mass added due to collisions with smaller droplets DO is = n+1, prt_count(k,j,i) rclass_s = particles(is)%class auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor IF ( particles(is)%radius < particles(n)%radius ) THEN sum1v(n) = sum1v(n) + particles(is)%radius**3 * auxs sum2v(is) = sum2v(is) + auxn ELSE sum2v(n) = sum2v(n) + auxs sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn ENDIF ENDDO ENDDO rclass_v = particles(1:prt_count(k,j,i))%class DO n = 1, prt_count(k,j,i) ck(n) = ckernel(rclass_v(n),rclass_v(n),eclass) ENDDO r3v = particles(1:prt_count(k,j,i))%radius**3 DO n = 1, prt_count(k,j,i) sum3 = 0.0_wp ddV = ddx * ddy / dz ! !-- Change of the current weighting factor sum3 = 1 - dt_3d * ddV * & ck(n) * & ( particles(n)%weight_factor - 1 ) * 0.5_wp - & dt_3d * ddV * sum2v(n) weight(n) = particles(n)%weight_factor * sum3 ! !-- Change of the current droplet radius rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/& sum3 )**0.33333333333333_wp ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) & * rad(n)**3 ENDDO 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(psi:pse)%radius = rad(1:prt_count(k,j,i)) particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) DEALLOCATE(rad, weight) ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & .NOT. use_kernel_tables ) THEN ! !-- Collision efficiencies 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:prt_count(k,j,i),1:prt_count(k,j,i),1:1) ) ! !-- Now calculate collision efficiencies for this box CALL recalculate_kernel( i, j, k ) ! !-- Droplet collision are calculated using collision-coalescence !-- formulation proposed by Wang (see PALM documentation) !-- Since new radii after collision are defined by radii of all !-- droplets before collision, temporary fields for new radii and !-- weighting factors are needed ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i))) rad = 0.0_wp weight = 0.0_wp DO n = psi, pse, 1 sum1 = 0.0_wp sum2 = 0.0_wp sum3 = 0.0_wp ! !-- Mass added due to collisions with smaller droplets DO is = psi, n-1 sum1 = sum1 + ( particles(is)%radius**3 * & ckernel(n,is,1) * & particles(is)%weight_factor ) ENDDO ! !-- Rate of collisions with larger droplets DO is = n+1, pse sum2 = sum2 + ( ckernel(n,is,1) * & particles(is)%weight_factor ) ENDDO r3 = particles(n)%radius**3 ddV = ddx * ddy / dz is = 1 ! !-- Change of the current weighting factor sum3 = 1 - dt_3d * ddV * & ckernel(n,n,1) * & ( particles(n)%weight_factor - 1 ) * 0.5_wp - & dt_3d * ddV * sum2 weight(n-is+1) = particles(n)%weight_factor * sum3 ! !-- Change of the current droplet radius rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/& sum3 )**0.33333333333333_wp IF ( weight(n-is+1) < 0.0_wp ) THEN WRITE( message_string, * ) 'negative weighting', & 'factor: ', weight(n-is+1) CALL message( 'lpm_droplet_collision', 'PA0037', & 2, 2, -1, 6, 1 ) ENDIF ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) & * rad(n-is+1)**3 ENDDO particles(psi:pse)%radius = rad(1:prt_count(k,j,i)) particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i)) DEALLOCATE( rad, weight, ckernel ) ELSEIF ( palm_kernel ) THEN ! !-- PALM collision kernel ! !-- Calculate the mean radius of all those particles which !-- are of smaller size than the current particle and !-- use this radius for calculating the collision efficiency DO n = psi+prt_count(k,j,i)-1, psi+1, -1 sl_r3 = 0.0_wp sl_r4 = 0.0_wp DO is = n-1, psi, -1 IF ( particles(is)%radius < particles(n)%radius ) & THEN sl_r3 = sl_r3 + particles(is)%weight_factor & * particles(is)%radius**3 sl_r4 = sl_r4 + particles(is)%weight_factor & * particles(is)%radius**4 ENDIF ENDDO IF ( ( sl_r3 ) > 0.0_wp ) THEN mean_r = ( sl_r4 ) / ( sl_r3 ) CALL collision_efficiency_rogers( mean_r, & particles(n)%radius, & effective_coll_efficiency ) ELSE effective_coll_efficiency = 0.0_wp ENDIF IF ( effective_coll_efficiency > 1.0_wp .OR. & effective_coll_efficiency < 0.0_wp ) & THEN WRITE( message_string, * ) 'collision_efficien' , & 'cy out of range:' ,effective_coll_efficiency CALL message( 'lpm_droplet_collision', 'PA0145', 2, & 2, -1, 6, 1 ) ENDIF ! !-- Interpolation of liquid water content ii = particles(n)%x * ddx jj = particles(n)%y * ddy kk = ( particles(n)%z + 0.5_wp * dz ) / dz x = particles(n)%x - ii * dx y = particles(n)%y - jj * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * & ql(kk,jj,ii+1) & + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * & ql(kk,jj+1,ii+1) & ) / ( 3.0_wp * gg ) ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * & ql(kk+1,jj,ii+1) & + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * & ql(kk+1,jj+1,ii+1) & ) / ( 3.0_wp * gg ) ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *& ( ql_int_u - ql_int_l ) ! !-- Interpolate u velocity-component ii = ( particles(n)%x + 0.5_wp * dx ) * ddx jj = particles(n)%y * ddy kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1 x = particles(n)%x + ( 0.5_wp - ii ) * dx y = particles(n)%y - jj * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * & u(kk,jj,ii+1) & + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * & u(kk,jj+1,ii+1) & ) / ( 3.0_wp * gg ) - u_gtrans IF ( kk+1 == nzt+1 ) THEN u_int = u_int_l ELSE u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * & u(kk+1,jj,ii+1) & + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * & u(kk+1,jj+1,ii+1) & ) / ( 3.0_wp * gg ) - u_gtrans u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz & * ( u_int_u - u_int_l ) ENDIF ! !-- Same procedure for interpolation of the v velocity-component !-- (adopt index k from u velocity-component) ii = particles(n)%x * ddx jj = ( particles(n)%y + 0.5_wp * dy ) * ddy x = particles(n)%x - ii * dx y = particles(n)%y + ( 0.5_wp - jj ) * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & v(kk,jj,ii+1) & + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & v(kk,jj+1,ii+1) & ) / ( 3.0_wp * gg ) - v_gtrans IF ( kk+1 == nzt+1 ) THEN v_int = v_int_l ELSE v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * & v(kk+1,jj,ii+1) & + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * & v(kk+1,jj+1,ii+1) & ) / ( 3.0_wp * gg ) - v_gtrans v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz & * ( v_int_u - v_int_l ) ENDIF ! !-- Same procedure for interpolation of the w velocity-component !-- (adopt index i from v velocity-component) jj = particles(n)%y * ddy kk = particles(n)%z / dz x = particles(n)%x - ii * dx y = particles(n)%y - jj * dy aa = x**2 + y**2 bb = ( dx - x )**2 + y**2 cc = x**2 + ( dy - y )**2 dd = ( dx - x )**2 + ( dy - y )**2 gg = aa + bb + cc + dd w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & w(kk,jj,ii+1) & + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & w(kk,jj+1,ii+1) & ) / ( 3.0_wp * gg ) IF ( kk+1 == nzt+1 ) THEN w_int = w_int_l ELSE w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * & w(kk+1,jj,ii+1) & + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * & w(kk+1,jj+1,ii+1) & ) / ( 3.0_wp * gg ) w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz & * ( w_int_u - w_int_l ) ENDIF ! !-- Change in radius due to collision delta_r = effective_coll_efficiency / 3.0_wp & * pi * sl_r3 * ddx * ddy / dz & * SQRT( ( u_int - particles(n)%speed_x )**2 & + ( v_int - particles(n)%speed_y )**2 & + ( w_int - particles(n)%speed_z )**2 & ) * dt_3d ! !-- Change in volume due to collision delta_v = particles(n)%weight_factor & * ( ( particles(n)%radius + delta_r )**3 & - particles(n)%radius**3 ) ! !-- Check if collected particles provide enough LWC for !-- volume change of collector particle IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0_wp ) THEN delta_r = ( ( sl_r3/particles(n)%weight_factor ) & + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp ) & - particles(n)%radius DO is = n-1, psi, -1 IF ( particles(is)%radius < particles(n)%radius ) THEN particles(is)%weight_factor = 0.0_wp particles(is)%particle_mask = .FALSE. deleted_particles = deleted_particles + 1 ENDIF ENDDO ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0_wp ) THEN DO is = n-1, psi, -1 IF ( particles(is)%radius < particles(n)%radius & .AND. sl_r3 > 0.0_wp ) THEN particles(is)%weight_factor = & ( ( particles(is)%weight_factor & * ( particles(is)%radius**3 ) ) & - ( delta_v & * particles(is)%weight_factor & * ( particles(is)%radius**3 ) & / sl_r3 ) ) & / ( particles(is)%radius**3 ) IF ( particles(is)%weight_factor < 0.0_wp ) THEN WRITE( message_string, * ) 'negative ', & 'weighting factor: ', & particles(is)%weight_factor CALL message( 'lpm_droplet_collision', & 'PA0039', & 2, 2, -1, 6, 1 ) ENDIF ENDIF ENDDO ENDIF particles(n)%radius = particles(n)%radius + delta_r ql_vp(k,j,i) = ql_vp(k,j,i) + & particles(n)%weight_factor * & ( particles(n)%radius**3 ) ENDDO ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & * particles(psi)%radius**3 ENDIF ! collision kernel ELSE IF ( prt_count(k,j,i) == 1 ) THEN psi = 1 ! !-- Calculate change of weighting factor due to self collision IF ( ( hall_kernel .OR. wang_kernel ) .AND. & use_kernel_tables ) THEN IF ( wang_kernel ) THEN eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.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 ddV = ddx * ddy / dz rclass_l = particles(psi)%class sum3 = 1 - dt_3d * ddV * & ( ckernel(rclass_l,rclass_l,eclass) * & ( particles(psi)%weight_factor-1 ) * 0.5_wp ) particles(psi)%radius = ( particles(psi)%radius**3 / & sum3 )**0.33333333333333_wp particles(psi)%weight_factor = particles(psi)%weight_factor & * sum3 ELSE IF ( ( hall_kernel .OR. wang_kernel ) .AND. & .NOT. use_kernel_tables ) THEN ! !-- Collision efficiencies 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(psi:psi, psi:psi, 1:1) ) ! !-- Now calculate collision efficiencies for this box CALL recalculate_kernel( i, j, k ) ddV = ddx * ddy / dz sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) * & ( particles(psi)%weight_factor - 1 ) * 0.5_wp ) particles(psi)%radius = ( particles(psi)%radius**3 / & sum3 )**0.33333333333333_wp particles(psi)%weight_factor = particles(psi)%weight_factor & * sum3 DEALLOCATE( ckernel ) ENDIF ql_vp(k,j,i) = particles(psi)%weight_factor * & particles(psi)%radius**3 ENDIF ! !-- Check if condensation of LWC was 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 CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' ) END SUBROUTINE lpm_droplet_collision