Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 824)
+++ palm/trunk/SOURCE/Makefile	(revision 825)
@@ -4,5 +4,5 @@
 # Current revisions:
 # ------------------
-# 
+# wang_kernel renamed lpm_collision_kernels
 #
 # Former revisions:
@@ -83,7 +83,7 @@
 	init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
 	init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
-	interaction_droplets_ptq.f90 local_flush.f90 \
-        local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
-	local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \
+	interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
+	local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90 \
+        lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \
 	package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
         plant_canopy_model.f90 poisfft.f90 \
@@ -108,5 +108,5 @@
 	user_particle_attributes.f90 user_read_restart_data.f90 \
 	user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
-	wang_kernel.f90 write_3d_binary.f90 write_compressed.f90 \
+	write_3d_binary.f90 write_compressed.f90 \
 	write_var_list.f90
 
@@ -124,11 +124,13 @@
 	data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \
 	diffusion_v.o diffusion_w.o diffusivities.o disturb_field.o \
-	disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o exchange_horiz_2d.o fft_xy.o \
-	flow_statistics.o global_min_max.o header.o impact_of_latent_heat.o \
-	inflow_turbulence.o init_1d_model.o init_3d_model.o init_advec.o init_cloud_physics.o \
-	init_coupling.o init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o init_pegrid.o \
-	init_pt_anomaly.o init_rankine.o init_slope.o \
+	disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o \
+	exchange_horiz_2d.o fft_xy.o flow_statistics.o global_min_max.o \
+	header.o impact_of_latent_heat.o inflow_turbulence.o init_1d_model.o \
+	init_3d_model.o init_advec.o init_cloud_physics.o init_coupling.o \
+	init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o \
+	init_pegrid.o init_pt_anomaly.o init_rankine.o init_slope.o \
 	interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
-	local_system.o local_tremain.o local_tremain_ini.o message.o modules.o \
+	local_system.o local_tremain.o local_tremain_ini.o \
+	lpm_collision_kernels.f90 message.o modules.o \
 	netcdf.o package_parin.o palm.o parin.o particle_boundary_conds.o \
         plant_canopy_model.o poisfft.o \
@@ -151,5 +153,5 @@
 	user_module.o user_parin.o user_particle_attributes.o \
 	user_read_restart_data.o user_spectra.o user_statistics.o \
-	wall_fluxes.o wang_kernel.o write_3d_binary.o write_compressed.o \
+	wall_fluxes.o write_3d_binary.o write_compressed.o \
 	write_var_list.o
 
@@ -181,5 +183,5 @@
 
 
-advec_particles.o: modules.o random_function.o wang_kernel.o
+advec_particles.o: modules.o lpm_collision_kernels.o random_function.o
 advec_s_bc.o: modules.o
 advec_s_pw.o: modules.o
@@ -257,4 +259,5 @@
 local_tremain.o: modules.o
 local_tremain_ini.o: modules.o
+lpm_collision_kernels.o: modules.o user_module.o
 message.o: modules.o
 modules.o: modules.f90
@@ -328,5 +331,4 @@
 user_statistics.o: modules.o user_module.o
 wall_fluxes.o: modules.o
-wang_kernel.o: modules.o user_module.o
 write_3d_binary.o: modules.o random_function.o
 write_compressed.o: modules.o
Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 824)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 825)
@@ -4,7 +4,9 @@
 ! Current revisions:
 ! ------------------
-! droplet growth by condensation may include curvature and solution effects.
-! initialisation of temporary particle array for resorting removed.
-! particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
+! droplet growth by condensation may include curvature and solution effects,
+! initialisation of temporary particle array for resorting removed,
+! particle attributes speed_x|y|z_sgs renamed rvar1|2|3,
+! module wang_kernel_mod renamed lpm_collision_kernels_mod,
+! wang_collision_kernel renamed wang_kernel
 !
 ! Former revisions:
@@ -130,4 +132,5 @@
     USE indices
     USE interfaces
+    USE lpm_collision_kernels_mod
     USE netcdf_control
     USE particle_attributes
@@ -135,5 +138,4 @@
     USE random_function_mod
     USE statistics
-    USE wang_kernel_mod
 
     IMPLICIT NONE
@@ -200,9 +202,5 @@
     CALL cpu_log( log_point(25), 'advec_particles', 'start' )
 
-!    IF ( number_of_particles /= number_of_tails )  THEN
-!       WRITE (9,*) '--- advec_particles: #1'
-!       WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
-!       CALL local_flush( 9 )
-!    ENDIF
+
 !
 !-- Write particle data on file for later analysis.
@@ -270,16 +268,4 @@
     CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
 
-!    WRITE ( 9, * ) '*** advec_particles: ##0.3'
-!    CALL local_flush( 9 )
-!    nd = 0
-!    DO  n = 1, number_of_particles
-!       IF ( .NOT. particle_mask(n) )  nd = nd + 1
-!    ENDDO
-!    IF ( nd /= deleted_particles ) THEN
-!       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
-!       CALL local_flush( 9 )
-!       CALL MPI_ABORT( comm2d, 9999, ierr )
-!    ENDIF
-
 !
 !-- Particle (droplet) growth by condensation/evaporation and collision
@@ -357,157 +343,6 @@
 !
 !--       Change in radius by condensation/evaporation
-          IF ( curvature_solution_effects .AND. particles(n)%radius < 1.0E-6 ) &
-          THEN
-!
-!--          Curvature and solutions effects are included in growth equation.
-!--          Change in Radius is calculated with 4th-order Rosenbrock method
-!--          for stiff o.d.e's with monitoring local truncation error to adjust
-!--          stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
-!--          For larger radii the simple analytic method (see ELSE) gives
-!--          almost the same results.
-!
-!--          Surface tension after (Straka, 2009)
-             sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
-
-             r_ros = particles(n)%radius
-             dt_ros_sum  = 0.0      ! internal integrated time (s)
-             internal_timestep_count = 0
-
-             ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
-                               ( l_v / ( r_v * t_int ) - 1.0 ) *               &
-                               rho_l * l_v / ( thermal_conductivity_l * t_int )&
-                             )
-              
-             afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
-
-             IF ( particles(n)%rvar3 == -9999999.9 )  THEN
-!
-!--             First particle timestep. Derivative has to be calculated.
-                drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
-                                             afactor / r_ros +    &
-                                             bfactor / r_ros**3 )
-             ELSE
-!
-!--             Take value from last PALM timestep
-                drdt_m = particles(n)%rvar3
-             ENDIF
-!
-!--          Take internal timestep values from the end of last PALM timestep
-             dt_ros_last = particles(n)%rvar1
-             dt_ros_next = particles(n)%rvar2
-!
-!--          Internal timestep must not be larger than PALM timestep
-             dt_ros = MIN( dt_ros_next, dt_3d )
-!
-!--          Integrate growth equation in time unless PALM timestep is reached
-             DO WHILE ( dt_ros_sum < dt_3d )
-
-                internal_timestep_count = internal_timestep_count + 1
-
-!
-!--             Derivative at starting value
-                drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
-                                          bfactor / r_ros**3 )
-                drdt_ini       = drdt
-                dt_ros_sum_ini = dt_ros_sum
-                r_ros_ini      = r_ros
-
-!
-!--             Calculate time derivative of dr/dt
-                d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 
-
-!
-!--             Calculate radial derivative of dr/dt
-                d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
-                                     2.0 * afactor / r_ros**3 -     &
-                                     4.0 * bfactor / r_ros**5 )
-!
-!--             Adjust stepsize unless required accuracy is reached
-                DO  jtry = 1, maxtry+1
-
-                   IF ( jtry == maxtry+1 )  THEN
-                      message_string = 'maxtry > 40 in Rosenbrock method'
-                      CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
-                   ENDIF
-
-                   aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
-                   g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
-                   r_ros = r_ros_ini + a21 * g1
-                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
-                                              afactor / r_ros + &
-                                              bfactor / r_ros**3 )
-
-                   g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
-                           / aa
-                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
-                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
-                                              afactor / r_ros + &
-                                              bfactor / r_ros**3 )
-
-                   g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
-                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
-                   g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
-                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
-                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
-
-                   dt_ros_sum = dt_ros_sum_ini + dt_ros
-
-                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
-                      message_string = 'zero stepsize in Rosenbrock method'
-                      CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
-                   ENDIF
-!
-!--                Calculate error
-                   err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 
-                   errmax  = 0.0
-                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
-!
-!--                Leave loop if accuracy is sufficient, otherwise try again
-!--                with a reduced stepsize
-                   IF ( errmax <= 1.0 )  THEN
-                      EXIT
-                   ELSE
-                      dt_ros = SIGN(                                           &
-                                    MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
-                                         shrnk * ABS( dt_ros ) ),              &
-                                    dt_ros                                     &
-                                   )
-                   ENDIF
-
-                ENDDO  ! loop for stepsize adjustment
-
-!
-!--             Calculate next internal timestep
-                IF ( errmax > errcon )  THEN
-                   dt_ros_next = 0.9 * dt_ros * errmax**pgrow
-                ELSE
-                   dt_ros_next = grow * dt_ros
-                ENDIF
-
-!
-!--             Estimated timestep is reduced if the PALM time step is exceeded
-                dt_ros_last = dt_ros
-                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
-                   dt_ros = dt_3d - dt_ros_sum
-                ELSE
-                   dt_ros = dt_ros_next
-                ENDIF
-
-                drdt_m = drdt
-
-             ENDDO
-!
-!--          Store derivative and internal timestep values for next PALM step
-             particles(n)%rvar1 = dt_ros_last
-             particles(n)%rvar2 = dt_ros_next
-             particles(n)%rvar3 = drdt_m
-
-             new_r = r_ros
-!
-!--          Radius should not fall below 1E-8 because Rosenbrock method may
-!--          lead to errors otherwise
-             new_r = MAX( new_r, 1.0E-8 )
-
-          ELSE
+          IF ( particles(n)%radius >= 1.0E-6  .OR.  &
+               .NOT. curvature_solution_effects )  THEN
 !
 !--          Approximation for large radii, where curvature and solution
@@ -523,4 +358,158 @@
                 new_r = SQRT( arg )
              ENDIF
+          ENDIF
+
+          IF ( curvature_solution_effects  .AND.                              &
+               ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
+          THEN
+!
+!--          Curvature and solutions effects are included in growth equation.
+!--          Change in Radius is calculated with 4th-order Rosenbrock method
+!--          for stiff o.d.e's with monitoring local truncation error to adjust
+!--          stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
+!--          For larger radii the simple analytic method (see ELSE) gives
+!--          almost the same results.
+!
+!--          Surface tension after (Straka, 2009)
+             sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
+
+             r_ros = particles(n)%radius
+             dt_ros_sum  = 0.0      ! internal integrated time (s)
+             internal_timestep_count = 0
+
+             ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
+                               ( l_v / ( r_v * t_int ) - 1.0 ) *               &
+                               rho_l * l_v / ( thermal_conductivity_l * t_int )&
+                             )
+              
+             afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
+
+             IF ( particles(n)%rvar3 == -9999999.9 )  THEN
+!
+!--             First particle timestep. Derivative has to be calculated.
+                drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
+                                             afactor / r_ros +    &
+                                             bfactor / r_ros**3 )
+             ELSE
+!
+!--             Take value from last PALM timestep
+                drdt_m = particles(n)%rvar3
+             ENDIF
+!
+!--          Take internal timestep values from the end of last PALM timestep
+             dt_ros_last = particles(n)%rvar1
+             dt_ros_next = particles(n)%rvar2
+!
+!--          Internal timestep must not be larger than PALM timestep
+             dt_ros = MIN( dt_ros_next, dt_3d )
+!
+!--          Integrate growth equation in time unless PALM timestep is reached
+             DO WHILE ( dt_ros_sum < dt_3d )
+
+                internal_timestep_count = internal_timestep_count + 1
+
+!
+!--             Derivative at starting value
+                drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
+                                          bfactor / r_ros**3 )
+                drdt_ini       = drdt
+                dt_ros_sum_ini = dt_ros_sum
+                r_ros_ini      = r_ros
+
+!
+!--             Calculate time derivative of dr/dt
+                d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 
+
+!
+!--             Calculate radial derivative of dr/dt
+                d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
+                                     2.0 * afactor / r_ros**3 -     &
+                                     4.0 * bfactor / r_ros**5 )
+!
+!--             Adjust stepsize unless required accuracy is reached
+                DO  jtry = 1, maxtry+1
+
+                   IF ( jtry == maxtry+1 )  THEN
+                      message_string = 'maxtry > 40 in Rosenbrock method'
+                      CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
+                   ENDIF
+
+                   aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
+                   g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
+                   r_ros = r_ros_ini + a21 * g1
+                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
+                                              afactor / r_ros + &
+                                              bfactor / r_ros**3 )
+
+                   g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
+                           / aa
+                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
+                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
+                                              afactor / r_ros + &
+                                              bfactor / r_ros**3 )
+
+                   g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
+                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
+                   g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
+                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
+                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
+
+                   dt_ros_sum = dt_ros_sum_ini + dt_ros
+
+                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
+                      message_string = 'zero stepsize in Rosenbrock method'
+                      CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
+                   ENDIF
+!
+!--                Calculate error
+                   err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 
+                   errmax  = 0.0
+                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
+!
+!--                Leave loop if accuracy is sufficient, otherwise try again
+!--                with a reduced stepsize
+                   IF ( errmax <= 1.0 )  THEN
+                      EXIT
+                   ELSE
+                      dt_ros = SIGN(                                           &
+                                    MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
+                                         shrnk * ABS( dt_ros ) ),              &
+                                    dt_ros                                     &
+                                   )
+                   ENDIF
+
+                ENDDO  ! loop for stepsize adjustment
+
+!
+!--             Calculate next internal timestep
+                IF ( errmax > errcon )  THEN
+                   dt_ros_next = 0.9 * dt_ros * errmax**pgrow
+                ELSE
+                   dt_ros_next = grow * dt_ros
+                ENDIF
+
+!
+!--             Estimated timestep is reduced if the PALM time step is exceeded
+                dt_ros_last = dt_ros
+                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
+                   dt_ros = dt_3d - dt_ros_sum
+                ELSE
+                   dt_ros = dt_ros_next
+                ENDIF
+
+                drdt_m = drdt
+
+             ENDDO
+!
+!--          Store derivative and internal timestep values for next PALM step
+             particles(n)%rvar1 = dt_ros_last
+             particles(n)%rvar2 = dt_ros_next
+             particles(n)%rvar3 = drdt_m
+
+             new_r = r_ros
+!
+!--          Radius should not fall below 1E-8 because Rosenbrock method may
+!--          lead to errors otherwise
+             new_r = MAX( new_r, 1.0E-8 )
 
           ENDIF
@@ -570,7 +559,7 @@
 !
 !--    Particle growth by collision
+       IF ( collision_kernel /= 'none' )  THEN
+
        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
-
-       IF ( wang_collision_kernel ) THEN
 
        DO  i = nxl, nxr
@@ -609,101 +598,338 @@
                    pse = psi + prt_count(k,j,i)-1
 
-                   ALLOCATE( kern(psi:pse,psi:pse) )
-
-!
-!--                Calculate collision kernel for all particles in grid box
-                   CALL colker( i, j, k, kern )
-!
-!--                collison kernel is calculated in cm**3/s but needed in m**3/s
-                   kern = kern * 1.0E-6
-
-                   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 * kern(n,is) * &
-                                      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 )
+                   IF ( hall_kernel  .OR.  wang_kernel )  THEN
+
+                      ALLOCATE( kern(psi:pse,psi:pse) )
+
+!
+!--                   Calculate collision kernel for all particles in grid box
+                      CALL colker( i, j, k, kern )
+!
+!--                   Collison kernel is calculated in cm**3/s but needed
+!--                   in m**3/s
+                      kern = kern * 1.0E-6   ! to be moved to colker
+
+                      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 * kern(n,is) *&
+                                         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
+
+                            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 - &
+                                  ( ( kern(n,is) * 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( 'advec_particles', '', 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
 
-!
-!--                   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
-
-                         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
-
-                         DO  is = psi, n-1
-
-                            particles(is)%weight_factor = 0.0
-                            particle_mask(is)  = .FALSE.
-                            deleted_particles = deleted_particles + 1
-
+                      DEALLOCATE( kern )
+
+                   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
 
-                      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 - &
-                               ( ( kern(n,is) * particles(is)%weight_factor / &
-                                   integral ) * delta_v )
-
-                            IF ( particles(is)%weight_factor < 0.0 )  THEN
-                                  WRITE( message_string, * ) 'negative ',      &
+                         IF ( ( sl_r3 ) > 0.0 )  THEN
+                            mean_r = ( sl_r4 ) / ( sl_r3 )
+
+                            CALL collision_efficiency( 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( 'advec_particles', '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( 'advec_particles', '', 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
+                                     CALL message( 'advec_particles', '', 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
 
-                   DEALLOCATE( kern )
 
                 ELSE IF ( prt_count(k,j,i) == 1 )  THEN
@@ -735,296 +961,9 @@
        ENDDO
 
-       ELSE
-
-       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
-
-!
-!--                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
-                   psi = prt_start_index(k,j,i)
-
-                   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( 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_efficiency ' , &
-                                   'out of range:' ,effective_coll_efficiency
-                         CALL message( 'advec_particles', '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 eq.dist
-
-                      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-compo-
-!--                   nent (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-compo-
-!--                   nent (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( 'advec_particles', '', 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 )
-
-                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.00001 .or.                &
-                       ql_vp(k,j,i)/ql_v(k,j,i) <= 0.99999  ) 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( 'advec_particles', '', 2, 2, -1, 6, 1 )
-                   ENDIF
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-       ENDIF
+       ENDIF   ! collision handling
 
        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
 
-    ENDIF
+    ENDIF   ! cloud droplet handling
 
 
@@ -2272,9 +2211,9 @@
              particles(n)%speed_x = particles(n)%speed_x * exp_term + &
                                     u_int * ( 1.0 - exp_term )
-              particles(n)%speed_y = particles(n)%speed_y * exp_term + &
+             particles(n)%speed_y = particles(n)%speed_y * exp_term + &
                                     v_int * ( 1.0 - exp_term )
-              particles(n)%speed_z = particles(n)%speed_z * exp_term +       &
-                              ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
-                                    * ( 1.0 - exp_term )
+             particles(n)%speed_z = particles(n)%speed_z * exp_term +       &
+                             ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
+                           * ( 1.0 - exp_term )
           ENDIF
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 824)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 825)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! check for curvature_solution_effects
+! check for collision_kernel and curvature_solution_effects
 !
 ! Former revisions:
@@ -709,4 +709,28 @@
     ENDIF
 
+!
+!-- Collision kernels:
+    SELECT CASE ( TRIM( collision_kernel ) )
+
+       CASE ( 'hall' )
+          hall_kernel = .TRUE.
+
+       CASE ( 'palm' )
+          palm_kernel = .TRUE.
+
+       CASE ( 'wang' )
+          wang_kernel = .TRUE.
+
+       CASE ( 'none' )
+
+
+       CASE DEFAULT
+          message_string = 'unknown collision kernel: collision_kernel = "' // &
+                           TRIM( collision_kernel ) // '"'
+          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
+
+    END SELECT
+
+
     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
Index: palm/trunk/SOURCE/data_output_ptseries.f90
===================================================================
--- palm/trunk/SOURCE/data_output_ptseries.f90	(revision 824)
+++ palm/trunk/SOURCE/data_output_ptseries.f90	(revision 825)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! mean particle radius added as output quantity.
+! mean/minimum/maximum particle radius added as output quantity.
 ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
 !
@@ -35,4 +35,5 @@
 !------------------------------------------------------------------------------!
 
+    USE cloud_parameters
     USE control_parameters
     USE cpulog
@@ -48,6 +49,5 @@
     INTEGER ::  i, inum, j, n
 
-    REAL, DIMENSION(0:number_of_particle_groups,30) ::  pts_value, pts_value_l
-
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  pts_value, pts_value_l
 
 
@@ -70,5 +70,9 @@
     ENDIF
 
+    ALLOCATE( pts_value(0:number_of_particle_groups,dopts_num), &
+              pts_value_l(0:number_of_particle_groups,dopts_num) )
+
     pts_value_l = 0.0
+    pts_value_l(:,16) = 9999999.9    ! for calculation of minimum radius
 
 !
@@ -88,7 +92,9 @@
        pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y     ! mean v
        pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z     ! mean w
-       pts_value_l(0,9)  = pts_value_l(0,9) + particles(n)%rvar1     ! mean sgsu
-       pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2    ! mean sgsv
-       pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3    ! mean sgsw
+       IF ( .NOT. curvature_solution_effects )  THEN
+          pts_value_l(0,9)  = pts_value_l(0,9) + particles(n)%rvar1  ! mean sgsu
+          pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
+          pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
+       ENDIF
        IF ( particles(n)%speed_z > 0.0 )  THEN
           pts_value_l(0,12) = pts_value_l(0,12) + 1.0  ! # of upward moving prts
@@ -99,7 +105,9 @@
                                             particles(n)%speed_z ! mean w down
        ENDIF
-       pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius
-       pts_value_l(0,16) = number_of_particles
-       pts_value_l(0,17) = number_of_particles
+       pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
+       pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
+       pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
+       pts_value_l(0,18) = number_of_particles
+       pts_value_l(0,19) = number_of_particles
 
 !
@@ -119,7 +127,9 @@
           pts_value_l(j,7)  = pts_value_l(j,7)  + particles(n)%speed_y
           pts_value_l(j,8)  = pts_value_l(j,8)  + particles(n)%speed_z
-          pts_value_l(j,9)  = pts_value_l(j,9)  + particles(n)%rvar1
-          pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2
-          pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3
+          IF ( .NOT. curvature_solution_effects )  THEN
+             pts_value_l(j,9)  = pts_value_l(j,9)  + particles(n)%rvar1
+             pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2
+             pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3
+          ENDIF
           IF ( particles(n)%speed_z > 0.0 )  THEN
              pts_value_l(j,12) = pts_value_l(j,12) + 1.0
@@ -129,6 +139,8 @@
           ENDIF
           pts_value_l(j,15) = pts_value_l(j,15) + particles(n)%radius
-          pts_value_l(j,16) = pts_value_l(j,16) + 1.0
-          pts_value_l(j,17) = pts_value_l(j,17) + 1.0
+          pts_value_l(j,16) = MIN( pts_value(j,16), particles(n)%radius )
+          pts_value_l(j,17) = MAX( pts_value(j,17), particles(n)%radius )
+          pts_value_l(j,18) = pts_value_l(j,18) + 1.0
+          pts_value_l(j,19) = pts_value_l(j,19) + 1.0
 
        ENDIF
@@ -146,15 +158,21 @@
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     CALL MPI_ALLREDUCE( pts_value_l(0,16), pts_value(0,16), inum, MPI_REAL, &
+                        MPI_MIN, comm2d, ierr )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
                         MPI_MAX, comm2d, ierr )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
+    CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum, MPI_REAL, &
+                        MPI_MAX, comm2d, ierr )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( pts_value_l(0,19), pts_value(0,19), inum, MPI_REAL, &
                         MPI_MIN, comm2d, ierr )
 #else
-    pts_value(:,1:17) = pts_value_l(:,1:17)
+    pts_value(:,1:19) = pts_value_l(:,1:19)
 #endif
 
 !
-!-- Normalize the above calculated quantities with the total number of
-!-- particles
+!-- Normalize the above calculated quantities (except min/max values) with the
+!-- total number of particles
     IF ( number_of_particle_groups > 1 )  THEN
        inum = number_of_particle_groups
@@ -186,22 +204,24 @@
     DO  n = 1, number_of_particles
 
-       pts_value_l(0,18) = pts_value_l(0,18) + ( particles(n)%x - &
+       pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
                            particles(n)%origin_x - pts_value(0,2) )**2 ! x*2
-       pts_value_l(0,19) = pts_value_l(0,19) + ( particles(n)%y - &
+       pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
                            particles(n)%origin_y - pts_value(0,3) )**2 ! y*2
-       pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%z - &
+       pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
                            particles(n)%origin_z - pts_value(0,4) )**2 ! z*2
-       pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%speed_x - &
-                                               pts_value(0,6) )**2     ! u*2
-       pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%speed_y - &
-                                               pts_value(0,7) )**2     ! v*2
-       pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_z - &
-                                               pts_value(0,8) )**2     ! w*2
-       pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%rvar1 - &
-                                               pts_value(0,9) )**2     ! u"2
-       pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%rvar2 - &
-                                               pts_value(0,10) )**2    ! v"2
-       pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar3 - &
-                                               pts_value(0,11) )**2    ! w"2
+       pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
+                                                pts_value(0,6) )**2   ! u*2
+       pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
+                                                 pts_value(0,7) )**2   ! v*2
+       pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
+                                                 pts_value(0,8) )**2   ! w*2
+       IF ( .NOT. curvature_solution_effects )  THEN
+          pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
+                                                    pts_value(0,9) )**2   ! u"2
+          pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
+                                                    pts_value(0,10) )**2  ! v"2
+          pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
+                                                    pts_value(0,11) )**2  ! w"2
+       ENDIF
 !
 !--    Repeat the same for the respective particle group
@@ -209,31 +229,33 @@
           j = particles(n)%group
 
-          pts_value_l(j,18) = pts_value_l(j,18) + ( particles(n)%x - &
+          pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%x - &
                               particles(n)%origin_x - pts_value(j,2) )**2
-          pts_value_l(j,19) = pts_value_l(j,19) + ( particles(n)%y - &
+          pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%y - &
                               particles(n)%origin_y - pts_value(j,3) )**2
-          pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%z - &
+          pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%z - &
                               particles(n)%origin_z - pts_value(j,4) )**2
-          pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%speed_x - &
-                                                  pts_value(j,6) )**2
-          pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%speed_y - &
-                                                  pts_value(j,7) )**2
-          pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_z - &
-                                                  pts_value(j,8) )**2
-          pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%rvar1 -   &
-                                                  pts_value(j,9) )**2
-          pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%rvar2 -   &
-                                                  pts_value(j,10) )**2
-          pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar3 -   &
-                                                  pts_value(j,11) )**2
+          pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_x - &
+                                                    pts_value(j,6) )**2
+          pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%speed_y - &
+                                                    pts_value(j,7) )**2
+          pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%speed_z - &
+                                                    pts_value(j,8) )**2
+          IF ( .NOT. curvature_solution_effects )  THEN
+             pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar1 - &
+                                                       pts_value(j,9) )**2
+             pts_value_l(j,27) = pts_value_l(j,27) + ( particles(n)%rvar2 - &
+                                                       pts_value(j,10) )**2
+             pts_value_l(j,28) = pts_value_l(j,28) + ( particles(n)%rvar3 - &
+                                                       pts_value(j,11) )**2
+          ENDIF
        ENDIF
 
     ENDDO
 
-    pts_value_l(0,27) = ( number_of_particles - pts_value(0,1) / numprocs )**2
+    pts_value_l(0,29) = ( number_of_particles - pts_value(0,1) / numprocs )**2
                                                  ! variance of particle numbers
     IF ( number_of_particle_groups > 1 )  THEN
        DO  j = 1, number_of_particle_groups
-          pts_value_l(j,27) = ( pts_value_l(j,1) - &
+          pts_value_l(j,29) = ( pts_value_l(j,1) - &
                                 pts_value(j,1) / numprocs )**2
        ENDDO
@@ -246,8 +268,8 @@
 
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum*10, MPI_REAL, &
+    CALL MPI_ALLREDUCE( pts_value_l(0,20), pts_value(0,20), inum*10, MPI_REAL, &
                         MPI_SUM, comm2d, ierr )
 #else
-    pts_value(:,18:27) = pts_value_l(:,18:27)
+    pts_value(:,20:29) = pts_value_l(:,20:29)
 #endif
 
@@ -264,7 +286,7 @@
 
        IF ( pts_value(j,1) > 0.0 )  THEN
-          pts_value(j,18:26) = pts_value(j,18:26) / pts_value(j,1)
-       ENDIF
-       pts_value(j,27) = pts_value(j,27) / numprocs
+          pts_value(j,20:28) = pts_value(j,20:28) / pts_value(j,1)
+       ENDIF
+       pts_value(j,29) = pts_value(j,29) / numprocs
 
     ENDDO
@@ -286,4 +308,6 @@
 #endif
 
+    DEALLOCATE( pts_value, pts_value_l )
+
     CALL cpu_log( log_point(36), 'data_output_ptseries','stop', 'nobarrier' )
 
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 824)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 825)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! wang_collision_kernel renamed wang_kernel
 !
 ! Former revisions:
@@ -164,5 +164,5 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -254,5 +254,5 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -268,5 +268,5 @@
 !
 !--    Boundary condition for dissipation
-       IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -372,5 +372,5 @@
 !
 !--    Store dissipation if needed for calculating the sgs particle velocities
-       IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
           DO  k = nzb_s_inner(j,i)+1, nzt
              diss(k,j,i) = dissipation(k)
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 824)
+++ palm/trunk/SOURCE/header.f90	(revision 825)
@@ -1478,4 +1478,9 @@
        WRITE ( io, 433 )
        IF ( curvature_solution_effects )  WRITE ( io, 434 )
+       IF ( collision_kernel /= 'none' )  THEN
+          WRITE ( io, 435 )  TRIM( collision_kernel )
+       ELSE
+          WRITE ( io, 436 )
+       ENDIF
     ENDIF
 
@@ -1955,4 +1960,6 @@
 434 FORMAT ('    Curvature and solution effecs are considered for growth of', &
                  ' droplets < 1.0E-6 m')
+435 FORMAT ('    Droplet collision is handled by ',A,'-kernel')
+436 FORMAT ('    Droplet collision is switched off')
 450 FORMAT (//' LES / Turbulence quantities:'/ &
               ' ---------------------------'/)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 824)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 825)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! wang_collision_kernel renamed wang_kernel
 !
 ! Former revisions:
@@ -334,5 +334,5 @@
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
 !-- particle velocities
-    IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
+    IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
        ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ELSE
Index: palm/trunk/SOURCE/lpm_collision_kernels.f90
===================================================================
--- palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 825)
+++ palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 825)
@@ -0,0 +1,759 @@
+MODULE lpm_collision_kernels_mod
+
+!------------------------------------------------------------------------------!
+! Current revisions:
+! -----------------
+! routine renamed from wang_kernel to lpm_collision_kernels,
+! turbulence_effects on collision replaced by wang_kernel
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+! 799 2011-12-21 17:48:03Z franke
+! speed optimizations and formatting
+! Bugfix: iq=1 is not allowed (routine effic) 
+! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) 
+!
+! 790 2011-11-29 03:11:20Z raasch
+! initial revision
+!
+! Description:
+! ------------
+! This routine calculates the effect of (SGS) turbulence on the collision
+! efficiency of droplets.
+! It is based on the original kernel developed by Wang (...)
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE cloud_parameters
+    USE constants
+    USE particle_attributes
+
+    IMPLICIT NONE
+
+    PRIVATE
+
+    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
+
+    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
+
+    REAL, SAVE :: epsilon, eps2, urms, urms2
+
+    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
+    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
+
+!
+!-- Public interfaces
+    INTERFACE colker
+       MODULE PROCEDURE colker
+    END INTERFACE colker
+
+    INTERFACE effic
+       MODULE PROCEDURE effic
+    END INTERFACE effic
+
+    INTERFACE fallg
+       MODULE PROCEDURE fallg
+    END INTERFACE fallg
+
+    INTERFACE phi
+       MODULE PROCEDURE phi
+    END INTERFACE phi
+
+    INTERFACE turbsd
+       MODULE PROCEDURE turbsd
+    END INTERFACE turbsd
+
+    INTERFACE turb_enhance_eff
+       MODULE PROCEDURE turb_enhance_eff
+    END INTERFACE turb_enhance_eff
+
+    INTERFACE zhi
+       MODULE PROCEDURE zhi
+    END INTERFACE zhi
+
+    CONTAINS
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of collision kernel
+!------------------------------------------------------------------------------!
+    SUBROUTINE colker( i1, j1, k1, kernel )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE cpulog
+       USE indices
+       USE interfaces
+       USE particle_attributes
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, i1, j, j1, k1
+
+       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
+
+!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
+
+       ip = i1
+       jp = j1
+       kp = k1
+
+       pstart = prt_start_index(kp,jp,ip)
+       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
+       psum   = prt_count(kp,jp,ip)
+
+       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
+
+       IF ( wang_kernel )  THEN
+
+          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
+
+          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
+          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
+
+          IF ( epsilon <= 0.001 )  THEN
+
+             CALL fallg
+             CALL effic
+
+             DO  j = pstart, pend
+                DO  i =  pstart, pend
+                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
+                                        particles(i)%radius * 100.0 )**2 &
+                                    * ec(i,j) * ABS( winf(j) - winf(i) )
+                ENDDO
+             ENDDO
+
+          ELSE
+
+             CALL turbsd
+             CALL turb_enhance_eff
+             CALL effic
+
+             DO  j = pstart, pend, 1
+                DO  i =  pstart, pend, 1
+                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
+                ENDDO
+             ENDDO
+
+          ENDIF
+
+          DEALLOCATE(gck, ecf)
+
+       ELSE
+
+!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
+          CALL fallg
+!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
+!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
+          CALL effic
+!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
+
+          DO  j = pstart, pend
+             DO  i =  pstart, pend
+                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
+                                     particles(i)%radius * 100.0 )**2 &
+                                 * ec(i,j) * ABS( winf(j) - winf(i) )
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+       DEALLOCATE( ec, winf )
+
+!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
+
+    END SUBROUTINE colker
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of w, g and gck
+!------------------------------------------------------------------------------!
+    SUBROUTINE turbsd
+! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
+! of droplets, air density, air viscosity, turbulent dissipation rate,
+! taylor microscale reynolds number, gravitational acceleration
+
+       USE constants
+       USE cloud_parameters
+       USE particle_attributes
+       USE arrays_3d
+
+       IMPLICIT NONE
+
+       REAL :: Relamda, &
+               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
+               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
+               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
+               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
+
+       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
+
+       REAL, DIMENSION(pstart:pend) :: St, tau
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       INTEGER :: i, j
+
+!
+!--   Initial assignment of constants
+       IF ( first )  THEN
+
+          first = .FALSE.
+          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
+          airdens   = 1.2250   !air density in mg/cm**3
+          waterdens = 1000.0   !water density in mg/cm**3
+          gravity   = 980.6650 !in cm/s**2
+          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
+
+       ENDIF   
+
+       Relamda = urms**2*sqrt(15.0/epsilon/anu)
+
+       Tl = urms**2/epsilon       !in s
+       Lf = 0.5 * (urms**3)/epsilon !in cm
+
+       tauk = (anu/epsilon)**0.5       !in s
+       eta  = (anu**3/epsilon)**0.25  !in cm
+       vk   = eta/tauk                   !in cm/s
+
+       ao = (11.+7.*Relamda)/(205.+Relamda)
+
+       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
+
+       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
+
+       CALL fallg !gives winf in cm/s
+
+       DO i = pstart, pend
+          tau(i) = winf(i)/gravity    !in s
+          St(i)  = tau(i)/tauk
+       ENDDO
+
+!***** TO CALCULATE wr ********************************
+!from Aayala 2008b, page 38f
+
+       z  = tt/Tl
+       be = sqrt(2.0)*lambda/Lf
+
+       bbb = sqrt(1.0-2.0*be**2)
+       d1 = (1.+bbb)/2.0/bbb
+       e1 = Lf*(1.0+bbb)/2.0   !in cm
+       d2 = (1.0-bbb)/2.0/bbb
+       e2 = Lf*(1.0-bbb)/2.0   !in cm
+
+       ccc = sqrt(1.0-2.0*z**2)
+       b1 = (1.+ccc)/2./ccc
+       c1 = Tl*(1.+ccc)/2.   !in s
+       b2 = (1.-ccc)/2./ccc
+       c2 = Tl*(1.-ccc)/2.   !in s
+
+       DO i = pstart, pend
+          v1 = winf(i)         !in cm/s
+          t1 = tau(i)         !in s
+
+          DO j = pstart,i
+             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
+             v2  = winf(j)         !in cm/s
+             t2  = tau(j)         !in s
+
+             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
+                     - b1*d2*PHI(c1,e2,v1,t1) &
+                     - b2*d1*PHI(c2,e1,v1,t1) &
+                     + b2*d2*PHI(c2,e2,v1,t1)
+             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
+
+             vrms1xy= sqrt(v1xysq)             !in cm/s
+
+             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
+                     - b1*d2*PHI(c1,e2,v2,t2) &
+                     - b2*d1*PHI(c2,e1,v2,t2) &
+                     + b2*d2*PHI(c2,e2,v2,t2)
+             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
+
+              vrms2xy= sqrt(v2xysq)             !in cm/s
+
+             IF(winf(i).ge.winf(j)) THEN
+                v1 = winf(i)
+                t1 = tau(i)
+                v2 = winf(j)
+                t2 = tau(j)
+             ELSE
+                v1 = winf(j)
+                t1 = tau(j)
+                v2 = winf(i)
+                t2 = tau(i)
+             ENDIF
+
+             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
+                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
+                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
+                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
+             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
+             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
+
+             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
+
+              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
+
+              wrgrav2=pi/8.*(winf(j)-winf(i))**2
+
+              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
+
+
+!***** TO CALCULATE gr ********************************
+
+             IF(St(j).gt.St(i)) THEN
+                SSt = St(j)
+             ELSE
+                SSt = St(i)
+             ENDIF
+
+             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
+                  -4.2942*SSt**2 + 5.3406*SSt
+
+             IF(XX.le.0.0) XX = 0.0
+
+             YY = 0.1886*exp(20.306/Relamda)
+
+             c1_gr = XX/(gravity/(vk/tauk))**YY
+
+             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
+             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
+             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
+
+             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
+             IF (grFIN.lt.1.0) grFIN = 1.0
+
+             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
+             gck(j,i) = gck(i,j)
+
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE TurbSD
+
+!------------------------------------------------------------------------------!
+! PHI as a function
+!------------------------------------------------------------------------------!
+  REAL FUNCTION PHI(a,b,vsett,tau0)
+
+       IMPLICIT NONE
+
+       REAL :: a, aa1, b, vsett, tau0
+
+       aa1 = 1./tau0 + 1./a + vsett/b
+
+       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
+
+    END FUNCTION PHI
+
+!------------------------------------------------------------------------------!
+! ZETA as a function
+!------------------------------------------------------------------------------!
+  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
+
+       IMPLICIT NONE
+
+       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
+
+        aa1 = vsett2/b - 1./tau2 - 1./a
+        aa2 = vsett1/b + 1./tau1 + 1./a
+        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
+        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
+        aa5 = vsett2/b + 1./tau2 + 1./a
+        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
+        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
+             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
+             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
+                                                           * 1./2./b/aa3             ! in s**2
+
+    END FUNCTION ZHI
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of terminal velocity winf
+!------------------------------------------------------------------------------!
+   SUBROUTINE fallg
+
+       USE constants
+       USE cloud_parameters
+       USE particle_attributes
+       USE arrays_3d
+
+      IMPLICIT NONE
+
+      INTEGER :: i, j
+
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
+
+      REAL :: bond, x, xrey, y
+
+      REAL, DIMENSION(1:7), SAVE  :: b
+      REAL, DIMENSION(1:6), SAVE  :: c
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
+             -0.578878e-3,0.855176e-4,-0.327815e-5/)
+         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
+             -0.542819e-1, 0.238449e-2/)
+
+         eta   = 1.818e-4     !in poise = g/(cm s)
+         xlamb = 6.62e-6    !in cm
+
+         rhow = 1.0       !in g/cm**3
+         rhoa = 1.225e-3    !in g/cm**3
+
+         grav = 980.665     !in cm/s**2
+         cunh = 1.257 * xlamb !in cm
+         t0   = 273.15        !in K
+         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
+         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
+         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
+         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
+         py   = phy**(1.0/6.0)
+
+      ENDIF
+
+!particle radius has to be in cm
+      DO j = pstart, pend
+
+         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
+
+            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
+
+         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
+
+            x = log(stb*(particles(j)%radius*100.0)**3)
+            y = 0.0
+
+            DO i = 1, 7
+               y = y + b(i) * (x**(i-1))
+            ENDDO
+
+            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
+            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
+
+         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
+
+            IF (particles(j)%radius*100.0.gt.0.35)  THEN
+               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
+            ELSE
+               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
+            ENDIF
+
+            x = log(16.0*bond*py/3.0)
+            y = 0.0
+
+            DO i = 1, 6
+               y = y + c(i) * (x**(i-1))
+            ENDDO
+
+            xrey = py*exp(y)
+
+            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
+              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
+            ELSE
+               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
+            ENDIF
+
+         ENDIF
+      ENDDO
+      RETURN
+   END SUBROUTINE fallg
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of collision efficencies
+!------------------------------------------------------------------------------!
+
+   SUBROUTINE effic
+
+      USE arrays_3d
+      USE constants
+      USE cloud_parameters
+      USE particle_attributes
+
+!collision efficiencies of hall kernel
+      IMPLICIT NONE
+
+      INTEGER :: i, ir, iq, j, k, kk
+
+      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL ::  ek, particle_radius, pp, qq, rq
+
+      REAL, DIMENSION(1:21), SAVE :: rat
+      REAL, DIMENSION(1:15), SAVE :: r0
+      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
+                 300./)
+         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
+                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
+
+         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
+                        0.001,0.001,0.001,0.001,0.001,0.001/)
+         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
+                        0.050,0.200,0.500,0.770,0.870,0.970/)
+         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
+                        0.430,0.580,0.790,0.930,0.960,1.000/)
+         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
+                        0.640,0.750,0.910,0.970,0.980,1.000/)
+         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
+                        0.770,0.840,0.950,0.970,1.000,1.000/)
+         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
+                        0.840,0.880,0.950,1.000,1.000,1.000/)
+         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
+                        0.870,0.900,0.950,1.000,1.000,1.000/)
+         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
+                        0.890,0.920,1.000,1.000,1.000,1.000/)
+         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
+                        0.900,0.940,1.000,1.000,1.000,1.000/)
+         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
+                        0.920,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
+                        0.930,0.970,1.000,1.000,1.000,1.000/)
+         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
+                        0.950,1.000,1.000,1.000,1.000,1.000/)
+         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
+                        1.000,1.020,1.020,1.020,1.020,1.020/)
+         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
+                        1.030,1.040,1.040,1.040,1.040,1.040/)
+         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
+                        1.700,2.300,2.300,2.300,2.300,2.300/)
+         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
+                        3.000,4.000,4.000,4.000,4.000,4.000/)
+      ENDIF
+
+!
+!--   Calculate the radius class index of particles with respect to array r
+      ALLOCATE( ira(pstart:pend) )
+      DO  j = pstart, pend
+         particle_radius = particles(j)%radius * 1.0E6
+         DO  k = 1, 15
+            IF ( particle_radius < r0(k) )  THEN
+               ira(j) = k
+               EXIT
+            ENDIF
+         ENDDO
+         IF ( particle_radius >= r0(15) )  ira(j) = 16
+      ENDDO
+
+!
+!--   Two-dimensional linear interpolation of the collision efficiency.
+!--   Radius has to be in Âµm
+      DO  j = pstart, pend
+         DO  i = pstart, j
+
+            ir = ira(j)
+
+            rq = particles(i)%radius / particles(j)%radius
+
+!            DO kk = 2, 21
+!               IF ( rq <= rat(kk) )  THEN
+!                  iq = kk
+!                  EXIT
+!               ENDIF
+!            ENDDO
+
+            iq = INT( rq * 20 ) + 1
+            iq = MAX(iq , 2)
+
+            IF ( ir < 16 )  THEN
+
+               IF ( ir >= 2 )  THEN
+                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
+                       ( r0(ir) - r0(ir-1) )
+                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
+                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
+                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
+                            + pp * qq * ecoll(ir,iq)
+               ELSE
+                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
+                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
+               ENDIF
+
+            ELSE
+               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
+               ec(j,i) = MIN( ek, 1.0 )
+            ENDIF
+
+            ec(i,j) = ec(j,i)
+            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
+
+         ENDDO
+      ENDDO
+
+      DEALLOCATE( ira )
+
+   END SUBROUTINE effic
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of enhancement factor collision efficencies
+!------------------------------------------------------------------------------!
+   SUBROUTINE turb_enhance_eff
+
+       USE constants
+       USE cloud_parameters
+       USE particle_attributes
+       USE arrays_3d
+
+      IMPLICIT NONE
+
+      INTEGER :: i, ik, ir, iq, j, k, kk
+
+      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
+
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL, DIMENSION(1:11), SAVE :: rat
+      REAL, DIMENSION(1:7), SAVE  :: r0
+      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+
+         r0  = (/10., 20., 30.,40., 50., 60.,100./)
+         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
+
+! 100 cm^2/s^3
+         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
+         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
+         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
+         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
+         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
+         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
+         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
+         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
+         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
+         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
+         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
+
+! 400 cm^2/s^3
+         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
+         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
+         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
+         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
+         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
+         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
+         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
+         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
+         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
+         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
+         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
+
+      ENDIF
+
+!
+!--   Calculate the radius class index of particles with respect to array r
+      ALLOCATE( ira(pstart:pend) )
+
+      DO  j = pstart, pend
+         particle_radius = particles(j)%radius * 1.0E6
+         DO  k = 1, 7
+            IF ( particle_radius < r0(k) )  THEN
+               ira(j) = k
+               EXIT
+            ENDIF
+         ENDDO
+         IF ( particle_radius >= r0(7) )  ira(j) = 8
+      ENDDO
+
+! two-dimensional linear interpolation of the collision efficiency
+      DO j =  pstart, pend
+         DO i = pstart, j
+
+            ir = ira(j)
+
+            rq = particles(i)%radius/particles(j)%radius
+
+            DO kk = 2, 11
+               IF ( rq <= rat(kk) )  THEN
+                  iq = kk
+                  EXIT
+               ENDIF
+            ENDDO
+
+! 0  cm2/s3
+            y1 = 1.0
+! 100 cm2/s3, 400 cm2/s3
+            IF (ir.lt.8) THEN
+               IF (ir.ge.2) THEN
+                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
+                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
+                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
+                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
+                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
+                            pp*qq*ecoll_100(ir,iq)
+                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
+                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
+                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
+                            pp*qq*ecoll_400(ir,iq)
+               ELSE
+                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
+                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
+                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
+               ENDIF
+            ELSE
+            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
+            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
+            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
+         ENDIF
+! linear interpolation
+! dissipation rate in cm2/s3
+         x1 = 0.0
+         x2 = 100.0
+         x3 = 400.0
+
+         IF (epsilon.le.100.) THEN
+            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
+                    +  (epsilon-0.)/(100.-0.) * y2
+         ELSE IF(epsilon.le.600.)THEN
+            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
+                    +  (epsilon-100.)/(400.-100.) * y3
+
+         ELSE
+            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
+                    +  (600.-100.)/(400.-100.) * y3
+         ENDIF
+
+         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
+
+         ecf(i,j)=ecf(j,i)
+      ENDDO
+   ENDDO
+
+   RETURN
+   END SUBROUTINE turb_enhance_eff
+
+ END MODULE lpm_collision_kernels_mod
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 824)
+++ palm/trunk/SOURCE/modules.f90	(revision 825)
@@ -6,6 +6,8 @@
 ! +bfactor, curvature_solution_effects, eps_ros, molecular_weight_of_water,
 ! vanthoff, -b_cond in cloud_parameters,
-! dopts_num increased to 27, particle attributes speed_x|y|z_sgs renamed
+! dopts_num increased to 29, particle attributes speed_x|y|z_sgs renamed
 ! rvar1|2|3
+! wang_collision_kernel and turbulence_effects_on_collision replaced by
+! collision_kernel, hall_kernel, palm_kernel, wang_kernel
 !
 ! Former revisions:
@@ -1047,5 +1049,5 @@
 #endif
 
-    INTEGER, PARAMETER ::  dopr_norm_num = 7, dopts_num = 27, dots_max = 100
+    INTEGER, PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
 
     INTEGER ::  dots_num = 23
@@ -1062,14 +1064,14 @@
           (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
              'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
-             'w_up   ', 'w_down ', 'radius ', 'npt_max', 'npt_min', 'x*2    ', &
-             'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', 'w*2    ', 'u"2    ', &
-             'v"2    ', 'w"2    ', 'npt*2  ' /)
+             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
+             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
+             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
 
     CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = &
           (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
              'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
-             'm/s    ', 'm/s    ', 'm      ', 'number ', 'number ', 'm2     ', &
-             'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', &
-             'm2/s2  ', 'm2/s2  ', 'number2' /)
+             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
+             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
+             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
 
     CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_label = &
@@ -1184,5 +1186,6 @@
 
     CHARACTER (LEN=15)  ::  bc_par_lr = 'cyclic',  bc_par_ns = 'cyclic', &
-                            bc_par_b  = 'reflect', bc_par_t  = 'absorb'
+                            bc_par_b  = 'reflect', bc_par_t  = 'absorb', &
+                            collision_kernel = 'none'
 
 #if defined( __parallel )
@@ -1206,11 +1209,10 @@
     INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  prt_count, prt_start_index
 
-    LOGICAL ::  turbulence_effects_on_collision = .FALSE.,                     &
+    LOGICAL ::  hall_kernel = .FALSE., palm_kernel = .FALSE.,                  &
                 particle_advection = .FALSE., random_start_position = .FALSE., &
                 read_particles_from_restartfile = .TRUE.,                      &
                 uniform_particles = .TRUE., use_particle_tails = .FALSE.,      &
                 use_sgs_for_particles = .FALSE.,                               &
-                wang_collision_kernel = .FALSE.,                               &
-                write_particle_statistics = .FALSE.
+                wang_kernel = .FALSE., write_particle_statistics = .FALSE.
 
     LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
Index: palm/trunk/SOURCE/package_parin.f90
===================================================================
--- palm/trunk/SOURCE/package_parin.f90	(revision 824)
+++ palm/trunk/SOURCE/package_parin.f90	(revision 825)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! wang_collision_kernel and turbulence_effects_on_collision in particles_par
+! replaced by collision_kernel
 !
 ! Former revisions:
@@ -71,4 +72,5 @@
 
     NAMELIST /particles_par/      bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,    &
+                                  collision_kernel,                            &
                                   density_ratio, radius, dt_dopts,             &
                                   dt_min_part, dt_prel, dt_sort_particles,     &
@@ -87,8 +89,6 @@
                                   read_particles_from_restartfile,             &
                                   skip_particles_for_tail,                     &
-                                  turbulence_effects_on_collision,             &
                                   use_particle_tails, use_sgs_for_particles,   &
                                   vertical_particle_advection,                 &
-                                  wang_collision_kernel,                       &
                                   write_particle_statistics
     NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 824)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 825)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! wang_collision_kernel renamed wang_kernel
 !
 ! Former revisions:
@@ -218,5 +218,5 @@
              CALL exchange_horiz( ql_vp, nbgp )
           ENDIF
-          IF ( wang_collision_kernel )  CALL exchange_horiz( diss, nbgp )
+          IF ( wang_kernel )  CALL exchange_horiz( diss, nbgp )
 
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
Index: palm/trunk/SOURCE/wang_kernel.f90
===================================================================
--- palm/trunk/SOURCE/wang_kernel.f90	(revision 824)
+++ 	(revision )
@@ -1,758 +1,0 @@
-MODULE wang_kernel_mod
-
-!------------------------------------------------------------------------------!
-! Current revisions:
-! -----------------
-!
-!
-! Former revisions:
-! -----------------
-! $Id$
-!
-! 799 2011-12-21 17:48:03Z franke
-! speed optimizations and formatting
-! Bugfix: iq=1 is not allowed (routine effic) 
-! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) 
-!
-! 790 2011-11-29 03:11:20Z raasch
-! initial revision
-!
-! Description:
-! ------------
-! This routine calculates the effect of (SGS) turbulence on the collision
-! efficiency of droplets.
-! It is based on the original kernel developed by Wang (...)
-!------------------------------------------------------------------------------!
-
-    USE arrays_3d
-    USE cloud_parameters
-    USE constants
-    USE particle_attributes
-
-    IMPLICIT NONE
-
-    PRIVATE
-
-    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
-
-    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
-
-    REAL, SAVE :: epsilon, eps2, urms, urms2
-
-    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
-    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
-
-!
-!-- Public interfaces
-    INTERFACE colker
-       MODULE PROCEDURE colker
-    END INTERFACE colker
-
-    INTERFACE effic
-       MODULE PROCEDURE effic
-    END INTERFACE effic
-
-    INTERFACE fallg
-       MODULE PROCEDURE fallg
-    END INTERFACE fallg
-
-    INTERFACE phi
-       MODULE PROCEDURE phi
-    END INTERFACE phi
-
-    INTERFACE turbsd
-       MODULE PROCEDURE turbsd
-    END INTERFACE turbsd
-
-    INTERFACE turb_enhance_eff
-       MODULE PROCEDURE turb_enhance_eff
-    END INTERFACE turb_enhance_eff
-
-    INTERFACE zhi
-       MODULE PROCEDURE zhi
-    END INTERFACE zhi
-
-    CONTAINS
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision kernel
-!------------------------------------------------------------------------------!
-    SUBROUTINE colker( i1, j1, k1, kernel )
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE constants
-       USE cpulog
-       USE indices
-       USE interfaces
-       USE particle_attributes
-
-       IMPLICIT NONE
-
-       INTEGER ::  i, i1, j, j1, k1
-
-       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
-
-       ip = i1
-       jp = j1
-       kp = k1
-
-       pstart = prt_start_index(kp,jp,ip)
-       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
-       psum   = prt_count(kp,jp,ip)
-
-       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
-
-       IF ( turbulence_effects_on_collision )  THEN
-
-          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
-
-          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
-          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
-
-          IF ( epsilon <= 0.001 )  THEN
-
-             CALL fallg
-             CALL effic
-
-             DO  j = pstart, pend
-                DO  i =  pstart, pend
-                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                        particles(i)%radius * 100.0 )**2 &
-                                    * ec(i,j) * ABS( winf(j) - winf(i) )
-                ENDDO
-             ENDDO
-
-          ELSE
-
-             CALL turbsd
-             CALL turb_enhance_eff
-             CALL effic
-
-             DO  j = pstart, pend, 1
-                DO  i =  pstart, pend, 1
-                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
-                ENDDO
-             ENDDO
-
-          ENDIF
-
-          DEALLOCATE(gck, ecf)
-
-       ELSE
-
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
-          CALL fallg
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
-          CALL effic
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
-
-          DO  j = pstart, pend
-             DO  i =  pstart, pend
-                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                     particles(i)%radius * 100.0 )**2 &
-                                 * ec(i,j) * ABS( winf(j) - winf(i) )
-             ENDDO
-          ENDDO
-
-       ENDIF
-
-       DEALLOCATE( ec, winf )
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
-
-    END SUBROUTINE colker
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of w, g and gck
-!------------------------------------------------------------------------------!
-    SUBROUTINE turbsd
-! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
-! of droplets, air density, air viscosity, turbulent dissipation rate,
-! taylor microscale reynolds number, gravitational acceleration
-
-       USE constants
-       USE cloud_parameters
-       USE particle_attributes
-       USE arrays_3d
-
-       IMPLICIT NONE
-
-       REAL :: Relamda, &
-               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
-               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
-               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
-               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
-
-       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
-
-       REAL, DIMENSION(pstart:pend) :: St, tau
-
-       LOGICAL, SAVE ::  first = .TRUE.
-
-       INTEGER :: i, j
-
-!
-!--   Initial assignment of constants
-       IF ( first )  THEN
-
-          first = .FALSE.
-          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
-          airdens   = 1.2250   !air density in mg/cm**3
-          waterdens = 1000.0   !water density in mg/cm**3
-          gravity   = 980.6650 !in cm/s**2
-          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
-
-       ENDIF   
-
-       Relamda = urms**2*sqrt(15.0/epsilon/anu)
-
-       Tl = urms**2/epsilon       !in s
-       Lf = 0.5 * (urms**3)/epsilon !in cm
-
-       tauk = (anu/epsilon)**0.5       !in s
-       eta  = (anu**3/epsilon)**0.25  !in cm
-       vk   = eta/tauk                   !in cm/s
-
-       ao = (11.+7.*Relamda)/(205.+Relamda)
-
-       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
-
-       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
-
-       CALL fallg !gives winf in cm/s
-
-       DO i = pstart, pend
-          tau(i) = winf(i)/gravity    !in s
-          St(i)  = tau(i)/tauk
-       ENDDO
-
-!***** TO CALCULATE wr ********************************
-!from Aayala 2008b, page 38f
-
-       z  = tt/Tl
-       be = sqrt(2.0)*lambda/Lf
-
-       bbb = sqrt(1.0-2.0*be**2)
-       d1 = (1.+bbb)/2.0/bbb
-       e1 = Lf*(1.0+bbb)/2.0   !in cm
-       d2 = (1.0-bbb)/2.0/bbb
-       e2 = Lf*(1.0-bbb)/2.0   !in cm
-
-       ccc = sqrt(1.0-2.0*z**2)
-       b1 = (1.+ccc)/2./ccc
-       c1 = Tl*(1.+ccc)/2.   !in s
-       b2 = (1.-ccc)/2./ccc
-       c2 = Tl*(1.-ccc)/2.   !in s
-
-       DO i = pstart, pend
-          v1 = winf(i)         !in cm/s
-          t1 = tau(i)         !in s
-
-          DO j = pstart,i
-             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
-             v2  = winf(j)         !in cm/s
-             t2  = tau(j)         !in s
-
-             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
-                     - b1*d2*PHI(c1,e2,v1,t1) &
-                     - b2*d1*PHI(c2,e1,v1,t1) &
-                     + b2*d2*PHI(c2,e2,v1,t1)
-             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
-
-             vrms1xy= sqrt(v1xysq)             !in cm/s
-
-             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
-                     - b1*d2*PHI(c1,e2,v2,t2) &
-                     - b2*d1*PHI(c2,e1,v2,t2) &
-                     + b2*d2*PHI(c2,e2,v2,t2)
-             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
-
-              vrms2xy= sqrt(v2xysq)             !in cm/s
-
-             IF(winf(i).ge.winf(j)) THEN
-                v1 = winf(i)
-                t1 = tau(i)
-                v2 = winf(j)
-                t2 = tau(j)
-             ELSE
-                v1 = winf(j)
-                t1 = tau(j)
-                v2 = winf(i)
-                t2 = tau(i)
-             ENDIF
-
-             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
-                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
-                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
-                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
-             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
-             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
-
-             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
-
-              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
-
-              wrgrav2=pi/8.*(winf(j)-winf(i))**2
-
-              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
-
-
-!***** TO CALCULATE gr ********************************
-
-             IF(St(j).gt.St(i)) THEN
-                SSt = St(j)
-             ELSE
-                SSt = St(i)
-             ENDIF
-
-             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
-                  -4.2942*SSt**2 + 5.3406*SSt
-
-             IF(XX.le.0.0) XX = 0.0
-
-             YY = 0.1886*exp(20.306/Relamda)
-
-             c1_gr = XX/(gravity/(vk/tauk))**YY
-
-             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
-             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
-             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
-
-             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
-             IF (grFIN.lt.1.0) grFIN = 1.0
-
-             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
-             gck(j,i) = gck(i,j)
-
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE TurbSD
-
-!------------------------------------------------------------------------------!
-! PHI as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION PHI(a,b,vsett,tau0)
-
-       IMPLICIT NONE
-
-       REAL :: a, aa1, b, vsett, tau0
-
-       aa1 = 1./tau0 + 1./a + vsett/b
-
-       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
-
-    END FUNCTION PHI
-
-!------------------------------------------------------------------------------!
-! ZETA as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
-
-       IMPLICIT NONE
-
-       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
-
-        aa1 = vsett2/b - 1./tau2 - 1./a
-        aa2 = vsett1/b + 1./tau1 + 1./a
-        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
-        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
-        aa5 = vsett2/b + 1./tau2 + 1./a
-        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
-        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
-             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
-             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
-                                                           * 1./2./b/aa3             ! in s**2
-
-    END FUNCTION ZHI
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of terminal velocity winf
-!------------------------------------------------------------------------------!
-   SUBROUTINE fallg
-
-       USE constants
-       USE cloud_parameters
-       USE particle_attributes
-       USE arrays_3d
-
-      IMPLICIT NONE
-
-      INTEGER :: i, j
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
-
-      REAL :: bond, x, xrey, y
-
-      REAL, DIMENSION(1:7), SAVE  :: b
-      REAL, DIMENSION(1:6), SAVE  :: c
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
-             -0.578878e-3,0.855176e-4,-0.327815e-5/)
-         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
-             -0.542819e-1, 0.238449e-2/)
-
-         eta   = 1.818e-4     !in poise = g/(cm s)
-         xlamb = 6.62e-6    !in cm
-
-         rhow = 1.0       !in g/cm**3
-         rhoa = 1.225e-3    !in g/cm**3
-
-         grav = 980.665     !in cm/s**2
-         cunh = 1.257 * xlamb !in cm
-         t0   = 273.15        !in K
-         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
-         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
-         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
-         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
-         py   = phy**(1.0/6.0)
-
-      ENDIF
-
-!particle radius has to be in cm
-      DO j = pstart, pend
-
-         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
-
-            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
-
-         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
-
-            x = log(stb*(particles(j)%radius*100.0)**3)
-            y = 0.0
-
-            DO i = 1, 7
-               y = y + b(i) * (x**(i-1))
-            ENDDO
-
-            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
-            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
-
-         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
-
-            IF (particles(j)%radius*100.0.gt.0.35)  THEN
-               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
-            ELSE
-               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
-            ENDIF
-
-            x = log(16.0*bond*py/3.0)
-            y = 0.0
-
-            DO i = 1, 6
-               y = y + c(i) * (x**(i-1))
-            ENDDO
-
-            xrey = py*exp(y)
-
-            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
-              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
-            ELSE
-               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
-            ENDIF
-
-         ENDIF
-      ENDDO
-      RETURN
-   END SUBROUTINE fallg
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision efficencies
-!------------------------------------------------------------------------------!
-
-   SUBROUTINE effic
-
-      USE arrays_3d
-      USE constants
-      USE cloud_parameters
-      USE particle_attributes
-
-!collision efficiencies of hall kernel
-      IMPLICIT NONE
-
-      INTEGER :: i, ir, iq, j, k, kk
-
-      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL ::  ek, particle_radius, pp, qq, rq
-
-      REAL, DIMENSION(1:21), SAVE :: rat
-      REAL, DIMENSION(1:15), SAVE :: r0
-      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
-                 300./)
-         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
-                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
-
-         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
-                        0.001,0.001,0.001,0.001,0.001,0.001/)
-         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
-                        0.050,0.200,0.500,0.770,0.870,0.970/)
-         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
-                        0.430,0.580,0.790,0.930,0.960,1.000/)
-         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
-                        0.640,0.750,0.910,0.970,0.980,1.000/)
-         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
-                        0.770,0.840,0.950,0.970,1.000,1.000/)
-         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
-                        0.840,0.880,0.950,1.000,1.000,1.000/)
-         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
-                        0.870,0.900,0.950,1.000,1.000,1.000/)
-         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
-                        0.890,0.920,1.000,1.000,1.000,1.000/)
-         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
-                        0.900,0.940,1.000,1.000,1.000,1.000/)
-         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
-                        0.920,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
-                        0.930,0.970,1.000,1.000,1.000,1.000/)
-         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
-                        0.950,1.000,1.000,1.000,1.000,1.000/)
-         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
-                        1.000,1.020,1.020,1.020,1.020,1.020/)
-         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
-                        1.030,1.040,1.040,1.040,1.040,1.040/)
-         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
-                        1.700,2.300,2.300,2.300,2.300,2.300/)
-         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
-                        3.000,4.000,4.000,4.000,4.000,4.000/)
-      ENDIF
-
-!
-!--   Calculate the radius class index of particles with respect to array r
-      ALLOCATE( ira(pstart:pend) )
-      DO  j = pstart, pend
-         particle_radius = particles(j)%radius * 1.0E6
-         DO  k = 1, 15
-            IF ( particle_radius < r0(k) )  THEN
-               ira(j) = k
-               EXIT
-            ENDIF
-         ENDDO
-         IF ( particle_radius >= r0(15) )  ira(j) = 16
-      ENDDO
-
-!
-!--   Two-dimensional linear interpolation of the collision efficiency.
-!--   Radius has to be in Âµm
-      DO  j = pstart, pend
-         DO  i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius / particles(j)%radius
-
-!            DO kk = 2, 21
-!               IF ( rq <= rat(kk) )  THEN
-!                  iq = kk
-!                  EXIT
-!               ENDIF
-!            ENDDO
-
-            iq = INT( rq * 20 ) + 1
-            iq = MAX(iq , 2)
-
-            IF ( ir < 16 )  THEN
-
-               IF ( ir >= 2 )  THEN
-                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
-                       ( r0(ir) - r0(ir-1) )
-                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
-                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
-                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
-                            + pp * qq * ecoll(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
-               ENDIF
-
-            ELSE
-               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
-               ec(j,i) = MIN( ek, 1.0 )
-            ENDIF
-
-            ec(i,j) = ec(j,i)
-            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
-
-         ENDDO
-      ENDDO
-
-      DEALLOCATE( ira )
-
-   END SUBROUTINE effic
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of enhancement factor collision efficencies
-!------------------------------------------------------------------------------!
-   SUBROUTINE turb_enhance_eff
-
-       USE constants
-       USE cloud_parameters
-       USE particle_attributes
-       USE arrays_3d
-
-      IMPLICIT NONE
-
-      INTEGER :: i, ik, ir, iq, j, k, kk
-
-      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
-
-      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL, DIMENSION(1:11), SAVE :: rat
-      REAL, DIMENSION(1:7), SAVE  :: r0
-      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-
-         r0  = (/10., 20., 30.,40., 50., 60.,100./)
-         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
-
-! 100 cm^2/s^3
-         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
-         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
-         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
-         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
-         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
-         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
-         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
-         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
-         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
-         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
-         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
-
-! 400 cm^2/s^3
-         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
-         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
-         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
-         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
-         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
-         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
-         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
-         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
-         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
-         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
-         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
-
-      ENDIF
-
-!
-!--   Calculate the radius class index of particles with respect to array r
-      ALLOCATE( ira(pstart:pend) )
-
-      DO  j = pstart, pend
-         particle_radius = particles(j)%radius * 1.0E6
-         DO  k = 1, 7
-            IF ( particle_radius < r0(k) )  THEN
-               ira(j) = k
-               EXIT
-            ENDIF
-         ENDDO
-         IF ( particle_radius >= r0(7) )  ira(j) = 8
-      ENDDO
-
-! two-dimensional linear interpolation of the collision efficiency
-      DO j =  pstart, pend
-         DO i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius/particles(j)%radius
-
-            DO kk = 2, 11
-               IF ( rq <= rat(kk) )  THEN
-                  iq = kk
-                  EXIT
-               ENDIF
-            ENDDO
-
-! 0  cm2/s3
-            y1 = 1.0
-! 100 cm2/s3, 400 cm2/s3
-            IF (ir.lt.8) THEN
-               IF (ir.ge.2) THEN
-                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
-                            pp*qq*ecoll_100(ir,iq)
-                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
-                            pp*qq*ecoll_400(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
-                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
-               ENDIF
-            ELSE
-            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
-            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
-         ENDIF
-! linear interpolation
-! dissipation rate in cm2/s3
-         x1 = 0.0
-         x2 = 100.0
-         x3 = 400.0
-
-         IF (epsilon.le.100.) THEN
-            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
-                    +  (epsilon-0.)/(100.-0.) * y2
-         ELSE IF(epsilon.le.600.)THEN
-            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
-                    +  (epsilon-100.)/(400.-100.) * y3
-
-         ELSE
-            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
-                    +  (600.-100.)/(400.-100.) * y3
-         ENDIF
-
-         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
-
-         ecf(i,j)=ecf(j,i)
-      ENDDO
-   ENDDO
-
-   RETURN
-   END SUBROUTINE turb_enhance_eff
-
- END MODULE wang_kernel_mod
