Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 789)
+++ palm/trunk/SOURCE/Makefile	(revision 790)
@@ -4,4 +4,5 @@
 # Current revisions:
 # ------------------
+# +wang_kernel
 #
 # Former revisions:
@@ -104,5 +105,6 @@
 	user_particle_attributes.f90 user_read_restart_data.f90 \
 	user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
-	write_3d_binary.f90 write_compressed.f90 write_var_list.f90
+	wang_kernel.f90 write_3d_binary.f90 write_compressed.f90 \
+	write_var_list.f90
 
 OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
@@ -146,5 +148,6 @@
 	user_module.o user_parin.o user_particle_attributes.o \
 	user_read_restart_data.o user_spectra.o user_statistics.o \
-	wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o
+	wall_fluxes.o wang_kernel.o write_3d_binary.o write_compressed.o \
+	write_var_list.o
 
 CC = cc
@@ -175,5 +178,5 @@
 
 
-advec_particles.o: modules.o random_function.o
+advec_particles.o: modules.o random_function.o wang_kernel.o
 advec_s_bc.o: modules.o
 advec_s_pw.o: modules.o
@@ -322,4 +325,5 @@
 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/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 789)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 790)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! bugfix: calculation of 'pr' must depend on the particle weighting factor
 !
 ! Former revisions:
@@ -303,6 +304,8 @@
                                s_r4 = 0.0
                                DO  n = psi, psi+prt_count(k,j,i)-1
-                                  s_r3 = s_r3 + particles(n)%radius**3
-                                  s_r4 = s_r4 + particles(n)%radius**4
+                                  s_r3 = s_r3 + particles(n)%radius**3 * &
+                                                particles(n)%weight_factor
+                                  s_r4 = s_r4 + particles(n)%radius**4 * &
+                                                particles(n)%weight_factor
                                ENDDO
                                IF ( s_r3 /= 0.0 )  THEN
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 789)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 790)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! bugfix: calculation of 'pr' must depend on the particle weighting factor,
+! nzt+1 replaced by nz_do3d for 'pr'
 !
 ! Former revisions:
@@ -207,11 +209,13 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb, nz_do3d
                          psi = prt_start_index(k,j,i)
                          s_r3 = 0.0
                          s_r4 = 0.0
                          DO  n = psi, psi+prt_count(k,j,i)-1
-                            s_r3 = s_r3 + particles(n)%radius**3
-                            s_r4 = s_r4 + particles(n)%radius**4
+                         s_r3 = s_r3 + particles(n)%radius**3 * &
+                                       particles(n)%weight_factor
+                         s_r4 = s_r4 + particles(n)%radius**4 * &
+                                       particles(n)%weight_factor
                          ENDDO
                          IF ( s_r3 /= 0.0 )  THEN
@@ -227,5 +231,5 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb, nz_do3d
                          local_pf(i,j,k) = tend(k,j,i)
                       ENDDO
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 789)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 790)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! diss is also calculated in case that the Wang kernel is used
 !
 ! Former revisions:
@@ -160,5 +161,5 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -250,5 +251,5 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -264,5 +265,5 @@
 !
 !--    Boundary condition for dissipation
-       IF ( use_sgs_for_particles )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -368,5 +369,5 @@
 !
 !--    Store dissipation if needed for calculating the sgs particle velocities
-       IF ( use_sgs_for_particles )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
           DO  k = nzb_s_inner(j,i)+1, nzt
              diss(k,j,i) = dissipation(k)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 789)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 790)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! diss is also allocated in case that the Wang kernel is used
 !
 ! Former revisions:
@@ -331,5 +331,5 @@
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
 !-- particle velocities
-    IF ( use_sgs_for_particles )  THEN
+    IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
        ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ELSE
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 789)
+++ palm/trunk/SOURCE/modules.f90	(revision 790)
@@ -5,5 +5,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +turbulence_effects_on_collision, wang_collision_kernel
 !
 ! Former revisions:
@@ -1170,8 +1170,10 @@
     INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  prt_count, prt_start_index
 
-    LOGICAL ::  particle_advection = .FALSE., random_start_position = .FALSE., &
+    LOGICAL ::  turbulence_effects_on_collision = .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.
 
Index: palm/trunk/SOURCE/package_parin.f90
===================================================================
--- palm/trunk/SOURCE/package_parin.f90	(revision 789)
+++ palm/trunk/SOURCE/package_parin.f90	(revision 790)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +turbulence_effects_on_collision, wang_collision_kernel in particles_par
 !
 ! Former revisions:
@@ -83,7 +83,9 @@
                                   random_start_position,                       &
                                   read_particles_from_restartfile,             &
-                                  skip_particles_for_tail, use_particle_tails, &
-                                  use_sgs_for_particles,                       &
+                                  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/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 789)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 790)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! bugfix: calculation of 'pr' must depend on the particle weighting factor
 !
 ! Former revisions:
@@ -295,6 +296,8 @@
                       s_r4 = 0.0
                       DO  n = psi, psi+prt_count(k,j,i)-1
-                         s_r3 = s_r3 + particles(n)%radius**3
-                         s_r4 = s_r4 + particles(n)%radius**4
+                         s_r3 = s_r3 + particles(n)%radius**3 * &
+                                       particles(n)%weight_factor
+                         s_r4 = s_r4 + particles(n)%radius**4 * &
+                                       particles(n)%weight_factor
                       ENDDO
                       IF ( s_r3 /= 0.0 )  THEN
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 789)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 790)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! exchange of ghostpoints for array diss
 !
 ! Former revisions:
@@ -208,5 +208,5 @@
              CALL exchange_horiz( prho, nbgp )
           ENDIF
-          IF (humidity  .OR.  passive_scalar) CALL exchange_horiz( q_p, nbgp )
+          IF (humidity  .OR.  passive_scalar)  CALL exchange_horiz( q_p, nbgp )
           IF ( cloud_droplets )  THEN
              CALL exchange_horiz( ql, nbgp )
@@ -215,4 +215,5 @@
              CALL exchange_horiz( ql_vp, nbgp )
           ENDIF
+          IF ( wang_collision_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 790)
+++ palm/trunk/SOURCE/wang_kernel.f90	(revision 790)
@@ -0,0 +1,649 @@
+MODULE wang_kernel_mod
+
+!------------------------------------------------------------------------------!
+! Current revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+!
+! 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 TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker
+
+    INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum
+    REAL, SAVE :: epsilon, urms
+    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf
+    REAL, DIMENSION(:),   ALLOCATABLE, SAVE :: winf
+    REAL, SAVE :: eps2, urms2
+
+!
+!-- Public interfaces
+    INTERFACE TurbSD
+       MODULE PROCEDURE TurbSD
+    END INTERFACE TurbSD
+
+    INTERFACE PHI
+       MODULE PROCEDURE PHI
+    END INTERFACE PHI
+
+    INTERFACE ZHI
+       MODULE PROCEDURE ZHI
+    END INTERFACE ZHI
+
+    INTERFACE turb_enhance_eff
+       MODULE PROCEDURE turb_enhance_eff
+    END INTERFACE turb_enhance_eff
+
+    INTERFACE effic
+       MODULE PROCEDURE effic
+    END INTERFACE effic
+
+    INTERFACE fallg
+       MODULE PROCEDURE fallg
+    END INTERFACE fallg
+
+    INTERFACE colker
+       MODULE PROCEDURE colker
+    END INTERFACE colker
+
+    CONTAINS
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of collision kernel
+!------------------------------------------------------------------------------!
+    SUBROUTINE colker(i1,j1,k1,kernel)
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       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
+
+       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(winf(pstart:pend), ec(pstart:pend,pstart:pend))
+
+       IF ( turbulence_effects_on_collision ) THEN
+
+          ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))
+
+          epsilon = diss(kp,jp,ip) * 10000.0  !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, 1
+                DO i =  pstart, pend, 1
+                   kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * 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 fallg
+          CALL effic
+
+          DO j = pstart, pend, 1
+             DO i =  pstart, pend, 1
+                kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+       DEALLOCATE(winf, ec)
+
+       RETURN
+
+    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 :: airvisc, airdens, waterdens, gravity, anu, 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, DIMENSION(pstart:pend) :: vST, tau, St
+
+       INTEGER :: i, j
+
+       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
+
+       Relamda = urms**2.0*sqrt(15.0/epsilon/anu)
+
+       Tl = urms**2.0/epsilon       !in s
+       Lf = 0.5 * (urms**3.0)/epsilon !in cm
+
+       tauk = (anu/epsilon)**0.5       !in s
+       eta  = (anu**3.0/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
+          vST(i) = winf(i)           !in cm/s
+          tau(i) = vST(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.0)
+       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.0)
+       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 = vST(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  = vST(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.0/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.0/t2    !in cm**2/s**2
+
+              vrms2xy= sqrt(v2xysq)             !in cm/s
+
+             IF(vST(i).ge.vST(j)) THEN
+                v1 = vST(i)
+                t1 = tau(i)
+                v2 = vST(j)
+                t2 = tau(j)
+             ELSE
+                v1 = vST(j)
+                t1 = tau(j)
+                v2 = vST(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.0/tau(i)/tau(j)    !in cm**2/s**2
+
+             wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy  !in cm**2/s**2
+
+              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
+
+              wrgrav2=pi/8.*(vST(j)-vST(i))**2.0
+
+              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.0 + 1.5275*SSt**3.0 &
+                  -4.2942*SSt**2.0 + 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.0
+             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
+             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
+
+             grFIN = ((eta**2.0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.)
+             IF (grFIN.lt.1.0) grFIN = 1.0
+
+             gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN       ! in cm**3/s
+             gck(j,i) = gck(i,j)
+
+          ENDDO
+       ENDDO
+
+
+       RETURN
+    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.0  !in s
+       RETURN
+    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.0 - (1./tau2 + 1./a)**2.0
+        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.0          &
+             + (4./aa4 - 1./aa5**2.0 - 1./aa1**2.0) * vsett2/2./b/aa6      &
+             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
+                                                           * 1./2./b/aa3             ! in s**2
+        RETURN
+    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
+
+      REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py,  &
+              x, y, xrey, bond
+
+      REAL, DIMENSION(1:7)  :: b
+      REAL, DIMENSION(1:6)  :: c
+      REAL, DIMENSION(1:20) :: rat
+      REAL, DIMENSION(1:15) :: r0
+      REAL, DIMENSION(1:15,1:20) :: ecoll
+
+      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.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa))
+      py   = phy**(1.0/6.0)
+
+!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.0)
+            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
+
+            bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma
+
+            IF (particles(j)%radius*100.0.gt.0.35)  THEN
+               bond = grav*(rhow-rhoa) * 0.35 * 0.35/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)
+            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
+
+            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
+              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
+            ENDIF
+
+         ENDIF
+      ENDDO
+      RETURN
+   END SUBROUTINE fallg
+
+!------------------------------------------------------------------------------!
+! SUBROUTINE for calculation of collision efficencies
+!------------------------------------------------------------------------------!
+
+   SUBROUTINE effic
+
+       USE constants
+       USE cloud_parameters
+       USE particle_attributes
+       USE arrays_3d
+
+!collision efficiencies of hall kernel
+      IMPLICIT NONE
+
+      INTEGER :: i, ir, iq, j, k, kk
+      REAL    :: ek, rq, pp, qq
+
+      REAL, DIMENSION(1:21) :: rat
+      REAL, DIMENSION(1:15) :: r0
+      REAL, DIMENSION(1:15,1:21) :: ecoll
+
+      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/)
+
+! two-dimensional linear interpolation of the collision efficiency
+! radius has to be in Âµm
+
+      DO j = pstart, pend
+         DO i = pstart, j
+
+            DO k = 2, 15
+               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
+                  ir = k
+               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
+                  ir = 16
+               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
+                  ir = 1
+               ENDIF
+            ENDDO
+
+            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
+
+            DO kk = 2, 21
+               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
+            ENDDO
+
+            IF (ir.lt.16) 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))
+                  ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) +  &
+                          pp * (1.-qq) * ecoll(ir,iq-1) +         &
+                          qq * (1.-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.-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).lt.1.e-20) stop 99
+         ENDDO
+      ENDDO
+      RETURN
+   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
+
+!collision efficiencies of hall kernel
+      IMPLICIT NONE
+
+      INTEGER :: i, ik, ir, iq, j, k, kk
+      REAL    :: rq, y1, pp, qq, y2, y3, x1, x2, x3
+
+      REAL, DIMENSION(1:11) :: rat
+      REAL, DIMENSION(1:7)  :: r0
+      REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400
+
+      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 /)
+
+! two-dimensional linear interpolation of the collision efficiency
+      DO j =  pstart, pend
+         DO i = pstart, j
+
+            DO k = 2, 7
+               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
+                  ir = k
+               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN
+                  ir = 8
+               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
+                  ir = 1
+               ENDIF
+            ENDDO
+
+            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
+
+            DO kk = 2, 11
+               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
+            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
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 789)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 790)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! bugfix: output of 'ql' also if cloud droplets are used
 !
 ! Former revisions:
@@ -141,5 +142,5 @@
        ENDIF
        WRITE ( 14 )  'q_m                 ';  WRITE ( 14 ) q_m
-       IF ( cloud_physics ) THEN
+       IF ( cloud_physics  .OR.  cloud_droplets )  THEN
           WRITE ( 14 )  'ql                  ';  WRITE ( 14 ) ql
           IF ( ALLOCATED( ql_av ) )  THEN
