Changeset 828 for palm/trunk

Timestamp:

Feb 21, 2012 12:00:36 PM (12 years ago)

Author:

raasch

Message:

New:
---

Changed:

---

Optimization of collision kernels. Collision tables can be calculated once at
simulation start for defined radius (and dissipation) classes instead of
re-calculating them at every timestep and for the particle ensemble in
every gridbox.
For this purpose the particle feature color is renamed class.
New parpar parameters radius_classes and dissipation_classes.
(Makefile, advec_particles, check_parameters, data_output_dvrp, header, init_particles, lpm_collision_kernels, modules, package_parin, set_particle_attributes)

Lower limit for droplet radius changed from 1E-7 to 1E-8.
(advec_particles)

Complete re-formatting of collision code (including changes in names of
variables, modules and subroutines).
(advec_particles, lpm_collision_kernels)

Errors:

---

Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
(advec_particles, lpm_collision_kernels)

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (3 diffs)
• advec_particles.f90 (modified) (17 diffs)
• check_parameters.f90 (modified) (4 diffs)
• data_output_dvrp.f90 (modified) (2 diffs)
• header.f90 (modified) (3 diffs)
• init_particles.f90 (modified) (7 diffs)
• lpm_collision_kernels.f90 (modified) (10 diffs)
• modules.f90 (modified) (5 diffs)
• package_parin.f90 (modified) (3 diffs)
• set_particle_attributes.f90 (modified) (4 diffs)

palm/trunk/SOURCE/Makefile

@@ -4 +4 @@
 # Current revisions:
 # ------------------
-# 
+# init_particles depends on lpm_collision_kernels
 #
 # Former revisions:
@@ -66 +66 @@
 
 RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
-        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
-        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
-        advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
-        average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
-        calc_liquid_water_content.f90 calc_precipitation.f90 \
-        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
-        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
-        coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
-        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
-        data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \
-        data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \
-        diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \
-        disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \
-        exchange_horiz.f90 exchange_horiz_2d.f90 \
-        fft_xy.f90 flow_statistics.f90 global_min_max.f90 \
-        header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \
-        init_1d_model.f90 init_3d_model.f90 init_advec.f90 \
-        init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 init_grid.f90 \
-        init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
-        init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
+        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
+        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
+        advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
+        average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
+        calc_liquid_water_content.f90 calc_precipitation.f90 \
+        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
+        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
+        coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
+        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
+        data_output_ptseries.f90 data_output_spectra.f90 \
+        data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \
+        diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 diffusion_v.f90 \
+        diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
+        disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
+        exchange_horiz_2d.f90 fft_xy.f90 flow_statistics.f90 \
+        global_min_max.f90 header.f90 impact_of_latent_heat.f90 \
+        inflow_turbulence.f90 init_1d_model.f90 init_3d_model.f90 \
+        init_advec.f90 init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 \
+        init_grid.f90 init_masks.f90 init_ocean.f90 init_particles.f90 \
+        init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
         interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
-        local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90 \
-        lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \
-        package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
-        plant_canopy_model.f90 poisfft.f90 \
+        local_stop.f90 local_system.f90 local_tremain.f90 \
+        local_tremain_ini.f90 lpm_collision_kernels.f90 message.f90 \
+        modules.f90 netcdf.f90 package_parin.f90 palm.f90 parin.f90 \
+        particle_boundary_conds.f90 plant_canopy_model.f90 poisfft.f90 \
         poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
         production_e.f90 prognostic_equations.f90 random_function.f90 \
         random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
         set_particle_attributes.f90 set_slicer_attributes_dvrp.f90 \
-        singleton.f90 sor.f90 spline_x.f90 \
+        singleton.f90 sor.f90 spline_x.f90 \
         spline_y.f90 spline_z.f90 subsidence.f90 \
-        sum_up_3d_data.f90 surface_coupler.f90 \
+        sum_up_3d_data.f90 surface_coupler.f90 \
         swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
-        time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
-        transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \
+        time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
+        transpose.f90 user_3d_data_averaging.f90 user_actions.f90 \
         user_additional_routines.f90 user_advec_particles.f90 \
         user_check_data_output.f90 user_check_data_output_pr.f90 \
         user_check_parameters.f90 user_data_output_2d.f90 \
-        user_data_output_3d.f90 user_data_output_dvrp.f90 user_data_output_mask.f90 \
-        user_define_netcdf_grid.f90 user_dvrp_coltab.f90 user_header.f90 \
-        user_init.f90 user_init_3d_model.f90 user_init_grid.f90 \
-        user_init_particles.f90 user_init_plant_canopy.f90 \
-        user_last_actions.f90 user_module.f90 user_parin.f90 \
-        user_particle_attributes.f90 user_read_restart_data.f90 \
+        user_data_output_3d.f90 user_data_output_dvrp.f90 \
+        user_data_output_mask.f90 user_define_netcdf_grid.f90 \
+        user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
+        user_init_3d_model.f90 user_init_grid.f90 user_init_particles.f90 \
+        user_init_plant_canopy.f90 user_last_actions.f90 user_module.f90 \
+        user_parin.f90 user_particle_attributes.f90 user_read_restart_data.f90 \
         user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
-        write_3d_binary.f90 write_compressed.f90 \
-        write_var_list.f90
+        write_3d_binary.f90 write_compressed.f90 write_var_list.f90
 
 OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
@@ -252 +251 @@
 init_masks.o: modules.o
 init_ocean.o: modules.o eqn_state_seawater.o
-init_particles.o: modules.o random_function.o
+init_particles.o: modules.o lpm_collision_kernels.o random_function.o
 init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
 init_pt_anomaly.o: modules.o

palm/trunk/SOURCE/advec_particles.f90

@@ -4 +4 @@
 ! Current revisions:
 ! ------------------
-! 
+! fast hall/wang kernels with fixed radius/dissipation classes added,
+! particle feature color renamed class, routine colker renamed
+! recalculate_kernel,
+! lower limit for droplet radius changed from 1E-7 to 1E-8
+!
+! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
 !
 ! Former revisions:
@@ -144 +149 @@
     IMPLICIT NONE
 
-    INTEGER ::  agp, deleted_particles, deleted_tails, i, ie, ii, inc,         &
-                internal_timestep_count, is, j,  &
-                jj, js, jtry, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi,   &
+    INTEGER ::  agp, deleted_particles, deleted_tails, eclass, i, ie, ii, inc, &
+                internal_timestep_count, is, j, jj, js, jtry, k, kk, kw, m, n, &
+                nc, nd, nn, num_gp, pse, psi, rclass_l, rclass_s,  &
                 tlength, trlp_count, trlp_count_sum, trlp_count_recv,          &
                 trlp_count_recv_sum, trlpt_count, trlpt_count_recv,            &
@@ -173 +178 @@
                 lagr_timescale, lw_max, 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, r_ros, r_ros_ini, sigma, 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
+                ql_int_u, random_gauss, r_ros, r_ros_ini, &
+                sigma, 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
 !
 !-- Parameters for Rosenbrock method
@@ -190 +195 @@
 
     REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
-    REAL, DIMENSION(:,:), ALLOCATABLE ::  kern
 
     REAL    ::  location(1:30,1:3)
@@ -356 +360 @@
                              thermal_conductivity_l +                         &
                              rho_l * r_v * t_int / diff_coeff_l / e_s )
-             IF ( arg < 1.0E-14 )  THEN
-                new_r = 1.0E-7
+             IF ( arg < 1.0E-16 )  THEN
+                new_r = 1.0E-8
              ELSE
                 new_r = SQRT( arg )
@@ -550 +554 @@
              CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 )
           ENDIF
-          particles(n)%radius = new_r
 
 !
 !--       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 *             &
-                                      particles(n)%radius**3
+          ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3
+
+          particles(n)%radius = new_r
+
+!
+!--       Determine radius class of the particle needed for collision
+          IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables ) &
+          THEN
+             particles(n)%class = ( LOG( new_r ) - rclass_lbound ) /  &
+                                  ( rclass_ubound - rclass_lbound ) * &
+                                  radius_classes
+             particles(n)%class = MIN( particles(n)%class, radius_classes )
+             particles(n)%class = MAX( particles(n)%class, 1 )
+          ENDIF
+
        ENDDO
        CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
@@ -601 +617 @@
                    pse = psi + prt_count(k,j,i)-1
 
-                   IF ( hall_kernel  .OR.  wang_kernel )  THEN
-
-                      ALLOCATE( kern(psi:pse,psi:pse) )
-
-!
-!--                   Calculate collision kernel for all particles in grid box
-                      CALL colker( i, j, k, kern )
-!
-!--                   Collison kernel is calculated in cm**3/s but needed
-!--                   in m**3/s
-                      kern = kern * 1.0E-6   ! to be moved to colker
+!
+!--                Now apply the different kernels
+                   IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
+                        use_kernel_tables )  THEN
+!
+!--                   Fast method with pre-calculated efficiencies for
+!--                   discrete radius- and dissipation-classes.
+!
+!--                   Determine dissipation class index of this gridbox
+                      eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
+                                    dissipation_classes ) + 1
+                      IF ( hall_kernel .OR. diss(k,j,i) * 1.0E4 < 0.001 )  THEN
+                         eclass = 0   ! Hall kernel is used
+                      ELSE
+                         eclass = MIN( dissipation_classes, eclass )
+                      ENDIF
 
                       DO  n = pse, psi+1, -1
@@ -617 +638 @@
                          integral = 0.0
                          lw_max   = 0.0
-
+                         rclass_l = particles(n)%class
 !
 !--                      Calculate growth of collector particle radius using all
@@ -623 +644 @@
                          DO  is = psi, n-1
 
-                            integral = integral +                              &
-                                       ( particles(is)%radius**3 * kern(n,is) *&
+                            rclass_s = particles(is)%class
+                            integral = integral +                             &
+                                       ( particles(is)%radius**3 *            &
+                                         ckernel(rclass_l,rclass_s,eclass) *  &
                                          particles(is)%weight_factor )
 !
@@ -653 +676 @@
 
                          ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
-
+!-- can this case really happen??
                             DO  is = psi, n-1
 
@@ -667 +690 @@
                             DO  is = psi, n-1
 
+                               rclass_s = particles(is)%class
                                particles(is)%weight_factor = &
                                   particles(is)%weight_factor - &
-                                  ( ( kern(n,is) * particles(is)%weight_factor &
+                                  ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor &
                                     / integral ) * delta_v )
 
@@ -699 +723 @@
                       ENDDO
 
-                      DEALLOCATE( kern )
+                   ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
+                            .NOT. use_kernel_tables )  THEN
+!
+!--                   Collision efficiencies are calculated for every new
+!--                   grid box. First, allocate memory for kernel table.
+!--                   Third dimension is 1, because table is re-calculated for
+!--                   every new dissipation value.
+                      ALLOCATE( ckernel(prt_start_index(k,j,i):                &
+                               prt_start_index(k,j,i)+prt_count(k,j,i)-1,      &
+                               prt_start_index(k,j,i):                         &
+                               prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
+!
+!--                   Now calculate collision efficiencies for this box
+                      CALL recalculate_kernel( i, j, k )
+
+                      DO  n = pse, psi+1, -1
+
+                         integral = 0.0
+                         lw_max   = 0.0
+!
+!--                      Calculate growth of collector particle radius using all
+!--                      droplets smaller than current droplet 
+                         DO  is = psi, n-1
+
+                            integral = integral + particles(is)%radius**3 * &
+                                                  ckernel(n,is,1) *         &
+                                                  particles(is)%weight_factor
+!
+!--                         Calculate volume of liquid water of the collected
+!--                         droplets which is the maximum liquid water available
+!--                         for droplet growth   
+                            lw_max = lw_max + ( particles(is)%radius**3 * &
+                                                particles(is)%weight_factor )
+
+                         ENDDO
+
+!
+!--                      Change in radius of collector droplet due to collision
+                         delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
+                                   integral * dt_3d * ddx * ddy / dz 
+
+!
+!--                      Change in volume of collector droplet due to collision
+                         delta_v = particles(n)%weight_factor                  &
+                                   * ( ( particles(n)%radius + delta_r )**3    &
+                                       - particles(n)%radius**3 )
+
+                         IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
+!-- replace by message call
+                            PRINT*, 'Particle has grown to much because the',  &
+                                    ' change of volume of particle is larger', &
+                                    ' than liquid water available!'
+
+                         ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
+!-- can this case really happen??
+                            DO  is = psi, n-1
+
+                               particles(is)%weight_factor = 0.0
+                               particle_mask(is)  = .FALSE.
+                               deleted_particles = deleted_particles + 1
+
+                            ENDDO
+
+                         ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
+!
+!--                         Calculate new weighting factor of collected droplets
+                            DO  is = psi, n-1
+
+                               particles(is)%weight_factor =                   &
+                                      particles(is)%weight_factor -            &
+                                      ( ckernel(n,is,1) / integral * delta_v * &
+                                        particles(is)%weight_factor )
+
+                               IF ( particles(is)%weight_factor < 0.0 )  THEN
+                                     WRITE( message_string, * ) 'negative ',   &
+                                                        'weighting factor: ',  &
+                                                     particles(is)%weight_factor
+                                     CALL message( 'advec_particles', '', 2, 2,&
+                                                   -1, 6, 1 )
+
+                               ELSEIF ( particles(is)%weight_factor == 0.0 )  &
+                               THEN
+
+                                   particles(is)%weight_factor = 0.0
+                                   particle_mask(is)  = .FALSE.
+                                   deleted_particles = deleted_particles + 1
+
+                               ENDIF
+
+                            ENDDO
+
+                         ENDIF
+
+                         particles(n)%radius = particles(n)%radius + delta_r
+                         ql_vp(k,j,i) = ql_vp(k,j,i) + &
+                                                    particles(n)%weight_factor &
+                                                    * particles(n)%radius**3
+
+                      ENDDO
+
+                      DEALLOCATE( ckernel )
 
                    ELSEIF ( palm_kernel )  THEN
@@ -2376 +2500 @@
                       particle_tail_coordinates(1,2,nn) = particles(n)%y
                       particle_tail_coordinates(1,3,nn) = particles(n)%z
-                      particle_tail_coordinates(1,4,nn) = particles(n)%color
+                      particle_tail_coordinates(1,4,nn) = particles(n)%class
                       particles(n)%tailpoints = 1
                       IF ( minimum_tailpoint_distance /= 0.0 )  THEN
@@ -2382 +2506 @@
                          particle_tail_coordinates(2,2,nn) = particles(n)%y
                          particle_tail_coordinates(2,3,nn) = particles(n)%z
-                         particle_tail_coordinates(2,4,nn) = particles(n)%color
+                         particle_tail_coordinates(2,4,nn) = particles(n)%class
                          particle_tail_coordinates(1:2,5,nn) = 0.0
                          particles(n)%tailpoints = 2
@@ -3743 +3867 @@
 
 !
-!-- Set particle attributes
-    CALL set_particle_attributes
+!-- Set particle attributes.
+!-- Feature is not available if collision is activated, because the respective
+!-- particle attribute (class) is then used for storing the particle radius
+!-- class.
+    IF ( collision_kernel == 'none' )  CALL set_particle_attributes
 
 !
 !-- Set particle attributes defined by the user
     CALL user_particle_attributes
-!    WRITE ( 9, * ) '*** advec_particles: ##10'
-!    CALL local_flush( 9 )
-!    DO  n = 1, number_of_particles
-!       IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
-!       THEN
-!          WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
-!          WRITE (9,*) '    id=',particles(n)%tail_id,' of (',number_of_tails,')'
-!          CALL local_flush( 9 )
-!          CALL MPI_ABORT( comm2d, 9999, ierr )
-!       ENDIF
-!    ENDDO
 
 !
@@ -3810 +3926 @@
              particle_tail_coordinates(1,2,nn) = particles(n)%y
              particle_tail_coordinates(1,3,nn) = particles(n)%z
-             particle_tail_coordinates(1,4,nn) = particles(n)%color
+             particle_tail_coordinates(1,4,nn) = particles(n)%class
 !             WRITE ( 9, * ) '*** advec_particles: ##10.4'
 !             CALL local_flush( 9 )
@@ -4058 +4174 @@
     CALL handle_netcdf_error( 'output_particles_netcdf', 15 )
 
-    nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%color,      &
+    nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class,      &
                             start = (/ 1, prt_time_count /),                  &
                             count = (/ maximum_number_of_particles /) )

palm/trunk/SOURCE/check_parameters.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! check of collision_kernel extended
 !
 ! Former revisions:
@@ -716 +716 @@
     SELECT CASE ( TRIM( collision_kernel ) )
 
-       CASE ( 'hall' )
+       CASE ( 'hall', 'hall_fast' )
           hall_kernel = .TRUE.
 
@@ -722 +722 @@
           palm_kernel = .TRUE.
 
-       CASE ( 'wang' )
+       CASE ( 'wang', 'wang_fast' )
           wang_kernel = .TRUE.
 
@@ -734 +734 @@
 
     END SELECT
+    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
 
 

palm/trunk/SOURCE/data_output_dvrp.f90

@@ -32 +32 @@
 ! Current revisions:
 ! -----------------
-! 
+! particle feature color renamed class
 !
 ! Former revisions:
@@ -226 +226 @@
                    p_y(k)   = particles(n)%y * superelevation_y
                    p_z(k)   = particles(n)%z * superelevation
-                   p_c(k)   = particles(n)%color
+                   p_c(k)   = particles(n)%class
                 ENDIF
              ENDDO

palm/trunk/SOURCE/header.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! output of dissipation_classes + radius_classes
 ! 
 ! Former revisions:
@@ -1483 +1483 @@
        IF ( collision_kernel /= 'none' )  THEN
           WRITE ( io, 435 )  TRIM( collision_kernel )
-       ELSE
-          WRITE ( io, 436 )
+          IF ( collision_kernel(6:9) == 'fast' )  THEN
+             WRITE ( io, 436 )  radius_classes, dissipation_classes
+          ENDIF
+       ELSE
+          WRITE ( io, 437 )
        ENDIF
     ENDIF
@@ -1964 +1967 @@
                  ' droplets < 1.0E-6 m')
 435 FORMAT ('    Droplet collision is handled by ',A,'-kernel')
-436 FORMAT ('    Droplet collision is switched off')
+436 FORMAT ('       Fast kernel with fixed radius- and dissipation classes ', &
+                    'are used'/ &
+            '          number of radius classes:       ',I3,'    interval ', &
+                       '[1.0E-6,2.0E-4] m'/ &
+            '          number of dissipation classes:   ',I2,'    interval ', &
+                       '[0,1000] cm**2/s**3')
+437 FORMAT ('    Droplet collision is switched off')
 450 FORMAT (//' LES / Turbulence quantities:'/ &
               ' ---------------------------'/)

palm/trunk/SOURCE/init_particles.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! call of init_kernels, particle feature color renamed class
 !
 ! Former revisions:
@@ -72 +72 @@
     USE grid_variables
     USE indices
+    USE lpm_collision_kernels_mod
     USE particle_attributes
     USE pegrid
@@ -148 +149 @@
        IF ( pdz(j) == 9999999.9  .OR.  pdz(j) == 0.0 )  pdz(j) = pdz(j-1)
     ENDDO
+
+!
+!-- Initialize collision kernels
+    IF ( collision_kernel /= 'none' )  CALL init_kernels
 
 !
@@ -343 +348 @@
                             particles(n)%radius      = particle_groups(i)%radius
                             particles(n)%weight_factor =initial_weighting_factor
-                            particles(n)%color         = 1
+                            particles(n)%class         = 1
                             particles(n)%group         = i
                             particles(n)%tailpoints    = 0
@@ -470 +475 @@
              IF ( MINVAL( particles(1:n)%dvrp_psize  ) ==     &
                   MAXVAL( particles(1:n)%dvrp_psize  )  .AND. &
-                  MINVAL( particles(1:n)%color ) ==     &
-                  MAXVAL( particles(1:n)%color ) )  THEN
+                  MINVAL( particles(1:n)%class ) ==     &
+                  MAXVAL( particles(1:n)%class ) )  THEN
                 uniform_particles_l = .TRUE.
              ELSE
@@ -532 +537 @@
                 particle_tail_coordinates(1,2,nn) = particles(n)%y
                 particle_tail_coordinates(1,3,nn) = particles(n)%z
-                particle_tail_coordinates(1,4,nn) = particles(n)%color
+                particle_tail_coordinates(1,4,nn) = particles(n)%class
                 particles(n)%tailpoints = 1
                 IF ( minimum_tailpoint_distance /= 0.0 )  THEN
@@ -538 +543 @@
                    particle_tail_coordinates(2,2,nn) = particles(n)%y
                    particle_tail_coordinates(2,3,nn) = particles(n)%z
-                   particle_tail_coordinates(2,4,nn) = particles(n)%color
+                   particle_tail_coordinates(2,4,nn) = particles(n)%class
                    particle_tail_coordinates(1:2,5,nn) = 0.0
                    particles(n)%tailpoints = 2

palm/trunk/SOURCE/lpm_collision_kernels.f90

@@ -1 @@
-MODULE lpm_collision_kernels_mod
+ MODULE lpm_collision_kernels_mod
 
 !------------------------------------------------------------------------------!
 ! Current revisions:
 ! -----------------
-! 
+! code has been completely reformatted, routine colker renamed
+! recalculate_kernel,
+! routine init_kernels added, radius is now communicated to the collision
+! routines by array radclass
+!
+! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
 !
 ! Former revisions:
@@ -24 +29 @@
 ! Description:
 ! ------------
-! This routine calculates the effect of (SGS) turbulence on the collision
-! efficiency of droplets.
-! It is based on the original kernel developed by Wang (...)
+! This module calculates collision efficiencies either due to pure gravitational
+! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
+! including the effects of (SGS) turbulence (Wang kernel, see Wang and
+! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
+! provided by L.-P. Wang but is substantially reformatted and speed optimized
+! here.
+!
+! ATTENTION:
+! Physical quantities (like g, densities, etc.) used in this module still
+! have to be adjusted to those values used in the main PALM code.
+! Also, quantities in CGS-units should be converted to SI-units eventually.
 !------------------------------------------------------------------------------!
 
@@ -33 +46 @@
     USE constants
     USE particle_attributes
+    USE pegrid
+
 
     IMPLICIT NONE
@@ -38 +53 @@
     PRIVATE
 
-    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
-
-    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
-
-    REAL, SAVE :: epsilon, eps2, urms, urms2
-
-    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
-    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
+    PUBLIC  ckernel, init_kernels,  rclass_lbound, rclass_ubound, &
+            recalculate_kernel            
+
+    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
+
+    REAL, DIMENSION(:),   ALLOCATABLE ::  epsclass, radclass, winf
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  ec, ecf, gck, hkernel, hwratio
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  ckernel
+
+    SAVE
 
 !
 !-- Public interfaces
-    INTERFACE colker
-       MODULE PROCEDURE colker
-    END INTERFACE colker
-
-    INTERFACE effic
-       MODULE PROCEDURE effic
-    END INTERFACE effic
-
-    INTERFACE fallg
-       MODULE PROCEDURE fallg
-    END INTERFACE fallg
-
-    INTERFACE phi
-       MODULE PROCEDURE phi
-    END INTERFACE phi
-
-    INTERFACE turbsd
-       MODULE PROCEDURE turbsd
-    END INTERFACE turbsd
-
-    INTERFACE turb_enhance_eff
-       MODULE PROCEDURE turb_enhance_eff
-    END INTERFACE turb_enhance_eff
-
-    INTERFACE zhi
-       MODULE PROCEDURE zhi
-    END INTERFACE zhi
+    INTERFACE init_kernels
+       MODULE PROCEDURE init_kernels
+    END INTERFACE init_kernels
+
+    INTERFACE recalculate_kernel
+       MODULE PROCEDURE recalculate_kernel
+    END INTERFACE recalculate_kernel
+
 
     CONTAINS
 
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision kernel
-!------------------------------------------------------------------------------!
-    SUBROUTINE colker( i1, j1, k1, kernel )
+
+    SUBROUTINE init_kernels
+!------------------------------------------------------------------------------!
+! Initialization of the collision efficiency matrix with fixed radius and
+! dissipation classes, calculated at simulation start only.
+!------------------------------------------------------------------------------!
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+
+!
+!--    Calculate collision efficiencies for fixed radius- and dissipation
+!--    classes
+       IF ( collision_kernel(6:9) == 'fast' )  THEN
+
+          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
+                    0:dissipation_classes), epsclass(1:dissipation_classes), &
+                    radclass(1:radius_classes) )
+
+!
+!--       Calculate the radius class bounds with logarithmic distances
+!--       in the interval [1.0E-6, 2.0E-4] m
+          rclass_lbound = LOG( 1.0E-6 )
+          rclass_ubound = LOG( 2.0E-4 )
+          radclass(1)   = 1.0E-6
+          DO  i = 2, radius_classes
+             radclass(i) = EXP( rclass_lbound +                                &
+                                ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
+                                ( radius_classes - 1.0 ) )
+!             IF ( myid == 0 )  THEN
+!                PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6
+!             ENDIF
+          ENDDO
+!
+!--       Collision routines expect radius to be in cm
+          radclass = radclass * 100.0
+
+!
+!--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
+          DO  i = 1, dissipation_classes
+             epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
+!             IF ( myid == 0 )  THEN
+!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
+!             ENDIF
+          ENDDO
+!
+!--       Calculate collision efficiencies of the Wang/ayala kernel
+          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
+                    ecf(1:radius_classes,1:radius_classes), &
+                    gck(1:radius_classes,1:radius_classes), &
+                    winf(1:radius_classes) )
+
+          DO  k = 1, dissipation_classes
+
+             epsilon = epsclass(k)
+             urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
+
+             CALL turbsd
+             CALL turb_enhance_eff
+             CALL effic
+
+             DO  j = 1, radius_classes
+                DO  i = 1, radius_classes
+                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
+                ENDDO
+             ENDDO
+
+          ENDDO
+
+!
+!--       Calculate collision efficiencies of the Hall kernel
+          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
+                    hwratio(1:radius_classes,1:radius_classes) )
+
+          CALL fallg
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
+                                  * ec(i,j) * ABS( winf(j) - winf(i) )
+                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
+              ENDDO
+          ENDDO
+
+!
+!--       Test output of efficiencies
+          IF ( j == -1 )  THEN
+
+             PRINT*, '*** Hall kernel'
+             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+             DO  j = 1, radius_classes
+                WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
+             ENDDO
+
+             DO  k = 1, dissipation_classes
+                DO  i = 1, radius_classes
+                   DO  j = 1, radius_classes
+                      IF ( hkernel(i,j) == 0.0 )  THEN
+                         hwratio(i,j) = 9999999.9
+                      ELSE
+                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
+                      ENDIF
+                   ENDDO
+                ENDDO
+
+                PRINT*, '*** epsilon = ', epsclass(k)
+                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+                DO  j = 1, radius_classes
+!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
+                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
+                ENDDO
+             ENDDO
+
+          ENDIF
+
+          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
+
+          ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
+
+       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
+       THEN
+!
+!--       Initial settings for Hall- and Wang-Kernel
+!--       To be done: move here parts from turbsd, fallg, ecoll, etc.
+       ENDIF
+
+    END SUBROUTINE init_kernels
+
+
+!------------------------------------------------------------------------------!
+! Calculation of collision kernels during each timestep and for each grid box
+!------------------------------------------------------------------------------!
+    SUBROUTINE recalculate_kernel( i1, j1, k1 )
 
        USE arrays_3d
@@ -94 +224 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, i1, j, j1, k1
-
-       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
-
-       ip = i1
-       jp = j1
-       kp = k1
-
-       pstart = prt_start_index(kp,jp,ip)
-       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
-       psum   = prt_count(kp,jp,ip)
-
-       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
-
-       IF ( wang_kernel )  THEN
-
-          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
-
-          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
-          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
-
-          IF ( epsilon <= 0.001 )  THEN
-
-             CALL fallg
-             CALL effic
-
-             DO  j = pstart, pend
-                DO  i =  pstart, pend
-                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                        particles(i)%radius * 100.0 )**2 &
-                                    * ec(i,j) * ABS( winf(j) - winf(i) )
-                ENDDO
+       INTEGER ::  i, i1, j, j1, k1, pend, pstart
+
+
+       pstart = prt_start_index(k1,j1,i1)
+       pend   = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
+       radius_classes = prt_count(k1,j1,i1)
+
+       ALLOCATE( ec(1:radius_classes,1:radius_classes), &
+                 radclass(1:radius_classes), winf(1:radius_classes) )
+
+!
+!--    Store particle radii on the radclass array. Collision routines
+!--    expect radii to be in cm.
+       radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
+
+       epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
+       urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
+
+       IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
+!
+!--       Call routines to calculate efficiencies for the Wang kernel
+          ALLOCATE( gck(1:radius_classes,1:radius_classes), &
+                    ecf(1:radius_classes,1:radius_classes) )
+
+          CALL turbsd
+          CALL turb_enhance_eff
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
              ENDDO
-
-          ELSE
-
-             CALL turbsd
-             CALL turb_enhance_eff
-             CALL effic
-
-             DO  j = pstart, pend, 1
-                DO  i =  pstart, pend, 1
-                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
-                ENDDO
+          ENDDO
+
+          DEALLOCATE( gck, ecf )
+
+       ELSE
+!
+!--       Call routines to calculate efficiencies for the Hall kernel
+          CALL fallg
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                ckernel(pstart+i-1,pstart+j-1,1) = pi *                       &
+                                          ( radclass(j) + radclass(i) )**2    &
+                                          * ec(i,j) * ABS( winf(j) - winf(i) )
              ENDDO
-
-          ENDIF
-
-          DEALLOCATE(gck, ecf)
-
-       ELSE
-
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
-          CALL fallg
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
-          CALL effic
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
-
-          DO  j = pstart, pend
-             DO  i =  pstart, pend
-                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                     particles(i)%radius * 100.0 )**2 &
-                                 * ec(i,j) * ABS( winf(j) - winf(i) )
-             ENDDO
           ENDDO
 
        ENDIF
 
-       DEALLOCATE( ec, winf )
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
-
-    END SUBROUTINE colker
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of w, g and gck
+       ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
+
+       DEALLOCATE( ec, radclass, winf )
+
+    END SUBROUTINE recalculate_kernel
+
+
+!------------------------------------------------------------------------------!
+! Calculation of gck
+! This is from Aayala 2008b, page 37ff.
+! Necessary input parameters: water density, radii of droplets, air density,
+! air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
+! gravitational acceleration  --> to be replaced by PALM parameters
 !------------------------------------------------------------------------------!
     SUBROUTINE turbsd
-! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
-! of droplets, air density, air viscosity, turbulent dissipation rate,
-! taylor microscale reynolds number, gravitational acceleration
 
        USE constants
@@ -186 +299 @@
        IMPLICIT NONE
 
-       REAL :: Relamda, &
-               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
-               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
-               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
-               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
-
-       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
-
-       REAL, DIMENSION(pstart:pend) :: St, tau
+       INTEGER ::  i, j
 
        LOGICAL, SAVE ::  first = .TRUE.
 
-       INTEGER :: i, j
-
-!
-!--   Initial assignment of constants
+       REAL ::  ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, &
+                e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp,   &
+                sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, &
+                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
+
+       REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
+
+       REAL, DIMENSION(1:radius_classes) ::  st, tau
+
+
+!
+!--    Initial assignment of constants
        IF ( first )  THEN
 
           first = .FALSE.
-          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
-          airdens   = 1.2250   !air density in mg/cm**3
-          waterdens = 1000.0   !water density in mg/cm**3
-          gravity   = 980.6650 !in cm/s**2
+          airvisc   = 0.1818     ! dynamic viscosity in mg/cm*s
+          airdens   = 1.2250     ! air density in mg/cm**3
+          waterdens = 1000.0     ! water density in mg/cm**3
+          gravity   = 980.6650   ! in cm/s**2
           anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
 
        ENDIF   
 
-       Relamda = urms**2*sqrt(15.0/epsilon/anu)
-
-       Tl = urms**2/epsilon       !in s
-       Lf = 0.5 * (urms**3)/epsilon !in cm
-
-       tauk = (anu/epsilon)**0.5       !in s
-       eta  = (anu**3/epsilon)**0.25  !in cm
-       vk   = eta/tauk                   !in cm/s
-
-       ao = (11.+7.*Relamda)/(205.+Relamda)
-
-       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
-
-       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
-
-       CALL fallg !gives winf in cm/s
-
-       DO i = pstart, pend
-          tau(i) = winf(i)/gravity    !in s
-          St(i)  = tau(i)/tauk
+       lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
+       lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
+       tl        = urms**2 / epsilon                       ! in s
+       lf        = 0.5 * urms**3 / epsilon                 ! in cm
+       tauk      = SQRT( anu / epsilon )                   ! in s
+       eta       = ( anu**3 / epsilon )**0.25              ! in cm
+       vk        = eta / tauk                              ! in cm/s
+
+       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
+       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
+
+       CALL fallg    ! gives winf in cm/s
+
+       DO  i = 1, radius_classes
+          tau(i) = winf(i) / gravity    ! in s
+          st(i)  = tau(i) / tauk
        ENDDO
 
-!***** TO CALCULATE wr ********************************
-!from Aayala 2008b, page 38f
-
-       z  = tt/Tl
-       be = sqrt(2.0)*lambda/Lf
-
-       bbb = sqrt(1.0-2.0*be**2)
-       d1 = (1.+bbb)/2.0/bbb
-       e1 = Lf*(1.0+bbb)/2.0   !in cm
-       d2 = (1.0-bbb)/2.0/bbb
-       e2 = Lf*(1.0-bbb)/2.0   !in cm
-
-       ccc = sqrt(1.0-2.0*z**2)
-       b1 = (1.+ccc)/2./ccc
-       c1 = Tl*(1.+ccc)/2.   !in s
-       b2 = (1.-ccc)/2./ccc
-       c2 = Tl*(1.-ccc)/2.   !in s
-
-       DO i = pstart, pend
-          v1 = winf(i)         !in cm/s
-          t1 = tau(i)         !in s
-
-          DO j = pstart,i
-             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
-             v2  = winf(j)         !in cm/s
-             t2  = tau(j)         !in s
-
-             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
-                     - b1*d2*PHI(c1,e2,v1,t1) &
-                     - b2*d1*PHI(c2,e1,v1,t1) &
-                     + b2*d2*PHI(c2,e2,v1,t1)
-             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
-
-             vrms1xy= sqrt(v1xysq)             !in cm/s
-
-             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
-                     - b1*d2*PHI(c1,e2,v2,t2) &
-                     - b2*d1*PHI(c2,e1,v2,t2) &
-                     + b2*d2*PHI(c2,e2,v2,t2)
-             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
-
-              vrms2xy= sqrt(v2xysq)             !in cm/s
-
-             IF(winf(i).ge.winf(j)) THEN
+!
+!--    Calculate wr (from Aayala 2008b, page 38f)
+       z   = tt / tl
+       be  = SQRT( 2.0 ) * lambda / lf
+       bbb = SQRT( 1.0 - 2.0 * be**2 )
+       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
+       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
+       d2  = ( 1.0 - bbb ) * 0.5 / bbb
+       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
+       ccc = SQRT( 1.0 - 2.0 * z**2 )
+       b1  = ( 1.0 + ccc ) * 0.5 / ccc
+       c1  = tl * ( 1.0 + ccc ) * 0.5   ! in s
+       b2  = ( 1.0 - ccc ) * 0.5 / ccc
+       c2  = tl * ( 1.0 - ccc ) * 0.5   ! in s
+
+       DO  i = 1, radius_classes
+
+          v1 = winf(i)        ! in cm/s
+          t1 = tau(i)         ! in s
+
+          DO  j = 1, i
+             rrp = radclass(i) + radclass(j)               ! radius in cm
+             v2  = winf(j)                                 ! in cm/s
+             t2  = tau(j)                                  ! in s
+
+             v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
+                     - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
+             v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
+             vrms1xy = SQRT( v1xysq )                       ! in cm/s
+
+             v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
+                     - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
+             v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
+             vrms2xy = SQRT( v2xysq )                       ! in cm/s
+
+             IF ( winf(i) >= winf(j) )  THEN
                 v1 = winf(i)
                 t1 = tau(i)
@@ -290 +391 @@
              ENDIF
 
-             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
-                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
-                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
-                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
-             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
-             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
-
-             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
-
-              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
-
-              wrgrav2=pi/8.*(winf(j)-winf(i))**2
-
-              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
-
-
-!***** TO CALCULATE gr ********************************
-
-             IF(St(j).gt.St(i)) THEN
-                SSt = St(j)
+             v1v2xy   =  b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
+                         b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
+                         b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
+                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
+             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
+             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
+             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
+             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
+             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
+             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
+
+!
+!--          Calculate gr
+             IF ( st(j) > st(i) )  THEN
+                sst = st(j)
              ELSE
-                SSt = St(i)
+                sst = st(i)
              ENDIF
 
-             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
-                  -4.2942*SSt**2 + 5.3406*SSt
-
-             IF(XX.le.0.0) XX = 0.0
-
-             YY = 0.1886*exp(20.306/Relamda)
-
-             c1_gr = XX/(gravity/(vk/tauk))**YY
-
-             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
-             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
-             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
-
-             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
-             IF (grFIN.lt.1.0) grFIN = 1.0
-
-             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
+             xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + &
+                   5.3406 * sst
+             IF ( xx < 0.0 )  xx = 0.0
+             yy = 0.1886 * EXP( 20.306 / lambda_re )
+
+             c1_gr  =  xx / ( gravity / vk * tauk )**yy
+
+             ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
+             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
+             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
+
+             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 )
+             IF ( grfin < 1.0 )  grfin = 1.0
+
+             gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin           ! in cm**3/s
              gck(j,i) = gck(i,j)
 
@@ -336 +430 @@
        ENDDO
 
-    END SUBROUTINE TurbSD
-
-!------------------------------------------------------------------------------!
-! PHI as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION PHI(a,b,vsett,tau0)
+    END SUBROUTINE turbsd
+
+
+!------------------------------------------------------------------------------!
+! phi as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION phi( a, b, vsett, tau0 )
 
        IMPLICIT NONE
 
-       REAL :: a, aa1, b, vsett, tau0
-
-       aa1 = 1./tau0 + 1./a + vsett/b
-
-       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
-
-    END FUNCTION PHI
-
-!------------------------------------------------------------------------------!
-! ZETA as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
+       REAL ::  a, aa1, b, tau0, vsett
+
+       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
+       phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
+
+    END FUNCTION phi
+
+
+!------------------------------------------------------------------------------!
+! zeta as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
 
        IMPLICIT NONE
 
-       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
-
-        aa1 = vsett2/b - 1./tau2 - 1./a
-        aa2 = vsett1/b + 1./tau1 + 1./a
-        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
-        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
-        aa5 = vsett2/b + 1./tau2 + 1./a
-        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
-        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
-             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
-             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
-                                                           * 1./2./b/aa3             ! in s**2
-
-    END FUNCTION ZHI
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of terminal velocity winf
-!------------------------------------------------------------------------------!
-   SUBROUTINE fallg
+       REAL ::  a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2
+
+       aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
+       aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
+       aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
+       aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2
+       aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
+       aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
+       zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
+           + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6&
+           + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )&
+           * 0.5 / b / aa3      ! in s**2
+
+    END FUNCTION zhi
+
+
+!------------------------------------------------------------------------------!
+! Calculation of terminal velocity winf
+!------------------------------------------------------------------------------!
+    SUBROUTINE fallg
 
        USE constants
@@ -385 +481 @@
        USE arrays_3d
 
-      IMPLICIT NONE
-
-      INTEGER :: i, j
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
-
-      REAL :: bond, x, xrey, y
-
-      REAL, DIMENSION(1:7), SAVE  :: b
-      REAL, DIMENSION(1:6), SAVE  :: c
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
+       IMPLICIT NONE
+
+       INTEGER ::  i, j
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
+                      t0, xlamb
+
+       REAL ::  bond, x, xrey, y
+
+       REAL, DIMENSION(1:7), SAVE  ::  b
+       REAL, DIMENSION(1:6), SAVE  ::  c
+
+!
+!--    Initial assignment of constants
+       IF ( first )  THEN
+
+          first = .FALSE.
+          b = (/  -0.318657E1,  0.992696E0, -0.153193E-2, -0.987059E-3, &
+                 -0.578878E-3, 0.855176E-4, -0.327815E-5 /)
+          c = (/  -0.500015E1,  0.523778E1,  -0.204914E1,   0.475294E0, &
+                 -0.542819E-1, 0.238449E-2 /)
+
+          eta   = 1.818E-4          ! in poise = g/(cm s)
+          xlamb = 6.62E-6           ! in cm
+          rhow  = 1.0               ! in g/cm**3
+          rhoa  = 1.225E-3          ! in g/cm**3
+          grav  = 980.665           ! in cm/s**2
+          cunh  = 1.257 * xlamb     ! in cm
+          t0    = 273.15            ! in K
+          sigma = 76.1 - 0.155 * ( 293.15 - t0 )              ! in N/m = g/s**2
+          stok  = 2.0  * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
+          stb   = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
+                                                                ! in 1/cm**3
+          phy   = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
+          py    = phy**( 1.0 / 6.0 )
+
+       ENDIF
+
+       DO  j = 1, radius_classes
+
+          IF ( radclass(j) <= 1.0E-3 ) THEN
+
+             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
+
+          ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
+
+             x = LOG( stb * radclass(j)**3 )
+             y = 0.0
+
+             DO  i = 1, 7
+                y = y + b(i) * x**(i-1)
+             ENDDO
+
+             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
+             winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
+
+          ELSEIF ( radclass(j) > 5.35E-2 )  THEN
+
+             IF ( radclass(j) > 0.35 )  THEN
+                bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
+             ELSE
+                bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
+             ENDIF
+
+             x = LOG( 16.0 * bond * py / 3.0 )
+             y = 0.0
+
+             DO  i = 1, 6
+                y = y + c(i) * x**(i-1)
+             ENDDO
+
+             xrey = py * EXP( y )
+
+             IF ( radclass(j) > 0.35 )  THEN
+                winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
+             ELSE
+                winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE fallg
+
+
+!------------------------------------------------------------------------------!
+! Calculation of collision efficencies for the Hall kernel
+!------------------------------------------------------------------------------!
+    SUBROUTINE effic
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE particle_attributes
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, iq, ir, j, k, kk
+
+       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       REAL ::  ek, particle_radius, pp, qq, rq
+
+       REAL, DIMENSION(1:21), SAVE ::  rat
+       REAL, DIMENSION(1:15), SAVE ::  r0
+       REAL, DIMENSION(1:15,1:21), SAVE ::  ecoll
+
+!
+!--    Initial assignment of constants
+       IF ( first )  THEN
 
          first = .FALSE.
-         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
-             -0.578878e-3,0.855176e-4,-0.327815e-5/)
-         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
-             -0.542819e-1, 0.238449e-2/)
-
-         eta   = 1.818e-4     !in poise = g/(cm s)
-         xlamb = 6.62e-6    !in cm
-
-         rhow = 1.0       !in g/cm**3
-         rhoa = 1.225e-3    !in g/cm**3
-
-         grav = 980.665     !in cm/s**2
-         cunh = 1.257 * xlamb !in cm
-         t0   = 273.15        !in K
-         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
-         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
-         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
-         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
-         py   = phy**(1.0/6.0)
-
-      ENDIF
-
-!particle radius has to be in cm
-      DO j = pstart, pend
-
-         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
-
-            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
-
-         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
-
-            x = log(stb*(particles(j)%radius*100.0)**3)
-            y = 0.0
-
-            DO i = 1, 7
-               y = y + b(i) * (x**(i-1))
-            ENDDO
-
-            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
-            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
-
-         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
-
-            IF (particles(j)%radius*100.0.gt.0.35)  THEN
-               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
-            ELSE
-               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
-            ENDIF
-
-            x = log(16.0*bond*py/3.0)
-            y = 0.0
-
-            DO i = 1, 6
-               y = y + c(i) * (x**(i-1))
-            ENDDO
-
-            xrey = py*exp(y)
-
-            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
-              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
-            ELSE
-               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
-            ENDIF
-
-         ENDIF
-      ENDDO
-      RETURN
-   END SUBROUTINE fallg
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision efficencies
-!------------------------------------------------------------------------------!
-
-   SUBROUTINE effic
-
-      USE arrays_3d
-      USE constants
-      USE cloud_parameters
-      USE particle_attributes
-
-!collision efficiencies of hall kernel
-      IMPLICIT NONE
-
-      INTEGER :: i, ir, iq, j, k, kk
-
-      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL ::  ek, particle_radius, pp, qq, rq
-
-      REAL, DIMENSION(1:21), SAVE :: rat
-      REAL, DIMENSION(1:15), SAVE :: r0
-      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
-                 300./)
-         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
-                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
-
-         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
-                        0.001,0.001,0.001,0.001,0.001,0.001/)
-         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
-                        0.050,0.200,0.500,0.770,0.870,0.970/)
-         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
-                        0.430,0.580,0.790,0.930,0.960,1.000/)
-         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
-                        0.640,0.750,0.910,0.970,0.980,1.000/)
-         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
-                        0.770,0.840,0.950,0.970,1.000,1.000/)
-         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
-                        0.840,0.880,0.950,1.000,1.000,1.000/)
-         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
-                        0.870,0.900,0.950,1.000,1.000,1.000/)
-         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
-                        0.890,0.920,1.000,1.000,1.000,1.000/)
-         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
-                        0.900,0.940,1.000,1.000,1.000,1.000/)
-         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
-                        0.920,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
-                        0.930,0.970,1.000,1.000,1.000,1.000/)
-         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
-                        0.950,1.000,1.000,1.000,1.000,1.000/)
-         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
-                        1.000,1.020,1.020,1.020,1.020,1.020/)
-         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
-                        1.030,1.040,1.040,1.040,1.040,1.040/)
-         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
-                        1.700,2.300,2.300,2.300,2.300,2.300/)
-         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
-                        3.000,4.000,4.000,4.000,4.000,4.000/)
-      ENDIF
-
-!
-!--   Calculate the radius class index of particles with respect to array r
-      ALLOCATE( ira(pstart:pend) )
-      DO  j = pstart, pend
-         particle_radius = particles(j)%radius * 1.0E6
-         DO  k = 1, 15
-            IF ( particle_radius < r0(k) )  THEN
-               ira(j) = k
-               EXIT
-            ENDIF
-         ENDDO
-         IF ( particle_radius >= r0(15) )  ira(j) = 16
-      ENDDO
-
-!
-!--   Two-dimensional linear interpolation of the collision efficiency.
-!--   Radius has to be in µm
-      DO  j = pstart, pend
-         DO  i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius / particles(j)%radius
-
-!            DO kk = 2, 21
-!               IF ( rq <= rat(kk) )  THEN
-!                  iq = kk
-!                  EXIT
-!               ENDIF
-!            ENDDO
-
-            iq = INT( rq * 20 ) + 1
-            iq = MAX(iq , 2)
-
-            IF ( ir < 16 )  THEN
-
-               IF ( ir >= 2 )  THEN
-                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
-                       ( r0(ir) - r0(ir-1) )
-                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
-                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
-                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
-                            + pp * qq * ecoll(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
-               ENDIF
-
-            ELSE
-               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
-               ec(j,i) = MIN( ek, 1.0 )
-            ENDIF
-
-            ec(i,j) = ec(j,i)
-            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
-
-         ENDDO
-      ENDDO
-
-      DEALLOCATE( ira )
-
-   END SUBROUTINE effic
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of enhancement factor collision efficencies
-!------------------------------------------------------------------------------!
-   SUBROUTINE turb_enhance_eff
+         r0  = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
+                  70.0, 100.0, 150.0, 200.0, 300.0 /)
+         rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
+                  0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
+                  1.00 /)
+
+         ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
+                        0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
+         ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
+                        0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
+         ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
+                        0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
+         ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
+                        0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
+         ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
+                        0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
+         ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
+                        0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
+         ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
+                        0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
+         ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
+                        0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
+                        0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
+                        0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
+                        0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
+                        0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
+                        0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
+                        0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
+                        0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
+         ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
+                        0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
+         ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
+                        0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
+         ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
+                        1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
+       ENDIF
+
+!
+!--    Calculate the radius class index of particles with respect to array r
+       ALLOCATE( ira(1:radius_classes) )
+       DO  j = 1, radius_classes
+          particle_radius = radclass(j) * 1.0E4   ! in microm
+          DO  k = 1, 15
+             IF ( particle_radius < r0(k) )  THEN
+                ira(j) = k
+                EXIT
+             ENDIF
+          ENDDO
+          IF ( particle_radius >= r0(15) )  ira(j) = 16
+       ENDDO
+
+!
+!--    Two-dimensional linear interpolation of the collision efficiency.
+!--    Radius has to be in µm
+       DO  j = 1, radius_classes
+          DO  i = 1, j
+
+             ir = ira(j)
+             rq = radclass(i) / radclass(j)
+             iq = INT( rq * 20 ) + 1
+             iq = MAX( iq , 2)
+
+             IF ( ir < 16 )  THEN
+                IF ( ir >= 2 )  THEN
+                   pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
+                        ( r0(ir) - r0(ir-1) )
+                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)  &
+                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)          &
+                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)          &
+                             + pp * qq * ecoll(ir,iq)
+                ELSE
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
+                ENDIF
+             ELSE
+                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
+                ec(j,i) = MIN( ek, 1.0 )
+             ENDIF 
+
+             ec(i,j) = ec(j,i)
+             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
+
+          ENDDO
+       ENDDO
+
+       DEALLOCATE( ira )
+
+    END SUBROUTINE effic
+
+
+!------------------------------------------------------------------------------!
+! Calculation of enhancement factor for collision efficencies due to turbulence
+!------------------------------------------------------------------------------!
+    SUBROUTINE turb_enhance_eff
 
        USE constants
@@ -627 +713 @@
        USE arrays_3d
 
-      IMPLICIT NONE
-
-      INTEGER :: i, ik, ir, iq, j, k, kk
-
-      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
-
-      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL, DIMENSION(1:11), SAVE :: rat
-      REAL, DIMENSION(1:7), SAVE  :: r0
-      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-
-         r0  = (/10., 20., 30.,40., 50., 60.,100./)
-         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
-
-! 100 cm^2/s^3
-         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
-         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
-         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
-         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
-         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
-         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
-         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
-         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
-         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
-         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
-         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
-
-! 400 cm^2/s^3
-         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
-         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
-         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
-         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
-         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
-         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
-         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
-         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
-         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
-         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
-         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
-
-      ENDIF
-
-!
-!--   Calculate the radius class index of particles with respect to array r
-      ALLOCATE( ira(pstart:pend) )
-
-      DO  j = pstart, pend
-         particle_radius = particles(j)%radius * 1.0E6
-         DO  k = 1, 7
-            IF ( particle_radius < r0(k) )  THEN
-               ira(j) = k
-               EXIT
-            ENDIF
-         ENDDO
-         IF ( particle_radius >= r0(7) )  ira(j) = 8
-      ENDDO
-
-! two-dimensional linear interpolation of the collision efficiency
-      DO j =  pstart, pend
-         DO i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius/particles(j)%radius
-
-            DO kk = 2, 11
-               IF ( rq <= rat(kk) )  THEN
-                  iq = kk
-                  EXIT
-               ENDIF
-            ENDDO
-
-! 0  cm2/s3
-            y1 = 1.0
-! 100 cm2/s3, 400 cm2/s3
-            IF (ir.lt.8) THEN
-               IF (ir.ge.2) THEN
-                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
-                            pp*qq*ecoll_100(ir,iq)
-                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
-                            pp*qq*ecoll_400(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
-                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
-               ENDIF
-            ELSE
-            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
-            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
-         ENDIF
-! linear interpolation
-! dissipation rate in cm2/s3
-         x1 = 0.0
-         x2 = 100.0
-         x3 = 400.0
-
-         IF (epsilon.le.100.) THEN
-            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
-                    +  (epsilon-0.)/(100.-0.) * y2
-         ELSE IF(epsilon.le.600.)THEN
-            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
-                    +  (epsilon-100.)/(400.-100.) * y3
-
-         ELSE
-            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
-                    +  (600.-100.)/(400.-100.) * y3
-         ENDIF
-
-         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
-
-         ecf(i,j)=ecf(j,i)
-      ENDDO
-   ENDDO
-
-   RETURN
-   END SUBROUTINE turb_enhance_eff
+       IMPLICIT NONE
+
+       INTEGER :: i, ik, iq, ir, j, k, kk
+
+       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+       REAL ::  particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       REAL, DIMENSION(1:11), SAVE ::  rat
+       REAL, DIMENSION(1:7), SAVE  ::  r0
+       REAL, DIMENSION(1:7,1:11), SAVE ::  ecoll_100, ecoll_400
+
+!
+!--    Initial assignment of constants
+       IF ( first )  THEN
+
+          first = .FALSE.
+
+          r0  = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
+          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
+!
+!--       In 100 cm^2/s^3
+          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
+          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
+          ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
+          ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
+          ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
+          ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
+          ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
+          ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
+          ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
+          ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
+          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
+!
+!--       In 400 cm^2/s^3
+          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
+          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
+          ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
+          ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
+          ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
+          ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
+          ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
+          ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
+          ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
+          ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
+          ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
+
+       ENDIF
+
+!
+!--    Calculate the radius class index of particles with respect to array r0
+       ALLOCATE( ira(1:radius_classes) )
+
+       DO  j = 1, radius_classes
+          particle_radius = radclass(j) * 1.0E4    ! in microm
+          DO  k = 1, 7
+             IF ( particle_radius < r0(k) )  THEN
+                ira(j) = k
+                EXIT
+             ENDIF
+          ENDDO
+          IF ( particle_radius >= r0(7) )  ira(j) = 8
+       ENDDO
+
+!
+!--    Two-dimensional linear interpolation of the collision efficiencies
+       DO  j =  1, radius_classes
+          DO  i = 1, j
+
+             ir = ira(j)
+             rq = radclass(i) / radclass(j)
+
+             DO  kk = 2, 11
+                IF ( rq <= rat(kk) )  THEN
+                   iq = kk
+                   EXIT
+                ENDIF
+             ENDDO
+
+             y1 = 1.0      ! 0  cm2/s3
+!
+!--          100 cm2/s3, 400 cm2/s3
+             IF ( ir < 8 )  THEN
+                IF ( ir >= 2 )  THEN
+                   pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
+                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
+                                qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
+                                pp * qq         * ecoll_100(ir,iq)
+                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
+                                pp * ( 1.0-qq ) * ecoll_400(ir,iq-1)   +  &
+                                qq * ( 1.0-pp ) * ecoll_400(ir-1,iq)   +  &
+                                pp * qq         * ecoll_400(ir,iq)
+                ELSE
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
+                   y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
+                ENDIF
+             ELSE
+                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
+                y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
+             ENDIF
+!
+!--          Linear interpolation of dissipation rate in cm2/s3
+             IF ( epsilon <= 100.0 )  THEN
+                ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
+                         + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
+             ELSEIF ( epsilon <= 600.0 )  THEN
+                ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
+                         + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
+             ELSE
+                ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
+                         + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
+             ENDIF
+
+             IF ( ecf(j,i) < 1.0 )  ecf(j,i) = 1.0
+
+             ecf(i,j) = ecf(j,i)
+
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE turb_enhance_eff
 
  END MODULE lpm_collision_kernels_mod

palm/trunk/SOURCE/modules.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! +dissipation_classes, radius_classes, use_kernel_tables,
+! particle feature color renamed class
 !
 ! Former revisions:
@@ -1195 +1196 @@
     INTEGER ::  mpi_particle_type
 #endif
-    INTEGER ::  ibc_par_lr, ibc_par_ns, ibc_par_b, ibc_par_t,                  &
-                iran_part = -1234567, maximum_number_of_particles = 1000,      &
+    INTEGER ::  dissipation_classes = 10, ibc_par_lr, ibc_par_ns, ibc_par_b,   &
+                ibc_par_t, iran_part = -1234567,                               &
+                maximum_number_of_particles = 1000,                            &
                 maximum_number_of_tailpoints = 100,                            &
                 maximum_number_of_tails = 0,                                   &
@@ -1203 +1205 @@
                 number_of_initial_tails = 0, offset_ocean_nzt = 0,             &
                 offset_ocean_nzt_m1 = 0, particles_per_point = 1,              &
-                particle_file_count = 0, skip_particles_for_tail = 100,        &
-                sort_count = 0, total_number_of_particles,                     &
-                total_number_of_tails = 0
+                particle_file_count = 0, radius_classes = 20,                  &
+                skip_particles_for_tail = 100, sort_count = 0,                 &
+                total_number_of_particles, total_number_of_tails = 0
 
     INTEGER, PARAMETER ::  max_number_of_particle_groups = 10
@@ -1215 +1217 @@
                 particle_advection = .FALSE., random_start_position = .FALSE., &
                 read_particles_from_restartfile = .TRUE.,                      &
-                uniform_particles = .TRUE., use_particle_tails = .FALSE.,      &
-                use_sgs_for_particles = .FALSE.,                               &
+                uniform_particles = .TRUE., use_kernel_tables = .FALSE.,       &
+                use_particle_tails = .FALSE., use_sgs_for_particles = .FALSE., &
                 wang_kernel = .FALSE., write_particle_statistics = .FALSE.
 
@@ -1247 +1249 @@
                    origin_z, radius, rvar1, rvar2, rvar3, speed_x, speed_y, &
                    speed_z, weight_factor, x, y, z
-       INTEGER ::  color, group, tailpoints, tail_id
+       INTEGER ::  class, group, tailpoints, tail_id
     END TYPE particle_type
 

palm/trunk/SOURCE/package_parin.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! +dissipation_classes, radius_classes in parpar
 !
 ! Former revisions:
@@ -75 +75 @@
 
     NAMELIST /particles_par/      bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,    &
-                                  collision_kernel,                            &
-                                  density_ratio, radius, dt_dopts,             &
+                                  collision_kernel, density_ratio,             &
+                                  dissipation_classes, dt_dopts,               &
                                   dt_min_part, dt_prel, dt_sort_particles,     &
                                   dt_write_particle_data,                      &
@@ -88 +88 @@
                                   particle_advection_start,                    &
                                   particle_maximum_age, pdx, pdy, pdz, psb,    &
-                                  psl, psn, psr, pss, pst,                     &
-                                  random_start_position,                       &
+                                  psl, psn, psr, pss, pst, radius,             &
+                                  radius_classes, random_start_position,       &
                                   read_particles_from_restartfile,             &
                                   skip_particles_for_tail,                     &

palm/trunk/SOURCE/set_particle_attributes.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! particle feature color renamed class
 !
 ! Former revisions:
@@ -121 +121 @@
 !
 !--       Number of available colors is defined in init_dvrp
-          particles(n)%color = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
+          particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
 
        ENDDO
@@ -192 +192 @@
 !
 !--       Number of available colors is defined in init_dvrp
-          particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
+          particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
 
        ENDDO
@@ -213 +213 @@
 !
 !--       Number of available colors is defined in init_dvrp
-          particles(n)%color = 1 + height * ( dvrp_colortable_entries_prt - 1 )
+          particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )
 
        ENDDO