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-2012 Leibniz University Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: lpm_droplet_collision.f90 1037 2012-10-22 14:10:22Z schwenkel $ ! ! 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 chang 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, integral, lw_max, & 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 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 / 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 DO n = pse, psi+1, -1 integral = 0.0 lw_max = 0.0 rclass_l = particles(n)%class ! !-- Calculate growth of collector particle radius using all !-- droplets smaller than current droplet DO is = psi, n-1 rclass_s = particles(is)%class integral = integral + & ( particles(is)%radius**3 * & ckernel(rclass_l,rclass_s,eclass) * & particles(is)%weight_factor ) ! !-- Calculate volume of liquid water of the collected !-- droplets which is the maximum liquid water available !-- for droplet growth lw_max = lw_max + ( particles(is)%radius**3 * & particles(is)%weight_factor ) ENDDO ! !-- Change in radius of collector droplet due to collision delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & integral * dt_3d * ddx * ddy / dz ! !-- Change in volume of collector droplet due to collision delta_v = particles(n)%weight_factor & * ( ( particles(n)%radius + delta_r )**3 & - particles(n)%radius**3 ) IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN !-- replace by message call PRINT*, 'Particle has grown to much because the', & ' change of volume of particle is larger', & ' than liquid water available!' ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN !-- can this case really happen?? DO is = psi, n-1 particles(is)%weight_factor = 0.0 particle_mask(is) = .FALSE. deleted_particles = deleted_particles + 1 ENDDO ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN ! !-- Calculate new weighting factor of collected droplets DO is = psi, n-1 rclass_s = particles(is)%class particles(is)%weight_factor = & particles(is)%weight_factor - & ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor & / integral ) * delta_v ) IF ( particles(is)%weight_factor < 0.0 ) THEN WRITE( message_string, * ) 'negative ', & 'weighting factor: ', & particles(is)%weight_factor CALL message( 'lpm_droplet_collision', '', & 2, 2, -1, 6, 1 ) ELSEIF ( particles(is)%weight_factor == 0.0 ) & THEN particles(is)%weight_factor = 0.0 particle_mask(is) = .FALSE. deleted_particles = deleted_particles + 1 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 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 ) DO n = pse, psi+1, -1 integral = 0.0 lw_max = 0.0 ! !-- Calculate growth of collector particle radius using all !-- droplets smaller than current droplet DO is = psi, n-1 integral = integral + particles(is)%radius**3 * & ckernel(n,is,1) * & particles(is)%weight_factor ! !-- Calculate volume of liquid water of the collected !-- droplets which is the maximum liquid water available !-- for droplet growth lw_max = lw_max + ( particles(is)%radius**3 * & particles(is)%weight_factor ) ENDDO ! !-- Change in radius of collector droplet due to collision delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & integral * dt_3d * ddx * ddy / dz ! !-- Change in volume of collector droplet due to collision delta_v = particles(n)%weight_factor & * ( ( particles(n)%radius + delta_r )**3 & - particles(n)%radius**3 ) IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN !-- replace by message call PRINT*, 'Particle has grown to much because the', & ' change of volume of particle is larger', & ' than liquid water available!' ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN !-- can this case really happen?? DO is = psi, n-1 particles(is)%weight_factor = 0.0 particle_mask(is) = .FALSE. deleted_particles = deleted_particles + 1 ENDDO ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN ! !-- Calculate new weighting factor of collected droplets DO is = psi, n-1 particles(is)%weight_factor = & particles(is)%weight_factor - & ( ckernel(n,is,1) / integral * delta_v * & particles(is)%weight_factor ) IF ( particles(is)%weight_factor < 0.0 ) THEN WRITE( message_string, * ) 'negative ', & 'weighting factor: ', & particles(is)%weight_factor CALL message( 'lpm_droplet_collision', '', & 2, 2, -1, 6, 1 ) ELSEIF ( particles(is)%weight_factor == 0.0 ) & THEN particles(is)%weight_factor = 0.0 particle_mask(is) = .FALSE. deleted_particles = deleted_particles + 1 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 DEALLOCATE( 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', '', & 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 ENDIF ! collision kernel ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & * particles(psi)%radius**3 ELSE IF ( prt_count(k,j,i) == 1 ) THEN psi = prt_start_index(k,j,i) 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', '', 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