SUBROUTINE lpm_droplet_collision !--------------------------------------------------------------------------------! ! 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 1323 2014-03-20 17:09:54Z maronga $ ! ! 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. !------------------------------------------------------------------------------! 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_mask, particle_type, & prt_count, prt_start_index, use_kernel_tables, wang_kernel 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 !: REAL(wp) :: aa !: 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 !: TYPE(particle_type) :: tmp_particle !: CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- 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! psi = prt_start_index(k,j,i) - 1 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 psi = prt_start_index(k,j,i) pse = psi + prt_count(k,j,i)-1 ! !-- 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 ENDIF IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 ) 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 weight = 0.0 DO n = psi, pse, 1 sum1 = 0.0 sum2 = 0.0 sum3 = 0.0 rclass_l = particles(n)%class ! !-- Mass added due to collisions with smaller droplets DO is = psi, n-1 rclass_s = particles(is)%class sum1 = sum1 + ( particles(is)%radius**3 * & ckernel(rclass_l,rclass_s,eclass) * & particles(is)%weight_factor ) ENDDO ! !-- Rate of collisions with larger droplets DO is = n+1, pse rclass_s = particles(is)%class sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) * & particles(is)%weight_factor ) ENDDO r3 = particles(n)%radius**3 ddV = ddx * ddy / dz is = prt_start_index(k,j,i) ! !-- Change of the current weighting factor sum3 = 1 - dt_3d * ddV * & ckernel(rclass_l,rclass_l,eclass) * & ( particles(n)%weight_factor - 1 ) * 0.5 - & 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 ) THEN WRITE( message_string, * ) 'negative weighting', & 'factor: ', weight(n-is+1) CALL message( 'lpm_droplet_collision', 'PA0028', & 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) 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(prt_start_index(k,j,i): & prt_start_index(k,j,i)+prt_count(k,j,i)-1, & prt_start_index(k,j,i): & prt_start_index(k,j,i)+prt_count(k,j,i)-1,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 weight = 0.0 DO n = psi, pse, 1 sum1 = 0.0 sum2 = 0.0 sum3 = 0.0 ! !-- 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 = prt_start_index(k,j,i) ! !-- Change of the current weighting factor sum3 = 1 - dt_3d * ddV * & ckernel(n,n,1) * & ( particles(n)%weight_factor - 1 ) * 0.5 - & 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 ) 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 sl_r4 = 0.0 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 ) 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 ENDIF IF ( effective_coll_efficiency > 1.0 .OR. & effective_coll_efficiency < 0.0 ) & 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 ... ii = particles(n)%x * ddx jj = particles(n)%y * ddy kk = ( particles(n)%z + 0.5 * 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 * 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 * 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 * dx ) * ddx jj = particles(n)%y * ddy kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1 x = particles(n)%x + ( 0.5 - 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 * 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 * 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-com- !-- ponent (adopt index k from u velocity-component) ii = particles(n)%x * ddx jj = ( particles(n)%y + 0.5 * dy ) * ddy x = particles(n)%x - ii * dx y = particles(n)%y + ( 0.5 - 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 * 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 * 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-com- !-- ponent (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 * 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 * 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 ) 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 particle_mask(is) = .FALSE. deleted_particles = deleted_particles + 1 ENDIF ENDDO ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN DO is = n-1, psi, -1 IF ( particles(is)%radius < particles(n)%radius & .AND. sl_r3 > 0.0 ) 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 ) 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 = prt_start_index(k,j,i) ! !-- 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 ENDIF IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 ) 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 ) 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 ) 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 ) THEN IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. & ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) 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 ENDDO ENDDO ENDDO CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' ) END SUBROUTINE lpm_droplet_collision