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 1310 2014-03-14 08:01:56Z heinze $ ! ! 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 USE cloud_parameters USE constants USE control_parameters USE cpulog USE grid_variables USE indices USE interfaces USE lpm_collision_kernels_mod USE particle_attributes IMPLICIT NONE INTEGER :: eclass, i, ii, inc, is, j, jj, js, k, kk, n, pse, psi, rclass_l, & rclass_s REAL :: aa, bb, cc, dd, delta_r, delta_v, gg, epsilon, mean_r, ql_int, & ql_int_l, ql_int_u, u_int, u_int_l, u_int_u, v_int, v_int_l, & v_int_u, w_int, w_int_l, w_int_u, sl_r3, sl_r4, x, y, sum1, sum2, & sum3, r3, ddV TYPE(particle_type) :: tmp_particle REAL, DIMENSION(:), ALLOCATABLE :: rad, weight 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 / 1000.0 * & dissipation_classes ) + 1 epsilon = diss(k,j,i) ELSE epsilon = 0.0 ENDIF IF ( hall_kernel .OR. epsilon * 1.0E4 < 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 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 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 & * 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./3. ) & - 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 / 1000.0 * & dissipation_classes ) + 1 epsilon = diss(k,j,i) ELSE epsilon = 0.0 ENDIF IF ( hall_kernel .OR. epsilon * 1.0E4 < 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 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 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