Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 415)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 420)
@@ -1,10 +1,10 @@
 New:
 ---
-humidity=.T. is now usable for runs with topography. wall_humidityflux and 
+humidity=.T. is now usable for runs with topography. wall_humidityflux and
 wall_scalarflux are the corresponding new parin arrays.
 (check_parameters, init_3d_model, parin)
 
-Large scale vertical motion (subsidence/ascent) can be applied to the 
-prognostic equation for the potential temperature. (check_parameters, header, 
+Large scale vertical motion (subsidence/ascent) can be applied to the
+prognostic equation for the potential temperature. (check_parameters, header,
 Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
 write_var_list)
@@ -18,5 +18,12 @@
 var_ts is replaced by dots_max (modules,init_3d_model)
 
-init_3d_model, modules
+Every cloud droplet has now an own weighting factor and can be deleted due to
+collisions. Condensation and collision of cloud droplets are adjusted
+accordingly. (advec_particles)
+
+Collision efficiency for large cloud droplets has changed according to table of
+Rogers and Yau. (collision_efficiency)
+
+advec_particles, collision_efficiency, init_3d_model, modules
 
 
@@ -25,4 +32,7 @@
 Dimension of array stat in cascade change to prevent type problems with
 mpi2 libraries (poisfft_hybrid)
+
+A loop has been splitted to make runs reproducible on HLRN
+systems. (disturb_field)
 
 Bugfix: exchange of ghost points for prho included (time_integration)
@@ -35,5 +45,11 @@
 messages from gfortran compiler (modules)
 
-calc_precipitation, modules, poisfft_hybrid, sum_up_3d_data, time_integration
+Bugfix: calculation of cloud droplet velocity (advec_particles)
+
+Bugfix: transfer of particles at south/left edge (advec_particles)
+
+Bugfix: calculation of collision_efficiency (collision_efficiency)
+
+advec_particles, calc_precipitation, collision_efficiency, disturb_field, modules, poisfft_hybrid, sum_up_3d_data, time_integration
 
 
Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 415)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 420)
@@ -5,4 +5,15 @@
 ! -----------------
 ! TEST: PRINT statements on unit 9 (commented out)
+! Own weighting factor for every cloud droplet is implemented and condensation
+! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3,
+! -s_r4
+! Initialization of variables for the (sub-) timestep is moved to the beginning
+! in order to enable deletion of cloud droplets due to collision processes.
+! Collision efficiency for large cloud droplets has changed according to table
+! of Rogers and Yau. (collision_efficiency)
+! Bugfix: calculation of cloud droplet velocity
+! Bugfix: transfer of particles at south/left edge when new position
+!         y>=(ny+0.5)*dy-1.e-12 or x>=(nx+0.5)*dx-1.e-12, very rare
+! Bugfix: calculation of y (collision_efficiency)
 !
 ! Former revisions:
@@ -117,15 +128,15 @@
     LOGICAL ::  dt_3d_reached, dt_3d_reached_l, prt_position
 
-    REAL    ::  aa, arg, bb, cc, dd, delta_r, dens_ratio, de_dt, de_dt_min,    &
-                de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, de_dy_int_l,   &
-                de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, diss_int,    &
-                diss_int_l, diss_int_u, distance, dt_gap, dt_particle,         &
-                dt_particle_m, d_radius, d_sum, e_a, e_int, e_int_l, e_int_u,  &
-                e_mean_int, e_s, exp_arg, exp_term, fs_int, gg,                &
-                lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,        &
+    REAL    ::  aa, arg, bb, cc, dd, delta_r, delta_v, dens_ratio, de_dt,      &
+                de_dt_min, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int,     &
+                de_dy_int_l, de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, &
+                diss_int, diss_int_l, diss_int_u, distance, dt_gap,            &
+                dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int,       &
+                e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int,  &
+                gg,lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,     &
                 pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, &
-                random_gauss, sl_r3, sl_r4, s_r3, s_r4, t_int, u_int, u_int_l, &
-                u_int_u, vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l,      &
-                w_int_u, x, y
+                random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,    &
+                vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u,      &
+                x, y
 
     REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
@@ -180,4 +191,22 @@
        CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
     ENDIF
+
+!
+!--    Initialize variables for the (sub-) timestep, i.e. for
+!--    marking those particles to be deleted after the timestep
+       particle_mask     = .TRUE.
+       deleted_particles = 0
+       trlp_count_recv   = 0
+       trnp_count_recv   = 0
+       trrp_count_recv   = 0
+       trsp_count_recv   = 0
+       trlpt_count_recv  = 0
+       trnpt_count_recv  = 0
+       trrpt_count_recv  = 0
+       trspt_count_recv  = 0
+       IF ( use_particle_tails )  THEN
+          tail_mask     = .TRUE.
+       ENDIF
+       deleted_tails = 0
 
 !
@@ -215,4 +244,5 @@
 !--    Reset summation arrays
        ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
+       delta_r=0.0;  delta_v=0.0
 
 !
@@ -343,5 +373,5 @@
 !--       Sum up the total volume of liquid water (needed below for
 !--       re-calculating the weighting factors)
-          ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
+          ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor *             &
                                       particles(n)%radius**3
        ENDDO
@@ -351,5 +381,5 @@
 !--    Particle growth by collision
        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
-       
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -374,5 +404,5 @@
                          tmp_particle = particles(psi+is)
                          js = is
-                         DO WHILE ( particles(psi+js-inc)%radius > &
+                         DO WHILE ( particles(psi+js-inc)%radius >             &
                                     tmp_particle%radius )
                             particles(psi+js) = particles(psi+js-inc)
@@ -386,31 +416,26 @@
 !
 !--                Calculate the mean radius of all those particles which
-!--                are of smaller or equal size than the current particle
+!--                are of smaller size than the current particle
 !--                and use this radius for calculating the collision efficiency
                    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
-!
-!--                   There may be some particles of size equal to the
-!--                   current particle but with larger index
+
+                   DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
                       sl_r3 = 0.0
                       sl_r4 = 0.0
-		      DO is = n, psi+prt_count(k,j,i)-2
-			 IF ( particles(is+1)%radius == &
-			      particles(is)%radius ) THEN
-                            sl_r3 = sl_r3 + particles(is+1)%radius**3
-                            sl_r4 = sl_r4 + particles(is+1)%radius**4
-		         ELSE
-			    EXIT
-			 ENDIF
-		      ENDDO
-
-                      IF ( ( s_r3 + sl_r3 ) > 0.0 )  THEN
-
-                         mean_r = ( s_r4 + sl_r4 ) / ( s_r3 + sl_r3 )
-
-                         CALL collision_efficiency( mean_r,              &
-                                                    particles(n)%radius, &
+
+                      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 )
 
@@ -418,12 +443,7 @@
                          effective_coll_efficiency = 0.0
                       ENDIF
-       
-!
-!--                   Contribution of the current particle to the next one
-                      s_r3 = s_r3 + particles(n)%radius**3
-                      s_r4 = s_r4 + particles(n)%radius**4
-
-                      IF ( effective_coll_efficiency > 1.0  .OR. &
-                           effective_coll_efficiency < 0.0 )     &
+
+                      IF ( effective_coll_efficiency > 1.0  .OR.               &
+                           effective_coll_efficiency < 0.0 )                   &
                       THEN
                          WRITE( message_string, * )  'collision_efficiency ' , &
@@ -431,5 +451,5 @@
                          CALL message( 'advec_particles', 'PA0145', 2, 2, -1,  &
                                        6, 1 )
-                      END IF
+                      ENDIF
 
 !
@@ -459,5 +479,5 @@
                                  ) / ( 3.0 * gg )
 
-                      ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * &
+                      ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *   &
                                           ( ql_int_u - ql_int_l )
 
@@ -478,8 +498,8 @@
                       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) &
+                      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
@@ -491,5 +511,5 @@
                                                              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_int_l + ( particles(n)%z - zu(kk) ) / dz *  &
                                            ( u_int_u - u_int_l )
                       ENDIF
@@ -558,35 +578,87 @@
 
 !
-!--                   Change in radius due to collision
-                      delta_r = effective_coll_efficiency * &
-                                ql_int * rho_surface / ( 1.0 - ql_int ) *   &
-                                0.25 / rho_l *                              &
-                                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 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)%radius**3
-
+                      ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
+                                      *( particles(n)%radius**3 )
                    ENDDO
 
-                ENDIF
-
-!
-!--             Re-calculate the weighting factor (total liquid water content
-!--             must be conserved during collision)
-                IF ( ql_vp(k,j,i) /= 0.0 )  THEN
-
-                   ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
-!
-!--                Re-assign this weighting factor to the particles of the
-!--                current gridbox
+                   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)
-                   DO  n = psi, psi + prt_count(k,j,i)-1
-                      particles(n)%weight_factor = ql_vp(k,j,i)
-                   ENDDO
-
+                   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
 
@@ -726,5 +798,5 @@
 !
 !--       Final values are obtained by division by the total number of grid
-!--       points used for the summation. 
+!--       points used for the summation.
           hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
           hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
@@ -761,5 +833,5 @@
           CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
-                  
+
 #else
           sums(:,8)  = sums_l(:,8,0)
@@ -779,5 +851,5 @@
        ENDIF
 
-    ENDIF 
+    ENDIF
 
 !
@@ -815,23 +887,4 @@
        dt_3d_reached_l = .TRUE.
 
-!
-!--    Initialize variables for the (sub-) timestep, i.e. for marking those
-!--    particles to be deleted after the timestep
-       particle_mask     = .TRUE.
-       deleted_particles = 0
-       trlp_count_recv   = 0
-       trnp_count_recv   = 0
-       trrp_count_recv   = 0
-       trsp_count_recv   = 0
-       trlpt_count_recv  = 0
-       trnpt_count_recv  = 0
-       trrpt_count_recv  = 0
-       trspt_count_recv  = 0
-       IF ( use_particle_tails )  THEN
-          tail_mask     = .TRUE.
-       ENDIF
-       deleted_tails = 0
-
-
        DO  n = 1, number_of_particles
 !
@@ -847,5 +900,5 @@
           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
               + offset_ocean_nzt                     ! only exact if equidistant
-                                               
+
 !
 !--       Interpolation of the velocity components in the xy-plane
@@ -942,5 +995,5 @@
                  + offset_ocean_nzt                      ! only exact if eq.dist
 
-             IF ( topography == 'flat' )  THEN        
+             IF ( topography == 'flat' )  THEN 
 
                 x  = particles(n)%x - i * dx
@@ -1066,5 +1119,5 @@
 !--             gp_outside_of_building(7): i+1,j,k+1
 !--             gp_outside_of_building(8): i+1,j+1,k+1
-           
+
                 gp_outside_of_building = 0
                 location = 0.0
@@ -1092,5 +1145,5 @@
                    location(num_gp,3) = k * dz - 0.5 * dz
                    ei(num_gp)     = e(k,j+1,i)
-                   dissi(num_gp)  = diss(k,j+1,i)        
+                   dissi(num_gp)  = diss(k,j+1,i)
                    de_dxi(num_gp) = de_dx(k,j+1,i)
                    de_dyi(num_gp) = de_dy(k,j+1,i)
@@ -1105,5 +1158,5 @@
                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
                    ei(num_gp)     = e(k+1,j,i)
-                   dissi(num_gp)  = diss(k+1,j,i)        
+                   dissi(num_gp)  = diss(k+1,j,i)
                    de_dxi(num_gp) = de_dx(k+1,j,i)
                    de_dyi(num_gp) = de_dy(k+1,j,i)
@@ -1119,5 +1172,5 @@
                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
                    ei(num_gp)     = e(k+1,j+1,i)
-                   dissi(num_gp)  = diss(k+1,j+1,i)        
+                   dissi(num_gp)  = diss(k+1,j+1,i)
                    de_dxi(num_gp) = de_dx(k+1,j+1,i)
                    de_dyi(num_gp) = de_dy(k+1,j+1,i)
@@ -1133,5 +1186,5 @@
                    location(num_gp,3) = k * dz - 0.5 * dz
                    ei(num_gp)     = e(k,j,i+1)
-                   dissi(num_gp)  = diss(k,j,i+1) 
+                   dissi(num_gp)  = diss(k,j,i+1)
                    de_dxi(num_gp) = de_dx(k,j,i+1)
                    de_dyi(num_gp) = de_dy(k,j,i+1)
@@ -1147,5 +1200,5 @@
                    location(num_gp,3) = k * dz - 0.5 * dz
                    ei(num_gp)     = e(k,j+1,i+1)
-                   dissi(num_gp)  = diss(k,j+1,i+1) 
+                   dissi(num_gp)  = diss(k,j+1,i+1)
                    de_dxi(num_gp) = de_dx(k,j+1,i+1)
                    de_dyi(num_gp) = de_dy(k,j+1,i+1)
@@ -1161,5 +1214,5 @@
                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
                    ei(num_gp)     = e(k+1,j,i+1)
-                   dissi(num_gp)  = diss(k+1,j,i+1) 
+                   dissi(num_gp)  = diss(k+1,j,i+1)
                    de_dxi(num_gp) = de_dx(k+1,j,i+1)
                    de_dyi(num_gp) = de_dy(k+1,j,i+1)
@@ -1175,5 +1228,5 @@
                    location(num_gp,3) = (k+1) * dz - 0.5 * dz
                    ei(num_gp)     = e(k+1,j+1,i+1)
-                   dissi(num_gp)  = diss(k+1,j+1,i+1) 
+                   dissi(num_gp)  = diss(k+1,j+1,i+1)
                    de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
                    de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
@@ -1182,5 +1235,5 @@
 
 !
-!--             If all gridpoints are situated outside of a building, then the 
+!--             If all gridpoints are situated outside of a building, then the
 !--             ordinary interpolation scheme can be used.
                 IF ( num_gp == 8 )  THEN
@@ -1309,7 +1362,7 @@
                       de_dxi(num_gp) = 0.0
                       de_dyi(num_gp) = de_dy(k,j,i)
-                      de_dzi(num_gp) = de_dz(k,j,i)                     
-                   ENDIF 
-   
+                      de_dzi(num_gp) = de_dz(k,j,i)
+                   ENDIF
+
                    IF ( gp_outside_of_building(5) == 1  .AND. &
                       gp_outside_of_building(1) == 0 )  THEN
@@ -1367,7 +1420,7 @@
                       de_dxi(num_gp) = 0.0
                       de_dyi(num_gp) = de_dy(k,j+1,i)
-                      de_dzi(num_gp) = de_dz(k,j+1,i)                     
-                   ENDIF 
-   
+                      de_dzi(num_gp) = de_dz(k,j+1,i)
+                   ENDIF
+
                    IF ( gp_outside_of_building(6) == 1  .AND. &
                         gp_outside_of_building(2) == 0 )  THEN
@@ -1425,7 +1478,7 @@
                       de_dxi(num_gp) = 0.0
                       de_dyi(num_gp) = de_dy(k+1,j,i)
-                      de_dzi(num_gp) = de_dz(k+1,j,i)                     
-                   ENDIF 
-   
+                      de_dzi(num_gp) = de_dz(k+1,j,i) 
+                   ENDIF
+
                    IF ( gp_outside_of_building(7) == 1  .AND. &
                         gp_outside_of_building(3) == 0 )  THEN
@@ -1483,7 +1536,7 @@
                       de_dxi(num_gp) = 0.0
                       de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                      de_dzi(num_gp) = de_dz(k+1,j+1,i)                     
-                   ENDIF 
-   
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                   ENDIF
+
                    IF ( gp_outside_of_building(8) == 1  .AND. &
                         gp_outside_of_building(4) == 0 )  THEN
@@ -1600,6 +1653,6 @@
                    IF ( num_gp == 1 )  THEN
 !
-!--                   If only one of the gridpoints is situated outside of the 
-!--                   building, it follows that the values at the particle 
+!--                   If only one of the gridpoints is situated outside of the
+!--                   building, it follows that the values at the particle
 !--                   location are the same as the gridpoint values
                       e_int     = ei(num_gp)
@@ -1612,6 +1665,6 @@
                       d_sum = 0.0
 !
-!--                   Evaluation of the distances between the gridpoints 
-!--                   contributing to the interpolated values, and the particle 
+!--                   Evaluation of the distances between the gridpoints
+!--                   contributing to the interpolated values, and the particle
 !--                   location
                       DO agp = 1, num_gp
@@ -1839,5 +1892,5 @@
              IF ( cloud_droplets )  THEN
                 exp_arg  = 4.5 * dens_ratio * molecular_viscosity /        &
-                           ( particles(n)%radius )**2 /                    &
+                           ( particles(n)%radius )**2 *                    &
                            ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius *    &
                              SQRT( ( u_int - particles(n)%speed_x )**2 +   &
@@ -2178,4 +2231,10 @@
                       deleted_particles = deleted_particles + 1
 
+                      IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN
+                         trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10
+                         trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x &
+                                                     - 1
+                      ENDIF
+
                       IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
                          trlpt_count = trlpt_count + 1
@@ -2617,4 +2676,10 @@
                          deleted_particles = deleted_particles + 1
 
+                         IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN
+                            trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10
+                            trsp(trsp_count)%origin_y =                        &
+                                                    trsp(trsp_count)%origin_y - 1
+                         ENDIF
+
                          IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
                             trspt_count = trspt_count + 1
@@ -3296,6 +3361,23 @@
        IF ( dt_3d_reached )  EXIT
 
+!
+!--    Initialize variables for the next (sub-) timestep, i.e. for marking those
+!--    particles to be deleted after the timestep
+       particle_mask     = .TRUE.
+       deleted_particles = 0
+       trlp_count_recv   = 0
+       trnp_count_recv   = 0
+       trrp_count_recv   = 0
+       trsp_count_recv   = 0
+       trlpt_count_recv  = 0
+       trnpt_count_recv  = 0
+       trrpt_count_recv  = 0
+       trspt_count_recv  = 0
+       IF ( use_particle_tails )  THEN
+          tail_mask     = .TRUE.
+       ENDIF
+       deleted_tails = 0
+
     ENDDO   ! timestep loop
-
 
 !
@@ -3308,11 +3390,4 @@
     ENDIF
 
-
-!
-!-- Re-evaluate the weighting factors. After advection, particles within a
-!-- grid box may have different weighting factors if some have been advected
-!-- from a neighbouring box. The factors are re-evaluated so that they are
-!-- the same for all particles of one box. This procedure must conserve the
-!-- liquid water content within one box.
     IF ( cloud_droplets )  THEN
 
@@ -3322,5 +3397,5 @@
 
 !
-!--    Re-calculate the weighting factors and calculate the liquid water content
+!--    Calculate the liquid water content
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -3328,32 +3403,28 @@
 
 !
-!--             Calculate the total volume of particles in the boxes (ql_vp) as
-!--             well as the real volume (ql_v, weighting factor has to be
-!--             included)
+!--             Calculate the total volume in the boxes (ql_v, weighting factor
+!--             has to beincluded)
                 psi = prt_start_index(k,j,i)
                 DO  n = psi, psi+prt_count(k,j,i)-1
-                   ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
-
-                   ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor * &
+                   ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor *  &
                                                  particles(n)%radius**3
                 ENDDO
 
 !
-!--             Re-calculate the weighting factors and calculate the liquid
-!--             water content
-                IF ( ql_vp(k,j,i) /= 0.0 )  THEN
-                   ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
-                   ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * &
-                                           ql_v(k,j,i) /             &
+!--             Calculate the liquid water content
+                IF ( ql_v(k,j,i) /= 0.0 )  THEN
+                   ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi *           &
+                                           ql_v(k,j,i) /                       &
                                            ( rho_surface * dx * dy * dz )
+
+                   IF ( ql(k,j,i) < 0.0 ) THEN
+                      WRITE( message_string, * )  'LWC out of range: ' , &
+                                                  ql(k,j,i)
+                         CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
+                   ENDIF
+
                 ELSE
                    ql(k,j,i) = 0.0
                 ENDIF
-
-!
-!--             Re-assign the weighting factor to the particles
-                DO  n = psi, psi+prt_count(k,j,i)-1
-                   particles(n)%weight_factor = ql_vp(k,j,i)
-                ENDDO
 
              ENDDO
@@ -3833,15 +3904,23 @@
        e = 0.001
     ELSEIF ( mean_rm >= 25.0 )  THEN
-       IF( j <= 3 )  e = 0.55
+       IF( j <= 2 )  e = 0.0
+       IF( j == 3 )  e = 0.47
        IF( j == 4 )  e = 0.8
        IF( j == 5 )  e = 0.9
        IF( j >=6  )  e = 1.0
     ELSEIF ( rm >= 3000.0 )  THEN
-       e = 1.0
+       IF( i == 1 )  e = 0.02
+       IF( i == 2 )  e = 0.16
+       IF( i == 3 )  e = 0.33
+       IF( i == 4 )  e = 0.55
+       IF( i == 5 )  e = 0.71
+       IF( i == 6 )  e = 0.81
+       IF( i == 7 )  e = 0.90
+       IF( i >= 8 )  e = 0.94
     ELSE
        x  = mean_rm - collected_r(i)
-       y  = rm - collected_r(j)
-       dx = collected_r(i+1) - collected_r(i) 
-       dy = collector_r(j+1) - collector_r(j) 
+       y  = rm - collector_r(j)
+       dx = collected_r(i+1) - collected_r(i)
+       dy = collector_r(j+1) - collector_r(j)
        aa = x**2 + y**2
        bb = ( dx - x )**2 + y**2
@@ -3854,5 +3933,5 @@
     ENDIF
 
- END SUBROUTINE collision_efficiency 
+ END SUBROUTINE collision_efficiency
 
 
@@ -3902,5 +3981,5 @@
           WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', &
                           j, ' k=', k,                                       &
-                          ' &nxl=', nxl, ' nxr=', nxr,                       &
+                          ' nxl=', nxl, ' nxr=', nxr,                        &
                           ' nys=', nys, ' nyn=', nyn,                        &
                           ' nzb=', nzb, ' nzt=', nzt
Index: palm/trunk/SOURCE/disturb_field.f90
===================================================================
--- palm/trunk/SOURCE/disturb_field.f90	(revision 415)
+++ palm/trunk/SOURCE/disturb_field.f90	(revision 420)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! A loop has been splitted to make runs reproducible on HLRN systems
 !
 ! Former revisions:
@@ -98,4 +98,5 @@
 !
 !-- Exchange of ghost points for the random perturbation
+
     CALL exchange_horiz( dist1 )
 
@@ -104,14 +105,20 @@
 !-- Neighboured grid points in all three directions are used for the
 !-- filter operation.
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
+!-- Loop has been splitted to make runs reproducible on HLRN systems using
+!-- compiler option -O3
+     DO  i = nxl, nxr
+        DO  j = nys, nyn
           DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
-             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) + dist1(k,j-1,i) &
-                            + dist1(k,j+1,i) + dist1(k+1,j,i) + dist1(k-1,j,i) &
+             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) &
+                            + dist1(k,j+1,i) + dist1(k+1,j,i) &
+                            ) / 12.0
+          ENDDO
+          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
+              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
                             + 6.0 * dist1(k,j,i)                               &
                             ) / 12.0
           ENDDO
-       ENDDO
-    ENDDO
+        ENDDO
+     ENDDO
 
 !
@@ -119,4 +126,5 @@
 !-- Afterwards, filter operation and exchange of ghost points are repeated.
     CALL exchange_horiz( dist2 )
+
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -129,4 +137,5 @@
        ENDDO
     ENDDO
+
     CALL exchange_horiz( dist1 )
 
