Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 1006)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 1007)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix: missing calculation of ql_vp added
 !
 ! Former revisions:
@@ -428,7 +428,32 @@
              CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
                 IF ( av == 0 )  THEN
+                   IF ( simulated_time >= particle_advection_start )  THEN
+                      DO  i = nxl, nxr
+                         DO  j = nys, nyn
+                            DO  k = nzb, nzt+1
+                               psi = prt_start_index(k,j,i)
+                               DO  n = psi, psi+prt_count(k,j,i)-1
+                                  tend(k,j,i) =  tend(k,j,i) + &
+                                                 particles(n)%weight_factor / &
+                                                 prt_count(k,j,i)
+                               ENDDO
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                      CALL exchange_horiz( tend, nbgp )
+                   ELSE
+                      tend = 0.0
+                   END IF
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = tend(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   CALL exchange_horiz( ql_vp_av, nbgp )
                    to_be_resorted => ql_vp
-                ELSE
-                   to_be_resorted => ql_vp_av
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 1006)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 1007)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: missing calculation of ql_vp added
 !
 ! Former revisions:
@@ -295,6 +295,27 @@
           CASE ( 'ql_vp' )
              IF ( av == 0 )  THEN
-                to_be_resorted => ql_vp
-             ELSE
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nz_do3d
+                         psi = prt_start_index(k,j,i)
+                         DO  n = psi, psi+prt_count(k,j,i)-1
+                            tend(k,j,i) = tend(k,j,i) + &
+                                          particles(n)%weight_factor / &
+                                          prt_count(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz( tend, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nz_do3d
+                         local_pf(i,j,k) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+                resorted = .TRUE.
+             ELSE
+                CALL exchange_horiz( ql_vp_av, nbgp )
                 to_be_resorted => ql_vp_av
              ENDIF
Index: /palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_mask.f90	(revision 1006)
+++ /palm/trunk/SOURCE/data_output_mask.f90	(revision 1007)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! Bugfix: calculation of pr must depend on the particle weighting factor,
+! missing calculation of ql_vp added
 !
 ! Former revisions:
@@ -159,6 +161,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
@@ -213,5 +217,5 @@
                 to_be_resorted => q_av
              ENDIF
-             
+
           CASE ( 'ql' )
              IF ( av == 0 )  THEN
@@ -237,6 +241,28 @@
           CASE ( 'ql_vp' )
              IF ( av == 0 )  THEN
-                to_be_resorted => ql_vp
-             ELSE
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nz_do3d
+                         psi = prt_start_index(k,j,i)
+                         DO  n = psi, psi+prt_count(k,j,i)-1
+                            tend(k,j,i) = tend(k,j,i) + &
+                                          particles(n)%weight_factor / &
+                                          prt_count(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz( tend, nbgp )
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                   mask_j(mid,j),mask_i(mid,i))
+                      ENDDO
+                   ENDDO
+                ENDDO
+                resorted = .TRUE.
+             ELSE
+                CALL exchange_horiz( ql_vp_av, nbgp )
                 to_be_resorted => ql_vp_av
              ENDIF
@@ -264,5 +290,5 @@
                 to_be_resorted => rho_av
              ENDIF
-             
+
           CASE ( 's' )
              IF ( av == 0 )  THEN
@@ -271,5 +297,5 @@
                 to_be_resorted => s_av
              ENDIF
-             
+
           CASE ( 'sa' )
              IF ( av == 0 )  THEN
@@ -278,5 +304,5 @@
                 to_be_resorted => sa_av
              ENDIF
-             
+
           CASE ( 'u' )
              IF ( av == 0 )  THEN
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 1006)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 1007)
@@ -4,4 +4,8 @@
 ! Current revisions:
 ! -----------------
+! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
+! turbulent fluxes of WS-scheme
+! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
+! droplets at nzb and nzt added
 !
 ! Former revisions:
@@ -151,5 +155,5 @@
                         'timestep'
        CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
-      
+
     ENDIF
 
@@ -176,14 +180,14 @@
        IF ( ws_scheme_mom .AND. sr == 0 )  THEN
 
-!       
+!
 !--       According to the Neumann bc for the horizontal velocity components,
 !--       the corresponding fluxes has to satisfiy the same bc.
           IF ( ocean )  THEN
              sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
-             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)    
+             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
           ENDIF
 
           DO  i = 0, threads_per_task-1
-!          
+!
 !--          Swap the turbulent quantities evaluated in advec_ws.
              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u* 
@@ -453,5 +457,5 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
-                
+
 !
 !--             Higher moments
@@ -578,4 +582,5 @@
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
+
                    IF ( cloud_physics ) THEN
                       sums_l(k,51,tn) = sums_l(k,51,tn)                        &
@@ -618,10 +623,15 @@
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
+                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
+                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
+                                                  qsws(j,i) )
+                   IF ( cloud_droplets )  THEN
+                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
+                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
+                                           ql(nzb,j,i) ) * shf(j,i) +  &
+                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
+                   ENDIF
                    IF ( cloud_physics )  THEN
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
-                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
-                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
-                                                     qsws(j,i)            &
-                                                              )
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
@@ -657,10 +667,15 @@
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
+                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
+                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
+                                                  qswst(j,i) )
+                   IF ( cloud_droplets )  THEN
+                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
+                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
+                                            ql(nzt,j,i) ) * tswst(j,i) +  &
+                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
+                   ENDIF
                    IF ( cloud_physics )  THEN
-                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
-                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
-                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
-                                                     qsws(j,i)            &
-                                                              )
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
@@ -714,20 +729,24 @@
 !--             content
                 IF ( humidity )  THEN
-                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
-                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
-                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
-                                                       rmask(j,i,sr)
-
-                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-                      pts = 0.5 *                                            &
-                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
-                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
-                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
+                   IF ( cloud_physics .OR. cloud_droplets )  THEN
+                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
+                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                           rmask(j,i,sr)
                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
                                                           rmask(j,i,sr)
-                   ENDIF
-                ENDIF
-
+                   ELSE
+                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
+                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
+                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
+                                                             rmask(j,i,sr)
+                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
+                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
+                                             sums_l(k,17,tn) +      &
+                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
+                      END IF
+                   END IF
+                ENDIF
 !
 !--             Passive scalar flux
@@ -762,5 +781,5 @@
                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
                               v(k+1,j,i) - hom(k+1,1,2,sr) )
-!                              
+!
 !--               Momentum flux w*u*
                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
@@ -1002,4 +1021,5 @@
           ENDDO
        ENDIF
+
 !
 !--    Collect horizontal average in hom.
Index: /palm/trunk/SOURCE/lpm_collision_kernels.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1006)
+++ /palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1007)
@@ -4,5 +4,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! converted all units to SI units and replaced some parameters by corresponding
+! PALM parameters
+! Bugfix: factor in calculation of enhancement factor for collision efficencies
+! changed from 10. to 1.0
 !
 ! Former revisions:
@@ -64,5 +67,5 @@
 
     PUBLIC  ckernel, collision_efficiency_rogers, init_kernels, &
-            rclass_lbound, rclass_ubound, recalculate_kernel            
+            rclass_lbound, rclass_ubound, recalculate_kernel
 
     REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
@@ -126,12 +129,9 @@
 !             ENDIF
           ENDDO
-!
-!--       Collision routines expect radius to be in cm
-          radclass = radclass * 100.0
-
-!
-!--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
+
+!
+!--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
           DO  i = 1, dissipation_classes
-             epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
+             epsclass(i) = 0.1 * REAL( i ) / dissipation_classes
 !             IF ( myid == 0 )  THEN
 !                PRINT*, 'i=', i, ' eps = ', epsclass(i)
@@ -148,5 +148,5 @@
 
              epsilon = epsclass(k)
-             urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
+             urms    = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 )
 
              CALL turbsd
@@ -183,7 +183,9 @@
 
              PRINT*, '*** Hall kernel'
-             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
+                                              i = 1,radius_classes )
              DO  j = 1, radius_classes
-                WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
+                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
+                                          ( hkernel(i,j), i = 1,radius_classes )
              ENDDO
 
@@ -200,8 +202,11 @@
 
                 PRINT*, '*** epsilon = ', epsclass(k)
-                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
+                                                 i = 1,radius_classes )
                 DO  j = 1, radius_classes
-!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
-                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
+!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E6, &
+!                                       ( ckernel(i,j,k), i = 1,radius_classes )
+                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j)*1.0E6, &
+                                          ( hwratio(i,j), i = 1,radius_classes )
                 ENDDO
              ENDDO
@@ -210,6 +215,4 @@
 
           DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
-
-          ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
 
        ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
@@ -249,16 +252,15 @@
 
 !
-!--    Store particle radii on the radclass array. Collision routines
-!--    expect radii to be in cm.
-       radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
+!--    Store particle radii on the radclass array
+       radclass(1:radius_classes) = particles(pstart:pend)%radius
 
        IF ( wang_kernel )  THEN
-          epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
+          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
        ELSE
           epsilon = 0.0
        ENDIF
-       urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
-
-       IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
+       urms    = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 )
+
+       IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
 !
 !--       Call routines to calculate efficiencies for the Wang kernel
@@ -294,6 +296,4 @@
        ENDIF
 
-       ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
-
        DEALLOCATE( ec, radclass, winf )
 
@@ -314,4 +314,5 @@
        USE particle_attributes
        USE arrays_3d
+       USE control_parameters
 
        IMPLICIT NONE
@@ -326,6 +327,4 @@
                 v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
 
-       REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
-
        REAL, DIMENSION(1:radius_classes) ::  st, tau
 
@@ -336,27 +335,22 @@
 
           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   
-
-       lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
-       lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
+
+       ENDIF
+
+       lambda    = urms * SQRT( 15.0 * molecular_viscosity / epsilon )    ! in m
+       lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity )
        tl        = urms**2 / epsilon                       ! in s
-       lf        = 0.5 * urms**3 / epsilon                 ! in cm
-       tauk      = SQRT( anu / epsilon )                   ! in s
-       eta       = ( anu**3 / epsilon )**0.25              ! in cm
-       vk        = eta / tauk                              ! in cm/s
+       lf        = 0.5 * urms**3 / epsilon                 ! in m
+       tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
+       eta       = ( molecular_viscosity**3 / epsilon )**0.25             ! in m
+       vk        = eta / tauk
 
        ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
        tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
 
-       CALL fallg    ! gives winf in cm/s
+       CALL fallg    ! gives winf in m/s
 
        DO  i = 1, radius_classes
-          tau(i) = winf(i) / gravity    ! in s
+          tau(i) = winf(i) / g    ! in s
           st(i)  = tau(i) / tauk
        ENDDO
@@ -368,7 +362,7 @@
        bbb = SQRT( 1.0 - 2.0 * be**2 )
        d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
-       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
+       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in m
        d2  = ( 1.0 - bbb ) * 0.5 / bbb
-       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
+       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in m
        ccc = SQRT( 1.0 - 2.0 * z**2 )
        b1  = ( 1.0 + ccc ) * 0.5 / ccc
@@ -379,21 +373,21 @@
        DO  i = 1, radius_classes
 
-          v1 = winf(i)        ! in cm/s
+          v1 = winf(i)        ! in m/s
           t1 = tau(i)         ! in s
 
           DO  j = 1, i
-             rrp = radclass(i) + radclass(j)               ! radius in cm
-             v2  = winf(j)                                 ! in cm/s
+             rrp = radclass(i) + radclass(j)
+             v2  = winf(j)                                 ! in m/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
+             v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
+                     - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
+             v1xysq  = v1xysq * urms**2 / t1                ! in m**2/s**2
+             vrms1xy = SQRT( v1xysq )                       ! in m/s
+
+             v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
+                     - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
+             v2xysq  = v2xysq * urms**2 / t2                ! in m**2/s**2
+             vrms2xy = SQRT( v2xysq )                       ! in m/s
 
              IF ( winf(i) >= winf(j) )  THEN
@@ -414,9 +408,9 @@
                          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.0 * v1v2xy   ! in cm**2/s**2
+             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in m**2/s**2
+             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in m**2/s**2
              IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
              wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
-             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
+             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in m/s
 
 !
@@ -433,7 +427,7 @@
              yy = 0.1886 * EXP( 20.306 / lambda_re )
 
-             c1_gr  =  xx / ( gravity / vk * tauk )**yy
-
-             ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
+             c1_gr  =  xx / ( g / vk * tauk )**yy
+
+             ao_gr  = ao + ( pi / 8.0) * ( g / vk * tauk )**2
              fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
              rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
@@ -452,7 +446,7 @@
 
 !------------------------------------------------------------------------------!
-! phi as a function
-!------------------------------------------------------------------------------!
-    REAL FUNCTION phi( a, b, vsett, tau0 )
+! phi_w as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION phi_w( a, b, vsett, tau0 )
 
        IMPLICIT NONE
@@ -461,11 +455,11 @@
 
        aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
-       phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
-
-    END FUNCTION phi
-
-
-!------------------------------------------------------------------------------!
-! zeta as a function
+       phi_w = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
+
+    END FUNCTION phi_w
+
+
+!------------------------------------------------------------------------------!
+! zhi as a function
 !------------------------------------------------------------------------------!
     REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
@@ -490,5 +484,6 @@
 
 !------------------------------------------------------------------------------!
-! Calculation of terminal velocity winf
+! Calculation of terminal velocity winf following Equations 10-138 to 10-145
+! from (Pruppacher and Klett, 1997)
 !------------------------------------------------------------------------------!
     SUBROUTINE fallg
@@ -498,4 +493,5 @@
        USE particle_attributes
        USE arrays_3d
+       USE control_parameters
 
        IMPLICIT NONE
@@ -505,5 +501,5 @@
        LOGICAL, SAVE ::  first = .TRUE.
 
-       REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
+       REAL, SAVE ::  cunh, eta, phy, py, rho_a, sigma, stb, stok, &
                       t0, xlamb
 
@@ -523,16 +519,14 @@
                  -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 ) )
+!
+!--       Parameter values for p = 1013,25 hPa and T = 293,15 K
+          eta   = 1.818E-5         ! in kg/(m s)
+          xlamb = 6.6E-8           ! in m
+          rho_a = 1.204            ! in kg/m**3
+          cunh  = 1.26 * xlamb     ! in m
+          sigma = 0.07363          ! in kg/s**2
+          stok  = 2.0  * g * ( rho_l - rho_a ) / ( 9.0 * eta ) ! in 1/(m s)
+          stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
+          phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
           py    = phy**( 1.0 / 6.0 )
 
@@ -541,9 +535,9 @@
        DO  j = 1, radius_classes
 
-          IF ( radclass(j) <= 1.0E-3 ) THEN
-
-             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
-
-          ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
+          IF ( radclass(j) <= 1.0E-5 ) THEN
+
+             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
+
+          ELSEIF ( radclass(j) > 1.0E-5  .AND.  radclass(j) <= 5.35E-4 )  THEN
 
              x = LOG( stb * radclass(j)**3 )
@@ -553,14 +547,17 @@
                 y = y + b(i) * x**(i-1)
              ENDDO
-
-             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
-             winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
-
-          ELSEIF ( radclass(j) > 5.35E-2 )  THEN
-
-             IF ( radclass(j) > 0.35 )  THEN
-                bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
+!
+!--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
+!--          for correct version see (Beard, 1976)
+             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y ) 
+
+             winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
+
+          ELSEIF ( radclass(j) > 5.35E-4 )  THEN
+
+             IF ( radclass(j) > 0.0035 )  THEN
+                bond = g * ( rho_l - rho_a ) * 0.0035**2 / sigma
              ELSE
-                bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
+               bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
              ENDIF
 
@@ -574,8 +571,8 @@
              xrey = py * EXP( y )
 
-             IF ( radclass(j) > 0.35 )  THEN
-                winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
+             IF ( radclass(j) > 0.0035 )  THEN
+                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 )
              ELSE
-                winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
+                winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
              ENDIF
 
@@ -668,7 +665,8 @@
 !
 !--    Calculate the radius class index of particles with respect to array r
+!--    Radius has to be in Âµm
        ALLOCATE( ira(1:radius_classes) )
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E4   ! in microm
+          particle_radius = radclass(j) * 1.0E6
           DO  k = 1, 15
              IF ( particle_radius < r0(k) )  THEN
@@ -693,5 +691,5 @@
              IF ( ir < 16 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
+                   pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / &
                         ( r0(ir) - r0(ir-1) )
                    qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
@@ -754,5 +752,5 @@
           rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
 !
-!--       In 100 cm^2/s^3
+!--       for 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 /)
@@ -767,5 +765,5 @@
           ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
 !
-!--       In 400 cm^2/s^3
+!--       for 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 /)
@@ -784,8 +782,9 @@
 !
 !--    Calculate the radius class index of particles with respect to array r0
+!--    Radius has to be in Âµm
        ALLOCATE( ira(1:radius_classes) )
 
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E4    ! in microm
+          particle_radius = radclass(j) * 1.0E6
           DO  k = 1, 7
              IF ( particle_radius < r0(k) )  THEN
@@ -799,4 +798,5 @@
 !
 !--    Two-dimensional linear interpolation of the collision efficiencies
+!--    Radius has to be in Âµm
        DO  j =  1, radius_classes
           DO  i = 1, j
@@ -812,14 +812,13 @@
              ENDDO
 
-             y1 = 1.0      ! 0  cm2/s3
-!
-!--          100 cm2/s3, 400 cm2/s3
+             y1 = 0.0001      ! for 0 m**2/s**3
+
              IF ( ir < 8 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
+                   pp = ( radclass(j)*1.0E6 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
                    qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
                    y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
                                 pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
-                                qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
+                                qq * ( 1.0-pp ) * ecoll_100(ir-1,iq)   +  &
                                 pp * qq         * ecoll_100(ir,iq)
                    y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
@@ -838,14 +837,14 @@
              ENDIF
 !
-!--          Linear interpolation of dissipation rate in cm2/s3
-             IF ( epsilon <= 100.0 )  THEN
-                ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
-                         + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
-             ELSEIF ( epsilon <= 600.0 )  THEN
-                ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
-                         + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
+!--          Linear interpolation of dissipation rate in m**2/s**3
+             IF ( epsilon <= 0.01 )  THEN
+                ecf(j,i) = ( epsilon - 0.01 ) / (   0.0 - 0.01 ) * y1 &
+                         + ( epsilon -   0.0 ) / ( 0.01 -   0.0 ) * y2
+             ELSEIF ( epsilon <= 0.06 )  THEN
+                ecf(j,i) = ( epsilon - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
+                         + ( epsilon - 0.01 ) / ( 0.04 - 0.01 ) * y3
              ELSE
-                ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
-                         + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
+                ecf(j,i) = (   0.06 - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
+                         + (   0.06 - 0.01 ) / ( 0.04 - 0.01 ) * y3
              ENDIF
 
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 1006)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 1007)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix: calculation of buoyancy production has to consider the liquid water
+! mixing ratio in case of cloud droplets
 !
 ! Former revisions:
@@ -164,5 +165,5 @@
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-                
+
              ENDDO
           ENDDO
@@ -189,5 +190,5 @@
 
                    IF ( wall_e_y(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                   Inconsistency removed: as the thermal stratification is 
 !--                   not taken into account for the evaluation of the wall 
@@ -219,5 +220,5 @@
 
                    IF ( wall_e_x(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                   Inconsistency removed: as the thermal stratification is 
 !--                   not taken into account for the evaluation of the wall 
@@ -275,5 +276,5 @@
 
                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                      Inconsistency removed: as the thermal stratification 
 !--                      is not taken into account for the evaluation of the 
@@ -301,5 +302,5 @@
 
                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                      Inconsistency removed: as the thermal stratification 
 !--                      is not taken into account for the evaluation of the 
@@ -416,5 +417,5 @@
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-                
+
                 ENDIF
 
@@ -551,8 +552,13 @@
                    DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
-                      IF ( .NOT. cloud_physics )  THEN
+                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
-                      ELSE
+                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
+                                         g / vpt(k,j,i) *                      &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+                                         ) * dd2zu(k)
+                      ELSE IF ( cloud_physics )  THEN
                          IF ( ql(k,j,i) == 0.0 )  THEN
                             k1 = 1.0 + 0.61 * q(k,j,i)
@@ -568,10 +574,20 @@
                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                          ENDIF
-                      ENDIF
-
-                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
+                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
+                                         g / vpt(k,j,i) *                      &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
                                          ) * dd2zu(k)
+                      ELSE IF ( cloud_droplets )  THEN
+                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61 * pt(k,j,i)
+                         tend(k,j,i) = tend(k,j,i) -                          & 
+                                       kh(k,j,i) * g / vpt(k,j,i) *           &
+                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +  &
+                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -  &
+                                         pt(k,j,i) * ( ql(k+1,j,i) -          &
+                                         ql(k-1,j,i) ) ) * dd2zu(k)
+                      ENDIF
+
                    ENDDO
 
@@ -584,8 +600,8 @@
                       k = nzb_diff_s_inner(j,i)-1
 
-                      IF ( .NOT. cloud_physics )  THEN
+                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
-                      ELSE
+                      ELSE IF ( cloud_physics )  THEN
                          IF ( ql(k,j,i) == 0.0 )  THEN
                             k1 = 1.0 + 0.61 * q(k,j,i)
@@ -601,4 +617,7 @@
                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                          ENDIF
+                      ELSE IF ( cloud_droplets )  THEN
+                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61 * pt(k,j,i)
                       ENDIF
 
@@ -615,8 +634,8 @@
                       k = nzt
 
-                      IF ( .NOT. cloud_physics )  THEN
+                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
-                      ELSE
+                      ELSE IF ( cloud_physics )  THEN
                          IF ( ql(k,j,i) == 0.0 )  THEN
                             k1 = 1.0 + 0.61 * q(k,j,i)
@@ -632,4 +651,7 @@
                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                          ENDIF
+                      ELSE IF ( cloud_droplets )  THEN
+                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61 * pt(k,j,i)
                       ENDIF
 
@@ -699,5 +721,5 @@
 
           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-                
+
        ENDDO
 
@@ -721,5 +743,5 @@
 
              IF ( wall_e_y(j,i) /= 0.0 )  THEN
-!                     
+!
 !--             Inconsistency removed: as the thermal stratification 
 !--             is not taken into account for the evaluation of the 
@@ -751,5 +773,5 @@
 
              IF ( wall_e_x(j,i) /= 0.0 )  THEN
-!                     
+!
 !--             Inconsistency removed: as the thermal stratification 
 !--             is not taken into account for the evaluation of the 
@@ -805,5 +827,5 @@
 
                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                Inconsistency removed: as the thermal stratification 
 !--                is not taken into account for the evaluation of the 
@@ -831,5 +853,5 @@
 
                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
-!                     
+!
 !--                Inconsistency removed: as the thermal stratification 
 !--                is not taken into account for the evaluation of the 
@@ -1041,8 +1063,12 @@
              DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
-                IF ( .NOT. cloud_physics )  THEN
+                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
-                ELSE
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
+                                         ) * dd2zu(k)
+                ELSE IF ( cloud_physics )  THEN
                    IF ( ql(k,j,i) == 0.0 )  THEN
                       k1 = 1.0 + 0.61 * q(k,j,i)
@@ -1058,10 +1084,17 @@
                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                    ENDIF
-                ENDIF
-
-                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
                                            k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
                                          ) * dd2zu(k)
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61 * pt(k,j,i)
+                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
+                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
+                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -    &
+                                       pt(k,j,i) * ( ql(k+1,j,i) -           &
+                                                     ql(k-1,j,i) ) ) * dd2zu(k)
+                ENDIF
              ENDDO
 
@@ -1069,8 +1102,8 @@
                 k = nzb_diff_s_inner(j,i)-1
 
-                IF ( .NOT. cloud_physics ) THEN
+                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
-                ELSE
+                ELSE IF ( cloud_physics )  THEN
                    IF ( ql(k,j,i) == 0.0 )  THEN
                       k1 = 1.0 + 0.61 * q(k,j,i)
@@ -1086,4 +1119,7 @@
                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                    ENDIF
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61 * pt(k,j,i)
                 ENDIF
 
@@ -1095,8 +1131,8 @@
                 k = nzt
 
-                IF ( .NOT. cloud_physics ) THEN
+                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
-                ELSE
+                ELSE IF ( cloud_physics )  THEN
                    IF ( ql(k,j,i) == 0.0 )  THEN
                       k1 = 1.0 + 0.61 * q(k,j,i)
@@ -1112,4 +1148,7 @@
                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
                    ENDIF
+                ELSE IF ( cloud_droplets )  THEN
+                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61 * pt(k,j,i)
                 ENDIF
 
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1006)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1007)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix in calculation of ql_vp
 !
 ! Former revisions:
@@ -388,8 +388,13 @@
 
           CASE ( 'ql_vp' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + ql_vp(k,j,i)
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nzt+1
+                      psi = prt_start_index(k,j,i)
+                      DO  n = psi, psi+prt_count(k,j,i)-1
+                         ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + &
+                                           particles(n)%weight_factor / &
+                                           prt_count(k,j,i)
+                      ENDDO
                    ENDDO
                 ENDDO
