Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 2262)
+++ /palm/trunk/SOURCE/Makefile	(revision 2263)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Implemented splitting and merging algorithm
+# 
+# 2259 2017-06-08 09:09:11Z gronemeier
 # Implemented synthetic turbulence generator
 #
@@ -354,7 +357,7 @@
 	lpm_data_output_particles.f90 lpm_droplet_collision.f90 \
 	lpm_droplet_condensation.f90 lpm_exchange_horiz.f90 \
-	lpm_init.f90 lpm_init_sgs_tke.f90 \
+	lpm_init.f90 lpm_init_sgs_tke.f90 lpm_merging.f90 \
 	lpm_pack_arrays.f90 lpm_read_restart_file.f90 \
-	lpm_set_attributes.f90 \
+	lpm_set_attributes.f90 lpm_splitting.f90 \
 	lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
 	ls_forcing_mod.f90 message.f90 microphysics_mod.f90 modules.f90 mod_kinds.f90 \
@@ -509,6 +512,6 @@
 local_tremain.o: modules.o cpulog_mod.o mod_kinds.o
 local_tremain_ini.o: modules.o cpulog_mod.o mod_kinds.o
-lpm.o: modules.o cpulog_mod.o lpm_exchange_horiz.o lpm_init.o lpm_pack_arrays.o \
-	mod_kinds.o mod_particle_attributes.o
+lpm.o: modules.o cpulog_mod.o lpm_exchange_horiz.o lpm_init.o lpm_merging.o lpm_pack_arrays.o \
+	lpm_splitting.o mod_kinds.o mod_particle_attributes.o
 lpm_advec.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o surface_mod.o
 lpm_boundary_conds.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
@@ -530,8 +533,11 @@
         surface_mod.o
 lpm_init_sgs_tke.o: modules.o mod_kinds.o mod_particle_attributes.o surface_mod.o
+lpm_merging.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
 lpm_pack_arrays.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
 lpm_read_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o \
 	lpm_pack_arrays.o
 lpm_set_attributes.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
+lpm_splitting.o: modules.o cpulog_mod.o mod_kinds.o \
+	mod_particle_attributes.o lpm_exchange_horiz.o
 lpm_write_exchange_statistics.o: modules.o mod_kinds.o mod_particle_attributes.o
 lpm_write_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o
Index: /palm/trunk/SOURCE/lpm.f90
===================================================================
--- /palm/trunk/SOURCE/lpm.f90	(revision 2262)
+++ /palm/trunk/SOURCE/lpm.f90	(revision 2263)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implemented splitting and merging algorithm
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -142,12 +145,12 @@
         ONLY:  collision_kernel, deleted_particles, deallocate_memory,         &
                dt_write_particle_data, dt_prel, end_time_prel,                 &
-               grid_particles, number_of_particles, number_of_particle_groups, &
-               particles, particle_groups, prt_count, step_dealloc,            &
-               time_prel, time_write_particle_data, trlp_count_sum,            &
-               trlp_count_recv_sum, trnp_count_sum,                            &
-               trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
-               trsp_count_sum, trsp_count_recv_sum,                            &
+               grid_particles, merging,  number_of_particles,                  &
+               number_of_particle_groups, particles, particle_groups,          &
+               prt_count, splitting, step_dealloc, time_prel,                  &
+               time_write_particle_data, trlp_count_sum, trlp_count_recv_sum,  &
+               trnp_count_sum, trnp_count_recv_sum, trrp_count_sum,            & 
+               trrp_count_recv_sum, trsp_count_sum, trsp_count_recv_sum,       &
                use_sgs_for_particles, write_particle_statistics
-
+                              
     USE pegrid
 
@@ -369,4 +372,14 @@
        IF ( dt_3d_reached )  time_prel = time_prel + dt_3d
 !
+!--    Apply splitting and merging algorithm
+       IF ( cloud_droplets )  THEN
+          IF ( splitting ) THEN
+             CALL lpm_splitting
+          ENDIF
+          IF ( merging ) THEN
+             CALL lpm_merging
+          ENDIF
+       ENDIF
+!
 !--    Move Particles local to PE to a different grid cell
        CALL lpm_move_particle
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 2262)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 2263)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implemented splitting and merging algorithm
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -175,5 +178,5 @@
                 block_offset, block_offset_def, collision_kernel,              &
                 curvature_solution_effects,                                    &
-                density_ratio, grid_particles,                                 &
+                density_ratio, grid_particles, isf,i_splitting_mode,           &
                 initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns,   &
                 ibc_par_t, iran_part, log_z_z0,                                &
@@ -190,5 +193,5 @@
                 radius, random_start_position, read_particles_from_restartfile,&
                 seed_follows_topography, sgs_wf_part, sort_count,              &
-                total_number_of_particles,                                     &
+                splitting_function, splitting_mode, total_number_of_particles,  &
                 use_sgs_for_particles,                                         &
                 write_particle_statistics, uniform_particles, zero_particle,   &
@@ -473,4 +476,39 @@
           
     END SELECT
+    SELECT CASE ( splitting_mode )
+    
+       CASE ( 'const' )
+          i_splitting_mode = 1
+
+       CASE ( 'cl_av' )
+          i_splitting_mode = 2
+
+       CASE ( 'gb_av' )
+          i_splitting_mode = 3
+          
+       CASE DEFAULT
+          WRITE( message_string, * )  'unknown splitting condition ',   &
+                                       'splitting_mode = "', TRIM( splitting_mode ), '"'
+          CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 )
+          
+    END SELECT
+    SELECT CASE ( splitting_function )
+    
+       CASE ( 'gamma' )
+          isf = 1
+
+       CASE ( 'log' )
+          isf = 2
+
+       CASE ( 'exp' )
+          isf = 3
+          
+       CASE DEFAULT
+          WRITE( message_string, * )  'unknown splitting function ',   &
+                                       'splitting_function = "', TRIM( splitting_function ), '"'
+          CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 )
+          
+    END SELECT
+    
 
 !
Index: /palm/trunk/SOURCE/lpm_merging.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_merging.f90	(revision 2263)
+++ /palm/trunk/SOURCE/lpm_merging.f90	(revision 2263)
@@ -0,0 +1,123 @@
+!> @file lpm_merging.f90
+!------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+!
+! Inital revision
+! 
+! 
+!
+! Description:
+! ------------
+! This routine is a part of the Lagrangian particle model. Two Super droplets 
+! which fulfill certain criterion's (e.g. a big weighting factor and a small
+! radius) can be merged into one super droplet with a increased number of 
+! represented particles of the super droplet. This mechanism ensures an
+! improved a feasible amount of computational costs. The limits of particle 
+! creation should be chosen carefully! The idea of this algorithm is based on 
+! Unterstrasser and Soelch, 2014. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_merging
+
+
+    USE arrays_3d,                                                             &
+        ONLY:  ql
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
+
+    USE kinds
+
+    USE particle_attributes,                                                   &
+        ONLY:  deleted_particles, grid_particles, initial_weighting_factor,    &     
+               merge_drp, merging, number_of_particles, particles, prt_count,  &
+               radius_merge, sum_merge_drp, weight_factor_merge
+
+    USE pegrid
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i         !<      
+    INTEGER(iwp) ::  j         !<       
+    INTEGER(iwp) ::  k         !<       
+    INTEGER(iwp) ::  n         !<   
+    
+    REAL(wp) ::  ql_crit = 1.0E-5_wp  !< threshold lwc for cloudy grid cells 
+                                      !< (e.g. Siebesma et al 2003, JAS, 60)
+    
+    CALL cpu_log( log_point_s(81), 'lpm_merging', 'start' )
+
+    merge_drp  = 0
+    
+    IF ( weight_factor_merge == -1.0_wp )  THEN
+       weight_factor_merge = 0.5_wp * initial_weighting_factor 
+    ENDIF
+
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          DO  k = nzb+1, nzt
+   
+             number_of_particles = prt_count(k,j,i)
+             IF ( number_of_particles <= 0  .OR.                               &
+                   ql(k,j,i) >= ql_crit )  CYCLE
+             particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+!
+!--          Start merging operations: This routine delete super droplets with
+!--          a small radius (radius <= radius_merge) and a low weighting 
+!--          factor (weight_factor  <= weight_factor_merge). The number of 
+!--          represented particles will be added to the next particle of the 
+!--          particle array. Tests showed that this simplified method can be 
+!--          used because it will only take place outside of cloudy grid 
+!--          boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned   
+!--          and subsequent evaporated super droplets will be merged.
+             DO  n = 1, number_of_particles-1
+                IF ( particles(n)%particle_mask                    .AND.       &
+                     particles(n+1)%particle_mask                  .AND.       &
+                     particles(n)%radius        <= radius_merge    .AND.       &
+                     particles(n)%weight_factor <= weight_factor_merge )       &   
+                THEN
+                   particles(n+1)%weight_factor  =                             &
+                                       particles(n+1)%weight_factor +          &
+                                       ( particles(n)%radius**3     /          &
+                                         particles(n+1)%radius**3   *          &
+                                         particles(n)%weight_factor            &
+                                       )
+                   particles(n)%particle_mask = .FALSE.
+                   deleted_particles          = deleted_particles + 1 
+                   merge_drp                  = merge_drp + 1
+                
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
+       
+    sum_merge_drp = sum_merge_drp + merge_drp
+
+    CALL cpu_log( log_point_s(81), 'lpm_merging', 'stop' )
+
+ END SUBROUTINE lpm_merging
Index: /palm/trunk/SOURCE/lpm_splitting.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_splitting.f90	(revision 2263)
+++ /palm/trunk/SOURCE/lpm_splitting.f90	(revision 2263)
@@ -0,0 +1,595 @@
+!> @file lpm_splitting.f90
+!------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+!
+! Initial revision
+! 
+! 
+!
+! Description:
+! ------------
+! This routine is a part of the Lagrangian particle model. Super droplets which 
+! fulfill certain criterion's (e.g. a big weighting factor and a large radius) 
+! can be split into several super droplets with a reduced number of 
+! represented particles of every super droplet. This mechanism ensures an
+! improved representation of the right tail of the drop size distribution with 
+! a feasible amount of computational costs. The limits of particle creation 
+! should be chosen carefully! The idea of this algorithm is based on 
+! Unterstrasser and Soelch, 2014. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_splitting
+
+
+    USE arrays_3d,                                                             &
+        ONLY:  ql
+
+    USE cloud_parameters,                                                      &
+        ONLY:  rho_l
+
+    USE constants,                                                             &
+        ONLY:  pi
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
+
+    USE kinds
+
+    USE lpm_exchange_horiz_mod,                                                &
+        ONLY:  realloc_particles_array
+
+    USE particle_attributes,                                                   &
+        ONLY:  grid_particles, iran_part, initial_weighting_factor, isf,       &
+               i_splitting_mode, max_number_particles_per_gridbox,             &  
+               new_particles, n_max, number_concentration,                     &
+               number_of_particles, number_particles_per_gridbox, particles,   &
+               particle_type, prt_count, radius_split, splitting,              &
+               splitting_factor, splitting_factor_max, splitting_mode,         &
+               sum_new_particles, weight_factor_split                        
+
+    USE pegrid
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i                !< 
+    INTEGER(iwp) ::  j                !<
+    INTEGER(iwp) ::  jpp              !<
+    INTEGER(iwp) ::  k                !<
+    INTEGER(iwp) ::  n                !<
+    INTEGER(iwp) ::  new_particles_gb !< counter of created particles within one grid box
+    INTEGER(iwp) ::  new_size         !< new particle array size
+    INTEGER(iwp) ::  np               !<
+    INTEGER(iwp) ::  old_size         !< old particle array size
+    
+    LOGICAL ::  first_loop_stride = .TRUE. !< flag to calculate constants only once
+
+    REAL(wp) ::  diameter                 !< diameter of droplet
+    REAL(wp) ::  dlog                     !< factor for DSD calculation
+    REAL(wp) ::  factor_volume_to_mass    !< pre calculate factor volume to mass
+    REAL(wp) ::  lambda                   !< slope parameter of gamma-distribution
+    REAL(wp) ::  lwc                      !< liquid water content of grid box
+    REAL(wp) ::  lwc_total                !< average liquid water content of cloud
+    REAL(wp) ::  m1                       !< first moment of DSD
+    REAL(wp) ::  m1_total                 !< average over all PEs of first moment of DSD
+    REAL(wp) ::  m2                       !< second moment of DSD
+    REAL(wp) ::  m2_total                 !< average average over all PEs second moment of DSD
+    REAL(wp) ::  m3                       !< third moment of DSD
+    REAL(wp) ::  m3_total                 !< average average over all PEs third moment of DSD
+    REAL(wp) ::  mu                       !< spectral shape parameter of gamma distribution
+    REAL(wp) ::  nrclgb                   !< number of cloudy grid boxes (ql >= 1.0E-5 kg/kg) 
+    REAL(wp) ::  nrclgb_total             !< average over all PEs of number of cloudy grid boxes
+    REAL(wp) ::  nr                       !< number concentration of cloud droplets
+    REAL(wp) ::  nr_total                 !< average over all PEs of number of cloudy grid boxes
+    REAL(wp) ::  nr0                      !< intercept parameter of gamma distribution
+    REAL(wp) ::  pirho_l                  !< pi * rho_l / 6.0
+    REAL(wp) ::  ql_crit = 1.0E-5_wp      !< threshold lwc for cloudy grid cells 
+                                          !< (Siebesma et al 2003, JAS, 60)
+    REAL(wp) ::  rm                       !< volume averaged mean radius
+    REAL(wp) ::  rm_total                 !< average over all PEs of volume averaged mean radius
+    REAL(wp) ::  r_min = 1.0E-6_wp        !< minimum radius of approximated spectra 
+    REAL(wp) ::  r_max = 1.0E-3_wp        !< maximum radius of approximated spectra
+    REAL(wp) ::  sigma_log = 1.5_wp       !< standard deviation of the LOG-distribution
+    REAL(wp) ::  zeta                     !< Parameter for DSD calculation of Seifert
+
+    REAL(wp), DIMENSION(0:n_max-1) ::  an_spl     !< size dependent critical weight factor
+    REAL(wp), DIMENSION(0:n_max-1) ::  r_bin_mid  !< mass weighted mean radius of a bin
+    REAL(wp), DIMENSION(0:n_max)   ::  r_bin      !< boundaries of a radius bin
+    
+    TYPE(particle_type) ::  tmp_particle   !< temporary particle TYPE
+
+    CALL cpu_log( log_point_s(80), 'lpm_splitting', 'start' )
+
+    IF ( first_loop_stride )  THEN
+       IF ( i_splitting_mode == 2  .OR.  i_splitting_mode == 3 )  THEN
+          dlog   = ( LOG10(r_max) - LOG10(r_min) ) / ( n_max - 1 )
+          DO  i = 0, n_max-1
+             r_bin(i) = 10.0_wp**( LOG10(r_min) + i * dlog - 0.5_wp * dlog )
+             r_bin_mid(i) = 10.0_wp**( LOG10(r_min) + i * dlog )
+          ENDDO
+          r_bin(n_max) = 10.0_wp**( LOG10(r_min) + n_max * dlog - 0.5_wp * dlog )
+       ENDIF   
+       factor_volume_to_mass =  4.0_wp / 3.0_wp * pi * rho_l
+       pirho_l  = pi * rho_l / 6.0_wp
+       IF ( weight_factor_split == -1.0_wp )  THEN
+          weight_factor_split = 0.1_wp * initial_weighting_factor 
+       ENDIF
+    ENDIF
+
+    new_particles  = 0
+
+    IF ( i_splitting_mode == 1 )  THEN
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+
+                new_particles_gb = 0
+                number_of_particles = prt_count(k,j,i)
+                IF ( number_of_particles <= 0  .OR.                            &  
+                     ql(k,j,i) < ql_crit )  CYCLE
+                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+!                
+!--             Start splitting operations. Each particle is checked if it
+!--             fulfilled the splitting criterion's. In splitting mode 'const'   
+!--             a critical radius  (radius_split) a critical weighting factor
+!--             (weight_factor_split) and a splitting factor (splitting_factor)
+!--             must  be prescribed (see particles_par). Super droplets which 
+!--             have a larger radius and larger weighting factor are split into
+!--             'splitting_factor' super droplets. Therefore, the weighting 
+!--             factor of  the super droplet and all created clones is reduced 
+!--             by the factor of 'splitting_factor'.
+                DO  n = 1, number_of_particles
+                   IF ( particles(n)%particle_mask  .AND.                      &
+                        particles(n)%radius >= radius_split  .AND.             & 
+                        particles(n)%weight_factor >= weight_factor_split )    &
+                   THEN          
+!
+!--                   Calculate the new number of particles.
+                      new_size = prt_count(k,j,i) + splitting_factor - 1                                    
+!
+!--                   Cycle if maximum number of particles per grid box
+!--                   is greater than the allowed maximum number.
+                      IF ( new_size >= max_number_particles_per_gridbox )  CYCLE                      
+!
+!--                   Reallocate particle array if necessary. 
+                      IF ( new_size > SIZE(particles) )  THEN 
+                         CALL realloc_particles_array(i,j,k,new_size)
+                      ENDIF
+                      old_size = prt_count(k,j,i)
+!
+!--                   Calculate new weighting factor.
+                      particles(n)%weight_factor =  & 
+                         particles(n)%weight_factor / splitting_factor
+                      tmp_particle = particles(n)
+!
+!--                   Create splitting_factor-1 new particles.
+                      DO  jpp = 1, splitting_factor-1
+                         grid_particles(k,j,i)%particles(jpp+old_size) =       & 
+                            tmp_particle         
+                      ENDDO  
+                      new_particles_gb = new_particles_gb + splitting_factor - 1
+!   
+!--                   Save the new number of super droplets for every grid box.
+                      prt_count(k,j,i) = prt_count(k,j,i) +                    &
+                                         splitting_factor - 1         
+                   ENDIF
+                ENDDO
+                
+                new_particles       = new_particles     + new_particles_gb
+                sum_new_particles   = sum_new_particles + new_particles_gb 
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSEIF ( i_splitting_mode == 2 )  THEN 
+!
+!--    Initialize summing variables.
+       lwc          = 0.0_wp
+       lwc_total    = 0.0_wp 
+       m1           = 0.0_wp
+       m1_total     = 0.0_wp
+       m2           = 0.0_wp
+       m2_total     = 0.0_wp
+       m3           = 0.0_wp
+       m3_total     = 0.0_wp
+       nr           = 0.0_wp   
+       nrclgb       = 0.0_wp
+       nrclgb_total = 0.0_wp 
+       nr_total     = 0.0_wp
+       rm           = 0.0_wp
+       rm_total     = 0.0_wp
+       
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                number_of_particles = prt_count(k,j,i)
+                IF ( number_of_particles <= 0  .OR.                            & 
+                     ql(k,j,i) < ql_crit )  CYCLE
+                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+                nrclgb = nrclgb + 1.0_wp
+!                
+!--             Calculate moments of DSD.                
+                DO  n = 1, number_of_particles
+                   IF ( particles(n)%particle_mask  .AND.                      &
+                        particles(n)%radius >= r_min )                         &
+                   THEN
+                      nr  = nr  + particles(n)%weight_factor
+                      rm  = rm  + factor_volume_to_mass  *                     &
+                                 particles(n)%radius**3  *                     &
+                                 particles(n)%weight_factor
+                      IF ( isf == 1 )  THEN           
+                         diameter   = particles(n)%radius * 2.0_wp            
+                         lwc = lwc + factor_volume_to_mass *                   &
+                                     particles(n)%radius**3 *                  & 
+                                     particles(n)%weight_factor 
+                         m1  = m1  + particles(n)%weight_factor * diameter                                                
+                         m2  = m2  + particles(n)%weight_factor * diameter**2           
+                         m3  = m3  + particles(n)%weight_factor * diameter**3
+                      ENDIF   
+                   ENDIF                        
+                ENDDO 
+             ENDDO
+          ENDDO
+       ENDDO
+
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( nr, nr_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       CALL MPI_ALLREDUCE( rm, rm_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( nrclgb, nrclgb_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( lwc, lwc_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( m1, m1_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( m2, m2_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( m3, m3_total, 1 , &
+       MPI_REAL, MPI_SUM, comm2d, ierr )
+#endif 
+
+!
+!--    Calculate number concentration and mean volume averaged radius.
+       nr_total = MERGE( nr_total / nrclgb_total,                              &
+                         0.0_wp, nrclgb_total > 0.0_wp                         &
+                       )
+       rm_total = MERGE( ( rm_total /                                          &
+                            ( nr_total * factor_volume_to_mass )               &
+                          )**0.3333333_wp, 0.0_wp, nrclgb_total > 0.0_wp       &
+                       )                         
+!
+!--    Check which function should be used to approximate the DSD.
+       IF ( isf == 1 )  THEN
+          lwc_total = MERGE( lwc_total / nrclgb_total,                         &
+                             0.0_wp, nrclgb_total > 0.0_wp                     &
+                           )
+          m1_total  = MERGE( m1_total / nrclgb_total,                          &
+                             0.0_wp, nrclgb_total > 0.0_wp                     &
+                           )
+          m2_total  = MERGE( m2_total / nrclgb_total,                          &
+                             0.0_wp, nrclgb_total > 0.0_wp                     &
+                           )
+          m3_total  = MERGE( m3_total / nrclgb_total,                          &
+                             0.0_wp, nrclgb_total > 0.0_wp                     &
+                           )
+          zeta = m1_total * m3_total / m2_total**2                             
+          mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /                &
+                        ( zeta - 1.0_wp ), 0.0_wp                              &
+                    )
+
+          lambda = ( pirho_l * nr_total / lwc_total *                          &
+                     ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &                                          
+                   )**0.3333333_wp
+          nr0 = nr_total / gamma( mu + 1.0_wp ) * lambda**( mu + 1.0_wp ) 
+          
+          DO  n = 0, n_max-1
+             diameter  = r_bin_mid(n) * 2.0_wp           
+             an_spl(n) = nr0 * diameter**mu * EXP( -lambda * diameter ) *      & 
+                         ( r_bin(n+1) - r_bin(n) ) * 2.0_wp 
+          ENDDO
+       ELSEIF ( isf == 2 )  THEN
+          DO  n = 0, n_max-1
+             an_spl(n) = nr_total / ( SQRT( 2.0_wp * pi ) *                    &
+                                     LOG(sigma_log) * r_bin_mid(n)             &
+                                     ) *                                       &
+                         EXP( -( LOG( r_bin_mid(n) / rm_total )**2 ) /         &
+                               ( 2.0_wp * LOG(sigma_log)**2 )                  & 
+                             ) *                                               & 
+                         ( r_bin(n+1) - r_bin(n) )
+          ENDDO
+       ELSEIF( isf == 3 )  THEN
+          DO  n = 0, n_max-1 
+             an_spl(n) = 3.0_wp * nr_total * r_bin_mid(n)**2 / rm_total**3  *  &
+                         EXP( - ( r_bin_mid(n)**3 / rm_total**3 ) )         *  &
+                         ( r_bin(n+1) - r_bin(n) )
+          ENDDO
+       ENDIF
+!
+!--    Criterion to avoid super droplets with a weighting factor < 1.0.
+       an_spl = MAX(an_spl, 1.0_wp)
+  
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                number_of_particles = prt_count(k,j,i)
+                IF ( number_of_particles <= 0  .OR.                            &  
+                     ql(k,j,i) < ql_crit )  CYCLE
+                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+                new_particles_gb = 0         
+!                
+!--             Start splitting operations. Each particle is checked if it
+!--             fulfilled the splitting criterion's. In splitting mode 'cl_av'   
+!--             a critical radius (radius_split) and a splitting function must 
+!--             be prescribed (see particles_par). The critical weighting factor
+!--             is calculated while approximating a 'gamma', 'log' or 'exp'-
+!--             drop size distribution. In this mode the DSD is calculated as 
+!--             an average over all cloudy grid boxes. Super droplets which 
+!--             have a larger radius and larger weighting factor are split into
+!--             'splitting_factor' super droplets. In this case the splitting  
+!--             factor is calculated of weighting factor of the super droplet  
+!--             and the approximated number concentration for droplet of such 
+!--             a size. Due to the splitting, the weighting factor of the  
+!--             super droplet and all created clones is reduced by the factor  
+!--             of 'splitting_facor'.
+                DO  n = 1, number_of_particles
+                   DO  np = 0, n_max-1
+                      IF ( r_bin(np) >= radius_split  .AND.                    &
+                           particles(n)%particle_mask  .AND.                   &
+                           particles(n)%radius >= r_bin(np)  .AND.             &
+                           particles(n)%radius < r_bin(np+1)  .AND.            &
+                           particles(n)%weight_factor >= an_spl(np)  )         &
+                      THEN
+!
+!--                      Calculate splitting factor
+                         splitting_factor =                                    & 
+                             MIN( INT( particles(n)%weight_factor /            &
+                                        an_spl(np)                             &
+                                     ), splitting_factor_max                   &
+                                )
+                         IF ( splitting_factor < 2 )  CYCLE
+!
+!--                      Calculate the new number of particles.                                                           
+                         new_size = prt_count(k,j,i) + splitting_factor - 1
+!
+!--                      Cycle if maximum number of particles per grid box
+!--                      is greater than the allowed maximum number.                         
+                         IF ( new_size >= max_number_particles_per_gridbox )   & 
+                         CYCLE
+!
+!--                      Reallocate particle array if necessary. 
+                         IF ( new_size > SIZE(particles) )  THEN 
+                            CALL realloc_particles_array(i,j,k,new_size)
+                         ENDIF
+                         old_size  = prt_count(k,j,i)                             
+                         new_particles_gb = new_particles_gb +                 &
+                                            splitting_factor - 1
+!
+!--                      Calculate new weighting factor.
+                         particles(n)%weight_factor =                          & 
+                            particles(n)%weight_factor / splitting_factor
+                         tmp_particle = particles(n)
+!
+!--                      Create splitting_factor-1 new particles.                                                  
+                         DO  jpp = 1, splitting_factor-1
+                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
+                                                                    tmp_particle         
+                         ENDDO
+!   
+!--                      Save the new number of super droplets. 
+                         prt_count(k,j,i) = prt_count(k,j,i) +                 &
+                                            splitting_factor - 1
+                      ENDIF
+                   ENDDO
+                ENDDO 
+                
+                new_particles       = new_particles     + new_particles_gb
+                sum_new_particles   = sum_new_particles + new_particles_gb                    
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSEIF ( i_splitting_mode == 3 )  THEN 
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+             
+!
+!--             Initialize summing variables.             
+                lwc = 0.0_wp
+                m1  = 0.0_wp
+                m2  = 0.0_wp
+                m3  = 0.0_wp
+                nr  = 0.0_wp
+                rm  = 0.0_wp  
+                
+                new_particles_gb = 0
+                number_of_particles = prt_count(k,j,i)
+                IF ( number_of_particles <= 0  .OR.                            & 
+                     ql(k,j,i) < ql_crit )  CYCLE
+                particles => grid_particles(k,j,i)%particles
+!                
+!--             Calculate moments of DSD.                
+                DO  n = 1, number_of_particles
+                   IF ( particles(n)%particle_mask  .AND.                      &
+                        particles(n)%radius >= r_min )                         &
+                   THEN
+                      nr  = nr + particles(n)%weight_factor
+                      rm  = rm + factor_volume_to_mass  *                      &
+                                 particles(n)%radius**3  *                     &
+                                 particles(n)%weight_factor
+                      IF ( isf == 1 )  THEN            
+                         diameter   = particles(n)%radius * 2.0_wp           
+                         lwc = lwc + factor_volume_to_mass *                   &
+                                     particles(n)%radius**3 *                  &  
+                                     particles(n)%weight_factor 
+                         m1  = m1 + particles(n)%weight_factor * diameter
+                         m2  = m2 + particles(n)%weight_factor * diameter**2
+                         m3  = m3 + particles(n)%weight_factor * diameter**3
+                      ENDIF      
+                   ENDIF                                            
+                ENDDO  
+
+                IF ( nr <= 0.0  .OR.  rm <= 0.0_wp )  CYCLE 
+!
+!--             Calculate mean volume averaged radius.                
+                rm = ( rm / ( nr * factor_volume_to_mass ) )**0.3333333_wp
+!
+!--             Check which function should be used to approximate the DSD.              
+                IF ( isf == 1 )  THEN
+!
+!--                Gamma size distribution to calculate  
+!--                critical weight_factor (e.g. Marshall + Palmer, 1948).
+                   zeta = m1 * m3 / m2**2
+                   mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /       &
+                                ( zeta - 1.0_wp ), 0.0_wp                      &
+                             )   
+                   lambda = ( pirho_l * nr / lwc *                             &
+                              ( mu + 3.0_wp ) * ( mu + 2.0_wp ) *              &
+                              ( mu + 1.0_wp )                                  &                                  
+                            )**0.3333333_wp
+                   nr0 =  ( nr / (gamma( mu + 1.0_wp ) ) ) *                   &
+                          lambda**( mu + 1.0_wp ) 
+
+                   DO  n = 0, n_max-1
+                      diameter         = r_bin_mid(n) * 2.0_wp           
+                      an_spl(n) = nr0 * diameter**mu *                         &
+                                  EXP( -lambda * diameter ) *                  & 
+                                  ( r_bin(n+1) - r_bin(n) ) * 2.0_wp 
+                   ENDDO
+                ELSEIF ( isf == 2 )  THEN
+!
+!--                Lognormal size distribution to calculate critical 
+!--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).   
+                   DO  n = 0, n_max-1
+                      an_spl(n) = nr / ( SQRT( 2.0_wp * pi ) *                 &
+                                              LOG(sigma_log) * r_bin_mid(n)    &
+                                        ) *                                    &
+                                  EXP( -( LOG( r_bin_mid(n) / rm )**2 ) /      &
+                                        ( 2.0_wp * LOG(sigma_log)**2 )         & 
+                                      ) *                                      & 
+                                  ( r_bin(n+1) - r_bin(n) )
+                   ENDDO
+                ELSEIF ( isf == 3 )  THEN
+!
+!--                Exponential size distribution to calculate critical 
+!--                weight_factor (e.g. Berry + Reinhardt, 1974).  
+                   DO  n = 0, n_max-1
+                      an_spl(n) = 3.0_wp * nr * r_bin_mid(n)**2 / rm**3 *     &
+                                  EXP( - ( r_bin_mid(n)**3 / rm**3 ) ) *      &
+                                  ( r_bin(n+1) - r_bin(n) )
+                   ENDDO
+                ENDIF
+                
+!
+!--             Criterion to avoid super droplets with a weighting factor < 1.0.                                   
+                an_spl = MAX(an_spl, 1.0_wp)
+!                
+!--             Start splitting operations. Each particle is checked if it
+!--             fulfilled the splitting criterion's. In splitting mode 'gb_av'   
+!--             a critical radius (radius_split) and a splitting function must 
+!--             be prescribed (see particles_par). The critical weighting factor
+!--             is calculated while appoximating a 'gamma', 'log' or 'exp'-
+!--             drop size distribution. In this mode a DSD is calculated for 
+!--             every cloudy grid box. Super droplets which have a larger 
+!--             radius and larger weighting factor are split into
+!--             'splitting_factor' super droplets. In this case the splitting  
+!--             factor is calculated of weighting factor of the super droplet  
+!--             and theapproximated number concentration for droplet of such 
+!--             a size. Due to the splitting, the weighting factor of the  
+!--             super droplet and all created clones is reduced by the factor  
+!--             of 'splitting_facor'.
+                DO  n = 1, number_of_particles
+                   DO  np = 0, n_max-1
+                      IF ( r_bin(np) >= radius_split  .AND.                    &
+                           particles(n)%particle_mask  .AND.                   &
+                           particles(n)%radius >= r_bin(np)    .AND.           &
+                           particles(n)%radius < r_bin(np+1)   .AND.           &
+                           particles(n)%weight_factor >= an_spl(np) )          &
+                      THEN
+!
+!--                      Calculate splitting factor.
+                         splitting_factor =                                    & 
+                             MIN( INT( particles(n)%weight_factor /            &
+                                        an_spl(np)                             &
+                                     ), splitting_factor_max                   &
+                                )
+                         IF ( splitting_factor < 2 )  CYCLE
+
+!
+!--                      Calculate the new number of particles.                                                                                         
+                         new_size = prt_count(k,j,i) + splitting_factor - 1
+!
+!--                      Cycle if maximum number of particles per grid box
+!--                      is greater than the allowed maximum number.                                                  
+                         IF ( new_size >= max_number_particles_per_gridbox )   &
+                         CYCLE
+!
+!--                      Reallocate particle array if necessary.                          
+                         IF ( new_size > SIZE(particles) )  THEN 
+                            CALL realloc_particles_array(i,j,k,new_size)
+                         ENDIF
+!
+!--                      Calculate new weighting factor.
+                         particles(n)%weight_factor = & 
+                            particles(n)%weight_factor / splitting_factor
+                         tmp_particle               = particles(n)
+                         old_size                   = prt_count(k,j,i)
+!
+!--                      Create splitting_factor-1 new particles.
+                         DO jpp = 1, splitting_factor-1
+                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
+                               tmp_particle                  
+                         ENDDO
+!
+!--                      Save the new number of droplets for every grid box.
+                         prt_count(k,j,i)    = prt_count(k,j,i) +              &
+                                               splitting_factor - 1
+                         new_particles_gb    = new_particles_gb +              &
+                                               splitting_factor - 1                                       
+                      ENDIF
+                   ENDDO  
+                ENDDO
+
+                new_particles       = new_particles + new_particles_gb
+                sum_new_particles   = sum_new_particles + new_particles_gb                                                 
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+        
+    CALL cpu_log( log_point_s(80), 'lpm_splitting', 'stop' )
+
+ END SUBROUTINE lpm_splitting
+ 
Index: /palm/trunk/SOURCE/mod_particle_attributes.f90
===================================================================
--- /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 2262)
+++ /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 2263)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implemented splitting and merging algorithm
+! 
+! 2183 2017-03-17 14:29:15Z schwenkel
 !
 ! 2182 2017-03-17 14:27:40Z schwenkel
@@ -86,60 +89,122 @@
     CHARACTER(LEN=15) ::  collision_algorithm = 'all_or_nothing'  !< collision algorithm
     CHARACTER(LEN=15) ::  collision_kernel = 'none'               !< collision kernel
-
-    INTEGER(iwp) ::  deleted_particles = 0,                                    &
-                     dissipation_classes = 10, ibc_par_lr,                     &
-                     ibc_par_ns, ibc_par_b, ibc_par_t, iran_part = -1234567,   &
-                     maximum_number_of_particles = 0,                          &
-                     min_nr_particle = 50,                                     &
-                     mpi_particle_type,                                        &
-                     number_of_particles = 0,                                  &
-                     number_of_particle_groups = 1,                            &
-                     number_of_sublayers = 20,                                 &
-                     number_particles_per_gridbox = -1,                        &
-                     offset_ocean_nzt = 0,                                     &
-                     offset_ocean_nzt_m1 = 0, particles_per_point = 1,         &
-                     particle_file_count = 0, radius_classes = 20,             &
-                     sort_count = 0, step_dealloc = 100,                       &
-                     total_number_of_particles,                                &
-                     trlp_count_sum, trlp_count_recv_sum, trrp_count_sum,      &
-                     trrp_count_recv_sum, trsp_count_sum, trsp_count_recv_sum, &
-                     trnp_count_sum, trnp_count_recv_sum
-
-    INTEGER(iwp), PARAMETER ::  max_number_of_particle_groups = 10
-
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  prt_count
-
-    LOGICAL ::  all_or_nothing = .FALSE., average_impact = .FALSE.,            &
-                curvature_solution_effects = .FALSE.,                          &
-                deallocate_memory = .TRUE.,                                    &
-                hall_kernel = .FALSE., particle_advection = .FALSE.,           &
-                random_start_position = .FALSE.,                               &
-                read_particles_from_restartfile = .TRUE.,                      &
-                seed_follows_topography = .FALSE.,                             &
-                uniform_particles = .TRUE., use_kernel_tables = .FALSE.,       &
-                use_sgs_for_particles = .FALSE., wang_kernel = .FALSE.,        &
-                write_particle_statistics = .FALSE.
-
+    CHARACTER(LEN=5) ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
+    CHARACTER(LEN=5) ::  splitting_mode = 'const'                 !< splitting mode
+
+
+    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step 
+    INTEGER(iwp) ::  dissipation_classes = 10                     !< number of dissipation classes
+    INTEGER(iwp) ::  ibc_par_b                                    !< particle bottom boundary condition dummy
+    INTEGER(iwp) ::  ibc_par_lr                                   !< particle left/right boundary condition dummy
+    INTEGER(iwp) ::  ibc_par_ns                                   !< particle north/south boundary condition dummy
+    INTEGER(iwp) ::  ibc_par_t                                    !< particle top boundary condition dummy
+    INTEGER(iwp) ::  iran_part = -1234567                         !< number for random generator
+    INTEGER(iwp) ::  isf                                          !< dummy for splitting function
+    INTEGER(iwp) ::  i_splitting_mode                             !< dummy for splitting mode
+    INTEGER(iwp) ::  maximum_number_of_particles = 0              !< maximum number of particles on a PE (can be removed?)
+    INTEGER(iwp) ::  max_number_particles_per_gridbox = 30        !< maximum number of particles per gridbox (used in splitting algorithm)
+    INTEGER(iwp) ::  merge_drp = 0                                !< number of merged droplets
+    INTEGER(iwp) ::  min_nr_particle = 50                         !< minimum number of particles for which memory is allocated at every grid cell         
+    INTEGER(iwp) ::  mpi_particle_type                            !< parameter for particle PE particle exchange
+    INTEGER(iwp) ::  new_particles = 0                            !< number of new particles
+    INTEGER(iwp) ::  n_max = 100                                  !< number of radii bin for splitting functions     
+    INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)           
+    INTEGER(iwp) ::  number_of_particle_groups = 1                !< number of particle groups to be used            
+    INTEGER(iwp) ::  number_of_sublayers = 20                     !< number of sublayers for particle velocities betwenn surface and first grid level         
+    INTEGER(iwp) ::  number_particles_per_gridbox = -1            !< number of particles which are created in every grid box            
+    INTEGER(iwp) ::  offset_ocean_nzt = 0                         !< in case of oceans runs, the vertical index calculations need an offset         
+    INTEGER(iwp) ::  offset_ocean_nzt_m1 = 0                      !< in case of oceans runs, the vertical index calculations need an offset
+    INTEGER(iwp) ::  particles_per_point = 1                      !< number of particles to be started per point
+    INTEGER(iwp) ::  particle_file_count = 0                      !< can be removed ?
+    INTEGER(iwp) ::  radius_classes = 20                          !< number of radius classes to be used in the collision efficiency table
+    INTEGER(iwp) ::  sort_count = 0                               !< counter for sorting particles
+    INTEGER(iwp) ::  splitting_factor = 2                         !< splitting factor 
+    INTEGER(iwp) ::  splitting_factor_max = 5                     !< maximum splittig factor
+    INTEGER(iwp) ::  step_dealloc = 100                           !< number of timesteps after which particle memory is deallocated
+    INTEGER(iwp) ::  sum_merge_drp = 0                            !< sum of merged super droplets
+    INTEGER(iwp) ::  sum_new_particles = 0                        !< sum of created particles (in splitting algorithm)
+    INTEGER(iwp) ::  total_number_of_particles                    !< total number of particles in the whole model domain           
+    INTEGER(iwp) ::  trlp_count_sum                               !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trlp_count_recv_sum                          !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trrp_count_sum                               !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trrp_count_recv_sum                          !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trsp_count_sum                               !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trsp_count_recv_sum                          !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trnp_count_sum                               !< parameter for particle exchange of PEs
+    INTEGER(iwp) ::  trnp_count_recv_sum                          !< parameter for particle exchange of PEs
+
+    INTEGER(iwp), PARAMETER ::  max_number_of_particle_groups = 10 !< maximum allowed number of particle groups
+
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  prt_count  !< 3d array of number of particles of every grid box
+
+    LOGICAL ::  all_or_nothing = .FALSE.                  !< flag for collision algorithm
+    LOGICAL ::  average_impact = .FALSE.                  !< flag for collision algortihm
+    LOGICAL ::  curvature_solution_effects = .FALSE.      !< parameter to consider solution and curvature effects on the equilibrium vapor pressure of cloud droplets                      
+    LOGICAL ::  deallocate_memory = .TRUE.                !< parameter to enable deallocation of unused memory                   
+    LOGICAL ::  hall_kernel = .FALSE.                     !< flag for collision kernel
+    LOGICAL ::  init_aerosol_probabilistic = .FALSE.      !< parameter to steer the initialization of the aerosol spectrum
+    LOGICAL ::  merging = .FALSE.                         !< parameter to enable merging of super droplets
+    LOGICAL ::  monodisperse_aerosols = .FALSE.           !< parameter to steer the initialization of the aerosol spectrum
+    LOGICAL ::  particle_advection = .FALSE.              !< parameter to steer the advection of particles
+    LOGICAL ::  random_start_position = .FALSE.           !< parameter to initialize particles on random positon (within one grid box)                   
+    LOGICAL ::  read_particles_from_restartfile = .TRUE.  !< read particle data from the previous run                   
+    LOGICAL ::  seed_follows_topography = .FALSE.         !< heights of initial particles are interpreted relative to the given topography
+    LOGICAL ::  splitting = .FALSE.                       !< parameter to enable the splitting of super droplets
+    LOGICAL ::  uniform_particles = .TRUE.                !< can be removed?
+    LOGICAL ::  use_kernel_tables = .FALSE.               !< parameter, which turns on the use of precalculated collision kernels
+    LOGICAL ::  use_sgs_for_particles = .FALSE.           !< parameter to use sgs velocities for particles 
+    LOGICAL ::  wang_kernel = .FALSE.                     !< flag for collision kernel
+    LOGICAL ::  write_particle_statistics = .FALSE.       !< switch on/off output of particle information (statistical information)
+                
     LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
-                vertical_particle_advection = .TRUE.
-
-    REAL(wp) ::  alloc_factor = 20.0_wp, c_0 = 3.0_wp,                         &
-                 dt_min_part = 0.0002_wp, dt_prel = 9999999.9_wp,              &
-                 dt_write_particle_data = 9999999.9_wp,                        &
-                 end_time_prel = 9999999.9_wp,                                 &
-                 initial_weighting_factor = 1.0_wp,                            &
-                 number_concentration     = -1.0_wp,                           &
-                 particle_advection_start = 0.0_wp,                            &
-                 sgs_wf_part, time_prel = 0.0_wp, time_sort_particles = 0.0_wp,&
-                 time_write_particle_data = 0.0_wp, z0_av_global
-
-    REAL(wp), DIMENSION(max_number_of_particle_groups) ::                      &
-                 density_ratio = 9999999.9_wp, pdx = 9999999.9_wp,             &
-                 pdy = 9999999.9_wp, pdz = 9999999.9_wp, psb = 9999999.9_wp,   &
-                 psl = 9999999.9_wp, psn = 9999999.9_wp, psr = 9999999.9_wp,   &
-                 pss = 9999999.9_wp, pst = 9999999.9_wp, radius = 9999999.9_wp
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0
-
+                vertical_particle_advection = .TRUE.              !< Switch on/off vertical particle transport
+
+    REAL(wp) ::  alloc_factor = 20.0_wp                           !< memory allocated additionally to the memory needed for initial particles at a given grid cell
+    REAL(wp) ::  c_0 = 3.0_wp                                     !< parameter for lagrangian timescale
+    REAL(wp) ::  dt_min_part = 0.0002_wp                          !< minimum particle time step when SGS velocities are used (s)
+    REAL(wp) ::  dt_prel = 9999999.9_wp                           !< temporal interval at which particles are to be released from a particle source (in s)
+    REAL(wp) ::  dt_write_particle_data = 9999999.9_wp            !< temporal interval for output of particle data (in s)             
+    REAL(wp) ::  end_time_prel = 9999999.9_wp                     !< time of the last release of particles (in s).            
+    REAL(wp) ::  initial_weighting_factor = 1.0_wp                !< weighting factor used for initialization
+    REAL(wp) ::  molecular_weight_of_solute = 0.05844_wp          !< mol. m. NaCl (kg mol-1)
+    REAL(wp) ::  molecular_weight_of_water = 0.01801528_wp        !< mol. m. H2O (kg mol-1)
+    REAL(wp) ::  n1 = 100.0_wp                                    !< number concentration of the first log-normal distribution (dry aerosol initialization)
+    REAL(wp) ::  n2 = 0.0_wp                                      !< see n1
+    REAL(wp) ::  n3 = 0.0_wp                                      !< see n1 
+    REAL(wp) ::  number_concentration = -1.0_wp                   !< initial particle number concentration (in units of 1/cm^-3)
+    REAL(wp) ::  particle_advection_start = 0.0_wp                !< time of the first release of particles (in s).
+    REAL(wp) ::  radius_merge = 1.0E-7_wp                         !< particles with a smaller radius are merged if merging is set
+    REAL(wp) ::  radius_split = 40.0E-6_wp                        !< particles with a bigger radius are splitted if splittig is set
+    REAL(wp) ::  rho_s = 2165.0_wp                                !< density of NaCl (kg m-3) 
+    REAL(wp) ::  rm1 = 0.05E-6_wp                                 !< mode radius of the first log-normal distribution (dry aerosol initialization)
+    REAL(wp) ::  rm2 = 0.05E-6_wp                                 !< see rm1
+    REAL(wp) ::  rm3 = 0.05E-6_wp                                 !< see rm2
+    REAL(wp) ::  s1 = 2.0_wp                                      !< geometric standard deviation of the first log-normal distribution (dry aerosol initialization)
+    REAL(wp) ::  s2 = 2.0_wp                                      !< see s1 
+    REAL(wp) ::  s3 = 2.0_wp                                      !< see s1
+    REAL(wp) ::  sgs_wf_part                                      !< parameter for sgs 
+    REAL(wp) ::  time_prel = 0.0_wp                               !< time for particle release
+    REAL(wp) ::  time_sort_particles = 0.0_wp                     !< can be removed?
+    REAL(wp) ::  time_write_particle_data = 0.0_wp                !< write particle data at current time on file
+    REAL(wp) ::  vanthoff = 2.0_wp                                !< van't Hoff factor for NaCl
+    REAL(wp) ::  weight_factor_merge = -1.0_wp                    !< critical weighting factor for merging
+    REAL(wp) ::  weight_factor_split = -1.0_wp                    !< critical weighting factor for splitting
+    REAL(wp) ::  z0_av_global                                     !< horizontal mean value of z0
+
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  density_ratio = 9999999.9_wp  !< ratio of the density of the fluid and the density of the particles   
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  pdx = 9999999.9_wp            !< distance along x between particles within a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  pdy = 9999999.9_wp            !< distance along y between particles within a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  pdz = 9999999.9_wp            !< distance along z between particles within a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  psb = 9999999.9_wp            !< bottom edge of a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  psl = 9999999.9_wp            !< left edge of a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  psn = 9999999.9_wp            !< rear ("north") edge of a particle source (in m)
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  psr = 9999999.9_wp            !< right edge of a particle source (in m). 
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  pss = 9999999.9_wp            !< front ("south") edge of a particle source (in m). 
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  pst = 9999999.9_wp            !< top edge of a particle source (in m).
+    REAL(wp), DIMENSION(max_number_of_particle_groups) ::  radius = 9999999.9_wp         !< particle radius (in m).
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)
+
+    
     TYPE particle_type
         SEQUENCE
@@ -181,18 +246,4 @@
     TYPE(block_offset_def), DIMENSION(0:7)         ::  block_offset
 
-!
-!-- Lagrangian cloud model constants (especially for steering aerosols)
-    REAL(wp) ::  molecular_weight_of_solute = 0.05844_wp    !< mol. m. NaCl (kg mol-1)
-    REAL(wp) ::  molecular_weight_of_water = 0.01801528_wp  !< mol. m. H2O (kg mol-1)
-    REAL(wp) ::  rho_s = 2165.0_wp                          !< density of NaCl (kg m-3) 
-    REAL(wp) ::  vanthoff = 2.0_wp                          !< van't Hoff factor for NaCl
-
-    REAL(wp) ::  n1 = 100.0_wp, s1 = 2.0_wp, rm1 = 0.05E-6_wp, &
-                 n2 =   0.0_wp, s2 = 2.0_wp, rm2 = 0.05E-6_wp, &
-                 n3 =   0.0_wp, s3 = 2.0_wp, rm3 = 0.05E-6_wp
-
-    LOGICAL  ::  monodisperse_aerosols      = .FALSE.       !< initialize monodisperse aerosols
-    LOGICAL  ::  init_aerosol_probabilistic = .FALSE.       !< how to initialize aerosol spectra
-
     SAVE
 
Index: /palm/trunk/SOURCE/package_parin.f90
===================================================================
--- /palm/trunk/SOURCE/package_parin.f90	(revision 2262)
+++ /palm/trunk/SOURCE/package_parin.f90	(revision 2263)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implemented splitting and merging algorithm
+! 
+! 2183 2017-03-17 14:29:15Z schwenkel
 !
 ! 2182 2017-03-17 14:27:40Z schwenkel
@@ -167,15 +170,17 @@
                dissipation_classes, dt_min_part, dt_prel,                      &
                dt_write_particle_data, end_time_prel, initial_weighting_factor,&
-               init_aerosol_probabilistic, min_nr_particle,                    &
-               monodisperse_aerosols, n1, n2, n3, number_concentration,        & 
-               number_of_particle_groups, number_particles_per_gridbox,        &
-               particles_per_point, particle_advection,                        &
-               particle_advection_start,                                       &
-               pdx, pdy, pdz, psb, psl, psn, psr, pss, pst, radius,            &
-               radius_classes, random_start_position,                          &
+               init_aerosol_probabilistic, max_number_particles_per_gridbox,   & 
+               merging, min_nr_particle, monodisperse_aerosols, n1, n2, n3,    &
+               number_concentration, number_of_particle_groups,                &
+               number_particles_per_gridbox, particles_per_point,              &
+               particle_advection, particle_advection_start, pdx, pdy, pdz,    &
+               psb, psl, psn, psr, pss, pst, radius, radius_classes,           &
+               radius_merge, radius_split, random_start_position,              &
                read_particles_from_restartfile, rm1, rm2, rm3,                 &
-               seed_follows_topography, step_dealloc, s1, s2, s3,              &
-               use_sgs_for_particles, vertical_particle_advection,             &
-               write_particle_statistics
+               seed_follows_topography, splitting, splitting_factor,           & 
+               splitting_factor_max, splitting_function, splitting_mode,       &
+               step_dealloc, s1, s2, s3, use_sgs_for_particles,                &
+               vertical_particle_advection, weight_factor_merge,               &
+               weight_factor_split, write_particle_statistics                
 
     IMPLICIT NONE
@@ -210,4 +215,5 @@
                                   end_time_prel, initial_weighting_factor,     &
                                   init_aerosol_probabilistic,                  &
+                                  max_number_particles_per_gridbox, merging,   &
                                   min_nr_particle, monodisperse_aerosols,      &
                                   n1, n2, n3, number_concentration,            &
@@ -218,10 +224,15 @@
                                   particle_maximum_age, pdx, pdy, pdz, psb,    &
                                   psl, psn, psr, pss, pst, radius,             &
-                                  radius_classes, random_start_position,       &
+                                  radius_classes, radius_merge, radius_split,  &   
+                                  random_start_position,                       &
                                   read_particles_from_restartfile,             &
                                   rm1, rm2, rm3,                               &
-                                  seed_follows_topography, step_dealloc,       &
-                                  s1, s2, s3, use_sgs_for_particles,           &
+                                  seed_follows_topography,                     &
+                                  splitting, splitting_factor,                 &
+                                  splitting_factor_max, splitting_function,    &
+                                  splitting_mode, step_dealloc, s1, s2, s3,    &
+                                  use_sgs_for_particles,                       &
                                   vertical_particle_advection,                 &
+                                  weight_factor_merge, weight_factor_split,    &
                                   write_particle_statistics
 
